# -*- coding: utf-8 -*- """ “Mesh”模块包含了有限元网格的不同类型. See the following implementations: * MeshTri,三角剖分网格 * MeshTet, 四面体剖分网格 * MeshQuad, 矩形剖分网格 * MeshHex, 六面体剖分网格 """ import matplotlib.pyplot as plt import matplotlib.tri as mtri import numpy as np import warnings from scipy.sparse import coo_matrix import skfem.mapping import skfem.element
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
class Mesh(): """A finite element mesh. This is an abstract superclass. See the following implementations: * MeshTri, triangular mesh * MeshTet, tetrahedral mesh * MeshQuad, quadrilateral mesh * MeshHex, hexahedral mesh """ refdom = "none" brefdom = "none" meshio_type = "none" p = np.array([]) t = np.array([])
1.1
def __init__(self): """Check that p and t are C_CONTIGUOUS as this leads to better performance.""" if self.p is not None: if self.p.flags['F_CONTIGUOUS']: if self.p.shape[1]>1000: warnings.warn("Mesh.__init__(): Transforming " + "over 100 vertices to C_CONTIGUOUS.") self.p = np.ascontiguousarray(self.p) if self.t is not None: if self.t.flags['F_CONTIGUOUS']: if self.t.shape[1]>1000: warnings.warn("Mesh.__init__(): Transforming " + "over 100 elements to C_CONTIGUOUS.") self.t = np.ascontiguousarray(self.t)
1.2
def __str__(self): return self.__repr__()
1.3
def __repr__(self): return "Mesh of type '" + str(type(self)) + "' " "with " + str(self.p.shape) + " vertices " "and " + str(self.t.shape) + " elements."
1.4
def show(self): """A wrapper包装纸 for matplotlib.pyplot.show().""" plt.show()
1.5
def dim(self): """Return the spatial空间的 dimension of the mesh.""" return int(self.p.shape[0])
1.6
def mapping(self): """Default不履行 local-to-global mapping for the mesh.""" raise NotImplementedError("Default mapping not implemented!")
1.7
def _uniform_refine(self): """Perform执行 a single uniform mesh refinement一致网格加密.""" raise NotImplementedError("Single refine not implemented " + "for this mesh type!")
1.8
def refine(self, N=None): """Refine the mesh. Parameters ---------- N : int (optional) Perform N refinements. """ if N is None: return self._uniform_refine() else: for itr in range(N): self._uniform_refine()
1.9
def remove_elements(self, element_indices): """Construct new mesh with elements removed based on their indices. Parameters ---------- element_indices : numpy array List of element indices to remove. Returns ------- skfem.Mesh A new mesh object with elements removed as per requested. """ keep = np.setdiff1d(np.arange(self.t.shape[1]), element_indices) newt = self.t[:, keep] ptix = np.unique(newt) reverse = np.zeros(self.p.shape[1]) reverse[ptix] = np.arange(len(ptix)) newt = reverse[newt] newp = self.p[:, ptix] meshclass = type(self) return meshclass(newp, newt.astype(np.intp))
1.10
def scale(self, scale): """Scale the mesh. Parameters ---------- scale : float OR tuple of size dim Scale each dimension by a factor. If a floating point number is provided, same scale is used for each dimension. """ for itr in range(int(self.dim())): if isinstance(scale, tuple): self.p[itr, :] *= scale[itr] else: self.p[itr, :] *= scale
1.11
def translate(self, vec): """Translate the mesh. Parameters ---------- vec : tuple of size dim Translate the mesh by a vector. """ for itr in range(int(self.dim())): self.p[itr, :] += vec[itr]
1.12
def _validate(self): """Perform mesh validity checks.""" # check that element connectivity contains integers # NOTE: this is neccessary for some plotting functionality if not np.issubdtype(self.t[0, 0], np.signedinteger): msg = ("Mesh._validate(): Element connectivity " "must consist of integers.") raise Exception(msg) # check that vertex matrix has "correct" size if self.p.shape[0] > 3: msg = ("Mesh._validate(): We do not allow meshes " "embedded into larger than 3-dimensional " "Euclidean space! Please check that " "the given vertex matrix is of size Ndim x Nvertices.") raise Exception(msg) # check that element connectivity matrix has correct size nvertices = {'line': 2, 'tri': 3, 'quad': 4, 'tet': 4, 'hex': 8} if self.t.shape[0] != nvertices[self.refdom]: msg = ("Mesh._validate(): The given connectivity " "matrix has wrong shape!") raise Exception(msg) # check that there are no duplicate points tmp = np.ascontiguousarray(self.p.T) if self.p.shape[1] != np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1])).shape[0]: msg = "Mesh._validate(): Mesh contains duplicate vertices." warnings.warn(msg) # check that all points are at least in some element if len(np.setdiff1d(np.arange(self.p.shape[1]), np.unique(self.t))): msg = ("Mesh._validate(): Mesh contains a vertex " "not belonging to any element.") raise Exception(msg)
1.13
def save(self, filename, pointData=None, cellData=None): """Export the mesh and fields using meshio. Parameters ---------- filename : string The filename for vtk-file. pointData : (optional) numpy array for one output or dict for multiple cellData : (optional) numpy array for one output or dict for multiple """ import meshio if pointData is not None: if type(pointData) != dict: pointData = {'0':pointData} if cellData is not None: if type(cellData) != dict: cellData = {'0':cellData} cells = { self.meshio_type : self.t.T } mesh = meshio.Mesh(self.p.T, cells, pointData, cellData) meshio.write(filename, mesh)
1.14
@classmethod def load(cls, filename): """Load a mesh from file using meshio. Parameters ---------- filename : string The filename for the mesh. """ import meshio mesh = meshio.read(filename) if issubclass(cls, Mesh2D): return cls(mesh.points[:, :2].T, mesh.cells[cls.meshio_type].T) else: return cls(mesh.points.T, mesh.cells[cls.meshio_type].T)
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
三维网格
class Mesh3D(Mesh): """Three dimensional meshes, common methods. See the following implementations: * MeshTet, tetrahedral mesh * MeshHex, hexahedral mesh """
2.1
def nodes_satisfying(self, test): """Return nodes that satisfy some condition. Parameters ---------- test : lambda function (3 params) Evaluates to 1 or True for nodes belonging to the output set. """ return np.nonzero(test(self.p[0, :], self.p[1, :], self.p[2, :]))[0]
2.2
def facets_satisfying(self, test): """Return facets whose midpoints satisfy some condition. Parameters ---------- test : lambda functions (3 params) Evaluates to 1 or True for facet midpoints of the facets belonging to the output set. """ mx = np.sum(self.p[0, self.facets], axis=0)/self.facets.shape[0] my = np.sum(self.p[1, self.facets], axis=0)/self.facets.shape[0] mz = np.sum(self.p[2, self.facets], axis=0)/self.facets.shape[0] return np.nonzero(test(mx, my, mz))[0]
2.3
def edges_satisfying(self, test): """Return edges whose midpoints satisfy some condition. Parameters ---------- test : lambda functions (3 params) Evaluates to 1 or True for edge midpoints of the edges belonging to the output set. """ mx = 0.5*(self.p[0, self.edges[0, :]] + self.p[0, self.edges[1, :]]) my = 0.5*(self.p[1, self.edges[0, :]] + self.p[1, self.edges[1, :]]) mz = 0.5*(self.p[2, self.edges[0, :]] + self.p[2, self.edges[1, :]]) return np.nonzero(test(mx, my, mz))[0]
2.4
def boundary_nodes(self): """Return an array of boundary node indices.""" return np.unique(self.facets[:, self.boundary_facets()])
2.5
def boundary_facets(self): """Return an array of boundary facet indices.""" return np.nonzero(self.f2t[1, :] == -1)[0]
2.6
def interior_facets(self): """Return an array of interior facet indices.""" return np.nonzero(self.f2t[1, :] >= 0)[0]
2.7
def boundary_edges(self): """Return an array of boundary edge indices.""" bnodes = self.boundary_nodes()[:, None] return np.nonzero(np.sum(self.edges[0, :] == bnodes, axis=0) * np.sum(self.edges[1, :] == bnodes, axis=0))[0]
2.8
def interior_nodes(self): """Return an array of interior node indices.""" return np.setdiff1d(np.arange(0, self.p.shape[1]), self.boundary_nodes())
2.9
def draw_vertices(self): """Draw all vertices using mplot3d.""" fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.scatter(self.p[0, :], self.p[1, :], self.p[2, :]) return fig
2.10
def param(self): """Return (maximum) mesh parameter.""" return np.max(np.sqrt(np.sum((self.p[:, self.edges[0, :]] - self.p[:, self.edges[1, :]])**2, axis=0)))
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
二维网格
class Mesh2D(Mesh): """Two dimensional meshes, common methods. See the following implementations: * MeshTri, triangular mesh * MeshQuad, quadrilateral mesh """
3.1
def jiggle(self, z=0.2): """Jiggle the interior nodes of the mesh. Parameters ---------- z : (optional, default=0.2) float Mesh parameter is multiplied by this number. The resulting number corresponds to the standard deviation of the jiggle. """ y = z*self.param() I = self.interior_nodes() self.p[0, I] = self.p[0, I] + y*np.random.rand(len(I)) self.p[1, I] = self.p[1, I] + y*np.random.rand(len(I))
3.2
def boundary_nodes(self): """Return an array of boundary node indices.""" return np.unique(self.facets[:, self.boundary_facets()])
3.3
def nodes_satisfying(self, test): """Return nodes that satisfy some condition. Parameters ---------- test : lambda An anonymous function with two parameters (x and y) and which returns True for the set of nodes that are to be included in the return set. Returns ------- numpy array An array of node indices. """ return np.nonzero(test(self.p[0, :], self.p[1, :]))[0]
3.4
def draw_nodes(self, nodes, mark='bo'): """Highlight some nodes. Parameters ---------- nodes : numpy array The indices of the nodes to highlight. mark : (optional, default='bo') string A standard matplotlib string to define the highlight style. """ plt.plot(self.p[0, nodes], self.p[1, nodes], mark)
3.5
def param(self): """Return mesh parameter.""" return np.max(np.sqrt(np.sum((self.p[:, self.facets[0, :]] - self.p[:, self.facets[1, :]])**2, axis=0)))
3.6
def interior_nodes(self): """Return an array of interior node indices.""" return np.setdiff1d(np.arange(0, self.p.shape[1]), self.boundary_nodes())
3.7
def facets_satisfying(self, test): """Return facets whose midpoints satisfy some condition. Parameters ---------- test : lambda function (2 params) An anonymous function with two parameters (x and y) and which returns True for the midpoints of the set of facets that are to be included in the return set. """ mx = 0.5*(self.p[0, self.facets[0, :]] + self.p[0, self.facets[1, :]]) my = 0.5*(self.p[1, self.facets[0, :]] + self.p[1, self.facets[1, :]]) return np.nonzero(test(mx, my))[0]
3.8
def elements_satisfying(self, test): """Return elements whose midpoints satisfy some condition. Parameters ---------- test : lambda function (2 params) An anonymous function with two parameters (x and y) and which returns True for the midpoints of the set of elements that are to be included in the return set. """ mx = np.sum(self.p[0, self.t], axis=0)/self.t.shape[0] my = np.sum(self.p[1, self.t], axis=0)/self.t.shape[0] return np.nonzero(test(mx, my))[0]
3.9
def interior_facets(self): """Return an array of interior facet indices.""" return np.nonzero(self.f2t[1, :] >= 0)[0]
3.10
def boundary_facets(self): """Return an array of boundary facet indices.""" return np.nonzero(self.f2t[1, :] == -1)[0]
3.11
def draw(self, ax=None, node_numbering=False, facet_numbering=False, element_numbering=False): """Draw the mesh. Parameters ---------- ax : (optional, default=None) matplotlib axis Use a predefined axis for plotting. node_numbering : (optional, default=False) Draw node numbering. facet_numbering: (optional, default=False) Draw facet numbering. element_numbering : (optional, default=False) Draw element numbering. """ if ax is None: # create new figure fig = plt.figure() ax = fig.add_subplot(111) # visualize the mesh faster plotting is achieved through # None insertion trick. xs = [] ys = [] for s, t, u, v in zip(self.p[0, self.facets[0, :]], self.p[1, self.facets[0, :]], self.p[0, self.facets[1, :]], self.p[1, self.facets[1, :]]): xs.append(s) xs.append(u) xs.append(None) ys.append(t) ys.append(v) ys.append(None) ax.plot(xs, ys, 'k', linewidth='0.5') if node_numbering: for itr in range(self.p.shape[1]): ax.text(self.p[0, itr], self.p[1, itr], str(itr)) if facet_numbering: mx = .5*(self.p[0, self.facets[0, :]] + self.p[0, self.facets[1, :]]) my = .5*(self.p[1, self.facets[0, :]] + self.p[1, self.facets[1, :]]) for itr in range(self.facets.shape[1]): ax.text(mx[itr], my[itr], str(itr)) if element_numbering: mx = np.sum(self.p[0, self.t], axis=0)/self.t.shape[0] my = np.sum(self.p[1, self.t], axis=0)/self.t.shape[0] for itr in range(self.t.shape[1]): ax.text(mx[itr], my[itr], str(itr)) return ax
3.12
def mirror_mesh(self, a, b, c): """Mirror a mesh by the line ax + by + c = 0. Parameters ---------- a : float b : float c : float """ tmp = -2.0*(a*self.p[0, :] + b*self.p[1, :] + c)/(a**2 + b**2) newx = a*tmp + self.p[0, :] newy = b*tmp + self.p[1, :] newpoints = np.vstack((newx, newy)) points = np.hstack((self.p, newpoints)) tris = np.hstack((self.t, self.t + self.p.shape[1])) # remove duplicates tmp = np.ascontiguousarray(points.T) tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)]*tmp.shape[1]), return_index=True, return_inverse=True) points = points[:, ixa] tris = ixb[tris] meshclass = type(self) return meshclass(points, tris)
3.13
def save(self, filename, pointData=None, cellData=None): """Export the mesh and fields using meshio. (2D version.) Parameters ---------- filename : string The filename for vtk-file. pointData : (optional) numpy array for one output or dict for multiple cellData : (optional) numpy array for one output or dict for multiple """ import meshio # a row of zeros required in 2D p = np.vstack((np.zeros(self.p.shape[1]), self.p)) if pointData is not None: if type(pointData) != dict: pointData = {'0':pointData} if cellData is not None: if type(cellData) != dict: cellData = {'0':cellData} cells = { self.meshio_type : self.t.T } mesh = meshio.Mesh(p.T, cells, pointData, cellData) meshio.write(filename, mesh)
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
class InterfaceMesh1D(Mesh): """An interface mesh for mortar methods."""
4.1
def __init__(self, mesh1, mesh2, rule, param, debug_plot=False): self.brefdom = mesh1.brefdom p1_ix = mesh1.nodes_satisfying(rule) p2_ix = mesh2.nodes_satisfying(rule) p1 = mesh1.p[:, p1_ix] p2 = mesh2.p[:, p2_ix] _, ix = np.unique(np.concatenate((param(p1[0, :], p1[1, :]), param(p2[0, :], p2[1, :]))), return_index=True) np1 = mesh1.p.shape[1] nt1 = mesh1.t.shape[1] ixorig = np.concatenate((p1_ix, p2_ix + np1))[ix] self.p = np.hstack((mesh1.p, mesh2.p)) self.t = np.hstack((mesh1.t, mesh2.t + np1)) self.facets = np.array([ixorig[:-1], ixorig[1:]]) self.t2f = -1 + 0*np.hstack((mesh1.t2f, mesh2.t2f)) # construct normals tangent_x = self.p[0, self.facets[0, :]] - self.p[0, self.facets[1, :]] tangent_y = self.p[1, self.facets[0, :]] - self.p[1, self.facets[1, :]] tangent_lengths = np.sqrt(tangent_x**2 + tangent_y**2) self.normals = np.array([-tangent_y/tangent_lengths, tangent_x/tangent_lengths]) if debug_plot: ax = mesh1.draw() mesh2.draw(ax=ax) xs = np.array([self.p[0, self.facets[0, :]], self.p[0, self.facets[1, :]]]) midx = np.sum(xs, axis=0)/2.0 ys = np.array([self.p[1, self.facets[0, :]], self.p[1, self.facets[1, :]]]) midy = np.sum(ys, axis=0)/2.0 xs = 0.9*(xs - midx) + midx ys = 0.9*(ys - midy) + midy ax.plot(xs, ys, 'x-') # mappings from facets to the original triangles # TODO vectorize self.f2t = self.facets*0-1 for itr in range(self.facets.shape[1]): mx = .5*(self.p[0, self.facets[0, itr]] + self.p[0, self.facets[1, itr]]) my = .5*(self.p[1, self.facets[0, itr]] + self.p[1, self.facets[1, itr]]) val = param(mx, my) for jtr in mesh1.boundary_facets(): fix1 = mesh1.facets[0, jtr] x1 = mesh1.p[0, fix1] y1 = mesh1.p[1, fix1] fix2 = mesh1.facets[1, jtr] x2 = mesh1.p[0, fix2] y2 = mesh1.p[1, fix2] if rule(x1, y1) > 0 or rule(x2, y2) > 0: if val > param(x1, y1) and val < param(x2, y2): # OK self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val < param(x1, y1) and val > param(x2, y2): # ye olde # OK self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val >= param(x1, y1) and val < param(x2, y2): # OK self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val > param(x1, y1) and val <= param(x2, y2): # OK self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val <= param(x1, y1) and val > param(x2, y2): # OK self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val < param(x1, y1) and val >= param(x2, y2): # OK self.f2t[0, itr] = mesh1.f2t[0, jtr] break for jtr in mesh2.boundary_facets(): fix1 = mesh2.facets[0, jtr] x1 = mesh2.p[0, fix1] y1 = mesh2.p[1, fix1] fix2 = mesh2.facets[1, jtr] x2 = mesh2.p[0, fix2] y2 = mesh2.p[1, fix2] if rule(x1, y1) > 0 or rule(x2, y2) > 0: if val > param(x1, y1) and val < param(x2, y2): # OK self.f2t[1, itr] = mesh2.f2t[0, jtr] + nt1 break elif val < param(x1, y1) and val > param(x2, y2): # OK self.f2t[1, itr] = mesh2.f2t[0, jtr] + nt1 break elif val >= param(x1, y1) and val < param(x2, y2): # OK self.f2t[1, itr] = mesh2.f2t[0, jtr] + nt1 break elif val > param(x1, y1) and val <= param(x2, y2): # OK self.f2t[1, itr] = mesh2.f2t[0, jtr] + nt1 break elif val <= param(x1, y1) and val > param(x2, y2): # OK self.f2t[1, itr] = mesh2.f2t[0, jtr] + nt1 break elif val < param(x1, y1) and val >= param(x2, y2): # OK self.f2t[1, itr] = mesh2.f2t[0, jtr] + nt1 break if (self.f2t>-1).all(): self.f2t[0, :] return else: print(self.f2t) raise Exception("All mesh facets corresponding to mortar facets not found!")
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
一维网格
class MeshLine(Mesh): """One-dimensional mesh.""" refdom = "line" brefdom = "point" p = np.array([]) t = np.array([])
5.1
def __init__(self, p=None, t=None, validate=True): if p is None and t is None: p = np.array([[0, 1]]) t = np.array([[0], [1]]) elif p is not None and t is None: t = np.array([np.arange(np.max(p.shape)-1), np.arange(np.max(p.shape)-1)+1]) if len(p.shape)==1: p = np.array([p]) self.p = p self.t = t if validate: self._validate() super(MeshLine, self).__init__()
5.2
@classmethod def init_refdom(cls): """Initialize a mesh constisting of the reference interval [0,1].""" return cls()
5.3
def adaptive_refine(self, marked): """Perform an adaptive refine which splits each marked element into two.""" t = self.t p = self.p mid = range(len(marked)) + np.max(t) + 1 nonmarked = np.setdiff1d(np.arange(t.shape[1]), marked) newp = np.hstack((p, 0.5*(p[:, self.t[0, marked]] + p[:, self.t[1, marked]]))) newt = np.vstack((t[0, marked], mid)) newt = np.hstack((t[:, nonmarked], newt, np.vstack((mid, t[1, marked])))) # update fields self.p = newp self.t = newt
5.4
def _uniform_refine(self): """Perform a single mesh refine that halves 'h'.""" # rename variables t = self.t p = self.p mid = range(self.t.shape[1]) + np.max(t) + 1 # new vertices and elements newp = np.hstack((p, 0.5*(p[:, self.t[0, :]] + p[:, self.t[1, :]]))) newt = np.vstack((t[0, :], mid)) newt = np.hstack((newt, np.vstack((mid, t[1, :])))) # update fields self.p = newp self.t = newt
5.5
# TODO implement prolongation def boundary_nodes(self): """Find the boundary nodes of the mesh.""" _, counts = np.unique(self.t.flatten(), return_counts=True) return np.nonzero(counts == 1)[0]
5.6
def interior_nodes(self): """Find the interior nodes of the mesh.""" _, counts = np.unique(self.t.flatten(), return_counts=True) return np.nonzero(counts == 2)[0]
5.7
def plot(self, u, ax=None, color='ko-'): """Plot a function defined at the nodes of the mesh.""" if ax is None: # create new figure fig = plt.figure() ax = fig.add_subplot(111) xs = [] ys = [] for y1, y2, s, t in zip(u[self.t[0, :]], u[self.t[1, :]], self.p[0, self.t[0, :]], self.p[0, self.t[1, :]]): xs.append(s) xs.append(t) xs.append(None) ys.append(y1) ys.append(y2) ys.append(None) ax.plot(xs, ys, color) return ax
5.8
def __mul__(self, other): """Tensor product mesh.""" return MeshQuad.init_tensor(self.p[0, :], other.p[0, :])
5.9
def param(self): return np.max(np.abs(self.p[0, self.t[1, :]] - self.p[0, self.t[0, :]]))
5.10
def mapping(self): return skfem.mapping.MappingAffine(self)
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
矩形剖分
class MeshQuad(Mesh2D): """A mesh consisting of quadrilateral elements. Attributes ---------- p : numpy array of size 2 x Nvertices The vertices of the mesh t : numpy array of size 4 x Nelements The element connectivity facets : numpy array of size 2 x Nfacets Each column contains a pair of indices to p. f2t : numpy array of size 2 x Nfacets Each column contains a pair of indices to t or -1 on the second row if the facet is on the boundary. t2f : numpy array of size 4 x Nelements Each column contains four indices to facets. """ refdom = "quad" brefdom = "line" meshio_type = "quad" p = np.array([]) t = np.array([]) facets = np.array([]) f2t = np.array([]) t2f = np.array([])
6.1
def __init__(self, p=None, t=None, validate=True): """Initialize a quadrilateral mesh. Parameters ---------- p : (optional) numpy array of size 2 x Nvertices The points of the mesh. t : (optional) numpy array of size 4 x Nelements The element connectivity, i.e. indices to p. These should be in counter-clockwise order. """ if p is None and t is None: p = np.array([[0., 0.], [1., 0.], [1., 1.], [0., 1.]]).T t = np.array([[0, 1, 2, 3]]).T elif p is None or t is None: raise Exception("Must provide p AND t or neither") self.p = p self.t = t if validate: self._validate() self._build_mappings() super(MeshQuad, self).__init__()
6.2
@classmethod def init_tensor(cls, x, y): """Initialize a tensor product mesh. Parameters ---------- x : numpy array (1d) The nodal coordinates in dimension x y : numpy array (1d) The nodal coordinates in dimension y """ npx = len(x) npy = len(y) X, Y = np.meshgrid(np.sort(x), np.sort(y)) p = np.vstack((X.flatten('F'), Y.flatten('F'))) ix = np.arange(npx*npy) ne = (npx-1)*(npy-1) t = np.zeros((4, ne)) ix = ix.reshape(npy, npx, order='F').copy() t[0, :] = ix[0:(npy-1), 0:(npx-1)].reshape(ne, 1, order='F').copy().flatten() t[1, :] = ix[1:npy, 0:(npx-1)].reshape(ne, 1, order='F').copy().flatten() t[2, :] = ix[1:npy, 1:npx].reshape(ne, 1, order='F').copy().flatten() t[3, :] = ix[0:(npy-1), 1:npx].reshape(ne, 1, order='F').copy().flatten() return cls(p, t.astype(np.int64))
6.3
@classmethod def init_refdom(cls): """Initialize a mesh that includes only the reference quad. The mesh topology is as follows: (-1,1) *-------------* (1,1) | | | | | | | | | | | | | | (-1,-1) *-------------* (1,-1) """ p = np.array([[-1., -1.], [1., -1.], [1., 1.], [-1., 1.]]).T t = np.array([[0, 1, 2, 3]]).T return cls(p, t)
6.4
def _build_mappings(self): # do not sort since order defines counterclockwise order # self.t=np.sort(self.t,axis=0) # define facets: in the order (0,1) (1,2) (2,3) (0,3) self.facets = np.sort(np.vstack((self.t[0, :], self.t[1, :])), axis=0) self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[1, :], self.t[2, :])), axis=0))) self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[2, :], self.t[3, :])), axis=0))) self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[0, :], self.t[3, :])), axis=0))) # get unique facets and build quad-to-facet mapping: 4 (edges) x Nquads tmp = np.ascontiguousarray(self.facets.T) tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)]*tmp.shape[1]), return_index=True, return_inverse=True) self.facets = self.facets[:, ixa] self.t2f = ixb.reshape((4, self.t.shape[1])) # build facet-to-quadrilateral mapping: 2 (quads) x Nedges e_tmp = np.hstack((self.t2f[0, :], self.t2f[1, :], self.t2f[2, :], self.t2f[3, :])) t_tmp = np.tile(np.arange(self.t.shape[1]), (1, 4))[0] e_first, ix_first = np.unique(e_tmp, return_index=True) # this emulates matlab unique(e_tmp,'last') e_last, ix_last = np.unique(e_tmp[::-1], return_index=True) ix_last = e_tmp.shape[0] - ix_last - 1 self.f2t = np.zeros((2, self.facets.shape[1]), dtype=np.int64) self.f2t[0, e_first] = t_tmp[ix_first] self.f2t[1, e_last] = t_tmp[ix_last] # second row to -1 if repeated (i.e., on boundary) self.f2t[1, np.nonzero(self.f2t[0, :] == self.f2t[1, :])[0]] = -1
6.5
def _uniform_refine(self): """Perform a single mesh refine that halves 'h'. Each quadrilateral is split into four. """ # rename variables t = self.t p = self.p e = self.facets sz = p.shape[1] t2f = self.t2f + sz # quadrilateral middle point mid = range(self.t.shape[1]) + np.max(t2f) + 1 # new vertices are the midpoints of edges ... newp1 = 0.5*np.vstack((p[0, e[0, :]] + p[0, e[1, :]], p[1, e[0, :]] + p[1, e[1, :]])) # ... and element middle points newp2 = 0.25*np.vstack((p[0, t[0, :]] + p[0, t[1, :]] + p[0, t[2, :]] + p[0, t[3, :]], p[1, t[0, :]] + p[1, t[1, :]] + p[1, t[2, :]] + p[1, t[3, :]])) newp = np.hstack((p, newp1, newp2)) # build new quadrilateral definitions newt = np.vstack((t[0, :], t2f[0, :], mid, t2f[3, :])) newt = np.hstack((newt, np.vstack((t2f[0, :], t[1, :], t2f[1, :], mid)))) newt = np.hstack((newt, np.vstack((mid, t2f[1, :], t[2, :], t2f[2, :])))) newt = np.hstack((newt, np.vstack((t2f[3, :], mid, t2f[2, :], t[3, :])))) # update fields self.p = newp self.t = newt self._build_mappings()
6.6
# TODO implement prolongation def _splitquads(self, x=None): """Split each quad into two triangles and return MeshTri.""" t = self.t[[0, 1, 3], :] t = np.hstack((t, self.t[[1, 2, 3]])) if x is not None: if len(x) == self.t.shape[1]: # preserve elemental constant functions X = np.concatenate((x, x)) else: raise Exception("The parameter x must have one value per element.") return MeshTri(self.p, t, validate=False), X else: return MeshTri(self.p, t, validate=False)
6.7
def _splitquads_symmetric(self): """Split quads into four triangles.""" t = np.vstack((self.t, np.arange(self.t.shape[1]) + self.p.shape[1])) newt = t[[0, 1, 4], :] newt = np.hstack((newt, t[[1, 2, 4], :])) newt = np.hstack((newt, t[[2, 3, 4], :])) newt = np.hstack((newt, t[[3, 0, 4], :])) mx = np.sum(self.p[0, self.t], axis=0)/self.t.shape[0] my = np.sum(self.p[1, self.t], axis=0)/self.t.shape[0] return MeshTri(np.hstack((self.p, np.vstack((mx, my)))), newt, validate=False)
6.8
def plot(self, z, smooth=False, edgecolors=None, ax=None, zlim=None): """Visualize nodal or elemental function (2d). The quadrilateral mesh is split into triangular mesh (MeshTri) and the respective plotting function for the triangular mesh is used. """ m, z = self._splitquads(z) return m.plot(z, smooth, ax=ax, zlim=zlim, edgecolors=edgecolors)
6.9
def plot3(self, z, smooth=False, ax=None): """Visualize nodal function (3d i.e. three axes). The quadrilateral mesh is split into triangular mesh (MeshTri) and the respective plotting function for the triangular mesh is used. """ m, z = self._splitquads(z) return m.plot3(z, smooth, ax=ax)
6.10
def mapping(self): return skfem.mapping.MappingIsoparametric(self, skfem.element.ElementQ1())
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
六面体剖分
class MeshHex(Mesh3D): """A mesh consisting of hexahedral elements. Attributes ---------- p : numpy array of size 3 x Nvertices The vertices of the mesh. Each column corresponds to a point. t : numpy array of size 8 x Nelements The element connectivity. Each column corresponds to a element and contains eight column indices to MeshHex.p. facets : numpy array of size 4 x Nfacets Each column contains four column indices to MeshHex.p. f2t : numpy array of size 2 x Nfacets Each column contains a pair of column indices to MeshHex.t or -1 on the second row if the corresponding facet is located on the boundary. t2f : numpy array of size 6 x Nelements Each column contains four indices to MeshHex.facets. edges : numpy array of size 2 x Nedges Each column corresponds to an edge and contains two indices to MeshHex.p. t2e : numpy array of size 12 x Nelements Each column contains twelve column indices of MeshHex.edges. """ refdom = "hex" brefdom = "quad" meshio_type = "hexahedron"
7.1
def __init__(self, p=None, t=None, validate=True): if p is None and t is None: p = np.array([[0., 0., 0.], [0., 0., 1.], [0., 1., 0.], [1., 0., 0.], [0., 1., 1.], [1., 0., 1.], [1., 1., 0.], [1., 1., 1.]]).T t = np.array([[0, 1, 2, 3, 4, 5, 6, 7]]).T elif p is None or t is None: raise Exception("Must provide p AND t or neither") # # TODO fix orientation if p and t is provided. refer to # the default mesh for correct orientation # # 2---6 # / /| # 4---7 3 # | |/ # 1---5 # # The hidden node is 0. # self.p = p self.t = t if validate: self._validate() self._build_mappings() super(MeshHex, self).__init__()
7.2
@classmethod def init_tensor(cls, x, y, z): """Initialize a tensor product mesh. Parameters ---------- x : numpy array (1d) The nodal coordinates in dimension x y : numpy array (1d) The nodal coordinates in dimension y z : numpy array (1d) The nodal coordinates in dimension z """ npx = len(x) npy = len(y) npz = len(z) X, Y, Z = np.meshgrid(np.sort(x), np.sort(y), np.sort(z)) p = np.vstack((X.flatten('F'), Y.flatten('F'), Z.flatten('F'))) ix = np.arange(npx*npy*npz) ne = (npx-1)*(npy-1)*(npz-1) t = np.zeros((8, ne)) ix = ix.reshape(npy, npx, npz, order='F').copy() t[0, :] = ix[0:(npy-1), 0:(npx-1), 0:(npz-1)].reshape(ne, 1, order='F').copy().flatten() t[1, :] = ix[1:npy, 0:(npx-1), 0:(npz-1)].reshape(ne, 1, order='F').copy().flatten() t[2, :] = ix[0:(npy-1), 1:npx, 0:(npz-1)].reshape(ne, 1, order='F').copy().flatten() t[3, :] = ix[0:(npy-1), 0:(npx-1), 1:npz].reshape(ne, 1, order='F').copy().flatten() t[4, :] = ix[1:npy, 1:npx, 0:(npz-1)].reshape(ne, 1, order='F').copy().flatten() t[5, :] = ix[1:npy, 0:(npx-1), 1:npz].reshape(ne, 1, order='F').copy().flatten() t[6, :] = ix[0:(npy-1), 1:npx, 1:npz].reshape(ne, 1, order='F').copy().flatten() t[7, :] = ix[1:npy, 1:npx, 1:npz].reshape(ne, 1, order='F').copy().flatten() return cls(p, t.astype(np.int64))
7.3
def _build_mappings(self): """Build element-to-facet, element-to-edges, etc. mappings.""" self.edges = np.sort(np.vstack((self.t[0, :], self.t[1, :])), axis=0) e = np.array([0, 2, 0, 3, 1, 4, 1, 5, 2, 4, 2, 6, 3, 5, 3, 6, 4, 7, 5, 7, 6, 7]) # see the picture in init for i in range(11): self.edges = np.hstack((self.edges, np.sort(np.vstack((self.t[e[2*i], :], self.t[e[2*i+1], :])), axis=0))) # unique edges self.edges, ixa, ixb = np.unique(self.edges, axis=1, return_index=True, return_inverse=True) self.edges = np.ascontiguousarray(self.edges) self.t2e = ixb.reshape((12, self.t.shape[1])) # define facets self.facets = np.sort(np.vstack((self.t[0, :], self.t[1, :], self.t[4, :], self.t[2, :])), axis=0) f = np.array([0, 3, 6, 2, 0, 1, 5, 3, 2, 6, 7, 4, 1, 5, 7, 4, 3, 5, 7, 6]) for i in range(5): self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[f[4*i], :], self.t[f[4*i+1], :], self.t[f[4*i+2], :], self.t[f[4*i+3], :])), axis=0))) # unique facets self.facets, ixa, ixb = np.unique(self.facets, axis=1, return_index=True, return_inverse=True) self.facets = np.ascontiguousarray(self.facets) self.t2f = ixb.reshape((6, self.t.shape[1])) # build facet-to-hexa mapping: 2 (hexes) x Nfacets e_tmp = np.hstack((self.t2f[0, :], self.t2f[1, :], self.t2f[2, :], self.t2f[3, :], self.t2f[4, :], self.t2f[5, :])) t_tmp = np.tile(np.arange(self.t.shape[1]), (1, 6))[0] e_first, ix_first = np.unique(e_tmp, return_index=True) e_last, ix_last = np.unique(e_tmp[::-1], return_index=True) ix_last = e_tmp.shape[0] - ix_last - 1 self.f2t = np.zeros((2, self.facets.shape[1]), dtype=np.int64) self.f2t[0, e_first] = t_tmp[ix_first] self.f2t[1, e_last] = t_tmp[ix_last] ## second row to zero if repeated (i.e., on boundary) self.f2t[1, np.nonzero(self.f2t[0, :] == self.f2t[1, :])[0]] = -1
7.4
def _uniform_refine(self): """Perform a single mesh refine that halves 'h'. Each hex is split into 8.""" # rename variables t = self.t p = self.p e = self.edges f = self.facets sz = p.shape[1] t2e = self.t2e + sz t2f = self.t2f + np.max(t2e) + 1 # hex middle point mid = range(self.t.shape[1]) + np.max(t2f) + 1 # new vertices are the midpoints of edges ... newp1 = 0.5*np.vstack((p[0, e[0, :]] + p[0, e[1, :]], p[1, e[0, :]] + p[1, e[1, :]], p[2, e[0, :]] + p[2, e[1, :]])) # ... midpoints of facets ... newp2 = 0.25*np.vstack((p[0, f[0, :]] + p[0, f[1, :]] + p[0, f[2, :]] + p[0, f[3, :]], p[1, f[0, :]] + p[1, f[1, :]] + p[1, f[2, :]] + p[1, f[3, :]], p[2, f[0, :]] + p[2, f[1, :]] + p[2, f[2, :]] + p[2, f[3, :]])) # ... and element middle points newp3 = 0.125*np.vstack((p[0, t[0, :]] + p[0, t[1, :]] + p[0, t[2, :]] + p[0, t[3, :]] + p[0, t[4, :]] + p[0, t[5, :]] + p[0, t[6, :]] + p[0, t[7, :]], p[1, t[0, :]] + p[1, t[1, :]] + p[1, t[2, :]] + p[1, t[3, :]] + p[1, t[4, :]] + p[1, t[5, :]] + p[1, t[6, :]] + p[1, t[7, :]], p[2, t[0, :]] + p[2, t[1, :]] + p[2, t[2, :]] + p[2, t[3, :]] + p[2, t[4, :]] + p[2, t[5, :]] + p[2, t[6, :]] + p[2, t[7, :]])) newp = np.hstack((p, newp1, newp2, newp3)) # build new hex indexing (this requires some serious meditation) newt = np.vstack((t[0, :], t2e[0, :], t2e[1, :], t2e[2, :], t2f[0, :], t2f[2, :], t2f[1, :], mid)) newt = np.hstack((newt, np.vstack((t2e[0, :], t[1, :], t2f[0, :], t2f[2, :], t2e[3, :], t2e[4, :], mid, t2f[4, :])))) newt = np.hstack((newt, np.vstack((t2e[1, :], t2f[0, :], t[2, :], t2f[1, :], t2e[5, :], mid, t2e[6, :], t2f[3, :] )))) newt = np.hstack((newt, np.vstack((t2e[2, :], t2f[2, :], t2f[1, :], t[3, :], mid, t2e[7, :], t2e[8, :], t2f[5, :] )))) newt = np.hstack((newt, np.vstack((t2f[0, :], t2e[3, :], t2e[5, :], mid, t[4, :], t2f[4, :], t2f[3, :], t2e[9, :] )))) newt = np.hstack((newt, np.vstack((t2f[2, :], t2e[4, :], mid, t2e[7, :], t2f[4, :], t[5, :], t2f[5, :], t2e[10, :], )))) newt = np.hstack((newt, np.vstack((t2f[1, :], mid, t2e[6, :], t2e[8, :], t2f[3, :], t2f[5, :], t[6, :], t2e[11, :] )))) newt = np.hstack((newt, np.vstack((mid, t2f[4, :], t2f[3, :], t2f[5, :], t2e[9, :], t2e[10, :], t2e[11, :], t[7, :] )))) # update fields self.p = newp self.t = newt self._build_mappings()
7.5
# TODO implement prolongation def save(self, filename, pointData=None, cellData=None): """Export the mesh and fields using meshio. (Hexahedron version.) Parameters ---------- filename : string The filename for vtk-file. pointData : (optional) numpy array for one output or dict for multiple cellData : (optional) numpy array for one output or dict for multiple """ import meshio # vtk requires a different ordering t = self.t[[0, 3, 6, 2, 1, 5, 7, 4], :] if pointData is not None: if type(pointData) != dict: pointData = {'0':pointData} if cellData is not None: if type(cellData) != dict: cellData = {'0':cellData} cells = { 'hexahedron' : t.T } mesh = meshio.Mesh(self.p.T, cells, pointData, cellData) meshio.write(filename, mesh)
7.6
def mapping(self): return skfem.mapping.MappingIsoparametric(self, skfem.element.ElementHex1())
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
四边体剖分
class MeshTet(Mesh3D): """A mesh consisting of tetrahedral elements. Attributes ---------- p : numpy array of size 3 x Nvertices The vertices of the mesh. Each column corresponds to a point. t : numpy array of size 4 x Nelements The element connectivity. Each column corresponds to a element and contains four column indices to MeshTet.p. facets : numpy array of size 3 x Nfacets Each column contains a triplet of column indices to MeshTet.p. Order: (0, 1, 2) (0, 1, 3) (0, 2, 3) (1, 2, 3) f2t : numpy array of size 2 x Nfacets Each column contains a pair of column indices to MeshTet.t or -1 on the second row if the facet is located on the boundary. t2f : numpy array of size 4 x Nelements Each column contains four indices to MeshTet.facets. edges : numpy array of size 2 x Nedges Each column corresponds to an edge and contains two indices to MeshTet.p. Order: (0, 1) (1, 2) (0, 2) (0, 3) (1, 3) (2, 3) t2e : numpy array of size 6 x Nelements Each column contains six indices to MeshTet.edges. Examples -------- Read a tetrahedral mesh from gmsh-format using meshio. >>> m = MeshTet.load('examples/box.msh') >>> m.p.shape (3, 358) """ refdom = "tet" brefdom = "tri" meshio_type = "tetra"
8.1
def __init__(self, p=None, t=None, validate=True): if p is None and t is None: p = np.array([[0., 0., 0.], [0., 0., 1.], [0., 1., 0.], [1., 0., 0.], [0., 1., 1.], [1., 0., 1.], [1., 1., 0.], [1., 1., 1.]]).T t = np.array([[0, 1, 2, 3], [3, 5, 1, 7], [2, 3, 6, 7], [2, 3, 1, 7], [1, 2, 4, 7]]).T elif p is None or t is None: raise Exception("Must provide p AND t or neither") self.p = p self.t = t if validate: self._validate() self.enable_facets = True self._build_mappings() super(MeshTet, self).__init__()
8.2
def _build_mappings(self): """Build element-to-facet, element-to-edges, etc. mappings.""" # define edges: in the order (0,1) (1,2) (0,2) (0,3) (1,3) (2,3) self.edges = np.sort(np.vstack((self.t[0, :], self.t[1, :])), axis=0) e = np.array([1, 2, 0, 2, 0, 3, 1, 3, 2, 3]) for i in range(5): self.edges = np.hstack((self.edges, np.sort(np.vstack((self.t[e[2*i], :], self.t[e[2*i+1], :])), axis=0))) # unique edges self.edges, ixa, ixb = np.unique(self.edges, axis=1, return_index=True, return_inverse=True) self.edges = np.ascontiguousarray(self.edges) self.t2e = ixb.reshape((6, self.t.shape[1])) # define facets if self.enable_facets: self.facets = np.sort(np.vstack((self.t[0, :], self.t[1, :], self.t[2, :])), axis=0) f = np.array([0, 1, 3, 0, 2, 3, 1, 2, 3]) for i in range(3): self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[f[2*i], :], self.t[f[2*i+1], :], self.t[f[2*i+2]])), axis=0))) # unique facets self.facets, ixa, ixb = np.unique(self.facets, axis=1, return_index=True, return_inverse=True) self.facets = np.ascontiguousarray(self.facets) self.t2f = ixb.reshape((4, self.t.shape[1])) # build facet-to-tetra mapping: 2 (tets) x Nfacets e_tmp = np.hstack((self.t2f[0, :], self.t2f[1, :], self.t2f[2, :], self.t2f[3, :])) t_tmp = np.tile(np.arange(self.t.shape[1]), (1, 4))[0] e_first, ix_first = np.unique(e_tmp, return_index=True) # this emulates matlab unique(e_tmp,'last') e_last, ix_last = np.unique(e_tmp[::-1], return_index=True) ix_last = e_tmp.shape[0] - ix_last-1 self.f2t = np.zeros((2, self.facets.shape[1]), dtype=np.int64) self.f2t[0, e_first] = t_tmp[ix_first] self.f2t[1, e_last] = t_tmp[ix_last] # second row to zero if repeated (i.e., on boundary) self.f2t[1, np.nonzero(self.f2t[0, :] == self.f2t[1, :])[0]] = -1
8.3
def refine(self, N=None): """Refine the mesh, tetrahedral optimization. Parameters ---------- N : (optional) int Perform N refinements. """ if N is None: return self._uniform_refine() else: self.enable_facets = False for itr in range(N-1): self._uniform_refine() self.enable_facets = True self._uniform_refine()
8.4
def _uniform_refine(self): """Perform a single mesh refine. Let the nodes of a tetrahedron be numbered as 0, 1, 2 and 3. It is assumed that the edges in self.t2e are given in the order I=(0,1), II=(1,2), III=(0,2), IV=(0,3), V=(1,3), VI=(2,3) by self._build_mappings(). Let I denote the midpoint of the edge (0,1), II denote the midpoint of the edge (1,2), etc. Then each tetrahedron is split into eight smaller subtetrahedra as follows. The first four subtetrahedra have the following nodes: 1. (0,I,III,IV) 2. (1,I,II,V) 3. (2,II,III,VI) 4. (3,IV,V,VI) The remaining middle-portion of the original tetrahedron consists of a union of two mirrored pyramids. This bi-pyramid can be splitted into four tetrahedra in a three different ways by connecting the midpoints of two opposing edges (there are three different pairs of opposite edges). For each tetrahedra in the original mesh, we split the bi-pyramid in such a way that the connection between the opposite edges is shortest. This minimizes the shape-regularity constant of the resulting mesh family. """ # rename variables t = self.t p = self.p e = self.edges sz = p.shape[1] t2e = self.t2e + sz # new vertices are the midpoints of edges newp = 0.5*np.vstack((p[0, e[0, :]] + p[0, e[1, :]], p[1, e[0, :]] + p[1, e[1, :]], p[2, e[0, :]] + p[2, e[1, :]])) newp = np.hstack((p, newp)) # new tets newt = np.vstack((t[0, :], t2e[0, :], t2e[2, :], t2e[3, :])) newt = np.hstack((newt, np.vstack((t[1, :], t2e[0, :], t2e[1, :], t2e[4, :])))) newt = np.hstack((newt, np.vstack((t[2, :], t2e[1, :], t2e[2, :], t2e[5, :])))) newt = np.hstack((newt, np.vstack((t[3, :], t2e[3, :], t2e[4, :], t2e[5, :])))) # compute middle pyramid diagonal lengths and choose shortest d1 = ((newp[0, t2e[2, :]] - newp[0, t2e[4, :]])**2 + (newp[1, t2e[2, :]] - newp[1, t2e[4, :]])**2) d2 = ((newp[0, t2e[1, :]] - newp[0, t2e[3, :]])**2 + (newp[1, t2e[1, :]] - newp[1, t2e[3, :]])**2) d3 = ((newp[0, t2e[0, :]] - newp[0, t2e[5, :]])**2 + (newp[1, t2e[0, :]] - newp[1, t2e[5, :]])**2) I1 = d1 < d2 I2 = d1 < d3 I3 = d2 < d3 c1 = I1*I2 c2 = (~I1)*I3 c3 = (~I2)*(~I3) # splitting the pyramid in the middle. # diagonals are [2,4], [1,3] and [0,5] # CASE 1: diagonal [2,4] newt = np.hstack((newt, np.vstack((t2e[2, c1], t2e[4, c1], t2e[0, c1], t2e[1, c1])))) newt = np.hstack((newt, np.vstack((t2e[2, c1], t2e[4, c1], t2e[0, c1], t2e[3, c1])))) newt = np.hstack((newt, np.vstack((t2e[2, c1], t2e[4, c1], t2e[1, c1], t2e[5, c1])))) newt = np.hstack((newt, np.vstack((t2e[2, c1], t2e[4, c1], t2e[3, c1], t2e[5, c1])))) # CASE 2: diagonal [1,3] newt = np.hstack((newt, np.vstack((t2e[1, c2], t2e[3, c2], t2e[0, c2], t2e[4, c2])))) newt = np.hstack((newt, np.vstack((t2e[1, c2], t2e[3, c2], t2e[4, c2], t2e[5, c2])))) newt = np.hstack((newt, np.vstack((t2e[1, c2], t2e[3, c2], t2e[5, c2], t2e[2, c2])))) newt = np.hstack((newt, np.vstack((t2e[1, c2], t2e[3, c2], t2e[2, c2], t2e[0, c2])))) # CASE 3: diagonal [0,5] newt = np.hstack((newt, np.vstack((t2e[0, c3], t2e[5, c3], t2e[1, c3], t2e[4, c3])))) newt = np.hstack((newt, np.vstack((t2e[0, c3], t2e[5, c3], t2e[4, c3], t2e[3, c3])))) newt = np.hstack((newt, np.vstack((t2e[0, c3], t2e[5, c3], t2e[3, c3], t2e[2, c3])))) newt = np.hstack((newt, np.vstack((t2e[0, c3], t2e[5, c3], t2e[2, c3], t2e[1, c3])))) # update fields self.p = newp self.t = newt self._build_mappings()
8.5
# TODO implement prolongation matrix def shapereg(self): """Return the largest shape-regularity constant.""" def edgelen(n): return np.sqrt(np.sum((self.p[:, self.edges[0, self.t2e[n, :]]] - self.p[:, self.edges[1, self.t2e[n, :]]])**2, axis=0)) edgelenmat = np.vstack(tuple(edgelen(i) for i in range(6))) return np.max(np.max(edgelenmat, axis=0)/np.min(edgelenmat, axis=0))
8.6
def mapping(self): return skfem.mapping.MappingAffine(self)
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
三角剖分
class MeshTri(Mesh2D): """A mesh consisting of triangular elements. Attributes ---------- p : numpy array of size 2 x Nvertices The vertices of the mesh t : numpy array of size 3 x Nelements The element connectivity facets : numpy array of size 2 x Nfacets Each column contains a pair of indices to p. f2t : numpy array of size 2 x Nfacets Each column contains a pair of indices to t or -1 on the second row if the facet is on the boundary. t2f : numpy array of size 3 x Nelements Each column contains three indices to facets. Examples -------- Initialize a symmetric mesh of the unit square. >>> m = MeshTri.init_sqsymmetric() >>> m.t.shape (3, 8) The different constructors are: * load (requires meshio) * init_symmetric * init_sqsymmetric * init_refdom Facets (edges) and mappings from triangles to facets and vice versa are automatically constructed. In the following example we have 5 facets (edges). >>> m = MeshTri() >>> m.facets array([[0, 0, 1, 1, 2], [1, 2, 2, 3, 3]]) >>> m.t2f array([[0, 2], [2, 4], [1, 3]]) >>> m.f2t array([[ 0, 0, 1, 1, 1], [-1, -1, 0, -1, -1]]) The value -1 implies that the facet (the edge) is on the boundary. Refine the triangular mesh of the unit square three times. >>> m = MeshTri() >>> m.refine(3) >>> m.p.shape (2, 81) """ refdom = "tri" brefdom = "line" meshio_type = "triangle" p = np.array([]) t = np.array([]) facets = np.array([]) f2t = np.array([]) t2f = np.array([])
9.1
def __init__(self, p=None, t=None, validate=True, sort_t=True): """Initialize a triangular mesh. Parameters ---------- p : (optional) numpy array of size 2 x Nvertices The points of the mesh. t : (optional) numpy array of size 3 x Nelements The element connectivity, i.e. indices to p. validate : (optional, default=True) bool Whether to run mesh validity checks or not. sort_t : (optional, default=True) bool Sort the element connectivity matrix before building mappings. """ if p is None and t is None: p = np.array([[0., 1., 0., 1.], [0., 0., 1., 1.]], dtype=np.float_) t = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.intp).T elif p is None or t is None: raise Exception("Must provide p AND t or neither") self.p = p self.t = t if validate: self._validate() self._build_mappings(sort_t=sort_t) super(MeshTri, self).__init__()
9.2
@classmethod def init_tensor(cls, x, y): """Initialize a tensor product mesh. Parameters ---------- x : numpy array (1d) The nodal coordinates in dimension x y : numpy array (1d) The nodal coordinates in dimension y """ return MeshQuad.init_tensor(x, y)._splitquads()
9.3
@classmethod def init_symmetric(cls): """Initialize a symmetric mesh of the unit square. The mesh topology is as follows: *------------* | /| | / | | / | | * | | / | | / | |/ | *------------* """ p = np.array([[0, 1, 1, 0, 0.5], [0, 0, 1, 1, 0.5]], dtype=np.float_) t = np.array([[0, 1, 4], [1, 2, 4], [2, 3, 4], [0, 3, 4]], dtype=np.intp).T return cls(p, t)
9.4
@classmethod def init_sqsymmetric(cls): """Initialize a symmetric mesh of the unit square. The mesh topology is as follows: *------------* | | /| | | / | | | / | |-----*------| | /| | | / | | |/ | | *------------* """ p = np.array([[0, 0.5, 1, 0, 0.5, 1, 0, 0.5, 1], [0, 0, 0, 0.5, 0.5, 0.5, 1, 1, 1]], dtype=np.float_) t = np.array([[0, 1, 4], [1, 2, 4], [2, 4, 5], [0, 3, 4], [3, 4, 6], [4, 6, 7], [4, 7, 8], [4, 5, 8]], dtype=np.intp).T return cls(p, t)
9.5
@classmethod def init_refdom(cls): """Initialize a mesh that includes only the reference triangle. The mesh topology is as follows: * | | | | | | | *-------------* """ p = np.array([[0., 1., 0.], [0., 0., 1.]], dtype=np.float_) t = np.array([[0, 1, 2]], dtype=np.intp).T return cls(p, t)
9.6
def _build_mappings(self, sort_t=True): # sort to preserve orientations etc. if sort_t: self.t = np.sort(self.t, axis=0) # define facets: in the order (0,1) (1,2) (0,2) self.facets = np.sort(np.vstack((self.t[0, :], self.t[1, :])), axis=0) self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[1, :], self.t[2, :])), axis=0))) self.facets = np.hstack((self.facets, np.sort(np.vstack((self.t[0, :], self.t[2, :])), axis=0))) # get unique facets and build triangle-to-facet # mapping: 3 (edges) x Ntris tmp = np.ascontiguousarray(self.facets.T) tmp, ixa, ixb = np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1]), return_index=True, return_inverse=True) self.facets = self.facets[:, ixa] self.t2f = ixb.reshape((3, self.t.shape[1])) # build facet-to-triangle mapping: 2 (triangles) x Nedges e_tmp = np.hstack((self.t2f[0, :], self.t2f[1, :], self.t2f[2, :])) t_tmp = np.tile(np.arange(self.t.shape[1]), (1, 3))[0] e_first, ix_first = np.unique(e_tmp, return_index=True) # this emulates matlab unique(e_tmp,'last') e_last, ix_last = np.unique(e_tmp[::-1], return_index=True) ix_last = e_tmp.shape[0] - ix_last - 1 self.f2t = np.zeros((2, self.facets.shape[1]), dtype=np.int64) self.f2t[0, e_first] = t_tmp[ix_first] self.f2t[1, e_last] = t_tmp[ix_last] # second row to zero if repeated (i.e., on boundary) self.f2t[1, np.nonzero(self.f2t[0, :] == self.f2t[1, :])[0]] = -1
9.7
def interpolator(self, x): """Return a function which interpolates values with P1 basis.""" triang = mtri.Triangulation(self.p[0, :], self.p[1, :], self.t.T) interpf = mtri.LinearTriInterpolator(triang, x) # contruct an interpolator handle def handle(X, Y): return interpf(X, Y).data return handle
9.8
def const_interpolator(self, x): """Return a function which interpolates values with P0 basis.""" triang = mtri.Triangulation(self.p[0, :], self.p[1, :], self.t.T) finder = triang.get_trifinder() # construct an interpolator handle def handle(X, Y): return x[finder(X, Y)] return handle
9.9
def smooth(self, k=1.0): """Apply smoothing to interior nodes.""" from skfem.assembly import InteriorBasis, asm from skfem.element import ElementTriP1 from skfem.utils import solve from skfem.models.poisson import laplace, mass e = ElementTriP1() ib = InteriorBasis(self, e) K = asm(laplace, ib) M = asm(mass, ib) I = self.interior_nodes() dx = - k*solve(M, K.dot(self.p[0, :])) dy = - k*solve(M, K.dot(self.p[1, :])) self.p[0, I] += dx[I] self.p[1, I] += dy[I]
9.10
def draw_debug(self): """Draw without mesh.facets. For debugging self.draw().""" fig = plt.figure() plt.hold('on') for itr in range(self.t.shape[1]): plt.plot(self.p[0,self.t[[0,1],itr]], self.p[1,self.t[[0,1],itr]], 'k-') plt.plot(self.p[0,self.t[[1,2],itr]], self.p[1,self.t[[1,2],itr]], 'k-') plt.plot(self.p[0,self.t[[0,2],itr]], self.p[1,self.t[[0,2],itr]], 'k-') return fig
9.11
def plot(self, z, smooth=False, ax=None, zlim=None, edgecolors=None): """Visualize nodal or elemental function (2d).""" if ax is None: fig = plt.figure() ax = fig.add_subplot(111) if edgecolors is None: edgecolors = 'k' if zlim == None: if smooth: ax.tripcolor(self.p[0, :], self.p[1, :], self.t.T, z, shading='gouraud', edgecolors=edgecolors) else: ax.tripcolor(self.p[0, :], self.p[1, :], self.t.T, z, edgecolors=edgecolors) else: if smooth: ax.tripcolor(self.p[0, :], self.p[1, :], self.t.T, z, shading='gouraud', vmin=zlim[0], vmax=zlim[1], edgecolors=edgecolors) else: ax.tripcolor(self.p[0, :], self.p[1, :], self.t.T, z, vmin=zlim[0], vmax=zlim[1], edgecolors=edgecolors) return ax
9.12
def plot3(self, z, smooth=False, ax=None): """Visualize nodal function (3d i.e. three axes).""" from mpl_toolkits.mplot3d import Axes3D if ax is None: fig = plt.figure() ax = Axes3D(fig) if len(z) == self.p.shape[1]: # use matplotlib ts = mtri.Triangulation(self.p[0, :], self.p[1, :], self.t.T) ax.plot_trisurf(self.p[0, :], self.p[1, :], z, triangles=ts.triangles, cmap=plt.cm.Spectral) elif len(z) == self.t.shape[1]: # one value per element (piecewise const) nt = self.t.shape[1] newt = np.arange(3*nt, dtype=np.int64).reshape((nt, 3)) newpx = self.p[0, self.t].flatten(order='F') newpy = self.p[1, self.t].flatten(order='F') newz = np.vstack((z, z, z)).flatten(order='F') ts = mtri.Triangulation(newpx, newpx, newt) ax.plot_trisurf(newpx, newpy, newz, triangles=ts.triangles, cmap=plt.cm.Spectral) elif len(z) == 3*self.t.shape[1]: # three values per element (piecewise linear) nt = self.t.shape[1] newt = np.arange(3*nt, dtype=np.int64).reshape((nt, 3)) newpx = self.p[0, self.t].flatten(order='F') newpy = self.p[1, self.t].flatten(order='F') ts = mtri.Triangulation(newpx, newpx, newt) ax.plot_trisurf(newpx, newpy, z, triangles=ts.triangles, cmap=plt.cm.Spectral) else: raise NotImplementedError("MeshTri.plot3: not implemented for " "the given shape of input vector!") return ax
9.13
def _uniform_refine(self): """Perform a single mesh refine.""" # rename variables t = self.t p = self.p e = self.facets sz = p.shape[1] t2f = self.t2f + sz # new vertices are the midpoints of edges newp = 0.5*np.vstack((p[0, e[0, :]] + p[0, e[1, :]], p[1, e[0, :]] + p[1, e[1, :]])) newp = np.hstack((p, newp)) # build new triangle definitions newt = np.vstack((t[0, :], t2f[0, :], t2f[2, :])) newt = np.hstack((newt, np.vstack((t[1, :], t2f[0, :], t2f[1, :])))) newt = np.hstack((newt, np.vstack((t[2, :], t2f[2, :], t2f[1, :])))) newt = np.hstack((newt, np.vstack((t2f[0, :], t2f[1, :], t2f[2, :])))) # update fields self.p = newp self.t = newt self._build_mappings() # prolongation matrix nsz = newp.shape[1] return coo_matrix( (np.hstack((np.ones(sz), .5 * np.ones(2 * (nsz - sz)))), (np.hstack((np.arange(sz), np.arange(nsz - sz) + sz, np.arange(nsz - sz) + sz)), np.hstack((np.arange(sz), e[0, :], e[1, :])))), shape=(nsz, sz)).tocsr()
9.14
def _uniform_refine(self): """Perform a single mesh refine.""" # rename variables t = self.t p = self.p e = self.facets sz = p.shape[1] t2f = self.t2f + sz # new vertices are the midpoints of edges newp = 0.5*np.vstack((p[0, e[0, :]] + p[0, e[1, :]], p[1, e[0, :]] + p[1, e[1, :]])) newp = np.hstack((p, newp)) # build new triangle definitions newt = np.vstack((t[0, :], t2f[0, :], t2f[2, :])) newt = np.hstack((newt, np.vstack((t[1, :], t2f[0, :], t2f[1, :])))) newt = np.hstack((newt, np.vstack((t[2, :], t2f[2, :], t2f[1, :])))) newt = np.hstack((newt, np.vstack((t2f[0, :], t2f[1, :], t2f[2, :])))) # update fields self.p = newp self.t = newt self._build_mappings() # prolongation matrix nsz = newp.shape[1] return coo_matrix( (np.hstack((np.ones(sz), .5 * np.ones(2 * (nsz - sz)))), (np.hstack((np.arange(sz), np.arange(nsz - sz) + sz, np.arange(nsz - sz) + sz)), np.hstack((np.arange(sz), e[0, :], e[1, :])))), shape=(nsz, sz)).tocsr()
9.15
def mapping(self): return skfem.mapping.MappingAffine(self)
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
if __name__ == "__main__": import doctest doctest.testmod()