C++ EBIndexSpace::define方法代码示例

//-----------------------------------------------------------------------
void
EBAMRGodunovGeomFactory::
createRegularCellsFromParms(Box& a_coarsestDomain,
                            RealVect& a_coarsestDx,
                            Vector<int>& a_refRatios)
{
  // Retrieve grid data.
  Box finestDomain;
  int maxGridSize, maxCoarsenings;
  EBAMRGodunovGeomFactory::createRefinedDomainsFromParms(a_coarsestDomain,
                                                         finestDomain,
                                                         maxGridSize,
                                                         maxCoarsenings,
                                                         a_refRatios);

  // Retrieve basic geometric information.
  RealVect domainExtents;
  RealVect origin, finestDx;
  EBAMRGodunovGeomFactory::getBasicGeometry(domainExtents, origin,
                                            a_coarsestDx, finestDx,
                                            a_refRatios);

  // now define the regular geometry.
  AllRegularService regServ;
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
  ebisPtr->define(finestDomain, origin, finestDx[0], regServ, maxGridSize,
                  maxCoarsenings);
}

void do_geo(DisjointBoxLayout&              a_dbl,
            EBISLayout       &              a_ebisl,
            Box              &              a_domain,
            Real             &              a_dx )
{

  //define grids
  int len = 32;
  int mid = len/2;
  Real xdom = 1.0;
  a_dx = xdom/len;
  const IntVect hi = (len-1)*IntVect::Unit;
  a_domain=  Box(IntVect::Zero, hi);
  Vector<Box> boxes(4);
  Vector<int> procs(4);
  boxes[0] = Box(IntVect(D_DECL(0  ,  0,   0)), IntVect(D_DECL(mid-1, mid-1, len-1)));
  boxes[1] = Box(IntVect(D_DECL(0  ,mid,   0)), IntVect(D_DECL(mid-1, len-1, len-1)));
  boxes[2] = Box(IntVect(D_DECL(mid,  0,   0)), IntVect(D_DECL(len-1, mid-1, len-1)));
  boxes[3] = Box(IntVect(D_DECL(mid,mid,   0)), IntVect(D_DECL(len-1, len-1, len-1)));

  LoadBalance(procs, boxes);
  /**
  int loProc = 0;
  int hiProc = numProc() -1;
  procs[0] = loProc;
  procs[1] = hiProc;
  procs[2] = loProc;
  procs[3] = hiProc;
  **/

  a_dbl = DisjointBoxLayout(boxes, procs);

  //define geometry
  RealVect rampNormal = RealVect::Zero;
  rampNormal[0] = -0.5;
  rampNormal[1] = 0.866025404;

  Real rampAlpha = -0.0625;

  RealVect rampPoint = RealVect::Zero;
  rampPoint[0] = rampAlpha / rampNormal[0];

  bool inside = true;

  PlaneIF ramp(rampNormal,rampPoint,inside);

  RealVect vectDx = RealVect::Unit;
  vectDx *= a_dx;

  GeometryShop mygeom( ramp, 0, vectDx );
  RealVect origin = RealVect::Zero;

  EBIndexSpace *ebisPtr = Chombo_EBIS::instance();
  ebisPtr->define(a_domain, origin, a_dx, mygeom);

  //fill layout
  ebisPtr->fillEBISLayout(a_ebisl, a_dbl, a_domain, 4 );
}

int
EBRestart::makeGeometry(EBIndexSpace&                a_ebIndexSpace,
                        const EBRestart::GeomParams& a_geomParams)
{
  int eekflag = 0;
  SphereIF outside(a_geomParams.radius, a_geomParams.center, false);

  RealVect vectDx = RealVect::Unit;
  vectDx *= a_geomParams.dx;

  GeometryShop workshop(outside,0,vectDx);

  workshop.m_phase = 0;

  ProblemDomain pdomain( a_geomParams.domain );
  a_ebIndexSpace.define(pdomain,
                        a_geomParams.origin,
                        a_geomParams.dx,
                        workshop);

  return eekflag;
}

