diff --git a/.github/workflows/test_pvade.yaml b/.github/workflows/test_pvade.yaml
index 3c9e05a4..897d46c5 100644
--- a/.github/workflows/test_pvade.yaml
+++ b/.github/workflows/test_pvade.yaml
@@ -32,6 +32,9 @@ jobs:
       - name: Run pytest
         shell: bash -l {0}
         run: PYTHONPATH=. pytest -sv pvade/tests/
+      - name: test all inputs
+        shell: bash -l {0}
+        run: pytest -sv  test_all_inputs.py 
 
   # Job 2 of 2 - enforce Black formatting
   formatting:
diff --git a/.gitignore b/.gitignore
index 5b92f102..65bed47d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,4 +8,5 @@
 *.h5
 
 docs/_build/
-results/
\ No newline at end of file
+results/input/
+input/
diff --git a/conftest.py b/conftest.py
new file mode 100644
index 00000000..05a95be0
--- /dev/null
+++ b/conftest.py
@@ -0,0 +1,29 @@
+from pathlib import Path
+
+
+def pytest_addoption(parser):
+    parser.addoption(
+        "--input-file",
+        action="store",
+        default=None,
+        help="Run test only for a specific input YAML file",
+    )
+
+
+def pytest_generate_tests(metafunc):
+    if "input_file" not in metafunc.fixturenames:
+        return
+
+    input_file_arg = metafunc.config.getoption("input_file")
+
+    if input_file_arg:
+        metafunc.parametrize("input_file", [Path(input_file_arg)])
+    else:
+        # 🔧 Always resolve input/ relative to location of conftest.py (PVade/)
+        this_dir = Path(__file__).resolve().parent  # PVade/
+        input_dir = this_dir / "input"
+
+        all_files = sorted(input_dir.glob("*.yaml"))
+        if not all_files:
+            raise RuntimeError(f"No input files found in {input_dir}")
+        metafunc.parametrize("input_file", all_files)
diff --git a/input/3d_cyld.yaml b/input/3d_cyld.yaml
index 32ad6386..8e56c2a6 100644
--- a/input/3d_cyld.yaml
+++ b/input/3d_cyld.yaml
@@ -11,26 +11,26 @@ domain:
   y_max: 0.41
   z_min: 0.0
   z_max: 0.41
-  l_char: 2.0 #.02
-pv_array:
-  stream_rows: 1 # not used
-  elevation: 1.5 # not used
-  stream_spacing: 7.0 # not used
-  panel_chord: 2.0 # not used
-  panel_span: 7.0 # not used
-  panel_thickness: 0.1 # not used
-  tracker_angle: 30.0 # not used
+  l_char: 5.0 #.02
+# pv_array:
+#   stream_rows: 1 # not used
+#   elevation: 1.5 # not used
+#   stream_spacing: 7.0 # not used
+#   panel_chord: 2.0 # not used
+#   panel_span: 7.0 # not used
+#   panel_thickness: 0.1 # not used
+#   tracker_angle: 30.0 # not used
 solver:
-  dt: 0.025
-  t_final: 0.1
-  solver1_ksp: cg
+  dt: 0.001
+  t_final: 0.01
+  solver1_ksp: gmres
   solver2_ksp: cg
   solver3_ksp: cg
   solver1_pc: jacobi
   solver2_pc: jacobi
   solver3_pc: jacobi
-  save_text_interval: 0.1
-  save_xdmf_interval: 0.1
+  save_text_interval: 0.001
+  save_xdmf_interval: 0.001
 fluid:
   velocity_profile_type: parabolic
   inflow_coeff: 16 # only necessary/used when applying a parabolic vel profile
diff --git a/input/heated_panels2d.yaml b/input/heated_panels2d.yaml
index a73abe9b..5c2e85ec 100644
--- a/input/heated_panels2d.yaml
+++ b/input/heated_panels2d.yaml
@@ -14,6 +14,7 @@ domain:
   l_char: 0.1
 pv_array:
   stream_rows: 8
+  span_rows: 1
   elevation: 1.5
   stream_spacing: 5.0
   panel_chord: 2.0
diff --git a/input/inflow_input.yaml b/input/inflow_input.yaml
index 7b384866..6f3baf81 100644
--- a/input/inflow_input.yaml
+++ b/input/inflow_input.yaml
@@ -24,7 +24,7 @@ solver:
   solver1_ksp: gmres
   solver2_ksp: gmres
   solver3_ksp: gmres
-  solver4_ksp: gmres
+  # solver4_ksp: gmres
   solver1_pc: hypre
   solver2_pc: hypre
   solver3_pc: hypre
diff --git a/input/panels2d.yaml b/input/panels2d.yaml
index cee63d41..7482b212 100644
--- a/input/panels2d.yaml
+++ b/input/panels2d.yaml
@@ -9,9 +9,10 @@ domain:
   x_max: 50
   y_min: 0
   y_max: 20
-  l_char: 0.9
+  l_char: 0.5
 pv_array:
   stream_rows: 1
+  span_rows: 1
   elevation: 1.5
   stream_spacing: 7.0
   panel_chord: 2.0
diff --git a/input/panels3d.yaml b/input/panels3d.yaml
index 93e4004b..8535cc61 100644
--- a/input/panels3d.yaml
+++ b/input/panels3d.yaml
@@ -56,4 +56,4 @@ structure:
   body_force_y: 0
   body_force_z: -1 #100
   bc_list: [left ]
-  tube_connection: False
\ No newline at end of file
+  tube_connection: False
diff --git a/input/single_heliostat.yaml b/input/single_heliostat.yaml
index 7bb4c39c..c0bd5bf4 100644
--- a/input/single_heliostat.yaml
+++ b/input/single_heliostat.yaml
@@ -11,7 +11,7 @@ domain:
   y_max: 20 # 20+39 39 is panel to panel 
   z_min: 0
   z_max: 40
-  l_char: 1.25 # 1.0
+  l_char: 4 #1.25 # 1.0
 pv_array:
   stream_rows: 1
   elevation: 5.5
@@ -38,7 +38,7 @@ solver:
   save_text_interval: .02 #0.01
   save_xdmf_interval: .02 #0.01
 fluid:
-  time_varying_inflow_bc: false # true
+  # time_varying_inflow_bc: false # true
   u_ref: 2.0
   rho: 1.0
   nu: 1.8e-05
diff --git a/input/test_heatedpanels2d.yaml b/input/test_heatedpanels2d.yaml
deleted file mode 100644
index 1219793e..00000000
--- a/input/test_heatedpanels2d.yaml
+++ /dev/null
@@ -1,73 +0,0 @@
-general:
-  geometry_module: panels2d
-  output_dir: output/test_heatedpanels2d
-  mesh_only: false
-  structural_analysis: True
-  fluid_analysis: True
-  thermal_analysis: True
-  debug_flag: True
-domain:
-  x_min: -0.3
-  x_max: 1.0 #0.7 #1.0 # 50
-  y_min: 0
-  y_max: 0.4 # 20
-  l_char: 0.006 #83 #0.9
-pv_array:
-  stream_rows: 3
-  elevation: 0.15
-  stream_spacing: 0.2 #
-  panel_chord: 0.1 # 2.0
-  panel_span: 7.0
-  panel_thickness: 0.03 # 0.1
-  tracker_angle: 30. #30.0
-solver:
-  dt: .01 # .001
-  t_final: 2.0
-  solver1_ksp: gmres
-  solver2_ksp: cg
-  solver3_ksp: cg
-  solver4_ksp: gmres
-  solver1_pc: jacobi
-  solver2_pc: jacobi
-  solver3_pc: jacobi
-  solver4_pc: lu
-  save_text_interval: 0.01 # must be same as or bigger than dt
-  save_xdmf_interval: 0.01
-fluid:
-  velocity_profile_type: uniform # loglaw
-  u_ref: 0.5 # 0.2 1.0 2.0
-  z0: 0.005 #m
-  d0: 0.05 #m
-  time_varying_inflow_window: 0.0
-  initialize_with_inflow_bc: true
-  # u_ref: 0.5 #0.5
-  nu: 15.89e-5 #0.001
-  g: 9.81
-  beta: 0.00333
-  alpha: 2.25e-05 # high Pe # 0.00225 # low Pe (no stability needed) # 
-  turbulence_model: 
-  periodic: false 
-  bc_y_max: slip # slip noslip free
-  bc_y_min: noslip # slip noslip free
-  T_ambient: 300.0
-  T_bottom: 320.0
-  T0_panel: 340.0
-structure:
-    dt : .01 # 0.002
-    # rho : 10000.0
-    # rho : 10000.0 # 10000.0
-    poissons_ratio: 0.3
-    elasticity_modulus: 1.0e+05
-    body_force_x: 0
-    body_force_y: -1
-    body_force_z: 0 #100
-    bc_list: ["top"]
-    motor_connection: False
-    tube_connection: False
-    beta_relaxation: 0.005
-
-    # elasticity_modulus: 1.0e+09
-    # poissons_ratio: 0.3
-    # body_force_x: 0
-    # body_force_y: 0
-    # bc_list: ["left"]
diff --git a/pvade/IO/DataStream.py b/pvade/IO/DataStream.py
index 8ebacc51..d602af7a 100644
--- a/pvade/IO/DataStream.py
+++ b/pvade/IO/DataStream.py
@@ -1,9 +1,14 @@
-from dolfinx.io import XDMFFile
-
-# import logging
+import dolfinx
+import ufl
 import sys
+
+import numpy as np
+
 from datetime import datetime
 
+# from dolfinx.fem import create_nonmatching_meshes_interpolation_data
+# import logging
+
 
 def start_print_and_log(rank, logfile_name):
 
@@ -74,7 +79,7 @@ class DataStream:
 
     """
 
-    def __init__(self, domain, flow, structure, params):
+    def __init__(self, domain, flow, structure=None, params=None):
         """Initialize the DataStream object
 
         This initializes an object that manages the I/O for all of PVade.
@@ -98,7 +103,9 @@ def __init__(self, domain, flow, structure, params):
                 f"{params.general.output_dir_sol}/solution_fluid.xdmf"
             )
 
-            with XDMFFile(self.comm, self.results_filename_fluid, "w") as xdmf_file:
+            with dolfinx.io.XDMFFile(
+                self.comm, self.results_filename_fluid, "w"
+            ) as xdmf_file:
                 tt = 0.0
                 xdmf_file.write_mesh(domain.fluid.msh)
                 xdmf_file.write_function(flow.u_k, 0.0)
