Here is an example of the evaluation of the Laplacian on a single level EBIndexSpace in three dimensions. There are three pieces to the code. First, we construct the IrregGeom on which the irregular differencing of Laplacian is to be performed. This consists of all of the irregular Node s, as well as the set of all regular Node s which have one or more irregular Node in its stencil.
// accrete code fragment.
EBIndexSpace isp;
IrregGeom ig();
...
ig.setIndexSpace(isp);
NodeIterator<EBIndexSpace>(isp) itn;
for (itn.begin(),itn.ok(),itn.next())
{
// add Nodes corresponding to irregular Nodes in isp.
IntVect iv = itn.getNode().gridIndex();
ig.addNodes(iv);
// add Nodes corresponding to cells within one cell of the current
grid index.
for (int ll = 0,ll < 3,ll++)
{
for (int s = -1, s < 2, s++)
{
IntVect is = BASISV(ll)*s;
ig.addNodes(iv + is);
}
}
}
The next step is to build the stencils used to apply the operator.
The strategy here is to use the 
flux calculation from
Johansen and Colella for edges for which it can be reliably computed,
and drop down to the simple centered difference stencil otherwise.
Since the latter is also for regular edges, the only place where
we are dropping order is for cases in which the geometry is seriously
underresolved.
void buildFluxStencils(
EBIndexSpace& isp,
IrregGeom& ig,
IEFab<NodeStencil>* fluxStencils,
IEFab<Vector<REAL>>* fluxCoefs
)
{
// Build stencils for computing higher order fluxes in 3D.
for (ll = 0,ll < 3,ll++)
{
int kPerp1 = (k + 1)%3;
int kPerp2 = (k + 2)%3;
fluxStencils[ll].define(ig,ll,SurroundingEdges,1);
fluxCoefs[ll].define(ig,ll,SurroundingEdges,1);
EdgeIterator<IrregGeom> ite(ig,ll,SurroundingEdges);
for (ite.begin(),ite.ok(),ite.next())
{
// Check conditions for using a nontrivial O(h**2) flux stencil.
// The conditions are:
// (1) Is this edge irregular ?
// (2) If so, are the components of the boundary normal in the directions
// tangent to the face of the same sign ?
// (3) If (1) or (2) are true, are the edges that will be used to
// computed the centered flues regular, and are the nodes making up
// those edges connected appropriately to the nodes making up
// edge at which the flux is being evaluated ?
// If any of these conditions fail, then we use the centered difference
// flux. This corresponds to dropping order for some of the irregular edges,
// but is still second order accurate for regular edges.
const Edge& ee = ite();
IntVect ive = ee.getIndex();
NodeStencil& fs = fluxStencil[ll](ee);
Vector<REAL>& fc = fluxCoefs[ll](ee);
bool makeSimpleStencil = 1;
if (!ee.isRegular())
{
// Condition (1) is met.
Node& nLow = ee.lowNode();
Node& nHigh = ee.highNode();
REAL normLow1 = nLow.getNormal(kPerp1);
REAL normHigh1 = nHigh.getNormal(kPerp1);
REAL normLow2 = nLow.getNormal(kPerp2);
REAL normHigh2 = nHigh.getNormal(kPerp2);
if ((normLow1*normHigh1 > 0) && (normLow2*normHigh2 > 0)
{
// Condition (2) is met.
REAL norm1 = .5*(normLow1 + normHigh1);
REAL norm2 = .5*(normLow2 + normHigh2);
int sign1;
int sign2;
if (norm1 > 0) {sign1 = 1} else {sign1 = -1};
if (norm2 > 0) {sign2 = 1} else {sign2 = -1};
IntVect ivPerp1;
ivPerp1.setVal(kPerp1,sign1);
IntVect ivPerp2;
ivPerp2.setVal(kPerp2,sign2);
ivPerp1 = ivPerp1 + ive;
ivPerp2 = ivPerp2 + ive;
bool edgesOk = (isp.isRegular(ivPerp1) &&
isp.isRegular(ivPerp1))
if (edgesOk)
// Get nodes for edges at ivPerp1 and ivPerp2 to check connectivity.
{
Edge e1 = getEdges(ivPerp1)[0];
Edge e2 = getEdges(ivPerp2)[0];
edgesOk = edgesOk &&
isConnected(e1.lowNode(),nLow);
edgesOk = edgesOk &&
isConnected(e2.lowNode(),nLow);
edgesOk = edgesOk &&
isConnected(e1.highNode(),nHigh);
edgesOk = edgesOk &&
isConnected(e2.highNode(),nHigh);
}
if (edgesOk)
{
// Condition (3) is met - construct the stencil.
// Compute the center of area for the face.
REAL xx1;
REAL xx2;
getCtrArea(ee.areaFrac(),norm1,norm2,xx1,xx2);
fs.addNode(ee.lowNode());
fc[0] = -(one-xx1-xx2);
fs.addNode(ee.highNode());
fc[1] = (one-xx1-xx2);
fs.addNode(ivPerp1.lowNode());
fc[2] = -xx1;
fs.addNode(ivPerp1.highNode());
fc[3] = xx1;
fs.addNode(ivPerp2.lowNode());
fc[4] = -xx2;
fs.addNode(ivPerp2.highNode());
fc[5] = xx2;
MakeSimpleStencil = 0;
}
}
}
if (MakeSimpleStencil)
{
// We don't need / have the more elaborate stencil
// - instead use simple centered differences.
fs.AddNode(ee.lowNode())
fc[0] = -one;
fs.AddNode(ee.highNode())
fc[1] = one;
} // End stencil construction for given edge.
} // End loop over edges for a given direction.
} // End loop over directions.
}
Finally, we specify how the operator is applied. We note that most of the lines of code are used to perform the irregular calculation. The regular calculation is performed using a call to a Fortran routine on the rectangular grid data.
void applyL(
EBIndexSpace& isp,
IrregGeom& ig,
FArrayBox& phi,
FArrayBox& LOfPhi,
INFab& iPhi,
INFab& iLOfPhi,
IEFab<NodeStencil>* fluxStencils,
IEFab<Vector<REAL>>* fluxCoefs,
Box D,
REAL h
)
{
// Declare and define irregular Flux fabs.
IEFab<REAL> Flux[3];
for (int ll = 0,ll < 3,ll++)
{
Flux[ll].define(ig,ll,SurroundingEdges,1);
EdgeIterator<IrregGeom> ite(ig,ll,SurroundingEdges);
for (ite.begin(),ite.ok(),ite.next())
{
NodeStencil& ns = fluxStencils[ll](ite());
Vector<REAL>& coefs = fluxCoefs[ll](ite());
REAL& f = Flux[ll](ite());
f = 0;
NodeIterator<NodeStencil> itnLoc;
for (itnLoc.begin(),itnLoc.ok(),itnLoc.next())
{
f = f + coefs[itnLoc.getInt()]*iPhi(itnLoc(),phi);
}
f = f/h;
}
}
// Compute flux divergence on ig.
NodeIterator<IrregGeom> itn(ig);
for (itn.begin(),itn.ok(),itn.next())
{
Node n = itn();
REAL& r = iLOfPhi(n);
for (int ll = 0,ll<3,ll++)
{
Vector<Edge> elo=n.lowEdge(ll);
Vector<Edge> ehi=n.highEdge(ll);
r = 0;
for (int le = 0,le < elo.length(),le++)
{
r = r - Flux[ll](elo[le])*elo[le].areaFrac();
}
for (int le = 0,le < ehi.length(),le++)
{
r = r + Flux[ll](ehi[le])*ehi[le].areaFrac();
}
}
r = r/h/n.volFrac();
}
// Copy values in regular cells of iPhi into corresponding locations in phi.
for (itn.begin(),itn.ok(),itn.next())
{
if (itn().isRegular()) phi(itn().gridIndex()) = iPhi(itn());
}
// Apply rectangular grid operator.
FORT_LAPLACIAN3D(
FAA(phi),
FAA(LOfPhi),
BAA(D),
&h);
}