//-----------------------------------------------------------------------
void
EBAMRGodunovGeomFactory::
createSphericalCavityFromParms(Box& a_coarsestDomain,
                               RealVect& a_coarsestDx,
                               Vector<int>& a_refRatios)
{
  // Retrieve grid data.
  Box finestDomain;
  int maxGridSize, maxCoarsenings;
  EBAMRGodunovGeomFactory::createRefinedDomainsFromParms(a_coarsestDomain,
                                                         finestDomain,
                                                         maxGridSize,
                                                         maxCoarsenings,
                                                         a_refRatios);

  // Retrieve basic geometric information.
  RealVect domainExtents;
  RealVect origin, finestDx;
  EBAMRGodunovGeomFactory::getBasicGeometry(domainExtents, origin,
                                            a_coarsestDx, finestDx,
                                            a_refRatios);

  // now define the cavity.
  ParmParse pp;
  vector<Real> x0(SpaceDim);
  pp.getarr("sphere_center",x0, 0, SpaceDim);
  RealVect sphereCenter;
  for (int i = 0; i < SpaceDim; ++i)
    sphereCenter[i] = x0[i];
  Real sphereRadius;
  pp.get("sphere_radius", sphereRadius);
  SphereIF sphereIF(sphereRadius, sphereCenter, true);
  int verbosity = 0;
  GeometryShop workshop(sphereIF,verbosity,finestDx);
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
  ebisPtr->define(finestDomain, origin, finestDx[0], workshop, maxGridSize,
                  maxCoarsenings);
}

void do_geo( Box &domain, Real &dx )
{
  int N = 10;
  Real xdom = 1.0;

  dx = xdom/N;
  const IntVect hi = (N-1)*IntVect::Unit;
  domain=  Box(IntVect::Zero, hi);

  int updir=1;
  int indv=0;
  Real start = 0.5;
  Real slope = tan( M_PI/6.0 );

  RealVect normal = RealVect::Zero;
  normal[updir] = 1.0;
  normal[indv] = -slope;

  RealVect point = RealVect::Zero;
  point[updir] = -slope*start;

  bool normalInside = true;

  PlaneIF myramp(normal,point,normalInside);

  RealVect vectDx = RealVect::Unit;
  vectDx *= dx;

  GeometryShop mygeom( myramp, 0, vectDx );

  RealVect origin = RealVect::Zero;

  EBIndexSpace *ebisPtr = Chombo_EBIS::instance();

  ebisPtr->define(domain, origin, dx, mygeom);

  return;
}

