Thermoelastic3d

engibench/problems/thermoelastic3d/README.md

@@ -0,0 +1,33 @@
+# ThermoElastic3D
+
+**Lead**: Gabriel Apaza @gapaza
+
+This is 3D topology optimization problem for minimizing weakly coupled thermo-elastic compliance subject to boundary conditions and a volume fraction constraint.
+
+
+## Design space
+This multi-physics topology optimization problem is governed by linear elasticity and steady-state heat conduction with a one-way coupling from the thermal domain to the elastic domain.
+The problem is defined over a cube 3D domain, where load elements and support elements are placed along the boundary to define a unique elastic condition.
+Similarly, heatsink elements are placed along the boundary to define a unique thermal condition.
+The design space is then defined by a 3D array representing density values (parameterized by DesignSpace = [0,1]^{nelx x nely x nelz}, where nelx, nely, and nelz denote the x, y, and z dimensions respectively).
+
+## Objectives
+The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material.
+Total compliance is defined as the sum of thermal compliance and structural compliance.
+
+## Simulator
+The simulation code is based on a Python adaptation of the popular 88-line topology optimization code, modified to handle the thermal domain in addition to thermal-elastic coupling.
+Optimization is conducted by reformulating the integer optimization problem as a continuous one (leveraging a SIMP approach), where a density filtering approach is used to prevent checkerboard-like artifacts.
+The optimization process itself operates by calculating the sensitivities of the design variables with respect to total compliance (done efficiently using the Adjoint method), calculating the sensitivities of the design variables with respect to the constraint value, and then updating the design variables by solving a convex-linear subproblem and taking a small step (using the method of moving asymptotes).
+The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold.
+
+## Conditions
+Problem conditions are defined by creating a python dict with the following info:
+- `fixed_elements`: Encodes a binary NxNxN matrix of the structurally fixed elements in the domain.
+- `force_elements_x`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the x-direction.
+- `force_elements_y`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the y-direction.
+- `force_elements_z`: Encodes a binary NxNxN matrix specifying elements that have a structural load in the z-direction.
+- `heatsink_elements`: Encodes a binary NxNxN matrix specifying elements that have a heat sink.
+- `volfrac`: Encodes the target volume fraction for the volume fraction constraint.
+- `rmin`: Encodes the filter size used in the optimization routine.
+- `weight`: Allows one to control which objective is optimized for. 1.0 Is pure structural optimization, while 0.0 is pure thermal optimization.

engibench/problems/thermoelastic3d/__init__.py

@@ -0,0 +1,5 @@
+"""thermoelastic3d problem module."""
+
+from engibench.problems.thermoelastic3d.v0 import ThermoElastic3D
+
+__all__ = ["ThermoElastic3D"]

engibench/problems/thermoelastic3d/model/fem_matrix_builder.py

