Add nutils solid participant to perpendicular flap (#431)

* add nutils solid participant
* Add new nutils solver to README list

---------

Co-authored-by: Benjamin Uekermann <benjamin.uekermann@ipvs.uni-stuttgart.de>

perpendicular-flap/README.md

@@ -37,6 +37,8 @@ Solid participant:
 
 * DUNE. For more information, have a look at the [experimental DUNE adapter](https://github.com/precice/dune-adapter) and send us your feedback.
 
+* Nutils. The structural model is currently limited to linear elasticity. For more information, have a look at the [Nutils adapter documentation](https://www.precice.org/adapter-nutils.html). This Nutils solver requires at least Nutils v8.0.
+
 * OpenFOAM (solidDisplacementFoam). For more information, have a look at the [OpenFOAM plateHole tutorial](https://www.openfoam.com/documentation/tutorial-guide/5-stress-analysis/5.1-stress-analysis-of-a-plate-with-a-hole). The solidDisplacementFoam solver only supports linear geometry. For general solid mechanics procedures in OpenFOAM, see solids4foam.
 
 * solids4foam. Like for CalculuX, the geometrically linear solver is used by default. For more information, see the [solids4foam documentation](https://bitbucket.org/philip_cardiff/solids4foam-release/src/master/documentation/overview.md) This case currently (November 2022) only works with the `nextRelease` branch of solids4foam, which is only compatible with up to OpenFOAM v2012 and OpenFOAM 9 (as well as foam-extend, with which the OpenFOAM-preCICE adapter is not compatible), as well as the OpenFOAM-preCICE adapter v1.2.0 or later. Note that, since both solids4foam and preCICE rely on Eigen, and due to some [solids4foam-preCICE compatibility issue](https://github.com/precice/openfoam-adapter/issues/238), you need to build preCICE and solids4foam with the same Eigen version.

perpendicular-flap/solid-nutils/clean.sh

@@ -0,0 +1,6 @@
+#!/bin/sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_nutils .

perpendicular-flap/solid-nutils/run.sh

@@ -0,0 +1,4 @@
+#!/bin/sh
+set -e -u
+
+python3 solid.py richoutput=no

perpendicular-flap/solid-nutils/solid.py

@@ -0,0 +1,102 @@
+#! /usr/bin/env python3
+
+from nutils import mesh, function, solver, export, cli
+from nutils.expression_v2 import Namespace
+import numpy
+import treelog
+import precice
+
+
+def main(young=4e6, density=3e3, poisson=.3, nelems=2, timestepsize=0.01, npoints_per_elem=3):
+
+    topo, geom = mesh.rectilinear([numpy.linspace(-.05, .05, nelems+1), numpy.linspace(0, 1, 10*nelems+1)])
+    wall = topo.boundary['left,top,right'].sample('uniform', npoints_per_elem)
+
+    ns = Namespace()
+    ns.X = geom
+    ns.δ = function.eye(2)
+    ns.define_for('X', gradient='∇', normal='n', jacobians=('dV', 'dS'))
+    ns.add_field('dt')
+    ns.add_field(('u', 'u0', 'testu', 'v', 'v0', 'testv'), topo.basis('std', degree=1), shape=(2,))
+    ns.add_field('F', wall.basis(), shape=(2,))
+    ns.qw = 1 / npoints_per_elem # quadrature weight
+    ns.t_i = 'F_i / qw dS'
+    ns.dudt_i = '(u_i - u0_i) / dt'
+    ns.dvdt_i = '(v_i - v0_i) / dt'
+    ns.ρ = density
+    ns.λ = young * poisson / ((1 + poisson) * (1 - 2 * poisson))
+    ns.μ = young / (2 * (1 + poisson))
+    ns.σ_ij = 'λ ∇_k(u_k) δ_ij + μ (∇_j(u_i) + ∇_i(u_j))'
+
+    # make sure we correctly scale point forces to tractions
+    testforce = numpy.random.normal(size=(wall.npoints, 2))
+    numpy.testing.assert_almost_equal(
+        actual=wall.integrate('t_i dS' @ ns, F=testforce),
+        desired=testforce.sum(0),
+        decimal=10,
+        err_msg='nutils error: failed to recover net force')
+
+    # continuum equations: ρ v' = ∇·σ + F, u' = v
+    res = topo.integral('testv_i (dudt_i - v_i) dV' @ ns, degree=2)
+    res += topo.integral('(testu_i ρ dvdt_i + ∇_j(testu_i) σ_ij) dV' @ ns, degree=2)
+    res -= wall.integral('testu_i t_i dS' @ ns)
+
+    # boundary conditions: fully constrained at y=0
+    sqr = topo.boundary['bottom'].integral('u_k u_k' @ ns, degree=2)
+    cons = solver.optimize('u,', sqr, droptol=1e-10)
+
+    # initial conditions: undeformed and unmoving
+    sqr = topo.integral('u_k u_k + v_k v_k' @ ns, degree=2)
+    arguments = solver.optimize('u,v', sqr, constrain=cons)
+
+    # preCICE setup
+    solverProcessIndex = 0
+    solverProcessSize = 1
+    interface = precice.Interface("Solid", "../precice-config.xml", solverProcessIndex, solverProcessSize)
+    meshID = interface.get_mesh_id("Solid-Mesh")
+    dataIndices = interface.set_mesh_vertices(meshID, wall.eval(ns.X))
+    writedataID = interface.get_data_id("Displacement", meshID)
+    readdataID = interface.get_data_id("Force", meshID)
+
+    # initialize preCICE
+    precice_dt = interface.initialize()
+
+    timestep = 0
+    force = numpy.zeros((wall.npoints, 2))
+
+    while interface.is_coupling_ongoing():
+      with treelog.context(f'timestep {timestep}'):
+
+        # read displacements from interface
+        if interface.is_read_data_available():
+            force = interface.read_block_vector_data(readdataID, dataIndices)
+
+        # save checkpoint
+        if interface.is_action_required(precice.action_write_iteration_checkpoint()):
+            checkpoint = timestep, arguments
+            interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
+
+        # advance variables
+        timestep += 1
+        dt = min(precice_dt, timestepsize)
+        arguments = dict(dt=dt, u0=arguments['u'], v0=arguments['v'], F=force)
+        arguments = solver.solve_linear('u:testu,v:testv', res, arguments=arguments, constrain=cons)
+
+        # write forces to interface
+        if interface.is_write_data_required(dt):
+            writedata = wall.eval(ns.u, **arguments)
+            interface.write_block_vector_data(writedataID, dataIndices, writedata)
+
+        # do the coupling
+        precice_dt = interface.advance(dt)
+
+        # read checkpoint if required
+        if interface.is_action_required(precice.action_read_iteration_checkpoint()):
+            timestep, arguments = checkpoint
+            interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
+
+    interface.finalize()
+
+
+if __name__ == '__main__':
+    cli.run(main)