moab
moab::Intx2Mesh Class Reference

#include <Intx2Mesh.hpp>

Inheritance diagram for moab::Intx2Mesh:
moab::Intx2MeshInPlane moab::Intx2MeshOnSphere

List of all members.

Public Member Functions

 Intx2Mesh (Interface *mbimpl)
virtual ~Intx2Mesh ()
ErrorCode intersect_meshes (EntityHandle mbs1, EntityHandle mbs2, EntityHandle &outputSet)
virtual int computeIntersectionBetweenRedAndBlue (EntityHandle red, EntityHandle blue, double *P, int &nP, double &area, int markb[MAXEDGES], int markr[MAXEDGES], int &nsidesBlue, int &nsidesRed, bool check_boxes_first=false)=0
virtual int findNodes (EntityHandle red, int nsRed, EntityHandle blue, int nsBlue, double *iP, int nP)=0
virtual void createTags ()
ErrorCode GetOrderedNeighbors (EntityHandle set, EntityHandle quad, EntityHandle neighbors[MAXEDGES])
void SetErrorTolerance (double eps)
void clean ()
virtual bool is_inside_element (double xyz[3], EntityHandle eh)=0
void set_box_error (double berror)
ErrorCode create_departure_mesh_2nd_alg (EntityHandle &euler_set, EntityHandle &covering_lagr_set)
ErrorCode create_departure_mesh_3rd_alg (EntityHandle &lagr_set, EntityHandle &covering_set)
ErrorCode build_processor_euler_boxes (EntityHandle euler_set, Range &local_verts)
void correct_polygon (EntityHandle *foundIds, int &nP)
void enable_debug ()
void disable_debug ()

Protected Attributes

Interfacemb
EntityHandle mbs1
EntityHandle mbs2
Range rs1
Range rs2
EntityHandle outSet
Tag BlueFlagTag
Tag RedFlagTag
Range RedEdges
Tag redParentTag
Tag blueParentTag
Tag countTag
const EntityHandleredConn
const EntityHandleblueConn
CartVect redCoords [MAXEDGES]
CartVect blueCoords [MAXEDGES]
double redCoords2D [MAXEDGES2]
double blueCoords2D [MAXEDGES2]
std::ofstream mout_1 [6]
int dbg_1
std::vector< std::vector
< EntityHandle > * > 
extraNodesVec
double epsilon_1
double epsilon_area
std::vector< double > allBoxes
double box_error
EntityHandle localRoot
Range localEnts
unsigned int my_rank
int max_edges

Detailed Description

Definition at line 41 of file Intx2Mesh.hpp.


Constructor & Destructor Documentation

Definition at line 22 of file Intx2Mesh.cpp.

                                      : mb(mbimpl)
#ifdef USE_MPI
   , parcomm(NULL), remote_cells(NULL)
#endif
{
  dbg_1=0;
  box_error=0;
  my_rank=0;
  BlueFlagTag=0;
  RedFlagTag=0;
  redParentTag =0;
  blueParentTag = 0;
  countTag = 0;
}

Definition at line 37 of file Intx2Mesh.cpp.

{
  // TODO Auto-generated destructor stub
#ifdef USE_MPI
  if (remote_cells)
  {
    delete remote_cells;
    remote_cells=NULL;
  }
#endif
}

Member Function Documentation

Definition at line 425 of file Intx2Mesh.cpp.

{
  //
  int indx = 0;
  for (Range::iterator eit = RedEdges.begin(); eit != RedEdges.end();
      eit++, indx++)
  {
    //EntityHandle edge = *eit;
    //delete extraNodesMap[edge];
    delete extraNodesVec[indx];
  }
  //extraNodesMap.clear();
  extraNodesVec.clear();
  // also, delete some bit tags, used to mark processed reds and blues
  mb->tag_delete(RedFlagTag);// to mark blue quads already considered
  mb->tag_delete(BlueFlagTag);

  /*mb->tag_delete(redParentTag);
  mb->tag_delete(blueParentTag);
  mb->tag_delete(countTag);
  RedEdges.clear();
  localEnts.clear();*/

}
virtual int moab::Intx2Mesh::computeIntersectionBetweenRedAndBlue ( EntityHandle  red,
EntityHandle  blue,
double *  P,
int &  nP,
double &  area,
int  markb[MAXEDGES],
int  markr[MAXEDGES],
int &  nsidesBlue,
int &  nsidesRed,
bool  check_boxes_first = false 
) [pure virtual]
void moab::Intx2Mesh::correct_polygon ( EntityHandle foundIds,
int &  nP 
)

