DataAnalyticsEngineering · sanathkeshav · May 20, 2025 · May 19, 2025 · May 19, 2025 · May 19, 2025
diff --git a/.gitignore b/.gitignore
@@ -205,6 +205,7 @@ test/input_files/**/*.json
 !test_LinearThermal.json
 !test_PseudoPlastic.json
 !test_J2Plasticity.json
+!test_MixedBCs.json
 
 # pixi environments
 .pixi

diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## latest
 
+- Add support for macroscale mixed stress-strain boundary conditions https://github.com/DataAnalyticsEngineering/FANS/pull/58
 - Add grain boundary diffusion material model for polycrystals https://github.com/DataAnalyticsEngineering/FANS/pull/52
 - Add a pixi environment for the FANS_dashboard and some tests https://github.com/DataAnalyticsEngineering/FANS/pull/55
 - Remove MPI initialization from pyFANS and add an integration test for it https://github.com/DataAnalyticsEngineering/FANS/pull/46

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -148,6 +148,7 @@ set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/solverFP.h
         include/solver.h
         include/setup.h
+        include/mixedBCs.h
 
         include/material_models/LinearThermal.h
         include/material_models/GBDiffusion.h

diff --git a/include/matmodel.h b/include/matmodel.h
@@ -41,7 +41,8 @@ class Matmodel {
     virtual void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) {}
     virtual void updateInternalVariables() {}
 
-    vector<double> macroscale_loading;
+    vector<double>               macroscale_loading;
+    Matrix<double, n_str, n_str> kapparef_mat; // Reference conductivity matrix
 
     virtual ~Matmodel() = default;
 
