Merge pull request #305 from kinnala/mortar-methods-fix

Mortar methods generalizations and fixes

CHANGELOG.md

@@ -4,6 +4,16 @@ 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

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

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'
 
 

docs/learning.rst → 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

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 \

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)

docs/examples/ex04.py

@@ -1,133 +1,213 @@
-"""
-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
 import numpy as np
-
-m = MeshTri()
-m.refine(3)
-M = MeshLine(np.linspace(0, 1, 6))
-M = M*M
-M = M._splitquads()
+from skfem.io import from_meshio
+from skfem.io.json import from_file, to_file
+
+
+# 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 = []
+    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))
+    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))
+M.refine()
 
-map = MappingAffine(m)
-Map = MappingAffine(M)
-e1 = ElementTriP1()
-e = ElementVectorH1(e1)
 
-ib = InteriorBasis(m, e, map, 2)
-Ib = InteriorBasis(M, e, Map, 2)
+# define elements and bases
+e1 = ElementTriP2()
+e = ElementVectorH1(e1)
 
-def rule(x, y):
-    return (x==1.0)
+E1 = ElementQuad2()
+E = ElementVectorH1(E1)
 
-def param(x, y):
-    return y
+ib = InteriorBasis(m, e, intorder=4)
+Ib = InteriorBasis(M, E, intorder=4)
 
-mortar = InterfaceMesh1D(m, M, rule, param)
+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 = {}
-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(m, e, mapping = mappings[0], intorder=4),
+    MortarBasis(M, E, mapping = mappings[1], intorder=4),
+]
 
+# define bilinear forms
 E1 = 1000.0
 E2 = 1000.0
 
 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=Mu1)
-weakform2 = linear_elasticity(Lambda=Lambda2, Mu=Mu2)
+weakform1 = linear_elasticity(Lambda=Lambda, Mu=Mu)
+weakform2 = linear_elasticity(Lambda=Lambda, Mu=Mu)
+
 
-alpha = 1
+alpha = 1000
+limit = 0.3
+
+# assemble the stiffness matrices
 K1 = asm(weakform1, ib)
 K2 = asm(weakform2, Ib)
-L = 0
+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)]])
+
+def Eps(dw):
+    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
         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])
-
-            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])
+            ju = (-1.) ** i * (u[0] * n[0] + u[1] * n[1])
+            jv = (-1.) ** j * (v[0] * n[0] + v[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 = asm(bilin_penalty, mb[i], mb[j]) + L
-
-@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])
 
 import scipy.sparse
-K = (scipy.sparse.bmat([[K1, None],[None, K2]]) + L).tocsr()
+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]
 
-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])
 
-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))
 I = np.setdiff1d(np.arange(K.shape[0]), D)
 
-x[I] = 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))
+
+
+# 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, :]]
 
-if __name__ == "__main__":
 
-    from skfem.visuals.matplotlib import draw, show
+# post processing
+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. * 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__":
+    from os.path import splitext
+    from sys import argv
+    from skfem.visuals.matplotlib import *
 
-    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] + '_solution.png')