[MOAB-dev] r2961 - in MOAB/branches/parallel_ghosting: . parallel

tautges at mcs.anl.gov tautges at mcs.anl.gov
Tue Jun 23 09:25:59 CDT 2009


Author: tautges
Date: 2009-06-23 09:25:58 -0500 (Tue, 23 Jun 2009)
New Revision: 2961

Modified:
   MOAB/branches/parallel_ghosting/AEntityFactory.cpp
   MOAB/branches/parallel_ghosting/HigherOrderFactory.cpp
   MOAB/branches/parallel_ghosting/MBCN.cpp
   MOAB/branches/parallel_ghosting/MBCN.hpp
   MOAB/branches/parallel_ghosting/MBCNArrays.hpp
   MOAB/branches/parallel_ghosting/MBSkinner.cpp
   MOAB/branches/parallel_ghosting/parallel/pcomm_unit.cpp
Log:
Back-porting MBCN and skinner changes from trunk, to fix problem 
calling skinner.

pcomm_unit works on 3d test now, and so do parallel ghost tests.  Still
dying in parallel_unit_tests.



Modified: MOAB/branches/parallel_ghosting/AEntityFactory.cpp
===================================================================
--- MOAB/branches/parallel_ghosting/AEntityFactory.cpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/AEntityFactory.cpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -806,7 +806,7 @@
     const MBCN::ConnMap &cmap = 
       MBCN::mConnectivityMap[source_type][target_dimension-1];
   
-    int verts_per_sub = cmap.num_nodes_per_sub_element[j];
+    int verts_per_sub = cmap.num_corners_per_sub_element[j];
 
       // get the corner vertices
     for (int i = 0; i < verts_per_sub; i++)

Modified: MOAB/branches/parallel_ghosting/HigherOrderFactory.cpp
===================================================================
--- MOAB/branches/parallel_ghosting/HigherOrderFactory.cpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/HigherOrderFactory.cpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -343,7 +343,7 @@
       tmp_face_conn[0] = element[entity_faces.conn[i][0]];
       tmp_face_conn[1] = element[entity_faces.conn[i][1]];
       tmp_face_conn[2] = element[entity_faces.conn[i][2]];
-      if(entity_faces.num_nodes_per_sub_element[i] == 4)
+      if(entity_faces.num_corners_per_sub_element[i] == 4)
         tmp_face_conn[3] = element[entity_faces.conn[i][3]];
       else
         tmp_face_conn[3] = 0;
@@ -359,7 +359,7 @@
       {
         EntitySequence* tmp_sequence = NULL;
         double sum_coords[3] = {0,0,0};
-        int max_nodes = entity_faces.num_nodes_per_sub_element[i];
+        int max_nodes = entity_faces.num_corners_per_sub_element[i];
         for(int k=0; k<max_nodes; k++)
         {
           seq_manager->find(tmp_face_conn[k], tmp_sequence);

Modified: MOAB/branches/parallel_ghosting/MBCN.cpp
===================================================================
--- MOAB/branches/parallel_ghosting/MBCN.cpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/MBCN.cpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -108,7 +108,7 @@
          ((source_dim > 0 && 
            *it1 < mConnectivityMap[this_type][source_dim-1].num_sub_elements) ||
           (source_dim == 0 && 
-           *it1 < mConnectivityMap[this_type][Dimension(this_type)-1].num_nodes_per_sub_element[0])) && 
+           *it1 < mConnectivityMap[this_type][Dimension(this_type)-1].num_corners_per_sub_element[0])) && 
          *it1 >= 0);
 
 
@@ -122,7 +122,7 @@
       // less than source_dim, which should give the connectivity of that sub element
     const ConnMap &cm = mConnectivityMap[this_type][source_dim-1];
     std::copy(cm.conn[source_indices[0]],
-              cm.conn[source_indices[0]]+cm.num_nodes_per_sub_element[source_indices[0]],
+              cm.conn[source_indices[0]]+cm.num_corners_per_sub_element[source_indices[0]],
               std::back_inserter(index_list));
     return 0;
   }
@@ -652,6 +652,49 @@
                                sub_index, sub_entity_conn);
 }
 
