# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import dolfin as df
import os
import numpy as np
from custEM.misc import logger_print as lp
import xml.etree.ElementTree
[docs]
def build_line_mesh(IB, axis, start, stop, n_segs, line_name, on_topo):
"""
Generate a 1D line-mesh for interpolation.
Required arguments
------------------
see *create_line_mesh* docs
"""
if line_name is None:
if axis == 'x':
out_name = 'x0_' + str(start) + '_x1_' + str(stop) +\
'_y_' + str(IB.y) + '_z_' + str(IB.z) + \
'_n_' + str(n_segs)
if axis == 'y':
out_name = 'y0_' + str(start) + '_y1_' + str(stop) +\
'_x_' + str(IB.x) + '_z_' + str(IB.z) + \
'_n_' + str(n_segs)
if axis == 'z':
out_name = 'z0_' + str(start) + '_z1_' + str(stop) +\
'_x_' + str(IB.x) + '_y_' + str(IB.y) + \
'_n_' + str(n_segs)
if on_topo:
out_name += '_on_topo'
else:
out_name = line_name
out_name = IB.MP.m_dir + '/lines/' + out_name + '_line_' + axis + '.xml'
if os.path.isfile(out_name) and not IB.force_create:
lp(IB.MP.logger, 10,
'... re-using existing "' + axis + '_line_mesh" ...',
barrier=False, pre_dash=False)
else:
lp(IB.MP.logger, 20,
'... building ' + axis + '-directed line_mesh ...',
barrier=False, pre_dash=False)
# generate coordinates
points = np.zeros((n_segs + 1, 3))
if axis == 'x':
points[:, 0] = np.linspace(start, stop, n_segs + 1)
points[:, 1] = IB.y
points[:, 2] = IB.z
if axis == 'y':
points[:, 0] = IB.x
points[:, 1] = np.linspace(start, stop, n_segs + 1)
points[:, 2] = IB.z
if axis == 'z':
points[:, 0] = IB.x
points[:, 1] = IB.y
points[:, 2] = np.linspace(start, stop, n_segs + 1)
# apply topography to horizontal lines if specified
if (axis == 'x' or axis == 'y') and on_topo:
points = np.array(points)
if os.path.isfile(IB.t_dir + '/' +
IB.MP.mesh_name + '_surf.xyz'):
from custEM.meshgen.dem_interpolator import DEM
dem = DEM(IB.t_dir + '/' + IB.MP.mesh_name +
'_surf.xyz', centering=False,
easting_shift=None, northing_shift=None)
points[:, 2] = dem(points[:, 0], points[:, 1], rotation=None)
else:
import custEM.misc.synthetic_definitions as sd
topo_func = getattr(sd, str(IB.MP.topo))
points[:, 2] = topo_func(points[:, 0], points[:, 1])
if IB.on_water_surface:
points[:, 2][points[:, 2] < 0.] = 0.
if IB.z != 0.:
points[:, 2] += IB.z
# add vertices and pseudo-elements to an empty mesh instance
mesh = df.Mesh(IB.MP.mpi_cs)
editor = df.MeshEditor()
editor.open(mesh, 'tetrahedron', 3, 3)
editor.init_vertices(n_segs + 1)
editor.init_cells(n_segs)
for pi, point in enumerate(points):
editor.add_vertex(pi, point)
if pi != 0:
editor.add_cell(pi - 1, [pi - 1, pi])
editor.close()
df.File(IB.MP.mpi_cs, out_name) << mesh
# hack some values in the xml file for a conformal view with Paraview
tree = xml.etree.ElementTree.parse(out_name)
for elem in tree.iter():
if elem.tag == 'mesh':
elem.attrib['celltype'] = 'interval'
if elem.tag == 'tetrahedron':
elem.tag = 'interval'
elem.attrib['v0'] = elem.attrib['v2']
elem.attrib['v1'] = elem.attrib['v3']
elem.attrib.pop('v2')
elem.attrib.pop('v3')
tree.write(out_name)
[docs]
def build_slice_mesh(IB, axis, dim, n_segs, slice_name, on_topo):
"""
Generate a 2D slice-mesh for interpolation.
Required arguments
------------------
see *create_slice_mesh* docs
"""
# reshape dim and n_segs if required
if type(dim) is not list:
dim = [dim, dim]
if type(n_segs) is not list:
n_segs = [n_segs, n_segs]
if slice_name is None:
out_name = 'dim_' + str(dim) +\
'_x_' + str(IB.x) + '_y_' + str(IB.y) +\
'_z_' + str(IB.z) + '_n_' + str(n_segs)
if on_topo:
out_name += '_on_topo'
else:
out_name = slice_name
out_name = IB.MP.m_dir + '/slices/' + out_name + '_slice_' + axis + '.xml'
if os.path.isfile(out_name) and not IB.force_create:
lp(IB.MP.logger, 10,
'... re-using existing "' + axis + '_slice_mesh" ...',
barrier=False, pre_dash=False)
else:
lp(IB.MP.logger, 20,
'... building ' + axis + '-normal slice_mesh ...',
barrier=False, pre_dash=False)
# generate coordinates
points = np.zeros(((n_segs[0] + 1) * (n_segs[1] + 1), 3))
vec_i = np.linspace(-dim[0]/2., dim[0]/2., n_segs[0] + 1)
vec_j = np.linspace(-dim[1]/2., dim[1]/2., n_segs[1] + 1)
counter = 0
for vi in vec_i:
for vj in vec_j:
if axis == 'x':
points[counter, 0] = IB.x
points[counter, 1] = vi
points[counter, 2] = vj
if axis == 'y':
points[counter, 0] = vi
points[counter, 1] = IB.y
points[counter, 2] = vj
if axis == 'z':
points[counter, 0] = vi
points[counter, 1] = vj
points[counter, 2] = IB.z
counter += 1
# apply topography to horizontal lines if specified
if axis == 'z' and on_topo:
points = np.array(points)
if os.path.isfile(IB.t_dir + '/' +
IB.MP.mesh_name + '_surf.xyz'):
from custEM.meshgen.dem_interpolator import DEM
dem = DEM(IB.t_dir + '/' + IB.MP.mesh_name +
'_surf.xyz', centering=False,
easting_shift=None, northing_shift=None)
points[:, 2] = dem(points[:, 0], points[:, 1], rotation=None)
else:
import custEM.misc.synthetic_definitions as sd
topo_func = getattr(sd, str(IB.MP.topo))
points[:, 2] = topo_func(points[:, 0], points[:, 1])
if IB.on_water_surface:
points[:, 2][points[:, 2] < 0.] = 0.
if IB.z != 0.:
points[:, 2] += IB.z
# shift to other origin than [0., 0., 0.] if specified
if axis == 'x':
points[:, 1] += IB.y
points[:, 2] += IB.z
if axis == 'y':
points[:, 0] += IB.x
points[:, 2] += IB.z
if axis == 'z':
points[:, 0] += IB.x
points[:, 1] += IB.y
# add vertices and pseudo-elements to an empty mesh instance
mesh = df.Mesh(IB.MP.mpi_cs)
editor = df.MeshEditor()
editor.open(mesh, 'tetrahedron', 3, 3)
editor.init_vertices((n_segs[0] + 1) * (n_segs[1] + 1))
editor.init_cells(n_segs[0] * n_segs[1] * 2)
counter = 0
cell_counter = 0
for i in range(len(vec_i)):
for j in range(len(vec_j)):
editor.add_vertex(counter, points[counter])
if i != 0 and j != 0:
editor.add_cell(cell_counter,
[counter - n_segs[1] - 2, counter, counter -1])
editor.add_cell(cell_counter + 1,
[counter - n_segs[1] - 2,
counter - n_segs[1] - 1, counter])
cell_counter += 2
counter += 1
editor.close()
df.File(IB.MP.mpi_cs, out_name) << mesh
# hack some values in the xml file for a conformal view with Paraview
tree = xml.etree.ElementTree.parse(out_name)
for elem in tree.iter():
if elem.tag == 'mesh':
elem.attrib['celltype'] = 'triangle'
if elem.tag == 'tetrahedron':
elem.tag = 'triangle'
elem.attrib['v0'] = elem.attrib['v1']
elem.attrib['v1'] = elem.attrib['v2']
elem.attrib['v2'] = elem.attrib['v3']
elem.attrib.pop('v3')
tree.write(out_name)
[docs]
def build_path_mesh(IB, points, path_name, on_topo, suffix):
"""
Generate an interpolation mesh for interpolation with an arbitrary list
of coordinates.
Required arguments
------------------
see *create_path_mesh* docs
"""
if path_name is None:
out_name = 'NOT_DEFINED'
lp(IB.MP.logger, 30,
'Warning! No export name for path interpolation mesh set.\n'
'Using *NOT_DEFINED* as dummy name. Continuing ...',
barrier=False, pre_dash=False, root_only=False)
else:
out_name = path_name
if suffix is None:
out_name = IB.MP.m_dir + '/paths/' + out_name + '_path.xml'
else:
if 'line' in suffix:
directory = '/lines/'
elif 'slice' in suffix:
directory = '/slices/'
else:
directory = '/paths/'
out_name = IB.MP.m_dir + directory + out_name + '_' + suffix + '.xml'
if os.path.isfile(out_name) and not IB.force_create:
lp(IB.MP.logger, 10,
'... re-using existing path mesh ...',
barrier=False, pre_dash=False, root_only=False)
else:
lp(IB.MP.logger, 20,
'... building path mesh from specified points ...',
barrier=False, pre_dash=False, root_only=False)
# apply topography to horizontal lines if specified
if on_topo:
points = np.array(points)
if os.path.isfile(IB.t_dir + '/' +
IB.MP.mesh_name + '_surf.xyz'):
from custEM.meshgen.dem_interpolator import DEM
dem = DEM(IB.t_dir + '/' + IB.MP.mesh_name +
'_surf.xyz', centering=False,
easting_shift=None, northing_shift=None)
points[:, 2] = dem(points[:, 0], points[:, 1], rotation=None)
else:
import custEM.misc.synthetic_definitions as sd
topo_func = getattr(sd, str(IB.MP.topo))
points[:, 2] = topo_func(points[:, 0], points[:, 1])
if IB.on_water_surface:
points[:, 2][points[:, 2] < 0.] = 0.
if IB.z != 0.:
points[:, 2] += IB.z
# add vertices and pseudo-elements to an empty mesh instance
mesh = df.Mesh(IB.MP.mpi_cs)
editor = df.MeshEditor()
editor.open(mesh, 'tetrahedron', 3, 3)
editor.init_vertices(len(points))
editor.init_cells(len(points) - 1)
for pi, point in enumerate(points):
editor.add_vertex(pi, point)
if pi != 0:
editor.add_cell(pi - 1, [pi - 1, pi])
editor.close()
df.File(IB.MP.mpi_cs, out_name) << mesh
# hack some values in the xml file for a conformal view with Paraview
tree = xml.etree.ElementTree.parse(out_name)
for elem in tree.iter():
if elem.tag == 'mesh':
elem.attrib['celltype'] = 'interval'
if elem.tag == 'tetrahedron':
elem.tag = 'interval'
elem.attrib['v0'] = elem.attrib['v2']
elem.attrib['v1'] = elem.attrib['v3']
elem.attrib.pop('v2')
elem.attrib.pop('v3')
tree.write(out_name)