@@ -0,0 +1,130 @@
+"""This module assembles the local stiffness matrices the elastic, thermal, and thermoelastic domains."""
+
+import numpy as np
+from numpy import typing as npt
+
+SHAPE_NORM = 0.125  # Hex8 shape function normalization factor
+
+
+def fe_melthm_3d(
+    nu: float, e: float, k: float, alpha: float
+) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+    """Build 3D Hex8 element matrices for thermo-elasticity.
+
+    Args:
+        nu (float): Poisson's ratio
+        e (float): Young's modulus
+        k (float): Thermal conductivity (isotropic)
+        alpha (float): Coefficient of thermal expansion
+
+    Returns:
+        ke    : 24x24 mechanical stiffness (Hex8, 3 dof/node)
+        k_eth :  8x8 thermal conductivity (Hex8, 1 dof/node)
+        c_ethm: 24x8 coupling mapping nodal temperatures to equivalent mechanical forces
+    """
+    # --- Gauss points (2x2x2) and weights ---
+    gp = np.array([-1 / np.sqrt(3), 1 / np.sqrt(3)])
+    w = np.array([1.0, 1.0])
+
+    # --- Hex8 shape functions and derivatives in natural coordinates ---
+    # Node order: (-,-,-), (+,-,-), (+,+,-), (-,+,-), (-,-,+), (+,-,+), (+,+,+), (-,+,+)
+    xi_nodes = np.array([-1, 1, 1, -1, -1, 1, 1, -1])
+    et_nodes = np.array([-1, -1, 1, 1, -1, -1, 1, 1])
+    ze_nodes = np.array([-1, -1, -1, -1, 1, 1, 1, 1])
+
+    def shape_fun_and_derivs(xi: float, eta: float, zeta: float) -> tuple[np.ndarray, np.ndarray]:
+        """Builds the shape function for the local elements.
+
+        Args:
+            xi (float): Gaussian x coordinate
+            eta (float): Gaussian y coordinate
+            zeta (float): Gaussian z coordinate
+
+        Returns:
+            n: (8,) shape functions
+            d_n_dxi: (8,3) derivatives w.r.t. [xi, eta, zeta]
+        """
+        n = SHAPE_NORM * (1 + xi_nodes * xi) * (1 + et_nodes * eta) * (1 + ze_nodes * zeta)
+
+        # Derivatives with respect to natural coords
+        d_n_dxi = np.zeros((8, 3))
+        d_n_dxi[:, 0] = SHAPE_NORM * xi_nodes * (1 + et_nodes * eta) * (1 + ze_nodes * zeta)  # ∂N/∂xi
+        d_n_dxi[:, 1] = SHAPE_NORM * et_nodes * (1 + xi_nodes * xi) * (1 + ze_nodes * zeta)  # ∂N/∂eta
+        d_n_dxi[:, 2] = SHAPE_NORM * ze_nodes * (1 + xi_nodes * xi) * (1 + et_nodes * eta)  # ∂N/∂zeta
+        return n, d_n_dxi
+
+    # --- Geometry & Jacobian ---
+    # For a unit cube in physical space mapped from [-1,1]^3:
+    # x = (xi+1)/2, y = (eta+1)/2, z = (zeta+1)/2 -> J = diag(0.5, 0.5, 0.5)
+    j = np.diag([0.5, 0.5, 0.5])
+    det_j = np.linalg.det(j)  # = 0.125
+    inv_j = np.linalg.inv(j)  # = diag(2,2,2)
+
+    # --- Elasticity matrix (Voigt 6x6) for 3D isotropic ---
+    lam = e * nu / ((1 + nu) * (1 - 2 * nu))
+    mu = e / (2 * (1 + nu))
+    d = np.array(
+        [
+            [lam + 2 * mu, lam, lam, 0, 0, 0],
+            [lam, lam + 2 * mu, lam, 0, 0, 0],
+            [lam, lam + 2 * mu, lam + 2 * mu, 0, 0, 0],
+            [0, 0, 0, mu, 0, 0],
+            [0, 0, 0, 0, mu, 0],
+            [0, 0, 0, 0, 0, mu],
+        ],
+        dtype=float,
+    )
+    d[2, 1] = lam
+
+    # Thermal "volumetric" strain direction in Voigt
+    e_th = np.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0])
+
+    # --- Allocate element matrices ---
+    ke = np.zeros((24, 24))
+    k_eth = np.zeros((8, 8))
+    c_ethm = np.zeros((24, 8))
+
+    # --- Integration loop ---
+    for i, xi in enumerate(gp):
+        for j_idx, eta in enumerate(gp):
+            for kq, zeta in enumerate(gp):
+                n, d_n_dxi = shape_fun_and_derivs(xi, eta, zeta)
+
+                # Gradients in physical coords: dN_dx = dN_dxi * invJ
+                d_n_dx = d_n_dxi @ inv_j  # (8,3)
+
+                # Build B-matrix (6 x 24) for Hex8, 3 dof/node (u,v,w)
+                b = np.zeros((6, 24))
+                for a in range(8):
+                    ix = 3 * a
+                    dy, dx_, dz = d_n_dx[a, 1], d_n_dx[a, 0], d_n_dx[a, 2]  # clarity
+                    # normal strains
+                    b[0, ix + 0] = dx_
+                    b[1, ix + 1] = dy
+                    b[2, ix + 2] = dz
+                    # shear strains (engineering)
+                    b[3, ix + 0] = dy
+                    b[3, ix + 1] = dx_
+                    b[4, ix + 1] = dz
+                    b[4, ix + 2] = dy
+                    b[5, ix + 0] = dz
+                    b[5, ix + 2] = dx_
+
+                # Thermal gradient matrix for conduction: G = grad(N) (3 x 8)
+                g = d_n_dx.T  # rows: [dN/dx; dN/dy; dN/dz]
+
+                # Weight factor
+                wt = w[i] * w[j_idx] * w[kq] * det_j
+
+                # --- Accumulate ---
+                # Mechanical stiffness
+                ke += (b.T @ d @ b) * wt
+
+                # Thermal conductivity (isotropic k * grad N^T grad N)
+                k_eth += (g.T @ g) * (k * wt)
+
+                # Thermo-mech coupling: columns correspond to nodal temperatures (via N_j)
+                de_th = d @ (alpha * e_th)  # (6,)
+                c_ethm += (b.T @ de_th[:, None] @ n[None, :]) * wt
+
+    return ke, k_eth, c_ethm