int makeGeometry(Box& a_domain)
{
  Real dx;
  RealVect origin;
  int eekflag =  0;
  //parse input file
  ParmParse pp;

  Vector<int> n_cell(SpaceDim);
  pp.getarr("n_cell",n_cell,0,SpaceDim);

  CH_assert(n_cell.size() == SpaceDim);
  IntVect lo = IntVect::Zero;
  IntVect hi;
  for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
      if (n_cell[ivec] <= 0)
        {
          pout() << " bogus number of cells input = " << n_cell[ivec];
          return(-1);
        }
      hi[ivec] = n_cell[ivec] - 1;
    }

  a_domain.setSmall(lo);
  a_domain.setBig(hi);

  Vector<Real> prob_lo(SpaceDim, 1.0);
  Real prob_hi;
  pp.getarr("prob_lo",prob_lo,0,SpaceDim);
  pp.get("prob_hi",prob_hi);
  dx = (prob_hi-prob_lo[0])/n_cell[0];
  RealVect dxVect = dx*RealVect::Unit;
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      origin[idir] = prob_lo[idir];
    }
  int verbosity = 0;
  int whichgeom;
  pp.get("which_geom",whichgeom);
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();

  if (whichgeom == 0)
    {
      //allregular
      pout() << "all regular geometry" << endl;
      AllRegularService regserv;
      ebisPtr->define(a_domain, origin, dx, regserv);
    }
  else if (whichgeom == 1)
    {
      pout() << "ramp geometry" << endl;
      int upDir;
      int indepVar;
      Real startPt;
      Real slope;
      pp.get("up_dir",upDir);
      pp.get("indep_var",indepVar);
      pp.get("start_pt", startPt);
      pp.get("ramp_slope", slope);

      RealVect normal = RealVect::Zero;
      normal[upDir] = 1.0;
      normal[indepVar] = -slope;

      RealVect point = RealVect::Zero;
      point[upDir] = -slope*startPt;

      bool normalInside = true;

      PlaneIF ramp(normal,point,normalInside);

      GeometryShop workshop(ramp,verbosity,dxVect);
      //this generates the new EBIS
      ebisPtr->define(a_domain, origin, dx, workshop);
    }
  else if (whichgeom == 5)
    {
      pout() << "sphere geometry" << endl;
      vector<Real> sphere_center(SpaceDim);
      pp.getarr("sphere_center",sphere_center, 0, SpaceDim);
      RealVect sphereCenter;
      for (int idir = 0; idir < SpaceDim; idir++)
        {
          sphereCenter[idir] = sphere_center[idir];
        }
      Real sphereRadius;
      pp.get("sphere_radius", sphereRadius);

      bool     insideRegular = false;
      SphereIF implicit(sphereRadius,sphereCenter,insideRegular);
      GeometryShop workshop(implicit,verbosity,dxVect);
      //this generates the new EBIS
      ebisPtr->define(a_domain, origin, dx, workshop);
    }
  else if (whichgeom == 13)
    {
      pout() << "rhodonea geometry" << endl;
      vector<Real> tmp(SpaceDim);
      pp.getarr("rhodonea_center", tmp, 0, SpaceDim);
//.........这里部分代码省略.........

BaseIF* makeGeometry(Box&      a_domain,
                     RealVect& a_origin,
                     Real&     a_dx)
{
  RealVect point;
  RealVect normal;
  bool     normalRegular;

  // parse input file
  ParmParse pp;

  Vector<int> n_cell(SpaceDim);
  pp.getarr("n_cell",n_cell,0,SpaceDim);

 CH_assert(n_cell.size() == SpaceDim);

  IntVect lo = IntVect::Zero;
  IntVect hi;

  for (int ivec = 0; ivec < SpaceDim; ivec++)
  {
    if (n_cell[ivec] <= 0)
    {
      pout() << "Bogus number of cells input = " << n_cell[ivec];
      exit(1);
    }

    hi[ivec] = n_cell[ivec] - 1;
  }

  a_domain.setSmall(lo);
  a_domain.setBig(hi);

  Vector<Real> prob_lo(SpaceDim,1.0);
  Real prob_hi;

  pp.getarr("prob_lo",prob_lo,0,SpaceDim);
  pp.get("prob_hi",prob_hi);

  a_dx = (prob_hi-prob_lo[0])/n_cell[0];

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    a_origin[idir] = prob_lo[idir];
  }

  // ParmParse doesn't get RealVects, so work-around with Vector<Real>
  Vector<Real> vectorCenter;
  pp.getarr("point",vectorCenter,0,SpaceDim);

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    point[idir] = vectorCenter[idir];
  }

  Vector<Real> vectorNormal;
  pp.getarr("normal",vectorNormal,0,SpaceDim);

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    normal[idir] = vectorNormal[idir];
  }

  // Parm Parse doesn't get bools, so work-around with int
  int intInsideRegular;
  pp.get("normalRegular",intInsideRegular);

  if (intInsideRegular != 0) normalRegular = true;
  if (intInsideRegular == 0) normalRegular = false;

  PlaneIF implicit(normal,point,normalRegular);

  RealVect vectDx = RealVect::Unit;
  vectDx *= a_dx;

  GeometryShop workshop(implicit,0,vectDx);

  // This generates the new EBIS
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
  ebisPtr->define(a_domain,a_origin,a_dx,workshop);

  return implicit.newImplicitFunction();
}

