# -*- 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 mesh2xml2h5(mesh_name, m_dir, rot):
"""
Converts TetGen meshes in Medit (*.mesh*) 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 .mesh and .xml, parser implemented as a
state machine:
0 = read 'Dimension'
1 = read dimension
2 = read 'Vertices'
3 = read number of vertices
4 = read next vertex
5 = read 'Triangles' or 'Tetrahedra'
6 = read number of cells
7 = read next cell
8 = read next domains
9 = done
"""
logger = logging.getLogger('custEM')
logger.info('... converting from Medit format (.mesh) to *xml* and '
'*h5* format ...')
ifile = open(m_dir + '/_mesh/' + mesh_name + '.mesh', "r")
ofile = open(m_dir + '/_xml/' + mesh_name + '.xml', "w")
cell_type = "tetrahedron"
dim = 3
# Step to beginning of file
ifile.seek(0)
# Write header
write_header_mesh(ofile, cell_type, dim)
# Current state
state = 0
# Write data
num_vertices_read = 0
num_cells_read = 0
marker = []
while 1:
# Read next line
line = ifile.readline()
try:
line[1]
except IndexError:
continue
if not line:
break
# Skip comments
if line[0] == '#':
continue
# Remove newline
if line[-1] == "\n":
line = line[:-1]
if state == 0:
if line[0] == 'D' and line[1] == 'i':
state += 1
elif state == 1:
num_dims = int(line)
state += 1
elif state == 2:
if line[0] == 'V' and line[1] == 'e':
state += 1
elif state == 3:
num_vertices = int(line)
write_header_vertices(ofile, num_vertices)
state += 1
elif state == 4:
if num_dims == 1:
(x, tmp1, tmp2) = line.split()
x = float(x)
y = 0.0
z = 0.0
elif num_dims == 2:
(x, y, tmp) = line.split()
x = float(x)
y = float(y)
z = 0.0
elif num_dims == 3:
(x, y, z, tmp) = line.split()
x = float(x)
y = float(y)
z = float(z)
if rot is not None:
# rotate pre-rotated mesh back to original rotation if
# post_rotate flag during mesh generation is not None
xyzrot = mu.rotate_around_point([[x, y, z]], [0., 0., 0.],
[np.deg2rad(rot)], ['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 == 5:
if line == "Intervals" and num_dims == 1:
state += 1
if line == "Triangles" and num_dims == 2:
state += 1
if line[0] == 'T' and line[1] == 'e' and num_dims == 3:
state += 1
elif state == 6:
num_cells = int(line)
num_domains = num_cells
write_header_cells(ofile, num_cells)
state += 1
elif state == 7:
if num_dims == 1:
(n0, n1, tmp1, tmp2) = line.split()
n0 = int(n0) - 1
n1 = int(n1) - 1
write_cell_interval(ofile, num_cells_read, n0, n1)
elif num_dims == 2:
(n0, n1, n2, tmp) = line.split()
n0 = int(n0) - 1
n1 = int(n1) - 1
n2 = int(n2) - 1
write_cell_triangle(ofile, num_cells_read, n0, n1, n2)
elif num_dims == 3:
(n0, n1, n2, n3, tmp) = line.split()
n0 = int(n0) - 1
n1 = int(n1) - 1
n2 = int(n2) - 1
n3 = int(n3) - 1
marker.append(tmp)
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 == 8:
write_header_domains(ofile, num_domains)
state += 1
elif state == 9:
for j in range(num_domains):
write_domain_marker(ofile, j, marker[j])
write_footer_domains(ofile)
state += 1
elif state == 10:
break
# Check that we got all data
if state == 10:
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_domains, marker)
# 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_graph_vertex(ofile, vertex, num_edges, weight=1):
"""
Write vertex graph to file.
"""
ofile.write(' <vertex index=\"%d\" num_edges=\"%d\" '
'weight=\"%d\"/>\n' % (vertex, num_edges, weight))
[docs]
def write_graph_edge(ofile, v1, v2, weight=1):
"""
Write edge graph to file.
"""
ofile.write(' <edge v1=\"%d\" v2=\"%d\" '
'weight=\"%d\"/>\n' % (v1, v2, weight))
[docs]
def write_cell_interval(ofile, cell, n0, n1):
"""
Write interval (1D) information to domain file.
"""
ofile.write(' <interval index=\"%d\" v0=\"%d\" '
'v1=\"%d\"/>\n' % (cell, n0, n1))
[docs]
def write_cell_triangle(ofile, cell, n0, n1, n2):
"""
Write triangle 2D information to domain file.
"""
ofile.write(' <triangle index=\"%d\" v0=\"%d\" v1=\"%d\" '
'v2=\"%d\"/>\n' % (cell, n0, n1, n2))
[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, marker):
"""
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, marker[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")