Add 3D orthogonal structure mesh type for mesh generators.

src/mesh/python/README.md

@@ -0,0 +1,23 @@
+# Mesh File Generation Python Subdirectory (src/mesh/python)
+
+This subdirectory, `src/mesh/python/`, is intended to serve as a location for
+different methods of generating X3D or RTT mesh files.
+The X3D and RTT mesh file formats are supported by C++ parsers in
+`src/RTT_Format_Reader/` and `src/mesh/`
+There is a python module called `mesh_types` with mesh classes, where each
+class has its own `__init__` method for calculating data needed to create X3D
+files (some of which should be usable for RTT file creation).
+For instance, one could add a mesh class that stochastically samples points and
+triangulates them to form a mesh of triangles.
+Currently, there is a script called `x3d_generator.py`, which takes command-line
+input to instantiate a mesh object from one of the mesh_type classes and outputs
+a set of X3D files (main mesh file and boundary node files) of the mesh object.
+
+## Example Usage
+
+To create X3D files for a 4x4x4 3D orthogonal structured mesh on domain (x,y,z)
+in [0,1]x[0,2]x[0,4]:
+
+```bash
+./x3d_generator.py --mesh_type orth_3d_mesh --num_per_dim 4 4 4 --bnd_per_dim 0 1 0 2 0 4
+```

src/mesh/python/mesh_types.py

@@ -121,12 +121,170 @@ def __init__(self, bounds_per_dim, num_cells_per_dim):
 # orthogonal 3D mesh type
 class orth_3d_mesh(base_mesh):
     '''
-    Class for orthogonally structured 2D mesh data.
+    Class for orthogonally structured 3D mesh data.
     This class generates an orthogonally structured mesh in an unstructured
     format suitable for creating unstructured-mesh input files.
     '''
     def __init__(self, bounds_per_dim, num_cells_per_dim):
-        assert (False), 'orth_3d_mesh type not yet implemented...'
+
+        # -- short-cuts
+        nx = num_cells_per_dim[0]
+        ny = num_cells_per_dim[1]
+        nz = num_cells_per_dim[2]
+
+        # -- number of dimensions
+        ndim = len(num_cells_per_dim)
+        self.ndim = ndim
+        assert (ndim == 3), 'ndim != 3, exiting...'
+        assert (len(bounds_per_dim) == ndim), 'len(bounds_per_dim) != ndim, exiting...'
+
+        # -- create grid arrays along each dimension
+        grid_per_dim = [np.linspace(bounds_per_dim[i][0], bounds_per_dim[i][1],
+                                    num_cells_per_dim[i] + 1) for i in range(ndim)]
+
+        # -- create node indices
+        num_nodes = (nx + 1) * (ny + 1) * (nz + 1)
+        self.num_nodes = num_nodes
+        self.coordinates_per_node = np.zeros((num_nodes, ndim))
+        for k in range(nz + 1):
+            k_offset = (nx + 1) * (ny + 1) * k
+            for j in range(ny + 1):
+                j_offset = (nx + 1) * j
+                for i in range(nx + 1):
+                    node = i + j_offset + k_offset
+                    self.coordinates_per_node[node, 0] = grid_per_dim[0][i]
+                    self.coordinates_per_node[node, 1] = grid_per_dim[1][j]
+                    self.coordinates_per_node[node, 2] = grid_per_dim[2][k]
+
+        # -- set total number of cells and faces
+        num_cells = nx * ny * nz
+        self.num_cells = num_cells
+        self.num_faces = 6 * num_cells
+
+        # -- constants for this mesh
+        self.num_faces_per_cell = 6 * np.ones((num_cells), dtype=int)
+        self.num_nodes_per_face = 4 * np.ones((6 * num_cells), dtype=int)
+
+        # -- set nodes per face and faces per cell
+        self.faces_per_cell = np.zeros((num_cells, 6), dtype=int)
+        self.nodes_per_face = np.zeros((6 * num_cells, 4), dtype=int)
+        for k in range(nz):
+            kno = k * (ny + 1) * (nx + 1)
+            knop1 = (k + 1) * (ny + 1) * (nx + 1)
+            for j in range(ny):
+                jno = j * (nx + 1)
+                jnop1 = (j + 1) * (nx + 1)
+                for i in range(nx):
+                    # -- cell index
+                    cell = i + nx * j + nx * ny * k
+                    # -- faces per cell
+                    self.faces_per_cell[cell, 0] = 6 * cell
+                    self.faces_per_cell[cell, 1] = 6 * cell + 1
+                    self.faces_per_cell[cell, 2] = 6 * cell + 2
+                    self.faces_per_cell[cell, 3] = 6 * cell + 3
+                    self.faces_per_cell[cell, 4] = 6 * cell + 4
+                    self.faces_per_cell[cell, 5] = 6 * cell + 5
+                    # -- nodes per face (per cell, counterclockwise around face normals) ...
+                    # ... bottom face (min z)
+                    f = 0
+                    self.nodes_per_face[6 * cell + f, 0] = i + jno + kno
+                    self.nodes_per_face[6 * cell + f, 1] = i + jnop1 + kno
+                    self.nodes_per_face[6 * cell + f, 2] = i + 1 + jnop1 + kno
+                    self.nodes_per_face[6 * cell + f, 3] = i + 1 + jno + kno
+                    # ... top face (max z)
+                    f = 1
+                    self.nodes_per_face[6 * cell + f, 0] = i + 1 + jno + knop1
+                    self.nodes_per_face[6 * cell + f, 1] = i + 1 + jnop1 + knop1
+                    self.nodes_per_face[6 * cell + f, 2] = i + jnop1 + knop1
+                    self.nodes_per_face[6 * cell + f, 3] = i + jno + knop1
+                    # ... down face (min y)
+                    f = 2
+                    self.nodes_per_face[6 * cell + f, 0] = i + jno + kno
+                    self.nodes_per_face[6 * cell + f, 1] = i + 1 + jno + kno
+                    self.nodes_per_face[6 * cell + f, 2] = i + 1 + jno + knop1
+                    self.nodes_per_face[6 * cell + f, 3] = i + jno + knop1
+                    # ... up face (max y)
+                    f = 3
+                    self.nodes_per_face[6 * cell + f, 0] = i + jnop1 + knop1
+                    self.nodes_per_face[6 * cell + f, 1] = i + 1 + jnop1 + knop1
+                    self.nodes_per_face[6 * cell + f, 2] = i + 1 + jnop1 + kno
+                    self.nodes_per_face[6 * cell + f, 3] = i + jnop1 + kno
+                    # ... strange face (min x)
+                    f = 4
+                    self.nodes_per_face[6 * cell + f, 0] = i + jno + kno
+                    self.nodes_per_face[6 * cell + f, 1] = i + jno + knop1
+                    self.nodes_per_face[6 * cell + f, 2] = i + jnop1 + knop1
+                    self.nodes_per_face[6 * cell + f, 3] = i + jnop1 + kno
+                    # ... charm face (max x)
+                    f = 5
+                    self.nodes_per_face[6 * cell + f, 0] = i + 1 + jnop1 + kno
+                    self.nodes_per_face[6 * cell + f, 1] = i + 1 + jnop1 + knop1
+                    self.nodes_per_face[6 * cell + f, 2] = i + 1 + jno + knop1
+                    self.nodes_per_face[6 * cell + f, 3] = i + 1 + jno + kno
+
+        # -- enumerate boundary nodes per face ("side")
+
+        # -- faces at x extrema
+        nodes_per_side_xlow = np.zeros((nz * ny, 4), dtype=int)
+        nodes_per_side_xhig = np.zeros((nz * ny, 4), dtype=int)
+        for k in range(nz):
+            kno = k * (ny + 1) * (nx + 1)
+            knop1 = (k + 1) * (ny + 1) * (nx + 1)
+            for j in range(ny):
+                f = j + ny * k
+                # -- low boundary
+                nodes_per_side_xlow[f, 0] = (nx + 1) * j + kno
+                nodes_per_side_xlow[f, 1] = (nx + 1) * j + knop1
+                nodes_per_side_xlow[f, 2] = (nx + 1) * (j + 1) + knop1
+                nodes_per_side_xlow[f, 3] = (nx + 1) * (j + 1) + kno
+                # -- high boundary
+                nodes_per_side_xhig[f, 0] = nx + (nx + 1) * (j + 1) + kno
+                nodes_per_side_xhig[f, 1] = nx + (nx + 1) * (j + 1) + knop1
+                nodes_per_side_xhig[f, 2] = nx + (nx + 1) * j + knop1
+                nodes_per_side_xhig[f, 3] = nx + (nx + 1) * j + kno
+
+        # -- faces at y extrema
+        nodes_per_side_ylow = np.zeros((nx * nz, 4), dtype=int)
+        nodes_per_side_yhig = np.zeros((nx * nz, 4), dtype=int)
+        for k in range(nz):
+            kno = k * (ny + 1) * (nx + 1)
+            knop1 = (k + 1) * (ny + 1) * (nx + 1)
+            for i in range(nx):
+                f = i + nx * k
+                # -- low boundary
+                nodes_per_side_ylow[f, 0] = i + kno
+                nodes_per_side_ylow[f, 1] = i + 1 + kno
+                nodes_per_side_ylow[f, 2] = i + 1 + knop1
+                nodes_per_side_ylow[f, 3] = i + knop1
+                # -- high boundary
+                nodes_per_side_yhig[f, 0] = i + (nx + 1) * ny + knop1
+                nodes_per_side_yhig[f, 1] = i + 1 + (nx + 1) * ny + knop1
+                nodes_per_side_yhig[f, 2] = i + 1 + (nx + 1) * ny + kno
+                nodes_per_side_yhig[f, 3] = i + (nx + 1) * ny + kno
+
+        # -- faces at z extrema
+        nodes_per_side_zlow = np.zeros((ny * nx, 4), dtype=int)
+        nodes_per_side_zhig = np.zeros((ny * nx, 4), dtype=int)
+        for j in range(ny):
+            jno = j * (nx + 1)
+            jnop1 = (j + 1) * (nx + 1)
+            for i in range(nx):
+                f = i + nx * j
+                # -- low boundary
+                nodes_per_side_zlow[f, 0] = i + jno
+                nodes_per_side_zlow[f, 1] = i + jnop1
+                nodes_per_side_zlow[f, 2] = i + 1 + jnop1
+                nodes_per_side_zlow[f, 3] = i + 1 + jno
+                # -- high boundary
+                nodes_per_side_zhig[f, 0] = i + 1 + jno + nz * (ny + 1) * (nx + 1)
+                nodes_per_side_zhig[f, 1] = i + 1 + jnop1 + nz * (ny + 1) * (nx + 1)
+                nodes_per_side_zhig[f, 2] = i + jnop1 + nz * (ny + 1) * (nx + 1)
+                nodes_per_side_zhig[f, 3] = i + jno + nz * (ny + 1) * (nx + 1)
+
+        # -- compile into one side face array list
+        self.nodes_per_side = [nodes_per_side_xlow, nodes_per_side_xhig,
+                               nodes_per_side_ylow, nodes_per_side_yhig,
+                               nodes_per_side_zlow, nodes_per_side_zhig]
 
 
 # ------------------------------------------------------------------------------------------------ #

src/mesh/python/x3d_generator.py

@@ -10,7 +10,7 @@
 import argparse
 
 # -- mesh class dictionary
-mesh_type_dict = {'orth_2d_mesh': mesh_types.orth_2d_mesh}
+mesh_type_dict = {'orth_2d_mesh': mesh_types.orth_2d_mesh, 'orth_3d_mesh': mesh_types.orth_3d_mesh}
 
 # ------------------------------------------------------------------------------------------------ #
 # -- create argument parser
@@ -22,6 +22,8 @@
                     help='Number of cells per dimension.')
 parser.add_argument('-bd', '--bnd_per_dim', type=float, nargs='+', default=[0.0, 1.0, 0.0, 1.0],
                     help='Length per dimension.')
+parser.add_argument('--name', type=str, default='mesh',
+                    help='Select file name (will be prefixed with x3d. and sufficed with .in).')
 
 # -- parse arguments from command line
 args = parser.parse_args()
@@ -47,7 +49,8 @@
 # -- write out the main mesh file in x3d format
 
 # -- open writable file
-fo = open('x3d.mesh.in', 'w')
+fname = 'x3d.' + args.name
+fo = open(fname + '.in', 'w')
 
 # -- write header block
 # -- for x3d:
@@ -174,7 +177,7 @@
 # -- assume two boundaries per dimension
 # -- x3d only needs a unique node list
 for i in range(2 * mesh.ndim):
-    np.savetxt('x3d.mesh.bdy' + str(i + 1) + '.in', np.unique(mesh.nodes_per_side[i] + 1), fmt='%d')
+    np.savetxt(fname + '.bdy' + str(i + 1) + '.in', np.unique(mesh.nodes_per_side[i] + 1), fmt='%d')
 
 # ------------------------------------------------------------------------------------------------ #
 # end of x3d_generator.py

src/mesh/test/tstX3D_Draco_Mesh_Reader.cc

@@ -235,6 +235,90 @@ void read_voronoi_mesh(rtt_c4::ParallelUnitTest &ut) {
   return;
 }
 
