From 7a473ae063cd52b5f2d1324ce16086dd0e6c9302 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Mon, 8 Sep 2025 16:38:13 -0600
Subject: [PATCH 01/26] adding new options for torque tube radius, sep.
 distance, and modules_per_span

---
 pvade/IO/input_schema.yaml | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index 1c99362..335dacb 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -136,7 +136,7 @@ properties:
         minimum: 0.0
         maximum: 10.0
         type: "number"
-        description: "The vertical distance between the center of the panel and the ground."
+        description: "The vertical distance between the center of the torque tube and the ground."
         units: "meter"
       stream_spacing:
         default: 7.0
@@ -198,6 +198,23 @@ properties:
         be unrolled such that tracking_angle_0 is applied to the panel at (x_min, y_min), tracking_angle_1 is applied to the panel at
         (x_min + x_spacing, y_min), etc., i.e., x-major direction. If argument is a single value, that tracking angle will be applied to the every panel."
         units: "degree"
+      torque_tube_separation:
+        default: 0.4
+        minimum: 0.0
+        maximum: 1.0
+        type: "number"
+        description: "The distance between the bottom of the panel structure and the centerline of the torque tube, in meters."
+      torque_tube_radius:
+        default: 0.1
+        minimum: 0.0
+        maximum: 1.0
+        type: "number"
+        description: "The radius of the torque tube, in meters."
+      modules_per_span:
+        default: 1
+        minimum: 1
+        type: "integer"
+        description: "The number of modules/panels per every panel_span distance. This can be thought of as the number of modules/panels per each table, counted in the spanwise direction only (i.e., 2-in-portait systems should report only the number in the spanwise direction, vs 2x that value)."
   solver:
     additionalProperties: false
     type: "object"

From e5e8bb7220d7334dce4086f2f55527c917b81114 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Mon, 8 Sep 2025 16:38:43 -0600
Subject: [PATCH 02/26] making panels differently, including torque tube
 standoff and discrete modules per span

---
 pvade/geometry/panels3d/DomainCreation.py | 393 +++++++++++++---------
 1 file changed, 225 insertions(+), 168 deletions(-)

diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py
index 5ead1e4..920eda6 100644
--- a/pvade/geometry/panels3d/DomainCreation.py
+++ b/pvade/geometry/panels3d/DomainCreation.py
@@ -152,7 +152,11 @@ def Rz(theta):
         panel_tag_list = []
         panel_ct = 0
 
-        prev_surf_tag = []
+        module_distances = np.linspace(
+            -half_span, half_span, params.pv_array.modules_per_span + 1
+        )
+
+        transformed_com = {}
 
         for panel_id_y, yy in enumerate(y_centers):
             for panel_id_x, xx in enumerate(x_centers):
@@ -169,29 +173,48 @@ def Rz(theta):
                 else:
                     tracker_angle_rad = np.radians(params.pv_array.tracker_angle)
 
-                # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0)
-                panel_id = self.gmsh_model.occ.addBox(
-                    -half_chord,
-                    -half_span,
-                    -half_thickness,
-                    params.pv_array.panel_chord,
-                    params.pv_array.panel_span,
-                    params.pv_array.panel_thickness,
-                )
+                this_panel_tag_list = []
+                this_panel_transformed_com = {}
 
-                panel_tag = (self.ndim, panel_id)
-                panel_tag_list.append(panel_tag)
+                for module_id in range(params.pv_array.modules_per_span):
+
+                    module_span = (
+                        module_distances[module_id + 1] - module_distances[module_id]
+                    )
+
+                    # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0)
+                    this_module = self.gmsh_model.occ.addBox(
+                        -half_chord,
+                        module_distances[module_id],
+                        params.pv_array.torque_tube_separation,
+                        params.pv_array.panel_chord,
+                        module_span,
+                        params.pv_array.panel_thickness,
+                    )
+
+                    this_standoff = self.gmsh_model.occ.addBox(
+                        -params.pv_array.torque_tube_radius,
+                        module_distances[module_id],
+                        0.0,
+                        2.0 * params.pv_array.torque_tube_radius,
+                        module_span,
+                        params.pv_array.torque_tube_separation,
+                    )
+
+                    panel_tag_list.append((self.ndim, this_module))
+                    panel_tag_list.append((self.ndim, this_standoff))
+
+                    this_panel_tag_list.append((self.ndim, this_module))
+                    this_panel_tag_list.append((self.ndim, this_standoff))
 
                 numpy_pt_list = []
                 embedded_lines_tag_list = []
 
                 # Add a bisecting line to the bottom of the panel in the spanwise direction
-                pt_1 = self.gmsh_model.occ.addPoint(0, -half_span, -half_thickness)
-                pt_2 = self.gmsh_model.occ.addPoint(0, half_span, -half_thickness)
+                pt_1 = self.gmsh_model.occ.addPoint(0, -half_span, 0.0)
+                pt_2 = self.gmsh_model.occ.addPoint(0, half_span, 0.0)
 
-                numpy_pt_list.append(
-                    [0, -half_span, -half_thickness, 0, half_span, -half_thickness]
-                )
+                numpy_pt_list.append([0, -half_span, 0.0, 0, half_span, 0.0])
 
                 torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
                 torque_tube_tag = (1, torque_tube_id)
@@ -224,150 +247,149 @@ def Rz(theta):
                         fixation_pts_list = params.pv_array.span_fixation_pts
 
                     for fp in fixation_pts_list:
-                        pt_1 = self.gmsh_model.occ.addPoint(
-                            -half_chord, -half_span + fp, -half_thickness
-                        )
-                        pt_2 = self.gmsh_model.occ.addPoint(
-                            half_chord, -half_span + fp, -half_thickness
-                        )
-
                         # FIXME: don't add the fixation points into the numpy tagging for now
                         numpy_pt_list.append(
                             [
-                                -half_chord,
+                                -params.pv_array.torque_tube_radius,
                                 -half_span + fp,
-                                -half_thickness,
-                                half_chord,
+                                0.0,
+                                params.pv_array.torque_tube_radius,
                                 -half_span + fp,
-                                -half_thickness,
+                                0.0,
                             ]
                         )
 
-                        fixed_pt_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
-                        fixed_pt_tag = (1, fixed_pt_id)
-
-                        embedded_lines_tag_list.append(fixed_pt_tag)
-
-                # Store the result of fragmentation, it holds all the small surfaces we need to tag
-                panel_frags = self.gmsh_model.occ.fragment(
-                    [panel_tag], embedded_lines_tag_list
+                # Fragment the lines into the surfaces (equivalent to embedding these lines)
+                self.gmsh_model.occ.fragment(
+                    this_panel_tag_list, embedded_lines_tag_list
                 )
 
-                # extract just the first entry, and remove the 3d entry in position 0
-                panel_surfs = panel_frags[0]
-                panel_surfs.pop(0)
-                panel_surfs = [k[1] for k in panel_surfs]
-
-                # TODO: USE THESE UNAMBIGUOUS NAMES IN A FUTURE REFACTOR
-                # self._add_to_domain_markers(f"x_min_{panel_ct:.0f}", [panel_surfs[0]], "facet")
-                # self._add_to_domain_markers(f"x_max_{panel_ct:.0f}", [panel_surfs[1]], "facet")
-                # self._add_to_domain_markers(f"y_min_{panel_ct:.0f}", [panel_surfs[2]], "facet")
-                # self._add_to_domain_markers(f"y_max_{panel_ct:.0f}", [panel_surfs[3]], "facet")
-                # self._add_to_domain_markers(f"z_min_{panel_ct:.0f}", panel_surfs[4:-1], "facet")
-                # self._add_to_domain_markers(f"z_max_{panel_ct:.0f}", [panel_surfs[-1]], "facet")
-
-                # Translate the panel by (x_center, y_center, elev)
-                self.gmsh_model.occ.translate(
-                    [panel_tag],
-                    xx,
-                    yy,
-                    params.pv_array.elevation,
-                )
+                for panel_tag in this_panel_tag_list:
+                    self.gmsh_model.occ.synchronize()
 
-                # Get the bounding box for this panel
-                panel_bounding_box = self.gmsh_model.occ.getBoundingBox(
-                    panel_tag[0], panel_tag[1]
-                )
+                    # Get the list of 2D surfaces (surfaces) that make up this panel
+                    surf_tags_for_this_panel = self.gmsh_model.getBoundary(
+                        [panel_tag], oriented=False
+                    )
 
-                # Store the x_min, y_min, x_max, y_max points
-                # corresponding to the un-rotated (tracking_angle = 0) configuration
-                x_min_panel = panel_bounding_box[0]
-                y_min_panel = panel_bounding_box[1]
-                z_min_panel = panel_bounding_box[2]
-                x_max_panel = panel_bounding_box[3]
-                y_max_panel = panel_bounding_box[4]
-                z_max_panel = panel_bounding_box[5]
+                    for surf_tag in surf_tags_for_this_panel:
+                        surf_dim = surf_tag[0]
+                        surf_id = surf_tag[1]
+                        com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id)
+
+                        # sturctures tagging
+                        if np.isclose(com[0], -half_chord) or np.isclose(
+                            com[0], -params.pv_array.torque_tube_radius
+                        ):
+                            key = f"left_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+
+                        elif np.isclose(com[0], half_chord) or np.isclose(
+                            com[0], params.pv_array.torque_tube_radius
+                        ):
+                            key = f"right_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+
+                        elif np.isclose(com[1], -half_span):
+                            key = f"front_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+
+                        elif np.isclose(com[1], half_span):
+                            key = f"back_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+
+                        elif np.isclose(
+                            com[2], params.pv_array.torque_tube_separation
+                        ) and not np.isclose(com[0], 0.0):
+                            key = f"bottom_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+
+                        elif np.isclose(com[2], 0.0):
+                            key = f"torque_tube_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+
+                        elif np.isclose(
+                            com[2],
+                            params.pv_array.torque_tube_separation
+                            + params.pv_array.panel_thickness,
+                        ):
+                            key = f"top_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
+                        else:
+                            self._add_to_domain_markers(
+                                f"trash_{panel_ct:.0f}", [surf_id], "facet"
+                            )
 
-                self.gmsh_model.occ.synchronize()
+                for key, val in this_panel_transformed_com.items():
+                    for row_num, com in enumerate(val):
+                        com_array = np.array(com)
 
-                # Get the list of 1D surfaces (edges) that make up this panel
-                surf_tags_for_this_panel = self.gmsh_model.getBoundary(
-                    [panel_tag], oriented=False
-                )
+                        com_array = np.dot(com_array, Ry(tracker_angle_rad).T)
 
-                for surf_tag in surf_tags_for_this_panel:
-                    surf_dim = surf_tag[0]
-                    surf_id = surf_tag[1]
-                    com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id)
-                    # sturctures tagging
-                    if np.isclose(com[0], x_min_panel):
-                        self._add_to_domain_markers(
-                            f"left_{panel_ct:.0f}", [surf_id], "facet"
-                        )
+                        com_array[0] += xx
+                        com_array[1] += yy
+                        com_array[2] += params.pv_array.elevation
 
-                    elif np.isclose(com[0], x_max_panel):
-                        self._add_to_domain_markers(
-                            f"right_{panel_ct:.0f}", [surf_id], "facet"
-                        )
+                        com_array[0] -= x_center_of_mass
+                        com_array[1] -= y_center_of_mass
 
-                    elif np.isclose(com[1], y_min_panel):
-                        self._add_to_domain_markers(
-                            f"front_{panel_ct:.0f}", [surf_id], "facet"
-                        )
+                        com_array = np.dot(com_array, Rz(array_rotation_rad).T)
 
-                    elif np.isclose(com[1], y_max_panel):
-                        self._add_to_domain_markers(
-                            f"back_{panel_ct:.0f}", [surf_id], "facet"
-                        )
+                        com_array[0] += x_center_of_mass
+                        com_array[1] += y_center_of_mass
 
-                    elif np.isclose(com[2], z_min_panel):
-                        self._add_to_domain_markers(
-                            f"bottom_{panel_ct:.0f}", [surf_id], "facet"
-                        )
-
-                    elif np.isclose(com[2], z_max_panel):
-                        self._add_to_domain_markers(
-                            f"top_{panel_ct:.0f}", [surf_id], "facet"
-                        )
+                        if key in transformed_com:
+                            transformed_com[key].append(com_array)
+                        else:
+                            transformed_com[key] = [com_array]
 
                 panel_ct += 1
 
                 # Rotate the panel by its tracking angle along the y-axis
-                # (currently centered at (xx, yy, params.pv_array.elevation))
+                # (currently centered at (0.0, 0.0, 0.0))
                 self.gmsh_model.occ.rotate(
-                    [panel_tag],
-                    xx,
-                    yy,
-                    params.pv_array.elevation,
+                    this_panel_tag_list,
+                    0.0,
+                    0.0,
+                    0.0,
                     0,
                     1,
                     0,
                     tracker_angle_rad,
                 )
 
-                # Note that the numpy pt panel array has NOT been shifted, i.e.,
-                # it still represents a panel centered at (0, 0, 0), so the rotation
-                # can be applied directly about the point (0, 0, 0) without any
-                # shifting like that which needs to be done above
-                numpy_pt_panel_array = np.array(numpy_pt_list)
-                numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim))
-
-                numpy_pt_panel_array = np.dot(
-                    numpy_pt_panel_array, Ry(tracker_angle_rad).T
+                # Translate the panel by (x_center, y_center, elev)
+                self.gmsh_model.occ.translate(
+                    this_panel_tag_list,
+                    xx,
+                    yy,
+                    params.pv_array.elevation,
                 )
 
-                # if not hasattr(self, "numpy_pt_array"):
-                #     numpy_pt_array = np.array(numpy_pt_list)
-                # else:
-                #     numpy_pt_array = np.vcat(numpy_pt_array, np.array(numpy_pt_list))
-
-                numpy_pt_panel_array[:, 0] += xx
-                numpy_pt_panel_array[:, 1] += yy
-                numpy_pt_panel_array[:, 2] += params.pv_array.elevation
-
                 # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0)
                 self.gmsh_model.occ.rotate(
-                    [panel_tag],
+                    this_panel_tag_list,
                     x_center_of_mass,
                     y_center_of_mass,
                     0,
@@ -377,6 +399,22 @@ def Rz(theta):
                     array_rotation_rad,
                 )
 
+                # Now, apply the same transformations to the numpy representation
+                # Rotate the panel by its tracking angle along the y-axis
+                # (currently centered at (0.0, 0.0, 0.0))
+                numpy_pt_panel_array = np.array(numpy_pt_list)
+                numpy_pt_panel_array = np.reshape(numpy_pt_panel_array, (-1, self.ndim))
+
+                numpy_pt_panel_array = np.dot(
+                    numpy_pt_panel_array, Ry(tracker_angle_rad).T
+                )
+
+                # Translate the panel by (x_center, y_center, elev)
+                numpy_pt_panel_array[:, 0] += xx
+                numpy_pt_panel_array[:, 1] += yy
+                numpy_pt_panel_array[:, 2] += params.pv_array.elevation
+
+                # Rotate the panel about the center of the full array as a proxy for changing wind direction (x_center, y_center, 0)
                 numpy_pt_panel_array[:, 0] -= x_center_of_mass
                 numpy_pt_panel_array[:, 1] -= y_center_of_mass
 
@@ -394,25 +432,26 @@ def Rz(theta):
                 else:
                     self.numpy_pt_total_array = np.copy(numpy_pt_panel_array)
 
-                # Check that this panel still exists in the confines of the domain
-                bbox = self.gmsh_model.occ.get_bounding_box(panel_tag[0], panel_tag[1])
-
-                if bbox[0] < params.domain.x_min:
-                    raise ValueError(
-                        f"Panel with location (x, y) = ({xx}, {yy}) extends past x_min wall."
-                    )
-                if bbox[1] < params.domain.y_min:
-                    raise ValueError(
-                        f"Panel with location (x, y) = ({xx}, {yy}) extends past y_min wall."
-                    )
-                if bbox[3] > params.domain.x_max:
-                    raise ValueError(
-                        f"Panel with location (x, y) = ({xx}, {yy}) extends past x_max wall."
-                    )
-                if bbox[4] > params.domain.y_max:
-                    raise ValueError(
-                        f"Panel with location (x, y) = ({xx}, {yy}) extends past y_max wall."
-                    )
+                # # Check that this panel still exists in the confines of the domain
+                # bbox = self.gmsh_model.occ.get_bounding_box(panel_tag_list)
+
+                # if bbox[0] < params.domain.x_min:
+                #     raise ValueError(
+                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past x_min wall."
+                #     )
+                # if bbox[1] < params.domain.y_min:
+                #     raise ValueError(
+                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past y_min wall."
+                #     )
+                # if bbox[3] > params.domain.x_max:
+                #     raise ValueError(
+                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past x_max wall."
+                #     )
+                # if bbox[4] > params.domain.y_max:
+                #     raise ValueError(
+                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past y_max wall."
+                #     )
+            # self.gmsh_model.occ.synchronize()
 
         # Fragment all panels from the overall domain
         self.gmsh_model.occ.fragment(domain_tag_list, panel_tag_list)
@@ -423,49 +462,69 @@ def Rz(theta):
             self.numpy_pt_total_array, (-1, int(2 * self.ndim))
         )
 
-        # import matplotlib.pyplot as plt
-        # for k in self.numpy_pt_total_array:
-        #     plt.plot([k[0], k[3]], [k[1], k[4]])
-        # plt.show()
-
-        # it is not necessary to loop over all surfaces, since the panel
-        # surfaces have been tagged already, but this ensures any ordering change
-        # doesn't cause problems in the future
+        # Loop over all the finalized surfaces after fragmentation and tag everything
         all_surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim - 1)
 
         for surf_tag in all_surf_tag_list:
             surf_id = surf_tag[1]
             com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id)
 
+            tagged_this_surface = False
+
             # sturctures tagging
             if np.isclose(com[0], params.domain.x_min):
                 self._add_to_domain_markers("x_min", [surf_id], "facet")
+                tagged_this_surface = True
 
             elif np.allclose(com[0], params.domain.x_max):
                 self._add_to_domain_markers("x_max", [surf_id], "facet")
+                tagged_this_surface = True
 
             elif np.allclose(com[1], params.domain.y_min):
                 self._add_to_domain_markers("y_min", [surf_id], "facet")
+                tagged_this_surface = True
 
             elif np.allclose(com[1], params.domain.y_max):
                 self._add_to_domain_markers("y_max", [surf_id], "facet")
+                tagged_this_surface = True
 
             elif np.allclose(com[2], params.domain.z_min):
                 self._add_to_domain_markers("z_min", [surf_id], "facet")
+                tagged_this_surface = True
 
             elif np.allclose(com[2], params.domain.z_max):
                 self._add_to_domain_markers("z_max", [surf_id], "facet")
+                tagged_this_surface = True
+
+            else:
+                for key, val in transformed_com.items():
+                    for target_com in val:
+                        # print(target_com)
+                        if np.allclose(np.array(com), target_com):
+                            self._add_to_domain_markers(key, [surf_id], "facet")
+                            tagged_this_surface = True
+
+            # if not tagged_this_surface:
+            #     print(f"Warning: Surface {surf_tag} has not been added to domain markers")
+
+        # print(self.domain_markers)
 
         # Volumes are the entities with dimension equal to the mesh dimension
         vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim)
         structure_vol_list = []
         fluid_vol_list = []
 
+        num_solids = len(vol_tag_list)
+
         for vol_tag in vol_tag_list:
             vol_id = vol_tag[1]
 
-            if vol_id <= params.pv_array.stream_rows * params.pv_array.span_rows:
+            if vol_id <= num_solids - 1:
                 # Solid Cell
+                # Since all panel/table structures were created after the box domain
+                # the final fragmentation removes the original vol_tag of the box and appends
+                # it to the end, thus, everything with id < num_solids is structure, and
+                # vol_tag = num_solids (the last-added, fragmented box) is fluid
                 structure_vol_list.append(vol_id)
             else:
                 # Fluid Cell
@@ -1144,23 +1203,21 @@ def set_length_scales(self, params, domain_markers):
         internal_surface_tags = []
 
         for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows):
