WIP: aste Replay Mode tutorial on elastic flap case

FSI/flap_perp/ASTE/.gitignore

@@ -0,0 +1,5 @@
+preCICE-output/
+Solid-fenics*
+*.log
+vtk*
+Solid/FSI-S

FSI/flap_perp/ASTE/Allclean

@@ -0,0 +1,64 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+echo "Cleaning..."
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+# Participant 1: Fluid
+Participant1="Fluid"
+cd ${Participant1}
+    # Clean the case
+    cleanCase
+    rm -rfv 0
+    # Create an empty .foam file for ParaView
+    # Note: ".foam" triggers the native OpenFOAM reader of ParaView.
+    # Change to ".OpenFOAM" to use the OpenFOAM reader provided with OpenFOAM.
+    touch ${Participant1}.foam
+cd ..
+# Remove the log files
+rm -fv ${Participant1}_blockMesh.log
+rm -fv ${Participant1}_checkMesh.log
+rm -fv ${Participant1}_potentialFoam.log
+rm -fv ${Participant1}_decomposePar.log
+rm -fv ${Participant1}.log
+rm -fv ${Participant1}_reconstructPar.log
+
+# Participant 2: Solid
+Participant2="Solid"
+cd ${Participant2}
+cd FSI-S
+    # Clean the case
+	rm -fv *.vtu
+	rm -fv *.pvd
+    rm -fv *.log
+cd ..
+cd ..
+# Remove the log files
+rm -fv spooles.out
+rm -fv ${Participant2}.log
+
+# Remove the precice-events-files
+rm -fv precice-fenics-events.json
+rm -fv precice-Fluid-events.json
+
+# Remove the preCICE-related log files
+echo "Deleting the preCICE log files..."
+rm -fv \
+    precice-*.log \
+    precice-postProcessingInfo.log
+
+# Output files for preCICE versions before 1.2:
+rm -fv \
+    iterations-${Participant1}.txt iterations-${Participant2}.txt \
+    convergence-${Participant1}.txt convergence-${Participant2}.txt \
+    Events-${Participant1}.log Events-${Participant2}.log \
+    EventTimings-${Participant1}.log EventTimings-${Participant2}.log
+
+# Remove the preCICE address files
+rm -rfv precice-run
+rm -fv .*.address
+
+echo "Cleaning complete!"
+#------------------------------------------------------------------------------

FSI/flap_perp/ASTE/README.md

@@ -0,0 +1,89 @@
+## Overview
+
+This tutorial is an example how the Artificial Solver Testing Environment (aste) can be used. The idea is that you run a coupled simulation with two regular solvers and save the coupling data in every timestep. Then, this data is converted to aste-format. Finally, aste replaces one of the solvers and the simulation can be "replayed" using only one solver and aste with the previously computed results.
+
+This tutorial uses the results from the [OpenFOAM-FEniCS perpendicular flap tutorial](https://github.com/precice/tutorials/tree/master/FSI/flap_perp/OpenFOAM-FEniCS) as a basis. aste then replaces OpenFOAM.
+
+## Requirements
+
+To run this tutorial you need to install the following components:
+- [preCICE](https://github.com/precice/precice/wiki/Get-preCICE)
+- [aste](https://github.com/precice/aste/tree/develop)
+- [FEniCS](https://fenicsproject.org/)
+- [FEniCS-Adapter](https://github.com/precice/fenics-adapter)
+- OpenFOAM, e.g. [OpenFOAM 7](https://openfoam.org/version/7/)
+- [OpenFOAM Adapter](https://github.com/precice/openfoam-adapter/wiki/Building) matching the OpenFOAM version.
+
+Make sure to add `aste/build` to the `PATH` such that the python scripts and `preciceMap` can be found from anywhere on your system.
+
+## Step-by-Step explanations
+
+### Generating vtk output during a simulation
+
+The base case for this tutorial is the OpenFOAM-FEniCS perpendicular flap tutorial. So let's start in the the root-directory of that simulation [`tutorials/FSI/flap_perp/OpenFOAM-FEniCS`](https://github.com/precice/tutorials/tree/master/FSI/flap_perp/OpenFOAM-FEniCS).
+
+To generate vtk output in preCICE, add the following line to the `precice-config.xml` in `tutorials/FSI/flap_perp/OpenFOAM-FEniCS` after [line 41](https://github.com/precice/tutorials/blob/develop/FSI/flap_perp/OpenFOAM-FEniCS/precice-config.xml#L41) in the Solid participant:
+```
+<export:vtk directory="preCICE-output" />
+```
+Then, run the simulation as explained in the [`README.md`](https://github.com/precice/tutorials/blob/develop/FSI/flap_perp/OpenFOAM-FEniCS/README.md) of the case.
+
+The exports can be found in the `preCICE-output` directory.
+
+### Converting the output to aste format.
+
+Copy the `preCICE-output` folder to the root directory of the aste tutorial `tutorials/FSI/flap_perp/ASTE`.
+To convert the files to the correct format, open the `preCICE-output` folder and run
+
+`precice_to_aste.py Solid-fenics -n 500 -t Forces0 --datadim 3`
+
+### Replay of the simulation with aste
+
+#### The quick way to run 
+
+1. Copy `perp_flap.py` from `tutorials/FSI/flap_perp/OpenFOAM-FEniCS/Solid` into `tutorials/FSI/flap_perp/ASTE/Solid`.
+
+2. Run `python3 Solid/perp-flap.py` and `preciceMap -v -c precice-config.xml -p A --mesh preCICE-output/Solid-fenics --vectordata` in two terminals.
+
+Read on if you want to know what to change in the configuration files starting from the OpenFOAM-FEniCS tutorial.
+
+#### The way to do it yourself
+
+The aste-FEniCS tutorial provides its own `precice-config.xml` and `precice-adapter-config-fsi-s.json` for getting started quickly.
+
+However, these can also be generated by modifiying the files from the OpenFOAM-FEniCS case. The main purpose of this tutorial is to explain how the replay mode with aste can be used for arbirtrary simulation setups.
+
+1. Copy the files from the OpenFOAM-FEniCS tutorial into the respective directory of the aste-tutorial. 
+2. Remove the line `"write_data_name": "Displacements0" ` from `precice-adapter-config-fsi-s.json`, because we don't write data to aste.
+3. Change the `read_data_name` from `Forces0` to `Data` since aste works with the general data name `Data`.
+4. In ```precice-config.xml``` we require some more changes:
+    1. Rename `Forces0` to `Data` throughout the whole file.
+    2. Delete all lines where `Displacements0` occurs, since we are only coupling in one direction (from aste to FEniCS).
+5. Then we process the Fluid meshes to fake the Fluid participant using aste:
+    1. Rename `Fluid-Mesh-Faces` to `MeshA` throughout the whole file.
+    2. Delete the `mesh` `"Fluid-Mesh-Nodes"` and all lines where `"Fluid-Mesh-Nodes"` occurs.
+    3. Rename the participant `Fluid` to `A`
+6. Change the coupling scheme to an explicit scheme:
+    1. Change `serial-implicit` to `serial-explicit`.
+    2. Explicit schemes don't iterate such that they don't need convergence measures, extrapolation or acceleration. Delete everything in the coulping scheme except for 
+    ```
+      <participants first="A" second="fenics" />
+      <max-time value="5.0" />
+      <time-window-size value="1.e-2" />
+      <exchange data="Data" mesh="Solid" from="A" to="fenics" />
+    ```
+
+#### Run
+
+After applying these preparations, we are now able to run the tutorial in two terminals with `python3 Solid/perp-flap.py` and `preciceMap -v -c precice-config.xml -p A --mesh preCICE-output/Solid-fenics --vectordata`.
+
+### Play Around
+
+If you want to explore more possibilities of the Replay-Mode here are some ideas:
+
+- Export the output at the Fluid participant. Instead of `preCICE-output/Solid-fenics.*` the files containing the forces will be `preCICE-output/Fluid-Mesh-Faces-Fluid.*`
+
+- Export the forces in the [Calculix-Openfoam-tutorial](https://github.com/precice/tutorials/tree/master/FSI/flap_perp/OpenFOAM-CalculiX) with the same geometry and feed them to FEniCS with aste. 
+
+- If you are familiar with the OpenFOAM-Adapter you can also substitute FEniCS with aste such that aste writes previously exported displacements to OpenFOAM.
+

FSI/flap_perp/ASTE/Solid/precice-adapter-config-fsi-s.json

@@ -0,0 +1,8 @@
+{
+  "solver_name": "fenics",
+  "config_file_name": "../precice-config.xml",
+  "interface": {
+      "coupling_mesh_name": "Solid",
+      "read_data_name": "Data"
+    }
+}

FSI/flap_perp/ASTE/precice-config.xml

@@ -0,0 +1,52 @@
+<?xml version="1.0"?>
+
+<precice-configuration>
+
+  <log enabled="true">
+    <sink type="stream" output="stdout" filter="%Severity% > debug" enabled="true" />
+    <sink type="file" output="filtered.log"  filter= "(%Severity% > debug) or (%Severity% >= trace and %Module% contains SolverInterfaceImpl)" enabled="true" />
+    <sink type="file" output="debug.log" filter="" enabled="true"/>
+  </log>
+
+  <solver-interface dimensions="3">
+
+    <!-- Data fields that are exchanged between the solvers -->
+    <data:vector name="Data" />
+
+    <!-- A common mesh that uses these data fields -->
+    <mesh name="MeshA" flip-normals="no">
+      <use-data name="Data" />
+    </mesh>
+
+    <mesh name="Solid" flip-normals="no">
+      <use-data name="Data" />
+    </mesh>
+
+    <m2n:sockets from="A" to="fenics"/> 
+
+    <participant name="A">
+      <use-mesh name="MeshA" provide="yes" />
+      <use-mesh name="Solid" provide="no" from="fenics" />
+      <write-data name="Data" mesh="MeshA"/>
+      <mapping:nearest-neighbor constraint="conservative" direction="write" from="MeshA" to="Solid" />
+      <export:vtk every-n-time-windows="1" directory="vtkA/" normals="0"/>
+    </participant>
+
+    <participant name="fenics">
+      <use-mesh name="Solid" provide="yes" />
+      <read-data name="Data" mesh="Solid" />
+
+      <export:vtk every-n-time-windows="1" directory="vtkB/" normals="0"/>
+    </participant>
+
+
+    <coupling-scheme:serial-explicit>
+      <participants first="A" second="fenics" />
+      <max-time value="5.0" />
+      <time-window-size value="1.e-2" />
+      <exchange data="Data" mesh="Solid" from="A" to="fenics" />
+    </coupling-scheme:serial-explicit>
+
+  </solver-interface>
+
+</precice-configuration>

FSI/flap_perp/OpenFOAM-FEniCS/Solid/perp-flap.py

@@ -1,18 +1,16 @@
 #Import required libs
-from fenics import *
 from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \
-TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad,\
-Identity, inner,dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource
+    TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, Identity, inner,dx, ds, sym, grad, lhs, \
+    rhs, dot, File, solve, PointSource, XDMFFile, TensorFunctionSpace, LocalSolver, assemble_system
 import dolfin
-
 from ufl import nabla_div
 import numpy as np
 import matplotlib.pyplot as plt
 from fenicsadapter import Adapter
 from enum import Enum
 
 
-#define the two kinds of boundary: clamped and coupling Neumann Boundary
+# define the two kinds of boundary: clamped and coupling Neumann Boundary
 def clamped_boundary(x, on_boundary):
     return on_boundary and abs(x[1])<tol
 
@@ -69,11 +67,10 @@ def Neumann_Boundary(x, on_boundary):
 
 coupling_boundary = AutoSubDomain(Neumann_Boundary)
 
-# create subdomain that resembles the 
 
-## get the adapter ready
+# get the adapter ready
 
-#read fenics-adapter json-config-file)
+# read fenics-adapter json-config-file)
 
 adapter_config_filename = "precice-adapter-config-fsi-s.json"
 
@@ -95,10 +92,10 @@ def Neumann_Boundary(x, on_boundary):
 # clamp the beam at the bottom
 bc = DirichletBC(V, Constant((0,0)), clamped_boundary)
 
-#alpha method parameters
+# alpha method parameters
 
-alpha_m = Constant(0.2)
-alpha_f = Constant(0.4)
+alpha_m = Constant(0.0)
+alpha_f = Constant(0.0)
 gamma   = Constant(0.5+alpha_f-alpha_m)
 beta    = Constant((gamma+0.5)**2/4.)
 
@@ -180,7 +177,7 @@ def avg(x_old, x_new, alpha):
 
 # Prepare for time-stepping
 
-#parameters for Time-Stepping
+# parameters for Time-Stepping
 t=0.0
 n=0
 time = []
@@ -196,8 +193,32 @@ def avg(x_old, x_new, alpha):
 u_np1.rename("Displacement", "")
 displacement_out << u_n
 
+# stress computation
+def local_project(v, V, u=None):
+    """Element-wise projection using LocalSolver"""
+    dv = TrialFunction(V)
+    v_ = TestFunction(V)
+    a_proj = inner(dv, v_)*dx
+    b_proj = inner(v, v_)*dx
+    solver = LocalSolver(a_proj, b_proj)
+    solver.factorize()
+    if u is None:
+        u = Function(V)
+        solver.solve_local_rhs(u)
+        return u
+    else:
+        solver.solve_local_rhs(u)
+        return
+
+Vsig = TensorFunctionSpace(mesh, "CG", 1)
+sig = Function(Vsig, name="sigma")
+
+xdmf_file = XDMFFile("elastodynamics-results.xdmf")
+xdmf_file.parameters["flush_output"] = True
+xdmf_file.parameters["functions_share_mesh"] = True
+xdmf_file.parameters["rewrite_function_mesh"] = False
 
-#time loop for coupling
+# time loop for coupling
 
 
 while precice.is_coupling_ongoing():
@@ -221,18 +242,21 @@ def avg(x_old, x_new, alpha):
         update_fields(u_np1, saved_u_old, v_n, a_n)
 
         if n % 10==0:
+            local_project(sigma(u_n), Vsig, sig)
+
             displacement_out << (u_n,t)
-
+            xdmf_file.write(u_n,t)
+            xdmf_file.write(sig,t)
         u_tip.append(u_n(0.,1.)[0])
         time.append(t)
 
-# Plot tip displacement evolution
-
+
+precice.finalize()
+
+# Plot tip displacement evolution    
 displacement_out << u_n 
 plt.figure()
 plt.plot(time, u_tip)
 plt.xlabel("Time")
 plt.ylabel("Tip displacement")
 plt.show()
-
-