int makeGeometry(Box& a_domain,
                 Real& a_dx,
                 RealVect& a_origin,
                 int& a_upDir,
                 int& a_indepVar,
                 Real& a_startPt,
                 Real& a_slope)
{
    int eekflag =  0;
    //parse input file
    ParmParse pp;

    Vector<int> n_cell(SpaceDim);
    pp.getarr("n_cell",n_cell,0,SpaceDim);

    CH_assert(n_cell.size() == SpaceDim);
    IntVect lo = IntVect::Zero;
    IntVect hi;
    for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
        if (n_cell[ivec] <= 0)
        {
            pout() << " bogus number of cells input = " << n_cell[ivec];
            return(-1);
        }
        hi[ivec] = n_cell[ivec] - 1;
    }

    a_domain.setSmall(lo);
    a_domain.setBig(hi);

    Vector<Real> prob_lo(SpaceDim, 1.0);
    Real prob_hi;
    pp.getarr("prob_lo",prob_lo,0,SpaceDim);
    pp.get("prob_hi",prob_hi);
    a_dx = (prob_hi-prob_lo[0])/n_cell[0];
    for (int idir = 0; idir < SpaceDim; idir++)
    {
        a_origin[idir] = prob_lo[idir];
    }
    pp.get("up_dir",a_upDir);
    a_upDir = SpaceDim-1;
    pp.get("indep_var",a_indepVar);
    pp.get("start_pt", a_startPt);
    pp.get("ramp_slope", a_slope);

    RealVect normal = RealVect::Zero;
    normal[a_upDir] = 1.0;
    normal[a_indepVar] = -a_slope;

    RealVect point = RealVect::Zero;
    point[a_upDir] = -a_slope*a_startPt;

    bool normalInside = true;

    PlaneIF ramp(normal,point,normalInside);

    RealVect vectDx = RealVect::Unit;
    vectDx *= a_dx;

    GeometryShop workshop(ramp,0,vectDx);
    //this generates the new EBIS
    EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
    ebisPtr->define(a_domain, a_origin, a_dx, workshop);

    return eekflag;
}

//.........这里部分代码省略.........
  Real outerRadius;
  pp.get("outer_radius",outerRadius);

  Real outerThick;
  pp.get("outer_thick",outerThick);

  Real outerHeight;
  pp.get("outer_height",outerHeight);

  Real outerOffset = ((prob_hi - prob_lo[0]) - outerHeight) / 2.0 + prob_lo[0];

  Real innerRadius;
  pp.get("inner_radius",innerRadius);

  Real innerThick;
  pp.get("inner_thick",innerThick);

  Real innerOffset = 0.0;
  innerOffset += outerOffset;

  Real innerHeight;
  pp.get("inner_height",innerHeight);

  Real plateHeight;
  pp.get("plate_height",plateHeight);
  plateHeight += outerOffset;

  Real plateThick;
  pp.get("plate_thick",plateThick);

  int doHoles;
  pp.get("do_holes",doHoles);

  Real holeRadius;
  pp.get("hole_radius",holeRadius);

  Real holeSpace;
  pp.get("hole_space",holeSpace);

  int vaneNum;
  pp.get("vane_num",vaneNum);

  Real vaneThick;
  pp.get("vane_thick",vaneThick);

  RealVect vaneNorm;

  Vector<Real> vectVaneNorm;
  pp.getarr("vane_norm",vectVaneNorm,0,SpaceDim);

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    vaneNorm[idir] = vectVaneNorm[idir];
  }

  Real vaneOffset;
  pp.get("vane_offset",vaneOffset);

  Real vaneHeight = innerHeight - 2*vaneOffset;

  vaneOffset += outerOffset;

  // Make the outer chamber
  BaseIF* outerChamber = makeChamber(outerRadius,outerThick,
                                     outerOffset,outerHeight);

  // Make the inner chamber
  BaseIF* innerChamber = makeChamber(innerRadius,innerThick,
                                     innerOffset,innerHeight);

  // Make the inner plate with holes
  BaseIF* holyPlate = makePlate(plateHeight,plateThick,innerRadius,
                                doHoles,holeRadius,holeSpace);

  // Make the vanes
  BaseIF* vanes = makeVanes(vaneNum,vaneThick,vaneNorm,innerRadius,outerRadius,
                            vaneOffset,vaneHeight);

  // Union all the pieces together
  Vector<BaseIF*> pieces;

  pieces.push_back(outerChamber);
  pieces.push_back(innerChamber);
  pieces.push_back(holyPlate);
  pieces.push_back(vanes);

  UnionIF insideSwirl(pieces);
  ComplementIF outsideSwirl(insideSwirl,true);

  RealVect vectDx = RealVect::Unit;
  vectDx *= a_dx;

  GeometryShop workshop(outsideSwirl,0,vectDx);

  // This generates the new EBIS
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
  ebisPtr->define(a_domain,a_origin,a_dx,workshop);

  return outsideSwirl.newImplicitFunction();
}