+void MBCN::SubEntityNodeIndices( const MBEntityType this_topo, 
+                                 const int num_nodes,
+                                 const int sub_dimension,
+                                 const int sub_index,
+                                 MBEntityType& subentity_topo,
+                                 int& num_sub_entity_nodes,
+                                 int sub_entity_conn[] )
+{
+    // If asked for a node, the special case...
+  if (sub_dimension == 0) {
+    assert( sub_index < num_nodes );
+    subentity_topo = MBVERTEX;
+    num_sub_entity_nodes = 1;
+    sub_entity_conn[0] = sub_index;
+    return;
+  }
+  
+  const int ho_bits = HasMidNodes( this_topo, num_nodes );
+  subentity_topo = SubEntityType(this_topo, sub_dimension, sub_index);
+  num_sub_entity_nodes = VerticesPerEntity(subentity_topo);
+  const short* corners = mConnectivityMap[this_topo][sub_dimension-1].conn[sub_index];
+  std::copy( corners, corners+num_sub_entity_nodes, sub_entity_conn );
+  
+  int sub_sub_corners[MB_MAX_SUB_ENTITY_VERTICES];
+  int side, sense, offset;
+  for (int dim = 1; dim <= sub_dimension; ++dim) {
+    if (!(ho_bits & (1<<dim)))
+      continue;
+    
+    const short num_mid = NumSubEntities( subentity_topo, dim );
+    for (int i = 0; i < num_mid; ++i) {
+      const MBEntityType sub_sub_topo = SubEntityType( subentity_topo, dim, i );
+      const int sub_sub_num_vert = VerticesPerEntity( sub_sub_topo );
+      SubEntityVertexIndices( subentity_topo, dim, i, sub_sub_corners );
+
+      for (int j = 0; j < sub_sub_num_vert; ++j)
+        sub_sub_corners[j] = corners[sub_sub_corners[j]];
+      SideNumber( this_topo, sub_sub_corners, sub_sub_num_vert, dim, side, sense, offset );
+      sub_entity_conn[num_sub_entity_nodes++] = HONodeIndex( this_topo, num_nodes, dim, side );
+    }
+  }
+}
+
   //! return the vertices of the specified sub entity
   //! \param parent_conn Connectivity of parent entity
   //! \param parent_type Entity type of parent entity

Modified: MOAB/branches/parallel_ghosting/MBCN.hpp
===================================================================
--- MOAB/branches/parallel_ghosting/MBCN.hpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/MBCN.hpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -38,6 +38,7 @@
 
 #include <vector>
 #include <algorithm>
+#include <cassert>
 
 #include "MBEntityType.h"
 
@@ -67,6 +68,11 @@
   
 public:
 
+  enum { MAX_NODES_PER_ELEMENT = 27 };
+  enum { MID_EDGE_BIT   = 1<<1,
+         MID_FACE_BIT   = 1<<2,
+         MID_REGION_BIT = 1<<3 };
+
     //! enum used to specify operation type
   enum {INTERSECT = 0, UNION};
 
@@ -82,7 +88,7 @@
     short int num_sub_elements;
     
       // Number of nodes in each sub-element of this dimension
-    short int num_nodes_per_sub_element[MB_MAX_SUB_ENTITIES];
+    short int num_corners_per_sub_element[MB_MAX_SUB_ENTITIES];
 
       // Type of each sub-element
     MBEntityType target_type[MB_MAX_SUB_ENTITIES];
@@ -93,7 +99,7 @@
 
     // mConnectivityMap[i=entity type][j=0,1,2]:
     //  num_sub_elements = # bounding edges(j=0) or faces(j=1) for entity type i, or self (j=2)
-    //  num_nodes_per_sub_element[k] (k=0..num_sub_elements-1) = number of nodes in sub-facet k
+    //  num_corners_per_sub_element[k] (k=0..num_sub_elements-1) = number of nodes in sub-facet k
     //    (can vary over sub-facets, e.g. faces bounding a pyramid) or self (j=2)
     //  target_type[k] = entity type of sub-facet k (e.g. MBTRI or MBQUAD bounding a pyramid) or self (j=2)
     //  conn[k][l] (l=0..MBCN::VerticesPerEntity[target_type[k]]) = vertex connectivity of sub-facet k,
@@ -118,6 +124,9 @@
     // (not documented with Doxygen)
   static const UpConnMap mUpConnMap[MBMAXTYPE][4][4];
 
+    // Mid-node bits indexed by number of nodes in element
+  static const unsigned char midNodesPerType[MBMAXTYPE][MAX_NODES_PER_ELEMENT+1];
+
     //! Permutation and reverse permutation vectors
   static short int permuteVec[MBMAXTYPE][3][MB_MAX_SUB_ENTITIES+1];
   static short int revPermuteVec[MBMAXTYPE][3][MB_MAX_SUB_ENTITIES+1];
@@ -167,6 +176,22 @@
                                      const int sub_index,
                                      int sub_entity_conn[]);
 
