diff --git a/src/mesh/python/mesh_types.py b/src/mesh/python/mesh_types.py new file mode 100644 index 000000000..9477d8779 --- /dev/null +++ b/src/mesh/python/mesh_types.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python +# -------------------------------------------*-python-*------------------------------------------- # +# file src/mesh/python/mesh_types.py +# date Monday, Jul 19, 2021, 12:14 pm +# brief This script provides mesh classes that calculate and contain unstructred mesh data. +# note Copyright (C) 2021, Triad National Security, LLC., All rights reserved. +# ------------------------------------------------------------------------------------------------ # +import numpy as np + + +# ------------------------------------------------------------------------------------------------ # +# base class (to show required member data) +class base_mesh: + ''' + Base class with data required by mesh file generators. + This contains no methods for creating the data members. + ''' + def __init__(self): + # -- required data + self.ndim = 0 # number of dimensions + self.num_nodes = 0 # total number of nodes + self.coordinates_per_node = np.array([]) # coordinate array indexed by node + self.num_cells = 0 # total number of cells + self.num_faces = 0 # total number of oriented faces + self.num_faces_per_cell = np.array([], dtype=int) # number of faces per cell + self.num_nodes_per_face = np.array([], dtype=int) # number of nodes per face + self.faces_per_cell = np.array([[]], dtype=int) # face indexes per cell + self.nodes_per_face = np.array([[]], dtype=int) # node indexes per face + self.nodes_per_side = [np.array([[]], dtype=int)] # list of arrays of node per bdy face + + +# ------------------------------------------------------------------------------------------------ # +# orthogonal 2D mesh type +class orth_2d_mesh(base_mesh): + ''' + Class for orthogonally structured 2D 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): + + # -- short-cuts + nx = num_cells_per_dim[0] + ny = num_cells_per_dim[1] + + # -- number of dimensions + ndim = len(num_cells_per_dim) + self.ndim = ndim + assert (ndim == 2), 'ndim != 2, 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) + self.num_nodes = num_nodes + self.coordinates_per_node = np.zeros((num_nodes, ndim)) + for j in range(ny + 1): + j_offset = (nx + 1) * j + for i in range(nx + 1): + node = i + j_offset + self.coordinates_per_node[node, 0] = grid_per_dim[0][i] + self.coordinates_per_node[node, 1] = grid_per_dim[1][j] + + # -- set total number of cells and faces + num_cells = nx * ny + self.num_cells = num_cells + self.num_faces = 4 * num_cells + + # -- constants for this mesh + self.num_faces_per_cell = 4 * np.ones((num_cells), dtype=int) + self.num_nodes_per_face = 2 * np.ones((4 * num_cells), dtype=int) + + # -- set nodes per face and faces per cell + self.faces_per_cell = np.zeros((num_cells, 4), dtype=int) + self.nodes_per_face = np.zeros((4 * num_cells, 2), dtype=int) + for j in range(ny): + for i in range(nx): + # -- cell index + cell = i + nx * j + # -- faces per cell + self.faces_per_cell[cell, 0] = 4 * cell + self.faces_per_cell[cell, 1] = 4 * cell + 1 + self.faces_per_cell[cell, 2] = 4 * cell + 2 + self.faces_per_cell[cell, 3] = 4 * cell + 3 + # -- nodes per face (per cell) + self.nodes_per_face[4 * cell, 0] = i + j * (nx + 1) + self.nodes_per_face[4 * cell, 1] = i + 1 + j * (nx + 1) + self.nodes_per_face[4 * cell + 1, 0] = i + 1 + j * (nx + 1) + self.nodes_per_face[4 * cell + 1, 1] = i + 1 + (j + 1) * (nx + 1) + self.nodes_per_face[4 * cell + 2, 0] = i + 1 + (j + 1) * (nx + 1) + self.nodes_per_face[4 * cell + 2, 1] = i + (j + 1) * (nx + 1) + self.nodes_per_face[4 * cell + 3, 0] = i + (j + 1) * (nx + 1) + self.nodes_per_face[4 * cell + 3, 1] = i + j * (nx + 1) + + # -- enumerate boundary nodes per face ("side") + # -- faces at x extrema + nodes_per_side_xlow = np.zeros((ny, 2), dtype=int) + nodes_per_side_xhig = np.zeros((ny, 2), dtype=int) + for j in range(ny): + nodes_per_side_xlow[ny - 1 - j, 0] = (j + 1) * (nx + 1) + nodes_per_side_xlow[ny - 1 - j, 1] = j * (nx + 1) + nodes_per_side_xhig[j, 0] = j * (nx + 1) + nx + nodes_per_side_xhig[j, 1] = (j + 1) * (nx + 1) + nx + # -- faces at y extrema + nodes_per_side_ylow = np.zeros((nx, 2), dtype=int) + nodes_per_side_yhig = np.zeros((nx, 2), dtype=int) + for i in range(nx): + nodes_per_side_ylow[i, 0] = i + nodes_per_side_ylow[i, 1] = i + 1 + nodes_per_side_yhig[nx - 1 - i, 0] = i + 1 + ny * (nx + 1) + nodes_per_side_yhig[nx - 1 - i, 1] = i + ny * (nx + 1) + # -- compile into one side face array list (expeced as counter-clockwise in 2D) + self.nodes_per_side = [nodes_per_side_ylow, nodes_per_side_xhig, + nodes_per_side_yhig, nodes_per_side_xlow] + + +# ------------------------------------------------------------------------------------------------ # +# orthogonal 3D mesh type +class orth_3d_mesh(base_mesh): + ''' + Class for orthogonally structured 2D 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...' + + +# ------------------------------------------------------------------------------------------------ # +# end of mesh_types.py +# ------------------------------------------------------------------------------------------------ # diff --git a/src/mesh/python/x3d_generator.py b/src/mesh/python/x3d_generator.py new file mode 100755 index 000000000..aae0e2720 --- /dev/null +++ b/src/mesh/python/x3d_generator.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python +# -------------------------------------------*-python-*------------------------------------------- # +# file src/mesh/python/x3d_generator.py +# date Monday, Jul 19, 2021, 12:14 pm +# brief This script generates X3D mesh files from a mesh object (assumed to have certain data). +# note Copyright (C) 2021, Triad National Security, LLC., All rights reserved. +# ------------------------------------------------------------------------------------------------ # +import mesh_types +import numpy as np +import argparse + +# -- mesh class dictionary +mesh_type_dict = {'orth_2d_mesh': mesh_types.orth_2d_mesh} + +# ------------------------------------------------------------------------------------------------ # +# -- create argument parser + +parser = argparse.ArgumentParser(description='Generate X3D mesh file.') +parser.add_argument('-mt', '--mesh_type', type=str, default='orth_2d_mesh', + help='Select mesh type.') +parser.add_argument('-nd', '--num_per_dim', type=int, nargs='+', default=[1, 1], + 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.') + +# -- parse arguments from command line +args = parser.parse_args() + +# ------------------------------------------------------------------------------------------------ # +# -- create the relevant mesh object + +# -- number of spatial dimensions should be length of num_per_dim +ndim = len(args.num_per_dim) + +# -- sanity check input +assert (ndim > 0), 'len(args.num_per_dim) <= 0' +assert (ndim < 4), 'len(args.num_per_dim) >= 4' +assert (len(args.bnd_per_dim) == 2 * ndim), 'len(args.bnd_per_dim) != 2 * ndim' + +# -- de-serialize spatial bound list +bnd_per_dim = [[args.bnd_per_dim[2 * i], args.bnd_per_dim[2 * i + 1]] for i in range(ndim)] + +# -- instantiate the class for the mesh type selected +mesh = mesh_type_dict[args.mesh_type](bnd_per_dim, args.num_per_dim) + +# ------------------------------------------------------------------------------------------------ # +# -- write out the main mesh file in x3d format + +# -- open writable file +fo = open('x3d.mesh.in', 'w') + +# -- write header block +# -- for x3d: +# -- nodes_per_face = max number of nodes per face +# -- faces_per_cell = max number of faces per cell +fo.write('ascii\n') +fo.write('header\n') +fo.write(' process 1\n') +fo.write(' numdim {0}\n'.format(mesh.ndim)) +fo.write(' materials 1\n') +fo.write(' nodes {0}\n'.format(mesh.num_nodes)) +fo.write(' faces {0}\n'.format(mesh.num_faces)) +fo.write(' elements {0}\n'.format(mesh.num_cells)) +fo.write(' ghost_nodes 0\n') +fo.write(' slaved_nodes 0\n') +fo.write(' nodes_per_slave 2\n') +fo.write(' nodes_per_face {0}\n'.format(np.max(mesh.num_nodes_per_face))) +fo.write(' faces_per_cell {0}\n'.format(np.max(mesh.num_faces_per_cell))) +fo.write(' node_data_fields 0\n') +fo.write(' cell_data_fields 2\n') +fo.write('end_header\n') + +# -- space +fo.write('\n') + +# -- write trivial material block +fo.write('matnames\n') +fo.write(' 1 0\n') +fo.write('end_matnames\n') + +# -- space +fo.write('\n') + +# -- write trivial material eos block +fo.write('mateos\n') +fo.write(' 1 -1\n') +fo.write('end_mateos\n') + +# -- space +fo.write('\n') + +# -- write trival material opc block +fo.write('matopc\n') +fo.write(' 1 -1\n') +fo.write('end_matopc\n') + +# -- space +fo.write('\n') + +# -- write node coordinate block +fo.write('nodes\n') +for node in range(mesh.num_nodes): + crds = [0.0, 0.0, 0.0] + for idim in range(mesh.ndim): + crds[idim] = mesh.coordinates_per_node[node, idim] + fo.write('{0:10d} {1:.15e} {2:.15e} {3:.15e}\n'.format(node + 1, crds[0], crds[1], crds[2])) +fo.write('end_nodes\n') + +# -- space +fo.write('\n') + +# -- write face block (excludes extra columns used for parallel meshes, for now) +fo.write('faces\n') +for cell in range(mesh.num_cells): + for j in range(mesh.num_faces_per_cell[cell]): + # -- get the full face index for the cell and local face index + face = mesh.faces_per_cell[cell, j] + # -- get the number of nodes for this face + nnpf = mesh.num_nodes_per_face[face] + # -- build face string + face_str = '{0:10d}{1:10d}'.format(face + 1, nnpf) + # -- append node indices to face string + for node in mesh.nodes_per_face[face, :]: + face_str += '{0:10d}'.format(node + 1) + # -- write face line + fo.write(face_str + '\n') +fo.write('end_faces\n') + +# -- space +fo.write('\n') + +# -- write cell block +fo.write('cells\n') +for cell in range(mesh.num_cells): + # -- get the number of faces for this cell + nfpc = mesh.num_faces_per_cell[cell] + # -- build cell string + cell_str = '{0:10d}{1:10d}'.format(cell + 1, nfpc) + # -- append face indices to cell string + for face in mesh.faces_per_cell[cell, :]: + cell_str += '{0:10d}'.format(face + 1) + # -- write cell line + fo.write(cell_str + '\n') +fo.write('end_cells\n') + +# -- space +fo.write('\n') + +# -- unused AMR-type/parallel data +fo.write('slaved_nodes 0\n') +fo.write('end_slaved_nodes\n') +fo.write('ghost_nodes 0\n') +fo.write('end_ghost_nodes\n') + +# -- space +fo.write('\n') + +# -- unused data per cell +fo.write('cell_data\n') +fo.write('matid\n') +fo.write(' 0\n') +fo.write('end_matid\n') +fo.write('partelm\n') +fo.write(' 1\n') +fo.write('end_partelm\n') +fo.write('end_cell_data\n') + +# -- close main mesh file +fo.close() + +# ------------------------------------------------------------------------------------------------ # +# -- write out the mesh boundary node lists to file (for boundary conditions) + +# -- 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') + +# ------------------------------------------------------------------------------------------------ # +# end of x3d_generator.py +# ------------------------------------------------------------------------------------------------ # diff --git a/src/mesh/test/x3d.mesh.bdy3.in b/src/mesh/test/x3d.mesh.bdy3.in index 81450f3d3..b94473479 100644 --- a/src/mesh/test/x3d.mesh.bdy3.in +++ b/src/mesh/test/x3d.mesh.bdy3.in @@ -1,2 +1,2 @@ -4 3 +4 diff --git a/src/mesh/test/x3d.mesh.bdy4.in b/src/mesh/test/x3d.mesh.bdy4.in index f00580c40..2b2f2e1b9 100644 --- a/src/mesh/test/x3d.mesh.bdy4.in +++ b/src/mesh/test/x3d.mesh.bdy4.in @@ -1,2 +1,2 @@ -3 1 +3 diff --git a/src/mesh/test/x3d.mesh.in b/src/mesh/test/x3d.mesh.in index 21f33e8ad..6d36c7553 100644 --- a/src/mesh/test/x3d.mesh.in +++ b/src/mesh/test/x3d.mesh.in @@ -28,17 +28,17 @@ matopc 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 + 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 end_nodes faces - 1 2 1 2 1 0 0 1 1 1 1 1 - 2 2 2 4 1 0 3 1 1 1 1 1 - 3 2 4 3 1 0 2 1 1 1 1 1 - 4 2 3 1 1 0 5 1 1 1 1 1 + 1 2 1 2 + 2 2 2 4 + 3 2 4 3 + 4 2 3 1 end_faces cells