-            internal_surface_tags.append(
-                domain_markers[f"bottom_{panel_id}"]["gmsh_tags"][0]
+            internal_surface_tags.extend(
+                domain_markers[f"bottom_{panel_id}"]["gmsh_tags"]
             )
-            internal_surface_tags.append(
-                domain_markers[f"top_{panel_id}"]["gmsh_tags"][0]
+            internal_surface_tags.extend(domain_markers[f"top_{panel_id}"]["gmsh_tags"])
+            internal_surface_tags.extend(
+                domain_markers[f"left_{panel_id}"]["gmsh_tags"]
             )
-            internal_surface_tags.append(
-                domain_markers[f"left_{panel_id}"]["gmsh_tags"][0]
-            )
-            internal_surface_tags.append(
-                domain_markers[f"right_{panel_id}"]["gmsh_tags"][0]
+            internal_surface_tags.extend(
+                domain_markers[f"right_{panel_id}"]["gmsh_tags"]
             )
-            internal_surface_tags.append(
-                domain_markers[f"front_{panel_id}"]["gmsh_tags"][0]
+            internal_surface_tags.extend(
+                domain_markers[f"front_{panel_id}"]["gmsh_tags"]
             )
-            internal_surface_tags.append(
-                domain_markers[f"back_{panel_id}"]["gmsh_tags"][0]
+            internal_surface_tags.extend(
+                domain_markers[f"back_{panel_id}"]["gmsh_tags"]
             )
 
         min_dist = []

From b1f4f5e26fec8e4d1b9e8c2e5f6861bc3bbc3a7a Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Mon, 8 Sep 2025 17:25:43 -0600
Subject: [PATCH 03/26] ignore panel surfaces marked trash, add back in checks
 on bbox, add new check on z_min/max clipping

---
 pvade/geometry/panels3d/DomainCreation.py | 70 +++++++++++------------
 1 file changed, 34 insertions(+), 36 deletions(-)

diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py
index 920eda6..ae793e9 100644
--- a/pvade/geometry/panels3d/DomainCreation.py
+++ b/pvade/geometry/panels3d/DomainCreation.py
@@ -337,9 +337,11 @@ def Rz(theta):
                             else:
                                 this_panel_transformed_com[key] = [com]
                         else:
-                            self._add_to_domain_markers(
-                                f"trash_{panel_ct:.0f}", [surf_id], "facet"
-                            )
+                            key = f"trash_{panel_ct:.0f}"
+                            if key in this_panel_transformed_com:
+                                this_panel_transformed_com[key].append(com)
+                            else:
+                                this_panel_transformed_com[key] = [com]
 
                 for key, val in this_panel_transformed_com.items():
                     for row_num, com in enumerate(val):
@@ -432,27 +434,6 @@ def Rz(theta):
                 else:
                     self.numpy_pt_total_array = np.copy(numpy_pt_panel_array)
 
-                # # Check that this panel still exists in the confines of the domain
-                # bbox = self.gmsh_model.occ.get_bounding_box(panel_tag_list)
-
-                # if bbox[0] < params.domain.x_min:
-                #     raise ValueError(
-                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past x_min wall."
-                #     )
-                # if bbox[1] < params.domain.y_min:
-                #     raise ValueError(
-                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past y_min wall."
-                #     )
-                # if bbox[3] > params.domain.x_max:
-                #     raise ValueError(
-                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past x_max wall."
-                #     )
-                # if bbox[4] > params.domain.y_max:
-                #     raise ValueError(
-                #         f"Panel with location (x, y) = ({xx}, {yy}) extends past y_max wall."
-                #     )
-            # self.gmsh_model.occ.synchronize()
-
         # Fragment all panels from the overall domain
         self.gmsh_model.occ.fragment(domain_tag_list, panel_tag_list)
 
@@ -469,45 +450,62 @@ def Rz(theta):
             surf_id = surf_tag[1]
             com = self.gmsh_model.occ.getCenterOfMass(self.ndim - 1, surf_id)
 
-            tagged_this_surface = False
+            located_this_surface = False
 
             # sturctures tagging
             if np.isclose(com[0], params.domain.x_min):
                 self._add_to_domain_markers("x_min", [surf_id], "facet")
-                tagged_this_surface = True
 
             elif np.allclose(com[0], params.domain.x_max):
                 self._add_to_domain_markers("x_max", [surf_id], "facet")
-                tagged_this_surface = True
 
             elif np.allclose(com[1], params.domain.y_min):
                 self._add_to_domain_markers("y_min", [surf_id], "facet")
-                tagged_this_surface = True
 
             elif np.allclose(com[1], params.domain.y_max):
                 self._add_to_domain_markers("y_max", [surf_id], "facet")
-                tagged_this_surface = True
 
             elif np.allclose(com[2], params.domain.z_min):
                 self._add_to_domain_markers("z_min", [surf_id], "facet")
-                tagged_this_surface = True
 
             elif np.allclose(com[2], params.domain.z_max):
                 self._add_to_domain_markers("z_max", [surf_id], "facet")
-                tagged_this_surface = True
 
             else:
                 for key, val in transformed_com.items():
                     for target_com in val:
                         # print(target_com)
                         if np.allclose(np.array(com), target_com):
-                            self._add_to_domain_markers(key, [surf_id], "facet")
-                            tagged_this_surface = True
+                            located_this_surface = True
+                            if "trash" not in key:
+                                self._add_to_domain_markers(key, [surf_id], "facet")
 
-            # if not tagged_this_surface:
-            #     print(f"Warning: Surface {surf_tag} has not been added to domain markers")
+                if not located_this_surface:
+                    print(
+                        f"Warning: Surface {surf_tag} has not been added to domain markers"
+                    )
+
+                # Since this is not one of the exterior walls, we should check if it extends
+                # past the boundaries x_min, x_max, ...
+                this_surf_bbox = self.gmsh_model.occ.get_bounding_box(
+                    self.ndim - 1, surf_id
+                )
 
-        # print(self.domain_markers)
+                # Test that the rotated point still exists in the box domain
+                if this_surf_bbox[0] < params.domain.x_min:
+                    raise ValueError(f"A panel extends past the x_min wall.")
+                if this_surf_bbox[0] > params.domain.x_max:
+                    raise ValueError(f"A panel extends past the x_max wall.")
+                if this_surf_bbox[1] < params.domain.y_min:
+                    raise ValueError(f"A panel extends past the y_min wall.")
+                if this_surf_bbox[1] > params.domain.y_max:
+                    raise ValueError(f"A panel extends past the y_max wall.")
+                if this_surf_bbox[2] < 0.0:
+                    raise ValueError(
+                        f"A panel extends past the z_min wall (ground level = 0.0)."
+                    )
+                if this_surf_bbox[2] > params.domain.z_max:
+                    raise ValueError(f"A panel extends past the z_max wall.")
 
         # Volumes are the entities with dimension equal to the mesh dimension
         vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim)

From 185b1fe6e1713c4e2e903b43538d02f961d573dd Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 12:40:42 -0600
Subject: [PATCH 04/26] change defaults for torque tube to zero to preserve old
 mesh behavior by default

---
 pvade/IO/input_schema.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index 335dacb..ddae4bb 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -199,13 +199,13 @@ properties:
         (x_min + x_spacing, y_min), etc., i.e., x-major direction. If argument is a single value, that tracking angle will be applied to the every panel."
         units: "degree"
       torque_tube_separation:
-        default: 0.4
+        default: 0.0
         minimum: 0.0
         maximum: 1.0
         type: "number"
         description: "The distance between the bottom of the panel structure and the centerline of the torque tube, in meters."
       torque_tube_radius:
-        default: 0.1
+        default: 0.0
         minimum: 0.0
         maximum: 1.0
         type: "number"

From 96672ad83dbbc461372232e347b2563a30ba1957 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 12:41:25 -0600
Subject: [PATCH 05/26] clean up facet marking tests and revert to old meshing
 behavior if torque tube is not specified (length zero)

---
 pvade/geometry/panels3d/DomainCreation.py | 228 ++++++++++++++--------
 1 file changed, 145 insertions(+), 83 deletions(-)

diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py
index ae793e9..31f4d42 100644
--- a/pvade/geometry/panels3d/DomainCreation.py
+++ b/pvade/geometry/panels3d/DomainCreation.py
@@ -158,6 +158,14 @@ def Rz(theta):
 
         transformed_com = {}
 
+        if (
+            params.pv_array.torque_tube_separation > 0.0
+            and params.pv_array.torque_tube_radius > 0.0
+        ):
+            modeling_torque_tube = True
+        else:
+            modeling_torque_tube = False
+
         for panel_id_y, yy in enumerate(y_centers):
             for panel_id_x, xx in enumerate(x_centers):
                 if isinstance(params.pv_array.tracker_angle, list):
@@ -182,30 +190,45 @@ def Rz(theta):
                         module_distances[module_id + 1] - module_distances[module_id]
                     )
 
-                    # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0)
-                    this_module = self.gmsh_model.occ.addBox(
-                        -half_chord,
-                        module_distances[module_id],
-                        params.pv_array.torque_tube_separation,
-                        params.pv_array.panel_chord,
-                        module_span,
-                        params.pv_array.panel_thickness,
-                    )
+                    if modeling_torque_tube:
+                        # Create an 0-tracking-degree panel centered at (x, y, z) = (0, 0, 0)
+                        this_module = self.gmsh_model.occ.addBox(
+                            -half_chord,
+                            module_distances[module_id],
+                            params.pv_array.torque_tube_separation,
+                            params.pv_array.panel_chord,
+                            module_span,
+                            params.pv_array.panel_thickness,
+                        )
 
-                    this_standoff = self.gmsh_model.occ.addBox(
-                        -params.pv_array.torque_tube_radius,
-                        module_distances[module_id],
-                        0.0,
-                        2.0 * params.pv_array.torque_tube_radius,
-                        module_span,
-                        params.pv_array.torque_tube_separation,
-                    )
+                        this_standoff = self.gmsh_model.occ.addBox(
+                            -params.pv_array.torque_tube_radius,
+                            module_distances[module_id],
+                            0.0,
+                            2.0 * params.pv_array.torque_tube_radius,
+                            module_span,
+                            params.pv_array.torque_tube_separation,
+                        )
+
+                        panel_tag_list.append((self.ndim, this_module))
+                        panel_tag_list.append((self.ndim, this_standoff))
+
+                        this_panel_tag_list.append((self.ndim, this_module))
+                        this_panel_tag_list.append((self.ndim, this_standoff))
+
+                    else:
+                        this_module = self.gmsh_model.occ.addBox(
+                            -half_chord,
+                            module_distances[module_id],
+                            0.0,
+                            params.pv_array.panel_chord,
+                            module_span,
+                            params.pv_array.panel_thickness,
+                        )
 
-                    panel_tag_list.append((self.ndim, this_module))
-                    panel_tag_list.append((self.ndim, this_standoff))
+                        panel_tag_list.append((self.ndim, this_module))
 
-                    this_panel_tag_list.append((self.ndim, this_module))
-                    this_panel_tag_list.append((self.ndim, this_standoff))
+                        this_panel_tag_list.append((self.ndim, this_module))
 
                 numpy_pt_list = []
                 embedded_lines_tag_list = []
@@ -248,16 +271,65 @@ def Rz(theta):
 
                     for fp in fixation_pts_list:
                         # FIXME: don't add the fixation points into the numpy tagging for now
-                        numpy_pt_list.append(
-                            [
-                                -params.pv_array.torque_tube_radius,
-                                -half_span + fp,
-                                0.0,
-                                params.pv_array.torque_tube_radius,
-                                -half_span + fp,
-                                0.0,
-                            ]
-                        )
+                        if modeling_torque_tube:
+                            numpy_pt_list.append(
+                                [
+                                    -params.pv_array.torque_tube_radius,
+                                    -half_span + fp,
+                                    0.0,
+                                    params.pv_array.torque_tube_radius,
+                                    -half_span + fp,
+                                    0.0,
+                                ]
+                            )
+
+                        else:
+                            numpy_pt_list.append(
+                                [
+                                    -half_chord,
+                                    -half_span + fp,
+                                    0.0,
+                                    half_chord,
+                                    -half_span + fp,
+                                    0.0,
+                                ]
+                            )
+
+                        # If the separation line already exists at a module division,
+                        # there's no need to redraw it, but otherwise, add the gmsh points
+                        # to the model for embedding
+                        if not np.any(np.isclose(-half_span + fp, module_distances)):
+                            print(
+                                f"Embedding a line for a no-deformation boundary condition."
+                            )
+
+                            if modeling_torque_tube:
+                                pt_1 = self.gmsh_model.occ.addPoint(
+                                    -params.pv_array.torque_tube_radius,
+                                    -half_span + fp,
+                                    0.0,
+                                )
+                                pt_2 = self.gmsh_model.occ.addPoint(
+                                    params.pv_array.torque_tube_radius,
+                                    -half_span + fp,
+                                    0.0,
+                                )
+
+                            else:
+                                pt_1 = self.gmsh_model.occ.addPoint(
+                                    -half_chord, -half_span + fp, 0.0
+                                )
+                                pt_2 = self.gmsh_model.occ.addPoint(
+                                    half_chord, -half_span + fp, 0.0
+                                )
+
+                            fixation_line_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
+                            fixation_line_tag = (1, fixation_line_id)
+                            embedded_lines_tag_list.append(fixation_line_tag)
+                        else:
+                            print(
+                                f"Applying no-deformation boundary condition at {fp}."
+                            )
 
                 # Fragment the lines into the surfaces (equivalent to embedding these lines)
                 self.gmsh_model.occ.fragment(
@@ -277,71 +349,61 @@ def Rz(theta):
                         surf_id = surf_tag[1]
                         com = self.gmsh_model.occ.getCenterOfMass(surf_dim, surf_id)
 
+                        target_key = None
+
                         # sturctures tagging
-                        if np.isclose(com[0], -half_chord) or np.isclose(
-                            com[0], -params.pv_array.torque_tube_radius
+                        if (
+                            np.isclose(com[0], -half_chord)
+                            or np.isclose(com[0], -params.pv_array.torque_tube_radius)
+                            and modeling_torque_tube
                         ):
-                            key = f"left_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
+                            target_key = f"left_{panel_ct:.0f}"
 
-                        elif np.isclose(com[0], half_chord) or np.isclose(
-                            com[0], params.pv_array.torque_tube_radius
+                        elif (
+                            np.isclose(com[0], half_chord)
+                            or np.isclose(com[0], params.pv_array.torque_tube_radius)
+                            and modeling_torque_tube
                         ):
-                            key = f"right_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
+                            target_key = f"right_{panel_ct:.0f}"
 
                         elif np.isclose(com[1], -half_span):
-                            key = f"front_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
+                            target_key = f"front_{panel_ct:.0f}"
 
                         elif np.isclose(com[1], half_span):
-                            key = f"back_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
-
-                        elif np.isclose(
-                            com[2], params.pv_array.torque_tube_separation
-                        ) and not np.isclose(com[0], 0.0):
-                            key = f"bottom_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
+                            target_key = f"back_{panel_ct:.0f}"
+
+                        elif (
+                            np.isclose(com[2], params.pv_array.torque_tube_separation)
+                            and not np.isclose(com[0], 0.0)
+                            and modeling_torque_tube
+                            or np.isclose(com[2], 0.0)
+                            and not modeling_torque_tube
+                        ):
+                            target_key = f"bottom_{panel_ct:.0f}"
 
-                        elif np.isclose(com[2], 0.0):
-                            key = f"torque_tube_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
+                        elif np.isclose(com[2], 0.0) and modeling_torque_tube:
+                            target_key = f"torque_tube_{panel_ct:.0f}"
 
-                        elif np.isclose(
-                            com[2],
-                            params.pv_array.torque_tube_separation
-                            + params.pv_array.panel_thickness,
+                        elif (
+                            np.isclose(
+                                com[2],
+                                params.pv_array.torque_tube_separation
+                                + params.pv_array.panel_thickness,
+                            )
+                            and modeling_torque_tube
+                            or np.isclose(com[2], params.pv_array.panel_thickness)
+                            and not modeling_torque_tube
                         ):
-                            key = f"top_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
-                            else:
-                                this_panel_transformed_com[key] = [com]
+                            target_key = f"top_{panel_ct:.0f}"
+
                         else:
-                            key = f"trash_{panel_ct:.0f}"
-                            if key in this_panel_transformed_com:
-                                this_panel_transformed_com[key].append(com)
+                            target_key = f"trash_{panel_ct:.0f}"
+
+                        if target_key is not None:
+                            if target_key in this_panel_transformed_com:
+                                this_panel_transformed_com[target_key].append(com)
                             else:
-                                this_panel_transformed_com[key] = [com]
+                                this_panel_transformed_com[target_key] = [com]
 
                 for key, val in this_panel_transformed_com.items():
                     for row_num, com in enumerate(val):

From bb62814666b81785342c93cee1caf8996a230a25 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 12:42:19 -0600
Subject: [PATCH 06/26] lower elevation by half panel thickness to coincide
 with new meshing strategy (panel bottom surface = 0 (z))

---
 pvade/tests/input/yaml/sim_params.yaml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvade/tests/input/yaml/sim_params.yaml b/pvade/tests/input/yaml/sim_params.yaml
index 7b6e4a6..2e63349 100644
--- a/pvade/tests/input/yaml/sim_params.yaml
+++ b/pvade/tests/input/yaml/sim_params.yaml
@@ -16,11 +16,11 @@ domain:
 pv_array:
   stream_rows: 7
   span_rows: 1
-  elevation: 1.5
+  elevation: 1.45 # 1.5 - 0.5*0.1 = 1.4
   stream_spacing: 7.0
   panel_chord: 2.0
   panel_span: 7.0
-  panel_thickness: 0.1
+  panel_thickness: 0.1 
   tracker_angle: -30.0
 solver:
   dt: 0.005

From 90a7db61356b38174f311e03ee2cd2b555141cb8 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 13:36:45 -0600
Subject: [PATCH 07/26] change the box to span [0,2] in the z direction vs
 being 0 centered

---
 pvade/tests/input/yaml/embedded_box.yaml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pvade/tests/input/yaml/embedded_box.yaml b/pvade/tests/input/yaml/embedded_box.yaml
index 0bd12d5..ad13c66 100644
--- a/pvade/tests/input/yaml/embedded_box.yaml
+++ b/pvade/tests/input/yaml/embedded_box.yaml
@@ -9,15 +9,15 @@ domain:
   x_max: 1.0
   y_min: -1.0
   y_max: 1.0
-  z_min: -1.0
-  z_max: 1.0
+  z_min: 0.0
+  z_max: 2.0
   l_char: 0.05
 pv_array:
   stream_rows: 1
   span_rows: 1
   stream_spacing: 1.0
   span_spacing: 1.0
-  elevation: 0.0
+  elevation: 0.5 # 0.0
   panel_chord: 1.2
   panel_span: 1.1
   panel_thickness: 1.0

From 251e117391c9030366aeadb7c34709ac24705dc3 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 13:37:58 -0600
Subject: [PATCH 08/26] update so that a panel's elevation is defined as the
 bottom of the panel to the ground

---
 pvade/tests/test_fsi_mesh.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pvade/tests/test_fsi_mesh.py b/pvade/tests/test_fsi_mesh.py
index 2d4f475..eba906e 100644
--- a/pvade/tests/test_fsi_mesh.py
+++ b/pvade/tests/test_fsi_mesh.py
@@ -148,10 +148,10 @@ def test_meshing_3dpanels_rotations(wind_direction, num_stream_rows, num_span_ro
             # Create the 4 corners of this table corresponding to the *top* surface (+0.5*thickness)
             top_surface_corners = np.array(
                 [
-                    [-0.5 * chord, -0.5 * span, 0.5 * thickness],
-                    [0.5 * chord, -0.5 * span, 0.5 * thickness],
-                    [0.5 * chord, 0.5 * span, 0.5 * thickness],
-                    [-0.5 * chord, 0.5 * span, 0.5 * thickness],
+                    [-0.5 * chord, -0.5 * span, thickness],
+                    [0.5 * chord, -0.5 * span, thickness],
+                    [0.5 * chord, 0.5 * span, thickness],
+                    [-0.5 * chord, 0.5 * span, thickness],
                 ]
             )
 

From 9d77229ed174560e523c5b71f7540762d387b93c Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 13:38:58 -0600
Subject: [PATCH 09/26] change dz to account for new panel position, change
 truth values to account for new positioning conventions

---
 pvade/tests/test_mesh_movement.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/pvade/tests/test_mesh_movement.py b/pvade/tests/test_mesh_movement.py
index 4de08c8..ec0798d 100644
--- a/pvade/tests/test_mesh_movement.py
+++ b/pvade/tests/test_mesh_movement.py
@@ -49,7 +49,7 @@ def test_calc_distance_to_panel_surface():
 
     dx = params.domain.x_max - 0.5 * params.pv_array.panel_chord
     dy = params.domain.y_max - 0.5 * params.pv_array.panel_span
-    dz = params.domain.z_max - 0.5 * params.pv_array.panel_thickness
+    dz = params.pv_array.elevation
     truth_max_dist = np.sqrt(dx * dx + dy * dy + dz * dz)
 
     assert np.isclose(max_dist, truth_max_dist)
@@ -113,9 +113,7 @@ def __init__(self, domain, x_shift, y_shift, z_shift):
     assert np.isclose(
         np.amin(structure_coords_before[:, 1]), -0.5 * params.pv_array.panel_span
     )
-    assert np.isclose(
-        np.amin(structure_coords_before[:, 2]), -0.5 * params.pv_array.panel_thickness
-    )
+    assert np.isclose(np.amin(structure_coords_before[:, 2]), params.pv_array.elevation)
 
     assert np.isclose(
         np.amax(structure_coords_before[:, 0]), 0.5 * params.pv_array.panel_chord
@@ -124,7 +122,8 @@ def __init__(self, domain, x_shift, y_shift, z_shift):
         np.amax(structure_coords_before[:, 1]), 0.5 * params.pv_array.panel_span
     )
     assert np.isclose(
-        np.amax(structure_coords_before[:, 2]), 0.5 * params.pv_array.panel_thickness
+        np.amax(structure_coords_before[:, 2]),
+        params.pv_array.elevation + params.pv_array.panel_thickness,
     )
 
     # Move the mesh by the amount prescribed in u_delta

From ca9df0043520112fac50c1b39f8bd9c16b0812f6 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 16:42:10 -0600
Subject: [PATCH 10/26] change the location of the tracked corner acceleration
 consistent with new meshing convention

---
 pvade/structure/StructureMain.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/pvade/structure/StructureMain.py b/pvade/structure/StructureMain.py
index 6e55265..67b5cba 100644
--- a/pvade/structure/StructureMain.py
+++ b/pvade/structure/StructureMain.py
@@ -110,9 +110,21 @@ def _north_east_corner(x):
                 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]
+                """
+                y
+
+                ^
+                | 
+                |---o NE corner (measured on bottom of panel)
+                |   |
+                |   |
+                |-----> x
+                """
+
+                corner = [
+                    0.5 * params.pv_array.panel_chord * np.cos(tracker_angle_rad),
+                    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)

From b6a5fd3471209a4a093b0f1796573f156b9a64da Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 16:44:22 -0600
Subject: [PATCH 11/26] changing variable names from north_west to north_east

---
 pvade/structure/ElasticityAnalysis.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 4a84d98..331bd1a 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -485,25 +485,27 @@ def solve(self, params, dataIO, structure):
         try:
             idx = structure.north_east_corner_dofs[0]
             # idx = self.north_east_corner_dofs[0]
-            nw_corner_accel = self.u.x.array[
+            north_east_corner_acccel = self.u.x.array[
                 structure.ndim * idx : structure.ndim * idx + structure.ndim
             ].astype(np.float64)
-            print(nw_corner_accel)
+            print(north_east_corner_acccel)
         except:
-            nw_corner_accel = np.zeros(structure.ndim, dtype=np.float64)
+            north_east_corner_acccel = np.zeros(structure.ndim, dtype=np.float64)
 
-        nw_corner_accel_global = np.zeros(
+        north_east_corner_accel_global = np.zeros(
             (self.num_procs, structure.ndim), dtype=np.float64
         )
 
-        self.comm.Gather(nw_corner_accel, nw_corner_accel_global, root=0)
+        self.comm.Gather(
+            north_east_corner_acccel, north_east_corner_accel_global, root=0
+        )
 
-        # print(f"Acceleration at North West corner = {nw_corner_accel}")
+        # print(f"Acceleration at North West corner = {north_east_corner_acccel}")
 
         if self.rank == 0:
-            norm2 = np.sum(nw_corner_accel_global**2, axis=1)
+            norm2 = np.sum(north_east_corner_accel_global**2, axis=1)
             max_norm2_idx = np.argmax(norm2)
-            np_accel = nw_corner_accel_global[max_norm2_idx, :]
+            np_accel = north_east_corner_accel_global[max_norm2_idx, :]
 
             accel_pos_filename = os.path.join(
                 params.general.output_dir_sol, "accel_pos.csv"

From 1929dc97be90bdcd6e41e276f6b7528fee13a377 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Tue, 9 Sep 2025 16:46:01 -0600
Subject: [PATCH 12/26] remove misleading comment

---
 pvade/structure/ElasticityAnalysis.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 331bd1a..92c43c0 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -500,8 +500,6 @@ def solve(self, params, dataIO, structure):
             north_east_corner_acccel, north_east_corner_accel_global, root=0
         )
 
-        # print(f"Acceleration at North West corner = {north_east_corner_acccel}")
-
         if self.rank == 0:
             norm2 = np.sum(north_east_corner_accel_global**2, axis=1)
             max_norm2_idx = np.argmax(norm2)

From 9b14d506369cd5423c10a50c879b286c1609ed47 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Wed, 10 Sep 2025 09:18:59 -0600
Subject: [PATCH 13/26] naive parsing to avoid issue of gha node name having
 underscore

---
 .github/workflows/test_pvade.yaml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.github/workflows/test_pvade.yaml b/.github/workflows/test_pvade.yaml
index 897d46c..477ec3d 100644
--- a/.github/workflows/test_pvade.yaml
+++ b/.github/workflows/test_pvade.yaml
@@ -33,6 +33,8 @@ jobs:
         shell: bash -l {0}
         run: PYTHONPATH=. pytest -sv pvade/tests/
       - name: test all inputs
+        env:
+          OMPI_MCA_regx: "naive"
         shell: bash -l {0}
         run: pytest -sv  test_all_inputs.py 
 

From d10aec8719d80acf4a5587f36db2e210c732c6e2 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Wed, 10 Sep 2025 09:52:35 -0600
Subject: [PATCH 14/26] save both acceleration and deformation in csv file,
 clean up unclear variable names

---
 pvade/structure/ElasticityAnalysis.py | 69 ++++++++++++++++++++++-----
 1 file changed, 56 insertions(+), 13 deletions(-)

diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 92c43c0..85ea43e 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -484,26 +484,52 @@ def solve(self, params, dataIO, structure):
 
         try:
             idx = structure.north_east_corner_dofs[0]
-            # idx = self.north_east_corner_dofs[0]
-            north_east_corner_acccel = self.u.x.array[
+            north_east_corner_deformation = self.u.x.array[
                 structure.ndim * idx : structure.ndim * idx + structure.ndim
             ].astype(np.float64)
-            print(north_east_corner_acccel)
+            print(f"Deformation: {north_east_corner_deformation}")
+
+        except:
+            north_east_corner_deformation = np.zeros(structure.ndim, dtype=np.float64)
+
+        try:
+            idx = structure.north_east_corner_dofs[0]
+            north_east_corner_acceleration = self.a_old.x.array[
+                structure.ndim * idx : structure.ndim * idx + structure.ndim
+            ].astype(np.float64)
+            print(f"Acceleration: {north_east_corner_acceleration}")
+
         except:
-            north_east_corner_acccel = np.zeros(structure.ndim, dtype=np.float64)
+            north_east_corner_acceleration = np.zeros(structure.ndim, dtype=np.float64)
 
-        north_east_corner_accel_global = np.zeros(
+        # Initialize a buffer to collect everything into
+        north_east_corner_deformation_global = np.zeros(
             (self.num_procs, structure.ndim), dtype=np.float64
         )
 
+        north_east_corner_acceleration_global = np.zeros(
+            (self.num_procs, structure.ndim), dtype=np.float64
+        )
+
+        # Gather all points (many of which are zeros) to rank 0
+        self.comm.Gather(
+            north_east_corner_deformation, north_east_corner_deformation_global, root=0
+        )
+
         self.comm.Gather(
-            north_east_corner_acccel, north_east_corner_accel_global, root=0
+            north_east_corner_acceleration,
+            north_east_corner_acceleration_global,
+            root=0,
         )
 
         if self.rank == 0:
-            norm2 = np.sum(north_east_corner_accel_global**2, axis=1)
+            norm2 = np.sum(north_east_corner_deformation_global**2, axis=1)
+            max_norm2_idx = np.argmax(norm2)
+            np_deformation = north_east_corner_deformation_global[max_norm2_idx, :]
+
+            norm2 = np.sum(north_east_corner_acceleration_global**2, axis=1)
             max_norm2_idx = np.argmax(norm2)
-            np_accel = north_east_corner_accel_global[max_norm2_idx, :]
+            np_acceleration = north_east_corner_acceleration_global[max_norm2_idx, :]
 
             accel_pos_filename = os.path.join(
                 params.general.output_dir_sol, "accel_pos.csv"
@@ -512,18 +538,35 @@ def solve(self, params, dataIO, structure):
             if self.first_call_to_solver:
 
                 with open(accel_pos_filename, "w") as fp:
-                    fp.write("#x-pos,y-pos,z-pos\n")
                     if structure.ndim == 3:
-                        fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n")
+                        fp.write(
+                            "#x-deformation,y-deformation,z-deformation,x-acceleration,y-acceleration,z-acceleration\n"
+                        )
+                        fp.write(
+                            f"{np_deformation[0]},{np_deformation[1]},{np_deformation[2]},"
+                        )
+                        fp.write(
+                            f"{np_acceleration[0]},{np_acceleration[1]},{np_acceleration[2]}\n"
+                        )
                     elif structure.ndim == 2:
-                        fp.write(f"{np_accel[0]},{np_accel[1]}\n")
+                        fp.write(
+                            "#x-deformation,y-deformation,x-acceleration,y-acceleration\n"
+                        )
+                        fp.write(f"{np_deformation[0]},{np_deformation[1]},")
+                        fp.write(f"{np_acceleration[0]},{np_acceleration[1]}\n")
 
             else:
                 with open(accel_pos_filename, "a") as fp:
                     if structure.ndim == 3:
-                        fp.write(f"{np_accel[0]},{np_accel[1]},{np_accel[2]}\n")
+                        fp.write(
+                            f"{np_deformation[0]},{np_deformation[1]},{np_deformation[2]},"
+                        )
+                        fp.write(
+                            f"{np_acceleration[0]},{np_acceleration[1]},{np_acceleration[2]}\n"
+                        )
                     elif structure.ndim == 2:
-                        fp.write(f"{np_accel[0]},{np_accel[1]}\n")
+                        fp.write(f"{np_deformation[0]},{np_deformation[1]},")
+                        fp.write(f"{np_acceleration[0]},{np_acceleration[1]}\n")
 
         if self.first_call_to_solver:
             self.first_call_to_solver = False

From 2a0a73ab4aff10ee8ae49f7e1244358ea61b4931 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Wed, 10 Sep 2025 09:52:56 -0600
Subject: [PATCH 15/26] change facets to vertices in variable name for clarity

---
 pvade/structure/StructureMain.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pvade/structure/StructureMain.py b/pvade/structure/StructureMain.py
index 67b5cba..ed3a756 100644
--- a/pvade/structure/StructureMain.py
+++ b/pvade/structure/StructureMain.py
@@ -133,12 +133,12 @@ def _north_east_corner(x):
 
             return north_east_corner
 
-        north_east_corner_facets = dolfinx.mesh.locate_entities_boundary(
+        north_east_corner_vertices = dolfinx.mesh.locate_entities_boundary(
             domain.structure.msh, 0, _north_east_corner
         )
 
         self.north_east_corner_dofs = dolfinx.fem.locate_dofs_topological(
-            self.elasticity.V, 0, north_east_corner_facets
+            self.elasticity.V, 0, north_east_corner_vertices
         )
 
     def build_boundary_conditions(self, domain, params):

From a2033edcb4706ecdafa9b42e39910b789a1081a2 Mon Sep 17 00:00:00 2001
From: Young <eyoung55@gmail.com>
Date: Wed, 10 Sep 2025 09:53:35 -0600
Subject: [PATCH 16/26] only compare the deformation outputs (first two cols)
 during the check to avoid dimension mismatch

---
 pvade/tests/test_solve.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pvade/tests/test_solve.py b/pvade/tests/test_solve.py
index 8734207..77d5457 100644
--- a/pvade/tests/test_solve.py
+++ b/pvade/tests/test_solve.py
@@ -310,7 +310,7 @@ def test_fsi2():
 
     # print('pos_data = ', pos_data)
 
-    assert np.allclose(pos_data, pos_data_truth)
+    assert np.allclose(pos_data[:, 0:2], pos_data_truth)
     print(lift_and_drag_data)
 
     # assert np.allclose(lift_and_drag_data[:, 0:3], lift_and_drag_data_truth[:, 0:3]) # needs new truth values to pass, mesh has changed

From 2a18d4b366f5c57a8fe609c80a1854e01ee3d46f Mon Sep 17 00:00:00 2001
From: xinhe <xh2205@columbia.edu>
Date: Thu, 23 Oct 2025 13:25:53 -0600
Subject: [PATCH 17/26] update the case study input

---
 input/duramat_case_study.yaml | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/input/duramat_case_study.yaml b/input/duramat_case_study.yaml
index fb5939e..d590670 100644
--- a/input/duramat_case_study.yaml
+++ b/input/duramat_case_study.yaml
@@ -23,6 +23,11 @@ pv_array:
   elevation: 2.1
   tracker_angle: 0.0
   span_fixation_pts: [13.2]
+  torque_tube_separation: 0.2 # gap between panel and tube center
+  torque_tube_radius: 0.1 # radius of the torque tube
+  modules_per_span: 10
+  block_chord_div_by_panel_chord: 0.02
+  
 solver:
   dt: 0.01
   t_final: 20.0

From 5e9a2f98d72d2ffd41269082ee7bec1c354bcbe4 Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Fri, 7 Nov 2025 13:59:01 -0700
Subject: [PATCH 18/26] implement structure simplification

---
 input/duramat_case_study.yaml             |   8 +-
 pvade/IO/input_schema.yaml                |  28 +-
 pvade/fluid/FlowManager.py                | 238 +++++++----
 pvade/geometry/MeshManager.py             |   2 +-
 pvade/geometry/panels3d/DomainCreation.py | 495 +++++++++++++++-------
 pvade/structure/ElasticityAnalysis.py     | 180 +++++++-
 pvade/structure/StructureMain.py          |  18 +-
 pvade/structure/boundary_conditions.py    |  34 +-
 pvade_main.py                             |   8 +-
 9 files changed, 749 insertions(+), 262 deletions(-)

diff --git a/input/duramat_case_study.yaml b/input/duramat_case_study.yaml
index d590670..2221b4e 100644
--- a/input/duramat_case_study.yaml
+++ b/input/duramat_case_study.yaml
@@ -24,7 +24,8 @@ pv_array:
   tracker_angle: 0.0
   span_fixation_pts: [13.2]
   torque_tube_separation: 0.2 # gap between panel and tube center
-  torque_tube_radius: 0.1 # radius of the torque tube
+  torque_tube_outer_radius: 0.1 # radius of the torque tube
+  torque_tube_inner_radius: 0.09 # radius of the torque tube
   modules_per_span: 10
   block_chord_div_by_panel_chord: 0.02
   
@@ -63,4 +64,7 @@ structure:
   bc_list: []
   motor_connection: true
   tube_connection: true
-  beta_relaxation: 0.5
\ No newline at end of file
+  beta_relaxation: 0.5
+  elasticity_modulus_tube: 2.0e+11
+  poissons_ratio_tube: 0.3
+  density_tube: 7800.0
\ No newline at end of file
diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index c127d0c..9822740 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -204,7 +204,13 @@ properties:
         maximum: 1.0
         type: "number"
         description: "The distance between the bottom of the panel structure and the centerline of the torque tube, in meters."
-      torque_tube_radius:
+      torque_tube_outer_radius:
+        default: 0.0
+        minimum: 0.0
+        maximum: 1.0
+        type: "number"
+        description: "The radius of the torque tube, in meters."
+      torque_tube_inner_radius:
         default: 0.0
         minimum: 0.0
         maximum: 1.0
@@ -215,6 +221,12 @@ properties:
         minimum: 1
         type: "integer"
         description: "The number of modules/panels per every panel_span distance. This can be thought of as the number of modules/panels per each table, counted in the spanwise direction only (i.e., 2-in-portait systems should report only the number in the spanwise direction, vs 2x that value)."
+      block_chord_div_by_panel_chord:
+        default: 0.1
+        minimum: 0.01
+        type: "number"
+        description: "block length or width divided by panel chord"
+      
   solver:
     additionalProperties: false
     type: "object"
@@ -665,6 +677,20 @@ properties:
         type: "number"
         description: "poissons ratio of the panel structure."
         units: "None"
+      elasticity_modulus_tube:
+        default: 200.0e+9
+        minimum: 1.0e+2
+        maximum: 500.0e+9
+        type: "number"
+        description: "The effective Young's modulus of the torque tube."
+        units: "Pascal"  
+      poissons_ratio_tube: 
+        default: 0.3
+        minimum: 0.1
+        maximum: 0.9
+        type: "number"
+        description: "poissons ratio of the torque tube."
+        units: "None"
       body_force_x: 
         default:  100
         type: "number"
diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index bd2c5c4..0f2a1c8 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -242,6 +242,54 @@ def build_forms(self, domain, params):
         self.p_k = dolfinx.fem.Function(self.Q, name="pressure")
         self.p_k1 = dolfinx.fem.Function(self.Q)
 
+        # define function to save x, y, z coordinates
+        self.spatial_X_coords = dolfinx.fem.Function(self.Q, name="X_coords")
+        self.spatial_Y_coords = dolfinx.fem.Function(self.Q, name="Y_coords")
+        self.spatial_Z_coords = dolfinx.fem.Function(self.Q, name="Z_coords")
+
+        class FillFunctionWithXCoords:
+            def __init__(self):
+                pass
+
+            def __call__(self, x):
+                coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
+
+                
+                coords[0] = x[0]
+
+                return coords
+
+        self.spatial_X_coords.interpolate(FillFunctionWithXCoords())
+        
+        class FillFunctionWithYCoords:
+            def __init__(self):
+                pass
+
+            def __call__(self, x):
+                coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
+
+                
+                coords[0] = x[1]
+
+                return coords
+
+        self.spatial_Y_coords.interpolate(FillFunctionWithYCoords())
+
+        class FillFunctionWithZCoords:
+            def __init__(self):
+                pass
+
+            def __call__(self, x):
+                coords = np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
+
+                
+                coords[0] = x[2]
+                
+
+                return coords
+
+        self.spatial_Z_coords.interpolate(FillFunctionWithZCoords())
+
         # initial conditions
         if params.fluid.initialize_with_inflow_bc:
             # self.inflow_profile = dolfinx.fem.Function(self.V)
@@ -534,94 +582,132 @@ def _all_interior_surfaces(x):
         ds_fluid = ufl.Measure(
             "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
         )
+    
+    def compute_panel_torques(self, domain, params):
 
-        self.integrated_force_x_form = []
-        self.integrated_force_y_form = []
-        self.integrated_force_z_form = []
+        ds_fluid = ufl.Measure(
+            "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
+        )
 
         for panel_id in range(
             int(params.pv_array.stream_rows * params.pv_array.span_rows)
         ):
-            # for loc in ["top_0", "bottom_0", "left_0", "right_0"]:
-            #     idx = domain.domain_markers[loc]["idx"]
-            #     s = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1.0*ds_fluid(idx)))
-            #     print(f"loc = {loc}, idx = {idx}, s = {s}")
-
-            self.integrated_force_x_form.append(0)
-            self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
-                domain.domain_markers[f"left_{panel_id:.0f}"]["idx"]
-            )
-            self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
-                domain.domain_markers[f"top_{panel_id:.0f}"]["idx"]
-            )
-            self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
-                domain.domain_markers[f"right_{panel_id:.0f}"]["idx"]
-            )
-            self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
-                domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"]
-            )
-            if self.ndim == 3:
-                self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
-                    domain.domain_markers[f"front_{panel_id:.0f}"]["idx"]
-                )
-                self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
-                    domain.domain_markers[f"back_{panel_id:.0f}"]["idx"]
-                )
-
-            self.integrated_force_x_form[-1] = dolfinx.fem.form(
-                self.integrated_force_x_form[-1]
-            )
-
-            self.integrated_force_y_form.append(0)
-            self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
-                domain.domain_markers[f"left_{panel_id:.0f}"]["idx"]
-            )
-            self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
-                domain.domain_markers[f"top_{panel_id:.0f}"]["idx"]
-            )
-            self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
-                domain.domain_markers[f"right_{panel_id:.0f}"]["idx"]
-            )
-            self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
-                domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"]
-            )
-            if self.ndim == 3:
-                self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
-                    domain.domain_markers[f"front_{panel_id:.0f}"]["idx"]
-                )
-                self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
-                    domain.domain_markers[f"back_{panel_id:.0f}"]["idx"]
-                )
-
-            self.integrated_force_y_form[-1] = dolfinx.fem.form(
-                self.integrated_force_y_form[-1]
-            )
+            for module_id in range(params.pv_array.modules_per_span):
+                
+                total_torque = 0
 
-            self.integrated_force_z_form.append(0)
-            if self.ndim == 3:
-                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
-                    domain.domain_markers[f"left_{panel_id:.0f}"]["idx"]
-                )
-                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
-                    domain.domain_markers[f"top_{panel_id:.0f}"]["idx"]
+                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
                 )
