diff --git a/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict b/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict
index a9e847699..85cf083c2 100644
--- a/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict
+++ b/turek-hron-fsi3/fluid-openfoam/system/decomposeParDict
@@ -6,12 +6,12 @@ FoamFile
object decomposeParDict;
}
-numberOfSubdomains 25;
+numberOfSubdomains 3;
method hierarchical;
hierarchicalCoeffs
{
- n (5 5 1);
+ n (3 1 1);
delta 0.001;
order xyz;
}
diff --git a/turek-hron-fsi3/fluid-openfoam/system/preciceDict b/turek-hron-fsi3/fluid-openfoam/system/preciceDict
index 872566d38..2a2b63d9a 100644
--- a/turek-hron-fsi3/fluid-openfoam/system/preciceDict
+++ b/turek-hron-fsi3/fluid-openfoam/system/preciceDict
@@ -26,10 +26,9 @@ interfaces
writeData
(
- Stress
+ Force
);
};
-
Interface2
{
mesh Fluid-Mesh-Nodes;
diff --git a/turek-hron-fsi3/plot-displacement.sh b/turek-hron-fsi3/plot-displacement.sh
index 141f54ec5..b44e3b759 100755
--- a/turek-hron-fsi3/plot-displacement.sh
+++ b/turek-hron-fsi3/plot-displacement.sh
@@ -5,5 +5,5 @@ gnuplot -p << EOF
set xlabel 'time [s]'
set ylabel 'y-displacement [m]'
set linestyle 1 lt 2 lc 1 # red-dashed
- plot "solid-dealii/precice-Solid-watchpoint-Flap-Tip.log" using 1:5 with lines notitle
+ plot "solid-fenics/precice-Solid-watchpoint-Flap-Tip.log" using 1:5 with lines notitle
EOF
diff --git a/turek-hron-fsi3/precice-config.xml b/turek-hron-fsi3/precice-config.xml
index ae004917b..26f3f77a5 100644
--- a/turek-hron-fsi3/precice-config.xml
+++ b/turek-hron-fsi3/precice-config.xml
@@ -7,11 +7,11 @@
enabled="true" />
-
+
-
+
@@ -20,45 +20,44 @@
-
+
-
+
-
+
+
+
+
-
-
+
-
-
-
-
+
-
+
-
+
-
+
diff --git a/turek-hron-fsi3/solid-fenics/clean.sh b/turek-hron-fsi3/solid-fenics/clean.sh
new file mode 100755
index 000000000..d3f3318b4
--- /dev/null
+++ b/turek-hron-fsi3/solid-fenics/clean.sh
@@ -0,0 +1,6 @@
+#!/usr/bin/env sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+clean_fenics .
diff --git a/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json b/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json
new file mode 100644
index 000000000..91a042b17
--- /dev/null
+++ b/turek-hron-fsi3/solid-fenics/precice-adapter-config-fsi-s.json
@@ -0,0 +1,9 @@
+{
+ "participant_name": "Solid",
+ "config_file_name": "../precice-config.xml",
+ "interface": {
+ "coupling_mesh_name": "Solid-Mesh",
+ "write_data_name": "Displacement",
+ "read_data_name": "Force"
+ }
+}
diff --git a/turek-hron-fsi3/solid-fenics/requirements.txt b/turek-hron-fsi3/solid-fenics/requirements.txt
new file mode 100644
index 000000000..e7a847a98
--- /dev/null
+++ b/turek-hron-fsi3/solid-fenics/requirements.txt
@@ -0,0 +1,11 @@
+fenicsprecice==2.2.0rc1
+numpy >1, <2
+matplotlib
+
+# Assuming FEniCS from ppa:fenics-packages/fenics was installed https://fenicsproject.org/download/archive/
+# Use --system-site-packages in venv
+fenics-dijitso==2019.2.0.dev0
+fenics-dolfin==2019.2.0.13.dev0
+fenics-ffc==2019.2.0.dev0
+fenics-fiat==2019.2.0.dev0
+fenics-ufl-legacy==2022.3.0
\ No newline at end of file
diff --git a/turek-hron-fsi3/solid-fenics/run.sh b/turek-hron-fsi3/solid-fenics/run.sh
new file mode 100755
index 000000000..867c1db13
--- /dev/null
+++ b/turek-hron-fsi3/solid-fenics/run.sh
@@ -0,0 +1,13 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 -m venv --system-site-packages .venv
+. .venv/bin/activate
+pip install -r requirements.txt
+
+python3 solid.py
+
+close_log
diff --git a/turek-hron-fsi3/solid-fenics/solid.py b/turek-hron-fsi3/solid-fenics/solid.py
new file mode 100644
index 000000000..eeae08ceb
--- /dev/null
+++ b/turek-hron-fsi3/solid-fenics/solid.py
@@ -0,0 +1,252 @@
+# Import required libs
+from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \
+ TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \
+ Identity, inner, dx, sym, grad, lhs, rhs, File, solve, assemble_system, div
+
+import numpy as np
+from fenicsprecice import Adapter
+from enum import Enum
+
+# Beam geometry
+dim = 2 # number of dimensions
+L = 0.35 # length
+H = 0.02 # height
+
+y_bottom = 0.2 - 0.5 * H # y coordinate of bottom surface of beam
+y_top = y_bottom + H # y coordinate of top surface of beam
+x_left = 0.25 # x coordinate of left surface of beam
+x_right = x_left + L # x coordinate of right surface of beam
+
+
+# define the two inside functions to determine the boundaries: clamped Dirichlet and coupling Neumann Boundary
+def left_boundary(x, on_boundary):
+ """
+ inside-function for the clamped Dirichlet Boundary.
+
+ Apply Dirichlet boundary on left part of the beam.
+ """
+ return on_boundary and abs(x[0] - x_left) < tol # left boundary of beam
+
+
+def remaining_boundary(x, on_boundary):
+ """
+ inside-function for the coupling Neumann Boundary.
+
+ Apply Neumann boundary on top, bottom and right (=remaining) part of the beam.
+ """
+ return on_boundary and (abs(x[1] - y_top) < tol or # top boundary of beam
+ abs(x[1] - y_bottom) < tol or # bottom boundary of beam
+ abs(x[0] - x_right) < tol) # right boundary of beam
+
+
+# Numerical properties
+tol = 1E-14
+
+# Beam material properties
+rho = 1000 # density
+E = 5600000.0 # Young's modulus
+nu = 0.4 # Poisson's ratio
+lambda_ = Constant(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) # first Lame constant
+mu = Constant(E / (2.0 * (1.0 + nu))) # second Lame constant
+
+# create Mesh
+n_x_Direction = 20 # DoFs in x direction
+n_y_Direction = 4 # DoFs in y direction
+mesh = RectangleMesh(Point(x_left, y_bottom), Point(x_right, y_top), n_x_Direction, n_y_Direction)
+
+# create Function Space
+V = VectorFunctionSpace(mesh, 'P', 2)
+
+# Trial and Test Functions
+du = TrialFunction(V)
+v = TestFunction(V)
+
+# displacement fields
+u_np1 = Function(V)
+saved_u_old = Function(V)
+
+# function known from previous timestep
+u_n = Function(V)
+v_n = Function(V)
+a_n = Function(V)
+
+# initial value for force and displacement field
+f_N_function = interpolate(Expression(("1", "0"), degree=1), V)
+u_function = interpolate(Expression(("0", "0"), degree=1), V)
+
+# define coupling boundary
+coupling_boundary = AutoSubDomain(remaining_boundary)
+
+# create Adapter
+precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")
+
+# create subdomains used by the adapter
+clamped_boundary_domain = AutoSubDomain(left_boundary)
+force_boundary = AutoSubDomain(remaining_boundary)
+
+# Initialize the coupling interface
+# Function space V is passed twice as both read and write functions are defined using the same space
+precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=clamped_boundary_domain)
+
+precice_dt = precice.get_max_time_step_size()
+fenics_dt = precice_dt
+dt = Constant(np.min([precice_dt, fenics_dt]))
+
+# alpha method parameters
+alpha_m = Constant(0.2)
+alpha_f = Constant(0.4)
+# alpha_m = Constant(0)
+# alpha_f = Constant(0)
+
+"""
+Check requirements for alpha_m and alpha_f from
+ Chung, J., and Hulbert, G. M. (June 1, 1993). "A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation:
+ The Generalized-α Method." ASME. J. Appl. Mech. June 1993; 60(2): 371–375. https://doi.org/10.1115/1.2900803
+"""
+assert (float(alpha_m) <= float(alpha_f))
+assert (float(alpha_f) <= 0.5)
+
+gamma = Constant(0.5 + alpha_f - alpha_m)
+beta = Constant((gamma + 0.5) ** 2 * 0.25)
+
+# clamp (u == 0) the beam at the left
+bc = DirichletBC(V, Constant((0, 0)), left_boundary)
+
+
+# Define strain
+def epsilon(u):
+ return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
+
+
+# Define Stress tensor
+def sigma(u):
+ return lambda_ * div(u) * Identity(dim) + 2 * mu * epsilon(u)
+
+
+# Define Mass form
+def m(u, v):
+ return rho * inner(u, v) * dx
+
+
+# Elastic stiffness form
+def k(u, v):
+ return inner(sigma(u), sym(grad(v))) * dx
+
+
+def update_a(u, u_old, v_old, a_old, ufl=True):
+ if ufl:
+ dt_ = dt
+ beta_ = beta
+ else:
+ dt_ = float(dt)
+ beta_ = float(beta)
+
+ return (u - u_old - dt_ * v_old) / beta / dt_ ** 2 - .5 * (1 - 2 * beta_) / beta_ * a_old
+
+
+def update_v(a, u_old, v_old, a_old, ufl=True):
+ if ufl:
+ dt_ = dt
+ gamma_ = gamma
+ else:
+ dt_ = float(dt)
+ gamma_ = float(gamma)
+
+ return v_old + dt_ * ((1 - gamma_) * a_old + gamma_ * a)
+
+
+def update_fields(u, u_old, v_old, a_old):
+ """Update all fields at the end of a timestep."""
+ u_vec, u0_vec = u.vector(), u_old.vector()
+ v0_vec, a0_vec = v_old.vector(), a_old.vector()
+
+ # call update functions
+ a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
+ v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)
+
+ # assign u->u_old
+ v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
+ u_old.vector()[:] = u.vector()
+
+
+def avg(x_old, x_new, alpha):
+ return alpha * x_old + (1 - alpha) * x_new
+
+
+# residual
+a_np1 = update_a(du, u_n, v_n, a_n, ufl=True)
+v_np1 = update_v(a_np1, u_n, v_n, a_n, ufl=True)
+
+res = m(avg(a_n, a_np1, alpha_m), v) + k(avg(u_n, du, alpha_f), v)
+
+a_form = lhs(res)
+L_form = rhs(res)
+
+# Prepare for time-stepping
+t = 0.0
+n = 0
+E_ext = 0
+
+displacement_out = File("output/u_fsi.pvd")
+
+u_n.rename("Displacement", "")
+u_np1.rename("Displacement", "")
+displacement_out << u_n
+
+# time loop for coupling
+while precice.is_coupling_ongoing():
+
+ if precice.requires_writing_checkpoint(): # write checkpoint
+ precice.store_checkpoint((u_n, v_n, a_n), t, n)
+
+ precice_dt = precice.get_max_time_step_size()
+ dt = Constant(np.min([precice_dt, fenics_dt]))
+
+ # read data from preCICE and get a new coupling expression
+ # sample force F at $F(t_{n+1-\alpha_f})$ (see generalized alpha paper)
+ read_data = precice.read_data((1 - float(alpha_f)) * dt)
+ # read_data = precice.read_data(dt)
+
+ # Update the point sources on the coupling boundary with the new read data
+ Forces_x, Forces_y = precice.get_point_sources(read_data)
+
+ A, b = assemble_system(a_form, L_form, bc)
+
+ b_forces = b.copy() # b is the same for every iteration, only forces change
+
+ for ps in Forces_x:
+ ps.apply(b_forces)
+ for ps in Forces_y:
+ ps.apply(b_forces)
+
+ assert (b is not b_forces)
+ solve(A, u_np1.vector(), b_forces)
+
+ # Write new displacements to preCICE
+ precice.write_data(u_np1)
+
+ # Call to advance coupling, also returns the optimum time step value
+ precice.advance(float(dt))
+
+ # Either revert to old step if timestep has not converged or move to next timestep
+ if precice.requires_reading_checkpoint(): # roll back to checkpoint
+ uva_cp, t_cp, n_cp = precice.retrieve_checkpoint()
+ u_cp, v_cp, a_cp = uva_cp
+ u_n.assign(u_cp)
+ v_n.assign(v_cp)
+ a_n.assign(a_cp)
+ t = t_cp
+ n = n_cp
+ else:
+ update_fields(u_np1, u_n, v_n, a_n)
+ u_n.assign(u_np1)
+ t += float(dt)
+ n += 1
+
+ if precice.is_time_window_complete():
+ if n % 20 == 0:
+ displacement_out << (u_n, t)
+
+displacement_out << u_n
+
+precice.finalize()