Added API to get homogenized stress and homogenized tangent (#31)

* worked on getting the algorithmic tangent via finite differences

* linear elastic triclinic material model

* refactored material properties from type map<string, vector<double>> to type json

* J2Plasticity time_step is double not a vector

* refactored compute_error to accept different error measure and type

* added additional error control when get_homogenized_tangent is called

* in get_homogenized_tangent linear materials are handled seperately

CMakeLists.txt

@@ -134,7 +134,7 @@ set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/setup.h
 
         include/material_models/LinearThermalIsotropic.h
-        include/material_models/LinearElasticIsotropic.h
+        include/material_models/LinearElastic.h
         include/material_models/PseudoPlastic.h
         include/material_models/J2Plasticity.h
 )

FANS_Dashboard/PlotYoungsModulus.py

@@ -1,76 +1,89 @@
 import numpy as np
 import plotly.graph_objs as go
+import meshio
 
 
-def compute_3d_youngs_modulus(C):
+def compute_YoungsModulus3D(C_batch):
     """
-    Compute Young's modulus for all directions in 3D.
+    Compute Young's modulus for all directions in 3D for a batch of stiffness tensors.
 
-    Parameters:
-    C : ndarray
-        Stiffness tensor in Mandel notation.
+    Args:
+        C_batch (ndarray): Batch of stiffness tensors in Mandel notation, shape (n, 6, 6).
 
     Returns:
-    E: ndarray
-        Young's modulus in all directions.
-    X, Y, Z: ndarrays
-        Coordinates for plotting the modulus surface.
+        tuple: A tuple containing:
+            - X_batch (ndarray): X-coordinates for plotting the modulus surface, shape (n, n_theta, n_phi).
+            - Y_batch (ndarray): Y-coordinates for plotting the modulus surface, shape (n, n_theta, n_phi).
+            - Z_batch (ndarray): Z-coordinates for plotting the modulus surface, shape (n, n_theta, n_phi).
+            - E_batch (ndarray): Young's modulus in all directions, shape (n, n_theta, n_phi).
     """
-
+    n = C_batch.shape[0]
     n_theta = 180
     n_phi = 360
-    theta = np.linspace(0, np.pi, n_theta)  # Polar angle
-    phi = np.linspace(0, 2 * np.pi, n_phi)  # Azimuthal angle
 
-    S = np.linalg.inv(C)
+    theta = np.linspace(0, np.pi, n_theta)
+    phi = np.linspace(0, 2 * np.pi, n_phi)
+    theta_grid, phi_grid = np.meshgrid(theta, phi, indexing="ij")
+
+    d_x = np.sin(theta_grid) * np.cos(phi_grid)  # Shape (n_theta, n_phi)
+    d_y = np.sin(theta_grid) * np.sin(phi_grid)
+    d_z = np.cos(theta_grid)
+
+    N = np.stack(
+        (
+            d_x**2,
+            d_y**2,
+            d_z**2,
+            np.sqrt(2) * d_x * d_y,
+            np.sqrt(2) * d_x * d_z,
+            np.sqrt(2) * d_y * d_z,
+        ),
+        axis=-1,
+    )  # Shape (n_theta, n_phi, 6)
 
-    E = np.zeros((n_theta, n_phi))
+    N_flat = N.reshape(-1, 6)  # Shape (n_points, 6)
 
-    for i in range(n_theta):
-        for j in range(n_phi):
-            d = np.array(
-                [
-                    np.sin(theta[i]) * np.cos(phi[j]),
-                    np.sin(theta[i]) * np.sin(phi[j]),
-                    np.cos(theta[i]),
-                ]
-            )
+    # Invert stiffness tensors to get compliance tensors
+    S_batch = np.linalg.inv(C_batch)  # Shape (n, 6, 6)
 
-            N = np.array(
-                [
-                    d[0] ** 2,
-                    d[1] ** 2,
-                    d[2] ** 2,
-                    np.sqrt(2.0) * d[0] * d[1],
-                    np.sqrt(2.0) * d[0] * d[2],
-                    np.sqrt(2.0) * d[2] * d[1],
-                ]
-            )
+    # Compute E for each tensor in the batch
+    NSN = np.einsum("pi,nij,pj->np", N_flat, S_batch, N_flat)  # Shape (n, n_points)
+    E_batch = 1.0 / NSN  # Shape (n, n_points)
 
