From f025c945ebdb15c0dae8e1a222030ded428cfe5d Mon Sep 17 00:00:00 2001
From: Niklas Vinnitchenko <niklas.vinnitchenko@tum.de>
Date: Mon, 17 Jun 2024 20:59:00 +0200
Subject: [PATCH 1/7] Added solid-fake for perpendicular flap tutorial

---
 perpendicular-flap/solid-fake/fake.py | 148 ++++++++++++++++++++++++++
 1 file changed, 148 insertions(+)
 create mode 100644 perpendicular-flap/solid-fake/fake.py

diff --git a/perpendicular-flap/solid-fake/fake.py b/perpendicular-flap/solid-fake/fake.py
new file mode 100644
index 000000000..b4716b50f
--- /dev/null
+++ b/perpendicular-flap/solid-fake/fake.py
@@ -0,0 +1,148 @@
+from __future__ import division
+
+import numpy as np
+import precice
+def displacement_param(t):
+    # this parameter exists to define a time-dependent displacement of the flap's tip
+    if t > 1:
+        return 0.5
+    else:
+        return 0.5*t
+
+
+# computes the flap displacements
+def compute_flap_displacements(t, height_flap, width_flap, vertices_per_side, dimensions=2):
+    a = displacement_param(t)
+    # ASSUMPTION: the middle of the tip must still have a length of 1
+    # For that, one can use the formula height_flap=int_0^y(sqrt(1+f'(z)^2) dz
+    #
+    # In the following, the function 1/a cosh(ay) - 1/a is used because the analytical solution can be determined for
+    # the height_flap
+    #
+    # We get for the y-coordinate of the tip's middle: y=arcsinh(a*height_flap)/a
+    # the y-displacement is therefore: y_displ =height_flap - y
+    # the x-coordinates can hence be computed by x(y)= 1/a cosh(ay) - 1/a in [0, y]
+
+    # for a = 0, we just use a straight line
+    if a == 0:
+        # if a=0 the flap is straight -> displacements are all 0
+        return np.zeros((2 * vertices_per_side, dimensions))
+
+    fun_x_coord = lambda y: (np.cosh(a * y) - 1) / a
+
+    ##########################
+    # IN THE FOLLOWING: a!=0 #
+    ##########################
+
+    #uniformly distributed vertices
+    h_vert = height_flap / (vertices_per_side - 1)
+    # with the function above, we know, that the middle flap's length can be computed with len=sinh(a*x)/a
+    # now, we will compute the points of the flap_shape, were the len is not the height_flap, but n*h_vert
+    # the reason: we guarantee that we compute the points of the borders uniformly distributed on the basis that
+    # the flap is straight => n*h_vert = sinh(a*x)/a <=> x = arcsinh(n*h_vert*a)/a
+    displacements_right = np.zeros((vertices_per_side, dimensions))
+    displacements_left = np.zeros((vertices_per_side, dimensions))
+
+    # get uniformly distanced (x,y) points of the flap
+    y_values_middle = np.array([np.arcsinh(i * h_vert * a) / a for i in range(vertices_per_side)])
+    x_values_middle = np.array([fun_x_coord(y_values_middle[i]) for i in range(len(y_values_middle))])
+    # dummy values for [0] entries to avoid warnings
+    x_values_middle[0] = y_values_middle[0] = 1
+
+    # now we compute the (x,y) values of the border according to our (x,y) values of the flap's middle
+
+    # step 1: compute the tangent with slope normal to the flap's middle shape
+    # derivative of fun_x_coord is sinh(ay)
+    slopes = -1 / np.sinh(a * y_values_middle)
+
+    # step 2: compute the x and y coordinates that are on the tangent of step 1 and distanced width_flap/2 from
+    # (x_values_middle,y_values_middle)
+    # use slope = dx/dy and (width/2)^2 = dy^2 + dx^2
+    # -> dy^2 = (width/2)^2/(slope^2+1)
+    # we get two results from that because the result is the absolute value!
+    dy = np.divide(width_flap / 2, np.sqrt(slopes * slopes + 1))
+    # -> dx = slope*dy (second dx only negative value of result)
+    dx = slopes * dy
+
+    # step 3: now, we can compute with dx, dy, x_values_middle and y_values_middle the coordinates
+    # of the border's vertices and eventually the displacements!
+    displacements_right[:, 0] = x_values_middle - dx
+    displacements_right[:, 1] = y_values_middle - dy
+    displacements_left[:, 0] = x_values_middle + dx  # second dx has inverted sign
+    displacements_left[:, 1] = y_values_middle + dy
+
+    # displacements: remember that y-coordinates of the vertices are normally uniformly distributed,
+    # i.e. y_0=0, y_1=h_vert, y_2=2*h_vert,...
+    # and the x-coordinates are for the left border -width/2 and for the right border width/2
+    default_x_left = -width_flap / 2
+    default_x_right = width_flap / 2
+    for i in range(vertices_per_side):
+        #displacements right
+        displacements_right[i][0] = displacements_right[i][0] - default_x_right
+        displacements_right[i][1] = displacements_right[i][1] - h_vert * i
+        #displacements left
+        displacements_left[i][0] = displacements_left[i][0] - default_x_left
+        displacements_left[i][1] = displacements_left[i][1] - h_vert * i
+    # bottom is fixed -> x & y displacement 0
+    displacements_left[0][:] = 0
+    displacements_right[0][:] = 0
+
+    return np.concatenate((displacements_left, displacements_right), axis=0)
+
+configuration_file_name = "../precice-config.xml"
+participant_name = "Solid"
+mesh_name = "Solid-Mesh"
+write_data_name = 'Displacement'
+
+solver_process_index = 0
+solver_process_size = 1
+
+# define mesh
+H = 1
+W = 0.1
+
+interface = precice.Participant(participant_name, configuration_file_name, solver_process_index, solver_process_size)
+dimensions = interface.get_mesh_dimensions(mesh_name)
+assert (dimensions == 2)
+
+x_left = 0.0 - 0.5 * W  # left boundary of the flap
+x_right = 0.5 * W  # right boundary of the flap
+y_bottom = 0.0  # bottom of the flap
+y_top = y_bottom + H  # top of the flap
+
+n = 24  # Number of vertices per side
+t = 0
+
+vertices = np.zeros((2 * n, dimensions))
+vertices[:n, 0] = x_left  # vertices at left side of flap
+vertices[n:, 0] = x_right  # vertices at right side of flap
+vertices[:n, 1] = np.linspace(y_bottom, y_top, n)  # equally disrtibuted vertices left
+vertices[n:, 1] = np.linspace(y_bottom, y_top, n)  # equally disrtibuted vertices right
+vertex_ids = interface.set_mesh_vertices(mesh_name, vertices)
+
+interface.initialize()
+# change if necessary
+solver_dt = np.inf
+# for checkpointing
+t_cp = 0
+
+while interface.is_coupling_ongoing():
+    if interface.requires_writing_checkpoint():
+        t_cp = t
+
+    precice_dt = interface.get_max_time_step_size()
+    dt = min([solver_dt, precice_dt])
+    # wiggle the flap
+    write_data = compute_flap_displacements(t, H, W, n)
+
+    interface.write_data(mesh_name, write_data_name, vertex_ids, write_data)
+    interface.advance(dt)
+
+    if interface.requires_reading_checkpoint():
+        t = t_cp
+    else:
+        # update t
+        t += dt
+
+interface.finalize()
+

From a94a44a0eac9526619c906ae861640acfcf58dc2 Mon Sep 17 00:00:00 2001
From: Niklas Vinnitchenko <niklas.vinnitchenko@tum.de>
Date: Mon, 17 Jun 2024 22:05:19 +0200
Subject: [PATCH 2/7] Formatting

---
 perpendicular-flap/solid-fake/fake.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/perpendicular-flap/solid-fake/fake.py b/perpendicular-flap/solid-fake/fake.py
index b4716b50f..3e82077d2 100644
--- a/perpendicular-flap/solid-fake/fake.py
+++ b/perpendicular-flap/solid-fake/fake.py
@@ -2,12 +2,14 @@
 
 import numpy as np
 import precice
+
+
 def displacement_param(t):
     # this parameter exists to define a time-dependent displacement of the flap's tip
     if t > 1:
         return 0.5
     else:
-        return 0.5*t
+        return 0.5 * t
 
 
 # computes the flap displacements
@@ -28,13 +30,13 @@ def compute_flap_displacements(t, height_flap, width_flap, vertices_per_side, di
         # if a=0 the flap is straight -> displacements are all 0
         return np.zeros((2 * vertices_per_side, dimensions))
 
-    fun_x_coord = lambda y: (np.cosh(a * y) - 1) / a
+    def fun_x_coord(y): return (np.cosh(a * y) - 1) / a
 
     ##########################
     # IN THE FOLLOWING: a!=0 #
     ##########################
 
-    #uniformly distributed vertices
+    # uniformly distributed vertices
     h_vert = height_flap / (vertices_per_side - 1)
     # with the function above, we know, that the middle flap's length can be computed with len=sinh(a*x)/a
     # now, we will compute the points of the flap_shape, were the len is not the height_flap, but n*h_vert
@@ -77,10 +79,10 @@ def compute_flap_displacements(t, height_flap, width_flap, vertices_per_side, di
     default_x_left = -width_flap / 2
     default_x_right = width_flap / 2
     for i in range(vertices_per_side):
-        #displacements right
+        # displacements right
         displacements_right[i][0] = displacements_right[i][0] - default_x_right
         displacements_right[i][1] = displacements_right[i][1] - h_vert * i
-        #displacements left
+        # displacements left
         displacements_left[i][0] = displacements_left[i][0] - default_x_left
         displacements_left[i][1] = displacements_left[i][1] - h_vert * i
     # bottom is fixed -> x & y displacement 0
@@ -89,6 +91,7 @@ def compute_flap_displacements(t, height_flap, width_flap, vertices_per_side, di
 
     return np.concatenate((displacements_left, displacements_right), axis=0)
 
+
 configuration_file_name = "../precice-config.xml"
 participant_name = "Solid"
 mesh_name = "Solid-Mesh"
@@ -145,4 +148,3 @@ def compute_flap_displacements(t, height_flap, width_flap, vertices_per_side, di
         t += dt
 
 interface.finalize()
-

From 806d99fc2e1500cb8b60bd1e2fd0617bbb21525c Mon Sep 17 00:00:00 2001
From: Niklas Vinnitchenko <niklas.vinnitchenko@tum.de>
Date: Sun, 23 Jun 2024 13:59:45 +0200
Subject: [PATCH 3/7] Simplify fake model

---
 perpendicular-flap/solid-fake/fake.py | 113 +++++---------------------
 1 file changed, 20 insertions(+), 93 deletions(-)

diff --git a/perpendicular-flap/solid-fake/fake.py b/perpendicular-flap/solid-fake/fake.py
index 3e82077d2..4f8f161df 100644
--- a/perpendicular-flap/solid-fake/fake.py
+++ b/perpendicular-flap/solid-fake/fake.py
@@ -4,92 +4,20 @@
 import precice
 
 
-def displacement_param(t):
-    # this parameter exists to define a time-dependent displacement of the flap's tip
-    if t > 1:
-        return 0.5
-    else:
-        return 0.5 * t
-
-
-# computes the flap displacements
-def compute_flap_displacements(t, height_flap, width_flap, vertices_per_side, dimensions=2):
-    a = displacement_param(t)
-    # ASSUMPTION: the middle of the tip must still have a length of 1
-    # For that, one can use the formula height_flap=int_0^y(sqrt(1+f'(z)^2) dz
-    #
-    # In the following, the function 1/a cosh(ay) - 1/a is used because the analytical solution can be determined for
-    # the height_flap
-    #
-    # We get for the y-coordinate of the tip's middle: y=arcsinh(a*height_flap)/a
-    # the y-displacement is therefore: y_displ =height_flap - y
-    # the x-coordinates can hence be computed by x(y)= 1/a cosh(ay) - 1/a in [0, y]
-
-    # for a = 0, we just use a straight line
-    if a == 0:
-        # if a=0 the flap is straight -> displacements are all 0
-        return np.zeros((2 * vertices_per_side, dimensions))
-
-    def fun_x_coord(y): return (np.cosh(a * y) - 1) / a
-
-    ##########################
-    # IN THE FOLLOWING: a!=0 #
-    ##########################
-
-    # uniformly distributed vertices
-    h_vert = height_flap / (vertices_per_side - 1)
-    # with the function above, we know, that the middle flap's length can be computed with len=sinh(a*x)/a
-    # now, we will compute the points of the flap_shape, were the len is not the height_flap, but n*h_vert
-    # the reason: we guarantee that we compute the points of the borders uniformly distributed on the basis that
-    # the flap is straight => n*h_vert = sinh(a*x)/a <=> x = arcsinh(n*h_vert*a)/a
-    displacements_right = np.zeros((vertices_per_side, dimensions))
-    displacements_left = np.zeros((vertices_per_side, dimensions))
-
-    # get uniformly distanced (x,y) points of the flap
-    y_values_middle = np.array([np.arcsinh(i * h_vert * a) / a for i in range(vertices_per_side)])
-    x_values_middle = np.array([fun_x_coord(y_values_middle[i]) for i in range(len(y_values_middle))])
-    # dummy values for [0] entries to avoid warnings
-    x_values_middle[0] = y_values_middle[0] = 1
-
-    # now we compute the (x,y) values of the border according to our (x,y) values of the flap's middle
-
-    # step 1: compute the tangent with slope normal to the flap's middle shape
-    # derivative of fun_x_coord is sinh(ay)
-    slopes = -1 / np.sinh(a * y_values_middle)
-
-    # step 2: compute the x and y coordinates that are on the tangent of step 1 and distanced width_flap/2 from
-    # (x_values_middle,y_values_middle)
-    # use slope = dx/dy and (width/2)^2 = dy^2 + dx^2
-    # -> dy^2 = (width/2)^2/(slope^2+1)
-    # we get two results from that because the result is the absolute value!
-    dy = np.divide(width_flap / 2, np.sqrt(slopes * slopes + 1))
-    # -> dx = slope*dy (second dx only negative value of result)
-    dx = slopes * dy
-
-    # step 3: now, we can compute with dx, dy, x_values_middle and y_values_middle the coordinates
-    # of the border's vertices and eventually the displacements!
-    displacements_right[:, 0] = x_values_middle - dx
-    displacements_right[:, 1] = y_values_middle - dy
-    displacements_left[:, 0] = x_values_middle + dx  # second dx has inverted sign
-    displacements_left[:, 1] = y_values_middle + dy
-
-    # displacements: remember that y-coordinates of the vertices are normally uniformly distributed,
-    # i.e. y_0=0, y_1=h_vert, y_2=2*h_vert,...
-    # and the x-coordinates are for the left border -width/2 and for the right border width/2
-    default_x_left = -width_flap / 2
-    default_x_right = width_flap / 2
-    for i in range(vertices_per_side):
-        # displacements right
-        displacements_right[i][0] = displacements_right[i][0] - default_x_right
-        displacements_right[i][1] = displacements_right[i][1] - h_vert * i
-        # displacements left
-        displacements_left[i][0] = displacements_left[i][0] - default_x_left
-        displacements_left[i][1] = displacements_left[i][1] - h_vert * i
-    # bottom is fixed -> x & y displacement 0
-    displacements_left[0][:] = 0
-    displacements_right[0][:] = 0
-
-    return np.concatenate((displacements_left, displacements_right), axis=0)
+def displace_flap(x, y, t, flap_tip_y):
+    x_displ = np.zeros_like(x)
+    y_displ = np.zeros_like(y)
+    # first, get displacement independent of x, only dependent on y and t
+    # maximal displacement at flap tip should be 0.5
+    # initially, flap's displacement is 0
+    x_displ = np.minimum(((np.sin(3*np.pi*t + np.arcsin(-0.95)) +0.95)/ 8) * y / flap_tip_y, 0.5 * y / flap_tip_y)
+
+    displ = np.zeros((len(x), 2))
+    displ[:, 0] = x_displ
+    displ[:, 1] = y_displ
+
+    return displ
+
 
 
 configuration_file_name = "../precice-config.xml"
@@ -116,12 +44,11 @@ def fun_x_coord(y): return (np.cosh(a * y) - 1) / a
 n = 24  # Number of vertices per side
 t = 0
 
-vertices = np.zeros((2 * n, dimensions))
-vertices[:n, 0] = x_left  # vertices at left side of flap
-vertices[n:, 0] = x_right  # vertices at right side of flap
-vertices[:n, 1] = np.linspace(y_bottom, y_top, n)  # equally disrtibuted vertices left
-vertices[n:, 1] = np.linspace(y_bottom, y_top, n)  # equally disrtibuted vertices right
-vertex_ids = interface.set_mesh_vertices(mesh_name, vertices)
+vertices_mid = np.zeros((2*n, dimensions))
+vertices_mid[:n, 1] = np.linspace(y_bottom, y_top, n)
+vertices_mid[n:, 1] = np.linspace(y_bottom, y_top, n)
+
+vertex_ids = interface.set_mesh_vertices(mesh_name, vertices_mid)
 
 interface.initialize()
 # change if necessary
@@ -136,7 +63,7 @@ def fun_x_coord(y): return (np.cosh(a * y) - 1) / a
     precice_dt = interface.get_max_time_step_size()
     dt = min([solver_dt, precice_dt])
     # wiggle the flap
-    write_data = compute_flap_displacements(t, H, W, n)
+    write_data = displace_flap(vertices_mid[:, 0], vertices_mid[:, 1], t, H)
 
     interface.write_data(mesh_name, write_data_name, vertex_ids, write_data)
     interface.advance(dt)

From 8686c5e4d17f4934d0c87473e72e4ef1474c862d Mon Sep 17 00:00:00 2001
From: Niklas Vinnitchenko <niklas.vinnitchenko@tum.de>
Date: Sun, 23 Jun 2024 14:02:39 +0200
Subject: [PATCH 4/7] Format

---
 perpendicular-flap/solid-fake/fake.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/perpendicular-flap/solid-fake/fake.py b/perpendicular-flap/solid-fake/fake.py
index 4f8f161df..5b610f886 100644
--- a/perpendicular-flap/solid-fake/fake.py
+++ b/perpendicular-flap/solid-fake/fake.py
@@ -10,7 +10,7 @@ def displace_flap(x, y, t, flap_tip_y):
     # first, get displacement independent of x, only dependent on y and t
     # maximal displacement at flap tip should be 0.5
     # initially, flap's displacement is 0
-    x_displ = np.minimum(((np.sin(3*np.pi*t + np.arcsin(-0.95)) +0.95)/ 8) * y / flap_tip_y, 0.5 * y / flap_tip_y)
+    x_displ = np.minimum(((np.sin(3 * np.pi * t + np.arcsin(-0.95)) + 0.95) / 8) * y / flap_tip_y, 0.5 * y / flap_tip_y)
 
     displ = np.zeros((len(x), 2))
     displ[:, 0] = x_displ
@@ -19,7 +19,6 @@ def displace_flap(x, y, t, flap_tip_y):
     return displ
 
 
-
 configuration_file_name = "../precice-config.xml"
 participant_name = "Solid"
 mesh_name = "Solid-Mesh"
@@ -44,7 +43,7 @@ def displace_flap(x, y, t, flap_tip_y):
 n = 24  # Number of vertices per side
 t = 0
 
-vertices_mid = np.zeros((2*n, dimensions))
+vertices_mid = np.zeros((2 * n, dimensions))
 vertices_mid[:n, 1] = np.linspace(y_bottom, y_top, n)
 vertices_mid[n:, 1] = np.linspace(y_bottom, y_top, n)
 

From e93232eefbe2f438f3888237e5e37217e4adcd8e Mon Sep 17 00:00:00 2001
From: Niklas <niklas.vinnitchenko@tum.de>
Date: Sat, 29 Jun 2024 00:28:10 +0200
Subject: [PATCH 5/7] Add solid fake to list of solid solvers in README

---
 perpendicular-flap/README.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/perpendicular-flap/README.md b/perpendicular-flap/README.md
index 9a6ae600b..c48a78fb9 100644
--- a/perpendicular-flap/README.md
+++ b/perpendicular-flap/README.md
@@ -51,6 +51,8 @@ Solid participant:
 
 * 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 and this case is only provided for quick testing purposes, leading to outlier results. For general solid mechanics procedures in OpenFOAM, see solids4foam.
 
+* Fake. A simple Python script that acts as a fake solver and provides an arbitrary time-dependent flap displacement in the x-direction, i.e., it performs a shear mapping on the resting flap. This solver can be used for debugging of the fluid participant and its adapter. It also technically works with implicit coupling, thus no changes to the preCICE configuration are necessary. Note that [ASTE's replay mode](https://precice.org/tooling-aste.html#replay-mode) has a similar use case and could also feed artificial or previously recorded real data, replacing an actual solver.
+
 ## Running the Simulation
 
 All listed solvers can be used in order to run the simulation. OpenFOAM can be executed in parallel using `run.sh -parallel`. The default setting uses 4 MPI ranks. Open two separate terminals and start the desired fluid and solid participant by calling the respective run script `run.sh` located in the participant directory. For example:

From 40ca2ab9dfd60fdc351aaeaa2b38ac8cf7bf76e7 Mon Sep 17 00:00:00 2001
From: Niklas <niklas.vinnitchenko@tum.de>
Date: Tue, 13 Aug 2024 13:50:52 +0200
Subject: [PATCH 6/7] Add suggested changes, run and requirements file

---
 perpendicular-flap/solid-fake/fake.py          | 12 +++++++-----
 perpendicular-flap/solid-fake/requirements.txt |  2 ++
 perpendicular-flap/solid-fake/run.sh           |  7 +++++++
 3 files changed, 16 insertions(+), 5 deletions(-)
 create mode 100644 perpendicular-flap/solid-fake/requirements.txt
 create mode 100755 perpendicular-flap/solid-fake/run.sh

diff --git a/perpendicular-flap/solid-fake/fake.py b/perpendicular-flap/solid-fake/fake.py
index 5b610f886..a641a85a9 100644
--- a/perpendicular-flap/solid-fake/fake.py
+++ b/perpendicular-flap/solid-fake/fake.py
@@ -43,11 +43,13 @@ def displace_flap(x, y, t, flap_tip_y):
 n = 24  # Number of vertices per side
 t = 0
 
-vertices_mid = np.zeros((2 * n, dimensions))
-vertices_mid[:n, 1] = np.linspace(y_bottom, y_top, n)
-vertices_mid[n:, 1] = np.linspace(y_bottom, y_top, n)
+vertices = np.zeros((2 * n, dimensions))
+vertices[:n, 1] = np.linspace(y_bottom, y_top, n)
+vertices[n:, 1] = np.linspace(y_bottom, y_top, n)
+vertices[:n, 0] = x_left
+vertices[n:, 0] = x_right
 
-vertex_ids = interface.set_mesh_vertices(mesh_name, vertices_mid)
+vertex_ids = interface.set_mesh_vertices(mesh_name, vertices)
 
 interface.initialize()
 # change if necessary
@@ -62,7 +64,7 @@ def displace_flap(x, y, t, flap_tip_y):
     precice_dt = interface.get_max_time_step_size()
     dt = min([solver_dt, precice_dt])
     # wiggle the flap
-    write_data = displace_flap(vertices_mid[:, 0], vertices_mid[:, 1], t, H)
+    write_data = displace_flap(vertices[:, 0], vertices[:, 1], t, H)
 
     interface.write_data(mesh_name, write_data_name, vertex_ids, write_data)
     interface.advance(dt)
diff --git a/perpendicular-flap/solid-fake/requirements.txt b/perpendicular-flap/solid-fake/requirements.txt
new file mode 100644
index 000000000..c815166a5
--- /dev/null
+++ b/perpendicular-flap/solid-fake/requirements.txt
@@ -0,0 +1,2 @@
+numpy >1, <2
+pyprecice==3
diff --git a/perpendicular-flap/solid-fake/run.sh b/perpendicular-flap/solid-fake/run.sh
new file mode 100755
index 000000000..487ea548d
--- /dev/null
+++ b/perpendicular-flap/solid-fake/run.sh
@@ -0,0 +1,7 @@
+#!/usr/bin/env sh
+set -e -u
+
+python3 -m venv .venv
+. .venv/bin/activate
+pip install -r requirements.txt
+python3 fake.py

From 81f679775ff7801c580f141476da130715b65c36 Mon Sep 17 00:00:00 2001
From: Niklas <niklas.vinnitchenko@tum.de>
Date: Mon, 19 Aug 2024 13:08:06 +0200
Subject: [PATCH 7/7] Add suggested changes

---
 perpendicular-flap/solid-fake/fake.py | 28 ++++++++++++++++++++-------
 1 file changed, 21 insertions(+), 7 deletions(-)

diff --git a/perpendicular-flap/solid-fake/fake.py b/perpendicular-flap/solid-fake/fake.py
index a641a85a9..fd5581292 100644
--- a/perpendicular-flap/solid-fake/fake.py
+++ b/perpendicular-flap/solid-fake/fake.py
@@ -7,12 +7,20 @@
 def displace_flap(x, y, t, flap_tip_y):
     x_displ = np.zeros_like(x)
     y_displ = np.zeros_like(y)
-    # first, get displacement independent of x, only dependent on y and t
-    # maximal displacement at flap tip should be 0.5
-    # initially, flap's displacement is 0
-    x_displ = np.minimum(((np.sin(3 * np.pi * t + np.arcsin(-0.95)) + 0.95) / 8) * y / flap_tip_y, 0.5 * y / flap_tip_y)
-
-    displ = np.zeros((len(x), 2))
+    # get displacement independent of x, only dependent on y and t
+    max_x_displ = 0.5
+    period_fac = 3 * np.pi
+    damping_fac = 8  # damps the amplitude of the sine
+    # defines how much the sine is shifted in y-direction
+    shift = 0.95
+    # wiggles the flap periodically.
+    # the arcsin(-shift) in the sine evaluation is necessary to start at a flap displacement of 0 at t=0
+    # (i.e. sin(arcsin(-shift))+shift = 0)
+    x_displ = np.minimum(((np.sin(period_fac * t + np.arcsin(-shift)) + shift) /
+                         damping_fac), max_x_displ) * y / flap_tip_y
+
+    dimensions = 2
+    displ = np.zeros((len(x), dimensions))
     displ[:, 0] = x_displ
     displ[:, 1] = y_displ
 
@@ -44,13 +52,19 @@ def displace_flap(x, y, t, flap_tip_y):
 t = 0
 
 vertices = np.zeros((2 * n, dimensions))
+# define vertices of flap's left side
 vertices[:n, 1] = np.linspace(y_bottom, y_top, n)
-vertices[n:, 1] = np.linspace(y_bottom, y_top, n)
 vertices[:n, 0] = x_left
+# define vertices of flap's right side
+vertices[n:, 1] = np.linspace(y_bottom, y_top, n)
 vertices[n:, 0] = x_right
 
 vertex_ids = interface.set_mesh_vertices(mesh_name, vertices)
 
+if interface.requires_initial_data():
+    # initially, there should be no displacement
+    interface.write_data(np.zeros_like(vertices))
+
 interface.initialize()
 # change if necessary
 solver_dt = np.inf