1 /*! 2 * \file mfront/src/NonLinearSystemSolverBase.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \brief 22 août 2014 6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 7 * reserved. 8 * This project is publicly released under either the GNU GPL Licence 9 * or the CECILL-A licence. A copy of thoses licences are delivered 10 * with the sources of TFEL. CEA or EDF may also distribute this 11 * project under specific licensing conditions. 12 */ 13 14 #include <cmath> 15 #include <limits> 16 #include <sstream> 17 #include "TFEL/Raise.hxx" 18 #include "MFront/BehaviourDescription.hxx" 19 #include "MFront/NonLinearSystemSolverBase.hxx" 20 21 namespace mfront { 22 getReservedNames() const23 std::vector<std::string> NonLinearSystemSolverBase::getReservedNames() const { 24 return {"jacobian_error", "error", "iter", "iterMax", "converged"}; 25 } // end of NonLinearSystemSolverBase::getReservedNames 26 getJacobianPart(const BehaviourDescription & mb,const VariableDescription & v1,const VariableDescription & v2,const SupportedTypes::TypeSize & n,const SupportedTypes::TypeSize & n2,const SupportedTypes::TypeSize & n3,const std::string & j,const std::string & p)27 std::string NonLinearSystemSolverBase::getJacobianPart( 28 const BehaviourDescription& mb, 29 const VariableDescription& v1, 30 const VariableDescription& v2, 31 const SupportedTypes::TypeSize& n, 32 const SupportedTypes::TypeSize& n2, 33 const SupportedTypes::TypeSize& n3, 34 const std::string& j, 35 const std::string& p) { 36 auto lthrow = [](const std::string& m) { 37 tfel::raise("NonLinearSystemSolverBase::getJacobianPart: " + m); 38 }; 39 std::ostringstream d; 40 if (mb.getTypeFlag(v1.type) == SupportedTypes::STENSOR) { 41 if (mb.getTypeFlag(v2.type) == SupportedTypes::STENSOR) { 42 d << "typename tfel::math::ST2toST2FromTinyMatrixView<N," << n2 << "," 43 << n2 << ",\n" 44 << n << "," << n3 << ",real>::type " + p + "df" << v1.name << "_dd" 45 << v2.name << "(" + j + ");\n"; 46 } else if (mb.getTypeFlag(v2.type) == SupportedTypes::SCALAR) { 47 d << "typename tfel::math::StensorFromTinyMatrixColumnView<N," << n2 48 << "," << n2 << ",\n" 49 << n << "," << n3 << ",real>::type " + p + "df" << v1.name << "_dd" 50 << v2.name << "(" + j + ");\n"; 51 } else { 52 lthrow("unsupported type for integration variable '" + v2.name + "'"); 53 } 54 } else if (mb.getTypeFlag(v1.type) == SupportedTypes::SCALAR) { 55 if (mb.getTypeFlag(v2.type) == SupportedTypes::STENSOR) { 56 d << "typename tfel::math::StensorFromTinyMatrixRowView<N," << n2 << "," 57 << n2 << ",\n" 58 << n << "," << n3 << ",real>::type " + p + "df" << v1.name << "_dd" 59 << v2.name << "(" + j + ");\n"; 60 } else if (mb.getTypeFlag(v2.type) == SupportedTypes::SCALAR) { 61 d << "real& " + p + "df" << v1.name << "_dd" << v2.name 62 << " = " + j + "(" << n << "," << n3 << ");\n"; 63 } else { 64 lthrow("unsupported type for integration variable '" + v2.name + "'"); 65 } 66 } else { 67 lthrow("unsupported type for integration variable '" + v2.name + "'"); 68 } 69 return d.str(); 70 } // end of NonLinearSystemSolverBase::getJacobianPart 71 writeEvaluateNumericallyComputedBlocks(std::ostream & out,const BehaviourDescription & mb,const Hypothesis h)72 void NonLinearSystemSolverBase::writeEvaluateNumericallyComputedBlocks( 73 std::ostream& out, const BehaviourDescription& mb, const Hypothesis h) { 74 auto throw_if = [](const bool c, const std::string& m) { 75 tfel::raise_if(c, 76 "NonLinearSystemSolverBase::" 77 "writeEvaluateNumericallyComputedBlocks: " + 78 m); 79 }; 80 const auto& d = mb.getBehaviourData(h); 81 if (!d.hasAttribute(BehaviourData::numericallyComputedJacobianBlocks)) { 82 return; 83 } 84 // First sort blocks by column 85 const auto blocs = [&d, throw_if] { 86 auto r = std::map<std::string, std::vector<std::string>>{}; 87 const auto& jbs = d.getAttribute<std::vector<std::string>>( 88 BehaviourData::numericallyComputedJacobianBlocks); 89 for (const auto& jb : jbs) { 90 const auto nd = [throw_if, 91 &jb]() -> std::pair<std::string, std::string> { 92 throw_if(jb.empty(), "empty jacobian block"); 93 throw_if(jb.size() < 6, "invalid jacobian block '" + jb + "'"); 94 throw_if(jb[0] != 'd', "invalid jacobian block '" + jb + "'"); 95 throw_if(jb[1] != 'f', "invalid jacobian block '" + jb + "'"); 96 const auto p = jb.find('_'); 97 throw_if(p == std::string::npos, 98 "invalid jacobian block '" + jb + "'"); 99 throw_if(p + 2 >= jb.size(), "invalid jacobian block '" + jb + "'"); 100 throw_if(jb[p + 1] != 'd', "invalid jacobian block '" + jb + "'"); 101 throw_if(jb[p + 2] != 'd', "invalid jacobian block '" + jb + "'"); 102 const auto nn = jb.substr(2, p - 2); 103 const auto dn = jb.substr(p + 3); 104 throw_if(nn.empty(), "invalid jacobian block '" + jb + "'"); 105 throw_if(dn.empty(), "invalid jacobian block '" + jb + "'"); 106 return {nn, dn}; 107 }(); 108 throw_if(!d.isIntegrationVariableName(nd.first), 109 "invalid jacobian block '" + jb + 110 "', " 111 "'" + 112 nd.first + "' is not an integration variable"); 113 throw_if(!d.isIntegrationVariableName(nd.second), 114 "invalid jacobian block '" + jb + 115 "', " 116 "'" + 117 nd.second + "' is not an integration variable"); 118 r[nd.second].push_back(nd.first); 119 } 120 return r; 121 }(); 122 const auto& ivs = d.getIntegrationVariables(); 123 const auto n = ivs.getTypeSize(); 124 out << "tmatrix<" << n << "," << n << ",real> tjacobian(this->jacobian);\n" 125 << "tvector<" << n << ",real> tfzeros(this->fzeros);\n"; 126 bool first = true; 127 for (const auto& b : blocs) { 128 auto getPositionAndSize = [&ivs, throw_if](const std::string& v) 129 -> std::pair<SupportedTypes::TypeSize, SupportedTypes::TypeSize> { 130 auto ns = SupportedTypes::TypeSize{}; 131 for (const auto& iv : ivs) { 132 const auto s = SupportedTypes::getTypeSize(iv.type, iv.arraySize); 133 if (iv.name == v) { 134 return {ns, s}; 135 } 136 ns += s; 137 } 138 throw_if(true, "no integration variable named '" + v + "'"); 139 }; 140 const auto pd = getPositionAndSize(b.first); 141 auto update_jacobian = [&out, &b, 142 getPositionAndSize](const std::string& j) { 143 // fzeros contain the new jacobian value 144 for (const auto& v : b.second) { 145 const auto pn = getPositionAndSize(v); 146 if (pn.second.isOne()) { 147 out << "tjacobian(" << pn.first << "," << j << ") = " 148 << "this->fzeros[" << pn.first << "];\n"; 149 } else { 150 const auto i = [&pn]() -> std::string { 151 if (!pn.first.isNull()) { 152 return to_string(pn.first) + "+idx2"; 153 } 154 return "idx2"; 155 }(); 156 out << "for(unsigned short idx2 = 0; idx2!= " << pn.second 157 << ";++idx2){\n" 158 << "tjacobian(" << i << "," << j << ") = " 159 << "this->fzeros[" << i << "];\n" 160 << "}\n"; 161 } 162 } 163 }; 164 auto compute_perturbation = [&out, &update_jacobian, &d, 165 &n](const std::string& j) { 166 out << "this->zeros(" << j 167 << ") -= this->numerical_jacobian_epsilon;\n"; 168 if (d.hasCode(BehaviourData::ComputeThermodynamicForces)) { 169 out << "this->computeThermodynamicForces();\n"; 170 } 171 out << "this->computeFdF(true);\n" 172 << "this->zeros = this->zeros_1;\n" 173 << "tvector<" << n << ",real> tfzeros2(this->fzeros);\n" 174 << "this->zeros(" << j 175 << ") += this->numerical_jacobian_epsilon;\n"; 176 if (d.hasCode(BehaviourData::ComputeThermodynamicForces)) { 177 out << "this->computeThermodynamicForces();\n"; 178 } 179 out << "this->computeFdF(true);\n" 180 << "this->zeros = this->zeros_1;\n" 181 << "this->fzeros = " 182 "(this->fzeros-tfzeros2)/" 183 "(2*(this->numerical_jacobian_epsilon));\n" 184 << "// update jacobian\n"; 185 update_jacobian(j); 186 }; 187 if (!first) { 188 out << "// restore residual values\n" 189 << "this->fzeros = tfzeros;\n"; 190 } 191 if (pd.second.isOne()) { 192 compute_perturbation(to_string(pd.first)); 193 } else { 194 out << "for(unsigned short idx = 0; idx!= " << pd.second 195 << ";++idx){\n"; 196 if (!pd.first.isNull()) { 197 compute_perturbation(to_string(pd.first) + "+idx"); 198 } else { 199 compute_perturbation("idx"); 200 } 201 out << "}\n"; 202 } 203 first = false; 204 } 205 out << "// update jacobian\n" 206 << "this->jacobian = tjacobian;\n" 207 << "// restore residual values\n" 208 << "this->fzeros = tfzeros;\n"; 209 } // end of 210 // NonLinearSystemSolverBase::writeEvaluateNumericallyComputedBlocks 211 writeComparisonToNumericalJacobian(std::ostream & out,const BehaviourDescription & mb,const Hypothesis h,const std::string & nj)212 void NonLinearSystemSolverBase::writeComparisonToNumericalJacobian( 213 std::ostream& out, 214 const BehaviourDescription& mb, 215 const Hypothesis h, 216 const std::string& nj) { 217 if (!mb.getAttribute(h, BehaviourData::compareToNumericalJacobian, false)) { 218 return; 219 } 220 const auto& d = mb.getBehaviourData(h); 221 SupportedTypes::TypeSize n; 222 SupportedTypes::TypeSize n3; 223 const auto n2 = d.getIntegrationVariables().getTypeSize(); 224 out << "auto jacobian_error = real{};" 225 << "this->computeNumericalJacobian(" << nj << ");\n"; 226 n = SupportedTypes::TypeSize(); 227 for (const auto& v : d.getIntegrationVariables()) { 228 if (v.arraySize == 1) { 229 n3 = SupportedTypes::TypeSize(); 230 for (const auto& v2 : d.getIntegrationVariables()) { 231 if ((v.arraySize == 1) && (v2.arraySize == 1)) { 232 out << "// derivative of variable f" << v.name << " by variable " 233 << v2.name << "\n" 234 << NonLinearSystemSolverBase::getJacobianPart(mb, v, v2, n, n2, 235 n3) 236 << "// numerical derivative of variable f" << v.name 237 << " by variable " << v2.name << "\n" 238 << NonLinearSystemSolverBase::getJacobianPart(mb, v, v2, n, n2, 239 n3, nj, "n"); 240 n3 += mb.getTypeSize(v2.type, v2.arraySize); 241 } 242 } 243 } 244 n += mb.getTypeSize(v.type, v.arraySize); 245 } 246 for (const auto& v1 : d.getIntegrationVariables()) { 247 for (const auto& v2 : d.getIntegrationVariables()) { 248 const auto nv1 = mb.getTypeSize(v1.type, 1u); 249 const auto nv2 = mb.getTypeSize(v2.type, 1u); 250 out << "jacobian_error=" << nv1 << "*" << nv2 << "*" 251 << "(this->jacobianComparisonCriterion)" 252 << ";\n"; 253 if ((v1.arraySize == 1u) && (v2.arraySize == 1u)) { 254 out << "if(abs(" 255 << "df" << v1.name << "_dd" << v2.name << "-" 256 << "ndf" << v1.name << "_dd" << v2.name << ") > jacobian_error)\n" 257 << "{\n" 258 << "cout << abs(" 259 << "df" << v1.name << "_dd" << v2.name << "-" 260 << "ndf" << v1.name << "_dd" << v2.name 261 << ") << \" \" << jacobian_error << endl;\n" 262 << "cout << \"df" << v1.name << "_dd" << v2.name << " :\\n\" << " 263 << "df" << v1.name << "_dd" << v2.name << " << endl;\n" 264 << "cout << \"ndf" << v1.name << "_dd" << v2.name << " :\\n\" << " 265 << "ndf" << v1.name << "_dd" << v2.name << " << endl;\n" 266 << "cout << \"df" << v1.name << "_dd" << v2.name << " - ndf" 267 << v1.name << "_dd" << v2.name << " :\\n\" << " 268 << "df" << v1.name << "_dd" << v2.name << "-" 269 << "ndf" << v1.name << "_dd" << v2.name << " << endl;\n" 270 << "cout << endl;\n" 271 << "}\n"; 272 } else if (((v1.arraySize != 1u) && (v2.arraySize == 1u)) || 273 ((v2.arraySize != 1u) && (v1.arraySize == 1u))) { 274 const auto asize = (v1.arraySize != 1u) ? v1.arraySize : v2.arraySize; 275 out << "for(unsigned short idx=0;idx!=" << asize << ";++idx){\n" 276 << "if(abs(" 277 << "df" << v1.name << "_dd" << v2.name << "(idx)-" 278 << "df" << v1.name << "_dd" << v2.name << "(" << nj 279 << ",idx)) > jacobian_error)\n" 280 << "{\n" 281 << "cout << abs(" 282 << "df" << v1.name << "_dd" << v2.name << "(idx)-" 283 << "df" << v1.name << "_dd" << v2.name << "(" << nj 284 << ",idx)) << \" \" << jacobian_error << endl;\n" 285 << "cout << \"df" << v1.name << "_dd" << v2.name 286 << "(\" << idx << \") :\\n\" << " 287 << "df" << v1.name << "_dd" << v2.name << "(idx) << endl;\n" 288 << "cout << \"ndf" << v1.name << "_dd" << v2.name 289 << "(\" << idx << \") :\\n\" << " 290 << "df" << v1.name << "_dd" << v2.name << "(" << nj 291 << ",idx) << endl;\n" 292 << "cout << \"df" << v1.name << "_dd" << v2.name << " - ndf" 293 << v1.name << "_dd" << v2.name << "(\" << idx << \") :\\n\" << " 294 << "df" << v1.name << "_dd" << v2.name << "(idx) -" 295 << "df" << v1.name << "_dd" << v2.name << "(" << nj 296 << ",idx) << endl;\n" 297 << "cout << endl;\n" 298 << "}\n" 299 << "}\n"; 300 } else { 301 out << "for(unsigned short idx=0;idx!=" << v1.arraySize 302 << ";++idx){\n" 303 << "for(unsigned short idx2=0;idx2!=" << v2.arraySize 304 << ";++idx2){\n" 305 << "if(abs(" 306 << "df" << v1.name << "_dd" << v2.name << "(idx,idx2)-" 307 << "df" << v1.name << "_dd" << v2.name << "(" << nj 308 << ",idx,idx2)) > jacobian_error)\n" 309 << "{\n" 310 << "cout << abs(" 311 << "df" << v1.name << "_dd" << v2.name << "(idx,idx2)-" 312 << "df" << v1.name << "_dd" << v2.name << "(" << nj 313 << ",idx,idx2)) << \" \" << jacobian_error << endl;\n" 314 << "cout << \"df" << v1.name << "_dd" << v2.name 315 << "(\" << idx << \",\" << idx2 << \") :\\n\" << " 316 << "df" << v1.name << "_dd" << v2.name << "(idx,idx2) << endl;\n" 317 << "cout << \"ndf" << v1.name << "_dd" << v2.name 318 << "(\" << idx << \",\" << idx2 << \") :\\n\" << " 319 << "df" << v1.name << "_dd" << v2.name << "(" << nj 320 << ",idx,idx2) << endl;\n" 321 << "cout << \"df" << v1.name << "_dd" << v2.name << " - ndf" 322 << v1.name << "_dd" << v2.name 323 << "(\" << idx << \",\" << idx2 << \") :\\n\" << " 324 << "df" << v1.name << "_dd" << v2.name << "(idx,idx2) -" 325 << "df" << v1.name << "_dd" << v2.name << "(" << nj 326 << ",idx,idx2) << endl;\n" 327 << "cout << endl;\n" 328 << "}\n" 329 << "}\n" 330 << "}\n"; 331 } 332 } 333 } 334 } // end of 335 // NonLinearSystemSolverBase::writeComparisonToNumericalJacobian 336 writeLimitsOnIncrementValues(std::ostream & out,const BehaviourDescription & mb,const Hypothesis h,const std::string & v)337 void NonLinearSystemSolverBase::writeLimitsOnIncrementValues( 338 std::ostream& out, 339 const BehaviourDescription& mb, 340 const Hypothesis h, 341 const std::string& v) { 342 const auto& d = mb.getBehaviourData(h); 343 SupportedTypes::TypeSize n; 344 for (const auto& iv : d.getIntegrationVariables()) { 345 const auto nv = mb.getTypeSize(iv.type, iv.arraySize); 346 if (mb.hasParameter(h, 347 iv.name + "_maximum_increment_value_per_iteration")) { 348 out << "for(unsigned short idx = 0; idx!=" << nv << ";++idx){\n"; 349 if (mb.hasAttribute(h, iv.name + "_normalisation_factor")) { 350 const auto& nf = mb.getAttribute<std::string>( 351 h, iv.name + "_normalisation_factor"); 352 out << "if(std::abs(" << v << "[" << n << "+idx])>" << nf 353 << "*(this->" << iv.name 354 << "_maximum_increment_value_per_iteration)){\n"; 355 } else { 356 out << "if(std::abs(" << v << "[" << n << "+idx])>this->" << iv.name 357 << "_maximum_increment_value_per_iteration){\n"; 358 } 359 out << "if(" << v << "[" << n << "+idx]<0){\n"; 360 if (mb.hasAttribute(h, iv.name + "_normalisation_factor")) { 361 const auto& nf = mb.getAttribute<std::string>( 362 h, iv.name + "_normalisation_factor"); 363 out << "" << v << "[" << n << "+idx] = -" << nf << "*(this->" 364 << iv.name << "_maximum_increment_value_per_iteration);\n"; 365 } else { 366 out << "" << v << "[" << n << "+idx] = -this->" << iv.name 367 << "_maximum_increment_value_per_iteration;\n"; 368 } 369 out << "} else {\n"; 370 if (mb.hasAttribute(h, iv.name + "_normalisation_factor")) { 371 const auto& nf = mb.getAttribute<std::string>( 372 h, iv.name + "_normalisation_factor"); 373 out << "" << v << "[" << n << "+idx] = " << nf << "*(this->" 374 << iv.name << "_maximum_increment_value_per_iteration);\n"; 375 } else { 376 out << "" << v << "[" << n << "+idx] = this->" << iv.name 377 << "_maximum_increment_value_per_iteration;\n"; 378 } 379 out << "}\n"; 380 out << "}\n"; 381 out << "}\n"; 382 } else if (mb.hasParameter(h, "maximum_increment_value_per_iteration")) { 383 out << "for(unsigned short idx = 0; idx!=" << nv << ";++idx){\n"; 384 out << "if(std::abs(" << v << "[" << n 385 << "+idx])>this->maximum_increment_value_per_iteration){\n"; 386 out << "if(" << v << "[" << n << "+idx]<0){\n"; 387 out << "" << v << "[" << n << "+idx] = " 388 "-this->maximum_increment_value_" 389 "per_iteration;\n"; 390 out << "} else {\n"; 391 out << "" << v << "[" << n << "+idx] = " 392 "this->maximum_increment_value_" 393 "per_iteration;\n"; 394 out << "}\n"; 395 out << "}\n"; 396 out << "}\n"; 397 } 398 n += nv; 399 } 400 } // end of NonLinearSystemSolverBase::writeLimitsOnIncrementValues 401 402 void NonLinearSystemSolverBase:: writeLimitsOnIncrementValuesBasedOnStateVariablesPhysicalBounds(std::ostream & out,const BehaviourDescription & mb,const Hypothesis h)403 writeLimitsOnIncrementValuesBasedOnStateVariablesPhysicalBounds( 404 std::ostream& out, 405 const BehaviourDescription& mb, 406 const Hypothesis h) { 407 const auto& d = mb.getBehaviourData(h); 408 auto n = SupportedTypes::TypeSize{}; 409 for (const auto& v : d.getStateVariables()) { 410 if (!v.hasPhysicalBounds()) { 411 n += mb.getTypeSize(v.type, v.arraySize); 412 continue; 413 } 414 const auto& b = v.getPhysicalBounds(); 415 // treating lower bounds 416 if ((b.boundsType == VariableBoundsDescription::LOWER) || 417 (b.boundsType == VariableBoundsDescription::LOWERANDUPPER)) { 418 if ((mb.getTypeFlag(v.type) == SupportedTypes::SCALAR) && 419 (v.arraySize == 1u)) { 420 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 421 const auto& nf = mb.getAttribute<std::string>( 422 h, v.name + "_normalisation_factor"); 423 out << "if(this->" << v.name << "+ " << nf << "*(this->zeros[" << n 424 << "]) <" << b.lowerBound << "){\n"; 425 } else { 426 out << "if(this->" << v.name << "+this->zeros[" << n << "]<" 427 << b.lowerBound << "){\n"; 428 } 429 if (std::abs(b.lowerBound) > 430 std::numeric_limits<long double>::min()) { 431 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 432 const auto& nf = mb.getAttribute<std::string>( 433 h, v.name + "_normalisation_factor"); 434 out << "this->zeros[" << n << "] = (" << b.lowerBound 435 << "- (this->" << v.name << "))/(" << nf << ");\n"; 436 } else { 437 out << "this->zeros[" << n << "] = " << b.lowerBound 438 << "- (this->" << v.name << ");\n"; 439 } 440 } else { 441 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 442 const auto& nf = mb.getAttribute<std::string>( 443 h, v.name + "_normalisation_factor"); 444 out << "this->zeros[" << n << "] = - (this->" << v.name << ")/(" 445 << nf << ");\n"; 446 } else { 447 out << "this->zeros[" << n << "] = - (this->" << v.name << ");\n"; 448 } 449 } 450 out << "}\n"; 451 } 452 if ((mb.getTypeFlag(v.type) == SupportedTypes::SCALAR) && 453 (v.arraySize != 1u)) { 454 out << "for(unsigned short idx=0;idx!=" << v.arraySize 455 << ";++idx){\n"; 456 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 457 const auto& nf = mb.getAttribute<std::string>( 458 h, v.name + "_normalisation_factor"); 459 out << "if(this->" << v.name << "[idx]+(" << nf << ")*(this->zeros[" 460 << n << "+idx])<" << b.lowerBound << "){\n"; 461 } else { 462 out << "if(this->" << v.name << "[idx]+this->zeros[" << n 463 << "+idx]<" << b.lowerBound << "){\n"; 464 } 465 if (std::abs(b.lowerBound) > 466 std::numeric_limits<long double>::min()) { 467 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 468 const auto& nf = mb.getAttribute<std::string>( 469 h, v.name + "_normalisation_factor"); 470 out << "this->zeros[" << n << "+idx] = (" << b.lowerBound 471 << "- (this->" << v.name << "[idx]))/(" << nf << ");\n"; 472 } else { 473 out << "this->zeros[" << n << "+idx] = " << b.lowerBound 474 << "- (this->" << v.name << "[idx]);\n"; 475 } 476 } else { 477 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 478 const auto& nf = mb.getAttribute<std::string>( 479 h, v.name + "_normalisation_factor"); 480 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 481 << "[idx])/(" << nf << ");\n"; 482 } else { 483 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 484 << "[idx]);\n"; 485 } 486 } 487 out << "}\n" 488 << "}\n"; 489 } 490 if ((mb.getTypeFlag(v.type) != SupportedTypes::SCALAR) && 491 (v.arraySize == 1u)) { 492 if (b.component == -1) { 493 const auto n2 = mb.getTypeSize(v.type, 1u); 494 out << "for(unsigned short idx=0;idx!=" << n2 << ";++idx){\n"; 495 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 496 const auto& nf = mb.getAttribute<std::string>( 497 h, v.name + "_normalisation_factor"); 498 out << "if(this->" << v.name << "[idx]+(" << nf 499 << ")*(this->zeros[" << n << "+idx])<" << b.lowerBound 500 << "){\n"; 501 } else { 502 out << "if(this->" << v.name << "[idx]+this->zeros[" << n 503 << "+idx]<" << b.lowerBound << "){\n"; 504 } 505 if (std::abs(b.lowerBound) > 506 std::numeric_limits<long double>::min()) { 507 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 508 const auto& nf = mb.getAttribute<std::string>( 509 h, v.name + "_normalisation_factor"); 510 out << "this->zeros[" << n << "+idx] = (" << b.lowerBound 511 << "- (this->" << v.name << "[idx]))/(" << nf << ");\n"; 512 } else { 513 out << "this->zeros[" << n << "+idx] = " << b.lowerBound 514 << "- (this->" << v.name << "[idx]);\n"; 515 } 516 } else { 517 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 518 const auto& nf = mb.getAttribute<std::string>( 519 h, v.name + "_normalisation_factor"); 520 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 521 << "[idx])/(" << nf << ");\n"; 522 } else { 523 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 524 << "[idx]);\n"; 525 } 526 } 527 out << "}\n" 528 << "}\n"; 529 } else { 530 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 531 const auto& nf = mb.getAttribute<std::string>( 532 h, v.name + "_normalisation_factor"); 533 out << "if(this->" << v.name << "[" << b.component << "]+(" << nf 534 << ")*(this->zeros[" << n << "+" << b.component << "])<" 535 << b.lowerBound << "){\n"; 536 } else { 537 out << "if(this->" << v.name << "[" << b.component 538 << "]+this->zeros[" << n << "+" << b.component << "]<" 539 << b.lowerBound << "){\n"; 540 } 541 if (std::abs(b.lowerBound) > 542 std::numeric_limits<long double>::min()) { 543 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 544 const auto& nf = mb.getAttribute<std::string>( 545 h, v.name + "_normalisation_factor"); 546 out << "this->zeros[" << n << "+" << b.component << "] = (" 547 << b.lowerBound << "- (this->" << v.name << "[" 548 << b.component << "]))/(" << nf << ");\n"; 549 } else { 550 out << "this->zeros[" << n << "+" << b.component 551 << "] = " << b.lowerBound << "- (this->" << v.name << "[" 552 << b.component << "]);\n"; 553 } 554 } else { 555 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 556 const auto& nf = mb.getAttribute<std::string>( 557 h, v.name + "_normalisation_factor"); 558 out << "this->zeros[" << n << "+" << b.component 559 << "] = - (this->" << v.name << "[" << b.component << "])/(" 560 << nf << ");\n"; 561 } else { 562 out << "this->zeros[" << n << "+" << b.component 563 << "] = - (this->" << v.name << "[" << b.component 564 << "]);\n"; 565 } 566 } 567 out << "}\n"; 568 } 569 } 570 if ((mb.getTypeFlag(v.type) != SupportedTypes::SCALAR) && 571 (v.arraySize != 1u)) { 572 const auto n2 = mb.getTypeSize(v.type, 1u); 573 out << "for(unsigned short idx=0;idx!=" << v.arraySize 574 << ";++idx){\n"; 575 if (b.component == -1) { 576 out << "for(unsigned short idx2=0;idx2!=" << n2 << ";++idx2){\n"; 577 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 578 const auto& nf = mb.getAttribute<std::string>( 579 h, v.name + "_normalisation_factor"); 580 out << "if(this->" << v.name << "[idx][idx2]+(" << nf 581 << ")*(this->zeros[" << n << "+" << n2 << "*idx+idx2])<" 582 << b.lowerBound << "){\n"; 583 } else { 584 out << "if(this->" << v.name << "[idx][idx2]+this->zeros[" << n 585 << "+" << n2 << "*idx+idx2]<" << b.lowerBound << "){\n"; 586 } 587 if (std::abs(b.lowerBound) > 588 std::numeric_limits<long double>::min()) { 589 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 590 const auto& nf = mb.getAttribute<std::string>( 591 h, v.name + "_normalisation_factor"); 592 out << "this->zeros[" << n << "+" << n2 << "*idx+idx2] = (" 593 << b.lowerBound << "- (this->" << v.name 594 << "[idx][idx2]))/(" << nf << ");\n"; 595 } else { 596 out << "this->zeros[" << n << "+" << n2 597 << "*idx+idx2] = " << b.lowerBound << "- (this->" << v.name 598 << "[idx][idx2]);\n"; 599 } 600 } else { 601 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 602 const auto& nf = mb.getAttribute<std::string>( 603 h, v.name + "_normalisation_factor"); 604 out << "this->zeros[" << n << "+" << n2 605 << "*idx+idx2] = - (this->" << v.name << "[idx][idx2])/(" 606 << nf << ");\n"; 607 } else { 608 out << "this->zeros[" << n << "+" << n2 609 << "*idx+idx2] = - (this->" << v.name << "[idx][idx2]);\n"; 610 } 611 } 612 out << "}\n" 613 << "}\n"; 614 } else { 615 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 616 const auto& nf = mb.getAttribute<std::string>( 617 h, v.name + "_normalisation_factor"); 618 out << "if(this->" << v.name << "[idx][" << b.component << "]+(" 619 << nf << ")*(this->zeros[" << n << "+" << n2 << "*idx+" 620 << b.component << "])<" << b.lowerBound << "){\n"; 621 } else { 622 out << "if(this->" << v.name << "[idx][" << b.component 623 << "]+this->zeros[" << n << "+" << n2 << "*idx+" 624 << b.component << "]<" << b.lowerBound << "){\n"; 625 } 626 if (std::abs(b.lowerBound) > 627 std::numeric_limits<long double>::min()) { 628 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 629 const auto& nf = mb.getAttribute<std::string>( 630 h, v.name + "_normalisation_factor"); 631 out << "this->zeros[" << n << "+" << n2 << "*idx+" 632 << b.component << "] = (" << b.lowerBound << "- (this->" 633 << v.name << "[idx][" << b.component << "]))/(" << nf 634 << ");\n"; 635 } else { 636 out << "this->zeros[" << n << "+" << n2 << "*idx+" 637 << b.component << "] = " << b.lowerBound << "- (this->" 638 << v.name << "[idx][" << b.component << "]);\n"; 639 } 640 } else { 641 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 642 const auto& nf = mb.getAttribute<std::string>( 643 h, v.name + "_normalisation_factor"); 644 out << "this->zeros[" << n << "+" << n2 << "*idx+" 645 << b.component << "] = - (this->" << v.name << "[idx][" 646 << b.component << "])/(" << nf << ");\n"; 647 } else { 648 out << "this->zeros[" << n << "+" << n2 << "*idx+" 649 << b.component << "] = - (this->" << v.name << "[idx][" 650 << b.component << "]);\n"; 651 } 652 } 653 out << "}\n"; 654 } 655 out << "}\n"; 656 } 657 } 658 // treating upper bounds 659 if ((b.boundsType == VariableBoundsDescription::UPPER) || 660 (b.boundsType == VariableBoundsDescription::LOWERANDUPPER)) { 661 if ((mb.getTypeFlag(v.type) == SupportedTypes::SCALAR) && 662 (v.arraySize == 1u)) { 663 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 664 const auto& nf = mb.getAttribute<std::string>( 665 h, v.name + "_normalisation_factor"); 666 out << "if(this->" << v.name << "+ " << nf << "*(this->zeros[" << n 667 << "]) >" << b.upperBound << "){\n"; 668 } else { 669 out << "if(this->" << v.name << "+this->zeros[" << n << "]>" 670 << b.upperBound << "){\n"; 671 } 672 if (std::abs(b.upperBound) > 673 std::numeric_limits<long double>::min()) { 674 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 675 const auto& nf = mb.getAttribute<std::string>( 676 h, v.name + "_normalisation_factor"); 677 out << "this->zeros[" << n << "] = (" << b.upperBound 678 << "- (this->" << v.name << "))/(" << nf << ");\n"; 679 } else { 680 out << "this->zeros[" << n << "] = " << b.upperBound 681 << "- (this->" << v.name << ");\n"; 682 } 683 } else { 684 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 685 const auto& nf = mb.getAttribute<std::string>( 686 h, v.name + "_normalisation_factor"); 687 out << "this->zeros[" << n << "] = - (this->" << v.name << ")/(" 688 << nf << ");\n"; 689 } else { 690 out << "this->zeros[" << n << "] = - (this->" << v.name << ");\n"; 691 } 692 } 693 out << "}\n"; 694 } 695 if ((mb.getTypeFlag(v.type) == SupportedTypes::SCALAR) && 696 (v.arraySize != 1u)) { 697 out << "for(unsigned short idx=0;idx!=" << v.arraySize 698 << ";++idx){\n"; 699 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 700 const auto& nf = mb.getAttribute<std::string>( 701 h, v.name + "_normalisation_factor"); 702 out << "if(this->" << v.name << "[idx]+(" << nf << ")*(this->zeros[" 703 << n << "+idx])>" << b.upperBound << "){\n"; 704 } else { 705 out << "if(this->" << v.name << "[idx]+this->zeros[" << n 706 << "+idx]>" << b.upperBound << "){\n"; 707 } 708 if (std::abs(b.upperBound) > 709 std::numeric_limits<long double>::min()) { 710 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 711 const auto& nf = mb.getAttribute<std::string>( 712 h, v.name + "_normalisation_factor"); 713 out << "this->zeros[" << n << "+idx] = (" << b.upperBound 714 << "- (this->" << v.name << "[idx]))/(" << nf << ");\n"; 715 } else { 716 out << "this->zeros[" << n << "+idx] = " << b.upperBound 717 << "- (this->" << v.name << "[idx]);\n"; 718 } 719 } else { 720 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 721 const auto& nf = mb.getAttribute<std::string>( 722 h, v.name + "_normalisation_factor"); 723 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 724 << "[idx])/(" << nf << ");\n"; 725 } else { 726 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 727 << "[idx]);\n"; 728 } 729 } 730 out << "}\n"; 731 out << "}\n"; 732 } 733 if ((mb.getTypeFlag(v.type) != SupportedTypes::SCALAR) && 734 (v.arraySize == 1u)) { 735 if (b.component == -1) { 736 const auto n2 = mb.getTypeSize(v.type, 1u); 737 out << "for(unsigned short idx=0;idx!=" << n2 << ";++idx){\n"; 738 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 739 const auto& nf = mb.getAttribute<std::string>( 740 h, v.name + "_normalisation_factor"); 741 out << "if(this->" << v.name << "[idx]+(" << nf 742 << ")*(this->zeros[" << n << "+idx])>" << b.upperBound 743 << "){\n"; 744 } else { 745 out << "if(this->" << v.name << "[idx]+this->zeros[" << n 746 << "+idx]>" << b.upperBound << "){\n"; 747 } 748 if (std::abs(b.upperBound) > 749 std::numeric_limits<long double>::min()) { 750 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 751 const auto& nf = mb.getAttribute<std::string>( 752 h, v.name + "_normalisation_factor"); 753 out << "this->zeros[" << n << "+idx] = (" << b.upperBound 754 << "- (this->" << v.name << "[idx]))/(" << nf << ");\n"; 755 } else { 756 out << "this->zeros[" << n << "+idx] = " << b.upperBound 757 << "- (this->" << v.name << "[idx]);\n"; 758 } 759 } else { 760 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 761 const auto& nf = mb.getAttribute<std::string>( 762 h, v.name + "_normalisation_factor"); 763 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 764 << "[idx])/(" << nf << ");\n"; 765 } else { 766 out << "this->zeros[" << n << "+idx] = - (this->" << v.name 767 << "[idx]);\n"; 768 } 769 } 770 out << "}\n" 771 << "}\n"; 772 } else { 773 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 774 const auto& nf = mb.getAttribute<std::string>( 775 h, v.name + "_normalisation_factor"); 776 out << "if(this->" << v.name << "[" << b.component << "]+(" << nf 777 << ")*(this->zeros[" << n << "+" << b.component << "])>" 778 << b.upperBound << "){\n"; 779 } else { 780 out << "if(this->" << v.name << "[" << b.component 781 << "]+this->zeros[" << n << "+" << b.component << "]>" 782 << b.upperBound << "){\n"; 783 } 784 if (std::abs(b.upperBound) > 785 std::numeric_limits<long double>::min()) { 786 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 787 const auto& nf = mb.getAttribute<std::string>( 788 h, v.name + "_normalisation_factor"); 789 out << "this->zeros[" << n << "+" << b.component << "] = (" 790 << b.upperBound << "- (this->" << v.name << "[" 791 << b.component << "]))/(" << nf << ");\n"; 792 } else { 793 out << "this->zeros[" << n << "+" << b.component 794 << "] = " << b.upperBound << "- (this->" << v.name << "[" 795 << b.component << "]);\n"; 796 } 797 } else { 798 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 799 const auto& nf = mb.getAttribute<std::string>( 800 h, v.name + "_normalisation_factor"); 801 out << "this->zeros[" << n << "+" << b.component 802 << "] = - (this->" << v.name << "[" << b.component << "])/(" 803 << nf << ");\n"; 804 } else { 805 out << "this->zeros[" << n << "+" << b.component 806 << "] = - (this->" << v.name << "[" << b.component 807 << "]);\n"; 808 } 809 } 810 out << "}\n"; 811 } 812 } 813 if ((mb.getTypeFlag(v.type) != SupportedTypes::SCALAR) && 814 (v.arraySize != 1u)) { 815 const auto n2 = mb.getTypeSize(v.type, 1u); 816 out << "for(unsigned short idx=0;idx!=" << v.arraySize 817 << ";++idx){\n"; 818 if (b.component == -1) { 819 out << "for(unsigned short idx2=0;idx2!=" << n2 << ";++idx2){\n"; 820 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 821 const auto& nf = mb.getAttribute<std::string>( 822 h, v.name + "_normalisation_factor"); 823 out << "if(this->" << v.name << "[idx][idx2]+(" << nf 824 << ")*(this->zeros[" << n << "+" << n2 << "*idx+idx2])>" 825 << b.upperBound << "){\n"; 826 } else { 827 out << "if(this->" << v.name << "[idx][idx2]+this->zeros[" << n 828 << "+" << n2 << "*idx+idx2]>" << b.upperBound << "){\n"; 829 } 830 if (std::abs(b.upperBound) > 831 std::numeric_limits<long double>::min()) { 832 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 833 const auto& nf = mb.getAttribute<std::string>( 834 h, v.name + "_normalisation_factor"); 835 out << "this->zeros[" << n << "+" << n2 << "*idx+idx2] = (" 836 << b.upperBound << "- (this->" << v.name 837 << "[idx][idx2]))/(" << nf << ");\n"; 838 } else { 839 out << "this->zeros[" << n << "+" << n2 840 << "*idx+idx2] = " << b.upperBound << "- (this->" << v.name 841 << "[idx][idx2]);\n"; 842 } 843 } else { 844 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 845 const auto& nf = mb.getAttribute<std::string>( 846 h, v.name + "_normalisation_factor"); 847 out << "this->zeros[" << n << "+" << n2 848 << "*idx+idx2] = - (this->" << v.name << "[idx][idx2])/(" 849 << nf << ");\n"; 850 } else { 851 out << "this->zeros[" << n << "+" << n2 852 << "*idx+idx2] = - (this->" << v.name << "[idx][idx2]);\n"; 853 } 854 } 855 out << "}\n" 856 << "}\n"; 857 } else { 858 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 859 const auto& nf = mb.getAttribute<std::string>( 860 h, v.name + "_normalisation_factor"); 861 out << "if(this->" << v.name << "[idx][" << b.component << "]+(" 862 << nf << ")*(this->zeros[" << n << "+" << n2 << "*idx+" 863 << b.component << "])>" << b.upperBound << "){\n"; 864 } else { 865 out << "if(this->" << v.name << "[idx][" << b.component 866 << "]+this->zeros[" << n << "+" << n2 << "*idx+" 867 << b.component << "]>" << b.upperBound << "){\n"; 868 } 869 if (std::abs(b.upperBound) > 870 std::numeric_limits<long double>::min()) { 871 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 872 const auto& nf = mb.getAttribute<std::string>( 873 h, v.name + "_normalisation_factor"); 874 out << "this->zeros[" << n << "+" << n2 << "*idx+" 875 << b.component << "] = (" << b.upperBound << "- (this->" 876 << v.name << "[idx][" << b.component << "]))/(" << nf 877 << ");\n"; 878 } else { 879 out << "this->zeros[" << n << "+" << n2 << "*idx+" 880 << b.component << "] = " << b.upperBound << "- (this->" 881 << v.name << "[idx][" << b.component << "]);\n"; 882 } 883 } else { 884 if (mb.hasAttribute(h, v.name + "_normalisation_factor")) { 885 const auto& nf = mb.getAttribute<std::string>( 886 h, v.name + "_normalisation_factor"); 887 out << "this->zeros[" << n << "+" << n2 << "*idx+" 888 << b.component << "] = - (this->" << v.name << "[idx][" 889 << b.component << "])/(" << nf << ");\n"; 890 } else { 891 out << "this->zeros[" << n << "+" << n2 << "*idx+" 892 << b.component << "] = - (this->" << v.name << "[idx][" 893 << b.component << "]);\n"; 894 } 895 } 896 out << "}\n"; 897 } 898 out << "}\n"; 899 } 900 } 901 n += mb.getTypeSize(v.type, v.arraySize); 902 } 903 } // end of 904 // NonLinearSystemSolverBase::writeLimitsOnIncrementValuesBasedOnStateVariablesPhysicalBounds 905 906 void NonLinearSystemSolverBase:: writeLimitsOnIncrementValuesBasedOnIntegrationVariablesIncrementsPhysicalBounds(std::ostream &,const BehaviourDescription &,const Hypothesis)907 writeLimitsOnIncrementValuesBasedOnIntegrationVariablesIncrementsPhysicalBounds( 908 std::ostream&, const BehaviourDescription&, const Hypothesis) { 909 } // end of 910 // NonLinearSystemSolverBase::writeLimitsOnIncrementValuesBasedOnIntegrationVariablesIncrementsPhysicalBounds 911 912 NonLinearSystemSolverBase::~NonLinearSystemSolverBase() = default; 913 914 } // end of namespace mfront 915