-                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
-                    domain.domain_markers[f"right_{panel_id:.0f}"]["idx"]
+                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid(
+                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
                 )
-                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
-                    domain.domain_markers[f"bottom_{panel_id:.0f}"]["idx"]
+                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid(
+                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
                 )
-                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
-                    domain.domain_markers[f"front_{panel_id:.0f}"]["idx"]
-                )
-                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
-                    domain.domain_markers[f"back_{panel_id:.0f}"]["idx"]
+                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
                 )
+                attr_name = f"total_torque_panel_{panel_id:.0f}_{module_id:.0f}"
+                setattr(self, attr_name, total_torque)
 
-                self.integrated_force_z_form[-1] = dolfinx.fem.form(
-                    self.integrated_force_z_form[-1]
-                )
+    def compute_double_integral_panel_torques(self, domain, params):
+        torque_function_on_panels = dolfinx.fem.Function(self.Q)
+
+        torque_function_on_panels.x.array[:] = self.spatial_X_coords.x.array[:]*self.traction[2].x.array[:]-self.spatial_Z_coords.x.array[:]*self.traction[0].x.array[:]
 
+        ds_fluid = ufl.Measure(
+            "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
+        )
+        for panel_id in range(
+            int(params.pv_array.stream_rows * params.pv_array.span_rows)
+        ):
+            for module_id in range(params.pv_array.modules_per_span):
+                
+                torque_double_integral = 0
+
+                whole_top_surface_facet_index = domain.fluid.facet_tags.find(domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"])
+                whole_bot_surface_facet_index = domain.fluid.facet_tags.find(domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"])
+                whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(self.fluid.mesh, 2, whole_top_surface_facet_index)
+                whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(self.fluid.mesh, 2, whole_bot_surface_facet_index)
+
+                dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh)
+                dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh)
+
+                whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1))
+                whole_top_submesh_func = dolfinx.fem.Function(whole_top_submesh_function_space)
+                whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1))
+                whole_bot_submesh_func = dolfinx.fem.Function(whole_bot_submesh_function_space)
+
+                coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates()
+                coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates()
+
+                top_integrated_func_along_x = np.zeros_like(whole_top_submesh_func.x.array[:])
+
+                bot_integrated_func_along_x = np.zeros_like(whole_bot_submesh_func.x.array[:])
+                
+                # Sort by x first, then y and z (group by x-levels)
+                sort_idx_top = np.lexsort((coords_top[:, 1], coords_top[:, 0]))  # x is sorted, from small to large
+                sort_idx_bot = np.lexsort((coords_bot[:, 1], coords_bot[:, 0]))  # x is sorted, from small to large
+                
+                top_coords_sorted = coords_top[sort_idx_top] #sorted from smallest x to largest x
+                bot_coords_sorted = coords_bot[sort_idx_bot] #sorted from smallest x to largest x
+                
+
+                ### integral of top surface:
+                # Unique x-values (up to numerical tolerance)
+                xs_top = np.unique(np.round(top_coords_sorted[:, 0], decimals=6))
+                
+                for x_value in xs_top:
+                    # Mask for points at this x-level
+                    x_mask_top = np.isclose(top_coords_sorted[:, 0], x_value, atol=1e-6)
+
+                    facet_centers_top = dolfinx.mesh.compute_midpoints(whole_top_surf_submesh, 2, np.array(np.arange(whole_top_surf_submesh.topology.index_map(2).size_local), dtype=np.int32))
+                    # Identify cells where the midpoint x-coordinate is less than 1
+                    marked_facet_top = np.where(facet_centers_top[:, 0] <= x_value)[0]
+                    
+                    # Create a submesh from those cells
+                    submesh_top_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_top_surf_submesh, 2, marked_facet_top.astype(np.int32))
+                    top_sub_function_space = dolfinx.fem.FunctionSpace(submesh_top_surface_part, ('CG', 1))
+                    top_sub_func = dolfinx.fem.Function(top_sub_function_space)
+                    top_sub_func.interpolate(torque_function_on_panels)
+
+                    dx_top = ufl.Measure("dx", domain=submesh_top_surface_part)
+                    integrated_top = dolfinx.fem.assemble_scalar(dolfinx.fem.form(top_sub_func * dx_top))
+                    top_integrated_func_along_x[sort_idx_top[x_mask_top]] = integrated_top
+
+                ### integral of bot surface:
+                # Unique x-values (up to numerical tolerance)
+                xs_bot = np.unique(np.round(bot_coords_sorted[:, 0], decimals=6))
+
+                for x_value in xs_bot:
+
+                    # Mask for points at this x-level
+                    x_mask_bot = np.isclose(bot_coords_sorted[:, 0], x_value, atol=1e-6)
+                    facet_centers_bot = dolfinx.mesh.compute_midpoints(whole_bot_surf_submesh, 2, np.array(np.arange(whole_bot_surf_submesh.topology.index_map(2).size_local), dtype=np.int32))
+                    # Identify cells where the midpoint x-coordinate is less than 1
+                    marked_facet_bot = np.where(facet_centers_bot[:, 0] <= x_value)[0]
+                    # Create a submesh from those cells
+                    submesh_bot_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_bot_surf_submesh, 2, marked_facet_bot.astype(np.int32))
+                    bot_sub_function_space = dolfinx.fem.FunctionSpace(submesh_bot_surface_part, ('CG', 1))
+                    bot_sub_func = dolfinx.fem.Function(bot_sub_function_space)
+                    bot_sub_func.interpolate(torque_function_on_panels)
+
+                    dx_bot = ufl.Measure("dx", domain=submesh_bot_surface_part)
+                    integrated_bot = dolfinx.fem.assemble_scalar(dolfinx.fem.form(bot_sub_func * dx_bot))
+                    bot_integrated_func_along_x[sort_idx_bot[x_mask_bot]] = integrated_bot
+
+                single_integral_function_top = dolfinx.fem.Function(whole_top_submesh_function_space, name='single_integral_function_top')  
+                single_integral_function_top.x.array[:] = top_integrated_func_along_x   
+
+                single_integral_function_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name='single_integral_function_bot')  
+                single_integral_function_bot.x.array[:] = bot_integrated_func_along_x 
+
+                torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_top * dx_top_whole))/self.pvarry.panel_span*params.pv_array.modules_per_span
+                torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_bot * dx_bot_whole))/self.pvarry.panel_span*params.pv_array.modules_per_span
+                attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{module_id:.0f}"
+                setattr(self, attr_name, torque_double_integral)
+           
     def _assemble_system(self, params):
         """Pre-assemble all LHS matrices and RHS vectors
 
diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py
index 13f1764..ebd8acd 100644
--- a/pvade/geometry/MeshManager.py
+++ b/pvade/geometry/MeshManager.py
@@ -309,7 +309,7 @@ def _transfer_mesh_tags_to_submeshes(self, params):
         num_facets = f_map.size_local + f_map.num_ghosts
         all_values = np.zeros(num_facets, dtype=np.int32)
 
-        # Assign non-zero facet tags using the facet tag indices
+        # Assign non-zero facet tags using the facet tag indices, save the facet marker to all_values
         all_values[self.facet_tags.indices] = self.facet_tags.values
 
         cell_to_facet = self.msh.topology.connectivity(self.ndim, facet_dim)
diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py
index 31f4d42..8ad4a7e 100644
--- a/pvade/geometry/panels3d/DomainCreation.py
+++ b/pvade/geometry/panels3d/DomainCreation.py
@@ -149,14 +149,20 @@ def Rz(theta):
         domain_tag_list = []
         domain_tag_list.append(domain_tag)
 
+        # panel_tag_list includes all structure, panel+connector
         panel_tag_list = []
         panel_ct = 0
+        
+        # only include panels
+        structure_panel_only_list = []
+        # only include connectors
+        structure_connector_only_list = []
 
         module_distances = np.linspace(
             -half_span, half_span, params.pv_array.modules_per_span + 1
         )
 
-        transformed_com = {}
+        transformed_com = {} # all panels
 
         if (
             params.pv_array.torque_tube_separation > 0.0
@@ -166,6 +172,8 @@ def Rz(theta):
         else:
             modeling_torque_tube = False
 
+
+        # start to add panel
         for panel_id_y, yy in enumerate(y_centers):
             for panel_id_x, xx in enumerate(x_centers):
                 if isinstance(params.pv_array.tracker_angle, list):
@@ -181,8 +189,11 @@ def Rz(theta):
                 else:
                     tracker_angle_rad = np.radians(params.pv_array.tracker_angle)
 
-                this_panel_tag_list = []
-                this_panel_transformed_com = {}
+                this_panel_tag_list = [] # for each row
+                this_panel_transformed_com = {} # for each row
+                
+                numpy_pt_list = [] # for each row
+                embedded_lines_tag_list = [] # for each row
 
                 for module_id in range(params.pv_array.modules_per_span):
 
@@ -200,22 +211,33 @@ def Rz(theta):
                             module_span,
                             params.pv_array.panel_thickness,
                         )
+                       
+                        panel_tag_list.append((self.ndim, this_module)) # all panels
+                        this_panel_tag_list.append((self.ndim, this_module)) # this panel array/row
 
+                        structure_panel_only_list.append((self.ndim, this_module))
+                        
                         this_standoff = self.gmsh_model.occ.addBox(
-                            -params.pv_array.torque_tube_radius,
+                            -params.pv_array.block_chord_div_by_panel_chord * half_chord,
                             module_distances[module_id],
                             0.0,
-                            2.0 * params.pv_array.torque_tube_radius,
-                            module_span,
+                            2.0 * params.pv_array.block_chord_div_by_panel_chord * half_chord,
+                            params.pv_array.block_chord_div_by_panel_chord * half_chord,
                             params.pv_array.torque_tube_separation,
                         )
-
-                        panel_tag_list.append((self.ndim, this_module))
                         panel_tag_list.append((self.ndim, this_standoff))
-
-                        this_panel_tag_list.append((self.ndim, this_module))
                         this_panel_tag_list.append((self.ndim, this_standoff))
 
+                        structure_connector_only_list.append((self.ndim, this_standoff))
+
+                        # Add a bisecting line to the bottom of the connector in the spanwise direction
+                        pt_1 = self.gmsh_model.occ.addPoint(0, module_distances[module_id], 0.0)
+                        pt_2 = self.gmsh_model.occ.addPoint(0, module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0)
+                        numpy_pt_list.append([0, module_distances[module_id], 0.0, 0, module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0]) # for this row
+                        torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
+                        torque_tube_tag = (1, torque_tube_id)
+                        embedded_lines_tag_list.append(torque_tube_tag) # for this panel row
+
                     else:
                         this_module = self.gmsh_model.occ.addBox(
                             -half_chord,
@@ -225,118 +247,134 @@ def Rz(theta):
                             module_span,
                             params.pv_array.panel_thickness,
                         )
-
                         panel_tag_list.append((self.ndim, this_module))
-
                         this_panel_tag_list.append((self.ndim, this_module))
 
-                numpy_pt_list = []
-                embedded_lines_tag_list = []
-
-                # Add a bisecting line to the bottom of the panel in the spanwise direction
-                pt_1 = self.gmsh_model.occ.addPoint(0, -half_span, 0.0)
-                pt_2 = self.gmsh_model.occ.addPoint(0, half_span, 0.0)
-
-                numpy_pt_list.append([0, -half_span, 0.0, 0, half_span, 0.0])
-
-                torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
-                torque_tube_tag = (1, torque_tube_id)
-                embedded_lines_tag_list.append(torque_tube_tag)
-
-                # Add lines in the streamwise direction to mimic sections of panel held rigid by motor
-                if params.pv_array.span_fixation_pts is not None:
-                    if not isinstance(params.pv_array.span_fixation_pts, list):
-                        num_fixation_pts = int(
-                            np.floor(
-                                params.pv_array.panel_span
-                                / params.pv_array.span_fixation_pts
-                            )
-                        )
-
-                        fixation_pts_list = []
-
-                        for k in range(1, num_fixation_pts + 1):
-                            next_pt = k * params.pv_array.span_fixation_pts
-
-                            eps = 1e-9
-
-                            if (
-                                next_pt > eps
-                                and next_pt < params.pv_array.panel_span - eps
-                            ):
-                                fixation_pts_list.append(next_pt)
-
-                    else:
-                        fixation_pts_list = params.pv_array.span_fixation_pts
-
-                    for fp in fixation_pts_list:
-                        # FIXME: don't add the fixation points into the numpy tagging for now
-                        if modeling_torque_tube:
-                            numpy_pt_list.append(
-                                [
-                                    -params.pv_array.torque_tube_radius,
-                                    -half_span + fp,
-                                    0.0,
-                                    params.pv_array.torque_tube_radius,
-                                    -half_span + fp,
-                                    0.0,
-                                ]
-                            )
-
-                        else:
-                            numpy_pt_list.append(
-                                [
-                                    -half_chord,
-                                    -half_span + fp,
-                                    0.0,
-                                    half_chord,
-                                    -half_span + fp,
-                                    0.0,
-                                ]
-                            )
-
-                        # If the separation line already exists at a module division,
-                        # there's no need to redraw it, but otherwise, add the gmsh points
-                        # to the model for embedding
-                        if not np.any(np.isclose(-half_span + fp, module_distances)):
-                            print(
-                                f"Embedding a line for a no-deformation boundary condition."
-                            )
-
-                            if modeling_torque_tube:
-                                pt_1 = self.gmsh_model.occ.addPoint(
-                                    -params.pv_array.torque_tube_radius,
-                                    -half_span + fp,
-                                    0.0,
-                                )
-                                pt_2 = self.gmsh_model.occ.addPoint(
-                                    params.pv_array.torque_tube_radius,
-                                    -half_span + fp,
+                        structure_panel_only_list.append((self.ndim, this_module))
+
+                        pt_1 = self.gmsh_model.occ.addPoint(0, module_distances[module_id], 0.0)
+                        pt_2 = self.gmsh_model.occ.addPoint(0, module_distances[module_id+1], 0.0)
+                        numpy_pt_list.append([0, module_distances[module_id], 0.0, 0, module_distances[module_id+1], 0.0])
+                        torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
+                        torque_tube_tag = (1, torque_tube_id)
+                        embedded_lines_tag_list.append(torque_tube_tag) # for this panel row
+                
+                if modeling_torque_tube: # add the last connector
+                    last_standoff = self.gmsh_model.occ.addBox(
+                                    -params.pv_array.block_chord_div_by_panel_chord * half_chord,
+                                    module_distances[module_id+1] - params.pv_array.block_chord_div_by_panel_chord * half_chord,
                                     0.0,
+                                    2.0 * params.pv_array.block_chord_div_by_panel_chord * half_chord,
+                                    params.pv_array.block_chord_div_by_panel_chord * half_chord,
+                                    params.pv_array.torque_tube_separation,
                                 )
-
-                            else:
-                                pt_1 = self.gmsh_model.occ.addPoint(
-                                    -half_chord, -half_span + fp, 0.0
-                                )
-                                pt_2 = self.gmsh_model.occ.addPoint(
-                                    half_chord, -half_span + fp, 0.0
-                                )
-
-                            fixation_line_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
-                            fixation_line_tag = (1, fixation_line_id)
-                            embedded_lines_tag_list.append(fixation_line_tag)
-                        else:
-                            print(
-                                f"Applying no-deformation boundary condition at {fp}."
-                            )
+                    this_panel_tag_list.append((self.ndim, last_standoff))
+                    panel_tag_list.append((self.ndim, last_standoff))
+
+                    structure_connector_only_list.append((self.ndim, last_standoff))
+
+                    pt_1 = self.gmsh_model.occ.addPoint(0, module_distances[module_id+1] - params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0)
+                    pt_2 = self.gmsh_model.occ.addPoint(0, module_distances[module_id+1], 0.0)
+                    numpy_pt_list.append([0, module_distances[module_id+1] - params.pv_array.block_chord_div_by_panel_chord * half_chord, 0.0, 0, module_distances[module_id+1], 0.0])
+                    torque_tube_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
+                    torque_tube_tag = (1, torque_tube_id)
+                    embedded_lines_tag_list.append(torque_tube_tag)
+                
+                """ this is not used as we are using block surface to simulate the fixed bc applied by motor """
+                # # Add lines in the streamwise direction to mimic sections of panel held rigid by motor
+                # if params.pv_array.span_fixation_pts is not None:
+                #     if not isinstance(params.pv_array.span_fixation_pts, list):
+                #         num_fixation_pts = int(
+                #             np.floor(
+                #                 params.pv_array.panel_span
+                #                 / params.pv_array.span_fixation_pts
+                #             )
+                #         )
+
+                #         fixation_pts_list = []
+
+                #         for k in range(1, num_fixation_pts + 1):
+                #             next_pt = k * params.pv_array.span_fixation_pts
+
+                #             eps = 1e-9
+
+                #             if (
+                #                 next_pt > eps
+                #                 and next_pt < params.pv_array.panel_span - eps
+                #             ):
+                #                 fixation_pts_list.append(next_pt)
+
+                #     else:
+                #         fixation_pts_list = params.pv_array.span_fixation_pts
+
+                #     for fp in fixation_pts_list:
+                #         # FIXME: don't add the fixation points into the numpy tagging for now
+                #         if modeling_torque_tube:
+                #             numpy_pt_list.append(
+                #                 [
+                #                     -params.pv_array.torque_tube_radius,
+                #                     -half_span + fp,
+                #                     0.0,
+                #                     params.pv_array.torque_tube_radius,
+                #                     -half_span + fp,
+                #                     0.0,
+                #                 ]
+                #             )
+
+                #         else:
+                #             numpy_pt_list.append(
+                #                 [
+                #                     -half_chord,
+                #                     -half_span + fp,
+                #                     0.0,
+                #                     half_chord,
+                #                     -half_span + fp,
+                #                     0.0,
+                #                 ]
+                #             )
+
+                #         # If the separation line already exists at a module division,
+                #         # there's no need to redraw it, but otherwise, add the gmsh points
+                #         # to the model for embedding
+                #         if not np.any(np.isclose(-half_span + fp, module_distances)):
+                #             print(
+                #                 f"Embedding a line for a no-deformation boundary condition."
+                #             )
+
+                #             if modeling_torque_tube:
+                #                 pt_1 = self.gmsh_model.occ.addPoint(
+                #                     -params.pv_array.torque_tube_radius,
+                #                     -half_span + fp,
+                #                     0.0,
+                #                 )
+                #                 pt_2 = self.gmsh_model.occ.addPoint(
+                #                     params.pv_array.torque_tube_radius,
+                #                     -half_span + fp,
+                #                     0.0,
+                #                 )
+
+                #             else:
+                #                 pt_1 = self.gmsh_model.occ.addPoint(
+                #                     -half_chord, -half_span + fp, 0.0
+                #                 )
+                #                 pt_2 = self.gmsh_model.occ.addPoint(
+                #                     half_chord, -half_span + fp, 0.0
+                #                 )
+
+                #             fixation_line_id = self.gmsh_model.occ.addLine(pt_1, pt_2)
+                #             fixation_line_tag = (1, fixation_line_id)
+                #             embedded_lines_tag_list.append(fixation_line_tag)
+                #         else:
+                #             print(
+                #                 f"Applying no-deformation boundary condition at {fp}."
+                #             )
 
                 # Fragment the lines into the surfaces (equivalent to embedding these lines)
                 self.gmsh_model.occ.fragment(
                     this_panel_tag_list, embedded_lines_tag_list
                 )
 
-                for panel_tag in this_panel_tag_list:
+                for panel_tag in this_panel_tag_list: # 3d cell domain
                     self.gmsh_model.occ.synchronize()
 
                     # Get the list of 2D surfaces (surfaces) that make up this panel
@@ -351,61 +389,189 @@ def Rz(theta):
 
                         target_key = None
 
+                        surface_located_or_not = False
+
                         # sturctures tagging
                         if (
-                            np.isclose(com[0], -half_chord)
-                            or np.isclose(com[0], -params.pv_array.torque_tube_radius)
-                            and modeling_torque_tube
+                            np.isclose(com[0], -half_chord) 
                         ):
-                            target_key = f"left_{panel_ct:.0f}"
+                            target_key = f"panel_left_{panel_ct:.0f}"
+                            surface_located_or_not = True
 
-                        elif (
+                        if (
                             np.isclose(com[0], half_chord)
-                            or np.isclose(com[0], params.pv_array.torque_tube_radius)
-                            and modeling_torque_tube
                         ):
-                            target_key = f"right_{panel_ct:.0f}"
-
-                        elif np.isclose(com[1], -half_span):
-                            target_key = f"front_{panel_ct:.0f}"
+                            target_key = f"panel_right_{panel_ct:.0f}"
+                            surface_located_or_not = True
 
-                        elif np.isclose(com[1], half_span):
-                            target_key = f"back_{panel_ct:.0f}"
+                        if (
+                            np.isclose(com[1], module_distances[0])
+                            and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0)
+                        ):
+                            target_key = f"panel_front_{panel_ct:.0f}"
+                            surface_located_or_not = True
 
-                        elif (
-                            np.isclose(com[2], params.pv_array.torque_tube_separation)
-                            and not np.isclose(com[0], 0.0)
-                            and modeling_torque_tube
-                            or np.isclose(com[2], 0.0)
-                            and not modeling_torque_tube
+                        if (
+                            np.isclose(com[1], module_distances[params.pv_array.modules_per_span])
+                            and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0)
                         ):
-                            target_key = f"bottom_{panel_ct:.0f}"
+                            target_key = f"panel_back_{panel_ct:.0f}"
+                            surface_located_or_not = True
 
-                        elif np.isclose(com[2], 0.0) and modeling_torque_tube:
-                            target_key = f"torque_tube_{panel_ct:.0f}"
+                        for module_id in range(params.pv_array.modules_per_span):
+                            if (
+                                com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord
+                                and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord
+                                and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation)
+                            ):
+                                target_key = f"panel_bottom_{panel_ct:.0f}_{module_id:.0f}"
+                                surface_located_or_not = True
+                                break
 
-                        elif (
-                            np.isclose(
-                                com[2],
-                                params.pv_array.torque_tube_separation
-                                + params.pv_array.panel_thickness,
-                            )
-                            and modeling_torque_tube
-                            or np.isclose(com[2], params.pv_array.panel_thickness)
-                            and not modeling_torque_tube
+                            if (
+                                com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord
+                                and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord
+                                and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness)
+                            ):
+                                target_key = f"panel_top_{panel_ct:.0f}_{module_id:.0f}"
+                                surface_located_or_not = True
+                                break
+                        # if (
+                        #         np.isclose(
+                        #             com[2],
+                        #             params.pv_array.torque_tube_separation
+                        #             + params.pv_array.panel_thickness
+                        #         )
+                                
+                        #     ):
+                        #         target_key = f"panel_top_{panel_ct:.0f}"
+                        #         surface_located_or_not = True
+                        
+                        # mark the last block in each row
+                        if (
+                                    np.isclose(com[0], params.pv_array.block_chord_div_by_panel_chord * half_chord)
+                                    and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                                ):
+                            target_key = f"block_right_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}"
+                            surface_located_or_not = True
+
+                        if (
+                            np.isclose(com[0], -params.pv_array.block_chord_div_by_panel_chord * half_chord)
+                            and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
                         ):
-                            target_key = f"top_{panel_ct:.0f}"
+                            target_key = f"block_left_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}"
+                            surface_located_or_not = True
 
-                        else:
-                            target_key = f"trash_{panel_ct:.0f}"
+                        if (
+                            np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord)
+                        ):
+                            target_key = f"block_front_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}"
+                            surface_located_or_not = True
+                        
+                        if (
+                            np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span])
+                        ):
+                            target_key = f"block_back_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}"
+                            surface_located_or_not = True
+                        
+                        if (
+                            np.isclose(com[2], 0.0) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                        ):
+                            target_key = f"block_bottom_{panel_ct:.0f}_{params.pv_array.modules_per_span:.0f}"
+                            surface_located_or_not = True
 
+                        if (
+                            np.isclose(com[2], params.pv_array.torque_tube_separation) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                        ):
+                            target_key = f"interior_surface_{panel_ct:.0f}"
+                            surface_located_or_not = True
+
+                        if not surface_located_or_not:
+                            for module_id in range(params.pv_array.modules_per_span):
+
+                                # sturctures tagging
+
+                                if (
+                                    np.isclose(com[1], module_distances[module_id])
+                                    and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0)
+                                ):
+                                    target_key = f"interior_surface_{panel_ct:.0f}"
+                                    surface_located_or_not = True
+                                    break
+
+                                if (
+                                    np.isclose(com[1], module_distances[module_id+1])
+                                    and np.isclose(com[0], 0) and np.isclose(com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness/2.0)
+                                ):
+                                    target_key = f"interior_surface_{panel_ct:.0f}"
+                                    surface_located_or_not = True
+                                    break
+
+                                # if (
+                                #     np.isclose(com[2], params.pv_array.torque_tube_separation)
+                                #     and np.isclose(com[0], 0.0) and com[1] >= module_distances[module_id]+module_span/2.0-params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0
+                                #     and com[1] <= module_distances[module_id]+module_span/2.0+params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0
+                                # ):
+                                #     target_key = f"panel_bottom_{panel_ct:.0f}"
+                                #     surface_located_or_not = True
+                                #     break
+                                
+                                if (
+                                    (np.isclose(com[2], params.pv_array.torque_tube_separation)
+                                    and np.isclose(com[0], 0.0) and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                                    and modeling_torque_tube)
+                                ):
+                                    target_key = f"interior_surface_{panel_ct:.0f}"
+                                    surface_located_or_not = True
+                                    break
+                                            
+                                if (
+                                    np.isclose(com[0], params.pv_array.block_chord_div_by_panel_chord * half_chord)
+                                    and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                                ):
+                                    target_key = f"block_right_{panel_ct:.0f}_{module_id:.0f}"
+                                    surface_located_or_not = True
+                                    break
+
+                                if (
+                                    np.isclose(com[0], -params.pv_array.block_chord_div_by_panel_chord * half_chord)
+                                    and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                                ):
+                                    target_key = f"block_left_{panel_ct:.0f}_{module_id:.0f}"
+                                    surface_located_or_not = True
+                                    break
+
+                                if (
+                                    np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[module_id])
+                                ):
+                                    target_key = f"block_front_{panel_ct:.0f}_{module_id:.0f}"
+                                    surface_located_or_not = True
+                                    break
+                                
+                                if (
+                                    np.isclose(com[2], params.pv_array.torque_tube_separation/2.0) and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord)
+                                ):
+                                    target_key = f"block_back_{panel_ct:.0f}_{module_id:.0f}"
+                                    surface_located_or_not = True
+                                    break
+                                
+                                if (
+                                    np.isclose(com[2], 0.0) and modeling_torque_tube and np.isclose(com[1], module_distances[module_id] + params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
+                                ):
+                                    target_key = f"block_bottom_{panel_ct:.0f}_{module_id:.0f}"
+                                    surface_located_or_not = True
+                                    break
+
+                        if not surface_located_or_not:
+                            target_key = f"trash_{panel_ct:.0f}"
+                        
                         if target_key is not None:
                             if target_key in this_panel_transformed_com:
                                 this_panel_transformed_com[target_key].append(com)
                             else:
                                 this_panel_transformed_com[target_key] = [com]
 
-                for key, val in this_panel_transformed_com.items():
+                for key, val in this_panel_transformed_com.items(): # all facets for each panel row.
                     for row_num, com in enumerate(val):
                         com_array = np.array(com)
 
@@ -426,8 +592,9 @@ def Rz(theta):
                         if key in transformed_com:
                             transformed_com[key].append(com_array)
                         else:
-                            transformed_com[key] = [com_array]
+                            transformed_com[key] = [com_array] # including all rows
 
+                # actually this is panel array count
                 panel_ct += 1
 
                 # Rotate the panel by its tracking angle along the y-axis
@@ -505,6 +672,10 @@ def Rz(theta):
             self.numpy_pt_total_array, (-1, int(2 * self.ndim))
         )
 
+        # print(transformed_com)
+        # exit()
+        
+        # for all panels
         # Loop over all the finalized surfaces after fragmentation and tag everything
         all_surf_tag_list = self.gmsh_model.occ.getEntities(self.ndim - 1)
 
@@ -590,6 +761,9 @@ def Rz(theta):
                 # Fluid Cell
                 fluid_vol_list.append(vol_id)
 
+        structure_panel_only_list = []
+        structure_connector_only_list = []
+        
         self._add_to_domain_markers("structure", structure_vol_list, "cell")
         self._add_to_domain_markers("fluid", fluid_vol_list, "cell")
 
@@ -1263,22 +1437,35 @@ def set_length_scales(self, params, domain_markers):
         internal_surface_tags = []
 
         for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows):
+            
             internal_surface_tags.extend(
-                domain_markers[f"bottom_{panel_id}"]["gmsh_tags"]
+                domain_markers[f"panel_left_{panel_id}"]["gmsh_tags"]
             )
-            internal_surface_tags.extend(domain_markers[f"top_{panel_id}"]["gmsh_tags"])
             internal_surface_tags.extend(
-                domain_markers[f"left_{panel_id}"]["gmsh_tags"]
+                domain_markers[f"panel_right_{panel_id}"]["gmsh_tags"]
             )
             internal_surface_tags.extend(
-                domain_markers[f"right_{panel_id}"]["gmsh_tags"]
+                domain_markers[f"panel_front_{panel_id}"]["gmsh_tags"]
             )
             internal_surface_tags.extend(
-                domain_markers[f"front_{panel_id}"]["gmsh_tags"]
+                domain_markers[f"panel_back_{panel_id}"]["gmsh_tags"]
             )
-            internal_surface_tags.extend(
-                domain_markers[f"back_{panel_id}"]["gmsh_tags"]
+
+            for module_id in range(params.pv_array.modules_per_span):
+                internal_surface_tags.extend(
+                    domain_markers[
+                        f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"
+                    ]["gmsh_tags"]
+            )
+                internal_surface_tags.extend(
+                    domain_markers[
+                        f"panel_top_{panel_id:.0f}_{module_id:.0f}"
+                    ]["gmsh_tags"]
             )
+        
+        # print(internal_surface_tags)
+        # print(len(internal_surface_tags))
+        # exit()
 
         min_dist = []
 
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 85ea43e..862dbeb 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -87,6 +87,143 @@ def build_boundary_conditions(self, domain, params):
 
         self.bc = build_structure_boundary_conditions(domain, params, self.V)
 
+    def calculate_K_for_Robin_BC(self, domain, flow, params):
+
+        # shear modulus of tube
+        Gt = params.structure.elasticity_modulus_tube / (2 * (1 + params.structure.poissons_ratio_tube))
+
+        Ipt = np.pi/2*(params.pv_array.torque_tube_outer_radius**4 - params.pv_array.torque_tube_inner_radius**4)
+
+        Gp = params.structure.elasticity_modulus / (2 * (1 + params.structure.poissons_ratio))
+
+        Ipp = 1/3.0*params.pv_array.panel_thickness * params.pv_array.panel_chord**3
+
+        # total number of pv rows
+        total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows
+
+        num_panel_right_fixed = params.pv_array.modules_per_span // 2
+        num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed
+
+        # for right fixed part:
+
+        for panel_id in range(total_num_panels):
+            total_torque_right_fixed = 0
+            total_torque_left_fixed = 0
+            
+            for i in range(num_panel_left_fixed, params.pv_array.modules_per_span):
+                name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}"
+                total_torque_right_fixed += getattr(flow, name)
+            
+            for i in range(num_panel_left_fixed):
+                name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}"
+                total_torque_left_fixed += getattr(flow, name)
+
+            theo_matrix_right_fixed = np.zeros((num_panel_right_fixed+1, num_panel_right_fixed+1))
+            theo_vector_right_fixed = np.zeros(num_panel_right_fixed+1)
+
+            theo_matrix_right_fixed[0, :] = 1.0
+            theo_vector_right_fixed[0] = total_torque_right_fixed
+            theo_matrix_right_fixed[1:, :-1] = params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
+            
+            # as the double integral of torque along the span from front to back, left fixed part to right fixed part, it need to be flipped
+            T_double_integral_right_fixed = []
+            T_right_fixed = []
+            for i in range(params.pv_array.modules_per_span):
+                name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}"
+                T_double_integral_right_fixed.append(getattr(flow, name))
+                name = f"total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}"
+                T_right_fixed.append(getattr(flow, name))
+            for i in range(num_panel_right_fixed):
+                for j in range(i+1):
+                    theo_vector_right_fixed[i+1] += T_double_integral_right_fixed[-1-j]/(Gp*Ipp)
+                    theo_vector_right_fixed[i+1] -= (i+1-j)*T_right_fixed[-1-j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span
+                for j in range(i):
+                    theo_matrix_right_fixed[i+1, :-1-j-1] += params.pv_array.panel_span/params.pv_array.modules_per_span/(Gt*Ipt)
+
+
+            for i in range(num_panel_right_fixed):
+                for j in range(i+1):
+                        theo_matrix_right_fixed[i+1, num_panel_right_fixed-j:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gp/Ipp
+
+            R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_right_fixed), theo_vector_right_fixed)
+            
+            tube_rotate_matrix = np.ones((num_panel_right_fixed, num_panel_right_fixed))*params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt)
+
+            for i in range(num_panel_right_fixed):
+                for j in range(i):
+                    tube_rotate_matrix[i, :num_panel_right_fixed-1-j] += params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt)
+            
+            # check the standalone code, why it need to be flipped.
+            phi = np.flip(np.dot(tube_rotate_matrix, R_reaction_torque[:-1]))
+
+            # rotation of connector
+            x_panel = np.arange(num_panel_right_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span) # location of panels
+
+            block_length = params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord
+            block_width = params.pv_array.block_chord_div_by_panel_chord*params.pv_array.panel_span/params.pv_array.modules_per_span/2
+
+            # To Do: check the angle, is this correct?
+            array_rotation = (params.fluid.wind_direction + 90.0) % 360.0
+            array_rotation_rad = np.radians(array_rotation)
+
+            for i in range(num_panel_right_fixed):
+
+                phi_i = phi[i]
+
+                K_i = 12*(R_reaction_torque[i])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width
+
+                name_K = f"spring_stiffness_{panel_id:.0f}_{params.pv_array.modules_per_span+1-i:.0f}"
+                
+                setattr(self, name_K, K_i)
+
+            # for left fixed part:
+            theo_matrix_left_fixed = np.zeros((num_panel_left_fixed+1, num_panel_left_fixed+1))
+            theo_vector_left_fixed = np.zeros(num_panel_left_fixed+1)
+            theo_matrix_left_fixed[0, :] = 1.0
+            theo_vector_left_fixed[0] = total_torque_left_fixed
+
+            T_double_integral_left_fixed = []
+            T_left_fixed = []
+
+            for i in range(num_panel_left_fixed):
+                name = f"double_integral_total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
+                T_double_integral_left_fixed.append(getattr(flow, name))
+                name = f"total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
+                T_left_fixed.append(getattr(flow, name))
+
+            for i in range(num_panel_left_fixed):
+                for j in range(i+1):
+                    theo_matrix_left_fixed[i+1, j] = (i+1-j)*params.pv_array.panel_span/params.pv_array.modules_per_span/(Gp*Ipp)
+                    theo_vector_left_fixed[i+1] += T_double_integral_left_fixed[j]/(Gp*Ipp)
+                for j in range(i):
+                    theo_vector_left_fixed[i+1] += (i-j)*T_left_fixed[j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span
+
+            for i in range(num_panel_left_fixed):    
+                theo_matrix_left_fixed[i+1:, i+1:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
+
+            R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_left_fixed), theo_vector_left_fixed)  # this is the torque applied to tube
+
+
+            tube_rotate_matrix = np.zeros((num_panel_left_fixed, num_panel_left_fixed))
+            for i in range(num_panel_left_fixed):    
+                tube_rotate_matrix[i:, i:] += params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
+
+            phi = np.dot(tube_rotate_matrix, R_reaction_torque[1:])
+
+            # rotation of connector
+            x_panel = np.arange(num_panel_left_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span)+params.pv_array.panel_span/params.pv_array.modules_per_span # location of panles
+
+
+            for i in range(num_panel_left_fixed):
+
+                phi_i = phi[i]
+
+                K_i = 12*(R_reaction_torque[i+1])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width
+
+                name_K = f"spring_stiffness_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
+                
+                setattr(self, name_K, K_i)
+
     def update_a(self, u, u_old, v_old, a_old, dt, beta, ufl=True):
         # Update formula for acceleration
         # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
@@ -126,7 +263,7 @@ def update_fields(self, u, u_old, v_old, a_old, dt, beta, gamma):
     def avg(self, x_old, x_new, alpha):
         return alpha * x_old + (1 - alpha) * x_new
 
-    def build_forms(self, domain, params, structure):
+    def build_forms(self, domain, params, structure, flow):
         """Builds all variational statements
 
         This method creates all the functions, expressions, and variational
