diff --git a/.github/workflows/clang-format-check.yml b/.github/workflows/clang-format-check.yml deleted file mode 100644 index fd6d7dd..0000000 --- a/.github/workflows/clang-format-check.yml +++ /dev/null @@ -1,21 +0,0 @@ -name: clang-format Check -on: - push: - branches: - - main - - develop - pull_request: - -jobs: - formatting-check: - name: Formatting Check - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - name: Run clang-format style check for C/C++/Protobuf programs. - uses: jidicula/clang-format-action@v4.13.0 - with: - clang-format-version: '13' - check-path: ${{ matrix.path }} - exclude-regex: 'include/json.hpp' - fallback-style: 'Mozilla' # optional diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 0240fee..6123ab2 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -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++] diff --git a/CHANGELOG.md b/CHANGELOG.md index f94b6f8..9d3656e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1918e68..c652e7b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 ) diff --git a/FANS_Dashboard/PlotYoungsModulus.py b/FANS_Dashboard/PlotYoungsModulus.py index 124e682..b6e5789 100644 --- a/FANS_Dashboard/PlotYoungsModulus.py +++ b/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") diff --git a/README.md b/README.md index 2c56764..b32f81d 100644 --- a/README.md +++ b/README.md @@ -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 diff --git a/docker/Dockerfile b/docker/Dockerfile index 5398805..4760e30 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -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 \ diff --git a/docker/Dockerfile_user_env_entrypoint.sh b/docker/Dockerfile_user_env_entrypoint.sh index 2f7ecc3..7fbdc79 100644 --- a/docker/Dockerfile_user_env_entrypoint.sh +++ b/docker/Dockerfile_user_env_entrypoint.sh @@ -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 diff --git a/docker/README.md b/docker/README.md index abb252d..ce27744 100644 --- a/docker/README.md +++ b/docker/README.md @@ -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 ``` diff --git a/include/general.h b/include/general.h index c164cff..395d089 100644 --- a/include/general.h +++ b/include/general.h @@ -16,6 +16,11 @@ using namespace std; +// JSON +#include +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 diff --git a/include/material_models/J2Plasticity.h b/include/material_models/J2Plasticity.h index 5551b76..ea918ed 100644 --- a/include/material_models/J2Plasticity.h +++ b/include/material_models/J2Plasticity.h @@ -6,17 +6,17 @@ class J2Plasticity : public MechModel { public: - J2Plasticity(vector l_e, map> materialProperties) + J2Plasticity(vector 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>(); + shear_modulus = materialProperties["shear_modulus"].get>(); + yield_stress = materialProperties["yield_stress"].get>(); // Initial yield stress + K = materialProperties["isotropic_hardening_parameter"].get>(); // Isotropic hardening parameter + H = materialProperties["kinematic_hardening_parameter"].get>(); // Kinematic hardening parameter + eta = materialProperties["viscosity"].get>(); // Viscosity parameter + dt = materialProperties["time_step"].get(); // 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 l_e, map> materialProperties) + J2ViscoPlastic_LinearIsotropicHardening(vector 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 l_e, map> materialProperties) + J2ViscoPlastic_NonLinearIsotropicHardening(vector 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>(); // Saturation stress + delta = materialProperties["saturation_exponent"].get>(); // Saturation exponent } catch (const std::out_of_range &e) { throw std::runtime_error("Missing material properties for the requested material model."); } diff --git a/include/material_models/LinearElastic.h b/include/material_models/LinearElastic.h new file mode 100644 index 0000000..94d6152 --- /dev/null +++ b/include/material_models/LinearElastic.h @@ -0,0 +1,150 @@ +#ifndef LINEARELASTIC_H +#define LINEARELASTIC_H + +#include "matmodel.h" +#include // For Eigen's aligned_allocator + +class LinearElasticIsotropic : public MechModel, public LinearModel<3> { + public: + LinearElasticIsotropic(vector l_e, json materialProperties) + : MechModel(l_e) + { + try { + bulk_modulus = materialProperties["bulk_modulus"].get>(); + mu = materialProperties["shear_modulus"].get>(); + } catch (const std::exception &e) { + throw std::runtime_error("Missing material properties for the requested material model."); + } + + n_mat = bulk_modulus.size(); + lambda.resize(n_mat); + mu.resize(n_mat); + + for (int i = 0; i < n_mat; ++i) { + lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i]; + } + + double lambda_ref = (*max_element(lambda.begin(), lambda.end()) + + *min_element(lambda.begin(), lambda.end())) / + 2; + double mu_ref = (*max_element(mu.begin(), mu.end()) + + *min_element(mu.begin(), mu.end())) / + 2; + + kapparef_mat = Matrix::Zero(); + kapparef_mat.topLeftCorner(3, 3).setConstant(lambda_ref); + kapparef_mat += 2 * mu_ref * Matrix::Identity(); + + phase_stiffness = new Matrix[n_mat]; + Matrix phase_kappa; + + for (int i = 0; i < n_mat; i++) { + phase_kappa.setZero(); + phase_stiffness[i] = Matrix::Zero(); + + phase_kappa.topLeftCorner(3, 3).setConstant(lambda[i]); + phase_kappa += 2 * mu[i] * Matrix::Identity(); + + for (int p = 0; p < 8; ++p) { + phase_stiffness[i] += B_int[p].transpose() * phase_kappa * B_int[p] * v_e * 0.1250; + } + } + } + + void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override + { + double buf1 = lambda[mat_index] * (eps(i, 0) + eps(i + 1, 0) + eps(i + 2, 0)); + double buf2 = 2 * mu[mat_index]; + sigma(i + 0, 0) = buf1 + buf2 * eps(i + 0, 0); + sigma(i + 1, 0) = buf1 + buf2 * eps(i + 1, 0); + sigma(i + 2, 0) = buf1 + buf2 * eps(i + 2, 0); + sigma(i + 3, 0) = buf2 * eps(i + 3, 0); + sigma(i + 4, 0) = buf2 * eps(i + 4, 0); + sigma(i + 5, 0) = buf2 * eps(i + 5, 0); + } + + private: + vector bulk_modulus; + vector lambda; + vector mu; +}; + +class LinearElasticTriclinic : public MechModel, public LinearModel<3> { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Ensure proper alignment for Eigen structures + + LinearElasticTriclinic(vector l_e, json materialProperties) + : MechModel(l_e) + { + vector C_keys = { + "C_11", "C_12", "C_13", "C_14", "C_15", "C_16", + "C_22", "C_23", "C_24", "C_25", "C_26", + "C_33", "C_34", "C_35", "C_36", + "C_44", "C_45", "C_46", + "C_55", "C_56", + "C_66"}; + + try { + n_mat = materialProperties.at("C_11").get>().size(); + size_t num_constants = C_keys.size(); + + // Initialize matrix to hold all constants + C_constants.resize(num_constants, n_mat); + + // Load material constants into matrix + for (size_t k = 0; k < num_constants; ++k) { + const auto &values = materialProperties.at(C_keys[k]).get>(); + if (values.size() != n_mat) { + throw std::runtime_error("Inconsistent size for material property: " + C_keys[k]); + } + C_constants.row(k) = Eigen::Map(values.data(), values.size()); + } + } catch (const std::exception &e) { + throw std::runtime_error("Missing or inconsistent material properties for the requested material model."); + } + + // Assemble stiffness matrices for each material + C_mats.resize(n_mat); + kapparef_mat = Matrix::Zero(); + + for (size_t i = 0; i < n_mat; ++i) { + Matrix C_i = Matrix::Zero(); + int k = 0; // Index for C_constants + + // Assign constants to the upper triangular part + for (int row = 0; row < 6; ++row) { + for (int col = row; col < 6; ++col) { + C_i(row, col) = C_constants(k++, i); + } + } + + // Symmetrize the matrix + C_i = C_i.selfadjointView(); + + C_mats[i] = C_i; + kapparef_mat += C_i; + } + + kapparef_mat /= n_mat; + + // Compute phase stiffness matrices + phase_stiffness = new Matrix[n_mat]; + for (size_t i = 0; i < n_mat; ++i) { + phase_stiffness[i] = Matrix::Zero(); + for (int p = 0; p < 8; ++p) { + phase_stiffness[i] += B_int[p].transpose() * C_mats[i] * B_int[p] * v_e * 0.1250; + } + } + } + + void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override + { + sigma.segment<6>(i) = C_mats[mat_index] * eps.segment<6>(i); + } + + private: + std::vector, Eigen::aligned_allocator>> C_mats; + MatrixXd C_constants; +}; + +#endif // LINEARELASTIC_H diff --git a/include/material_models/LinearElasticIsotropic.h b/include/material_models/LinearElasticIsotropic.h deleted file mode 100644 index 684ba79..0000000 --- a/include/material_models/LinearElasticIsotropic.h +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef LINEARELASTICISOTROPIC_H -#define LINEARELASTICISOTROPIC_H - -#include "matmodel.h" - -class LinearElasticIsotropic : public MechModel, public LinearModel<3> { - public: - LinearElasticIsotropic(vector l_e, map> materialProperties) - : MechModel(l_e) - { - try { - bulk_modulus = materialProperties["bulk_modulus"]; - mu = materialProperties["shear_modulus"]; - } catch (const std::exception &e) { - throw std::runtime_error("Missing material properties for the requested material model."); - } - - n_mat = bulk_modulus.size(); - lambda.resize(n_mat); - mu.resize(n_mat); - - for (int i = 0; i < n_mat; ++i) { - lambda[i] = bulk_modulus[i] - (2.0 / 3.0) * mu[i]; - } - - double lambda_ref = (*max_element(lambda.begin(), lambda.end()) + - *min_element(lambda.begin(), lambda.end())) / - 2; - double mu_ref = (*max_element(mu.begin(), mu.end()) + - *min_element(mu.begin(), mu.end())) / - 2; - - kapparef_mat = Matrix::Zero(); - kapparef_mat.topLeftCorner(3, 3).setConstant(lambda_ref); - kapparef_mat += 2 * mu_ref * Matrix::Identity(); - - phase_stiffness = new Matrix[n_mat]; - Matrix phase_kappa; - - for (int i = 0; i < n_mat; i++) { - phase_kappa.setZero(); - phase_stiffness[i] = Matrix::Zero(); - - phase_kappa.topLeftCorner(3, 3).setConstant(lambda[i]); - phase_kappa += 2 * mu[i] * Matrix::Identity(); - - for (int p = 0; p < 8; ++p) { - phase_stiffness[i] += B_int[p].transpose() * phase_kappa * B_int[p] * v_e * 0.1250; - } - } - } - - void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override - { - double buf1 = lambda[mat_index] * (eps(i, 0) + eps(i + 1, 0) + eps(i + 2, 0)); - double buf2 = 2 * mu[mat_index]; - sigma(i + 0, 0) = buf1 + buf2 * eps(i + 0, 0); - sigma(i + 1, 0) = buf1 + buf2 * eps(i + 1, 0); - sigma(i + 2, 0) = buf1 + buf2 * eps(i + 2, 0); - sigma(i + 3, 0) = buf2 * eps(i + 3, 0); - sigma(i + 4, 0) = buf2 * eps(i + 4, 0); - sigma(i + 5, 0) = buf2 * eps(i + 5, 0); - } - - private: - vector bulk_modulus; - vector lambda; - vector mu; -}; - -#endif // LINEARELASTICISOTROPIC_H diff --git a/include/material_models/LinearThermal.h b/include/material_models/LinearThermal.h new file mode 100644 index 0000000..bb8fc47 --- /dev/null +++ b/include/material_models/LinearThermal.h @@ -0,0 +1,115 @@ +#ifndef LINEARTHERMAL_H +#define LINEARTHERMAL_H + +#include "matmodel.h" +#include // For Eigen's aligned_allocator + +class LinearThermalIsotropic : public ThermalModel, public LinearModel<1> { + public: + LinearThermalIsotropic(vector l_e, json materialProperties) + : ThermalModel(l_e) + { + try { + conductivity = materialProperties["conductivity"].get>(); + } catch (const std::exception &e) { + throw std::runtime_error("Missing material properties for the requested material model."); + } + n_mat = conductivity.size(); + + double kappa_ref = (*max_element(conductivity.begin(), conductivity.end()) + + *min_element(conductivity.begin(), conductivity.end())) / + 2; + kapparef_mat = kappa_ref * Matrix3d::Identity(); + + Matrix3d phase_kappa; + phase_stiffness = new Matrix[n_mat]; + + for (size_t i = 0; i < n_mat; ++i) { + phase_stiffness[i] = Matrix::Zero(); + phase_kappa = conductivity[i] * Matrix3d::Identity(); + + for (int p = 0; p < 8; ++p) { + phase_stiffness[i] += B_int[p].transpose() * phase_kappa * B_int[p] * v_e * 0.1250; + } + } + } + + void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override + { + sigma(i + 0, 0) = conductivity[mat_index] * eps(i + 0, 0); + sigma(i + 1, 0) = conductivity[mat_index] * eps(i + 1, 0); + sigma(i + 2, 0) = conductivity[mat_index] * eps(i + 2, 0); + } + + private: + vector conductivity; +}; + +class LinearThermalTriclinic : public ThermalModel, public LinearModel<1> { + public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW // Ensure proper alignment for Eigen structures + + LinearThermalTriclinic(vector l_e, json materialProperties) + : ThermalModel(l_e) + { + vector K_keys = { + "K_11", "K_12", "K_13", + "K_22", "K_23", + "K_33"}; + + try { + n_mat = materialProperties.at("K_11").get>().size(); + size_t num_constants = K_keys.size(); + + // Initialize matrix to hold all constants + K_constants.resize(num_constants, n_mat); + + // Load material constants into matrix + for (size_t k = 0; k < num_constants; ++k) { + const auto &values = materialProperties.at(K_keys[k]).get>(); + if (values.size() != n_mat) { + throw std::runtime_error("Inconsistent size for material property: " + K_keys[k]); + } + K_constants.row(k) = Eigen::Map(values.data(), values.size()); + } + } catch (const std::exception &e) { + throw std::runtime_error("Missing or inconsistent material properties for the requested material model."); + } + + // Assemble conductivity matrices for each material + K_mats.resize(n_mat); + kapparef_mat = Matrix3d::Zero(); + + for (size_t i = 0; i < n_mat; ++i) { + Matrix3d K_i; + K_i << K_constants(0, i), K_constants(1, i), K_constants(2, i), + K_constants(1, i), K_constants(3, i), K_constants(4, i), + K_constants(2, i), K_constants(4, i), K_constants(5, i); + + K_mats[i] = K_i; + kapparef_mat += K_i; + } + + kapparef_mat /= n_mat; + + // Compute phase stiffness matrices + phase_stiffness = new Matrix[n_mat]; + for (size_t i = 0; i < n_mat; ++i) { + phase_stiffness[i] = Matrix::Zero(); + for (int p = 0; p < 8; ++p) { + phase_stiffness[i] += B_int[p].transpose() * K_mats[i] * B_int[p] * v_e * 0.1250; + } + } + } + + void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override + { + sigma.segment<3>(i) = K_mats[mat_index] * eps.segment<3>(i); + } + + private: + std::vector> K_mats; + MatrixXd K_constants; +}; + +#endif // LINEARTHERMAL_H diff --git a/include/material_models/LinearThermalIsotropic.h b/include/material_models/LinearThermalIsotropic.h deleted file mode 100644 index 0951608..0000000 --- a/include/material_models/LinearThermalIsotropic.h +++ /dev/null @@ -1,43 +0,0 @@ -#ifndef LINEARTHERMALISOTROPIC_H -#define LINEARTHERMALISOTROPIC_H - -#include "matmodel.h" - -class LinearThermalIsotropic : public ThermalModel, public LinearModel<1> { - public: - LinearThermalIsotropic(vector l_e, map> materialProperties) - : ThermalModel(l_e) - { - conductivity = materialProperties["conductivity"]; - n_mat = conductivity.size(); - - double kappa_ref = (*max_element(conductivity.begin(), conductivity.end()) + - *min_element(conductivity.begin(), conductivity.end())) / - 2; - kapparef_mat = kappa_ref * Matrix3d::Identity(); - - Matrix3d phase_kappa; - phase_stiffness = new Matrix[n_mat]; - - for (size_t i = 0; i < n_mat; ++i) { - phase_stiffness[i] = Matrix::Zero(); - phase_kappa = conductivity[i] * Matrix3d::Identity(); - - for (int p = 0; p < 8; ++p) { - phase_stiffness[i] += B_int[p].transpose() * phase_kappa * B_int[p] * v_e * 0.1250; - } - } - } - - void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override - { - sigma(i + 0, 0) = conductivity[mat_index] * eps(i + 0, 0); - sigma(i + 1, 0) = conductivity[mat_index] * eps(i + 1, 0); - sigma(i + 2, 0) = conductivity[mat_index] * eps(i + 2, 0); - } - - private: - vector conductivity; -}; - -#endif // LINEARTHERMALISOTROPIC_H diff --git a/include/material_models/PseudoPlastic.h b/include/material_models/PseudoPlastic.h index f3f8d21..b47f8ac 100644 --- a/include/material_models/PseudoPlastic.h +++ b/include/material_models/PseudoPlastic.h @@ -20,13 +20,13 @@ class PseudoPlastic : public MechModel { public: - PseudoPlastic(vector l_e, map> materialProperties) + PseudoPlastic(vector l_e, json materialProperties) : MechModel(l_e) { try { - bulk_modulus = materialProperties["bulk_modulus"]; - shear_modulus = materialProperties["shear_modulus"]; - yield_stress = materialProperties["yield_stress"]; + bulk_modulus = materialProperties["bulk_modulus"].get>(); + shear_modulus = materialProperties["shear_modulus"].get>(); + yield_stress = materialProperties["yield_stress"].get>(); } catch (const std::exception &e) { throw std::runtime_error("Missing material properties for the requested material model."); } @@ -84,11 +84,11 @@ class PseudoPlastic : public MechModel { class PseudoPlasticLinearHardening : public PseudoPlastic { public: - PseudoPlasticLinearHardening(vector l_e, map> materialProperties) + PseudoPlasticLinearHardening(vector l_e, json materialProperties) : PseudoPlastic(l_e, materialProperties) { try { - hardening_parameter = materialProperties["hardening_parameter"]; + hardening_parameter = materialProperties["hardening_parameter"].get>(); } catch (const std::exception &e) { throw std::runtime_error("Missing material properties for the requested material model."); } @@ -134,12 +134,12 @@ class PseudoPlasticLinearHardening : public PseudoPlastic { class PseudoPlasticNonLinearHardening : public PseudoPlastic { public: - PseudoPlasticNonLinearHardening(vector l_e, map> materialProperties) + PseudoPlasticNonLinearHardening(vector l_e, json materialProperties) : PseudoPlastic(l_e, materialProperties) { try { - hardening_exponent = materialProperties["hardening_exponent"]; - eps_0 = materialProperties["eps_0"]; // ε0 parameter + hardening_exponent = materialProperties["hardening_exponent"].get>(); + eps_0 = materialProperties["eps_0"].get>(); // ε0 parameter } catch (const std::exception &e) { throw std::runtime_error("Missing material properties for the requested material model."); } diff --git a/include/matmodel.h b/include/matmodel.h index f99519d..4bee51c 100644 --- a/include/matmodel.h +++ b/include/matmodel.h @@ -41,6 +41,8 @@ class Matmodel { virtual void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) {} virtual void updateInternalVariables() {} + vector macroscale_loading; + protected: double l_e_x; double l_e_y; @@ -152,6 +154,7 @@ void Matmodel::getStrainStress(double *strain, double *stress, Matrix void Matmodel::setGradient(vector _g0) { + macroscale_loading = _g0; for (int i = 0; i < n_str; i++) { for (int j = 0; j < 8; ++j) { g0(n_str * j + i, 0) = _g0[i]; diff --git a/include/reader.h b/include/reader.h index ca07ea0..9a0c7ba 100644 --- a/include/reader.h +++ b/include/reader.h @@ -14,8 +14,9 @@ class Reader { char ms_filename[4096]; // Name of Micro-structure hdf5 file char ms_datasetname[4096]; // Absolute path of Micro-structure in hdf5 file int n_mat; - map> materialProperties; + json materialProperties; double TOL; + json errorParameters; int n_it; vector>> g0; string problemType; diff --git a/include/setup.h b/include/setup.h index 79d3581..a9fec2f 100644 --- a/include/setup.h +++ b/include/setup.h @@ -2,10 +2,10 @@ #include "solverFP.h" // Thermal models -#include "material_models/LinearThermalIsotropic.h" +#include "material_models/LinearThermal.h" // Mechanical models -#include "material_models/LinearElasticIsotropic.h" +#include "material_models/LinearElastic.h" #include "material_models/PseudoPlastic.h" #include "material_models/J2Plasticity.h" @@ -17,6 +17,8 @@ Matmodel<1> *createMatmodel(const Reader &reader) { if (reader.matmodel == "LinearThermalIsotropic") { return new LinearThermalIsotropic(reader.l_e, reader.materialProperties); + } else if (reader.matmodel == "LinearThermalTriclinic") { + return new LinearThermalTriclinic(reader.l_e, reader.materialProperties); } else { throw std::invalid_argument(reader.matmodel + " is not a valid matmodel for thermal problem"); } @@ -28,6 +30,8 @@ Matmodel<3> *createMatmodel(const Reader &reader) // Linear Elastic models if (reader.matmodel == "LinearElasticIsotropic") { return new LinearElasticIsotropic(reader.l_e, reader.materialProperties); + } else if (reader.matmodel == "LinearElasticTriclinic") { + return new LinearElasticTriclinic(reader.l_e, reader.materialProperties); // Pseudo Plastic models } else if (reader.matmodel == "PseudoPlasticLinearHardening") { diff --git a/include/solver.h b/include/solver.h index 17531db..d1d2111 100644 --- a/include/solver.h +++ b/include/solver.h @@ -24,7 +24,7 @@ class Solver { const ptrdiff_t local_1_start; const int n_it; //!< Max number of FANS iterations - const double TOL; //!< Tolerance on relative error norm + double TOL; //!< Tolerance on relative error norm Matmodel *matmodel; //!< Material Model unsigned char *ms; // Micro-structure Binary @@ -43,7 +43,7 @@ class Solver { void iterateCubes(F f); void solve(); - virtual void internalSolve(){}; // important to have "{}" here, otherwise we get an error about undefined reference to vtable + virtual void internalSolve() {}; // important to have "{}" here, otherwise we get an error about undefined reference to vtable template void compute_residual_basic(RealArray &r_matrix, RealArray &u_matrix, F f); @@ -56,6 +56,12 @@ class Solver { double compute_error(RealArray &r); void CreateFFTWPlans(double *in, fftw_complex *transformed, double *out); + VectorXd homogenized_stress; + VectorXd get_homogenized_stress(); + + MatrixXd homogenized_tangent; + MatrixXd get_homogenized_tangent(double pert_param); + protected: fftw_plan planfft, planifft; clock_t fft_time, buftime; @@ -354,27 +360,41 @@ void Solver::convolution() template double Solver::compute_error(RealArray &r) { - double err, err0, err_local; + double err_local; + const std::string &measure = reader.errorParameters["measure"].get(); + if (measure == "L1") { + err_local = r.matrix().lpNorm<1>(); + } else if (measure == "L2") { + err_local = r.matrix().lpNorm<2>(); + } else if (measure == "Linfinity") { + err_local = r.matrix().lpNorm(); + } else { + throw std::runtime_error("Unknown measure type: " + measure); + } - err_local = r.matrix().lpNorm(); + double err; MPI_Allreduce(&err_local, &err, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); err_all[iter] = err; - err0 = err_all[0]; + double err0 = err_all[0]; double err_rel = (iter == 0 ? 100 : err / err0); if (world_rank == 0) { if (iter == 0) { printf("Before 1st iteration: %16.8e\n", err0); - } else if (iter == 1) { - printf("it %3lu .... err %16.8e / %8.4e, ratio: ------------- , FFT time: %2.6f sec\n", iter, err, err / err0, double(buftime) / CLOCKS_PER_SEC); } else { - printf("it %3lu .... err %16.8e / %8.4e, ratio: %4.8e, FFT time: %2.6f sec\n", iter, err, err / err0, err / err_all[iter - 1], double(buftime) / CLOCKS_PER_SEC); + printf("it %3lu .... err %16.8e / %8.4e, ratio: %4.8e, FFT time: %2.6f sec\n", iter, err, err / err0, (iter == 1 ? 0.0 : err / err_all[iter - 1]), double(buftime) / CLOCKS_PER_SEC); } } - return err; // returns absolute error - // return err_rel; // returns relative error + const std::string &error_type = reader.errorParameters["type"].get(); + if (error_type == "absolute") { + return err; + } else if (error_type == "relative") { + return err_rel; + } else { + throw std::runtime_error("Unknown error type: " + error_type); + } } template @@ -435,21 +455,20 @@ void Solver::postprocess(Reader reader, const char resultsFileName[], i if (world_rank == 0) { printf("# Effective Stress .. ("); for (int i = 0; i < n_str; ++i) - printf("%+f ", stress_average[i]); + printf("%+.12f ", stress_average[i]); printf(") \n"); printf("# Effective Strain .. ("); for (int i = 0; i < n_str; ++i) - printf("%+f ", strain_average[i]); + printf("%+.12f ", strain_average[i]); printf(") \n\n"); } // Write results to results h5 file - auto writeData = [&](const char *resultName, const char *resultPrefix, auto *data, hsize_t size) { + auto writeData = [&](const char *resultName, const char *resultPrefix, auto *data, hsize_t *dims, int ndims) { if (std::find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), resultName) != reader.resultsToWrite.end()) { char name[5096]; sprintf(name, "%s/load%i/time_step%i/%s", reader.ms_datasetname, load_idx, time_idx, resultPrefix); - hsize_t dims[1] = {size}; - reader.WriteData(data, resultsFileName, name, dims, 1); + reader.WriteData(data, resultsFileName, name, dims, ndims); } }; @@ -464,18 +483,20 @@ void Solver::postprocess(Reader reader, const char resultsFileName[], i for (int i = 0; i < world_size; ++i) { if (i == world_rank) { if (world_rank == 0) { - writeData("stress_average", "stress_average", stress_average.data(), n_str); - writeData("strain_average", "strain_average", strain_average.data(), n_str); - writeData("absolute_error", "absolute_error", err_all.data(), iter + 1); + hsize_t dims[1] = {static_cast(n_str)}; + writeData("stress_average", "stress_average", stress_average.data(), dims, 1); + writeData("strain_average", "strain_average", strain_average.data(), dims, 1); for (int mat_index = 0; mat_index < n_mat; ++mat_index) { char stress_name[512]; char strain_name[512]; sprintf(stress_name, "phase_stress_average_phase%d", mat_index); sprintf(strain_name, "phase_strain_average_phase%d", mat_index); - writeData("phase_stress_average", stress_name, phase_stress_average[mat_index].data(), n_str); - writeData("phase_strain_average", strain_name, phase_strain_average[mat_index].data(), n_str); + writeData("phase_stress_average", stress_name, phase_stress_average[mat_index].data(), dims, 1); + writeData("phase_strain_average", strain_name, phase_strain_average[mat_index].data(), dims, 1); } + dims[0] = iter + 1; + writeData("absolute_error", "absolute_error", err_all.data(), dims, 1); } writeSlab("microstructure", "microstructure", ms, 1); writeSlab("displacement", "displacement", v_u, howmany); @@ -486,6 +507,86 @@ void Solver::postprocess(Reader reader, const char resultsFileName[], i MPI_Barrier(MPI_COMM_WORLD); } matmodel->postprocess(*this, reader, resultsFileName, load_idx, time_idx); + + // Compute homogenized tangent + if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "homogenized_tangent") != reader.resultsToWrite.end()) { + homogenized_tangent = get_homogenized_tangent(1e-6); + hsize_t dims[2] = {static_cast(n_str), static_cast(n_str)}; + if (world_rank == 0) { + cout << "# Homogenized tangent: " << endl + << setprecision(12) << homogenized_tangent << endl + << endl; + writeData("homogenized_tangent", "homogenized_tangent", homogenized_tangent.data(), dims, 2); + } + } +} + +template +VectorXd Solver::get_homogenized_stress() +{ + + int n_str = matmodel->n_str; + VectorXd strain = VectorXd::Zero(local_n0 * n_y * n_z * n_str); + VectorXd stress = VectorXd::Zero(local_n0 * n_y * n_z * n_str); + homogenized_stress = VectorXd::Zero(n_str); + + MPI_Sendrecv(v_u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0, + v_u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); + + Matrix ue; + int mat_index; + iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *idxPadding) { + for (int i = 0; i < 8; ++i) { + for (int j = 0; j < howmany; ++j) { + ue(howmany * i + j, 0) = v_u[howmany * idx[i] + j]; + } + } + mat_index = ms[idx[0]]; + + matmodel->getStrainStress(strain.segment(n_str * idx[0], n_str).data(), stress.segment(n_str * idx[0], n_str).data(), ue, mat_index, idx[0]); + homogenized_stress += stress.segment(n_str * idx[0], n_str); + }); + + MPI_Allreduce(MPI_IN_PLACE, homogenized_stress.data(), n_str, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + homogenized_stress /= (n_x * n_y * n_z); + + return homogenized_stress; +} + +template +MatrixXd Solver::get_homogenized_tangent(double pert_param) +{ + int n_str = matmodel->n_str; + homogenized_tangent = MatrixXd::Zero(n_str, n_str); + VectorXd unperturbed_stress = get_homogenized_stress(); + VectorXd perturbed_stress; + vector pert_strain(n_str, 0.0); + vector g0 = this->matmodel->macroscale_loading; + bool islinear = dynamic_cast *>(this->matmodel) != nullptr; + + this->reader.errorParameters["type"] = "relative"; + this->TOL = max(1e-6, this->TOL); + + // TODO: a deep copy of the solver object is needed here to avoid modifying the history of the solver object + + for (int i = 0; i < n_str; i++) { + if (islinear) { + pert_strain = vector(n_str, 0.0); + pert_strain[i] = 1.0; + } else { + pert_strain = g0; + pert_strain[i] += pert_param; + } + + matmodel->setGradient(pert_strain); + solve(); + perturbed_stress = get_homogenized_stress(); + + homogenized_tangent.col(i) = islinear ? perturbed_stress : (perturbed_stress - unperturbed_stress) / pert_param; + } + + homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval(); + return homogenized_tangent; } #endif diff --git a/src/reader.cpp b/src/reader.cpp index 5a2e2b1..973bc4f 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -10,16 +10,28 @@ #include "H5FDmpi.h" #include "H5FDmpio.h" -using namespace std; - -#include -using nlohmann::json; -using namespace nlohmann; - void Reader::ComputeVolumeFractions() { - if (world_rank == 0) + // TODO: this is not the best way to compute the number of materials + // Determine n_mat as the maximum material index plus one + unsigned char local_max = 0; + size_t local_size = local_n0 * dims[1] * dims[2]; + + // Find the local maximum material index + for (size_t i = 0; i < local_size; i++) { + if (ms[i] > local_max) { + local_max = ms[i]; + } + } + // Find the global maximum material index + unsigned char global_max; + MPI_Allreduce(&local_max, &global_max, 1, MPI_UNSIGNED_CHAR, MPI_MAX, MPI_COMM_WORLD); + n_mat = global_max + 1; // Set n_mat to the maximum material index plus one + + if (world_rank == 0) { + printf("# Number of materials: %i\n", n_mat); printf("# Volume fractions\n"); + } long vol_frac[n_mat]; double v_frac[n_mat]; for (int i = 0; i < n_mat; i++) { @@ -53,9 +65,10 @@ void Reader ::ReadInputFile(char fn[]) L = j["ms_L"].get>(); - TOL = j["TOL"].get(); - n_it = j["n_it"].get(); - g0 = j["macroscale_loading"].get>>>(); + errorParameters = j["error_parameters"]; + TOL = errorParameters["tolerance"].get(); + n_it = j["n_it"].get(); + g0 = j["macroscale_loading"].get>>>(); problemType = j["problem_type"].get(); matmodel = j["matmodel"].get(); @@ -67,17 +80,35 @@ void Reader ::ReadInputFile(char fn[]) if (world_rank == 0) { printf("# microstructure file name: \t '%s'\n", ms_filename); printf("# microstructure dataset name: \t '%s'\n", ms_datasetname); - printf("# FANS Tolerance: \t %10.5e\n# Max iterations: \t %6i\n", TOL, n_it); + printf( + "# FANS error measure: \t %s %s error \n", + errorParameters["type"].get().c_str(), + errorParameters["measure"].get().c_str()); + printf("# FANS Tolerance: \t %10.5e\n", errorParameters["tolerance"].get()); + printf("# Max iterations: \t %6i\n", n_it); } for (auto it = j_mat.begin(); it != j_mat.end(); ++it) { - materialProperties[it.key()] = it.value().get>(); - n_mat = materialProperties[it.key()].size(); + materialProperties[it.key()] = it.value(); if (world_rank == 0) { cout << "# " << it.key() << ":\t "; - for (double d : materialProperties[it.key()]) { - printf(" %10.5f", d); + if (it.value().is_array()) { + for (const auto &elem : it.value()) { + if (elem.is_number()) { + printf(" %10.5f", elem.get()); + } else if (elem.is_string()) { + cout << " " << elem.get(); + } else { + cout << " " << elem; + } + } + } else if (it.value().is_number()) { + printf(" %10.5f", it.value().get()); + } else if (it.value().is_string()) { + cout << " " << it.value().get(); + } else { + cout << " " << it.value(); } printf("\n"); } diff --git a/test/input_files/test_J2Plasticity.json b/test/input_files/test_J2Plasticity.json index ef770ee..8b5e46c 100644 --- a/test/input_files/test_J2Plasticity.json +++ b/test/input_files/test_J2Plasticity.json @@ -12,14 +12,18 @@ "isotropic_hardening_parameter": [0.0, 0.0], "kinematic_hardening_parameter": [0.0, 0.0], "viscosity": [1, 1], - "time_step": [0.01], + "time_step": 0.01, "saturation_stress": [0.15, 10000], "saturation_exponent": [1000, 1000] }, "method": "cg", - "TOL": 1e-10, + "error_parameters":{ + "measure": "Linfinity", + "type": "absolute", + "tolerance": 1e-10 + }, "n_it": 100, "macroscale_loading": [ [ [0.0000, 0, 0, 0, 0, 0], [0.0001, 0, 0, 0, 0, 0], diff --git a/test/input_files/test_LinearElastic.json b/test/input_files/test_LinearElastic.json index 61f9f7a..6f7f448 100644 --- a/test/input_files/test_LinearElastic.json +++ b/test/input_files/test_LinearElastic.json @@ -11,7 +11,11 @@ }, "method": "cg", - "TOL": 1e-10, + "error_parameters":{ + "measure": "Linfinity", + "type": "absolute", + "tolerance": 1e-10 + }, "n_it": 100, "macroscale_loading": [ [[1, 0, 0, 0, 0, 0]], diff --git a/test/input_files/test_LinearThermal.json b/test/input_files/test_LinearThermal.json index d1c069f..91ecb65 100644 --- a/test/input_files/test_LinearThermal.json +++ b/test/input_files/test_LinearThermal.json @@ -10,7 +10,11 @@ }, "method": "cg", - "TOL": 1e-10, + "error_parameters":{ + "measure": "Linfinity", + "type": "absolute", + "tolerance": 1e-10 + }, "n_it": 100, "macroscale_loading": [ [[1, 0, 0]], diff --git a/test/input_files/test_PseudoPlastic.json b/test/input_files/test_PseudoPlastic.json index 12deb08..713a15e 100644 --- a/test/input_files/test_PseudoPlastic.json +++ b/test/input_files/test_PseudoPlastic.json @@ -15,7 +15,11 @@ }, "method": "cg", - "TOL": 1e-10, + "error_parameters":{ + "measure": "Linfinity", + "type": "absolute", + "tolerance": 1e-10 + }, "n_it": 100, "macroscale_loading": [ [