Add single-plane reflection and periodicity capability to native mesh…

… interface (OpenFUSIONToolkit#9)

Add single-plane reflection and periodicity capability to native mesh interface
 - Add tests for single-plane reflection and periodicity with linear and quadratic Hex elements from CUBIT grids
 - Update CUBIT conversion script to handle nodesets for periodicity
 - Add support for conversion of CUBIT files with lower-dimensional information (eg. CAD linkage)

src/grid/multigrid_build.F90

@@ -26,7 +26,7 @@ module multigrid_build
   multigrid_level, trimesh_mg_globals, quadmesh_mg_globals, multigrid_reffix_ho, &
   multigrid_reffix_ho_surf
 use oft_mesh_native, only: native_load_mesh, native_load_smesh, mesh_native_id, &
-  native_hobase
+  native_hobase, native_set_periodic
 use oft_mesh_t3d, only: mesh_t3d_load, mesh_t3d_cadsync, mesh_t3d_cadlink, &
   mesh_t3d_add_quad, mesh_t3d_reffix, mesh_t3d_add_quad, &
   mesh_t3d_set_periodic, smesh_t3d_load, mesh_t3d_id
@@ -68,6 +68,7 @@ subroutine multigrid_load(cad_type)
     CALL native_load_mesh
     CALL mesh_global_init(mesh)
     CALL native_hobase(mesh)
+    CALL native_set_periodic
   case(mesh_t3d_id) ! T3D Mesh
     CALL mesh_t3d_load
     CALL mesh_global_init(mesh)
@@ -750,6 +751,7 @@ subroutine multigrid_load_surf(cad_type)
     CALL native_load_smesh
     CALL smesh_global_init(smesh)
     CALL native_hobase(smesh)
+    ! CALL native_set_periodic
   case(mesh_t3d_id) ! T3D Mesh
     CALL smesh_t3d_load
     CALL smesh_global_init(smesh)

src/tests/grid/CMakeLists.txt

@@ -23,7 +23,13 @@ if( OFT_BUILD_TESTS )
   file( COPY test_cubit.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
   file( COPY ref_test.3dm DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
   file( COPY ref_tet4_test.in.e DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+  file( COPY ref_tet4_test.h5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
   file( COPY ref_tet10_test.in.e DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+  file( COPY ref_tet10_test.h5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+  file( COPY ref_hex8_test.in.e DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+  file( COPY ref_hex8_test.h5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+  file( COPY ref_hex27_test.in.e DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+  file( COPY ref_hex27_test.h5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
   file( COPY sphere_test.3dm DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
   file( COPY sphere_tet4_test.in.e DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
   file( COPY sphere_tet4_test.h5 DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )

src/tests/grid/ref_test.jou

@@ -12,10 +12,33 @@ mesh volume 1
 #
 nodeset 1 surface 3
 #
-refine parallel fileroot 'ref_test' overwrite no_execute
+#refine parallel fileroot 'ref_test' overwrite no_execute
+export Genesis "ref_tet4_test.g" dimension 3 overwrite block all
 #
 block 1 volume 1
 block 1 element type tetra10
 set large exodus file on
-export Genesis  "ref_test2.in.e" block 1 overwrite
+export Genesis "ref_tet10_test.g" dimension 3 overwrite block all
 
+reset
+#
+create Cylinder height 0.5 radius 1 
+move Volume 1  x 0 y 0 z 0.25 include_merged
+#
+surface 2  scheme circle
+surface 2 size .5
+mesh surface 2
+volume 1  redistribute nodes off 
+volume 1  scheme Sweep sweep transform least squares 
+volume 1  autosmooth target on  fixed imprints off  smart smooth off 
+volume 1 size .5
+mesh volume 1
+#
+nodeset 1 surface 3
+#
+export Genesis "ref_hex8_test.g" dimension 3 overwrite block all
+#
+block 1 volume 1
+block 1 element type hex27
+set large exodus file on
+export Genesis  "ref_hex27_test.g" dimension 3 overwrite block all

src/tests/grid/test_cubit.py

@@ -38,6 +38,8 @@
 
 &native_mesh_options
  filename='{2}.h5'
+ reflect={5}
+ ref_periodic={10}
 /
 """
 
@@ -51,11 +53,14 @@ def cubit_setup(nbase, nlevels, prefix, inp_prefix=None, grid_order=1,
     nproc = 1
     if nbase != nlevels:
         nproc = 2
+    ref_periodic = 'F'
+    if per_ns > 0:
+        ref_periodic = 'T'
     #
     os.chdir(test_dir)
     with open('oft.in', 'w+') as fid:
         fid.write(oft_in_template.format(nbase, nlevels, prefix, prefix2,
-                                       grid_order, reflect, per_ns, zstretch, test_2d, cad_type))
+                                       grid_order, reflect, per_ns, zstretch, test_2d, cad_type, ref_periodic))
     return run_OFT("./test_cubit", nproc, 60)
 
 # Validate results against expected values
@@ -161,10 +166,18 @@ def test_cut_1ref(top_lev):
 #============================================================================
 # Test runners for reflection
 @pytest.mark.parametrize("top_lev", (1, 2))
-def test_reflect_base(top_lev):
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_reflect_base(top_lev,cad_type):
     volume_cubit = 3.02070
     area_cubit =  12.26361
-    assert cubit_setup(1,top_lev,'ref_tet4_test',reflect='T')
+    assert cubit_setup(1,top_lev,'ref_tet4_test',reflect='T',cad_type=cad_type)
+    assert check_result(volume_cubit, area_cubit)
+@pytest.mark.parametrize("top_lev", (1, 2))
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_reflect_hex_base(top_lev,cad_type):
+    volume_cubit = 3.00000
+    area_cubit =  12.21166
+    assert cubit_setup(1,top_lev,'ref_hex8_test',reflect='T',cad_type=cad_type)
     assert check_result(volume_cubit, area_cubit)
 @pytest.mark.parametrize("top_lev", (2, 3))
 def test_reflect_1ref(top_lev):
@@ -183,10 +196,18 @@ def test_reflect_quad(top_lev):
 #============================================================================
 # Test runners for periodic reflection
 @pytest.mark.parametrize("top_lev", (1, 2))
-def test_perreflect_base(top_lev):
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_perreflect_base(top_lev,cad_type):
     volume_cubit = 3.02070
     area_cubit =   6.22221
-    assert cubit_setup(1,top_lev,'ref_tet4_test',reflect='T',per_ns=1)
+    assert cubit_setup(1,top_lev,'ref_tet4_test',reflect='T',per_ns=1,cad_type=cad_type)
+    assert check_result(volume_cubit, area_cubit)
+@pytest.mark.parametrize("top_lev", (1, 2))
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_perreflect_hex_base(top_lev,cad_type):
+    volume_cubit = 3.00000
+    area_cubit =   6.21166
+    assert cubit_setup(1,top_lev,'ref_hex8_test',reflect='T',per_ns=1,cad_type=cad_type)
     assert check_result(volume_cubit, area_cubit)
 @pytest.mark.parametrize("top_lev", (2, 3))
 def test_perreflect_1ref(top_lev):
@@ -205,17 +226,34 @@ def test_perreflect_quad(top_lev):
 #============================================================================
 # Test runners for high order reflected meshes
 @pytest.mark.parametrize("top_lev", (1, 2))
-def test_tet10reflect_quad(top_lev):
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_tet10reflect_quad(top_lev,cad_type):
     volume_cubit = 3.14123
     area_cubit =  12.56531
-    assert cubit_setup(1,top_lev,'ref_tet10_test',grid_order=2,reflect='T')
+    assert cubit_setup(1,top_lev,'ref_tet10_test',grid_order=2,reflect='T',cad_type=cad_type)
     assert check_result(volume_cubit, area_cubit)
 @pytest.mark.parametrize("top_lev", (2, 3))
-def test_tet10reflect_1ref(top_lev):
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_tet10reflect_1ref(top_lev,cad_type):
     volume_cubit = 3.11110
     area_cubit =  12.49011
     minlev = 4 - top_lev
-    assert cubit_setup(minlev,top_lev,'ref_tet10_test',reflect='T')
+    assert cubit_setup(minlev,top_lev,'ref_tet10_test',reflect='T',cad_type=cad_type)
+    assert check_result(volume_cubit, area_cubit)
+@pytest.mark.parametrize("top_lev", (1, 2))
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_hex27reflect_quad(top_lev,cad_type):
+    volume_cubit = 3.14110
+    area_cubit =  12.56491
+    assert cubit_setup(1,top_lev,'ref_hex27_test',grid_order=2,reflect='T',cad_type=cad_type)
+    assert check_result(volume_cubit, area_cubit)
+@pytest.mark.parametrize("top_lev", (2, 3))
+@pytest.mark.parametrize("cad_type", (0, 2))
+def test_hex27reflect_1ref(top_lev,cad_type):
+    volume_cubit = 3.10583
+    area_cubit =  12.47691
+    minlev = 4 - top_lev
+    assert cubit_setup(minlev,top_lev,'ref_hex27_test',reflect='T',cad_type=cad_type)
     assert check_result(volume_cubit, area_cubit)
 
 #============================================================================

src/utilities/convert_cubit.py

@@ -47,14 +47,14 @@
     [2,3,7,6]
 ])
 element_type_map = {
-    'TRI3': {'ncp': 3, 'ncp_lin': 3, 'type': 'TRI_p1', 'ed_map': tri_ed_map},
-    'TRI6': {'ncp': 6, 'ncp_lin': 3, 'type': 'TRI_p2', 'ed_map': tri_ed_map},
-    'QUAD4': {'ncp': 4, 'ncp_lin': 4, 'type': 'QUAD_p1', 'ed_map': quad_ed_map},
-    'QUAD9': {'ncp': 9, 'ncp_lin': 4, 'type': 'QUAD_p2', 'ed_map': quad_ed_map, 'cellpoint': True},
-    'TETRA4': {'ncp': 4, 'ncp_lin': 4, 'type': 'TET_p1', 'ed_map': tet_ed_map},
-    'TETRA10': {'ncp': 10, 'ncp_lin': 4, 'type': 'TET_p2', 'ed_map': tet_ed_map},
-    'HEX8': {'ncp': 8, 'ncp_lin': 8, 'type': 'HEX_p1', 'ed_map': hex_ed_map, 'face_map': hex_face_map},
-    'HEX27': {'ncp': 27, 'ncp_lin': 8, 'type': 'HEX_p2', 'ed_map': hex_ed_map, 'face_map': hex_face_map, 'cellpoint': True}
+    'TRI3': {'dim': 2, 'ncp': 3, 'ncp_lin': 3, 'type': 'TRI_p1', 'ed_map': tri_ed_map},
+    'TRI6': {'dim': 2, 'ncp': 6, 'ncp_lin': 3, 'type': 'TRI_p2', 'ed_map': tri_ed_map},
+    'QUAD4': {'dim': 2, 'ncp': 4, 'ncp_lin': 4, 'type': 'QUAD_p1', 'ed_map': quad_ed_map},
+    'QUAD9': {'dim': 2, 'ncp': 9, 'ncp_lin': 4, 'type': 'QUAD_p2', 'ed_map': quad_ed_map, 'cellpoint': True},
+    'TETRA4': {'dim': 3, 'ncp': 4, 'ncp_lin': 4, 'type': 'TET_p1', 'ed_map': tet_ed_map},
+    'TETRA10': {'dim': 3, 'ncp': 10, 'ncp_lin': 4, 'type': 'TET_p2', 'ed_map': tet_ed_map},
+    'HEX8': {'dim': 3, 'ncp': 8, 'ncp_lin': 8, 'type': 'HEX_p1', 'ed_map': hex_ed_map, 'face_map': hex_face_map},
+    'HEX27': {'dim': 3, 'ncp': 27, 'ncp_lin': 8, 'type': 'HEX_p2', 'ed_map': hex_ed_map, 'face_map': hex_face_map, 'cellpoint': True}
 }
 element_type_map['TRI'] = element_type_map['TRI3']
 element_type_map['TETRA'] = element_type_map['TETRA4']
@@ -79,24 +79,36 @@ def read_mesh(filename):
         r = np.transpose(np.vstack(r)).copy()
     # Read regions
     lc = []
-    reg = []
     node_sets = []
     side_sets = []
     block_types = []
+    max_logical_dim = 0
     for varname, variable in ncdf_file.variables.items():
         if varname.startswith('connect'):
             for attrname in variable.ncattrs():
                 if attrname.startswith('elem_type'):
-                    block_types.append(element_type_map.get(getattr(variable, attrname),{'type': 'unknown'}))
+                    block_type = getattr(variable, attrname)
+                    block_type_info = element_type_map.get(block_type,{'type': block_type, 'dim': -1})
+                    max_logical_dim = max(max_logical_dim,block_type_info['dim'])
+                    block_types.append(block_type_info)
             lc_tmp = np.asarray(variable)
             lc.append(lc_tmp)
-            nReg = len(lc)
-            reg.append(nReg*np.ones((lc_tmp.shape[0],), dtype=np.int32))
         elif varname.startswith('node_ns'):
             node_sets.append(np.asarray(variable))
         elif varname.startswith('elem_ss'):
             side_sets.append(np.asarray(variable))
-    nReg = len(lc)
+    # Remove lower level geometry
+    keep_inds = []
+    reg = []
+    for i, block_type in enumerate(block_types):
+        if block_type['dim'] < max_logical_dim:
+            print("  Note: Removing block {0} of type {1}".format(i+1,block_type['type']))
+            continue
+        keep_inds.append(i)
+        nReg = len(keep_inds)
+        reg.append(nReg*np.ones((lc_tmp.shape[0],), dtype=np.int32))
+    lc = [lc[i] for i in keep_inds]
+    block_types = [block_types[i] for i in keep_inds]
     lc = np.vstack(lc)
     reg = np.hstack(reg)
     mesh_order = 1
@@ -167,7 +179,7 @@ def read_mesh(filename):
     #
     return r_new, lc_new, reg, node_sets, side_sets, ho_info
 
-def write_file(filename, r, lc, reg, node_sets=[], side_sets=[], ho_info=None):
+def write_file(filename, r, lc, reg, node_sets=[], side_sets=[], ho_info=None, periodic_info=None):
     print()
     print("Saving mesh: {0}".format(filename))
     h5_file = h5py.File(filename, 'w')
@@ -187,17 +199,27 @@ def write_file(filename, r, lc, reg, node_sets=[], side_sets=[], ho_info=None):
         h5_file.create_dataset('mesh/NUM_SIDESETS', data=[len(side_sets),], dtype='i4')
         for i, side_set in enumerate(side_sets):
             h5_file.create_dataset('mesh/SIDESET{0:04d}'.format(i+1), data=side_set, dtype='i4')
+    if periodic_info is not None:
+        h5_file.create_dataset('mesh/periodicity/nodes', data=periodic_info, dtype='i4')
 
 
 parser = argparse.ArgumentParser()
 parser.description = "Pre-processing script for mesh files"
 parser.add_argument("--in_file", type=str, required=True, help="Input mesh file")
 parser.add_argument("--out_file", type=str, default=None, help="Ouput mesh file")
+parser.add_argument("--periodic_nodeset", type=int, default=None, help="Index of perioidic nodeset")
 options = parser.parse_args()
 
 out_file = options.out_file
 if out_file is None:
     out_file = os.path.splitext(options.in_file)[0] + ".h5"
 
 r, lc, reg, node_sets, side_sets, ho_info = read_mesh(options.in_file)
-write_file(out_file, r, lc, reg, node_sets, side_sets, ho_info)
+
+periodic_info = None
+if options.periodic_nodeset is not None:
+    if options.periodic_nodeset > len(node_sets):
+        raise ValueError("Periodic nodeset ({0}) is out of bounds ({1})".format(options.periodic_nodeset, len(node_sets)))
+    periodic_info = node_sets.pop(options.periodic_nodeset-1)
+
+write_file(out_file, r, lc, reg, node_sets, side_sets, ho_info, periodic_info)