+  //! return the node indices of the specified sub-entity.
+  //! \param this_topo            The topology of the queried element type
+  //! \param num_nodes            The number of nodes in the queried element type.
+  //! \param sub_dimension        Dimension of sub-entity
+  //! \param sub_index            Index of sub-entity
+  //! \param sub_entity_topo      (Output) Topology of requested sub-entity.
+  //! \param num_sub_entity_nodes (Output) Number of nodes in the requested sub-entity.
+  //! \param sub_entity_conn      (Output) Connectivity of sub-entity
+  static void SubEntityNodeIndices(const MBEntityType this_topo, 
+                                   const int num_nodes,
+                                   const int sub_dimension,
+                                   const int sub_index,
+                                   MBEntityType& sub_entity_topo,
+                                   int& num_sub_entity_nodes,
+                                   int sub_entity_conn[]);
+
   //! return the vertices of the specified sub entity
   //! \param parent_conn Connectivity of parent entity
   //! \param parent_type Entity type of parent entity
@@ -376,6 +401,12 @@
                           const int num_verts, 
                           int mid_nodes[4]);
 
+  //! Same as above, except returns a single integer with the bits, from
+  //! least significant to most significant set to one if the corresponding
+  //! mid nodes on sub entities of the least dimension (0) to the highest
+  //! dimension (3) are present in the elment type.
+  static int HasMidNodes( const MBEntityType this_type, const int num_verts );
+
   //! given data about an element and a vertex in that element, return the dimension
   //! and index of the sub-entity that the vertex resolves.  If it does not resolve a
   //! sub-entity, either because it's a corner node or it's not in the element, -1 is
@@ -418,7 +449,7 @@
 
 inline short int MBCN::VerticesPerEntity(const MBEntityType t) 
 {
-  return (MBVERTEX == t ? 1 : mConnectivityMap[t][mConnectivityMap[t][0].topo_dimension-1].num_nodes_per_sub_element[0]);
+  return (MBVERTEX == t ? 1 : mConnectivityMap[t][mConnectivityMap[t][0].topo_dimension-1].num_corners_per_sub_element[0]);
 }
 
 inline short int MBCN::NumSubEntities(const MBEntityType t, const int d)
@@ -506,23 +537,20 @@
     return false;
 }
 
+inline int MBCN::HasMidNodes( const MBEntityType this_type, const int num_nodes )
+{
+  assert( (unsigned)num_nodes <= (unsigned)MAX_NODES_PER_ELEMENT );
+  return midNodesPerType[this_type][num_nodes];
+}
+
 inline void MBCN::HasMidNodes(const MBEntityType this_type, const int num_nodes,
                               int mid_nodes[4])
 {
-  if (num_nodes <= VerticesPerEntity(this_type) ||
-    // poly elements never have mid nodes as far as canonical ordering is concerned
-      MBPOLYGON == this_type || MBPOLYHEDRON == this_type) {
-    mid_nodes[0] = 0;
-    mid_nodes[1] = 0;
-    mid_nodes[2] = 0;
-    mid_nodes[3] = 0;
-    return;
-  }
-  
+  const int bits = HasMidNodes( this_type, num_nodes );
   mid_nodes[0] = 0;
-  mid_nodes[1] = HasMidEdgeNodes(this_type, num_nodes);
-  mid_nodes[2] = HasMidFaceNodes(this_type, num_nodes);
-  mid_nodes[3] = HasMidRegionNodes(this_type, num_nodes);
+  mid_nodes[1] = (bits & (1<<1)) >> 1;
+  mid_nodes[2] = (bits & (1<<2)) >> 2;
+  mid_nodes[3] = (bits & (1<<3)) >> 3;
 }
 
 //! Set permutation or reverse permutation vector

Modified: MOAB/branches/parallel_ghosting/MBCNArrays.hpp
===================================================================
--- MOAB/branches/parallel_ghosting/MBCNArrays.hpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/MBCNArrays.hpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -687,4 +687,33 @@
     {{{0}, {{0}} }, {{0}, {{0}} }, {{0}, {{0}} }} // source dim 3
   } // end type MBENTITYSET
 };