@@ -115,13 +122,31 @@ def __init__(self, domain, flow, structure, params):
                 f"{params.general.output_dir_sol}/solution_structure.xdmf"
             )
 
-            with XDMFFile(self.comm, self.results_filename_structure, "w") as xdmf_file:
+            # Since this is the first call, build solvers that project forms on an as-saved basis
+            # Set up the linear problem used for the projection, cg solver and jacobi pc
+            petsc_options = {"ksp_type": "cg", "pc_type": "jacobi"}
+            self.prob_k_nominal = dolfinx.fem.petsc.LinearProblem(
+                ufl.lhs(structure.elasticity.k_nominal_proj),
+                ufl.rhs(structure.elasticity.k_nominal_proj),
+                bcs=[],
+                petsc_options=petsc_options,
+            )
+
+            # Solve for the stress tensor
+            sol_k_nominal = self.prob_k_nominal.solve()
+            structure.elasticity.internal_stress.x.array[:] = sol_k_nominal.x.array[:]
+            structure.elasticity.internal_stress.x.scatter_forward()
+
+            with dolfinx.io.XDMFFile(
+                self.comm, self.results_filename_structure, "w"
+            ) as xdmf_file:
                 tt = 0.0
                 xdmf_file.write_mesh(domain.structure.msh)
                 xdmf_file.write_function(structure.elasticity.u, 0.0)
-                xdmf_file.write_function(structure.elasticity.stress, 0.0)
                 xdmf_file.write_function(structure.elasticity.v_old, 0.0)
-                xdmf_file.write_function(structure.elasticity.sigma_vm_h, 0.0)
+                xdmf_file.write_function(structure.elasticity.a_old, 0.0)
+                xdmf_file.write_function(structure.elasticity.stress, 0.0)
+                xdmf_file.write_function(structure.elasticity.internal_stress, 0.0)
 
         if self.comm.rank == 0 and self.comm.size > 1 and params.general.test:
             self.log_filename_structure = f"{params.general.output_dir_sol}/log_str.txt"
@@ -188,7 +213,9 @@ def save_XDMF_files(self, fsi_object, domain, tt):
         """
 
         if fsi_object.name == "fluid":
-            with XDMFFile(self.comm, self.results_filename_fluid, "a") as xdmf_file:
+            with dolfinx.io.XDMFFile(
+                self.comm, self.results_filename_fluid, "a"
+            ) as xdmf_file:
                 xdmf_file.write_function(fsi_object.u_k, tt)
                 xdmf_file.write_function(fsi_object.p_k, tt)
                 xdmf_file.write_function(fsi_object.panel_stress, tt)
@@ -197,11 +224,19 @@ def save_XDMF_files(self, fsi_object, domain, tt):
                     xdmf_file.write_function(fsi_object.theta_k, tt)
 
         elif fsi_object.name == "structure":
-            with XDMFFile(self.comm, self.results_filename_structure, "a") as xdmf_file:
+            # Solve for the stress tensor
+            sol_k_nominal = self.prob_k_nominal.solve()
+            fsi_object.elasticity.internal_stress.x.array[:] = sol_k_nominal.x.array[:]
+            fsi_object.elasticity.internal_stress.x.scatter_forward()
+
+            with dolfinx.io.XDMFFile(
+                self.comm, self.results_filename_structure, "a"
+            ) as xdmf_file:
                 xdmf_file.write_function(fsi_object.elasticity.u, tt)
-                xdmf_file.write_function(fsi_object.elasticity.stress, tt)
                 xdmf_file.write_function(fsi_object.elasticity.v_old, tt)
-                xdmf_file.write_function(fsi_object.elasticity.sigma_vm_h, tt)
+                xdmf_file.write_function(fsi_object.elasticity.a_old, tt)
+                xdmf_file.write_function(fsi_object.elasticity.stress, tt)
+                xdmf_file.write_function(fsi_object.elasticity.internal_stress, tt)
 
         else:
             raise ValueError(
diff --git a/pvade/IO/Utilities.py b/pvade/IO/Utilities.py
index 0fbb35d4..250dbcfd 100644
--- a/pvade/IO/Utilities.py
+++ b/pvade/IO/Utilities.py
@@ -23,7 +23,7 @@ def get_input_file():
     return input_file_path
 
 
-def write_metrics(flow, elasticity, prof_filename="profiling.txt"):
+def write_metrics(flow, prof_filename="profiling.txt"):
     with open(prof_filename, "r") as output_file:
         if flow.fluid_analysis == True:
             # solver_line = [line for line in output_file if "(solve)" in line]
diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index 9496ea50..1c993621 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -118,15 +118,15 @@ properties:
     type: "object"
     properties:
       stream_rows:
-        default: 3
-        minimum: 1
+        default: 0
+        minimum: 0
         maximum: 10
         type: "integer"
         description: "The number of panel rows in the streamwise direction."
         units: "unitless"
       span_rows:
-        default: 1
-        minimum: 1
+        default: 0
+        minimum: 0
         maximum: 10
         type: "integer"
         description: "The number of panel rows in the spanwise direction."
@@ -233,6 +233,15 @@ properties:
         enum:
           - "cg"
           - "gmres"
+          - "bicg"
+          - "bicgstab"
+          - "minres"
+          - "preonly"
+          - "richardson"
+          - "chebyshev"
+          - "tfqmr"
+          - "pipecg"
+          - "pipegmres"
       solver2_ksp:
         default: "cg"
         type: "string"
@@ -240,6 +249,15 @@ properties:
         enum:
           - "cg"
           - "gmres"
+          - "bicg"
+          - "bicgstab"
+          - "minres"
+          - "preonly"
+          - "richardson"
+          - "chebyshev"
+          - "tfqmr"
+          - "pipecg"
+          - "pipegmres"
       solver3_ksp:
         default: "cg"
         type: "string"
@@ -247,55 +265,132 @@ properties:
         enum:
           - "cg"
           - "gmres"
+          - "bicg"
+          - "bicgstab"
+          - "minres"
+          - "preonly"
+          - "richardson"
+          - "chebyshev"
+          - "tfqmr"
+          - "pipecg"
+          - "pipegmres"
       solver4_ksp:
         default: "gmres"
         type: "string"
         description: "The solver to use for fluid solve 4 of 5"
         enum:
-          - "cg" # not tested
+          - "cg"
           - "gmres"
-          - "preonly" # for debugging
+          - "bicg"
+          - "bicgstab"
+          - "minres"
+          - "preonly"
+          - "richardson"
+          - "chebyshev"
+          - "tfqmr"
+          - "pipecg"
+          - "pipegmres"
       solver5_ksp:
         default: "cg"
         type: "string"
         description: "The solver to use for fluid solve 5 of 5"
         enum:
           - "cg"
-          - "gmres"    
+          - "gmres"
+          - "bicg"
+          - "bicgstab"
+          - "minres"
+          - "preonly"
+          - "richardson"
+          - "chebyshev"
+          - "tfqmr"
+          - "pipecg"
+          - "pipegmres" 
       solver1_pc:
         default: "hypre"
         type: "string"
         description: "The preconditioner to use for fluid solve 1 of 5"
         enum:
+          - "none"
           - "jacobi"
+          - "sor"
+          - "ilu"
+          - "icc"
+          - "asm"
+          - "gamg"
           - "hypre"
+          - "ml"
+          - "lu"
+          - "cholesky"
+          - "python"
       solver2_pc:
         default: "hypre"
         type: "string"
         description: "The preconditioner to use for fluid solve 2 of 5"
         enum:
+          - "none"
           - "jacobi"
+          - "sor"
+          - "ilu"
+          - "icc"
+          - "asm"
+          - "gamg"
           - "hypre"
+          - "ml"
+          - "lu"
+          - "cholesky"
+          - "python"
       solver3_pc:
         default: "hypre"
         type: "string"
         description: "The preconditioner to use for fluid solve 3 of 5"
         enum:
+          - "none"
           - "jacobi"
+          - "sor"
+          - "ilu"
+          - "icc"
+          - "asm"
+          - "gamg"
           - "hypre"
+          - "ml"
+          - "lu"
+          - "cholesky"
+          - "python"
       solver4_pc:
         default: "lu"
         type: "string"
         description: "The preconditioner to use for fluid solve 4 of 5"
         enum:
+          - "none"
+          - "jacobi"
+          - "sor"
+          - "ilu"
+          - "icc"
+          - "asm"
+          - "gamg"
+          - "hypre"
+          - "ml"
           - "lu"
+          - "cholesky"
+          - "python"
       solver5_pc:
         default: "hypre"
         type: "string"
         description: "The preconditioner to use for fluid solve 5 of 5"
         enum:
+          - "none"
           - "jacobi"
-          - "hypre"    
+          - "sor"
+          - "ilu"
+          - "icc"
+          - "asm"
+          - "gamg"
+          - "hypre"
+          - "ml"
+          - "lu"
+          - "cholesky"
+          - "python" 
   fluid:
     additionalProperties: false
     type: "object"
diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index 5cd3f364..f31df144 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -80,15 +80,32 @@ def __init__(self, domain, params):
 
             # find hmin in mesh
             num_cells = domain.fluid.msh.topology.index_map(self.ndim).size_local
+
             h = dolfinx.cpp.mesh.h(domain.fluid.msh, self.ndim, range(num_cells))
 
-            # This value of hmin is local to the mesh portion owned by the process
-            hmin_local = np.amin(h)
+            # # This value of hmin is local to the mesh portion owned by the process
+            # hmin_local = np.amin(h)
+
+            if len(h) > 0:
+                hmin_local = np.amin(h)
+            else:
+                hmin_local = np.inf
 
-            # collect the minimum hmin from all processes
+            print(hmin_local)
             self.hmin = np.zeros(1)
-            self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN)
-            self.hmin = self.hmin[0]
+            self.hmin = self.comm.allreduce(hmin_local, op=MPI.MIN)
+
+            # ------------------------------------------------------------------------------
+
+            # h = dolfinx.cpp.mesh.h(domain.fluid.msh, self.ndim, range(num_cells))
+
+            # # This value of hmin is local to the mesh portion owned by the process
+            # hmin_local = np.amin(h)
+
+            # # collect the minimum hmin from all processes
+            # self.hmin = np.zeros(1)
+            # self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN)
+            # self.hmin = self.hmin[0]
 
             self.num_Q_dofs = (
                 self.Q.dofmap.index_map_bs * self.Q.dofmap.index_map.size_global
@@ -98,7 +115,6 @@ def __init__(self, domain, params):
             )
 
             if self.rank == 0:
-                print(f"hmin on fluid  = {self.hmin}")
                 print(f"Total num dofs on fluid= {self.num_Q_dofs + self.num_V_dofs}")
 
     def build_boundary_conditions(self, domain, params):
@@ -734,7 +750,9 @@ def solve(self, domain, params, current_time):
 
         self.compute_lift_and_drag(params, current_time)
 
-        self.compute_pressure_drop_between_points(domain, params)
+        # Compute the pressure drop between the inlet and outlet
+        if params.pv_array.stream_rows > 0:
+            self.compute_pressure_drop_between_points(domain, params)
 
         # Update new -> old variables
         self.u_k2.x.array[:] = self.u_k1.x.array
@@ -946,11 +964,24 @@ def compute_cfl(self):
         self.solver_6.solve(self.b6, self.cfl_vec.vector)
         self.cfl_vec.x.scatter_forward()
 
-        cfl_max_local = np.amax(self.cfl_vec.vector.array)
+        # cfl_max_local = np.amax(self.cfl_vec.vector.array)
+
+        # # collect the minimum hmin from all processes
+        # self.cfl_max = np.zeros(1)
+        # self.comm.Allreduce(cfl_max_local, self.cfl_max, op=MPI.MAX)
+        # self.cfl_max = self.cfl_max[0]
 
-        # collect the minimum hmin from all processes
-        self.cfl_max = np.zeros(1)
-        self.comm.Allreduce(cfl_max_local, self.cfl_max, op=MPI.MAX)
+        # Compute local max only if array has values
+        if self.cfl_vec.vector.array.size > 0:
+            cfl_max_local = np.amax(self.cfl_vec.vector.array)
+        else:
+            cfl_max_local = -np.inf  # So it won't affect global max
+
+        # Prepare buffer and reduce across all ranks
+        self.cfl_max = np.zeros(1, dtype=np.float64)
+        self.comm.Allreduce(
+            np.array(cfl_max_local, dtype=np.float64), self.cfl_max, op=MPI.MAX
+        )
         self.cfl_max = self.cfl_max[0]
 
     def compute_lift_and_drag(self, params, current_time):
diff --git a/pvade/fluid/boundary_conditions.py b/pvade/fluid/boundary_conditions.py
index 1d82be98..f4e68193 100644
--- a/pvade/fluid/boundary_conditions.py
+++ b/pvade/fluid/boundary_conditions.py
@@ -309,8 +309,8 @@ def build_velocity_boundary_conditions(domain, params, functionspace, current_ti
         domain, params, functionspace, current_time
     )
 
-    if domain.rank == 0:
-        print("inflow_function = ", inflow_function.x.array[:])
+    # if domain.rank == 0:
+    #     print("inflow_function = ", inflow_function.x.array[:])
 
     dofs = get_facet_dofs_by_gmsh_tag(domain, functionspace, "x_min")
     bcu.append(dolfinx.fem.dirichletbc(inflow_function, dofs))
diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py
index be30aab3..13f1764d 100644
--- a/pvade/geometry/MeshManager.py
+++ b/pvade/geometry/MeshManager.py
@@ -214,15 +214,17 @@ def build(self, params):
         self.cell_tags.name = "cell_tags"
         self.facet_tags.name = "facet_tags"
 
-        if (
-            params.general.geometry_module == "panels3d"
-            or params.general.geometry_module == "heliostats3d"
-            or params.general.geometry_module == "flag2d"
-            or params.general.geometry_module == "panels2d"
-        ):
-            self._create_submeshes_from_parent(params)
-            self._transfer_mesh_tags_to_submeshes(params)
-
+        # if (
+        #     params.general.geometry_module == "panels3d"
+        #     or params.general.geometry_module == "heliostats3d"
+        #     or params.general.geometry_module == "flag2d"
+        #     or params.general.geometry_module == "panels2d"
+        # ):
+        #     self._create_submeshes_from_parent(params)
+        #     self._transfer_mesh_tags_to_submeshes(params)
+
+        self._create_submeshes_from_parent(params)
+        self._transfer_mesh_tags_to_submeshes(params)
         self._save_submeshes_for_reload_hack(params)
 
         if params.general.fluid_analysis == True:
@@ -253,67 +255,6 @@ def _save_submeshes_for_reload_hack(self, params):
         self.write_mesh_files(params)
         self.read_mesh_files(params.general.output_dir_mesh, params)
 
-        # if params.general.fluid_analysis == True and params.general.structural_analysis == True :
-        #     sub_domain_list = ["fluid", "structure"]
-        #     for sub_domain_name in sub_domain_list:
-        #         mesh_filename = os.path.join(params.general.output_dir_mesh, f"temp_{sub_domain_name}.xdmf")
-        #         sub_domain = getattr(self, sub_domain_name)
-        #         sub_domain.msh.name = "temp_mesh"
-
-        #         # Write this submesh to immediately read it
-        #         with dolfinx.io.XDMFFile(self.comm, mesh_filename, "w") as xdmf:
-        #             xdmf.write_mesh(sub_domain.msh)
-        #             xdmf.write_meshtags(sub_domain.cell_tags)
-        #             sub_domain.msh.topology.create_connectivity(self.ndim-1, self.ndim)
-        #             xdmf.write_meshtags(sub_domain.facet_tags)
-
-        #         # Read the just-written mesh
-        #         with dolfinx.io.XDMFFile(self.comm, mesh_filename, "r") as xdmf:
-        #             submesh = xdmf.read_mesh(name=sub_domain.msh.name)
-        #             cell_tags = xdmf.read_meshtags(submesh, name="cell_tags")
-        #             ndim = submesh.topology.dim
-        #             submesh.topology.create_connectivity(ndim - 1, ndim)
-        #             facet_tags = xdmf.read_meshtags(submesh, name="facet_tags")
-
-        #         sub_domain.msh = submesh
-        #         sub_domain.facet_tags = facet_tags
-        #         sub_domain.cell_tags = cell_tags
-
-        # elif params.general.geometry_module == "panels3d" and params.general.fluid_analysis == False and params.general.structural_analysis == True :
-        #         sub_domain_name = "structure"
-
-        #         mesh_filename = os.path.join(params.general.output_dir_mesh, f"temp_{sub_domain_name}.xdmf")
-        #         # sub_domain = getattr(self, sub_domain_list[0])
-        #         # sub_domain.msh.name = "temp_mesh"
-
-        #         # self.msh.name  = "temp_mesh"
-
-        #         # Write this submesh to immediately read it
-        #         with dolfinx.io.XDMFFile(self.comm, mesh_filename, "w") as xdmf:
-        #             xdmf.write_mesh(self.msh)
-        #             xdmf.write_meshtags(self.cell_tags)
-        #             self.msh.topology.create_connectivity(self.ndim-1, self.ndim)
-        #             xdmf.write_meshtags(self.facet_tags)
-
-        #         # Read the just-written mesh
-        #         with dolfinx.io.XDMFFile(self.comm, mesh_filename, "r") as xdmf:
-        #             mesh = xdmf.read_mesh(name=self.msh.name)
-        #             cell_tags = xdmf.read_meshtags(mesh, name="cell_tags")
-        #             ndim = mesh.topology.dim
-        #             mesh.topology.create_connectivity(ndim - 1, ndim)
-        #             facet_tags = xdmf.read_meshtags(mesh, name="facet_tags")
-
-        #         class FSISubDomain:
-        #             pass
-
-        #         sub_domain = FSISubDomain()
-
-        #         sub_domain.msh = mesh
-        #         sub_domain.facet_tags = facet_tags
-        #         sub_domain.cell_tags = cell_tags
-
-        #         setattr(self, sub_domain_name, sub_domain)
-
     def _create_submeshes_from_parent(self, params):
         """Create submeshes from a parent mesh by cell tags.
 
