diff --git a/amigo/expressions.py b/amigo/expressions.py index e2e5a98..75f74a7 100644 --- a/amigo/expressions.py +++ b/amigo/expressions.py @@ -455,15 +455,16 @@ def __init__(self, consts=[], data=[], inputs=[], vars=[], rhs=[], lhs=[]): self.counter = 0 self.key_to_id = {} self.nodes = {} + self.expr_to_id = {} # Serialize the constants, data, inputs and vars - self.const_list = list(self._serialize_expr(c) for c in self.consts) - self.data_list = list(self._serialize_expr(d) for d in self.data) - self.input_list = list(self._serialize_expr(x) for x in self.inputs) - self.var_list = list(self._serialize_expr(v) for v in self.vars) + self.const_list = [self._serialize_expr(c) for c in self.consts] + self.data_list = [self._serialize_expr(d) for d in self.data] + self.input_list = [self._serialize_expr(x) for x in self.inputs] + self.var_list = [self._serialize_expr(v) for v in self.vars] # Serialize the expressions - self.rhs_list = list(self._serialize_expr(e) for e in self.rhs) + self.rhs_list = [self._serialize_expr(e) for e in self.rhs] self.counts = [0] * len(self.nodes) self.node_exprs = [None] * len(self.nodes) @@ -524,12 +525,14 @@ def __init__(self, consts=[], data=[], inputs=[], vars=[], rhs=[], lhs=[]): return def _serialize_expr(self, e: Expr): - key = e.serialize() - if key in self.key_to_id: - return self.key_to_id[key] + # Fast path: same Expr object already seen + eid = id(e) + if eid in self.expr_to_id: + return self.expr_to_id[eid] node = e.node if isinstance(node, ConstNode): + key = ("const", node.name, node.value, node.type) info = { "type": "const", "name": node.name, @@ -537,6 +540,7 @@ def _serialize_expr(self, e: Expr): "vartype": _type_to_str[node.type], } elif isinstance(node, VarNode): + key = ("var", node.name, node.shape, node.type, node.active) info = { "type": "var", "name": node.name, @@ -545,36 +549,53 @@ def _serialize_expr(self, e: Expr): "active": node.active, } elif isinstance(node, IndexNode): + expr_id = self._serialize_expr(node.expr) + key = ("index", expr_id, node.index) info = { "type": "index", - "expr": self._serialize_expr(node.expr), + "expr": expr_id, "index": node.index, } elif isinstance(node, PassiveNode): + expr_id = self._serialize_expr(node.expr) + key = ("passive", expr_id) info = { "type": "passive", - "expr": self._serialize_expr(node.expr), + "expr": expr_id, } elif isinstance(node, UnaryNode): + expr_id = self._serialize_expr(node.expr) + key = ("unary", node.op, expr_id) info = { "type": "unary", "op": node.op, - "expr": self._serialize_expr(node.expr), + "expr": expr_id, } elif isinstance(node, BinaryNode): + left_id = self._serialize_expr(node.left) + right_id = self._serialize_expr(node.right) + key = ("binary", node.op, left_id, right_id) info = { "type": "binary", "op": node.op, - "left": self._serialize_expr(node.left), - "right": self._serialize_expr(node.right), + "left": left_id, + "right": right_id, } else: raise TypeError(type(node)) + # Structural CSE: catches distinct Expr objects with identical structure + if key in self.key_to_id: + count = self.key_to_id[key] + self.expr_to_id[eid] = count + return count + count = self.counter + self.counter += 1 + self.key_to_id[key] = count + self.expr_to_id[eid] = count self.nodes[count] = info - self.counter += 1 return count diff --git a/amigo/fem/fem.py b/amigo/fem/fem.py index 6598493..ab53b0f 100644 --- a/amigo/fem/fem.py +++ b/amigo/fem/fem.py @@ -593,8 +593,10 @@ def create_model(self, module_name: str): model.link(f"{comp_name}.{name}", f"outputs.{name}[0]") # Set the node locations directly + spatial_names = ["x", "y", "z"][: self.mesh.X.shape[1]] for k, name in enumerate(self.geo_space.get_names("H1")): - model.set_data(f"geo.{name}", self.mesh.X[:, k]) + if name in spatial_names: + model.set_data(f"geo.{name}", self.mesh.X[:, k]) # Link the output to the finite element class return model diff --git a/amigo/model.py b/amigo/model.py index 2f806fc..89995a6 100644 --- a/amigo/model.py +++ b/amigo/model.py @@ -1755,7 +1755,7 @@ def _create_opt_problem(self, comm=COMM_WORLD): else: fixed_vars.append(self.get_indices(expr)[indices]) if fixed_vars: - fixed_vars = np.concatenate(fixed_vars) + fixed_vars = np.concatenate([np.atleast_1d(x) for x in fixed_vars]) vtypes[fixed_vars] = OptVarType.FIXED # Create the optimization problem diff --git a/examples/fem/mitc_cylinder/cylinder.py b/examples/fem/mitc_cylinder/cylinder.py new file mode 100644 index 0000000..d19c1fc --- /dev/null +++ b/examples/fem/mitc_cylinder/cylinder.py @@ -0,0 +1,666 @@ +import argparse +import numpy as np +import amigo as am +from amigo.fem import MITCTyingStrain, MITCElement, SolutionSpace, Mesh, Problem +from amigo.fem.basis import QuadLagrangeBasis, LagrangeBasis2D +from scipy.sparse.linalg import spsolve +from utils import write_vtu, get_exact_solution +import time + +from amigo.fem.basis import dot_product, eval_2d_monomials, eval_2d_monomial_grad + + +class NaturalShellGeoBasis(QuadLagrangeBasis): + """Shell Geometry Basis --> Computes Jacobian""" + + def __init__(self, degree, names, kind="data"): + super().__init__(degree, names, kind=kind) + + def compute_transform(self, geo): + x1, x2 = geo["x"]["grad"] + y1, y2 = geo["y"]["grad"] + z1, z2 = geo["z"]["grad"] + + nx0 = geo["nx"]["value"] + ny0 = geo["ny"]["value"] + nz0 = geo["nz"]["value"] + + nx1, nx2 = geo["nx"]["grad"] + ny1, ny2 = geo["ny"]["grad"] + nz1, nz2 = geo["nz"]["grad"] + + # Compute the normalized direction from the interpolation of the + # nodal normal vectors. This is used to compute the transformation + n0_norm = am.sqrt(nx0**2 + ny0**2 + nz0**2) + n0_inv = 1.0 / n0_norm + nx = n0_inv * nx0 + ny = n0_inv * ny0 + nz = n0_inv * nz0 + + # Project x,1 into the tangent plane of n + x1_dot_n = x1 * nx + y1 * ny + z1 * nz + + t1x0 = x1 - x1_dot_n * nx + t1y0 = y1 - x1_dot_n * ny + t1z0 = z1 - x1_dot_n * nz + t1_inv = 1.0 / am.sqrt(t1x0**2 + t1y0**2 + t1z0**2) + + t1x = t1_inv * t1x0 + t1y = t1_inv * t1y0 + t1z = t1_inv * t1z0 + + # Find the second in-plane direction t2 = n x t1 + t2x = ny * t1z - nz * t1y + t2y = nz * t1x - nx * t1z + t2z = nx * t1y - ny * t1x + + # Form the transformation matrix + T = [ + [t1x, t2x, nx], + [t1y, t2y, ny], + [t1z, t2z, nz], + ] + + # Compute the matrix: Jinv = X,xi^{-1} * T + # [x,1 x,2 n0x]^{-1} [t1x t2x nx] [ * * 0 ] + # [y,1 y,2 n0y] [t1y t2y ny] = [ * * 0 ] + # [z,1 z,2 n0z] [t1z t2z nz] [ 0 0 n0_inv ] + # Note that since n0 is parallel with n, the resulting matrix is + + # Covariant metric tensor + g11 = x1 * x1 + y1 * y1 + z1 * z1 + g12 = x1 * x2 + y1 * y2 + z1 * z2 + g22 = x2 * x2 + y2 * y2 + z2 * z2 + detg = g11 * g22 - g12 * g12 + detg_inv = 1.0 / detg + + # Determinant of the Jacobian transformation + detJ = am.sqrt(detg) + + # Contravariant tensor components + g11c = detg_inv * g22 + g12c = -detg_inv * g12 + g22c = detg_inv * g11 + + d11 = x1 * t1x + y1 * t1y + z1 * t1z + d12 = x1 * t2x + y1 * t2y + z1 * t2z + d21 = x2 * t1x + y2 * t1y + z2 * t1z + d22 = x2 * t2x + y2 * t2y + z2 * t2z + + # Components of the in-plane tensor + Jinv = [ + [g11c * d11 + g12c * d21, g11c * d12 + g12c * d22], + [g12c * d11 + g22c * d21, g12c * d12 + g22c * d22], + ] + + # Compute zJinv = X,xi^{-1} * n0,xi * X,xi^{-1} * T + def apply_Xinv_to_vec(qx, qy, qz): + a1_dot_q = x1 * qx + y1 * qy + z1 * qz + a2_dot_q = x2 * qx + y2 * qy + z2 * qz + + r1 = g11c * a1_dot_q + g12c * a2_dot_q + r2 = g12c * a1_dot_q + g22c * a2_dot_q + + return r1, r2 + + q1x = nx1 * Jinv[0][0] + nx2 * Jinv[1][0] + q1y = ny1 * Jinv[0][0] + ny2 * Jinv[1][0] + q1z = nz1 * Jinv[0][0] + nz2 * Jinv[1][0] + + zJ11, zJ21 = apply_Xinv_to_vec(q1x, q1y, q1z) + + q2x = nx1 * Jinv[0][1] + nx2 * Jinv[1][1] + q2y = ny1 * Jinv[0][1] + ny2 * Jinv[1][1] + q2z = nz1 * Jinv[0][1] + nz2 * Jinv[1][1] + + zJ12, zJ22 = apply_Xinv_to_vec(q2x, q2y, q2z) + + zJ31 = n0_inv * n0_inv * (nx0 * q1x + ny0 * q1y + nz0 * q1z) + zJ32 = n0_inv * n0_inv * (nx0 * q2x + ny0 * q2y + nz0 * q2z) + + # Set the rate terms + zJinv = [ + [zJ11, zJ12], + [zJ21, zJ22], + [zJ31, zJ32], + ] + + # To transform derivatives into the local coordinates + Jdict = {"Jinv": Jinv, "zJinv": zJinv, "J33": n0_inv, "T": T} + + return detJ, Jdict + + +class ShellSolnBasis(QuadLagrangeBasis): + def __init__(self, degree, kind="input"): + names = ["u", "v", "w", "rx", "ry", "rz"] + super().__init__(degree, names, kind=kind) + + def eval(self, comp, pt): + xi = pt[0] + eta = pt[1] + + # Evaluate the monomials + m = eval_2d_monomials(self.p, xi, eta, self.exps) + N = m @ self.C + + # Evaluate the derivatives of the monomials + mgrad = eval_2d_monomial_grad(self.p, xi, eta, self.exps) + Nxi = mgrad[:, 0] @ self.C + Neta = mgrad[:, 1] @ self.C + + # Get the rotations at the nodes + ry = comp.inputs["ry"] + rx = comp.inputs["rx"] + rz = comp.inputs["rz"] + + # Get the normals at the nodes + nx = comp.data["nx"] + ny = comp.data["ny"] + nz = comp.data["nz"] + + # Set up dx, dy, dz to interpolate the directors from the nodes + dx = [None] * self.nnodes + dy = [None] * self.nnodes + dz = [None] * self.nnodes + for i in range(self.nnodes): + dx[i] = ry[i] * nz[i] - rz[i] * ny[i] + dy[i] = rz[i] * nx[i] - rx[i] * nz[i] + dz[i] = rx[i] * ny[i] - ry[i] * nx[i] + + soln = {} + for name in self.names: + if self.kind == "input": + u = comp.inputs[name] + elif self.kind == "data": + u = comp.data[name] + elif self.kind == "multiplier": + u = comp.constraints.get_multipliers()[f"res_{name}"] + + soln[name] = { + "value": dot_product(u, N, n=self.nnodes), + "grad": [ + dot_product(u, Nxi, n=self.nnodes), + dot_product(u, Neta, n=self.nnodes), + ], + } + + for name, u in zip(["dx", "dy", "dz"], [dx, dy, dz]): + soln[name] = { + "value": dot_product(u, N, n=self.nnodes), + "grad": [ + dot_product(u, Nxi, n=self.nnodes), + dot_product(u, Neta, n=self.nnodes), + ], + } + + return soln + + def transform(self, detJ, Jdict, orig): + # Get the transformation information + Jinv = Jdict["Jinv"] + zJinv = Jdict["zJinv"] + T = Jdict["T"] + + # Compute the the values + soln = {} + + # The interpolated values from the solution field + soln["u"] = {"value": orig["u"]["value"]} + soln["v"] = {"value": orig["v"]["value"]} + soln["w"] = {"value": orig["w"]["value"]} + + soln["rx"] = {"value": orig["rx"]["value"]} + soln["ry"] = {"value": orig["ry"]["value"]} + soln["rz"] = {"value": orig["rz"]["value"]} + + u1, u2 = orig["u"]["grad"] + v1, v2 = orig["v"]["grad"] + w1, w2 = orig["w"]["grad"] + + dx = orig["dx"]["value"] + dy = orig["dy"]["value"] + dz = orig["dz"]["value"] + + dx1, dx2 = orig["dx"]["grad"] + dy1, dy2 = orig["dy"]["grad"] + dz1, dz2 = orig["dz"]["grad"] + + # Compute c = T^{T} * [d,1 d,2 0] + c11 = T[0][0] * dx1 + T[1][0] * dy1 + T[2][0] * dz1 + c21 = T[0][1] * dx1 + T[1][1] * dy1 + T[2][1] * dz1 + + c12 = T[0][0] * dx2 + T[1][0] * dy2 + T[2][0] * dz2 + c22 = T[0][1] * dx2 + T[1][1] * dy2 + T[2][1] * dz2 + + c31 = T[0][2] * dx1 + T[1][2] * dy1 + T[2][2] * dz1 + c32 = T[0][2] * dx2 + T[1][2] * dy2 + T[2][2] * dz2 + + # Compute b = T^{T} * [u,1 u,2 d] + b11 = T[0][0] * u1 + T[1][0] * v1 + T[2][0] * w1 + b12 = T[0][0] * u2 + T[1][0] * v2 + T[2][0] * w2 + b13 = T[0][0] * dx + T[1][0] * dy + T[2][0] * dz + + b21 = T[0][1] * u1 + T[1][1] * v1 + T[2][1] * w1 + b22 = T[0][1] * u2 + T[1][1] * v2 + T[2][1] * w2 + b23 = T[0][1] * dx + T[1][1] * dy + T[2][1] * dz + + b31 = T[0][2] * u1 + T[1][2] * v1 + T[2][2] * w1 + b32 = T[0][2] * u2 + T[1][2] * v2 + T[2][2] * w2 + b33 = T[0][2] * dx + T[1][2] * dy + T[2][2] * dz + + def row_transform(c1, c2, b1, b2, b3): + r1 = ( + c1 * Jinv[0][0] + + c2 * Jinv[1][0] + - b1 * zJinv[0][0] + - b2 * zJinv[1][0] + - b3 * zJinv[2][0] + ) + + r2 = ( + c1 * Jinv[0][1] + + c2 * Jinv[1][1] + - b1 * zJinv[0][1] + - b2 * zJinv[1][1] + - b3 * zJinv[2][1] + ) + + return [r1, r2] + + # Compute c * Jinv - b * zJinv + soln["u1"] = {"grad": row_transform(c11, c12, b11, b12, b13)} + soln["v1"] = {"grad": row_transform(c21, c22, b21, b22, b23)} + soln["w1"] = {"grad": row_transform(c31, c32, b31, b32, b33)} + + return soln + + +class MITC4ShellTying(MITCTyingStrain): + """ + MITC4 tying strains for a 3D curved shell + DOFS: u,v,w, rx, ry, rz + + """ + + def __init__(self): + super().__init__() + + def get_tying_points(self): + pts = [] + for index in range(9): + pts.append(self._get_tying_point(index)) + return pts + + def _get_tying_field_and_offset(self, idx): + if idx == 0 or idx == 1: + return "g11", 0 + if idx == 2 or idx == 3: + return "g22", 2 + if idx == 4: + return "g12", 4 + if idx == 5 or idx == 6: + return "g23", 5 + if idx == 7 or idx == 8: + return "g13", 7 + + def _get_tying_point(self, idx): + field, offset = self._get_tying_field_and_offset(idx) + + if field == "g11" or field == "g13": + return [(0, -1), (0, 1)][idx - offset] + if field == "g12": + return (0, 0) + if field == "g22" or field == "g23": + return [(-1, 0), (1, 0)][idx - offset] + + def eval_tying_strain(self, idx, geo, soln): + nx0 = geo["nx"]["value"] + ny0 = geo["ny"]["value"] + nz0 = geo["nz"]["value"] + + x1, x2 = geo["x"]["grad"] + y1, y2 = geo["y"]["grad"] + z1, z2 = geo["z"]["grad"] + + dx = soln["dx"]["value"] + dy = soln["dy"]["value"] + dz = soln["dz"]["value"] + + u1, u2 = soln["u"]["grad"] + v1, v2 = soln["v"]["grad"] + w1, w2 = soln["w"]["grad"] + + field, _ = self._get_tying_field_and_offset(idx) + + if field == "g11": + # G11 component + g11 = u1 * x1 + v1 * y1 + w1 * z1 + return g11 + elif field == "g22": + # G22 component + g22 = u2 * x2 + v2 * y2 + w2 * z2 + return g22 + elif field == "g12": + # G12 component + g12 = 0.5 * ((u1 * x2 + v1 * y2 + w1 * z2) + (u2 * x1 + v2 * y1 + w2 * z1)) + return g12 + elif field == "g23": + # G23 component + g23 = 0.5 * ( + (x2 * dx + y2 * dy + z2 * dz) + (nx0 * u2 + ny0 * v2 + nz0 * w2) + ) + return g23 + elif field == "g13": + # G13 component + g13 = 0.5 * ( + (x1 * dx + y1 * dy + z1 * dz) + (nx0 * u1 + ny0 * v1 + nz0 * w1) + ) + return g13 + + def interp_and_transform(self, pt, Jdict, e): + # Interpolate the tensorial components of the tying strains + g11 = 0.5 * ((1.0 - pt[1]) * e[0] + (1.0 + pt[1]) * e[1]) + g22 = 0.5 * ((1.0 - pt[0]) * e[2] + (1.0 + pt[0]) * e[3]) + g12 = e[4] + g23 = 0.5 * ((1.0 - pt[0]) * e[5] + (1.0 + pt[0]) * e[6]) + g13 = 0.5 * ((1.0 - pt[1]) * e[7] + (1.0 + pt[1]) * e[8]) + + # Extract the transformation matrix + Jinv = Jdict["Jinv"] + J33 = Jdict.get("J33", 1.0) + + J11 = Jinv[0][0] + J12 = Jinv[0][1] + J21 = Jinv[1][0] + J22 = Jinv[1][1] + + # First compute H = G * Jinv + h11 = g11 * J11 + g12 * J21 + h12 = g11 * J12 + g12 * J22 + + h21 = g12 * J11 + g22 * J21 + h22 = g12 * J12 + g22 * J22 + + # Then compute Gbar = Jinv^T * H + ex = J11 * h11 + J21 * h21 + ey = J12 * h12 + J22 * h22 + + # Transform and convert to engineering shear strain + gxy = 2.0 * (J11 * h12 + J21 * h22) + gxz = 2.0 * J33 * (J11 * g13 + J21 * g23) + gyz = 2.0 * J33 * (J12 * g13 + J22 * g23) + + out = {} + out["ex"] = {"value": ex} + out["ey"] = {"value": ey} + out["gxy"] = {"value": gxy} + + out["gxz"] = {"value": gxz} + out["gyz"] = {"value": gyz} + + return out + + +def integrand(soln, data=None, geo=None): + x = geo["x"]["value"] + y = geo["y"]["value"] + z = geo["z"]["value"] + + nx0 = geo["nx"]["value"] + ny0 = geo["ny"]["value"] + nz0 = geo["nz"]["value"] + + rx = soln["rx"]["value"] + ry = soln["ry"]["value"] + rz = soln["rz"]["value"] + + u = soln["u"]["value"] + v = soln["v"]["value"] + w = soln["w"]["value"] + + # Gradients for the bending terms + u1x, u1y = soln["u1"]["grad"] + v1x, v1y = soln["v1"]["grad"] + + # In-plane strains from MITC interpolation + ex = soln["ex"]["value"] + ey = soln["ey"]["value"] + gxy = soln["gxy"]["value"] + + # Shear strains from MITC interpolation + gxz = soln["gxz"]["value"] + gyz = soln["gyz"]["value"] + + kx = u1x + ky = v1y + kxy = u1y + v1x + + E = 70e9 + nu = 0.3 + ks = 5.0 / 6.0 + + # Set up the cylinder parameters + t = 0.01 + R = 1.0 + L = 2.0 + + M = 4 + N = 3 + + alpha = 4.0 / R + beta = 3 * np.pi / L + + theta = am.atan2(x, y) + + # Make the pressure dependent on the coordinate + pressure = am.sin(alpha * R * theta) * am.sin(beta * z) + + # Compute the in-plane energy + A = E * t / (1.0 - nu**2) + Um = 0.5 * A * ((ex**2 + ey**2 + 2.0 * nu * ex * ey) + 0.5 * (1.0 - nu) * gxy**2) + + # Compute the bending energy + D = E * t**3 / (12 * (1.0 - nu**2)) + Ub = 0.5 * D * ((kx**2 + ky**2 + 2.0 * nu * kx * ky) + 0.5 * (1.0 - nu) * kxy**2) + + # Compute the shear energy + G = 0.5 * E / (1.0 + nu) + Us = 0.5 * ks * G * t * (gxz**2 + gyz**2) + + # Compute the drill penalty + rot = rx * nx0 + ry * ny0 + rz * nz0 + + k_drill = 1e-4 * E * t + Ud = 0.5 * k_drill * rot**2 + + U = Um + Ub + Us + Ud + + # Set the work term + W = pressure * (u * nx0 + v * ny0 + w * nz0) + + return U - W + + +parser = argparse.ArgumentParser() +parser.add_argument("--build", action="store_true", default=False) +parser.add_argument( + "--solver", + dest="solver", + choices=["cholesky", "cholesky_left", "ldl", "scipy", "cuda"], + default="cholesky", +) +args = parser.parse_args() + +# Load the mesh +mesh = Mesh("cylinder.inp") +domains = mesh.get_domains() + +lateral_surfaces = ["SURFACE1"] +bottom_line = "LINE3" +lateral_line = "LINE2" +top_line = "LINE1" +print(f"Lateral: {lateral_surfaces}, bottom: {bottom_line}, top: {top_line}") + +# 6 DOF/node: u, v, w translations + rx, ry, rz global rotations +soln_space = SolutionSpace( + {"u": "H1", "v": "H1", "w": "H1", "rx": "H1", "ry": "H1", "rz": "H1"} +) +geo_space = SolutionSpace( + {"x": "H1", "y": "H1", "z": "H1", "nx": "H1", "ny": "H1", "nz": "H1"} +) +data_space = SolutionSpace({}) + +etype = "CPS4" + +degree = 1 +soln_basis = ShellSolnBasis(degree, kind="input") +geo_basis = NaturalShellGeoBasis(degree, ["x", "y", "z", "nx", "ny", "nz"], kind="data") +quadrature = mesh.get_quadrature(etype) +data_basis = mesh.get_basis(data_space, etype, kind="data") +mitc = MITC4ShellTying() + +shell_elem = MITCElement( + "Shell", soln_basis, data_basis, geo_basis, quadrature, mitc, integrand +) + +integrand_map = { + "shell": { + "target": lateral_surfaces, + "integrand": integrand, + }, +} +bc_map = { + "pinned_bottom": { + "type": "dirichlet", + "input": ["u", "v", "rz"], + "target": [bottom_line], + }, + "pinned_top": { + "type": "dirichlet", + "input": ["u", "v", "rz"], + "target": [top_line], + }, +} + +problem = Problem( + mesh, + soln_space, + data_space, + geo_space, + integrand_map=integrand_map, + bc_map=bc_map, + element_objs={("shell", etype): shell_elem}, +) + +model = problem.create_model("cylinder_shell") + +# Add a fixed boundary condition for the zeroth node. +model.add_fixed("soln.w[0]") + +if args.build: + model.build_module() +model.initialize() + +R = 1.0 +data = model.get_data_vector() +data["geo.nx"] = data["geo.x"] / R +data["geo.ny"] = data["geo.y"] / R + +# Create the vectors and matrices for the model +x = model.create_vector() +g = model.create_vector() +mat = model.create_matrix() + +# Copy the data over to the GPU +data = model.get_data_vector() +data.copy_host_to_device() + +print("Evaluating the Hessian...") +model.eval_gradient(x, g) +model.eval_hessian(x, mat) + +num_factors = 1 +if args.solver == "cuda": + from amigo.amigo import CSRMatFactorCuda + + # Duplicate the matrix + mat_copy = mat.duplicate() + mat_copy.copy(mat) + + pivot_eps = 1e-12 + solver = CSRMatFactorCuda(mat_copy, pivot_eps) + solver.factor() + + start_time = time.perf_counter() + for i in range(num_factors): + mat_copy.copy(mat) + solver.factor() + end_time = time.perf_counter() + tfactor = (end_time - start_time) / num_factors + + solver.solve(g.get_vector(), x.get_vector()) + + x.copy_device_to_host() +else: + g.copy_device_to_host() + mat.copy_data_device_to_host() + + if args.solver == "cholesky" or args.solver == "ldl": + stype = am.SolverType.CHOLESKY + if args.solver == "ldl": + stype = am.SolverType.LDL + + ldl = am.SparseLDL(mat, stype, ustab=0.05) + flag = ldl.factor() + + start_time = time.perf_counter() + for i in range(num_factors): + flag = ldl.factor() + end_time = time.perf_counter() + if flag != 0: + print(f"LDL factor flag {flag}") + + x[:] = g[:] + ldl.solve(x.get_vector()) + if stype == am.SolverType.LDL: + print("Inertia: ", ldl.get_inertia()) + + tfactor = (end_time - start_time) / num_factors + elif args.solver == "cholesky_left": + chol = am.SparseCholesky(mat) + start_time = time.perf_counter() + for i in range(num_factors): + flag = chol.factor() + end_time = time.perf_counter() + if flag != 0: + print(f"Cholesky factor flag {flag}") + + x[:] = g[:] + chol.solve(x.get_vector()) + + tfactor = (end_time - start_time) / num_factors + elif args.solver == "scipy": + csr = am.tocsr(mat) + + # This isn't a completely fair comparison + start_time = time.perf_counter() + x[:] = spsolve(csr, g[:]) + end_time = time.perf_counter() + + tfactor = end_time - start_time + +print(f"Factor time... {tfactor:.6f} seconds") + +u = x["soln.u"] +v = x["soln.v"] +w = x["soln.w"] + +u_ex, v_ex, w_ex = get_exact_solution(data["geo.x"], data["geo.y"], data["geo.z"]) + +w_diff = w - w_ex +print("Error = ", np.max(np.absolute(w_diff)) / np.max(w)) + +conn = np.vstack([mesh.get_conn(s, "CPS4") for s in lateral_surfaces]) +write_vtu(mesh, conn, u, v, w, filename="cylinder_shell.vtu") +write_vtu(mesh, conn, u_ex, v_ex, w_ex, filename="cylinder_shell_exact.vtu") diff --git a/examples/fem/mitc_cylinder/cylinder_mesh.py b/examples/fem/mitc_cylinder/cylinder_mesh.py index b01767c..8fea0e2 100644 --- a/examples/fem/mitc_cylinder/cylinder_mesh.py +++ b/examples/fem/mitc_cylinder/cylinder_mesh.py @@ -1,25 +1,69 @@ import gmsh - -R = 0.5 # radius -L = 2.0 # length -lc = 0.05 # mesh size +import math gmsh.initialize() gmsh.model.add("cylinder") -gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, L, R) -gmsh.model.occ.synchronize() +# Parameters +R = 1.0 +H = 2.0 +n_circ = 512 # number of quads around circumference +n_axial = 256 # number of quads along height + +# Use the built-in geo kernel +geo = gmsh.model.geo + +# Create a rectangular parametric domain: +# s in [0, 2*pi*R], z in [0, H] +p1 = geo.addPoint(0.0, 0.0, 0.0) +p2 = geo.addPoint(2.0 * math.pi * R, 0.0, 0.0) +p3 = geo.addPoint(2.0 * math.pi * R, H, 0.0) +p4 = geo.addPoint(0.0, H, 0.0) + +l1 = geo.addLine(p1, p2) +l2 = geo.addLine(p2, p3) +l3 = geo.addLine(p3, p4) +l4 = geo.addLine(p4, p1) + +loop = geo.addCurveLoop([l1, l2, l3, l4]) +surf = geo.addPlaneSurface([loop]) + +# Structured mesh constraints +# Number of points = number of elements + 1 +geo.mesh.setTransfiniteCurve(l1, n_circ + 1) +geo.mesh.setTransfiniteCurve(l3, n_circ + 1) +geo.mesh.setTransfiniteCurve(l2, n_axial + 1) +geo.mesh.setTransfiniteCurve(l4, n_axial + 1) + +geo.mesh.setTransfiniteSurface(surf) +geo.mesh.setRecombine(2, surf) -gmsh.option.setNumber("Mesh.Algorithm", 8) -gmsh.option.setNumber("Mesh.RecombineAll", 1) -gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) -gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc) -gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc) +geo.synchronize() +# Generate the flat quad mesh gmsh.model.mesh.generate(2) -gmsh.model.mesh.recombine() +node_tags, coords, _ = gmsh.model.mesh.getNodes() + +for i, tag in enumerate(node_tags): + s = coords[3 * i + 0] + z = coords[3 * i + 1] + + theta = s / R + x = R * math.cos(theta) + y = R * math.sin(theta) + + gmsh.model.mesh.setNode(int(tag), [x, y, z], []) + +# Optional: remove duplicate nodes on the seam +# This merges the theta=0 and theta=2*pi node columns. +gmsh.model.mesh.removeDuplicateNodes() + +# Save gmsh.write("cylinder.inp") gmsh.write("cylinder.msh") + +# Launch GUI if desired +# gmsh.fltk.run() + gmsh.finalize() -print("Done") diff --git a/examples/fem/mitc_cylinder/cylinder_shell.py b/examples/fem/mitc_cylinder/cylinder_shell.py deleted file mode 100644 index bcef3c9..0000000 --- a/examples/fem/mitc_cylinder/cylinder_shell.py +++ /dev/null @@ -1,339 +0,0 @@ -""" -MITC4 shell FEM on a cylinder. -- 3D curved surface (x, y, z nodes from cylinder.inp) -- CPS4 bilinear quad elements -- 5 DOF/node: u, v, w (translations) + t1, t2 (rotations in local shell frame) -- Pinned bottom (z=0), unit axial pressure load on top ring -""" - -import argparse -import numpy as np -import amigo as am -from amigo.fem import MITCTyingStrain, MITCElement, SolutionSpace, Mesh, Problem -from amigo.fem.basis import QuadLagrangeBasis, LagrangeBasis2D -from scipy.sparse.linalg import spsolve -from utils import write_vtu - - -class ShellGeoBasis(QuadLagrangeBasis): - """Shell Geometry Basis --> Computes Jacobian""" - - def __init__(self, names, kind="data"): - super().__init__(1, names, kind=kind) - - def compute_transform(self, geo): - x1, x2 = geo["x"]["grad"] - y1, y2 = geo["y"]["grad"] - z1, z2 = geo["z"]["grad"] - - # Surface normal via cross product a1 x a2 - nx = y1 * z2 - z1 * y2 - ny = z1 * x2 - x1 * z2 - nz = x1 * y2 - y1 * x2 - detJ = am.sqrt(nx**2 + ny**2 + nz**2) - - # Covariant metric tensor - g11 = x1 * x1 + y1 * y1 + z1 * z1 - g12 = x1 * x2 + y1 * y2 + z1 * z2 - g22 = x2 * x2 + y2 * y2 + z2 * z2 - inv_det_g = 1.0 / (g11 * g22 - g12 * g12) - - # Contravariant metric - g11c = g22 * inv_det_g - g12c = -g12 * inv_det_g - g22c = g11 * inv_det_g - - # 2x3 left-pseudoinverse: Jinv = G^{-1} * J^T - Jinv = [ - [g11c * x1 + g12c * x2, g11c * y1 + g12c * y2, g11c * z1 + g12c * z2], - [g12c * x1 + g22c * x2, g12c * y1 + g22c * y2, g12c * z1 + g22c * z2], - ] - return detJ, Jinv - - def transform(self, detJ, Jinv, orig): - # Pass through parametric grads unchanged — integrand uses covariant components - return { - name: {"value": orig[name]["value"], "grad": orig[name]["grad"]} - for name in orig - } - - -class ShellSolnBasis(QuadLagrangeBasis): - def __init__(self, names, kind="input"): - super().__init__(1, names, kind=kind) - - def transform(self, detJ, Jinv, orig): - # Keep parametric grads; integrand and tying code work in covariant coords - return { - name: {"value": orig[name]["value"], "grad": orig[name]["grad"]} - for name in orig - } - - -class MITC4ShellTying(MITCTyingStrain): - """ - MITC4 tying strains for a 3D curved shell - DOFS: u,v,w, rx, ry, rz - - """ - - def get_tying_points(self): - return [(-1, 0), (1, 0), (0, -1), (0, 1)] - - def eval_tying_strain(self, idx, geo, soln): - x1, x2 = geo["x"]["grad"] - y1, y2 = geo["y"]["grad"] - z1, z2 = geo["z"]["grad"] - - # Unit normal - nx = y1 * z2 - z1 * y2 - ny = z1 * x2 - x1 * z2 - nz = x1 * y2 - y1 * x2 - nmag = am.sqrt(nx * nx + ny * ny + nz * nz) - nx = nx / nmag - ny = ny / nmag - nz = nz / nmag - - rx = soln["rx"]["value"] - ry = soln["ry"]["value"] - rz = soln["rz"]["value"] - u1, u2 = soln["u"]["grad"] - v1, v2 = soln["v"]["grad"] - w1, w2 = soln["w"]["grad"] - - # d = theta x n = (rx,ry,rz) x (nx,ny,nz) - dx = ry * nz - rz * ny - dy = rz * nx - rx * nz - dz = rx * ny - ry * nx - - if idx < 2: - # gamma_23: u,2 · n + d · a2 - return u2 * nx + v2 * ny + w2 * nz + dx * x2 + dy * y2 + dz * z2 - else: - # gamma_13: u,1 · n + d · a1 - return u1 * nx + v1 * ny + w1 * nz + dx * x1 + dy * y1 + dz * z1 - - def interp_and_transform(self, pt, Jinv, e): - xi, eta = pt - gamma_2 = 0.5 * ((1.0 - xi) * e[0] + (1.0 + xi) * e[1]) - gamma_1 = 0.5 * ((1.0 - eta) * e[2] + (1.0 + eta) * e[3]) - return {"gs1": {"value": gamma_1}, "gs2": {"value": gamma_2}} - - -def shell_integrand(soln, data=None, geo=None): - """ - Shell Potential Energy Integrand - Membrane energy + Bending energy + Shear energy + Drill Energy - """ - E = 70e9 - nu = 0.3 - t = 0.002 - ks = 5.0 / 6.0 - G = E / (2.0 * (1.0 + nu)) - - x1, x2 = geo["x"]["grad"] - y1, y2 = geo["y"]["grad"] - z1, z2 = geo["z"]["grad"] - - g11 = x1 * x1 + y1 * y1 + z1 * z1 - g12 = x1 * x2 + y1 * y2 + z1 * z2 - g22 = x2 * x2 + y2 * y2 + z2 * z2 - inv_det_g = 1.0 / (g11 * g22 - g12 * g12) - g11c = g22 * inv_det_g - g12c = -g12 * inv_det_g - g22c = g11 * inv_det_g - - u1, u2 = soln["u"]["grad"] - v1, v2 = soln["v"]["grad"] - w1, w2 = soln["w"]["grad"] - rx1, rx2 = soln["rx"]["grad"] - ry1, ry2 = soln["ry"]["grad"] - rz1, rz2 = soln["rz"]["grad"] - - # Unit normal - nx = y1 * z2 - z1 * y2 - ny = z1 * x2 - x1 * z2 - nz = x1 * y2 - y1 * x2 - nmag = am.sqrt(nx * nx + ny * ny + nz * nz) - nx = nx / nmag - ny = ny / nmag - nz = nz / nmag - - # Membrane strains: eps_ab = u_,a · a_b (symmetrized) - eps11 = u1 * x1 + v1 * y1 + w1 * z1 - eps22 = u2 * x2 + v2 * y2 + w2 * z2 - eps12 = 0.5 * (u1 * x2 + v1 * y2 + w1 * z2 + u2 * x1 + v2 * y1 + w2 * z1) - - tr_eps = g11c * eps11 + 2.0 * g12c * eps12 + g22c * eps22 - eps11c = g11c * g11c * eps11 + 2.0 * g11c * g12c * eps12 + g12c * g12c * eps22 - eps12c = ( - g11c * g12c * eps11 + (g11c * g22c + g12c * g12c) * eps12 + g12c * g22c * eps22 - ) - eps22c = g12c * g12c * eps11 + 2.0 * g12c * g22c * eps12 + g22c * g22c * eps22 - - fac_m = E * t / (1.0 - nu**2) - U_membrane = ( - 0.5 - * fac_m - * ( - nu * tr_eps * tr_eps - + (1.0 - nu) * (eps11c * eps11 + 2.0 * eps12c * eps12 + eps22c * eps22) - ) - ) - - # Bending strains: kappa_ab = 0.5*(d,a · a_b + d,b · a_a) - # d = theta x n - # d,1 = (rx1, ry1, rz1) x (nx, ny, nz) - d1x = ry1 * nz - rz1 * ny - d1y = rz1 * nx - rx1 * nz - d1z = rx1 * ny - ry1 * nx - # d,2 = (rx2, ry2, rz2) x (nx, ny, nz) - d2x = ry2 * nz - rz2 * ny - d2y = rz2 * nx - rx2 * nz - d2z = rx2 * ny - ry2 * nx - - kap11 = d1x * x1 + d1y * y1 + d1z * z1 - kap22 = d2x * x2 + d2y * y2 + d2z * z2 - kap12 = 0.5 * (d1x * x2 + d1y * y2 + d1z * z2 + d2x * x1 + d2y * y1 + d2z * z1) - - tr_kap = g11c * kap11 + 2.0 * g12c * kap12 + g22c * kap22 - kap11c = g11c * g11c * kap11 + 2.0 * g11c * g12c * kap12 + g12c * g12c * kap22 - kap12c = ( - g11c * g12c * kap11 + (g11c * g22c + g12c * g12c) * kap12 + g12c * g22c * kap22 - ) - kap22c = g12c * g12c * kap11 + 2.0 * g12c * g22c * kap12 + g22c * g22c * kap22 - - fac_b = E * t**3 / (12.0 * (1.0 - nu**2)) - U_bending = ( - 0.5 - * fac_b - * ( - nu * tr_kap * tr_kap - + (1.0 - nu) * (kap11c * kap11 + 2.0 * kap12c * kap12 + kap22c * kap22) - ) - ) - - # Transverse shear (MITC tying, covariant) - gs1 = soln["gs1"]["value"] - gs2 = soln["gs2"]["value"] - gs1c = g11c * gs1 + g12c * gs2 - gs2c = g12c * gs1 + g22c * gs2 - U_shear = 0.5 * ks * G * t * (gs1c * gs1 + gs2c * gs2) - - # Drilling penalty - # Drilling rotation = theta · n = rx * nx + ry * ny + rz * nz - rx = soln["rx"]["value"] - ry = soln["ry"]["value"] - rz = soln["rz"]["value"] - drill = rx * nx + ry * ny + rz * nz - - # Also penalize the antisymmetric in-plane displacement gradient - omega = 0.5 * (u1 * x2 + v1 * y2 + w1 * z2 - u2 * x1 - v2 * y1 - w2 * z1) - drill_alpha = 1.0e-3 * E * t / (1.0 - nu**2) - U_drill = 0.5 * drill_alpha * (drill * drill + omega * omega) - - # Distributed load q [N/m^2] in +z: W_ext = q * w - q = -1e6 - W_ext = q * soln["w"]["value"] - - return U_membrane + U_bending + U_shear + U_drill - W_ext - - -parser = argparse.ArgumentParser() -parser.add_argument("--build", action="store_true", default=False) -args = parser.parse_args() - -mesh = Mesh("cylinder.inp") -domains = mesh.get_domains() - -lateral_surfaces = ["SURFACE1"] -bottom_line = "LINE3" -lateral_line = "LINE2" -top_line = "LINE1" -print(f"Lateral: {lateral_surfaces}, bottom: {bottom_line}, top: {top_line}") - -# 6 DOF/node: u, v, w translations + rx, ry, rz global rotations -soln_space = SolutionSpace( - {"u": "H1", "v": "H1", "w": "H1", "rx": "H1", "ry": "H1", "rz": "H1"} -) -geo_space = SolutionSpace({"x": "H1", "y": "H1", "z": "H1"}) -data_space = SolutionSpace({}) - -etype = "CPS4" - -soln_basis = ShellSolnBasis(["u", "v", "w", "rx", "ry", "rz"], kind="input") -data_basis = mesh.get_basis(data_space, etype, kind="data") -geo_basis = ShellGeoBasis(["x", "y", "z"], kind="data") -quadrature = mesh.get_quadrature(etype) -mitc = MITC4ShellTying() - -shell_elem = MITCElement( - "Shell", soln_basis, data_basis, geo_basis, quadrature, mitc, shell_integrand -) - -integrand_map = { - "shell": { - "target": lateral_surfaces, - "integrand": shell_integrand, - }, -} -bc_map = { - "pinned_bottom": { - "type": "dirichlet", - "input": ["u", "v", "w", "rx", "ry", "rz"], - "target": [bottom_line], - }, -} - -problem = Problem( - mesh, - soln_space, - data_space, - geo_space, - integrand_map=integrand_map, - bc_map=bc_map, - element_objs={("shell", etype): shell_elem}, -) - -model = problem.create_model("cylinder_shell") - -if args.build: - model.build_module() - -model.initialize() - -x = model.create_vector() -g = model.create_vector() -mat = model.create_matrix() - -model.eval_gradient(x, g) -model.eval_hessian(x, mat) - -K = am.tocsr(mat) -# Penalize any zero-stiffness DOFs (cap nodes not in lateral surface) -diag = K.diagonal() -from scipy.sparse import diags - -K = K + diags(np.where(np.abs(diag) < 1e-10, 1.0, 0.0)) - -x.get_vector().get_array()[:] = spsolve(K, g.get_vector().get_array()) - -# ldl = am.SparseLDL(mat, am.SolverType.LDL, ustab=0.05) -# flag = ldl.factor() - -# x.copy(g) -# ldl.solve(x.get_vector()) - -u = x["soln.u"] -v = x["soln.v"] -w = x["soln.w"] - -disp_mag = np.sqrt(u**2 + v**2 + w**2) -print(f"Max displacement magnitude: {disp_mag.max():.6e}") -print(f"u: [{u.min():.3e}, {u.max():.3e}]") -print(f"v: [{v.min():.3e}, {v.max():.3e}]") -print(f"w: [{w.min():.3e}, {w.max():.3e}]") - -conn = np.vstack([mesh.get_conn(s, "CPS4") for s in lateral_surfaces]) -write_vtu(mesh, conn, u, v, w) -print("Wrote cylinder_shell.vtu") diff --git a/examples/fem/mitc_cylinder/cylinder_shell_endtraction.py b/examples/fem/mitc_cylinder/cylinder_shell_endtraction.py deleted file mode 100644 index 825a693..0000000 --- a/examples/fem/mitc_cylinder/cylinder_shell_endtraction.py +++ /dev/null @@ -1,318 +0,0 @@ -""" -MITC4 shell FEM on a cylinder — end traction variant. -Load applied as axial traction integrated along top ring (LINE1, T3D2), -instead of a distributed surface pressure. -Author: Jack Turbush -""" - -import argparse -import numpy as np -import amigo as am -from amigo.fem import MITCTyingStrain, MITCElement, SolutionSpace, Mesh, Problem -from amigo.fem.basis import QuadLagrangeBasis -from scipy.sparse.linalg import spsolve -from scipy.sparse import diags -from utils import write_vtu - - -class ShellGeoBasis(QuadLagrangeBasis): - def __init__(self, names, kind="data"): - super().__init__(1, names, kind=kind) - - def compute_transform(self, geo): - x1, x2 = geo["x"]["grad"] - y1, y2 = geo["y"]["grad"] - z1, z2 = geo["z"]["grad"] - - # Surface normal via cross product a1 x a2 - nx = y1 * z2 - z1 * y2 - ny = z1 * x2 - x1 * z2 - nz = x1 * y2 - y1 * x2 - detJ = am.sqrt(nx**2 + ny**2 + nz**2) - - # Covariant metric tensor - g11 = x1 * x1 + y1 * y1 + z1 * z1 - g12 = x1 * x2 + y1 * y2 + z1 * z2 - g22 = x2 * x2 + y2 * y2 + z2 * z2 - inv_det_g = 1.0 / (g11 * g22 - g12 * g12) - - # Contravariant metric - g11c = g22 * inv_det_g - g12c = -g12 * inv_det_g - g22c = g11 * inv_det_g - - # 2x3 left-pseudoinverse: Jinv = G^{-1} * J^T - Jinv = [ - [g11c * x1 + g12c * x2, g11c * y1 + g12c * y2, g11c * z1 + g12c * z2], - [g12c * x1 + g22c * x2, g12c * y1 + g22c * y2, g12c * z1 + g22c * z2], - ] - return detJ, Jinv - - def transform(self, detJ, Jinv, orig): - # Pass through parametric grads unchanged — integrand uses covariant components - return { - name: {"value": orig[name]["value"], "grad": orig[name]["grad"]} - for name in orig - } - - -class ShellSolnBasis(QuadLagrangeBasis): - def __init__(self, names, kind="input"): - super().__init__(1, names, kind=kind) - - def transform(self, detJ, Jinv, orig): - return { - name: {"value": orig[name]["value"], "grad": orig[name]["grad"]} - for name in orig - } - - -class MITC4ShellTying(MITCTyingStrain): - def get_tying_points(self): - return [(-1, 0), (1, 0), (0, -1), (0, 1)] - - def eval_tying_strain(self, idx, geo, soln): - x1, x2 = geo["x"]["grad"] - y1, y2 = geo["y"]["grad"] - z1, z2 = geo["z"]["grad"] - - # Unit normal - nx = y1 * z2 - z1 * y2 - ny = z1 * x2 - x1 * z2 - nz = x1 * y2 - y1 * x2 - nmag = am.sqrt(nx * nx + ny * ny + nz * nz) - nx, ny, nz = nx / nmag, ny / nmag, nz / nmag - - rx, ry, rz = soln["rx"]["value"], soln["ry"]["value"], soln["rz"]["value"] - u1, u2 = soln["u"]["grad"] - v1, v2 = soln["v"]["grad"] - w1, w2 = soln["w"]["grad"] - - # d = theta x n = (rx,ry,rz) x (nx,ny,nz) - dx = ry * nz - rz * ny - dy = rz * nx - rx * nz - dz = rx * ny - ry * nx - - if idx < 2: - # gamma_23 - return u2 * nx + v2 * ny + w2 * nz + dx * x2 + dy * y2 + dz * z2 - else: - # gamma_13 - return u1 * nx + v1 * ny + w1 * nz + dx * x1 + dy * y1 + dz * z1 - - def interp_and_transform(self, pt, Jinv, e): - xi, eta = pt - gamma_2 = 0.5 * ((1.0 - xi) * e[0] + (1.0 + xi) * e[1]) - gamma_1 = 0.5 * ((1.0 - eta) * e[2] + (1.0 + eta) * e[3]) - return {"gs1": {"value": gamma_1}, "gs2": {"value": gamma_2}} - - -def shell_integrand(soln, data=None, geo=None): - """ - Shell Potential Energy Integrand - Membrane energy + Bending energy + Shear energy + Drill Energy - """ - - E = 70e9 - nu = 0.3 - t = 0.002 - ks = 5.0 / 6.0 - G = E / (2.0 * (1.0 + nu)) - - x1, x2 = geo["x"]["grad"] - y1, y2 = geo["y"]["grad"] - z1, z2 = geo["z"]["grad"] - - g11 = x1 * x1 + y1 * y1 + z1 * z1 - g12 = x1 * x2 + y1 * y2 + z1 * z2 - g22 = x2 * x2 + y2 * y2 + z2 * z2 - inv_det_g = 1.0 / (g11 * g22 - g12 * g12) - g11c = g22 * inv_det_g - g12c = -g12 * inv_det_g - g22c = g11 * inv_det_g - - u1, u2 = soln["u"]["grad"] - v1, v2 = soln["v"]["grad"] - w1, w2 = soln["w"]["grad"] - rx1, rx2 = soln["rx"]["grad"] - ry1, ry2 = soln["ry"]["grad"] - rz1, rz2 = soln["rz"]["grad"] - - # Unit normal - nx = y1 * z2 - z1 * y2 - ny = z1 * x2 - x1 * z2 - nz = x1 * y2 - y1 * x2 - nmag = am.sqrt(nx * nx + ny * ny + nz * nz) - nx, ny, nz = nx / nmag, ny / nmag, nz / nmag - - eps11 = u1 * x1 + v1 * y1 + w1 * z1 - eps22 = u2 * x2 + v2 * y2 + w2 * z2 - eps12 = 0.5 * (u1 * x2 + v1 * y2 + w1 * z2 + u2 * x1 + v2 * y1 + w2 * z1) - - tr_eps = g11c * eps11 + 2.0 * g12c * eps12 + g22c * eps22 - eps11c = g11c * g11c * eps11 + 2.0 * g11c * g12c * eps12 + g12c * g12c * eps22 - eps12c = ( - g11c * g12c * eps11 + (g11c * g22c + g12c * g12c) * eps12 + g12c * g22c * eps22 - ) - eps22c = g12c * g12c * eps11 + 2.0 * g12c * g22c * eps12 + g22c * g22c * eps22 - - fac_m = E * t / (1.0 - nu**2) - U_membrane = ( - 0.5 - * fac_m - * ( - nu * tr_eps * tr_eps - + (1.0 - nu) * (eps11c * eps11 + 2.0 * eps12c * eps12 + eps22c * eps22) - ) - ) - - # Bending strains: kappa_ab = 0.5*(d,a · a_b + d,b · a_a) - # d = theta x n - # d,1 = (rx1, ry1, rz1) x (nx, ny, nz) - d1x = ry1 * nz - rz1 * ny - d1y = rz1 * nx - rx1 * nz - d1z = rx1 * ny - ry1 * nx - # d,2 = (rx2, ry2, rz2) x (nx, ny, nz) - d2x = ry2 * nz - rz2 * ny - d2y = rz2 * nx - rx2 * nz - d2z = rx2 * ny - ry2 * nx - - kap11 = d1x * x1 + d1y * y1 + d1z * z1 - kap22 = d2x * x2 + d2y * y2 + d2z * z2 - kap12 = 0.5 * (d1x * x2 + d1y * y2 + d1z * z2 + d2x * x1 + d2y * y1 + d2z * z1) - - tr_kap = g11c * kap11 + 2.0 * g12c * kap12 + g22c * kap22 - kap11c = g11c * g11c * kap11 + 2.0 * g11c * g12c * kap12 + g12c * g12c * kap22 - kap12c = ( - g11c * g12c * kap11 + (g11c * g22c + g12c * g12c) * kap12 + g12c * g22c * kap22 - ) - kap22c = g12c * g12c * kap11 + 2.0 * g12c * g22c * kap12 + g22c * g22c * kap22 - - fac_b = E * t**3 / (12.0 * (1.0 - nu**2)) - U_bending = ( - 0.5 - * fac_b - * ( - nu * tr_kap * tr_kap - + (1.0 - nu) * (kap11c * kap11 + 2.0 * kap12c * kap12 + kap22c * kap22) - ) - ) - # Transverse shear (MITC tying, covariant) - gs1, gs2 = soln["gs1"]["value"], soln["gs2"]["value"] - gs1c = g11c * gs1 + g12c * gs2 - gs2c = g12c * gs1 + g22c * gs2 - U_shear = 0.5 * ks * G * t * (gs1c * gs1 + gs2c * gs2) - - # Drilling penalty - # Drilling rotation = dot(theta, n) = rx*nx + ry*ny + rz*nz - rx, ry, rz = soln["rx"]["value"], soln["ry"]["value"], soln["rz"]["value"] - drill = rx * nx + ry * ny + rz * nz - - # Also penalize the antisymmetric in-plane displacement gradient - omega = 0.5 * (u1 * x2 + v1 * y2 + w1 * z2 - u2 * x1 - v2 * y1 - w2 * z1) - drill_alpha = 1.0e-3 * E * t / (1.0 - nu**2) - U_drill = 0.5 * drill_alpha * (drill * drill + omega * omega) - - return U_membrane + U_bending + U_shear + U_drill - - -def end_traction_integrand(soln, data=None, geo=None): - """Axial traction on top ring: W_ext = tz * w, integrated along LINE1 (T3D2).""" - tz = 1e6 # N/m — axial load per unit arc length (negative = compression) - return -tz * soln["w"]["value"] - - -parser = argparse.ArgumentParser() -parser.add_argument("--build", action="store_true", default=False) -args = parser.parse_args() - -mesh = Mesh("cylinder_coarse.inp") -domains = mesh.get_domains() - - -# groups made with cylinder_mesh.py -lateral_surfaces = ["SURFACE1"] -bottom_line = "Line3" -lateral_line = "Line2" -top_line = "Line1" -print(f"Lateral: {lateral_surfaces}, bottom: {bottom_line}, top: {top_line}") - -soln_space = SolutionSpace( - {"u": "H1", "v": "H1", "w": "H1", "rx": "H1", "ry": "H1", "rz": "H1"} -) -geo_space = SolutionSpace({"x": "H1", "y": "H1", "z": "H1"}) -data_space = SolutionSpace({}) - -etype = "CPS4" - -soln_basis = ShellSolnBasis(["u", "v", "w", "rx", "ry", "rz"], kind="input") -data_basis = mesh.get_basis(data_space, etype, kind="data") -geo_basis = ShellGeoBasis(["x", "y", "z"], kind="data") -quadrature = mesh.get_quadrature(etype) -mitc = MITC4ShellTying() - -shell_elem = MITCElement( - "Shell", soln_basis, data_basis, geo_basis, quadrature, mitc, shell_integrand -) - -integrand_map = { - "shell": { - "target": lateral_surfaces, - "integrand": shell_integrand, - }, - "end_traction": { - "target": [top_line], # LINE1: top ring, T3D2 elements - "integrand": end_traction_integrand, - }, -} -bc_map = { - "pinned_bottom": { - "type": "dirichlet", - "input": ["u", "v", "w", "rx", "ry", "rz"], - "target": [bottom_line], - }, -} - -problem = Problem( - mesh, - soln_space, - data_space, - geo_space, - integrand_map=integrand_map, - bc_map=bc_map, - element_objs={("shell", etype): shell_elem}, -) - -model = problem.create_model("cylinder_shell") - -if args.build: - model.build_module() - -model.initialize() - -x = model.create_vector() -g = model.create_vector() -mat = model.create_matrix() - -model.eval_gradient(x, g) -model.eval_hessian(x, mat) - -K = am.tocsr(mat) -diag = K.diagonal() -K = K + diags(np.where(np.abs(diag) < 1e-10, 1.0, 0.0)) - -x.get_vector().get_array()[:] = spsolve(K, g.get_vector().get_array()) - -u = x["soln.u"] -v = x["soln.v"] -w = x["soln.w"] - -disp_mag = np.sqrt(u**2 + v**2 + w**2) -print(f"Max displacement magnitude: {disp_mag.max():.6e}") -print(f"u: [{u.min():.3e}, {u.max():.3e}]") -print(f"v: [{v.min():.3e}, {v.max():.3e}]") -print(f"w: [{w.min():.3e}, {w.max():.3e}]") - -conn = np.vstack([mesh.get_conn(s, "CPS4") for s in lateral_surfaces]) -write_vtu(mesh, conn, u, v, w, filename="cylinder_shell_endtraction.vtu") -print("Wrote cylinder_shell_endtraction.vtu") diff --git a/examples/fem/mitc_cylinder/utils.py b/examples/fem/mitc_cylinder/utils.py index 12b0ffd..c9d1952 100644 --- a/examples/fem/mitc_cylinder/utils.py +++ b/examples/fem/mitc_cylinder/utils.py @@ -1,9 +1,110 @@ -# write VTU function for amigo fem results +import numpy as np -def write_vtu(mesh, conns, u, v, w, filename="cylinder_shell.vtu"): - import numpy as np +def compute_coefficients(E=70e9, nu=0.3, R=1.0, L=2.0, t=0.01, ks=5.0 / 6.0, M=4, N=3): + """ + Compute the scalar coefficients for a cylindrical mesh + + Returns + ------- + U, V, W, theta, phi : scalars + Solution coefficients. + """ + + alpha = 4.0 / R + beta = 3 * np.pi / L + ainv = 1.0 / R + + G = 0.5 * E / (1.0 + nu) + A11 = A22 = E * t / (1.0 - nu**2) + A12 = nu * A11 + A33 = G * t + D11 = D22 = E * t**3 / (12.0 * (1.0 - nu**2)) + D12 = nu * D11 + D33 = G * t**3 / 12.0 + bA11 = bA22 = ks * G * t + + rhs = np.zeros(5, dtype=float) + rhs[2] = 1.0 + + # Explicit matrix A, mainly for comparison/debugging + A = np.zeros((5, 5), dtype=float) + + # The first equation for U + A[0, 0] = ( + -(A11 * beta * beta + A33 * alpha * alpha) - D33 * ainv * ainv * alpha * alpha + ) + A[0, 1] = -(A33 + A12) * alpha * beta + A[0, 2] = A12 * beta * ainv + A[0, 3] = D33 * ainv * alpha * alpha + A[0, 4] = D33 * ainv * alpha * beta + + # The second equation for V + A[1, 0] = -(A12 + A33) * alpha * beta + A[1, 1] = ( + -(A33 * beta * beta + A22 * alpha * alpha) + - ainv * ainv * bA11 + - D22 * ainv * ainv * alpha * alpha + ) + A[1, 2] = (A22 + bA11) * ainv * alpha + D22 * alpha * ainv**3 + A[1, 3] = D12 * ainv * alpha * beta + A[1, 4] = bA11 * ainv + D22 * ainv * alpha * alpha + + # The third equation for W + A[2, 0] = A12 * beta * ainv + A[2, 1] = (bA11 + A22) * alpha * ainv + D22 * alpha * ainv**3 + A[2, 2] = ( + -(bA11 * alpha * alpha + bA22 * beta * beta) - A22 * ainv * ainv - D22 * ainv**4 + ) + A[2, 3] = -bA22 * beta - D12 * beta * ainv * ainv + A[2, 4] = -bA11 * alpha - D22 * alpha * ainv * ainv + + # Fourth equation for theta + A[3, 0] = D33 * ainv * alpha * alpha + A[3, 1] = D12 * ainv * alpha * beta + A[3, 2] = -bA22 * beta - D12 * beta * ainv * ainv + A[3, 3] = -(D11 * beta * beta + D33 * alpha * alpha) - bA22 + A[3, 4] = -(D12 + D33) * alpha * beta + + # Fifth equation for phi + A[4, 0] = D33 * ainv * alpha * beta + A[4, 1] = bA11 * ainv + D22 * ainv * alpha * alpha + A[4, 2] = -bA11 * alpha - D22 * alpha * ainv * ainv + A[4, 3] = -(D33 + D12) * alpha * beta + A[4, 4] = -(D33 * beta * beta + D22 * alpha * alpha) - bA11 + sol = np.linalg.solve(A, rhs) + + U, V, W, theta, phi = sol + return U, V, W, theta, phi + + +def get_exact_solution(x, y, z, R=1.0, L=2.0, M=4, N=3): + theta = np.atan2(x, y) + s = R * theta + + alpha = 4.0 / R + beta = 3 * np.pi / L + + U, V, W, _, _ = compute_coefficients( + E=70e9, nu=0.3, R=1.0, L=2.0, t=0.01, ks=5.0 / 6.0, M=M, N=N + ) + + # This is in the coordinate system aligned with x = axial direction + # v = tangential, w = circumferential directions + axial = U * np.sin(alpha * s) * np.cos(beta * z) + tangent = V * np.cos(alpha * s) * np.sin(beta * z) + circum = W * np.sin(alpha * s) * np.sin(beta * z) + + # Use these values to find the local coordinates + u = circum * (x / R) + tangent * (y / R) + v = circum * (y / R) - tangent * (x / R) + w = axial + + return u, v, w + + +def write_vtu(mesh, conns, u, v, w, filename="cylinder_shell.vtu"): X = mesh.X nnodes = X.shape[0] conn = np.vstack(conns) if isinstance(conns, list) else conns