"""
Computational domain, consisting of the mesh and regions.
"""
import time
import numpy as nm
from sfepy.base.base import output, Struct
from geometry_element import GeometryElement
from sfepy.discrete.common.domain import Domain
from sfepy.discrete.fem.refine import refine_2_3, refine_2_4, refine_3_4, refine_3_8
from sfepy.discrete.fem.fe_surface import FESurface
from sfepy.discrete.fem.mesh import make_inverse_connectivity
import fea
[docs]class FEDomain(Domain):
"""
Domain is divided into groups, whose purpose is to have homogeneous
data shapes.
"""
def __init__(self, name, mesh, verbose=False, **kwargs):
"""Create a Domain.
Parameters
----------
name : str
Object name.
mesh : Mesh
A mesh defining the domain.
"""
Domain.__init__(self, name, mesh=mesh, verbose=verbose, **kwargs)
self.geom_els = geom_els = {}
for ig, desc in enumerate(mesh.descs):
gel = GeometryElement(desc)
if gel.dim > 1:
# Create geometry elements of dimension - 1.
gel.create_surface_facet()
geom_els[desc] = gel
self.geom_interps = interps = {}
for gel in geom_els.itervalues():
key = gel.get_interpolation_name()
gel.interp = interps.setdefault(key, fea.Interpolant(key, gel))
gel = gel.surface_facet
if gel is not None:
key = gel.get_interpolation_name()
gel.interp = interps.setdefault(key, fea.Interpolant(key, gel))
self.vertex_set_bcs = self.mesh.nodal_bcs
self.mat_ids_to_i_gs = {}
for ig, mat_id in enumerate(mesh.mat_ids):
self.mat_ids_to_i_gs[mat_id[0]] = ig
n_nod, dim = self.mesh.coors.shape
self.shape = Struct(n_nod=n_nod, dim=dim, tdim=0,
n_el=0,
n_gr=len(self.mesh.conns))
self.setup_groups()
self.fix_element_orientation()
self.reset_regions()
self.clear_surface_groups()
from sfepy.discrete.fem.geometry_element import create_geometry_elements
from sfepy.discrete.fem.extmods.cmesh import CMesh
self.cmesh = CMesh.from_mesh(mesh)
gels = create_geometry_elements()
self.cmesh.set_local_entities(gels)
self.cmesh.setup_entities()
self.shape.tdim = self.cmesh.tdim
self.cell_offsets = self.mesh.el_offsets
[docs] def setup_groups(self):
self.groups = {}
for ii in range(self.shape.n_gr):
gel = self.geom_els[self.mesh.descs[ii]] # Shortcut.
conn = self.mesh.conns[ii]
vertices = nm.unique(conn)
n_vertex = vertices.shape[0]
n_el, n_ep = conn.shape
n_edge = gel.n_edge
n_edge_total = n_edge * n_el
if gel.dim == 3:
n_face = gel.n_face
n_face_total = n_face * n_el
else:
n_face = n_face_total = 0
shape = Struct(n_vertex=n_vertex, n_el=n_el, n_ep=n_ep,
n_edge=n_edge, n_edge_total=n_edge_total,
n_face=n_face, n_face_total=n_face_total,
dim=self.mesh.dims[ii])
self.groups[ii] = Struct(ig=ii, vertices=vertices, conn=conn,
gel=gel, shape=shape)
self.shape.n_el += n_el
[docs] def iter_groups(self, igs=None):
if igs is None:
for ig in xrange(self.shape.n_gr): # sorted by ig.
yield self.groups[ig]
else:
for ig in igs:
yield ig, self.groups[ig]
[docs] def get_cell_offsets(self):
offs = {}
for ig in range(self.shape.n_gr):
offs[ig] = self.cell_offsets[ig]
return offs
[docs] def get_mesh_coors(self, actual=False):
"""
Return the coordinates of the underlying mesh vertices.
"""
if actual and hasattr(self.mesh, 'coors_act'):
return self.mesh.coors_act
else:
return self.mesh.coors
[docs] def get_mesh_bounding_box(self):
"""
Return the bounding box of the underlying mesh.
Returns
-------
bbox : ndarray (2, dim)
The bounding box with min. values in the first row and max. values
in the second row.
"""
return self.mesh.get_bounding_box()
[docs] def get_diameter(self):
"""
Return the diameter of the domain.
Notes
-----
The diameter corresponds to the Friedrichs constant.
"""
bbox = self.get_mesh_bounding_box()
return (bbox[1,:] - bbox[0,:]).max()
[docs] def fix_element_orientation(self):
"""
Ensure element nodes ordering giving positive element volume.
The groups with elements of lower dimension than the space dimension
are skipped.
"""
from extmods.cmesh import orient_elements
coors = self.mesh.coors
for ii, group in self.groups.iteritems():
if group.shape.dim < self.shape.dim: continue
ori, conn = group.gel.orientation, group.conn
itry = 0
while itry < 2:
flag = -nm.ones(conn.shape[0], dtype=nm.int32)
# Changes orientation if it is wrong according to swap*!
# Changes are indicated by positive flag.
orient_elements(flag, conn, coors,
ori.roots, ori.vecs,
ori.swap_from, ori.swap_to)
if nm.alltrue(flag == 0):
if itry > 0: output('...corrected')
itry = -1
break
output('warning: bad element orientation, trying to correct...')
itry += 1
if itry == 2 and flag[0] != -1:
raise RuntimeError('elements cannot be oriented! (%d, %s)'
% (ii, self.mesh.descs[ii]))
elif flag[0] == -1:
output('warning: element orienation not checked')
[docs] def get_element_diameters(self, ig, cells, vg, mode, square=True):
group = self.groups[ig]
diameters = nm.empty((len(cells), 1, 1, 1), dtype=nm.float64)
if vg is None:
diameters.fill(1.0)
else:
vg.get_element_diameters(diameters, group.gel.edges,
self.get_mesh_coors().copy(), group.conn,
cells.astype(nm.int32), mode)
if square:
out = diameters.squeeze()
else:
out = nm.sqrt(diameters.squeeze())
return out
[docs] def get_evaluate_cache(self, cache=None, share_geometry=False):
"""
Get the evaluate cache for :func:`Variable.evaluate_at()
<sfepy.discrete.variables.Variable.evaluate_at()>`.
Parameters
----------
cache : Struct instance, optional
Optionally, use the provided instance to store the cache data.
share_geometry : bool
Set to True to indicate that all the probes will work on the same
domain. Certain data are then computed only for the first probe and
cached.
Returns
-------
cache : Struct instance
The evaluate cache.
"""
try:
from scipy.spatial import cKDTree as KDTree
except ImportError:
from scipy.spatial import KDTree
if cache is None:
cache = Struct(name='evaluate_cache')
tt = time.clock()
if (cache.get('iconn', None) is None) or not share_geometry:
mesh = self.mesh
offsets, iconn = make_inverse_connectivity(mesh.conns, mesh.n_nod,
ret_offsets=True)
ii = nm.where(offsets[1:] == offsets[:-1])[0]
if len(ii):
raise ValueError('some vertices not in any element! (%s)' % ii)
cache.offsets = offsets
cache.iconn = iconn
output('iconn: %f s' % (time.clock()-tt))
tt = time.clock()
if (cache.get('kdtree', None) is None) or not share_geometry:
cache.kdtree = KDTree(mesh.coors)
output('kdtree: %f s' % (time.clock()-tt))
return cache
[docs] def clear_surface_groups(self):
"""
Remove surface group data.
"""
self.surface_groups = {}
[docs] def create_surface_group(self, region):
"""
Create a new surface group corresponding to `region` if it does
not exist yet.
Notes
-----
Surface groups define surface facet connectivity that is needed
for :class:`sfepy.discrete.fem.mappings.SurfaceMapping`.
"""
for ig in region.igs:
groups = self.surface_groups.setdefault(ig, {})
if region.name not in groups:
group = self.groups[ig]
gel_faces = group.gel.get_surface_entities()
name = 'surface_group_%s_%d' % (region.name, ig)
surface_group = FESurface(name, region, gel_faces,
group.conn, ig)
groups[region.name] = surface_group
[docs] def refine(self):
"""
Uniformly refine the domain mesh.
Returns
-------
domain : FEDomain instance
The new domain with the refined mesh.
Notes
-----
Works only for meshes with single element type! Does not
preserve node groups!
"""
names = set()
for group in self.groups.itervalues():
names.add(group.gel.name)
if len(names) != 1:
msg = 'refine() works only for meshes with single element type!'
raise NotImplementedError(msg)
el_type = names.pop()
if el_type == '2_3':
mesh = refine_2_3(self.mesh, self.cmesh)
elif el_type == '2_4':
mesh = refine_2_4(self.mesh, self.cmesh)
elif el_type == '3_4':
mesh = refine_3_4(self.mesh, self.cmesh)
elif el_type == '3_8':
mesh = refine_3_8(self.mesh, self.cmesh)
else:
msg = 'unsupported element type! (%s)' % el_type
raise NotImplementedError(msg)
domain = FEDomain(self.name + '_r', mesh)
return domain