1 /*! 2 * \file mfront/src/CalculiXInterface.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 17 Jan 2007 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 <fstream> 15 #include <sstream> 16 #include <cstdlib> 17 #include <stdexcept> 18 19 #include "TFEL/Raise.hxx" 20 #include "TFEL/Config/GetInstallPath.hxx" 21 #include "TFEL/Utilities/StringAlgorithms.hxx" 22 #include "TFEL/System/System.hxx" 23 24 #include "MFront/DSLUtilities.hxx" 25 #include "MFront/MFrontLock.hxx" 26 #include "MFront/MFrontUtilities.hxx" 27 #include "MFront/MFrontLogStream.hxx" 28 #include "MFront/MFrontDebugMode.hxx" 29 #include "MFront/FileDescription.hxx" 30 #include "MFront/TargetsDescription.hxx" 31 #include "MFront/CalculiXSymbolsGenerator.hxx" 32 #include "MFront/CalculiXInterface.hxx" 33 34 #ifndef _MSC_VER 35 static const char* const constexpr_c = "constexpr"; 36 #else 37 static const char* const constexpr_c = "const"; 38 #endif 39 40 namespace mfront { 41 checkFiniteStrainStrategy(const std::string & fs)42 static void checkFiniteStrainStrategy(const std::string& fs) { 43 tfel::raise_if((fs != "FiniteRotationSmallStrain") && 44 (fs != "MieheApelLambrechtLogarithmicStrain") && 45 (fs != "MieheApelLambrechtLogarithmicStrainII"), 46 "checkFiniteStrainStrategy: " 47 "unsupported strategy '" + 48 fs + 49 "'\n" 50 "The only supported strategies are " 51 "'FiniteRotationSmallStrain', " 52 "'MieheApelLambrechtLogarithmicStrain' and " 53 "'MieheApelLambrechtLogarithmicStrainII'."); 54 } // end of checkFiniteStrainStrategy 55 checkFiniteStrainStrategyDefinitionConsistency(const BehaviourDescription & bd,const std::string & fs)56 static void checkFiniteStrainStrategyDefinitionConsistency( 57 const BehaviourDescription& bd, const std::string& fs) { 58 auto throw_if = [](const bool c, const std::string& msg) { 59 tfel::raise_if(c, 60 "checkFiniteStrainStrategyDefinitionConsistency: " + msg); 61 }; 62 checkFiniteStrainStrategy(fs); 63 if (bd.isStrainMeasureDefined()) { 64 const auto ms = bd.getStrainMeasure(); 65 if (ms == BehaviourDescription::LINEARISED) { 66 throw_if(fs != "Native", 67 "incompatible finite strain strategy " 68 "'" + 69 fs + "' (only `Native` accepted)"); 70 } else if (ms == BehaviourDescription::GREENLAGRANGE) { 71 throw_if(fs != "FiniteRotationSmallStrain", 72 "incompatible finite strain strategy " 73 "'" + 74 fs + "' (only `FiniteRotationSmallStrain` accepted)"); 75 } else if (ms == BehaviourDescription::HENCKY) { 76 throw_if(fs != "MieheApelLambrechtLogarithmicStrain", 77 "incompatible finite strain strategy '" + fs + 78 "' " 79 "(only `MieheApelLambrechtLogarithmicStrain` accepted)"); 80 } else { 81 throw_if(true, "unsupported finite strain strategy"); 82 } 83 } 84 } // end of checkFiniteStrainStrategyDefinitionConsistency 85 checkFiniteStrainStrategyDefinitionConsistency(const BehaviourDescription & bd)86 static void checkFiniteStrainStrategyDefinitionConsistency( 87 const BehaviourDescription& bd) { 88 auto throw_if = [](const bool c, const std::string& msg) { 89 tfel::raise_if(c, 90 "checkFiniteStrainStrategyDefinitionConsistency: " + msg); 91 }; 92 if (bd.getBehaviourType() != 93 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) { 94 throw_if(bd.hasAttribute(CalculiXInterface::finiteStrainStrategy), 95 "finite strain strategy is only supported for strain based " 96 "behaviours"); 97 } else { 98 if (bd.hasAttribute(CalculiXInterface::finiteStrainStrategy)) { 99 const auto fs = bd.getAttribute<std::string>( 100 CalculiXInterface::finiteStrainStrategy); 101 checkFiniteStrainStrategyDefinitionConsistency(bd, fs); 102 } 103 } 104 } // end of checkFiniteStrainStrategyDefinitionConsistency 105 hasFiniteStrainStrategy(const BehaviourDescription & bd)106 bool CalculiXInterface::hasFiniteStrainStrategy( 107 const BehaviourDescription& bd) { 108 checkFiniteStrainStrategyDefinitionConsistency(bd); 109 if (bd.isStrainMeasureDefined()) { 110 return bd.getStrainMeasure() != BehaviourDescription::LINEARISED; 111 } 112 return bd.hasAttribute(CalculiXInterface::finiteStrainStrategy); 113 } // end of CalculiXInterface::hasFiniteStrainStrategy 114 getFiniteStrainStrategy(const BehaviourDescription & bd)115 static std::string getFiniteStrainStrategy(const BehaviourDescription& bd) { 116 checkFiniteStrainStrategyDefinitionConsistency(bd); 117 auto throw_if = [](const bool c, const std::string& msg) { 118 tfel::raise_if(c, "getFiniteStrainStrategy: " + msg); 119 }; 120 if (bd.isStrainMeasureDefined()) { 121 const auto ms = bd.getStrainMeasure(); 122 if (ms == BehaviourDescription::GREENLAGRANGE) { 123 return "FiniteRotationSmallStrain"; 124 } else if (ms == BehaviourDescription::HENCKY) { 125 return "MieheApelLambrechtLogarithmicStrain"; 126 } else { 127 throw_if(true, "unsupported strain measure"); 128 } 129 } 130 throw_if(!bd.hasAttribute(CalculiXInterface::finiteStrainStrategy), 131 "no finite strain strategy defined"); 132 return bd.getAttribute<std::string>( 133 CalculiXInterface::finiteStrainStrategy); 134 } // end of getFiniteStrainStrategy 135 writeArguments(std::ostream & out,const BehaviourDescription & mb,const bool base)136 static void writeArguments(std::ostream& out, 137 const BehaviourDescription& mb, 138 const bool base) { 139 if (!base) { 140 const auto requires_strain = [&mb] { 141 if (mb.getBehaviourType() == 142 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) { 143 if (CalculiXInterface::hasFiniteStrainStrategy(mb)) { 144 return getFiniteStrainStrategy(mb) == "FiniteRotationSmallStrain"; 145 } 146 } 147 return true; 148 }(); 149 out << "(const char * const amat,\n" 150 << " const calculix::CalculiXInt* const iel,\n" 151 << " const calculix::CalculiXInt* const iint,\n" 152 << " const calculix::CalculiXInt* const NPROPS,\n" 153 << " const calculix::CalculiXReal* const MPROPS,\n"; 154 if (requires_strain) { 155 out << " const calculix::CalculiXReal* const STRAN1,\n" 156 << " const calculix::CalculiXReal* const STRAN0,\n"; 157 } else { 158 out << " const calculix::CalculiXReal* const,\n" 159 << " const calculix::CalculiXReal* const,\n"; 160 } 161 out << " const calculix::CalculiXReal* const beta,\n" 162 << " const calculix::CalculiXReal* const XOKL,\n" 163 << " const calculix::CalculiXReal* const voj,\n" 164 << " const calculix::CalculiXReal* const XKL,\n" 165 << " const calculix::CalculiXReal* const vj,\n" 166 << " const calculix::CalculiXInt* const ithermal,\n" 167 << " const calculix::CalculiXReal* const TEMP1,\n" 168 << " const calculix::CalculiXReal* const DTIME,\n" 169 << " const calculix::CalculiXReal* const time,\n" 170 << " const calculix::CalculiXReal* const ttime,\n" 171 << " const calculix::CalculiXInt* const icmd,\n" 172 << " const calculix::CalculiXInt* const ielas,\n" 173 << " const calculix::CalculiXInt* const mi,\n" 174 << " const calculix::CalculiXInt* const NSTATV,\n" 175 << " const calculix::CalculiXReal* const STATEV0,\n" 176 << " calculix::CalculiXReal* const STATEV1,\n" 177 << " calculix::CalculiXReal* const STRESS,\n" 178 << " calculix::CalculiXReal* const DDSDDE,\n" 179 << " const calculix::CalculiXInt* const iorien,\n" 180 << " const calculix::CalculiXReal* const pgauss,\n" 181 << " const calculix::CalculiXReal* const orab,\n" 182 << " calculix::CalculiXReal* const PNEWDT,\n" 183 << " const calculix::CalculiXInt* const ipkon,\n" 184 << " const int size)"; 185 } else { 186 const auto requires_strain = 187 (mb.getBehaviourType() == 188 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR); 189 const auto requires_F = 190 (mb.getBehaviourType() == 191 BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR); 192 out << "(const char * const,\n" 193 << " const calculix::CalculiXInt* const iel,\n" 194 << " const calculix::CalculiXInt* const iint,\n" 195 << " const calculix::CalculiXInt* const,\n" // NPROPS 196 << " const calculix::CalculiXReal* const MPROPS,\n"; 197 if (requires_strain) { 198 out << " const calculix::CalculiXReal *const DSTRAN,\n" 199 << " const calculix::CalculiXReal *const STRAN0,\n"; 200 } else { 201 out << " const calculix::CalculiXReal *const,\n" 202 << " const calculix::CalculiXReal *const,\n"; 203 } 204 out << " const calculix::CalculiXReal* const,\n"; 205 if (requires_F) { 206 out << " const calculix::CalculiXReal* const XOKL,\n" 207 << " const calculix::CalculiXReal* const ,\n" 208 << " const calculix::CalculiXReal* const XKL,\n" 209 << " const calculix::CalculiXReal* const,\n"; 210 } else { 211 out << " const calculix::CalculiXReal* const,\n" 212 << " const calculix::CalculiXReal* const ,\n" 213 << " const calculix::CalculiXReal* const,\n" 214 << " const calculix::CalculiXReal* const,\n"; 215 } 216 out << " const calculix::CalculiXInt* const,\n" // ithermal 217 << " const calculix::CalculiXReal* const TEMP1,\n" 218 << " const calculix::CalculiXReal* const DTIME,\n" 219 << " const calculix::CalculiXReal* const,\n" 220 << " const calculix::CalculiXReal* const,\n" 221 << " const calculix::CalculiXInt* const,\n" // icmd 222 << " const calculix::CalculiXInt* const,\n" // ielas 223 << " const calculix::CalculiXInt* const mi,\n" 224 << " const calculix::CalculiXInt* const NSTATV,\n" 225 << " const calculix::CalculiXReal* const STATEV0,\n" 226 << " calculix::CalculiXReal* const STATEV1,\n" 227 << " calculix::CalculiXReal* const STRESS,\n" 228 << " calculix::CalculiXReal* const DDSDDE,\n" 229 << " const calculix::CalculiXInt* const,\n" 230 << " const calculix::CalculiXReal* const,\n" 231 << " const calculix::CalculiXReal* const,\n" 232 << " calculix::CalculiXReal* const PNEWDT,\n" 233 << " const calculix::CalculiXInt* const,\n" 234 << " const int)"; 235 } 236 } // end of writeArguments 237 writeArguments(std::ostream & out)238 static void writeArguments(std::ostream& out) { 239 out << "(const char * const,\n" 240 << " const calculix::CalculiXInt* const,\n" 241 << " const calculix::CalculiXInt* const,\n" 242 << " const calculix::CalculiXInt* const,\n" 243 << " const calculix::CalculiXReal* const,\n" 244 << " const calculix::CalculiXReal* const,\n" 245 << " const calculix::CalculiXReal* const,\n" 246 << " const calculix::CalculiXReal* const,\n" 247 << " const calculix::CalculiXReal* const,\n" 248 << " const calculix::CalculiXReal* const,\n" 249 << " const calculix::CalculiXReal* const,\n" 250 << " const calculix::CalculiXReal* const,\n" 251 << " const calculix::CalculiXInt* const,\n" 252 << " const calculix::CalculiXReal* const,\n" 253 << " const calculix::CalculiXReal* const,\n" 254 << " const calculix::CalculiXReal* const,\n" 255 << " const calculix::CalculiXReal* const,\n" 256 << " const calculix::CalculiXInt* const,\n" 257 << " const calculix::CalculiXInt* const,\n" 258 << " const calculix::CalculiXInt* const,\n" 259 << " const calculix::CalculiXInt* const,\n" 260 << " const calculix::CalculiXReal* const,\n" 261 << " calculix::CalculiXReal* const,\n" 262 << " calculix::CalculiXReal* const,\n" 263 << " calculix::CalculiXReal* const,\n" 264 << " const calculix::CalculiXInt* const,\n" 265 << " const calculix::CalculiXReal* const,\n" 266 << " const calculix::CalculiXReal* const,\n" 267 << " calculix::CalculiXReal* const,\n" 268 << " const calculix::CalculiXInt* const,\n" 269 << " const int)"; 270 } // end of writeArguments 271 272 const char* const CalculiXInterface::finiteStrainStrategy = 273 "calculix::finiteStrainStrategy"; 274 getName()275 std::string CalculiXInterface::getName() { return "calculix"; } 276 getInterfaceName() const277 std::string CalculiXInterface::getInterfaceName() const { 278 return "CalculiX"; 279 } // end of CalculiXInterface::getInterfaceName 280 281 std::pair<bool, CalculiXInterface::tokens_iterator> treatKeyword(BehaviourDescription & bd,const std::string & k,const std::vector<std::string> & i,tokens_iterator current,const tokens_iterator end)282 CalculiXInterface::treatKeyword(BehaviourDescription& bd, 283 const std::string& k, 284 const std::vector<std::string>& i, 285 tokens_iterator current, 286 const tokens_iterator end) { 287 using tfel::utilities::CxxTokenizer; 288 auto throw_if = [](const bool b, const std::string& m) { 289 tfel::raise_if(b, "CalculiXInterface::treatKeyword: " + m); 290 }; 291 if (!i.empty()) { 292 if (std::find(i.begin(), i.end(), this->getName()) != i.end()) { 293 const auto keys = 294 std::vector<std::string>{"@CalculiXFiniteStrainStrategy", 295 "@CalculiXGenerateMTestFileOnFailure", 296 "@CalculiXStrainPerturbationValue"}; 297 throw_if(std::find(keys.begin(), keys.end(), k) == keys.end(), 298 "unsupported key '" + k + "'"); 299 } else { 300 return {false, current}; 301 } 302 } 303 if (k == "@CalculiXFiniteStrainStrategy") { 304 throw_if(bd.hasAttribute(CalculiXInterface::finiteStrainStrategy), 305 "a finite strain strategy has already been defined"); 306 throw_if(current == end, "unexpected end of file"); 307 const auto fs = current->value; 308 throw_if(++current == end, "unexpected end of file"); 309 throw_if(current->value != ";", 310 "expected ';', read '" + current->value + '\''); 311 ++(current); 312 checkFiniteStrainStrategyDefinitionConsistency(bd, fs); 313 bd.setAttribute(CalculiXInterface::finiteStrainStrategy, fs, false); 314 return {true, current}; 315 } else if (k == "@CalculiXGenerateMTestFileOnFailure") { 316 this->setGenerateMTestFileOnFailureAttribute( 317 bd, this->readBooleanValue(k, current, end)); 318 return {true, current}; 319 } 320 return {false, current}; 321 } // end of CalculiXInterface::treatKeyword 322 endTreatment(const BehaviourDescription & mb,const FileDescription & fd) const323 void CalculiXInterface::endTreatment(const BehaviourDescription& mb, 324 const FileDescription& fd) const { 325 using namespace tfel::system; 326 auto throw_if = [](const bool b, const std::string& m) { 327 tfel::raise_if(b, "CalculiXInterface::endTreatment: " + m); 328 }; 329 throw_if(!((mb.getBehaviourType() == 330 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) || 331 (mb.getBehaviourType() == 332 BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR)), 333 "the calculix interface only supports small and " 334 "finite strain behaviours"); 335 checkFiniteStrainStrategyDefinitionConsistency(mb); 336 // the only supported modelling hypothesis 337 constexpr const auto h = ModellingHypothesis::TRIDIMENSIONAL; 338 const auto& d = mb.getBehaviourData(h); 339 throw_if(d.getExternalStateVariables().size() != 1u, 340 "external state variables are not supported " 341 "by CalculiX's native interface"); 342 // get the modelling hypotheses to be treated 343 const auto name = mb.getLibrary() + mb.getClassName(); 344 // output directories 345 systemCall::mkdir("include/MFront"); 346 systemCall::mkdir("include/MFront/CalculiX"); 347 systemCall::mkdir("calculix"); 348 349 std::ofstream out; 350 351 // header 352 auto fname = "calculix" + name + ".hxx"; 353 out.open("include/MFront/CalculiX/" + fname); 354 throw_if(!out, "could not open file '" + fname + "'"); 355 356 out << "/*!\n" 357 << "* \\file " << fname << '\n' 358 << "* \\brief This file declares the calculix interface for the " 359 << mb.getClassName() << " behaviour law\n" 360 << "* \\author " << fd.authorName << '\n' 361 << "* \\date " << fd.date << '\n' 362 << "*/\n\n"; 363 364 const auto header = this->getHeaderGuard(mb); 365 out << "#ifndef " << header << "\n" 366 << "#define " << header << "\n\n" 367 << "#include\"TFEL/Config/TFELConfig.hxx\"\n\n"; 368 if ((mb.getBehaviourType() == 369 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) && 370 (CalculiXInterface::hasFiniteStrainStrategy(mb))) { 371 const auto fs = getFiniteStrainStrategy(mb); 372 if ((fs == "MieheApelLambrechtLogarithmicStrain") || 373 (fs == "MieheApelLambrechtLogarithmicStrainII")) { 374 out << "#include\"TFEL/Material/LogarithmicStrainHandler.hxx\"\n\n"; 375 } 376 } 377 out << "#include\"MFront/CalculiX/CalculiX.hxx\"\n" 378 << "#include\"MFront/CalculiX/CalculiXData.hxx\"\n\n" 379 << "#ifdef __cplusplus\n" 380 << "#include\"MFront/CalculiX/CalculiXTraits.hxx\"\n" 381 << "#include\"TFEL/Material/" << mb.getClassName() << ".hxx\"\n\n" 382 << "#endif /* __cplusplus */\n\n"; 383 384 this->writeVisibilityDefines(out); 385 386 out << "#ifdef __cplusplus\n\n" 387 << "namespace calculix{\n\n"; 388 389 this->writeCalculiXBehaviourTraits(out, mb); 390 391 out << "} // end of namespace calculix\n\n" 392 << "#endif /* __cplusplus */\n\n" 393 << "#ifdef __cplusplus\n" 394 << "extern \"C\"{\n" 395 << "#endif /* __cplusplus */\n\n"; 396 397 this->writeSetOutOfBoundsPolicyFunctionDeclaration(out, name); 398 this->writeSetParametersFunctionsDeclarations(out, mb, name); 399 400 out << "MFRONT_SHAREDOBJ void\n" << this->getFunctionNameBasis(name); 401 writeArguments(out); 402 out << ";\n\n"; 403 404 out << "#ifdef __cplusplus\n" 405 << "}\n" 406 << "#endif /* __cplusplus */\n\n" 407 << "#endif /* " << header << " */\n"; 408 409 out.close(); 410 411 fname = "calculix" + name + ".cxx"; 412 out.open("src/" + fname); 413 throw_if(!out, "could not open file '" + fname + "'"); 414 415 out << "/*!\n" 416 << "* \\file " << fname << '\n' 417 << "* \\brief This file implements the calculix interface for the " 418 << mb.getClassName() << " behaviour law\n" 419 << "* \\author " << fd.authorName << '\n' 420 << "* \\date " << fd.date << '\n' 421 << "*/\n\n"; 422 423 this->getExtraSrcIncludes(out, mb); 424 425 out << "#include\"TFEL/Material/OutOfBoundsPolicy.hxx\"\n" 426 << "#include\"TFEL/Material/" << mb.getClassName() << ".hxx\"\n"; 427 if (mb.getAttribute(BehaviourData::profiling, false)) { 428 out << "#include\"MFront/BehaviourProfiler.hxx\"\n\n"; 429 } 430 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 431 out << "#include\"MFront/CalculiX/CalculiXRotationMatrix.hxx\"\n"; 432 } 433 out << "#include\"MFront/CalculiX/" 434 "CalculiXStressFreeExpansionHandler.hxx\"\n\n" 435 << "#include\"MFront/CalculiX/CalculiXInterface.hxx\"\n\n" 436 << "#include\"MFront/CalculiX/calculix" << name << ".hxx\"\n\n"; 437 438 this->writeGetOutOfBoundsPolicyFunctionImplementation(out, name); 439 440 out << "extern \"C\"{\n\n"; 441 442 CalculiXSymbolsGenerator sg; 443 sg.generateGeneralSymbols(out, *this, mb, fd, {h}, name); 444 sg.generateSymbols(out, *this, mb, fd, name, h); 445 446 this->writeSetParametersFunctionsImplementations(out, mb, name); 447 this->writeSetOutOfBoundsPolicyFunctionImplementation(out, name); 448 449 if (mb.getBehaviourType() == 450 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) { 451 if (CalculiXInterface::hasFiniteStrainStrategy(mb)) { 452 const auto fs = getFiniteStrainStrategy(mb); 453 if (fs == "FiniteRotationSmallStrain") { 454 this->writeFiniteRotationSmallStrainFunction(out, mb, name); 455 } else if ((fs == "MieheApelLambrechtLogarithmicStrain") || 456 (fs == "MieheApelLambrechtLogarithmicStrainII")) { 457 this->writeMieheApelLambrechtLogarithmicStrainFunction(out, mb, name); 458 } else { 459 throw_if(true, "unsupported finite strain strategy !"); 460 } 461 } else { 462 this->writeSmallStrainFunction(out, mb, name); 463 } 464 } else if (mb.getBehaviourType() == 465 BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR) { 466 this->writeFiniteStrainFunction(out, mb, name); 467 } else { 468 throw_if(true, 469 "the calculix interface only supports small " 470 "and finite strain behaviours"); 471 } 472 out << "} // end of extern \"C\"\n"; 473 out.close(); 474 this->writeInputFileExample(mb, fd, true); 475 } // end of CalculiXInterface::endTreatment 476 writeFunctionBase(std::ostream & out,const BehaviourDescription & mb,const std::string & name,const std::string & sfeh) const477 void CalculiXInterface::writeFunctionBase(std::ostream& out, 478 const BehaviourDescription& mb, 479 const std::string& name, 480 const std::string& sfeh) const { 481 auto throw_if = [](const bool b, const std::string& m) { 482 tfel::raise_if(b, "CalculiXInterface::writeFunctionBase: " + m); 483 }; 484 std::string dv0, dv1, sig, statev, nstatev; 485 const auto btype = mb.getBehaviourType(); 486 out << "static void\n" << name << "_base"; 487 writeArguments(out, mb, true); 488 out << "{\n"; 489 if (btype == BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) { 490 dv0 = "STRAN0"; 491 dv1 = "DSTRAN"; 492 } else if (btype == BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR) { 493 dv0 = "XOKL"; 494 dv1 = "XKL"; 495 } else { 496 throw_if(true, 497 "the calculix interface only supports small " 498 "and finite strain behaviours"); 499 } 500 out << "using calculix::CalculiXData;\n" 501 << "const auto ivs0 = STATEV0+(*NSTATV)*((*iint-1)+(*mi)*(*iel-1));\n" 502 << "const auto ivs1 = STATEV1+(*NSTATV)*((*iint-1)+(*mi)*(*iel-1));\n" 503 << "CalculiXData d = {STRESS,PNEWDT,DDSDDE,ivs1,*DTIME,ivs0,\n" 504 << " " << dv0 << "," << dv1 << ",TEMP1,MPROPS,\n" 505 << this->getFunctionNameBasis(name) << "_getOutOfBoundsPolicy()," 506 << sfeh << "};\n" 507 << "if(calculix::CalculiXInterface<tfel::material::" 508 << mb.getClassName() << ">::exe(d)!=0){\n" 509 << "*PNEWDT = 0.2;\n" 510 << "return;\n" 511 << "}\n" 512 << "}\n\n"; 513 } // end of CalculiXInterface::writeFunctionBase 514 writeFiniteStrainFunction(std::ostream & out,const BehaviourDescription & mb,const std::string & name) const515 void CalculiXInterface::writeFiniteStrainFunction( 516 std::ostream& out, 517 const BehaviourDescription& mb, 518 const std::string& name) const { 519 auto throw_if = [](const bool b, const std::string& m) { 520 tfel::raise_if(b, "CalculiXInterface::writeFiniteStrainFunction: " + m); 521 }; 522 const std::string sfeh = "nullptr"; 523 this->writeFunctionBase(out, mb, name, sfeh); 524 out << "MFRONT_SHAREDOBJ void\n" << this->getFunctionNameBasis(name); 525 writeArguments(out, mb, false); 526 out << "{\n" 527 << "using namespace tfel::math;\n" 528 << "using real = calculix::CalculiXReal;\n"; 529 if (mb.getAttribute(BehaviourData::profiling, false)) { 530 out << "using mfront::BehaviourProfiler;\n" 531 << "using tfel::material::" << mb.getClassName() << "Profiler;\n" 532 << "BehaviourProfiler::Timer total_timer(" << mb.getClassName() 533 << "Profiler::getProfiler(),\n" 534 << "BehaviourProfiler::TOTALTIME);\n"; 535 } 536 out << "st2tost2<3u,real> D = {0,0,0,0,0,0,\n" 537 << " 0,0,0,0,0,0,\n" 538 << " 0,0,0,0,0,0,\n" 539 << " 0,0,0,0,0,0,\n" 540 << " 0,0,0,0,0,0,\n" 541 << " 0,0,0,0,0,0};\n"; 542 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 543 out << "if(*iorien==0){\n" 544 << " std::cerr << \"" << this->getFunctionNameBasis(name) << ":\"\n" 545 << " << \"no orientation defined for an orthotropic " 546 "behaviour\\n\";\n" 547 << " std::exit(-1);\n" 548 << "}\n" 549 << "const auto r = " 550 "calculix::getRotationMatrix(orab+7*(*iorien-1),pgauss);\n" 551 << "const auto rb = transpose(r);\n"; 552 } else { 553 throw_if(mb.getSymmetryType() != mfront::ISOTROPIC, 554 "unsupported symmetry type"); 555 out << "if(*iorien!=0){\n" 556 << " std::cerr << \"" << this->getFunctionNameBasis(name) << ":\"\n" 557 << " << \"no orientation shall be defined for an istropic " 558 "behaviour\\n\";\n" 559 << " std::exit(-1);\n" 560 << "}\n"; 561 } 562 out << "const auto F0 = tensor<3u,real>::buildFromFortranMatrix(XOKL);\n" 563 << "const auto F1 = tensor<3u,real>::buildFromFortranMatrix(XKL);\n" 564 << "stensor<3u,real> pk2;\n" 565 << "pk2.importTab(STRESS);\n" 566 << "auto s = " 567 "convertSecondPiolaKirchhoffStressToCauchyStress(pk2,F1);\n"; 568 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 569 out << "const auto rF0 = change_basis(F0,r);\n" 570 << "const auto rF1 = change_basis(F1,r);\n" 571 << "s.changeBasis(r);\n" 572 << name << "_base" 573 << "(amat,iel,iint,NPROPS,MPROPS,STRAN1,STRAN0,beta,rF0.begin()," 574 << " voj,rF1.begin(),vj,ithermal,TEMP1,DTIME,time,ttime,icmd," 575 << " ielas,mi,NSTATV,STATEV0,STATEV1,s.begin(),D.begin()," 576 << "iorien,pgauss,orab,PNEWDT,ipkon,size);\n"; 577 } else { 578 out << name << "_base" 579 << "(amat,iel,iint,NPROPS,MPROPS,STRAN1,STRAN0,beta,F0.begin()," 580 << " voj,F1.begin(),vj,ithermal,TEMP1,DTIME,time,ttime,icmd," 581 << " ielas,mi,NSTATV,STATEV0,STATEV1,s.begin(),D.begin()," 582 << "iorien,pgauss,orab,PNEWDT,ipkon,size);\n"; 583 } 584 out << "if(*PNEWDT>=1){\n" 585 << "*PNEWDT=-1;\n"; 586 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 587 out << "s.changeBasis(rb);\n" 588 << "D = change_basis(st2tost2<3u,real>(D),rb);\n"; 589 } 590 out << "// turning Cauchy stress to pk2\n" 591 << "pk2 = convertCauchyStressToSecondPiolaKirchhoffStress(s,F1);\n" 592 << "pk2.exportTab(STRESS);\n" 593 << "// converting the consistent tangent operator\n" 594 << "calculix::ConvertUnsymmetricTangentOperator::exe(DDSDDE,D.begin());" 595 "\n"; 596 if (getDebugMode()) { 597 out << "std::cout << \"Dt :\" << std::endl;\n" 598 << "const calculix::CalculiXReal *p = DDSDDE;\n" 599 << "for(calculix::CalculiXInt i=0;i!=6;++i){\n" 600 << "for(calculix::CalculiXInt j=0;j!=i+1;++j,++p){\n" 601 << "std::cout << *p << \" \";\n" 602 << "}\n" 603 << "std::cout << std::endl;\n" 604 << "}\n" 605 << "std::cout << std::endl;\n"; 606 } 607 out << "}\n" 608 << "} // end of " << this->getFunctionNameBasis(name) << "\n\n"; 609 } 610 writeSmallStrainFunction(std::ostream & out,const BehaviourDescription & mb,const std::string & name) const611 void CalculiXInterface::writeSmallStrainFunction( 612 std::ostream& out, 613 const BehaviourDescription& mb, 614 const std::string& name) const { 615 auto throw_if = [](const bool b, const std::string& m) { 616 tfel::raise_if(b, "CalculiXInterface::writeSmallStrainFunction: " + m); 617 }; 618 const std::string sfeh = 619 "calculix::CalculiXStandardSmallStrainStressFreeExpansionHandler"; 620 this->writeFunctionBase(out, mb, name, sfeh); 621 out << "MFRONT_SHAREDOBJ void\n" << this->getFunctionNameBasis(name); 622 writeArguments(out, mb, false); 623 out << "{\n" 624 << "using namespace tfel::math;\n" 625 << "using real = calculix::CalculiXReal;\n" 626 << "TFEL_CONSTEXPR const real sqrt2 = Cste<real>::sqrt2;\n"; 627 if (mb.getAttribute(BehaviourData::profiling, false)) { 628 out << "using mfront::BehaviourProfiler;\n" 629 << "using tfel::material::" << mb.getClassName() << "Profiler;\n" 630 << "BehaviourProfiler::Timer total_timer(" << mb.getClassName() 631 << "Profiler::getProfiler(),\n" 632 << "BehaviourProfiler::TOTALTIME);\n"; 633 } 634 if (this->shallGenerateMTestFileOnFailure(mb)) { 635 this->generateMTestFile1(out, mb); 636 } 637 out << "st2tost2<3u,real> D = {0,0,0,0,0,0,\n" 638 << " 0,0,0,0,0,0,\n" 639 << " 0,0,0,0,0,0,\n" 640 << " 0,0,0,0,0,0,\n" 641 << " 0,0,0,0,0,0,\n" 642 << " 0,0,0,0,0,0};\n"; 643 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 644 out << "stensor<3u,real> eto = {STRAN0[0],\n" 645 << " STRAN0[1],\n" 646 << " STRAN0[2],\n" 647 << " STRAN0[3]*sqrt2,\n" 648 << " STRAN0[4]*sqrt2,\n" 649 << " STRAN0[5]*sqrt2};\n" 650 << "stensor<3u,real> deto = {STRAN1[0]-STRAN0[0],\n" 651 << " STRAN1[1]-STRAN0[1],\n" 652 << " STRAN1[2]-STRAN0[2],\n" 653 << " (STRAN1[3]-STRAN0[3])*sqrt2,\n" 654 << " (STRAN1[4]-STRAN0[4])*sqrt2,\n" 655 << " (STRAN1[5]-STRAN0[5])*sqrt2};\n" 656 << "stensor<3u,real> s;\n" 657 << "s.importTab(STRESS);\n" 658 << "if(*iorien==0){\n" 659 << " std::cerr << \"" << this->getFunctionNameBasis(name) << ":\"\n" 660 << " << \"no orientation defined for an orthotropic " 661 "behaviour\\n\";\n" 662 << " std::exit(-1);\n" 663 << "}\n" 664 << "const auto r = " 665 "calculix::getRotationMatrix(orab+7*(*iorien-1),pgauss);\n" 666 << "eto.changeBasis(r);\n" 667 << "deto.changeBasis(r);\n" 668 << "s.changeBasis(r);\n"; 669 } else { 670 throw_if(mb.getSymmetryType() != mfront::ISOTROPIC, 671 "unsupported symmetry type"); 672 out << "const stensor<3u,real> eto = {STRAN0[0],\n" 673 << " STRAN0[1],\n" 674 << " STRAN0[2],\n" 675 << " STRAN0[3]*sqrt2,\n" 676 << " STRAN0[4]*sqrt2,\n" 677 << " STRAN0[5]*sqrt2};\n" 678 << "const stensor<3u,real> deto = {STRAN1[0]-STRAN0[0],\n" 679 << " STRAN1[1]-STRAN0[1],\n" 680 << " STRAN1[2]-STRAN0[2],\n" 681 << " (STRAN1[3]-STRAN0[3])*sqrt2,\n" 682 << " (STRAN1[4]-STRAN0[4])*sqrt2,\n" 683 << " (STRAN1[5]-STRAN0[5])*sqrt2};\n" 684 << "stensor<3u,real> s;\n" 685 << "s.importTab(STRESS);\n" 686 << "if(*iorien!=0){\n" 687 << " std::cerr << \"" << this->getFunctionNameBasis(name) << ":\"\n" 688 << " << \"no orientation shall be defined for an istropic " 689 "behaviour\\n\";\n" 690 << " std::exit(-1);\n" 691 << "}\n"; 692 } 693 out << name << "_base" 694 << "(amat,iel,iint,NPROPS,MPROPS,deto.begin(),eto.begin(),beta,XOKL," 695 << " voj,XKL,vj,ithermal,TEMP1,DTIME,time,ttime,icmd," 696 << " ielas,mi,NSTATV,STATEV0,STATEV1,s.begin(),D.begin()," 697 << "iorien,pgauss,orab,PNEWDT,ipkon,size);\n" 698 << "if(*PNEWDT>=1){\n" 699 << "*PNEWDT=-1;\n"; 700 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 701 out << "const auto rb = transpose(r);\n" 702 << "s.changeBasis(rb);\n" 703 << "D = change_basis(st2tost2<3u,real>(D),rb);\n"; 704 } 705 out << "s.exportTab(STRESS);\n" 706 << "// converting the consistent tangent operator\n" 707 << "calculix::ConvertUnsymmetricTangentOperator::exe(DDSDDE,D.begin());" 708 "\n"; 709 if (getDebugMode()) { 710 out << "std::cout << \"Dt :\" << std::endl;\n" 711 << "const calculix::CalculiXReal *p = DDSDDE;\n" 712 << "for(calculix::CalculiXInt i=0;i!=6;++i){\n" 713 << "for(calculix::CalculiXInt j=0;j!=i+1;++j,++p){\n" 714 << "std::cout << *p << \" \";\n" 715 << "}\n" 716 << "std::cout << std::endl;\n" 717 << "}\n" 718 << "std::cout << std::endl;\n"; 719 } 720 if (this->shallGenerateMTestFileOnFailure(mb)) { 721 out << "} else {\n"; 722 this->generateMTestFile2(out, mb, mb.getBehaviourType(), name, ""); 723 } 724 out << "}\n"; 725 out << "} // end of " << this->getFunctionNameBasis(name) << "\n\n"; 726 } 727 writeFiniteRotationSmallStrainFunction(std::ostream & out,const BehaviourDescription & mb,const std::string & name) const728 void CalculiXInterface::writeFiniteRotationSmallStrainFunction( 729 std::ostream& out, 730 const BehaviourDescription& mb, 731 const std::string& name) const { 732 this->writeSmallStrainFunction(out, mb, name); 733 } // end of CalculiXInterface::writeFiniteRotationSmallStrainFunction 734 writeMieheApelLambrechtLogarithmicStrainFunction(std::ostream & out,const BehaviourDescription & mb,const std::string & name) const735 void CalculiXInterface::writeMieheApelLambrechtLogarithmicStrainFunction( 736 std::ostream& out, 737 const BehaviourDescription& mb, 738 const std::string& name) const { 739 auto throw_if = [](const bool b, const std::string& m) { 740 tfel::raise_if(b, 741 "CalculiXInterface::writeMieheApelLambrecht" 742 "LogarithmicStrainFunction: " + 743 m); 744 }; 745 constexpr const auto h = ModellingHypothesis::TRIDIMENSIONAL; 746 const auto& ivs = mb.getBehaviourData(h).getPersistentVariables(); 747 const auto nivs = ivs.getTypeSize().getValueForModellingHypothesis(h); 748 const std::string sfeh = 749 "calculix::CalculiXLogarithmicStrainStressFreeExpansionHandler"; 750 const auto variant = [&mb, throw_if] { 751 const auto fs = getFiniteStrainStrategy(mb); 752 throw_if((fs != "MieheApelLambrechtLogarithmicStrain") && 753 (fs != "MieheApelLambrechtLogarithmicStrainII"), 754 "invalid finite strain strategy (internal error)"); 755 return fs == "MieheApelLambrechtLogarithmicStrain"; 756 }(); 757 this->writeFunctionBase(out, mb, name, sfeh); 758 out << "MFRONT_SHAREDOBJ void\n" << this->getFunctionNameBasis(name); 759 writeArguments(out, mb, false); 760 out << "{\n" 761 << "using namespace tfel::math;\n" 762 << "using namespace tfel::material;\n" 763 << "using real = calculix::CalculiXReal;\n"; 764 if (this->shallGenerateMTestFileOnFailure(mb)) { 765 this->generateMTestFile1(out, mb); 766 } 767 out << "st2tost2<3u,real> D = {0,0,0,0,0,0,\n" 768 << " 0,0,0,0,0,0,\n" 769 << " 0,0,0,0,0,0,\n" 770 << " 0,0,0,0,0,0,\n" 771 << " 0,0,0,0,0,0,\n" 772 << " 0,0,0,0,0,0};\n"; 773 if (variant) { 774 out << "LogarithmicStrainHandler<3u,real> " 775 << "lsh0(LogarithmicStrainHandlerBase::LAGRANGIAN,\n" 776 << " tensor<3u,real>::buildFromFortranMatrix(XOKL));\n" 777 << "LogarithmicStrainHandler<3u,real> " 778 << "lsh1(LogarithmicStrainHandlerBase::LAGRANGIAN,\n" 779 << " tensor<3u,real>::buildFromFortranMatrix(XKL));\n" 780 << "auto eto0 = lsh0.getHenckyLogarithmicStrain();\n" 781 << "tfel::math::stensor<3u,real> pk2;\n" 782 << "pk2.importTab(STRESS);\n" 783 << "auto T = lsh0.convertFromSecondPiolaKirchhoffStress(pk2);\n"; 784 } else { 785 out << "LogarithmicStrainHandler<3u,real> " 786 << "lsh1(LogarithmicStrainHandlerBase::LAGRANGIAN,\n" 787 << " tensor<3u,real>::buildFromFortranMatrix(XKL));\n" 788 << "const auto ivs0 = " 789 "STATEV0+(*NSTATV)*((*iint-1)+(*mi)*(*iel-1));\n" 790 << "const auto ivs1 = " 791 "STATEV1+(*NSTATV)*((*iint-1)+(*mi)*(*iel-1));\n" 792 << "stensor<3u,real> eto0(ivs0+" << nivs << ");\n" 793 << "stensor<3u,real> T;\n" 794 << "tfel::fsalgo::copy<6u>::exe(ivs1+" << nivs + 6 795 << ",T.begin());\n"; 796 } 797 out << "auto eto1 = lsh1.getHenckyLogarithmicStrain();\n"; 798 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 799 out << "if(*iorien==0){\n" 800 << " std::cerr << \"" << this->getFunctionNameBasis(name) << ":\"\n" 801 << " << \"no orientation defined for an orthotropic " 802 "behaviour\\n\";\n" 803 << " std::exit(-1);\n" 804 << "}\n" 805 << "const auto r = " 806 "calculix::getRotationMatrix(orab+7*(*iorien-1),pgauss);\n" 807 << "const auto rb = transpose(r);\n"; 808 if (variant) { 809 out << "eto0.changeBasis(r);\n"; 810 } 811 out << "eto1.changeBasis(r);\n"; 812 if (variant) { 813 out << "T.changeBasis(r);\n"; 814 } 815 } else { 816 throw_if(mb.getSymmetryType() != mfront::ISOTROPIC, 817 "unsupported symmetry type"); 818 out << "if(*iorien!=0){\n" 819 << " std::cerr << \"" << this->getFunctionNameBasis(name) << ":\"\n" 820 << " << \"no orientation shall be defined for an istropic " 821 "behaviour\\n\";\n" 822 << " std::exit(-1);\n" 823 << "}\n"; 824 } 825 out << "auto deto = eval(eto1-eto0);\n" 826 << "// behaviour integration\n" 827 << name << "_base" 828 << "(amat,iel,iint,NPROPS,MPROPS,deto.begin(),eto0.begin(),beta,XOKL," 829 << " voj,XKL,vj,ithermal,TEMP1,DTIME,time,ttime,icmd," 830 << " ielas,mi,NSTATV,STATEV0,STATEV1,T.begin(),D.begin()," 831 << "iorien,pgauss,orab,PNEWDT,ipkon,size);\n" 832 << "if(*PNEWDT>=1){\n" 833 << "*PNEWDT=-1;\n"; 834 if (!variant) { 835 // saving the stresses in the material frame 836 out << "tfel::fsalgo::copy<6u>::exe(eto1.begin(),ivs1+" << nivs << ");\n" 837 << "tfel::fsalgo::copy<6u>::exe(T.begin(),ivs1+" << nivs + 6 838 << ");\n"; 839 } 840 out << "// stress at the end of the time step\n"; 841 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 842 out << "auto s = lsh1.convertToSecondPiolaKirchhoffStress(T);\n" 843 << "s.changeBasis(rb);\n"; 844 } else { 845 out << "const auto s = lsh1.convertToSecondPiolaKirchhoffStress(T);\n"; 846 } 847 out << "D = lsh1.convertToMaterialTangentModuli(st2tost2<3u,real>(D),T);\n"; 848 if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 849 out << "D = change_basis(st2tost2<3u,real>(D),rb);\n"; 850 } 851 out << "// converting the stress\n" 852 << "s.exportTab(STRESS);\n" 853 << "// converting the consistent tangent operator\n" 854 << "calculix::ConvertUnsymmetricTangentOperator::exe(DDSDDE,D.begin());" 855 "\n"; 856 if (getDebugMode()) { 857 out << "std::cout << \"Dt :\" << std::endl;\n" 858 << "const calculix::CalculiXReal *p = DDSDDE;\n" 859 << "for(calculix::CalculiXInt i=0;i!=6;++i){\n" 860 << "for(calculix::CalculiXInt j=0;j!=i+1;++j,++p){\n" 861 << "std::cout << *p << \" \";\n" 862 << "}\n" 863 << "std::cout << std::endl;\n" 864 << "}\n" 865 << "std::cout << std::endl;\n"; 866 } 867 if (this->shallGenerateMTestFileOnFailure(mb)) { 868 out << "} else {\n"; 869 this->generateMTestFile2(out, mb, mb.getBehaviourType(), name, ""); 870 } 871 out << "}\n" 872 << "} // end of " << this->getFunctionNameBasis(name) << "\n\n"; 873 } // end of 874 // CalculiXInterface::writeMieheApelLambrechtLogarithmicStrainFunction 875 writeInterfaceSpecificIncludes(std::ostream & out,const BehaviourDescription &) const876 void CalculiXInterface::writeInterfaceSpecificIncludes( 877 std::ostream& out, const BehaviourDescription&) const { 878 out << "#include\"MFront/CalculiX/CalculiX.hxx\"\n" 879 << "#include\"MFront/CalculiX/CalculiXConvert.hxx\"\n\n"; 880 } // end of CalculiXInterface::writeInterfaceSpecificIncludes 881 writeBehaviourDataGradientSetter(std::ostream & os,const Gradient & v,const SupportedTypes::TypeSize o) const882 void CalculiXInterface::writeBehaviourDataGradientSetter( 883 std::ostream& os, 884 const Gradient& v, 885 const SupportedTypes::TypeSize o) const { 886 const auto iprefix = makeUpperCase(this->getInterfaceName()); 887 tfel::raise_if(!o.isNull(), 888 "CalculiXInterface::writeBehaviourDataMainVariablesSetter : " 889 "only one driving variable supported"); 890 if (Gradient::isIncrementKnown(v)) { 891 os << "calculix::ImportGradients::exe(this->" << v.name << "," 892 << iprefix << "stran);\n"; 893 } else { 894 os << "calculix::ImportGradients::exe(this->" << v.name << "0," 895 << iprefix << "stran);\n"; 896 } 897 } // end of CalculiXInterface::writeBehaviourDataGradientSetter 898 writeIntegrationDataGradientSetter(std::ostream & os,const Gradient & v,const SupportedTypes::TypeSize o) const899 void CalculiXInterface::writeIntegrationDataGradientSetter( 900 std::ostream& os, 901 const Gradient& v, 902 const SupportedTypes::TypeSize o) const { 903 const auto iprefix = makeUpperCase(this->getInterfaceName()); 904 tfel::raise_if( 905 !o.isNull(), 906 "CalculiXInterface::writeIntegrationDataMainVariablesSetter : " 907 "only one driving variable supported"); 908 if (Gradient::isIncrementKnown(v)) { 909 os << "calculix::ImportGradients::exe(this->d" << v.name << "," 910 << iprefix << "dstran);\n"; 911 } else { 912 os << "calculix::ImportGradients::exe(this->" << v.name << "1," 913 << iprefix << "dstran);\n"; 914 } 915 } // end of CalculiXInterface::writeIntegrationDataGradientSetter 916 writeBehaviourDataThermodynamicForceSetter(std::ostream & os,const ThermodynamicForce & f,const SupportedTypes::TypeSize o) const917 void CalculiXInterface::writeBehaviourDataThermodynamicForceSetter( 918 std::ostream& os, 919 const ThermodynamicForce& f, 920 const SupportedTypes::TypeSize o) const { 921 const auto iprefix = makeUpperCase(this->getInterfaceName()); 922 if (SupportedTypes::getTypeFlag(f.type) == SupportedTypes::STENSOR) { 923 os << "calculix::ImportThermodynamicForces::exe(this->" << f.name << ","; 924 if (!o.isNull()) { 925 os << iprefix << "stress_+" << o << ");\n"; 926 } else { 927 os << iprefix << "stress_);\n"; 928 } 929 } else { 930 tfel::raise( 931 "CalculiXInterface::writeBehaviourDataMainVariablesSetters : " 932 "unsupported forces type"); 933 } 934 } // end of CalculiXInterface::writeBehaviourDataThermodynamicForceSetter 935 exportThermodynamicForce(std::ostream & out,const std::string & a,const ThermodynamicForce & f,const SupportedTypes::TypeSize o) const936 void CalculiXInterface::exportThermodynamicForce( 937 std::ostream& out, 938 const std::string& a, 939 const ThermodynamicForce& f, 940 const SupportedTypes::TypeSize o) const { 941 const auto iprefix = makeUpperCase(this->getInterfaceName()); 942 const auto flag = SupportedTypes::getTypeFlag(f.type); 943 if (flag == SupportedTypes::STENSOR) { 944 if (!o.isNull()) { 945 out << "calculix::ExportThermodynamicForces::exe(" << a << "+" << o 946 << ",this->sig);\n"; 947 } else { 948 out << "calculix::ExportThermodynamicForces::exe(" << a 949 << ",this->sig);\n"; 950 } 951 } else { 952 tfel::raise( 953 "CalculiXInterface::exportThermodynamicForce: " 954 "unsupported forces type"); 955 } 956 } // end of CalculiXInterface::exportThermodynamicForce 957 getTargetsDescription(TargetsDescription & d,const BehaviourDescription & bd)958 void CalculiXInterface::getTargetsDescription( 959 TargetsDescription& d, const BehaviourDescription& bd) { 960 const auto lib = this->getLibraryName(bd); 961 const auto name = bd.getLibrary() + bd.getClassName(); 962 const auto tfel_config = tfel::getTFELConfigExecutableName(); 963 auto& l = d.getLibrary(lib); 964 insert_if(l.cppflags, 965 "$(shell " + tfel_config + " --cppflags --compiler-flags)"); 966 insert_if(l.include_directories, 967 "$(shell " + tfel_config + " --include-path)"); 968 insert_if(l.sources, "calculix" + name + ".cxx"); 969 d.headers.push_back("MFront/CalculiX/calculix" + name + ".hxx"); 970 insert_if(l.link_directories, 971 "$(shell " + tfel_config + " --library-path)"); 972 insert_if(l.link_libraries, 973 tfel::getLibraryInstallName("CalculiXInterface")); 974 if (this->shallGenerateMTestFileOnFailure(bd)) { 975 insert_if(l.link_libraries, 976 tfel::getLibraryInstallName("MTestFileGenerator")); 977 } 978 #if __cplusplus >= 201703L 979 insert_if(l.link_libraries, "$(shell " + tfel_config + 980 " --library-dependency " 981 "--material --mfront-profiling)"); 982 #else /* __cplusplus < 201703L */ 983 insert_if(l.link_libraries, 984 "$(shell " + tfel_config + 985 " --library-dependency " 986 "--material --mfront-profiling --physical-constants)"); 987 #endif /* __cplusplus < 201703L */ 988 insert_if(l.epts, this->getFunctionNameBasis(name)); 989 } // end of CalculiXInterface::getTargetsDescription 990 getLibraryName(const BehaviourDescription & mb) const991 std::string CalculiXInterface::getLibraryName( 992 const BehaviourDescription& mb) const { 993 auto lib = std::string{}; 994 if (mb.getLibrary().empty()) { 995 if (!mb.getMaterialName().empty()) { 996 lib = this->getInterfaceName() + mb.getMaterialName(); 997 } else { 998 lib = this->getInterfaceName() + "Behaviour"; 999 } 1000 } else { 1001 lib = this->getInterfaceName() + mb.getLibrary(); 1002 } 1003 return makeUpperCase(lib); 1004 } // end of CalculiXInterface::getLibraryName 1005 getFunctionNameBasis(const std::string & name) const1006 std::string CalculiXInterface::getFunctionNameBasis( 1007 const std::string& name) const { 1008 return makeUpperCase(name); 1009 } // end of CalculiXInterface::getFunctionName 1010 1011 std::set<CalculiXInterface::Hypothesis> getModellingHypothesesToBeTreated(const BehaviourDescription & mb) const1012 CalculiXInterface::getModellingHypothesesToBeTreated( 1013 const BehaviourDescription& mb) const { 1014 const auto& bh = mb.getModellingHypotheses(); 1015 tfel::raise_if(bh.find(ModellingHypothesis::TRIDIMENSIONAL) == bh.end(), 1016 "CalculiXInterface::getModellingHypothesesToBeTreated : " 1017 "the 'Tridimensional' hypothesis is not supported, " 1018 "which is required for the CalculiX interface"); 1019 return {ModellingHypothesis::TRIDIMENSIONAL}; 1020 } // end of CalculiXInterface::getModellingHypothesesToBeTreated 1021 writeCalculiXBehaviourTraits(std::ostream & out,const BehaviourDescription & mb) const1022 void CalculiXInterface::writeCalculiXBehaviourTraits( 1023 std::ostream& out, const BehaviourDescription& mb) const { 1024 constexpr const auto h = ModellingHypothesis::TRIDIMENSIONAL; 1025 const auto mvs = mb.getMainVariablesSize(); 1026 const auto mprops = this->buildMaterialPropertiesList(mb, h); 1027 out << "template<typename Type"; 1028 if (mb.useQt()) { 1029 out << ",bool use_qt"; 1030 } 1031 out << ">\n" 1032 << "struct CalculiXTraits<tfel::material::" << mb.getClassName() 1033 << "<tfel::material::ModellingHypothesis::TRIDIMENSIONAL,"; 1034 out << "Type,"; 1035 if (mb.useQt()) { 1036 out << "use_qt"; 1037 } else { 1038 out << "false"; 1039 } 1040 out << ">>\n{\n" 1041 << "//! behaviour type\n"; 1042 if (mb.getBehaviourType() == 1043 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) { 1044 out << "static " << constexpr_c << " CalculiXBehaviourType btype = " 1045 "calculix::" 1046 "STANDARDSTRAINBASEDBEHAVIOUR;\n"; 1047 } else if (mb.getBehaviourType() == 1048 BehaviourDescription::STANDARDFINITESTRAINBEHAVIOUR) { 1049 out << "static " << constexpr_c << " CalculiXBehaviourType btype = " 1050 "calculix::" 1051 "STANDARDFINITESTRAINBEHAVIOUR;\n"; 1052 } else { 1053 tfel::raise( 1054 "CalculiXInterface::writeCalculiXBehaviourTraits : " 1055 "unsupported behaviour type"); 1056 } 1057 out << "//! space dimension\n" 1058 << "static " << constexpr_c << " unsigned short N " 1059 << "= " 1060 "tfel::material::ModellingHypothesisToSpaceDimension<tfel::material:" 1061 ":ModellingHypothesis::TRIDIMENSIONAL>::value;\n" 1062 << "// tiny vector size\n" 1063 << "static " << constexpr_c << " unsigned short TVectorSize = N;\n" 1064 << "// symmetric tensor size\n" 1065 << "static " << constexpr_c 1066 << " unsigned short StensorSize = " 1067 "tfel::math::StensorDimeToSize<N>::value;\n" 1068 << "// tensor size\n" 1069 << "static " << constexpr_c 1070 << " unsigned short TensorSize = " 1071 "tfel::math::TensorDimeToSize<N>::value;\n" 1072 << "// size of the driving variable array\n" 1073 << "static " << constexpr_c 1074 << " unsigned short GradientSize = " << mvs.first << ";\n" 1075 << "// size of the thermodynamic force variable array (STRESS)\n" 1076 << "static " << constexpr_c 1077 << " unsigned short ThermodynamicForceVariableSize = " << mvs.second 1078 << ";\n"; 1079 if (mb.getAttribute(BehaviourDescription::requiresUnAlteredStiffnessTensor, 1080 false)) { 1081 out << "static " << constexpr_c 1082 << " bool requiresUnAlteredStiffnessTensor = true;\n"; 1083 } else { 1084 out << "static " << constexpr_c 1085 << " bool requiresUnAlteredStiffnessTensor = false;\n"; 1086 } 1087 if (mb.getAttribute(BehaviourDescription::requiresStiffnessTensor, false)) { 1088 out << "static " << constexpr_c 1089 << " bool requiresStiffnessTensor = true;\n"; 1090 } else { 1091 out << "static " << constexpr_c 1092 << " bool requiresStiffnessTensor = false;\n"; 1093 } 1094 if (mb.getAttribute( 1095 BehaviourDescription::requiresThermalExpansionCoefficientTensor, 1096 false)) { 1097 out << "static " << constexpr_c 1098 << " bool requiresThermalExpansionCoefficientTensor = true;\n"; 1099 } else { 1100 out << "static " << constexpr_c 1101 << " bool requiresThermalExpansionCoefficientTensor = false;\n"; 1102 } 1103 if (mb.getSymmetryType() == mfront::ISOTROPIC) { 1104 out << "static " << constexpr_c 1105 << " CalculiXSymmetryType type = calculix::ISOTROPIC;\n"; 1106 } else if (mb.getSymmetryType() == mfront::ORTHOTROPIC) { 1107 out << "static " << constexpr_c 1108 << " CalculiXSymmetryType type = calculix::ORTHOTROPIC;\n"; 1109 } else { 1110 tfel::raise( 1111 "CalculiXInterface::writeCalculiXBehaviourTraits: " 1112 "unsupported behaviour type.\n" 1113 "The calculix interface only support isotropic or orthotropic " 1114 "behaviour at this time."); 1115 } 1116 // computing material properties size 1117 auto msize = SupportedTypes::TypeSize{}; 1118 if (!mprops.first.empty()) { 1119 const auto& m = mprops.first.back(); 1120 msize = m.offset; 1121 msize += SupportedTypes::getTypeSize(m.type, m.arraySize); 1122 msize -= mprops.second; 1123 } 1124 out << "static " << constexpr_c 1125 << " unsigned short material_properties_nb = " << msize << ";\n"; 1126 if (mb.getElasticSymmetryType() == mfront::ISOTROPIC) { 1127 out << "static " << constexpr_c 1128 << " CalculiXSymmetryType etype = calculix::ISOTROPIC;\n"; 1129 if (mb.getAttribute(BehaviourDescription::requiresStiffnessTensor, 1130 false)) { 1131 out << "static " << constexpr_c 1132 << " unsigned short elasticPropertiesOffset = 2u;\n"; 1133 } else { 1134 out << "static " << constexpr_c 1135 << " unsigned short elasticPropertiesOffset = 0u;\n"; 1136 } 1137 if (mb.getAttribute( 1138 BehaviourDescription::requiresThermalExpansionCoefficientTensor, 1139 false)) { 1140 out << "static " << constexpr_c 1141 << " unsigned short thermalExpansionPropertiesOffset = 1u;\n"; 1142 } else { 1143 out << "static " << constexpr_c 1144 << " unsigned short thermalExpansionPropertiesOffset = 0u;\n"; 1145 } 1146 } else if (mb.getElasticSymmetryType() == mfront::ORTHOTROPIC) { 1147 out << "static " << constexpr_c 1148 << " CalculiXSymmetryType etype = calculix::ORTHOTROPIC;\n"; 1149 if (mb.getAttribute(BehaviourDescription::requiresStiffnessTensor, 1150 false)) { 1151 out << "static " << constexpr_c 1152 << " unsigned short elasticPropertiesOffset " 1153 << "= 9u;\n"; 1154 } else { 1155 out << "static " << constexpr_c 1156 << " unsigned short elasticPropertiesOffset = 0u;\n"; 1157 } 1158 if (mb.getAttribute( 1159 BehaviourDescription::requiresThermalExpansionCoefficientTensor, 1160 false)) { 1161 out << "static " << constexpr_c 1162 << " unsigned short thermalExpansionPropertiesOffset = 3u;\n"; 1163 } else { 1164 out << "static " << constexpr_c 1165 << " unsigned short thermalExpansionPropertiesOffset = 0u;\n"; 1166 } 1167 } else { 1168 tfel::raise( 1169 "CalculiXInterface::writeCalculiXBehaviourTraits: " 1170 "unsupported behaviour type.\n" 1171 "The calculix interface only support isotropic or " 1172 "orthotropic behaviour at this time."); 1173 } 1174 out << "}; // end of class CalculiXTraits\n\n"; 1175 } 1176 1177 std::map<UMATInterfaceBase::Hypothesis, std::string> gatherModellingHypothesesAndTests(const BehaviourDescription & mb) const1178 CalculiXInterface::gatherModellingHypothesesAndTests( 1179 const BehaviourDescription& mb) const { 1180 auto res = std::map<Hypothesis, std::string>{}; 1181 if ((mb.getSymmetryType() == mfront::ORTHOTROPIC) && 1182 ((mb.getAttribute(BehaviourDescription::requiresStiffnessTensor, 1183 false)) || 1184 (mb.getAttribute( 1185 BehaviourDescription::requiresThermalExpansionCoefficientTensor, 1186 false)))) { 1187 for (const auto& h : this->getModellingHypothesesToBeTreated(mb)) { 1188 res.insert({h, this->getModellingHypothesisTest(h)}); 1189 } 1190 return res; 1191 } 1192 return UMATInterfaceBase::gatherModellingHypothesesAndTests(mb); 1193 } // end of CalculiXInterface::gatherModellingHypothesesAndTests 1194 getModellingHypothesisTest(const Hypothesis h) const1195 std::string CalculiXInterface::getModellingHypothesisTest( 1196 const Hypothesis h) const { 1197 if (h == ModellingHypothesis::TRIDIMENSIONAL) { 1198 return "true"; 1199 } 1200 tfel::raise( 1201 "CalculiXInterface::getModellingHypothesisTest : " 1202 "unsupported modelling hypothesis"); 1203 } // end of CalculiXInterface::gatherModellingHypothesesAndTests 1204 areExternalStateVariablesSupported() const1205 bool CalculiXInterface::areExternalStateVariablesSupported() const { 1206 return false; 1207 } // end of CalculiXInterface::areExternalStateVariablesSupported() 1208 isTemperatureIncrementSupported() const1209 bool CalculiXInterface::isTemperatureIncrementSupported() const { 1210 return false; 1211 } // end of CalculiXInterface::isTemperatureIncrementSupported() 1212 writeMTestFileGeneratorSetModellingHypothesis(std::ostream & out) const1213 void CalculiXInterface::writeMTestFileGeneratorSetModellingHypothesis( 1214 std::ostream& out) const { 1215 out << "mg.setModellingHypothesis(ModellingHypothesis::TRIDIMENSIONAL);\n"; 1216 } // end of CalculiXInterface::writeMTestFileGeneratorSetModellingHypothesis 1217 writeInputFileExample(const BehaviourDescription & mb,const FileDescription & fd,const bool b) const1218 void CalculiXInterface::writeInputFileExample(const BehaviourDescription& mb, 1219 const FileDescription& fd, 1220 const bool b) const { 1221 auto throw_if = [](const bool c, const std::string& m) { 1222 tfel::raise_if(c, "CalculiXInterface::writeInputFileExample: " + m); 1223 }; 1224 const auto name = mb.getLibrary() + mb.getClassName(); 1225 const auto mn = this->getLibraryName(mb) + "_" + mb.getClassName(); 1226 const auto fn = (b ? "calculix/" : "calculix-explicit/") + name + ".inp"; 1227 std::ofstream out{fn}; 1228 throw_if(!out, "could not open file '" + fn + "'"); 1229 // header 1230 out << "** \n" 1231 << "** File generated by MFront from the " << fd.fileName << " source\n" 1232 << "** Example of how to use the " << mb.getClassName() 1233 << " behaviour law\n" 1234 << "** Author " << fd.authorName << '\n' 1235 << "** Date " << fd.date << '\n' 1236 << "**\n\n"; 1237 // loop over hypothesis 1238 const auto h = ModellingHypothesis::TRIDIMENSIONAL; 1239 const auto& d = mb.getBehaviourData(h); 1240 const auto mps = this->buildMaterialPropertiesList(mb, h); 1241 auto msize = SupportedTypes::TypeSize{}; 1242 if (!mps.first.empty()) { 1243 const auto& m = mps.first.back(); 1244 msize = m.offset; 1245 msize += SupportedTypes::getTypeSize(m.type, m.arraySize); 1246 } 1247 const auto& persistentVarsHolder = d.getPersistentVariables(); 1248 auto vs = SupportedTypes::TypeSize{}; 1249 for (const auto& v : persistentVarsHolder) { 1250 vs += SupportedTypes::getTypeSize(v.type, v.arraySize); 1251 } 1252 // this is included for gcc 4.7.2 1253 const auto vsize = [mb, &vs, h, this]() -> unsigned int { 1254 if ((mb.getBehaviourType() == 1255 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR) && 1256 (CalculiXInterface::hasFiniteStrainStrategy(mb)) && 1257 (getFiniteStrainStrategy(mb) == 1258 "MieheApelLambrechtLogarithmicStrainII")) { 1259 return vs.getValueForModellingHypothesis(h) + 12u; 1260 } 1261 return vs.getValueForModellingHypothesis(h); 1262 }(); 1263 out << "*Material, name=@" << this->getFunctionNameBasis(mn) << '\n'; 1264 if (!b) { 1265 out << "*DENSITY\n<density>\n"; 1266 } 1267 if (vsize != 0) { 1268 out << "*Depvar\n" << vsize << "\n"; 1269 } 1270 if (!mps.first.empty()) { 1271 out << "** The material properties are given as if we used parameters to " 1272 "explicitly\n" 1273 << "** display their names. Users shall replace those declaration " 1274 "by\n" 1275 << "** theirs values\n"; 1276 } 1277 out << "*User Material, constants=" 1278 << msize.getValueForModellingHypothesis(h); 1279 out << '\n'; 1280 if (!mps.first.empty()) { 1281 int i = 1; 1282 auto write = [&i, &out](const std::string& n) { 1283 if (i % 9 == 0) { 1284 out << "\n"; 1285 i = 1; 1286 } 1287 out << '<' << n << '>'; 1288 ++i; 1289 }; 1290 for (auto pm = mps.first.begin(); pm != mps.first.end();) { 1291 if (pm->arraySize == 1u) { 1292 write(pm->name); 1293 } else { 1294 for (unsigned short a = 0; a != pm->arraySize;) { 1295 write(pm->name + "_" + std::to_string(a)); 1296 if (++a != pm->arraySize) { 1297 out << ", "; 1298 } 1299 } 1300 } 1301 if (++pm != mps.first.end()) { 1302 out << ", "; 1303 } 1304 } 1305 } 1306 out << "\n\n"; 1307 } // end of CalculiXInterface::writeInputFileExample 1308 1309 CalculiXInterface::~CalculiXInterface() = default; 1310 1311 } // end of namespace mfront 1312