@@ -212,6 +349,9 @@ def c(u, u_):
 
         def k_nominal(u, u_):
             return ufl.inner(P_(u), ufl.grad(u_))
+        
+        def k_nominal_connector(u, u_):
+            return ufl.inner(P_connector(u), ufl.grad(u_))
 
         # The deformation gradient, F = I + dy/dX
         def F_(u):
@@ -241,6 +381,17 @@ def S_(u):
 
             S_svk = structure.lame_lambda * ufl.tr(E) * I + 2.0 * structure.lame_mu * E
             return S_svk
+        
+        # The second Piola–Kirchhoff stress, S
+        def S_connector(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 = structure.lame_lambda_connector * ufl.tr(E) * I + 2.0 * structure.lame_mu_connector * E
+            return S_svk
 
         # The first Piola–Kirchhoff stress tensor, P = F*S
         def P_(u):
@@ -248,6 +399,12 @@ def P_(u):
             S = S_(u)
             # return ufl.inv(F) * S
             return F * S
+        
+        def P_connector(u):
+            F = F_(u)
+            S = S_connector(u)
+            # return ufl.inv(F) * S
+            return F * S
 
         # self.uh_exp = dolfinx.fem.Function(self.V,  name="Deformation")
 
@@ -302,15 +459,32 @@ def P_(u):
 
         F = ufl.grad(self.u) + ufl.Identity(len(self.u))
         J = ufl.det(F)
+
+        self.z_unit_vector = dolfinx.fem.Constant(domain.structure.msh, [0.0,0.0,1.0])  # surface traction, N/m^2
+        
+        self.calculate_K_for_Robin_BC(domain, flow, params)
+
+        # To Do: how to differentiate the connector part and the panel part, dx_connector and dx_panel
         self.res = (
             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
+            + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * ufl.dx_panel
+            + k_nominal_connector(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * ufl.dx_connector
+            - structure.rho * ufl.inner(self.f, self.u_) * ufl.dx_panel
+            - structure.rho_connector * ufl.inner(self.f, self.u_) * ufl.dx_connector
             - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_)
             * self.ds
         )  # - Wext(self.u)
 
+        # Robin boundary condition terms
+        for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows):
+            for i in range(params.pv_array.modules_per_span):
+                name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
+                K_springs = dolfinx.fem.Constant(domain.structure.msh, float(getattr(self, name_K)))
+                self.res += ufl.dot(K_springs * self.u_, self.z_unit_vector)*ds_bottom.panel_connector_markers[panel_id][i]
+        for i in range(params.pv_array.modules_per_span):
+            self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*ds_bottom....
+
         # self.a = dolfinx.fem.form(ufl.lhs(res))
         # self.L = dolfinx.fem.form(ufl.rhs(res))
 
