Skip to content

Commit

Permalink
Re-structure boundary structure in fieldreader for particles.
Browse files Browse the repository at this point in the history 
This also affected the fluxfunction (and made that it didn't compile),
because it had some hardcoded particleParameters in the field
loading code. Now that has been factored out.

Also, the fluxfunction tool supports equatorial plane files now
(although the support is quite hacky)
  • Loading branch information
Urs Ganse committed Jul 20, 2016
1 parent 42947a8 commit 5b54291
Show file tree
Hide file tree
Showing 6 changed files with 368 additions and 242 deletions.
256 changes: 124 additions & 132 deletions particles/boundaries.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,166 +2,158 @@

#include "particles.h"

class Boundary
struct Boundary
{
public:
// Handle a particle's boundary behaviour.
// returns "true" if the particle is still part of the simulation
// afterwards, or "false" if it is to be removed.
virtual bool handleParticle(Particle& p) = 0;

// Handle cell coordinate in the spatial dimension this boundary
// object cares about (for example: wrap in a periodic direction)
virtual int cellCoordinate(int c) = 0;

Boundary(int _dimension) : dimension(_dimension) {};
virtual void setExtent(double _min, double _max, int _cells) {
min=_min;
max=_max;
cells=_cells;
}
// Handle a particle's boundary behaviour.
// returns "true" if the particle is still part of the simulation
// afterwards, or "false" if it is to be removed.
virtual bool handleParticle(Particle& p) = 0;

// Handle cell coordinate in the spatial dimension this boundary
// object cares about (for example: wrap in a periodic direction)
virtual int cellCoordinate(int c) = 0;

Boundary(int _dimension) : dimension(_dimension) {};
virtual void setExtent(double _min, double _max, int _cells) {
min=_min;
max=_max;
cells=_cells;
}

virtual ~Boundary(){};
virtual ~Boundary(){};

protected:
// Which spatial dimension to handle
int dimension;
// Which spatial dimension to handle
int dimension;

// Minimum and maximum spatial extents in this dimension
double min, max;
// Minimum and maximum spatial extents in this dimension
double min, max;

// Number of cells in this dimension
int cells;
// Number of cells in this dimension
int cells;
};

// Boundary for a spatial dimension that is only 1 cell thick (pseudo-periodic)
class CompactSpatialDimension : public Boundary
struct CompactSpatialDimension : public Boundary
{
public:
virtual bool handleParticle(Particle& p) {
// This boundary does not affect particles
return true;
}
virtual int cellCoordinate(int c) {
// Actual cell coordinates in this direction are
// always mapped to 0.
return 0;
}
CompactSpatialDimension(int _dimension) : Boundary(_dimension){};
virtual bool handleParticle(Particle& p) {
// This boundary does not affect particles
return true;
}
virtual int cellCoordinate(int c) {
// Actual cell coordinates in this direction are
// always mapped to 0.
return 0;
}
CompactSpatialDimension(int _dimension) : Boundary(_dimension){};
};

// Open boundary, which removes particles if they fly out
class OpenBoundary : public Boundary
struct OpenBoundary : public Boundary
{
public:
virtual bool handleParticle(Particle& p) {

// Delete particles that are outside our boundaries.
if(p.x[dimension] <= min || p.x[dimension] >= max) {
return false;
} else {
// Leave all others be.
return true;
}
}
virtual bool handleParticle(Particle& p) {

virtual int cellCoordinate(int c) {
// Cell coordinates are clamped
// TODO: Should this print warnings?
if(c < 0) {
return 0;
} else if(c >= cells) {
return cells-1;
} else {
return c;
}
// Delete particles that are outside our boundaries.
if(p.x[dimension] <= min || p.x[dimension] >= max) {
return false;
} else {
// Leave all others be.
return true;
}
}

virtual void setExtent(double _min, double _max, int _cells) {
double dx = (_max-_min)/((double)_cells);
min=_min+2*dx; // 2 Cells border.
max=_max-2*dx;
cells=_cells;
virtual int cellCoordinate(int c) {
// Cell coordinates are clamped
// TODO: Should this print warnings?
if(c < 0) {
return 0;
} else if(c >= cells) {
return cells-1;
} else {
return c;
}
OpenBoundary(int _dimension) : Boundary(_dimension){};
}

virtual void setExtent(double _min, double _max, int _cells) {
double dx = (_max-_min)/((double)_cells);
min=_min+2*dx; // 2 Cells border.
max=_max-2*dx;
cells=_cells;
}
OpenBoundary(int _dimension) : Boundary(_dimension){};
};

class ReflectBoundary : public Boundary
struct ReflectBoundary : public Boundary
{
public:
virtual bool handleParticle(Particle& p) {
// Particles outside of bounds get their velocities flipped
if(p.x[dimension] <= min || p.x[dimension] >= max) {
p.v *= flip_v;
}
return true;
virtual bool handleParticle(Particle& p) {
// Particles outside of bounds get their velocities flipped
if(p.x[dimension] <= min || p.x[dimension] >= max) {
p.v *= flip_v;
}
return true;
}

virtual int cellCoordinate(int c) {
// Cell coordinates are clamped
// TODO: Should this print warnings?
if(c < 0) {
return 0;
} else if(c >= cells) {
return cells-1;
} else {
return c;
}
}

// Constructor
ReflectBoundary(int _dimension) : Boundary(_dimension) {
double flip[3] = {1.,1.,1.};
flip[dimension] = -1.;
flip_v.load(flip);
}
virtual void setExtent(double _min, double _max, int _cells) {
double dx = (_max-_min)/((double)_cells);
min=_min+2*dx; // 2 Cells border.
max=_max-2*dx;
cells=_cells;
virtual int cellCoordinate(int c) {
// Cell coordinates are clamped
// TODO: Should this print warnings?
if(c < 0) {
return 0;
} else if(c >= cells) {
return cells-1;
} else {
return c;
}

protected:
// Vector to multiply with in order to flip velocity
// vectors for our dimension
Vec3d flip_v;
}

// Constructor
ReflectBoundary(int _dimension) : Boundary(_dimension) {
double flip[3] = {1.,1.,1.};
flip[dimension] = -1.;
flip_v.load(flip);
}
virtual void setExtent(double _min, double _max, int _cells) {
double dx = (_max-_min)/((double)_cells);
min=_min+2*dx; // 2 Cells border.
max=_max-2*dx;
cells=_cells;
}

// Vector to multiply with in order to flip velocity
// vectors for our dimension
Vec3d flip_v;

};

class PeriodicBoundary : public Boundary
struct PeriodicBoundary : public Boundary
{
public:
virtual bool handleParticle(Particle& p) {
if(p.x[dimension] < min) {
p.x += offset_p;
} else if(p.x[dimension] >= max) {
p.x -= offset_p;
}
return true;
virtual bool handleParticle(Particle& p) {
if(p.x[dimension] < min) {
p.x += offset_p;
} else if(p.x[dimension] >= max) {
p.x -= offset_p;
}

virtual int cellCoordinate(int c) {
return c % cells;
}

// Constructor
PeriodicBoundary(int _dimension) : Boundary(_dimension) {
}
virtual void setExtent(double _min, double _max, int _cells) {
min=_min;
max=_max;
cells=_cells;

double offset[3] = {0.,0.,0.};
offset[dimension] = max-min;
offset_p.load(offset);
}

protected:
// Vector to offset particle positions that leave through
// one boundary with, to come out the other end
Vec3d offset_p;
return true;
}

virtual int cellCoordinate(int c) {
return c % cells;
}

// Constructor
PeriodicBoundary(int _dimension) : Boundary(_dimension) {
}
virtual void setExtent(double _min, double _max, int _cells) {
min=_min;
max=_max;
cells=_cells;

double offset[3] = {0.,0.,0.};
offset[dimension] = max-min;
offset_p.load(offset);
}

// Vector to offset particle positions that leave through
// one boundary with, to come out the other end
Vec3d offset_p;

};

Expand Down
36 changes: 20 additions & 16 deletions particles/field.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,36 +12,39 @@ struct Field
// Time at which this field is "valid"
double time;

// Coordinate boundaries
double min[3];
double max[3];

// Mesh spacing
double dx[3];

// Mesh cells
int cells[3];

bool periodic[3];
// Information about spatial dimensions (like extent, boundaries, etc)
Boundary* dimension[3];

// The actual field data
std::vector<double> data;

// Constructor (primarily here to make sure boundaries are properly initialized as zero)
Field() {
for(int i=0; i<3; i++) {
dimension[i] = nullptr;
}
}

double* getCellRef(int x, int y, int z) {

if(cells[2] == 1) {
return &(data[4*(y*cells[0]+x)]);
if(dimension[2]->cells == 1) {
// Equatorial plane
return &(data[4*(y*dimension[0]->cells+x)]);
} else {
return &(data[4*(z*cells[0]*cells[1] + y*cells[0] + x)]);
// General 3d case
return &(data[4*(z*dimension[0]->cells*dimension[1]->cells + y*dimension[0]->cells + x)]);
}
}

Vec3d getCell(int x, int y, int z) {

// Map these cell coordinates using the boundaries
x = ParticleParameters::boundary_behaviour_x->cellCoordinate(x);
y = ParticleParameters::boundary_behaviour_y->cellCoordinate(y);
z = ParticleParameters::boundary_behaviour_z->cellCoordinate(z);
x = dimension[0]->cellCoordinate(x);
y = dimension[1]->cellCoordinate(y);
z = dimension[2]->cellCoordinate(z);

double* cell = getCellRef(x,y,z);
Vec3d retval;
Expand All @@ -52,6 +55,7 @@ struct Field
// Round-Brace indexing: indexing by physical location, with interpolation
Vec3d operator()(Vec3d v) {
Vec3d vmin,vdx;
double min[3] = { min[0] = dimension[0]->min, dimension[1]->min, dimension[2]->min};
vmin.load(min);
vdx.load(dx);

Expand All @@ -63,7 +67,7 @@ struct Field
truncate_to_int(v).store(index);
(v-truncate(v)).store(fract);

if(cells[2] <= 1) {
if(dimension[2]->cells <= 1) {
// Equatorial plane
Vec3d interp[4];
interp[0] = getCell(index[0],index[1],index[2]);
Expand All @@ -73,7 +77,7 @@ struct Field

return fract[0]*(fract[1]*interp[3]+(1.-fract[1])*interp[1])
+ (1.-fract[0])*(fract[1]*interp[2]+(1.-fract[1])*interp[0]);
} else if (cells[1] <= 1) {
} else if (dimension[1]->cells <= 1) {
// Polar plane
Vec3d interp[4];

Expand Down
Loading

0 comments on commit 5b54291

Please sign in to comment.