Added inertial residual in the shell model for modal analysis.

examples/cantilever_plate/cantilever_plate.py

@@ -28,7 +28,7 @@
         "plate_2_10_quad_40_200.xdmf",
         "plate_2_10_quad_80_400.xdmf",]
 
-filename = "../../mesh/mesh-examples/clamped-RM-plate/"+beam[2]
+filename = "./clamped-RM-plate/"+beam[2]
 with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r") as xdmf:
     mesh = xdmf.read_mesh(name="Grid")
 

examples/pegasus/pegasus_wing.py

@@ -4,6 +4,9 @@
 Note: to run the example with the mesh files associated, you need to
 have `git lfs` installed to download the actual mesh files. Please
 refer to instructions on their official website at https://git-lfs.github.com/
+
+The concentrated forces are from AVL simulation
+https://web.mit.edu/drela/Public/web/avl/
 -----------------------------------------------------------
 """
 from dolfinx.io import XDMFFile
@@ -27,32 +30,34 @@
 nel = shell_mesh.topology.index_map(shell_mesh.topology.dim).size_local
 nn = shell_mesh.topology.index_map(0).size_local
 
-E_val = 6.8E10 # unit: Pa (N/m^2)
-nu_val = 0.35
+# Aluminum 7050
+E_val = 6.9E10 # unit: Pa (N/m^2)
+nu_val = 0.327
 h_val = 3E-3 # overall thickness (unit: m)
-
+rho = 2700
 # Scaled body force
-f_d = 10. # force per unit area (unit: N/m^2)
-
+# f_d = Constant(shell_mesh,(0.,0.,-rho*9.81)) # force per unit area (unit: N/m^2)
+f_d = ufl.as_vector([0.,0.,-rho*h_val*9.81])
 E = Constant(shell_mesh,E_val) # Young's modulus
 nu = Constant(shell_mesh,nu_val) # Poisson ratio
 
 # ################## Constant thickness  ###################
 # h = Constant(shell_mesh,h_val) # Shell thickness
 
-################### Varying thickness distribution ###################
-hrcs = np.reshape(np.loadtxt(path+'pegasus_t_med_SI.csv'),(nn,1))
-h_nodal = np.arange(nn)
-node_indices = shell_mesh.geometry.input_global_indices
-h_array = sortIndex(hrcs, h_nodal, node_indices)
-## Apply element-wise thickness ###
-VT = FunctionSpace(shell_mesh, ("CG", 1))
-h_cg1 = Function(VT)
-h_cg1.vector.setArray(h_array)
-
+# ################### Varying thickness distribution ###################
+# hrcs = np.reshape(np.loadtxt(path+'pegasus_t_med_SI.csv'),(nn,1))
+# h_nodal = np.arange(nn)
+# node_indices = shell_mesh.geometry.input_global_indices
+# h_array = sortIndex(hrcs, h_nodal, node_indices)
+# ## Apply element-wise thickness ###
+# VT = FunctionSpace(shell_mesh, ("CG", 1))
+# h_cg1 = Function(VT)
+# h_cg1.vector.setArray(h_array)
+#
 V0 = FunctionSpace(shell_mesh, ('DG', 0))
 h = Function(V0)
-project(h_cg1, h, lump_mass=False)
+h.vector.set(h_val)
+# project(h_cg1, h, lump_mass=False)
 # f = as_vector([0,0,f_d]) # Body force per unit area
 
 
@@ -67,18 +72,18 @@
 #                shear_deg=3
                 )
 
-# VE1 = VectorElement("Lagrange",shell_mesh.ufl_cell(),1)
-# WE = MixedElement([VE1,VE1])
-# W = FunctionSpace(shell_mesh,WE)
-
 W = element.W
 w = Function(W)
 dx_inplane, dx_shear = element.dx_inplane, element.dx_shear
 
 ################# Apply concentrated forces #################
+#
+# x0 = np.load("coords_SI.npy")
+# f_c = np.load("loads_SI.npy")
 
-x0 = np.load("coords_SI.npy")
-f_c = np.load("loads_SI.npy")
+x0 = np.load("coords_beam.npy")
+# f_c = np.load("loads_1g_beam.npy")
+f_c = np.load("loads_2p5g_n1g_beam.npy")
 
 from FSI_coupling.VLM_sim_handling import *
 from FSI_coupling.shellmodule_utils import *
@@ -95,11 +100,13 @@
 
 
 # Define force functions and aero-elastic coupling object ########
-print("x0 shape:", np.shape(x0))
-coupling_obj = FEniCSx_concentrated_load_coupling(shell_mesh, x0[:,0,:],
-                    W, RBF_width_par=.5, RBF_func=RadialBasisFunctions.Gaussian)
+# print("x0 shape:", np.shape(x0))
+coupling_obj = FEniCSx_concentrated_load_coupling(shell_mesh, x0,
+                    W, RBF_width_par=2., RBF_func=RadialBasisFunctions.Gaussian)
 # print("G mat shape:", np.shape(coupling_obj.G_mat.map))
-f_array = coupling_obj.compute_dist_solid_force_from_point_load(f_c)
+# f_c_reordered = f_c.reshape(17,3).T.flatten()
+f_c_reordered = f_c.T.flatten()
+f_array = coupling_obj.compute_dist_solid_force_from_point_load(f_c_reordered)
 
 # print(f_array.shape)
 # apply array in function space
@@ -110,36 +117,40 @@
 fx_sum = 0.
 fy_sum = 0.
 fz_sum = 0.
-f_c_sum = np.sum(f_c.reshape(17,3),axis=0)
+f_c_sum = np.sum(f_c,axis=0)
 for i in range(nn):
     fx_sum += f_array[3*i]
     fy_sum += f_array[3*i+1]
     fz_sum += f_array[3*i+2]
 print("-"*60)
-print("                               Projected      "+"     Original     ")
+print("                               Projected CG1  "+"     Original     ")
 print("-"*60)
 print("Sum of forces in x-direction:", fx_sum, f_c_sum[0])
 print("Sum of forces in y-direction:", fy_sum, f_c_sum[1])
 print("Sum of forces in z-direction:", fz_sum, f_c_sum[2])
 print("-"*60)
 #################################################################
-exit()
+# exit()
 
+f.vector.setArray(f_array) # Nodal force in Newton
 
-f.vector.setArray(0.001*f_array) # Nodal force in Newton
-
+dims = shell_mesh.topology.dim
+# h_mesh = dolfinx.cpp.mesh.h(shell_mesh, dims, range(nel))
+h_mesh = CellDiameter(shell_mesh)
 #### Compute the CLT model from the material properties (for single-layer material)
 material_model = MaterialModel(E=E,nu=nu,h=h)
 elastic_model = ElasticModel(shell_mesh,w,material_model.CLT)
 elastic_energy = elastic_model.elasticEnergy(E, h, dx_inplane,dx_shear)
-F = elastic_model.weakFormResidual(elastic_energy, f)
+F = elastic_model.weakFormResidual(elastic_energy, f_d)
 
 ############ Set the BCs for the airplane model ###################
 
 u0 = Function(W)
 u0.vector.set(0.0)
-
-
+# locate_BC1 = locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]),
+#                                     lambda x: np.less(x[1], 0.9858135))
+# locate_BC2 = locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]),
+#                                     lambda x: np.less(x[1], 0.9858135))
 locate_BC1 = locate_dofs_geometrical((W.sub(0), W.sub(0).collapse()[0]),
                                     lambda x: np.less(x[1], 1e-6))
 locate_BC2 = locate_dofs_geometrical((W.sub(1), W.sub(1).collapse()[0]),
@@ -153,25 +164,43 @@
        ]
 
 ########### Apply the point load #############################
+#
+
 
 f1 = Function(W)
 f1_0,_ = f1.split()
-f1_0.interpolate(f)
+
+
+delta_mpt = Delta_mpt(x0=x0, f_p=f_c)
+f1_0.interpolate(delta_mpt.eval)
+
+
+f1_x = computeNodalDisp(f1.sub(0))[0]
+f1_y = computeNodalDisp(f1.sub(0))[1]
+f1_z = computeNodalDisp(f1.sub(0))[2]
+print("-"*60)
+print("                               Projected CG2   "+"     Original     ")
+print("-"*60)
+print("Sum of forces in x-direction:", np.sum(f1_x), f_c_sum[0])
+print("Sum of forces in y-direction:", np.sum(f1_y), f_c_sum[1])
+print("Sum of forces in z-direction:", np.sum(f1_z), f_c_sum[2])
+print("-"*60)
 
 # Assemble linear system
 a = derivative(F,w)
 L = -F
 A = assemble_matrix(form(a), bcs)
 A.assemble()
 b = assemble_vector(form(L))
+b_array = b.getArray()
 b.setArray(f1.vector.getArray())
 b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
 dolfinx.fem.set_bc(b, bcs)
 
 
 
 ######### Solve the linear system with KSP solver ############
-# solveKSP_mumps(A, b, w.vector)
+solveKSP_mumps(A, b, w.vector)
 
 
 u_mid, _ = w.split()
@@ -268,6 +297,8 @@ def dmax_vmdw(w,vm_stress,dx,rho=200,alpha=None,m=1e-6):
 with XDMFFile(MPI.COMM_WORLD, "solutions_varied_thickness/von_Mises_top.xdmf", "w") as xdmf:
     xdmf.write_mesh(shell_mesh)
     xdmf.write_function(von_Mises_top_func)
+
+project(f1_0,f)
 with XDMFFile(MPI.COMM_WORLD, "solutions_varied_thickness/distributed_force.xdmf", "w") as xdmf:
     xdmf.write_mesh(shell_mesh)
     xdmf.write_function(f)
@@ -363,6 +394,32 @@ def dmax_vmdw(w,vm_stress,dx,rho=200,alpha=None,m=1e-6):
 # np.save("coords_SI.npy", x0_)
 # np.save("loads_SI.npy", f_c_)
 
+# beam results (1g)
+# displacement:  [[ 0.00000000e+00  0.00000000e+00  0.00000000e+00]
+#  [-3.51460786e-06 -1.88959688e-06  7.39187952e-03]
+#  [-7.20382954e-06 -1.17740054e-04  1.32997525e-02]
+#  [-1.82763701e-05 -4.32212137e-04  2.91158859e-02]
+#  [-3.40535986e-05 -8.46548973e-04  4.95745879e-02]
+#  [-5.41226536e-05 -1.34601723e-03  7.37860405e-02]
+#  [-7.80528722e-05 -1.91819116e-03  1.01011492e-01]
+#  [-1.05401197e-04 -2.55236919e-03  1.30636053e-01]
+#  [-1.33567953e-04 -3.07832847e-03  1.62184107e-01]
+#  [-1.55046782e-04 -2.96455543e-03  1.95412848e-01]
+#  [-1.78141362e-04 -2.85014706e-03  2.30185681e-01]
+#  [-2.02571174e-04 -2.73225775e-03  2.66352544e-01]
+#  [-2.28078788e-04 -2.61025630e-03  3.03742686e-01]
+#  [-2.54419451e-04 -2.48472959e-03  3.42163681e-01]
+#  [-2.81352704e-04 -2.35631956e-03  3.81400949e-01]
+#  [-3.08655082e-04 -2.22569927e-03  4.21222111e-01]
+#  [-3.36140400e-04 -2.09350284e-03  4.61392829e-01]
+#  [-3.63682498e-04 -1.96017659e-03  5.01712560e-01]
+#  [-3.77638757e-04 -1.89222007e-03  5.22156462e-01]]
+# stress [5.04713118e+08 4.23850490e+08 3.61016326e+08 2.92575438e+08
+#  2.35824175e+08 1.88632942e+08 1.49302768e+08 1.21143565e+08
+#  1.13182922e+08 1.03400960e+08 9.21621299e+07 7.94862092e+07
+#  6.52547663e+07 4.95601508e+07 3.29864561e+07 1.72764582e+07
+#  5.34804478e+06 8.04720005e+01]
+
 """
 Observations from numerical experiments:
 >> CG1 space for thickness would lead to strange stress concentration at

