Merge pull request #33 from DataAnalyticsEngineering/FANS-v0.3.0

* Mention dependency on compiler support for OpenMP in README

* Rename container user also in fans_dev image to fix broken entrypoint script (#27)

* Add python3-dev to fans-ci stage already

* Remove clang-format check via GitHub Action (#29)

* Bump clang-format version

* Formatting

* Use a different Action workflow to run clang-format

* Remove the Action to run clang-format, as it will go out of sync with the pre-commit run

* Bump mirror-clang-format to latest version in the pre-commit-config

* Update clang-format config file to also format JSON

* Manually specify which extensions to format in the mirror-clang-format hooks in the pre-commit config

* Exclude JSON files from clang-format

* Exclude JSON files from clang-format

* Correctly exclude JSON files from clang-format checking

* Rather than excluding json, specifically suggest c++

* 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

* added linear thermal triclinic model  (#32)

* added linear thermal triclinic model similar to linear elastic triclinic model

* fixed readme.md indent causing linting error

* Update CHANGELOG

* Bump version

---------

Co-authored-by: Moritz Sigg <[email protected]>
Co-authored-by: Moritz Sigg <[email protected]>
Co-authored-by: Ishaan Desai <[email protected]>

Author: 3 people

.github/workflows/clang-format-check.yml

@@ -16,8 +16,9 @@ repos:
     - id: black
 # clang-format for C/C++ formatting
 -   repo: https://github.com/pre-commit/mirrors-clang-format
-    rev: v8.0.1
+    rev: v19.1.2
     hooks:
     -   id: clang-format
         args: ['--style=file']
         exclude: "include/json.hpp"
+        types_or: [c++]

.pre-commit-config.yaml

@@ -1,5 +1,10 @@
 # FANS Changelog
 
+## v0.3.0
+
+- Added Linear thermal and mechanical triclinic material models https://github.com/DataAnalyticsEngineering/FANS/pull/32
+- Added API to get homogenized stress and homogenized tangent https://github.com/DataAnalyticsEngineering/FANS/pull/31
+
 ## v0.2.0
 
 - Add integration tests https://github.com/DataAnalyticsEngineering/FANS/pull/20

CHANGELOG.md

@@ -5,7 +5,7 @@ cmake_minimum_required(VERSION 3.0...3.28)
 # ##############################################################################
 
 project(FANS
-        VERSION 0.2.0
+        VERSION 0.3.0
         LANGUAGES C CXX
 )
 
@@ -133,8 +133,8 @@ set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/solver.h
         include/setup.h
 
-        include/material_models/LinearThermalIsotropic.h
-        include/material_models/LinearElasticIsotropic.h
+        include/material_models/LinearThermal.h
+        include/material_models/LinearElastic.h
         include/material_models/PseudoPlastic.h
         include/material_models/J2Plasticity.h
 )

CMakeLists.txt

@@ -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")

FANS_Dashboard/PlotYoungsModulus.py

@@ -16,7 +16,7 @@ Fourier Accelerated Nodal Solvers (FANS) is an FFT-based homogenization solver d
 
 FANS has the following dependencies:
 
-- A C++ compiler (e.g. GCC)
+- A C++ compiler with OpenMP support (e.g. GCC, or Clang with OpenMP libraries installed)
 - CMake (version 3.0 or higher) (+ Unix file utility for creating .deb packages)
 - Git (for cloning this repo)
 - MPI (mpicc and mpic++)
@@ -178,13 +178,20 @@ FANS requires a JSON input file specifying the problem parameters. Example input
 
 ```json
 "method": "cg",
-"TOL": 1e-10,
-"n_it": 100
+"error_parameters":{
+    "measure": "Linfinity",
+    "type": "absolute",
+    "tolerance": 1e-10
+},
+"n_it": 100,
 ```
 
 - `method`: This indicates the numerical method to be used for solving the system of equations. `cg` stands for the Conjugate Gradient method, and `fp` stands for the Fixed Point method.
-- `TOL`: This sets the tolerance level for the solver. It defines the convergence criterion which is based on the L-infinity norm of the nodal finite element residual; the solver iterates until the solution meets this tolerance.
-- `n_it`: This specifies the maximum number of iterations allowed for the FANS solver.
+- `error_parameters`: This section defines the error parameters for the solver. Error control is applied on the finite element nodal residual of the problem.
+  - `measure`: Specifies the norm used to measure the error. Options include `Linfinity`, `L1`, or `L2`.
+  - `type`: Defines the type of error measurement. Options are `absolute` or `relative`.
+  - `tolerance`: Sets the tolerance level for the solver, defining the convergence criterion based on the chosen error measure. The solver iterates until the solution meets this tolerance.
+- `n_it`: Specifies the maximum number of iterations allowed for the FANS solver.
 
 ### Macroscale Loading Conditions
 

README.md

@@ -40,6 +40,8 @@ RUN apt-get update -qq && apt-get install -y --no-install-recommends \
     libeigen3-dev \
     libfftw3-dev \
     libfftw3-mpi-dev \
+# Required for preCICE Micro Manager Python bindings
+    python3-dev \
 # Clean up
     && apt-get clean \
     && apt-get autoremove --purge -y \
@@ -60,11 +62,9 @@ RUN apt-get update -qq && apt-get install -y --no-install-recommends \
     time \
     htop \
     vim \
-    python3 \
     python3-pip \
     python3-venv \
     python-is-python3 \
-    python3-dev \
 # Clean up
     && apt-get clean \
     && apt-get autoremove --purge -y \

docker/Dockerfile

@@ -8,7 +8,7 @@ set -e
 # USAGE: docker run -e HOST_UID=$(id -u) -e HOST_GID=$(id -g) ...
 # open issue on this topic: https://github.com/docker/roadmap/issues/398
 hostgroup="hostgroup"
-container_user="develop"
+container_user="fans"
 
 if [ "$(id -u -n)" = "root" ]; then
     if [ -n "$HOST_UID" ] && [ -n "$HOST_GID" ]; then

docker/Dockerfile_user_env_entrypoint.sh

@@ -3,7 +3,7 @@
 We provide a set of docker images for different use cases on our [Dockerhub profile](https://hub.docker.com/u/unistuttgartdae):
 
 - **fans-ci**: Contains the minimum tools to build FANS (including dev packages of dependencies with the required headers), but does not include FANS itself. Meant for a CI workflow.
-- **fans-dev**: Based upon fans-ci, but offers a non-root user (`develop`) and handling of UID and GID to not mess up permissions when volume mounting into the container. Meant as an quick to setup build environment for FANS.
+- **fans-dev**: Based upon fans-ci, but offers a non-root user (`fans`) and handling of UID and GID to not mess up permissions when volume mounting into the container. Meant as an quick to setup build environment for FANS.
 
 Both images are built for linux/amd64 and linux/arm64 as well as for the three most recent Ubuntu LTS versions (focal, jammy, noble). The Ubuntu version can be selected through tags, e.g. `fans-dev:focal`; `noble` is equivalent to the `latest` tag. The architecture is selected automatically depending on your host platform.
 
@@ -75,7 +75,7 @@ You can attach VS Code to the newly created container in order to actually work
 
 To attach VS Code you need to install the `Remote Development Extension Pack` and the `Docker` Extension. Then open the Docker menu, right click our newly created `fans-dev` container and select "Start" (if not running already) and then "Attach Visual Studio Code".
 
-After attaching VS Code you unfortunately are user `root` in VS Code due to the way the UID and GID mapping is implemented: The container starts as root, executes the entrypoint script which changes UID and GID and only then drops privileges using `gosu`. VS Code though skips the entrypoint script and thus doesn't switch to the non-root user `develop`. You however can do so manually by typing `gosu develop bash` in your terminal sessions inside VS Code.
+After attaching VS Code you unfortunately are user `root` in VS Code due to the way the UID and GID mapping is implemented: The container starts as root, executes the entrypoint script which changes UID and GID and only then drops privileges using `gosu`. VS Code though skips the entrypoint script and thus doesn't switch to the non-root user `fans`. You however can do so manually by typing `gosu fans bash` in your terminal sessions inside VS Code.
 
 For further reading and alternative approaches like a full DevContainer setup have a look at
 
@@ -87,10 +87,10 @@ For further reading and alternative approaches like a full DevContainer setup ha
 
 By building inside the container, FANS is linked against the container's libs and therefore must run inside the container. After attaching to the container you can then continue to use FANS as described in the main [README](../README.md#usage). Just remember that any input and output files need to visible to the container and thus must lie somewhere inside the mounted volumes.
 
-Special care has to be taken if you need to use FANS within scripts on the host, as Docker's interactive mode (`-i`) is not suitable in this case. Instead you need to use `docker exec`. One basically replaces the original `FANS` call by `docker exec -u develop -w /FANS/test fans-dev [original call]`. For example in conjunction with nohup:
+Special care has to be taken if you need to use FANS within scripts on the host, as Docker's interactive mode (`-i`) is not suitable in this case. Instead you need to use `docker exec`. One basically replaces the original `FANS` call by `docker exec -u fans -w /FANS/test fans-dev [original call]`. For example in conjunction with nohup:
 
 ```bash
 docker start fans-dev
-nohup /usr/bin/time -v docker exec -u develop -w /FANS/test fans-dev [original call] &
+nohup /usr/bin/time -v docker exec -u fans -w /FANS/test fans-dev [original call] &
 docker stop fans-dev
 ```

docker/README.md

@@ -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
@@ -49,4 +54,4 @@ inline V *FANS_malloc(size_t n)
 
 #define VERBOSITY 0
 
-//#define EIGEN_RUNTIME_NO_MALLOC
+// #define EIGEN_RUNTIME_NO_MALLOC