diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py index 9477d8779..28ce63b66 100644 --- a/src/mesh/python/mesh_types.py +++ b/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] # ------------------------------------------------------------------------------------------------ # diff --git a/src/mesh/python/x3d_generator.py b/src/mesh/python/x3d_generator.py index aae0e2720..b9488b67f 100755 --- a/src/mesh/python/x3d_generator.py +++ b/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 diff --git a/src/mesh/test/tstX3D_Draco_Mesh_Reader.cc b/src/mesh/test/tstX3D_Draco_Mesh_Reader.cc index 6b13d519a..345ff2e27 100644 --- a/src/mesh/test/tstX3D_Draco_Mesh_Reader.cc +++ b/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 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 bdy_flags = {3, 1, 0, 2, 1, 5}; + + // construct reader + std::shared_ptr 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 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> 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> 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> &bc_node_map = x3d_reader->get_bc_node_map(); + + FAIL_IF_NOT(bc_node_map.size() == 6); + + std::vector> 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); } diff --git a/src/mesh/test/x3d.mesh3d_1cell.bdy1.in b/src/mesh/test/x3d.mesh3d_1cell.bdy1.in new file mode 100644 index 000000000..9b4237ab1 --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.bdy1.in @@ -0,0 +1,4 @@ +1 +3 +5 +7 diff --git a/src/mesh/test/x3d.mesh3d_1cell.bdy2.in b/src/mesh/test/x3d.mesh3d_1cell.bdy2.in new file mode 100644 index 000000000..628d2a8fd --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.bdy2.in @@ -0,0 +1,4 @@ +2 +4 +6 +8 diff --git a/src/mesh/test/x3d.mesh3d_1cell.bdy3.in b/src/mesh/test/x3d.mesh3d_1cell.bdy3.in new file mode 100644 index 000000000..785c59f29 --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.bdy3.in @@ -0,0 +1,4 @@ +1 +2 +5 +6 diff --git a/src/mesh/test/x3d.mesh3d_1cell.bdy4.in b/src/mesh/test/x3d.mesh3d_1cell.bdy4.in new file mode 100644 index 000000000..145411447 --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.bdy4.in @@ -0,0 +1,4 @@ +3 +4 +7 +8 diff --git a/src/mesh/test/x3d.mesh3d_1cell.bdy5.in b/src/mesh/test/x3d.mesh3d_1cell.bdy5.in new file mode 100644 index 000000000..94ebaf900 --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.bdy5.in @@ -0,0 +1,4 @@ +1 +2 +3 +4 diff --git a/src/mesh/test/x3d.mesh3d_1cell.bdy6.in b/src/mesh/test/x3d.mesh3d_1cell.bdy6.in new file mode 100644 index 000000000..3aac70f84 --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.bdy6.in @@ -0,0 +1,4 @@ +5 +6 +7 +8 diff --git a/src/mesh/test/x3d.mesh3d_1cell.in b/src/mesh/test/x3d.mesh3d_1cell.in new file mode 100644 index 000000000..75fa81422 --- /dev/null +++ b/src/mesh/test/x3d.mesh3d_1cell.in @@ -0,0 +1,66 @@ +ascii +header + process 1 + numdim 3 + materials 1 + nodes 8 + faces 6 + elements 1 + ghost_nodes 0 + slaved_nodes 0 + nodes_per_slave 2 + nodes_per_face 4 + faces_per_cell 6 + node_data_fields 0 + cell_data_fields 2 +end_header + +matnames + 1 0 +end_matnames + +mateos + 1 -1 +end_mateos + +matopc + 1 -1 +end_matopc + +nodes + 1 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 2 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 3 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 4 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 5 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 + 6 1.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 + 7 0.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 + 8 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 +end_nodes + +faces + 1 4 1 3 4 2 + 2 4 6 8 7 5 + 3 4 1 2 6 5 + 4 4 7 8 4 3 + 5 4 1 5 7 3 + 6 4 4 8 6 2 +end_faces + +cells + 1 6 1 2 3 4 5 6 +end_cells + +slaved_nodes 0 +end_slaved_nodes +ghost_nodes 0 +end_ghost_nodes + +cell_data +matid + 0 +end_matid +partelm + 1 +end_partelm +end_cell_data