Definition at line 452 of file Intx2Mesh.cpp.

{
  int i = 0;
  while(i<nP)
  {
    int nextIndex = (i+1)%nP;
    if (nodes[i]==nodes[nextIndex])
    {
      // we need to reduce nP, and collapse nodes
      if (dbg_1)
      {
        std::cout<<" nodes duplicated in list: " ;
        for (int j=0; j<nP; j++)
          std::cout<<nodes[j] << " " ;
        std::cout<<"\n";
        std::cout<<" node " << nodes[i] << " at index " << i << " is duplicated" << "\n";
      }
      // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0, then next thing does nothing
      //  (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
      for (int k=i; k<nP-1; k++)
        nodes[k] = nodes[k+1];
      nP--; // decrease the number of nodes; also, decrease i, just if we may need to check again
      i--;
    }
    i++;
  }
  return;
}
void moab::Intx2Mesh::createTags ( ) [virtual]

Definition at line 48 of file Intx2Mesh.cpp.

{
  if (redParentTag)
    mb->tag_delete(redParentTag);
  if(blueParentTag)
    mb->tag_delete(blueParentTag);
  if (countTag)
    mb->tag_delete(countTag);
    /*RedEdges.clear();
    localEnts.clear()*/

  unsigned char def_data_bit = 0; // unused by default
  ErrorCode rval = mb->tag_get_handle("blueFlag", 1, MB_TYPE_BIT, BlueFlagTag,
      MB_TAG_CREAT, &def_data_bit);
  if (MB_SUCCESS != rval)
    return;
  // maybe the red tag is better to be deleted every time, and recreated;
  // or is it easy to set all values to something again? like 0?
  rval = mb->tag_get_handle("redFlag", 1, MB_TYPE_BIT, RedFlagTag, MB_TAG_CREAT,
      &def_data_bit);
  ERRORV(rval, "can't get red flag tag");

  // assume that the edges are on the red triangles
  Range redElements;
  //Range redEdges;
  rval = mb->get_entities_by_dimension(mbs2, 2, redElements, false); // so all tri, quad, poly
  ERRORV(rval, "can't get ents by dimension");

  // create red edges if they do not exist yet; so when they are looked upon, they are found
  // this is the only call that is potentially NlogN, in the whole method
  rval = mb->get_adjacencies(redElements, 1, true, RedEdges, Interface::UNION);
  ERRORV(rval, "can't get adjacent red edges");

  // now, create a map from each edge to a list of potential new nodes on a red edge
  // this memory has to be cleaned up
  // change it to a vector, and use the index in range of red edges
  int indx = 0;
  extraNodesVec.reserve(RedEdges.size());
  for (Range::iterator eit = RedEdges.begin(); eit != RedEdges.end();
      eit++, indx++)
  {
    //EntityHandle edge = *eit;
    //extraNodesMap[edge] = new std::vector<EntityHandle>;
    std::vector<EntityHandle> * nv = new std::vector<EntityHandle>;
    extraNodesVec.push_back(nv);
  }

  int defaultInt = 0;

  rval = mb->tag_get_handle("RedParent", 1, MB_TYPE_INTEGER, redParentTag,
      MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt);
  ERRORV(rval, "can't create positive tag");

  rval = mb->tag_get_handle("BlueParent", 1, MB_TYPE_INTEGER, blueParentTag,
      MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt);
  ERRORV(rval, "can't create negative tag");

  rval = mb->tag_get_handle("Counting", 1, MB_TYPE_INTEGER, countTag,
        MB_TAG_SPARSE | MB_TAG_CREAT, &defaultInt);
  ERRORV(rval, "can't create Counting tag");

  return;
}
void moab::Intx2Mesh::disable_debug ( ) [inline]

Definition at line 102 of file Intx2Mesh.hpp.

{dbg_1 = 0;}
void moab::Intx2Mesh::enable_debug ( ) [inline]

Definition at line 101 of file Intx2Mesh.hpp.

{dbg_1 = 1;}
virtual int moab::Intx2Mesh::findNodes ( EntityHandle  red,
int  nsRed,
EntityHandle  blue,
int  nsBlue,
double *  iP,
int  nP 
) [pure virtual]

Definition at line 116 of file Intx2Mesh.cpp.

