from __future__ import print_function
import numpy as nm
from sfepy.base.base import output, assert_, OneTypeList, Struct
from sfepy.base.timing import Timer
from sfepy.discrete.common.region import (Region, get_dependency_graph,
sort_by_dependency, get_parents)
from sfepy.discrete.parse_regions import create_bnf, visit_stack, ParseException
[docs]def region_leaf(domain, regions, rdef, functions, tdim):
"""
Create/setup a region instance according to rdef.
"""
n_coor = domain.shape.n_nod
dim = domain.shape.dim
cmesh = domain.cmesh_tdim[tdim]
def _region_leaf(level, op):
token, details = op['token'], op['orig']
if token != 'KW_Region':
parse_def = token + '<' + ' '.join(details) + '>'
if token != 'E_COG':
region = Region('leaf', rdef, domain, parse_def=parse_def,
tdim=tdim)
if token == 'KW_Region':
details = details[1][2:]
aux = regions.find(details)
if not aux:
raise ValueError('region %s does not exist' % details)
else:
if rdef[:4] == 'copy':
region = aux.copy()
else:
region = aux
elif token == 'KW_All':
region.vertices = nm.arange(n_coor, dtype=nm.uint32)
elif token == 'E_VIR':
where = details[2]
if where[0] == '[':
vertices = nm.array(eval(where), dtype=nm.uint32)
assert_(nm.amin(vertices) >= 0)
assert_(nm.amax(vertices) < n_coor)
else:
coors = cmesh.coors
y = z = None
x = coors[:,0]
if dim > 1:
y = coors[:,1]
if dim > 2:
z = coors[:,2]
coor_dict = {'x' : x, 'y' : y, 'z': z}
vertices = nm.where(eval(where, {}, coor_dict))[0]
region.vertices = vertices
elif token == 'E_VOS':
facets = cmesh.get_surface_facets()
region.set_kind('facet')
region.facets = facets
elif token == 'E_VBF':
where = details[2]
coors = cmesh.coors
fun = functions[where]
vertices = fun(coors, domain=domain)
region.vertices = vertices
elif token == 'E_CBF':
where = details[2]
coors = domain.get_centroids(dim)
fun = functions[where]
cells = fun(coors, domain=domain)
region.cells = cells
elif token == 'E_COG':
group = int(details[3])
td = 0
for k in range(4):
if domain.cmesh_tdim[k] is not None:
cg = domain.cmesh_tdim[k].cell_groups
if nm.any(cg == group):
td = k
break
if td == 0:
raise ValueError('cell group %s does not exist' % group)
region = Region('leaf', rdef, domain, parse_def=parse_def, tdim=td)
region.cells = \
nm.where(domain.cmesh_tdim[td].cell_groups == group)[0]
elif token == 'E_COSET':
raise NotImplementedError('element sets not implemented!')
elif token == 'E_VOG':
group = int(details[3])
region.vertices = nm.where(cmesh.vertex_groups == group)[0]
elif token == 'E_VOSET':
try:
vertices = domain.vertex_set_bcs[details[3]]
except KeyError:
msg = 'undefined vertex set! (%s)' % details[3]
raise ValueError(msg)
region.vertices = vertices
elif token == 'E_OVIR':
aux = regions[details[3][2:]]
region.vertices = aux.vertices[0:1]
elif token == 'E_VI':
region.vertices = nm.array([int(ii) for ii in details[1:]],
dtype=nm.uint32)
elif token == 'E_CI':
region.cells = nm.array([int(ii) for ii in details[1:]],
dtype=nm.uint32)
else:
output('token "%s" unkown - check regions!' % token)
raise NotImplementedError
return region
return _region_leaf
[docs]def region_op(level, op_code, item1, item2):
token = op_code['token']
op = {'S' : '-', 'A' : '+', 'I' : '*'}[token[3]]
if token[-1] == 'V':
return item1.eval_op_vertices(item2, op)
elif token[-1] == 'E':
return item1.eval_op_edges(item2, op)
elif token[-1] == 'F':
return item1.eval_op_faces(item2, op)
elif token[-1] == 'S':
return item1.eval_op_facets(item2, op)
elif token[-1] == 'C':
return item1.eval_op_cells(item2, op)
else:
raise ValueError('unknown region operator token! (%s)' % token)
[docs]class Domain(Struct):
def __init__(self, name, mesh=None, nurbs=None, bmesh=None, regions=None,
verbose=False):
Struct.__init__(self, name=name, mesh=mesh, nurbs=nurbs, bmesh=bmesh,
regions=regions, verbose=verbose)
[docs] def get_centroids(self, dim):
"""
Return the coordinates of centroids of mesh entities with dimension
`dim`.
"""
return self.cmesh.get_centroids(dim)
[docs] def has_faces(self):
return self.shape.tdim == 3
[docs] def reset_regions(self):
"""
Reset the list of regions associated with the domain.
"""
self.regions = OneTypeList(Region)
self._region_stack = []
self._bnf = create_bnf(self._region_stack)
[docs] def create_extra_tdim_region(self, region, functions, tdim):
from sfepy.discrete.fem.geometry_element import (GeometryElement,
create_geometry_elements)
from sfepy.discrete import PolySpace
"""
Create a new region which has its own cmesh with
topological dimension tdim.
"""
mesh = self.mesh
if mesh.cmesh_tdim[tdim] is not None:
raise ValueError(f'cmesh of dimension {tdim} already exists!')
aux = mesh.from_region(region, mesh, tdim=tdim)
cmesh = aux.cmesh
new_mat_id = nm.max([nm.max(k.cell_groups) for k in mesh.cmesh_tdim
if k is not None]) + 1
cmesh.cell_groups[:] = new_mat_id
mesh.cmesh_tdim[tdim] = cmesh
mesh.descs += aux.descs
mesh.dims += aux.dims
mesh.n_el += aux.n_el
cmesh.set_local_entities(create_geometry_elements())
cmesh.setup_entities()
gel = GeometryElement(aux.descs[0])
if gel.dim > 0:
gel.create_surface_facet()
new_gel_entry = {aux.descs[0]: gel}
self.geom_els.update(new_gel_entry)
self.fix_element_orientation(geom_els=new_gel_entry, force_check=True)
key = gel.get_interpolation_name()
gel.poly_space = PolySpace.any_from_args(key, gel, 1)
gel = gel.surface_facet
if gel is not None:
key = gel.get_interpolation_name()
gel.poly_space = PolySpace.any_from_args(key, gel, 1)
select = f'cells of group {new_mat_id}'
self._bnf.parseString(select)
region = visit_stack(self._region_stack, region_op,
region_leaf(self, self.regions, select,
functions, tdim))
region.field_dim = tdim
return region
[docs] def create_region(self, name, select, kind='cell', parent=None,
check_parents=True, extra_options=None, functions=None,
add_to_regions=True, allow_empty=False):
"""
Region factory constructor. Append the new region to
self.regions list.
"""
if check_parents:
parents = get_parents(select)
for p in parents:
if p not in [region.name for region in self.regions]:
msg = 'parent region %s of %s not found!' % (p, name)
raise ValueError(msg)
stack = self._region_stack
try:
self._bnf.parseString(select)
except ParseException:
print('parsing failed:', select)
raise
tdim = self.shape.tdim if parent is None else self.regions[parent].tdim
region = visit_stack(stack, region_op,
region_leaf(self, self.regions, select,
functions, tdim))
finalize = True
eopts = extra_options
if eopts is not None:
if 'mesh_dim' in eopts:
region = self.create_extra_tdim_region(region, functions,
eopts['mesh_dim'])
if not(eopts.get('finalize', True)):
finalize = False
if 'vertices_from' in eopts:
vreg = eopts['vertices_from']
region.entities[0] = self.regions[vreg].vertices.copy()
finalize = False
region.name = name
region.definition = select
region.set_kind(kind)
if finalize:
region.finalize(allow_empty=allow_empty)
region.parent = parent
region.extra_options = extra_options
region.update_shape()
if add_to_regions:
self.regions.append(region)
return region
[docs] def create_regions(self, region_defs, functions=None, allow_empty=False):
output('creating regions...')
timer = Timer(start=True)
self.reset_regions()
##
# Sort region definitions by dependencies.
graph, name_to_sort_name = get_dependency_graph(region_defs)
sorted_regions = sort_by_dependency(graph)
##
# Define regions.
for name in sorted_regions:
sort_name = name_to_sort_name[name]
rdef = region_defs[sort_name]
region = self.create_region(name, rdef.select,
kind=rdef.get('kind', 'cell'),
parent=rdef.get('parent', None),
check_parents=False,
extra_options=rdef.get('extra_options',
None),
functions=functions,
allow_empty=allow_empty)
output(' ', region.name)
output('...done in %.2f s' % timer.stop())
return self.regions
[docs] def save_regions(self, filename, region_names=None):
"""
Save regions as individual meshes.
Parameters
----------
filename : str
The output filename.
region_names : list, optional
If given, only the listed regions are saved.
"""
import os
if region_names is None:
region_names = self.regions.get_names()
trunk, suffix = os.path.splitext(filename)
output('saving regions...')
for name in region_names:
region = self.regions[name]
output(name)
dim = region.tdim
is_surf = not region.can[dim] and region.can[dim - 1]
aux = self.mesh.from_region(region, self.mesh, is_surface=is_surf)
aux.write('%s_%s%s' % (trunk, region.name, suffix),
io='auto')
output('...done')
[docs] def save_regions_as_groups(self, filename, region_names=None):
"""
Save regions in a single mesh but mark them by using different
element/node group numbers.
If regions overlap, the result is undetermined, with exception of the
whole domain region, which is marked by group id 0.
Region masks are also saved as scalar point data for output formats
that support this.
Parameters
----------
filename : str
The output filename.
region_names : list, optional
If given, only the listed regions are saved.
"""
output('saving regions as groups...')
aux = self.mesh.copy()
n_ig = c_ig = 0
n_nod = self.shape.n_nod
# The whole domain region should go first.
names = (region_names if region_names is not None
else self.regions.get_names())
for name in names:
region = self.regions[name]
if region.vertices.shape[0] == n_nod:
names.remove(region.name)
names = [region.name] + names
break
out = {}
for name in names:
region = self.regions[name]
output(region.name)
aux.cmesh.vertex_groups[region.vertices] = n_ig
n_ig += 1
mask = nm.zeros((n_nod, 1), dtype=nm.float64)
mask[region.vertices] = 1.0
out[name] = Struct(name='region', mode='vertex', data=mask,
var_name=name, dofs=None)
if region.has_cells():
ii = region.get_cells()
aux.cmesh.cell_groups[ii] = c_ig
c_ig += 1
aux.write(filename, io='auto', out=out)
output('...done')