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