-            E[i, j] = 1.0 / (N.T @ S @ N)
+    # Reshape E_batch back to (n, n_theta, n_phi)
+    E_batch = E_batch.reshape(n, *d_x.shape)
 
-    X = E * np.sin(theta)[:, np.newaxis] * np.cos(phi)[np.newaxis, :]
-    Y = E * np.sin(theta)[:, np.newaxis] * np.sin(phi)[np.newaxis, :]
-    Z = E * np.cos(theta)[:, np.newaxis]
+    X_batch = E_batch * d_x  # Shape (n, n_theta, n_phi)
+    Y_batch = E_batch * d_y
+    Z_batch = E_batch * d_z
 
-    return X, Y, Z, E
+    return X_batch, Y_batch, Z_batch, E_batch
 
 
-def plot_3d_youngs_modulus_surface(C, title="Young's Modulus Surface"):
+def plot_YoungsModulus3D(C, title="Young's Modulus Surface"):
     """
     Plot a 3D surface of Young's modulus.
 
-    Parameters:
-    C : ndarray
-        Stiffness tensor in Mandel notation.
-    title : str
-        Title of the plot.
+    Args:
+        C (ndarray): Stiffness tensor in Mandel notation. Can be a single tensor of shape (6,6) or a batch of tensors of shape (n,6,6).
+        title (str): Title of the plot.
 
+    Raises:
+        ValueError: If C is not of shape (6,6) or (1,6,6).
     """
-    X, Y, Z, E = compute_3d_youngs_modulus(C)
+    if C.shape == (6, 6):
+        C_batch = C[np.newaxis, :, :]
+    elif C.shape == (1, 6, 6):
+        C_batch = C
+    else:
+        raise ValueError(
+            "C must be either a (6,6) tensor or a batch with one tensor of shape (1,6,6)."
+        )
+
+    X_batch, Y_batch, Z_batch, E_batch = compute_YoungsModulus3D(C_batch)
+    X, Y, Z, E = X_batch[0], Y_batch[0], Z_batch[0], E_batch[0]
 
     surface = go.Surface(x=X, y=Y, z=Z, surfacecolor=E, colorscale="Viridis")