diff --git a/pvade/structure/StructureMain.py b/pvade/structure/StructureMain.py
index ed3a756..8f81ea2 100644
--- a/pvade/structure/StructureMain.py
+++ b/pvade/structure/StructureMain.py
@@ -80,6 +80,10 @@ def __init__(self, domain, params):
             domain.structure.msh, params.structure.rho
         )  # Constant(0.)
 
+        self.rho_connector = dolfinx.fem.Constant(
+            domain.structure.msh, params.structure.rho_tube
+        )
+
         # Define structural properties
         self.E = params.structure.elasticity_modulus  # 1.0e9
         self.poissons_ratio = params.structure.poissons_ratio  # 0.3
@@ -89,6 +93,16 @@ def __init__(self, domain, params):
             * self.poissons_ratio
             / ((1.0 + self.poissons_ratio) * (1.0 - 2.0 * self.poissons_ratio))
         )
+        
+        # we are assuming the connectors has same mechanical properties as the tubes
+        self.E_connector = params.structure.elasticity_modulus_tube  # 1.0e9
+        self.poissons_ratio_connector = params.structure.poissons_ratio_tube  # 0.3
+        self.lame_mu_connector = self.E_connector / (2.0 * (1.0 + self.poissons_ratio_connector))
+        self.lame_lambda_connector = (
+            self.E_connector
+            * self.poissons_ratio_connector
+            / ((1.0 + self.poissons_ratio_connector) * (1.0 - 2.0 * self.poissons_ratio_connector))
+        )
 
         if self.rank == 0:
             print(
@@ -201,7 +215,7 @@ def update_fields(self, u, u_old, v_old, a_old, dt, beta, gamma):
     def avg(self, x_old, x_new, alpha):
         return alpha * x_old + (1 - alpha) * x_new
 
-    def build_forms(self, domain, params):
+    def build_forms(self, domain, params, flow):
         """Builds all variational statements
 
         This method creates all the functions, expressions, and variational
@@ -217,7 +231,7 @@ def build_forms(self, domain, params):
             params (:obj:`pvade.Parameters.SimParams`): A SimParams object
 
         """
-        self.elasticity.build_forms(domain, params, self)
+        self.elasticity.build_forms(domain, params, self, flow)
 
     def _assemble_system(self, params):
         """Pre-assemble all LHS matrices and RHS vectors
diff --git a/pvade/structure/boundary_conditions.py b/pvade/structure/boundary_conditions.py
index 0aa0ea0..ebfe5f7 100644
--- a/pvade/structure/boundary_conditions.py
+++ b/pvade/structure/boundary_conditions.py
@@ -344,20 +344,12 @@ def connection_point_up_helper(nodes_to_pin_between):
 
         return fn_handle
 
-    # Start pinning along the lines expressed byt numpy_pt_total_array
-    # First, determine the total number of rows (number of lines to pin)
-    num_nodes = np.shape(domain.numpy_pt_total_array)[0]
-
-    # Determine how many pinning lines exist per each panel (e.g., 24 lines distributed on 8 panels means 3 lines per panel)
-    nodes_per_panel = int(num_nodes / total_num_panels)
-
-    # The torque tube entry (oriented spanwise along the middle, divides panel into upstream and downstream rectangular halves)
-    # is always the first entry, e.g., [0, ..., ..., 3, ..., ..., 6, ..., ...], [0, 3, 6] are the torque tubes
-    tube_nodes_idx = np.arange(0, num_nodes, nodes_per_panel, dtype=np.int64)
+    # # Start pinning along the lines expressed byt numpy_pt_total_array
+    # The center line of connectors bottom surface is fixed to remove rigid body motion
 
     if params.structure.tube_connection == True:
         # If making torque tube connections, pass only those pinning lines to the BC identification function
-        tube_nodes = domain.numpy_pt_total_array[tube_nodes_idx, :]
+        tube_nodes = domain.numpy_pt_total_array[:, :]
 
         facet_uppoint = dolfinx.mesh.locate_entities(
             domain.structure.msh, 1, connection_point_up_helper(tube_nodes)
@@ -369,18 +361,18 @@ def connection_point_up_helper(nodes_to_pin_between):
         bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace))
 
     if params.structure.motor_connection == True:
-        # If making motor mount connections, pass only those pinning lines to the BC identification function
-        # this is done by making a copy of the numpy_pt_total_array with the torque tube lines *deleted*
-        # not done in place, so numpy_pt_total_array remains unaltered.
-        motor_nodes = np.delete(domain.numpy_pt_total_array, tube_nodes_idx, axis=0)
+        # The bottom surface of the center connector is fixed to represent the motor mount
+        motor_location = params.pv_array.modules_per_span // 2
+        for panel_id in range(total_num_panels):
 
-        facet_uppoint = dolfinx.mesh.locate_entities(
-            domain.structure.msh, 1, connection_point_up_helper(motor_nodes)
-        )
-        dofs_disp = dolfinx.fem.locate_dofs_topological(
-            functionspace, 1, [facet_uppoint]
+            mount_facet = domain.structure.facet_tags.find(
+                domain.domain_markers[f"block_bottom_{panel_id:.0f}_{motor_location:.0f}"]["idx"]
         )
+        
+            dofs_disp = dolfinx.fem.locate_dofs_topological(
+                functionspace, 2, [mount_facet]
+            )
 
-        bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace))
+            bc.append(dolfinx.fem.dirichletbc(zero_vec, dofs_disp, functionspace))
 
     return bc
diff --git a/pvade_main.py b/pvade_main.py
index 8bc322a..6f262f5 100644
--- a/pvade_main.py
+++ b/pvade_main.py
@@ -16,6 +16,10 @@
 import os
 from mpi4py import MPI
 
+"""
+to run: python -u pvade_main.py --input_file=input/duramat_case_study.yaml --domain.l_char=3.0 --general.mesh_only=true --pv_array.torque_tube_radius=0.1 --pv_array.torque_tube_separation=0.4 --pv_array.stream_rows=2 --pv_array.tracker_angle=[-40,40] --pv_array.modules_per_span=10
+"""
+
 
 def main(input_file=None):
     # Get the path to the input file from the command line
@@ -40,7 +44,7 @@ def main(input_file=None):
         domain.read_mesh_files(params.general.input_mesh_dir, params)
     else:
         domain.build(params)
-
+    
     # If we only want to create the mesh, we can stop here
     if params.general.mesh_only:
         list_timings(params.comm, [TimingType.wall])
@@ -60,7 +64,7 @@ def main(input_file=None):
         structure = Structure(domain, params)
         structure.build_boundary_conditions(domain, params)
         # # # Build the fluid forms
-        structure.build_forms(domain, params)
+        structure.build_forms(domain, params, flow)
     else:
         structure = None
 

From 64f4efccae199a39bd9f65cd4da9f0369889d69d Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Tue, 11 Nov 2025 17:14:25 -0700
Subject: [PATCH 19/26] dx and ds for modules and connectors

---
 pvade/IO/input_schema.yaml                |  7 ++++
 pvade/fluid/FlowManager.py                |  3 ++
 pvade/geometry/MeshManager.py             | 42 +++++++++++++++++++---
 pvade/geometry/panels3d/DomainCreation.py | 43 ++++++++++++++++++++---
 pvade/structure/ElasticityAnalysis.py     | 18 +++++-----
 5 files changed, 96 insertions(+), 17 deletions(-)

diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index 9822740..f9911da 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -691,6 +691,13 @@ properties:
         type: "number"
         description: "poissons ratio of the torque tube."
         units: "None"
+      density_tube:
+        default: 1.0
+        # minimum: 0.001
+        # maximum: 1000.0
+        type: "number"
+        description: "The density of the structure."
+        units: "kg/meter^3"
       body_force_x: 
         default:  100
         type: "number"
diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index 0f2a1c8..f52ec07 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -856,6 +856,9 @@ def solve(self, domain, params, current_time):
 
         self.compute_lift_and_drag(params, current_time)
 
+        self.compute_panel_torques(domain, params)
+        self.compute_double_integral_panel_torques(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)
diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py
index ebd8acd..618f064 100644
--- a/pvade/geometry/MeshManager.py
+++ b/pvade/geometry/MeshManager.py
@@ -207,6 +207,8 @@ def build(self, params):
             gdim=self.ndim,
         )
 
+
+
         self.msh.topology.create_connectivity(self.ndim, self.ndim - 1)
 
         # Specify names for the mesh elements
@@ -214,6 +216,11 @@ def build(self, params):
         self.cell_tags.name = "cell_tags"
         self.facet_tags.name = "facet_tags"
 
+        with dolfinx.io.XDMFFile(self.comm, f"parent_mesh_only.xdmf", "w") as fp:
+            fp.write_mesh(self.msh)
+            fp.write_meshtags(self.cell_tags)
+            fp.write_meshtags(self.facet_tags)
+
         # if (
         #     params.general.geometry_module == "panels3d"
         #     or params.general.geometry_module == "heliostats3d"
@@ -276,10 +283,17 @@ def _create_submeshes_from_parent(self, params):
                 print(f"Creating {sub_domain_name} submesh")
 
             # Get the idx associated with either "fluid" or "structure"
-            marker_id = self.domain_markers[sub_domain_name]["idx"]
-
-            # Find all cells where cell tag = marker_id
-            submesh_cells = self.cell_tags.find(marker_id)
+            if sub_domain_name == "structure" and "structure" not in self.domain_markers:
+                marker_id = self.domain_markers["modules"]["idx"] 
+                # Find all cells where cell tag = marker_id
+                submesh_cells_modules = self.cell_tags.find(marker_id)
+                marker_id = self.domain_markers["connectors"]["idx"]
+                submesh_cells = np.hstack((self.cell_tags.find(marker_id), submesh_cells_modules))
+                
+            else:
+                marker_id = self.domain_markers[sub_domain_name]["idx"]
+                # Find all cells where cell tag = marker_id
+                submesh_cells = self.cell_tags.find(marker_id)
 
             # Use those found cells to construct a new mesh
             submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(
@@ -345,14 +359,32 @@ def _transfer_mesh_tags_to_submeshes(self, params):
                 for child, parent in zip(child_facets, parent_facets):
                     sub_values[child] = all_values[parent]
 
+            # sub_cell_map = sub_domain.msh.topology.index_map(self.ndim)
+            f_map_cell = self.msh.topology.index_map(self.ndim)
+
+            # Get the total number of cells in the parent mesh
+            num_cells = f_map_cell.size_local + f_map_cell.num_ghosts
+            all_cell_values = np.zeros(num_cells, dtype=np.int32)
+            all_cell_values[self.cell_tags.indices] = self.cell_tags.values
+
+
             sub_cell_map = sub_domain.msh.topology.index_map(self.ndim)
             sub_num_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts
 
+            sub_cell_values = np.empty(sub_num_cells, dtype=np.int32)
+
+
+            for k, entity in enumerate(sub_domain.entity_map):
+                sub_cell_values[k] = all_cell_values[entity]
+
+ 
+            # sub_num_cells = sub_cell_map.size_local + sub_cell_map.num_ghosts
+
             sub_domain.cell_tags = dolfinx.mesh.meshtags(
                 sub_domain.msh,
                 sub_domain.msh.topology.dim,
                 np.arange(sub_num_cells, dtype=np.int32),
-                np.ones(sub_num_cells, dtype=np.int32),
+                sub_cell_values,
             )
             sub_domain.cell_tags.name = "cell_tags"
 
diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py
index 8ad4a7e..4709c9b 100644
--- a/pvade/geometry/panels3d/DomainCreation.py
+++ b/pvade/geometry/panels3d/DomainCreation.py
@@ -166,13 +166,14 @@ def Rz(theta):
 
         if (
             params.pv_array.torque_tube_separation > 0.0
-            and params.pv_array.torque_tube_radius > 0.0
+            and params.pv_array.torque_tube_outer_radius > 0.0
         ):
             modeling_torque_tube = True
         else:
             modeling_torque_tube = False
 
-
+        vol_tags_modules = []
+        vol_tags_connectors = []
         # start to add panel
         for panel_id_y, yy in enumerate(y_centers):
             for panel_id_x, xx in enumerate(x_centers):
@@ -382,6 +383,20 @@ def Rz(theta):
                         [panel_tag], oriented=False
                     )
 
+                    vol_com = self.gmsh_model.occ.getCenterOfMass(self.ndim, panel_tag[1])
+                    if np.isclose(vol_com[2], params.pv_array.torque_tube_separation + params.pv_array.panel_thickness / 2.0):
+                        target_key = f"modules"
+                    elif np.isclose(vol_com[2], params.pv_array.torque_tube_separation / 2.0):
+                        target_key = f"connectors"
+                    
+                    
+                        
+                    if target_key is not None:
+                        if target_key in this_panel_transformed_com:
+                            this_panel_transformed_com[target_key].append(vol_com)
+                        else:
+                            this_panel_transformed_com[target_key] = [vol_com]
+
                     for surf_tag in surf_tags_for_this_panel:
                         surf_dim = surf_tag[0]
                         surf_id = surf_tag[1]
@@ -739,6 +754,24 @@ def Rz(theta):
                     )
                 if this_surf_bbox[2] > params.domain.z_max:
                     raise ValueError(f"A panel extends past the z_max wall.")
+        
+
+        # Loop over all the finalized volumes after fragmentation and tag everything
+        all_vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim)
+
+        for vol_tag in all_vol_tag_list:
+            vol_id = vol_tag[1]
+            com = self.gmsh_model.occ.getCenterOfMass(self.ndim, vol_id)
+
+            located_this_volume = False
+
+            for key, val in transformed_com.items():
+                    for target_com in val:
+                        # print(target_com)
+                        if np.allclose(np.array(com), target_com):
+                            located_this_volume = True
+                            if "trash" not in key:
+                                self._add_to_domain_markers(key, [vol_id], "cell")
 
         # Volumes are the entities with dimension equal to the mesh dimension
         vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim)
@@ -763,8 +796,10 @@ def Rz(theta):
 
         structure_panel_only_list = []
         structure_connector_only_list = []
-        
-        self._add_to_domain_markers("structure", structure_vol_list, "cell")
+# 
+        # # self._add_to_domain_markers("structure", structure_vol_list, "cell")
+        # for individual_vol_id in structure_vol_list:
+        #     self._add_to_domain_markers(f"structure_{individual_vol_id:03.0f}", [individual_vol_id], "cell")
         self._add_to_domain_markers("fluid", fluid_vol_list, "cell")
 
         # Record all the data collected to domain_markers as physics groups with physical names
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 862dbeb..76d61cc 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -464,14 +464,18 @@ def P_connector(u):
         
         self.calculate_K_for_Robin_BC(domain, flow, params)
 
+        dx_structure = ufl.Measure(
+            "dx", domain=domain.structure.msh, subdomain_data=domain.structure.cell_tags
+        )
+
         # To Do: how to differentiate the connector part and the panel part, dx_connector and dx_panel
         self.res = (
             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_panel
-            + k_nominal_connector(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * ufl.dx_connector
-            - structure.rho * ufl.inner(self.f, self.u_) * ufl.dx_panel
-            - structure.rho_connector * ufl.inner(self.f, self.u_) * ufl.dx_connector
+            + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * dx_structure(domain.domain_markers["modules"]["idx"])
+            + k_nominal_connector(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * dx_structure(domain.domain_markers["connectors"]["idx"])
+            - structure.rho * ufl.inner(self.f, self.u_) * dx_structure(domain.domain_markers["modules"]["idx"])
+            - structure.rho_connector * ufl.inner(self.f, self.u_) * dx_structure(domain.domain_markers["connectors"]["idx"])
             - ufl.dot(ufl.dot(self.stress_predicted * J * ufl.inv(F.T), n), self.u_)
             * self.ds
         )  # - Wext(self.u)
@@ -481,10 +485,8 @@ def P_connector(u):
             for i in range(params.pv_array.modules_per_span):
                 name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
                 K_springs = dolfinx.fem.Constant(domain.structure.msh, float(getattr(self, name_K)))
-                self.res += ufl.dot(K_springs * self.u_, self.z_unit_vector)*ds_bottom.panel_connector_markers[panel_id][i]
-        for i in range(params.pv_array.modules_per_span):
-            self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*ds_bottom....
-
+                self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*self.ds(domain.domain_markers[f"block_bottom_{panel_id:.0f}_{i:.0f}"])
+       
         # self.a = dolfinx.fem.form(ufl.lhs(res))
         # self.L = dolfinx.fem.form(ufl.rhs(res))
 

From dc61252eb6af24e8501b385d37f2af5f3ccbc423 Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Tue, 16 Dec 2025 14:14:46 -0700
Subject: [PATCH 20/26] new derivations without split the left fixed and right
 fixed part

---
 input/duramat_case_study.yaml             |   1 +
 pvade/IO/input_schema.yaml                |   4 +
 pvade/fluid/FlowManager.py                | 231 +++++++++-------
 pvade/geometry/MeshManager.py             |   2 +
 pvade/geometry/panels3d/DomainCreation.py |   6 +-
 pvade/structure/ElasticityAnalysis.py     | 321 ++++++++++++++++------
 pvade/structure/boundary_conditions.py    |   6 +-
 7 files changed, 377 insertions(+), 194 deletions(-)

diff --git a/input/duramat_case_study.yaml b/input/duramat_case_study.yaml
index 2221b4e..bbb8277 100644
--- a/input/duramat_case_study.yaml
+++ b/input/duramat_case_study.yaml
@@ -27,6 +27,7 @@ pv_array:
   torque_tube_outer_radius: 0.1 # radius of the torque tube
   torque_tube_inner_radius: 0.09 # radius of the torque tube
   modules_per_span: 10
+  fixed_location: 5   # fixed location of the panel along the span (from 0 (fixed at left) to modules_per_span (fixed at right))
   block_chord_div_by_panel_chord: 0.02
   
 solver:
diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index f9911da..f9f3bf2 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -221,6 +221,10 @@ properties:
         minimum: 1
         type: "integer"
         description: "The number of modules/panels per every panel_span distance. This can be thought of as the number of modules/panels per each table, counted in the spanwise direction only (i.e., 2-in-portait systems should report only the number in the spanwise direction, vs 2x that value)."
+      fixed_location:
+        default: 5
+        type: "integer"
+        description: "If an integer between 0 and modules_per_span, indicates the fixed location
       block_chord_div_by_panel_chord:
         default: 0.1
         minimum: 0.01
diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index f52ec07..bf892f4 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -583,33 +583,35 @@ def _all_interior_surfaces(x):
             "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
         )
     
-    def compute_panel_torques(self, domain, params):
+    # def compute_total_torques_on_each_panel(self, domain, params):
 
-        ds_fluid = ufl.Measure(
-            "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
-        )
+    #     ds_fluid = ufl.Measure(
+    #         "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
+    #     )
 
-        for panel_id in range(
-            int(params.pv_array.stream_rows * params.pv_array.span_rows)
-        ):
-            for module_id in range(params.pv_array.modules_per_span):
-                
-                total_torque = 0
+    #     for panel_id in range(
+    #         int(params.pv_array.stream_rows * params.pv_array.span_rows)
+    #     ):
+    #         total_torque = 0
 
-                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid(
-                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
-                )
-                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid(
-                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
-                )
-                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid(
-                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
-                )
-                total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid(
-                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
-                )
-                attr_name = f"total_torque_panel_{panel_id:.0f}_{module_id:.0f}"
-                setattr(self, attr_name, total_torque)
+    #         for module_id in range(params.pv_array.modules_per_span):
+                
+    #             total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid(
+    #                 domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
+    #             )
+    #             total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_X_coords, self.traction[2]) * ds_fluid(
+    #                 domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
+    #             )
+    #             total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid(
+    #                 domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
+    #             )
+    #             total_torque += dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(self.spatial_Z_coords, self.traction[0]) * ds_fluid(
+    #                 domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]))
+    #             )
+
+    #         attr_name = f"total_torque_panel_{panel_id:.0f}"
+
+    #         setattr(self, attr_name, total_torque)
 
     def compute_double_integral_panel_torques(self, domain, params):
         torque_function_on_panels = dolfinx.fem.Function(self.Q)
@@ -619,94 +621,125 @@ def compute_double_integral_panel_torques(self, domain, params):
         ds_fluid = ufl.Measure(
             "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
         )
+        
         for panel_id in range(
             int(params.pv_array.stream_rows * params.pv_array.span_rows)
         ):
+            whole_top_surface_facet_index = [] # facet index of the whole panel top surface in a row
+            whole_bot_surface_facet_index = [] # facet index of the whole panel bottom surface in a row
+
             for module_id in range(params.pv_array.modules_per_span):
                 
-                torque_double_integral = 0
+                whole_top_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist()
+                whole_bot_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist()
 
-                whole_top_surface_facet_index = domain.fluid.facet_tags.find(domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"])
-                whole_bot_surface_facet_index = domain.fluid.facet_tags.find(domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"])
-                whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(self.fluid.mesh, 2, whole_top_surface_facet_index)
-                whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(self.fluid.mesh, 2, whole_bot_surface_facet_index)
+            whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_top_surface_facet_index))
+            whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_bot_surface_facet_index))
 
-                dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh)
-                dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh)
+            dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh)
+            dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh)
 
-                whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1))
-                whole_top_submesh_func = dolfinx.fem.Function(whole_top_submesh_function_space)
-                whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1))
-                whole_bot_submesh_func = dolfinx.fem.Function(whole_bot_submesh_function_space)
+            whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1))
+            whole_top_submesh_func = dolfinx.fem.Function(whole_top_submesh_function_space)
+            whole_top_submesh_func.interpolate()(torque_function_on_panels)
+            whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1))
+            whole_bot_submesh_func = dolfinx.fem.Function(whole_bot_submesh_function_space)
+            whole_bot_submesh_func.interpolate()(torque_function_on_panels)
 
-                coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates()
-                coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates()
+            total_torque_on_panel_array = dolfinx.fem.assemble_scalar(dolfinx.fem.form(whole_top_submesh_func * dx_top_whole)) \
+                                        + dolfinx.fem.assemble_scalar(dolfinx.fem.form(whole_bot_submesh_func * dx_bot_whole))
+            
+            attr_name = f"total_torque_panel_{panel_id:.0f}"
 
-                top_integrated_func_along_x = np.zeros_like(whole_top_submesh_func.x.array[:])
+            setattr(self, attr_name, total_torque_on_panel_array)
 
-                bot_integrated_func_along_x = np.zeros_like(whole_bot_submesh_func.x.array[:])
-                
-                # Sort by x first, then y and z (group by x-levels)
-                sort_idx_top = np.lexsort((coords_top[:, 1], coords_top[:, 0]))  # x is sorted, from small to large
-                sort_idx_bot = np.lexsort((coords_bot[:, 1], coords_bot[:, 0]))  # x is sorted, from small to large
-                
-                top_coords_sorted = coords_top[sort_idx_top] #sorted from smallest x to largest x
-                bot_coords_sorted = coords_bot[sort_idx_bot] #sorted from smallest x to largest x
-                
+            coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates()
+            coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates()
 
-                ### integral of top surface:
-                # Unique x-values (up to numerical tolerance)
-                xs_top = np.unique(np.round(top_coords_sorted[:, 0], decimals=6))
+            top_integrated_func_along_y = np.zeros_like(whole_top_submesh_func.x.array[:])
+            bot_integrated_func_along_y = np.zeros_like(whole_bot_submesh_func.x.array[:])
+
+            # in PVade, the front has the minimum y coordinates, if we integral from large y to small y, it is negative,
+            # to keep it consistant as the theoretical model, we need to flip the sign of the integral
                 
-                for x_value in xs_top:
-                    # Mask for points at this x-level
-                    x_mask_top = np.isclose(top_coords_sorted[:, 0], x_value, atol=1e-6)
-
-                    facet_centers_top = dolfinx.mesh.compute_midpoints(whole_top_surf_submesh, 2, np.array(np.arange(whole_top_surf_submesh.topology.index_map(2).size_local), dtype=np.int32))
-                    # Identify cells where the midpoint x-coordinate is less than 1
-                    marked_facet_top = np.where(facet_centers_top[:, 0] <= x_value)[0]
-                    
-                    # Create a submesh from those cells
-                    submesh_top_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_top_surf_submesh, 2, marked_facet_top.astype(np.int32))
-                    top_sub_function_space = dolfinx.fem.FunctionSpace(submesh_top_surface_part, ('CG', 1))
-                    top_sub_func = dolfinx.fem.Function(top_sub_function_space)
-                    top_sub_func.interpolate(torque_function_on_panels)
-
-                    dx_top = ufl.Measure("dx", domain=submesh_top_surface_part)
-                    integrated_top = dolfinx.fem.assemble_scalar(dolfinx.fem.form(top_sub_func * dx_top))
-                    top_integrated_func_along_x[sort_idx_top[x_mask_top]] = integrated_top
-
-                ### integral of bot surface:
-                # Unique x-values (up to numerical tolerance)
-                xs_bot = np.unique(np.round(bot_coords_sorted[:, 0], decimals=6))
-
-                for x_value in xs_bot:
-
-                    # Mask for points at this x-level
-                    x_mask_bot = np.isclose(bot_coords_sorted[:, 0], x_value, atol=1e-6)
-                    facet_centers_bot = dolfinx.mesh.compute_midpoints(whole_bot_surf_submesh, 2, np.array(np.arange(whole_bot_surf_submesh.topology.index_map(2).size_local), dtype=np.int32))
-                    # Identify cells where the midpoint x-coordinate is less than 1
-                    marked_facet_bot = np.where(facet_centers_bot[:, 0] <= x_value)[0]
-                    # Create a submesh from those cells
-                    submesh_bot_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_bot_surf_submesh, 2, marked_facet_bot.astype(np.int32))
-                    bot_sub_function_space = dolfinx.fem.FunctionSpace(submesh_bot_surface_part, ('CG', 1))
-                    bot_sub_func = dolfinx.fem.Function(bot_sub_function_space)
-                    bot_sub_func.interpolate(torque_function_on_panels)
-
-                    dx_bot = ufl.Measure("dx", domain=submesh_bot_surface_part)
-                    integrated_bot = dolfinx.fem.assemble_scalar(dolfinx.fem.form(bot_sub_func * dx_bot))
-                    bot_integrated_func_along_x[sort_idx_bot[x_mask_bot]] = integrated_bot
-
-                single_integral_function_top = dolfinx.fem.Function(whole_top_submesh_function_space, name='single_integral_function_top')  
-                single_integral_function_top.x.array[:] = top_integrated_func_along_x   
-
-                single_integral_function_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name='single_integral_function_bot')  
-                single_integral_function_bot.x.array[:] = bot_integrated_func_along_x 
-
-                torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_top * dx_top_whole))/self.pvarry.panel_span*params.pv_array.modules_per_span
-                torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_bot * dx_bot_whole))/self.pvarry.panel_span*params.pv_array.modules_per_span
-                attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{module_id:.0f}"
+            # Sort by y first, then z and z (group by y-levels)
+            sort_idx_top = np.lexsort((coords_top[:, 0], -coords_top[:, 1]))  # y is sorted, from large to small 
+            sort_idx_bot = np.lexsort((coords_bot[:, 0], -coords_bot[:, 1]))  # y is sorted, from large to small 
+
+            top_coords_sorted = coords_top[sort_idx_top] #sorted from largest y to smallest y
+            bot_coords_sorted = coords_bot[sort_idx_bot] #sorted from largest y to smallest y
+
+            ### integral of top surface:
+            # Unique x-values (up to numerical tolerance)
+            ys_top = np.unique(np.round(top_coords_sorted[:, 0], decimals=6))
+
+            # get f(x) = integral of torque on top surface from left end to x (x is all unique x coordinates on top surface)
+
+            for y_value in ys_top:
+                # Mask for points at this y-level
+                y_mask_top = np.isclose(top_coords_sorted[:, 1], y_value, atol=1e-6)
+
+                facet_centers_top = dolfinx.mesh.compute_midpoints(whole_top_surf_submesh, 2, np.array(np.arange(whole_top_surf_submesh.topology.index_map(2).size_local), dtype=np.int32))
+                # Identify cells where the midpoint y-coordinate is higher
+                marked_facet_top = np.where(facet_centers_top[:, 1] >= y_value)[0]
+
+                # Create a submesh from those cells
+                submesh_top_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_top_surf_submesh, 2, marked_facet_top.astype(np.int32))
+                top_sub_function_space = dolfinx.fem.FunctionSpace(submesh_top_surface_part, ('CG', 1))
+                top_sub_func = dolfinx.fem.Function(top_sub_function_space)
+                top_sub_func.interpolate(torque_function_on_panels)
+
+                dx_top = ufl.Measure("dx", domain=submesh_top_surface_part)
+                integrated_top = dolfinx.fem.assemble_scalar(dolfinx.fem.form(top_sub_func * dx_top))
+                top_integrated_func_along_y[sort_idx_top[y_mask_top]] = integrated_top
+
+            ### integral of bot surface:
+            # Unique y-values (up to numerical tolerance)
+            ys_bot = np.unique(np.round(bot_coords_sorted[:, 1], decimals=6))
+
+            for y_value in ys_bot:
+
+                # Mask for points at this y-level
+                y_mask_bot = np.isclose(bot_coords_sorted[:, 1], y_value, atol=1e-6)
+                facet_centers_bot = dolfinx.mesh.compute_midpoints(whole_bot_surf_submesh, 2, np.array(np.arange(whole_bot_surf_submesh.topology.index_map(2).size_local), dtype=np.int32))
+                # Identify cells where the midpoint y-coordinate is less than 1
+                marked_facet_bot = np.where(facet_centers_bot[:, 1] >= y_value)[0]
+                # Create a submesh from those cells
+                submesh_bot_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_bot_surf_submesh, 2, marked_facet_bot.astype(np.int32))
+                bot_sub_function_space = dolfinx.fem.FunctionSpace(submesh_bot_surface_part, ('CG', 1))
+                bot_sub_func = dolfinx.fem.Function(bot_sub_function_space)
+                bot_sub_func.interpolate(torque_function_on_panels)
+
+                dx_bot = ufl.Measure("dx", domain=submesh_bot_surface_part)
+                integrated_bot = dolfinx.fem.assemble_scalar(dolfinx.fem.form(bot_sub_func * dx_bot))
+                bot_integrated_func_along_y[sort_idx_bot[y_mask_bot]] = integrated_bot
+
+            single_integral_function_top = dolfinx.fem.Function(whole_top_submesh_function_space, name='single_integral_function_top')  
+            single_integral_function_top.x.array[:] = top_integrated_func_along_y
+
+            single_integral_function_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name='single_integral_function_bot')  
+            single_integral_function_bot.x.array[:] = bot_integrated_func_along_y
+
+            torque_double_integral = 0.0
+
+            attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_0"
+            setattr(self, attr_name, torque_double_integral)
+
+            single_integral_function_top_fluid = dolfinx.fem.Function(self.Q)
+            single_integral_function_top_fluid.interpolate(single_integral_function_top)
+            single_integral_function_bot_fluid = dolfinx.fem.Function(self.Q)
+            single_integral_function_bot_fluid.interpolate(single_integral_function_bot)
+
+            for connector_id in range(params.pv_array.modules_per_span):
+                torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_top_fluid * ds_fluid(
+                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{params.panel_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
+                torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_bot_fluid * ds_fluid(
+                    domain.domain_markers[f"panel_bot_{panel_id:.0f}_{params.panel_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
+                attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{connector_id+1:.0f}"
                 setattr(self, attr_name, torque_double_integral)
+
+                # double_integral_total_torque_panel_0_0 is the double integral of first panel array at the most back connector (maximum y), 
+                # double_integral_total_torque_panel_0_10 is the double integral of first panel array at the most front connector (smallest y)
            
     def _assemble_system(self, params):
         """Pre-assemble all LHS matrices and RHS vectors
@@ -856,7 +889,7 @@ def solve(self, domain, params, current_time):
 
         self.compute_lift_and_drag(params, current_time)
 
-        self.compute_panel_torques(domain, params)
+        # self.compute_panel_torques(domain, params)
         self.compute_double_integral_panel_torques(domain, params)
 
         # Compute the pressure drop between the inlet and outlet
diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py
index 618f064..11f54be 100644
--- a/pvade/geometry/MeshManager.py
+++ b/pvade/geometry/MeshManager.py
@@ -283,6 +283,8 @@ def _create_submeshes_from_parent(self, params):
                 print(f"Creating {sub_domain_name} submesh")
 
             # Get the idx associated with either "fluid" or "structure"
+
+            # if structure includes modules and connectors
             if sub_domain_name == "structure" and "structure" not in self.domain_markers:
                 marker_id = self.domain_markers["modules"]["idx"] 
                 # Find all cells where cell tag = marker_id
diff --git a/pvade/geometry/panels3d/DomainCreation.py b/pvade/geometry/panels3d/DomainCreation.py
index 4709c9b..fdd46af 100644
--- a/pvade/geometry/panels3d/DomainCreation.py
+++ b/pvade/geometry/panels3d/DomainCreation.py
@@ -498,9 +498,10 @@ def Rz(theta):
                         if (
                             np.isclose(com[2], params.pv_array.torque_tube_separation) and modeling_torque_tube and np.isclose(com[1], module_distances[params.pv_array.modules_per_span] - params.pv_array.block_chord_div_by_panel_chord * half_chord/2.0)
                         ):
-                            target_key = f"interior_surface_{panel_ct:.0f}"
+                            target_key = f"interior_surface_{panel_ct:.0f}" # block/panel interface
                             surface_located_or_not = True
 
+                        # if not the most left/right panel boundary, it is panel/panel interface
                         if not surface_located_or_not:
                             for module_id in range(params.pv_array.modules_per_span):
 
@@ -756,7 +757,7 @@ def Rz(theta):
                     raise ValueError(f"A panel extends past the z_max wall.")
         
 
-        # Loop over all the finalized volumes after fragmentation and tag everything
+        # mark the panel and connector volumes
         all_vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim)
 
         for vol_tag in all_vol_tag_list:
