add new example

docs/examples/ex51.py

@@ -0,0 +1,130 @@
+"""Contact problem."""
+from skfem import *
+from skfem.experimental.autodiff import *
+from skfem.experimental.autodiff.helpers import *
+from skfem.supermeshing import intersect, elementwise_quadrature
+import jax.numpy as jnp
+import numpy as np
+
+
+m1 = (MeshQuad
+      .init_tensor(np.linspace(0, 5, 60), np.linspace(0, 0.25, 5))
+      .with_defaults())
+m2 = (MeshQuad
+      .init_tensor(np.linspace(0, 5, 51), np.linspace(-0.5, -0.25, 4))
+      .with_defaults())
+
+e1 = ElementVector(ElementQuad1())
+e2 = ElementVector(ElementQuad1())
+
+basis1 = Basis(m1, e1)
+basis2 = Basis(m2, e2)
+
+# forms are based on the energy minimization problem
+
+@NonlinearForm(hessian=True)
+def g(u1, u2, w):
+    eps = 1e-2
+    return 1 / eps * jnp.minimum(0.25 + u1[1] - u2[1], 0) ** 2
+
+
+@NonlinearForm(hessian=True)
+def J1(u, w):
+    eps = 1/2 * (grad(u) + transpose(grad(u)) + mul(transpose(grad(u)), grad(u)))
+    sig = 2 * eps + eye(trace(eps), 2)
+    return 1/2 * ddot(sig, eps)
+
+
+@NonlinearForm(hessian=True)
+def J2(u, w):
+    eps = 1/2 * (grad(u) + transpose(grad(u)) + mul(transpose(grad(u)), grad(u)))
+    sig = 2 * eps + eye(trace(eps), 2)
+    # apply body force for domain 2
+    return 1/2 * ddot(sig, eps) - 1e-2 * u[1] - 5e-2 * u[0]
+
+
+x = np.zeros(basis1.N + basis2.N)
+m1defo = m1.copy()
+m2defo = m2.copy()
+
+for itr in range(40):
+
+    # use deformed mesh for mapping between domain 1 and 2
+    m1t, orig1 = m1defo.trace('bottom',
+                              mtype=MeshLine,
+                              project=lambda p: np.array(p[0]))
+    m2t, orig2 = m2defo.trace('top',
+                              mtype=MeshLine,
+                              project=lambda p: np.array(p[0]))
+    m12, t1, t2 = intersect(m1t, m2t)
+
+    fbasis1 = FacetBasis(m1, e1,
+                         quadrature=elementwise_quadrature(m1t, m12, t1),
+                         facets=orig1[t1])
+    fbasis2 = FacetBasis(m2, e2,
+                         quadrature=elementwise_quadrature(m2t, m12, t2),
+                         facets=orig2[t2])
+    fbasis = fbasis1 * fbasis2
+
+    # assemble jacobians for 1 and 2 separately and create a block matrix
+    jac1, rhs1 = J1.assemble(basis1, x=x[:basis1.N])
+    jac2, rhs2 = J2.assemble(basis2, x=x[basis1.N:])
+    jac = bmat([[jac1, None],
+                [None, jac2]])
+    rhs = np.concatenate((rhs1, rhs2))
+
+    # g assembled using fbasis1 * fbasis2 has automatically the correct shape
+    jacg, rhsg = g.assemble(fbasis, x=x)
+    jac = jac + jacg
+    rhs = rhs + rhsg
+
+    dx = solve(*enforce(jac, np.array(rhs), D=np.concatenate((
+        basis1.get_dofs({'left', 'right'}).all(),
+        basis2.get_dofs({'left', 'right'}).all() + basis1.N,
+    ))))
+    x += 0.95 * dx  # regularization
+    print(np.linalg.norm(dx))
+    if np.linalg.norm(dx) < 1e-8:
+        break
+    (x1, _), (x2, _) = fbasis.split(x)
+    m1defo = m1.translated(x1[basis1.nodal_dofs])
+    m2defo = m2.translated(x2[basis2.nodal_dofs])
+
+# postprocessing: calculate stress and visualize
+
+# calculate stress
+u1 = basis1.interpolate(x1)
+eps1 = 1/2 * (grad(u1) + transpose(grad(u1)) + mul(transpose(grad(u1)), grad(u1)))
+sig1 = 2 * eps1 + eye(trace(eps1), 2)
+
+u2 = basis2.interpolate(x2)
+eps2 = 1/2 * (grad(u2) + transpose(grad(u2)) + mul(transpose(grad(u2)), grad(u2)))
+sig2 = 2 * eps2 + eye(trace(eps2), 2)
+
+pbasis1 = basis1.with_element(ElementQuad0())
+pbasis2 = basis2.with_element(ElementQuad0())
+
+# basis with deformed mesh for convenient plotting
+pbasis1defo = Basis(m1defo, ElementQuad0())
+pbasis2defo = Basis(m2defo, ElementQuad0())
+
+sigyy1 = pbasis1.project(np.array(sig1[1, 1]))
+sigyy2 = pbasis2.project(np.array(sig2[1, 1]))
+
+if __name__ == "__main__":
+
+    ax = m1defo.draw()
+    m2defo.draw(ax=ax)
+    ax.set_xlim(0.5, 4.5)
+    pbasis1defo.plot(pbasis1.project(np.array(sig1[1, 1])),
+                     vmax=0.0,
+                     vmin=-0.01,
+                     ax=ax,
+                     cmap='viridis')
+    pbasis2defo.plot(pbasis2.project(np.array(sig2[1, 1])),
+                     vmax=0.0,
+                     vmin=-0.01,
+                     ax=ax,
+                     colorbar={'orientation': 'horizontal',
+                               'label': r'$\sigma_{yy}$'},
+                     cmap='viridis').show()

