Small update (#4)

* pushing latest

* fix a bug with facets

* expand test

opensg/__init__.py

@@ -2,7 +2,7 @@
 # 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
-from opensg.compute_utils import solve_ksp, solve_boun, compute_nullspace, \
+from opensg.compute_utils import solve_ksp, solve_eb_boundary, compute_nullspace, \
     create_gamma_e, R_sig, Dee, sigma, eps, local_boun, local_frame_1D, \
     local_frame, local_frame_1D_manual, local_grad, deri, ddot, gamma_d, \
     construct_gamma_e, gamma_h, gamma_l, A_mat, initialize_array, dof_mapping_quad, generate_boundary_markers

opensg/compute_utils.py

@@ -424,7 +424,23 @@ def A_mat(ABD, e_l, x_l, dx_l, nullspace_l, v_l, dvl, nphases):
     A_l.setNullSpace(nullspace_l) 
     return A_l
 
-def solve_boun(ABD, meshdata, nphases):
+def solve_eb_boundary(ABD, meshdata, nphases):
+    """_summary_
+
+    Parameters
+    ----------
+    ABD : List of ABD matrices for each phase
+        _description_
+    meshdata : _type_
+        _description_
+    nphases : _type_
+        _description_
+
+    Returns
+    -------
+    _type_
+        _description_
+    """
     mesh = meshdata["mesh"]
     # frame = meshdata["frame"]
     frame = local_frame_1D(mesh)

opensg/mesh.py

@@ -56,6 +56,36 @@ def _generate_material_database(self):
         self.material_database = material_database
         return
 
+
+    def _generate_layup_data(self):
+        layup_database = dict()
+
+        mat_names, thick, angle, nlay = [], [], [], []
+        for section in self.sections:
+            name_components = section['elementSet'].split('_')
+            if(len(name_components) > 2):
+                material_name, t, an = [], [], []
+                # if(int(name_components[1]) == self.segment_index):
+                layup = section['layup'] # layup = [[material_name: str, thickness: float, angle:]]
+                nlay.append(len(layup))
+                for layer in layup:
+                    material_name.append(layer[0])     
+                mat_names.append(material_name)
+                for layer in layup:
+                    t.append(layer[1])
+                thick.append(t)
+                for layer in layup:
+                    an.append(layer[2])
+                angle.append(an) 
+
+        layup_database["mat_names"] = mat_names
+        layup_database["thick"] = thick
+        layup_database["angle"] = angle
+        layup_database["nlay"] = nlay
+        self.layup_database = layup_database
+
+        return layup_database
+
 
     def generate_segment_mesh(self, segment_index, filename):
         segment_node_labels = -1 * np.ones(self.num_nodes, dtype=int)
@@ -254,35 +284,37 @@ def _build_local_orientations(self):
 
 
     def _build_boundary_submeshes(self):
-        # extract geometry
         pp = self.mesh.geometry.x
 
         is_left_boundary, is_right_boundary = opensg.generate_boundary_markers(
             min(pp[:,0]), max(pp[:,0]))
 
-        facets_left = dolfinx.mesh.locate_entities_boundary(
+        left_facets = dolfinx.mesh.locate_entities_boundary(
             self.mesh, dim=self.fdim, marker=is_left_boundary)
-        facets_right = dolfinx.mesh.locate_entities_boundary(
+        right_facets = dolfinx.mesh.locate_entities_boundary(
             self.mesh, dim=self.fdim, marker=is_right_boundary)
 
         left_mesh, left_entity_map, left_vertex_map, left_geom_map = dolfinx.mesh.create_submesh(
-            self.mesh, self.fdim, facets_left)
+            self.mesh, self.fdim, left_facets)
         right_mesh, right_entity_map, right_vertex_map, right_geom_map = dolfinx.mesh.create_submesh(
-            self.mesh, self.fdim, facets_right)
+            self.mesh, self.fdim, right_facets)
 
         self.left_submesh = {
             "mesh": left_mesh, 
             "entity_map": left_entity_map, 
             "vertex_map": left_vertex_map, 
-            "geom_map": left_geom_map}
+            "geom_map": left_geom_map,
+            "marker": is_left_boundary}
+            # "facets": left_facets}
 
         self.right_submesh = {
             "mesh": right_mesh, 
             "entity_map": right_entity_map, 
             "vertex_map": right_vertex_map, 
-            "geom_map": right_geom_map}
+            "geom_map": right_geom_map,
+            "marker": is_right_boundary}
+            # "facets": right_facets}
 
-        # generate subdomains
         self.mesh.topology.create_connectivity(2,1)  # (quad mesh topology, boundary(1D) mesh topology)
         cell_of_facet_mesh = self.mesh.topology.connectivity(2,1)
 
@@ -296,12 +328,11 @@ def _build_boundary_submeshes(self):
         #     cell_edge_map.append(c)
         # cell_edge_map = np.ndarray.flatten(np.array(cell_edge_map))
 
