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
@@ -177,6 +177,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,16 +41,15 @@ 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
 
   protected:
     double l_e_x;
     double l_e_y;
     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);
 
@@ -62,6 +62,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;
@@ -583,6 +600,7 @@ MatrixXd Solver<howmany>::get_homogenized_tangent(double pert_param)
         }
 
         matmodel->setGradient(pert_strain);
+        disableMixedBC();
         solve();
         perturbed_stress = get_homogenized_stress();
 

diff --git a/include/solverCG.h b/include/solverCG.h
@@ -93,6 +93,11 @@ void SolverCG<howmany>::internalSolve()
 
         if (islinear) {
             Matrix<double, howmany * 8, 1> res_e;
+            if (this->isMixedBCActive()) {
+                // if it is mixed and linear, we need to update the residual
+                this->updateMixedBC();
+                this->template compute_residual<2>(v_r_real, v_u_real);
+            }
             this->template compute_residual_basic<0>(rnew_real, d_real,
                                                      [&](Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx) -> Matrix<double, howmany * 8, 1> & {
                                                          res_e.noalias() = linearModel->phase_stiffness[mat_index] * ue;
@@ -129,6 +134,7 @@ void SolverCG<howmany>::LineSearchSecant()
     while (((_iter < MaxIter) && (err > tol))) {
 
         v_u_real += d_real * (alpha_new - alpha_old);
+        this->updateMixedBC();
         this->template compute_residual<0>(rnew_real, v_u_real);
         r1pd = dotProduct(rnew_real, d_real);
 

diff --git a/include/solverFP.h b/include/solverFP.h
@@ -44,6 +44,7 @@ void SolverFP<howmany>::internalSolve()
 
         this->convolution();
         v_u_real -= v_r_real;
+        this->updateMixedBC();
         this->template compute_residual<2>(v_r_real, v_u_real);
 
         iter++;

diff --git a/src/main.cpp b/src/main.cpp
@@ -11,23 +11,20 @@ void runSolver(Reader &reader, const char *output_file_basename)
 {
     reader.ReadMS(howmany);
 
-    for (size_t load_path_idx = 0; load_path_idx < reader.g0.size(); ++load_path_idx) {
+    for (size_t load_path_idx = 0; load_path_idx < reader.load_cases.size(); ++load_path_idx) {
         Matmodel<howmany> *matmodel = createMatmodel<howmany>(reader);
         Solver<howmany>   *solver   = createSolver(reader, matmodel);
 
-        const auto &load_path = reader.g0[load_path_idx];
-        for (size_t time_step_idx = 0; time_step_idx < load_path.size(); ++time_step_idx) {
-            const auto &g0 = load_path[time_step_idx];
-
-            if (g0.size() != matmodel->n_str) {
-                throw std::invalid_argument("Invalid length of loading g0");
+        for (size_t time_step_idx = 0; time_step_idx < reader.load_cases[load_path_idx].n_steps; ++time_step_idx) {
+            if (reader.load_cases[load_path_idx].mixed) {
+                solver->enableMixedBC(reader.load_cases[load_path_idx].mbc, time_step_idx);
+            } else {
+                const auto &g0 = reader.load_cases[load_path_idx].g0_path[time_step_idx];
+                matmodel->setGradient(g0);
             }
-
-            matmodel->setGradient(g0);
             solver->solve();
             solver->postprocess(reader, output_file_basename, load_path_idx, time_step_idx);
         }
-
         delete solver;
         delete matmodel;
     }