int makeGeometry(Box& a_domain)
{
  int eekflag =  0;
  //parse input file
  ParmParse pp;
  RealVect origin = RealVect::Zero;
  Vector<int> n_cell(SpaceDim);
  pp.getarr("n_cell",n_cell,0,SpaceDim);
  Real probhi;
  pp.get("prob_hi", probhi);
  Real dx = probhi/n_cell[0];

 CH_assert(n_cell.size() == SpaceDim);
  IntVect lo = IntVect::Zero;
  IntVect hi;
  for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
      if (n_cell[ivec] <= 0)
        {
          pout() << " bogus number of cells input = " << n_cell[ivec];
          return(-1);
        }
      hi[ivec] = n_cell[ivec] - 1;
    }

  a_domain.setSmall(lo);
  a_domain.setBig(hi);

  int iverbose;
  pp.get("verbose", iverbose);
  bool verbose = (iverbose==1);

  int whichgeom;
  pp.get("which_geom",whichgeom);
  if (whichgeom == 0)
    {
      //allregular
      if (verbose)
        pout() << "all regular geometry" << endl;
      AllRegularService regserv;
      EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
      ebisPtr->define(a_domain, origin, dx, regserv, 2048);
    }
  else if (whichgeom == 1)
    {
      if (verbose)
        pout() << "ramp geometry" << endl;
      int upDir;
      int indepVar;
      Real startPt;
      Real slope;
      pp.get("up_dir",upDir);
      pp.get("indep_var",indepVar);
      pp.get("start_pt", startPt);
      pp.get("ramp_slope", slope);

      RealVect normal = RealVect::Zero;
      normal[upDir] = 1.0;
      normal[indepVar] = -slope;

      RealVect point = RealVect::Zero;
      point[upDir] = -slope*startPt;

      bool normalInside = true;

      PlaneIF ramp(normal,point,normalInside);

      RealVect vectDx = RealVect::Unit;
      vectDx *= dx;

      GeometryShop workshop(ramp,0,vectDx);
      //this generates the new EBIS
      EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
      ebisPtr->define(a_domain, origin, dx, workshop, 2048);
    }
  else if (whichgeom == 2)
    {
      if (verbose)
        pout() << "slab geometry" << endl;
      vector<int> slab_lo(SpaceDim);
      pp.getarr("slab_lo",slab_lo,0,SpaceDim);
      vector<int> slab_hi(SpaceDim);
      pp.getarr("slab_hi",slab_hi,0,SpaceDim);
      IntVect lo, hi;
      for (int idir = 0; idir < SpaceDim; idir++)
        {
          lo[idir] = slab_lo[idir];
          hi[idir] = slab_hi[idir];
        }
      Box coveredBox(lo,hi);
      SlabService slab(coveredBox);
      //this generates the new EBIS
      EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
      RealVect origin = RealVect::Zero;
      ebisPtr->define(a_domain, origin, dx, slab,2048);
    }
  else
    {
      //bogus which_geom
      pout() << " bogus which_geom input = " << whichgeom;
//.........这里部分代码省略.........

void
sphereGeometry(Box& a_coarsestDomain,
               Real& a_dx)
{
    ParmParse ppgodunov;
    //parse input file
    int max_level = 0;
    ppgodunov.get("max_level",max_level);

    int num_read_levels = Max(max_level,1);
    std::vector<int> refRatios; // (num_read_levels,1);
    // note this requires a refRatio to be defined for the
    // finest level (even though it will never be used)

    ppgodunov.getarr("ref_ratio",refRatios,0,num_read_levels+1);
    ParmParse pp;
    RealVect origin = RealVect::Zero;
    Vector<int> n_cell(SpaceDim);
    pp.getarr("n_cell",n_cell,0,SpaceDim);

    CH_assert(n_cell.size() == SpaceDim);
    IntVect lo = IntVect::Zero;
    IntVect hi;
    for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
        if (n_cell[ivec] <= 0)
        {
            pout() << " bogus number of cells input = " << n_cell[ivec];
            MayDay::Error();
        }
        hi[ivec] = n_cell[ivec] - 1;
    }

    a_coarsestDomain = Box(lo, hi);
    Box finestDomain = a_coarsestDomain;
    for (int ilev = 0; ilev < max_level; ilev++)
    {
        finestDomain.refine(refRatios[ilev]);
    }

    Real prob_hi;
    int numOpen = n_cell[0];
    pp.get("domain_length",prob_hi);
    a_dx = prob_hi/numOpen;
    Real fineDx = a_dx;
    int ebMaxCoarsen = 2;
    for (int ilev = 0; ilev < max_level; ilev++)
    {
        fineDx /= refRatios[ilev];
        ebMaxCoarsen += refRatios[ilev]/2;
    }
    int ebMaxSize;
    ppgodunov.get("max_grid_size", ebMaxSize);
    EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
    pout() << "sphere geometry" << endl;
    vector<Real> sphere_center(SpaceDim);
    pp.getarr("sphere_center",sphere_center, 0, SpaceDim);
    RealVect sphereCenter;
    for (int idir = 0; idir < SpaceDim; idir++)
    {
        sphereCenter[idir] = sphere_center[idir];
    }
    Real sphereRadius;
    pp.get("sphere_radius", sphereRadius);
    int inside_fluid;
    pp.get("sphere_fluid_inside", inside_fluid);
    bool insideFluid = (inside_fluid==1);
    SphereIF sphereIF(sphereRadius, sphereCenter, insideFluid);
    RealVect vectFineDx(D_DECL(fineDx,fineDx,fineDx));
    GeometryShop workshop(sphereIF,0,vectFineDx);
    //this generates the new EBIS
    ebisPtr->define(finestDomain, origin, fineDx, workshop, ebMaxSize, ebMaxCoarsen);
}

