1 /*! 2 * \file OctaveMaterialPropertyInterface.cxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 06 mai 2008 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<sstream> 15 #include<stdexcept> 16 #include<vector> 17 #include<string> 18 #include<iterator> 19 #include<algorithm> 20 #include<cstdlib> 21 22 #include"TFEL/Raise.hxx" 23 #include"TFEL/System/System.hxx" 24 #include"MFront/MFrontHeader.hxx" 25 #include"MFront/DSLUtilities.hxx" 26 #include"MFront/FileDescription.hxx" 27 #include"MFront/TargetsDescription.hxx" 28 #include"MFront/MaterialPropertyDescription.hxx" 29 #include"MFront/MaterialPropertyParametersHandler.hxx" 30 #include"MFront/CMaterialPropertyInterfaceBase.hxx" 31 #include"MFront/OctaveMaterialPropertyInterface.hxx" 32 33 namespace mfront 34 { 35 36 static std::vector<std::string> tokenize(const std::string & s,const char c)37 tokenize(const std::string& s, 38 const char c) 39 { 40 using namespace std; 41 vector<string> res; 42 string::size_type b = 0u; 43 string::size_type e = s.find_first_of(c, b); 44 while (string::npos != e || string::npos != b){ 45 // Found a token, add it to the vector. 46 res.push_back(s.substr(b, e - b)); 47 b = s.find_first_not_of(c, e); 48 e = s.find_first_of(c, b); 49 } 50 return res; 51 } // end of tokenize 52 53 std::string toString(const unsigned short src)54 OctaveMaterialPropertyInterface::toString(const unsigned short src) 55 { 56 std::ostringstream os; 57 os << src; 58 return os.str(); 59 } 60 61 std::string getName()62 OctaveMaterialPropertyInterface::getName() 63 { 64 return "octave"; 65 } 66 67 OctaveMaterialPropertyInterface::OctaveMaterialPropertyInterface() = default; 68 69 std::pair<bool,OctaveMaterialPropertyInterface::tokens_iterator> treatKeyword(const std::string & k,const std::vector<std::string> & i,tokens_iterator current,const tokens_iterator)70 OctaveMaterialPropertyInterface::treatKeyword(const std::string& k, 71 const std::vector<std::string>& i, 72 tokens_iterator current, 73 const tokens_iterator) 74 { 75 tfel::raise_if(std::find(i.begin(),i.end(),"octave")!=i.end(), 76 "OctaveMaterialPropertyInterface::treatKeyword: " 77 "unsupported key '"+k+"'"); 78 return {false,current}; 79 } // end of treatKeyword 80 81 void getTargetsDescription(TargetsDescription & td,const MaterialPropertyDescription & mpd) const82 OctaveMaterialPropertyInterface::getTargetsDescription(TargetsDescription& td, 83 const MaterialPropertyDescription& mpd) const 84 { 85 const char * mkoctfile = ::getenv("MKOCTFILE"); 86 const auto name = (mpd.material.empty()) ? mpd.className : mpd.material+"_"+mpd.className; 87 const auto target="../octave/"+name+".oct"; 88 std::string cmd = "@cd ../octave/ && CXXFLAGS=\"$(CXXFLAGS) -std=c++11\" "; 89 if(mkoctfile==nullptr){ 90 cmd += "mkoctfile"; 91 } else { 92 cmd += mkoctfile; 93 } 94 cmd += " $(INCLUDES) -L../src/"; 95 cmd += " "+name+".cpp"; 96 auto& res = td.specific_targets; 97 res[target].first.push_back("../octave/"+name+".cpp"); 98 res[target].second.push_back(cmd); 99 res["all"].first.push_back("../octave/"+name+".oct"); 100 } // end of OctaveMaterialPropertyInterface::getSpecificTargets 101 writeOutputFiles(const MaterialPropertyDescription & mpd,const FileDescription & fd) const102 void OctaveMaterialPropertyInterface::writeOutputFiles(const MaterialPropertyDescription& mpd, 103 const FileDescription& fd) const 104 { 105 using namespace std; 106 tfel::system::systemCall::mkdir("octave"); 107 const auto name = (mpd.material.empty()) ? mpd.className : mpd.material+"_"+mpd.className; 108 const auto fn = "octave/"+ name+".cpp"; 109 std::ofstream out{fn}; 110 tfel::raise_if(!out,"OctaveMaterialPropertyInterface::writeOutputFiles: " 111 "unable to open file '"+fn+"'"); 112 out.exceptions(std::ios::badbit|std::ios::failbit); 113 out << "/*!\n" 114 << "* \\file " << fn << '\n' 115 << "* \\brief " << "this file implements the " 116 << name << " MaterialLaw.\n" 117 << "* File generated by " 118 << MFrontHeader::getVersionName() << " " 119 << "version " << MFrontHeader::getVersionNumber() 120 << '\n'; 121 if(!fd.authorName.empty()){ 122 out << "* \\author " << fd.authorName << '\n'; 123 } 124 if(!fd.date.empty()){ 125 out << "* \\date " << fd.date << '\n'; 126 } 127 if(!fd.description.empty()){ 128 out << fd.description << '\n'; 129 } 130 out << " */\n\n" 131 << "#include<algorithm>\n" 132 << "#include<iostream>\n" 133 << "#include<iterator>\n" 134 << "#include<sstream>\n" 135 << "#include<fstream>\n" 136 << "#include<cstring>\n" 137 << "#include<cstdlib>\n" 138 << "#include<string>\n" 139 << "#include<vector>\n" 140 << "#include<cmath>\n" 141 << "#include<octave/oct.h>\n\n"; 142 if(!mpd.includes.empty()){ 143 out << mpd.includes << "\n\n"; 144 } 145 writeMaterialPropertyParametersHandler(out,mpd,name,"double","octave"); 146 writeExportDirectives(out); 147 out << "#ifdef __cplusplus\n" 148 << "extern \"C\"{\n" 149 << "#endif /* __cplusplus */\n\n"; 150 // mfront metadata 151 writeEntryPointSymbol(out,name); 152 writeTFELVersionSymbol(out,name); 153 writeInterfaceSymbol(out,name,"Octave"); 154 writeMaterialSymbol(out,name,mpd.material); 155 writeMaterialKnowledgeTypeSymbol(out,name,MATERIALPROPERTY); 156 157 if(mpd.inputs.size()>1){ 158 out << "static double\n" 159 << "get_scalar_value(const octave_value& value,const octave_idx_type,const octave_idx_type){\n" 160 << "return value.scalar_value();\n" 161 << "} // end of get_scalar_value\n\n" 162 << "static double\n" 163 << "get_matrix_value(const octave_value& value,const octave_idx_type i,const octave_idx_type j){\n" 164 << "return (value.matrix_value())(i,j);\n" 165 << "} // end of get_matrix_value\n\n" 166 << "static double\n" 167 << "get_range_value(const octave_value& value,const octave_idx_type,const octave_idx_type j){\n" 168 << "return value.range_value().elem(j);\n" 169 << "} // end of get_range_value\n\n"; 170 } 171 out << "static double " << name <<"_compute("; 172 if(!mpd.inputs.empty()){ 173 for(auto p3=mpd.inputs.begin();p3!=mpd.inputs.end();){ 174 out << "const double " << p3->name; 175 if((++p3)!=mpd.inputs.end()){ 176 out << ","; 177 } 178 } 179 } else { 180 out << "void"; 181 } 182 out << ")\n{\n" 183 << "using namespace std;\n" 184 << "using real = double;\n"; 185 // material laws 186 writeMaterialLaws(out,mpd.materialLaws); 187 // static variables 188 writeStaticVariables(out,mpd.staticVars,fd.fileName); 189 // parameters 190 if(!mpd.parameters.empty()){ 191 writeAssignMaterialPropertyParameters(out,mpd,name, 192 "double","octave"); 193 } 194 out << "real " << mpd.output.name << ";\n" 195 << mpd.f.body 196 << "return " << mpd.output.name << ";\n" 197 << "} // end of " << name << "_compute\n\n"; 198 if((hasBounds(mpd.inputs))||(hasPhysicalBounds(mpd.inputs))){ 199 out << "static double " << name <<"_checkBounds("; 200 for(auto p3=mpd.inputs.begin();p3!=mpd.inputs.end();){ 201 out << "const double " << p3->name; 202 if((++p3)!=mpd.inputs.end()){ 203 out << ","; 204 } 205 } 206 out << ")\n{\n"; 207 out << "using namespace std;\n"; 208 if(hasPhysicalBounds(mpd.inputs)){ 209 out << "// treating physical bounds\n"; 210 for(const auto& i : mpd.inputs){ 211 if(!i.hasPhysicalBounds()){ 212 continue; 213 } 214 const auto& b = i.getPhysicalBounds(); 215 const auto nbr = CMaterialPropertyInterfaceBase::getVariableNumber(mpd,i.name); 216 if(b.boundsType==VariableBoundsDescription::LOWER){ 217 out << "if(" << i.name<< " < "<< b.lowerBound << "){\n"; 218 out << "error(\"" << name << " :" << i.name 219 << " is below its physical lower bound.\");\n"; 220 out << "return -" << nbr << ";\n"; 221 out << "}\n"; 222 } else if(b.boundsType==VariableBoundsDescription::UPPER){ 223 out << "if(" << i.name<< " > "<< b.upperBound << "){\n"; 224 out << "error(\"" << name << " : " << i.name 225 << " is over its physical upper bound.\");\n"; 226 out << "return -" << nbr << ";\n"; 227 out << "}\n"; 228 } else { 229 out << "if((" << i.name<< " < "<< b.lowerBound << ")||" 230 << "(" << i.name<< " > "<< b.upperBound << ")){\n"; 231 out << "if(" << i.name<< " < "<< b.lowerBound << "){\n"; 232 out << "error(\"" << name << " : " << i.name 233 << " is below its physical lower bound.\");\n"; 234 out << "}\n"; 235 out << "if(" << i.name<< " > "<< b.upperBound << "){\n"; 236 out << "error(\"" << name << " : " << i.name 237 << " is over its physical upper bound.\");\n"; 238 out << "}\n"; 239 out << "return -" << nbr << ";\n"; 240 out << "}\n"; 241 } 242 } 243 } 244 if(hasBounds(mpd.inputs)){ 245 out << "// treating standard bounds\n"; 246 for(const auto& i: mpd.inputs){ 247 if(!i.hasBounds()){ 248 continue; 249 } 250 const auto& b = i.getBounds(); 251 const auto nbr = CMaterialPropertyInterfaceBase::getVariableNumber(mpd,i.name); 252 if(b.boundsType==VariableBoundsDescription::LOWER){ 253 out << "if(" << i.name<< " < "<< b.lowerBound << "){\n"; 254 out << "const octave_value policy = get_global_value(" 255 << "\"OCTAVE_OUT_OF_BOUNDS_POLICY\", true);\n"; 256 out << "if(policy.is_defined()){\n"; 257 out << "if(policy.is_string()){\n"; 258 out << "string msg(\"" << name << " : " << i.name 259 << " is below its physical lower bound.\");\n"; 260 out << "if(policy.string_value()==\"STRICT\"){\n"; 261 out << "error(msg.c_str());\n"; 262 out << "return -" << nbr << ";\n"; 263 out << "} if(policy.string_value()==\"WARNING\"){\n"; 264 out << "octave_stdout << msg << \"\\n\";\n"; 265 out << "} else {\n"; 266 out << "return " << nbr << ";\n"; 267 out << "}\n"; 268 out << "} else {\n"; 269 out << "return " << nbr << ";\n"; 270 out << "}\n"; 271 out << "} else {\n"; 272 out << "return " << nbr << ";\n"; 273 out << "}\n"; 274 out << "return " << nbr << ";\n"; 275 out << "}\n"; 276 } else if(b.boundsType==VariableBoundsDescription::UPPER){ 277 out << "if(" << i.name<< " < "<< b.lowerBound << "){\n"; 278 out << "const octave_value policy = get_global_value(" 279 << "\"OCTAVE_OUT_OF_BOUNDS_POLICY\", true);\n"; 280 out << "if(policy.is_defined()){\n"; 281 out << "if(policy.is_string()){\n"; 282 out << "string msg(\"" << i.name 283 << " is over its physical upper bound.\");\n"; 284 out << "if(policy.string_value()==\"STRICT\"){\n"; 285 out << "error(msg.c_str());\n"; 286 out << "return -" << nbr << ";\n"; 287 out << "} if(policy.string_value()==\"WARNING\"){\n"; 288 out << "octave_stdout << msg << \"\\n\";\n"; 289 out << "} else {\n"; 290 out << "return " << nbr << ";\n"; 291 out << "}\n"; 292 out << "} else {\n"; 293 out << "return " << nbr << ";\n"; 294 out << "}\n"; 295 out << "} else {\n"; 296 out << "return " << nbr << ";\n"; 297 out << "}\n"; 298 out << "return " << nbr << ";\n"; 299 out << "}\n"; 300 } else { 301 out << "if((" << i.name<< " < "<< b.lowerBound << ")||" 302 << "(" << i.name<< " > "<< b.upperBound << ")){\n"; 303 out << "const octave_value policy = get_global_value(" 304 << "\"OCTAVE_OUT_OF_BOUNDS_POLICY\", true);\n"; 305 out << "if(policy.is_defined()){\n"; 306 out << "if(policy.is_string()){\n"; 307 out << "string msg(\"" << i.name 308 << " is out of its bounds.\");\n"; 309 out << "if(policy.string_value()==\"STRICT\"){\n"; 310 out << "error(msg.c_str());\n"; 311 out << "return -" << nbr << ";\n"; 312 out << "} if(policy.string_value()==\"WARNING\"){\n"; 313 out << "octave_stdout << msg << \"\\n\";\n"; 314 out << "} else {\n"; 315 out << "return " << nbr << ";\n"; 316 out << "}\n"; 317 out << "} else {\n"; 318 out << "return " << nbr << ";\n"; 319 out << "}\n"; 320 out << "} else {\n"; 321 out << "return " << nbr << ";\n"; 322 out << "}\n"; 323 out << "return " << nbr << ";\n"; 324 out << "}\n"; 325 } 326 } 327 } 328 out << "return 0;\n} // end of " << name << "_checkBounds\n\n"; 329 } 330 331 out << "DEFUN_DLD(" << name << ",args,nargout,\n" 332 << "\"this function implements the " << name << " material law.\\n\"\n"; 333 if(!fd.description.empty()){ 334 const auto desc = tokenize(fd.description,'\n'); 335 for(auto p=desc.begin();p!=desc.end();){ 336 out << "\"" << this->treatDescriptionString(*p) << "\\n\""; 337 if(++p!=desc.end()){ 338 out << "\n"; 339 } 340 } 341 } 342 out << ")\n"; 343 out << "{\n"; 344 if(mpd.inputs.size()>1){ 345 //* out << name << "OctaveVarType varTypes[" 346 //* << mpd.inputs.size() << "];\n"; 347 out << "typedef double (*PtrGetFunction)(const octave_value&,const octave_idx_type,const octave_idx_type);\n"; 348 out << "PtrGetFunction getfunction[" << mpd.inputs.size() << "];\n"; 349 } 350 if(mpd.inputs.size()>=1){ 351 out << "// local variables used to convert ranges to matrices\n"; 352 out << "octave_idx_type i,j;\n"; 353 out << "octave_idx_type row = -1;\n"; 354 out << "octave_idx_type col = -1;\n"; 355 out << "bool areAllVariablesMatrices = true;\n"; 356 out << "bool areAllVariablesScalars = true;\n"; 357 } 358 out << "octave_value retval;\n"; 359 if(!mpd.parameters.empty()){ 360 const auto hn = getMaterialPropertyParametersHandlerClassName(name); 361 out << "if(!octave::" << hn << "::get" << hn << "().ok){\n" 362 << "error(octave::"<< name << "MaterialPropertyHandler::get" 363 << name << "MaterialPropertyHandler().msg.c_str());\n" 364 << "return retval;\n" 365 << "}\n"; 366 } 367 out << "if(args.length()!=" << mpd.inputs.size() << "){\n"; 368 if(mpd.inputs.size()==0){ 369 out << "error(\"" << name << " : no argument required\");\n"; 370 } else if(mpd.inputs.size()==1){ 371 out << "error(\"" << name << " : one argument and only one required\");\n"; 372 } else { 373 out << "error(\"" << name << " : " << mpd.inputs.size() 374 << " arguments and only " << mpd.inputs.size() << " arguments required\");\n"; 375 } 376 out << "return retval;\n"; 377 out << "}\n"; 378 if(mpd.inputs.size()==0){ 379 out << "retval = " << name << "_compute();\n"; 380 } else if(mpd.inputs.size()==1){ 381 out << "if(args(0).is_real_scalar()){\n"; 382 if((hasBounds(mpd.inputs))|| 383 (hasPhysicalBounds(mpd.inputs))){ 384 out << "if("<< name << "_checkBounds("; 385 out << "args(0).scalar_value())<0){\n"; 386 out << "return retval;\n"; 387 out << "}\n"; 388 } 389 out << "retval = " << name << "_compute("; 390 out << "args(0).scalar_value());\n"; 391 out << "} else if(args(0).is_real_matrix()||args(0).is_range()){\n"; 392 out << "Matrix xin0(args(0).is_range() ? args(0).range_value().matrix_value() : args(0).matrix_value());\n"; 393 out << "Matrix xout(xin0.rows(),xin0.cols());\n"; 394 out << "for(i=0;i!=xin0.rows();++i){\n"; 395 out << "for(j=0;j!=xin0.cols();++j){\n"; 396 if((hasBounds(mpd.inputs))|| 397 (hasPhysicalBounds(mpd.inputs))){ 398 out << "if("<< name << "_checkBounds("; 399 out << "xin0(i,j))<0){\n"; 400 out << "return retval;\n"; 401 out << "}\n"; 402 } 403 out << "xout(i,j) = " << name << "_compute("; 404 out << "xin0(i,j));\n"; 405 out << "}\n"; 406 out << "}\n"; 407 out << "retval = xout;\n"; 408 out << "} else {"; 409 out << "error(\"" << name << " : argument must be either a matrix or scalar\");\n"; 410 out << "return retval;\n"; 411 out << "}\n"; 412 } else { 413 // scalars 414 out << "for(i=0;i!=" << mpd.inputs.size() << ";++i){\n" 415 << "if(args(i).is_real_scalar()){\n" 416 << "areAllVariablesMatrices = false;\n" 417 << "getfunction[i] = &get_scalar_value;\n"; 418 // matrices 419 out << "} else if(args(i).is_real_matrix()){\n" 420 << "areAllVariablesScalars = false;\n" 421 << "getfunction[i] = &get_matrix_value;\n" 422 << "if(row==-1){\n" 423 << "row = args(i).matrix_value().rows();\n" 424 << "col = args(i).matrix_value().cols();\n" 425 << "} else {\n" 426 << "if((row!=args(i).matrix_value().rows())||\n" 427 << "(col!=args(i).matrix_value().cols())){\n" 428 << "error(\"" << name << " : all arguments shall have the same size\");\n" 429 << "return retval;\n" 430 << "}\n" 431 << "}\n"; 432 //ranges 433 out << "} else if(args(i).is_range()){\n" 434 << "areAllVariablesScalars = false;\n" 435 << "areAllVariablesMatrices = false;\n" 436 << "getfunction[i] = &get_range_value;\n" 437 << "if(row==-1){\n" 438 << "row = 1;\n" 439 << "col = args(i).range_value().nelem();\n" 440 << "} else {\n" 441 << "if((row!=1)||(col!=args(i).range_value().nelem())){\n" 442 << "error(\"" << name << " : all arguments shall have the same size\");\n" 443 << "return retval;\n" 444 << "}\n" 445 << "}\n"; 446 // unsupported type 447 out << "} else {\n" 448 << "error(\"" << name << " : arguments must be either a matrix or scalar\");\n" 449 << "return retval;\n" 450 << "}\n" 451 << "}\n"; 452 // all scalar case 453 out << "if(areAllVariablesScalars){\n"; 454 decltype(mpd.inputs.size()) i; 455 if((hasBounds(mpd.inputs))|| 456 (hasPhysicalBounds(mpd.inputs))){ 457 out << "if("<< name << "_checkBounds("; 458 for(i=0;i!=mpd.inputs.size()-1;++i){ 459 out << "args(" << i << ").scalar_value(),"; 460 } 461 out << "args(" << i << ").scalar_value())<0){\n"; 462 out << "return retval;\n"; 463 out << "}\n"; 464 } 465 out << "retval = " << name << "_compute("; 466 for(i=0;i!=mpd.inputs.size()-1;++i){ 467 out << "args(" << i << ").scalar_value(),"; 468 } 469 out << "args(" << i << ").scalar_value());\n"; 470 // all matrices case 471 out << "} else if(areAllVariablesMatrices){\n"; 472 for(i=0;i!=mpd.inputs.size();++i){ 473 out << "Matrix xin" << i<< "(args(" << i << ").matrix_value());\n"; 474 } 475 out << "Matrix xout(row,col);\n"; 476 out << "for(i=0;i!=row;++i){\n"; 477 out << "for(j=0;j!=col;++j){\n"; 478 if((hasBounds(mpd.inputs))|| 479 (hasPhysicalBounds(mpd.inputs))){ 480 out << "if("<< name << "_checkBounds("; 481 for(i=0;i!=mpd.inputs.size()-1;++i){ 482 out << "xin" << i << "(i,j),"; 483 } 484 out << "xin" << i << "(i,j))<0){\n"; 485 out << "return retval;\n"; 486 out << "}\n"; 487 } 488 out << "xout(i,j) = " << name << "_compute("; 489 for(i=0;i!=mpd.inputs.size()-1;++i){ 490 out << "xin" << i << "(i,j),"; 491 } 492 out << "xin" << i << "(i,j));\n"; 493 out << "}\n"; 494 out << "}\n"; 495 out << "retval = xout;\n"; 496 // general case 497 out << "} else {\n"; 498 out << "Matrix xout(row,col);\n"; 499 out << "for(i=0;i!=row;++i){\n"; 500 out << "for(j=0;j!=col;++j){\n"; 501 if((hasBounds(mpd.inputs))||(hasPhysicalBounds(mpd.inputs))){ 502 out << "if("<< name << "_checkBounds("; 503 for(i=0;i!=mpd.inputs.size()-1;++i){ 504 out << "(*(getfunction[" << i << "]))(args(" << i << "),i,j),\n"; 505 } 506 out << "(*(getfunction[" << i << "]))(args(" << i << "),i,j))<0){\n"; 507 out << "return retval;\n"; 508 out << "}\n"; 509 } 510 out << "xout(i,j) = " << name << "_compute("; 511 for(i=0;i!=mpd.inputs.size()-1;++i){ 512 out << "(*(getfunction[" << i << "]))(args(" << i << "),i,j),\n"; 513 } 514 out << "(*(getfunction[" << i << "]))(args(" << i << "),i,j));\n"; 515 out << "}\n"; 516 out << "}\n"; 517 out << "retval = xout;\n"; 518 out << "}\n"; 519 } 520 out << "return retval;\n"; 521 out << "} // end of " << name << "\n\n"; 522 out << "#ifdef __cplusplus\n"; 523 out << "} // end of extern \"C\"\n"; 524 out << "#endif /* __cplusplus */\n\n"; 525 out.close(); 526 } // end of OctaveMaterialPropertyInterface::writeOutputFiles 527 528 void replace(std::string & s,const std::string::value_type a,const std::string::value_type b)529 OctaveMaterialPropertyInterface::replace(std::string& s, 530 const std::string::value_type a, 531 const std::string::value_type b) 532 { 533 std::string::size_type p = 0; 534 while((p=s.find(a))!=std::string::npos){ 535 s[p]=b; 536 } 537 } // end of MFrontCLawInterfaceBase::replace 538 539 std::string treatDescriptionString(const std::string & s) const540 OctaveMaterialPropertyInterface::treatDescriptionString(const std::string& s) const 541 { 542 std::string res(s); 543 OctaveMaterialPropertyInterface::replace(res,'"',' '); 544 return res; 545 } // end of OctaveMaterialPropertyInterface::treatDescriptionString 546 547 OctaveMaterialPropertyInterface::~OctaveMaterialPropertyInterface() = default; 548 549 } // end of namespace mfront 550