mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
Finished basic build of multibox region. Still bug-filled, needs fixing.
This commit is contained in:
parent
5f88186863
commit
3154f19ac0
1 changed files with 97 additions and 27 deletions
|
@ -18,42 +18,42 @@
|
||||||
#include <algorithm>
|
#include <algorithm>
|
||||||
|
|
||||||
#include "region_generator.hh"
|
#include "region_generator.hh"
|
||||||
|
|
||||||
#include "point_file_reader.hh"
|
#include "point_file_reader.hh"
|
||||||
|
|
||||||
|
|
||||||
typedef std::vector<int> col;
|
typedef std::vector<unsigned> col;
|
||||||
typedef std::vector<col> slice;
|
typedef std::vector<col> slice;
|
||||||
typedef std::vector<slice> grid;
|
typedef std::vector<slice> grid;
|
||||||
|
|
||||||
typedef struct {
|
typedef struct {
|
||||||
int x;
|
unsigned x;
|
||||||
int y;
|
unsigned y;
|
||||||
int z;
|
unsigned z;
|
||||||
} coord;
|
} coord;
|
||||||
|
|
||||||
typedef std::vector<coord> region;
|
typedef std::vector<coord> region;
|
||||||
|
|
||||||
class region_multibox_plugin : public region_generator_plugin{
|
class region_multibox_plugin : public region_generator_plugin{
|
||||||
private:
|
private:
|
||||||
int res, reflevel;
|
unsigned res;
|
||||||
|
double vfac_;
|
||||||
grid refgrid;
|
grid refgrid;
|
||||||
|
|
||||||
region find_equal(grid ingrid, int level)
|
region where(unsigned level)
|
||||||
{
|
{
|
||||||
slice curslice;
|
slice curslice;
|
||||||
col curcol;
|
col curcol;
|
||||||
region equal_region;
|
region equal_region;
|
||||||
coord cp;
|
coord cp;
|
||||||
for(int i=0;i<ingrid.size();i++)
|
for(unsigned i=0;i<refgrid.size();i++)
|
||||||
{
|
{
|
||||||
cp.x = i;
|
cp.x = i;
|
||||||
curslice = ingrid[i];
|
curslice = refgrid[i];
|
||||||
for(int j=0;j<curslice.size();j++)
|
for(unsigned j=0;j<curslice.size();j++)
|
||||||
{
|
{
|
||||||
cp.y = j;
|
cp.y = j;
|
||||||
curcol = curslice[j];
|
curcol = curslice[j];
|
||||||
for(int k=0;k<curcol.size();k++)
|
for(unsigned k=0;k<curcol.size();k++)
|
||||||
{
|
{
|
||||||
cp.z = k;
|
cp.z = k;
|
||||||
if(curcol[k] == level)
|
if(curcol[k] == level)
|
||||||
|
@ -66,32 +66,32 @@ private:
|
||||||
return equal_region;
|
return equal_region;
|
||||||
}
|
}
|
||||||
//This function takes the grid, a vector of doubles containing the particle
|
//This function takes the grid, a vector of doubles containing the particle
|
||||||
//coordinates, and sets each particle-containing gridcell to maxlevel.
|
//coordinates, and sets each particle-containing gridcell to levelmax_.
|
||||||
void deposit_particles(grid ingrid, std::vector<double> pp, int maxlevel, int res)
|
void deposit_particles(std::vector<double> pp, unsigned res)
|
||||||
{
|
{
|
||||||
int i,j,k;
|
unsigned i,j,k;
|
||||||
std::vector<double>::iterator cp = pp.begin();
|
std::vector<double>::iterator cp = pp.begin();
|
||||||
while(cp != pp.end())
|
while(cp != pp.end())
|
||||||
{
|
{
|
||||||
i = *(cp);
|
i = (*(cp)+0.5)*res;
|
||||||
cp++;
|
cp++;
|
||||||
j = *(cp);
|
j = (*(cp)+0.5)*res;
|
||||||
cp++;
|
cp++;
|
||||||
k = *(cp);
|
k = (*(cp)+0.5)*res;
|
||||||
cp++;
|
cp++;
|
||||||
ingrid[int((i+0.5)*res)][int((j+0.5)*res)][int((k+0.5)*res)] = maxlevel;
|
refgrid[i][j][k] = levelmax_;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
//This function takes the grid, which has been already had particles
|
//This function takes the grid, which has been already had particles
|
||||||
//deposited onto it, and set to the maximum refinement level. It then
|
//deposited onto it, and set to the maximum refinement level. It then
|
||||||
//fills the remaining refinement levels
|
//fills the remaining refinement levels
|
||||||
void build_refgrid(grid ingrid, int maxlevel)
|
void build_refgrid()
|
||||||
{
|
{
|
||||||
region curregion;
|
region curregion;
|
||||||
for(int curlevel=maxlevel; curlevel>2; curlevel--)
|
for(unsigned curlevel=levelmax_; curlevel>(levelmin_+1); curlevel--)
|
||||||
{
|
{
|
||||||
curregion = find_equal(ingrid, curlevel);
|
curregion = where(curlevel);
|
||||||
for(region::iterator cp= curregion.begin(); cp != curregion.end(); ++cp)
|
for(region::iterator cp= curregion.begin(); cp != curregion.end(); ++cp)
|
||||||
{
|
{
|
||||||
for(int i=-1; i<2; i++)
|
for(int i=-1; i<2; i++)
|
||||||
|
@ -100,9 +100,9 @@ private:
|
||||||
{
|
{
|
||||||
for(int k=-1; k<2; k++)
|
for(int k=-1; k<2; k++)
|
||||||
{
|
{
|
||||||
if(ingrid[cp->x+i][cp->y+j][cp->z+k] == 0)
|
if(refgrid[cp->x+i][cp->y+j][cp->z+k] == levelmin_)
|
||||||
{
|
{
|
||||||
ingrid[cp->x+i][cp->y+j][cp->z+k] = curlevel-1;
|
refgrid[cp->x+i][cp->y+j][cp->z+k] = curlevel-1;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -115,33 +115,103 @@ public:
|
||||||
explicit region_multibox_plugin( config_file& cf )
|
explicit region_multibox_plugin( config_file& cf )
|
||||||
: region_generator_plugin( cf )
|
: region_generator_plugin( cf )
|
||||||
{
|
{
|
||||||
|
//check parameters
|
||||||
|
if ( !cf.containsKey("setup", "region_point_file"))
|
||||||
|
{
|
||||||
|
LOGERR("Missing parameter \'region_point_file\' needed for region=multibox");
|
||||||
|
throw std::runtime_error("Missing parameter for region=multibox");
|
||||||
|
}
|
||||||
|
vfac_ = cf.getValue<double>("cosmology","vfact");
|
||||||
|
if (levelmin_ >= levelmax_)
|
||||||
|
{
|
||||||
|
LOGERR("Why are you specifying a region if you aren't refining?");
|
||||||
|
throw std::runtime_error("region=multibox needs levelmax>levelmin");
|
||||||
|
}
|
||||||
std::vector<double> pp;
|
std::vector<double> pp;
|
||||||
res = 128;
|
point_reader pfr;
|
||||||
|
pfr.read_points_from_file(cf.getValue<std::string>("setup", "region_point_file"), vfac_, pp);
|
||||||
|
res = 1<<(levelmin_-1);
|
||||||
|
LOGINFO("Multibox Resolution: %d", res);
|
||||||
//Initialize the grid with zeros, the base level
|
//Initialize the grid with zeros, the base level
|
||||||
refgrid = grid(res,slice(res,col(res,0)));
|
refgrid = grid(res,slice(res,col(res,levelmin_)));
|
||||||
deposit_particles(refgrid, pp, reflevel, res);
|
//Build the highest refinement region
|
||||||
build_refgrid(refgrid, reflevel);
|
deposit_particles(pp, res);
|
||||||
|
LOGINFO("Deposited Multigrid Particles");
|
||||||
|
//Build the intermediate refinement regions
|
||||||
|
build_refgrid();
|
||||||
|
LOGINFO("Built Multigrid");
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
void get_AABB( double *left, double *right, unsigned level )
|
void get_AABB( double *left, double *right, unsigned level )
|
||||||
{
|
{
|
||||||
|
left[0] = left[1] = left[2] = 0.5;
|
||||||
|
right[0] = right[1] = right[2] = -0.5;
|
||||||
|
if( level <= levelmin_ )
|
||||||
|
{
|
||||||
|
left[0] = left[1] = left[2] = 0.0;
|
||||||
|
right[0] = right[1] = right[2] = 1.0;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
region myres = where(level);
|
||||||
|
std::vector<coord>::iterator cp = myres.begin();
|
||||||
|
while (cp != myres.end())
|
||||||
|
{
|
||||||
|
//Check left side
|
||||||
|
if ((cp->x-0.5)/res-0.5 < left[0])
|
||||||
|
left[0] = (cp->x-0.5)/res-0.5;
|
||||||
|
if ((cp->y-0.5)/res-0.5 < left[1])
|
||||||
|
left[1] = (cp->y-0.5)/res-0.5;
|
||||||
|
if ((cp->z-0.5)/res-0.5 < left[2])
|
||||||
|
left[2] = (cp->z-0.5)/res-0.5;
|
||||||
|
//Check right side
|
||||||
|
if ((cp->x+0.5)/res-0.5 > right[0])
|
||||||
|
right[0] = (cp->x+0.5)/res-0.5;
|
||||||
|
if ((cp->y+0.5)/res-0.5 > right[1])
|
||||||
|
right[1] = (cp->y+0.5)/res-0.5;
|
||||||
|
if ((cp->z+0.5)/res-0.5 > right[2])
|
||||||
|
right[2] = (cp->z+0.5)/res-0.5;
|
||||||
|
cp++;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void update_AABB( double *left, double *right )
|
void update_AABB( double *left, double *right )
|
||||||
{
|
{
|
||||||
|
//no need for this I think?
|
||||||
}
|
}
|
||||||
|
|
||||||
bool query_point( double *x, int level )
|
bool query_point( double *x, int level )
|
||||||
{
|
{
|
||||||
|
if(fabs(x[0]) > 0.5 || fabs(x[1]) > 0.5 || fabs(x[2]) > 0.5)
|
||||||
|
{
|
||||||
|
printf("Outside point: %3.2e %3.2e %3.2e\n", x[0], x[1], x[2]);
|
||||||
|
return 0;
|
||||||
|
}
|
||||||
|
return (level == int(refgrid[(x[0]+0.5)*res][(x[1]+0.5)*res][(x[2]+0.5)*res]));
|
||||||
}
|
}
|
||||||
|
|
||||||
bool is_grid_dim_forced( size_t* ndims )
|
bool is_grid_dim_forced( size_t* ndims )
|
||||||
{
|
{
|
||||||
|
return false; //is this true?
|
||||||
}
|
}
|
||||||
|
|
||||||
void get_center( double *xc )
|
void get_center( double *xc )
|
||||||
{
|
{
|
||||||
|
region hires = where(levelmax_);
|
||||||
|
std::vector<coord>::iterator cp = hires.begin();
|
||||||
|
int n=0;
|
||||||
|
while(cp != hires.end())
|
||||||
|
{
|
||||||
|
xc[0] += cp->x-0.5;
|
||||||
|
xc[1] += cp->y-0.5;
|
||||||
|
xc[2] += cp->z-0.5;
|
||||||
|
n++;
|
||||||
|
cp++;
|
||||||
|
}
|
||||||
|
xc[0] /= n;
|
||||||
|
xc[1] /= n;
|
||||||
|
xc[2] /= n;
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
void get_center_unshifted( double *xc )
|
void get_center_unshifted( double *xc )
|
||||||
|
|
Loading…
Reference in a new issue