BaseIF* makeGeometry(Box&      a_domain,
                     RealVect& a_origin,
                     Real&     a_dx)
{
  RealVect center1;
  Real     radius1;

  RealVect center2;
  Real     radius2;

  bool     insideRegular;

  // parse input file
  ParmParse pp;

  Vector<int> n_cell(SpaceDim);
  pp.getarr("n_cell",n_cell,0,SpaceDim);

 CH_assert(n_cell.size() == SpaceDim);

  IntVect lo = IntVect::Zero;
  IntVect hi;

  for (int ivec = 0; ivec < SpaceDim; ivec++)
  {
    if (n_cell[ivec] <= 0)
    {
      pout() << "Bogus number of cells input = " << n_cell[ivec];
      exit(1);
    }

    hi[ivec] = n_cell[ivec] - 1;
  }

  a_domain.setSmall(lo);
  a_domain.setBig(hi);

  Vector<Real> prob_lo(SpaceDim,1.0);
  Real prob_hi;

  pp.getarr("prob_lo",prob_lo,0,SpaceDim);
  pp.get("prob_hi",prob_hi);

  a_dx = (prob_hi-prob_lo[0])/n_cell[0];

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    a_origin[idir] = prob_lo[idir];
  }

  // ParmParse doesn't get RealVects, so work-around with Vector<Real>
  Vector<Real> vectorCenter;
  pp.getarr("center1",vectorCenter,0,SpaceDim);

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    center1[idir] = vectorCenter[idir];
  }

  pp.get("radius1",radius1);

  pp.getarr("center2",vectorCenter,0,SpaceDim);

  for (int idir = 0; idir < SpaceDim; idir++)
  {
    center2[idir] = vectorCenter[idir];
  }

  pp.get("radius2",radius2);

  // Parm Parse doesn't get bools, so work-around with int
  int intInsideRegular;
  pp.get("insideRegular",intInsideRegular);

  if (intInsideRegular != 0) insideRegular = true;
  if (intInsideRegular == 0) insideRegular = false;

  // Stay inside until the end
  bool inside = true;

  // Sphere 1 is the given "radius1" and "center1"
  SphereIF sphere1(radius1,center1,inside);

  // Sphere 2 is the given "radius2" and "center2"
  SphereIF sphere2(radius2,center2,inside);

  // Take the union of the spheres
  UnionIF implicit(sphere1,sphere2);

  // Complement if necessary
  ComplementIF insideOut(implicit,!insideRegular);

  RealVect vectDx = RealVect::Unit;
  vectDx *= a_dx;

  GeometryShop workshop(insideOut,0,vectDx);

  // This generates the new EBIS
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
  ebisPtr->define(a_domain,a_origin,a_dx,workshop);