examples/t_junction/t_junction.py

@@ -28,6 +28,9 @@
 filename = "../../mesh/mesh-examples/t-junction/"+t_beam[2]
 with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r") as xdmf:
         mesh = xdmf.read_mesh(name="Grid")
+# sometimes it should be `name="mesh"` to avoid the error
+nel = mesh.topology.index_map(mesh.topology.dim).size_local
+nn = mesh.topology.index_map(0).size_local
 
 
 E_val = 1e7
@@ -71,11 +74,32 @@ def boundary(x):
 
 ########### Apply the point load #############################
 delta = Delta(x0=np.array([20.0, 0.0, 0.0]), f_p=(0.,0.,f_val))
-f1 = Function(W)
-f1_0,_ = f1.split()
-f1_0.interpolate(delta.eval)
 
+V1 = VectorFunctionSpace(mesh,("CG",1))
+f_1 = Function(V1)
+
+dofs = locate_dofs_geometrical(V1,lambda x: np.isclose(x.T,[20.0, 0.0, 0.0]).all(axis=1))
 
+f_1.x.array[dofs] = f_val
+print(f_1.vector.getArray()[dofs])
+# print(f_array)
+# delta = Delta_mp_1(x0=np.array([[20.0, 0.0, 0.0],[10.0, 0.0, 0.0]]), f_p=np.array([[0.,0.,f_val],[0.,0.,f_val]]))
+f1 = Function(W)
+f1_0,_ = f1.split()
+print("call evaluation...")
+# f1_0.interpolate(delta.eval)
+f1_0.interpolate(f_1)
+
+f1_x = computeNodalDisp(f1.sub(0))[0]
+f1_y = computeNodalDisp(f1.sub(0))[1]
+f1_z = computeNodalDisp(f1.sub(0))[2]
+print("-"*60)
+print("                               Projected CG2   "+"     Original     ")
+print("-"*60)
+print("Sum of forces in x-direction:", np.sum(f1_x))
+print("Sum of forces in y-direction:", np.sum(f1_y))
+print("Sum of forces in z-direction:", np.sum(f1_z))
+print("-"*60)
 # Assemble linear system
 a = derivative(F,w)
 L = -F