-        # 
-        def subdomains_boundary(
-            boundary_mesh, 
-            boundary_marker, 
-            boundary_entity_map
-            ):
+        # generate subdomains
+        def _build_boundary_subdomains(boundary_meshdata):
+            boundary_mesh = boundary_meshdata["mesh"]
+            boundary_entity_map = boundary_meshdata["entity_map"]
+            boundary_marker = boundary_meshdata["marker"]
             boundary_VV = dolfinx.fem.functionspace(
                 boundary_mesh, basix.ufl.element("DG", boundary_mesh.topology.cell_name(), 0, shape=(3, )))
 
@@ -310,13 +341,14 @@ def subdomains_boundary(
             boundary_n = dolfinx.fem.Function(boundary_VV)
 
             boundary_facets = dolfinx.mesh.locate_entities(boundary_mesh, self.fdim, boundary_marker)
+
             # TODO: review the subdomain assingments with akshat
             boundary_subdomains = []
             for i, xx in enumerate(boundary_entity_map):
                 # assign subdomain
                 # 4 is for number of nodes in quad element
                 # NOTE: we should find a different way to do this that doesn't assume quad elements if
-                #    we plan to expand to other elements.
+                #    we plan to expand to other elements. -klb
                 idx = int(np.where(cell_of_facet_mesh.array==xx)[0]/4) 
                 boundary_subdomains.append(self.subdomains.values[idx])
                 # assign orientation
@@ -335,10 +367,8 @@ def subdomains_boundary(
             # Mapping the orinetation data from quad mesh to boundary. The alternative is to use local_frame_1D(self.left_submesh["mesh"]).
             # Either of both can be used in local_boun subroutine 
 
-        self.left_submesh["subdomains"], self.left_submesh["frame"], self.left_submesh["facets"] = subdomains_boundary(
-            self.left_submesh["mesh"], is_left_boundary, self.left_submesh["entity_map"]) 
-        self.right_submesh["subdomains"], self.right_submesh["frame"], self.right_submesh["facets"] = subdomains_boundary(
-            self.right_submesh["mesh"], is_right_boundary, self.right_submesh["entity_map"])
+        self.left_submesh["subdomains"], self.left_submesh["frame"], self.left_submesh["facets"] = _build_boundary_subdomains(self.left_submesh) 
+        self.right_submesh["subdomains"], self.right_submesh["frame"], self.right_submesh["facets"] = _build_boundary_subdomains(self.right_submesh)
 
         return
 
@@ -347,11 +377,10 @@ def compute_ABD(self):
         ABD_ = []
         for i in range(nphases):
             ABD_.append(opensg.compute_ABD_matrix(
-                i, 
-                thick=self.layup_database["thick"], 
-                nlay=self.layup_database["nlay"], 
-                mat_names=self.layup_database["mat_names"],
-                angle=self.layup_database["angle"],
+                thick=self.layup_database["thick"][i], 
+                nlay=self.layup_database["nlay"][i], 
+                mat_names=self.layup_database["mat_names"][i],
+                angle=self.layup_database["angle"][i],
                 material_database=self.blade_mesh.material_database
                 ))
 

opensg/solve.py

@@ -18,7 +18,7 @@
 ### ABD matrix computation
 # @profile
 # NOTE can pass in thick[ii], nlay[ii], etc instead of the dictionaries
-def compute_ABD_matrix(ii, thick, nlay, angle, mat_names, material_database):
+def compute_ABD_matrix(thick, nlay, angle, mat_names, material_database):
     """Compute the ABD matrix for a composite layup structure
 
     Constructs a local stiffness matrix for a composite laminate
@@ -52,14 +52,14 @@ def compute_ABD_matrix(ii, thick, nlay, angle, mat_names, material_database):
 
     # nodes (1D SG)
     th, s = [0], 0  # Reference starting point
-    for k in thick[ii]:
+    for k in thick:
         s = s - k  # Add the thickness of each layer
         th.append(s)
     points = np.array(th)
 
     # elements
     cell = []
-    for k in range(nlay[ii]):
+    for k in range(nlay):
         cell.append([k, k + 1])
     cellss = np.array(cell)
 
@@ -87,9 +87,9 @@ def compute_ABD_matrix(ii, thick, nlay, angle, mat_names, material_database):
     # Weak form of energy
     F2 = 0
     for j in range(nphases):
-        mat_name = mat_names[ii][j]
+        mat_name = mat_names[j]
         material_props = material_database[mat_name]
-        theta = angle[ii][j]
+        theta = angle[j]
         sigma_val = opensg.compute_utils.sigma(u, material_props, theta, Eps = gamma_e[:,0])[0]
         inner_val = inner(sigma_val, opensg.compute_utils.eps(v)[0])
         F2 += inner_val * dx(j)
@@ -111,9 +111,9 @@ def compute_ABD_matrix(ii, thick, nlay, angle, mat_names, material_database):
         # weak form
         F2 = 0
         for j in range(nphases):
-            mat_name = mat_names[ii][j]
+            mat_name = mat_names[j]
             material_props = material_database[mat_name]
-            theta = angle[ii][j]
+            theta = angle[j]
             sigma_val = opensg.compute_utils.sigma(u, material_props, theta, Eps)[0]
             inner_val = inner(sigma_val, opensg.compute_utils.eps(v)[0])
             F2 += inner_val * dx(j)
@@ -137,9 +137,9 @@ def compute_ABD_matrix(ii, thick, nlay, angle, mat_names, material_database):
             # f=dolfinx.fem.form(sum([opensg.compute_utils.Dee(x, u, material_props, theta, Eps)[s,k]*dx(i) for i in range(nphases)])) # Scalar assembly
             f = 0
             for j in range(nphases):
-                mat_name = mat_names[ii][j]
+                mat_name = mat_names[j]
                 material_props = material_database[mat_name]
-                theta = angle[ii][j]
+                theta = angle[j]
                 dee_val = opensg.compute_utils.Dee(x, u, material_props, theta, Eps)[s, k]
                 f += dee_val * dx(j)
             f = dolfinx.fem.form(f)
@@ -189,20 +189,34 @@ def compute_stiffness_EB_blade_segment(
     pp = mesh.geometry.x # point data
     x_min, x_max=min(pp[:,0]), max(pp[:,0])
     # 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"])
+    # 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"])
+    # 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_r = dolfinx.fem.functionspace(mesh, basix.ufl.element(
+        "S", mesh.topology.cell_name(), 2, shape = (3, )))
 
     l_submesh["nullspace"] = opensg.compute_nullspace(V_l)
     r_submesh["nullspace"] = opensg.compute_nullspace(V_r)
 
-    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)
+    # 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.solve_boun(ABD, l_submesh, nphases)
-    V0_r = opensg.solve_boun(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
 
     # 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)
 
@@ -298,6 +312,28 @@ def compute_stiffness_EB_blade_segment(
     print(np.around(D_eff))
 
     return D_eff
+
+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
+    ):
+
+    # 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"])
+
+    l_submesh["nullspace"] = opensg.compute_nullspace(V_l)
+    r_submesh["nullspace"] = opensg.compute_nullspace(V_r)
+
+    V0_l = opensg.solve_eb_boundary(ABD, l_submesh, nphases)
+    V0_r = opensg.solve_eb_boundary(ABD, r_submesh, nphases)
+
+    return V0_l, V0_r
 
 
 def compute_timo_boun(ABD, mesh, subdomains, frame, nullspace, sub_nullspace, nphases):

opensg/tests/test_workflow.py

@@ -2,6 +2,8 @@
 from os.path import abspath, dirname, join
 import numpy as np
 import opensg
+import filecmp
+
 
 example_data_path = ""
 testdir = dirname(abspath(str(__file__)))
@@ -17,8 +19,30 @@ def setUp(self):
 
         self.blade_mesh = opensg.BladeMesh(mesh_data)
         return super().setUp()
+
+    def test_segment_creation(self):
+        mesh_filename = "section.msh"
+        expected_mesh_filename = join(datadir, 'test_section.msh')
+
+        _ = self.blade_mesh.generate_segment_mesh(segment_index=1, filename=mesh_filename)
+
+        assert filecmp.cmp(mesh_filename, expected_mesh_filename)
+
+    def test_ABD_matrix(self):
+        section_mesh = self.blade_mesh.generate_segment_mesh(segment_index=1, filename="section.msh")
+        abd = section_mesh.compute_ABD()
+
+        abd_concat = np.concat(abd)
+        # np.savetxt(join(datadir, 'abd_test.txt'), abd_concat) # use to reset ground truth
+
+        expected_abd = np.loadtxt(join(datadir, 'abd_test.txt'))
+        assert np.isclose(abd_concat, expected_abd).all()
+
+    def test_timo_stiffness(self):
+        pass
+
 
-    def test_example(self):
+    def test_EB_stiffness(self):
         # load expected values
         # expected_ABD = np.loadtxt("")
         # expected_stiffness = np.loadtxt("")
@@ -27,8 +51,22 @@ def test_example(self):
         ABD = section_mesh.compute_ABD()
         # assert np.isclsoe(ABD, expected_ABD).all()
 
-        stiffness_matrix = section_mesh.compute_stiffness_EB(ABD)
-        assert stiffness_matrix.shape == (4,4)
-        # assert np.isclsoe(stiffness_matrix, expected_stiffness).all()
+        computed_stiffness_matrix = section_mesh.compute_stiffness_EB(ABD)
+        # np.savetxt(join(datadir, 'stiffness_test.txt'), computed_stiffness_matrix) # use to reset ground truth
+
+        assert computed_stiffness_matrix.shape == (4,4)
+
+        expected_stiffness_matrix = np.loadtxt(join(datadir, 'stiffness_test.txt'))
+        assert np.isclose(computed_stiffness_matrix, expected_stiffness_matrix).all()
+
+
 
+"""
+[[ 5.6138e+09  7.8930e+07 -9.1787e+07 -6.2902e+08]
+ [ 7.8930e+07  2.2724e+10 -6.0954e+08  2.0795e+08]
+ [-9.1787e+07 -6.0954e+08  1.0064e+10  9.9959e+08]
+ [-6.2902e+08  2.0795e+08  9.9959e+08  1.2617e+10]]
+"""
 
+if __name__ == "__main__":
+    unittest.main()