From 7b998e43632c09b9b7438223c18995f4fcab744d Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Tue, 7 Jan 2020 17:04:08 +0200 Subject: [PATCH 01/32] Add MeshMortar and fix ex04 --- docs/examples/ex04.py | 49 ++++----- skfem/__init__.py | 1 + skfem/assembly/basis/basis.py | 25 +++-- skfem/assembly/basis/mortar_basis.py | 62 ++++++++---- skfem/mesh/__init__.py | 1 + skfem/mesh/mesh_mortar.py | 144 +++++++++++++++++++++++++++ 6 files changed, 223 insertions(+), 59 deletions(-) create mode 100644 skfem/mesh/mesh_mortar.py diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index 00045dd5c..a799fbd8c 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -1,41 +1,36 @@ -""" -Author: kinnala - -Solve a linearized contact problem between two linear elastic -bodies using the Nitsche's method. -""" - from skfem import * from skfem.models.elasticity import linear_elasticity +from skfem.visuals.matplotlib import * import numpy as np +# create meshes m = MeshTri() m.refine(3) M = MeshLine(np.linspace(0, 1, 6)) M = M*M -M = M._splitquads() M.translate((1.0, 0.0)) -map = MappingAffine(m) -Map = MappingAffine(M) e1 = ElementTriP1() e = ElementVectorH1(e1) -ib = InteriorBasis(m, e, map, 2) -Ib = InteriorBasis(M, e, Map, 2) +E1 = ElementQuad1() +E = ElementVectorH1(E1) + +ib = InteriorBasis(m, e) +Ib = InteriorBasis(M, E) -def rule(x, y): - return (x==1.0) +def rule(x): + return (x[0]==1.0) def param(x, y): return y -mortar = InterfaceMesh1D(m, M, rule, param) +mortar = MeshMortar(m, M, rule, param) -mb = {} -mortar_map = MappingAffine(mortar) -mb[0] = FacetBasis(mortar, e, mortar_map, 2, side=0) -mb[1] = FacetBasis(mortar, e, mortar_map, 2, side=1) +mb = [ + MortarBasis(mortar, e, intorder=2, side=0), + MortarBasis(mortar, E, intorder=2, side=1) +] E1 = 1000.0 E2 = 1000.0 @@ -52,13 +47,13 @@ def param(x, y): Mu = Mu1 Lambda = Lambda1 -weakform1 = linear_elasticity(Lambda=Lambda1, Mu=Mu1) -weakform2 = linear_elasticity(Lambda=Lambda2, Mu=Mu2) +weakform1 = linear_elasticity(Lambda=Lambda1, Mu=Mu) +weakform2 = linear_elasticity(Lambda=Lambda2, Mu=Mu) alpha = 1 K1 = asm(weakform1, ib) K2 = asm(weakform2, Ib) -L = 0 +L = [[None,None],[None,None]] for i in range(2): for j in range(2): @bilinear_form @@ -88,7 +83,7 @@ def Eps(dw): h = w.h return 1.0/(alpha*h)*ju*jv - mu*jv - mv*ju - L = asm(bilin_penalty, mb[i], mb[j]) + L + L[i][j] = asm(bilin_penalty, mb[i], mb[j]) @linear_form def load(v, dv, w): @@ -98,13 +93,13 @@ def load(v, dv, w): f2 = np.zeros(K2.shape[0]) import scipy.sparse -K = (scipy.sparse.bmat([[K1, None],[None, K2]]) + L).tocsr() +K = (scipy.sparse.bmat([[K1 + L[0][0], L[1][0]],[L[0][1], K2 + L[1][1]]])).tocsr() i1 = np.arange(K1.shape[0]) i2 = np.arange(K2.shape[0]) + K1.shape[0] -D1 = ib.get_dofs(lambda x, y: x==0.0).all() -D2 = Ib.get_dofs(lambda x, y: x==2.0).all() +D1 = ib.get_dofs(lambda x: x[0]==0.0).all() +D2 = Ib.get_dofs(lambda x: x[0]==2.0).all() x = np.zeros(K.shape[0]) @@ -114,7 +109,7 @@ def load(v, dv, w): D = np.concatenate((D1, D2 + ib.N)) I = np.setdiff1d(np.arange(K.shape[0]), D) -x[I] = solve(*condense(K, f, I=I)) +x = solve(*condense(K, f, I=I)) sf = 1 diff --git a/skfem/__init__.py b/skfem/__init__.py index 3e85b0db3..c3d94aa54 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -12,6 +12,7 @@ 'MeshQuad', 'MeshTet', 'MeshTri', + 'MeshMortar', 'FacetBasis', 'Basis', 'InteriorBasis', diff --git a/skfem/assembly/basis/basis.py b/skfem/assembly/basis/basis.py index 82f412993..85ba71bb4 100644 --- a/skfem/assembly/basis/basis.py +++ b/skfem/assembly/basis/basis.py @@ -1,3 +1,4 @@ +import warnings from typing import List, Optional, NamedTuple, Any import numpy as np @@ -28,6 +29,19 @@ def __init__(self, mesh, elem, mapping, intorder): self.mapping = mapping self._build_dofnum(mesh, elem) + + # human readable names + self.dofnames = elem.dofnames + + # global degree-of-freedom location + if hasattr(elem, 'doflocs'): + doflocs = self.mapping.F(elem.doflocs.T) + self.doflocs = np.zeros((doflocs.shape[0], self.N)) + for itr in range(doflocs.shape[0]): + for jtr in range(self.element_dofs.shape[0]): + self.doflocs[itr, self.element_dofs[jtr]] =\ + doflocs[itr, :, jtr] + self.mesh = mesh self.elem = elem @@ -104,17 +118,8 @@ def _build_dofnum(self, mesh, element): # interior dofs self.element_dofs = np.vstack((self.element_dofs, self.interior_dofs)) - # number-of-dofs and human readable names + # total dofs self.N = np.max(self.element_dofs) + 1 - self.dofnames = element.dofnames - - if hasattr(element, 'doflocs'): - doflocs = self.mapping.F(element.doflocs.T) - self.doflocs = np.zeros((doflocs.shape[0], self.N)) - for itr in range(doflocs.shape[0]): - for jtr in range(self.element_dofs.shape[0]): - self.doflocs[itr, self.element_dofs[jtr]] =\ - doflocs[itr, :, jtr] def complement_dofs(self, *D): if type(D[0]) is dict: diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index 9164e1369..6dff2997a 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -1,5 +1,3 @@ -from typing import Optional - import numpy as np from numpy import ndarray @@ -11,49 +9,69 @@ class MortarBasis(Basis): """Global basis functions evaluated at integration points on the mortar boundary. """ + def __init__(self, mesh, elem, - mapping, - intorder: Optional[int] = None, + intorder: int = None, side: int = 0): - super(MortarBasis, self).__init__(mesh, elem, mapping, intorder) - - self.ib1 = InteriorBasis(mesh.mesh1, elem) - self.ib2 = InteriorBasis(mesh.mesh2, elem) + """Initialize a basis for assembling mortar matrices. + + Parameters + ---------- + mesh + An object of type :class:`~skfem.mesh.MeshMortar`. + elem + An object of type :class:`~skfem.element.Element`. + intorder + Required integration order, i.e. the degree of polynomials that are + integrated exactly by the used quadrature. + side + The numbers 0 and 1 refer to the different sides of the mortar + interface. Side 0 corresponds to the indices mesh.f2t[0, :]. + + """ + if intorder is None: + raise ValueError("Please specify 'intorder' for MortarBasis " + "and use consistent orders for both sides.") + + # build mapping for the mortar mesh + super(MortarBasis, self).__init__(mesh, elem, None, intorder) + + # build dofnum and mappings for the target mesh + self._build_dofnum(mesh.target_mesh[side], elem) + target_mapping = mesh.target_mesh[side].mapping() self.X, self.W = get_quadrature(self.brefdom, self.intorder) - self.find = np.nonzero(self.mesh.f2t[1, :] != -1)[0] + self.find = np.arange(self.mesh.facets.shape[1]) self.tind = self.mesh.f2t[side, self.find] # boundary refdom to global facet x = self.mapping.G(self.X, find=self.find) # global facet to refdom facet - Y = self.mapping.invF(x, tind=self.tind) + Y = target_mapping.invF(x, tind=self.tind) - self.normals = np.repeat(mesh.normals[:, :, None], len(self.W), axis=2) + # normals are defined in the mortar mesh + self.normals = np.repeat(mesh.normals[:, :, None], + len(self.W), + axis=2) - self.nelems = len(self.find) - - self.basis = [self.elem.gbasis(self.mapping, Y, j, self.tind) + self.basis = [self.elem.gbasis(target_mapping, Y, j, self.tind) for j in range(self.Nbfun)] + self.nelems = len(self.find) self.dx = np.abs(self.mapping.detDG(self.X, find=self.find)) *\ np.tile(self.W, (self.nelems, 1)) - if side == 0: - self.element_dofs = self.ib1.element_dofs[:, self.tind] - elif side == 1: - self.element_dofs = self.ib2.element_dofs[:, self.tind - mesh.mesh1.t.shape[1]] + self.ib1.N + self.element_dofs = self.element_dofs[:, self.tind] - self.N = self.ib1.N + self.ib2.N def default_parameters(self): """Return default parameters for `~skfem.assembly.asm`.""" - return {'x':self.global_coordinates(), - 'h':self.mesh_parameters(), - 'n':self.normals} + return {'x': self.global_coordinates(), + 'h': self.mesh_parameters(), + 'n': self.normals} def global_coordinates(self) -> ndarray: return self.mapping.G(self.X, find=self.find) diff --git a/skfem/mesh/__init__.py b/skfem/mesh/__init__.py index 8777fd592..1f85ee0fd 100644 --- a/skfem/mesh/__init__.py +++ b/skfem/mesh/__init__.py @@ -17,3 +17,4 @@ from .mesh_line import MeshLine from .mesh2d import Mesh2D, MeshTri, MeshQuad from .mesh3d import Mesh3D, MeshTet, MeshHex +from .mesh_mortar import MeshMortar diff --git a/skfem/mesh/mesh_mortar.py b/skfem/mesh/mesh_mortar.py new file mode 100644 index 000000000..74cce0a15 --- /dev/null +++ b/skfem/mesh/mesh_mortar.py @@ -0,0 +1,144 @@ +import numpy as np + +from .mesh import Mesh +from ..mapping import MappingAffine + + +class MeshMortar(Mesh): + """An interface mesh for mortar methods.""" + + 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.array([[0, 1, 2]]).T + self.facets = np.array([ixorig[:-1], ixorig[1:]]) + self.t2f = -1 + 0 * self.t + + self.target_mesh = [None,None] + self.target_mesh[0] = mesh1 + self.target_mesh[1] = mesh2 + + # 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 + self.f2f = 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.f2f[0, itr] = jtr + self.f2t[0, itr] = mesh1.f2t[0, jtr] + break + elif val < param(x1, y1) and val > param(x2, y2): # ye olde + # OK + self.f2f[0, itr] = jtr + self.f2t[0, itr] = mesh1.f2t[0, jtr] + break + elif val >= param(x1, y1) and val < param(x2, y2): + # OK + self.f2f[0, itr] = jtr + self.f2t[0, itr] = mesh1.f2t[0, jtr] + break + elif val > param(x1, y1) and val <= param(x2, y2): + # OK + self.f2f[0, itr] = jtr + self.f2t[0, itr] = mesh1.f2t[0, jtr] + break + elif val <= param(x1, y1) and val > param(x2, y2): + # OK + self.f2f[0, itr] = jtr + self.f2t[0, itr] = mesh1.f2t[0, jtr] + break + elif val < param(x1, y1) and val >= param(x2, y2): + # OK + self.f2f[0, itr] = jtr + 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.f2f[1, itr] = jtr + self.f2t[1, itr] = mesh2.f2t[0, jtr] + break + elif val < param(x1, y1) and val > param(x2, y2): + # OK + self.f2f[1, itr] = jtr + self.f2t[1, itr] = mesh2.f2t[0, jtr] + break + elif val >= param(x1, y1) and val < param(x2, y2): + # OK + self.f2f[1, itr] = jtr + self.f2t[1, itr] = mesh2.f2t[0, jtr] + break + elif val > param(x1, y1) and val <= param(x2, y2): + # OK + self.f2f[1, itr] = jtr + self.f2t[1, itr] = mesh2.f2t[0, jtr] + break + elif val <= param(x1, y1) and val > param(x2, y2): + # OK + self.f2f[1, itr] = jtr + self.f2t[1, itr] = mesh2.f2t[0, jtr] + break + elif val < param(x1, y1) and val >= param(x2, y2): + # OK + self.f2f[1, itr] = jtr + self.f2t[1, itr] = mesh2.f2t[0, jtr] + 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!") + + def mapping(self): + return MappingAffine(self) From e5da7ce5fb18c9642a391b1e1ad7f96e4efcfe12 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Tue, 7 Jan 2020 17:04:19 +0200 Subject: [PATCH 02/32] Add doflocs to ElementHex1 --- skfem/element/element_hex/element_hex1.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/skfem/element/element_hex/element_hex1.py b/skfem/element/element_hex/element_hex1.py index af2d0f79a..874ef14e6 100644 --- a/skfem/element/element_hex/element_hex1.py +++ b/skfem/element/element_hex/element_hex1.py @@ -7,6 +7,14 @@ class ElementHex1(ElementH1): dim = 3 maxdeg = 3 dofnames = ['u'] + doflocs = np.array([[-1., -1., -1.], + [-1., -1., 1.], + [-1., 1., -1.], + [ 1., -1., -1.], + [-1., 1., 1.], + [ 1., -1., 1.], + [ 1., 1., -1.], + [ 1., 1., 1.]]) def lbasis(self, X, i): x, y, z = X From 326bc447f18ac18a08807e27475edce4bef7ed67 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Tue, 7 Jan 2020 17:09:23 +0200 Subject: [PATCH 03/32] Clean up --- skfem/mesh/mesh_mortar.py | 44 ++++----------------------------------- 1 file changed, 4 insertions(+), 40 deletions(-) diff --git a/skfem/mesh/mesh_mortar.py b/skfem/mesh/mesh_mortar.py index 74cce0a15..764820379 100644 --- a/skfem/mesh/mesh_mortar.py +++ b/skfem/mesh/mesh_mortar.py @@ -26,7 +26,7 @@ def __init__(self, mesh1, mesh2, rule, param, debug_plot=False): self.facets = np.array([ixorig[:-1], ixorig[1:]]) self.t2f = -1 + 0 * self.t - self.target_mesh = [None,None] + self.target_mesh = 2 * [None] self.target_mesh[0] = mesh1 self.target_mesh[1] = mesh2 @@ -37,21 +37,9 @@ def __init__(self, mesh1, mesh2, rule, param, debug_plot=False): 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 - self.f2f = self.facets*0-1 + 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]]) @@ -65,33 +53,21 @@ def __init__(self, mesh1, mesh2, rule, param, debug_plot=False): 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.f2f[0, itr] = jtr self.f2t[0, itr] = mesh1.f2t[0, jtr] break - elif val < param(x1, y1) and val > param(x2, y2): # ye olde - # OK - self.f2f[0, itr] = jtr + elif val < param(x1, y1) and val > param(x2, y2): self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val >= param(x1, y1) and val < param(x2, y2): - # OK - self.f2f[0, itr] = jtr self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val > param(x1, y1) and val <= param(x2, y2): - # OK - self.f2f[0, itr] = jtr self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val <= param(x1, y1) and val > param(x2, y2): - # OK - self.f2f[0, itr] = jtr self.f2t[0, itr] = mesh1.f2t[0, jtr] break elif val < param(x1, y1) and val >= param(x2, y2): - # OK - self.f2f[0, itr] = jtr self.f2t[0, itr] = mesh1.f2t[0, jtr] break for jtr in mesh2.boundary_facets(): @@ -103,37 +79,25 @@ def __init__(self, mesh1, mesh2, rule, param, debug_plot=False): 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.f2f[1, itr] = jtr self.f2t[1, itr] = mesh2.f2t[0, jtr] break elif val < param(x1, y1) and val > param(x2, y2): - # OK - self.f2f[1, itr] = jtr self.f2t[1, itr] = mesh2.f2t[0, jtr] break elif val >= param(x1, y1) and val < param(x2, y2): - # OK - self.f2f[1, itr] = jtr self.f2t[1, itr] = mesh2.f2t[0, jtr] break elif val > param(x1, y1) and val <= param(x2, y2): - # OK - self.f2f[1, itr] = jtr self.f2t[1, itr] = mesh2.f2t[0, jtr] break elif val <= param(x1, y1) and val > param(x2, y2): - # OK - self.f2f[1, itr] = jtr self.f2t[1, itr] = mesh2.f2t[0, jtr] break elif val < param(x1, y1) and val >= param(x2, y2): - # OK - self.f2f[1, itr] = jtr self.f2t[1, itr] = mesh2.f2t[0, jtr] break - if (self.f2t>-1).all(): + if (self.f2t > -1).all(): self.f2t[0, :] return else: From 77b7bed66a6ff37376bf7c4ddd939247c2bd0473 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 10 Jan 2020 15:00:09 +0200 Subject: [PATCH 04/32] Support more general directions and cases in MeshMortar --- docs/examples/ex04.py | 5 +- skfem/mesh/mesh_mortar.py | 177 ++++++++++++++++++++------------------ tests/test_mesh.py | 42 +++++++++ 3 files changed, 139 insertions(+), 85 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index a799fbd8c..c714b9a33 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -25,7 +25,10 @@ def rule(x): def param(x, y): return y -mortar = MeshMortar(m, M, rule, param) +mortar = MeshMortar.init_1D(m, M, + m.facets_satisfying(lambda x: x[0]==1.0), + M.facets_satisfying(lambda x: x[0]==1.0), + np.array([0.0, 1.0])) mb = [ MortarBasis(mortar, e, intorder=2, side=0), diff --git a/skfem/mesh/mesh_mortar.py b/skfem/mesh/mesh_mortar.py index 764820379..18c5ba085 100644 --- a/skfem/mesh/mesh_mortar.py +++ b/skfem/mesh/mesh_mortar.py @@ -7,102 +7,111 @@ class MeshMortar(Mesh): """An interface mesh for mortar methods.""" - def __init__(self, mesh1, mesh2, rule, param, debug_plot=False): + name = "Mortar" + + def __init__(self, mesh1, mesh2, p, facets, f2t, normals): + """Create an interface mesh for mortar methods. + + Parameters + ---------- + mesh1 + mesh2 + p + facets + f2t + normals + + """ + self.p = p + self.facets = facets + self.f2t = f2t + self.normals = normals + self.target_mesh = 2 * [None] + self.target_mesh[0] = mesh1 + self.target_mesh[1] = mesh2 + + # for quadrature rules self.brefdom = mesh1.brefdom - p1_ix = mesh1.nodes_satisfying(rule) - p2_ix = mesh2.nodes_satisfying(rule) + # dummy values so that initialising Mapping doesn't crash + self.t = mesh1.t[:, :1] + self.t2f = -1 + 0 * self.t + + @classmethod + def init_1D(cls, mesh1, mesh2, boundary1, boundary2, tangent): + """Create a mortar mesh between two 2D meshes via projection. + + Parameters + ---------- + mesh1 + mesh2 + boundary1 + A subset of facets from mesh1. + boundary2 + A subset of facets from mesh2. + tangent + A tangent vector defining the direction of the projection. + + """ + tangent /= np.linalg.norm(tangent) + + # find nodes on the two boundaries + p1_ix = np.unique(mesh1.facets[:, boundary1].flatten()) + p2_ix = np.unique(mesh2.facets[:, boundary2].flatten()) 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] + def proj(p): + """Project onto the line defined by 'tangent'.""" + return np.outer(tangent, tangent) @ p - self.p = np.hstack((mesh1.p, mesh2.p)) - self.t = np.array([[0, 1, 2]]).T - self.facets = np.array([ixorig[:-1], ixorig[1:]]) - self.t2f = -1 + 0 * self.t + def param(p): + """Calculate signed distances of projected points from origin.""" + y = proj(p) + return np.linalg.norm(y, axis=0) * np.dot(tangent, y) - self.target_mesh = 2 * [None] - self.target_mesh[0] = mesh1 - self.target_mesh[1] = mesh2 + # find unique facets by combining nodes from both sides + _, ix = np.unique(np.concatenate((param(p1), param(p2))), return_index=True) + ixorig = np.concatenate((p1_ix, p2_ix + mesh1.p.shape[1]))[ix] + p = np.hstack((mesh1.p, mesh2.p)) # TODO has duplicate nodes + facets = np.array([ixorig[:-1], ixorig[1:]]) # 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]) + normal = np.array([tangent[1], -tangent[0]]) + normals = normal[:, None].repeat(facets.shape[1], axis=1) # 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): - self.f2t[0, itr] = mesh1.f2t[0, jtr] - break - elif val < param(x1, y1) and val > param(x2, y2): - self.f2t[0, itr] = mesh1.f2t[0, jtr] - break - elif val >= param(x1, y1) and val < param(x2, y2): - self.f2t[0, itr] = mesh1.f2t[0, jtr] - break - elif val > param(x1, y1) and val <= param(x2, y2): - self.f2t[0, itr] = mesh1.f2t[0, jtr] - break - elif val <= param(x1, y1) and val > param(x2, y2): - self.f2t[0, itr] = mesh1.f2t[0, jtr] - break - elif val < param(x1, y1) and val >= param(x2, y2): - 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): - self.f2t[1, itr] = mesh2.f2t[0, jtr] - break - elif val < param(x1, y1) and val > param(x2, y2): - self.f2t[1, itr] = mesh2.f2t[0, jtr] - break - elif val >= param(x1, y1) and val < param(x2, y2): - self.f2t[1, itr] = mesh2.f2t[0, jtr] - break - elif val > param(x1, y1) and val <= param(x2, y2): - self.f2t[1, itr] = mesh2.f2t[0, jtr] - break - elif val <= param(x1, y1) and val > param(x2, y2): - self.f2t[1, itr] = mesh2.f2t[0, jtr] - break - elif val < param(x1, y1) and val >= param(x2, y2): - self.f2t[1, itr] = mesh2.f2t[0, jtr] - 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!") + f2t = facets * 0 - 1 + for itr in range(facets.shape[1]): + mp = .5 * (p[:, facets[0, itr]] + p[:, facets[1, itr]]) + val = param(mp) + for jtr in boundary1: + x1 = mesh1.p[:, mesh1.facets[0, jtr]] + x2 = mesh1.p[:, mesh1.facets[1, jtr]] + if (val > param(x1) and val < param(x2) or + val < param(x1) and val > param(x2) or + val >= param(x1) and val < param(x2) or + val > param(x1) and val <= param(x2) or + val <= param(x1) and val > param(x2) or + val < param(x1) and val >= param(x2)): + f2t[0, itr] = mesh1.f2t[0, jtr] + break + for jtr in boundary2: + x1 = mesh2.p[:, mesh2.facets[0, jtr]] + x2 = mesh2.p[:, mesh2.facets[1, jtr]] + if (val > param(x1) and val < param(x2) or + val < param(x1) and val > param(x2) or + val >= param(x1) and val < param(x2) or + val > param(x1) and val <= param(x2) or + val <= param(x1) and val > param(x2) or + val < param(x1) and val >= param(x2)): + f2t[1, itr] = mesh2.f2t[0, jtr] + break + if not (f2t > -1).all(): + raise Exception("All mesh facets corresponding to mortar facets not found! {}".format(f2t)) + + return cls(mesh1, mesh2, p, facets, f2t, normals) def mapping(self): return MappingAffine(self) diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 9f3753e14..afa950c8f 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -129,5 +129,47 @@ def runTest(self): self.assertTrue((m.subdomains[k] == M.subdomains[k]).all()) +class TestMortarMeshTriQuad(unittest.TestCase): + + def create_mesh1(self): + m1 = MeshTri() + m1.refine(3) + m1.translate((0.0, -0.5)) + return m1 + + def create_mesh2(self): + m2 = MeshQuad() + m2.refine(2) + m2.translate((1.0, -0.5)) + return m2 + + def create_meshes(self): + m1 = self.create_mesh1() + m2 = self.create_mesh2() + mortar = MeshMortar.init_1D(m1, m2, + m1.facets_satisfying(lambda x: x[0]==1.0), + m2.facets_satisfying(lambda x: x[0]==1.0), + np.array([0.0, 1.0])) + return m1, m2, mortar + + def runTest(self): + m1, m2, mortar = self.create_meshes() + + # check that f2t indices make sense + self.assertTrue((mortar.f2t[0, :] >= 0).all()) + self.assertTrue((mortar.f2t[0, :] <= np.max(m1.t)).all()) + self.assertTrue((mortar.f2t[1, :] >= 0).all()) + self.assertTrue((mortar.f2t[1, :] <= np.max(m2.t)).all()) + + +class TestMortarMeshTriTri(TestMortarMeshTriQuad): + + def create_mesh2(self): + m2 = MeshTri() + m2.refine(4) + m2.translate((1.0, -0.5)) + return m2 + + if __name__ == '__main__': unittest.main() From b7f5037f3efc1eaa80a9cf96ea2e22a2924c0d5c Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sat, 11 Jan 2020 23:05:44 +0200 Subject: [PATCH 05/32] Allow multiple inputs to MappingAffine.G --- skfem/mapping/mapping_affine.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py index 542abe436..45a4beed7 100644 --- a/skfem/mapping/mapping_affine.py +++ b/skfem/mapping/mapping_affine.py @@ -152,7 +152,12 @@ def G(self, X, find=None): else: B, c = self.B[:, :, find], self.c[:, find] - return (np.einsum('ijk,jl', B, X).T + c.T).T + if len(X.shape) == 2: + return (np.einsum('ijk,jl', B, X).T + c.T).T + elif len(X.shape) == 3: + return (np.einsum('ijk,jkl->ikl', B, X).T + c.T).T + else: + raise Exception("Wrong dimension of input.") def detDG(self, X, find=None): if find is None: From 0e7e567322e9226755ad874a3dda8b33d97ddc00 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sun, 12 Jan 2020 00:45:58 +0200 Subject: [PATCH 06/32] Refactor MeshMortar and MortarBasis Uses a completely new logic for evaluating the global basis functions in MortarBasis constructor. The idea is to use an explicit d-1 dimensional space for the mortaring mesh. In this d-1 dimensional space we have three meshes: two helper meshes created from the original facet meshes via projection and one additional supermesh which contains all nodes from both helper meshes. Quadrature rule is generated for the supermesh and then inverse mapped using the helper mesh reference mappings. Next, the reference (facet) mapping for the original mesh is used to get global quadrature points. Finally, these global quadrature points are inverse mapped to the reference (volume) element of the original mesh and basis functions evaluated. This elaborate technique hopefully allows enough flexibility for general use cases in the future. --- docs/examples/ex04.py | 13 +-- skfem/assembly/basis/mortar_basis.py | 28 ++++-- skfem/mesh/mesh_mortar.py | 123 ++++++++++++++++----------- 3 files changed, 98 insertions(+), 66 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index c714b9a33..a2cb07465 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -9,25 +9,20 @@ M = MeshLine(np.linspace(0, 1, 6)) M = M*M M.translate((1.0, 0.0)) +M = M._splitquads() e1 = ElementTriP1() e = ElementVectorH1(e1) -E1 = ElementQuad1() +E1 = ElementTriP1() E = ElementVectorH1(E1) ib = InteriorBasis(m, e) Ib = InteriorBasis(M, E) -def rule(x): - return (x[0]==1.0) - -def param(x, y): - return y - mortar = MeshMortar.init_1D(m, M, - m.facets_satisfying(lambda x: x[0]==1.0), - M.facets_satisfying(lambda x: x[0]==1.0), + m.facets_satisfying(lambda x: x[0] == 1.0), + M.facets_satisfying(lambda x: x[0] == 1.0), np.array([0.0, 1.0])) mb = [ diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index 6dff2997a..132d6a882 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -36,21 +36,37 @@ def __init__(self, "and use consistent orders for both sides.") # build mapping for the mortar mesh - super(MortarBasis, self).__init__(mesh, elem, None, intorder) + #super(MortarBasis, self).__init__(mesh, elem, None, intorder) + + self.mesh = mesh.target_mesh[side] + #self.brefdom = mesh.brefdom + # self.refdom = mesh.brefdom + # mesh.refdom = mesh.brefdom + #self.intorder = 3 + #self.mapping = mesh.mapping() + self.elem = elem # build dofnum and mappings for the target mesh self._build_dofnum(mesh.target_mesh[side], elem) + integ_mapping = mesh.supermesh.mapping() target_mapping = mesh.target_mesh[side].mapping() + helper_mapping = mesh.helper_mesh[side].mapping() + self.Nbfun = self.element_dofs.shape[0] - self.X, self.W = get_quadrature(self.brefdom, self.intorder) + self.X, self.W = get_quadrature('line', 3) - self.find = np.arange(self.mesh.facets.shape[1]) - self.tind = self.mesh.f2t[side, self.find] + # self.find = np.arange(self.mesh.facets.shape[1]) + self.find = mesh.I[side] + self.tind = mesh.target_mesh[side].f2t[0, mesh.I[side]] + self.mapping = target_mapping # boundary refdom to global facet - x = self.mapping.G(self.X, find=self.find) + x = integ_mapping.F(self.X) + x = helper_mapping.invF(x, tind=mesh.ix[side]) + x = target_mapping.G(x, find=mesh.I[side]) + # global facet to refdom facet - Y = target_mapping.invF(x, tind=self.tind) + Y = target_mapping.invF(x, tind=mesh.target_mesh[side].f2t[0, mesh.I[side]]) # normals are defined in the mortar mesh self.normals = np.repeat(mesh.normals[:, :, None], diff --git a/skfem/mesh/mesh_mortar.py b/skfem/mesh/mesh_mortar.py index 18c5ba085..89aedd8c8 100644 --- a/skfem/mesh/mesh_mortar.py +++ b/skfem/mesh/mesh_mortar.py @@ -1,6 +1,7 @@ import numpy as np from .mesh import Mesh +from .mesh_line import MeshLine from ..mapping import MappingAffine @@ -9,7 +10,7 @@ class MeshMortar(Mesh): name = "Mortar" - def __init__(self, mesh1, mesh2, p, facets, f2t, normals): + def __init__(self, m_super, m1, m2, mesh1, mesh2, normals, ix1,ix2,I1,I2): """Create an interface mesh for mortar methods. Parameters @@ -22,20 +23,22 @@ def __init__(self, mesh1, mesh2, p, facets, f2t, normals): normals """ - self.p = p - self.facets = facets - self.f2t = f2t + self.p = m_super.p + self.t = m_super.t self.normals = normals + self.supermesh = m_super + self.helper_mesh = 2 * [None] + self.helper_mesh[0] = m1 + self.helper_mesh[1] = m2 self.target_mesh = 2 * [None] self.target_mesh[0] = mesh1 self.target_mesh[1] = mesh2 - - # for quadrature rules - self.brefdom = mesh1.brefdom - - # dummy values so that initialising Mapping doesn't crash - self.t = mesh1.t[:, :1] - self.t2f = -1 + 0 * self.t + self.ix={} + self.ix[0]=ix1 + self.ix[1]=ix2 + self.I ={} + self.I[0]=I1 + self.I[1]=I2 @classmethod @@ -71,47 +74,65 @@ def param(p): y = proj(p) return np.linalg.norm(y, axis=0) * np.dot(tangent, y) + param_p1 = param(p1) + param_p2 = param(p2) + # find unique facets by combining nodes from both sides - _, ix = np.unique(np.concatenate((param(p1), param(p2))), return_index=True) + _, ix = np.unique(np.concatenate((param_p1, param_p2)), return_index=True) ixorig = np.concatenate((p1_ix, p2_ix + mesh1.p.shape[1]))[ix] - p = np.hstack((mesh1.p, mesh2.p)) # TODO has duplicate nodes - facets = np.array([ixorig[:-1], ixorig[1:]]) + p = np.array([np.hstack((param(mesh1.p), param(mesh2.p)))]) + t = np.array([ixorig[:-1], ixorig[1:]]) + + # supermesh + p = p[:, np.concatenate((t[0], np.array([t[1, -1]])))] + t = np.array([np.arange(p.shape[1] - 1), np.arange(1, p.shape[1])]) + m_super = MeshLine(p, t) + + # helper meshes + m1 = MeshLine(np.sort(param_p1), np.array([np.arange(p1.shape[1] - 1), + np.arange(1, p1.shape[1])])) + m2 = MeshLine(np.sort(param_p2), np.array([np.arange(p2.shape[1] - 1), + np.arange(1, p2.shape[1])])) - # construct normals + # construct normals by rotating 'tangent' normal = np.array([tangent[1], -tangent[0]]) - normals = normal[:, None].repeat(facets.shape[1], axis=1) - - # mappings from facets to the original triangles - f2t = facets * 0 - 1 - for itr in range(facets.shape[1]): - mp = .5 * (p[:, facets[0, itr]] + p[:, facets[1, itr]]) - val = param(mp) - for jtr in boundary1: - x1 = mesh1.p[:, mesh1.facets[0, jtr]] - x2 = mesh1.p[:, mesh1.facets[1, jtr]] - if (val > param(x1) and val < param(x2) or - val < param(x1) and val > param(x2) or - val >= param(x1) and val < param(x2) or - val > param(x1) and val <= param(x2) or - val <= param(x1) and val > param(x2) or - val < param(x1) and val >= param(x2)): - f2t[0, itr] = mesh1.f2t[0, jtr] - break - for jtr in boundary2: - x1 = mesh2.p[:, mesh2.facets[0, jtr]] - x2 = mesh2.p[:, mesh2.facets[1, jtr]] - if (val > param(x1) and val < param(x2) or - val < param(x1) and val > param(x2) or - val >= param(x1) and val < param(x2) or - val > param(x1) and val <= param(x2) or - val <= param(x1) and val > param(x2) or - val < param(x1) and val >= param(x2)): - f2t[1, itr] = mesh2.f2t[0, jtr] - break - if not (f2t > -1).all(): - raise Exception("All mesh facets corresponding to mortar facets not found! {}".format(f2t)) - - return cls(mesh1, mesh2, p, facets, f2t, normals) - - def mapping(self): - return MappingAffine(self) + normals = normal[:, None].repeat(t.shape[1], axis=1) + + # initialize mappings for orienting + map_super = m_super.mapping() + map_m1 = m1.mapping() + map_m2 = m2.mapping() + map_mesh1 = mesh1.mapping() + map_mesh2 = mesh2.mapping() + + # orient helper meshes + mps = map_super.F(np.array([[0.5]])) + ix1 = np.digitize(mps[0,:,0], m1.p[0]) - 1 + ix2 = np.digitize(mps[0,:,0], m2.p[0]) - 1 + + # for each element, map two points to global coordinates, reparametrize + # the points, and flip corresponding helper mesh element indices if + # sorting is wrong + f1mps = .5 * (mesh1.p[:, mesh1.facets[0, boundary1]] + + mesh1.p[:, mesh1.facets[1, boundary1]]) + sort_boundary1 = np.argsort(param(f1mps)) + z1 = map_mesh1.G(map_m1.invF(map_super.F(np.array([[0.25, 0.75]])), tind=ix1), + find=boundary1[sort_boundary1][ix1]) + ix1_flip = np.unique(ix1[param(z1[:,:,1]) < param(z1[:,:,0])]) + m1.t[:, ix1_flip] = np.flipud(m1.t[:, ix1_flip]) + + f2mps = .5 * (mesh2.p[:, mesh2.facets[0, boundary2]] + + mesh2.p[:, mesh2.facets[1, boundary2]]) + sort_boundary2 = np.argsort(param(f2mps)) + z2 = map_mesh2.G(map_m2.invF(map_super.F(np.array([[0.25, 0.75]])), tind=ix2), + find=boundary2[sort_boundary2][ix2]) + ix2_flip = np.unique(ix2[param(z2[:,:,1]) < param(z2[:,:,0])]) + m2.t[:, ix2_flip] = np.flipud(m2.t[:, ix2_flip]) + + return cls(m_super, + m1, m2, + mesh1, mesh2, + normals, + ix1, ix2, + boundary1[sort_boundary1][ix1], + boundary2[sort_boundary2][ix2]) From 055236b7c90132c8d8110aa756c61e1980e3bf06 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 13 Jan 2020 10:52:52 +0200 Subject: [PATCH 07/32] Allow element-wise inputs to detDG --- skfem/mapping/mapping_affine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py index 45a4beed7..b3aa5d9f2 100644 --- a/skfem/mapping/mapping_affine.py +++ b/skfem/mapping/mapping_affine.py @@ -165,7 +165,7 @@ def detDG(self, X, find=None): else: detDG = self.detB[find] - return np.tile(detDG, (X.shape[1], 1)).T + return np.tile(detDG, (X.shape[-1], 1)).T def normals(self, X, tind, find, t2f): if self.dim == 1: From 8498fac4f087d06aabd9955a487d040b74301dbc Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 13 Jan 2020 10:53:14 +0200 Subject: [PATCH 08/32] Further testing of the approach --- docs/examples/ex04.py | 9 ++++++--- skfem/assembly/basis/mortar_basis.py | 10 +++++----- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index a2cb07465..f14105420 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -11,6 +11,9 @@ M.translate((1.0, 0.0)) M = M._splitquads() +m.refine(3) +M.refine(3) + e1 = ElementTriP1() e = ElementVectorH1(e1) @@ -48,7 +51,7 @@ weakform1 = linear_elasticity(Lambda=Lambda1, Mu=Mu) weakform2 = linear_elasticity(Lambda=Lambda2, Mu=Mu) -alpha = 1 +alpha = 0.001 K1 = asm(weakform1, ib) K2 = asm(weakform2, Ib) L = [[None,None],[None,None]] @@ -79,13 +82,13 @@ def Eps(dw): n[1]*C(Eps(dv))[1, 0]*n[0] +\ n[1]*C(Eps(dv))[1, 1]*n[1]) h = w.h - return 1.0/(alpha*h)*ju*jv - mu*jv - mv*ju + return 1.0/(alpha*h)*ju*jv - 0*mu*jv - 0*mv*ju L[i][j] = asm(bilin_penalty, mb[i], mb[j]) @linear_form def load(v, dv, w): - return -50*v[1] + return 500*v[0] f1 = asm(load, ib) f2 = np.zeros(K2.shape[0]) diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index 132d6a882..f06f3ab2b 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -62,11 +62,11 @@ def __init__(self, # boundary refdom to global facet x = integ_mapping.F(self.X) - x = helper_mapping.invF(x, tind=mesh.ix[side]) - x = target_mapping.G(x, find=mesh.I[side]) - + X = helper_mapping.invF(x, tind=mesh.ix[side]) + x = target_mapping.G(X, find=self.find) + # global facet to refdom facet - Y = target_mapping.invF(x, tind=mesh.target_mesh[side].f2t[0, mesh.I[side]]) + Y = target_mapping.invF(x, tind=self.tind) # normals are defined in the mortar mesh self.normals = np.repeat(mesh.normals[:, :, None], @@ -77,7 +77,7 @@ def __init__(self, for j in range(self.Nbfun)] self.nelems = len(self.find) - self.dx = np.abs(self.mapping.detDG(self.X, find=self.find)) *\ + self.dx = np.abs(self.mapping.detDG(X, find=self.find)) *\ np.tile(self.W, (self.nelems, 1)) self.element_dofs = self.element_dofs[:, self.tind] From ac485a5cdb1456007d4f6858b60f9d3ab6fe2c1e Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 13 Jan 2020 17:36:26 +0200 Subject: [PATCH 09/32] Build mortar mappings instead of mortar meshes --- docs/examples/ex04.py | 41 ++++++++++------- skfem/__init__.py | 18 +++++++- skfem/assembly/basis/mortar_basis.py | 67 ++++++++++------------------ skfem/mesh/mesh_mortar.py | 39 +++++++++++----- 4 files changed, 92 insertions(+), 73 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index f14105420..2c59fe710 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -4,15 +4,19 @@ import numpy as np # create meshes -m = MeshTri() -m.refine(3) +m = MeshTri.init_symmetric() +#m = MeshTri() +m.refine(4) M = MeshLine(np.linspace(0, 1, 6)) M = M*M +#M = MeshTri.init_symmetric() +#M = MeshTri() +M.refine(1) M.translate((1.0, 0.0)) M = M._splitquads() -m.refine(3) -M.refine(3) +#m.refine(3) +#M.refine(3) e1 = ElementTriP1() e = ElementVectorH1(e1) @@ -29,8 +33,8 @@ np.array([0.0, 1.0])) mb = [ - MortarBasis(mortar, e, intorder=2, side=0), - MortarBasis(mortar, E, intorder=2, side=1) + MortarBasis(m, e, mapping = mortar[0], intorder=2), + MortarBasis(M, E, mapping = mortar[1], intorder=2) ] E1 = 1000.0 @@ -39,19 +43,20 @@ nu1 = 0.3 nu2 = 0.3 -Mu1 = E1/(2.0*(1.0 + nu1)) -Mu2 = E2/(2.0*(1.0 + nu2)) +Mu1 = E1 / (2. * (1. + nu1)) +Mu2 = E2 / (2. * (1. + nu2)) -Lambda1 = E1*nu1/((1.0 + nu1)*(1.0 - 2.0*nu1)) -Lambda2 = E2*nu2/((1.0 + nu2)*(1.0 - 2.0*nu2)) +Lambda1 = E1 * nu1 / ((1. + nu1) * (1. - 2. * nu1)) +Lambda2 = E2 * nu2 / ((1. + nu2) * (1. - 2. * nu2)) Mu = Mu1 Lambda = Lambda1 -weakform1 = linear_elasticity(Lambda=Lambda1, Mu=Mu) -weakform2 = linear_elasticity(Lambda=Lambda2, Mu=Mu) +weakform1 = linear_elasticity(Lambda=Lambda, Mu=Mu) +weakform2 = linear_elasticity(Lambda=Lambda, Mu=Mu) -alpha = 0.001 + +alpha = 1 K1 = asm(weakform1, ib) K2 = asm(weakform2, Ib) L = [[None,None],[None,None]] @@ -60,8 +65,8 @@ @bilinear_form def bilin_penalty(u, du, v, dv, w): n = w.n - ju = (-1.0)**i*(u[0]*n[0] + u[1]*n[1]) - jv = (-1.0)**j*(v[0]*n[0] + v[1]*n[1]) + ju = (-1.) ** i * (u[0] * n[0] + u[1] * n[1]) + jv = (-1.) ** j * (v[0] * n[0] + v[1] * n[1]) def tr(T): return T[0, 0] + T[1, 1] @@ -82,13 +87,15 @@ def Eps(dw): n[1]*C(Eps(dv))[1, 0]*n[0] +\ n[1]*C(Eps(dv))[1, 1]*n[1]) h = w.h - return 1.0/(alpha*h)*ju*jv - 0*mu*jv - 0*mv*ju + return 1.0/(alpha*h)*ju*jv - mu*jv - mv*ju L[i][j] = asm(bilin_penalty, mb[i], mb[j]) +import pdb; pdb.set_trace() + @linear_form def load(v, dv, w): - return 500*v[0] + return 50*v[1] f1 = asm(load, ib) f2 = np.zeros(K2.shape[0]) diff --git a/skfem/__init__.py b/skfem/__init__.py index c3d94aa54..2e9d35641 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -1,9 +1,25 @@ +"""Support for wildcard import.""" + + +def mark_unstable(fun): + def call(*args, **kwargs): + import warnings + warnings.warn("You are using an unstable feature '{}' that might " + "change in the near future. This is not yet part " + "of the publicly documented API and hence is not " + "governed by the usual assurance of semantic " + "versioning.".format(fun.__module__)) + return fun(*args, **kwargs) + return call + + from skfem.mesh import * from skfem.assembly import * from skfem.mapping import * from skfem.element import * from skfem.utils import * + __all__ = ['Mesh', 'Mesh2D', 'Mesh3D', @@ -58,4 +74,4 @@ 'ElementTriRT0', 'ElementVectorH1', 'ElementLineP1', - 'ElementLineP2',] + 'ElementLineP2'] diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index f06f3ab2b..b24ddf586 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -1,88 +1,69 @@ import numpy as np from numpy import ndarray -from skfem.quadrature import get_quadrature +from skfem import mark_unstable +from ...quadrature import get_quadrature from .basis import Basis from .interior_basis import InteriorBasis class MortarBasis(Basis): - """Global basis functions evaluated at integration points on the mortar - boundary. """ + """Global basis functions at integration points on the mortar boundary.""" + @mark_unstable def __init__(self, mesh, elem, - intorder: int = None, - side: int = 0): + mapping, + intorder): """Initialize a basis for assembling mortar matrices. Parameters ---------- mesh - An object of type :class:`~skfem.mesh.MeshMortar`. + An object of type :class:`~skfem.mesh.Mesh`. elem An object of type :class:`~skfem.element.Element`. + mapping + Mortar mapping. intorder - Required integration order, i.e. the degree of polynomials that are - integrated exactly by the used quadrature. - side - The numbers 0 and 1 refer to the different sides of the mortar - interface. Side 0 corresponds to the indices mesh.f2t[0, :]. + Integration order, i.e. the degree of polynomials that are + integrated exactly by the used quadrature. Please use equivalent + integration orders for the both sides of the mortar boundary. """ if intorder is None: raise ValueError("Please specify 'intorder' for MortarBasis " "and use consistent orders for both sides.") - # build mapping for the mortar mesh - #super(MortarBasis, self).__init__(mesh, elem, None, intorder) + super(MortarBasis, self).__init__(mesh, elem, mapping, intorder) - self.mesh = mesh.target_mesh[side] - #self.brefdom = mesh.brefdom - # self.refdom = mesh.brefdom - # mesh.refdom = mesh.brefdom - #self.intorder = 3 - #self.mapping = mesh.mapping() - self.elem = elem + self.X, self.W = get_quadrature(self.brefdom, self.intorder) - # build dofnum and mappings for the target mesh - self._build_dofnum(mesh.target_mesh[side], elem) - integ_mapping = mesh.supermesh.mapping() - target_mapping = mesh.target_mesh[side].mapping() - helper_mapping = mesh.helper_mesh[side].mapping() - self.Nbfun = self.element_dofs.shape[0] - - self.X, self.W = get_quadrature('line', 3) - - # self.find = np.arange(self.mesh.facets.shape[1]) - self.find = mesh.I[side] - self.tind = mesh.target_mesh[side].f2t[0, mesh.I[side]] - self.mapping = target_mapping + # global facets where basis is evaluated + self.find = self.mapping.find + self.tind = self.mesh.f2t[0, self.find] # boundary refdom to global facet - x = integ_mapping.F(self.X) - X = helper_mapping.invF(x, tind=mesh.ix[side]) - x = target_mapping.G(X, find=self.find) + x = self.mapping.G(self.X) # global facet to refdom facet - Y = target_mapping.invF(x, tind=self.tind) + Y = self.mapping.invF(x, tind=self.tind) - # normals are defined in the mortar mesh - self.normals = np.repeat(mesh.normals[:, :, None], + # normals are defined in the mortar mapping + self.normals = np.repeat(self.mapping.normals[:, :, None], len(self.W), axis=2) - self.basis = [self.elem.gbasis(target_mapping, Y, j, self.tind) + self.basis = [self.elem.gbasis(self.mapping, Y, j, self.tind) for j in range(self.Nbfun)] self.nelems = len(self.find) - self.dx = np.abs(self.mapping.detDG(X, find=self.find)) *\ + self.dx = np.abs(self.mapping.detDG(self.X, find=self.find)) *\ np.tile(self.W, (self.nelems, 1)) self.element_dofs = self.element_dofs[:, self.tind] - def default_parameters(self): """Return default parameters for `~skfem.assembly.asm`.""" return {'x': self.global_coordinates(), @@ -90,7 +71,7 @@ def default_parameters(self): 'n': self.normals} def global_coordinates(self) -> ndarray: - return self.mapping.G(self.X, find=self.find) + return self.mapping.G(self.X) def mesh_parameters(self) -> ndarray: if self.mesh.dim() == 1: diff --git a/skfem/mesh/mesh_mortar.py b/skfem/mesh/mesh_mortar.py index 89aedd8c8..9c99f9c94 100644 --- a/skfem/mesh/mesh_mortar.py +++ b/skfem/mesh/mesh_mortar.py @@ -72,7 +72,7 @@ def proj(p): def param(p): """Calculate signed distances of projected points from origin.""" y = proj(p) - return np.linalg.norm(y, axis=0) * np.dot(tangent, y) + return np.linalg.norm(y, axis=0) * np.sign(np.dot(tangent, y)) param_p1 = param(p1) param_p2 = param(p2) @@ -107,8 +107,8 @@ def param(p): # orient helper meshes mps = map_super.F(np.array([[0.5]])) - ix1 = np.digitize(mps[0,:,0], m1.p[0]) - 1 - ix2 = np.digitize(mps[0,:,0], m2.p[0]) - 1 + ix1 = np.digitize(mps[0, :, 0], m1.p[0]) - 1 + ix2 = np.digitize(mps[0, :, 0], m2.p[0]) - 1 # for each element, map two points to global coordinates, reparametrize # the points, and flip corresponding helper mesh element indices if @@ -118,7 +118,7 @@ def param(p): sort_boundary1 = np.argsort(param(f1mps)) z1 = map_mesh1.G(map_m1.invF(map_super.F(np.array([[0.25, 0.75]])), tind=ix1), find=boundary1[sort_boundary1][ix1]) - ix1_flip = np.unique(ix1[param(z1[:,:,1]) < param(z1[:,:,0])]) + ix1_flip = np.unique(ix1[param(z1[:, :, 1]) < param(z1[:, :, 0])]) m1.t[:, ix1_flip] = np.flipud(m1.t[:, ix1_flip]) f2mps = .5 * (mesh2.p[:, mesh2.facets[0, boundary2]] + @@ -126,13 +126,28 @@ def param(p): sort_boundary2 = np.argsort(param(f2mps)) z2 = map_mesh2.G(map_m2.invF(map_super.F(np.array([[0.25, 0.75]])), tind=ix2), find=boundary2[sort_boundary2][ix2]) - ix2_flip = np.unique(ix2[param(z2[:,:,1]) < param(z2[:,:,0])]) + ix2_flip = np.unique(ix2[param(z2[:, :, 1]) < param(z2[:, :, 0])]) m2.t[:, ix2_flip] = np.flipud(m2.t[:, ix2_flip]) - return cls(m_super, - m1, m2, - mesh1, mesh2, - normals, - ix1, ix2, - boundary1[sort_boundary1][ix1], - boundary2[sort_boundary2][ix2]) + + map_m1 = m1.mapping() + map_m2 = m2.mapping() + + nmap_mesh1 = mesh1.mapping() + nmap_mesh2 = mesh2.mapping() + + nmap_mesh1.G = lambda X: map_mesh1.G(map_m1.invF(map_super.F(X), tind=ix1), + find=boundary1[sort_boundary1][ix1]) + nmap_mesh2.G = lambda X: map_mesh2.G(map_m2.invF(map_super.F(X), tind=ix2), + find=boundary2[sort_boundary2][ix2]) + nmap_mesh1.find = boundary1[sort_boundary1][ix1] + nmap_mesh2.find = boundary2[sort_boundary2][ix2] + nmap_mesh1.normals = normals + nmap_mesh2.normals = normals + + segments1 = nmap_mesh1.G(np.array([[0.0, 1.0]])) + segments2 = nmap_mesh2.G(np.array([[0.0, 1.0]])) + nmap_mesh1.detDG = lambda *_,**__: np.sqrt(np.diff(segments1[0], axis=1)**2 + np.diff(segments1[1], axis=1)**2) + nmap_mesh2.detDG = lambda *_,**__: np.sqrt(np.diff(segments2[0], axis=1)**2 + np.diff(segments2[1], axis=1)**2) + + return nmap_mesh1, nmap_mesh2 From 40d259edc1561b5c361492d665812f4d6ca1f787 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 13 Jan 2020 22:47:37 +0200 Subject: [PATCH 10/32] Change warning message --- skfem/__init__.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/skfem/__init__.py b/skfem/__init__.py index 2e9d35641..7ea81a8f1 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -5,10 +5,8 @@ def mark_unstable(fun): def call(*args, **kwargs): import warnings warnings.warn("You are using an unstable feature '{}' that might " - "change in the near future. This is not yet part " - "of the publicly documented API and hence is not " - "governed by the usual assurance of semantic " - "versioning.".format(fun.__module__)) + "change in the near future. It is not part of the " + "publicly documented API.".format(fun.__module__)) return fun(*args, **kwargs) return call From 900a0f2683f2da7d7c93476ac9ea71d33e9f2220 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 13 Jan 2020 22:47:54 +0200 Subject: [PATCH 11/32] Remove ValueError on missing intorder --- skfem/assembly/basis/mortar_basis.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index b24ddf586..87bb734c0 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -29,13 +29,9 @@ def __init__(self, intorder Integration order, i.e. the degree of polynomials that are integrated exactly by the used quadrature. Please use equivalent - integration orders for the both sides of the mortar boundary. + integration orders on both sides of the mortar boundary. """ - if intorder is None: - raise ValueError("Please specify 'intorder' for MortarBasis " - "and use consistent orders for both sides.") - super(MortarBasis, self).__init__(mesh, elem, mapping, intorder) self.X, self.W = get_quadrature(self.brefdom, self.intorder) @@ -77,4 +73,5 @@ def mesh_parameters(self) -> ndarray: if self.mesh.dim() == 1: return np.array([0.0]) else: - return np.abs(self.mapping.detDG(self.X, self.find)) ** (1.0 / (self.mesh.dim() - 1)) + return (np.abs(self.mapping.detDG(self.X, find=self.find)) + ** (1.0 / (self.mesh.dim() - 1))) From f7a642cde66d44aa9a30fe4d5d3e94477fd773fd Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 15 Jan 2020 13:44:15 +0200 Subject: [PATCH 12/32] Fix shapes of the boundary mapping outputs --- skfem/mapping/mapping_isoparametric.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skfem/mapping/mapping_isoparametric.py b/skfem/mapping/mapping_isoparametric.py index 4102eab08..690fc58d8 100644 --- a/skfem/mapping/mapping_isoparametric.py +++ b/skfem/mapping/mapping_isoparametric.py @@ -81,7 +81,7 @@ def bndmap(i, X, find=None): out += p[i, facets[itr, :]][:, None] * phi return out else: - out = np.zeros((len(find), X.shape[1])) + out = np.zeros((len(find), X.shape[-1])) for itr in range(facets.shape[0]): phi, _ = bndelem.lbasis(X, itr) out += p[i, facets[itr, find]][:, None] * phi @@ -95,7 +95,7 @@ def bndJ(i, j, X, find=None): out += p[i, facets[itr, :]][:, None] * dphi[j] return out else: - out = np.zeros((len(find), X.shape[1])) + out = np.zeros((len(find), X.shape[-1])) for itr in range(facets.shape[0]): _, dphi = bndelem.lbasis(X, itr) out += p[i, facets[itr, find]][:, None] * dphi[j] From 178ca578107f69a039fc82cbe8c8f67f8e9fd3ad Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 15 Jan 2020 13:44:49 +0200 Subject: [PATCH 13/32] Improve docstring --- skfem/assembly/basis/basis.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/skfem/assembly/basis/basis.py b/skfem/assembly/basis/basis.py index 85ba71bb4..629c8e4f5 100644 --- a/skfem/assembly/basis/basis.py +++ b/skfem/assembly/basis/basis.py @@ -193,6 +193,12 @@ def get_dofs(self, facets: Optional[Any] = None): Mesh.facets_satisfying for each entry to get facets and then find dofs for those. + Returns + ------- + Dofs + A subset of degrees-of-freedom as :class:`skfem.assembly.dofs.Dofs`. + + """ if facets is None: facets = self.mesh.boundary_facets() From f6e626f07a3ea45f6904ede9c0cd8bc9fedffd29 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 15 Jan 2020 14:43:02 +0200 Subject: [PATCH 14/32] Add ElementQuadDG --- skfem/__init__.py | 1 + skfem/element/__init__.py | 4 ++-- skfem/element/element_quad/__init__.py | 1 + skfem/element/element_quad/element_quad_dg.py | 16 ++++++++++++++++ 4 files changed, 20 insertions(+), 2 deletions(-) create mode 100644 skfem/element/element_quad/element_quad_dg.py diff --git a/skfem/__init__.py b/skfem/__init__.py index 7ea81a8f1..a97b320a4 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -60,6 +60,7 @@ def call(*args, **kwargs): 'ElementQuad0', 'ElementQuad1', 'ElementQuad2', + 'ElementQuadDG', 'ElementTetN0', 'ElementTetP0', 'ElementTetP1', diff --git a/skfem/element/__init__.py b/skfem/element/__init__.py index 71f76c25c..2c63b6dce 100644 --- a/skfem/element/__init__.py +++ b/skfem/element/__init__.py @@ -21,9 +21,9 @@ from .element_tri import ElementTriP1, ElementTriP2, ElementTriDG,\ ElementTriP0, ElementTriRT0, ElementTriMorley,\ ElementTriArgyris -from .element_quad import ElementQuad0, ElementQuad1, ElementQuad2 +from .element_quad import ElementQuad0, ElementQuad1, ElementQuad2,\ + ElementQuadDG from .element_tet import ElementTetP0, ElementTetP1, ElementTetP2,\ ElementTetRT0, ElementTetN0 from .element_hex import ElementHex1 from .element_line import ElementLineP1, ElementLineP2 - diff --git a/skfem/element/element_quad/__init__.py b/skfem/element/element_quad/__init__.py index e3788daa7..58d9f2de4 100644 --- a/skfem/element/element_quad/__init__.py +++ b/skfem/element/element_quad/__init__.py @@ -1,3 +1,4 @@ from .element_quad0 import ElementQuad0 from .element_quad1 import ElementQuad1 from .element_quad2 import ElementQuad2 +from .element_quad_dg import ElementQuadDG diff --git a/skfem/element/element_quad/element_quad_dg.py b/skfem/element/element_quad/element_quad_dg.py new file mode 100644 index 000000000..01f1d023d --- /dev/null +++ b/skfem/element/element_quad/element_quad_dg.py @@ -0,0 +1,16 @@ +from ..element_h1 import ElementH1 + + +class ElementQuadDG(ElementH1): + dim = 2 + + def __init__(self, elem): + # change all dofs to interior dofs + self.elem = elem + self.maxdeg = elem.maxdeg + self.interior_dofs = (4 * elem.nodal_dofs + + 4 * elem.facet_dofs + + elem.interior_dofs) + + def lbasis(self, X, i): + return self.elem.lbasis(X, i) From 1aeb3e616f59222f4b3279a447841a429d37715b Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Wed, 15 Jan 2020 15:39:40 +0200 Subject: [PATCH 15/32] More generic example for ex04 --- docs/examples/ex04.py | 185 ++++++++++++++++++++++++++++++------------ 1 file changed, 131 insertions(+), 54 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index 2c59fe710..a5ecd1713 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -2,39 +2,57 @@ from skfem.models.elasticity import linear_elasticity from skfem.visuals.matplotlib import * import numpy as np +from pygmsh import generate_mesh +from pygmsh.built_in import Geometry +from skfem.importers import from_meshio -# create meshes -m = MeshTri.init_symmetric() -#m = MeshTri() -m.refine(4) -M = MeshLine(np.linspace(0, 1, 6)) -M = M*M -#M = MeshTri.init_symmetric() -#M = MeshTri() -M.refine(1) +geom = Geometry() + +points = [] +lines = [] +points.append(geom.add_point([0., 0., 0.], .1)) +points.append(geom.add_point([0., 1., 0.], .1)) +points.append(geom.add_point([0.,-1., 0.], .1)) +lines.append(geom.add_circle_arc(points[2], points[0], points[1])) +geom.add_physical(lines[-1], 'contact') +lines.append(geom.add_line(points[1], points[2])) +geom.add_physical(lines[-1], 'dirichlet') +geom.add_physical(geom.add_plane_surface(geom.add_line_loop(lines)), 'domain') + +m = from_meshio(generate_mesh(geom, dim=2)) + + +M = MeshLine(np.linspace(0, 1, 6)) * MeshLine(np.linspace(-1, 1, 10)) M.translate((1.0, 0.0)) -M = M._splitquads() +M.refine() +#M = M._splitquads() + +ax = draw(m) +draw(M, ax=ax) +show() #m.refine(3) #M.refine(3) -e1 = ElementTriP1() +e1 = ElementTriP2() e = ElementVectorH1(e1) -E1 = ElementTriP1() +E1 = ElementQuad2() E = ElementVectorH1(E1) -ib = InteriorBasis(m, e) -Ib = InteriorBasis(M, E) +ib = InteriorBasis(m, e, intorder=4) +Ib = InteriorBasis(M, E, intorder=4) +fb = FacetBasis(m, e) +Fb = FacetBasis(M, E) mortar = MeshMortar.init_1D(m, M, - m.facets_satisfying(lambda x: x[0] == 1.0), + m.boundaries['contact'], M.facets_satisfying(lambda x: x[0] == 1.0), np.array([0.0, 1.0])) mb = [ - MortarBasis(m, e, mapping = mortar[0], intorder=2), - MortarBasis(M, E, mapping = mortar[1], intorder=2) + MortarBasis(m, e, mapping = mortar[0], intorder=4), + MortarBasis(M, E, mapping = mortar[1], intorder=4) ] E1 = 1000.0 @@ -56,10 +74,24 @@ weakform2 = linear_elasticity(Lambda=Lambda, Mu=Mu) -alpha = 1 +alpha = 1000 +limit = 0.3 K1 = asm(weakform1, ib) K2 = asm(weakform2, Ib) -L = [[None,None],[None,None]] +K = [[None,None],[None,None]] +f = [None] * 2 + +def tr(T): + return T[0, 0] + T[1, 1] + +def C(T): + return np.array([[2*Mu*T[0, 0] + Lambda*tr(T), 2*Mu*T[0, 1]], + [2*Mu*T[1, 0], 2*Mu*T[1, 1] + Lambda*tr(T)]]) + +def Eps(dw): + return np.array([[dw[0][0], 0.5*(dw[0][1] + dw[1][0])], + [0.5*(dw[1][0] + dw[0][1]), dw[1][1]]]) + for i in range(2): for j in range(2): @bilinear_form @@ -67,47 +99,47 @@ def bilin_penalty(u, du, v, dv, w): n = w.n ju = (-1.) ** i * (u[0] * n[0] + u[1] * n[1]) jv = (-1.) ** j * (v[0] * n[0] + v[1] * n[1]) - - def tr(T): - return T[0, 0] + T[1, 1] - - def C(T): - return np.array([[2*Mu*T[0, 0] + Lambda*tr(T), 2*Mu*T[0, 1]], - [2*Mu*T[1, 0], 2*Mu*T[1, 1] + Lambda*tr(T)]]) - - def Eps(dw): - return np.array([[dw[0][0], 0.5*(dw[0][1] + dw[1][0])], - [0.5*(dw[1][0] + dw[0][1]), dw[1][1]]]) - mu = 0.5*(n[0]*C(Eps(du))[0, 0]*n[0] +\ - n[0]*C(Eps(du))[0, 1]*n[1] +\ - n[1]*C(Eps(du))[1, 0]*n[0] +\ - n[1]*C(Eps(du))[1, 1]*n[1]) - mv = 0.5*(n[0]*C(Eps(dv))[0, 0]*n[0] +\ - n[0]*C(Eps(dv))[0, 1]*n[1] +\ - n[1]*C(Eps(dv))[1, 0]*n[0] +\ - n[1]*C(Eps(dv))[1, 1]*n[1]) + mu = .5 * (n[0] * C(Eps(du))[0, 0] * n[0] + + n[0] * C(Eps(du))[0, 1] * n[1] + + n[1] * C(Eps(du))[1, 0] * n[0] + + n[1] * C(Eps(du))[1, 1] * n[1]) + mv = .5 * (n[0] * C(Eps(dv))[0, 0] * n[0] + + n[0] * C(Eps(dv))[0, 1] * n[1] + + n[1] * C(Eps(dv))[1, 0] * n[0] + + n[1] * C(Eps(dv))[1, 1] * n[1]) h = w.h - return 1.0/(alpha*h)*ju*jv - mu*jv - mv*ju - - L[i][j] = asm(bilin_penalty, mb[i], mb[j]) - -import pdb; pdb.set_trace() - -@linear_form -def load(v, dv, w): - return 50*v[1] - -f1 = asm(load, ib) -f2 = np.zeros(K2.shape[0]) + return (1. / (alpha * h) * ju * jv - mu * jv - mv * ju) *(np.abs(w.x[1]) <= limit) + + K[j][i] = asm(bilin_penalty, mb[i], mb[j]) + + @linear_form + def lin_penalty(v, dv, w): + n = w.n + jv = (-1.) ** i * (v[0] * n[0] + v[1] * n[1]) + mv = .5 * (n[0] * C(Eps(dv))[0, 0] * n[0] + + n[0] * C(Eps(dv))[0, 1] * n[1] + + n[1] * C(Eps(dv))[1, 0] * n[0] + + n[1] * C(Eps(dv))[1, 1] * n[1]) + h = w.h + def gap(x): + return (1. - np.sqrt(1. - x[1] ** 2)) + return (1. / (alpha * h) * gap(w.x) * jv - gap(w.x) * mv) * (np.abs(w.x[1]) <= limit) + f[i] = asm(lin_penalty, mb[i]) + +K[0][0] += K1 +K[1][1] += K2 + +f1 = f[0] +f2 = np.zeros(K2.shape[0]) + f[1] import scipy.sparse -K = (scipy.sparse.bmat([[K1 + L[0][0], L[1][0]],[L[0][1], K2 + L[1][1]]])).tocsr() +K = (scipy.sparse.bmat(K)).tocsr() i1 = np.arange(K1.shape[0]) i2 = np.arange(K2.shape[0]) + K1.shape[0] -D1 = ib.get_dofs(lambda x: x[0]==0.0).all() -D2 = Ib.get_dofs(lambda x: x[0]==2.0).all() +D1 = ib.get_dofs(lambda x: x[0] == 0.0).all() +D2 = Ib.get_dofs(lambda x: x[0] == 2.0).all() x = np.zeros(K.shape[0]) @@ -117,7 +149,10 @@ def load(v, dv, w): D = np.concatenate((D1, D2 + ib.N)) I = np.setdiff1d(np.arange(K.shape[0]), D) -x = solve(*condense(K, f, I=I)) +x[ib.get_dofs(lambda x: x[0] == 0.0).nodal['u^1']] = 0.1 +x[ib.get_dofs(lambda x: x[0] == 0.0).facet['u^1']] = 0.1 + +x = solve(*condense(K, f, I=I, x=x)) sf = 1 @@ -127,6 +162,48 @@ def load(v, dv, w): M.p[0, :] = M.p[0, :] + sf*x[i2][Ib.nodal_dofs[0, :]] M.p[1, :] = M.p[1, :] + sf*x[i2][Ib.nodal_dofs[1, :]] + +s, S = {}, {} +e_dg = ElementTriDG(ElementTriP1()) +E_dg = ElementQuadDG(ElementQuad1()) + +for itr in range(2): + for jtr in range(2): + @bilinear_form + def proj_cauchy(u, du, v, dv, w): + return C(Eps(du))[itr, jtr] * v + + @bilinear_form + def mass(u, du, v, dv, w): + return u * v + + ib_dg = InteriorBasis(m, e_dg, intorder=4) + Ib_dg = InteriorBasis(M, E_dg, intorder=4) + + s[itr, jtr] = solve(asm(mass, ib_dg), asm(proj_cauchy, ib, ib_dg) @ x[i1]) + S[itr, jtr] = solve(asm(mass, Ib_dg), asm(proj_cauchy, Ib, Ib_dg) @ x[i2]) + +s[2, 2] = nu1*(s[0, 0] + s[1, 1]) +S[2, 2] = nu2*(S[0, 0] + S[1, 1]) + +vonmises1 = np.sqrt(.5*((s[0, 0] - s[1, 1])**2 +\ + (s[1, 1] - s[2, 2])**2 +\ + (s[2, 2] - s[0, 0])**2 +\ + 6.0*s[0, 1]**2)) + +vonmises2 = np.sqrt(.5*((S[0, 0] - S[1, 1])**2 +\ + (S[1, 1] - S[2, 2])**2 +\ + (S[2, 2] - S[0, 0])**2 +\ + 6.0*S[0, 1]**2)) + +ax = plot(ib_dg, vonmises1, shading='gouraud') +draw(m, ax=ax) +plot(Ib_dg, vonmises2, ax=ax, Nrefs=3, shading='gouraud') +draw(M, ax=ax) +show() + +import pdb; pdb.set_trace() + if __name__ == "__main__": from skfem.visuals.matplotlib import draw, show From d8a8ecad787cead8c7519aeeafc988d780d2d826 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Thu, 16 Jan 2020 17:44:51 +0200 Subject: [PATCH 16/32] Rename MeshMortar to MortarPair --- docs/examples/ex04.py | 12 +-- skfem/__init__.py | 2 + skfem/mapping/__init__.py | 1 + .../mesh_mortar.py => mapping/mortar_pair.py} | 89 ++++++++----------- 4 files changed, 44 insertions(+), 60 deletions(-) rename skfem/{mesh/mesh_mortar.py => mapping/mortar_pair.py} (74%) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index a5ecd1713..277d6d4b4 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -45,14 +45,14 @@ fb = FacetBasis(m, e) Fb = FacetBasis(M, E) -mortar = MeshMortar.init_1D(m, M, - m.boundaries['contact'], - M.facets_satisfying(lambda x: x[0] == 1.0), - np.array([0.0, 1.0])) +mappings = MortarPair.init_2D(m, M, + m.boundaries['contact'], + M.facets_satisfying(lambda x: x[0] == 1.0), + np.array([0.0, 1.0])) mb = [ - MortarBasis(m, e, mapping = mortar[0], intorder=4), - MortarBasis(M, E, mapping = mortar[1], intorder=4) + MortarBasis(m, e, mapping = mappings[0], intorder=4), + MortarBasis(M, E, mapping = mappings[1], intorder=4) ] E1 = 1000.0 diff --git a/skfem/__init__.py b/skfem/__init__.py index a1f30edee..0616e83d9 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -32,8 +32,10 @@ def call(*args, **kwargs): 'InteriorBasis', 'MortarBasis', 'asm', + 'Mapping', 'MappingAffine', 'MappingIsoparametric', + 'MortarPair', 'adaptive_theta', 'bilinear_form', 'build_pc_ilu', diff --git a/skfem/mapping/__init__.py b/skfem/mapping/__init__.py index f40c9e142..c3087c8ea 100644 --- a/skfem/mapping/__init__.py +++ b/skfem/mapping/__init__.py @@ -8,3 +8,4 @@ from .mapping import Mapping from .mapping_affine import MappingAffine from .mapping_isoparametric import MappingIsoparametric +from .mortar_pair import MortarPair diff --git a/skfem/mesh/mesh_mortar.py b/skfem/mapping/mortar_pair.py similarity index 74% rename from skfem/mesh/mesh_mortar.py rename to skfem/mapping/mortar_pair.py index 9c99f9c94..83033788e 100644 --- a/skfem/mesh/mesh_mortar.py +++ b/skfem/mapping/mortar_pair.py @@ -1,49 +1,26 @@ -import numpy as np - -from .mesh import Mesh -from .mesh_line import MeshLine -from ..mapping import MappingAffine - +from typing import NamedTuple -class MeshMortar(Mesh): - """An interface mesh for mortar methods.""" - - name = "Mortar" +import numpy as np +from numpy import ndarray - def __init__(self, m_super, m1, m2, mesh1, mesh2, normals, ix1,ix2,I1,I2): - """Create an interface mesh for mortar methods. +from ..mesh.mesh2d import Mesh2D +from .mapping import Mapping - Parameters - ---------- - mesh1 - mesh2 - p - facets - f2t - normals - """ - self.p = m_super.p - self.t = m_super.t - self.normals = normals - self.supermesh = m_super - self.helper_mesh = 2 * [None] - self.helper_mesh[0] = m1 - self.helper_mesh[1] = m2 - self.target_mesh = 2 * [None] - self.target_mesh[0] = mesh1 - self.target_mesh[1] = mesh2 - self.ix={} - self.ix[0]=ix1 - self.ix[1]=ix2 - self.I ={} - self.I[0]=I1 - self.I[1]=I2 +class MortarPair(NamedTuple): + """A pair of mappings for mortar methods.""" + mapping1: Mapping = None + mapping2: Mapping = None @classmethod - def init_1D(cls, mesh1, mesh2, boundary1, boundary2, tangent): - """Create a mortar mesh between two 2D meshes via projection. + def init_2D(cls, + mesh1: Mesh2D, + mesh2: Mesh2D, + boundary1: ndarray, + boundary2: ndarray, + tangent: ndarray): + """Create mortar mappings for two 2D meshes via projection. Parameters ---------- @@ -57,9 +34,11 @@ def init_1D(cls, mesh1, mesh2, boundary1, boundary2, tangent): A tangent vector defining the direction of the projection. """ + from ..mesh import MeshLine + tangent /= np.linalg.norm(tangent) - # find nodes on the two boundaries + # find unique nodes on the two boundaries p1_ix = np.unique(mesh1.facets[:, boundary1].flatten()) p2_ix = np.unique(mesh2.facets[:, boundary2].flatten()) p1 = mesh1.p[:, p1_ix] @@ -74,21 +53,20 @@ def param(p): y = proj(p) return np.linalg.norm(y, axis=0) * np.sign(np.dot(tangent, y)) + # find unique supermesh facets by combining nodes from both sides param_p1 = param(p1) param_p2 = param(p2) - - # find unique facets by combining nodes from both sides _, ix = np.unique(np.concatenate((param_p1, param_p2)), return_index=True) ixorig = np.concatenate((p1_ix, p2_ix + mesh1.p.shape[1]))[ix] p = np.array([np.hstack((param(mesh1.p), param(mesh2.p)))]) t = np.array([ixorig[:-1], ixorig[1:]]) - # supermesh + # create 1-dimensional supermesh p = p[:, np.concatenate((t[0], np.array([t[1, -1]])))] t = np.array([np.arange(p.shape[1] - 1), np.arange(1, p.shape[1])]) m_super = MeshLine(p, t) - # helper meshes + # helper meshes for creating the mappings m1 = MeshLine(np.sort(param_p1), np.array([np.arange(p1.shape[1] - 1), np.arange(1, p1.shape[1])])) m2 = MeshLine(np.sort(param_p2), np.array([np.arange(p2.shape[1] - 1), @@ -98,14 +76,14 @@ def param(p): normal = np.array([tangent[1], -tangent[0]]) normals = normal[:, None].repeat(t.shape[1], axis=1) - # initialize mappings for orienting + # initialize mappings (for orienting) map_super = m_super.mapping() map_m1 = m1.mapping() map_m2 = m2.mapping() map_mesh1 = mesh1.mapping() map_mesh2 = mesh2.mapping() - # orient helper meshes + # matching of elements in the supermesh and the helper meshes mps = map_super.F(np.array([[0.5]])) ix1 = np.digitize(mps[0, :, 0], m1.p[0]) - 1 ix2 = np.digitize(mps[0, :, 0], m2.p[0]) - 1 @@ -129,25 +107,28 @@ def param(p): ix2_flip = np.unique(ix2[param(z2[:, :, 1]) < param(z2[:, :, 0])]) m2.t[:, ix2_flip] = np.flipud(m2.t[:, ix2_flip]) - - map_m1 = m1.mapping() - map_m2 = m2.mapping() - + # create new mapping objects with G replaced by the supermesh mapping nmap_mesh1 = mesh1.mapping() nmap_mesh2 = mesh2.mapping() - + map_m1 = m1.mapping() + map_m2 = m2.mapping() nmap_mesh1.G = lambda X: map_mesh1.G(map_m1.invF(map_super.F(X), tind=ix1), find=boundary1[sort_boundary1][ix1]) nmap_mesh2.G = lambda X: map_mesh2.G(map_m2.invF(map_super.F(X), tind=ix2), find=boundary2[sort_boundary2][ix2]) + + # these are used by :class:`skfem.assembly.basis.MortarBasis` nmap_mesh1.find = boundary1[sort_boundary1][ix1] nmap_mesh2.find = boundary2[sort_boundary2][ix2] nmap_mesh1.normals = normals nmap_mesh2.normals = normals + # method for calculating the lengths of the segments ('detDG') segments1 = nmap_mesh1.G(np.array([[0.0, 1.0]])) segments2 = nmap_mesh2.G(np.array([[0.0, 1.0]])) - nmap_mesh1.detDG = lambda *_,**__: np.sqrt(np.diff(segments1[0], axis=1)**2 + np.diff(segments1[1], axis=1)**2) - nmap_mesh2.detDG = lambda *_,**__: np.sqrt(np.diff(segments2[0], axis=1)**2 + np.diff(segments2[1], axis=1)**2) + nmap_mesh1.detDG = lambda *_,**__: np.sqrt(np.diff(segments1[0], axis=1)**2 + + np.diff(segments1[1], axis=1)**2) + nmap_mesh2.detDG = lambda *_,**__: np.sqrt(np.diff(segments2[0], axis=1)**2 + + np.diff(segments2[1], axis=1)**2) - return nmap_mesh1, nmap_mesh2 + return cls(mapping1=nmap_mesh1, mapping2=nmap_mesh2) From bd48e03419a6f606c0720b9b2d41897a862f1dfe Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 14:40:48 +0200 Subject: [PATCH 17/32] Remove references to MeshMortar --- skfem/__init__.py | 1 - skfem/mesh/__init__.py | 1 - tests/test_mesh.py | 42 ------------------------------------------ 3 files changed, 44 deletions(-) diff --git a/skfem/__init__.py b/skfem/__init__.py index 0616e83d9..417919490 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -26,7 +26,6 @@ def call(*args, **kwargs): 'MeshQuad', 'MeshTet', 'MeshTri', - 'MeshMortar', 'FacetBasis', 'Basis', 'InteriorBasis', diff --git a/skfem/mesh/__init__.py b/skfem/mesh/__init__.py index 1f85ee0fd..8777fd592 100644 --- a/skfem/mesh/__init__.py +++ b/skfem/mesh/__init__.py @@ -17,4 +17,3 @@ from .mesh_line import MeshLine from .mesh2d import Mesh2D, MeshTri, MeshQuad from .mesh3d import Mesh3D, MeshTet, MeshHex -from .mesh_mortar import MeshMortar diff --git a/tests/test_mesh.py b/tests/test_mesh.py index afa950c8f..9f3753e14 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -129,47 +129,5 @@ def runTest(self): self.assertTrue((m.subdomains[k] == M.subdomains[k]).all()) -class TestMortarMeshTriQuad(unittest.TestCase): - - def create_mesh1(self): - m1 = MeshTri() - m1.refine(3) - m1.translate((0.0, -0.5)) - return m1 - - def create_mesh2(self): - m2 = MeshQuad() - m2.refine(2) - m2.translate((1.0, -0.5)) - return m2 - - def create_meshes(self): - m1 = self.create_mesh1() - m2 = self.create_mesh2() - mortar = MeshMortar.init_1D(m1, m2, - m1.facets_satisfying(lambda x: x[0]==1.0), - m2.facets_satisfying(lambda x: x[0]==1.0), - np.array([0.0, 1.0])) - return m1, m2, mortar - - def runTest(self): - m1, m2, mortar = self.create_meshes() - - # check that f2t indices make sense - self.assertTrue((mortar.f2t[0, :] >= 0).all()) - self.assertTrue((mortar.f2t[0, :] <= np.max(m1.t)).all()) - self.assertTrue((mortar.f2t[1, :] >= 0).all()) - self.assertTrue((mortar.f2t[1, :] <= np.max(m2.t)).all()) - - -class TestMortarMeshTriTri(TestMortarMeshTriQuad): - - def create_mesh2(self): - m2 = MeshTri() - m2.refine(4) - m2.translate((1.0, -0.5)) - return m2 - - if __name__ == '__main__': unittest.main() From 2a839f13f07ac1a31d425c6499d86fde7ea92823 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 14:44:38 +0200 Subject: [PATCH 18/32] Create supermesh in the intersection of projected boundaries --- docs/examples/ex04.py | 77 +++++++++++++++++------------------- skfem/mapping/mortar_pair.py | 3 ++ tests/test_examples.py | 10 +++++ 3 files changed, 49 insertions(+), 41 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index 277d6d4b4..aa81b1f13 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -1,13 +1,13 @@ from skfem import * from skfem.models.elasticity import linear_elasticity -from skfem.visuals.matplotlib import * import numpy as np from pygmsh import generate_mesh from pygmsh.built_in import Geometry from skfem.importers import from_meshio -geom = Geometry() +# create meshes +geom = Geometry() points = [] lines = [] points.append(geom.add_point([0., 0., 0.], .1)) @@ -18,22 +18,14 @@ lines.append(geom.add_line(points[1], points[2])) geom.add_physical(lines[-1], 'dirichlet') geom.add_physical(geom.add_plane_surface(geom.add_line_loop(lines)), 'domain') - m = from_meshio(generate_mesh(geom, dim=2)) - M = MeshLine(np.linspace(0, 1, 6)) * MeshLine(np.linspace(-1, 1, 10)) M.translate((1.0, 0.0)) M.refine() -#M = M._splitquads() -ax = draw(m) -draw(M, ax=ax) -show() - -#m.refine(3) -#M.refine(3) +# define elements and bases e1 = ElementTriP2() e = ElementVectorH1(e1) @@ -55,6 +47,8 @@ MortarBasis(M, E, mapping = mappings[1], intorder=4) ] + +# define bilinear forms E1 = 1000.0 E2 = 1000.0 @@ -76,22 +70,27 @@ alpha = 1000 limit = 0.3 + + +# assemble the stiffness matrices K1 = asm(weakform1, ib) K2 = asm(weakform2, Ib) -K = [[None,None],[None,None]] +K = [[K1, 0.], [0., K2]] f = [None] * 2 def tr(T): return T[0, 0] + T[1, 1] def C(T): - return np.array([[2*Mu*T[0, 0] + Lambda*tr(T), 2*Mu*T[0, 1]], - [2*Mu*T[1, 0], 2*Mu*T[1, 1] + Lambda*tr(T)]]) + return np.array([[2. * Mu * T[0, 0] + Lambda * tr(T), 2. * Mu * T[0, 1]], + [2. * Mu * T[1, 0], 2. * Mu * T[1, 1] + Lambda * tr(T)]]) def Eps(dw): - return np.array([[dw[0][0], 0.5*(dw[0][1] + dw[1][0])], - [0.5*(dw[1][0] + dw[0][1]), dw[1][1]]]) + return np.array([[dw[0][0], .5 * (dw[0][1] + dw[1][0])], + [.5 * (dw[1][0] + dw[0][1]), dw[1][1]]]) + +# assemble the mortar matrices for i in range(2): for j in range(2): @bilinear_form @@ -110,7 +109,7 @@ def bilin_penalty(u, du, v, dv, w): h = w.h return (1. / (alpha * h) * ju * jv - mu * jv - mv * ju) *(np.abs(w.x[1]) <= limit) - K[j][i] = asm(bilin_penalty, mb[i], mb[j]) + K[j][i] += asm(bilin_penalty, mb[i], mb[j]) @linear_form def lin_penalty(v, dv, w): @@ -124,17 +123,14 @@ def lin_penalty(v, dv, w): def gap(x): return (1. - np.sqrt(1. - x[1] ** 2)) return (1. / (alpha * h) * gap(w.x) * jv - gap(w.x) * mv) * (np.abs(w.x[1]) <= limit) - f[i] = asm(lin_penalty, mb[i]) -K[0][0] += K1 -K[1][1] += K2 - -f1 = f[0] -f2 = np.zeros(K2.shape[0]) + f[1] + f[i] = asm(lin_penalty, mb[i]) import scipy.sparse K = (scipy.sparse.bmat(K)).tocsr() + +# set boundary conditions and solve i1 = np.arange(K1.shape[0]) i2 = np.arange(K2.shape[0]) + K1.shape[0] @@ -143,7 +139,7 @@ def gap(x): x = np.zeros(K.shape[0]) -f = np.hstack((f1, f2)) +f = np.hstack((f[0], f[1])) x = np.zeros(K.shape[0]) D = np.concatenate((D1, D2 + ib.N)) @@ -154,15 +150,18 @@ def gap(x): x = solve(*condense(K, f, I=I, x=x)) + +# create a displaced mesh sf = 1 -m.p[0, :] = m.p[0, :] + sf*x[i1][ib.nodal_dofs[0, :]] -m.p[1, :] = m.p[1, :] + sf*x[i1][ib.nodal_dofs[1, :]] +m.p[0, :] = m.p[0, :] + sf * x[i1][ib.nodal_dofs[0, :]] +m.p[1, :] = m.p[1, :] + sf * x[i1][ib.nodal_dofs[1, :]] -M.p[0, :] = M.p[0, :] + sf*x[i2][Ib.nodal_dofs[0, :]] -M.p[1, :] = M.p[1, :] + sf*x[i2][Ib.nodal_dofs[1, :]] +M.p[0, :] = M.p[0, :] + sf * x[i2][Ib.nodal_dofs[0, :]] +M.p[1, :] = M.p[1, :] + sf * x[i2][Ib.nodal_dofs[1, :]] +# post processing s, S = {}, {} e_dg = ElementTriDG(ElementTriP1()) E_dg = ElementQuadDG(ElementQuad1()) @@ -183,8 +182,8 @@ def mass(u, du, v, dv, w): s[itr, jtr] = solve(asm(mass, ib_dg), asm(proj_cauchy, ib, ib_dg) @ x[i1]) S[itr, jtr] = solve(asm(mass, Ib_dg), asm(proj_cauchy, Ib, Ib_dg) @ x[i2]) -s[2, 2] = nu1*(s[0, 0] + s[1, 1]) -S[2, 2] = nu2*(S[0, 0] + S[1, 1]) +s[2, 2] = nu1 * (s[0, 0] + s[1, 1]) +S[2, 2] = nu2 * (S[0, 0] + S[1, 1]) vonmises1 = np.sqrt(.5*((s[0, 0] - s[1, 1])**2 +\ (s[1, 1] - s[2, 2])**2 +\ @@ -196,18 +195,14 @@ def mass(u, du, v, dv, w): (S[2, 2] - S[0, 0])**2 +\ 6.0*S[0, 1]**2)) -ax = plot(ib_dg, vonmises1, shading='gouraud') -draw(m, ax=ax) -plot(Ib_dg, vonmises2, ax=ax, Nrefs=3, shading='gouraud') -draw(M, ax=ax) -show() - -import pdb; pdb.set_trace() if __name__ == "__main__": + from os.path import splitext + from sys import argv + from skfem.visuals.matplotlib import * - from skfem.visuals.matplotlib import draw, show - - ax = draw(m) + ax = plot(ib_dg, vonmises1, shading='gouraud') + draw(m, ax=ax) + plot(Ib_dg, vonmises2, ax=ax, Nrefs=3, shading='gouraud') draw(M, ax=ax) - show() + savefig(splitext(argv[0])[0] + '_vonmises.png') diff --git a/skfem/mapping/mortar_pair.py b/skfem/mapping/mortar_pair.py index 83033788e..ca490499b 100644 --- a/skfem/mapping/mortar_pair.py +++ b/skfem/mapping/mortar_pair.py @@ -63,6 +63,9 @@ def param(p): # create 1-dimensional supermesh p = p[:, np.concatenate((t[0], np.array([t[1, -1]])))] + range_max = np.min([np.max(param_p1), np.max(param_p2)]) + range_min = np.max([np.min(param_p1), np.min(param_p2)]) + p = np.array([p[0, (p[0] <= range_max) * (p[0] >= range_min)]]) t = np.array([np.arange(p.shape[1] - 1), np.arange(1, p.shape[1])]) m_super = MeshLine(p, t) diff --git a/tests/test_examples.py b/tests/test_examples.py index 15f0cc782..f745321f9 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -1,4 +1,5 @@ """These tests run the examples and check that their output stays constant.""" + import unittest import numpy as np @@ -28,6 +29,15 @@ def runTest(self): self.assertAlmostEqual(ex03.L[0], 0.00418289) +class TestEx04(unittest.TestCase): + """Run examples/ex04.py""" + + def runTest(self): + import docs.examples.ex04 as ex04 + self.assertAlmostEqual(np.max(ex04.vonmises1), 64.57919892367978) + self.assertAlmostEqual(np.max(ex04.vonmises2), 67.91419753783893) + + class TestEx05(unittest.TestCase): """Run examples/ex05.py""" From a63cd5fb5d1fdcb60c538f72bf3789886be3118b Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 15:49:30 +0200 Subject: [PATCH 19/32] Avoid regenerating mesh at test suite --- docs/examples/ex04.py | 31 +- docs/examples/ex04_mesh.off | 617 ++++++++++++++++++++++++++++++++++++ skfem/mesh/mesh.py | 24 +- 3 files changed, 655 insertions(+), 17 deletions(-) create mode 100644 docs/examples/ex04_mesh.off diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index aa81b1f13..b7b52b26d 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -4,21 +4,28 @@ from pygmsh import generate_mesh from pygmsh.built_in import Geometry from skfem.importers import from_meshio +import meshio # create meshes -geom = Geometry() -points = [] -lines = [] -points.append(geom.add_point([0., 0., 0.], .1)) -points.append(geom.add_point([0., 1., 0.], .1)) -points.append(geom.add_point([0.,-1., 0.], .1)) -lines.append(geom.add_circle_arc(points[2], points[0], points[1])) -geom.add_physical(lines[-1], 'contact') -lines.append(geom.add_line(points[1], points[2])) -geom.add_physical(lines[-1], 'dirichlet') -geom.add_physical(geom.add_plane_surface(geom.add_line_loop(lines)), 'domain') -m = from_meshio(generate_mesh(geom, dim=2)) +try: + m = MeshTri.load("docs/examples/ex04_mesh.off") + m.boundaries = {'contact': m.facets_satisfying(lambda x: x[0] > 0., + boundaries_only=True)} +except meshio.ReadError: + geom = Geometry() + points = [] + lines = [] + points.append(geom.add_point([0., 0., 0.], .1)) + points.append(geom.add_point([0., 1., 0.], .1)) + points.append(geom.add_point([0.,-1., 0.], .1)) + lines.append(geom.add_circle_arc(points[2], points[0], points[1])) + geom.add_physical(lines[-1], 'contact') + lines.append(geom.add_line(points[1], points[2])) + geom.add_physical(lines[-1], 'dirichlet') + geom.add_physical(geom.add_plane_surface(geom.add_line_loop(lines)), 'domain') + m = from_meshio(generate_mesh(geom, dim=2)) + m.save("docs/examples/ex04_mesh.off") M = MeshLine(np.linspace(0, 1, 6)) * MeshLine(np.linspace(-1, 1, 10)) M.translate((1.0, 0.0)) diff --git a/docs/examples/ex04_mesh.off b/docs/examples/ex04_mesh.off new file mode 100644 index 000000000..8fa111c23 --- /dev/null +++ b/docs/examples/ex04_mesh.off @@ -0,0 +1,617 @@ +OFF +# Created by meshio + +222 390 0 + +0.0 1.0 0.0 +0.0 -1.0 0.0 +0.09801714055591637 -0.9951847266499028 0.0 +0.19509032241883761 -0.9807852803231266 0.0 +0.2902846778931259 -0.9569403355384724 0.0 +0.38268343326549487 -0.9238795321383267 0.0 +0.47139673797238385 -0.8819212637355989 0.0 +0.5555702343496238 -0.8314696114138532 0.0 +0.6343932856511478 -0.7730104521419754 0.0 +0.7071067827656586 -0.7071067796074364 0.0 +0.773010454973133 -0.6343932822013723 0.0 +0.8314696138527624 -0.5555702306995381 0.0 +0.8819212658270285 -0.47139673405959426 0.0 +0.9238795338254099 -0.38268342919251586 0.0 +0.9569403368184075 -0.2902846736737454 0.0 +0.9807852811883714 -0.1950903180689581 0.0 +0.995184727094887 -0.09801713603791533 0.0 +1.0 4.5837337757057685e-09 0.0 +0.9951847262476915 0.09801714463963587 0.0 +0.980785279615548 0.19509032597607545 0.0 +0.9569403346479524 0.290284680828777 0.0 +0.9238795311905263 0.38268343555368756 0.0 +0.8819212628015901 0.4713967397197912 0.0 +0.8314696106349614 0.5555702355153178 0.0 +0.7730104516368775 0.6343932862666112 0.0 +0.7071067794805702 0.7071067828925248 0.0 +0.634393282512301 0.7730104547179604 0.0 +0.5555702315457777 0.8314696132873234 0.0 +0.4713967355549627 0.8819212650277374 0.0 +0.38268343125990784 0.9238795329690681 0.0 +0.2902846763567968 0.9569403360045128 0.0 +0.195090321428782 0.9807852805200609 0.0 +0.09801714007694617 0.9951847266970772 0.0 +0.0 0.9000000000002774 0.0 +0.0 0.8000000000005548 0.0 +0.0 0.7000000000008322 0.0 +0.0 0.6000000000011096 0.0 +0.0 0.500000000001387 0.0 +0.0 0.40000000000166436 0.0 +0.0 0.3000000000019418 0.0 +0.0 0.20000000000221918 0.0 +0.0 0.10000000000249643 0.0 +0.0 2.752242878045763e-12 0.0 +0.0 -0.09999999999750342 0.0 +0.0 -0.19999999999778084 0.0 +0.0 -0.29999999999805826 0.0 +0.0 -0.39999999999833546 0.0 +0.0 -0.4999999999986129 0.0 +0.0 -0.5999999999988903 0.0 +0.0 -0.699999999999168 0.0 +0.0 -0.7999999999994454 0.0 +0.0 -0.8999999999997226 0.0 +0.9159764153258458 0.05383987386802238 0.0 +0.8314138091584239 -0.39785805095354165 0.0 +0.8264043258240841 0.39254386345989056 0.0 +0.08657154386342972 -0.04993181939059333 0.0 +0.08825881533998217 0.24854065189332641 0.0 +0.0866467339337712 -0.350125666409955 0.0 +0.6151960220089083 0.6914850111403807 0.0 +0.611656851878296 -0.6767753615474275 0.0 +0.08792437843765298 0.546553776392746 0.0 +0.40214879730821057 -0.8329018889819535 0.0 +0.0864975927752666 -0.5501598487639185 0.0 +0.4031239663030456 0.8299929282327413 0.0 +0.08552443345288481 0.7412874524475859 0.0 +0.08387010326216801 -0.7452678551627052 0.0 +0.9064888004049059 -0.13475830751070036 0.0 +0.883006687779775 0.21591742664103838 0.0 +0.7318295262956729 0.5462376817273685 0.0 +0.7310857439949313 -0.5534884214698972 0.0 +0.22203717119110486 0.8864220786001916 0.0 +0.222037172228374 -0.88642207782408 0.0 +0.12735200151713993 0.8953724394140379 0.0 +0.1264872216201089 -0.8973542598638573 0.0 +0.08603181563097448 -0.6492168486242395 0.0 +0.1728569295463326 -0.5999642076514183 0.0 +0.1729345952627196 -0.500480145586956 0.0 +0.2591331121034306 -0.5504140700880753 0.0 +0.2592714087128541 -0.4510637781789997 0.0 +0.3451538888162647 -0.5008525665656047 0.0 +0.3453660173120208 -0.4016689274236115 0.0 +0.43112082579940325 -0.45140415758047114 0.0 +0.4313132869549403 -0.35228060068609146 0.0 +0.3456417826557642 -0.30202707863747263 0.0 +0.4316914876509484 -0.25313887750341396 0.0 +0.5172439124243171 -0.30297215045699327 0.0 +0.5175060234485399 -0.20398089749625326 0.0 +0.4320078005036593 -0.1541618871667654 0.0 +0.5178249360980043 -0.10497328859987322 0.0 +0.4323517746522692 -0.05513829699300258 0.0 +0.5181810627769347 -0.005926661553542319 0.0 +0.4326873030792038 0.0439068356495514 0.0 +0.5185129534634814 0.0931212502063908 0.0 +0.43303358988383006 0.14295033231019513 0.0 +0.5188407112939315 0.19216073999298827 0.0 +0.4333907822760996 0.2419806037163972 0.0 +0.5191570773658916 0.29117889053487644 0.0 +0.4336167327285912 0.34098097416379747 0.0 +0.3466297470498466 0.09502752487561306 0.0 +0.6036466835467159 -0.05588497882022823 0.0 +0.6043483291878866 0.14226361071690696 0.0 +0.6030083111892219 -0.253739769204638 0.0 +0.5194358558074463 0.3901995369490897 0.0 +0.4339131655358476 0.44001601280296204 0.0 +0.34799605056576544 0.3907547089718833 0.0 +0.3483736150587413 0.4897859002884928 0.0 +0.3461095808402675 -0.10314722872110887 0.0 +0.3451380056642011 -0.5999795967606039 0.0 +0.43423779016255415 0.5389955255642274 0.0 +0.3487855783748346 0.5887556881868931 0.0 +0.6027356100339158 -0.35253659554114325 0.0 +0.6048915795006872 0.34040851731280114 0.0 +0.26231771840506496 0.4405891472347052 0.0 +0.26200491468128767 0.34136888541754196 0.0 +0.4341509042932531 0.6378572118248705 0.0 +0.17613909944252154 0.39128549747516345 0.0 +0.6867867143117476 -0.3005776176954207 0.0 +0.6883938442262401 -0.20422841407407338 0.0 +0.7887504482032948 -0.24363163143334435 0.0 +0.6039194807805792 0.4371233138451532 0.0 +0.34717469685206775 0.6842928036187731 0.0 +0.2630163234652642 0.6380173006925967 0.0 +0.26231461886291657 0.7347372834745751 0.0 +0.6821568960369653 -0.38932236004593684 0.0 +0.6024613510918597 -0.4500576975192292 0.0 +0.7718417301410816 -0.14476793206076038 0.0 +0.6840585160728114 0.38136305112776203 0.0 +0.6872609976885783 0.28791528538360167 0.0 +0.5168363646187115 0.5833344105436389 0.0 +0.17550407052083822 0.6900813696897046 0.0 +0.17619557798382957 0.4914856785795687 0.0 +0.42016802198142106 0.7346497809142203 0.0 +0.17302511054291092 -0.40062265779056283 0.0 +0.17237477181374233 -0.6988080454878537 0.0 +0.1714952098378757 -0.8012321435395398 0.0 +0.25599371682296734 -0.7410015264285403 0.0 +0.5190057914509976 0.4879530306260381 0.0 +0.618290339602385 0.536875949487155 0.0 +0.4308956703486012 -0.5504811711457109 0.0 +0.43130650380920504 -0.6485484160152994 0.0 +0.5175687390320705 -0.6021217093581201 0.0 +0.5092206012248724 -0.7035085970951118 0.0 +0.620737462322611 -0.555765633253074 0.0 +0.3429556442908701 -0.6929402428430802 0.0 +0.25848406726552225 -0.6479093927842566 0.0 +0.6885382570987393 -0.10392056646537914 0.0 +0.769286205655087 -0.05057534362618522 0.0 +0.6040111157907584 0.043242214936212164 0.0 +0.6897389251565653 0.09246051544826223 0.0 +0.6895451688254398 0.1909341621919352 0.0 +0.7853297010802 0.13990094156251878 0.0 +0.7756810503792961 0.2413610566339871 0.0 +0.7701883492540708 0.03853734593193141 0.0 +0.5190243392045759 -0.5017121515376488 0.0 +0.6040073106437357 0.2408102010221849 0.0 +0.687568256083656 -0.006023468765897795 0.0 +0.08680748951240319 0.14978745087550502 0.0 +0.17485631474456298 0.19832056843973647 0.0 +0.1730517611561072 0.10005310843209428 0.0 +0.26087327445847697 0.14679457604766757 0.0 +0.08657886954653793 0.04998478457738213 0.0 +0.17291704371860217 -0.00010105186068302731 0.0 +0.1729970659840663 -0.10002826078482492 0.0 +0.2593646476016979 -0.04985929515148815 0.0 +0.08660511306657999 -0.1498281482654842 0.0 +0.17305207548825177 -0.20001247954407836 0.0 +0.08658138073658056 -0.24999061371589165 0.0 +0.17306526794695942 -0.3003972076357125 0.0 +0.2594694617813657 -0.25019468739155504 0.0 +0.517316554251502 -0.4018272255535962 0.0 +0.08648158353870185 -0.4502095033908465 0.0 +0.6031530092679103 -0.1544546524434075 0.0 +0.26271089280937426 0.539695638559989 0.0 +0.25930650815864587 -0.35099572284298564 0.0 +0.3481417000647074 0.29225483630730303 0.0 +0.08775490432998115 0.6445771191516094 0.0 +0.17312433591087684 0.7949143361266013 0.0 +0.08917340860778944 0.4445319633762133 0.0 +0.17610593477516273 0.5911079955205527 0.0 +0.34591484982492143 -0.20259666030335644 0.0 +0.2596845571558199 0.04809501117136811 0.0 +0.259197333421417 -0.15100131147531112 0.0 +0.34645554142539703 -0.004036758239607918 0.0 +0.3476513202248709 0.19411470667553155 0.0 +0.5165083519158323 0.6874008844465628 0.0 +0.26285614408485514 0.24520073416291943 0.0 +0.8702436047190465 0.12724803991148598 0.0 +0.1765715058218143 0.29520431385172685 0.0 +0.08940011612821333 0.34630648225343424 0.0 +0.4166298503101692 -0.7429585203863014 0.0 +0.8549693515455798 0.30930121195949895 0.0 +0.8950418208143169 -0.2241963012461981 0.0 +0.860025803539879 -0.3124749658989302 0.0 +0.6724940254177364 0.6208040077969013 0.0 +0.7854994399143052 0.4754049055547044 0.0 +0.9215823721068173 -0.045274434635522595 0.0 +0.785420238160406 -0.4789576989502909 0.0 +0.31501693368973693 0.863121801168385 0.0 +0.31408164224217594 -0.8662391964201028 0.0 +0.688719459186926 -0.6255058956158415 0.0 +0.844562532980393 -0.007366724481464558 0.0 +0.6969550069593997 0.46476187067423613 0.0 +0.7641662461091518 0.3352234361081791 0.0 +0.766062195642033 -0.3394573773342276 0.0 +0.6967061893284834 -0.4732097765718286 0.0 +0.48759115339420067 0.7864769275402721 0.0 +0.33363382533391917 0.7769765419565695 0.0 +0.4861251931921746 -0.7917065301079622 0.0 +0.5603470596005116 0.7537033611358295 0.0 +0.25400471431432753 0.8150646543307231 0.0 +0.3311957483405575 -0.7829906385661866 0.0 +0.8445851642438251 -0.08836500410921283 0.0 +0.5617849379877878 -0.7574792990349553 0.0 +0.07305726944797225 0.8238422029156591 0.0 +0.07358418293812316 -0.8256623870972735 0.0 +0.5883410243610833 0.6196215619927677 0.0 +0.2566752832473187 -0.8195973033555625 0.0 +0.841341592761485 -0.16714383527204318 0.0 +0.7514167069759504 0.40985942538495446 0.0 +0.7553021899259311 -0.4137557280795954 0.0 +0.8372601206719111 0.0704318953584988 0.0 +0.9360577467536026 0.13890639314725461 0.0 +3 32 33 72 +3 2 51 73 +3 67 150 151 +3 67 150 186 +3 0 32 33 +3 1 2 51 +3 124 142 153 +3 31 70 72 +3 3 71 73 +3 31 32 72 +3 2 3 73 +3 68 137 193 +3 70 72 176 +3 59 140 142 +3 118 191 192 +3 68 137 201 +3 140 142 153 +3 20 21 190 +3 59 140 141 +3 118 191 217 +3 124 142 204 +3 71 73 134 +3 72 176 213 +3 69 142 204 +3 131 184 205 +3 184 205 208 +3 37 60 177 +3 60 130 177 +3 69 142 199 +3 20 67 190 +3 21 54 190 +3 64 129 175 +3 64 129 176 +3 118 192 203 +3 129 175 178 +3 38 39 188 +3 60 175 178 +3 196 204 219 +3 37 38 177 +3 60 130 178 +3 98 159 180 +3 67 151 190 +3 167 168 173 +3 83 168 173 +3 106 163 181 +3 98 159 183 +3 98 180 182 +3 38 177 188 +3 163 180 182 +3 69 196 204 +3 57 132 167 +3 132 167 173 +3 162 163 181 +3 146 152 155 +3 158 159 180 +3 106 163 182 +3 127 149 151 +3 127 149 154 +3 135 143 144 +3 106 179 181 +3 75 133 144 +3 133 135 144 +3 57 132 170 +3 165 167 168 +3 99 145 155 +3 145 146 155 +3 45 57 166 +3 124 153 169 +3 57 166 167 +3 35 64 175 +3 165 166 167 +3 168 179 181 +3 65 133 134 +3 158 160 161 +3 83 168 179 +3 102 119 136 +3 128 136 137 +3 62 74 75 +3 99 145 171 +3 119 136 137 +3 133 134 135 +3 75 77 144 +3 65 74 133 +3 100 147 148 +3 96 111 154 +3 48 62 74 +3 107 138 139 +3 111 127 154 +3 148 152 155 +3 147 148 155 +3 74 75 133 +3 92 100 147 +3 44 45 166 +3 81 153 169 +3 102 103 136 +3 79 107 138 +3 88 99 171 +3 46 57 170 +3 156 158 160 +3 100 149 154 +3 138 140 153 +3 62 76 170 +3 41 156 160 +3 76 132 170 +3 138 139 140 +3 77 79 107 +3 139 140 141 +3 75 76 77 +3 62 75 76 +3 43 44 164 +3 81 138 153 +3 110 124 169 +3 49 65 74 +3 77 78 79 +3 114 120 131 +3 79 81 138 +3 100 148 149 +3 94 96 154 +3 108 114 128 +3 56 156 157 +3 85 110 169 +3 111 119 126 +3 161 162 163 +3 157 158 159 +3 148 149 150 +3 34 35 64 +3 101 110 116 +3 49 50 65 +3 55 160 161 +3 96 102 111 +3 102 111 119 +3 110 116 123 +3 108 109 114 +3 79 80 81 +3 99 147 155 +3 44 164 166 +3 85 101 110 +3 90 92 147 +3 92 94 100 +3 85 86 101 +3 158 161 180 +3 116 117 118 +3 47 48 62 +3 149 150 151 +3 88 90 99 +3 94 100 154 +3 156 157 158 +3 55 161 162 +3 103 108 136 +3 82 84 85 +3 76 77 78 +3 107 139 143 +3 90 99 147 +3 36 37 60 +3 40 56 156 +3 148 150 152 +3 42 55 160 +3 78 79 80 +3 109 114 120 +3 117 118 125 +3 101 116 117 +3 81 82 169 +3 48 49 74 +3 84 85 86 +3 39 40 56 +3 111 126 127 +3 96 97 102 +3 103 105 108 +3 45 46 57 +3 80 81 82 +3 42 43 55 +3 101 117 171 +3 94 95 96 +3 162 164 165 +3 117 145 171 +3 108 128 136 +3 41 42 160 +3 92 93 94 +3 82 85 169 +3 40 41 156 +3 90 91 92 +3 43 55 164 +3 105 108 109 +3 110 123 124 +3 97 102 103 +3 47 62 170 +3 88 89 90 +3 86 87 88 +3 125 145 146 +3 86 88 171 +3 77 107 144 +3 82 83 84 +3 46 47 170 +3 95 96 97 +3 162 165 181 +3 55 162 164 +3 86 101 171 +3 93 94 95 +3 164 165 166 +3 64 176 213 +3 120 121 122 +3 91 92 93 +3 117 125 145 +3 89 90 91 +3 84 86 87 +3 87 88 89 +3 15 16 66 +3 19 20 67 +3 103 104 105 +3 80 82 83 +3 23 24 68 +3 76 78 132 +3 21 22 54 +3 17 18 52 +3 10 11 69 +3 12 13 53 +3 25 26 58 +3 35 36 175 +3 97 103 104 +3 8 9 59 +3 109 120 121 +3 137 193 215 +3 107 143 144 +3 91 93 98 +3 112 115 130 +3 28 29 63 +3 87 89 106 +3 5 6 61 +3 104 112 113 +3 122 129 176 +3 112 130 172 +3 105 109 172 +3 78 80 173 +3 104 105 112 +3 78 132 173 +3 80 83 173 +3 95 97 174 +3 112 113 115 +3 36 60 175 +3 97 104 174 +3 121 122 129 +3 109 121 172 +3 30 31 70 +3 3 4 71 +3 104 113 174 +3 105 112 172 +3 83 84 179 +3 115 130 177 +3 39 56 188 +3 87 106 179 +3 91 98 182 +3 89 106 182 +3 56 157 187 +3 93 98 183 +3 121 172 178 +3 95 174 183 +3 165 168 181 +3 157 159 185 +3 84 87 179 +3 121 129 178 +3 130 172 178 +3 161 163 180 +3 33 72 213 +3 51 73 214 +3 89 91 182 +3 114 128 184 +3 93 95 183 +3 56 187 188 +3 157 185 187 +3 141 189 207 +3 73 134 214 +3 174 183 185 +3 114 131 184 +3 141 207 212 +3 139 143 189 +3 159 183 185 +3 68 194 201 +3 143 189 210 +3 194 201 218 +3 113 185 187 +3 113 174 185 +3 113 115 187 +3 53 192 203 +3 115 187 188 +3 115 177 188 +3 59 142 199 +3 197 206 209 +3 118 125 217 +3 119 137 201 +3 139 141 189 +3 128 184 215 +3 61 189 210 +3 63 131 206 +3 116 123 203 +3 58 184 215 +3 120 131 206 +3 65 134 214 +3 63 197 206 +3 33 34 213 +3 50 51 214 +3 125 146 211 +3 70 176 209 +3 14 15 191 +3 16 17 195 +3 22 23 194 +3 13 14 192 +3 24 25 193 +3 11 12 196 +3 15 66 191 +3 16 66 195 +3 23 68 194 +3 22 54 194 +3 24 68 193 +3 17 52 195 +3 11 69 196 +3 13 53 192 +3 12 53 196 +3 25 58 193 +3 29 30 197 +3 4 5 198 +3 150 152 220 +3 29 63 197 +3 5 61 198 +3 123 204 219 +3 30 70 197 +3 4 71 198 +3 127 151 202 +3 10 69 199 +3 9 59 199 +3 126 127 202 +3 150 186 220 +3 63 131 205 +3 59 141 212 +3 28 63 205 +3 146 152 200 +3 26 58 208 +3 6 61 207 +3 9 10 199 +3 19 67 221 +3 18 52 221 +3 8 59 212 +3 119 126 201 +3 14 191 192 +3 58 184 208 +3 125 211 217 +3 151 190 202 +3 67 186 221 +3 52 186 221 +3 61 189 207 +3 116 118 203 +3 123 124 204 +3 61 198 210 +3 71 134 216 +3 27 28 205 +3 120 122 206 +3 134 135 216 +3 52 195 200 +3 122 176 209 +3 135 143 210 +3 58 193 215 +3 126 201 218 +3 198 210 216 +3 123 203 219 +3 6 7 207 +3 26 27 208 +3 54 190 202 +3 122 206 209 +3 128 137 215 +3 7 8 212 +3 34 64 213 +3 50 65 214 +3 18 19 221 +3 66 195 211 +3 27 205 208 +3 7 207 212 +3 70 197 209 +3 52 186 220 +3 195 200 211 +3 126 202 218 +3 66 211 217 +3 135 210 216 +3 66 191 217 +3 71 198 216 +3 146 200 211 +3 54 194 218 +3 53 196 219 +3 54 202 218 +3 52 200 220 +3 152 200 220 +3 53 203 219 diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py index 00a00d36a..46eddc728 100644 --- a/skfem/mesh/mesh.py +++ b/skfem/mesh/mesh.py @@ -312,7 +312,9 @@ def element_finder(self) -> Callable[[ndarray], ndarray]: raise NotImplementedError("element_finder not implemented " "for the given Mesh type.") - def nodes_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray: + def nodes_satisfying(self, + test: Callable[[ndarray], bool], + boundaries_only: bool = False) -> ndarray: """Return nodes that satisfy some condition. Parameters @@ -320,11 +322,18 @@ def nodes_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray: test A function which returns True for the set of nodes that are to be included in the return set. + boundaries_only + If True, include only boundary facets. """ - return np.nonzero(test(self.p))[0] - - def facets_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray: + nodes = np.nonzero(test(self.p))[0] + if boundaries_only: + nodes = np.intersect1d(nodes, self.boundary_nodes()) + return nodes + + def facets_satisfying(self, + test: Callable[[ndarray], bool], + boundaries_only: bool = False) -> ndarray: """Return facets whose midpoints satisfy some condition. Parameters @@ -332,11 +341,16 @@ def facets_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray: test A function which returns True for the facet midpoints that are to be included in the return set. + boundaries_only + If True, include only boundary facets. """ midp = [np.sum(self.p[itr, self.facets], axis=0) / self.facets.shape[0] for itr in range(self.p.shape[0])] - return np.nonzero(test(np.array(midp)))[0] + facets = np.nonzero(test(np.array(midp)))[0] + if boundaries_only: + facets = np.intersect1d(facets, self.boundary_facets()) + return facets def elements_satisfying(self, test: Callable[[ndarray], bool]) -> ndarray: From f7329380d4b714588d1496f36e449aabc97a4a86 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 16:43:30 +0200 Subject: [PATCH 20/32] Add io for .json files --- docs/examples/ex04.py | 9 ++++----- skfem/importers/json.py | 34 ++++++++++++++++++++++++++++++++++ skfem/importers/meshio.py | 1 + skfem/mesh/mesh.py | 2 +- 4 files changed, 40 insertions(+), 6 deletions(-) create mode 100644 skfem/importers/json.py diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index b7b52b26d..7499c76f8 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -4,15 +4,14 @@ from pygmsh import generate_mesh from pygmsh.built_in import Geometry from skfem.importers import from_meshio +from skfem.importers.json import from_file, to_file import meshio # create meshes try: - m = MeshTri.load("docs/examples/ex04_mesh.off") - m.boundaries = {'contact': m.facets_satisfying(lambda x: x[0] > 0., - boundaries_only=True)} -except meshio.ReadError: + m = from_file("docs/examples/ex04_mesh.json") +except FileNotFoundError: geom = Geometry() points = [] lines = [] @@ -25,7 +24,7 @@ geom.add_physical(lines[-1], 'dirichlet') geom.add_physical(geom.add_plane_surface(geom.add_line_loop(lines)), 'domain') m = from_meshio(generate_mesh(geom, dim=2)) - m.save("docs/examples/ex04_mesh.off") + to_file(m, "docs/examples/ex04_mesh.json") M = MeshLine(np.linspace(0, 1, 6)) * MeshLine(np.linspace(-1, 1, 10)) M.translate((1.0, 0.0)) diff --git a/skfem/importers/json.py b/skfem/importers/json.py new file mode 100644 index 000000000..0d7fb337e --- /dev/null +++ b/skfem/importers/json.py @@ -0,0 +1,34 @@ +"""Import mesh from JSON as defined by :class:`skfem.mesh.to_dict`.""" + +import json + +from skfem.mesh import * + + +def from_file(filename: str): + with open(filename, 'r') as handle: + d = json.load(handle) + + # detect dimension and number of vertices + dim = len(d['p'][0]) + nverts = len(d['t'][0]) + + if dim == 1: + mesh_type = MeshLine + elif dim == 2 and nverts == 3: + mesh_type = MeshTri + elif dim == 2 and nverts == 4: + mesh_type = MeshQuad + elif dim == 3 and nverts == 4: + mesh_type = MeshTet + elif dim == 3 and nverts == 8: + mesh_type = MeshHex + else: + raise NotImplementedError("The given mesh is not supported.") + + return mesh_type.from_dict(d) + + +def to_file(mesh: Mesh, filename: str): + with open(filename, 'w') as handle: + d = json.dump(mesh.to_dict(), handle) diff --git a/skfem/importers/meshio.py b/skfem/importers/meshio.py index 032d350d6..977a13e5d 100644 --- a/skfem/importers/meshio.py +++ b/skfem/importers/meshio.py @@ -8,6 +8,7 @@ import skfem + MESH_TYPE_MAPPING = OrderedDict([ ('tetra', skfem.MeshTet), ('hexahedron', skfem.MeshHex), diff --git a/skfem/mesh/mesh.py b/skfem/mesh/mesh.py index 46eddc728..3db862b4c 100644 --- a/skfem/mesh/mesh.py +++ b/skfem/mesh/mesh.py @@ -371,7 +371,7 @@ def elements_satisfying(self, def from_dict(cls: Type[MeshType], d) -> MeshType: """Initialize a mesh from a dictionary.""" if 'p' not in d or 't' not in d: - raise ValueError("Dictionary must contain keys 'p' and 't'") + raise ValueError("Dictionary must contain keys 'p' and 't'.") else: d['p'] = np.array(d['p']).T d['t'] = np.array(d['t']).T From 5dfa746c4625f40e0e3ea0650f849ecf6a944147 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 16:43:57 +0200 Subject: [PATCH 21/32] Add the ex04_mesh.json --- docs/examples/ex04_mesh.json | 1 + docs/examples/ex04_mesh.off | 617 ----------------------------------- 2 files changed, 1 insertion(+), 617 deletions(-) create mode 100644 docs/examples/ex04_mesh.json delete mode 100644 docs/examples/ex04_mesh.off diff --git a/docs/examples/ex04_mesh.json b/docs/examples/ex04_mesh.json new file mode 100644 index 000000000..8064125bd --- /dev/null +++ b/docs/examples/ex04_mesh.json @@ -0,0 +1 @@ +{"p": [[0.0, 1.0], [0.0, -1.0], [0.09801714055591637, -0.9951847266499028], [0.19509032241883761, -0.9807852803231266], [0.2902846778931259, -0.9569403355384724], [0.38268343326549487, -0.9238795321383267], [0.47139673797238385, -0.8819212637355989], [0.5555702343496238, -0.8314696114138532], [0.6343932856511478, -0.7730104521419754], [0.7071067827656586, -0.7071067796074364], [0.773010454973133, -0.6343932822013723], [0.8314696138527624, -0.5555702306995381], [0.8819212658270285, -0.47139673405959426], [0.9238795338254099, -0.38268342919251586], [0.9569403368184075, -0.2902846736737454], [0.9807852811883714, -0.1950903180689581], [0.995184727094887, -0.09801713603791533], [1.0, 4.5837337757057685e-09], [0.9951847262476915, 0.09801714463963587], [0.980785279615548, 0.19509032597607545], [0.9569403346479524, 0.290284680828777], [0.9238795311905263, 0.38268343555368756], [0.8819212628015901, 0.4713967397197912], [0.8314696106349614, 0.5555702355153178], [0.7730104516368775, 0.6343932862666112], [0.7071067794805702, 0.7071067828925248], [0.634393282512301, 0.7730104547179604], [0.5555702315457777, 0.8314696132873234], [0.4713967355549627, 0.8819212650277374], [0.38268343125990784, 0.9238795329690681], [0.2902846763567968, 0.9569403360045128], [0.195090321428782, 0.9807852805200609], [0.09801714007694617, 0.9951847266970772], [0.0, 0.9000000000002774], [0.0, 0.8000000000005548], [0.0, 0.7000000000008322], [0.0, 0.6000000000011096], [0.0, 0.500000000001387], [0.0, 0.40000000000166436], [0.0, 0.3000000000019418], [0.0, 0.20000000000221918], [0.0, 0.10000000000249643], [0.0, 2.752242878045763e-12], [0.0, -0.09999999999750342], [0.0, -0.19999999999778084], [0.0, -0.29999999999805826], [0.0, -0.39999999999833546], [0.0, -0.4999999999986129], [0.0, -0.5999999999988903], [0.0, -0.699999999999168], [0.0, -0.7999999999994454], [0.0, -0.8999999999997226], [0.9159764153258458, 0.05383987386802238], [0.8314138091584239, -0.39785805095354165], [0.8264043258240841, 0.39254386345989056], [0.08657154386342972, -0.04993181939059333], [0.08825881533998217, 0.24854065189332641], [0.0866467339337712, -0.350125666409955], [0.6151960220089083, 0.6914850111403807], [0.611656851878296, -0.6767753615474275], [0.08792437843765298, 0.546553776392746], [0.40214879730821057, -0.8329018889819535], [0.0864975927752666, -0.5501598487639185], [0.4031239663030456, 0.8299929282327413], [0.08552443345288481, 0.7412874524475859], [0.08387010326216801, -0.7452678551627052], [0.9064888004049059, -0.13475830751070036], [0.883006687779775, 0.21591742664103838], [0.7318295262956729, 0.5462376817273685], [0.7310857439949313, -0.5534884214698972], [0.22203717119110486, 0.8864220786001916], [0.222037172228374, -0.88642207782408], [0.12735200151713993, 0.8953724394140379], [0.1264872216201089, -0.8973542598638573], [0.08603181563097448, -0.6492168486242395], [0.1728569295463326, -0.5999642076514183], [0.1729345952627196, -0.500480145586956], [0.2591331121034306, -0.5504140700880753], [0.2592714087128541, -0.4510637781789997], [0.3451538888162647, -0.5008525665656047], [0.3453660173120208, -0.4016689274236115], [0.43112082579940325, -0.45140415758047114], [0.4313132869549403, -0.35228060068609146], [0.3456417826557642, -0.30202707863747263], [0.4316914876509484, -0.25313887750341396], [0.5172439124243171, -0.30297215045699327], [0.5175060234485399, -0.20398089749625326], [0.4320078005036593, -0.1541618871667654], [0.5178249360980043, -0.10497328859987322], [0.4323517746522692, -0.05513829699300258], [0.5181810627769347, -0.005926661553542319], [0.4326873030792038, 0.0439068356495514], [0.5185129534634814, 0.0931212502063908], [0.43303358988383006, 0.14295033231019513], [0.5188407112939315, 0.19216073999298827], [0.4333907822760996, 0.2419806037163972], [0.5191570773658916, 0.29117889053487644], [0.4336167327285912, 0.34098097416379747], [0.3466297470498466, 0.09502752487561306], [0.6036466835467159, -0.05588497882022823], [0.6043483291878866, 0.14226361071690696], [0.6030083111892219, -0.253739769204638], [0.5194358558074463, 0.3901995369490897], [0.4339131655358476, 0.44001601280296204], [0.34799605056576544, 0.3907547089718833], [0.3483736150587413, 0.4897859002884928], [0.3461095808402675, -0.10314722872110887], [0.3451380056642011, -0.5999795967606039], [0.43423779016255415, 0.5389955255642274], [0.3487855783748346, 0.5887556881868931], [0.6027356100339158, -0.35253659554114325], [0.6048915795006872, 0.34040851731280114], [0.26231771840506496, 0.4405891472347052], [0.26200491468128767, 0.34136888541754196], [0.4341509042932531, 0.6378572118248705], [0.17613909944252154, 0.39128549747516345], [0.6867867143117476, -0.3005776176954207], [0.6883938442262401, -0.20422841407407338], [0.7887504482032948, -0.24363163143334435], [0.6039194807805792, 0.4371233138451532], [0.34717469685206775, 0.6842928036187731], [0.2630163234652642, 0.6380173006925967], [0.26231461886291657, 0.7347372834745751], [0.6821568960369653, -0.38932236004593684], [0.6024613510918597, -0.4500576975192292], [0.7718417301410816, -0.14476793206076038], [0.6840585160728114, 0.38136305112776203], [0.6872609976885783, 0.28791528538360167], [0.5168363646187115, 0.5833344105436389], [0.17550407052083822, 0.6900813696897046], [0.17619557798382957, 0.4914856785795687], [0.42016802198142106, 0.7346497809142203], [0.17302511054291092, -0.40062265779056283], [0.17237477181374233, -0.6988080454878537], [0.1714952098378757, -0.8012321435395398], [0.25599371682296734, -0.7410015264285403], [0.5190057914509976, 0.4879530306260381], [0.618290339602385, 0.536875949487155], [0.4308956703486012, -0.5504811711457109], [0.43130650380920504, -0.6485484160152994], [0.5175687390320705, -0.6021217093581201], [0.5092206012248724, -0.7035085970951118], [0.620737462322611, -0.555765633253074], [0.3429556442908701, -0.6929402428430802], [0.25848406726552225, -0.6479093927842566], [0.6885382570987393, -0.10392056646537914], [0.769286205655087, -0.05057534362618522], [0.6040111157907584, 0.043242214936212164], [0.6897389251565653, 0.09246051544826223], [0.6895451688254398, 0.1909341621919352], [0.7853297010802, 0.13990094156251878], [0.7756810503792961, 0.2413610566339871], [0.7701883492540708, 0.03853734593193141], [0.5190243392045759, -0.5017121515376488], [0.6040073106437357, 0.2408102010221849], [0.687568256083656, -0.006023468765897795], [0.08680748951240319, 0.14978745087550502], [0.17485631474456298, 0.19832056843973647], [0.1730517611561072, 0.10005310843209428], [0.26087327445847697, 0.14679457604766757], [0.08657886954653793, 0.04998478457738213], [0.17291704371860217, -0.00010105186068302731], [0.1729970659840663, -0.10002826078482492], [0.2593646476016979, -0.04985929515148815], [0.08660511306657999, -0.1498281482654842], [0.17305207548825177, -0.20001247954407836], [0.08658138073658056, -0.24999061371589165], [0.17306526794695942, -0.3003972076357125], [0.2594694617813657, -0.25019468739155504], [0.517316554251502, -0.4018272255535962], [0.08648158353870185, -0.4502095033908465], [0.6031530092679103, -0.1544546524434075], [0.26271089280937426, 0.539695638559989], [0.25930650815864587, -0.35099572284298564], [0.3481417000647074, 0.29225483630730303], [0.08775490432998115, 0.6445771191516094], [0.17312433591087684, 0.7949143361266013], [0.08917340860778944, 0.4445319633762133], [0.17610593477516273, 0.5911079955205527], [0.34591484982492143, -0.20259666030335644], [0.2596845571558199, 0.04809501117136811], [0.259197333421417, -0.15100131147531112], [0.34645554142539703, -0.004036758239607918], [0.3476513202248709, 0.19411470667553155], [0.5165083519158323, 0.6874008844465628], [0.26285614408485514, 0.24520073416291943], [0.8702436047190465, 0.12724803991148598], [0.1765715058218143, 0.29520431385172685], [0.08940011612821333, 0.34630648225343424], [0.4166298503101692, -0.7429585203863014], [0.8549693515455798, 0.30930121195949895], [0.8950418208143169, -0.2241963012461981], [0.860025803539879, -0.3124749658989302], [0.6724940254177364, 0.6208040077969013], [0.7854994399143052, 0.4754049055547044], [0.9215823721068173, -0.045274434635522595], [0.785420238160406, -0.4789576989502909], [0.31501693368973693, 0.863121801168385], [0.31408164224217594, -0.8662391964201028], [0.688719459186926, -0.6255058956158415], [0.844562532980393, -0.007366724481464558], [0.6969550069593997, 0.46476187067423613], [0.7641662461091518, 0.3352234361081791], [0.766062195642033, -0.3394573773342276], [0.6967061893284834, -0.4732097765718286], [0.48759115339420067, 0.7864769275402721], [0.33363382533391917, 0.7769765419565695], [0.4861251931921746, -0.7917065301079622], [0.5603470596005116, 0.7537033611358295], [0.25400471431432753, 0.8150646543307231], [0.3311957483405575, -0.7829906385661866], [0.8445851642438251, -0.08836500410921283], [0.5617849379877878, -0.7574792990349553], [0.07305726944797225, 0.8238422029156591], [0.07358418293812316, -0.8256623870972735], [0.5883410243610833, 0.6196215619927677], [0.2566752832473187, -0.8195973033555625], [0.841341592761485, -0.16714383527204318], [0.7514167069759504, 0.40985942538495446], [0.7553021899259311, -0.4137557280795954], [0.8372601206719111, 0.0704318953584988], [0.9360577467536026, 0.13890639314725461]], "t": [[32, 33, 72], [2, 51, 73], [67, 150, 151], [67, 150, 186], [0, 32, 33], [1, 2, 51], [124, 142, 153], [31, 70, 72], [3, 71, 73], [31, 32, 72], [2, 3, 73], [68, 137, 193], [70, 72, 176], [59, 140, 142], [118, 191, 192], [68, 137, 201], [140, 142, 153], [20, 21, 190], [59, 140, 141], [118, 191, 217], [124, 142, 204], [71, 73, 134], [72, 176, 213], [69, 142, 204], [131, 184, 205], [184, 205, 208], [37, 60, 177], [60, 130, 177], [69, 142, 199], [20, 67, 190], [21, 54, 190], [64, 129, 175], [64, 129, 176], [118, 192, 203], [129, 175, 178], [38, 39, 188], [60, 175, 178], [196, 204, 219], [37, 38, 177], [60, 130, 178], [98, 159, 180], [67, 151, 190], [167, 168, 173], [83, 168, 173], [106, 163, 181], [98, 159, 183], [98, 180, 182], [38, 177, 188], [163, 180, 182], [69, 196, 204], [57, 132, 167], [132, 167, 173], [162, 163, 181], [146, 152, 155], [158, 159, 180], [106, 163, 182], [127, 149, 151], [127, 149, 154], [135, 143, 144], [106, 179, 181], [75, 133, 144], [133, 135, 144], [57, 132, 170], [165, 167, 168], [99, 145, 155], [145, 146, 155], [45, 57, 166], [124, 153, 169], [57, 166, 167], [35, 64, 175], [165, 166, 167], [168, 179, 181], [65, 133, 134], [158, 160, 161], [83, 168, 179], [102, 119, 136], [128, 136, 137], [62, 74, 75], [99, 145, 171], [119, 136, 137], [133, 134, 135], [75, 77, 144], [65, 74, 133], [100, 147, 148], [96, 111, 154], [48, 62, 74], [107, 138, 139], [111, 127, 154], [148, 152, 155], [147, 148, 155], [74, 75, 133], [92, 100, 147], [44, 45, 166], [81, 153, 169], [102, 103, 136], [79, 107, 138], [88, 99, 171], [46, 57, 170], [156, 158, 160], [100, 149, 154], [138, 140, 153], [62, 76, 170], [41, 156, 160], [76, 132, 170], [138, 139, 140], [77, 79, 107], [139, 140, 141], [75, 76, 77], [62, 75, 76], [43, 44, 164], [81, 138, 153], [110, 124, 169], [49, 65, 74], [77, 78, 79], [114, 120, 131], [79, 81, 138], [100, 148, 149], [94, 96, 154], [108, 114, 128], [56, 156, 157], [85, 110, 169], [111, 119, 126], [161, 162, 163], [157, 158, 159], [148, 149, 150], [34, 35, 64], [101, 110, 116], [49, 50, 65], [55, 160, 161], [96, 102, 111], [102, 111, 119], [110, 116, 123], [108, 109, 114], [79, 80, 81], [99, 147, 155], [44, 164, 166], [85, 101, 110], [90, 92, 147], [92, 94, 100], [85, 86, 101], [158, 161, 180], [116, 117, 118], [47, 48, 62], [149, 150, 151], [88, 90, 99], [94, 100, 154], [156, 157, 158], [55, 161, 162], [103, 108, 136], [82, 84, 85], [76, 77, 78], [107, 139, 143], [90, 99, 147], [36, 37, 60], [40, 56, 156], [148, 150, 152], [42, 55, 160], [78, 79, 80], [109, 114, 120], [117, 118, 125], [101, 116, 117], [81, 82, 169], [48, 49, 74], [84, 85, 86], [39, 40, 56], [111, 126, 127], [96, 97, 102], [103, 105, 108], [45, 46, 57], [80, 81, 82], [42, 43, 55], [101, 117, 171], [94, 95, 96], [162, 164, 165], [117, 145, 171], [108, 128, 136], [41, 42, 160], [92, 93, 94], [82, 85, 169], [40, 41, 156], [90, 91, 92], [43, 55, 164], [105, 108, 109], [110, 123, 124], [97, 102, 103], [47, 62, 170], [88, 89, 90], [86, 87, 88], [125, 145, 146], [86, 88, 171], [77, 107, 144], [82, 83, 84], [46, 47, 170], [95, 96, 97], [162, 165, 181], [55, 162, 164], [86, 101, 171], [93, 94, 95], [164, 165, 166], [64, 176, 213], [120, 121, 122], [91, 92, 93], [117, 125, 145], [89, 90, 91], [84, 86, 87], [87, 88, 89], [15, 16, 66], [19, 20, 67], [103, 104, 105], [80, 82, 83], [23, 24, 68], [76, 78, 132], [21, 22, 54], [17, 18, 52], [10, 11, 69], [12, 13, 53], [25, 26, 58], [35, 36, 175], [97, 103, 104], [8, 9, 59], [109, 120, 121], [137, 193, 215], [107, 143, 144], [91, 93, 98], [112, 115, 130], [28, 29, 63], [87, 89, 106], [5, 6, 61], [104, 112, 113], [122, 129, 176], [112, 130, 172], [105, 109, 172], [78, 80, 173], [104, 105, 112], [78, 132, 173], [80, 83, 173], [95, 97, 174], [112, 113, 115], [36, 60, 175], [97, 104, 174], [121, 122, 129], [109, 121, 172], [30, 31, 70], [3, 4, 71], [104, 113, 174], [105, 112, 172], [83, 84, 179], [115, 130, 177], [39, 56, 188], [87, 106, 179], [91, 98, 182], [89, 106, 182], [56, 157, 187], [93, 98, 183], [121, 172, 178], [95, 174, 183], [165, 168, 181], [157, 159, 185], [84, 87, 179], [121, 129, 178], [130, 172, 178], [161, 163, 180], [33, 72, 213], [51, 73, 214], [89, 91, 182], [114, 128, 184], [93, 95, 183], [56, 187, 188], [157, 185, 187], [141, 189, 207], [73, 134, 214], [174, 183, 185], [114, 131, 184], [141, 207, 212], [139, 143, 189], [159, 183, 185], [68, 194, 201], [143, 189, 210], [194, 201, 218], [113, 185, 187], [113, 174, 185], [113, 115, 187], [53, 192, 203], [115, 187, 188], [115, 177, 188], [59, 142, 199], [197, 206, 209], [118, 125, 217], [119, 137, 201], [139, 141, 189], [128, 184, 215], [61, 189, 210], [63, 131, 206], [116, 123, 203], [58, 184, 215], [120, 131, 206], [65, 134, 214], [63, 197, 206], [33, 34, 213], [50, 51, 214], [125, 146, 211], [70, 176, 209], [14, 15, 191], [16, 17, 195], [22, 23, 194], [13, 14, 192], [24, 25, 193], [11, 12, 196], [15, 66, 191], [16, 66, 195], [23, 68, 194], [22, 54, 194], [24, 68, 193], [17, 52, 195], [11, 69, 196], [13, 53, 192], [12, 53, 196], [25, 58, 193], [29, 30, 197], [4, 5, 198], [150, 152, 220], [29, 63, 197], [5, 61, 198], [123, 204, 219], [30, 70, 197], [4, 71, 198], [127, 151, 202], [10, 69, 199], [9, 59, 199], [126, 127, 202], [150, 186, 220], [63, 131, 205], [59, 141, 212], [28, 63, 205], [146, 152, 200], [26, 58, 208], [6, 61, 207], [9, 10, 199], [19, 67, 221], [18, 52, 221], [8, 59, 212], [119, 126, 201], [14, 191, 192], [58, 184, 208], [125, 211, 217], [151, 190, 202], [67, 186, 221], [52, 186, 221], [61, 189, 207], [116, 118, 203], [123, 124, 204], [61, 198, 210], [71, 134, 216], [27, 28, 205], [120, 122, 206], [134, 135, 216], [52, 195, 200], [122, 176, 209], [135, 143, 210], [58, 193, 215], [126, 201, 218], [198, 210, 216], [123, 203, 219], [6, 7, 207], [26, 27, 208], [54, 190, 202], [122, 206, 209], [128, 137, 215], [7, 8, 212], [34, 64, 213], [50, 65, 214], [18, 19, 221], [66, 195, 211], [27, 205, 208], [7, 207, 212], [70, 197, 209], [52, 186, 220], [195, 200, 211], [126, 202, 218], [66, 211, 217], [135, 210, 216], [66, 191, 217], [71, 198, 216], [146, 200, 211], [54, 194, 218], [53, 196, 219], [54, 202, 218], [52, 200, 220], [152, 200, 220], [53, 203, 219]], "boundaries": {"contact": [0, 2, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31, 34, 37, 40, 43, 46, 49, 52, 55, 58, 61, 64, 67, 70, 73, 76, 79, 82, 85, 88, 91], "dirichlet": [1, 3, 96, 99, 102, 105, 108, 111, 114, 117, 120, 123, 126, 129, 132, 135, 138, 141, 144, 147]}, "subdomains": {"domain": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389]}} \ No newline at end of file diff --git a/docs/examples/ex04_mesh.off b/docs/examples/ex04_mesh.off deleted file mode 100644 index 8fa111c23..000000000 --- a/docs/examples/ex04_mesh.off +++ /dev/null @@ -1,617 +0,0 @@ -OFF -# Created by meshio - -222 390 0 - -0.0 1.0 0.0 -0.0 -1.0 0.0 -0.09801714055591637 -0.9951847266499028 0.0 -0.19509032241883761 -0.9807852803231266 0.0 -0.2902846778931259 -0.9569403355384724 0.0 -0.38268343326549487 -0.9238795321383267 0.0 -0.47139673797238385 -0.8819212637355989 0.0 -0.5555702343496238 -0.8314696114138532 0.0 -0.6343932856511478 -0.7730104521419754 0.0 -0.7071067827656586 -0.7071067796074364 0.0 -0.773010454973133 -0.6343932822013723 0.0 -0.8314696138527624 -0.5555702306995381 0.0 -0.8819212658270285 -0.47139673405959426 0.0 -0.9238795338254099 -0.38268342919251586 0.0 -0.9569403368184075 -0.2902846736737454 0.0 -0.9807852811883714 -0.1950903180689581 0.0 -0.995184727094887 -0.09801713603791533 0.0 -1.0 4.5837337757057685e-09 0.0 -0.9951847262476915 0.09801714463963587 0.0 -0.980785279615548 0.19509032597607545 0.0 -0.9569403346479524 0.290284680828777 0.0 -0.9238795311905263 0.38268343555368756 0.0 -0.8819212628015901 0.4713967397197912 0.0 -0.8314696106349614 0.5555702355153178 0.0 -0.7730104516368775 0.6343932862666112 0.0 -0.7071067794805702 0.7071067828925248 0.0 -0.634393282512301 0.7730104547179604 0.0 -0.5555702315457777 0.8314696132873234 0.0 -0.4713967355549627 0.8819212650277374 0.0 -0.38268343125990784 0.9238795329690681 0.0 -0.2902846763567968 0.9569403360045128 0.0 -0.195090321428782 0.9807852805200609 0.0 -0.09801714007694617 0.9951847266970772 0.0 -0.0 0.9000000000002774 0.0 -0.0 0.8000000000005548 0.0 -0.0 0.7000000000008322 0.0 -0.0 0.6000000000011096 0.0 -0.0 0.500000000001387 0.0 -0.0 0.40000000000166436 0.0 -0.0 0.3000000000019418 0.0 -0.0 0.20000000000221918 0.0 -0.0 0.10000000000249643 0.0 -0.0 2.752242878045763e-12 0.0 -0.0 -0.09999999999750342 0.0 -0.0 -0.19999999999778084 0.0 -0.0 -0.29999999999805826 0.0 -0.0 -0.39999999999833546 0.0 -0.0 -0.4999999999986129 0.0 -0.0 -0.5999999999988903 0.0 -0.0 -0.699999999999168 0.0 -0.0 -0.7999999999994454 0.0 -0.0 -0.8999999999997226 0.0 -0.9159764153258458 0.05383987386802238 0.0 -0.8314138091584239 -0.39785805095354165 0.0 -0.8264043258240841 0.39254386345989056 0.0 -0.08657154386342972 -0.04993181939059333 0.0 -0.08825881533998217 0.24854065189332641 0.0 -0.0866467339337712 -0.350125666409955 0.0 -0.6151960220089083 0.6914850111403807 0.0 -0.611656851878296 -0.6767753615474275 0.0 -0.08792437843765298 0.546553776392746 0.0 -0.40214879730821057 -0.8329018889819535 0.0 -0.0864975927752666 -0.5501598487639185 0.0 -0.4031239663030456 0.8299929282327413 0.0 -0.08552443345288481 0.7412874524475859 0.0 -0.08387010326216801 -0.7452678551627052 0.0 -0.9064888004049059 -0.13475830751070036 0.0 -0.883006687779775 0.21591742664103838 0.0 -0.7318295262956729 0.5462376817273685 0.0 -0.7310857439949313 -0.5534884214698972 0.0 -0.22203717119110486 0.8864220786001916 0.0 -0.222037172228374 -0.88642207782408 0.0 -0.12735200151713993 0.8953724394140379 0.0 -0.1264872216201089 -0.8973542598638573 0.0 -0.08603181563097448 -0.6492168486242395 0.0 -0.1728569295463326 -0.5999642076514183 0.0 -0.1729345952627196 -0.500480145586956 0.0 -0.2591331121034306 -0.5504140700880753 0.0 -0.2592714087128541 -0.4510637781789997 0.0 -0.3451538888162647 -0.5008525665656047 0.0 -0.3453660173120208 -0.4016689274236115 0.0 -0.43112082579940325 -0.45140415758047114 0.0 -0.4313132869549403 -0.35228060068609146 0.0 -0.3456417826557642 -0.30202707863747263 0.0 -0.4316914876509484 -0.25313887750341396 0.0 -0.5172439124243171 -0.30297215045699327 0.0 -0.5175060234485399 -0.20398089749625326 0.0 -0.4320078005036593 -0.1541618871667654 0.0 -0.5178249360980043 -0.10497328859987322 0.0 -0.4323517746522692 -0.05513829699300258 0.0 -0.5181810627769347 -0.005926661553542319 0.0 -0.4326873030792038 0.0439068356495514 0.0 -0.5185129534634814 0.0931212502063908 0.0 -0.43303358988383006 0.14295033231019513 0.0 -0.5188407112939315 0.19216073999298827 0.0 -0.4333907822760996 0.2419806037163972 0.0 -0.5191570773658916 0.29117889053487644 0.0 -0.4336167327285912 0.34098097416379747 0.0 -0.3466297470498466 0.09502752487561306 0.0 -0.6036466835467159 -0.05588497882022823 0.0 -0.6043483291878866 0.14226361071690696 0.0 -0.6030083111892219 -0.253739769204638 0.0 -0.5194358558074463 0.3901995369490897 0.0 -0.4339131655358476 0.44001601280296204 0.0 -0.34799605056576544 0.3907547089718833 0.0 -0.3483736150587413 0.4897859002884928 0.0 -0.3461095808402675 -0.10314722872110887 0.0 -0.3451380056642011 -0.5999795967606039 0.0 -0.43423779016255415 0.5389955255642274 0.0 -0.3487855783748346 0.5887556881868931 0.0 -0.6027356100339158 -0.35253659554114325 0.0 -0.6048915795006872 0.34040851731280114 0.0 -0.26231771840506496 0.4405891472347052 0.0 -0.26200491468128767 0.34136888541754196 0.0 -0.4341509042932531 0.6378572118248705 0.0 -0.17613909944252154 0.39128549747516345 0.0 -0.6867867143117476 -0.3005776176954207 0.0 -0.6883938442262401 -0.20422841407407338 0.0 -0.7887504482032948 -0.24363163143334435 0.0 -0.6039194807805792 0.4371233138451532 0.0 -0.34717469685206775 0.6842928036187731 0.0 -0.2630163234652642 0.6380173006925967 0.0 -0.26231461886291657 0.7347372834745751 0.0 -0.6821568960369653 -0.38932236004593684 0.0 -0.6024613510918597 -0.4500576975192292 0.0 -0.7718417301410816 -0.14476793206076038 0.0 -0.6840585160728114 0.38136305112776203 0.0 -0.6872609976885783 0.28791528538360167 0.0 -0.5168363646187115 0.5833344105436389 0.0 -0.17550407052083822 0.6900813696897046 0.0 -0.17619557798382957 0.4914856785795687 0.0 -0.42016802198142106 0.7346497809142203 0.0 -0.17302511054291092 -0.40062265779056283 0.0 -0.17237477181374233 -0.6988080454878537 0.0 -0.1714952098378757 -0.8012321435395398 0.0 -0.25599371682296734 -0.7410015264285403 0.0 -0.5190057914509976 0.4879530306260381 0.0 -0.618290339602385 0.536875949487155 0.0 -0.4308956703486012 -0.5504811711457109 0.0 -0.43130650380920504 -0.6485484160152994 0.0 -0.5175687390320705 -0.6021217093581201 0.0 -0.5092206012248724 -0.7035085970951118 0.0 -0.620737462322611 -0.555765633253074 0.0 -0.3429556442908701 -0.6929402428430802 0.0 -0.25848406726552225 -0.6479093927842566 0.0 -0.6885382570987393 -0.10392056646537914 0.0 -0.769286205655087 -0.05057534362618522 0.0 -0.6040111157907584 0.043242214936212164 0.0 -0.6897389251565653 0.09246051544826223 0.0 -0.6895451688254398 0.1909341621919352 0.0 -0.7853297010802 0.13990094156251878 0.0 -0.7756810503792961 0.2413610566339871 0.0 -0.7701883492540708 0.03853734593193141 0.0 -0.5190243392045759 -0.5017121515376488 0.0 -0.6040073106437357 0.2408102010221849 0.0 -0.687568256083656 -0.006023468765897795 0.0 -0.08680748951240319 0.14978745087550502 0.0 -0.17485631474456298 0.19832056843973647 0.0 -0.1730517611561072 0.10005310843209428 0.0 -0.26087327445847697 0.14679457604766757 0.0 -0.08657886954653793 0.04998478457738213 0.0 -0.17291704371860217 -0.00010105186068302731 0.0 -0.1729970659840663 -0.10002826078482492 0.0 -0.2593646476016979 -0.04985929515148815 0.0 -0.08660511306657999 -0.1498281482654842 0.0 -0.17305207548825177 -0.20001247954407836 0.0 -0.08658138073658056 -0.24999061371589165 0.0 -0.17306526794695942 -0.3003972076357125 0.0 -0.2594694617813657 -0.25019468739155504 0.0 -0.517316554251502 -0.4018272255535962 0.0 -0.08648158353870185 -0.4502095033908465 0.0 -0.6031530092679103 -0.1544546524434075 0.0 -0.26271089280937426 0.539695638559989 0.0 -0.25930650815864587 -0.35099572284298564 0.0 -0.3481417000647074 0.29225483630730303 0.0 -0.08775490432998115 0.6445771191516094 0.0 -0.17312433591087684 0.7949143361266013 0.0 -0.08917340860778944 0.4445319633762133 0.0 -0.17610593477516273 0.5911079955205527 0.0 -0.34591484982492143 -0.20259666030335644 0.0 -0.2596845571558199 0.04809501117136811 0.0 -0.259197333421417 -0.15100131147531112 0.0 -0.34645554142539703 -0.004036758239607918 0.0 -0.3476513202248709 0.19411470667553155 0.0 -0.5165083519158323 0.6874008844465628 0.0 -0.26285614408485514 0.24520073416291943 0.0 -0.8702436047190465 0.12724803991148598 0.0 -0.1765715058218143 0.29520431385172685 0.0 -0.08940011612821333 0.34630648225343424 0.0 -0.4166298503101692 -0.7429585203863014 0.0 -0.8549693515455798 0.30930121195949895 0.0 -0.8950418208143169 -0.2241963012461981 0.0 -0.860025803539879 -0.3124749658989302 0.0 -0.6724940254177364 0.6208040077969013 0.0 -0.7854994399143052 0.4754049055547044 0.0 -0.9215823721068173 -0.045274434635522595 0.0 -0.785420238160406 -0.4789576989502909 0.0 -0.31501693368973693 0.863121801168385 0.0 -0.31408164224217594 -0.8662391964201028 0.0 -0.688719459186926 -0.6255058956158415 0.0 -0.844562532980393 -0.007366724481464558 0.0 -0.6969550069593997 0.46476187067423613 0.0 -0.7641662461091518 0.3352234361081791 0.0 -0.766062195642033 -0.3394573773342276 0.0 -0.6967061893284834 -0.4732097765718286 0.0 -0.48759115339420067 0.7864769275402721 0.0 -0.33363382533391917 0.7769765419565695 0.0 -0.4861251931921746 -0.7917065301079622 0.0 -0.5603470596005116 0.7537033611358295 0.0 -0.25400471431432753 0.8150646543307231 0.0 -0.3311957483405575 -0.7829906385661866 0.0 -0.8445851642438251 -0.08836500410921283 0.0 -0.5617849379877878 -0.7574792990349553 0.0 -0.07305726944797225 0.8238422029156591 0.0 -0.07358418293812316 -0.8256623870972735 0.0 -0.5883410243610833 0.6196215619927677 0.0 -0.2566752832473187 -0.8195973033555625 0.0 -0.841341592761485 -0.16714383527204318 0.0 -0.7514167069759504 0.40985942538495446 0.0 -0.7553021899259311 -0.4137557280795954 0.0 -0.8372601206719111 0.0704318953584988 0.0 -0.9360577467536026 0.13890639314725461 0.0 -3 32 33 72 -3 2 51 73 -3 67 150 151 -3 67 150 186 -3 0 32 33 -3 1 2 51 -3 124 142 153 -3 31 70 72 -3 3 71 73 -3 31 32 72 -3 2 3 73 -3 68 137 193 -3 70 72 176 -3 59 140 142 -3 118 191 192 -3 68 137 201 -3 140 142 153 -3 20 21 190 -3 59 140 141 -3 118 191 217 -3 124 142 204 -3 71 73 134 -3 72 176 213 -3 69 142 204 -3 131 184 205 -3 184 205 208 -3 37 60 177 -3 60 130 177 -3 69 142 199 -3 20 67 190 -3 21 54 190 -3 64 129 175 -3 64 129 176 -3 118 192 203 -3 129 175 178 -3 38 39 188 -3 60 175 178 -3 196 204 219 -3 37 38 177 -3 60 130 178 -3 98 159 180 -3 67 151 190 -3 167 168 173 -3 83 168 173 -3 106 163 181 -3 98 159 183 -3 98 180 182 -3 38 177 188 -3 163 180 182 -3 69 196 204 -3 57 132 167 -3 132 167 173 -3 162 163 181 -3 146 152 155 -3 158 159 180 -3 106 163 182 -3 127 149 151 -3 127 149 154 -3 135 143 144 -3 106 179 181 -3 75 133 144 -3 133 135 144 -3 57 132 170 -3 165 167 168 -3 99 145 155 -3 145 146 155 -3 45 57 166 -3 124 153 169 -3 57 166 167 -3 35 64 175 -3 165 166 167 -3 168 179 181 -3 65 133 134 -3 158 160 161 -3 83 168 179 -3 102 119 136 -3 128 136 137 -3 62 74 75 -3 99 145 171 -3 119 136 137 -3 133 134 135 -3 75 77 144 -3 65 74 133 -3 100 147 148 -3 96 111 154 -3 48 62 74 -3 107 138 139 -3 111 127 154 -3 148 152 155 -3 147 148 155 -3 74 75 133 -3 92 100 147 -3 44 45 166 -3 81 153 169 -3 102 103 136 -3 79 107 138 -3 88 99 171 -3 46 57 170 -3 156 158 160 -3 100 149 154 -3 138 140 153 -3 62 76 170 -3 41 156 160 -3 76 132 170 -3 138 139 140 -3 77 79 107 -3 139 140 141 -3 75 76 77 -3 62 75 76 -3 43 44 164 -3 81 138 153 -3 110 124 169 -3 49 65 74 -3 77 78 79 -3 114 120 131 -3 79 81 138 -3 100 148 149 -3 94 96 154 -3 108 114 128 -3 56 156 157 -3 85 110 169 -3 111 119 126 -3 161 162 163 -3 157 158 159 -3 148 149 150 -3 34 35 64 -3 101 110 116 -3 49 50 65 -3 55 160 161 -3 96 102 111 -3 102 111 119 -3 110 116 123 -3 108 109 114 -3 79 80 81 -3 99 147 155 -3 44 164 166 -3 85 101 110 -3 90 92 147 -3 92 94 100 -3 85 86 101 -3 158 161 180 -3 116 117 118 -3 47 48 62 -3 149 150 151 -3 88 90 99 -3 94 100 154 -3 156 157 158 -3 55 161 162 -3 103 108 136 -3 82 84 85 -3 76 77 78 -3 107 139 143 -3 90 99 147 -3 36 37 60 -3 40 56 156 -3 148 150 152 -3 42 55 160 -3 78 79 80 -3 109 114 120 -3 117 118 125 -3 101 116 117 -3 81 82 169 -3 48 49 74 -3 84 85 86 -3 39 40 56 -3 111 126 127 -3 96 97 102 -3 103 105 108 -3 45 46 57 -3 80 81 82 -3 42 43 55 -3 101 117 171 -3 94 95 96 -3 162 164 165 -3 117 145 171 -3 108 128 136 -3 41 42 160 -3 92 93 94 -3 82 85 169 -3 40 41 156 -3 90 91 92 -3 43 55 164 -3 105 108 109 -3 110 123 124 -3 97 102 103 -3 47 62 170 -3 88 89 90 -3 86 87 88 -3 125 145 146 -3 86 88 171 -3 77 107 144 -3 82 83 84 -3 46 47 170 -3 95 96 97 -3 162 165 181 -3 55 162 164 -3 86 101 171 -3 93 94 95 -3 164 165 166 -3 64 176 213 -3 120 121 122 -3 91 92 93 -3 117 125 145 -3 89 90 91 -3 84 86 87 -3 87 88 89 -3 15 16 66 -3 19 20 67 -3 103 104 105 -3 80 82 83 -3 23 24 68 -3 76 78 132 -3 21 22 54 -3 17 18 52 -3 10 11 69 -3 12 13 53 -3 25 26 58 -3 35 36 175 -3 97 103 104 -3 8 9 59 -3 109 120 121 -3 137 193 215 -3 107 143 144 -3 91 93 98 -3 112 115 130 -3 28 29 63 -3 87 89 106 -3 5 6 61 -3 104 112 113 -3 122 129 176 -3 112 130 172 -3 105 109 172 -3 78 80 173 -3 104 105 112 -3 78 132 173 -3 80 83 173 -3 95 97 174 -3 112 113 115 -3 36 60 175 -3 97 104 174 -3 121 122 129 -3 109 121 172 -3 30 31 70 -3 3 4 71 -3 104 113 174 -3 105 112 172 -3 83 84 179 -3 115 130 177 -3 39 56 188 -3 87 106 179 -3 91 98 182 -3 89 106 182 -3 56 157 187 -3 93 98 183 -3 121 172 178 -3 95 174 183 -3 165 168 181 -3 157 159 185 -3 84 87 179 -3 121 129 178 -3 130 172 178 -3 161 163 180 -3 33 72 213 -3 51 73 214 -3 89 91 182 -3 114 128 184 -3 93 95 183 -3 56 187 188 -3 157 185 187 -3 141 189 207 -3 73 134 214 -3 174 183 185 -3 114 131 184 -3 141 207 212 -3 139 143 189 -3 159 183 185 -3 68 194 201 -3 143 189 210 -3 194 201 218 -3 113 185 187 -3 113 174 185 -3 113 115 187 -3 53 192 203 -3 115 187 188 -3 115 177 188 -3 59 142 199 -3 197 206 209 -3 118 125 217 -3 119 137 201 -3 139 141 189 -3 128 184 215 -3 61 189 210 -3 63 131 206 -3 116 123 203 -3 58 184 215 -3 120 131 206 -3 65 134 214 -3 63 197 206 -3 33 34 213 -3 50 51 214 -3 125 146 211 -3 70 176 209 -3 14 15 191 -3 16 17 195 -3 22 23 194 -3 13 14 192 -3 24 25 193 -3 11 12 196 -3 15 66 191 -3 16 66 195 -3 23 68 194 -3 22 54 194 -3 24 68 193 -3 17 52 195 -3 11 69 196 -3 13 53 192 -3 12 53 196 -3 25 58 193 -3 29 30 197 -3 4 5 198 -3 150 152 220 -3 29 63 197 -3 5 61 198 -3 123 204 219 -3 30 70 197 -3 4 71 198 -3 127 151 202 -3 10 69 199 -3 9 59 199 -3 126 127 202 -3 150 186 220 -3 63 131 205 -3 59 141 212 -3 28 63 205 -3 146 152 200 -3 26 58 208 -3 6 61 207 -3 9 10 199 -3 19 67 221 -3 18 52 221 -3 8 59 212 -3 119 126 201 -3 14 191 192 -3 58 184 208 -3 125 211 217 -3 151 190 202 -3 67 186 221 -3 52 186 221 -3 61 189 207 -3 116 118 203 -3 123 124 204 -3 61 198 210 -3 71 134 216 -3 27 28 205 -3 120 122 206 -3 134 135 216 -3 52 195 200 -3 122 176 209 -3 135 143 210 -3 58 193 215 -3 126 201 218 -3 198 210 216 -3 123 203 219 -3 6 7 207 -3 26 27 208 -3 54 190 202 -3 122 206 209 -3 128 137 215 -3 7 8 212 -3 34 64 213 -3 50 65 214 -3 18 19 221 -3 66 195 211 -3 27 205 208 -3 7 207 212 -3 70 197 209 -3 52 186 220 -3 195 200 211 -3 126 202 218 -3 66 211 217 -3 135 210 216 -3 66 191 217 -3 71 198 216 -3 146 200 211 -3 54 194 218 -3 53 196 219 -3 54 202 218 -3 52 200 220 -3 152 200 220 -3 53 203 219 From fb914104dae70557c9d7a66c372427a13be2ce35 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 16:53:48 +0200 Subject: [PATCH 22/32] Rename skfem.importers to skfem.io --- docs/examples/ex04.py | 20 +++--- skfem/importers/__init__.py | 7 +- skfem/importers/meshio.py | 121 +------------------------------- skfem/io/__init__.py | 1 + skfem/{importers => io}/json.py | 0 skfem/io/meshio.py | 121 ++++++++++++++++++++++++++++++++ 6 files changed, 141 insertions(+), 129 deletions(-) create mode 100644 skfem/io/__init__.py rename skfem/{importers => io}/json.py (100%) create mode 100644 skfem/io/meshio.py diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index 7499c76f8..0f1a59fb6 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -3,8 +3,8 @@ import numpy as np from pygmsh import generate_mesh from pygmsh.built_in import Geometry -from skfem.importers import from_meshio -from skfem.importers.json import from_file, to_file +from skfem.io import from_meshio +from skfem.io.json import from_file, to_file import meshio @@ -191,15 +191,15 @@ def mass(u, du, v, dv, w): s[2, 2] = nu1 * (s[0, 0] + s[1, 1]) S[2, 2] = nu2 * (S[0, 0] + S[1, 1]) -vonmises1 = np.sqrt(.5*((s[0, 0] - s[1, 1])**2 +\ - (s[1, 1] - s[2, 2])**2 +\ - (s[2, 2] - s[0, 0])**2 +\ - 6.0*s[0, 1]**2)) +vonmises1 = np.sqrt(.5 * ((s[0, 0] - s[1, 1]) ** 2 + + (s[1, 1] - s[2, 2]) ** 2 + + (s[2, 2] - s[0, 0]) ** 2 + + 6. * s[0, 1]**2)) -vonmises2 = np.sqrt(.5*((S[0, 0] - S[1, 1])**2 +\ - (S[1, 1] - S[2, 2])**2 +\ - (S[2, 2] - S[0, 0])**2 +\ - 6.0*S[0, 1]**2)) +vonmises2 = np.sqrt(.5 * ((S[0, 0] - S[1, 1]) ** 2 + + (S[1, 1] - S[2, 2]) ** 2 + + (S[2, 2] - S[0, 0]) ** 2 + + 6. * S[0, 1]**2)) if __name__ == "__main__": diff --git a/skfem/importers/__init__.py b/skfem/importers/__init__.py index d65017e61..8fb1158f2 100644 --- a/skfem/importers/__init__.py +++ b/skfem/importers/__init__.py @@ -1 +1,6 @@ -from .meshio import from_meshio +import warnings + +warnings.warn("DeprecationWarning: skfem.importers was renamed to " + "skfem.io and is removed in the next major release.") + +from ..io import * diff --git a/skfem/importers/meshio.py b/skfem/importers/meshio.py index 977a13e5d..60ecc391c 100644 --- a/skfem/importers/meshio.py +++ b/skfem/importers/meshio.py @@ -1,121 +1,6 @@ -"""Import any formats supported by meshio.""" - import warnings -from collections import OrderedDict - -import meshio -import numpy as np - -import skfem - - -MESH_TYPE_MAPPING = OrderedDict([ - ('tetra', skfem.MeshTet), - ('hexahedron', skfem.MeshHex), - ('triangle', skfem.MeshTri), - ('quad', skfem.MeshQuad), - ('line', skfem.MeshLine), -]) - - -def from_meshio(m, force_mesh_type=None): - """Convert meshio mesh into :class:`skfem.mesh.Mesh`. - - Parameters - ---------- - m - The mesh object from meshio. - force_mesh_type - An optional string forcing the mesh type if automatic detection - fails. See skfem.importers.meshio.MESH_TYPE_MAPPING for possible values. - - """ - - if force_mesh_type is None: - for k, v in MESH_TYPE_MAPPING.items(): - # find first if match - if k in m.cells: - meshio_type, mesh_type = (k, MESH_TYPE_MAPPING[k]) - break - else: - meshio_type, mesh_type = (force_mesh_type, - MESH_TYPE_MAPPING[force_mesh_type]) - - def strip_extra_coordinates(p): - if meshio_type == "line": - return p[:, :1] - if meshio_type in ("quad", "triangle"): - return p[:, :2] - return p - - # create p and t - p = np.ascontiguousarray(strip_extra_coordinates(m.points).T) - t = np.ascontiguousarray(m.cells[meshio_type].T) - - mtmp = mesh_type(p, - (t[[0, 4, 3, 1, 7, 5, 2, 6]] - if meshio_type == 'hexahedron' - # vtk requires a different ordering - else t)) - - try: - # element to boundary element type mapping - bnd_type = { - 'line': 'vertex', - 'triangle': 'line', - 'quad': 'line', - 'tetra': 'triangle', - 'hexahedron': 'quad', - }[meshio_type] - - def find_tagname(tag): - for key in m.field_data: - if m.field_data[key][0] == tag: - return key - return None - - # find subdomains - if meshio_type in m.cell_data and 'gmsh:physical' in m.cell_data[meshio_type]: - elements_tag = m.cell_data[meshio_type]['gmsh:physical'] - - subdomains = {} - tags = np.unique(elements_tag) - - for tag in tags: - t_set = np.nonzero(tag == elements_tag)[0] - subdomains[find_tagname(tag)] = t_set - - # find tagged boundaries - if bnd_type in m.cell_data and 'gmsh:physical' in m.cell_data[bnd_type]: - facets = m.cells[bnd_type] - facets_tag = m.cell_data[bnd_type]['gmsh:physical'] - bndfacets = mtmp.boundary_facets() - - # put meshio facets to dict - dic = {tuple(np.sort(facets[i])): facets_tag[i] - for i in range(facets.shape[0])} - - # get index of corresponding Mesh.facets for each meshio - # facet found in the dict - index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i] - for i in bndfacets - if tuple(np.sort(mtmp.facets[:, i])) in dic]) - - # read meshio tag numbers and names - tags = index[:, 0] - boundaries = {} - for tag in np.unique(tags): - tagindex = np.nonzero(tags == tag)[0] - boundaries[find_tagname(tag)] = index[tagindex, 1] - - mtmp.boundaries = boundaries - mtmp.subdomains = subdomains - - except Exception as e: - warnings.warn("Unable to load tagged boundaries/subdomains.") - - return mtmp +warnings.warn("DeprecationWarning: skfem.importers was renamed to " + "skfem.io and is removed in the next major release.") -def from_file(filename): - return from_meshio(meshio.read(filename)) +from ..io.meshio import * diff --git a/skfem/io/__init__.py b/skfem/io/__init__.py new file mode 100644 index 000000000..d65017e61 --- /dev/null +++ b/skfem/io/__init__.py @@ -0,0 +1 @@ +from .meshio import from_meshio diff --git a/skfem/importers/json.py b/skfem/io/json.py similarity index 100% rename from skfem/importers/json.py rename to skfem/io/json.py diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py new file mode 100644 index 000000000..977a13e5d --- /dev/null +++ b/skfem/io/meshio.py @@ -0,0 +1,121 @@ +"""Import any formats supported by meshio.""" + +import warnings +from collections import OrderedDict + +import meshio +import numpy as np + +import skfem + + +MESH_TYPE_MAPPING = OrderedDict([ + ('tetra', skfem.MeshTet), + ('hexahedron', skfem.MeshHex), + ('triangle', skfem.MeshTri), + ('quad', skfem.MeshQuad), + ('line', skfem.MeshLine), +]) + + +def from_meshio(m, force_mesh_type=None): + """Convert meshio mesh into :class:`skfem.mesh.Mesh`. + + Parameters + ---------- + m + The mesh object from meshio. + force_mesh_type + An optional string forcing the mesh type if automatic detection + fails. See skfem.importers.meshio.MESH_TYPE_MAPPING for possible values. + + """ + + if force_mesh_type is None: + for k, v in MESH_TYPE_MAPPING.items(): + # find first if match + if k in m.cells: + meshio_type, mesh_type = (k, MESH_TYPE_MAPPING[k]) + break + else: + meshio_type, mesh_type = (force_mesh_type, + MESH_TYPE_MAPPING[force_mesh_type]) + + def strip_extra_coordinates(p): + if meshio_type == "line": + return p[:, :1] + if meshio_type in ("quad", "triangle"): + return p[:, :2] + return p + + # create p and t + p = np.ascontiguousarray(strip_extra_coordinates(m.points).T) + t = np.ascontiguousarray(m.cells[meshio_type].T) + + mtmp = mesh_type(p, + (t[[0, 4, 3, 1, 7, 5, 2, 6]] + if meshio_type == 'hexahedron' + # vtk requires a different ordering + else t)) + + try: + # element to boundary element type mapping + bnd_type = { + 'line': 'vertex', + 'triangle': 'line', + 'quad': 'line', + 'tetra': 'triangle', + 'hexahedron': 'quad', + }[meshio_type] + + def find_tagname(tag): + for key in m.field_data: + if m.field_data[key][0] == tag: + return key + return None + + # find subdomains + if meshio_type in m.cell_data and 'gmsh:physical' in m.cell_data[meshio_type]: + elements_tag = m.cell_data[meshio_type]['gmsh:physical'] + + subdomains = {} + tags = np.unique(elements_tag) + + for tag in tags: + t_set = np.nonzero(tag == elements_tag)[0] + subdomains[find_tagname(tag)] = t_set + + # find tagged boundaries + if bnd_type in m.cell_data and 'gmsh:physical' in m.cell_data[bnd_type]: + facets = m.cells[bnd_type] + facets_tag = m.cell_data[bnd_type]['gmsh:physical'] + bndfacets = mtmp.boundary_facets() + + # put meshio facets to dict + dic = {tuple(np.sort(facets[i])): facets_tag[i] + for i in range(facets.shape[0])} + + # get index of corresponding Mesh.facets for each meshio + # facet found in the dict + index = np.array([[dic[tuple(np.sort(mtmp.facets[:, i]))], i] + for i in bndfacets + if tuple(np.sort(mtmp.facets[:, i])) in dic]) + + # read meshio tag numbers and names + tags = index[:, 0] + boundaries = {} + for tag in np.unique(tags): + tagindex = np.nonzero(tags == tag)[0] + boundaries[find_tagname(tag)] = index[tagindex, 1] + + mtmp.boundaries = boundaries + mtmp.subdomains = subdomains + + except Exception as e: + warnings.warn("Unable to load tagged boundaries/subdomains.") + + return mtmp + + +def from_file(filename): + return from_meshio(meshio.read(filename)) From 3a1afc45538790e411d69a25a5641dc5e600fc5e Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 16:55:24 +0200 Subject: [PATCH 23/32] Add changelog entry --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e66fec27c..9891c0530 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Additional keyword arguments to skfem.utils.solve get passed on to solvers. +### Changed +- Renamed skfem.importers to skfem.io. + ## [0.4.0] - 2020-01-03 ### Changed From 5f31aa9bd1ce18593e6a05a5583a7d91010e417e Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 17:03:14 +0200 Subject: [PATCH 24/32] Make deprecation message better --- docs/examples/ex04.py | 6 +----- skfem/importers/__init__.py | 4 ++-- skfem/importers/meshio.py | 4 ++-- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index 0f1a59fb6..05fb8931e 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -40,8 +40,6 @@ ib = InteriorBasis(m, e, intorder=4) Ib = InteriorBasis(M, E, intorder=4) -fb = FacetBasis(m, e) -Fb = FacetBasis(M, E) mappings = MortarPair.init_2D(m, M, m.boundaries['contact'], @@ -50,10 +48,9 @@ mb = [ MortarBasis(m, e, mapping = mappings[0], intorder=4), - MortarBasis(M, E, mapping = mappings[1], intorder=4) + MortarBasis(M, E, mapping = mappings[1], intorder=4), ] - # define bilinear forms E1 = 1000.0 E2 = 1000.0 @@ -77,7 +74,6 @@ alpha = 1000 limit = 0.3 - # assemble the stiffness matrices K1 = asm(weakform1, ib) K2 = asm(weakform2, Ib) diff --git a/skfem/importers/__init__.py b/skfem/importers/__init__.py index 8fb1158f2..a08941c99 100644 --- a/skfem/importers/__init__.py +++ b/skfem/importers/__init__.py @@ -1,6 +1,6 @@ import warnings -warnings.warn("DeprecationWarning: skfem.importers was renamed to " - "skfem.io and is removed in the next major release.") +warnings.warn("DeprecationWarning: skfem.importers is removed " + "in the next release. Use skfem.io instead.") from ..io import * diff --git a/skfem/importers/meshio.py b/skfem/importers/meshio.py index 60ecc391c..101abc08a 100644 --- a/skfem/importers/meshio.py +++ b/skfem/importers/meshio.py @@ -1,6 +1,6 @@ import warnings -warnings.warn("DeprecationWarning: skfem.importers was renamed to " - "skfem.io and is removed in the next major release.") +warnings.warn("DeprecationWarning: skfem.importers is removed " + "in the next release. Use skfem.io instead.") from ..io.meshio import * From e0ac052ff75495f5aaeec7ac8cbe39f5302c8c6a Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 17 Jan 2020 17:05:43 +0200 Subject: [PATCH 25/32] Add changelog entry about mortaring --- CHANGELOG.md | 1 + skfem/__init__.py | 11 ----------- skfem/assembly/basis/mortar_basis.py | 2 -- 3 files changed, 1 insertion(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9891c0530..67cdb2eb3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Additional keyword arguments to skfem.utils.solve get passed on to solvers. +- Add MortarBasis and MortarPair for mortaring in 2D. ### Changed - Renamed skfem.importers to skfem.io. diff --git a/skfem/__init__.py b/skfem/__init__.py index 417919490..7cb6ce084 100644 --- a/skfem/__init__.py +++ b/skfem/__init__.py @@ -1,16 +1,5 @@ """Support for wildcard import.""" - -def mark_unstable(fun): - def call(*args, **kwargs): - import warnings - warnings.warn("You are using an unstable feature '{}' that might " - "change in the near future. It is not part of the " - "publicly documented API.".format(fun.__module__)) - return fun(*args, **kwargs) - return call - - from skfem.mesh import * from skfem.assembly import * from skfem.mapping import * diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index 87bb734c0..ad9d1d535 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -1,7 +1,6 @@ import numpy as np from numpy import ndarray -from skfem import mark_unstable from ...quadrature import get_quadrature from .basis import Basis from .interior_basis import InteriorBasis @@ -10,7 +9,6 @@ class MortarBasis(Basis): """Global basis functions at integration points on the mortar boundary.""" - @mark_unstable def __init__(self, mesh, elem, From f66c484f9cf0a90054c75f93cc7105558383a37f Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sat, 18 Jan 2020 23:53:49 +0200 Subject: [PATCH 26/32] Use DeprecationWarning category --- skfem/importers/__init__.py | 3 --- skfem/importers/meshio.py | 4 ++-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/skfem/importers/__init__.py b/skfem/importers/__init__.py index a08941c99..6f8efa2b7 100644 --- a/skfem/importers/__init__.py +++ b/skfem/importers/__init__.py @@ -1,6 +1,3 @@ import warnings -warnings.warn("DeprecationWarning: skfem.importers is removed " - "in the next release. Use skfem.io instead.") - from ..io import * diff --git a/skfem/importers/meshio.py b/skfem/importers/meshio.py index 101abc08a..37d0bf473 100644 --- a/skfem/importers/meshio.py +++ b/skfem/importers/meshio.py @@ -1,6 +1,6 @@ import warnings -warnings.warn("DeprecationWarning: skfem.importers is removed " - "in the next release. Use skfem.io instead.") +warnings.warn("skfem.importers is removed in the next release. " + "Use skfem.io instead.", DeprecationWarning) from ..io.meshio import * From e42b20a9d8c94f634610b9600f0031aeaeb2a09e Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sat, 18 Jan 2020 23:58:15 +0200 Subject: [PATCH 27/32] Improve dosctring --- skfem/assembly/basis/mortar_basis.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index ad9d1d535..8e10a2123 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -23,7 +23,8 @@ def __init__(self, elem An object of type :class:`~skfem.element.Element`. mapping - Mortar mapping. + Mapping to the relevant facets of the mesh, see + e.g. :class:`~skfem.mapping.MortarPair`. intorder Integration order, i.e. the degree of polynomials that are integrated exactly by the used quadrature. Please use equivalent From ac5e2f340e4076fdd5c737453c9120c64190cdcc Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sun, 19 Jan 2020 00:07:15 +0200 Subject: [PATCH 28/32] Improve docstring --- skfem/mapping/mortar_pair.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/skfem/mapping/mortar_pair.py b/skfem/mapping/mortar_pair.py index ca490499b..35aa49bb5 100644 --- a/skfem/mapping/mortar_pair.py +++ b/skfem/mapping/mortar_pair.py @@ -8,7 +8,17 @@ class MortarPair(NamedTuple): - """A pair of mappings for mortar methods.""" + """A pair of mappings for mortar methods. + + In mortar methods, we are interested in interface conditions on non-matching + meshes. We need to generate matching quadrature points on both sides of the + interface. + + The pair of mappings in MortarPair is defined so that when we map any local + point in the reference element using both of the mappings, we get matching + global points on the facets of the original non-matching meshes. + + """ mapping1: Mapping = None mapping2: Mapping = None @@ -25,11 +35,13 @@ def init_2D(cls, Parameters ---------- mesh1 + An object of the type :class:`~skfem.mesh.mesh_2d.Mesh2D`. mesh2 + An object of the type :class:`~skfem.mesh.mesh_2d.Mesh2D`. boundary1 - A subset of facets from mesh1. + A subset of facets to use from mesh1. boundary2 - A subset of facets from mesh2. + A subset of facets to use from mesh2. tangent A tangent vector defining the direction of the projection. From dbdd766d07eeaa3715f4b24ca35b04eaba452ece Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sun, 19 Jan 2020 14:35:48 +0200 Subject: [PATCH 29/32] Add tests for MortarPair --- skfem/assembly/basis/mortar_basis.py | 4 +- skfem/mapping/mortar_pair.py | 54 ++++++++++++++++---------- tests/test_assembly.py | 1 - tests/test_convergence.py | 3 -- tests/test_examples.py | 57 ++++++++-------------------- tests/test_manufactured.py | 1 + tests/test_mapping.py | 49 ++++++++++++++++++++++++ tests/test_utils.py | 1 - tests/test_visuals.py | 1 - 9 files changed, 102 insertions(+), 69 deletions(-) diff --git a/skfem/assembly/basis/mortar_basis.py b/skfem/assembly/basis/mortar_basis.py index 8e10a2123..fdc600b31 100644 --- a/skfem/assembly/basis/mortar_basis.py +++ b/skfem/assembly/basis/mortar_basis.py @@ -24,7 +24,7 @@ def __init__(self, An object of type :class:`~skfem.element.Element`. mapping Mapping to the relevant facets of the mesh, see - e.g. :class:`~skfem.mapping.MortarPair`. + :class:`~skfem.mapping.MortarPair`. intorder Integration order, i.e. the degree of polynomials that are integrated exactly by the used quadrature. Please use equivalent @@ -64,7 +64,7 @@ def default_parameters(self): return {'x': self.global_coordinates(), 'h': self.mesh_parameters(), 'n': self.normals} - + def global_coordinates(self) -> ndarray: return self.mapping.G(self.X) diff --git a/skfem/mapping/mortar_pair.py b/skfem/mapping/mortar_pair.py index 35aa49bb5..1c6633fd0 100644 --- a/skfem/mapping/mortar_pair.py +++ b/skfem/mapping/mortar_pair.py @@ -10,13 +10,21 @@ class MortarPair(NamedTuple): """A pair of mappings for mortar methods. - In mortar methods, we are interested in interface conditions on non-matching - meshes. We need to generate matching quadrature points on both sides of the - interface. + In mortar methods, we are enforcing interface conditions on non-matching + meshes. To this end, we must generate matching quadrature points on both + sides of the interface. - The pair of mappings in MortarPair is defined so that when we map any local + :class:`~skfem.mapping.MortarPair` is a pair of + :class:`~skfem.mapping.Mapping` objects defined so that mapping any local point in the reference element using both of the mappings, we get matching - global points on the facets of the original non-matching meshes. + global points on the original non-matching meshes. + + Attributes + ---------- + mapping1 + Mapping to the facets of the first mesh. + mapping2 + Mapping to the facets of the second mesh. """ @@ -29,7 +37,7 @@ def init_2D(cls, mesh2: Mesh2D, boundary1: ndarray, boundary2: ndarray, - tangent: ndarray): + tangent: ndarray = None): """Create mortar mappings for two 2D meshes via projection. Parameters @@ -44,10 +52,15 @@ def init_2D(cls, A subset of facets to use from mesh2. tangent A tangent vector defining the direction of the projection. + If not given, use first facet of boundary1 as the vector. """ from ..mesh import MeshLine + if tangent is None: + tangent = (mesh1.p[:, mesh1.facets[0, boundary1[0]]] - + mesh1.p[:, mesh1.facets[1, boundary1[0]]]) + tangent /= np.linalg.norm(tangent) # find unique nodes on the two boundaries @@ -73,7 +86,8 @@ def param(p): p = np.array([np.hstack((param(mesh1.p), param(mesh2.p)))]) t = np.array([ixorig[:-1], ixorig[1:]]) - # create 1-dimensional supermesh + # create 1-dimensional supermesh from the intersections of the projected + # facet elements p = p[:, np.concatenate((t[0], np.array([t[1, -1]])))] range_max = np.min([np.max(param_p1), np.max(param_p2)]) range_min = np.max([np.min(param_p1), np.min(param_p2)]) @@ -123,27 +137,27 @@ def param(p): m2.t[:, ix2_flip] = np.flipud(m2.t[:, ix2_flip]) # create new mapping objects with G replaced by the supermesh mapping - nmap_mesh1 = mesh1.mapping() - nmap_mesh2 = mesh2.mapping() + new_map1 = mesh1.mapping() + new_map2 = mesh2.mapping() map_m1 = m1.mapping() map_m2 = m2.mapping() - nmap_mesh1.G = lambda X: map_mesh1.G(map_m1.invF(map_super.F(X), tind=ix1), + new_map1.G = lambda X: map_mesh1.G(map_m1.invF(map_super.F(X), tind=ix1), find=boundary1[sort_boundary1][ix1]) - nmap_mesh2.G = lambda X: map_mesh2.G(map_m2.invF(map_super.F(X), tind=ix2), + new_map2.G = lambda X: map_mesh2.G(map_m2.invF(map_super.F(X), tind=ix2), find=boundary2[sort_boundary2][ix2]) # these are used by :class:`skfem.assembly.basis.MortarBasis` - nmap_mesh1.find = boundary1[sort_boundary1][ix1] - nmap_mesh2.find = boundary2[sort_boundary2][ix2] - nmap_mesh1.normals = normals - nmap_mesh2.normals = normals + new_map1.find = boundary1[sort_boundary1][ix1] + new_map2.find = boundary2[sort_boundary2][ix2] + new_map1.normals = normals + new_map2.normals = normals # method for calculating the lengths of the segments ('detDG') - segments1 = nmap_mesh1.G(np.array([[0.0, 1.0]])) - segments2 = nmap_mesh2.G(np.array([[0.0, 1.0]])) - nmap_mesh1.detDG = lambda *_,**__: np.sqrt(np.diff(segments1[0], axis=1)**2 + + segments1 = new_map1.G(np.array([[0.0, 1.0]])) + segments2 = new_map2.G(np.array([[0.0, 1.0]])) + new_map1.detDG = lambda *_,**__: np.sqrt(np.diff(segments1[0], axis=1)**2 + np.diff(segments1[1], axis=1)**2) - nmap_mesh2.detDG = lambda *_,**__: np.sqrt(np.diff(segments2[0], axis=1)**2 + + new_map2.detDG = lambda *_,**__: np.sqrt(np.diff(segments2[0], axis=1)**2 + np.diff(segments2[1], axis=1)**2) - return cls(mapping1=nmap_mesh1, mapping2=nmap_mesh2) + return cls(mapping1=new_map1, mapping2=new_map2) diff --git a/tests/test_assembly.py b/tests/test_assembly.py index 544d9b3fb..49b525522 100644 --- a/tests/test_assembly.py +++ b/tests/test_assembly.py @@ -1,4 +1,3 @@ -"""Unit tests for assembly operations.""" import unittest import numpy as np diff --git a/tests/test_convergence.py b/tests/test_convergence.py index 39ff41683..75c716db1 100644 --- a/tests/test_convergence.py +++ b/tests/test_convergence.py @@ -1,6 +1,3 @@ -# -*- coding: utf-8 -*- -"""Check the convergence rates of the elements.""" - import unittest import numpy as np diff --git a/tests/test_examples.py b/tests/test_examples.py index 1d4781f2e..63c902059 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -1,37 +1,27 @@ -"""These tests run the examples and check that their output stays constant.""" - import unittest import numpy as np class TestEx01(unittest.TestCase): - """Run examples/ex01.py""" - def runTest(self): import docs.examples.ex01 as ex01 self.assertAlmostEqual(np.max(ex01.x), 0.07344576657) class TestEx02(unittest.TestCase): - """Run examples/ex02.py""" - def runTest(self): import docs.examples.ex02 as ex02 self.assertAlmostEqual(np.max(ex02.x), 0.001217973811129439) class TestEx03(unittest.TestCase): - """Run examples/ex03.py""" - def runTest(self): import docs.examples.ex03 as ex03 self.assertAlmostEqual(ex03.L[0], 0.00418289) class TestEx04(unittest.TestCase): - """Run examples/ex04.py""" - def runTest(self): import docs.examples.ex04 as ex04 self.assertAlmostEqual(np.max(ex04.vonmises1), 64.57919892367978) @@ -39,56 +29,42 @@ def runTest(self): class TestEx05(unittest.TestCase): - """Run examples/ex05.py""" - def runTest(self): import docs.examples.ex05 as ex05 self.assertAlmostEqual(np.max(ex05.x), 0.93570751751091152) class TestEx06(unittest.TestCase): - """Run examples/ex06.py""" - def runTest(self): import docs.examples.ex06 as ex06 self.assertAlmostEqual(np.max(ex06.x), 0.073651530833125131) class TestEx07(unittest.TestCase): - """Run examples/ex07.py""" - def runTest(self): import docs.examples.ex07 as ex07 self.assertAlmostEqual(np.max(ex07.x), 0.07869083767545548) class TestEx08(unittest.TestCase): - """Run examples/ex08.py""" - def runTest(self): - pass + import docs.examples.ex08 as ex08 # only run the initialization, nothing to test class TestEx09(unittest.TestCase): - """Run examples/ex09.py""" - def runTest(self): import docs.examples.ex09 as ex09 self.assertAlmostEqual(np.max(ex09.x), 0.055596791644282988) class TestEx10(unittest.TestCase): - """Run examples/ex10.py""" - def runTest(self): import docs.examples.ex10 as ex10 self.assertAlmostEqual(np.mean(ex10.x), 0.277931521728906) class TestEx11(unittest.TestCase): - """Run examples/ex11.py""" - def runTest(self): import docs.examples.ex11 as ex11 u = ex11.u @@ -115,8 +91,6 @@ def runTest(self): class TestEx14(unittest.TestCase): - """Run examples/ex14.py""" - def runTest(self): import docs.examples.ex14 u = docs.examples.ex14.u @@ -125,35 +99,35 @@ def runTest(self): class TestEx15(unittest.TestCase): - """Run examples/ex15.py""" - def runTest(self): - import docs.examples.ex15 - self.assertTrue(np.max(docs.examples.ex15.x) - 0.1234567 < 1e-5) + import docs.examples.ex15 as ex15 + self.assertTrue(np.max(ex15.x) - 0.1234567 < 1e-5) class TestEx16(unittest.TestCase): - """Run examples/ex16.py""" - def runTest(self): - import docs.examples.ex16 - self.assertTrue(np.linalg.norm(np.array([0, 2, 6, 12, 20, 30]) - docs.examples.ex16.ks) < 0.4) - self.assertTrue(docs.examples.ex16.ks[-1], 30.309720458315521) + import docs.examples.ex16 as ex16 + self.assertTrue(np.linalg.norm(np.array([0, 2, 6, 12, 20, 30]) + - ex16.ks) < 0.4) + self.assertTrue(ex16.ks[-1], 30.309720458315521) class TestEx17(unittest.TestCase): def runTest(self): import docs.examples.ex17 as ex + # TODO improve class TestEx18(unittest.TestCase): def runTest(self): import docs.examples.ex18 as ex + # TODO improve class TestEx19(unittest.TestCase): def runTest(self): import docs.examples.ex19 as ex + # TODO improve class TestEx20(unittest.TestCase): @@ -187,22 +161,23 @@ def runTest(self): self.assertAlmostEqual(max(ex.lmbda_list), ex.turning_point, delta=5e-5) + class TestEx24(unittest.TestCase): def runTest(self): - pass + import docs.examples.ex24 as ex24 class TestEx25(unittest.TestCase): def runTest(self): - import docs.examples.ex25 as ex - mu = np.mean(ex.t) + import docs.examples.ex25 as ex25 + mu = np.mean(ex25.t) self.assertAlmostEqual(mu, 0.4642600944590631, places=5) - self.assertAlmostEqual(np.mean(ex.t0), mu, places=2) + self.assertAlmostEqual(np.mean(ex25.t0), mu, places=2) class TestEx26(unittest.TestCase): def runTest(self): - pass + import docs.examples.ex26 as ex26 class TestEx27(unittest.TestCase): diff --git a/tests/test_manufactured.py b/tests/test_manufactured.py index d9c6e0cd5..6fe667b2b 100644 --- a/tests/test_manufactured.py +++ b/tests/test_manufactured.py @@ -1,4 +1,5 @@ """Solve problems that have manufactured solutions.""" + import unittest import numpy as np diff --git a/tests/test_mapping.py b/tests/test_mapping.py index 44f3d1035..6b2d807ce 100644 --- a/tests/test_mapping.py +++ b/tests/test_mapping.py @@ -73,5 +73,54 @@ def initialize_meshes(self): return m +class TestMortarPair(unittest.TestCase): + """Check that mapped points match.""" + + mesh1_type = MeshTri + mesh2_type = MeshTri + nrefs1 = 2 + nrefs2 = 3 + translate_y = 0.0 + + def init_meshes(self): + m1 = self.mesh1_type() + m1.refine(self.nrefs1) + m2 = self.mesh2_type() + m2.refine(self.nrefs2) + m2.translate((1.0, self.translate_y)) + return m1, m2 + + def runTest(self): + m1, m2 = self.init_meshes() + mp = MortarPair.init_2D(m1, m2, + m1.facets_satisfying(lambda x: x[0] == 1.), + m2.facets_satisfying(lambda x: x[0] == 1.)) + test_points = np.array([np.linspace(0., 1., 7)]) + self.assertTrue((mp.mapping1.G(test_points) - + mp.mapping1.G(test_points) < 1e-10).all()) + + +class TestMortarPairTriQuad(TestMortarPair): + mesh1_type = MeshTri + mesh2_type = MeshQuad + + +class TestMortarPairQuadQuad(TestMortarPair): + mesh1_type = MeshQuad + mesh2_type = MeshQuad + + +class TestMortarPairNoMatch1(TestMortarPair): + mesh1_type = MeshQuad + mesh2_type = MeshTri + translate_y = 0.1 + + +class TestMortarPairNoMatch2(TestMortarPair): + mesh1_type = MeshQuad + mesh2_type = MeshTri + translate_y = -np.pi / 10. + + if __name__ == '__main__': unittest.main() diff --git a/tests/test_utils.py b/tests/test_utils.py index afab46d2f..d482bddc3 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,4 +1,3 @@ -"""Unit tests for utils package.""" import unittest import numpy as np diff --git a/tests/test_visuals.py b/tests/test_visuals.py index e36868f64..cebd1873f 100644 --- a/tests/test_visuals.py +++ b/tests/test_visuals.py @@ -1,6 +1,5 @@ import unittest - from skfem.mesh import * from skfem.visuals.matplotlib import * From b8ed6939d3f0fb0d7ebe801848e7f35507ad7ea0 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Sun, 19 Jan 2020 14:45:11 +0200 Subject: [PATCH 30/32] Update changelog --- CHANGELOG.md | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e3261e0d5..0e9008434 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,14 +4,20 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Changed +- Renamed skfem.importers to skfem.io. + +### Added +- MortarBasis and MortarPair that support mortar methods in 2D. +- Serialization of meshes to json-files via skfem.io.json. +- ElementQuadDG for transforming H1 elements to DG elements. + ## [0.4.1] - 2020-01-19 ### Added - Additional keyword arguments to skfem.utils.solve get passed on to solvers. -- Add MortarBasis and MortarPair for mortaring in 2D. - -### Changed -- Renamed skfem.importers to skfem.io. ### Fixed - Made skfem.visuals.matplotlib Python 3.6 compatible. From 531aea884b033546dcb103076ee8cb40bdb89e8f Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 20 Jan 2020 12:31:51 +0200 Subject: [PATCH 31/32] Add ex04 to docs and lots of small improvements --- docs/api.rst | 2 + docs/conf.py | 2 +- docs/{learning.rst => examples.rst} | 11 +-- docs/examples/Makefile | 1 + docs/examples/ex01.py | 4 +- docs/examples/ex04.py | 19 ++--- docs/examples/ex04.rst | 69 +++++++++++++++++++ docs/examples/ex06.rst | 8 +++ docs/index.rst | 2 +- docs/tutorial.rst | 2 +- .../ex01.rst => tutorials/poisson.rst} | 16 ++--- docs/tutorials/walkthrough.rst | 6 +- 12 files changed, 113 insertions(+), 29 deletions(-) rename docs/{learning.rst => examples.rst} (97%) create mode 100644 docs/examples/ex04.rst rename docs/{examples/ex01.rst => tutorials/poisson.rst} (89%) diff --git a/docs/api.rst b/docs/api.rst index f9503afb3..f9983c6ca 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -169,6 +169,8 @@ Defining a global basis .. autoclass:: skfem.assembly.InteriorBasis +.. automethod:: skfem.assembly.InteriorBasis.refinterp + .. autoclass:: skfem.assembly.FacetBasis Assembling matrices diff --git a/docs/conf.py b/docs/conf.py index fc5b10195..e8e3b6326 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -19,7 +19,7 @@ # -- Project information ----------------------------------------------------- project = 'scikit-fem' -copyright = '2018-2019, Tom Gustafsson, Geordie McBain' +copyright = '2018-20, Tom Gustafsson, Geordie McBain' author = 'Tom Gustafsson, Geordie McBain' diff --git a/docs/learning.rst b/docs/examples.rst similarity index 97% rename from docs/learning.rst rename to docs/examples.rst index f7432a536..9dccffffa 100644 --- a/docs/learning.rst +++ b/docs/examples.rst @@ -15,13 +15,14 @@ Poisson equation examples/ex14.rst examples/ex12.rst examples/ex13.rst - examples/ex15.rst examples/ex26.rst examples/ex09.rst examples/ex07.rst examples/ex05.rst examples/ex10.rst examples/ex22.rst + examples/ex15.rst + examples/ex06.rst Heat transfer ############# @@ -51,11 +52,12 @@ Structural mechanics .. toctree:: - examples/ex03.rst - examples/ex11.rst - examples/ex08.rst examples/ex02.rst examples/ex21.rst + examples/ex04.rst + examples/ex11.rst + examples/ex03.rst + examples/ex08.rst Others ###### @@ -64,4 +66,3 @@ Others examples/ex23.rst examples/ex16.rst - examples/ex06.rst diff --git a/docs/examples/Makefile b/docs/examples/Makefile index d71a6dcba..4ff37034d 100644 --- a/docs/examples/Makefile +++ b/docs/examples/Makefile @@ -1,5 +1,6 @@ all: ex01_solution.png \ ex02_solution.png \ + ex04_solution.png \ ex05_solution.png \ ex06_solution.png \ ex17_solution.png \ diff --git a/docs/examples/ex01.py b/docs/examples/ex01.py index ff85d171f..a45a75668 100644 --- a/docs/examples/ex01.py +++ b/docs/examples/ex01.py @@ -8,11 +8,11 @@ @bilinear_form def laplace(u, du, v, dv, w): - return du[0]*dv[0] + du[1]*dv[1] + return du[0] * dv[0] + du[1] * dv[1] @linear_form def load(v, dv, w): - return 1.0*v + return 1. * v A = asm(laplace, basis) b = asm(load, basis) diff --git a/docs/examples/ex04.py b/docs/examples/ex04.py index 05fb8931e..474632caf 100644 --- a/docs/examples/ex04.py +++ b/docs/examples/ex04.py @@ -1,17 +1,16 @@ from skfem import * from skfem.models.elasticity import linear_elasticity import numpy as np -from pygmsh import generate_mesh -from pygmsh.built_in import Geometry from skfem.io import from_meshio from skfem.io.json import from_file, to_file -import meshio # create meshes try: m = from_file("docs/examples/ex04_mesh.json") except FileNotFoundError: + from pygmsh import generate_mesh + from pygmsh.built_in import Geometry geom = Geometry() points = [] lines = [] @@ -109,7 +108,8 @@ def bilin_penalty(u, du, v, dv, w): n[1] * C(Eps(dv))[1, 0] * n[0] + n[1] * C(Eps(dv))[1, 1] * n[1]) h = w.h - return (1. / (alpha * h) * ju * jv - mu * jv - mv * ju) *(np.abs(w.x[1]) <= limit) + return ((1. / (alpha * h) * ju * jv - mu * jv - mv * ju) + * (np.abs(w.x[1]) <= limit)) K[j][i] += asm(bilin_penalty, mb[i], mb[j]) @@ -124,7 +124,8 @@ def lin_penalty(v, dv, w): h = w.h def gap(x): return (1. - np.sqrt(1. - x[1] ** 2)) - return (1. / (alpha * h) * gap(w.x) * jv - gap(w.x) * mv) * (np.abs(w.x[1]) <= limit) + return ((1. / (alpha * h) * gap(w.x) * jv - gap(w.x) * mv) + * (np.abs(w.x[1]) <= limit)) f[i] = asm(lin_penalty, mb[i]) @@ -181,8 +182,10 @@ def mass(u, du, v, dv, w): ib_dg = InteriorBasis(m, e_dg, intorder=4) Ib_dg = InteriorBasis(M, E_dg, intorder=4) - s[itr, jtr] = solve(asm(mass, ib_dg), asm(proj_cauchy, ib, ib_dg) @ x[i1]) - S[itr, jtr] = solve(asm(mass, Ib_dg), asm(proj_cauchy, Ib, Ib_dg) @ x[i2]) + s[itr, jtr] = solve(asm(mass, ib_dg), + asm(proj_cauchy, ib, ib_dg) @ x[i1]) + S[itr, jtr] = solve(asm(mass, Ib_dg), + asm(proj_cauchy, Ib, Ib_dg) @ x[i2]) s[2, 2] = nu1 * (s[0, 0] + s[1, 1]) S[2, 2] = nu2 * (S[0, 0] + S[1, 1]) @@ -207,4 +210,4 @@ def mass(u, du, v, dv, w): draw(m, ax=ax) plot(Ib_dg, vonmises2, ax=ax, Nrefs=3, shading='gouraud') draw(M, ax=ax) - savefig(splitext(argv[0])[0] + '_vonmises.png') + savefig(splitext(argv[0])[0] + '_solution.png') diff --git a/docs/examples/ex04.rst b/docs/examples/ex04.rst new file mode 100644 index 000000000..add352e97 --- /dev/null +++ b/docs/examples/ex04.rst @@ -0,0 +1,69 @@ +.. _mortaring: + +Contact problem +--------------- + +.. note:: + + This example requires the external package `pygmsh + `_. + +Mortar methods allow setting interface conditions on non-matching meshes. +They are useful also when solving variational inequalities such as +`elastic contact problems `_. + +This example solves the first contact iteration for the following prototype +contact problem: find :math:`\boldsymbol{u}_i : \Omega_i \rightarrow +\mathbb{R}^2`, :math:`i = 1,2`, such that + + +.. math:: + + \begin{aligned} + \boldsymbol{\mathrm{div}}\,\boldsymbol{\sigma}_i(\boldsymbol{u}_i)&=\boldsymbol{0} \quad && \text{in $\Omega_i$,} \\ + \boldsymbol{u}_1&=(0.1, 0) \quad && \text{on $\Gamma_{D,1}$,} \\ + \boldsymbol{u}_2&=\boldsymbol{0} \quad && \text{on $\Gamma_{D,2}$,} \\ + \boldsymbol{\sigma}_2(\boldsymbol{u}_2) \boldsymbol{n}_2 &=\boldsymbol{0} \quad && \text{on $\Gamma_{N,2}$,} \\ + \boldsymbol{\sigma}_{i,t}(\boldsymbol{u}_i) &= \boldsymbol{0} && \text{on $\Gamma$,} \\ + \sigma_{1,n}(\boldsymbol{u}_1(\boldsymbol{\gamma}(\boldsymbol{x})) - \sigma_{2,n}(\boldsymbol{u}_2)&=0 && \text{on $\Gamma$,} \\ + [[u_n]] - g &\geq 0 && \text{on $\Gamma$,} \\ + \sigma_{i,n}(\boldsymbol{u}_i)&\leq 0 && \text{on $\Gamma$,} \\ + ([[u_n]] - g) \sigma_{i,n}(\boldsymbol{u}_i) &= 0 && \text{on $\Gamma$,} + \end{aligned} + + +where + +* :math:`\Omega_1 = \{ (x, y) : x^2 + y^2 < 1 \} \setminus \{ (x, y) : x < 0\}`, +* :math:`\Omega_2 = (1, 2) \times (-1, 1)`, +* :math:`\Gamma_{D,1} = \{ (x, y) \in \Omega_1 : x=0 \}`, +* :math:`\Gamma_{D,2} = \{ (x, y) \in \Omega_2 : x=2 \}`, +* :math:`\Gamma = \{ (x, y) \in \Omega_2 : x=1 \}`, +* :math:`g((x,y)) = 1 - \sqrt{1 - y^2}`, +* :math:`\boldsymbol{\gamma} : \Gamma \rightarrow \{ (x, y) \in \partial + \Omega_1 : x > 0 \}`, :math:`\boldsymbol{\gamma}((x,y)) = (g(x-1)+1, y)`, + +and the directions for evaluating :math:`\sigma_{1,n}` and +:math:`\boldsymbol{\sigma}_{1,t}` are defined as :math:`\boldsymbol{n}=(1,0)` +and :math:`\boldsymbol{t}=(0,1)`. +This is a nonlinear problem since we do not know a priori which subset +:math:`\Gamma_C \subset \Gamma` satisfies :math:`([[u_n]] - g)|_{\Gamma_C} = 0`. + +.. note:: + + The example solves a simplified prototype problem. + Instead of iterating for the true contact boundary, + we solve a single contact iteration (a linear problem) with the initial + guess :math:`\{ (x, y) \in \Gamma : |y| < 0.1 \}`. + Solving a real contact problem involves repeatedly solving and guessing a new + candidate boundary :math:`\Gamma_C` until convergence. + Extending this example should be straightforward. + +.. figure:: ex04_solution.png + + The von Mises stress in a displaced geometry. + +The complete source code reads as follows: + +.. literalinclude:: ex04.py + :linenos: diff --git a/docs/examples/ex06.rst b/docs/examples/ex06.rst index f4188d899..6976cbef0 100644 --- a/docs/examples/ex06.rst +++ b/docs/examples/ex06.rst @@ -15,6 +15,14 @@ interpolates any solution vector. The resulting mesh is non-conforming, i.e. the connectivity between neighboring elements is lost, and hence it can be used only for visualisation purposes. +.. note:: + + As of 0.4.0, this functionality is included in :func:`~skfem.visuals.matplotlib.plot`, + i.e. inputting :class:`~skfem.assembly.InteriorBasis` instead of :class:`~skfem.mesh.Mesh` + uses :meth:`~skfem.assembly.InteriorBasis.refinterp` automatically. + The steps in this example are still useful when, e.g., exporting to different + formats for visualization purposes. + As an example, we solve the Poisson equation in a unit square with zero boundary conditions and biquadratic basis on quadrilateral elements. The quadrilateral elements are defined using an isoparametric local-to-global mapping. diff --git a/docs/index.rst b/docs/index.rst index 67284fa28..ef75f5ce7 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -19,5 +19,5 @@ questions, do not hesitate to drop in and say hello at our `Gitter chat gettingstarted tutorial - learning + examples api diff --git a/docs/tutorial.rst b/docs/tutorial.rst index b0e1c1b87..79420ba13 100644 --- a/docs/tutorial.rst +++ b/docs/tutorial.rst @@ -6,4 +6,4 @@ Tutorials .. toctree:: tutorials/walkthrough.rst - examples/ex01.rst + tutorials/poisson.rst diff --git a/docs/examples/ex01.rst b/docs/tutorials/poisson.rst similarity index 89% rename from docs/examples/ex01.rst rename to docs/tutorials/poisson.rst index d80f6cfa6..02ba586f9 100644 --- a/docs/examples/ex01.rst +++ b/docs/tutorials/poisson.rst @@ -39,7 +39,7 @@ times: After creating the mesh, we evaluate the finite element basis at the global quadrature points. -.. literalinclude:: ex01.py +.. literalinclude:: ../examples/ex01.py :start-at: ElementTriP1 :end-at: InteriorBasis @@ -54,13 +54,13 @@ The bilinear and linear forms are defined using the decorators :func:`~skfem.assembly.linear_form`. It is important to have the order of the form arguments correct. -.. literalinclude:: ex01.py +.. literalinclude:: ../examples/ex01.py :start-at: bilinear_form - :end-at: 1.0*v + :end-at: 1. * v All assembly operations are performed using the function :func:`~skfem.assembly.asm`. -.. literalinclude:: ex01.py +.. literalinclude:: ../examples/ex01.py :start-at: asm(laplace, basis) :end-at: asm(load, basis) @@ -71,21 +71,21 @@ We are left with solving the assembled linear system. We eliminate the boundary degrees-of-freedom using :func:`~skfem.utils.condense` and call :func:`~skfem.utils.solve`. -.. literalinclude:: ex01.py +.. literalinclude:: ../examples/ex01.py :start-at: solve :end-at: solve The solution can now be visualised using :meth:`~skfem.visuals.matplotlib.plot`. -.. literalinclude:: ex01.py +.. literalinclude:: ../examples/ex01.py :start-at: main :end-at: solution.png -.. figure:: ex01_solution.png +.. figure:: ../examples/ex01_solution.png The solution of Poisson equation. The complete source code reads as follows: -.. literalinclude:: ex01.py +.. literalinclude:: ../examples/ex01.py :linenos: diff --git a/docs/tutorials/walkthrough.rst b/docs/tutorials/walkthrough.rst index 38c9fdcf9..8b1ea5da9 100644 --- a/docs/tutorials/walkthrough.rst +++ b/docs/tutorials/walkthrough.rst @@ -1,7 +1,7 @@ .. _basic-features: -Basic features --------------- +Overview of basic features +-------------------------- This is an overview of the basic features of scikit-fem that are needed in most finite element computations. @@ -136,7 +136,7 @@ See :ref:`learning` for more use cases and instructions. Setting essential boundary conditions ##################################### -.. note:: +.. warning:: Using the assembled matrices requires basic understanding of the finite element method. In particular, to understand From e556b8e80359c19a66b6c017e600b5945aab7484 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Mon, 20 Jan 2020 12:46:07 +0200 Subject: [PATCH 32/32] Add the definition of jump --- docs/examples/ex04.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/examples/ex04.rst b/docs/examples/ex04.rst index add352e97..0a7b3fb82 100644 --- a/docs/examples/ex04.rst +++ b/docs/examples/ex04.rst @@ -27,8 +27,8 @@ contact problem: find :math:`\boldsymbol{u}_i : \Omega_i \rightarrow \boldsymbol{\sigma}_{i,t}(\boldsymbol{u}_i) &= \boldsymbol{0} && \text{on $\Gamma$,} \\ \sigma_{1,n}(\boldsymbol{u}_1(\boldsymbol{\gamma}(\boldsymbol{x})) - \sigma_{2,n}(\boldsymbol{u}_2)&=0 && \text{on $\Gamma$,} \\ [[u_n]] - g &\geq 0 && \text{on $\Gamma$,} \\ - \sigma_{i,n}(\boldsymbol{u}_i)&\leq 0 && \text{on $\Gamma$,} \\ - ([[u_n]] - g) \sigma_{i,n}(\boldsymbol{u}_i) &= 0 && \text{on $\Gamma$,} + \sigma_{2,n}(\boldsymbol{u}_2)&\leq 0 && \text{on $\Gamma$,} \\ + ([[u_n]] - g) \sigma_{2,n}(\boldsymbol{u}_2) &= 0 && \text{on $\Gamma$,} \end{aligned} @@ -42,8 +42,9 @@ where * :math:`g((x,y)) = 1 - \sqrt{1 - y^2}`, * :math:`\boldsymbol{\gamma} : \Gamma \rightarrow \{ (x, y) \in \partial \Omega_1 : x > 0 \}`, :math:`\boldsymbol{\gamma}((x,y)) = (g(x-1)+1, y)`, +* :math:`[[u_n]] = \boldsymbol{u}_1(\boldsymbol{\gamma}(\boldsymbol{x})) \cdot \boldsymbol{n} - \boldsymbol{u}_2(\boldsymbol{x}) \cdot \boldsymbol{n}`, -and the directions for evaluating :math:`\sigma_{1,n}` and +and the directions for evaluating :math:`[[u_n]]`, :math:`\sigma_{1,n}` and :math:`\boldsymbol{\sigma}_{1,t}` are defined as :math:`\boldsymbol{n}=(1,0)` and :math:`\boldsymbol{t}=(0,1)`. This is a nonlinear problem since we do not know a priori which subset