Source code for custEM.meshgen.vtk_convert

# -*- coding: utf-8 -*-
"""
@author: modified by Rochlitz.R
"""

# ! /usr/bin/env python
#
# Copyright (C) 2006 Anders Logg, modified by Rochlitz.R
# Licensed under the GNU LGPL Version 2.1

import dolfin as df
import logging
import numpy as np
from custEM.meshgen import meshgen_utils as mu

"""
Modified dolfin-convert script for usage with domain markers in mpirun calls.
"""


[docs] def vtk2xml2h5(mesh_name, m_dir, rot, overwrite_markers): """ Converts TetGen meshes in VTK (*.vtk*) format to *xml* and *h5* files for the usage in FEniCS. This function is called automatically during the intialization of the **MOD** class and converts a mesh if it does not already exist as *xml* and *h5* file. Required arguments ------------------ - mesh_name, type str mesh name - m_dir, type str path to mesh directory Functionality ------------- Convert between .vtk and .xml, parser implemented as a state machine: 0 = read number of vertices 1 = read vertices 2 = read number of cells 3 = read cells 4 = read number of markers 5 = read cell markers 6 = done """ logger = logging.getLogger('custEM') logger.info('... converting from VTK format (.vtk) to *xml* and ' '*h5* ...') ifile = open(m_dir + '/_vtk/' + mesh_name + '.vtk', "r") ofile = open(m_dir + '/_xml/' + mesh_name + '.xml', "w") cell_type = "tetrahedron" num_dims = 3 # Step to beginning of file ifile.seek(0) # Write header write_header_mesh(ofile, cell_type, num_dims) # Current state state = 0 # Write data num_vertices_read = 0 num_cells_read = 0 num_markers_read = 0 markers = [] new_mcounter = 2 while 1: # Read next line line = ifile.readline() if state == 0: if line[0] == 'P' and (line[1] == 'O' or line[1] == 'o'): num_vertices = [int(s) for s in line.split() if s.isdigit()][0] write_header_vertices(ofile, num_vertices) state += 1 elif state == 1: coords = line.split() x = float(coords[0]) y = float(coords[1]) z = float(coords[2]) if rot is not None: # rotate pre-rotated mesh back to original rotation if # post_rotate flag during mesh generation is not None if len(rot) == 1: # rotate declination logger.info('... rotating declination ...') xyzrot = mu.rotate_around_point([[x, y, z]], [0., 0., 0.], [np.deg2rad(rot[0])], ['z']) # if len(rot) == 2: # # rotate declination # # logger.info('... rotating declination ...') # xyzrot = mu.rotate_around_point([[x, y, z]], # [0., 0., 0.], # [np.deg2rad(rot[1])], # ['z']) # # rotate inclination # # logger.info('... rotating inclination ...') # xyzrot = mu.rotate_around_point(xyzrot, # [0., 0., 0.], # [np.deg2rad(rot[0])], # ['y']) if len(rot) == 2: # rotate declination # logger.info('... rotating declination ...') xyzrot = mu.rotate_around_point([[x, y, z]], [0., 0., 0.], [np.deg2rad(rot[0])], ['x']) # rotate inclination # logger.info('... rotating inclination ...') xyzrot = mu.rotate_around_point(xyzrot, [0., 0., 0.], [np.deg2rad(rot[1])], ['z']) x, y, z = xyzrot[0, 0], xyzrot[0, 1], xyzrot[0, 2] write_vertex(ofile, num_vertices_read, x, y, z) num_vertices_read += 1 if num_vertices == num_vertices_read: write_footer_vertices(ofile) state += 1 elif state == 2: if line[0] == 'C' and (line[1] == 'E' or line[1] == 'e'): num_cells = [int(s) for s in line.split() if s.isdigit()][0] write_header_cells(ofile, num_cells) state += 1 elif state == 3: nodes = line.split() n0 = int(nodes[1]) n1 = int(nodes[2]) n2 = int(nodes[3]) n3 = int(nodes[4]) write_cell_tetrahedron(ofile, num_cells_read, n0, n1, n2, n3) num_cells_read += 1 if num_cells == num_cells_read: write_footer_cells(ofile) state += 1 elif state == 4: if line[0] == 'L' and (line[1] == 'O' or line[1] == 'o'): write_header_domains(ofile, num_cells) state += 1 elif state == 5: marker = int(line.split()[0]) if overwrite_markers is not None: if marker in overwrite_markers: marker = new_mcounter new_mcounter += 1 if marker == 99999999: logger.info('... Reached Marker 99999999 ...') marker = 1 markers.append(marker) write_domain_marker(ofile, num_markers_read, marker) num_markers_read += 1 if num_cells == num_markers_read: write_footer_domains(ofile) state += 1 elif state == 6: break # Check that we got all data if state == 6: pass else: logger.critical('Error! Missing mesh data, unable to convert. ' 'Aborting ...') raise SystemExit # Write footer write_footer_mesh(ofile) # Close files ifile.close() ofile.close() # Write xml domains file since MeshFunction xml output from one process # is not supported mpirun mode write_xml_domain_file(m_dir + '/_xml/' + mesh_name + '_domains.xml', num_cells, markers) # Convert xml files to h5 format xml_to_hdf5(m_dir + '/_xml/' + mesh_name + '.xml', m_dir + '/_h5/' + mesh_name + '.h5')
[docs] def write_header_vertices(ofile, num_vertices): """ Write header *vertices* to file. """ logger = logging.getLogger('custEM') logger.debug('... expecting %d vertices ...' % num_vertices) ofile.write(" <vertices size=\"%d\">\n" % num_vertices)
[docs] def write_vertex(ofile, vertex, x, y, z): """ Write vertex to file. """ ofile.write(' <vertex index=\"%d\" x=\"%.16g\" y=\"%.16g\" ' 'z=\"%.16g\"/>\n' % (vertex, x, y, z))
[docs] def write_header_mesh(ofile, cell_type, dim): """ Write header of mesh file. """ ofile.write("""\ <?xml version=\"1.0\" encoding=\"UTF-8\"?> <dolfin xmlns:dolfin=\"http://www.fenics.org/dolfin/\"> <mesh celltype="%s" dim="%d"> """ % (cell_type, dim))
[docs] def write_header_cells(ofile, num_cells): """ Write header *cells* to file. """ logger = logging.getLogger('custEM') ofile.write(" <cells size=\"%d\">\n" % num_cells) logger.debug('... expecting %d cells ...' % num_cells)
[docs] def write_cell_tetrahedron(ofile, cell, n0, n1, n2, n3): """ Write tetrahedron 3D information to file. """ ofile.write(' <tetrahedron index=\"%d\" v0=\"%d\" v1=\"%d\" ' 'v2=\"%d\" v3=\"%d\"/>\n' % (cell, n0, n1, n2, n3))
[docs] def write_domain_marker(ofile, cell, marker, dfile=False): """ Write domain markers to domain file. """ if not dfile: ofile.write(' <value cell_index=\"%d\" local_entity=\"%s\" ' 'value=\"%s\"/>\n' % (cell, 0, marker)) elif dfile: ofile.write(' <value cell_index=\"%d\" local_entity=\"%s\" ' 'value=\"%s\"/>\n' % (cell, 0, marker))
[docs] def write_header_domains(ofile, num_domains, dfile=False): """ Write header to domain file. """ logger = logging.getLogger('custEM') if not dfile: logger.debug('... expecting %d domains ...' % num_domains) ofile.write(' <domains>\n') ofile.write(' <mesh_value_collection name="m" type="uint" ' 'dim="2" size=\"%d\">\n' % 0) ofile.write(" </mesh_value_collection>\n") ofile.write(' <mesh_value_collection name="m" type="uint" ' 'dim="3" size=\"%d\">\n' % num_domains) elif dfile: ofile.write('<?xml version="1.0"?>\n') ofile.write('<dolfin xmlns:dolfin="http://fenicsproject.org">\n') ofile.write(' <mesh_function>\n') ofile.write(' <mesh_value_collection name="m" type="uint" ' 'dim="3" size=\"%d\">\n' % num_domains)
[docs] def write_xml_domain_file(ofile_name, num_domains, markers): """ Write domain marker to mesh function file. """ ofile = open(ofile_name, "w") write_header_domains(ofile, num_domains, dfile=True) for j in range(num_domains): write_domain_marker(ofile, j, markers[j], dfile=True) write_footer_domains(ofile, dfile=True) ofile.close()
[docs] def xml_to_hdf5(xml, h5): """ Convert *xml* to HDF5 (*h5*) file. Required arguments ------------------ - xml, type str name of *xml* file - h5, type str name of *h5* file """ mesh = df.Mesh(df.MPI.comm_self, xml) domains = df.MeshFunction("size_t", mesh, 3, mesh.domains()) mesh_file = df.HDF5File(df.MPI.comm_self, h5, 'w') mesh_file.write(mesh, '/mesh') mesh_file.write(domains, "/domains")