Add Neo-Hookean (NH) model and unit tests in FIM

fim/refactor/main_VFM.py

@@ -6,9 +6,15 @@
 import time
 
 import numpy as np
-from material_model import MaterialModel
 from scipy.optimize import least_squares
-from vws_models import central_differentiation, increase_matrix_size, map_elements_to_centraldiff, read_input_file
+
+from fim.refactor.material_model import MaterialModel
+from fim.refactor.vws_models import (
+    central_differentiation,
+    increase_matrix_size,
+    map_elements_to_centraldiff,
+    read_input_file,
+)
 
 # Setup logging
 logging.basicConfig(level=logging.INFO)
@@ -18,11 +24,12 @@
 DEFAULT_PATHS = {
     "linear": os.path.join(DATA_ROOT, "80um"),
     "hgo": os.path.join(DATA_ROOT, "HGO"),
+    "nh": os.path.join(DATA_ROOT, "NH"),
 }
 
 # CLI
 parser = argparse.ArgumentParser(description="Run FIM Material Model Evaluation")
-parser.add_argument("--model", type=str, choices=["linear", "hgo"], help="Material model type")
+parser.add_argument("--model", type=str, choices=["linear", "hgo", "nh"], help="Material model type")
 parser.add_argument("--data_path", type=str, help="Path to input data folder")
 args = parser.parse_args()
 
@@ -63,6 +70,16 @@ def residual(x):
                 displacement_field, X, Y, Z, C10, D1, k1, k2, kappa, volume_matrix, Force, L, H
             )
 
+    elif name == "nh":
+
+        def residual(x):
+            C10, D1 = x
+            L = material_model.get_parameter("L")
+            H = material_model.get_parameter("H")
+            Force = material_model.get_parameter("Force")
+            volume_matrix = material_model.get_parameter("volume_matrix")
+            return material_model.model_func(displacement_field, X, Y, Z, C10, D1, volume_matrix, Force, L, H)
+
     else:
         raise ValueError("Unknown material model type")
 
@@ -115,6 +132,19 @@ def load_hgo_fields(folder):
     return X, Y, Z, tensor_displacement_list, L, W, H, volume_matrix
 
 
+def load_nh_fields(folder):
+    """Loads NH-specific displacement, volume, and mesh dimensions."""
+    X, Y, Z, tensor_displacement_list, volume_matrix = load_common_fields(folder)
+
+    undeformed_nodes, connectivity = read_input_file(f"{folder}/335k_32um.inp")
+
+    L = abs(np.max(undeformed_nodes[:, 1]) - np.min(undeformed_nodes[:, 1]))
+    W = abs(np.max(undeformed_nodes[:, 2]) - np.min(undeformed_nodes[:, 2]))
+    H = abs(np.max(undeformed_nodes[:, 3]) - np.min(undeformed_nodes[:, 3]))
+
+    return X, Y, Z, tensor_displacement_list, L, W, H, volume_matrix
+
+
 if __name__ == "__main__":
     start_time = time.time()
 
@@ -150,7 +180,7 @@ def load_hgo_fields(folder):
 
         # Run sensitivity analysis
         deviation = 0.05
-        sens = linear_model.sensitivity_analysis(disp_tensor, X, Y, Z, volume_matrix, L, H, deviation)
+        sens = linear_model.sensitivity_analysis_linear(disp_tensor, X, Y, Z, volume_matrix, L, H, deviation)
 
     elif model_name == "hgo":
         # === HGO Model ===
@@ -178,4 +208,29 @@ def load_hgo_fields(folder):
         # logging.info(f"HGO model result: {result_hgo}")
         logging.info("HGO model result: C10 = %.2f, D1 = %.2e, kappa = %.3f", *result_hgo)
 
+    elif model_name == "nh":
+        # === NH Model ===
+        X, Y, Z, disp_tensor, L, W, H, volume_matrix = load_nh_fields(data_path)
+        nh_params = {
+            "C10": 267 * 0.95,
+            "D1": 8e-4 * 0.95,
+            # "k1": 2000,
+            # "k2": 5,
+            # "kappa": 0.05,
+            "L": L,
+            "W": W,
+            "H": H,
+            "volume_matrix": volume_matrix,
+            "Force": 1.05e-05,
+        }
+        nh_model = MaterialModel("nh", nh_params)
+
+        initial_guess = [nh_params["C10"], nh_params["D1"]]
+        bounds = ((100, 1e-6), (1000, 1e-3))
+
+        # Run optimization
+        result_nh = run_inverse_model(disp_tensor, X, Y, Z, volume_matrix, initial_guess, bounds, nh_model)
+        # logging.info(f"HGO model result: {result_hgo}")
+        logging.info("NH model result: C10 = %.2f, D1 = %.2e", *result_nh)
+
     logging.info(f"Total runtime: {time.time() - start_time}")

fim/refactor/material_model.py