@@ -773,6 +774,7 @@ def Rz(theta):
                             if "trash" not in key:
                                 self._add_to_domain_markers(key, [vol_id], "cell")
 
+        # Mark the fluid volume
         # Volumes are the entities with dimension equal to the mesh dimension
         vol_tag_list = self.gmsh_model.occ.getEntities(self.ndim)
         structure_vol_list = []
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 76d61cc..24d344c 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -101,128 +101,269 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
         # total number of pv rows
         total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows
 
-        num_panel_right_fixed = params.pv_array.modules_per_span // 2
-        num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed
+        for panel_id in range(total_num_panels):
+            # construct the matrix, size: modules_per_span+1 by modules_per_span+1
+            theo_matrix = np.zeros((params.pv_array.modules_per_span+1, params.pv_array.modules_per_span+1))
+            theo_vector = np.zeros(params.pv_array.modules_per_span+1)
+
+            connector_locations = np.linspace(
+                0, params.pv_array.panel_span, params.pv_array.modules_per_span + 1
+            )
 
-        # for right fixed part:
+            total_torque_on_this_panel_name = f"total_torque_panel_{panel_id:.0f}"
+            total_torque_on_this_panel = getattr(flow, total_torque_on_this_panel_name)
 
-        for panel_id in range(total_num_panels):
-            total_torque_right_fixed = 0
-            total_torque_left_fixed = 0
+            theo_matrix[params.pv_array.modules_per_span, :] = 1/Gp/Ipp
+            theo_vector[params.pv_array.modules_per_span] = -total_torque_on_this_panel/Gp/Ipp
+
+            # contribution from panel:
+
+                # contribution from C0 to matrix:
+
+            theo_matrix[:params.pv_array.modules_per_span, :params.pv_array.fixed_location] += connector_locations[params.pv_array.fixed_location]/(Gp*Ipp)
+            theo_matrix[:params.pv_array.modules_per_span, 1:params.pv_array.fixed_location] -= connector_locations[1:params.pv_array.fixed_location]/(Gp*Ipp)
+
+                # contribution from C0 to vector:
+            T_double_integral_at_fixed_location_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.fixed_location:.0f}"
+            T_double_integral_at_fixed_location = getattr(flow, T_double_integral_at_fixed_location_name)
+            theo_vector[:params.pv_array.modules_per_span] += -T_double_integral_at_fixed_location/Gp/Ipp
+            
+                # contribution of panels to matrix
+            theo_matrix[0, 0] += -connector_locations[0]/(Gp*Ipp)
+            theo_matrix[1, 0] += -connector_locations[1]/(Gp*Ipp)
+
+            for i in range(2, params.pv_array.fixed_location):
+                theo_matrix[i, :i] += -connector_locations[i]/(Gp*Ipp)
+                theo_matrix[i, 1:i] += connector_locations[1:i]/(Gp*Ipp)
+
+            for i in range(params.pv_array.fixed_location, params.pv_array.modules_per_span):
+                theo_matrix[i, :i+1] += -connector_locations[i+1]/(Gp*Ipp)
+                theo_matrix[i, 1:i+1] += connector_locations[1:i+1]/(Gp*Ipp)
+
+                # contribution of panels to vector
+            T_double_integral_array = []
+            for i in range(params.pv_array.modules_per_span+1):
+                name = f"double_integral_total_torque_panel_{panel_id:.0f}_{i:.0f}"
+                T_double_integral_array.append(getattr(flow, name))
+            T_double_integral_array = np.array(T_double_integral_array)
+            theo_vector[:params.pv_array.modules_per_span] += np.delete(T_double_integral_array, params.pv_array.fixed_location)/Gp/Ipp
+
+            # contribution from tube:
+
+                # contribution from C0 to matrix:
+
+            theo_matrix[:params.pv_array.modules_per_span, params.pv_array.fixed_location:params.pv_array.modules_per_span+1] += -connector_locations[params.pv_array.fixed_location]/(Gt*Ipt)
+            theo_matrix[:params.pv_array.modules_per_span, 1:params.pv_array.fixed_location] += -connector_locations[1:params.pv_array.fixed_location]/(Gt*Ipt)
+
+                # contribution from C0 to vector:
+            theo_vector[:params.pv_array.modules_per_span] += total_torque_on_this_panel*connector_locations[params.pv_array.fixed_location]/Gt/Ipt
+
+                # contribution of tube to matrix
+            theo_matrix[0, 1:params.pv_array.modules_per_span+1] += connector_locations[0]/(Gt*Ipt)
+            theo_matrix[1, 1:params.pv_array.modules_per_span+1] += connector_locations[1]/(Gt*Ipt)
+
+            for i in range(2, params.pv_array.fixed_location):
+                theo_matrix[i, i:params.pv_array.modules_per_span+1] += connector_locations[i]/(Gt*Ipt)
+                theo_matrix[i, 1:i] += connector_locations[1:i]/(Gt*Ipt)
+
+            for i in range(params.pv_array.fixed_location, params.pv_array.modules_per_span):
+                theo_matrix[i, i+1:params.pv_array.modules_per_span+1] += connector_locations[i+1]/(Gt*Ipt)
+                theo_matrix[i, 1:i+1] += connector_locations[1:i+1]/(Gt*Ipt)
+
+                # contribution of tube to vector
+            theo_vector[:params.pv_array.fixed_location] += -total_torque_on_this_panel/Gt/Ipt*connector_locations[:params.pv_array.fixed_location]
+            theo_vector[params.pv_array.fixed_location:params.pv_array.modules_per_span] += -total_torque_on_this_panel/Gt/Ipt*connector_locations[params.pv_array.fixed_location]
+
+            R_reaction_torque = np.dot(np.linalg.inv(theo_matrix), theo_vector)  # this is the torque applied to tube
+            
+            # rotation of tube at each connector location
+            tube_rotate_matrix = np.zeros((params.pv_array.modules_per_span, params.pv_array.modules_per_span+1))
+            tube_rotate_vector = np.zeros(params.pv_array.modules_per_span)
+
+            # contribution from panel:
+
+                # contribution from C0 to matrix:
+
+            tube_rotate_matrix[:params.pv_array.modules_per_span, :params.pv_array.fixed_location] += connector_locations[params.pv_array.fixed_location]/(Gp*Ipp)
+            tube_rotate_matrix[:params.pv_array.modules_per_span, 1:params.pv_array.fixed_location] -= connector_locations[1:params.pv_array.fixed_location]/(Gp*Ipp)
+
+                # contribution from C0 to vector:
+            tube_rotate_vector[:params.pv_array.modules_per_span] += T_double_integral_array[params.pv_array.fixed_location]/Gp/Ipp
+
+
+                # contribution of panels to matrix
+            tube_rotate_matrix[0, 0] += -connector_locations[0]/(Gp*Ipp)
+            tube_rotate_matrix[1, 0] += -connector_locations[1]/(Gp*Ipp)
+
+            for i in range(2, params.pv_array.fixed_location):
+                tube_rotate_matrix[i, :i] += -connector_locations[i]/(Gp*Ipp)
+                tube_rotate_matrix[i, 1:i] += connector_locations[1:i]/(Gp*Ipp)
+
+            for i in range(params.pv_array.fixed_location, params.pv_array.modules_per_span):
+                tube_rotate_matrix[i, :i+1] += -connector_locations[i+1]/(Gp*Ipp)
+                tube_rotate_matrix[i, 1:i+1] += connector_locations[1:i+1]/(Gp*Ipp)
+
+                # contribution of panels to vector
+            tube_rotate_vector[:params.pv_array.modules_per_span] += -np.delete(T_double_integral_array, params.pv_array.fixed_location)/Gp/Ipp
+
+
+            phi = np.dot(tube_rotate_matrix, R_reaction_torque) + tube_rotate_vector
+
+            if isinstance(params.pv_array.tracker_angle, list):
+                    if panel_id == 0:
+                        assert (
+                            len(params.pv_array.tracker_angle)
+                            == params.pv_array.stream_rows * params.pv_array.span_rows
+                        ), f"Length of tracker angle list ({len(params.pv_array.tracker_angle)}) not equal to total number of PV tables ({params.pv_array.stream_rows * params.pv_array.span_rows})."
+
+                    tracker_angle_rad = np.radians(
+                        params.pv_array.tracker_angle[panel_id]
+                    )
+            else:
+                tracker_angle_rad = np.radians(params.pv_array.tracker_angle)
+
+
+            K = np.abs(12*(np.delete(R_reaction_torque, params.pv_array.fixed_location))/((params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord)**3)/np.cos(tracker_angle_rad+phi)/(np.sin(tracker_angle_rad+phi)-np.sin(tracker_angle_rad))/(params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord/2))
+
+            # assume K is 0 at the fixed connector
+            K = np.flip(np.insert(K, params.pv_array.fixed_location, 0))
+
+            # K[0] is the stiffness at the most back connector (highest y)
+            # K[10] is the stiffness at the most front connector (lowest y)
+            # for block_bot_surface, it is numbered from lowest y to highest y, so flip K
             
-            for i in range(num_panel_left_fixed, params.pv_array.modules_per_span):
-                name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}"
-                total_torque_right_fixed += getattr(flow, name)
+            # if there are 10 modules per array, there are 11 connectors, K has shape of 10, the K at the fixed connector is not included.
+            for i in range(params.pv_array.modules_per_span+1):
+
+                name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
+                setattr(self, name_K, K[i])
+               
+
+                
+
+
+
+        # num_panel_right_fixed = params.pv_array.modules_per_span // 2
+        # num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed
+
+        # # for right fixed part:
+
+        # for panel_id in range(total_num_panels):
+        #     total_torque_right_fixed = 0
+        #     total_torque_left_fixed = 0
+            
+        #     for i in range(num_panel_left_fixed, params.pv_array.modules_per_span):
+        #         name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}"
+        #         total_torque_right_fixed += getattr(flow, name)
             
-            for i in range(num_panel_left_fixed):
-                name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}"
-                total_torque_left_fixed += getattr(flow, name)
+        #     for i in range(num_panel_left_fixed):
+        #         name = f"total_torque_panel_{panel_id:.0f}_{i:.0f}"
+        #         total_torque_left_fixed += getattr(flow, name)
 
-            theo_matrix_right_fixed = np.zeros((num_panel_right_fixed+1, num_panel_right_fixed+1))
-            theo_vector_right_fixed = np.zeros(num_panel_right_fixed+1)
+        #     theo_matrix_right_fixed = np.zeros((num_panel_right_fixed+1, num_panel_right_fixed+1))
+        #     theo_vector_right_fixed = np.zeros(num_panel_right_fixed+1)
 
-            theo_matrix_right_fixed[0, :] = 1.0
-            theo_vector_right_fixed[0] = total_torque_right_fixed
-            theo_matrix_right_fixed[1:, :-1] = params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
+        #     theo_matrix_right_fixed[0, :] = 1.0
+        #     theo_vector_right_fixed[0] = total_torque_right_fixed
+        #     theo_matrix_right_fixed[1:, :-1] = params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
             
-            # as the double integral of torque along the span from front to back, left fixed part to right fixed part, it need to be flipped
-            T_double_integral_right_fixed = []
-            T_right_fixed = []
-            for i in range(params.pv_array.modules_per_span):
-                name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}"
-                T_double_integral_right_fixed.append(getattr(flow, name))
-                name = f"total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}"
-                T_right_fixed.append(getattr(flow, name))
-            for i in range(num_panel_right_fixed):
-                for j in range(i+1):
-                    theo_vector_right_fixed[i+1] += T_double_integral_right_fixed[-1-j]/(Gp*Ipp)
-                    theo_vector_right_fixed[i+1] -= (i+1-j)*T_right_fixed[-1-j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span
-                for j in range(i):
-                    theo_matrix_right_fixed[i+1, :-1-j-1] += params.pv_array.panel_span/params.pv_array.modules_per_span/(Gt*Ipt)
-
-
-            for i in range(num_panel_right_fixed):
-                for j in range(i+1):
-                        theo_matrix_right_fixed[i+1, num_panel_right_fixed-j:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gp/Ipp
-
-            R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_right_fixed), theo_vector_right_fixed)
+        #     # as the double integral of torque along the span from front to back, left fixed part to right fixed part, it need to be flipped
+        #     T_double_integral_right_fixed = []
+        #     T_right_fixed = []
+        #     for i in range(params.pv_array.modules_per_span):
+        #         name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}"
+        #         T_double_integral_right_fixed.append(getattr(flow, name))
+        #         name = f"total_torque_panel_{panel_id:.0f}_{params.pv_array.modules_per_span-1-i:.0f}"
+        #         T_right_fixed.append(getattr(flow, name))
+        #     for i in range(num_panel_right_fixed):
+        #         for j in range(i+1):
+        #             theo_vector_right_fixed[i+1] += T_double_integral_right_fixed[-1-j]/(Gp*Ipp)
+        #             theo_vector_right_fixed[i+1] -= (i+1-j)*T_right_fixed[-1-j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span
+        #         for j in range(i):
+        #             theo_matrix_right_fixed[i+1, :-1-j-1] += params.pv_array.panel_span/params.pv_array.modules_per_span/(Gt*Ipt)
+
+
+        #     for i in range(num_panel_right_fixed):
+        #         for j in range(i+1):
+        #                 theo_matrix_right_fixed[i+1, num_panel_right_fixed-j:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gp/Ipp
+
+        #     R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_right_fixed), theo_vector_right_fixed)
             