docs/listofexamples.rst

@@ -341,6 +341,19 @@ of a hyperelastic Neo-Hookean solid.
 
 See the source code of :exlink:`43` for more information.
 
+Example 51: Contact problem
+---------------------------
+
+This example solves the so-called large deformation
+contact problem between
+two nonlinear elastic solids.
+
+.. figure:: https://private-user-images.githubusercontent.com/973268/349649283-7527098e-b511-4f41-a359-d590a4479039.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjEyNDUyNDUsIm5iZiI6MTcyMTI0NDk0NSwicGF0aCI6Ii85NzMyNjgvMzQ5NjQ5MjgzLTc1MjcwOThlLWI1MTEtNGY0MS1hMzU5LWQ1OTBhNDQ3OTAzOS5wbmc_WC1BbXotQWxnb3JpdGhtPUFXUzQtSE1BQy1TSEEyNTYmWC1BbXotQ3JlZGVudGlhbD1BS0lBVkNPRFlMU0E1M1BRSzRaQSUyRjIwMjQwNzE3JTJGdXMtZWFzdC0xJTJGczMlMkZhd3M0X3JlcXVlc3QmWC1BbXotRGF0ZT0yMDI0MDcxN1QxOTM1NDVaJlgtQW16LUV4cGlyZXM9MzAwJlgtQW16LVNpZ25hdHVyZT0yZDhlZGI5NmMzOGI1MTZhMDQ1OTUyZWUxNDg5N2MyYTUzMDAzNmU3Y2FkZWQ0ZWQ1NTM0NzkxYWMzOTEzZTg1JlgtQW16LVNpZ25lZEhlYWRlcnM9aG9zdCZhY3Rvcl9pZD0wJmtleV9pZD0wJnJlcG9faWQ9MCJ9.GVSzG1xLmxGgi9GxoE1IwGKAhuMEdcFo17WB7zpMZ1w
+
+   The deformed mesh and stress of Example 51.
+
+See the source code of :exlink:`51` for more information.
+
 Fluid mechanics
 ===============
 

skfem/mesh/mesh.py

@@ -845,6 +845,9 @@ def init_refdom(cls):
         """Initialize a mesh corresponding to the reference domain."""
         return cls(cls.elem.refdom.p, cls.elem.refdom.t, validate=False)
 
+    def copy(self):
+        return replace(self)
+
     def morphed(self, *args):
         """Morph the mesh using functions.
 

tests/test_examples.py

@@ -416,3 +416,13 @@ def runTest(self):
         import docs.examples.ex50 as ex50
         self.assertAlmostEqual(ex50.u.max(),
                                1.0026836382652844)
+
+
+class TestEx51(TestCase):
+
+    def runTest(self):
+        import docs.examples.ex51 as ex51
+        self.assertAlmostEqual(ex51.x1.max(),
+                               0.15983570088457588)
+        self.assertAlmostEqual(ex51.x2.max(),
+                               0.4062806125571157)