@@ -620,86 +561,92 @@ def read_mesh_files(self, read_mesh_dir, params):
             mesh_name = f"{sub_domain_name}_mesh.xdmf"
             mesh_filename = os.path.join(read_mesh_dir, mesh_name)
 
-            try:
-                if self.rank == 0:
-                    print(f"Reading {sub_domain_name} mesh.")
-
-                # Read the subdomain mesh
-                with dolfinx.io.XDMFFile(self.comm, mesh_filename, "r") as xdmf:
-                    submesh = xdmf.read_mesh(name=mesh_name)
-                    cell_tags = xdmf.read_meshtags(submesh, name="cell_tags")
-                    ndim = submesh.topology.dim
-                    submesh.topology.create_connectivity(ndim - 1, ndim)
-                    facet_tags = xdmf.read_meshtags(submesh, name="facet_tags")
-
-            except:
-                if self.rank == 0:
-                    print(
-                        f"Could not find subdomain {sub_domain_name} mesh file, not reading this mesh."
-                    )
+            if (
+                params.general.fluid_analysis == True
+                and sub_domain_name == "fluid"
+                or params.general.structural_analysis == True
+                and sub_domain_name == "structure"
+            ):
+                try:
+                    if self.rank == 0:
+                        print(f"Reading {sub_domain_name} mesh.")
+
+                    # Read the subdomain mesh
+                    with dolfinx.io.XDMFFile(self.comm, mesh_filename, "r") as xdmf:
+                        submesh = xdmf.read_mesh(name=mesh_name)
+                        cell_tags = xdmf.read_meshtags(submesh, name="cell_tags")
+                        ndim = submesh.topology.dim
+                        submesh.topology.create_connectivity(ndim - 1, ndim)
+                        facet_tags = xdmf.read_meshtags(submesh, name="facet_tags")
+
+                except:
+                    if self.rank == 0:
+                        print(
+                            f"Could not find subdomain {sub_domain_name} mesh file, not reading this mesh."
+                        )
 
-            else:
+                else:
 
-                class FSISubDomain:
-                    pass
+                    class FSISubDomain:
+                        pass
 
-                if self.ndim is None:
-                    self.ndim = submesh.topology.dim
-                else:
-                    assert self.ndim == submesh.topology.dim
+                    if self.ndim is None:
+                        self.ndim = submesh.topology.dim
+                    else:
+                        assert self.ndim == submesh.topology.dim
 
-                submesh.topology.create_connectivity(self.ndim, self.ndim - 1)
+                    submesh.topology.create_connectivity(self.ndim, self.ndim - 1)
 
-                sub_domain = FSISubDomain()
+                    sub_domain = FSISubDomain()
 