-
     layout = go.Layout(
         title=title,
         scene=dict(
@@ -85,14 +98,64 @@ def plot_3d_youngs_modulus_surface(C, title="Young's Modulus Surface"):
     fig.show()
 
 
+def export_YoungsModulus3D_to_vtk(C, prefix="youngs_modulus_surface"):
+    """
+    Export the computed Young's modulus surfaces to VTK files for Paraview visualization.
+
+    Args:
+        C (ndarray): Stiffness tensor in Mandel notation. Can be a single tensor of shape (6,6) or a batch of tensors of shape (n,6,6).
+        prefix (str): Prefix for the output files.
+
+    Returns:
+        None
+    """
+    X_batch, Y_batch, Z_batch, E_batch = compute_YoungsModulus3D(C)
+    n, n_theta, n_phi = X_batch.shape
+
+    for k in range(n):
+        points = np.vstack(
+            (X_batch[k].ravel(), Y_batch[k].ravel(), Z_batch[k].ravel())
+        ).T
+        cells = [
+            (
+                "quad",
+                np.array(
+                    [
+                        [
+                            i * n_phi + j,
+                            (i + 1) * n_phi + j,
+                            (i + 1) * n_phi + (j + 1),
+                            i * n_phi + (j + 1),
+                        ]
+                        for i in range(n_theta - 1)
+                        for j in range(n_phi - 1)
+                    ],
+                    dtype=np.int32,
+                ),
+            )
+        ]
+        mesh = meshio.Mesh(
+            points=points,
+            cells=cells,
+            point_data={"Youngs_Modulus": E_batch[k].ravel()},
+        )
+        filename = f"{prefix}_{k}.vtk"
+        meshio.write(filename, mesh)
+        print(f"Exported {filename}")
+
+
 def demoCubic():
     """
-    Demonstrates the Young's modulus surface plotting routine for a cubic material (Copper)
+    Demonstrates the Young's modulus surface plotting routine for a cubic material (Copper).
+
+    This function generates the stiffness tensor for a cubic material, specifically copper,
+    and then plots the 3D Young's modulus surface using the generated tensor.
 
-    Returns
-    -------
-    None.
+    Args:
+        None
 
+    Returns:
+        None
     """
     P1 = np.zeros((6, 6))
     P1[:3, :3] = 1.0 / 3.0
@@ -104,7 +167,5 @@ def demoCubic():
     l1, l2, l3 = 136.67, 46, 150
     C = 3 * l1 * P1 + l2 * P2 + l3 * P3
 
-    print(C)
-
     # show the 3D Young's modulus plot for copper
-    plot_3d_youngs_modulus_surface(C, title="Young's Modulus Surface for Copper")
+    plot_YoungsModulus3D(C, title="Young's Modulus Surface for Copper")

include/general.h

@@ -16,6 +16,11 @@
 
 using namespace std;
 
+// JSON
+#include <json.hpp>
+using nlohmann::json;
+using namespace nlohmann;
+
 // Packages
 #include "fftw3-mpi.h"
 #include "fftw3.h" // this includes the serial fftw as well as mpi header files! See https://fftw.org/doc/MPI-Files-and-Data-Types.html

include/material_models/J2Plasticity.h

@@ -6,17 +6,17 @@
 
 class J2Plasticity : public MechModel {
   public:
-    J2Plasticity(vector<double> l_e, map<string, vector<double>> materialProperties)
+    J2Plasticity(vector<double> l_e, json materialProperties)
         : MechModel(l_e)
     {
         try {
-            bulk_modulus  = materialProperties["bulk_modulus"];
-            shear_modulus = materialProperties["shear_modulus"];
-            yield_stress  = materialProperties["yield_stress"];                  // Initial yield stress
-            K             = materialProperties["isotropic_hardening_parameter"]; // Isotropic hardening parameter
-            H             = materialProperties["kinematic_hardening_parameter"]; // Kinematic hardening parameter
-            eta           = materialProperties["viscosity"];                     // Viscosity parameter
-            dt            = materialProperties["time_step"][0];                  // Time step
+            bulk_modulus  = materialProperties["bulk_modulus"].get<vector<double>>();
+            shear_modulus = materialProperties["shear_modulus"].get<vector<double>>();
+            yield_stress  = materialProperties["yield_stress"].get<vector<double>>();                  // Initial yield stress
+            K             = materialProperties["isotropic_hardening_parameter"].get<vector<double>>(); // Isotropic hardening parameter
+            H             = materialProperties["kinematic_hardening_parameter"].get<vector<double>>(); // Kinematic hardening parameter
+            eta           = materialProperties["viscosity"].get<vector<double>>();                     // Viscosity parameter
+            dt            = materialProperties["time_step"].get<double>();                             // Time step
         } catch (const std::out_of_range &e) {
             throw std::runtime_error("Missing material properties for the requested material model.");
         }
@@ -158,7 +158,7 @@ class J2Plasticity : public MechModel {
 // Derived Class Linear Isotropic Hardening
 class J2ViscoPlastic_LinearIsotropicHardening : public J2Plasticity {
   public:
-    J2ViscoPlastic_LinearIsotropicHardening(vector<double> l_e, map<string, vector<double>> materialProperties)
+    J2ViscoPlastic_LinearIsotropicHardening(vector<double> l_e, json materialProperties)
         : J2Plasticity(l_e, materialProperties)
     {
     }
@@ -177,12 +177,12 @@ class J2ViscoPlastic_LinearIsotropicHardening : public J2Plasticity {
 // Derived Class Non-Linear (Exponential law) Isotropic Hardening
 class J2ViscoPlastic_NonLinearIsotropicHardening : public J2Plasticity {
   public:
-    J2ViscoPlastic_NonLinearIsotropicHardening(vector<double> l_e, map<string, vector<double>> materialProperties)
+    J2ViscoPlastic_NonLinearIsotropicHardening(vector<double> l_e, json materialProperties)
         : J2Plasticity(l_e, materialProperties)
     {
         try {
-            sigma_inf = materialProperties["saturation_stress"];   // Saturation stress
-            delta     = materialProperties["saturation_exponent"]; // Saturation exponent
+            sigma_inf = materialProperties["saturation_stress"].get<vector<double>>();   // Saturation stress
+            delta     = materialProperties["saturation_exponent"].get<vector<double>>(); // Saturation exponent
         } catch (const std::out_of_range &e) {
             throw std::runtime_error("Missing material properties for the requested material model.");
         }