//.........这里部分代码省略.........

BaseIF* makeGeometry(Box&      a_domain,
                     RealVect& a_origin,
                     Real&     a_dx)
{
    RealVect center;
    RealVect radii;
    bool     insideRegular;

    // parse input file
    ParmParse pp;

    Vector<int> n_cell(SpaceDim);
    pp.getarr("n_cell",n_cell,0,SpaceDim);

    CH_assert(n_cell.size() == SpaceDim);

    IntVect lo = IntVect::Zero;
    IntVect hi;

    for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
        if (n_cell[ivec] <= 0)
        {
            pout() << "Bogus number of cells input = " << n_cell[ivec];
            exit(1);
        }

        hi[ivec] = n_cell[ivec] - 1;
    }

    a_domain.setSmall(lo);
    a_domain.setBig(hi);

    Vector<Real> prob_lo(SpaceDim,1.0);
    Real prob_hi;

    pp.getarr("prob_lo",prob_lo,0,SpaceDim);
    pp.get("prob_hi",prob_hi);

    a_dx = (prob_hi-prob_lo[0])/n_cell[0];

    for (int idir = 0; idir < SpaceDim; idir++)
    {
        a_origin[idir] = prob_lo[idir];
    }

    // ParmParse doesn't get RealVects, so work-around with Vector<Real>
    Vector<Real> vectorCenter;
    pp.getarr("center",vectorCenter,0,SpaceDim);

    for (int idir = 0; idir < SpaceDim; idir++)
    {
        center[idir] = vectorCenter[idir];
    }

    Vector<Real> vectorRadii;
    pp.getarr("radii",vectorRadii,0,SpaceDim);

    for (int idir = 0; idir < SpaceDim; idir++)
    {
        radii[idir] = vectorRadii[idir];
    }

    int halfTurns;
    pp.get("halfturns",halfTurns);

    // Parm Parse doesn't get bools, so work-around with int
    int intInsideRegular;
    pp.get("insideRegular",intInsideRegular);

    if (intInsideRegular != 0) insideRegular = true;
    if (intInsideRegular == 0) insideRegular = false;

    IntVect zero=IntVect::Zero;
    bool inside = true;
    EllipsoidIF ellipsoid(radii,center,inside);

    Real coef = halfTurns / 2.0;
    IntVect power(D_DECL(1,0,0));

    Vector<PolyTerm> term;
    term.resize(1);
    term[0].coef   = coef;
    term[0].powers = power;

    PolynomialIF rotate(term,inside);

    LatheIF implicit(ellipsoid,rotate,center,insideRegular);

    RealVect vectDx = RealVect::Unit;
    vectDx *= a_dx;

    GeometryShop workshop(implicit,0,vectDx);

    // This generates the new EBIS
    EBIndexSpace* ebisPtr = Chombo_EBIS::instance();
    ebisPtr->define(a_domain,a_origin,a_dx,workshop);

    return implicit.newImplicitFunction();
}