-                sub_domain.msh = submesh
-                sub_domain.cell_tags = cell_tags
-                sub_domain.facet_tags = facet_tags
+                    sub_domain.msh = submesh
+                    sub_domain.cell_tags = cell_tags
+                    sub_domain.facet_tags = facet_tags
 
-                # These elements do not need to be created when reading a mesh
-                # they are only used in the transfer of facet tags, and since
-                # those can be read directly from a file, we don't need these
-                sub_domain.entity_map = None
-                sub_domain.vertex_map = None
-                sub_domain.geom_map = None
+                    # These elements do not need to be created when reading a mesh
+                    # they are only used in the transfer of facet tags, and since
+                    # those can be read directly from a file, we don't need these
+                    sub_domain.entity_map = None
+                    sub_domain.vertex_map = None
+                    sub_domain.geom_map = None
 
-                setattr(self, sub_domain_name, sub_domain)
+                    setattr(self, sub_domain_name, sub_domain)
 
-                if sub_domain_name == "fluid":
-                    domain_ufl = ufl.Mesh(
-                        self.fluid.msh.ufl_domain().ufl_coordinate_element()
-                    )
-                    fluid_undeformed = dolfinx.mesh.Mesh(
-                        self.comm,
-                        self.fluid.msh.topology,
-                        self.fluid.msh.geometry,
-                        domain_ufl,
-                    )
+                    if sub_domain_name == "fluid":
+                        domain_ufl = ufl.Mesh(
+                            self.fluid.msh.ufl_domain().ufl_coordinate_element()
+                        )
+                        fluid_undeformed = dolfinx.mesh.Mesh(
+                            self.comm,
+                            self.fluid.msh.topology,
+                            self.fluid.msh.geometry,
+                            domain_ufl,
+                        )
 
-                    fluid_undeformed.topology.create_connectivity(
-                        self.ndim, self.ndim - 1
-                    )
+                        fluid_undeformed.topology.create_connectivity(
+                            self.ndim, self.ndim - 1
+                        )
 
-                    sub_domain_undeformed = FSISubDomain()
+                        sub_domain_undeformed = FSISubDomain()
 
-                    sub_domain_undeformed.msh = fluid_undeformed
-                    sub_domain_undeformed.cell_tags = cell_tags
-                    sub_domain_undeformed.facet_tags = facet_tags
+                        sub_domain_undeformed.msh = fluid_undeformed
+                        sub_domain_undeformed.cell_tags = cell_tags
+                        sub_domain_undeformed.facet_tags = facet_tags
 
-                    # These elements do not need to be created when reading a mesh
-                    # they are only used in the transfer of facet tags, and since
-                    # those can be read directly from a file, we don't need these
-                    sub_domain_undeformed.entity_map = None
-                    sub_domain_undeformed.vertex_map = None
-                    sub_domain_undeformed.geom_map = None
+                        # These elements do not need to be created when reading a mesh
+                        # they are only used in the transfer of facet tags, and since
+                        # those can be read directly from a file, we don't need these
+                        sub_domain_undeformed.entity_map = None
+                        sub_domain_undeformed.vertex_map = None
+                        sub_domain_undeformed.geom_map = None
 
-                    setattr(self, "fluid_undeformed", sub_domain_undeformed)
+                        setattr(self, "fluid_undeformed", sub_domain_undeformed)
 
-                    # assert np.all(self.fluid.msh.geometry.x[:] == self.fluid_undeformed.msh.geometry.x[:])
-                    # assert np.shape(self.fluid.msh.geometry.x[:]) == np.shape(self.fluid_undeformed.msh.geometry.x[:])
+                        # assert np.all(self.fluid.msh.geometry.x[:] == self.fluid_undeformed.msh.geometry.x[:])
+                        # assert np.shape(self.fluid.msh.geometry.x[:]) == np.shape(self.fluid_undeformed.msh.geometry.x[:])
 
-            if self.rank == 0:
-                print(f"Finished read {sub_domain_name} mesh.")
+                if self.rank == 0:
+                    print(f"Finished read {sub_domain_name} mesh.")
 
         self.ndim = submesh.topology.dim
 
diff --git a/pvade/geometry/cylinder2d/DomainCreation.py b/pvade/geometry/cylinder2d/DomainCreation.py
index 525820e7..cb17ddef 100644
--- a/pvade/geometry/cylinder2d/DomainCreation.py
+++ b/pvade/geometry/cylinder2d/DomainCreation.py
@@ -21,7 +21,7 @@ def __init__(self, params):
         """
         super().__init__(params)
 
-    def build(self, params):
+    def build_FSI(self, params):
         """This function creates the computational domain for a flow around a 2D cylinder.
 
         Returns:
@@ -45,3 +45,5 @@ def build(self, params):
         self.gmsh_model.occ.cut([(2, rectangle)], [(2, obstacle)])
 
         self.gmsh_model.occ.synchronize()
+
+        self.numpy_pt_total_array = np.zeros((3, 6))
diff --git a/pvade/geometry/cylinder3d/DomainCreation.py b/pvade/geometry/cylinder3d/DomainCreation.py
index b54cb8d5..e139c42b 100644
--- a/pvade/geometry/cylinder3d/DomainCreation.py
+++ b/pvade/geometry/cylinder3d/DomainCreation.py
@@ -20,7 +20,7 @@ def __init__(self, params):
         """
         super().__init__(params)
 
-    def build(self, params):
+    def build_FSI(self, params):
         """This function creates the computational domain for a flow around a 3D cylinder.
 
         Returns:
@@ -52,6 +52,7 @@ def build(self, params):
         self.gmsh_model.occ.cut([(3, domain)], [(3, cylinder)])
 
         self.gmsh_model.occ.synchronize()
+        self.numpy_pt_total_array = np.zeros((3, 6))
 
     def set_length_scales(self, params, domain_markers):
         res_min = params.domain.l_char
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 13ecba2d..4a84d988 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -40,9 +40,6 @@ def __init__(self, domain, structural_analysis, params):
         self.rank = domain.rank
         self.num_procs = domain.num_procs
 
-        # # domain.structure.msh = dolfinx.mesh.refine(domain.structure.msh,None)
-        # # domain.structure.msh = dolfinx.mesh.refine(domain.structure.msh)
-
         P1 = ufl.VectorElement("Lagrange", domain.structure.msh.ufl_cell(), 2)
         self.V = dolfinx.fem.FunctionSpace(domain.structure.msh, P1)
 
@@ -69,85 +66,6 @@ def __init__(self, domain, structural_analysis, params):
         # time step
         self.dt_st = dolfinx.fem.Constant(domain.structure.msh, (params.structure.dt))
 