shell_analysis_fenicsx/shell_model.py

@@ -242,7 +242,14 @@ def penaltyResidual(self,u,v,g,dss,dSS):
                 (beta/h_E*inner(u-g, v))("+")*dSS + \
                 (beta/h_E*inner(u-g, v))("-")*dSS
 
-
+    def inertialResidual(self, rho, h):
+        h_mesh = CellDiameter(self.mesh)
+        retval = 0
+        retval += rho*h*inner(self.u_mid, self.du_mid)*dx
+        retval += rho*h*h_mesh**2*inner(self.theta, self.dtheta)*dx
+        drilling_inertia = Constant(self.mesh,1.0)*h**3*dx
+        # retval += drilling_inertia
+        return retval
 
 class ShellStressRM:
     """

shell_analysis_fenicsx/utils.py

@@ -106,6 +106,63 @@ def eval(self, x):
                 values[0][i] = self.f_p[0]
                 values[1][i] = self.f_p[1]
                 values[2][i] = self.f_p[2]
+                print(i)
+        print(np.sum(values,axis=1))
+        return values
+
+class Delta_cpt:
+    """
+    Delta function on closest points
+    """
+    def __init__(self, x0, f_p, **kwargs):
+        self.x0 = x0
+        self.f_p = f_p
+
+    def eval(self, x):
+        dist = []
+        values = np.zeros((3, x.shape[1]))
+        closest_local = []
+        for i in range(x.shape[1]):
+            x_i = np.array([x[0][i], x[1][i], x[2][i]])
+            dist_i = np.linalg.norm(x_i-self.x0)
+            dist.append(dist_i)
+
+        closest_local = np.argsort(np.array(dist))[:4]
+        print(closest_local)
+        print("applying forces to the closest point...")
+        values[0][closest_local] = self.f_p[0]
+        values[1][closest_local] = self.f_p[1]
+        values[2][closest_local] = self.f_p[2]
+        return values
+
+class Delta_mpt:
+    """
+    Multi-point delta function applied on the closest points
+    """
+    def __init__(self, x0, f_p, **kwargs):
+        self.x0 = x0
+        self.f_p = f_p
+
+    def eval(self, x):
+        values = np.zeros((3, x.shape[1]))
+        for j in range(self.x0.shape[0]):
+            x0_j = self.x0[:][j]
+            f_p_j = self.f_p[:][j]
+            print(x0_j, f_p_j)
+
+            dist = []
+            closest_local = []
+            for i in range(x.shape[1]):
+                x_i = np.array([x[0][i], x[1][i], x[2][i]])
+                dist_i = np.linalg.norm(x_i-x0_j)
+                dist.append(dist_i)
+
+            closest_local = np.argsort(np.array(dist))[:4]
+            print(closest_local)
+            print("applying forces to the closest point...")
+            values[0][closest_local] = f_p_j[0]
+            values[1][closest_local] = f_p_j[1]
+            values[2][closest_local] = f_p_j[2]
         return values