{
  int nnodes = 3;
  // will get the nnodes ordered neighbors;
  // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
  const EntityHandle * conn4;
  ErrorCode rval = mb->get_connectivity(cell, conn4, nnodes);
  int nsides = nnodes;
  // account for possible padded polygons
  while (conn4[nsides-2]==conn4[nsides-1] && nsides>3)
    nsides--;
  ERRORR(rval, "can't get connectivity on an element");
  for (int i = 0; i < nsides; i++)
  {
    EntityHandle v[2];
    v[0] = conn4[i];
    v[1] = conn4[(i + 1) % nsides];
    // get quads adjacent to vertices
    std::vector<EntityHandle> cells;
    std::vector<EntityHandle> cellsInSet;
    rval = mb->get_adjacencies(v, 2, 2, false, cells, Interface::INTERSECT);
    ERRORR(rval, "can't get adjacencies on 2 nodes");
    size_t siz = cells.size();
    for (size_t j = 0; j < siz; j++)
      if (mb->contains_entities(set, &(cells[j]), 1))
        cellsInSet.push_back(cells[j]);
    siz = cellsInSet.size();

    if (siz > 2)
    {
      std::cout << "non manifold mesh, error"
          << mb->list_entities(&(cellsInSet[0]), cellsInSet.size()) << "\n";
      return MB_FAILURE; // non-manifold
    }
    if (siz == 1)
    {
      // it must be the border,
      neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border
      // borders do not appear for a sphere in serial, but they do appear for
      // parallel processing anyway
      continue;
    }
    // here siz ==2, it is either the first or second
    if (cell == cellsInSet[0])
      neighbors[i] = cellsInSet[1];
    else
      neighbors[i] = cellsInSet[0];
  }
  return MB_SUCCESS;
}

Definition at line 169 of file Intx2Mesh.cpp.