-        # # Store the dimension of the problem for convenience
-        # self.ndim = domain.structure.msh.topology.dim
-        # self.facet_dim = self.ndim - 1
-
-        # # find hmin in mesh
-        # num_cells = domain.structure.msh.topology.index_map(self.ndim).size_local
-        # h = dolfinx.cpp.mesh.h(domain.structure.msh, self.ndim, range(num_cells))
-
-        # # This value of hmin is local to the mesh portion owned by the process
-        # hmin_local = np.amin(h)
-
-        # # collect the minimum hmin from all processes
-        # self.hmin = np.zeros(1)
-        # self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN)
-        # self.hmin = self.hmin[0]
-
-        # if self.rank == 0:
-        #     print(f"hmin on structure = {self.hmin}")
-        #     print(f"Total num dofs on structure = {self.num_V_dofs}")
-
-        # # Mass density
-        # self.rho = dolfinx.fem.Constant(
-        #     domain.structure.msh, params.structure.rho
-        # )  # Constant(0.)
-
-        # # Rayleigh damping coefficients
-        # self.eta_m = dolfinx.fem.Constant(domain.structure.msh, 0.0)  # Constant(0.)
-        # self.eta_k = dolfinx.fem.Constant(domain.structure.msh, 0.0)  # Constant(0.)
-
-        # # Generalized-alpha method parameters
-        # self.alpha_m = dolfinx.fem.Constant(domain.structure.msh, 0.2)
-        # self.alpha_f = dolfinx.fem.Constant(domain.structure.msh, 0.4)
-        # self.gamma = 0.5 + self.alpha_f - self.alpha_m
-        # self.beta = (self.gamma + 0.5) ** 2 / 4.0
-
-        # # Define structural properties
-        # self.E = params.structure.elasticity_modulus  # 1.0e9
-        # self.ν = params.structure.poissons_ratio  # 0.3
-        # self.μ = self.E / (2.0 * (1.0 + self.ν))
-        # self.λ = self.E * self.ν / ((1.0 + self.ν) * (1.0 - 2.0 * self.ν))
-
-        # if self.rank == 0:
-        #     print(
-        #         f"mu = {self.μ} lambda = {self.λ} E = {self.E} nu = {self.ν} density = {self.rho.value}"
-        #     )
-
-        # # time step
-        # self.dt_st = dolfinx.fem.Constant(domain.structure.msh, (params.structure.dt))
-
-        # def _north_east_corner(x):
-        #     eps = 1.0e-6
-
-        #     # TODO: allow the probing of (x,y,z) points as something specified in the yaml file
-        #     if params.general.geometry_module == "flag2d":
-        #         # TEMP Hack for Turek and Hron Flag
-        #         x1 = 0.6
-        #         x2 = 0.2
-        #         corner = [x1, x2]
-        #     else:
-        #         tracker_angle_rad = np.radians(params.pv_array.tracker_angle)
-        #         x1 = 0.5 * params.pv_array.panel_chord * np.cos(tracker_angle_rad)
-        #         x2 = 0.5 * params.pv_array.panel_thickness * np.sin(tracker_angle_rad)
-        #         corner = [x1 - x2, 0.5 * params.pv_array.panel_span]
-
-        #     east_edge = np.logical_and(corner[0] - eps < x[0], x[0] < corner[0] + eps)
-        #     north_edge = np.logical_and(corner[1] - eps < x[1], x[1] < corner[1] + eps)
-
-        #     north_east_corner = np.logical_and(east_edge, north_edge)
-
-        #     return north_east_corner
-
-        # north_east_corner_facets = dolfinx.mesh.locate_entities_boundary(
-        #     domain.structure.msh, 0, _north_east_corner
-        # )
-
-        # self.north_east_corner_dofs = dolfinx.fem.locate_dofs_topological(
-        #     self.V, 0, north_east_corner_facets
-        # )
-
     def build_boundary_conditions(self, domain, params):
         """Build the boundary conditions
 
@@ -232,70 +150,68 @@ def build_forms(self, domain, params, structure):
         P3 = ufl.TensorElement("Lagrange", domain.structure.msh.ufl_cell(), 2)
         self.T = dolfinx.fem.FunctionSpace(domain.structure.msh, P3)
 
+        self.trial_tensor = ufl.TrialFunction(self.T)
+        self.test_tensor = ufl.TestFunction(self.T)
+        self.internal_stress = dolfinx.fem.Function(self.T, name="stress_structure")
+
         self.stress = dolfinx.fem.Function(self.T, name="stress_fluid")
         self.stress_old = dolfinx.fem.Function(self.T, name="stress_fluid_old")
         self.stress_predicted = dolfinx.fem.Function(
             self.T, name="stress_fluid_predicted"
         )
 
-        self.sigma_vm_h = dolfinx.fem.Function(self.W, name="Stress")
+        # self.sigma_vm_h = dolfinx.fem.Function(self.W, name="Stress")
 
         # discplacement
-        self.u = dolfinx.fem.Function(self.V, name="Deformation")
-        self.u_old = dolfinx.fem.Function(self.V, name="Deformation_old")
-        self.u_delta = dolfinx.fem.Function(self.V, name="Deformation_change")
+        self.u = dolfinx.fem.Function(self.V, name="deformation")
+        self.u_old = dolfinx.fem.Function(self.V, name="deformation_old")
+        self.u_delta = dolfinx.fem.Function(self.V, name="deformation_change")
 
         # velocity
-        self.v = dolfinx.fem.Function(self.V, name="Velocity")
-        self.v_old = dolfinx.fem.Function(self.V, name="Velocity_old")
+        self.v = dolfinx.fem.Function(self.V)
+        self.v_old = dolfinx.fem.Function(self.V, name="velocity")
 
         # acceleration
-        self.a = dolfinx.fem.Function(self.V, name="acceleration")
-        self.a_old = dolfinx.fem.Function(self.V, name="acceleration_old")
+        self.a = dolfinx.fem.Function(self.V)
+        self.a_old = dolfinx.fem.Function(self.V, name="acceleration")
 
         # dss = ufl.ds(subdomain_data=boundary_subdomains)
 
-        # # Stress tensor
         # def sigma(r):
-        #     return dolfinx.fem.form(2.0*self.μ*ufl.sym(ufl.grad(r)) + self.λ *ufl.tr(ufl.sym(ufl.grad(r)))*ufl.Identity(len(r)))
+        #     return dolfinx.fem.form(2.0*self.lame_mu*ufl.sym(ufl.grad(r)) + self.lame_lambda *ufl.tr(ufl.sym(ufl.grad(r)))*ufl.Identity(len(r)))
 
         # # Mass form
         # def m(u, u_):
         #     return dolfinx.fem.form(self.rho*ufl.inner(u, u_)*ufl.dx)
 
         # # Elastic stiffness form
-        # def k(u, u_):
+        # def k_nominal(u, u_):
         #     return dolfinx.fem.form(ufl.inner(sigma(u), ufl.sym(ufl.grad(u_)))*ufl.dx)
 
         # # Rayleigh damping form
         # def c(u, u_):
-        #     return dolfinx.fem.form(self.eta_m*m(u, u_) + self.eta_k*k(u, u_))
+        #     return dolfinx.fem.form(self.eta_m*m(u, u_) + self.eta_k*k_nominal(u, u_))
 
         # # Work of external forces
         # def Wext(u_):
         #     return ufl.dot(u_, self.f)*self.ds #dss(3)
 
-        def sigma(r):
-            return structure.λ * ufl.nabla_div(r) * ufl.Identity(
-                len(r)
-            ) + 2 * structure.μ * ufl.sym(ufl.grad(r))
+        # def sigma(r):
+        #     return structure.lame_lambda * ufl.nabla_div(r) * ufl.Identity(
+        #         len(r)
+        #     ) + 2 * structure.lame_mu * ufl.sym(ufl.grad(r))
 
         def m(u, u_):
-            return structure.rho * ufl.inner(u, u_) * ufl.dx
+            return structure.rho * ufl.inner(u, u_)
 
-        def k(u, u_):
-            return ufl.inner(sigma(u), ufl.grad(u_)) * ufl.dx
+        def c(u, u_):
+            return self.eta_m * m(u, u_) + self.eta_k * k_nominal(u, u_)
 
-        def k(u, u_):
-            # return ufl.inner(sigma(u),ufl.grad(u_))*ufl.dx
-            # updated lagrangian form
-            F = ufl.grad(u) + ufl.Identity(len(u))
-            E = 0.5 * (F.T * F - ufl.Identity(len(u)))
-            # S = self.λ *ufl.tr(E)*ufl.Identity(len(u))   + 2*self.μ * (E - ufl.tr(E)*ufl.Identity(len(u))  /3.0)
-            S = structure.λ * ufl.tr(E) * ufl.Identity(len(u)) + 2 * structure.μ * (E)
+        # def k_cauchy(u, u_):
+        #     return ufl.inner(sigma(u), ufl.grad(u_))
 
-            # return ufl.inner(F * S, ufl.grad(u_)) * ufl.dx
-            return ufl.inner(P_(u), ufl.grad(u_)) * ufl.dx
+        def k_nominal(u, u_):
+            return ufl.inner(P_(u), ufl.grad(u_))
 
         # The deformation gradient, F = I + dy/dX
         def F_(u):
@@ -323,7 +239,7 @@ def S_(u):
             # return lamda * ufl.tr(E) * I + 2.0 * mu * (E - ufl.tr(E) * I / 3.0)
             # TODO: Why does the above form give a better result and where does it come from?
 
-            S_svk = structure.λ * ufl.tr(E) * I + 2.0 * structure.μ * E
+            S_svk = structure.lame_lambda * ufl.tr(E) * I + 2.0 * structure.lame_mu * E
             return S_svk
 
         # The first Piola–Kirchhoff stress tensor, P = F*S
@@ -333,16 +249,13 @@ def P_(u):
             # return ufl.inv(F) * S
             return F * S
 
-        def c(u, u_):
-            return self.eta_m * m(u, u_) + self.eta_k * k(u, u_)
-
         # self.uh_exp = dolfinx.fem.Function(self.V,  name="Deformation")
 
-        def σ(v):
-            """Return an expression for the stress σ given a displacement field"""
-            return 2.0 * structure.μ * ufl.sym(ufl.grad(v)) + structure.λ * ufl.tr(
-                ufl.sym(ufl.grad(v))
-            ) * ufl.Identity(len(v))
+        # def σ(v):
+        #     """Return an expression for the stress σ given a displacement field"""
+        #     return 2.0 * structure.lame_mu * ufl.sym(ufl.grad(v)) + structure.lame_lambda * ufl.tr(
+        #         ufl.sym(ufl.grad(v))
+        #     ) * ufl.Identity(len(v))
 
         # source term ($f = \rho \omega^2 [x_0, \, x_1]$)
         # self.ω, self.ρ = 300.0, 10.0
@@ -390,9 +303,9 @@ def σ(v):
         F = ufl.grad(self.u) + ufl.Identity(len(self.u))
         J = ufl.det(F)
         self.res = (
-            m(self.avg(self.a_old, a_new, self.alpha_m), self.u_)
-            + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_)
-            + k(self.avg(self.u_old, self.u, self.alpha_f), self.u_)
+            m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * ufl.dx
+            + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * ufl.dx
+            + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * ufl.dx
             - structure.rho * ufl.inner(self.f, self.u_) * ufl.dx
             - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_)
             * self.ds
@@ -401,6 +314,16 @@ def σ(v):
         # self.a = dolfinx.fem.form(ufl.lhs(res))
         # self.L = dolfinx.fem.form(ufl.rhs(res))
 
+        # Save a form to project the first Piola–Kirchhoff, P_, stress tensor in the structure
+        # u * v * dx = P_ * v * dx, where u and v are trial and test functions on tensor function space
+        F_k_nominal_proj = ufl.inner(self.trial_tensor, self.test_tensor) * ufl.dx
+        F_k_nominal_proj -= (
+            ufl.inner(P_(self.avg(self.u_old, self.u, self.alpha_f)), self.test_tensor)
+            * ufl.dx
+        )
+
+        self.k_nominal_proj = F_k_nominal_proj
+
         # self.a = dolfinx.fem.form(ufl.inner(σ(self.u), ufl.grad(self.v)) * ufl.dx)
         # self.L = dolfinx.fem.form(
         #     ufl.dot(self.f, self.v) * ufl.dx
@@ -518,7 +441,7 @@ def build_nullspace(self, V):
     def solve(self, params, dataIO, structure):
         # def σ(v):
         #     """Return an expression for the stress σ given a displacement field"""
-        #     return 2.0 * self.μ * ufl.sym(ufl.grad(v)) + self.λ * ufl.tr(
+        #     return 2.0 * self.lame_mu * ufl.sym(ufl.grad(v)) + self.lame_lambda * ufl.tr(
         #         ufl.sym(ufl.grad(v))
         #     ) * ufl.Identity(len(v))
 
diff --git a/pvade/structure/ModalAnalysis.py b/pvade/structure/ModalAnalysis.py
index 732943c3..62163505 100644
--- a/pvade/structure/ModalAnalysis.py
+++ b/pvade/structure/ModalAnalysis.py
@@ -93,13 +93,17 @@ def __init__(self, domain, structural_analysis, params):
 
         # Define structural properties
         self.E = params.structure.elasticity_modulus  # 1.0e9
-        self.ν = params.structure.poissons_ratio  # 0.3
-        self.μ = self.E / (2.0 * (1.0 + self.ν))
-        self.λ = self.E * self.ν / ((1.0 + self.ν) * (1.0 - 2.0 * self.ν))
+        self.poissons_ratio = params.structure.poissons_ratio  # 0.3
+        self.lame_mu = self.E / (2.0 * (1.0 + self.poissons_ratio))
+        self.lame_lambda = (
+            self.E
+            * self.poissons_ratio
+            / ((1.0 + self.poissons_ratio) * (1.0 - 2.0 * self.poissons_ratio))
+        )
 
         if self.rank == 0:
             print(
-                f"mu = {self.μ} lambda = {self.λ} E = {self.E} nu = {self.ν} density = {self.rho.value}"
+                f"mu = {self.lame_mu} lambda = {self.lame_lambda} E = {self.E} nu = {self.poissons_ratio} density = {self.rho.value}"
             )
 
         # time step
@@ -334,42 +338,44 @@ def build_forms(self, domain, params):
 
         # # Stress tensor
         # def sigma(r):
-        #     return dolfinx.fem.form(2.0*self.μ*ufl.sym(ufl.grad(r)) + self.λ *ufl.tr(ufl.sym(ufl.grad(r)))*ufl.Identity(len(r)))
+        #     return dolfinx.fem.form(2.0*self.lame_mu*ufl.sym(ufl.grad(r)) + self.lame_lambda *ufl.tr(ufl.sym(ufl.grad(r)))*ufl.Identity(len(r)))
 
         # # Mass form
         # def m(u, u_):
         #     return dolfinx.fem.form(self.rho*ufl.inner(u, u_)*ufl.dx)
 
         # # Elastic stiffness form
-        # def k(u, u_):
+        # def k_nominal(u, u_):
         #     return dolfinx.fem.form(ufl.inner(sigma(u), ufl.sym(ufl.grad(u_)))*ufl.dx)
 
         # # Rayleigh damping form
         # def c(u, u_):
-        #     return dolfinx.fem.form(self.eta_m*m(u, u_) + self.eta_k*k(u, u_))
+        #     return dolfinx.fem.form(self.eta_m*m(u, u_) + self.eta_k*k_nominal(u, u_))
 
         # # Work of external forces
         # def Wext(u_):
         #     return ufl.dot(u_, self.f)*self.ds #dss(3)
 
         def sigma(r):
-            return self.λ * ufl.nabla_div(r) * ufl.Identity(
+            return self.lame_lambda * ufl.nabla_div(r) * ufl.Identity(
                 len(r)
-            ) + 2 * self.μ * ufl.sym(ufl.grad(r))
+            ) + 2 * self.lame_mu * ufl.sym(ufl.grad(r))
 
         def m(u, u_):
             return self.rho * ufl.inner(u, u_) * ufl.dx
 
-        def k(u, u_):
+        def k_cauchy(u, u_):
             return ufl.inner(sigma(u), ufl.grad(u_)) * ufl.dx
 
-        def k(u, u_):
+        def k_nominal(u, u_):
             # return ufl.inner(sigma(u),ufl.grad(u_))*ufl.dx
             # updated lagrangian form
             F = ufl.grad(u) + ufl.Identity(len(u))
             E = 0.5 * (F.T * F - ufl.Identity(len(u)))
-            # S = self.λ *ufl.tr(E)*ufl.Identity(len(u))   + 2*self.μ * (E - ufl.tr(E)*ufl.Identity(len(u))  /3.0)
-            S = self.λ * ufl.tr(E) * ufl.Identity(len(u)) + 2 * self.μ * (E)
+            # S = self.lame_lambda *ufl.tr(E)*ufl.Identity(len(u))   + 2*self.lame_mu * (E - ufl.tr(E)*ufl.Identity(len(u))  /3.0)
+            S = self.lame_lambda * ufl.tr(E) * ufl.Identity(
+                len(u)
+            ) + 2 * self.lame_mu * (E)
 
             # return ufl.inner(F * S, ufl.grad(u_)) * ufl.dx
             return ufl.inner(P_(u), ufl.grad(u_)) * ufl.dx
@@ -400,7 +406,7 @@ def S_(u):
             # return lamda * ufl.tr(E) * I + 2.0 * mu * (E - ufl.tr(E) * I / 3.0)
             # TODO: Why does the above form give a better result and where does it come from?
 
-            S_svk = self.λ * ufl.tr(E) * I + 2.0 * self.μ * E
+            S_svk = self.lame_lambda * ufl.tr(E) * I + 2.0 * self.lame_mu * E
             return S_svk
 
         # The first Piola–Kirchhoff stress tensor, P = F*S
@@ -411,15 +417,15 @@ def P_(u):
             return F * S
 
         def c(u, u_):
-            return self.eta_m * m(u, u_) + self.eta_k * k(u, u_)
+            return self.eta_m * m(u, u_) + self.eta_k * k_nominal(u, u_)
 
         # self.uh_exp = dolfinx.fem.Function(self.V,  name="Deformation")
 
-        def σ(v):
-            """Return an expression for the stress σ given a displacement field"""
-            return 2.0 * self.μ * ufl.sym(ufl.grad(v)) + self.λ * ufl.tr(
-                ufl.sym(ufl.grad(v))
-            ) * ufl.Identity(len(v))
+        # def σ(v):
+        #     """Return an expression for the stress σ given a displacement field"""
+        #     return 2.0 * self.lame_mu * ufl.sym(ufl.grad(v)) + self.lame_lambda * ufl.tr(
+        #         ufl.sym(ufl.grad(v))
+        #     ) * ufl.Identity(len(v))
 
         # source term ($f = \rho \omega^2 [x_0, \, x_1]$)
         # self.ω, self.ρ = 300.0, 10.0
@@ -458,7 +464,7 @@ def σ(v):
         self.res = (
             m(self.avg(self.a_old, a_new, self.alpha_m), self.u_)
             + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_)
-            + k(self.avg(self.u_old, self.u, self.alpha_f), self.u_)
+            + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_)
             - self.rho * ufl.inner(self.f, self.u_) * ufl.dx
             - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_)
             * self.ds
@@ -584,7 +590,7 @@ def build_nullspace(self, V):
     def solve(self, params, dataIO):
         # def σ(v):
         #     """Return an expression for the stress σ given a displacement field"""
-        #     return 2.0 * self.μ * ufl.sym(ufl.grad(v)) + self.λ * ufl.tr(
+        #     return 2.0 * self.lame_mu * ufl.sym(ufl.grad(v)) + self.lame_lambda * ufl.tr(
         #         ufl.sym(ufl.grad(v))
         #     ) * ufl.Identity(len(v))
 
diff --git a/pvade/structure/StructureMain.py b/pvade/structure/StructureMain.py
index 68fd33f2..6e55265d 100644
--- a/pvade/structure/StructureMain.py
+++ b/pvade/structure/StructureMain.py
@@ -42,37 +42,9 @@ def __init__(self, domain, params):
         self.rank = domain.rank
         self.num_procs = domain.num_procs
 
-        # domain.structure.msh = dolfinx.mesh.refine(domain.structure.msh,None)
-        # domain.structure.msh = dolfinx.mesh.refine(domain.structure.msh)
-
-        # physics = 'ElasticityAnalysis'
-
-        # domain_creation_module = (
-        #     # f"pvade.geometry.{params.general.geometry_module}.DomainCreation"
-        #     f"pvade.structure.{physics}.Elasticity"
-        # )
-        # try:
-        #     # This is equivalent to "import pvade.geometry.{params.general.geometry_module}.DomainCreation as dcm"
-        #     dcm = import_module(domain_creation_module)
-        # except:
-        #     raise ValueError(f"Could not import {domain_creation_module}")
-
         # init physics
         self.elasticity = Elasticity(domain, self.structural_analysis, params)
 
-        # P1 = ufl.VectorElement("Lagrange", domain.structure.msh.ufl_cell(), 2)
-        # self.V = dolfinx.fem.FunctionSpace(domain.structure.msh, P1)
-
-        # self.W = dolfinx.fem.FunctionSpace(
-        #     domain.structure.msh, ("Discontinuous Lagrange", 0)
-        # )
-
-        # self.first_call_to_solver = True
-
-        # self.num_V_dofs = (
-        #     self.V.dofmap.index_map_bs * self.V.dofmap.index_map.size_global
-        # )
-
         # Store the dimension of the problem for convenience
         self.ndim = domain.structure.msh.topology.dim
         self.facet_dim = self.ndim - 1
@@ -81,13 +53,22 @@ def __init__(self, domain, params):
         num_cells = domain.structure.msh.topology.index_map(self.ndim).size_local
         h = dolfinx.cpp.mesh.h(domain.structure.msh, self.ndim, range(num_cells))
 
-        # This value of hmin is local to the mesh portion owned by the process
-        hmin_local = np.amin(h)
+        # # This value of hmin is local to the mesh portion owned by the process
+        # hmin_local = np.amin(h)
+
+        if len(h) > 0:
+            hmin_local = np.amin(h)
+        else:
+            hmin_local = np.inf
 
-        # collect the minimum hmin from all processes
+        print(hmin_local)
         self.hmin = np.zeros(1)
-        self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN)
-        self.hmin = self.hmin[0]
+        self.hmin = self.comm.allreduce(hmin_local, op=MPI.MIN)
+
+        # # collect the minimum hmin from all processes
+        # self.hmin = np.zeros(1)
+        # self.comm.Allreduce(hmin_local, self.hmin, op=MPI.MIN)
+        # self.hmin = self.hmin[0]
 
         self.num_V_dofs = self.elasticity.num_V_dofs
         if self.rank == 0:
@@ -99,25 +80,19 @@ def __init__(self, domain, params):
             domain.structure.msh, params.structure.rho
         )  # Constant(0.)
 
-        # # Rayleigh damping coefficients
-        # self.eta_m = dolfinx.fem.Constant(domain.structure.msh, 0.0)  # Constant(0.)
-        # self.eta_k = dolfinx.fem.Constant(domain.structure.msh, 0.0)  # Constant(0.)
-
-        # # Generalized-alpha method parameters
-        # self.alpha_m = dolfinx.fem.Constant(domain.structure.msh, 0.2)
-        # self.alpha_f = dolfinx.fem.Constant(domain.structure.msh, 0.4)
-        # self.gamma = 0.5 + self.alpha_f - self.alpha_m
-        # self.beta = (self.gamma + 0.5) ** 2 / 4.0
-
         # Define structural properties
         self.E = params.structure.elasticity_modulus  # 1.0e9
-        self.ν = params.structure.poissons_ratio  # 0.3
-        self.μ = self.E / (2.0 * (1.0 + self.ν))
-        self.λ = self.E * self.ν / ((1.0 + self.ν) * (1.0 - 2.0 * self.ν))
+        self.poissons_ratio = params.structure.poissons_ratio  # 0.3
+        self.lame_mu = self.E / (2.0 * (1.0 + self.poissons_ratio))
+        self.lame_lambda = (
+            self.E
+            * self.poissons_ratio
+            / ((1.0 + self.poissons_ratio) * (1.0 - 2.0 * self.poissons_ratio))
+        )
 
         if self.rank == 0:
             print(
-                f"mu = {self.μ} lambda = {self.λ} E = {self.E} nu = {self.ν} density = {self.rho.value}"
+                f"mu = {self.lame_mu} lambda = {self.lame_lambda} E = {self.E} nu = {self.poissons_ratio} density = {self.rho.value}"
             )
 
         def _north_east_corner(x):
@@ -232,177 +207,6 @@ def build_forms(self, domain, params):
         """
         self.elasticity.build_forms(domain, params, self)
 
