# Source code for sfepy.discrete.fem.utils

```from __future__ import absolute_import
import numpy as nm

import sfepy.linalg as la
from sfepy.discrete.integrals import Integral
from sfepy.discrete import PolySpace
from six.moves import range

[docs]def prepare_remap(indices, n_full):
"""
Prepare vector for remapping range `[0, n_full]` to its subset given
by `indices`.
"""
remap = nm.empty((n_full,), dtype=nm.int32)
remap.fill(-1)
remap[indices] = nm.arange(indices.shape[0], dtype=nm.int32)

return remap

[docs]def invert_remap(remap):
"""
Return the inverse of `remap`, i.e. a mapping from a sub-range
indices to a full range, see :func:`prepare_remap()`.
"""
if remap is not None:
inverse = nm.where(remap >= 0)[0].astype(nm.int32)

else:
inverse = None

return inverse

[docs]def prepare_translate(old_indices, new_indices):
"""
Prepare vector for translating `old_indices` to `new_indices`.

Returns
-------
translate : array
The translation vector. Then `new_ar = translate[old_ar]`.
"""
old_indices = nm.asarray(old_indices)
new_indices = nm.asarray(new_indices)

translate = nm.zeros(old_indices.max() + 1, dtype=new_indices.dtype)
translate[old_indices] = new_indices

return translate

[docs]def compute_nodal_normals(nodes, region, field, return_imap=False):
"""
Nodal normals are computed by simple averaging of element normals of
elements every node is contained in.
"""
dim = region.dim

field.domain.create_surface_group(region)
field.setup_surface_data(region)

# Custom integral with quadrature points in nodes.
ps = PolySpace.any_from_args('', field.gel.surface_facet,
field.approx_order)
qp_coors = ps.node_coors
# Unit normals -> weights = ones.
qp_weights = nm.ones(qp_coors.shape[0], dtype=nm.float64)

integral = Integral('aux', coors=qp_coors, weights=qp_weights)

normals = nm.zeros((nodes.shape[0], dim), dtype=nm.float64)
mask = nm.zeros((nodes.max() + 1,), dtype=nm.int32)
imap.fill(nodes.shape[0]) # out-of-range index for normals.
imap[nodes] = nm.arange(nodes.shape[0], dtype=nm.int32)

cmap, _ = field.get_mapping(region, integral, 'surface')
e_normals = cmap.normal[..., 0]

sd = field.surface_data[region.name]
econn = sd.get_connectivity()

# normals[imap[econn]] += e_normals
im = imap[econn]
for ii, en in enumerate(e_normals):
normals[im[ii]] += en

# All nodes must have a normal.
raise ValueError('region %s has not complete faces!' % region.name)

norm = la.norm_l2_along_axis(normals)[:, nm.newaxis]
if (norm < 1e-15).any():
raise ValueError('zero nodal normal! (a node in volume?)')
normals /= norm

if return_imap:
return normals, imap

else:
return normals

"""
Get a path in an edge graph starting with seed. The mask is incremented by
one at positions of the path vertices.
"""
return []

path = [seed]

row = graph[seed].indices
nv = len(row)
while nv:
if nv == 2:
if cycle:
path.append(seed)
break

else:
vert = row[1]

else:
vert = row[0]

break

else:
vert = row[0]

path.append(vert)

row = graph[vert].indices
nv = len(row)

path = nm.array(path, dtype=nm.int32)

return path

"""
Get all edge paths in a graph with non-masked vertices. The mask is
updated.
"""
nodes = nm.unique(graph.indices)
npv = nm.diff(graph.indptr)

if npv.max() > 2:
raise ValueError('more than 2 edges sharing a vertex!')

seeds = nm.where(npv == 1)[0]

# 1. get paths.
paths = []
for seed in seeds:
if len(path):
paths.append(path)

# 2. get possible remaing cycles.
while 1:
if not len(ii):
break

path = _get_edge_path(graph, nodes[ii[0]], mask, cycle=True)
if len(path):
paths.append(path)

return paths

[docs]def compute_nodal_edge_dirs(nodes, region, field, return_imap=False):
"""
Nodal edge directions are computed by simple averaging of direction vectors
of edges a node is contained in. Edges are assumed to be straight and a
node must be on a single edge (a border node) or shared by exactly two
edges.
"""
coors = region.domain.mesh.coors
dim = coors.shape[1]

graph = region.get_edge_graph()

imap = prepare_remap(nodes, nodes.max() + 1)

try:

except ValueError:
raise ValueError('more than 2 edges sharing a vertex in region %s!'
% region.name)

# All nodes must have an edge direction.
raise ValueError('region %s has not complete edges!' % region.name)

edge_dirs = nm.zeros((nodes.shape[0], dim), dtype=nm.float64)
for path in paths:
pcoors = coors[path]

edirs = nm.diff(pcoors, axis=0)
la.normalize_vectors(edirs, eps=1e-12)

im = imap[nm.c_[path[:-1], path[1:]]]
for ii, edir in enumerate(edirs):
edge_dirs[im[ii]] += edir

la.normalize_vectors(edge_dirs, eps=1e-12)

if return_imap:
return edge_dirs, imap

else:
return edge_dirs

[docs]def get_min_value(dofs):
"""
Get a reasonable minimal value of DOFs suitable for extending over a
whole domain.
"""
if dofs.shape[1] > 1: # Vector.
val = 0.0

else: # Scalar.
val = dofs.min()

return val

[docs]def extend_cell_data(data, domain, rname, val=None, is_surface=False,
average_surface=True):
"""
Extend cell data defined in a region to the whole domain.

Parameters
----------
data : array
The data defined in the region.
domain : FEDomain instance
The FE domain.
rname : str
The region name.
val : float, optional
The value for filling cells not covered by the region. If not given,
the smallest value in data is used.
is_surface : bool
If True, the data are defined on a surface region. In that case the
values are averaged or summed into the cells containing the region
surface faces (a cell can have several faces of the surface), see
`average_surface`.
average_surface : bool
If True, the data defined on a surface region are averaged, otherwise
the data are summed.

Returns
-------
edata : array
The data extended to all domain elements.
"""
n_el = domain.shape.n_el
if data.shape[0] == n_el: return data

if val is None:
if data.shape[2] > 1: # Vector.
val = nm.amin(nm.abs(data))
else: # Scalar.
val = nm.amin(data)

edata = nm.empty((n_el,) + data.shape[1:], dtype=data.dtype)
edata.fill(val)

region = domain.regions[rname]

if not is_surface:
edata[region.get_cells()] = data

else:
cells = region.get_cells(true_cells_only=False)
ucells = nm.unique(cells)

if len(cells) != len(region.facets):
raise ValueError('region %s has an inner face!'
% region.name)

if average_surface:
avg = nm.bincount(cells, minlength=n_el)[ucells]

else:
avg = 1.0

for ic in range(data.shape[2]):
if nm.isrealobj(data):
evals = nm.bincount(cells, weights=data[:, 0, ic, 0],
minlength=n_el)[ucells]

else:
evals = (nm.bincount(cells, weights=data[:, 0, ic, 0].real,
minlength=n_el)[ucells]
+ 1j *
nm.bincount(cells, weights=data[:, 0, ic, 0].imag,
minlength=n_el)[ucells])

edata[ucells, 0, ic, 0] = evals / avg

return edata

[docs]def refine_mesh(filename, level):
"""
Uniformly refine `level`-times a mesh given by `filename`.

The refined mesh is saved to a file with name constructed from base
name of `filename` and `level`-times appended `'_r'` suffix.

Parameters
----------
filename : str
The mesh file name.
level : int
The refinement level.
"""
import os
from sfepy.base.base import output
from sfepy.discrete.fem import Mesh, FEDomain

if level > 0:
mesh = Mesh.from_file(filename)
domain = FEDomain(mesh.name, mesh)
for ii in range(level):
output('refine %d...' % ii)
domain = domain.refine()
output('... %d nodes %d elements'
% (domain.shape.n_nod, domain.shape.n_el))

suffix = os.path.splitext(filename)[1]
filename = domain.name + suffix

domain.mesh.write(filename, io='auto')

return filename
```