diff --git a/CMakeLists.txt b/CMakeLists.txt index 1918e68..d80ea8a 100644 --- a/CMakeLists.txt +++ b/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 ) 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/include/general.h b/include/general.h index 70d9fc4..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 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/LinearThermalIsotropic.h b/include/material_models/LinearThermalIsotropic.h index 0951608..35446b8 100644 --- a/include/material_models/LinearThermalIsotropic.h +++ b/include/material_models/LinearThermalIsotropic.h @@ -5,11 +5,15 @@ class LinearThermalIsotropic : public ThermalModel, public LinearModel<1> { public: - LinearThermalIsotropic(vector l_e, map> materialProperties) + LinearThermalIsotropic(vector l_e, json materialProperties) : ThermalModel(l_e) { - conductivity = materialProperties["conductivity"]; - n_mat = conductivity.size(); + 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())) / 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..beecd1c 100644 --- a/include/setup.h +++ b/include/setup.h @@ -5,7 +5,7 @@ #include "material_models/LinearThermalIsotropic.h" // Mechanical models -#include "material_models/LinearElasticIsotropic.h" +#include "material_models/LinearElastic.h" #include "material_models/PseudoPlastic.h" #include "material_models/J2Plasticity.h" @@ -28,6 +28,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 f2ccc1e..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 @@ -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": [ [