@@ -51,8 +52,6 @@ class Matmodel {
     double l_e_z;
     double v_e;
 
-    Matrix<double, n_str, n_str> kapparef_mat; // Reference conductivity matrix
-
     Matrix<double, n_str, howmany * 8>     B_el_mean; //!< precomputed mean B matrix over the element
     Matrix<double, n_str, howmany * 8>     B_int[8];  //!< precomputed B matrix at all integration points
     Matrix<double, 8 * n_str, howmany * 8> B;

diff --git a/include/mixedBCs.h b/include/mixedBCs.h
@@ -0,0 +1,202 @@
+#ifndef MIXED_BC_H
+#define MIXED_BC_H
+
+// ============================================================================
+//  mixedBCs.h
+//  --------------------------------------------------------------------------
+//  • Holds MixedBC + LoadCase structs
+//  • Provides:   MixedBC::from_json(), finalize()
+// ============================================================================
+
+#include <Eigen/Dense>
+using namespace Eigen;
+#include <json.hpp>
+using nlohmann::json;
+
+// ---------------------------------------------------------------------------
+struct MixedBC {
+    /*  Index sets (0‑based)  */
+    VectorXi idx_E; // strain‑controlled components
+    VectorXi idx_F; // stress‑controlled components
+
+    /*  Time paths  */
+    MatrixXd F_E_path; // (#steps × |E|)
+    MatrixXd P_F_path; // (#steps × |F|)
+
+    /*  Projectors & auxiliary matrix  */
+    MatrixXd Q_E, Q_F, M; // M = (Q_Fᵀ C0 Q_F)⁻¹
+
+    // ------------------------------------------------------------
+    // build Q_E, Q_F, M
+    // ------------------------------------------------------------
+    void finalize(const MatrixXd &C0)
+    {
+        const int n_str = static_cast<int>(C0.rows());
+
+        Q_E = MatrixXd::Zero(n_str, idx_E.size());
+        for (int c = 0; c < idx_E.size(); ++c)
+            Q_E(idx_E(c), c) = 1.0;
+
+        Q_F = MatrixXd::Zero(n_str, idx_F.size());
+        for (int c = 0; c < idx_F.size(); ++c)
+            Q_F(idx_F(c), c) = 1.0;
+
+        if (idx_F.size() > 0)
+            M = (Q_F.transpose() * C0 * Q_F).inverse();
+        else
+            M.resize(0, 0); // pure‑strain case
+    }
+
+    // ------------------------------------------------------------
+    // Factory: parse MixedBC from a JSON object
+    // ------------------------------------------------------------
+    static MixedBC from_json(const json &jc, int n_str)
+    {
+        auto toEigen = [](const vector<int> &v) {
+            VectorXi e(v.size());
+            for (size_t i = 0; i < v.size(); ++i)
+                e(static_cast<int>(i)) = v[i];
+            return e;
+        };
+
+        MixedBC bc;
+
+        if (!jc.contains("strain_indices") || !jc.contains("stress_indices"))
+            throw runtime_error("mixed BC: strain_indices or stress_indices missing");
+
+        bc.idx_E = toEigen(jc["strain_indices"].get<vector<int>>());
+        bc.idx_F = toEigen(jc["stress_indices"].get<vector<int>>());
+
+        // ---- sanity: disjoint + complementary -----
+        vector<char> present(n_str, 0);
+        for (int i = 0; i < bc.idx_E.size(); ++i) {
+            int k = bc.idx_E(i);
+            if (k < 0 || k >= n_str)
+                throw runtime_error("strain index out of range");
+            present[k] = 1;
+        }
+        for (int i = 0; i < bc.idx_F.size(); ++i) {
+            int k = bc.idx_F(i);
+            if (k < 0 || k >= n_str)
+                throw runtime_error("stress index out of range");
+            if (present[k])
+                throw runtime_error("index appears in both strain_indices and stress_indices");
+            present[k] = 1;
+        }
+        for (int k = 0; k < n_str; ++k)
+            if (!present[k])
+                throw runtime_error("each component must be either strain‑ or stress‑controlled");
+
+        // ---- parse 2‑D arrays (allow empty for |E|==0 etc.) -----
+        auto get2D = [](const json &arr) {
+            return arr.get<vector<vector<double>>>();
+        };
+
+        vector<vector<double>> strain_raw, stress_raw;
+        size_t                 n_steps = 0;
+
+        if (bc.idx_E.size() > 0) {
+            if (!jc.contains("strain"))
+                throw runtime_error("strain array missing");
+            strain_raw = get2D(jc["strain"]);
+            n_steps    = strain_raw.size();
+        }
+        if (bc.idx_F.size() > 0) {
+            if (!jc.contains("stress"))
+                throw runtime_error("stress array missing");
+            stress_raw = get2D(jc["stress"]);
+            n_steps    = max(n_steps, stress_raw.size());
+        }
+        if (n_steps == 0)
+            throw runtime_error("mixed BC: at least one of strain/stress must have timesteps");
+
+        // default‑fill for missing part (constant 0)
+        if (strain_raw.empty())
+            strain_raw.resize(n_steps, vector<double>(0));
+        if (stress_raw.empty())
+            stress_raw.resize(n_steps, vector<double>(0));
+
+        // consistency checks & build Eigen matrices
+        bc.F_E_path = MatrixXd::Zero(n_steps, bc.idx_E.size());
+        bc.P_F_path = MatrixXd::Zero(n_steps, bc.idx_F.size());
+
+        for (size_t t = 0; t < n_steps; ++t) {
+            if (strain_raw[t].size() != static_cast<size_t>(bc.idx_E.size()))
+                throw runtime_error("strain row length mismatch");
+            for (int c = 0; c < bc.idx_E.size(); ++c)
+                bc.F_E_path(static_cast<int>(t), c) = (bc.idx_E.size() ? strain_raw[t][c] : 0.0);
+
+            if (stress_raw[t].size() != static_cast<size_t>(bc.idx_F.size()))
+                throw runtime_error("stress row length mismatch");
+            for (int c = 0; c < bc.idx_F.size(); ++c)
+                bc.P_F_path(static_cast<int>(t), c) = (bc.idx_F.size() ? stress_raw[t][c] : 0.0);
+        }
+        return bc;
+    }
+};
+
+// ---------------------------------------------------------------------------
+//  LoadCase : covers both legacy and mixed formats (used by Reader)
+// ---------------------------------------------------------------------------
+struct LoadCase {
+    bool                   mixed = false;
+    vector<vector<double>> g0_path;     // legacy pure‑strain
+    MixedBC                mbc;         // mixed BC data
+    size_t                 n_steps = 0; // number of time steps
+};
+
+// ---------------------------------------------------------------------------
+//  Helper mix‑in with enable/update for Solver (header‑only)
+// ---------------------------------------------------------------------------
+template <int howmany>
+struct MixedBCController {
+    bool mixed_active = false;
+
+  protected:
+    MixedBC  mbc_local;
+    size_t   step_idx = 0;
+    VectorXd g0_vec; // current macro strain (size n_str)
+
+    // call from user code after v_u update each iteration
+    template <typename SolverType>
+    void update(SolverType &solver)
+    {
+        if (!mixed_active)
+            return;
+
+        VectorXd Pbar = solver.get_homogenized_stress();
+
+        VectorXd PF = (mbc_local.idx_F.size() ? mbc_local.P_F_path.row(step_idx).transpose() : VectorXd());
+
+        if (mbc_local.idx_F.size()) {
+            VectorXd rhs     = PF - mbc_local.Q_F.transpose() * Pbar;
+            VectorXd delta_E = mbc_local.M * rhs;
+            g0_vec += mbc_local.Q_F * delta_E;
+        }
+
+        vector<double> gvec(g0_vec.data(), g0_vec.data() + g0_vec.size());
+        solver.matmodel->setGradient(gvec);
+    }
+
+    template <typename SolverType>
+    void activate(SolverType &solver, const MixedBC &mbc_in, size_t t)
+    {
+        mixed_active = true;
+        mbc_local    = mbc_in;
+        step_idx     = t;
+        mbc_local.finalize(solver.matmodel->kapparef_mat);
+
+        int n_str = solver.matmodel->n_str;
+        g0_vec    = VectorXd::Zero(n_str);
+        if (mbc_local.idx_E.size())
+            g0_vec += mbc_local.Q_E * mbc_local.F_E_path.row(t).transpose();
+
+        vector<double> gvec(g0_vec.data(), g0_vec.data() + n_str);
+        solver.matmodel->setGradient(gvec);
+
+        // Do one update to set the initial strain
+        solver.updateMixedBC();
+    }
+};
+
+#endif // MIXED_BC_H
diff --git a/include/reader.h b/include/reader.h
@@ -5,25 +5,26 @@
 #include <map>
 #include <string>
 #include <vector>
+#include "mixedBCs.h"
 
 using namespace std;
 
 class Reader {
   public:
     // contents of input file:
-    char                           ms_filename[4096];    // Name of Micro-structure hdf5 file
-    char                           ms_datasetname[4096]; // Absolute path of Micro-structure in hdf5 file
-    char                           results_prefix[4096];
-    int                            n_mat;
-    json                           materialProperties;
-    double                         TOL;
-    json                           errorParameters;
-    json                           microstructure;
-    int                            n_it;
-    vector<vector<vector<double>>> g0;
-    string                         problemType;
-    string                         matmodel;
-    string                         method;
+    char             ms_filename[4096];    // Name of Micro-structure hdf5 file
+    char             ms_datasetname[4096]; // Absolute path of Micro-structure in hdf5 file
+    char             results_prefix[4096];
+    int              n_mat;
+    json             materialProperties;
+    double           TOL;
+    json             errorParameters;
+    json             microstructure;
+    int              n_it;
+    vector<LoadCase> load_cases;
+    string           problemType;
+    string           matmodel;
+    string           method;
 
     vector<string> resultsToWrite;
 

diff --git a/include/solver.h b/include/solver.h
@@ -6,7 +6,7 @@
 typedef Map<Array<double, Dynamic, Dynamic>, Unaligned, OuterStride<>> RealArray;
 
 template <int howmany>
-class Solver {
+class Solver : private MixedBCController<howmany> {
   public:
     Solver(Reader reader, Matmodel<howmany> *matmodel);
     virtual ~Solver() = default;
@@ -63,6 +63,23 @@ class Solver {
     MatrixXd homogenized_tangent;
     MatrixXd get_homogenized_tangent(double pert_param);
 
+    void enableMixedBC(const MixedBC &mbc, size_t step)
+    {
+        this->activate(*this, mbc, step);
+    }
+    void disableMixedBC()
+    {
+        this->mixed_active = false;
+    }
+    bool isMixedBCActive()
+    {
+        return this->mixed_active;
+    }
+    void updateMixedBC()
+    {
+        this->update(*this);
+    }
+
   protected:
     fftw_plan planfft, planifft;
     clock_t   fft_time, buftime;
@@ -464,6 +481,49 @@ void Solver<howmany>::postprocess(Reader reader, const char resultsFileName[], i
         printf(") \n\n");
     }
 
+    /* ====================================================================== *
+     *  u_total = g0·X  +  ũ          (vector or scalar, decided at compile time)
+     * ====================================================================== */
+    const double     dx  = reader.l_e[0];
+    const double     dy  = reader.l_e[1];
+    const double     dz  = reader.l_e[2];
+    const double     Lx2 = reader.L[0] / 2.0;
+    const double     Ly2 = reader.L[1] / 2.0;
+    const double     Lz2 = reader.L[2] / 2.0;
+    constexpr double rs2 = 1.0 / std::sqrt(2.0);
+    VectorXd         u_total(local_n0 * n_y * n_z * howmany);
+    /* ---------- single sweep ------------------------------------------------- */
+    ptrdiff_t n = 0;
+    for (ptrdiff_t ix = 0; ix < local_n0; ++ix) {
+        const double x = (local_0_start + ix) * dx - Lx2;
+        for (ptrdiff_t iy = 0; iy < n_y; ++iy) {
+            const double y = iy * dy - Ly2;
+            for (ptrdiff_t iz = 0; iz < n_z; ++iz, ++n) {
+                const double z = iz * dz - Lz2;
+                if (howmany == 3) { /* ===== mechanics (vector) ===== */
+                    const double    g11 = strain_average[0];
+                    const double    g22 = strain_average[1];
+                    const double    g33 = strain_average[2];
+                    const double    g12 = strain_average[3] * rs2;
+                    const double    g13 = strain_average[4] * rs2;
+                    const double    g23 = strain_average[5] * rs2;
+                    const double    ux  = g11 * x + g12 * y + g13 * z;
+                    const double    uy  = g12 * x + g22 * y + g23 * z;
+                    const double    uz  = g13 * x + g23 * y + g33 * z;
+                    const ptrdiff_t b   = 3 * n;
+                    u_total[b]          = v_u[b] + ux;
+                    u_total[b + 1]      = v_u[b + 1] + uy;
+                    u_total[b + 2]      = v_u[b + 2] + uz;
+                } else { /* ===== scalar (howmany==1) ==== */
+                    const double g1 = strain_average[0];
+                    const double g2 = strain_average[1];
+                    const double g3 = strain_average[2];
+                    u_total[n]      = v_u[n] + (g1 * x + g2 * y + g3 * z);
+                }
+            }
+        }
+    }
+
     // Concatenate reader.ms_datasetname and reader.results_prefix into reader.ms_datasetname
     strcat(reader.ms_datasetname, "_results/");
     strcat(reader.ms_datasetname, reader.results_prefix);
@@ -504,7 +564,8 @@ void Solver<howmany>::postprocess(Reader reader, const char resultsFileName[], i
                 writeData("absolute_error", "absolute_error", err_all.data(), dims, 1);
             }
             writeSlab("microstructure", "microstructure", ms, 1);
-            writeSlab("displacement", "displacement", v_u, howmany);
+            writeSlab("displacement_fluctuation", "displacement_fluctuation", v_u, howmany);
+            writeSlab("displacement", "displacement", u_total.data(), howmany);
             writeSlab("residual", "residual", v_r, howmany);
             writeSlab("strain", "strain", strain.data(), n_str);
             writeSlab("stress", "stress", stress.data(), n_str);
@@ -584,6 +645,7 @@ MatrixXd Solver<howmany>::get_homogenized_tangent(double pert_param)
         }
 
         matmodel->setGradient(pert_strain);
+        disableMixedBC();
         solve();
         perturbed_stress = get_homogenized_stress();