add nutils solid participant

Author: gertjanvanzwieten

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,96 @@
+#! /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):
+
+    topo, geom = mesh.rectilinear([numpy.linspace(-.05, .05, nelems+1), numpy.linspace(0, 1, 10*nelems+1)])
+    bezier = topo.sample('bezier', 2)
+    wall = topo.boundary['left,top,right']
+    wallsample = wall.sample('gauss', degree=2)
+
+    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('tdS', wallsample.basis(), shape=(2,))
+    ns.dudt_i = '(u_i - u0_i) / dt'
+    ns.dvdt_i = '(v_i - v0_i) / dt'
+    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))'
+    ns.ρ = density
+    ns.x_i = 'X_i + u_i'
+
+    # continuum equations: ρ v' = ∇·σ + F, u' = v
+    res = topo.integral('(testv_i (dudt_i - v_i) + testu_i ρ dvdt_i + ∇_j(testu_i) σ_ij) dV' @ ns, degree=2)
+    res -= wallsample.integral('testu_i tdS_i' @ 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 still
+    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, wallsample.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
+    traction = numpy.zeros((wallsample.npoints, 2))
+
+    while interface.is_coupling_ongoing():
+      with treelog.context(f'timestep {timestep}'):
+
+        # read displacements from interface
+        if interface.is_read_data_available():
+            traction = 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'], tdS=traction)
+        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 = wallsample.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())
+
+        if interface.is_time_window_complete():
+            treelog.user('tip displacement:', topo.boundary['top'].sample('uniform', 1).eval(ns.u, **arguments).mean(0))
+
+    interface.finalize()
+
+
+if __name__ == '__main__':
+    cli.run(main)