[MOAB-dev] commit/MOAB: 15 new changesets

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Mon Jul 14 17:56:06 CDT 2014


15 new commits in MOAB:

https://bitbucket.org/fathomteam/moab/commits/a6e546d1a491/
Changeset:   a6e546d1a491
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:57
Summary:     add large mesh example

Affected #:  2 files

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
new file mode 100644
index 0000000..6ef1dc9
--- /dev/null
+++ b/examples/GenLargeMesh.cpp
@@ -0,0 +1,105 @@
+/** @example GenLargeMesh.cpp \n
+ * \brief Create a large structured mesh, partitioned \n
+ * <b>To run</b>: mpiexec -np 2 GenLargeMesh \n
+ *
+ *  It shows how to load create a mesh on the fly, on multiple
+ *  processors, block wise
+ *  Each processor will create its version of a block mesh, partitioned
+ *  as AxBxC blocks. Each block will be with blockSize^3 hexahedrons, and will get a
+ *  different PARALLEL_PARTITION tag
+ *
+ *  The number of tasks will be MxNxK, and it must match the mpi size
+ *  Each task will generate its mesh at location (m,n,k)
+ *
+ *  By default M=2, N=1, K=1, so by default it should be launched on 2 procs
+ *  By default, blockSize is 4, and A=2, B=2, C=2, so each task will generate locally
+ *    blockSize^3 x A x B x C hexahedrons (value = 64x8 = 512 hexas, in 8 partitions)
+ *
+ *  The total number f partitions will be A*B*C*M*N*K (default 16)
+ *
+ *  Each part in partition will get a proper tag
+ *
+ *  The vertices will get a proper global id, which will be used to resolve the
+ *  shared entities
+ *  The output will be written in parallel, and we will try sizes as big as we can
+ *  (up to a billion vertices, because we use int for global ids)
+ *
+ *  Within each partition, the hexas will be numbered contiguously, and also the
+ *  vertices; The global id will be determined easily, for a vertex, but the entity
+ *  handle space will be more interesting to handle, within a partition (we want
+ *  contiguous handles within a partition)
+ *
+ *  We may or may not use ScdInterface, because we want control over global id and
+ *  entity handle within a partition
+ *
+ *
+ *  mpiexec -np 2 GenLargeMesh
+ */
+
+#include "moab/Core.hpp"
+#include "moab/ProgOptions.hpp"
+#include "moab/ParallelComm.hpp"
+#include "moab/CN.hpp"
+
+#include <iostream>
+#include <vector>
+
+using namespace moab;
+using namespace std;
+int main(int argc, char **argv)
+{
+  int A=2, B=2, C=2, M=2, N=1, K=1;
+  int blockSize = 4;
+
+  MPI_Init(&argc, &argv);
+
+  ProgOptions opts;
+
+  opts.addOpt<int>(std::string("blockSize,b"),
+      std::string("Block size of mesh (default=4)"), &blockSize);
+  opts.addOpt<int>(std::string("xproc,M"),
+      std::string("Number of processors in x dir (default=2)"), &M);
+  opts.addOpt<int>(std::string("yproc,N"),
+      std::string("Number of processors in y dir (default=1)"), &N);
+  opts.addOpt<int>(std::string("zproc,K"),
+      std::string("Number of processors in z dir (default=1)"), &K);
+
+  opts.addOpt<int>(std::string("xblocks,A"),
+      std::string("Number of blocks on a task in x dir (default=2)"), &A);
+  opts.addOpt<int>(std::string("yblocks,B"),
+      std::string("Number of blocks on a task in y dir (default=2)"), &B);
+  opts.addOpt<int>(std::string("zblocks,C"),
+      std::string("Number of blocks on a task in x dir (default=2)"), &C);
+
+
+  opts.parseCommandLine(argc, argv);
+
+  Interface* mb = new Core;
+  if (NULL == mb)
+  {
+    MPI_Finalize();
+    return 1;
+  }
+
+  MPI_Comm comm ;
+  int rank, size;
+  MPI_Comm_rank( MPI_COMM_WORLD, &rank );
+  MPI_Comm_rank( MPI_COMM_WORLD, &size );
+
+  if (M*N*K != size)
+  {
+    if (rank==0)
+    {
+      cout <<"M*N*K=" << M*N*K << "  != size = "<< " size\n";
+    }
+    MPI_Finalize();
+    return 1;
+  }
+
+  // the global id of each vertex will be, on task
+
+  MPI_Finalize();
+  return 0;
+}
+
+

diff --git a/examples/makefile b/examples/makefile
index 787c8d3..f41def7 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -71,6 +71,9 @@ VisTags: VisTags.o ${MOAB_LIBDIR}/libMOAB.la
 ReadWriteTest: ReadWriteTest.o ${MOAB_LIBDIR}/libMOAB.la
 	${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} 
 	
+GenLargeMesh: GenLargeMesh.o ${MOAB_LIBDIR}/libMOAB.la
+	${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} 
+	
 clean:
 	rm -rf *.o *.mod *.h5m ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES} ${F90EXAMPLES}
 


https://bitbucket.org/fathomteam/moab/commits/f6bbe84f4c2e/
Changeset:   f6bbe84f4c2e
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:57
Summary:     is_valid should be on the interface too

not only on the Core.hpp
How do we know if an entity handle is valid or not?
is there another way?

Affected #:  1 file

diff --git a/src/moab/Interface.hpp b/src/moab/Interface.hpp
index 09f4219..c8102b8 100644
--- a/src/moab/Interface.hpp
+++ b/src/moab/Interface.hpp
@@ -1981,6 +1981,8 @@ public:
                                         SetIterator *&set_iter) = 0;
     /**@}*/
 
+
+  virtual bool is_valid(EntityHandle entity) const = 0;
 };
 
 //! predicate for STL algorithms.  Returns true if the entity handle is


https://bitbucket.org/fathomteam/moab/commits/acb938397bc9/
Changeset:   acb938397bc9
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:57
Summary:     validate sharedEnts

it is becoming important if some entities are deleted after they are resolved
they are part of the sharedEnts vector, but they are not deleted
from there.
This is the advantage of keeping entities in lists in member data
we should hunt them down :(

Affected #:  2 files

diff --git a/src/parallel/ParallelComm.cpp b/src/parallel/ParallelComm.cpp
index 886ef77..4fb8afc 100644
--- a/src/parallel/ParallelComm.cpp
+++ b/src/parallel/ParallelComm.cpp
@@ -8791,7 +8791,22 @@ ErrorCode ParallelComm::settle_intersection_points(Range & edges, Range & shared
 
   return MB_SUCCESS;
   // end copy
+}
+ErrorCode ParallelComm::correct_shared_entities()
+{
+
+  std::vector<EntityHandle> good_ents;
+  for (size_t i=0; i<sharedEnts.size(); i++)
+  {
+    if (mbImpl->is_valid(sharedEnts[i]))
+    {
+      good_ents.push_back(sharedEnts[i]);
+    }
   }
+  sharedEnts = good_ents;
+
+  return MB_SUCCESS;
+}
 
 void ParallelComm::print_pstatus(unsigned char pstat, std::string &ostr) 
 {

diff --git a/src/parallel/moab/ParallelComm.hpp b/src/parallel/moab/ParallelComm.hpp
index 3424c51..2a16ac2 100644
--- a/src/parallel/moab/ParallelComm.hpp
+++ b/src/parallel/moab/ParallelComm.hpp
@@ -932,6 +932,12 @@ namespace moab {
     ErrorCode settle_intersection_points(Range & edges, Range & shared_edges_owned,
         std::vector<std::vector<EntityHandle> *> & extraNodesVec, double tolerance);
 
+    /* \brief check the shared entities if they exist anymore
+     * will check the shared ents array, and clean it if necessary
+     *
+     */
+    ErrorCode correct_shared_entities();
+
   private:
 
     ErrorCode reduce_void(int tag_data_type, const MPI_Op mpi_op, int num_ents, void *old_vals, void *new_vals);


https://bitbucket.org/fathomteam/moab/commits/834ce0f7daa2/
Changeset:   834ce0f7daa2
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:57
Summary:     improve the large mesh example

this example pushes the limits of moab
It will be used later to create meshes larger than 1 billion vertices

Each processor creates a structure of AxBxC blocks, each block with
blockSize ^ 3 hexas
The example should be run on MxNxK array of processors.
After resolving and sharing the vertices, mesh is written in
parallel.

each block will represent a partition
It is different than test paralle_write_test because there the number
of partitions created is equal to the number of processors.
Here, each processor creates AxBxC partitions
For a total of AxBxC x MxNxK partitions

block size, A, B, .., K are options, with default values

The total number of hexas will be

 AxBxC x MxNxK x block^3

edges and faces are removed after resolve sharing

the file will have just vertices and hexas, and it will be already
partitioned

We will add some tags of interest, to make it better.

The writer failed if I just deleted the edges and faces.
I had to introduce a method to correct sharedEnts in PC

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 6ef1dc9..1e0fa73 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -40,12 +40,16 @@
 #include "moab/ProgOptions.hpp"
 #include "moab/ParallelComm.hpp"
 #include "moab/CN.hpp"
+#include "moab/ReadUtilIface.hpp"
+#include "moab/MergeMesh.hpp"
 
 #include <iostream>
 #include <vector>
 
 using namespace moab;
 using namespace std;
+
+#define CHECKE(message) if(rval != MB_SUCCESS) { cout << (message) << "\n"; MPI_Finalize(); return 1; }
 int main(int argc, char **argv)
 {
   int A=2, B=2, C=2, M=2, N=1, K=1;
@@ -58,7 +62,7 @@ int main(int argc, char **argv)
   opts.addOpt<int>(std::string("blockSize,b"),
       std::string("Block size of mesh (default=4)"), &blockSize);
   opts.addOpt<int>(std::string("xproc,M"),
-      std::string("Number of processors in x dir (default=2)"), &M);
+      std::string("Number of processors in x dir (default=1)"), &M);
   opts.addOpt<int>(std::string("yproc,N"),
       std::string("Number of processors in y dir (default=1)"), &N);
   opts.addOpt<int>(std::string("zproc,K"),
@@ -80,23 +84,175 @@ int main(int argc, char **argv)
     MPI_Finalize();
     return 1;
   }
+  ReadUtilIface *iface;
+  ErrorCode rval = mb->query_interface(iface);
+  CHECKE("Can't get reader interface");
 
-  MPI_Comm comm ;
   int rank, size;
   MPI_Comm_rank( MPI_COMM_WORLD, &rank );
-  MPI_Comm_rank( MPI_COMM_WORLD, &size );
+  MPI_Comm_size( MPI_COMM_WORLD, &size );
 
   if (M*N*K != size)
   {
     if (rank==0)
     {
-      cout <<"M*N*K=" << M*N*K << "  != size = "<< " size\n";
+      cout <<"M*N*K=" << M*N*K << "  != size = "<< size << "\n";
     }
     MPI_Finalize();
     return 1;
   }
 
-  // the global id of each vertex will be, on task
+  // determine m, n, k for processor rank
+  int m,n,k;
+  k = rank/(M*N);
+  int leftover = rank%(M*N);
+  n = leftover/M;
+  m = leftover%M;
+  if (rank==size-1)
+  {
+    cout << "m, n, k for last rank:" << m << " " << n << " " << k << "\n";
+  }
+
+  // so there are a total of M * A * blockSize elements in x direction (so M * A * blockSize + 1 verts in x direction)
+  // so there are a total of N * B * blockSize elements in y direction (so N * B * blockSize + 1 verts in y direction)
+  // so there are a total of K * C * blockSize elements in z direction (so K * C * blockSize + 1 verts in z direction)
+
+  // there are ( M * A blockSize )      *  ( N * B * blockSize)      * (K * C * blockSize )    hexas
+  // there are ( M * A * blockSize + 1) *  ( N * B * blockSize + 1 ) * (K * C * blockSize + 1) vertices
+  // x is the first dimension that varies
+
+  int NX = ( M * A * blockSize + 1);
+  int NY = ( N * B * blockSize + 1);
+  int NZ = ( K * C * blockSize + 1);
+  int blockSize1 = blockSize + 1;// used for vertices
+  int bsq = blockSize * blockSize;
+  int b1sq = blockSize1 * blockSize1;
+  // generate the block at (a, b, c); it will represent a partition , it will get a partition tag
+
+  Tag global_id_tag;
+  mb->tag_get_handle("GLOBAL_ID", 1, MB_TYPE_INTEGER,
+               global_id_tag);
+
+  Tag part_tag;
+  int dum_id=-1;
+  mb->tag_get_handle("PARALLEL_PARTITION", 1, MB_TYPE_INTEGER,
+                 part_tag, MB_TAG_CREAT|MB_TAG_SPARSE, &dum_id);
+
+  for (int a=0; a<A; a++)
+  {
+    for (int b=0; b<B; b++)
+    {
+      for (int c=0; c<C; c++)
+      {
+        // we will generate (block+1)^3 vertices, and block^3 hexas
+        // the global id of the vertices will come from m, n, k, a, b, c
+        // x will vary from  m*A*block + a*block to m*A*block+(a+1)*block etc;
+        int num_nodes = (blockSize+1)*(blockSize+1)*(blockSize+1);
+
+        vector<double*> arrays;
+        EntityHandle startv;
+        rval = iface->get_node_coords(3, num_nodes, 0, startv, arrays);
+        CHECKE("Can't get node coords.");
+
+        int x = m*A*blockSize + a*blockSize;
+        int y = n*B*blockSize + b*blockSize;
+        int z = k*C*blockSize + c*blockSize;
+        int ix=0;
+        vector<int>  gids(num_nodes);
+        Range verts(startv, startv+num_nodes-1);
+        for (int kk=0; kk<blockSize+1; kk++)
+        {
+          for (int jj=0; jj<blockSize+1; jj++)
+          {
+            for (int ii=0; ii<blockSize+1; ii++)
+            {
+              arrays[0][ix] = (double)x+ii;
+              arrays[1][ix] = (double)y+jj;
+              arrays[2][ix] = (double)z+kk;
+              gids[ix] = 1 + (z+kk) * (NX*NY) + (y+jj) * NX + (x+ii);
+              ix++;
+            }
+          }
+        }
+        mb->tag_set_data(global_id_tag, verts, &gids[0]);
+        int num_hexas = (blockSize)*(blockSize)*(blockSize);
+        EntityHandle starte; // connectivity
+        EntityHandle * conn;
+        rval = iface->get_element_connect(num_hexas, 8, MBHEX, 0, starte, conn);
+        CHECKE("Can't get hexa  connectivity.");
+
+        Range hexas(starte, starte+num_hexas-1);
+        // fill  hexas
+        ix=0;
+        for (int kk=0; kk<blockSize; kk++)
+        {
+          for (int jj=0; jj<blockSize; jj++)
+          {
+            for (int ii=0; ii<blockSize; ii++)
+            {
+              EntityHandle corner=startv + ii + (jj*blockSize1) + kk * b1sq;
+              conn[ix]=corner;
+              conn[ix+1]=corner+1;
+              conn[ix+2]=corner+ 1 + blockSize1;
+              conn[ix+3]=corner+ blockSize1;
+              conn[ix+4]=corner + b1sq;
+              conn[ix+5]=corner+1 + b1sq;
+              conn[ix+6]=corner+ 1 + blockSize1 + b1sq;
+              conn[ix+7]=corner+ blockSize1 + b1sq;
+              ix+=8;
+            }
+          }
+        }
+        EntityHandle part_set;
+        rval = mb->create_meshset(MESHSET_SET, part_set);
+        CHECKE("Can't create mesh set.");
+        rval = mb->add_entities(part_set, hexas);
+        CHECKE("Can't add entities to set.");
+        int part_num= a +  m*A + (b + n*B)*(M*A) + (c+k*C)*(M*A * N*B);
+        rval = mb->tag_set_data(part_tag, &part_set, 1, &part_num);
+        CHECKE("Can't set part tag on set");
+
+
+
+      }
+    }
+  }
+  // after the mesh is generated on each proc, merge the vertices
+  MergeMesh mm(mb);
+  Range allhexas;
+  rval = mb->get_entities_by_type(0, MBHEX, allhexas);
+
+  CHECKE("Can't get all hexa  elements.");
+
+  rval = mm.merge_entities( allhexas, 0.0001);
+  CHECKE("Can't merge");
+  ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0);
+  if (pcomm==NULL)
+  {
+    pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
+  }
+  rval = pcomm->resolve_shared_ents( 0, allhexas, 3, 0 );
+  CHECKE("Can't resolve shared ents");
+
+  // delete all quads and edges
+  Range toDelete;
+  rval =  mb->get_entities_by_dimension(0, 1, toDelete);
+  CHECKE("Can't get edges");
+  rval = mb->delete_entities(toDelete);
+  CHECKE("Can't delete edges");
+
+  toDelete.clear();
+
+  rval = mb->get_entities_by_dimension(0, 2, toDelete);
+  CHECKE("Can't get faces");
+  rval = mb->delete_entities(toDelete);
+  CHECKE("Can't delete faces");
+
+  rval = pcomm->correct_shared_entities() ;
+  CHECKE("Can't correct local shared")
+
+  rval = mb->write_file("test1.h5m", 0, ";;PARALLEL=WRITE_PART");
+  CHECKE("Can't write in parallel");
 
   MPI_Finalize();
   return 0;


https://bitbucket.org/fathomteam/moab/commits/f916f2f0c3fd/
Changeset:   f916f2f0c3fd
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:57
Summary:     changes to the large structured mesh example

fun project
see how big of a model we can write in moab
use a different merging strategy, in which the global id tag is
used for determining what should be merged. (instead of a
adaptive KD tree)
It took 21 minutes on a task to merge a mesh with 500k vertices
This is ridiculous; the new method takes 3 minutes
I think that there is still room for improvement.

In the new method, we get the tag on each vertex and then order
The vertices that get the same tag will be merged (duh..)
It is away to generate a big model on a small machine

This new merging method might need a new branch/pull request

Maybe this will be accepted as is :) I will not hold my breath.

Still, right now it took 21 minutes to merge locally an 80x80x80
mesh

this is way too much. With the new method (-w option) it is taking less,
Still need to quantify

Affected #:  3 files

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 1e0fa73..52cca34 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -42,6 +42,7 @@
 #include "moab/CN.hpp"
 #include "moab/ReadUtilIface.hpp"
 #include "moab/MergeMesh.hpp"
+#include <time.h>
 
 #include <iostream>
 #include <vector>
@@ -52,9 +53,11 @@ using namespace std;
 #define CHECKE(message) if(rval != MB_SUCCESS) { cout << (message) << "\n"; MPI_Finalize(); return 1; }
 int main(int argc, char **argv)
 {
-  int A=2, B=2, C=2, M=2, N=1, K=1;
+  int A=2, B=2, C=2, M=1, N=1, K=1;
   int blockSize = 4;
 
+  bool newMergeMethod=false;
+
   MPI_Init(&argc, &argv);
 
   ProgOptions opts;
@@ -75,6 +78,8 @@ int main(int argc, char **argv)
   opts.addOpt<int>(std::string("zblocks,C"),
       std::string("Number of blocks on a task in x dir (default=2)"), &C);
 
+  opts.addOpt<void>("newMerge,w", "use new merging method",
+          &newMergeMethod);
 
   opts.parseCommandLine(argc, argv);
 
@@ -121,11 +126,13 @@ int main(int argc, char **argv)
   // there are ( M * A * blockSize + 1) *  ( N * B * blockSize + 1 ) * (K * C * blockSize + 1) vertices
   // x is the first dimension that varies
 
+  clock_t tt = clock();
+
   int NX = ( M * A * blockSize + 1);
   int NY = ( N * B * blockSize + 1);
-  int NZ = ( K * C * blockSize + 1);
+  // int NZ = ( K * C * blockSize + 1); not used
   int blockSize1 = blockSize + 1;// used for vertices
-  int bsq = blockSize * blockSize;
+  // int bsq = blockSize * blockSize;
   int b1sq = blockSize1 * blockSize1;
   // generate the block at (a, b, c); it will represent a partition , it will get a partition tag
 
@@ -217,6 +224,14 @@ int main(int argc, char **argv)
       }
     }
   }
+
+  if (0==rank)
+  {
+      std::cout << "generate local mesh:  "
+            << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+      tt = clock();
+  }
+
   // after the mesh is generated on each proc, merge the vertices
   MergeMesh mm(mb);
   Range allhexas;
@@ -224,8 +239,23 @@ int main(int argc, char **argv)
 
   CHECKE("Can't get all hexa  elements.");
 
-  rval = mm.merge_entities( allhexas, 0.0001);
+  Range verts;
+  rval = mb->get_entities_by_dimension(0, 0, verts);
+
+  CHECKE("Can't get all vertices.");
+
+  if (newMergeMethod)
+    rval = mm.merge_using_integer_tag( verts, global_id_tag);
+  else
+    rval = mm.merge_entities(allhexas, 0.0001);
   CHECKE("Can't merge");
+  if (0==rank)
+  {
+     std::cout << "merge locally:  "
+           << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+     tt = clock();
+  }
+
   ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0);
   if (pcomm==NULL)
   {
@@ -234,6 +264,13 @@ int main(int argc, char **argv)
   rval = pcomm->resolve_shared_ents( 0, allhexas, 3, 0 );
   CHECKE("Can't resolve shared ents");
 
+  if (0==rank)
+  {
+     std::cout << "resolve shared entities:  "
+           << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+     tt = clock();
+  }
+
   // delete all quads and edges
   Range toDelete;
   rval =  mb->get_entities_by_dimension(0, 1, toDelete);
@@ -251,9 +288,23 @@ int main(int argc, char **argv)
   rval = pcomm->correct_shared_entities() ;
   CHECKE("Can't correct local shared")
 
+  if (0==rank)
+  {
+    std::cout << "delete edges and faces, and correct sharedEnts:  "
+          << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+    tt = clock();
+  }
+
   rval = mb->write_file("test1.h5m", 0, ";;PARALLEL=WRITE_PART");
   CHECKE("Can't write in parallel");
 
+  if (0==rank)
+  {
+    std::cout << "write file test1.h5m in parallel  "
+          << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+    tt = clock();
+  }
+
   MPI_Finalize();
   return 0;
 }

diff --git a/src/MergeMesh.cpp b/src/MergeMesh.cpp
index cb43b85..62b0022 100644
--- a/src/MergeMesh.cpp
+++ b/src/MergeMesh.cpp
@@ -139,7 +139,79 @@ ErrorCode MergeMesh::perform_merge(Tag merge_tag)
   result = mbImpl->delete_entities(deadEnts);
   return result;
 }
+// merge vertices according to an input tag
+// merge them if the tags are equal
+struct handle_id
+{
+  EntityHandle eh;
+  int val;
+};
+
+// handle structure comparison function for qsort
+// if the id is the same , compare the handle.
+bool compare_handle_id(handle_id ia, handle_id ib) {
+
+  if(ia.val == ib.val) {
+    return ia.eh<ia.eh;
+  } else {
+    return ia.val<ib.val;
+  }
+}
+
+ErrorCode MergeMesh::merge_using_integer_tag(Range & verts, Tag user_tag, Tag merge_tag)
+{
+  ErrorCode rval;
+  DataType tag_type;
+  rval = mbImpl->tag_get_data_type(user_tag, tag_type);
+  if (rval!=MB_SUCCESS || tag_type!=MB_TYPE_INTEGER)
+    return MB_FAILURE;
 
+  std::vector<int> vals(verts.size());
+  rval = mbImpl->tag_get_data(user_tag, verts, &vals[0]);
+  if (rval!=MB_SUCCESS)
+    return rval;
+
+  if (0 == merge_tag)
+  {
+    EntityHandle def_val = 0;
+    rval = mbImpl->tag_get_handle("__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag,
+        MB_TAG_DENSE | MB_TAG_EXCL, &def_val);
+    if (MB_SUCCESS != rval)
+      return rval;
+  }
+  else
+    mbMergeTag = merge_tag;
+
+  std::vector<handle_id>  handles(verts.size());
+  int i=0;
+  for (Range::iterator vit = verts.begin(); vit!= verts.end(); vit++ )
+  {
+    handles[i].eh=*vit;
+    handles[i].val = vals[i];
+    i++;
+  }
+  std::sort(handles.begin(), handles.end(), compare_handle_id);
+
+  i=0;
+  while (i<(int)verts.size()-1)
+  {
+    handle_id  first = handles[i];
+    int j=i+1;
+    while (handles[j].val == first.val && j<(int)verts.size())
+    {
+      rval= mbImpl->tag_set_data(mbMergeTag, &(handles[j].eh), 1, &(first.eh));
+      if (rval!=MB_SUCCESS)
+        return rval;
+      deadEnts.insert(handles[j].eh);
+      j++;
+    }
+    i=j;
+  }
+
+  rval = perform_merge(mbMergeTag);
+
+  return rval;
+}
 ErrorCode MergeMesh::find_merged_to(EntityHandle &tree_root,
     AdaptiveKDTree &tree, Tag merge_tag)
 {

diff --git a/src/moab/MergeMesh.hpp b/src/moab/MergeMesh.hpp
index 8eef9a8..8f4d5bc 100644
--- a/src/moab/MergeMesh.hpp
+++ b/src/moab/MergeMesh.hpp
@@ -32,6 +32,9 @@ public:
   //Identify higher dimension to be merged
   ErrorCode merge_higher_dimensions(Range &elems);
 
+  // merge vertices according to an input tag
+  ErrorCode merge_using_integer_tag(Range & verts, Tag user_tag, Tag merge_tag=0);
+
   //- perform the actual merge
   ErrorCode perform_merge(Tag merged_to);
 private:
@@ -39,8 +42,6 @@ private:
 
   double mergeTol, mergeTolSq;
 
-  Tag mergeTag;
-
   //- given a kdtree, set tag on vertices in leaf nodes with vertices
   //- to which they should be merged
   ErrorCode find_merged_to(EntityHandle &tree_root,


https://bitbucket.org/fathomteam/moab/commits/e162e9ae5eee/
Changeset:   e162e9ae5eee
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:57
Summary:     add quadratic option

with -q option, the hexas generated will be quadratic hex 27 elements
it will allow for more vertices compared to elements (~ 8 times more vertices)

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 52cca34..7de1a8b 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -57,6 +57,7 @@ int main(int argc, char **argv)
   int blockSize = 4;
 
   bool newMergeMethod=false;
+  bool quadratic=false;
 
   MPI_Init(&argc, &argv);
 
@@ -81,6 +82,9 @@ int main(int argc, char **argv)
   opts.addOpt<void>("newMerge,w", "use new merging method",
           &newMergeMethod);
 
+  opts.addOpt<void>("quadratic,q", "use hex 27 elements",
+      &quadratic);
+
   opts.parseCommandLine(argc, argv);
 
   Interface* mb = new Core;
@@ -128,10 +132,16 @@ int main(int argc, char **argv)
 
   clock_t tt = clock();
 
-  int NX = ( M * A * blockSize + 1);
-  int NY = ( N * B * blockSize + 1);
+  int q = 1;
+  if (quadratic)
+  {
+    q = 2;
+  }
+  int NX = (q * M * A * blockSize + 1);
+  int NY = (q * N * B * blockSize + 1);
   // int NZ = ( K * C * blockSize + 1); not used
-  int blockSize1 = blockSize + 1;// used for vertices
+  int blockSize1 = q*blockSize + 1;// used for vertices
+
   // int bsq = blockSize * blockSize;
   int b1sq = blockSize1 * blockSize1;
   // generate the block at (a, b, c); it will represent a partition , it will get a partition tag
@@ -151,27 +161,28 @@ int main(int argc, char **argv)
     {
       for (int c=0; c<C; c++)
       {
-        // we will generate (block+1)^3 vertices, and block^3 hexas
+        // we will generate (q*block+1)^3 vertices, and block^3 hexas; q is 1 for linear, 2 for quadratic
         // the global id of the vertices will come from m, n, k, a, b, c
-        // x will vary from  m*A*block + a*block to m*A*block+(a+1)*block etc;
-        int num_nodes = (blockSize+1)*(blockSize+1)*(blockSize+1);
+        // x will vary from  m*A*q*block + a*q*block to m*A*q*block+(a+1)*q*block etc;
+        int num_nodes = blockSize1*blockSize1*blockSize1;
 
         vector<double*> arrays;
         EntityHandle startv;
         rval = iface->get_node_coords(3, num_nodes, 0, startv, arrays);
         CHECKE("Can't get node coords.");
 
-        int x = m*A*blockSize + a*blockSize;
-        int y = n*B*blockSize + b*blockSize;
-        int z = k*C*blockSize + c*blockSize;
+        // will start with the lower corner:
+        int x = m*A*q*blockSize + a*q*blockSize;
+        int y = n*B*q*blockSize + b*q*blockSize;
+        int z = k*C*q*blockSize + c*q*blockSize;
         int ix=0;
         vector<int>  gids(num_nodes);
         Range verts(startv, startv+num_nodes-1);
-        for (int kk=0; kk<blockSize+1; kk++)
+        for (int kk=0; kk<blockSize1; kk++)
         {
-          for (int jj=0; jj<blockSize+1; jj++)
+          for (int jj=0; jj<blockSize1; jj++)
           {
-            for (int ii=0; ii<blockSize+1; ii++)
+            for (int ii=0; ii<blockSize1; ii++)
             {
               arrays[0][ix] = (double)x+ii;
               arrays[1][ix] = (double)y+jj;
@@ -185,7 +196,12 @@ int main(int argc, char **argv)
         int num_hexas = (blockSize)*(blockSize)*(blockSize);
         EntityHandle starte; // connectivity
         EntityHandle * conn;
-        rval = iface->get_element_connect(num_hexas, 8, MBHEX, 0, starte, conn);
+        if (quadratic)
+        {
+          rval = iface->get_element_connect(num_hexas, 27, MBHEX, 0, starte, conn);
+        }
+        else
+          rval = iface->get_element_connect(num_hexas, 8, MBHEX, 0, starte, conn);
         CHECKE("Can't get hexa  connectivity.");
 
         Range hexas(starte, starte+num_hexas-1);
@@ -197,16 +213,50 @@ int main(int argc, char **argv)
           {
             for (int ii=0; ii<blockSize; ii++)
             {
-              EntityHandle corner=startv + ii + (jj*blockSize1) + kk * b1sq;
-              conn[ix]=corner;
-              conn[ix+1]=corner+1;
-              conn[ix+2]=corner+ 1 + blockSize1;
-              conn[ix+3]=corner+ blockSize1;
-              conn[ix+4]=corner + b1sq;
-              conn[ix+5]=corner+1 + b1sq;
-              conn[ix+6]=corner+ 1 + blockSize1 + b1sq;
-              conn[ix+7]=corner+ blockSize1 + b1sq;
-              ix+=8;
+              EntityHandle corner=startv + q*ii + (q*jj*blockSize1) + q*kk * b1sq;
+              if (quadratic)
+              {
+                conn[ix]=corner;
+                conn[ix+1]=corner+2;
+                conn[ix+2]=corner+2 + 2 * blockSize1;
+                conn[ix+3]=corner+   2 * blockSize1;
+                conn[ix+4]=corner +                   2 * b1sq;
+                conn[ix+5]=corner+ 2 +                2 * b1sq;
+                conn[ix+6]=corner+ 2 + 2 * blockSize1 + 2 * b1sq;
+                conn[ix+7]=corner+       2*blockSize1 +  2*b1sq;
+                conn[ix+8]=corner + 1  ;   // 0-1
+                conn[ix+9]=  corner+ 2 + blockSize1; //1-2
+                conn[ix+10]= corner + 1 + 2 * blockSize1;// 2-3
+                conn[ix+11] =corner + blockSize1; // 3-0
+                conn[ix+12]= corner + b1sq; // 0-4
+                conn[ix+13]= corner+2 + b1sq; // 1-5
+                conn[ix+14]=  corner + 2+ 2*blockSize1 +   b1sq; // 2-6
+                conn[ix+15] = corner +    2* blockSize1 +   b1sq; // 3-7
+                conn[ix+16] = corner + 1 +                 2*b1sq; // 4-5
+                conn[ix+17] = corner + 2 + blockSize1 +   2* b1sq; //5-6
+                conn[ix+18] = corner + 1 + 2*blockSize1 + 2* b1sq; // 6-7
+                conn[ix+19] = corner +     blockSize1 +   2* b1sq; // 4-7
+                conn[ix+20] = corner + 1 +                 b1sq; // 0154
+                conn[ix+21] = corner + 2 + blockSize1 +    b1sq; // 1265
+                conn[ix+22] = corner + 1 + 2*blockSize1 +  b1sq ; // 2376
+                conn[ix+23] = corner +     blockSize1 +    b1sq; // 0374
+                conn[ix+24] = corner + 1 + blockSize1          ; // 0123
+                conn[ix+25] = corner + 1 + blockSize1 +   2*b1sq; // 4567
+                conn[ix+26] = corner + 1 + blockSize1 +    b1sq; // center
+                ix+=27;
+              }
+              else
+              {
+                conn[ix]=corner;
+                conn[ix+1]=corner+1;
+                conn[ix+2]=corner+ 1 + blockSize1;
+                conn[ix+3]=corner+ blockSize1;
+                conn[ix+4]=corner + b1sq;
+                conn[ix+5]=corner+1 + b1sq;
+                conn[ix+6]=corner+ 1 + blockSize1 + b1sq;
+                conn[ix+7]=corner+ blockSize1 + b1sq;
+                ix+=8;
+              }
             }
           }
         }
@@ -232,6 +282,9 @@ int main(int argc, char **argv)
       tt = clock();
   }
 
+  // before merge locally
+  rval = mb->write_file("test0.h5m", 0, ";;PARALLEL=WRITE_PART");
+  CHECKE("Can't write in parallel, before merging");
   // after the mesh is generated on each proc, merge the vertices
   MergeMesh mm(mb);
   Range allhexas;


https://bitbucket.org/fathomteam/moab/commits/461454043992/
Changeset:   461454043992
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     add global id tag on elements too

also, add ystride and zstride variables, for clarity, and
a drawing with hex27 connectivity

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 7de1a8b..16cd9a0 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -11,11 +11,11 @@
  *  The number of tasks will be MxNxK, and it must match the mpi size
  *  Each task will generate its mesh at location (m,n,k)
  *
- *  By default M=2, N=1, K=1, so by default it should be launched on 2 procs
+ *  By default M=1, N=1, K=1, so by default it should be launched on 1 proc
  *  By default, blockSize is 4, and A=2, B=2, C=2, so each task will generate locally
  *    blockSize^3 x A x B x C hexahedrons (value = 64x8 = 512 hexas, in 8 partitions)
  *
- *  The total number f partitions will be A*B*C*M*N*K (default 16)
+ *  The total number of partitions will be A*B*C*M*N*K (default 8)
  *
  *  Each part in partition will get a proper tag
  *
@@ -33,7 +33,14 @@
  *  entity handle within a partition
  *
  *
- *  mpiexec -np 2 GenLargeMesh
+ *  to run: ./GenLargeMesh
+ *
+ *  when launched on more procs, you have to make sure
+ *  num procs = M*N*K
+ *
+ *  so you can launch with
+ *  mpiexec -np 8 ./GenLargeMesh -M 2 -N 2 -K 2
+ *
  */
 
 #include "moab/Core.hpp"
@@ -139,11 +146,15 @@ int main(int argc, char **argv)
   }
   int NX = (q * M * A * blockSize + 1);
   int NY = (q * N * B * blockSize + 1);
+  int nex = M * A * blockSize; // number of elements in x direction, used for global id on element
+  int ney = N * B * blockSize; // number of elements in y direction  ....
   // int NZ = ( K * C * blockSize + 1); not used
   int blockSize1 = q*blockSize + 1;// used for vertices
 
-  // int bsq = blockSize * blockSize;
-  int b1sq = blockSize1 * blockSize1;
+  //int xstride = 1;
+  int ystride = blockSize1;
+
+  int zstride = blockSize1 * blockSize1;
   // generate the block at (a, b, c); it will represent a partition , it will get a partition tag
 
   Tag global_id_tag;
@@ -187,7 +198,7 @@ int main(int argc, char **argv)
               arrays[0][ix] = (double)x+ii;
               arrays[1][ix] = (double)y+jj;
               arrays[2][ix] = (double)z+kk;
-              gids[ix] = 1 + (z+kk) * (NX*NY) + (y+jj) * NX + (x+ii);
+              gids[ix] = 1 + (x+ii) + (y+jj) * NX + (z+kk) * (NX*NY) ;
               ix++;
             }
           }
@@ -207,54 +218,79 @@ int main(int argc, char **argv)
         Range hexas(starte, starte+num_hexas-1);
         // fill  hexas
         ix=0;
+        // identify the elements at the lower corner, for their global ids
+        int xe = m*A*blockSize + a*blockSize;
+        int ye = n*B*blockSize + b*blockSize;
+        int ze = k*C*blockSize + c*blockSize;
+        gids.resize(num_hexas);
+        int ie=0; // index now in the elements, for global ids
         for (int kk=0; kk<blockSize; kk++)
         {
           for (int jj=0; jj<blockSize; jj++)
           {
             for (int ii=0; ii<blockSize; ii++)
             {
-              EntityHandle corner=startv + q*ii + (q*jj*blockSize1) + q*kk * b1sq;
+              EntityHandle corner=startv + q * ii + q * jj * ystride + q * kk * zstride;
+              gids[ie] = 1 + (xe+ii) + (ye+jj) * nex + (ze+kk) * (nex*ney) ;
+              ie++;
               if (quadratic)
               {
-                conn[ix]=corner;
-                conn[ix+1]=corner+2;
-                conn[ix+2]=corner+2 + 2 * blockSize1;
-                conn[ix+3]=corner+   2 * blockSize1;
-                conn[ix+4]=corner +                   2 * b1sq;
-                conn[ix+5]=corner+ 2 +                2 * b1sq;
-                conn[ix+6]=corner+ 2 + 2 * blockSize1 + 2 * b1sq;
-                conn[ix+7]=corner+       2*blockSize1 +  2*b1sq;
-                conn[ix+8]=corner + 1  ;   // 0-1
-                conn[ix+9]=  corner+ 2 + blockSize1; //1-2
-                conn[ix+10]= corner + 1 + 2 * blockSize1;// 2-3
-                conn[ix+11] =corner + blockSize1; // 3-0
-                conn[ix+12]= corner + b1sq; // 0-4
-                conn[ix+13]= corner+2 + b1sq; // 1-5
-                conn[ix+14]=  corner + 2+ 2*blockSize1 +   b1sq; // 2-6
-                conn[ix+15] = corner +    2* blockSize1 +   b1sq; // 3-7
-                conn[ix+16] = corner + 1 +                 2*b1sq; // 4-5
-                conn[ix+17] = corner + 2 + blockSize1 +   2* b1sq; //5-6
-                conn[ix+18] = corner + 1 + 2*blockSize1 + 2* b1sq; // 6-7
-                conn[ix+19] = corner +     blockSize1 +   2* b1sq; // 4-7
-                conn[ix+20] = corner + 1 +                 b1sq; // 0154
-                conn[ix+21] = corner + 2 + blockSize1 +    b1sq; // 1265
-                conn[ix+22] = corner + 1 + 2*blockSize1 +  b1sq ; // 2376
-                conn[ix+23] = corner +     blockSize1 +    b1sq; // 0374
-                conn[ix+24] = corner + 1 + blockSize1          ; // 0123
-                conn[ix+25] = corner + 1 + blockSize1 +   2*b1sq; // 4567
-                conn[ix+26] = corner + 1 + blockSize1 +    b1sq; // center
+
+  //                    4   ----- 19   -----  7
+  //                .   |                 .   |
+  //            16         25         18      |
+  //         .          |          .          |
+  //      5   ----- 17   -----  6             |
+  //      |            12       | 23         15
+  //      |                     |             |
+  //      |     20      |  26   |     22      |
+  //      |                     |             |
+  //     13         21  |      14             |
+  //      |             0   ----- 11   -----  3
+  //      |         .           |         .
+  //      |      8         24   |     10
+  //      |  .                  |  .
+  //      1   -----  9   -----  2
+  //
+                conn[ix]=    corner;
+                conn[ix+1]=  corner + 2;
+                conn[ix+2]=  corner + 2 + 2 * ystride;
+                conn[ix+3]=  corner +     2 * ystride;
+                conn[ix+4]=  corner                   + 2 * zstride;
+                conn[ix+5]=  corner + 2               + 2 * zstride;
+                conn[ix+6]=  corner + 2 + 2 * ystride + 2 * zstride;
+                conn[ix+7]=  corner +     2 * ystride + 2 * zstride;
+                conn[ix+8]=  corner + 1  ;                                         // 0-1
+                conn[ix+9]=  corner + 2 +     ystride ;                            // 1-2
+                conn[ix+10]= corner + 1 + 2 * ystride ;                            // 2-3
+                conn[ix+11]= corner +         ystride;                             // 3-0
+                conn[ix+12]= corner +                       zstride;               // 0-4
+                conn[ix+13]= corner + 2 +                   zstride;               // 1-5
+                conn[ix+14]= corner + 2 + 2 * ystride +     zstride;               // 2-6
+                conn[ix+15]= corner +     2 * ystride +     zstride;               // 3-7
+                conn[ix+16]= corner + 1 +               2 * zstride;               // 4-5
+                conn[ix+17]= corner + 2 +     ystride + 2 * zstride;               //5-6
+                conn[ix+18]= corner + 1 + 2 * ystride + 2 * zstride;               // 6-7
+                conn[ix+19]= corner +         ystride + 2 * zstride;               // 4-7
+                conn[ix+20]= corner + 1 +                   zstride;               // 0154
+                conn[ix+21]= corner + 2 +     ystride +     zstride;               // 1265
+                conn[ix+22]= corner + 1 + 2 * ystride +     zstride;               // 2376
+                conn[ix+23]= corner +         ystride +     zstride;               // 0374
+                conn[ix+24]= corner + 1 +     ystride              ;               // 0123
+                conn[ix+25]= corner + 1 +     ystride + 2 * zstride;               // 4567
+                conn[ix+26]= corner + 1 +     ystride +     zstride;               // center
                 ix+=27;
               }
               else
               {
-                conn[ix]=corner;
-                conn[ix+1]=corner+1;
-                conn[ix+2]=corner+ 1 + blockSize1;
-                conn[ix+3]=corner+ blockSize1;
-                conn[ix+4]=corner + b1sq;
-                conn[ix+5]=corner+1 + b1sq;
-                conn[ix+6]=corner+ 1 + blockSize1 + b1sq;
-                conn[ix+7]=corner+ blockSize1 + b1sq;
+                conn[ix]=  corner;
+                conn[ix+1]=corner + 1;
+                conn[ix+2]=corner + 1 + ystride;
+                conn[ix+3]=corner +     ystride;
+                conn[ix+4]=corner +               zstride;
+                conn[ix+5]=corner + 1 +           zstride;
+                conn[ix+6]=corner + 1 + ystride + zstride;
+                conn[ix+7]=corner +     ystride + zstride;
                 ix+=8;
               }
             }
@@ -265,6 +301,8 @@ int main(int argc, char **argv)
         CHECKE("Can't create mesh set.");
         rval = mb->add_entities(part_set, hexas);
         CHECKE("Can't add entities to set.");
+        rval = mb->tag_set_data(global_id_tag, hexas, &gids[0]);
+        CHECKE("Can't set global ids to elements.");
         int part_num= a +  m*A + (b + n*B)*(M*A) + (c+k*C)*(M*A * N*B);
         rval = mb->tag_set_data(part_tag, &part_set, 1, &part_num);
         CHECKE("Can't set part tag on set");


https://bitbucket.org/fathomteam/moab/commits/649779ff57a1/
Changeset:   649779ff57a1
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     add tags to the large mesh example

use -i option for integer tags name on vertices,
 and -d for double tags on cells
we can have multiple integer and double tags, with something like
GenLargeMesh -i tag1 -i tag2 -d tag3 -d tag4

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 16cd9a0..08dad08 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -92,8 +92,19 @@ int main(int argc, char **argv)
   opts.addOpt<void>("quadratic,q", "use hex 27 elements",
       &quadratic);
 
+  vector<string> intTagNames;
+  string firstIntTag;
+  opts.addOpt<string>(std::string("int_tag_vert,i"), string("add integer tag on vertices"), &firstIntTag);
+
+  vector<string> doubleTagNames;
+  string firstDoubleTag;
+  opts.addOpt<string>(std::string("double_tag_cell,d"), string("add double tag on cells"), &firstDoubleTag);
+
   opts.parseCommandLine(argc, argv);
 
+  opts.getOptAllArgs("int_tag_vert,i", intTagNames);
+  opts.getOptAllArgs("double_tag_cell,d", doubleTagNames);
+
   Interface* mb = new Core;
   if (NULL == mb)
   {
@@ -166,6 +177,22 @@ int main(int argc, char **argv)
   mb->tag_get_handle("PARALLEL_PARTITION", 1, MB_TYPE_INTEGER,
                  part_tag, MB_TAG_CREAT|MB_TAG_SPARSE, &dum_id);
 
+  // create tags on vertices and cells, look in the list of options
+  vector<Tag> intTags(intTagNames.size());
+  vector<Tag> doubleTags(doubleTagNames.size());
+  for (size_t i=0; i<intTagNames.size(); i++)
+  {
+    rval = mb->tag_get_handle(intTagNames[i].c_str(), 1, MB_TYPE_INTEGER, intTags[i],
+        MB_TAG_CREAT|MB_TAG_DENSE, &dum_id);
+    CHECKE("Can't create integer tag.");
+  }
+  double defval=0.;
+  for (size_t i=0; i<doubleTagNames.size(); i++)
+  {
+    rval = mb->tag_get_handle(doubleTagNames[i].c_str(), 1, MB_TYPE_DOUBLE, doubleTags[i],
+        MB_TAG_CREAT|MB_TAG_DENSE, &defval);
+    CHECKE("Can't create double tag.");
+  }
   for (int a=0; a<A; a++)
   {
     for (int b=0; b<B; b++)
@@ -199,6 +226,13 @@ int main(int argc, char **argv)
               arrays[1][ix] = (double)y+jj;
               arrays[2][ix] = (double)z+kk;
               gids[ix] = 1 + (x+ii) + (y+jj) * NX + (z+kk) * (NX*NY) ;
+              // set int tags, some nice values?
+              EntityHandle v = startv + ix;
+              for (size_t i=0; i<intTags.size(); i++)
+              {
+                int valv=gids[ix]/2+3 + i*1000;
+                mb->tag_set_data(intTags[i], &v, 1, &valv);
+              }
               ix++;
             }
           }
@@ -232,7 +266,15 @@ int main(int argc, char **argv)
             {
               EntityHandle corner=startv + q * ii + q * jj * ystride + q * kk * zstride;
               gids[ie] = 1 + (xe+ii) + (ye+jj) * nex + (ze+kk) * (nex*ney) ;
+
+              EntityHandle eh = starte + ie;
               ie++;
+              for (size_t i=0; i<doubleTags.size(); i++)
+              {
+                double valv=gids[ie]/30. + i*5000.;
+                mb->tag_set_data(doubleTags[i], &eh, 1, &valv);
+              }
+
               if (quadratic)
               {
 
@@ -320,9 +362,9 @@ int main(int argc, char **argv)
       tt = clock();
   }
 
-  // before merge locally
+  /*// before merge locally
   rval = mb->write_file("test0.h5m", 0, ";;PARALLEL=WRITE_PART");
-  CHECKE("Can't write in parallel, before merging");
+  CHECKE("Can't write in parallel, before merging");*/
   // after the mesh is generated on each proc, merge the vertices
   MergeMesh mm(mb);
   Range allhexas;


https://bitbucket.org/fathomteam/moab/commits/1e6784293d7f/
Changeset:   1e6784293d7f
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     introduce ParallelComm::delete_entities(Range &)

when a user deletes an entity, in a parallel environment,
it needs to make sure it is deleted from all tasks that might
share it
Also, update the sharedEnts local vector on each

Affected #:  3 files

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 08dad08..8e5e029 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -408,18 +408,12 @@ int main(int argc, char **argv)
   Range toDelete;
   rval =  mb->get_entities_by_dimension(0, 1, toDelete);
   CHECKE("Can't get edges");
-  rval = mb->delete_entities(toDelete);
-  CHECKE("Can't delete edges");
-
-  toDelete.clear();
 
   rval = mb->get_entities_by_dimension(0, 2, toDelete);
   CHECKE("Can't get faces");
-  rval = mb->delete_entities(toDelete);
-  CHECKE("Can't delete faces");
 
-  rval = pcomm->correct_shared_entities() ;
-  CHECKE("Can't correct local shared")
+  rval = pcomm->delete_entities(toDelete) ;
+  CHECKE("Can't delete entities")
 
   if (0==rank)
   {

diff --git a/src/parallel/ParallelComm.cpp b/src/parallel/ParallelComm.cpp
index 4fb8afc..49c1606 100644
--- a/src/parallel/ParallelComm.cpp
+++ b/src/parallel/ParallelComm.cpp
@@ -8792,19 +8792,75 @@ ErrorCode ParallelComm::settle_intersection_points(Range & edges, Range & shared
   return MB_SUCCESS;
   // end copy
 }
-ErrorCode ParallelComm::correct_shared_entities()
+ErrorCode ParallelComm::delete_entities(Range & to_delete)
 {
 
+  // will not look at shared sets yet, but maybe we should
+  // first, see if any of the entities to delete is shared; then inform the other processors
+  // about their fate (to be deleted), using a crystal router transfer
+  ErrorCode rval=MB_SUCCESS;
+  unsigned char pstat;
+  EntityHandle tmp_handles[MAX_SHARING_PROCS];
+  int tmp_procs[MAX_SHARING_PROCS];
+  unsigned int num_ps;
+  TupleList ents_to_delete;
+  ents_to_delete.initialize(1, 0, 1, 0, to_delete.size() * (MAX_SHARING_PROCS+1) );// a little bit of overkill
+  ents_to_delete.enableWriteAccess();
+  unsigned int i = 0;
+  for (Range::iterator it=to_delete.begin(); it!=to_delete.end(); it++)
+  {
+    EntityHandle eh=*it; // entity to be deleted
+
+    rval = get_sharing_data(eh, tmp_procs, tmp_handles,
+                                   pstat, num_ps);
+    if (rval!=MB_SUCCESS || num_ps==0)
+      continue;
+    // add to the tuple list the information to be sent (to the remote procs)
+
+    for (unsigned int p = 0; p < num_ps; p++)
+    {
+      ents_to_delete.vi_wr[i] = tmp_procs[p];
+      ents_to_delete.vul_wr[i] = (unsigned long)tmp_handles[p];
+      i++;
+      ents_to_delete.inc_n();
+    }
+  }
+
+  gs_data::crystal_data *cd = this->procConfig.crystal_router();
+  // all communication happens here; no other mpi calls
+  // also, this is a collective call
+  rval = cd->gs_transfer(1,ents_to_delete,0);
+
+  if (MB_SUCCESS!= rval)
+  {
+    std::cout << "error in tuple transfer\n";
+    return rval;
+  }
+  // add to the range of ents to delete the new ones that were sent from other procs
+  unsigned int received = ents_to_delete.get_n();
+  for (i=0; i< received; i++)
+  {
+    //int from = ents_to_delete.vi_rd[i];
+    unsigned long valrec = ents_to_delete.vul_rd[i];
+    to_delete.insert((EntityHandle)valrec);
+  }
+  rval = mbImpl->delete_entities(to_delete);
+  if (MB_SUCCESS!= rval)
+  {
+    std::cout << "error in deleting actual entities\n";
+    return rval;
+  }
   std::vector<EntityHandle> good_ents;
-  for (size_t i=0; i<sharedEnts.size(); i++)
+  for (size_t j=0; j<sharedEnts.size(); j++)
   {
-    if (mbImpl->is_valid(sharedEnts[i]))
+    if (mbImpl->is_valid(sharedEnts[j]))
     {
-      good_ents.push_back(sharedEnts[i]);
+      good_ents.push_back(sharedEnts[j]);
     }
   }
   sharedEnts = good_ents;
 
+  // what about shared sets? who is updating them?
   return MB_SUCCESS;
 }
 

diff --git a/src/parallel/moab/ParallelComm.hpp b/src/parallel/moab/ParallelComm.hpp
index 2a16ac2..976236e 100644
--- a/src/parallel/moab/ParallelComm.hpp
+++ b/src/parallel/moab/ParallelComm.hpp
@@ -932,11 +932,11 @@ namespace moab {
     ErrorCode settle_intersection_points(Range & edges, Range & shared_edges_owned,
         std::vector<std::vector<EntityHandle> *> & extraNodesVec, double tolerance);
 
-    /* \brief check the shared entities if they exist anymore
+    /* \brief delete entities from moab database
      * will check the shared ents array, and clean it if necessary
      *
      */
-    ErrorCode correct_shared_entities();
+    ErrorCode delete_entities(Range & to_delete);
 
   private:
 


https://bitbucket.org/fathomteam/moab/commits/2484581d45f1/
Changeset:   2484581d45f1
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     add the option to keep the skins

these skins result from resolving shared entities
we may want to keep them

By default, we delete them

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 8e5e029..ce67ac9 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -65,6 +65,7 @@ int main(int argc, char **argv)
 
   bool newMergeMethod=false;
   bool quadratic=false;
+  bool keep_skins=false;
 
   MPI_Init(&argc, &argv);
 
@@ -92,6 +93,9 @@ int main(int argc, char **argv)
   opts.addOpt<void>("quadratic,q", "use hex 27 elements",
       &quadratic);
 
+  opts.addOpt<void>("keep_skins,k", "keep skins with shared entities",
+        &keep_skins);
+
   vector<string> intTagNames;
   string firstIntTag;
   opts.addOpt<string>(std::string("int_tag_vert,i"), string("add integer tag on vertices"), &firstIntTag);
@@ -404,22 +408,25 @@ int main(int argc, char **argv)
      tt = clock();
   }
 
-  // delete all quads and edges
-  Range toDelete;
-  rval =  mb->get_entities_by_dimension(0, 1, toDelete);
-  CHECKE("Can't get edges");
+  if (!keep_skins) // default is to delete the 1- and 2-dimensional entities
+  {
+    // delete all quads and edges
+    Range toDelete;
+    rval =  mb->get_entities_by_dimension(0, 1, toDelete);
+    CHECKE("Can't get edges");
 
-  rval = mb->get_entities_by_dimension(0, 2, toDelete);
-  CHECKE("Can't get faces");
+    rval = mb->get_entities_by_dimension(0, 2, toDelete);
+    CHECKE("Can't get faces");
 
-  rval = pcomm->delete_entities(toDelete) ;
-  CHECKE("Can't delete entities")
+    rval = pcomm->delete_entities(toDelete) ;
+    CHECKE("Can't delete entities")
 
-  if (0==rank)
-  {
-    std::cout << "delete edges and faces, and correct sharedEnts:  "
-          << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
-    tt = clock();
+    if (0==rank)
+    {
+      std::cout << "delete edges and faces, and correct sharedEnts:  "
+            << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+      tt = clock();
+    }
   }
 
   rval = mb->write_file("test1.h5m", 0, ";;PARALLEL=WRITE_PART");


https://bitbucket.org/fathomteam/moab/commits/b11de15ea2b2/
Changeset:   b11de15ea2b2
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     add tetra option

also, add option for output file, default was test1.h5m
tetra does not work with quadratic yet, but it will someday

I mean, the tetras are all linear right now.
if -q and -t are used together, quadratic hexas will be generated

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index ce67ac9..5918a09 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -7,6 +7,7 @@
  *  Each processor will create its version of a block mesh, partitioned
  *  as AxBxC blocks. Each block will be with blockSize^3 hexahedrons, and will get a
  *  different PARALLEL_PARTITION tag
+ *  When -t option is used, instead of each hex, we are creating 6 tetrahedrons
  *
  *  The number of tasks will be MxNxK, and it must match the mpi size
  *  Each task will generate its mesh at location (m,n,k)
@@ -14,7 +15,7 @@
  *  By default M=1, N=1, K=1, so by default it should be launched on 1 proc
  *  By default, blockSize is 4, and A=2, B=2, C=2, so each task will generate locally
  *    blockSize^3 x A x B x C hexahedrons (value = 64x8 = 512 hexas, in 8 partitions)
- *
+ *   (if -t, multiple by 6 for total number of cells/tets)
  *  The total number of partitions will be A*B*C*M*N*K (default 8)
  *
  *  Each part in partition will get a proper tag
@@ -41,6 +42,8 @@
  *  so you can launch with
  *  mpiexec -np 8 ./GenLargeMesh -M 2 -N 2 -K 2
  *
+ *  We also added -q option; it works now only for hexa mesh, it will generate
+ *  quadratic hex27 elements
  */
 
 #include "moab/Core.hpp"
@@ -66,6 +69,7 @@ int main(int argc, char **argv)
   bool newMergeMethod=false;
   bool quadratic=false;
   bool keep_skins=false;
+  bool tetra = false;
 
   MPI_Init(&argc, &argv);
 
@@ -96,6 +100,9 @@ int main(int argc, char **argv)
   opts.addOpt<void>("keep_skins,k", "keep skins with shared entities",
         &keep_skins);
 
+  opts.addOpt<void>("tetrahedrons,t", "generate tetrahedrons ",
+        &tetra);
+
   vector<string> intTagNames;
   string firstIntTag;
   opts.addOpt<string>(std::string("int_tag_vert,i"), string("add integer tag on vertices"), &firstIntTag);
@@ -104,6 +111,9 @@ int main(int argc, char **argv)
   string firstDoubleTag;
   opts.addOpt<string>(std::string("double_tag_cell,d"), string("add double tag on cells"), &firstDoubleTag);
 
+  string outFileName="test1.h5m";
+  opts.addOpt<std::string> ("outFile,o", "Specify the output file name string (default test1.h5m)", &outFileName );
+
   opts.parseCommandLine(argc, argv);
 
   opts.getOptAllArgs("int_tag_vert,i", intTagNames);
@@ -159,6 +169,9 @@ int main(int argc, char **argv)
   {
     q = 2;
   }
+  int factor =1;
+  if (tetra)
+    factor =6;
   int NX = (q * M * A * blockSize + 1);
   int NY = (q * N * B * blockSize + 1);
   int nex = M * A * blockSize; // number of elements in x direction, used for global id on element
@@ -243,24 +256,26 @@ int main(int argc, char **argv)
         }
         mb->tag_set_data(global_id_tag, verts, &gids[0]);
         int num_hexas = (blockSize)*(blockSize)*(blockSize);
+        int num_el = num_hexas * factor;
+
         EntityHandle starte; // connectivity
         EntityHandle * conn;
         if (quadratic)
-        {
-          rval = iface->get_element_connect(num_hexas, 27, MBHEX, 0, starte, conn);
-        }
+          rval = iface->get_element_connect(num_el, 27, MBHEX, 0, starte, conn);
+        else if (tetra)
+          rval = iface->get_element_connect(num_el, 4, MBTET, 0, starte, conn);
         else
-          rval = iface->get_element_connect(num_hexas, 8, MBHEX, 0, starte, conn);
-        CHECKE("Can't get hexa  connectivity.");
+          rval = iface->get_element_connect(num_el, 8, MBHEX, 0, starte, conn);
+        CHECKE("Can't get element  connectivity.");
 
-        Range hexas(starte, starte+num_hexas-1);
+        Range hexas(starte, starte+num_el-1); // should be elements
         // fill  hexas
         ix=0;
         // identify the elements at the lower corner, for their global ids
         int xe = m*A*blockSize + a*blockSize;
         int ye = n*B*blockSize + b*blockSize;
         int ze = k*C*blockSize + c*blockSize;
-        gids.resize(num_hexas);
+        gids.resize(num_el);
         int ie=0; // index now in the elements, for global ids
         for (int kk=0; kk<blockSize; kk++)
         {
@@ -269,16 +284,15 @@ int main(int argc, char **argv)
             for (int ii=0; ii<blockSize; ii++)
             {
               EntityHandle corner=startv + q * ii + q * jj * ystride + q * kk * zstride;
-              gids[ie] = 1 + (xe+ii) + (ye+jj) * nex + (ze+kk) * (nex*ney) ;
+              gids[ie] = 1 + ((xe+ii) + (ye+jj) * nex + (ze+kk) * (nex*ney))*factor ; // 6 more for tetra
 
               EntityHandle eh = starte + ie;
-              ie++;
               for (size_t i=0; i<doubleTags.size(); i++)
               {
                 double valv=gids[ie]/30. + i*5000.;
                 mb->tag_set_data(doubleTags[i], &eh, 1, &valv);
               }
-
+              ie++;
               if (quadratic)
               {
 
@@ -327,7 +341,73 @@ int main(int argc, char **argv)
                 conn[ix+26]= corner + 1 +     ystride +     zstride;               // center
                 ix+=27;
               }
-              else
+              else if (tetra)
+              {
+
+                //        E      H
+                //     F     G
+                //
+                //        A     D
+                //     B     C
+                EntityHandle AA = corner;
+                EntityHandle BB = corner + 1;
+                EntityHandle CC = corner + 1 + ystride;
+                EntityHandle D = corner +     ystride;
+                EntityHandle E = corner +               zstride;
+                EntityHandle F = corner + 1 +           zstride;
+                EntityHandle G = corner + 1 + ystride + zstride;
+                EntityHandle H = corner +     ystride + zstride;
+
+                // tet EDHG
+                conn[ix]    = E;
+                conn[ix+1]  = D;
+                conn[ix+2]  = H;
+                conn[ix+3]  = G;
+
+                // tet ABCF
+                conn[ix+4]  = AA;
+                conn[ix+5]  = BB;
+                conn[ix+6]  = CC;
+                conn[ix+7]  = F;
+
+                // tet ADEF
+                conn[ix+8]  = AA;
+                conn[ix+9]  = D;
+                conn[ix+10] = E;
+                conn[ix+11] = F;
+
+                // tet CGDF
+                conn[ix+12] = CC;
+                conn[ix+13] = G;
+                conn[ix+14] = D;
+                conn[ix+15] = F;
+
+                // tet ACDF
+                conn[ix+16] = AA;
+                conn[ix+17] = CC;
+                conn[ix+18] = D;
+                conn[ix+19] = F;
+
+                // tet DGEF
+                conn[ix+20] = D;
+                conn[ix+21] = G;
+                conn[ix+22] = E;
+                conn[ix+23] = F;
+                ix+=24;
+                for (int ff=0; ff<factor-1; ff++)
+                {
+                  gids[ie] = gids[ie-1]+1 ; // 6 more for tetra
+
+                  eh = starte + ie;
+                  for (size_t i=0; i<doubleTags.size(); i++)
+                  {
+                    double valv=gids[ie]/30. + i*5000.;
+                    mb->tag_set_data(doubleTags[i], &eh, 1, &valv);
+                  }
+                  ie++;
+                }
+              }
+              else // linear hex
               {
                 conn[ix]=  corner;
                 conn[ix+1]=corner + 1;
@@ -371,10 +451,10 @@ int main(int argc, char **argv)
   CHECKE("Can't write in parallel, before merging");*/
   // after the mesh is generated on each proc, merge the vertices
   MergeMesh mm(mb);
-  Range allhexas;
-  rval = mb->get_entities_by_type(0, MBHEX, allhexas);
+  Range all3dcells;
+  rval = mb->get_entities_by_dimension(0, 3, all3dcells);
 
-  CHECKE("Can't get all hexa  elements.");
+  CHECKE("Can't get all 3d cells  elements.");
 
   Range verts;
   rval = mb->get_entities_by_dimension(0, 0, verts);
@@ -384,7 +464,7 @@ int main(int argc, char **argv)
   if (newMergeMethod)
     rval = mm.merge_using_integer_tag( verts, global_id_tag);
   else
-    rval = mm.merge_entities(allhexas, 0.0001);
+    rval = mm.merge_entities(all3dcells, 0.0001);
   CHECKE("Can't merge");
   if (0==rank)
   {
@@ -398,7 +478,7 @@ int main(int argc, char **argv)
   {
     pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
   }
-  rval = pcomm->resolve_shared_ents( 0, allhexas, 3, 0 );
+  rval = pcomm->resolve_shared_ents( 0, all3dcells, 3, 0 );
   CHECKE("Can't resolve shared ents");
 
   if (0==rank)
@@ -429,12 +509,12 @@ int main(int argc, char **argv)
     }
   }
 
-  rval = mb->write_file("test1.h5m", 0, ";;PARALLEL=WRITE_PART");
+  rval = mb->write_file(outFileName.c_str(), 0, ";;PARALLEL=WRITE_PART");
   CHECKE("Can't write in parallel");
 
   if (0==rank)
   {
-    std::cout << "write file test1.h5m in parallel  "
+    std::cout << "write file " << outFileName << " in "
           << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
     tt = clock();
   }


https://bitbucket.org/fathomteam/moab/commits/2ab0e45260af/
Changeset:   2ab0e45260af
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     add option for faces and edges

it will create faces and edges too, and output them to the file
They will not be part of parallel partition

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 5918a09..25cbaa0 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -70,6 +70,7 @@ int main(int argc, char **argv)
   bool quadratic=false;
   bool keep_skins=false;
   bool tetra = false;
+  bool adjEnts = false;
 
   MPI_Init(&argc, &argv);
 
@@ -103,6 +104,8 @@ int main(int argc, char **argv)
   opts.addOpt<void>("tetrahedrons,t", "generate tetrahedrons ",
         &tetra);
 
+  opts.addOpt<void>("faces_edges,f", "create all faces and edges", &adjEnts);
+
   vector<string> intTagNames;
   string firstIntTag;
   opts.addOpt<string>(std::string("int_tag_vert,i"), string("add integer tag on vertices"), &firstIntTag);
@@ -143,6 +146,10 @@ int main(int argc, char **argv)
     return 1;
   }
 
+  if (adjEnts)
+  {
+    keep_skins = true; // do not delete anything
+  }
   // determine m, n, k for processor rank
   int m,n,k;
   k = rank/(M*N);
@@ -473,6 +480,15 @@ int main(int argc, char **argv)
      tt = clock();
   }
 
+  if (adjEnts)
+  {
+    // generate all adj entities dimension 1 and 2 (edges and faces/ tri or quads)
+    Range edges, faces;
+    rval = mb->get_adjacencies(all3dcells, 1, true, edges, Interface::UNION );
+    CHECKE("Can't get edges");
+    rval = mb->get_adjacencies(all3dcells, 2, true, faces, Interface::UNION );
+    CHECKE("Can't get faces");
+  }
   ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0);
   if (pcomm==NULL)
   {


https://bitbucket.org/fathomteam/moab/commits/7b1158dc7db4/
Changeset:   7b1158dc7db4
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     the edges and faces should be added to the partition

if created, it would be better to put the edges and faces to their partitions
for this, create them early, in the loop
Need to update adjacencies, otherwise next time in the loop it does not work
(to create edges and faces based on adjacencies)
it is the n-th time I made this mistake :(

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 25cbaa0..4ec1e1b 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -267,16 +267,23 @@ int main(int argc, char **argv)
 
         EntityHandle starte; // connectivity
         EntityHandle * conn;
+        int num_v_per_elem = 8;
         if (quadratic)
+        {
+          num_v_per_elem = 27;
           rval = iface->get_element_connect(num_el, 27, MBHEX, 0, starte, conn);
+        }
         else if (tetra)
+        {
+          num_v_per_elem = 4;
           rval = iface->get_element_connect(num_el, 4, MBTET, 0, starte, conn);
+        }
         else
           rval = iface->get_element_connect(num_el, 8, MBHEX, 0, starte, conn);
         CHECKE("Can't get element  connectivity.");
 
-        Range hexas(starte, starte+num_el-1); // should be elements
-        // fill  hexas
+        Range cells(starte, starte+num_el-1); // should be elements
+        // fill  cells
         ix=0;
         // identify the elements at the lower corner, for their global ids
         int xe = m*A*blockSize + a*blockSize;
@@ -432,9 +439,28 @@ int main(int argc, char **argv)
         EntityHandle part_set;
         rval = mb->create_meshset(MESHSET_SET, part_set);
         CHECKE("Can't create mesh set.");
-        rval = mb->add_entities(part_set, hexas);
+        rval = mb->add_entities(part_set, cells);
         CHECKE("Can't add entities to set.");
-        rval = mb->tag_set_data(global_id_tag, hexas, &gids[0]);
+        // if needed, add all edges and faces
+        if (adjEnts)
+        {
+          // we need to update adjacencies now, because some elements are new
+          rval = iface->update_adjacencies(starte, num_el, num_v_per_elem, conn);
+          CHECKE("Can't update adjacencies");
+          // generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua
+          Range edges, faces;
+          rval = mb->get_adjacencies(cells, 1, true, edges,
+              Interface::UNION);
+          CHECKE("Can't get edges");
+          rval = mb->get_adjacencies(cells, 2, true, faces,
+              Interface::UNION);
+          CHECKE("Can't get faces");
+          rval = mb->add_entities(part_set, edges);
+          CHECKE("Can't add edges to partition set.");
+          rval = mb->add_entities(part_set, faces);
+          CHECKE("Can't add faces to partition set.");
+        }
+        rval = mb->tag_set_data(global_id_tag, cells, &gids[0]);
         CHECKE("Can't set global ids to elements.");
         int part_num= a +  m*A + (b + n*B)*(M*A) + (c+k*C)*(M*A * N*B);
         rval = mb->tag_set_data(part_tag, &part_set, 1, &part_num);
@@ -480,15 +506,6 @@ int main(int argc, char **argv)
      tt = clock();
   }
 
-  if (adjEnts)
-  {
-    // generate all adj entities dimension 1 and 2 (edges and faces/ tri or quads)
-    Range edges, faces;
-    rval = mb->get_adjacencies(all3dcells, 1, true, edges, Interface::UNION );
-    CHECKE("Can't get edges");
-    rval = mb->get_adjacencies(all3dcells, 2, true, faces, Interface::UNION );
-    CHECKE("Can't get faces");
-  }
   ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0);
   if (pcomm==NULL)
   {


https://bitbucket.org/fathomteam/moab/commits/0f126314d793/
Changeset:   0f126314d793
Branch:      None
User:        iulian07
Date:        2014-07-14 21:19:58
Summary:     skip merge if not needed and skip PCOMM if on one proc

will be able to create meshes faster, if the partitioning
is not needed

Affected #:  1 file

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 4ec1e1b..4338601 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -493,55 +493,58 @@ int main(int argc, char **argv)
   rval = mb->get_entities_by_dimension(0, 0, verts);
 
   CHECKE("Can't get all vertices.");
-
-  if (newMergeMethod)
-    rval = mm.merge_using_integer_tag( verts, global_id_tag);
-  else
-    rval = mm.merge_entities(all3dcells, 0.0001);
-  CHECKE("Can't merge");
-  if (0==rank)
+  if (A*B*C!=1) //  merge needed
   {
-     std::cout << "merge locally:  "
-           << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
-     tt = clock();
+    if (newMergeMethod)
+      rval = mm.merge_using_integer_tag( verts, global_id_tag);
+    else
+      rval = mm.merge_entities(all3dcells, 0.0001);
+    CHECKE("Can't merge");
+    if (0==rank)
+    {
+       std::cout << "merge locally:  "
+             << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+       tt = clock();
+    }
   }
-
-  ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0);
-  if (pcomm==NULL)
+  if (size>1)
   {
-    pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
-  }
-  rval = pcomm->resolve_shared_ents( 0, all3dcells, 3, 0 );
-  CHECKE("Can't resolve shared ents");
+    ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0);
+    if (pcomm==NULL)
+    {
+      pcomm = new ParallelComm( mb, MPI_COMM_WORLD );
+    }
+    rval = pcomm->resolve_shared_ents( 0, all3dcells, 3, 0 );
+    CHECKE("Can't resolve shared ents");
 
-  if (0==rank)
-  {
-     std::cout << "resolve shared entities:  "
-           << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
-     tt = clock();
-  }
+    if (0==rank)
+    {
+       std::cout << "resolve shared entities:  "
+             << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+       tt = clock();
+    }
 
-  if (!keep_skins) // default is to delete the 1- and 2-dimensional entities
-  {
-    // delete all quads and edges
-    Range toDelete;
-    rval =  mb->get_entities_by_dimension(0, 1, toDelete);
-    CHECKE("Can't get edges");
+    if (!keep_skins) // default is to delete the 1- and 2-dimensional entities
+    {
+      // delete all quads and edges
+      Range toDelete;
+      rval =  mb->get_entities_by_dimension(0, 1, toDelete);
+      CHECKE("Can't get edges");
 
-    rval = mb->get_entities_by_dimension(0, 2, toDelete);
-    CHECKE("Can't get faces");
+      rval = mb->get_entities_by_dimension(0, 2, toDelete);
+      CHECKE("Can't get faces");
 
-    rval = pcomm->delete_entities(toDelete) ;
-    CHECKE("Can't delete entities")
+      rval = pcomm->delete_entities(toDelete) ;
+      CHECKE("Can't delete entities")
 
-    if (0==rank)
-    {
-      std::cout << "delete edges and faces, and correct sharedEnts:  "
-            << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
-      tt = clock();
+      if (0==rank)
+      {
+        std::cout << "delete edges and faces, and correct sharedEnts:  "
+              << (clock() - tt) / (double) CLOCKS_PER_SEC << " seconds" << std::endl;
+        tt = clock();
+      }
     }
   }
-
   rval = mb->write_file(outFileName.c_str(), 0, ";;PARALLEL=WRITE_PART");
   CHECKE("Can't write in parallel");
 


https://bitbucket.org/fathomteam/moab/commits/fdfef1b01fe7/
Changeset:   fdfef1b01fe7
Branch:      iulian07/largemesh
User:        iulian07
Date:        2014-07-14 22:39:22
Summary:     remove the is_valid from interface
cosmetic changes to example (too many braces, better definition for q, factor)

Affected #:  3 files

diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp
index 4338601..eb26805 100644
--- a/examples/GenLargeMesh.cpp
+++ b/examples/GenLargeMesh.cpp
@@ -139,17 +139,14 @@ int main(int argc, char **argv)
   if (M*N*K != size)
   {
     if (rank==0)
-    {
       cout <<"M*N*K=" << M*N*K << "  != size = "<< size << "\n";
-    }
     MPI_Finalize();
     return 1;
   }
 
   if (adjEnts)
-  {
     keep_skins = true; // do not delete anything
-  }
+
   // determine m, n, k for processor rank
   int m,n,k;
   k = rank/(M*N);
@@ -157,9 +154,7 @@ int main(int argc, char **argv)
   n = leftover/M;
   m = leftover%M;
   if (rank==size-1)
-  {
     cout << "m, n, k for last rank:" << m << " " << n << " " << k << "\n";
-  }
 
   // so there are a total of M * A * blockSize elements in x direction (so M * A * blockSize + 1 verts in x direction)
   // so there are a total of N * B * blockSize elements in y direction (so N * B * blockSize + 1 verts in y direction)
@@ -171,14 +166,11 @@ int main(int argc, char **argv)
 
   clock_t tt = clock();
 
-  int q = 1;
-  if (quadratic)
-  {
-    q = 2;
-  }
-  int factor =1;
-  if (tetra)
-    factor =6;
+  // used for nodes increments
+  int q = (quadratic)? 2 : 1;
+  // used for element increments
+  int factor = (tetra)? 6 : 1;
+
   int NX = (q * M * A * blockSize + 1);
   int NY = (q * N * B * blockSize + 1);
   int nex = M * A * blockSize; // number of elements in x direction, used for global id on element

diff --git a/src/moab/Interface.hpp b/src/moab/Interface.hpp
index c8102b8..09f4219 100644
--- a/src/moab/Interface.hpp
+++ b/src/moab/Interface.hpp
@@ -1981,8 +1981,6 @@ public:
                                         SetIterator *&set_iter) = 0;
     /**@}*/
 
-
-  virtual bool is_valid(EntityHandle entity) const = 0;
 };
 
 //! predicate for STL algorithms.  Returns true if the entity handle is

diff --git a/src/parallel/ParallelComm.cpp b/src/parallel/ParallelComm.cpp
index 49c1606..f660c0a 100644
--- a/src/parallel/ParallelComm.cpp
+++ b/src/parallel/ParallelComm.cpp
@@ -8853,10 +8853,9 @@ ErrorCode ParallelComm::delete_entities(Range & to_delete)
   std::vector<EntityHandle> good_ents;
   for (size_t j=0; j<sharedEnts.size(); j++)
   {
-    if (mbImpl->is_valid(sharedEnts[j]))
-    {
+    int index=to_delete.index(sharedEnts[j]);
+    if (-1==index)
       good_ents.push_back(sharedEnts[j]);
-    }
   }
   sharedEnts = good_ents;

Repository URL: https://bitbucket.org/fathomteam/moab/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.


More information about the moab-dev mailing list