-            tube_rotate_matrix = np.ones((num_panel_right_fixed, num_panel_right_fixed))*params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt)
+        #     tube_rotate_matrix = np.ones((num_panel_right_fixed, num_panel_right_fixed))*params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt)
 
-            for i in range(num_panel_right_fixed):
-                for j in range(i):
-                    tube_rotate_matrix[i, :num_panel_right_fixed-1-j] += params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt)
+        #     for i in range(num_panel_right_fixed):
+        #         for j in range(i):
+        #             tube_rotate_matrix[i, :num_panel_right_fixed-1-j] += params.pv_array.panel_span/params.pv_array.modules_per_span*(1.0/Gt/Ipt)
             
-            # check the standalone code, why it need to be flipped.
-            phi = np.flip(np.dot(tube_rotate_matrix, R_reaction_torque[:-1]))
+        #     # check the standalone code, why it need to be flipped.
+        #     phi = np.flip(np.dot(tube_rotate_matrix, R_reaction_torque[:-1]))
 
-            # rotation of connector
-            x_panel = np.arange(num_panel_right_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span) # location of panels
+        #     # rotation of connector
+        #     x_panel = np.arange(num_panel_right_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span) # location of panels
 
-            block_length = params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord
-            block_width = params.pv_array.block_chord_div_by_panel_chord*params.pv_array.panel_span/params.pv_array.modules_per_span/2
+        #     block_length = params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord
+        #     block_width = params.pv_array.block_chord_div_by_panel_chord*params.pv_array.panel_span/params.pv_array.modules_per_span/2
 
-            # To Do: check the angle, is this correct?
-            array_rotation = (params.fluid.wind_direction + 90.0) % 360.0
-            array_rotation_rad = np.radians(array_rotation)
+        #     # To Do: check the angle, is this correct?
+        #     array_rotation = (params.fluid.wind_direction + 90.0) % 360.0
+        #     array_rotation_rad = np.radians(array_rotation)
 
-            for i in range(num_panel_right_fixed):
+        #     for i in range(num_panel_right_fixed):
 
-                phi_i = phi[i]
+        #         phi_i = phi[i]
 
-                K_i = 12*(R_reaction_torque[i])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width
+        #         K_i = 12*(R_reaction_torque[i])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width
 
-                name_K = f"spring_stiffness_{panel_id:.0f}_{params.pv_array.modules_per_span+1-i:.0f}"
-                
-                setattr(self, name_K, K_i)
+        #         name_K = f"spring_stiffness_{panel_id:.0f}_{params.pv_array.modules_per_span+1-i:.0f}"
+        #         setattr(self, name_K, K_i)
 
-            # for left fixed part:
-            theo_matrix_left_fixed = np.zeros((num_panel_left_fixed+1, num_panel_left_fixed+1))
-            theo_vector_left_fixed = np.zeros(num_panel_left_fixed+1)
-            theo_matrix_left_fixed[0, :] = 1.0
-            theo_vector_left_fixed[0] = total_torque_left_fixed
+        #     # for left fixed part:
+        #     theo_matrix_left_fixed = np.zeros((num_panel_left_fixed+1, num_panel_left_fixed+1))
+        #     theo_vector_left_fixed = np.zeros(num_panel_left_fixed+1)
+        #     theo_matrix_left_fixed[0, :] = 1.0
+        #     theo_vector_left_fixed[0] = total_torque_left_fixed
 
-            T_double_integral_left_fixed = []
-            T_left_fixed = []
+        #     T_double_integral_left_fixed = []
+        #     T_left_fixed = []
 
-            for i in range(num_panel_left_fixed):
-                name = f"double_integral_total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
-                T_double_integral_left_fixed.append(getattr(flow, name))
-                name = f"total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
-                T_left_fixed.append(getattr(flow, name))
+        #     for i in range(num_panel_left_fixed):
+        #         name = f"double_integral_total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
+        #         T_double_integral_left_fixed.append(getattr(flow, name))
+        #         name = f"total_torque_panel_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
+        #         T_left_fixed.append(getattr(flow, name))
 
-            for i in range(num_panel_left_fixed):
-                for j in range(i+1):
-                    theo_matrix_left_fixed[i+1, j] = (i+1-j)*params.pv_array.panel_span/params.pv_array.modules_per_span/(Gp*Ipp)
-                    theo_vector_left_fixed[i+1] += T_double_integral_left_fixed[j]/(Gp*Ipp)
-                for j in range(i):
-                    theo_vector_left_fixed[i+1] += (i-j)*T_left_fixed[j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span
+        #     for i in range(num_panel_left_fixed):
+        #         for j in range(i+1):
+        #             theo_matrix_left_fixed[i+1, j] = (i+1-j)*params.pv_array.panel_span/params.pv_array.modules_per_span/(Gp*Ipp)
+        #             theo_vector_left_fixed[i+1] += T_double_integral_left_fixed[j]/(Gp*Ipp)
+        #         for j in range(i):
+        #             theo_vector_left_fixed[i+1] += (i-j)*T_left_fixed[j]/(Gp*Ipp)*params.pv_array.panel_span/params.pv_array.modules_per_span
 
-            for i in range(num_panel_left_fixed):    
-                theo_matrix_left_fixed[i+1:, i+1:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
+        #     for i in range(num_panel_left_fixed):    
+        #         theo_matrix_left_fixed[i+1:, i+1:] -= params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
 
-            R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_left_fixed), theo_vector_left_fixed)  # this is the torque applied to tube
+        #     R_reaction_torque = np.dot(np.linalg.inv(theo_matrix_left_fixed), theo_vector_left_fixed)  # this is the torque applied to tube
 
 
-            tube_rotate_matrix = np.zeros((num_panel_left_fixed, num_panel_left_fixed))
-            for i in range(num_panel_left_fixed):    
-                tube_rotate_matrix[i:, i:] += params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
+        #     tube_rotate_matrix = np.zeros((num_panel_left_fixed, num_panel_left_fixed))
+        #     for i in range(num_panel_left_fixed):    
+        #         tube_rotate_matrix[i:, i:] += params.pv_array.panel_span/params.pv_array.modules_per_span/Gt/Ipt
 
-            phi = np.dot(tube_rotate_matrix, R_reaction_torque[1:])
+        #     phi = np.dot(tube_rotate_matrix, R_reaction_torque[1:])
 
-            # rotation of connector
-            x_panel = np.arange(num_panel_left_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span)+params.pv_array.panel_span/params.pv_array.modules_per_span # location of panles
+        #     # rotation of connector
+        #     x_panel = np.arange(num_panel_left_fixed)*(params.pv_array.panel_span/params.pv_array.modules_per_span)+params.pv_array.panel_span/params.pv_array.modules_per_span # location of panles
 
 
-            for i in range(num_panel_left_fixed):
+        #     for i in range(num_panel_left_fixed):
 
-                phi_i = phi[i]
+        #         phi_i = phi[i]
 
-                K_i = 12*(R_reaction_torque[i+1])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width
+        #         K_i = 12*(R_reaction_torque[i+1])/((block_length)**3)/np.cos(array_rotation_rad+phi_i)/(np.sin(array_rotation_rad+phi_i)-np.sin(array_rotation_rad))/block_width
 
-                name_K = f"spring_stiffness_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
+        #         name_K = f"spring_stiffness_{panel_id:.0f}_{num_panel_left_fixed-1-i:.0f}"
                 
-                setattr(self, name_K, K_i)
+        #         setattr(self, name_K, K_i)
 
     def update_a(self, u, u_old, v_old, a_old, dt, beta, ufl=True):
         # Update formula for acceleration
@@ -482,7 +623,7 @@ def P_connector(u):
 
         # Robin boundary condition terms
         for panel_id in range(params.pv_array.stream_rows * params.pv_array.span_rows):
-            for i in range(params.pv_array.modules_per_span):
+            for i in range(params.pv_array.modules_per_span+1):
                 name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
                 K_springs = dolfinx.fem.Constant(domain.structure.msh, float(getattr(self, name_K)))
                 self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*self.ds(domain.domain_markers[f"block_bottom_{panel_id:.0f}_{i:.0f}"])
diff --git a/pvade/structure/boundary_conditions.py b/pvade/structure/boundary_conditions.py
index ebfe5f7..956b857 100644
--- a/pvade/structure/boundary_conditions.py
+++ b/pvade/structure/boundary_conditions.py
@@ -272,7 +272,7 @@ def build_structure_boundary_conditions(domain, params, functionspace):
     total_num_panels = params.pv_array.stream_rows * params.pv_array.span_rows
 
     for num_panel in range(total_num_panels):
-        for location in params.structure.bc_list:
+        for location in params.structure.bc_list: # it is empty
             location_panel = f"{location}_{num_panel}"
             # f"front_{num_panel}" , f"back_{num_panel}":
             # for location in  [f"left_{num_panel}"]:# , f"right_{num_panel}":
@@ -362,11 +362,11 @@ def connection_point_up_helper(nodes_to_pin_between):
 
     if params.structure.motor_connection == True:
         # The bottom surface of the center connector is fixed to represent the motor mount
-        motor_location = params.pv_array.modules_per_span // 2
+        motor_location = params.pv_array.fixed_location  # fixed location along the span, if fixed_location=5, it means the left boundary of the 6th connector is fixed
         for panel_id in range(total_num_panels):
 
             mount_facet = domain.structure.facet_tags.find(
-                domain.domain_markers[f"block_bottom_{panel_id:.0f}_{motor_location:.0f}"]["idx"]
+                domain.domain_markers[f"block_left_{panel_id:.0f}_{motor_location:.0f}"]["idx"]
         )
         
             dofs_disp = dolfinx.fem.locate_dofs_topological(

From 3f764ece8771b49d2482f0de98e64e0ac0179edd Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Tue, 16 Dec 2025 14:55:02 -0700
Subject: [PATCH 21/26] fix bugs in the input file

---
 input/duramat_case_study.yaml         | 2 +-
 pvade/IO/input_schema.yaml            | 4 ++--
 pvade/fluid/FlowManager.py            | 2 ++
 pvade/geometry/MeshManager.py         | 2 +-
 pvade/structure/ElasticityAnalysis.py | 4 ++++
 pvade_main.py                         | 2 +-
 6 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/input/duramat_case_study.yaml b/input/duramat_case_study.yaml
index bbb8277..1f40a19 100644
--- a/input/duramat_case_study.yaml
+++ b/input/duramat_case_study.yaml
@@ -68,4 +68,4 @@ structure:
   beta_relaxation: 0.5
   elasticity_modulus_tube: 2.0e+11
   poissons_ratio_tube: 0.3
-  density_tube: 7800.0
\ No newline at end of file
+  rho_tube: 7800.0
\ No newline at end of file
diff --git a/pvade/IO/input_schema.yaml b/pvade/IO/input_schema.yaml
index f9f3bf2..18f1beb 100644
--- a/pvade/IO/input_schema.yaml
+++ b/pvade/IO/input_schema.yaml
@@ -224,7 +224,7 @@ properties:
       fixed_location:
         default: 5
         type: "integer"
-        description: "If an integer between 0 and modules_per_span, indicates the fixed location
+        description: "If an integer between 0 and modules_per_span, indicates the fixed location"
       block_chord_div_by_panel_chord:
         default: 0.1
         minimum: 0.01
@@ -695,7 +695,7 @@ properties:
         type: "number"
         description: "poissons ratio of the torque tube."
         units: "None"
-      density_tube:
+      rho_tube:
         default: 1.0
         # minimum: 0.001
         # maximum: 1000.0
diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index bf892f4..197a2d4 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -653,6 +653,8 @@ def compute_double_integral_panel_torques(self, domain, params):
 
             setattr(self, attr_name, total_torque_on_panel_array)
 
+            print('check total torque definiition', self.total_torque_panel_0)
+
             coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates()
             coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates()
 
diff --git a/pvade/geometry/MeshManager.py b/pvade/geometry/MeshManager.py
index 11f54be..f35a7aa 100644
--- a/pvade/geometry/MeshManager.py
+++ b/pvade/geometry/MeshManager.py
@@ -747,7 +747,7 @@ def test_mesh_functionspace(self):
 
         print(f"Rank {self.rank} owns {num_nodes_owned_by_proc} nodes\n{coords}")
 
-    def test_submesh_transfer(self, params):
+    def test_submesh_transfer(self, params, domain, elasticity):
         P2 = ufl.VectorElement("Lagrange", self.msh.ufl_cell(), 2)
         # P2 = ufl.FiniteElement("Lagrange", self.fluid.msh.ufl_cell(), 1)
 
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 24d344c..38fecea 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -110,6 +110,10 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
                 0, params.pv_array.panel_span, params.pv_array.modules_per_span + 1
             )
 
+            # ?? double integral is not available until the second time step. Should make it be zero for the 
+            # first time step then it will be updated once the first step is finished? but the assemble is out of 
+            # the time loop. Should we solve fluid then structure?
+            
             total_torque_on_this_panel_name = f"total_torque_panel_{panel_id:.0f}"
             total_torque_on_this_panel = getattr(flow, total_torque_on_this_panel_name)
 
diff --git a/pvade_main.py b/pvade_main.py
index 6f262f5..2af05fe 100644
--- a/pvade_main.py
+++ b/pvade_main.py
@@ -17,7 +17,7 @@
 from mpi4py import MPI
 
 """
-to run: python -u pvade_main.py --input_file=input/duramat_case_study.yaml --domain.l_char=3.0 --general.mesh_only=true --pv_array.torque_tube_radius=0.1 --pv_array.torque_tube_separation=0.4 --pv_array.stream_rows=2 --pv_array.tracker_angle=[-40,40] --pv_array.modules_per_span=10
+to run: python -u pvade_main.py --input_file=input/duramat_case_study.yaml --domain.l_char=3.0 --pv_array.torque_tube_radius=0.1 --pv_array.torque_tube_separation=0.2 --pv_array.stream_rows=2 --pv_array.tracker_angle=[20,20] --pv_array.modules_per_span=10
 """
 
 

From 5bbb1dda6ce7ed0576d8e6e6617ba8c0fb1dde49 Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Fri, 19 Dec 2025 16:45:40 -0700
Subject: [PATCH 22/26] fix initialization of total torque

---
 pvade/fluid/FlowManager.py            | 32 ++++++++++++++++++++++-----
 pvade/structure/ElasticityAnalysis.py | 26 ++++++++++++++++++----
 2 files changed, 48 insertions(+), 10 deletions(-)

diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index 197a2d4..accb25c 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -165,6 +165,17 @@ def build_forms(self, domain, params):
 
         """
 
+        # initialize total torque on each panel to zero
+        for panel_id in range(
+            int(params.pv_array.stream_rows * params.pv_array.span_rows)
+        ): 
+            attr_name = f"total_torque_panel_{panel_id:.0f}"
+            setattr(self, attr_name, dolfinx.fem.Constant(domain.fluid.msh, 0.0))
+
+            for module_id in range(params.pv_array.modules_per_span+1):
+                attr_d_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{module_id:.0f}"
+                setattr(self, attr_d_name, dolfinx.fem.Constant(domain.fluid.msh, 0.0))
+
         # Define fluid properties
         if self.ndim == 2:
             self.dpdx = dolfinx.fem.Constant(domain.fluid.msh, (0.0, 0.0))
@@ -624,7 +635,8 @@ def compute_double_integral_panel_torques(self, domain, params):
         
         for panel_id in range(
             int(params.pv_array.stream_rows * params.pv_array.span_rows)
-        ):
+        ): 
+
             whole_top_surface_facet_index = [] # facet index of the whole panel top surface in a row
             whole_bot_surface_facet_index = [] # facet index of the whole panel bottom surface in a row
 
@@ -651,7 +663,9 @@ def compute_double_integral_panel_torques(self, domain, params):
             
             attr_name = f"total_torque_panel_{panel_id:.0f}"
 
-            setattr(self, attr_name, total_torque_on_panel_array)
+            total_torque_constant = getattr(self, attr_name)
+
+            total_torque_constant.value = total_torque_on_panel_array
 
             print('check total torque definiition', self.total_torque_panel_0)
 
@@ -724,8 +738,12 @@ def compute_double_integral_panel_torques(self, domain, params):
 
             torque_double_integral = 0.0
 
+
             attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_0"
-            setattr(self, attr_name, torque_double_integral)
+            # setattr(self, attr_name, torque_double_integral)
+
+            double_integral_total_torque_panel_constant = getattr(self, attr_name)
+            double_integral_total_torque_panel_constant.value = torque_double_integral
 
             single_integral_function_top_fluid = dolfinx.fem.Function(self.Q)
             single_integral_function_top_fluid.interpolate(single_integral_function_top)
@@ -734,11 +752,13 @@ def compute_double_integral_panel_torques(self, domain, params):
 
             for connector_id in range(params.pv_array.modules_per_span):
                 torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_top_fluid * ds_fluid(
-                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{params.panel_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
+                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
                 torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_bot_fluid * ds_fluid(
-                    domain.domain_markers[f"panel_bot_{panel_id:.0f}_{params.panel_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
+                    domain.domain_markers[f"panel_bot_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
                 attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{connector_id+1:.0f}"
-                setattr(self, attr_name, torque_double_integral)
+                
+                double_integral_total_torque_panel_constant = getattr(self, attr_name)
+                double_integral_total_torque_panel_constant.value = torque_double_integral
 
                 # double_integral_total_torque_panel_0_0 is the double integral of first panel array at the most back connector (maximum y), 
                 # double_integral_total_torque_panel_0_10 is the double integral of first panel array at the most front connector (smallest y)
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 38fecea..0e80a82 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -115,7 +115,7 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
             # the time loop. Should we solve fluid then structure?
             
             total_torque_on_this_panel_name = f"total_torque_panel_{panel_id:.0f}"
-            total_torque_on_this_panel = getattr(flow, total_torque_on_this_panel_name)
+            total_torque_on_this_panel = getattr(flow, total_torque_on_this_panel_name).value
 
             theo_matrix[params.pv_array.modules_per_span, :] = 1/Gp/Ipp
             theo_vector[params.pv_array.modules_per_span] = -total_torque_on_this_panel/Gp/Ipp
@@ -129,7 +129,7 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
 
                 # contribution from C0 to vector:
             T_double_integral_at_fixed_location_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{params.pv_array.fixed_location:.0f}"
-            T_double_integral_at_fixed_location = getattr(flow, T_double_integral_at_fixed_location_name)
+            T_double_integral_at_fixed_location = getattr(flow, T_double_integral_at_fixed_location_name).value
             theo_vector[:params.pv_array.modules_per_span] += -T_double_integral_at_fixed_location/Gp/Ipp
             
                 # contribution of panels to matrix
@@ -148,7 +148,7 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
             T_double_integral_array = []
             for i in range(params.pv_array.modules_per_span+1):
                 name = f"double_integral_total_torque_panel_{panel_id:.0f}_{i:.0f}"
-                T_double_integral_array.append(getattr(flow, name))
+                T_double_integral_array.append(getattr(flow, name).value)
             T_double_integral_array = np.array(T_double_integral_array)
             theo_vector[:params.pv_array.modules_per_span] += np.delete(T_double_integral_array, params.pv_array.fixed_location)/Gp/Ipp
 
@@ -613,6 +613,24 @@ def P_connector(u):
             "dx", domain=domain.structure.msh, subdomain_data=domain.structure.cell_tags
         )
 
+        # print(domain.domain_markers["modules"])
+        # print(np.unique(domain.structure.cell_tags.values))
+
+        print(domain.domain_markers["modules"]["idx"])
+
+        print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain.structure.msh, 1.0) * dx_structure(167))))
+        print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain.structure.msh, 1.0) * dx_structure(168))))
+        
+        print(id(domain.structure.msh))
+        print(id(self.V.mesh))
+        print(id(domain.structure.cell_tags.mesh))
+
+        print(np.unique(domain.structure.cell_tags.values))
+        print(domain.domain_markers["modules"]["idx"])
+        print(domain.domain_markers["connectors"]["idx"])
+        
+        exit()
+
         # To Do: how to differentiate the connector part and the panel part, dx_connector and dx_panel
         self.res = (
             m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * ufl.dx
@@ -630,7 +648,7 @@ def P_connector(u):
             for i in range(params.pv_array.modules_per_span+1):
                 name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
                 K_springs = dolfinx.fem.Constant(domain.structure.msh, float(getattr(self, name_K)))
-                self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*self.ds(domain.domain_markers[f"block_bottom_{panel_id:.0f}_{i:.0f}"])
+                self.res -= ufl.dot(K_springs * self.u_, self.z_unit_vector)*self.ds(domain.domain_markers[f"block_bottom_{panel_id:.0f}_{i:.0f}"]["idx"])
        
         # self.a = dolfinx.fem.form(ufl.lhs(res))
         # self.L = dolfinx.fem.form(ufl.rhs(res))

From 527e0820f4c765c292e3d6d40911e6d8a18fbb7d Mon Sep 17 00:00:00 2001
From: xinhe2205 <45882470+xinhe2205@users.noreply.github.com>
Date: Tue, 23 Dec 2025 16:18:47 -0700
Subject: [PATCH 23/26] fix the Assertion Error

The variational form should include only ufl.dx or dx_structure. It cause assertion error if we use ufl.dx and dx_structure at the same time.
---
 pvade/structure/ElasticityAnalysis.py | 21 ++-------------------
 1 file changed, 2 insertions(+), 19 deletions(-)

diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 0e80a82..e382e0c 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -613,28 +613,11 @@ def P_connector(u):
             "dx", domain=domain.structure.msh, subdomain_data=domain.structure.cell_tags
         )
 
-        # print(domain.domain_markers["modules"])
-        # print(np.unique(domain.structure.cell_tags.values))
-
-        print(domain.domain_markers["modules"]["idx"])
-
-        print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain.structure.msh, 1.0) * dx_structure(167))))
-        print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain.structure.msh, 1.0) * dx_structure(168))))
         
-        print(id(domain.structure.msh))
-        print(id(self.V.mesh))
-        print(id(domain.structure.cell_tags.mesh))
-
-        print(np.unique(domain.structure.cell_tags.values))
-        print(domain.domain_markers["modules"]["idx"])
-        print(domain.domain_markers["connectors"]["idx"])
-        
-        exit()
-
         # To Do: how to differentiate the connector part and the panel part, dx_connector and dx_panel
         self.res = (
-            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
+            m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) * dx_structure
+            + c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) * dx_structure
             + k_nominal(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * dx_structure(domain.domain_markers["modules"]["idx"])
             + k_nominal_connector(self.avg(self.u_old, self.u, self.alpha_f), self.u_) * dx_structure(domain.domain_markers["connectors"]["idx"])
             - structure.rho * ufl.inner(self.f, self.u_) * dx_structure(domain.domain_markers["modules"]["idx"])

From 0681dd5cccb1b6f6768706f24f2b4fa373005182 Mon Sep 17 00:00:00 2001
From: xinhe2205 <45882470+xinhe2205@users.noreply.github.com>
Date: Tue, 23 Dec 2025 16:37:49 -0700
Subject: [PATCH 24/26] Initialize spring stiffness

we are solving the structure first. For the first step, no stress is applied to the structure, so the rotation is zero from first time step, which return nan stiffness value. Initialize the spring stiffness to zero at the first time step (or phi=0).
---
 pvade/structure/ElasticityAnalysis.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index e382e0c..038ee3e 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -226,8 +226,11 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
             else:
                 tracker_angle_rad = np.radians(params.pv_array.tracker_angle)
 
-
-            K = np.abs(12*(np.delete(R_reaction_torque, params.pv_array.fixed_location))/((params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord)**3)/np.cos(tracker_angle_rad+phi)/(np.sin(tracker_angle_rad+phi)-np.sin(tracker_angle_rad))/(params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord/2))
+            if np.linalg.norm(phi) == 0:
+                print('yes 0')
+                K = np.zeros(params.pv_array.modules_per_span)
+            else:
+                K = np.abs(12*(np.delete(R_reaction_torque, params.pv_array.fixed_location))/((params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord)**3)/np.cos(tracker_angle_rad+phi)/(np.sin(tracker_angle_rad+phi)-np.sin(tracker_angle_rad))/(params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord/2))
 
             # assume K is 0 at the fixed connector
             K = np.flip(np.insert(K, params.pv_array.fixed_location, 0))

From e5423e2e20b2b0869534aa063387f1efa6c3515f Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Wed, 7 Jan 2026 09:58:41 -0700
Subject: [PATCH 25/26] verified

---
 pvade/fluid/FlowManager.py            | 98 +++++++++++++++++++--------
 pvade/structure/ElasticityAnalysis.py |  5 --
 2 files changed, 69 insertions(+), 34 deletions(-)

diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index accb25c..6969054 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -625,55 +625,93 @@ def _all_interior_surfaces(x):
     #         setattr(self, attr_name, total_torque)
 
     def compute_double_integral_panel_torques(self, domain, params):
-        torque_function_on_panels = dolfinx.fem.Function(self.Q)
-
-        torque_function_on_panels.x.array[:] = self.spatial_X_coords.x.array[:]*self.traction[2].x.array[:]-self.spatial_Z_coords.x.array[:]*self.traction[0].x.array[:]
 
         ds_fluid = ufl.Measure(
             "ds", domain=domain.fluid.msh, subdomain_data=domain.fluid.facet_tags
         )
-        
+
         for panel_id in range(
             int(params.pv_array.stream_rows * params.pv_array.span_rows)
-        ): 
-
-            whole_top_surface_facet_index = [] # facet index of the whole panel top surface in a row
-            whole_bot_surface_facet_index = [] # facet index of the whole panel bottom surface in a row
+        ):
+            whole_top_surface_facet_index = []  # facet index of the whole panel top surface in a row
+            whole_bot_surface_facet_index = []  # facet index of the whole panel bottom surface in a row
 
             for module_id in range(params.pv_array.modules_per_span):
-                
                 whole_top_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist()
                 whole_bot_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist()
 
-            whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_top_surface_facet_index))
-            whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_bot_surface_facet_index))
+            whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_top_surface_facet_index, dtype=np.int32))
+            whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_bot_surface_facet_index, dtype=np.int32))
+            
+            spatial_X_coords_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="X_coords_top")
+            spatial_X_coords_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="X_coords_bot")
+            spatial_Z_coords_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="Z_coords_top")
+            spatial_Z_coords_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="Z_coords_bot")
+            spatial_X_coords_top.interpolate(self.spatial_X_coords)
+            spatial_X_coords_bot.interpolate(self.spatial_X_coords)
+            spatial_Z_coords_top.interpolate(self.spatial_Z_coords)
+            spatial_Z_coords_bot.interpolate(self.spatial_Z_coords)
+            """
+            the torque should be relative to central of tube, traction x times z (z=0 at center tube) + traction z times x, (x=0 at center tube)
+            for z, how we should put the tube center at zero? It is not minus elevation
+            """
+            traction_z_array_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="traction_z_top")
+            traction_x_array_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="traction_x_top")
+
+            traction_z_array_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="traction_z_bot")
+            traction_x_array_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="traction_x_bot")
+            
+            traction_z_array_top.interpolate(self.traction[2])
+            traction_x_array_top.interpolate(self.traction[0])
+            traction_z_array_bot.interpolate(self.traction[2])
+            traction_x_array_bot.interpolate(self.traction[0])
+
+            # traction_z_array_top.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180)
+            # traction_z_array_bot.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180)
+            # traction_x_array_top.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180)
+            # traction_x_array_bot.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180)
+
+            torque_function_on_panels_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)))
+            torque_function_on_panels_top.x.array[:] = (
+                (spatial_X_coords_top.x.array[:] - np.linspace(
+            0.0,
+            params.pv_array.stream_spacing * (params.pv_array.stream_rows - 1),
+            params.pv_array.stream_rows,
+        )[panel_id]) * traction_z_array_top.x.array[:]
+                - (spatial_Z_coords_top.x.array[:] - params.pv_array.elevation) * traction_x_array_top.x.array[:]
+            )
+
+            torque_function_on_panels_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)))
+            torque_function_on_panels_bot.x.array[:] = (
+                (spatial_X_coords_bot.x.array[:] - np.linspace(
+            0.0,
+            params.pv_array.stream_spacing * (params.pv_array.stream_rows - 1),
+            params.pv_array.stream_rows,
+        )[panel_id]) * traction_z_array_bot.x.array[:]
+                - (spatial_Z_coords_bot.x.array[:] - params.pv_array.elevation) * traction_x_array_bot.x.array[:]
+            )
 
             dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh)
             dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh)
 
+            total_torque_on_panel_array = (
+                dolfinx.fem.assemble_scalar(dolfinx.fem.form(torque_function_on_panels_top * dx_top_whole))
+                + dolfinx.fem.assemble_scalar(dolfinx.fem.form(torque_function_on_panels_bot * dx_bot_whole))
+            )
+
             whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1))
-            whole_top_submesh_func = dolfinx.fem.Function(whole_top_submesh_function_space)
-            whole_top_submesh_func.interpolate()(torque_function_on_panels)
             whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1))
-            whole_bot_submesh_func = dolfinx.fem.Function(whole_bot_submesh_function_space)
-            whole_bot_submesh_func.interpolate()(torque_function_on_panels)
 
-            total_torque_on_panel_array = dolfinx.fem.assemble_scalar(dolfinx.fem.form(whole_top_submesh_func * dx_top_whole)) \
-                                        + dolfinx.fem.assemble_scalar(dolfinx.fem.form(whole_bot_submesh_func * dx_bot_whole))
-            
             attr_name = f"total_torque_panel_{panel_id:.0f}"
-
             total_torque_constant = getattr(self, attr_name)
 
             total_torque_constant.value = total_torque_on_panel_array
 
-            print('check total torque definiition', self.total_torque_panel_0)
-
             coords_top = whole_top_submesh_function_space.tabulate_dof_coordinates()
             coords_bot = whole_bot_submesh_function_space.tabulate_dof_coordinates()
 
-            top_integrated_func_along_y = np.zeros_like(whole_top_submesh_func.x.array[:])
-            bot_integrated_func_along_y = np.zeros_like(whole_bot_submesh_func.x.array[:])
+            top_integrated_func_along_y = np.zeros_like(torque_function_on_panels_top.x.array[:])
+            bot_integrated_func_along_y = np.zeros_like(torque_function_on_panels_bot.x.array[:])
 
             # in PVade, the front has the minimum y coordinates, if we integral from large y to small y, it is negative,
             # to keep it consistant as the theoretical model, we need to flip the sign of the integral
@@ -687,7 +725,7 @@ def compute_double_integral_panel_torques(self, domain, params):
 
             ### integral of top surface:
             # Unique x-values (up to numerical tolerance)
-            ys_top = np.unique(np.round(top_coords_sorted[:, 0], decimals=6))
+            ys_top = np.unique(np.round(top_coords_sorted[:, 1], decimals=6))
 
             # get f(x) = integral of torque on top surface from left end to x (x is all unique x coordinates on top surface)
 
@@ -703,7 +741,7 @@ def compute_double_integral_panel_torques(self, domain, params):
                 submesh_top_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_top_surf_submesh, 2, marked_facet_top.astype(np.int32))
                 top_sub_function_space = dolfinx.fem.FunctionSpace(submesh_top_surface_part, ('CG', 1))
                 top_sub_func = dolfinx.fem.Function(top_sub_function_space)
-                top_sub_func.interpolate(torque_function_on_panels)
+                top_sub_func.interpolate(torque_function_on_panels_top)
 
                 dx_top = ufl.Measure("dx", domain=submesh_top_surface_part)
                 integrated_top = dolfinx.fem.assemble_scalar(dolfinx.fem.form(top_sub_func * dx_top))
@@ -724,7 +762,7 @@ def compute_double_integral_panel_torques(self, domain, params):
                 submesh_bot_surface_part, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(whole_bot_surf_submesh, 2, marked_facet_bot.astype(np.int32))
                 bot_sub_function_space = dolfinx.fem.FunctionSpace(submesh_bot_surface_part, ('CG', 1))
                 bot_sub_func = dolfinx.fem.Function(bot_sub_function_space)
-                bot_sub_func.interpolate(torque_function_on_panels)
+                bot_sub_func.interpolate(torque_function_on_panels_bot)
 
                 dx_bot = ufl.Measure("dx", domain=submesh_bot_surface_part)
                 integrated_bot = dolfinx.fem.assemble_scalar(dolfinx.fem.form(bot_sub_func * dx_bot))
@@ -752,9 +790,9 @@ def compute_double_integral_panel_torques(self, domain, params):
 
             for connector_id in range(params.pv_array.modules_per_span):
                 torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_top_fluid * ds_fluid(
-                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
+                    domain.domain_markers[f"panel_top_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/params.pv_array.panel_chord
                 torque_double_integral += dolfinx.fem.assemble_scalar(dolfinx.fem.form(single_integral_function_bot_fluid * ds_fluid(
-                    domain.domain_markers[f"panel_bot_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/self.pvarry.panel_chord
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{params.pv_array.modules_per_span-1-connector_id:.0f}"]["idx"])))/params.pv_array.panel_chord
                 attr_name = f"double_integral_total_torque_panel_{panel_id:.0f}_{connector_id+1:.0f}"
                 
                 double_integral_total_torque_panel_constant = getattr(self, attr_name)
@@ -762,7 +800,7 @@ def compute_double_integral_panel_torques(self, domain, params):
 
                 # double_integral_total_torque_panel_0_0 is the double integral of first panel array at the most back connector (maximum y), 
                 # double_integral_total_torque_panel_0_10 is the double integral of first panel array at the most front connector (smallest y)
-           
+
     def _assemble_system(self, params):
         """Pre-assemble all LHS matrices and RHS vectors
 
