src/Boundary/SphericalBoundary.cc
//------------------------------------------------------------------------------
// Set the ghost nodes for the given NodeList.
//------------------------------------------------------------------------------
Set the ghost nodes for the given NodeList.
Set the ghost nodes for the given NodeList .
void
SphericalBoundary::
setGhostNodes(NodeList<Dim<3> >& nodeList) {
  REQUIRE(mGhostPositions.fieldForNodeList(nodeList) < mGhostPositions.end());
  addNodeList(nodeList);
  BoundaryNodes& boundaryNodes = accessBoundaryNodes(nodeList);
  vector<int>& controlNodes = boundaryNodes.controlNodes;
  vector<int>& ghostNodes = boundaryNodes.ghostNodes;
  controlNodes = vector<int>();
  ghostNodes = vector<int>();
  Field<Dim<3>, vector<Vector> >& ghostPositions = **(mGhostPositions.fieldForNodeList(nodeList));
  const int currentNumGhostNodes = nodeList.numGhostNodes();
  int ghostNodeIndex = nodeList.numNodes();
  Field<Dimension, Vector>& position = nodeList.positions();
  Field<Dimension, SymTensor>& H = nodeList.Hfield();
  const double nPerh = nodeList.nodesPerSmoothingScale();
  const double kernelExtent = nodeList.neighbor().kernelExtent();
  CHECK(distinctlyGreaterThan(nPerh, 0.0));
  CHECK(distinctlyGreaterThan(kernelExtent, 0.0));
  const double phicrit = 0.25*M_PI/(nPerh*kernelExtent);
  for (auto i = 0u; i != nodeList.numInternalNodes(); ++i) {
    if (ghostPositions(i).size() == 0) {
      const double ri = position(i).x();
      const SymTensor& Hi = H(i);
      const double hyz = sqrt((Hi*Vector(1.0, 0.0, 0.0)).magnitude()/Hi.Determinant());
      const double delta = hyz/nPerh;
      const double dphi = 2.0*atan2(delta, 2.0*ri);
      if (dphi < phicrit) {
        const int nbins = int(kernelExtent*nPerh + 0.5);
        CHECK(nbins >= 1);
        for (int ibin = 0; ibin != nbins; ++ibin) {
          const double phij = (ibin + 1)*dphi;
          CHECK(phij < 0.5*M_PI);
          const double c = ri*sqrt(2.0*(1.0 - cos(phij)));
          const double ryz = c*sin(0.5*(M_PI - phij));
          const double circ = 2.0*M_PI*ryz;
          const int nTheta = max(1, int(circ/delta + 0.5));
          const double dTheta = 2.0*M_PI/nTheta;
          for (int j = 0; j != nTheta; ++j) {
            const double thetaj = (j + 0.5)*dTheta;
            ghostPositions(i).push_back(Vector(ri*cos(phij),
                                               ri*cos(thetaj)*sin(phij),
                                               ri*sin(thetaj)*sin(phij)));
          }
        }
      } else {
        const int nbins = max(1, int(0.5*M_PI/dphi + 0.5));
        CHECK(nbins >= 1);
        const double dphieff = 0.5*M_PI/nbins;
        for (int ibin = 0; ibin != nbins; ++ibin) {
          const double phij = (ibin + 1)*dphieff;
          CHECK(phij <= 0.5*M_PI);
          const double c = ri*sqrt(2.0*(1.0 - cos(phij)));
          const double ryz = c*sin(0.5*(M_PI - phij));
          const double circ = 2.0*M_PI*ryz;
          const int nTheta = max(1, int(circ/delta + 0.5));
          const double dTheta = 2.0*M_PI/nTheta;
          for (int j = 0; j != nTheta; ++j) {
            const double thetaj = (j + 0.5)*dTheta;
            ghostPositions(i).push_back(Vector(ri*cos(phij),
                                               ri*cos(thetaj)*sin(phij),
                                               ri*sin(thetaj)*sin(phij)));
            ghostPositions(i).push_back(Vector(ri*cos(phij + M_PI),
                                               ri*cos(thetaj)*sin(phij + M_PI),
                                               ri*sin(thetaj)*sin(phij + M_PI)));
          }
        }
        ghostPositions(i).push_back(Vector(-ri, 0.0, 0.0));
      }
    }
    const int nghosti = ghostPositions(i).size();
    CHECK(nghosti > 0);
    controlNodes.reserve(controlNodes.size() + nghosti);
    ghostNodes.reserve(ghostNodes.size() + nghosti);
    for (int j = 0; j != nghosti; ++j, ++ghostNodeIndex) {
      controlNodes.push_back(i);
      ghostNodes.push_back(ghostNodeIndex);
    }
  }
  const int numNewGhostNodes = ghostNodes.size();
  CHECK((int)controlNodes.size() == numNewGhostNodes);
  nodeList.numGhostNodes(currentNumGhostNodes + numNewGhostNodes);
  this->updateGhostNodes(nodeList);
}
void SphericalBoundary :: setGhostNodes ( NodeList < Dim < 3 > > & nodeList ) { REQUIRE ( mGhostPositions . fieldForNodeList ( nodeList ) < mGhostPositions . end ( ) ) ; addNodeList ( nodeList ) ; BoundaryNodes & boundaryNodes = accessBoundaryNodes ( nodeList ) ; vector < int > & controlNodes = boundaryNodes . controlNodes ; vector < int > & ghostNodes = boundaryNodes . ghostNodes ; controlNodes = vector < int > ( ) ; ghostNodes = vector < int > ( ) ; Field < Dim < 3 > , vector < Vector > > & ghostPositions = * * ( mGhostPositions . fieldForNodeList ( nodeList ) ) ; const int currentNumGhostNodes = nodeList . numGhostNodes ( ) ; int ghostNodeIndex = nodeList . numNodes ( ) ; Field < Dimension , Vector > & position = nodeList . positions ( ) ; Field < Dimension , SymTensor > & H = nodeList . Hfield ( ) ; const double nPerh = nodeList . nodesPerSmoothingScale ( ) ; const double kernelExtent = nodeList . neighbor ( ) . kernelExtent ( ) ; CHECK ( distinctlyGreaterThan ( nPerh , 0.0 ) ) ; CHECK ( distinctlyGreaterThan ( kernelExtent , 0.0 ) ) ; const double phicrit = 0.25 * M_PI / ( nPerh * kernelExtent ) ; for ( auto i = 0u ; i != nodeList . numInternalNodes ( ) ; ++ i ) { if ( ghostPositions ( i ) . size ( ) == 0 ) { const double ri = position ( i ) . x ( ) ; const SymTensor & Hi = H ( i ) ; const double hyz = sqrt ( ( Hi * Vector ( 1.0 , 0.0 , 0.0 ) ) . magnitude ( ) / Hi . Determinant ( ) ) ; const double delta = hyz / nPerh ; const double dphi = 2.0 * atan2 ( delta , 2.0 * ri ) ; if ( dphi < phicrit ) { const int nbins = int ( kernelExtent * nPerh + 0.5 ) ; CHECK ( nbins >= 1 ) ; for ( int ibin = 0 ; ibin != nbins ; ++ ibin ) { const double phij = ( ibin + 1 ) * dphi ; CHECK ( phij < 0.5 * M_PI ) ; const double c = ri * sqrt ( 2.0 * ( 1.0 - cos ( phij ) ) ) ; const double ryz = c * sin ( 0.5 * ( M_PI - phij ) ) ; const double circ = 2.0 * M_PI * ryz ; const int nTheta = max ( 1 , int ( circ / delta + 0.5 ) ) ; const double dTheta = 2.0 * M_PI / nTheta ; for ( int j = 0 ; j != nTheta ; ++ j ) { const double thetaj = ( j + 0.5 ) * dTheta ; ghostPositions ( i ) . push_back ( Vector ( ri * cos ( phij ) , ri * cos ( thetaj ) * sin ( phij ) , ri * sin ( thetaj ) * sin ( phij ) ) ) ; } } } else { const int nbins = max ( 1 , int ( 0.5 * M_PI / dphi + 0.5 ) ) ; CHECK ( nbins >= 1 ) ; const double dphieff = 0.5 * M_PI / nbins ; for ( int ibin = 0 ; ibin != nbins ; ++ ibin ) { const double phij = ( ibin + 1 ) * dphieff ; CHECK ( phij <= 0.5 * M_PI ) ; const double c = ri * sqrt ( 2.0 * ( 1.0 - cos ( phij ) ) ) ; const double ryz = c * sin ( 0.5 * ( M_PI - phij ) ) ; const double circ = 2.0 * M_PI * ryz ; const int nTheta = max ( 1 , int ( circ / delta + 0.5 ) ) ; const double dTheta = 2.0 * M_PI / nTheta ; for ( int j = 0 ; j != nTheta ; ++ j ) { const double thetaj = ( j + 0.5 ) * dTheta ; ghostPositions ( i ) . push_back ( Vector ( ri * cos ( phij ) , ri * cos ( thetaj ) * sin ( phij ) , ri * sin ( thetaj ) * sin ( phij ) ) ) ; ghostPositions ( i ) . push_back ( Vector ( ri * cos ( phij + M_PI ) , ri * cos ( thetaj ) * sin ( phij + M_PI ) , ri * sin ( thetaj ) * sin ( phij + M_PI ) ) ) ; } } ghostPositions ( i ) . push_back ( Vector ( - ri , 0.0 , 0.0 ) ) ; } } const int nghosti = ghostPositions ( i ) . size ( ) ; CHECK ( nghosti > 0 ) ; controlNodes . reserve ( controlNodes . size ( ) + nghosti ) ; ghostNodes . reserve ( ghostNodes . size ( ) + nghosti ) ; for ( int j = 0 ; j != nghosti ; ++ j , ++ ghostNodeIndex ) { controlNodes . push_back ( i ) ; ghostNodes . push_back ( ghostNodeIndex ) ; } } const int numNewGhostNodes = ghostNodes . size ( ) ; CHECK ( ( int ) controlNodes . size ( ) == numNewGhostNodes ) ; nodeList . numGhostNodes ( currentNumGhostNodes + numNewGhostNodes ) ; this -> updateGhostNodes ( nodeList ) ; }
Set the ghost nodes for the given NodeList.
Set the ghost nodes for the given NodeList .
// Add this NodeList, creating space for control & ghost nodes. // The field of ghost positions for this NodeList. // The ghost node indices. // Get the state for this NodeList. // Get the number of nodes per h, and the kernel extent. // The critical phi angle, at which the surface of ghost nodes will  // extend for 45 degrees around the full sphere.  At this point we // decide it is worth simply filling in the entire spherical surface // with ghost nodes. // Iterate over the internal nodes. // Have we computed the set of ghost nodes for this node yet? // This is a new node we have not computed ghost nodes for yet. // Go ahead and create the appropriate set of ghost nodes for // it. // Determine H extent in the yz plane. // The delta expected between nodes, and corresponding angle. // Do we need to create a full sphere of ghost nodes? // We're outside the critical point for doing the entire sphere, so  // just create a section of the surface with enough ghost nodes that  // this one shouldn't see beyond them. // We're near enough the origin that we need to create the full sphere // of ghost points at this radius. // Create the mirror image of the control node across the origin. // Set the control and ghost node indices for this internal node. // Create the set of ghost nodes on the NodeList. // Use the updateGhostNodes method to set the node positions and H's.