+const unsigned char E = MBCN::MID_EDGE_BIT;
+const unsigned char F = MBCN::MID_FACE_BIT;
+const unsigned char R = MBCN::MID_REGION_BIT;
+const unsigned char MBCN::midNodesPerType[MBMAXTYPE][MAX_NODES_PER_ELEMENT+1] = {
+// vertex
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// edge
+  { 0, 0, 0, E, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// tri
+  { 0, 0, 0, 0, F, 0, E, E|F, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// quad
+  { 0, 0, 0, 0, 0, F, 0, 0, E, E|F, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// polygon
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// tet 1, 2, 3, 4, 5, 6, 7, 8, 9,  10, 11, 12,13, 14,  15
+  { 0, 0, 0, 0, 0, R, 0, 0, F, F|R, E, E|R, 0, 0, E|F, E|F|R, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// pyramid   3, 4, 5, 6, 7, 8, 9,10, 11, 12,13, 14, 15,16,17, 18,  19
+  { 0, 0, 0, 0, 0, 0, R, 0, 0, 0, F, F|R, 0, E, E|R, 0, 0, 0, E|F, E|F|R, 0, 0, 0, 0, 0, 0, 0, 0 },
+// prism  2, 3, 4, 5, 6, 7, 8, 9,10,11, 12, 13,14,15, 16, 17,18,19, 20,  21
+  { 0, 0, 0, 0, 0, 0, 0, R, 0, 0, 0, F, F|R, 0, 0, E, E|R, 0, 0, 0, E|F, E|F|R, 0, 0, 0, 0, 0, 0 },
+// knife  2, 3, 4, 5, 6, 7, 8, 9,10,11, 12,13, 14,15,16,17, 18, 19,20,21, 22,  23
+  { 0, 0, 0, 0, 0, 0, 0, 0, R, 0, 0, 0, F, F|R, 0, 0, 0, E, E|R, 0, 0, 0, E|F, E|F|R, 0, 0, 0, 0 },
+// hex 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15, 16,17,18,19,20, 21, 22,23,24,25, 26,  27
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, R, 0, 0, 0, 0, F, F|R, 0, 0, 0, 0, E, E|R, 0, 0, 0, 0, E|F, E|F|R },
+// polyhedron
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+// set
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
+  };
 #endif

Modified: MOAB/branches/parallel_ghosting/MBSkinner.cpp
===================================================================
--- MOAB/branches/parallel_ghosting/MBSkinner.cpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/MBSkinner.cpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -280,7 +280,7 @@
   for(iter = source_entities.begin(); iter != end_iter; ++iter)
   {
     // get the connectivity of this entity
-    result = thisMB->get_connectivity(*iter, tmp_conn, num_nodes, true);
+    result = thisMB->get_connectivity(*iter, tmp_conn, num_nodes, false);
     if (MB_SUCCESS == result)
       conn = tmp_conn;
     else {
@@ -288,7 +288,7 @@
         // which doesn't store connectivity explicitly; use a connect
         // vector instead
       tmp_conn_vec.clear();
-      result = thisMB->get_connectivity(&(*iter), 1, tmp_conn_vec, true);
+      result = thisMB->get_connectivity(&(*iter), 1, tmp_conn_vec, false);
       if (MB_SUCCESS != result) return MB_FAILURE;
       conn = &tmp_conn_vec[0];
       num_nodes = tmp_conn_vec.size();
@@ -302,7 +302,7 @@
     const struct MBCN::ConnMap* conn_map = &(MBCN::mConnectivityMap[type][mTargetDim-1]);
     for(int i=0; i<conn_map->num_sub_elements; i++)
     {
-      int num_sub_nodes = conn_map->num_nodes_per_sub_element[i];
+      int num_sub_nodes = conn_map->num_corners_per_sub_element[i];
       assert(num_sub_nodes <= 32);
       for(int j=0; j<num_sub_nodes; j++)
         sub_conn[j] = conn[conn_map->conn[i][j]];
@@ -312,23 +312,36 @@
         result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, source_dim, false,
                                          dum_elems);
         if (MB_SUCCESS != result) return result;
-        assert(0 < dum_elems.size() && 3 > dum_elems.size());
-        if (1 == dum_elems.size() ||
-            (2 == dum_elems.size() && 
-             (source_entities.find(*dum_elems.begin()) == source_entities.end() ||
-              source_entities.find(*dum_elems.rbegin()) == source_entities.end()))) {
+        dum_elems = dum_elems.intersect( source_entities );
+        if (dum_elems.empty()) {
+          assert(false);  // should never happen
+          return MB_FAILURE;
+        }
+        // if (dum_elems.size() > 2) { 
+          // source entities do not form a valid source-dim patch (t-junction).
+          // do we care?
+        // }
+        
+        if (1 == dum_elems.size()) {
             // this sub_element is on the skin
 
             // check for existing entity
           dum_sub_elems.clear();
-          result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, source_dim-1, false,
+          result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, mTargetDim, false,
                                            dum_sub_elems);
           if (MB_SUCCESS != result) return result;
           if (dum_sub_elems.empty()) {
               // need to create one
             MBEntityHandle tmphndl=0;
-            result = thisMB->create_element(conn_map->target_type[i], sub_conn, num_sub_nodes,
-                                            tmphndl);
+            int indices[MB_MAX_SUB_ENTITY_VERTICES];
+            MBEntityType new_type;
+            int num_new_nodes;
+            MBCN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
+            for(int j=0; j<num_new_nodes; j++)
+              sub_conn[j] = conn[indices[j]];
+        
+            result = thisMB->create_element(new_type, sub_conn,  
+                                            num_new_nodes, tmphndl);
             forward_target_entities.insert(tmphndl);
           }
           else {
@@ -358,7 +371,13 @@
         if(match == 0)
         {
           MBEntityHandle tmphndl=0;
-          result = thisMB->create_element(conn_map->target_type[i], sub_conn, num_sub_nodes,
+          int indices[MB_MAX_SUB_ENTITY_VERTICES];
+          MBEntityType new_type;
+          int num_new_nodes;
+          MBCN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
+          for(int j=0; j<num_new_nodes; j++)
+            sub_conn[j] = conn[indices[j]];
+          result = thisMB->create_element(new_type, sub_conn, num_new_nodes,
                                           tmphndl);
           assert(MB_SUCCESS == result);
           add_adjacency(tmphndl, sub_conn, num_sub_nodes);
@@ -652,7 +671,7 @@
     const struct MBCN::ConnMap* conn_map = &(MBCN::mConnectivityMap[type][0]);
     for(int i=0; i<conn_map->num_sub_elements; i++)
     {
-      int num_sub_nodes = conn_map->num_nodes_per_sub_element[i];
+      int num_sub_nodes = conn_map->num_corners_per_sub_element[i];
       assert(num_sub_nodes <= 32);
       for(int j=0; j<num_sub_nodes; j++)
         sub_conn[j] = conn[conn_map->conn[i][j]];

Modified: MOAB/branches/parallel_ghosting/parallel/pcomm_unit.cpp
===================================================================
--- MOAB/branches/parallel_ghosting/parallel/pcomm_unit.cpp	2009-06-22 18:07:19 UTC (rev 2960)
+++ MOAB/branches/parallel_ghosting/parallel/pcomm_unit.cpp	2009-06-23 14:25:58 UTC (rev 2961)
@@ -7,6 +7,7 @@
 #include <algorithm>
 #include <vector>
 #include <set>
+#include <sstream>
 
 #ifdef USE_MPI
 #  include <mpi.h>
@@ -224,66 +225,6 @@
   delete [] elems;
 }
 
-MBErrorCode set_owners(unsigned char pstatus,
-                       MBParallelComm *pc0, MBEntityHandle ent0, 
-                       MBParallelComm *pc1, MBEntityHandle ent1, 
-                       MBParallelComm *pc2 = NULL, MBEntityHandle ent2 = 0, 
-                       MBParallelComm *pc3 = NULL, MBEntityHandle ent3 = 0)
-{
-  int owners[MAX_SHARING_PROCS];
-  MBEntityHandle tmp_handles[MAX_SHARING_PROCS];
-  
-  MBErrorCode result = MB_SUCCESS;
-  
-  std::fill(owners, owners+MAX_SHARING_PROCS, -1);
-  std::fill(tmp_handles, tmp_handles+MAX_SHARING_PROCS, 0);
-  owners[0] = pc0->proc_config().proc_rank(); owners[1] = pc1->proc_config().proc_rank();
-  tmp_handles[0] = ent0; tmp_handles[1] = ent1;
-  int np = 2;
-  if (pc2) {
-    owners[2] = pc2->proc_config().proc_rank();
-    tmp_handles[2] = ent2;
-    np++;
-  }
-  if (pc3) {
-    owners[3] = pc3->proc_config().proc_rank();
-    tmp_handles[3] = ent3;
-    np++;
-  }
-  if (np > 2) {
-    result = pc0->get_moab()->tag_set_data(pc0->sharedps_tag(), &ent0, 1, owners);
-    result = pc0->get_moab()->tag_set_data(pc0->sharedhs_tag(), &ent0, 1, tmp_handles);
-    result = pc0->get_moab()->tag_set_data(pc0->pstatus_tag(), &ent0, 1, &pstatus);
-      // 2nd and greater don't own them
-    pstatus |= PSTATUS_NOT_OWNED;
-    result = pc1->get_moab()->tag_set_data(pc1->sharedps_tag(), &ent1, 1, owners);
-    result = pc1->get_moab()->tag_set_data(pc1->sharedhs_tag(), &ent1, 1, tmp_handles);
-    result = pc1->get_moab()->tag_set_data(pc1->pstatus_tag(), &ent1, 1, &pstatus);
-    if (pc2) {
-      result = pc2->get_moab()->tag_set_data(pc2->sharedps_tag(), &ent2, 1, owners);
-      result = pc2->get_moab()->tag_set_data(pc2->sharedhs_tag(), &ent2, 1, tmp_handles);
-      result = pc2->get_moab()->tag_set_data(pc2->pstatus_tag(), &ent2, 1, &pstatus);
-    }
-    if (pc3) {
-      result = pc3->get_moab()->tag_set_data(pc3->sharedps_tag(), &ent3, 1, owners);
-      result = pc3->get_moab()->tag_set_data(pc3->sharedhs_tag(), &ent3, 1, tmp_handles);
-      result = pc3->get_moab()->tag_set_data(pc3->pstatus_tag(), &ent3, 1, &pstatus);
-    }
-  }
-  else {
-    result = pc0->get_moab()->tag_set_data(pc0->sharedp_tag(), &ent0, 1, &owners[1]);
-    result = pc0->get_moab()->tag_set_data(pc0->sharedh_tag(), &ent0, 1, &tmp_handles[1]);
-    result = pc0->get_moab()->tag_set_data(pc0->pstatus_tag(), &ent0, 1, &pstatus);
-      // 2nd and greater don't own them
-    pstatus |= PSTATUS_NOT_OWNED;
-    result = pc1->get_moab()->tag_set_data(pc1->sharedp_tag(), &ent1, 1, owners);
-    result = pc1->get_moab()->tag_set_data(pc1->sharedh_tag(), &ent1, 1, tmp_handles);
-    result = pc1->get_moab()->tag_set_data(pc1->pstatus_tag(), &ent1, 1, &pstatus);
-  }
-  
-  return result;
-}
-
 MBErrorCode create_patch(MBInterface *moab, MBRange &verts, MBRange &quads,
                          unsigned int n, double *xyz, int *gids) 
 {
@@ -320,80 +261,6 @@
   return result;
 }
 
-MBErrorCode set_owners(MBParallelComm **pc, MBRange *verts, MBRange *quads,
-                       const bool edges_too) 
-{
-    // P0-P2
-  unsigned char pstat = PSTATUS_SHARED | PSTATUS_INTERFACE;
-  MBErrorCode rval = set_owners(pstat, pc[0], verts[0][6], pc[2], verts[2][0]);
-  rval = set_owners(pstat, pc[0], verts[0][7], pc[2], verts[2][1]); CHECK_ERR(rval);
-
-    // P1-P2
-  rval = set_owners(pstat, pc[1], verts[1][7], pc[2], verts[2][5]); CHECK_ERR(rval);
-  rval = set_owners(pstat, pc[1], verts[1][8], pc[2], verts[2][8]); CHECK_ERR(rval);
-    // P0-P1-P3
-  pstat |= PSTATUS_MULTISHARED;
-  rval = set_owners(pstat, pc[0], verts[0][2], pc[1], verts[1][0], pc[3], verts[3][0]); CHECK_ERR(rval);
-  rval = set_owners(pstat, pc[0], verts[0][5], pc[1], verts[1][3], pc[3], verts[3][3]); CHECK_ERR(rval);
-    // P0-P1-P2-P3
-  rval = set_owners(pstat, pc[0], verts[0][8], pc[1], verts[1][6], 
-                    pc[2], verts[2][2], pc[3], verts[3][6]); CHECK_ERR(rval);
-
-  if (edges_too) {
-    MeshTopoUtil *mtu[4];
-    MBInterface *mb[4];
-    MBRange dum_range;
-      // create mtu's and explicit edges
-    for (unsigned int i = 0; i < 4; i++) {
-      mb[i] = pc[i]->get_moab();
-      assert(mb[i]);
-      mtu[i] = new MeshTopoUtil(mb[i]);
-      rval = mb[i]->get_adjacencies(quads[i], 1, true, dum_range, 
-                                    MBInterface::UNION);
-      CHECK_ERR(rval);
-      dum_range.clear();
-    }
-    
-    MBEntityHandle edge1, edge2, edge3;
-    pstat = PSTATUS_SHARED | PSTATUS_INTERFACE;
-      // P0-P2
-    edge1 = mtu[0]->common_entity(verts[0][6], verts[0][7], 1);
-    edge2 = mtu[2]->common_entity(verts[2][0], verts[2][1], 1);
-    assert(edge1 && edge2);
-    rval = set_owners(pstat, pc[0], edge1, pc[2], edge2); CHECK_ERR(rval);
-    edge1 = mtu[0]->common_entity(verts[0][7], verts[0][8], 1);
-    edge2 = mtu[2]->common_entity(verts[2][1], verts[2][2], 1);
-    assert(edge1 && edge2);
-    rval = set_owners(pstat, pc[0], edge1, pc[2], edge2); CHECK_ERR(rval);
-      // P1-P2
-    edge1 = mtu[1]->common_entity(verts[1][6], verts[1][7], 1);
-    edge2 = mtu[2]->common_entity(verts[2][2], verts[2][5], 1);
-    assert(edge1 && edge2);
-    rval = set_owners(pstat, pc[1], edge1, pc[2], edge2); CHECK_ERR(rval);
-    edge1 = mtu[1]->common_entity(verts[1][7], verts[1][8], 1);
-    edge2 = mtu[2]->common_entity(verts[2][5], verts[2][8], 1);
-    assert(edge1 && edge2);
-    rval = set_owners(pstat, pc[1], edge1, pc[2], edge2); CHECK_ERR(rval);
-      // P0-P1-P3
-    pstat |= PSTATUS_MULTISHARED;
-    edge1 = mtu[0]->common_entity(verts[0][2], verts[0][5], 1);
-    edge2 = mtu[1]->common_entity(verts[1][0], verts[1][3], 1);
-    edge3 = mtu[3]->common_entity(verts[3][0], verts[3][3], 1);
-    assert(edge1 && edge2 && edge3);
-    rval = set_owners(pstat, pc[0], edge1, pc[1], edge2, pc[3], edge3); CHECK_ERR(rval);
-    edge1 = mtu[0]->common_entity(verts[0][5], verts[0][8], 1);
-    edge2 = mtu[1]->common_entity(verts[1][3], verts[1][6], 1);
-    edge3 = mtu[3]->common_entity(verts[3][3], verts[3][6], 1);
-    assert(edge1 && edge2 && edge3);
-    rval = set_owners(pstat, pc[0], edge1, pc[1], edge2, pc[3], edge3); CHECK_ERR(rval);
-
-    for (unsigned int i = 0; i < 4; i++)
-      delete mtu[i];
-  }
-    
-  return MB_SUCCESS;
-}
-
 MBErrorCode create_shared_grid_2d(MBParallelComm **pc, MBRange *verts, MBRange *quads) 
 {
 //          
@@ -456,15 +323,17 @@
 //    
 //   4 _____   _____ 
 //    |  |  | |  |  |
-//   3|__|__| |__|__|
+//   3|__|__| |__|__|          GIDS               HANDLES
 //    |P1|  | |P2|  |
-//    |__|__| |__|__|
-//    2 ___________         P3 - k = 2..4
-//     |  |  |  |  | 
+//    |__|__| |__|__|   20 21 22   22 23 24      7  8  9    7  8  9
+//    2 ___________     15 16 17   17 18 19      4  5  6    4  5  6
+//     |  |  |  |  |    10 11 12   12 13 14      1  2  3    1  2  3
 // /  1|__|__|__|__|         
-// |   |  |P0|  |  |         
-// J  0|__|__|__|__|        
-// I-> 0  1  2  3  4
+// |   |  |P0|  |  |     10 11 12 13 14           10 11 12 13 14     
+// J  0|__|__|__|__|      5  6  7  8  9            5  6  7  8  9
+// I-> 0  1  2  3  4      0  1  2  3  4            0  1  2  3  4
+//
+//  P3 - k = 2..4
   
     // create structured meshes
     // ijkmin[p][ijk], ijkmax[p][ijk]
@@ -474,7 +343,9 @@
 
   int nijk[P][3];
   int NIJK[3] = {0, 0, 0};
-#define INDEX(i, j, k) (k * NIJK[1] * NIJK[0] + j * NIJK[0] + i)
+#define INDEXG(i, j, k) (k * NIJK[1] * NIJK[0] + j * NIJK[0] + i)
+#define INDEXL(i, j, k) ((k-ijkmin[p][2])*nijk[p][1]*nijk[p][0] + \
+                         (j-ijkmin[p][1])*nijk[p][0] + (i - ijkmin[p][0]))
 
   int p, i, j, k;
   for (int p = 0; p < P; p++) {
@@ -489,7 +360,6 @@
   MBErrorCode rval;
   MBTag gid_tag;
   int dum_default = -1;
-  if (MB_SUCCESS != rval) return rval;
 
   for (p = 0; p < P; p++) {
     rval = pc[p]->get_moab()->tag_create(GLOBAL_ID_TAG_NAME, 
@@ -502,44 +372,56 @@
     int nverts = nijk[p][0] * nijk[p][1] * nijk[p][2];
     xyz.resize(3*nverts);
     gids.resize(nverts);
-    rval = pc[p]->get_moab()->create_vertices(&xyz[0], nverts, verts[p]);
-    CHECK_ERR(rval);
 
       // set vertex gids
     int nv = 0;
-    for (k = ijkmin[p][2]; k < ijkmax[p][2]; k++) 
-      for (j = ijkmin[p][1]; j < ijkmax[p][1]; j++) 
-        for (i = ijkmin[p][0]; i < ijkmax[p][0]; i++)
+    for (k = ijkmin[p][2]; k <= ijkmax[p][2]; k++) 
+      for (j = ijkmin[p][1]; j <= ijkmax[p][1]; j++) 
+        for (i = ijkmin[p][0]; i <= ijkmax[p][0]; i++) {
+            // xyz
+          xyz[3*nv] = i;
+          xyz[3*nv+1] = j;
+          xyz[3*nv+2] = k;
+          
             // gid
-          gids[nv++] = INDEX(i, j, k);
+          gids[nv++] = INDEXG(i, j, k);
+        }
+    
 
+    rval = pc[p]->get_moab()->create_vertices(&xyz[0], nverts, verts[p]);
+    CHECK_ERR(rval);
+
     rval = pc[p]->get_moab()->tag_set_data(gid_tag, verts[p], &gids[0]);
     if (MB_SUCCESS != rval) return rval;
 
       // make elements
     nv = 0;
     MBEntityHandle connect[8], dum_hex;
-    for (k = ijkmin[p][2]; k < ijkmax[p][2]-1; k++) 
-      for (j = ijkmin[p][1]; j < ijkmax[p][1]-1; j++) 
-        for (i = ijkmin[p][0]; i < ijkmax[p][0]-1; i++) {
+    for (k = ijkmin[p][2]; k < ijkmax[p][2]; k++) 
+      for (j = ijkmin[p][1]; j < ijkmax[p][1]; j++) 
+        for (i = ijkmin[p][0]; i < ijkmax[p][0]; i++) {
             // gid
-          connect[0] = verts[p][INDEX(i, j, k)];
-          connect[1] = verts[p][INDEX(i+1, j, k)];
-          connect[2] = verts[p][INDEX(i+1, j+1, k)];
-          connect[3] = verts[p][INDEX(i, j+1, k)];
-          connect[4] = verts[p][INDEX(i, j, k+1)];
-          connect[5] = verts[p][INDEX(i+1, j, k+1)];
-          connect[6] = verts[p][INDEX(i+1, j+1, k+1)];
-          connect[7] = verts[p][INDEX(i, j+1, k+1)];
+          connect[0] = verts[p][INDEXL(i, j, k)];
+          connect[1] = verts[p][INDEXL(i+1, j, k)];
+          connect[2] = verts[p][INDEXL(i+1, j+1, k)];
+          connect[3] = verts[p][INDEXL(i, j+1, k)];
+          connect[4] = verts[p][INDEXL(i, j, k+1)];
+          connect[5] = verts[p][INDEXL(i+1, j, k+1)];
+          connect[6] = verts[p][INDEXL(i+1, j+1, k+1)];
+          connect[7] = verts[p][INDEXL(i, j+1, k+1)];
           rval = pc[p]->get_moab()->create_element(MBHEX, connect, 8, dum_hex);
           hexes[p].insert(dum_hex);
-          gids[nv++] = INDEX(i, j, k);
+          gids[nv++] = INDEXG(i, j, k);
         }
     rval = pc[p]->get_moab()->tag_set_data(gid_tag, hexes[p], &gids[0]);
     if (MB_SUCCESS != rval) return rval;
+
+    std::ostringstream fname;
+    fname << "tmp" << p << ".h5m";
+    rval = pc[p]->get_moab()->write_file(fname.str().c_str());
+    if (MB_SUCCESS != rval) return rval;
   }
-  
-  rval = MBParallelComm::resolve_shared_ents(pc, 3, 3);
+  rval = MBParallelComm::resolve_shared_ents(pc, 4, 3);
   CHECK_ERR(rval);
   return rval;
 }



More information about the moab-dev mailing list