@@ -891,6 +929,8 @@ def solve(self, domain, params, current_time):
         self.mesh_vel.vector.array[:] = self.mesh_vel.vector.array / params.solver.dt
         self.mesh_vel.x.scatter_forward()
 
+        self.compute_double_integral_panel_torques(domain, params)
+
         # Calculate the tentative velocity
         self._solver_step_1(params)
 
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index 038ee3e..dde7b95 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -244,12 +244,7 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
 
                 name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
                 setattr(self, name_K, K[i])
-               
-
                 
-
-
-
         # num_panel_right_fixed = params.pv_array.modules_per_span // 2
         # num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed
 

From 25c1fd78718ecaf5789a5844705fdd246d4a9d29 Mon Sep 17 00:00:00 2001
From: He <xhe@xhe-37925s.nrel.gov>
Date: Thu, 8 Jan 2026 14:01:01 -0700
Subject: [PATCH 26/26] add lift and drag force calculation

---
 pvade/fluid/FlowManager.py            | 172 ++++++++++++++++++++++----
 pvade/structure/ElasticityAnalysis.py |   3 +-
 pvade_main.py                         |   1 +
 3 files changed, 149 insertions(+), 27 deletions(-)

diff --git a/pvade/fluid/FlowManager.py b/pvade/fluid/FlowManager.py
index 6969054..cdb97ce 100644
--- a/pvade/fluid/FlowManager.py
+++ b/pvade/fluid/FlowManager.py
@@ -623,6 +623,93 @@ def _all_interior_surfaces(x):
     #         attr_name = f"total_torque_panel_{panel_id:.0f}"
 
     #         setattr(self, attr_name, total_torque)
+        
+        # define integrated force forms
+        self.integrated_force_x_form = []
+        self.integrated_force_y_form = []
+        self.integrated_force_z_form = []
+
+        for panel_id in range(
+            int(params.pv_array.stream_rows * params.pv_array.span_rows)
+        ):
+            # for loc in ["top_0", "bottom_0", "left_0", "right_0"]:
+            #     idx = domain.domain_markers[loc]["idx"]
+            #     s = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1.0*ds_fluid(idx)))
+            #     print(f"loc = {loc}, idx = {idx}, s = {s}")
+
+            self.integrated_force_x_form.append(0)
+            self.integrated_force_y_form.append(0)
+            self.integrated_force_z_form.append(0)
+
+            self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
+                    domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"]
+                )
+            self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
+                    domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"]
+                )
+            self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
+                    domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"]
+                )
+            self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
+                    domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"]
+                )
+            self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
+                    domain.domain_markers[f"panel_left_{panel_id:.0f}"]["idx"]
+                )
+            self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
+                    domain.domain_markers[f"panel_right_{panel_id:.0f}"]["idx"]
+                )
+            
+            if self.ndim == 3:
+                    self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
+                        domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"]
+                    )
+                    self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
+                        domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"]
+                    )
+                    self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
+                        domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"]
+                    )
+                    self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
+                        domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"]
+                    )
+                    self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
+                        domain.domain_markers[f"panel_front_{panel_id:.0f}"]["idx"]
+                    )
+                    self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
+                        domain.domain_markers[f"panel_back_{panel_id:.0f}"]["idx"]
+                    )
+            for module_id in range(params.pv_array.modules_per_span):
+                
+                self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]
+                )
+                self.integrated_force_x_form[-1] += self.traction[0] * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]
+                )
+                self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]
+                )
+                self.integrated_force_y_form[-1] += self.traction[1] * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]
+                )
+                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]
+                )
+                self.integrated_force_z_form[-1] += self.traction[2] * ds_fluid(
+                    domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]
+                )
+                
+
+            self.integrated_force_x_form[-1] = dolfinx.fem.form(
+                self.integrated_force_x_form[-1]
+            )
+            self.integrated_force_y_form[-1] = dolfinx.fem.form(
+                self.integrated_force_y_form[-1]
+            )
+            self.integrated_force_z_form[-1] = dolfinx.fem.form(
+                self.integrated_force_z_form[-1]
+            )
 
     def compute_double_integral_panel_torques(self, domain, params):
 
@@ -635,14 +722,14 @@ def compute_double_integral_panel_torques(self, domain, params):
         ):
             whole_top_surface_facet_index = []  # facet index of the whole panel top surface in a row
             whole_bot_surface_facet_index = []  # facet index of the whole panel bottom surface in a row
-
+            
             for module_id in range(params.pv_array.modules_per_span):
                 whole_top_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_top_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist()
                 whole_bot_surface_facet_index += domain.fluid.facet_tags.find(domain.domain_markers[f"panel_bottom_{panel_id:.0f}_{module_id:.0f}"]["idx"]).tolist()
 
             whole_top_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_top_surface_facet_index, dtype=np.int32))
             whole_bot_surf_submesh, entity_map, vertex_map, geom_map = dolfinx.mesh.create_submesh(domain.fluid.msh, 2, np.array(whole_bot_surface_facet_index, dtype=np.int32))
-            
+
             spatial_X_coords_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="X_coords_top")
             spatial_X_coords_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="X_coords_bot")
             spatial_Z_coords_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="Z_coords_top")
@@ -655,23 +742,62 @@ def compute_double_integral_panel_torques(self, domain, params):
             the torque should be relative to central of tube, traction x times z (z=0 at center tube) + traction z times x, (x=0 at center tube)
             for z, how we should put the tube center at zero? It is not minus elevation
             """
-            traction_z_array_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="traction_z_top")
-            traction_x_array_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)), name="traction_x_top")
 
-            traction_z_array_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="traction_z_bot")
-            traction_x_array_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)), name="traction_x_bot")
+            whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1))
+            whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1))
+
+            traction_z_array_top = dolfinx.fem.Function(whole_top_submesh_function_space, name="traction_z_top")
+            traction_x_array_top = dolfinx.fem.Function(whole_top_submesh_function_space, name="traction_x_top")
+            traction_z_array_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name="traction_z_bot")
+            traction_x_array_bot = dolfinx.fem.Function(whole_bot_submesh_function_space, name="traction_x_bot")
             
-            traction_z_array_top.interpolate(self.traction[2])
-            traction_x_array_top.interpolate(self.traction[0])
-            traction_z_array_bot.interpolate(self.traction[2])
-            traction_x_array_bot.interpolate(self.traction[0])
+            traction_z_fluid = dolfinx.fem.Function(self.Q, name="traction_z_fluid")
+            traction_x_fluid = dolfinx.fem.Function(self.Q, name="traction_x_fluid")
+
+            # traction is expression on whole fluid mesh, we need to
+            # project it to a function defined on whole fluid mesh, then
+            # interpolate it to the top and bot submesh function space.
+            
+            u_inter_fluid = ufl.TrialFunction(self.Q)
+            v_inter_fluid = ufl.TestFunction(self.Q)
+           
+            a_inter_fluid = ufl.inner(u_inter_fluid, v_inter_fluid) * ufl.ds
+            L_inter_fluid_z = ufl.inner(self.traction[2], v_inter_fluid) * ufl.ds
+            L_inter_fluid_x = ufl.inner(self.traction[0], v_inter_fluid) * ufl.ds
+
+            # Assemble and solve
+            A_inter_fluid = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a_inter_fluid))
+            A_inter_fluid.assemble()
+            
+            b_inter_fluid_z = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(L_inter_fluid_z))
+            b_inter_fluid_x = dolfinx.fem.petsc.assemble_vector(dolfinx.fem.form(L_inter_fluid_x))
+
+            solver = PETSc.KSP().create(traction_z_array_top.function_space.mesh.comm)
+            solver.setOperators(A_inter_fluid)
+            solver.setType("cg")
+            solver.getPC().setType("jacobi")
+            solver.setFromOptions()
+            solver.solve(b_inter_fluid_z, traction_z_fluid.vector)
+            traction_z_fluid.x.scatter_forward()
+
+            solver = PETSc.KSP().create(traction_x_array_top.function_space.mesh.comm)
+            solver.setOperators(A_inter_fluid)
+            solver.setType("cg")
+            solver.getPC().setType("jacobi")
+            solver.setFromOptions()
+            solver.solve(b_inter_fluid_x, traction_x_fluid.vector)
+            traction_x_fluid.x.scatter_forward()
+
+            traction_z_array_top.interpolate(traction_z_fluid)
+            traction_x_array_top.interpolate(traction_x_fluid)
+            traction_z_array_bot.interpolate(traction_z_fluid)
+            traction_x_array_bot.interpolate(traction_x_fluid)
 
             # traction_z_array_top.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180)
             # traction_z_array_bot.x.array[:] = 10*np.sin(params.pv_array.tracker_angle[0]*np.pi/180)
             # traction_x_array_top.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180)
             # traction_x_array_bot.x.array[:] = 10*np.cos(params.pv_array.tracker_angle[0]*np.pi/180)
-
-            torque_function_on_panels_top = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1)))
+            torque_function_on_panels_top = dolfinx.fem.Function(whole_top_submesh_function_space)
             torque_function_on_panels_top.x.array[:] = (
                 (spatial_X_coords_top.x.array[:] - np.linspace(
             0.0,
@@ -681,7 +807,7 @@ def compute_double_integral_panel_torques(self, domain, params):
                 - (spatial_Z_coords_top.x.array[:] - params.pv_array.elevation) * traction_x_array_top.x.array[:]
             )
 
-            torque_function_on_panels_bot = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1)))
+            torque_function_on_panels_bot = dolfinx.fem.Function(whole_bot_submesh_function_space)
             torque_function_on_panels_bot.x.array[:] = (
                 (spatial_X_coords_bot.x.array[:] - np.linspace(
             0.0,
@@ -690,7 +816,6 @@ def compute_double_integral_panel_torques(self, domain, params):
         )[panel_id]) * traction_z_array_bot.x.array[:]
                 - (spatial_Z_coords_bot.x.array[:] - params.pv_array.elevation) * traction_x_array_bot.x.array[:]
             )
-
             dx_top_whole = ufl.Measure("dx", domain=whole_top_surf_submesh)
             dx_bot_whole = ufl.Measure("dx", domain=whole_bot_surf_submesh)
 
@@ -698,10 +823,7 @@ def compute_double_integral_panel_torques(self, domain, params):
                 dolfinx.fem.assemble_scalar(dolfinx.fem.form(torque_function_on_panels_top * dx_top_whole))
                 + dolfinx.fem.assemble_scalar(dolfinx.fem.form(torque_function_on_panels_bot * dx_bot_whole))
             )
-
-            whole_top_submesh_function_space = dolfinx.fem.FunctionSpace(whole_top_surf_submesh, ('CG', 1))
-            whole_bot_submesh_function_space = dolfinx.fem.FunctionSpace(whole_bot_surf_submesh, ('CG', 1))
-
+            
             attr_name = f"total_torque_panel_{panel_id:.0f}"
             total_torque_constant = getattr(self, attr_name)
 
@@ -715,7 +837,7 @@ def compute_double_integral_panel_torques(self, domain, params):
 
             # in PVade, the front has the minimum y coordinates, if we integral from large y to small y, it is negative,
             # to keep it consistant as the theoretical model, we need to flip the sign of the integral
-                
+
             # Sort by y first, then z and z (group by y-levels)
             sort_idx_top = np.lexsort((coords_top[:, 0], -coords_top[:, 1]))  # y is sorted, from large to small 
             sort_idx_bot = np.lexsort((coords_bot[:, 0], -coords_bot[:, 1]))  # y is sorted, from large to small 
@@ -728,7 +850,7 @@ def compute_double_integral_panel_torques(self, domain, params):
             ys_top = np.unique(np.round(top_coords_sorted[:, 1], decimals=6))
 
             # get f(x) = integral of torque on top surface from left end to x (x is all unique x coordinates on top surface)
-
+            indicator_top = 0
             for y_value in ys_top:
                 # Mask for points at this y-level
                 y_mask_top = np.isclose(top_coords_sorted[:, 1], y_value, atol=1e-6)
@@ -746,11 +868,11 @@ def compute_double_integral_panel_torques(self, domain, params):
                 dx_top = ufl.Measure("dx", domain=submesh_top_surface_part)
                 integrated_top = dolfinx.fem.assemble_scalar(dolfinx.fem.form(top_sub_func * dx_top))
                 top_integrated_func_along_y[sort_idx_top[y_mask_top]] = integrated_top
-
+                indicator_top += 1
             ### integral of bot surface:
             # Unique y-values (up to numerical tolerance)
             ys_bot = np.unique(np.round(bot_coords_sorted[:, 1], decimals=6))
-
+            indicator_bot = 0
             for y_value in ys_bot:
 
                 # Mask for points at this y-level
@@ -768,6 +890,7 @@ def compute_double_integral_panel_torques(self, domain, params):
                 integrated_bot = dolfinx.fem.assemble_scalar(dolfinx.fem.form(bot_sub_func * dx_bot))
                 bot_integrated_func_along_y[sort_idx_bot[y_mask_bot]] = integrated_bot
 
+                indicator_bot += 1
             single_integral_function_top = dolfinx.fem.Function(whole_top_submesh_function_space, name='single_integral_function_top')  
             single_integral_function_top.x.array[:] = top_integrated_func_along_y
 
@@ -929,11 +1052,10 @@ def solve(self, domain, params, current_time):
         self.mesh_vel.vector.array[:] = self.mesh_vel.vector.array / params.solver.dt
         self.mesh_vel.x.scatter_forward()
 
-        self.compute_double_integral_panel_torques(domain, params)
+        # self.compute_double_integral_panel_torques(domain, params)
 
         # Calculate the tentative velocity
         self._solver_step_1(params)
-
         # Calculate the change in pressure to enforce incompressibility
         self._solver_step_2(params)
 
@@ -957,7 +1079,7 @@ def solve(self, domain, params, current_time):
         # 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
         self.u_k1.x.array[:] = self.u_k.x.array
diff --git a/pvade/structure/ElasticityAnalysis.py b/pvade/structure/ElasticityAnalysis.py
index dde7b95..e585054 100644
--- a/pvade/structure/ElasticityAnalysis.py
+++ b/pvade/structure/ElasticityAnalysis.py
@@ -227,7 +227,6 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
                 tracker_angle_rad = np.radians(params.pv_array.tracker_angle)
 
             if np.linalg.norm(phi) == 0:
-                print('yes 0')
                 K = np.zeros(params.pv_array.modules_per_span)
             else:
                 K = np.abs(12*(np.delete(R_reaction_torque, params.pv_array.fixed_location))/((params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord)**3)/np.cos(tracker_angle_rad+phi)/(np.sin(tracker_angle_rad+phi)-np.sin(tracker_angle_rad))/(params.pv_array.block_chord_div_by_panel_chord * params.pv_array.panel_chord/2))
@@ -244,7 +243,7 @@ def calculate_K_for_Robin_BC(self, domain, flow, params):
 
                 name_K = f"spring_stiffness_{panel_id:.0f}_{i:.0f}"
                 setattr(self, name_K, K[i])
-                
+
         # num_panel_right_fixed = params.pv_array.modules_per_span // 2
         # num_panel_left_fixed = params.pv_array.modules_per_span - num_panel_right_fixed
 
diff --git a/pvade_main.py b/pvade_main.py
index 2af05fe..5e52d56 100644
--- a/pvade_main.py
+++ b/pvade_main.py
@@ -79,6 +79,7 @@ def main(input_file=None):
     for k in range(params.solver.t_steps):
         current_time = (k + 1) * params.solver.dt
 
+
         if (
             structural_analysis
             and (k + 1) % solve_structure_interval_n == 0