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.fem.poly_spaces 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 = nm.empty_like(mask) 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() mask[econn] += 1 # 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. if not nm.all(mask[nodes] > 0): 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
def _get_edge_path(graph, seed, mask, cycle=False): """ Get a path in an edge graph starting with seed. The mask is incremented by one at positions of the path vertices. """ if mask[seed]: return [] path = [seed] mask[seed] = 1 row = graph[seed].indices nv = len(row) while nv: if nv == 2: if mask[row[0]]: if mask[row[1]]: if cycle: path.append(seed) break else: vert = row[1] else: vert = row[0] elif mask[row[0]]: break else: vert = row[0] path.append(vert) mask[vert] = 1 row = graph[vert].indices nv = len(row) path = nm.array(path, dtype=nm.int32) return path
[docs]def get_edge_paths(graph, mask): """ 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: path = _get_edge_path(graph, seed, mask) if len(path): paths.append(path) # 2. get possible remaing cycles. while 1: ii = nm.where(mask[nodes] == 0)[0] 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) mask = nm.zeros_like(imap) try: paths = get_edge_paths(graph, mask) except ValueError: raise ValueError('more than 2 edges sharing a vertex in region %s!' % region.name) # All nodes must have an edge direction. if not nm.all(mask[nodes]): 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