int makeGeometry(Box& a_domain)
{
  Real dx;
  RealVect origin;
  int eekflag =  0;
  //parse input file
  ParmParse pp;

  Vector<int> n_cell(SpaceDim);
  pp.getarr("n_cell",n_cell,0,SpaceDim);

  CH_assert(n_cell.size() == SpaceDim);
  IntVect lo = IntVect::Zero;
  IntVect hi;
  for (int ivec = 0; ivec < SpaceDim; ivec++)
    {
      if (n_cell[ivec] <= 0)
        {
          pout() << " bogus number of cells input = " << n_cell[ivec];
          return(-1);
        }
      hi[ivec] = n_cell[ivec] - 1;
    }

  a_domain.setSmall(lo);
  a_domain.setBig(hi);

  Vector<Real> prob_lo(SpaceDim, 1.0);
  Real prob_hi;
  pp.getarr("prob_lo",prob_lo,0,SpaceDim);
  pp.get("prob_hi",prob_hi);
  dx = (prob_hi-prob_lo[0])/n_cell[0];
  RealVect dxVect = dx*RealVect::Unit;
  for (int idir = 0; idir < SpaceDim; idir++)
    {
      origin[idir] = prob_lo[idir];
    }
  int verbosity = 0;
  int whichgeom;
  pp.get("which_geom",whichgeom);
  EBIndexSpace* ebisPtr = Chombo_EBIS::instance();

  if (whichgeom == 0)
    {
      //allregular
      pout() << "all regular geometry" << endl;
      AllRegularService regserv;
      ebisPtr->define(a_domain, origin, dx, regserv);
    }
  else if (whichgeom == 1)
    {
      pout() << "ramp geometry" << endl;
      int upDir;
      int indepVar;
      Real startPt;
      Real slope;
      pp.get("up_dir",upDir);
      pp.get("indep_var",indepVar);
      pp.get("start_pt", startPt);
      pp.get("ramp_slope", slope);

      RealVect normal = RealVect::Zero;
      normal[upDir] = 1.0;
      normal[indepVar] = -slope;

      RealVect point = RealVect::Zero;
      point[upDir] = -slope*startPt;

      bool normalInside = true;

      PlaneIF ramp(normal,point,normalInside);

      GeometryShop workshop(ramp,verbosity,dxVect);
      //this generates the new EBIS
      ebisPtr->define(a_domain, origin, dx, workshop);
    }
  else if (whichgeom == 5)
    {
      pout() << "sphere geometry" << endl;
      vector<Real> sphere_center(SpaceDim);
      pp.getarr("sphere_center",sphere_center, 0, SpaceDim);
      RealVect sphereCenter;
      for (int idir = 0; idir < SpaceDim; idir++)
        {
          sphereCenter[idir] = sphere_center[idir];
        }
      Real sphereRadius;
      pp.get("sphere_radius", sphereRadius);

      bool     insideRegular = false;
      SphereIF implicit(sphereRadius,sphereCenter,insideRegular);
      GeometryShop workshop(implicit,verbosity,dxVect);
      //this generates the new EBIS
      ebisPtr->define(a_domain, origin, dx, workshop);
    }
  else if (whichgeom == 13)
    {
      pout() << "rhodonea geometry" << endl;
      vector<Real> tmp(SpaceDim);
      pp.getarr("rhodonea_center", tmp, 0, SpaceDim);
//.........这里部分代码省略.........