# -*- 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_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_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_footer_domains(ofile, dfile=False):
"""
Write footer to domain file.
"""
logger = logging.getLogger('custEM')
if not dfile:
ofile.write(" </mesh_value_collection>\n")
ofile.write(" </domains>\n")
logger.debug('... found all domains ...')
elif dfile:
ofile.write(' </mesh_value_collection>\n')
ofile.write(' </mesh_function>\n')
ofile.write('</dolfin>\n')
[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")