{

  ErrorCode rval;
  mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
  mbs2 = mbset2;
  outSet = outputSet;

  // really, should be something from t1 and t2; blue is 1 (lagrange), red is 2 (euler)
  createTags(); //
  EntityHandle startBlue, startRed;

  mb->get_entities_by_dimension(mbs1, 2, rs1);
  mb->get_entities_by_dimension(mbs2, 2, rs2);
  Range rs22=rs2; // a copy of the initial range; we will remove from it elements as we
                 // advance ; rs2 is needed for marking the polygon to the red parent
  while (!rs22.empty())
  {
    for (Range::iterator it = rs1.begin(); it != rs1.end(); it++)
    {
      startBlue = *it;
      int found = 0;
      for (Range::iterator it2 = rs22.begin(); it2 != rs22.end() && !found; it2++)
      {
        startRed = *it2;
        double area = 0;
        // if area is > 0 , we have intersections
        double P[10*MAXEDGES]; // max 8 intx points + 8 more in the polygon
        //
        int nP = 0;
        int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
        int nsRed, nsBlue;
        computeIntersectionBetweenRedAndBlue(startRed, startBlue, P, nP, area, nb, nr,
            nsBlue, nsRed, true);
        if (area > 0)
        {
          found = 1;
          break; // found 2 elements that intersect; these will be the seeds
        }
      }
      if (found)
        break;
    }

    std::queue<EntityHandle> blueQueue; // these are corresponding to Ta,
    blueQueue.push(startBlue);
    std::queue<EntityHandle> redQueue;
    redQueue.push(startRed);

    Range toResetBlues; // will be used to reset blue flags for every red quad
    // processed

    /*if (my_rank==0)
      dbg_1 = 1;*/
    unsigned char used = 1;
    unsigned char unused = 0; // for red flags
    // mark the start blue quad as used, so it will not come back again
    mb->tag_set_data(RedFlagTag, &startRed, 1, &used);
    while (!redQueue.empty())
    {
      // flags for the side : 0 means a blue quad not found on side
      // a paired blue not found yet for the neighbors of red
      EntityHandle n[MAXEDGES] = { EntityHandle(0) };

      int nsidesRed; // will be initialized later, when we compute intersection
      EntityHandle currentRed = redQueue.front();

      redQueue.pop();
      //        for (k=0; k<m_numPos; k++)
      //          redFlag[k] = 0;
      //        redFlag[m_numPos] = 1; // to guard for the boundary
      // all reds that were tagged, are now cleared
      for (Range::iterator itr = toResetBlues.begin(); itr != toResetBlues.end();
          itr++)
      {
        EntityHandle ttt = *itr;
        rval = mb->tag_set_data(BlueFlagTag, &ttt, 1, &unused);
        ERRORR(rval, "can't set blue unused tag");
      }
      //rval = mb2->tag_set_data(RedFlagTag, toResetReds, &unused);
      /*if (dbg_1)
      {
        std::cout << "reset blues: ";
        mb->list_entities(toResetBlues);
      }*/
      //rval = mb2->tag_set_data(RedFlagTag, toResetReds, &unused);
      if (dbg_1)
      {
        std::cout << "reset blues: ";
        for (Range::iterator itr = toResetBlues.begin(); itr != toResetBlues.end();
            itr++)
          std::cout << mb->id_from_handle(*itr) << " ";
        std::cout << std::endl;
      }
      EntityHandle currentBlue = blueQueue.front(); // where do we check for redQueue????
      // red and blue queues are parallel
      blueQueue.pop(); // mark the current red
      //redFlag[currentRed] = 1; //
      toResetBlues.clear(); // empty the range of used blues, will have to be set unused again,
      // at the end of red element processing
      toResetBlues.insert(currentBlue);
      rval = mb->tag_set_data(BlueFlagTag, &currentBlue, 1, &used);
      ERRORR(rval, "can't set blue tag");
      //mb2->set_tag_data
      std::queue<EntityHandle> localBlue;
      localBlue.push(currentBlue);
      while (!localBlue.empty())
      {
        //
        EntityHandle blueT = localBlue.front();
        localBlue.pop();
        double P[10*MAXEDGES], area; //
        int nP = 0;
        int nb[MAXEDGES] = {0};
        int nr[MAXEDGES] = {0};

        int nsidesBlue; 
        // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
        // intersection points could include the vertices of initial elements
        // nb [j] = 0 means no intersection on the side k for element blue (markers)
        // nb [j] = 1 means that the side j (from j to j+1) of blue poly intersects the
        // red poly.  A potential next poly is the red poly that is adjacent to this side
        computeIntersectionBetweenRedAndBlue(/* red */currentRed, blueT, P, nP,
            area, nb, nr, nsidesBlue, nsidesRed);
        if (nP > 0)
        {
          if (dbg_1)
          {
            for (int k=0; k<3; k++)
            {
              std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
            }
          }
          // intersection found: output P and original triangles if nP > 2

          EntityHandle neighbors[MAXEDGES];
          rval = GetOrderedNeighbors(mbs1, blueT, neighbors);
          if (rval != MB_SUCCESS)
          {
            std::cout << " can't get the neighbors for blue element "
                << mb->id_from_handle(blueT);
            return MB_FAILURE;
          }

          // add neighbors to the localBlue queue, if they are not marked
          for (int nn = 0; nn < nsidesBlue; nn++)
          {
            EntityHandle neighbor = neighbors[nn];
            if (neighbor > 0 && nb[nn]>0) // advance across blue boundary n
            {
              //n[nn] = redT; // start from 0!!
              unsigned char status = 0;
              mb->tag_get_data(BlueFlagTag, &neighbor, 1, &status);
              if (status == 0)
              {
                localBlue.push(neighbor);
                /*if (dbg_1)
                {
                  std::cout << " local blue elem "
                      << mb->list_entities(&neighbor, 1) << "\n  for red:"
                      << mb->list_entities(&currentRed, 1) << "\n";
                }*/
                if (dbg_1)
                {
                  std::cout << " local blue elem " << mb->id_from_handle(neighbor)
                      << " for red:" << mb->id_from_handle(currentRed) << "\n";
                  mb->list_entities(&neighbor, 1);
                }
                rval = mb->tag_set_data(BlueFlagTag, &neighbor, 1, &used);
                //redFlag[neighbor] = 1; // flag it to not be added anymore
                toResetBlues.insert(neighbor); // this is used to reset the red flag
              }
            }
            // n(find(nc>0))=ac;        % ac is starting candidate for neighbor
            if (nr[nn] > 0)
              n[nn] = blueT;

          }
          if (nP > 1) // this will also construct triangles/polygons in the new mesh, if needed
            findNodes(currentRed, nsidesRed, blueT, nsidesBlue, P, nP);
        }
        else if (dbg_1)
        /*{
          std::cout << " red, blue, do not intersect: "
              << mb->list_entities(&currentRed, 1) << " "
              << mb->list_entities(&blueT, 1) << "\n";
        }*/
        {
          std::cout << " red, blue, do not intersect: "
              << mb->id_from_handle(currentRed) << " "
              << mb->id_from_handle(blueT) << "\n";
        }

      } // end while (!localBlue.empty())
      // here, we are finished with redCurrent, take it out of the rs22 range (red, arrival mesh)
      rs22.erase(currentRed);
      // also, look at its neighbors, and add to the seeds a next one

      EntityHandle redNeighbors[MAXEDGES];
      rval = GetOrderedNeighbors(mbs2, currentRed, redNeighbors);
      ERRORR(rval, "can't get neighbors");
      if (dbg_1)
      {
        std::cout << "Next: neighbors for current red ";
        for (int kk = 0; kk < nsidesRed; kk++)
        {
          if (redNeighbors[kk] > 0)
            std::cout << mb->id_from_handle(redNeighbors[kk]) << " ";
          else
            std::cout << 0 << " ";
        }
        std::cout << std::endl;
      }
      for (int j = 0; j < nsidesRed; j++)
      {
        EntityHandle redNeigh = redNeighbors[j];
        unsigned char status = 1;
        if (redNeigh == 0)
          continue;
        mb->tag_get_data(RedFlagTag, &redNeigh, 1, &status); // status 0 is unused
        if (status == 0 && n[j] > 0) // not treated yet and marked as a neighbor
        {
          // we identified red quad n[j] as intersecting with neighbor j of the blue quad
          redQueue.push(redNeigh);
          blueQueue.push(n[j]);
          if (dbg_1)
            std::cout << "new polys pushed: blue, red:"
                << mb->id_from_handle(redNeigh) << " "
                << mb->id_from_handle(n[j]) << std::endl;
          mb->tag_set_data(RedFlagTag, &redNeigh, 1, &used);
        }
      }

    } // end while (!redQueue.empty())
  }
  if (dbg_1)
  {
    for (int k = 0; k < 6; k++)
      mout_1[k].close();
  }
  // before cleaning up , we need to settle the position of the intersection points
  // on the boundary edges
  // this needs to be collective, so we should maybe wait something
#ifdef USE_MPI
  rval = correct_intersection_points_positions();
  if (rval!=MB_SUCCESS)
  {
    std::cout << "can't correct position, Intx2Mesh.cpp \n";
  }
#endif
  clean();
  return MB_SUCCESS;
}
virtual bool moab::Intx2Mesh::is_inside_element ( double  xyz[3],
EntityHandle  eh 
) [pure virtual]
void moab::Intx2Mesh::set_box_error ( double  berror) [inline]