@@ -5,7 +5,7 @@
 
 import logging
 
-from vws_models import calculate_VWS_hgo, calculate_VWS_linear, sensitivity_full
+from fim.refactor.vws_models import calculate_VWS_hgo, calculate_VWS_linear, calculate_VWS_nh, sensitivity_full_linear
 
 
 class MaterialModel:
@@ -22,12 +22,14 @@ def _get_material_model(self, model_name):
             return calculate_VWS_linear
         if model_name == "hgo":
             return calculate_VWS_hgo
+        if model_name == "nh":
+            return calculate_VWS_nh
         raise ValueError(f"Unsupported model name: {model_name}")
 
     def get_parameter(self, key, default=None):
         return self.params.get(key, default)
 
-    def sensitivity_analysis(self, tensor_displacement_list, X, Y, Z, volume_matrix, L, H, deviation=0.05):
+    def sensitivity_analysis_linear(self, tensor_displacement_list, X, Y, Z, volume_matrix, L, H, deviation=0.05):
         if self.name != "linear":
             raise NotImplementedError("Sensitivity analysis is only available for the 'linear' model.")
 
@@ -38,7 +40,7 @@ def sensitivity_analysis(self, tensor_displacement_list, X, Y, Z, volume_matrix,
         Gt = self.get_parameter("Gt")
         Force = self.get_parameter("Force")
 
-        return sensitivity_full(
+        return sensitivity_full_linear(
             tensor_displacement_list, E1, E2, v12, v23, Gt, X, Y, Z, Force, volume_matrix, L, H, deviation
         )
 

fim/refactor/vws_models.py

@@ -587,6 +587,26 @@ def calculate_VWS_hgo(tensor_displacement_list, X, Y, Z, C10, D1, k1, k2, kappa,
     return phi * 1e10
 
 
+def calculate_VWS_nh(tensor_displacement_list, X, Y, Z, C10, D1, volume_matrix, Force, L, H):
+    I3 = np.eye(3)
+    F = tensor_displacement_list + I3
+    J = np.linalg.det(F)
+    Finv = np.linalg.inv(F)
+
+    b = np.einsum("...ik,...jk->...ij", F, F)
+    c = np.einsum("...ji,...jk->...ik", F, F)
+    I1 = np.trace(c, axis1=-2, axis2=-1)[..., None, None]
+
+    sigma_iso = (2 * C10 * J ** (-5 / 3))[..., None, None] * (b - (I3 * I1 / 3))
+    sigma_vol = ((1 / D1) * (J - 1 / J))[..., None, None] * I3
+    sigma = sigma_iso + sigma_vol
+
+    pk1 = J[..., None, None] * np.einsum("...ij,...kj->...ik", sigma, Finv)
+
+    phi = calculate_VWS_virtual_work(pk1, X, Y, Z, volume_matrix, Force, L, H, mode="nh")
+    return phi * 1e10
+
+
 def calculate_VWS_virtual_work(pk1, X, Y, Z, volume_element, Force, L, H, mode):
     """Compute internal (IVW) and external (EVW) virtual work contributions for five
     predefined virtual fields, and return the residual vector phi for the given mode.
@@ -665,14 +685,26 @@ def calculate_VWS_virtual_work(pk1, X, Y, Z, volume_element, Force, L, H, mode):
         # HGO mode: use virtual fields 5 and 3
         phi = np.array([total_IVW_5 - evw_5, total_IVW_3])
 
+    elif mode == "nh":
+        du_star_5 = np.zeros_like(pk1)
+        du_star_5[..., 2, 0] = U_star_z_pw_vol_devX(X1, X2, X3, L, H)
+        du_star_5[..., 2, 1] = U_star_z_pw_vol_devY(X1, X2, X3, L, H)
+        du_star_5[..., 2, 2] = U_star_z_pw_vol_devZ(X1, X2, X3, L, H)
+        ivw_5 = np.sum(pk1 * du_star_5, axis=(-2, -1)) * volume_element
+        total_IVW_5 = np.sum(ivw_5)
+        evw_5 = -Force * U_star_z_pw_vol(0, 0, H, L, H)
+        phi = np.array([total_IVW_3, total_IVW_5 - evw_5])
+        phi = np.sqrt(phi[0] ** 2 + phi[1] ** 2)
+
     else:
         raise ValueError(f"Unsupported mode '{mode}', choose 'linear' or 'hgo'.")
 
-    # Apply scaling for numerical stability
     return phi
 
 
-def sensitivity_full(tensor_displacement_list, E1, E2, v12, v23, Gt, X, Y, Z, Force, volume_matrix, L, H, deviation):
+def sensitivity_full_linear(
+    tensor_displacement_list, E1, E2, v12, v23, Gt, X, Y, Z, Force, volume_matrix, L, H, deviation
+):
     """Perform sensitivity analysis for the linear material model using finite difference.
 
     Parameters: