1 /*! 2 * \file mfront/src/ComsolInterface.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 30/07/2019 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 "TFEL/Raise.hxx" 16 #include "TFEL/Config/GetInstallPath.hxx" 17 #include "TFEL/System/System.hxx" 18 #include "MFront/MFrontUtilities.hxx" 19 #include "MFront/FileDescription.hxx" 20 #include "MFront/TargetsDescription.hxx" 21 #include "MFront/ComsolInterface.hxx" 22 23 namespace mfront { 24 getName()25 std::string ComsolInterface::getName() { 26 return "comsol"; 27 } // end of getName 28 getInterfaceName() const29 std::string ComsolInterface::getInterfaceName() const { 30 return "Comsol"; 31 } // end of ComsolInterface::getInterfaceName 32 33 std::pair<bool, ComsolInterface::tokens_iterator> treatKeyword(BehaviourDescription &,const std::string &,const std::vector<std::string> &,tokens_iterator current,const tokens_iterator)34 ComsolInterface::treatKeyword(BehaviourDescription&, 35 const std::string&, 36 const std::vector<std::string>&, 37 tokens_iterator current, 38 const tokens_iterator) { 39 return {false, current}; 40 } // end of ComsolInterface::treatKeyword 41 42 std::set<ComsolInterface::Hypothesis> getModellingHypothesesToBeTreated(const BehaviourDescription &) const43 ComsolInterface::getModellingHypothesesToBeTreated( 44 const BehaviourDescription&) const { 45 return {ModellingHypothesis::TRIDIMENSIONAL}; 46 } // end of ComsolInterface::getModellingHypothesesToBeTreated 47 writeInterfaceSpecificIncludes(std::ostream &,const BehaviourDescription &) const48 void ComsolInterface::writeInterfaceSpecificIncludes( 49 std::ostream&, const BehaviourDescription&) const { 50 } // end of ComsolInterface::writeInterfaceSpecificIncludes 51 getTargetsDescription(TargetsDescription & d,const BehaviourDescription & bd)52 void ComsolInterface::getTargetsDescription( 53 TargetsDescription& d, const BehaviourDescription& bd) { 54 const auto lib = this->getLibraryName(bd); 55 const auto name = bd.getLibrary() + bd.getClassName(); 56 const auto tfel_config = tfel::getTFELConfigExecutableName(); 57 auto& l = d.getLibrary(lib); 58 insert_if(l.cppflags, 59 "$(shell " + tfel_config + " --cppflags --compiler-flags)"); 60 insert_if(l.include_directories, 61 "$(shell " + tfel_config + " --include-path)"); 62 insert_if(l.sources, "comsol" + name + ".cxx"); 63 d.headers.push_back("MFront/Comsol/comsol" + name + ".hxx"); 64 insert_if(l.link_directories, 65 "$(shell " + tfel_config + " --library-path)"); 66 if (this->shallGenerateMTestFileOnFailure(bd)) { 67 insert_if(l.link_libraries, 68 tfel::getLibraryInstallName("MTestFileGenerator")); 69 } 70 #if __cplusplus >= 201703L 71 insert_if(l.link_libraries, "$(shell " + tfel_config + 72 " --library-dependency " 73 "--material --mfront-profiling)"); 74 #else /* __cplusplus < 201703L */ 75 insert_if(l.link_libraries, 76 "$(shell " + tfel_config + 77 " --library-dependency " 78 "--material --mfront-profiling --physical-constants)"); 79 #endif /* __cplusplus < 201703L */ 80 for (const auto h : this->getModellingHypothesesToBeTreated(bd)) { 81 if (h != ModellingHypothesis::TRIDIMENSIONAL) { 82 tfel::raise( 83 "ComsolInterface::getTargetsDescription: " 84 "unsupported modelling hypothesis"); 85 } 86 insert_if(l.epts, this->getFunctionNameBasis(name)); 87 } 88 } // end of ComsolInterface::getTargetsDescription 89 endTreatment(const BehaviourDescription & bd,const FileDescription & fd) const90 void ComsolInterface::endTreatment(const BehaviourDescription& bd, 91 const FileDescription& fd) const { 92 using namespace tfel::system; 93 auto throw_if = [](const bool b, const std::string& m) { 94 tfel::raise_if(b, "ComsolInterface::endTreatment: " + m); 95 }; 96 throw_if(bd.getBehaviourType() != 97 BehaviourDescription::STANDARDSTRAINBASEDBEHAVIOUR, 98 "the comsol interface only supports strain based behaviours"); 99 if (bd.isStrainMeasureDefined()) { 100 const auto ms = bd.getStrainMeasure(); 101 throw_if(!((ms == BehaviourDescription::LINEARISED) || 102 (ms == BehaviourDescription::GREENLAGRANGE)), 103 "the comsol interface only supports the linearised " 104 "and Green-Lagrange strain measures"); 105 } 106 throw_if(bd.getSymmetryType() != mfront::ISOTROPIC, 107 "the comsol interface only supports isotropic behaviours"); 108 // the only supported modelling hypothesis 109 constexpr const auto h = ModellingHypothesis::TRIDIMENSIONAL; 110 const auto& d = bd.getBehaviourData(h); 111 throw_if(d.getExternalStateVariables().size() != 1u, 112 "external state variables are not supported " 113 "by comsol interface"); 114 // get the modelling hypotheses to be treated 115 const auto name = bd.getLibrary() + bd.getClassName(); 116 // output directories 117 systemCall::mkdir("include/MFront"); 118 systemCall::mkdir("include/MFront/Comsol"); 119 systemCall::mkdir("comsol"); 120 121 std::ofstream out; 122 123 // header 124 auto fname = "comsol" + name + ".hxx"; 125 out.open("include/MFront/Comsol/" + fname); 126 throw_if(!out, "could not open file '" + fname + "'"); 127 128 out << "/*!\n" 129 << "* \\file " << fname << '\n' 130 << "* \\brief This file declares the comsol interface for the " 131 << bd.getClassName() << " behaviour law\n" 132 << "* \\author " << fd.authorName << '\n' 133 << "* \\date " << fd.date << '\n' 134 << "*/\n\n"; 135 136 const auto header = this->getHeaderGuard(bd); 137 out << "#ifndef " << header << "\n" 138 << "#define " << header << "\n\n" 139 << "#ifdef __cplusplus\n" 140 << "#include\"TFEL/Config/TFELConfig.hxx\"\n" 141 << "#include\"TFEL/Material/" << bd.getClassName() << ".hxx\"\n" 142 << "#endif /* __cplusplus */\n\n"; 143 144 this->writeVisibilityDefines(out); 145 146 out << "#ifdef __cplusplus\n" 147 << "extern \"C\"{\n" 148 << "#endif /* __cplusplus */\n\n"; 149 150 this->writeSetOutOfBoundsPolicyFunctionDeclaration(out, name); 151 this->writeSetParametersFunctionsDeclarations(out, bd, name); 152 153 if (!d.getPersistentVariables().empty()) { 154 out << "MFRONT_SHAREDOBJ int\n" 155 << this->getFunctionNameBasis(name) 156 << "(const double * const, // strain\n" 157 << "double * const, // stress\n" 158 << "double * const, // tangent operator\n" 159 << "const int * const, // number of material properties\n" 160 << "const double *const, // material properties\n" 161 << "const int *const, // number of strain components\n" 162 << "double *const, // strain at the beginning of the time step\n" 163 << "const int *const, // number of state variables\n" 164 << "double *const);// state variables\n\n"; 165 } else { 166 out << "MFRONT_SHAREDOBJ int\n" 167 << this->getFunctionNameBasis(name) 168 << "(const double * const, // strain\n" 169 << "double * const, // stress\n" 170 << "double * const, // tangent operator\n" 171 << "const int * const, // number of material properties\n" 172 << "const double *const, // material properties\n" 173 << "const int *const, // number of strain components\n" 174 << "double *const); // strain at the beginning of the time step\n\n"; 175 } 176 177 out << "#ifdef __cplusplus\n" 178 << "}\n" 179 << "#endif /* __cplusplus */\n\n" 180 << "#endif /* " << header << " */\n"; 181 182 out.close(); 183 184 fname = "comsol" + name + ".cxx"; 185 out.open("src/" + fname); 186 throw_if(!out, "could not open file '" + fname + "'"); 187 188 out << "/*!\n" 189 << "* \\file " << fname << '\n' 190 << "* \\brief This file implements the comsol interface for the " 191 << bd.getClassName() << " behaviour law\n" 192 << "* \\author " << fd.authorName << '\n' 193 << "* \\date " << fd.date << '\n' 194 << "*/\n\n"; 195 196 // this->getExtraSrcIncludes(out, mb); 197 // 198 out << "#include\"TFEL/Math/General/MathConstants.hxx\"\n" 199 << "#include\"TFEL/Math/stensor.hxx\"\n" 200 << "#include\"TFEL/Material/OutOfBoundsPolicy.hxx\"\n" 201 << "#include\"TFEL/Material/" << bd.getClassName() << ".hxx\"\n"; 202 if (bd.getAttribute(BehaviourData::profiling, false)) { 203 out << "#include\"MFront/BehaviourProfiler.hxx\"\n\n"; 204 } 205 // out << "#include\"MFront/Comsol/" 206 // "ComsolStressFreeExpansionHandler.hxx\"\n\n" 207 // << "#include\"MFront/Comsol/ComsolInterface.hxx\"\n\n" 208 out << "#include\"MFront/Comsol/comsol" << name << ".hxx\"\n\n"; 209 // 210 this->writeGetOutOfBoundsPolicyFunctionImplementation(out, name); 211 212 out << "extern \"C\"{\n\n"; 213 // 214 // ComsolSymbolsGenerator sg; 215 // sg.generateGeneralSymbols(out, *this, mb, fd, {h}, name); 216 // sg.generateSymbols(out, *this, mb, fd, name, h); 217 // 218 this->writeSetParametersFunctionsImplementations(out, bd, name); 219 this->writeSetOutOfBoundsPolicyFunctionImplementation(out, name); 220 221 if (!d.getPersistentVariables().empty()) { 222 out << "MFRONT_SHAREDOBJ int\n" 223 << this->getFunctionNameBasis(name) // 224 << "(const double * const e,\n" 225 << "double * const s,\n" 226 << "double * const D,\n" 227 << "const int * const nprops,\n" 228 << "const double *const mprops,\n" 229 << "const int *const nstran,\n" 230 << "double *const stran,\n" 231 << "const int *const nstatv,\n" 232 << "double *const statev){\n"; 233 } else { 234 out << "MFRONT_SHAREDOBJ int\n" 235 << this->getFunctionNameBasis(name) // 236 << "(const double * const e,\n" 237 << "double * const s,\n" 238 << "double * const D,\n" 239 << "const int * const nprops,\n" 240 << "const double *const mprops,\n" 241 << "const int *const nstran,\n" 242 << "double *const stran){\n"; 243 } 244 out << "using tfel::material::ModellingHypothesis;\n" 245 << "using Behaviour = tfel::material::" << bd.getClassName() 246 << "<ModellingHypothesis::TRIDIMENSIONAL, double, " 247 "false>;\n"; 248 if (bd.requiresStressFreeExpansionTreatment(h)) { 249 out << "using StressFreeExpansionType = " 250 "Behaviour::StressFreeExpansionType;\n"; 251 } 252 out << "constexpr const auto cste = tfel::math::Cste<double>::sqrt2;\n"; 253 // inputs checks 254 const auto nprops = 255 d.getMaterialProperties().getTypeSize().getValueForDimension(3) + 1; 256 out << "if (*nprops != " << nprops << " ){\n" 257 << "return 1;\n" 258 << "}\n" 259 << "if (*nstran!= 6){\n" 260 << "return 2;\n" 261 << "}\n"; 262 if (!d.getPersistentVariables().empty()) { 263 const auto nstatv = 264 d.getPersistentVariables().getTypeSize().getValueForDimension(3); 265 out << "if (*nstatv!= " << nstatv << " ){\n" 266 << "return 2;\n" 267 << "}\n"; 268 } 269 out << "auto rdt = double{1};\n" 270 << "try{\n" 271 << "const auto p = " << this->getFunctionNameBasis(name) 272 << "_getOutOfBoundsPolicy();\n" 273 << "const auto dt = double{0};\n" 274 << "const auto T = mprops[" << nprops << "];\n" 275 << "const tfel::math::stensor<3u,double> e0(stran);\n" 276 << "tfel::math::stensor<3u,double> e1(e);\n" 277 << "tfel::math::stensor<3u,double> eto;\n" 278 << "tfel::math::stensor<3u,double> deto;\n"; 279 if (!d.getPersistentVariables().empty()) { 280 out << "Behaviour b(&dt, &T, mprops, statev);\n"; 281 } else { 282 out << "Behaviour b(&dt, &T, mprops, nullptr);\n"; 283 } 284 out << "// convert e1 in TFEL conventions\n" 285 << "e1[3] *= cste;\n" 286 << "e1[4] *= cste;\n" 287 << "e1[5] *= cste;\n" 288 << "std::swap(e1[3],e1[5]);\n"; 289 if (bd.requiresStressFreeExpansionTreatment(h)) { 290 out << "// stress free expansion\n" 291 << "auto s = StressFreeExpansionType{};\n" 292 << "b.computeStressFreeExpansion(s);\n" 293 << "const auto& s0 = s.first;\n" 294 << "const auto& s1 = s.second;\n" 295 << "eto = e0-s0;\n" 296 << "deto = e1-e0-(s1-s0);\n"; 297 } else { 298 out << "eto = e0;\n" 299 << "deto = e1-e0;\n"; 300 } 301 out << "// to Voigt conventions\n" 302 << "eto[3] *= cste;\n" 303 << "eto[4] *= cste;\n" 304 << "eto[5] *= cste;\n" 305 << "deto[3] *= cste;\n" 306 << "deto[4] *= cste;\n" 307 << "deto[5] *= cste;\n" 308 << "b.setCOMSOLBehaviourDataGradients(eto.begin());\n" 309 << "b.setCOMSOLIntegrationDataGradients(deto.begin());\n" 310 << "b.setOutOfBoundsPolicy(p);\n" 311 << "b.initialize();\n" 312 << "b.checkBounds();\n" 313 << "const auto smflag = Behaviour::STANDARDTANGENTOPERATOR;\n" 314 << "const auto smt = Behaviour::CONSISTENTTANGENTOPERATOR;\n" 315 << "auto tsf = b.computeAPrioriTimeStepScalingFactor(rdt);\n" 316 << "rdt = tsf.second;\n" 317 << "if (!tsf.first) {\n" 318 << " return -1;\n" 319 << "}\n" 320 << "const auto r = b.integrate(smflag, smt);\n" 321 << "if (r == Behaviour::FAILURE) {\n" 322 << " return -1;\n" 323 << "}\n" 324 << "const auto atsf = " 325 "b.computeAPosterioriTimeStepScalingFactor(rdt);\n" 326 << "if (!atsf.first) {\n" 327 << " return -1;\n" 328 << "}\n" 329 << "rdt = std::min(atsf.second, rdt);\n"; 330 if (!d.getPersistentVariables().empty()) { 331 out << "b.COMSOLexportStateData(s,statev);\n"; 332 } else { 333 out << "b.COMSOLexportStateData(s,nullptr);\n"; 334 } 335 out << "std::swap(s[3],s[5]);\n" 336 << "const auto& K = b.getTangentOperator();\n" 337 << "std::copy(K.begin(),K.end(), D);\n" 338 << "tfel::fsalgo::copy<6>::exe(e1.begin(),stran);\n" 339 << "} catch(...){\n" 340 << "return -1;\n" 341 << "}\n" 342 << "return rdt < 0.99 ? -1 : 0;\n" 343 << "} // end of " << this->getFunctionNameBasis(name) << "\n\n"; 344 345 out << "} // end of extern \"C\"\n"; 346 out.close(); 347 } // end of ComsolInterface::endTreatment 348 areExternalStateVariablesSupported() const349 bool ComsolInterface::areExternalStateVariablesSupported() const { 350 return false; 351 } // end of ComsolInterface::areExternalStateVariablesSupported() 352 isTemperatureIncrementSupported() const353 bool ComsolInterface::isTemperatureIncrementSupported() const { 354 return false; 355 } // end of ComsolInterface::isTemperatureIncrementSupported() 356 getLibraryName(const BehaviourDescription & bd) const357 std::string ComsolInterface::getLibraryName( 358 const BehaviourDescription& bd) const { 359 if (bd.getLibrary().empty()) { 360 if (!bd.getMaterialName().empty()) { 361 return this->getInterfaceName() + bd.getMaterialName(); 362 } else { 363 return this->getInterfaceName() + "Behaviour"; 364 } 365 } 366 return this->getInterfaceName() + bd.getLibrary(); 367 } // end of ComsolInterface::getLibraryName 368 getFunctionNameBasis(const std::string & n) const369 std::string ComsolInterface::getFunctionNameBasis( 370 const std::string& n) const { 371 return n; 372 } // end of ComsolInterface::getFunctionName 373 writeMTestFileGeneratorSetModellingHypothesis(std::ostream & out) const374 void ComsolInterface::writeMTestFileGeneratorSetModellingHypothesis( 375 std::ostream& out) const { 376 out << "mg.setModellingHypothesis(ModellingHypothesis::TRIDIMENSIONAL);\n"; 377 } // end of ComsolInterface::writeMTestFileGeneratorSetModellingHypothesis 378 379 ComsolInterface::~ComsolInterface() = default; 380 381 } // end of namespace mfront 382