-        # # Define trial and test functions for deformation
-        # # self.du = ufl.TrialFunction(self.V)
-        # self.u_ = ufl.TestFunction(self.V)
-
-        # P3 = ufl.TensorElement("Lagrange", domain.structure.msh.ufl_cell(), 2)
-        # self.T = dolfinx.fem.FunctionSpace(domain.structure.msh, P3)
-
-        # self.stress = dolfinx.fem.Function(self.T, name="stress_fluid")
-        # self.stress_old = dolfinx.fem.Function(self.T, name="stress_fluid_old")
-        # self.stress_predicted = dolfinx.fem.Function(
-        #     self.T, name="stress_fluid_predicted"
-        # )
-
-        # self.sigma_vm_h = dolfinx.fem.Function(self.W, name="Stress")
-
-        # # discplacement
-        # self.u = dolfinx.fem.Function(self.V, name="Deformation")
-        # self.u_old = dolfinx.fem.Function(self.V, name="Deformation_old")
-        # self.u_delta = dolfinx.fem.Function(self.V, name="Deformation_change")
-
-        # # velocity
-        # self.v = dolfinx.fem.Function(self.V, name="Velocity")
-        # self.v_old = dolfinx.fem.Function(self.V, name="Velocity_old")
-
-        # # acceleration
-        # self.a = dolfinx.fem.Function(self.V, name="acceleration")
-        # self.a_old = dolfinx.fem.Function(self.V, name="acceleration_old")
-
-        # # dss = ufl.ds(subdomain_data=boundary_subdomains)
-
-        # # # Stress tensor
-        # # def sigma(r):
-        # #     return dolfinx.fem.form(2.0*self.μ*ufl.sym(ufl.grad(r)) + self.λ *ufl.tr(ufl.sym(ufl.grad(r)))*ufl.Identity(len(r)))
-
-        # # # Mass form
-        # # def m(u, u_):
-        # #     return dolfinx.fem.form(self.rho*ufl.inner(u, u_)*ufl.dx)
-
-        # # # Elastic stiffness form
-        # # def k(u, u_):
-        # #     return dolfinx.fem.form(ufl.inner(sigma(u), ufl.sym(ufl.grad(u_)))*ufl.dx)
-
-        # # # Rayleigh damping form
-        # # def c(u, u_):
-        # #     return dolfinx.fem.form(self.eta_m*m(u, u_) + self.eta_k*k(u, u_))
-
-        # # # Work of external forces
-        # # def Wext(u_):
-        # #     return ufl.dot(u_, self.f)*self.ds #dss(3)
-
-        # def sigma(r):
-        #     return self.λ * ufl.nabla_div(r) * ufl.Identity(
-        #         len(r)
-        #     ) + 2 * self.μ * ufl.sym(ufl.grad(r))
-
-        # def m(u, u_):
-        #     return self.rho * ufl.inner(u, u_) * ufl.dx
-
-        # def k(u, u_):
-        #     return ufl.inner(sigma(u), ufl.grad(u_)) * ufl.dx
-
-        # def k(u, u_):
-        #     # return ufl.inner(sigma(u),ufl.grad(u_))*ufl.dx
-        #     # updated lagrangian form
-        #     F = ufl.grad(u) + ufl.Identity(len(u))
-        #     E = 0.5 * (F.T * F - ufl.Identity(len(u)))
-        #     # S = self.λ *ufl.tr(E)*ufl.Identity(len(u))   + 2*self.μ * (E - ufl.tr(E)*ufl.Identity(len(u))  /3.0)
-        #     S = self.λ * ufl.tr(E) * ufl.Identity(len(u)) + 2 * self.μ * (E)
-
-        #     # return ufl.inner(F * S, ufl.grad(u_)) * ufl.dx
-        #     return ufl.inner(P_(u), ufl.grad(u_)) * ufl.dx
-
-        # # The deformation gradient, F = I + dy/dX
-        # def F_(u):
-        #     I = ufl.Identity(len(u))
-        #     return I + ufl.grad(u)
-
-        # # The Cauchy-Green deformation tensor, C = F.T * F
-        # def C_(u):
-        #     F = F_(u)
-        #     return F.T * F
-
-        # # Green–Lagrange strain tensor, E = 0.5*(C - I)
-        # def E_(u):
-        #     I = ufl.Identity(len(u))
-        #     C = C_(u)
-
-        #     return 0.5 * (C - I)
-        #     # return 0.5 * (ufl.grad(u) + ufl.grad(u).T)
-
-        # # The second Piola–Kirchhoff stress, S
-        # def S_(u):
-        #     E = E_(u)
-        #     I = ufl.Identity(len(u))
-
-        #     # return lamda * ufl.tr(E) * I + 2.0 * mu * (E - ufl.tr(E) * I / 3.0)
-        #     # TODO: Why does the above form give a better result and where does it come from?
-
-        #     S_svk = self.λ * ufl.tr(E) * I + 2.0 * self.μ * E
-        #     return S_svk
-
-        # # The first Piola–Kirchhoff stress tensor, P = F*S
-        # def P_(u):
-        #     F = F_(u)
-        #     S = S_(u)
-        #     # return ufl.inv(F) * S
-        #     return F * S
-
-        # def c(u, u_):
-        #     return self.eta_m * m(u, u_) + self.eta_k * k(u, u_)
-
-        # # self.uh_exp = dolfinx.fem.Function(self.V,  name="Deformation")
-
-        # def σ(v):
-        #     """Return an expression for the stress σ given a displacement field"""
-        #     return 2.0 * self.μ * ufl.sym(ufl.grad(v)) + self.λ * ufl.tr(
-        #         ufl.sym(ufl.grad(v))
-        #     ) * ufl.Identity(len(v))
-
-        # # source term ($f = \rho \omega^2 [x_0, \, x_1]$)
-        # # self.ω, self.ρ = 300.0, 10.0
-        # # x = ufl.SpatialCoordinate(domain.structure.msh)
-        # # self.f = ufl.as_vector((0*self.ρ * self.ω**2 * x[0], self.ρ * self.ω**2 * x[1], 0.0))
-        # # self.f_structure = dolfinx.fem.Constant(
-        # #     domain.structure.msh,
-        # #     (PETSc.ScalarType(0), PETSc.ScalarType(0), PETSc.ScalarType(0)),
-        # # )
-        # # self.f = ufl.as_vector((0*self.ρ * self.ω**2 * x[0], self.ρ * self.ω**2 * x[1], 0.0))
-        # # self.T = dolfinx.fem.Constant(domain.structure.msh, PETSc.ScalarType((0, 1.e-3, 0)))
-        # # self.f = dolfinx.fem.Constant(domain.structure.msh, PETSc.ScalarType((0,100,100)))
-        # self.f = dolfinx.fem.Constant(
-        #     domain.structure.msh,
-        #     PETSc.ScalarType(
-        #         (
-        #             params.structure.body_force_x,
-        #             params.structure.body_force_y,
-        #             params.structure.body_force_z,
-        #         )
-        #     ),
-        # )
-        # self.ds = ufl.Measure("ds", domain=domain.structure.msh)
-        # n = ufl.FacetNormal(domain.structure.msh)
-
-        # # Residual
-        # a_new = self.update_a(
-        #     self.u, self.u_old, self.v_old, self.a_old, self.dt_st, self.beta, ufl=True
-        # )
-        # v_new = self.update_v(
-        #     a_new, self.u_old, self.v_old, self.a_old, self.dt_st, self.gamma, ufl=True
-        # )
-
-        # F = ufl.grad(self.u) + ufl.Identity(len(self.u))
-        # J = ufl.det(F)
-        # self.res = (
-        #     m(self.avg(self.a_old, a_new, self.alpha_m), self.u_)
-        #     + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_)
-        #     + k(self.avg(self.u_old, self.u, self.alpha_f), self.u_)
-        #     - self.rho * ufl.inner(self.f, self.u_) * ufl.dx
-        #     - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_)
-        #     * self.ds
-        # )  # - Wext(self.u)
-
-        # # self.a = dolfinx.fem.form(ufl.lhs(res))
-        # # self.L = dolfinx.fem.form(ufl.rhs(res))
-
-        # # self.a = dolfinx.fem.form(ufl.inner(σ(self.u), ufl.grad(self.v)) * ufl.dx)
-        # # self.L = dolfinx.fem.form(
-        # #     ufl.dot(self.f, self.v) * ufl.dx
-        # #     + ufl.dot(ufl.dot(self.stress, n), self.v) * self.ds
-        # # )
-
     def _assemble_system(self, params):
         """Pre-assemble all LHS matrices and RHS vectors
 
@@ -514,88 +318,3 @@ def build_nullspace(self, V):
     def solve(self, params, dataIO):
 
         self.elasticity.solve(params, dataIO, self)
-        # # def σ(v):
-        # #     """Return an expression for the stress σ given a displacement field"""
-        # #     return 2.0 * self.μ * ufl.sym(ufl.grad(v)) + self.λ * ufl.tr(
-        # #         ufl.sym(ufl.grad(v))
-        # #     ) * ufl.Identity(len(v))
-
-        # if self.first_call_to_solver:
-        #     if self.rank == 0:
-        #         print("Starting Strutural Solution")
-
-        #     self._assemble_system(params)
-
-        # num_its, converged = self.solver.solve(self.u)  # solve the current time step
-        # assert converged
-        # self.u.x.scatter_forward()
-
-        # # Calculate the change in the displacement (new - old) this is what moves the mesh
-        # self.u_delta.vector.array[:] = (
-        #     self.u.vector.array[:] - self.u_old.vector.array[:]
-        # )
-        # self.u_delta.x.scatter_forward()
-
-        # # Update old fields with new quantities
-        # self.update_fields(
-        #     self.u,
-        #     self.u_old,
-        #     self.v_old,
-        #     self.a_old,
-        #     self.dt_st,
-        #     self.beta,
-        #     self.gamma,
-        # )
-
-        # # sigma_dev = σ(self.u) - (1 / 3) * ufl.tr(σ(self.u)) * ufl.Identity(len(self.u))
-        # # sigma_vm = ufl.sqrt((3 / 2) * ufl.inner(sigma_dev, sigma_dev))
-
-        # # sigma_vm_expr = dolfinx.fem.Expression(
-        # #     sigma_vm, self.W.element.interpolation_points()
-        # # )
-        # # self.sigma_vm_h.interpolate(sigma_vm_expr)
-
-        # # self.unorm = self.u.x.norm()
-
-        # try:
-        #     idx = self.north_east_corner_dofs[0]
-        #     nw_corner_accel = self.u.vector.array[3 * idx : 3 * idx + 3].astype(
-        #         np.float64
-        #     )
-        #     print(nw_corner_accel)
-        # except:
-        #     nw_corner_accel = np.zeros(3, dtype=np.float64)
-
-        # nw_corner_accel_global = np.zeros((self.num_procs, 3), dtype=np.float64)
-
-        # self.comm.Gather(nw_corner_accel, nw_corner_accel_global, root=0)
-
-        # # print(f"Acceleration at North West corner = {nw_corner_accel}")
-
-        # if self.rank == 0:
-        #     norm2 = np.sum(nw_corner_accel_global**2, axis=1)
-        #     max_norm2_idx = np.argmax(norm2)
-        #     np_accel = nw_corner_accel_global[max_norm2_idx, :]
-
-        #     accel_pos_filename = os.path.join(
-        #         params.general.output_dir_sol, "accel_pos.csv"
-        #     )
-
-        #     if self.first_call_to_solver:
-
-        #         with open(accel_pos_filename, "w") as fp:
-        #             fp.write("#x-pos,y-pos,z-pos\n")
-        #             if self.ndim == 3:
-        #                 fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n")
-        #             elif self.ndim == 2:
-        #                 fp.write(f"{np_accel[0]},{np_accel[1]}\n")
-
-        #     else:
-        #         with open(accel_pos_filename, "a") as fp:
-        #             if self.ndim == 3:
-        #                 fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n")
-        #             elif self.ndim == 2:
-        #                 fp.write(f"{np_accel[0]},{np_accel[1]}\n")
-
-        # if self.first_call_to_solver:
-        #     self.first_call_to_solver = False
diff --git a/pvade/tests/input/yaml/sim_params.yaml b/pvade/tests/input/yaml/sim_params.yaml
index fd216a06..7b6e4a6b 100644
--- a/pvade/tests/input/yaml/sim_params.yaml
+++ b/pvade/tests/input/yaml/sim_params.yaml
@@ -15,6 +15,7 @@ domain:
   l_char: 20
 pv_array:
   stream_rows: 7
+  span_rows: 1
   elevation: 1.5
   stream_spacing: 7.0
   panel_chord: 2.0
diff --git a/pvade/tests/test_input_files.py b/pvade/tests/test_input_files.py
index 5354530e..01c13c86 100644
--- a/pvade/tests/test_input_files.py
+++ b/pvade/tests/test_input_files.py
@@ -19,7 +19,7 @@ def launch_sim(test_file):
     tf = dt * 10  # ten timesteps
 
     command = (
-        f"mpirun -n 2 python "
+        f"mpirun -n 8 python "
         + rootdir
         + "/pvade_main.py --input_file "
         + test_file["path_to_file"]
diff --git a/pvade/tests/test_solve.py b/pvade/tests/test_solve.py
index 41a8eb1f..8734207d 100644
--- a/pvade/tests/test_solve.py
+++ b/pvade/tests/test_solve.py
@@ -22,7 +22,7 @@
 
 solve_iter = 10
 
-rtol = 1.0e-4
+rtol = 1.0e-1
 
 
 @pytest.mark.unit
diff --git a/pvade_main.py b/pvade_main.py
index 1c2e7689..8bc322aa 100644
--- a/pvade_main.py
+++ b/pvade_main.py
@@ -14,6 +14,7 @@
 
 from pvade.structure.StructureMain import Structure
 import os
+from mpi4py import MPI
 
 
 def main(input_file=None):
@@ -25,6 +26,7 @@ def main(input_file=None):
     # Load the parameters object specified by the input file
     params = SimParams(input_file)
 
+    # Save the parameters object to a logfile
     logfile_name = os.path.join(params.general.output_dir, "logfile.log")
     start_print_and_log(params.rank, logfile_name)
 
@@ -34,64 +36,45 @@ def main(input_file=None):
 
     # Initialize the domain and construct the initial mesh
     domain = FSIDomain(params)
-
     if params.general.input_mesh_dir is not None:
         domain.read_mesh_files(params.general.input_mesh_dir, params)
     else:
         domain.build(params)
-        # domain.write_mesh_files(params)
 
+    # If we only want to create the mesh, we can stop here
     if params.general.mesh_only:
         list_timings(params.comm, [TimingType.wall])
         structure, flow = [], []
         return params, structure, flow
 
-    # Check to ensure mesh node matching for periodic simulations
-    # if domain.periodic_simulation:
-    # domain.check_mesh_periodicity(params)
-    # sys.exit()
-
-    # if params.general.debug_flag:
-    #     print('before Flow init')
-
-    flow = Flow(domain, params)
-    structure = Structure(domain, params)
-
     if fluid_analysis:
         flow = Flow(domain, params)
         # # # Specify the boundary conditions
         flow.build_boundary_conditions(domain, params)
         # # # Build the fluid forms
         flow.build_forms(domain, params)
+    else:
+        flow = None
 
     if structural_analysis:
+        structure = Structure(domain, params)
         structure.build_boundary_conditions(domain, params)
         # # # Build the fluid forms
         structure.build_forms(domain, params)
+    else:
+        structure = None
 
     if structural_analysis and fluid_analysis:
-        # pass
         domain.move_mesh(structure, params)
 
     dataIO = DataStream(domain, flow, structure, params)
     FSI_inter = FSI(domain, flow, structure, params)
-    # if domain.rank == 0:
-    #     progress = tqdm.autonotebook.tqdm(
-    #         desc="Solving PDE", total=params.solver.t_steps
-    #     )
 
     solve_structure_interval_n = int(params.structure.dt / params.solver.dt)
 
     for k in range(params.solver.t_steps):
         current_time = (k + 1) * params.solver.dt
 
-        # if domain.rank == 0:
-        #     progress.update(1)
-        #     print()
-        # Solve the fluid problem at each timestep
-
-        # print("time step is : ", current_time)
-        # print("reaminder from modulo ",current_time % params.structure.dt )
         if (
             structural_analysis
             and (k + 1) % solve_structure_interval_n == 0
@@ -139,12 +122,13 @@ def main(input_file=None):
                         print(f"| f_x (drag) = {fx:.4f}")
                         print(f"| f_y (lift) = {fy:.4f}")
                     else:
-                        # still print info, but just for first row
-                        fx = flow.integrated_force_x[0]
-                        fy = flow.integrated_force_y[0]
+                        if structure is not None:
+                            # still print info, but just for first row
+                            fx = flow.integrated_force_x[0]
+                            fy = flow.integrated_force_y[0]
 
-                        print(f"| f_x (drag) of 1st row = {fx:.4f}")
-                        print(f"| f_y (lift) of 1st row = {fy:.4f}")
+                            print(f"| f_x (drag) of 1st row = {fx:.4f}")
+                            print(f"| f_y (lift) of 1st row = {fy:.4f}")
                     if thermal_analysis:
                         print(f"| T = {flow.theta_max:.4f}")
 
@@ -156,15 +140,34 @@ def main(input_file=None):
                 and current_time > params.fluid.warm_up_time
             ):
 
-                local_def_max = np.amax(
-                    np.sum(
-                        structure.elasticity.u.vector.array.reshape(-1, domain.ndim)
-                        ** 2,
-                        axis=1,
-                    )
+                vec = structure.elasticity.u.vector
+                u_local = vec.array  # This gives a NumPy array on local rank
+
+                if u_local.size > 0:
+                    u_reshaped = u_local.reshape(-1, domain.ndim)
+                    local_def_max = np.amax(np.sum(u_reshaped**2, axis=1))
+                else:
+                    local_def_max = -np.inf  # So it doesn't interfere in max
+                print(
+                    f"Rank {domain.comm.rank}: u_vec.size = {u_local.size}, local_def_max = {local_def_max}"
                 )
-                global_def_max_list = np.zeros(params.num_procs, dtype=np.float64)
-                params.comm.Gather(local_def_max, global_def_max_list, root=0)
+
+                # Wrap scalar in NumPy array (required for Gather)
+                sendbuf = np.array([local_def_max], dtype=np.float64)
+
+                # Preallocate only on root
+                if params.comm.rank == 0:
+                    global_def_max_list = np.empty(params.num_procs, dtype=np.float64)
+                else:
+                    global_def_max_list = None
+
+                # Perform the gather
+                params.comm.Gather(sendbuf, global_def_max_list, root=0)
+
+                # Handle on root
+                if params.comm.rank == 0:
+                    print("Per-rank def max values:", global_def_max_list)
+                    print("Global def max:", np.max(global_def_max_list))
 
                 if domain.rank == 0:
                     # print("Structural time is : ", current_time)
@@ -198,4 +201,4 @@ def main(input_file=None):
             sys.stdout = sys.__stdout__
 
         if not params.general.mesh_only:
-            write_metrics(flow, structure.elasticity, prof_filename=prof_txt_filename)
+            write_metrics(flow, prof_filename=prof_txt_filename)
diff --git a/test_all_inputs.py b/test_all_inputs.py
new file mode 100644
index 00000000..3e16def7
--- /dev/null
+++ b/test_all_inputs.py
@@ -0,0 +1,46 @@
+import subprocess
+from pathlib import Path
+import pytest
+
+
+@pytest.mark.parametrize("mesh_only", [True, False])
+@pytest.mark.parametrize("nprocs", [1])
+def test_pvade_run(input_file, mesh_only, nprocs):
+    cmd = [
+        "mpirun",
+        "-n",
+        str(nprocs),
+        "python",
+        "pvade_main.py",
+        "--input",
+        str(input_file.resolve()),
+        "--solver.dt",
+        "0.001",
+        "--solver.t_final",
+        "0.005",
+        "--general.mesh_only",
+        str(mesh_only).lower(),
+    ]
+
+    # Add special argument for duramat_case_study.yaml
+    if input_file.name == "duramat_case_study.yaml":
+        cmd += ["--domain.l_char", "4"]
+
+    print(
+        f"\n=== Running: {input_file.name} | mesh_only={mesh_only} | nprocs={nprocs} ==="
+    )
+    print("Command:", " ".join(cmd))
+
+    proc = subprocess.Popen(
+        cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True
+    )
+
+    for line in proc.stdout:
+        print(line, end="")
+
+    proc.wait()
+
+    assert proc.returncode == 0, (
+        f"\nFAILED: {input_file.name} | mesh_only={mesh_only} | nprocs={nprocs}\n"
+        f"Exit code: {proc.returncode}"
+    )