Source code for sfepy.discrete.fem.mesh

from __future__ import absolute_import
import numpy as nm
import scipy.sparse as sp

from sfepy.base.base import Struct, invert_dict, get_default, output,\
     assert_, is_sequence
from sfepy.base.timing import Timer
from .meshio import MeshIO
from scipy.spatial import cKDTree

eps = 1e-9

[docs]def set_accuracy(eps): globals()['eps'] = eps
[docs]def find_map(x1, x2, allow_double=False, join=True): """ Find a mapping between common coordinates in x1 and x2, such that x1[cmap[:,0]] == x2[cmap[:,1]] """ kdtree = cKDTree(nm.vstack([x1, x2])) cmap = kdtree.query_pairs(eps, output_type='ndarray') dns1 = nm.where(cmap[:, 1] < x1.shape[0])[0] dns2 = nm.where(cmap[:, 0] >= x1.shape[0])[0] if (dns1.size + dns2.size): output('double node(s) in:') for dn in dns1: idxs = cmap[dn, :] output('x1: %d %d -> %s %s' % (idxs[0], idxs[1], x1[idxs[0], :], x1[idxs[1], :])) for dn in dns2: idxs = cmap[dn, :] output('x2: %d %d -> %s %s' % (idxs[0], idxs[1], x2[idxs[0], :], x2[idxs[1], :])) if not allow_double: raise ValueError('double node(s)! (see above)') cmap[:, 1] -= x1.shape[0] return cmap if join else (cmap[:, 0], cmap[:, 1])
[docs]def merge_mesh(x1, ngroups1, conn1, mat_ids1, x2, ngroups2, conn2, mat_ids2, cmap): """ Merge two meshes in common coordinates found in x1, x2. Notes ----- Assumes the same number and kind of element groups in both meshes! """ n1 = x1.shape[0] n2 = x2.shape[0] err = nm.sum(nm.sum(nm.abs(x1[cmap[:,0],:-1] - x2[cmap[:,1],:-1]))) if abs(err) > (10.0 * eps): raise ValueError('nonmatching meshes! (error: %e)' % err) mask = nm.ones((n2,), dtype=nm.int32) mask[cmap[:,1]] = 0 remap = nm.cumsum(mask) + n1 - 1 remap[cmap[:,1]] = cmap[:,0] i2 = nm.setdiff1d(nm.arange( n2, dtype=nm.int32), cmap[:,1]) xx = nm.r_[x1, x2[i2]] ngroups = nm.r_[ngroups1, ngroups2[i2]] conn = nm.vstack((conn1, remap[conn2])) mat_ids = None if (mat_ids1 is not None) and (mat_ids2 is not None): mat_ids = nm.concatenate((mat_ids1, mat_ids2)) return xx, ngroups, conn, mat_ids
[docs]def fix_double_nodes(coor, ngroups, conns): """ Detect and attempt fixing double nodes in a mesh. The double nodes are nodes having the same coordinates w.r.t. precision given by `eps`. """ n_nod, dim = coor.shape cmap = find_map(coor, nm.zeros((0,dim)), allow_double=True) if cmap.size: output('double nodes in input mesh!') output('trying to fix...') while cmap.size: # Just like in Variable.equation_mapping()... ii = nm.argsort(cmap[:,1]) scmap = cmap[ii] eq = nm.arange(n_nod) eq[scmap[:,1]] = -1 eqi = eq[eq >= 0] eq[eqi] = nm.arange(eqi.shape[0]) remap = eq.copy() remap[scmap[:,1]] = eq[scmap[:,0]] output(coor.shape) coor = coor[eqi] ngroups = ngroups[eqi] output(coor.shape) ccs = [] for conn in conns: ccs.append(remap[conn]) conns = ccs cmap = find_map(coor, nm.zeros((0,dim)), allow_double=True) output('...done') return coor, ngroups, conns
[docs]def get_min_vertex_distance(coor, guess): """Can miss the minimum, but is enough for our purposes.""" # Sort by x. ix = nm.argsort(coor[:,0]) scoor = coor[ix] mvd = 1e16 # Get mvd in chunks potentially smaller than guess. n_coor = coor.shape[0] i0 = i1 = 0 x0 = scoor[i0,0] while 1: while ((scoor[i1,0] - x0) < guess) and (i1 < (n_coor - 1)): i1 += 1 ## print i0, i1, x0, scoor[i1,0] aim, aa1, aa2, aux = get_min_vertex_distance_naive(scoor[i0:i1+1]) if aux < mvd: im, a1, a2 = aim, aa1 + i0, aa2 + i0 mvd = min(mvd, aux) i0 = i1 = int(0.5 * (i1 + i0)) + 1 ## i0 += 1 x0 = scoor[i0,0] ## print '-', i0 if i1 == n_coor - 1: break ## print im, ix[a1], ix[a2], a1, a2, scoor[a1], scoor[a2] return mvd
[docs]def get_min_vertex_distance_naive(coor): ii = nm.arange(coor.shape[0]) i1, i2 = nm.meshgrid(ii, ii) i1 = i1.flatten() i2 = i2.flatten() ii = nm.where(i1 < i2) aux = coor[i1[ii]] - coor[i2[ii]] aux = nm.sum(aux**2.0, axis=1) im = aux.argmin() return im, i1[ii][im], i2[ii][im], nm.sqrt(aux[im])
[docs]def make_mesh(coor, ngroups, conns, mesh_in): """Create a mesh reusing mat_ids and descs of mesh_in.""" mat_ids = [] for ii, conn in enumerate(conns): mat_id = nm.empty((conn.shape[0],), dtype=nm.int32) mat_id.fill(mesh_in.mat_ids[ii][0]) mat_ids.append(mat_id) mesh_out = Mesh.from_data('merged mesh', coor, ngroups, conns, mat_ids, mesh_in.descs) return mesh_out
[docs]class Mesh(Struct): """ The Mesh class is a light proxy to CMesh. Input and output is handled by the MeshIO class and subclasses. """
[docs] @staticmethod def from_file(filename=None, io='auto', prefix_dir=None, omit_facets=False, file_format=None): """ Read a mesh from a file. Parameters ---------- filename : string or function or MeshIO instance or Mesh instance The name of file to read the mesh from. For convenience, a mesh creation function or a MeshIO instance or directly a Mesh instance can be passed in place of the file name. io : *MeshIO instance Passing *MeshIO instance has precedence over filename. prefix_dir : str If not None, the filename is relative to that directory. omit_facets : bool If True, do not read cells of lower dimension than the space dimension (faces and/or edges). Only some MeshIO subclasses support this! """ if isinstance(filename, Mesh): return filename if io == 'auto': if filename is None: output('filename or io must be specified!') raise ValueError else: io = MeshIO.any_from_filename(filename, prefix_dir=prefix_dir, file_format=file_format) output('reading mesh (%s)...' % io.filename) timer = Timer(start=True) trunk = io.get_filename_trunk() mesh = Mesh(trunk) mesh =, omit_facets=omit_facets) output('...done in %.2f s' % timer.stop()) mesh._set_shape_info() return mesh
[docs] @staticmethod def from_region(region, mesh_in, localize=False, is_surface=False, tdim=None): """ Create a mesh corresponding to cells, or lower dimensional entities according tdim parameter, of a given region. If `is_surface` is True, then tdim = dim - 1. """ cmesh_in = region.cmesh if is_surface: tdim = region.tdim - 1 if tdim is None: if not region.shape.n_cell: raise ValueError('region %s has no cells!' % cmesh = cmesh_in.create_new(region.cells, region.tdim, localize=localize) else: entity = region.entities[tdim] if len(region.entities[tdim]) < 1: raise ValueError(f'region {} has no entities' f' of dimension {tdim}!') cmesh = cmesh_in.create_new(entity, tdim, localize=localize) mesh = Mesh( + "_reg", cmesh=cmesh) if localize: remap = nm.empty(cmesh_in.n_coor, dtype=nm.int32) remap[:] = -1 remap[region.vertices] = nm.arange(region.vertices.shape[0]) mesh.nodal_bcs = {} for key, val in mesh_in.nodal_bcs.items(): new_val = remap[val] mesh.nodal_bcs[key] = new_val[new_val >= 0] else: mesh.nodal_bcs = mesh_in.nodal_bcs.copy() return mesh
[docs] @staticmethod def from_data(name, coors, ngroups, conns, mat_ids, descs, nodal_bcs=None): """ Create a mesh from mesh IO data. """ mesh = Mesh(name) mesh._set_io_data(coors=coors, ngroups=ngroups, conns=conns, mat_ids=mat_ids, descs=descs, nodal_bcs=nodal_bcs) mesh._set_shape_info() return mesh
def __init__(self, name='mesh', cmesh=None): """ Create a Mesh. By default, the mesh is empty, see the `cmesh` argument. Parameters ---------- name : str Object name. cmesh : CMesh, optional If given, use this as the cmesh. """ Struct.__init__(self, name=name, nodal_bcs={}, io=None) if cmesh is not None: self.cmesh_tdim = [None] * 4 self.cmesh = self.cmesh_tdim[cmesh.tdim] = cmesh self._collect_descs() self._coors = self.cmesh.coors self._set_shape_info()
[docs] def copy(self, name=None): """ Make a deep copy of the mesh. Parameters ---------- name : str Name of the copied mesh. """ if name is None: name = cmesh = self.cmesh.create_new() return Mesh(name=name, cmesh=cmesh)
def __add__(self, other): """ Merge the two meshes, assuming they have the same kind of the single element group. """ cmap = find_map(self.coors, other.coors) desc = self.descs[0] aux = merge_mesh(self.coors, self.cmesh.vertex_groups, self.get_conn(desc), self.cmesh.cell_groups, other.coors, other.cmesh.vertex_groups, other.get_conn(desc), other.cmesh.cell_groups, cmap) coors, ngroups, conn, mat_ids = aux mesh = Mesh.from_data( + ' + ' +, coors, ngroups, [conn], [mat_ids], [desc]) return mesh def _collect_descs(self): cmesh = self.cmesh i2k = invert_dict(cmesh.key_to_index) cts = nm.unique(cmesh.cell_types) self.descs = [i2k[ct] for ct in cts] def _set_shape_info(self): self.n_nod, self.dim = self.coors.shape self.n_el = nm.sum([cm.n_el for cm in self.cmesh_tdim if cm]) self.dims = [int(ii[0]) for ii in self.descs] def _set_io_data(self, coors, ngroups, conns, mat_ids, descs, nodal_bcs=None): """ Set mesh data. Parameters ---------- coors : array Coordinates of mesh nodes. ngroups : array Node groups. conns : list of arrays The array of mesh elements (connectivities) for each element group. mat_ids : list of arrays The array of material ids for each element group. descs: list of strings The element type for each element group. nodal_bcs : dict of arrays, optional The nodes defining regions for boundary conditions referred to by the dict keys in problem description files. """ ac = nm.ascontiguousarray coors = ac(coors, dtype=nm.float64) if ngroups is None: ngroups = nm.zeros((coors.shape[0],), dtype=nm.int32) self.descs = descs self.nodal_bcs = get_default(nodal_bcs, {}) max_tdim = nm.max([int(k[0]) for k in descs]) self.cmesh_tdim = [None] * 4 copy_coors = True from sfepy.discrete.common.extmods.cmesh import CMesh for idim in range(1, max_tdim + 1): conns_, mat_ids_, descs_ = [], [], [] for k, dsc in enumerate(descs): if int(dsc[0]) == idim: conns_.append(ac(conns[k], dtype=nm.int32)) mat_ids_.append(mat_ids[k]) descs_.append(dsc) if len(conns_) > 0: mat_ids_ = ac(nm.concatenate(mat_ids_)) self.cmesh_tdim[idim] = CMesh.from_data(coors, ac(ngroups), conns_, mat_ids_, descs_, copy_coors=copy_coors) if copy_coors: copy_coors = False coors = self.cmesh_tdim[idim].coors self.cmesh = self.cmesh_tdim[max_tdim] self._coors = self.cmesh.coors
[docs] def get_cmesh(self, desc): return self.cmesh_tdim[self.dims[self.descs.index(desc)]]
def _get_io_data(self, cell_dim_only=None): """ Return data to be used by `MeshIO`. """ conns, mat_ids = [], [] if cell_dim_only is not None: if not is_sequence(cell_dim_only): cell_dim_only = [cell_dim_only] descs = [ii for ii in self.descs if int(ii[0]) in cell_dim_only] else: descs = self.descs for desc in descs: cmesh = self.get_cmesh(desc) conn, cells = self.get_conn(desc, ret_cells=True, tdim=cmesh.tdim) conns.append(conn) mat_ids.append(cmesh.cell_groups[cells]) return cmesh.coors, cmesh.vertex_groups, conns, mat_ids, descs @property def coors(self): return self._coors
[docs] def write(self, filename=None, io=None, out=None, float_format=None, file_format=None, **kwargs): """ Write mesh + optional results in `out` to a file. Parameters ---------- filename : str, optional The file name. If None, the mesh name is used instead. io : MeshIO instance or 'auto', optional Passing 'auto' respects the extension of `filename`. out : dict, optional The output data attached to the mesh vertices and/or cells. float_format : str, optional The format string used to print floats in case of a text file format. **kwargs : dict, optional Additional arguments that can be passed to the `MeshIO` instance. """ if filename is None: filename = + '.mesh' if io is None: io = if io is None: io = 'auto' if io == 'auto': io = MeshIO.any_from_filename(filename, file_format=file_format, mode='w') io.set_float_format(float_format) io.write(filename, self, out, **kwargs)
[docs] def get_bounding_box(self): return nm.vstack((nm.amin(self.coors, 0), nm.amax(self.coors, 0)))
[docs] def get_conn(self, desc, ret_cells=False, tdim=None): """ Get the rectangular cell-vertex connectivity corresponding to `desc`. If `ret_cells` is True, the corresponding cells are returned as well. """ if desc not in self.descs: raise ValueError("'%s' not in %s!" % (desc, self.descs)) cmesh = self.cmesh if tdim is None else self.cmesh_tdim[tdim] cells = nm.where(cmesh.cell_types == cmesh.key_to_index[desc]) cells = cells[0].astype(nm.uint32) cdim = int(desc[0]) conn = cmesh.get_incident(0, cells, cdim) conn = conn.reshape((cells.shape[0], -1)).astype(nm.int32) if ret_cells: return conn, cells else: return conn
[docs] def transform_coors(self, mtx_t, ref_coors=None): """ Transform coordinates of the mesh by the given transformation matrix. Parameters ---------- mtx_t : array The transformation matrix `T` (2D array). It is applied depending on its shape: - `(dim, dim): x = T * x` - `(dim, dim + 1): x = T[:, :-1] * x + T[:, -1]` ref_coors : array, optional Alternative coordinates to use for the transformation instead of the mesh coordinates, with the same shape as `self.coors`. """ if ref_coors is None: ref_coors = self.coors if mtx_t.shape[1] > self.coors.shape[1]: self.coors[:] =, mtx_t[:,:-1].T) + mtx_t[:,-1] else: self.coors[:] =, mtx_t.T)
[docs] def create_conn_graph(self, verbose=True): """ Create a graph of mesh connectivity. Returns ------- graph : csr_matrix The mesh connectivity graph as a SciPy CSR matrix. """ from sfepy.discrete.common.extmods.cmesh import create_mesh_graph shape = (self.n_nod, self.n_nod) output('graph shape:', shape, verbose=verbose) if == 0: output('no graph (zero size)!', verbose=verbose) return None output('assembling mesh graph...', verbose=verbose) timer = Timer(start=True) conn = [self.get_conn(self.descs[0], tdim=k) for k in self.dims] nnz, prow, icol = create_mesh_graph(shape[0], shape[1], 1, conn, conn) output('...done in %.2f s' % timer.stop(), verbose=verbose) output('graph nonzeros: %d (%.2e%% fill)' \ % (nnz, 100.0 * float(nnz) /, verbose=verbose) data = nm.ones((nnz,), dtype=bool) graph = sp.csr_matrix((data, icol, prow), shape) return graph