Commit 06fdb909 authored by frey_m's avatar frey_m
Browse files

co-moving mesh for EllipticDomain only

parent 960ea5a6
......@@ -129,8 +129,7 @@ void PartBunch::computeSelfFields(int binNumber) {
if(fs_m->hasValidSolver()) {
/// Mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// Scatter charge onto space charge grid.
this->Q *= this->dt;
......@@ -411,12 +410,14 @@ void PartBunch::computeSelfFields(int binNumber) {
}
void PartBunch::resizeMesh() {
if (fs_m->getFieldSolverType() != "SAAMG") {
return;
}
double xmin = fs_m->solver_m->getXRangeMin();
double xmax = fs_m->solver_m->getXRangeMax();
double ymin = fs_m->solver_m->getYRangeMin();
double ymax = fs_m->solver_m->getYRangeMax();
double zmin = fs_m->solver_m->getZRangeMin();
double zmax = fs_m->solver_m->getZRangeMax();
if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
ymin > rmin_m[1] || ymax < rmax_m[1]) {
......@@ -437,14 +438,14 @@ void PartBunch::resizeMesh() {
boundp();
get_bounds(rmin_m, rmax_m);
}
Vector_t mymin = Vector_t(xmin, ymin , zmin);
Vector_t mymax = Vector_t(xmax, ymax , zmax);
for (int i = 0; i < 3; i++)
hr_m[i] = (mymax[i] - mymin[i])/nr_m[i];
Vector_t origin = Vector_t(0.0, 0.0, 0.0);
// update the mesh origin and mesh spacing hr_m
fs_m->solver_m->resizeMesh(origin, hr_m, this);
getMesh().set_meshSpacing(&(hr_m[0]));
getMesh().set_origin(mymin);
getMesh().set_origin(origin);
rho_m.initialize(getMesh(),
getFieldLayout(),
......@@ -467,8 +468,7 @@ void PartBunch::computeSelfFields() {
if(fs_m->hasValidSolver()) {
//mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
//scatter charges onto grid
this->Q *= this->dt;
......@@ -652,8 +652,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......@@ -844,8 +843,7 @@ void PartBunch::computeSelfFields_cycl(int bin) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......
......@@ -242,6 +242,8 @@ public:
virtual void set_meshEnlargement(double dh);
double getMeshEnlargement() const;
void gatherLoadBalanceStatistics();
size_t getLoadBalance(int p) const;
......@@ -394,6 +396,8 @@ public:
// virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
virtual FieldLayout_t &getFieldLayout() = 0;
virtual void resizeMesh() { };
/*
* Wrapped member functions of IpplParticleBase
*/
......@@ -507,6 +511,10 @@ public:
// get 2nd order momentum matrix
FMatrix<double, 2 * Dim, 2 * Dim> getSigmaMatrix();
const Vector_t& getLowerBound() const;
const Vector_t& getUpperBound() const;
private:
// save particles in case of one core
std::unique_ptr<Inform> pmsg_m;
......
......@@ -1194,6 +1194,12 @@ void PartBunchBase<T, Dim>::set_meshEnlargement(double dh) {
}
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getMeshEnlargement() const {
return dh_m;
}
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::gatherLoadBalanceStatistics() {
......@@ -2277,4 +2283,15 @@ FMatrix<double, 2 * Dim, 2 * Dim> PartBunchBase<T, Dim>::getSigmaMatrix() {
return sigmaMatrix;
}
template <class T, unsigned Dim>
const Vector_t& PartBunchBase<T, Dim>::getLowerBound() const {
return rmin_m;
}
template <class T, unsigned Dim>
const Vector_t& PartBunchBase<T, Dim>::getUpperBound() const {
return rmax_m;
}
#endif
\ No newline at end of file
......@@ -159,6 +159,24 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
}
}
void EllipticDomain::resizeMesh(Vector_t& origin, Vector_t& hr, PartBunchBase<double, 3>* bunch) {
Vector_t rmin = bunch->getLowerBound();
Vector_t rmax = bunch->getUpperBound();
Vector_t mymax = Vector_t(0.0, 0.0, 0.0);
double dh = bunch->getMeshEnlargement();
// apply bounding box increment, i.e., "BBOXINCR" input argument
double zmin = std::signbit(rmin[2]) ? rmin[2] * (1.0 + dh) : rmin[2] * (1.0 - dh);
double zmax = std::signbit(rmax[2]) ? rmax[2] * (1.0 - dh) : rmax[2] * (1.0 + dh);
origin = Vector_t(-semiMajor_m, -semiMinor_m, zmin);
mymax = Vector_t( semiMajor_m, semiMinor_m, zmax);
for (int i = 0; i < 3; ++i)
hr[i] = (mymax[i] - origin[i]) / nr[i];
}
void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W,
double &E, double &S, double &N,
double &F, double &B, double &C,
......
......@@ -97,6 +97,9 @@ public:
bool hasGeometryChanged() { return hasGeometryChanged_m; }
void resizeMesh(Vector_t& origin, Vector_t& hr, PartBunchBase<double, 3>* bunch);
private:
/// Map from a single coordinate (x or y) to a list of intersection values with
......
......@@ -27,6 +27,7 @@
#include <string>
#include "Algorithms/PBunchDefs.h"
#include "Algorithms/Quaternion.h"
#include "Algorithms/PartBunchBase.h"
/// enumeration corresponding to different interpolation methods at the boundary
enum {
......@@ -126,6 +127,22 @@ public:
virtual bool hasGeometryChanged() = 0;
virtual ~IrregularDomain() {};
// FIXME The function body should be implemented in derived classes.
virtual void resizeMesh(Vector_t& origin, Vector_t& hr, PartBunchBase<double, 3>* /*bunch*/) {
double xmin = getXRangeMin();
double xmax = getXRangeMax();
double ymin = getYRangeMin();
double ymax = getYRangeMax();
double zmin = getZRangeMin();
double zmax = getZRangeMax();
origin = Vector_t(xmin, ymin , zmin);
Vector_t mymax = Vector_t(xmax, ymax , zmax);
for (int i = 0; i < 3; i++)
hr[i] = (mymax[i] - origin[i]) / nr[i];
};
protected:
// a irregular domain is always defined on a grid
......
......@@ -58,6 +58,8 @@ public:
virtual void test(PartBunchBase<double, 3> *bunch) = 0 ;
virtual ~PoissonSolver(){};
virtual void resizeMesh(Vector_t& /*origin*/, Vector_t& /*hr*/, PartBunchBase<double, 3>* /*bunch*/) { };
};
inline Inform &operator<<(Inform &os, const PoissonSolver &/*fs*/) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment