Two ellipsoid muscles

two_muscle_case/ellipsoid_muscles/README.md

@@ -0,0 +1,13 @@
+# Bayesian Optimization for two coupled cuboid muscles
+
+## Setup
+- Two ellipsoid muscle geometries. The 3d geometries are created in ```ellipsoid_mesh_generation.py```, and the fibers' geometries are given by ```mesh_left.json``` and ```mesh_right.json```.
+- The solvers for the contraction are coupled mechanics solver and fastmonodomain solver. In the contraction process we set dynamic to `True`, fix the outer end of the muscles and apply a Neumann boundary condition at the inner end of the muscles. The force in the Neumann boundary condition is updated at every timestep and it is different for each muscle. We compute the force such that it mimics the effect of an immaginary tendon connecting the two muscle. We assume the tendon behaves like a spring and compute the force using Hooke's law: $F_{tendon} = k_{tendon} (l − l_0 )$. The spring's constant is $500N/cm$ by default, but can be modified in ```variables/variables.py``` by changing ```tendon_spring_constant``` to the desired value.
+- It uses the electrophysiology CellML model "hodgkin_huxley-razumova" and the incompressible mechanics model "Mooney-Rivlin".
+- No preCICE involved. 
+
+## How to run
+To run a single simulation of contracting the muscles, go to build_release and run:
+```
+./two_muscles_contraction ../settings_contraction.py
+```

two_muscle_case/ellipsoid_muscles/ellipsoid_mesh_generation.py

@@ -0,0 +1,140 @@
+import math
+import matplotlib.pyplot as plt
+
+import numpy as np
+
+def plot_2d_projection(points):
+    x_vals = [p[0] for p in points]
+    y_vals = [p[1] for p in points]
+
+    plt.figure()
+    plt.scatter(x_vals, y_vals, c='red', marker='o')
+
+    # Add index labels
+    for i, (x, y) in enumerate(zip(x_vals, y_vals)):
+        plt.text(x, y, str(i), fontsize=9, ha='right', va='bottom', color='blue')
+
+    plt.xlabel('X')
+    plt.ylabel('Y')
+    plt.title('2D Projection (Z constant)')
+    plt.axis('equal')
+    plt.grid(True)
+    plt.show()
+
+def plot_3d_points(points):
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    # Unpack x, y, z coordinates
+    x_vals = [p[0] for p in points]
+    y_vals = [p[1] for p in points]
+    z_vals = [p[2] for p in points]
+
+    ax.scatter(x_vals, y_vals, z_vals, c='blue', marker='o')
+
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Z')
+    plt.title('3D Point Visualization')
+    plt.show()
+
+# x^2/c^2 + y^2/c^2 + z^2/a^2 = 1
+# 2*r^2/c^2 + z^2/a^2 = 1
+
+# assume el_x is even, and el_x = el_y
+def mesh_nodes(el_x,el_y, el_z,c,a,amin,amax):
+
+    if el_x != el_y:
+        print("Only el_x = el_y supported")
+
+    nodes = []
+    z_points = [(amax-amin)/el_z*i+amin for i in range(el_z+1)]
+    for z in z_points:
+        # r^2 = c²/2*(1-z^2/a^2)
+        rz=np.sqrt(c**2/2*(1-z**2/a**2))
+        circular_nodes = create_points(rz,el_x,z)
+        for node in circular_nodes:
+            nodes.append(node)
+
+
+    return nodes
+
+def compute_outer_points(el):    
+    return (el + 1)*4-4
+
+def create_points(r,el,z):
+
+    nodes = []
+    circular_matrix = []
+    circular_matrix = circle_points(r,(el+1)*4-4,z)
+
+    xpoints = el +1
+
+    if el == 2:
+        # y = 0
+        for i in range(xpoints):
+            nodes.append(circular_matrix[i])
+        # y = 1
+        nodes.append(circular_matrix[-1])
+        nodes.append([0,0,z])
+        nodes.append(circular_matrix[i+1])
+        # y = 2
+        for i in range(xpoints):
+            nodes.append(circular_matrix[-2-i])     
+
+    if el == 4:
+        xpoints_small = el-1
+        circular_matrix_small = circle_points(r/2,(el-1)*4-4,z)
+
+        # y = 0
+        for i in range(xpoints):
+            nodes.append(circular_matrix[i])
+        # y = 1
+        nodes.append(circular_matrix[-1])
+        for j in range(xpoints_small):
+            nodes.append(circular_matrix_small[j])
+        nodes.append(circular_matrix[i+1])
+        # y = 2
+        nodes.append(circular_matrix[-2])
+        nodes.append(circular_matrix_small[-1])
+        nodes.append([0,0,z])
+        nodes.append(circular_matrix_small[j+1])
+        nodes.append(circular_matrix[i+2])
+
+        # y = 3
+        nodes.append(circular_matrix[-3])
+        for j in range(xpoints_small):
+            nodes.append(circular_matrix_small[-2-j])
+        nodes.append(circular_matrix[i+3])
+        # y = 4
+        for i in range(xpoints):
+            nodes.append(circular_matrix[-4-i]) 
+
+    return nodes 
+
+def circle_points(r, points, z):
+    result = []
+    for k in range(points):
+        theta = math.pi + 2 * math.pi * k / points  # Start at π (180°)
+        x = r * math.cos(theta)
+        y = r * math.sin(theta)
+        result.append([x, y, z])
+    return result
+
+
+
+
+
+def apply_offset(nodes,offset):
+    nodes_with_offset = []
+    new_node = [0,0,0]
+    for node in nodes:
+        new_node = [node[0]+offset[0],node[1]+offset[1],node[2]+offset[2]]
+
+        nodes_with_offset.append(new_node)
+    return nodes_with_offset
+
+
+################## 
+
+#plot_2d_projection(create_points(2,4,0))