Boundaries (#5)

* pushing latest

* fix a bug with facets

* removing unneccesary computations in segment stiffness

* add function for boundary stiffness

* compute nphases from length of ABD

* add function to API

* add new script for boundary computations

* add some docstring information

* update test

* add back in missing nullspaces

* cleanup

* expand test

* fix bug where wrong mesh was used to initialize nullspaces

* add testing files

examples/compute_eb_blade_boundaries.py

@@ -0,0 +1,23 @@
+from os.path import join
+
+import opensg
+import numpy as np
+import time
+
+
+mesh_yaml = join("data", "bar_urc_shell_mesh.yaml")
+mesh_data = opensg.load_yaml(mesh_yaml)
+
+blade_mesh = opensg.BladeMesh(mesh_data)
+stiffness_matrices = []
+for i in range(1, 10):
+
+    segment_mesh = blade_mesh.generate_segment_mesh(segment_index=i, filename="section.msh")
+
+    ABD = segment_mesh.compute_ABD()
+
+    stiffness_matrix_l, stiffness_matrix_r = segment_mesh.compute_stiffness_EB_boundary(ABD)
+
+    stiffness_matrices.append((stiffness_matrix_l, stiffness_matrix_r))
+
+pause

opensg/__init__.py

@@ -1,4 +1,4 @@
-from opensg.solve import compute_ABD_matrix, compute_timo_boun, compute_stiffness_EB_blade_segment
+from opensg.solve import compute_ABD_matrix, compute_timo_boun, compute_stiffness_EB_blade_segment, compute_eb_blade_segment_boundary
 # from kirklocal.timo import local_frame_1D, directional_derivative, local_grad, ddot
 from opensg.io import load_yaml, write_yaml
 from opensg.mesh import BladeMesh

opensg/compute_utils.py

@@ -415,16 +415,41 @@ def local_boun(mesh, frame, subdomains):
 
     return frame, V, dv, v, x, dx
 
+def A_mat(ABD, e_l, x_l, dx_l, nullspace_l, v_l, dvl, nphases):
+    """Assembly matrix
 
+    Parameters
+    ----------
+    ABD : _type_
+        _description_
+    e_l : _type_
+        _description_
+    x_l : _type_
+        _description_
+    dx_l : _type_
+        _description_
+    nullspace_l : _type_
+        _description_
+    v_l : _type_
+        _description_
+    dvl : _type_
+        _description_
+    nphases : _type_
+        _description_
 
-def A_mat(ABD, e_l, x_l, dx_l, nullspace_l, v_l, dvl, nphases):
+    Returns
+    -------
+    _type_
+        _description_
+    """
     F2 = sum([dot(dot(as_tensor(ABD[i]),gamma_h(e_l,x_l,dvl)), gamma_h(e_l,x_l,v_l))*dx_l(i) for i in range(nphases)])   
     A_l = assemble_matrix(form(F2))
     A_l.assemble()
     A_l.setNullSpace(nullspace_l) 
     return A_l
 
-def solve_eb_boundary(ABD, meshdata, nphases):
+def solve_eb_boundary(ABD, meshdata):
+
     """_summary_
 
     Parameters
@@ -441,6 +466,10 @@ def solve_eb_boundary(ABD, meshdata, nphases):
     _type_
         _description_
     """
+
+    # assert len(ABD) == nphases
+    nphases = len(ABD)
+
     mesh = meshdata["mesh"]
     # frame = meshdata["frame"]
     frame = local_frame_1D(mesh)
@@ -475,12 +504,35 @@ def initialize_array(V_l):
     V1s = np.zeros((xxx,4))
     return V0,Dle,Dhe,Dhd,Dld,D_ed,D_dd,D_ee,V1s
 
-# dof mapping makes solved unknown value w_l(Function(V_l)) assigned to v2a (Function(V)). 
-# The boundary of wind blade mesh is a 1D curve. The facet/edge number is obtained from cell to edge connectivity (conn3) showed in subdomain subroutine.
-# The same facet/edge number of extracted mesh_l (submesh) is obtaine din entity_mapl (gloabl mesh number). refer how submesh was generated.
-#Therefore, once identifying the edge number being same for global(mesh)&boundary mesh(mesh_l), we equate the dofs and store w_l to v2a.
-# The dofs can be verified by comparing the coordinates of local and global dofs if required. 
+
 def dof_mapping_quad(V, v2a, V_l, w_ll, boundary_facets_left, entity_mapl):
+    """dof mapping makes solved unknown value w_l(Function(V_l)) assigned to v2a (Function(V)). 
+    The boundary of wind blade mesh is a 1D curve. The facet/edge number is obtained from cell to edge connectivity (conn3) showed in subdomain subroutine.
+    The same facet/edge number of extracted mesh_l (submesh) is obtaine din entity_mapl (gloabl mesh number). refer how submesh was generated.
+    Therefore, once identifying the edge number being same for global(mesh)&boundary mesh(mesh_l), we equate the dofs and store w_l to v2a.
+    The dofs can be verified by comparing the coordinates of local and global dofs if required. 
+
+    Parameters
+    ----------
+    V : _type_
+        _description_
+    v2a : _type_
+        _description_
+    V_l : _type_
+        _description_
+    w_ll : 1D array (len(4))
+        Fluctuating function data for case p
+    boundary_facets_left : _type_
+        _description_
+    entity_mapl : _type_
+        _description_
+
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
     dof_S2L = []
     deg = 2
     for i,xx in enumerate(entity_mapl):

opensg/mesh.py

@@ -417,6 +417,12 @@ def compute_stiffness_EB(self, ABD):
 
         return D_eff
 
+    def compute_stiffness_EB_boundary(self, ABD):
+        m_l = opensg.compute_eb_blade_segment_boundary(ABD, self.left_submesh)
+        m_r = opensg.compute_eb_blade_segment_boundary(ABD, self.right_submesh)
+
+        return m_l, m_r
+
     def compute_boundary_stiffness_timo(self, ABD):
 
         left_stiffness = opensg.compute_timo_boun(

opensg/solve.py

@@ -164,8 +164,8 @@ def compute_stiffness_EB_blade_segment(
 
     Parameters
     ----------
-    ABD : _type_
-        _description_
+    ABD : list[array]
+        list of ABD matrices for each phase
     mesh : _type_
         _description_
     frame : _type_
@@ -185,39 +185,23 @@ def compute_stiffness_EB_blade_segment(
     tdim = mesh.topology.dim
     fdim = tdim - 1
     nphases = max(subdomains.values[:]) + 1
-
+    assert nphases == len(ABD)
     pp = mesh.geometry.x # point data
     x_min, x_max=min(pp[:,0]), max(pp[:,0])
-    # Initialize terms
-    # NOTE: only V_l/V_r used so replacing for now
-    # e_l, V_l, dvl, v_l, x_l, dx_l = opensg.local_boun(
-    #     l_submesh["mesh"], l_submesh["frame"], l_submesh["subdomains"])
-
-    # e_r, V_r, dvr, v_r, x_r, dx_r = opensg.local_boun(
-    #     r_submesh["mesh"], r_submesh["frame"], r_submesh["subdomains"])
-
+
     # Initialize nullspaces
-    V_l = dolfinx.fem.functionspace(mesh, basix.ufl.element(
-        "S", mesh.topology.cell_name(), 2, shape = (3, )))
+    V_l = dolfinx.fem.functionspace(l_submesh["mesh"], basix.ufl.element(
+        "S", l_submesh["mesh"].topology.cell_name(), 2, shape = (3, )))
 
-    V_r = dolfinx.fem.functionspace(mesh, basix.ufl.element(
-        "S", mesh.topology.cell_name(), 2, shape = (3, )))
+    V_r = dolfinx.fem.functionspace(r_submesh["mesh"], basix.ufl.element(
+        "S", r_submesh["mesh"].topology.cell_name(), 2, shape = (3, )))
 
-    l_submesh["nullspace"] = opensg.compute_nullspace(V_l)
-    r_submesh["nullspace"] = opensg.compute_nullspace(V_r)
+    # Compute boundaries
+    # NOTE: not the stiffness matrix
 
-    # NOTE: unused code
-    # A_l = opensg.A_mat(ABD, e_l, x_l, dx_l, l_submesh["nullspace"],v_l, dvl, nphases)
-    # A_r = opensg.A_mat(ABD, e_r, x_r, dx_r, r_submesh["nullspace"],v_r, dvr, nphases)
-
-    V0_l = opensg.solve_eb_boundary(ABD, l_submesh, nphases)
-    V0_r = opensg.solve_eb_boundary(ABD, r_submesh, nphases)
-
-    # V0_l = opensg.compute_eb_blade_segment_boundary(ABD, l_submesh, nphases)
-    # V0_r = opensg.compute_eb_blade_segment_boundary(ABD, r_submesh, nphases)
-
-    # compute_eb_blade_segment_boundary
-
+    V0_l = compute_eb_blade_segment_boundary(ABD, l_submesh)
+    V0_r = compute_eb_blade_segment_boundary(ABD, r_submesh)
+
     # The local_frame_l(submesh["mesh"]) can be replaced with frame_l, if we want to use mapped orientation from given direction cosine matrix (orien mesh data-yaml)
 
     # Quad mesh
@@ -264,7 +248,8 @@ def compute_stiffness_EB_blade_segment(
     for p in range(4): # 4 load cases meaning 
         # Boundary 
         v2a = Function(V)
-        v2a = opensg.dof_mapping_quad(V, v2a,V_l,V0_l[:,p], l_submesh["facets"], l_submesh["entity_map"])
+        # NOTE: V0 first dimension is degrees of freedom and second dimension is load cases
+        v2a = opensg.dof_mapping_quad(V, v2a, V_l, V0_l[:,p], l_submesh["facets"], l_submesh["entity_map"])
         # NOTE: does this second call overwrite previous, or is the data combined? -klb
         v2a = opensg.dof_mapping_quad(V, v2a,V_r,V0_r[:,p], r_submesh["facets"], r_submesh["entity_map"])  
 
@@ -315,25 +300,22 @@ def compute_stiffness_EB_blade_segment(
 
 def compute_eb_blade_segment_boundary(
     ABD, # array
-    nphases, 
-    l_submesh, # dictionary with mesh data for l boundary
-    r_submesh # dictionary with mesh data for r boundary
+    submeshdata, # dictionary with mesh data for l boundary
+    # nphases
+    # r_submesh # dictionary with mesh data for r boundary
     ):
 
     # Initialize terms
-    e_l, V_l, dvl, v_l, x_l, dx_l = opensg.local_boun(
-        l_submesh["mesh"], l_submesh["frame"], l_submesh["subdomains"])
-
-    e_r, V_r, dvr, v_r, x_r, dx_r = opensg.local_boun(
-        r_submesh["mesh"], r_submesh["frame"], r_submesh["subdomains"])
+    # e, V, dv, v, x, dx = opensg.local_boun(
+    #     submeshdata["mesh"], submeshdata["frame"], submeshdata["subdomains"])
+    V = dolfinx.fem.functionspace(submeshdata["mesh"], basix.ufl.element(
+        "S", submeshdata["mesh"].topology.cell_name(), 2, shape = (3, )))
 
-    l_submesh["nullspace"] = opensg.compute_nullspace(V_l)
-    r_submesh["nullspace"] = opensg.compute_nullspace(V_r)
+    submeshdata["nullspace"] = opensg.compute_nullspace(V)
 
-    V0_l = opensg.solve_eb_boundary(ABD, l_submesh, nphases)
-    V0_r = opensg.solve_eb_boundary(ABD, r_submesh, nphases)
+    V0 = opensg.solve_eb_boundary(ABD, submeshdata)
 
-    return V0_l, V0_r
+    return V0
 
 
 def compute_timo_boun(ABD, mesh, subdomains, frame, nullspace, sub_nullspace, nphases):
@@ -353,7 +335,7 @@ def compute_timo_boun(ABD, mesh, subdomains, frame, nullspace, sub_nullspace, np
         sub_nullspace.remove(F_l)
         w_l = opensg.compute_utils.solve_ksp(A_l,F_l,V_l)
         Dhe[:,p]=  F_l[:]
-        V0[:,p]= w_l.vector[:]  
+        V0[:,p]= w_l.vector[:]
     D1 = np.matmul(V0.T,-Dhe)   
     for s in range(4):
         for k in range(4):