Definition at line 81 of file Intx2Mesh.hpp.

   {box_error = berror;}
void moab::Intx2Mesh::SetErrorTolerance ( double  eps) [inline]

Definition at line 72 of file Intx2Mesh.hpp.

{ epsilon_1=eps; epsilon_area = eps*eps/2;}

Member Data Documentation

std::vector<double> moab::Intx2Mesh::allBoxes [protected]

Definition at line 146 of file Intx2Mesh.hpp.

Definition at line 129 of file Intx2Mesh.hpp.

Definition at line 131 of file Intx2Mesh.hpp.

Definition at line 133 of file Intx2Mesh.hpp.

Definition at line 114 of file Intx2Mesh.hpp.

Definition at line 123 of file Intx2Mesh.hpp.

double moab::Intx2Mesh::box_error [protected]

Definition at line 147 of file Intx2Mesh.hpp.

Definition at line 124 of file Intx2Mesh.hpp.

int moab::Intx2Mesh::dbg_1 [protected]

Definition at line 136 of file Intx2Mesh.hpp.

double moab::Intx2Mesh::epsilon_1 [protected]

Definition at line 143 of file Intx2Mesh.hpp.

double moab::Intx2Mesh::epsilon_area [protected]

Definition at line 144 of file Intx2Mesh.hpp.

std::vector<std::vector<EntityHandle> *> moab::Intx2Mesh::extraNodesVec [protected]

Definition at line 141 of file Intx2Mesh.hpp.

Definition at line 150 of file Intx2Mesh.hpp.

Definition at line 149 of file Intx2Mesh.hpp.

int moab::Intx2Mesh::max_edges [protected]

Definition at line 158 of file Intx2Mesh.hpp.

Definition at line 104 of file Intx2Mesh.hpp.

Definition at line 106 of file Intx2Mesh.hpp.

Definition at line 107 of file Intx2Mesh.hpp.

std::ofstream moab::Intx2Mesh::mout_1[6] [protected]

Definition at line 135 of file Intx2Mesh.hpp.

unsigned int moab::Intx2Mesh::my_rank [protected]

Definition at line 151 of file Intx2Mesh.hpp.

Definition at line 111 of file Intx2Mesh.hpp.

Definition at line 128 of file Intx2Mesh.hpp.

Definition at line 130 of file Intx2Mesh.hpp.

Definition at line 132 of file Intx2Mesh.hpp.

Definition at line 118 of file Intx2Mesh.hpp.

Definition at line 116 of file Intx2Mesh.hpp.

Definition at line 122 of file Intx2Mesh.hpp.

Definition at line 108 of file Intx2Mesh.hpp.

Definition at line 109 of file Intx2Mesh.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines