[MOAB-dev] commit/MOAB: iulian07: add quadratic option

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Mon Jul 7 00:02:32 CDT 2014


1 new commit in MOAB:

https://bitbucket.org/fathomteam/moab/commits/6dab2094d347/
Changeset:   6dab2094d347
Branch:      iulian07/largemesh
User:        iulian07
Date:        2014-07-07 07:00:20
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;

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