+// Parse an X3D file format and compare to reference
+void read_x3d_mesh_3d(rtt_c4::ParallelUnitTest &ut) {
+
+  // >>> PARSE MESH
+
+  const std::string inputpath = ut.getTestSourcePath();
+  const std::string filename = inputpath + "x3d.mesh3d_1cell.in";
+  const std::vector<std::string> bdy_filenames = {
+      inputpath + "x3d.mesh3d_1cell.bdy1.in", inputpath + "x3d.mesh3d_1cell.bdy2.in",
+      inputpath + "x3d.mesh3d_1cell.bdy3.in", inputpath + "x3d.mesh3d_1cell.bdy4.in",
+      inputpath + "x3d.mesh3d_1cell.bdy5.in", inputpath + "x3d.mesh3d_1cell.bdy6.in"};
+  const std::vector<unsigned> bdy_flags = {3, 1, 0, 2, 1, 5};
+
+  // construct reader
+  std::shared_ptr<X3D_Draco_Mesh_Reader> x3d_reader(
+      new X3D_Draco_Mesh_Reader(filename, bdy_filenames, bdy_flags));
+
+  // read mesh
+  x3d_reader->read_mesh();
+
+  // >>> CHECK HEADER DATA
+
+  FAIL_IF_NOT(x3d_reader->get_process() == 0);
+  FAIL_IF_NOT(x3d_reader->get_numdim() == 3);
+  FAIL_IF_NOT(x3d_reader->get_numcells() == 1);
+  FAIL_IF_NOT(x3d_reader->get_numnodes() == 8);
+
+  // >>> CHECK CELL-NODE DATA
+
+  FAIL_IF_NOT(x3d_reader->get_celltype(0) == 6);
+
+  std::vector<unsigned> test_cellnodes = {0, 2, 3, 1, 5, 7, 6, 4, 0, 1, 5, 4,
+                                          6, 7, 3, 2, 0, 4, 6, 2, 3, 7, 5, 1};
+  FAIL_IF_NOT(x3d_reader->get_cellnodes(0) == test_cellnodes);
+
+  // >>> CHECK NODE-COORD DATA
+
+  std::vector<std::vector<double>> test_coords = {{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0},
+                                                  {1.0, 1.0, 0.0}, {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0},
+                                                  {0.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
+  for (int node = 0; node < 8; ++node)
+    FAIL_IF_NOT(x3d_reader->get_nodecoord(node) == test_coords[node]);
+
+  // >>> CHECK SIDE DATA
+
+  FAIL_IF_NOT(x3d_reader->get_numsides() == 6);
+
+  // set test boundary nodes
+  std::vector<std::vector<unsigned>> test_sidenodes = {{0, 4, 6, 2}, {3, 7, 5, 1}, {0, 1, 5, 4},
+                                                       {6, 7, 3, 2}, {0, 2, 3, 1}, {5, 7, 6, 4}};
+
+  // check each side's data
+  for (int side = 0; side < 6; ++side) {
+
+    // sides must always give 2 nodes per face in X3D
+    FAIL_IF_NOT(x3d_reader->get_sidetype(side) == 4);
+
+    // boundary conditions are not supplied in X3D (note this check is specialized for the 1-cell
+    // mesh)
+    FAIL_IF_NOT(x3d_reader->get_sideflag(side) == bdy_flags[side]);
+
+    // check node indices
+    FAIL_IF_NOT(x3d_reader->get_sidenodes(side) == test_sidenodes[side]);
+  }
+
+  // >>> CHECK BC-NODE MAP
+
+  const std::map<size_t, std::vector<unsigned>> &bc_node_map = x3d_reader->get_bc_node_map();
+
+  FAIL_IF_NOT(bc_node_map.size() == 6);
+
+  std::vector<std::vector<unsigned>> test_bc_nodes = {{0, 2, 4, 6}, {1, 3, 5, 7}, {0, 1, 4, 5},
+                                                      {2, 3, 6, 7}, {0, 1, 2, 3}, {4, 5, 6, 7}};
+
+  for (size_t ibc = 0; ibc < 6; ++ibc) {
+    FAIL_IF_NOT(bc_node_map.at(ibc) == test_bc_nodes[ibc]);
+  }
+
+  // successful test output
+  if (ut.numFails == 0)
+    PASSMSG("3D X3D_Draco_Mesh_Reader parsing tests ok.");
+  return;
+}
+
 //------------------------------------------------------------------------------------------------//
 int main(int argc, char *argv[]) {
   rtt_c4::ParallelUnitTest ut(argc, argv, rtt_dsxx::release);
@@ -243,6 +327,7 @@ int main(int argc, char *argv[]) {
     read_x3d_mesh_2d(ut);
     build_x3d_mesh_2d(ut);
     read_voronoi_mesh(ut);
+    read_x3d_mesh_3d(ut);
   }
   UT_EPILOG(ut);
 }

src/mesh/test/x3d.mesh3d_1cell.bdy1.in

@@ -0,0 +1,4 @@
+1
+3
+5
+7

src/mesh/test/x3d.mesh3d_1cell.bdy2.in

@@ -0,0 +1,4 @@
+2
+4
+6
+8

src/mesh/test/x3d.mesh3d_1cell.bdy3.in

@@ -0,0 +1,4 @@
+1
+2
+5
+6

src/mesh/test/x3d.mesh3d_1cell.bdy4.in

@@ -0,0 +1,4 @@
+3
+4
+7
+8

src/mesh/test/x3d.mesh3d_1cell.bdy5.in

@@ -0,0 +1,4 @@
+1
+2
+3
+4

src/mesh/test/x3d.mesh3d_1cell.bdy6.in

@@ -0,0 +1,4 @@
+5
+6
+7
+8