|
| 1 | +#ifndef MIXED_BC_H |
| 2 | +#define MIXED_BC_H |
| 3 | + |
| 4 | +// ============================================================================ |
| 5 | +// mixedBCs.h |
| 6 | +// -------------------------------------------------------------------------- |
| 7 | +// • Holds MixedBC + LoadCase structs |
| 8 | +// • Provides: MixedBC::from_json(), finalize() |
| 9 | +// ============================================================================ |
| 10 | + |
| 11 | +#include <Eigen/Dense> |
| 12 | +using namespace Eigen; |
| 13 | +#include <json.hpp> |
| 14 | +using nlohmann::json; |
| 15 | + |
| 16 | +// --------------------------------------------------------------------------- |
| 17 | +struct MixedBC { |
| 18 | + /* Index sets (0‑based) */ |
| 19 | + VectorXi idx_E; // strain‑controlled components |
| 20 | + VectorXi idx_F; // stress‑controlled components |
| 21 | + |
| 22 | + /* Time paths */ |
| 23 | + MatrixXd F_E_path; // (#steps × |E|) |
| 24 | + MatrixXd P_F_path; // (#steps × |F|) |
| 25 | + |
| 26 | + /* Projectors & auxiliary matrix */ |
| 27 | + MatrixXd Q_E, Q_F, M; // M = (Q_Fᵀ C0 Q_F)⁻¹ |
| 28 | + |
| 29 | + // ------------------------------------------------------------ |
| 30 | + // build Q_E, Q_F, M |
| 31 | + // ------------------------------------------------------------ |
| 32 | + void finalize(const MatrixXd &C0) |
| 33 | + { |
| 34 | + const int n_str = static_cast<int>(C0.rows()); |
| 35 | + |
| 36 | + Q_E = MatrixXd::Zero(n_str, idx_E.size()); |
| 37 | + for (int c = 0; c < idx_E.size(); ++c) |
| 38 | + Q_E(idx_E(c), c) = 1.0; |
| 39 | + |
| 40 | + Q_F = MatrixXd::Zero(n_str, idx_F.size()); |
| 41 | + for (int c = 0; c < idx_F.size(); ++c) |
| 42 | + Q_F(idx_F(c), c) = 1.0; |
| 43 | + |
| 44 | + if (idx_F.size() > 0) |
| 45 | + M = (Q_F.transpose() * C0 * Q_F).inverse(); |
| 46 | + else |
| 47 | + M.resize(0, 0); // pure‑strain case |
| 48 | + } |
| 49 | + |
| 50 | + // ------------------------------------------------------------ |
| 51 | + // Factory: parse MixedBC from a JSON object |
| 52 | + // ------------------------------------------------------------ |
| 53 | + static MixedBC from_json(const json &jc, int n_str) |
| 54 | + { |
| 55 | + auto toEigen = [](const vector<int> &v) { |
| 56 | + VectorXi e(v.size()); |
| 57 | + for (size_t i = 0; i < v.size(); ++i) |
| 58 | + e(static_cast<int>(i)) = v[i]; |
| 59 | + return e; |
| 60 | + }; |
| 61 | + |
| 62 | + MixedBC bc; |
| 63 | + |
| 64 | + if (!jc.contains("strain_indices") || !jc.contains("stress_indices")) |
| 65 | + throw runtime_error("mixed BC: strain_indices or stress_indices missing"); |
| 66 | + |
| 67 | + bc.idx_E = toEigen(jc["strain_indices"].get<vector<int>>()); |
| 68 | + bc.idx_F = toEigen(jc["stress_indices"].get<vector<int>>()); |
| 69 | + |
| 70 | + // ---- sanity: disjoint + complementary ----- |
| 71 | + vector<char> present(n_str, 0); |
| 72 | + for (int i = 0; i < bc.idx_E.size(); ++i) { |
| 73 | + int k = bc.idx_E(i); |
| 74 | + if (k < 0 || k >= n_str) |
| 75 | + throw runtime_error("strain index out of range"); |
| 76 | + present[k] = 1; |
| 77 | + } |
| 78 | + for (int i = 0; i < bc.idx_F.size(); ++i) { |
| 79 | + int k = bc.idx_F(i); |
| 80 | + if (k < 0 || k >= n_str) |
| 81 | + throw runtime_error("stress index out of range"); |
| 82 | + if (present[k]) |
| 83 | + throw runtime_error("index appears in both strain_indices and stress_indices"); |
| 84 | + present[k] = 1; |
| 85 | + } |
| 86 | + for (int k = 0; k < n_str; ++k) |
| 87 | + if (!present[k]) |
| 88 | + throw runtime_error("each component must be either strain‑ or stress‑controlled"); |
| 89 | + |
| 90 | + // ---- parse 2‑D arrays (allow empty for |E|==0 etc.) ----- |
| 91 | + auto get2D = [](const json &arr) { |
| 92 | + return arr.get<vector<vector<double>>>(); |
| 93 | + }; |
| 94 | + |
| 95 | + vector<vector<double>> strain_raw, stress_raw; |
| 96 | + size_t n_steps = 0; |
| 97 | + |
| 98 | + if (bc.idx_E.size() > 0) { |
| 99 | + if (!jc.contains("strain")) |
| 100 | + throw runtime_error("strain array missing"); |
| 101 | + strain_raw = get2D(jc["strain"]); |
| 102 | + n_steps = strain_raw.size(); |
| 103 | + } |
| 104 | + if (bc.idx_F.size() > 0) { |
| 105 | + if (!jc.contains("stress")) |
| 106 | + throw runtime_error("stress array missing"); |
| 107 | + stress_raw = get2D(jc["stress"]); |
| 108 | + n_steps = max(n_steps, stress_raw.size()); |
| 109 | + } |
| 110 | + if (n_steps == 0) |
| 111 | + throw runtime_error("mixed BC: at least one of strain/stress must have timesteps"); |
| 112 | + |
| 113 | + // default‑fill for missing part (constant 0) |
| 114 | + if (strain_raw.empty()) |
| 115 | + strain_raw.resize(n_steps, vector<double>(0)); |
| 116 | + if (stress_raw.empty()) |
| 117 | + stress_raw.resize(n_steps, vector<double>(0)); |
| 118 | + |
| 119 | + // consistency checks & build Eigen matrices |
| 120 | + bc.F_E_path = MatrixXd::Zero(n_steps, bc.idx_E.size()); |
| 121 | + bc.P_F_path = MatrixXd::Zero(n_steps, bc.idx_F.size()); |
| 122 | + |
| 123 | + for (size_t t = 0; t < n_steps; ++t) { |
| 124 | + if (strain_raw[t].size() != static_cast<size_t>(bc.idx_E.size())) |
| 125 | + throw runtime_error("strain row length mismatch"); |
| 126 | + for (int c = 0; c < bc.idx_E.size(); ++c) |
| 127 | + bc.F_E_path(static_cast<int>(t), c) = (bc.idx_E.size() ? strain_raw[t][c] : 0.0); |
| 128 | + |
| 129 | + if (stress_raw[t].size() != static_cast<size_t>(bc.idx_F.size())) |
| 130 | + throw runtime_error("stress row length mismatch"); |
| 131 | + for (int c = 0; c < bc.idx_F.size(); ++c) |
| 132 | + bc.P_F_path(static_cast<int>(t), c) = (bc.idx_F.size() ? stress_raw[t][c] : 0.0); |
| 133 | + } |
| 134 | + return bc; |
| 135 | + } |
| 136 | +}; |
| 137 | + |
| 138 | +// --------------------------------------------------------------------------- |
| 139 | +// LoadCase : covers both legacy and mixed formats (used by Reader) |
| 140 | +// --------------------------------------------------------------------------- |
| 141 | +struct LoadCase { |
| 142 | + bool mixed = false; |
| 143 | + vector<vector<double>> g0_path; // legacy pure‑strain |
| 144 | + MixedBC mbc; // mixed BC data |
| 145 | + size_t n_steps = 0; // number of time steps |
| 146 | +}; |
| 147 | + |
| 148 | +// --------------------------------------------------------------------------- |
| 149 | +// Helper mix‑in with enable/update for Solver (header‑only) |
| 150 | +// --------------------------------------------------------------------------- |
| 151 | +template <int howmany> |
| 152 | +struct MixedBCController { |
| 153 | + bool mixed_active = false; |
| 154 | + |
| 155 | + protected: |
| 156 | + MixedBC mbc_local; |
| 157 | + size_t step_idx = 0; |
| 158 | + VectorXd g0_vec; // current macro strain (size n_str) |
| 159 | + |
| 160 | + // call from user code after v_u update each iteration |
| 161 | + template <typename SolverType> |
| 162 | + void update(SolverType &solver) |
| 163 | + { |
| 164 | + if (!mixed_active) |
| 165 | + return; |
| 166 | + |
| 167 | + VectorXd Pbar = solver.get_homogenized_stress(); |
| 168 | + |
| 169 | + VectorXd PF = (mbc_local.idx_F.size() ? mbc_local.P_F_path.row(step_idx).transpose() : VectorXd()); |
| 170 | + |
| 171 | + if (mbc_local.idx_F.size()) { |
| 172 | + VectorXd rhs = PF - mbc_local.Q_F.transpose() * Pbar; |
| 173 | + VectorXd delta_E = mbc_local.M * rhs; |
| 174 | + g0_vec += mbc_local.Q_F * delta_E; |
| 175 | + } |
| 176 | + |
| 177 | + vector<double> gvec(g0_vec.data(), g0_vec.data() + g0_vec.size()); |
| 178 | + solver.matmodel->setGradient(gvec); |
| 179 | + } |
| 180 | + |
| 181 | + template <typename SolverType> |
| 182 | + void activate(SolverType &solver, const MixedBC &mbc_in, size_t t) |
| 183 | + { |
| 184 | + mixed_active = true; |
| 185 | + mbc_local = mbc_in; |
| 186 | + step_idx = t; |
| 187 | + mbc_local.finalize(solver.matmodel->kapparef_mat); |
| 188 | + |
| 189 | + int n_str = solver.matmodel->n_str; |
| 190 | + g0_vec = VectorXd::Zero(n_str); |
| 191 | + if (mbc_local.idx_E.size()) |
| 192 | + g0_vec += mbc_local.Q_E * mbc_local.F_E_path.row(t).transpose(); |
| 193 | + |
| 194 | + vector<double> gvec(g0_vec.data(), g0_vec.data() + n_str); |
| 195 | + solver.matmodel->setGradient(gvec); |
| 196 | + |
| 197 | + // Do one update to set the initial strain |
| 198 | + solver.updateMixedBC(); |
| 199 | + } |
| 200 | +}; |
| 201 | + |
| 202 | +#endif // MIXED_BC_H |
0 commit comments