diff --git a/perpendicular-flap/README.md b/perpendicular-flap/README.md index 2d92ce849..e1e20c8cb 100644 --- a/perpendicular-flap/README.md +++ b/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. diff --git a/perpendicular-flap/solid-nutils/clean.sh b/perpendicular-flap/solid-nutils/clean.sh new file mode 100755 index 000000000..6893c1ea5 --- /dev/null +++ b/perpendicular-flap/solid-nutils/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_nutils . diff --git a/perpendicular-flap/solid-nutils/run.sh b/perpendicular-flap/solid-nutils/run.sh new file mode 100755 index 000000000..2f48dbf3c --- /dev/null +++ b/perpendicular-flap/solid-nutils/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +python3 solid.py richoutput=no diff --git a/perpendicular-flap/solid-nutils/solid.py b/perpendicular-flap/solid-nutils/solid.py new file mode 100644 index 000000000..829cfc0de --- /dev/null +++ b/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)