diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp index 3e2202f19c1fe14c935077bf393e17376278817d..c15a5df6bd94174801ab8b1c4e098e2a45f23421 100644 --- a/dG3D/src/dG3DIPVariable.cpp +++ b/dG3D/src/dG3DIPVariable.cpp @@ -2033,47 +2033,62 @@ double torchANNBasedDG3DIPVariable::get(const int comp) const { #if defined(HAVE_TORCH) auto Vars_a = _internalVars.accessor<float,3>(); - if (comp == IPField::USER0) + if (comp == IPField::USER0) { - if (_n > 1) return Vars_a[0][0][0]; + if (_extraVariables.size2() > 0) return _extraVariables(0,0); else return 0.; } else if (comp == IPField::USER1) { - if (_n > 2) return Vars_a[0][0][1]; + if (_extraVariables.size2() > 1) return _extraVariables(0,1); else return 0.; } else if (comp == IPField::USER2) + { + if (_extraVariables.size2() > 2) return _extraVariables(0,2); + else return 0.; + } + else if (comp == IPField::USER3) + { + if (_n > 1) return Vars_a[0][0][0]; + else return 0.; + } + else if (comp == IPField::USER4) + { + if (_n > 2) return Vars_a[0][0][1]; + else return 0.; + } + else if (comp == IPField::USER5) { if (_n > 3) return Vars_a[0][0][2]; else return 0.; } - else if (comp == IPField::USER3) + else if (comp == IPField::USER6) { if (_n > 4) return Vars_a[0][0][3]; else return 0.; } - else if (comp == IPField::USER4) + else if (comp == IPField::USER7) { if (_n > 5) return Vars_a[0][0][4]; else return 0.; } - else if (comp == IPField::USER5) + else if (comp == IPField::USER8) { if (_n > 6) return Vars_a[0][0][5]; else return 0.; } - else if (comp == IPField::USER6) + else if (comp == IPField::USER9) { if (_n > 7) return Vars_a[0][0][6]; else return 0.; } - else if (comp == IPField::USER7) + else if (comp == IPField::USER10) { if (_n > 8) return Vars_a[0][0][7]; else return 0.; } - else if (comp == IPField::USER8) + /*else if (comp == IPField::USER8) { if (_n > 9) return Vars_a[0][0][8]; else return 0.; @@ -2087,7 +2102,7 @@ double torchANNBasedDG3DIPVariable::get(const int comp) const { if (_n > 11) return Vars_a[0][0][10]; else return 0.; - } + }*/ else { return dG3DIPVariable::get(comp); diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp index ff0d1b378d6d1935f76c733eee44fc688f519871..3d63ef77f3e7d5d76d7b55eef377f18632e5c564 100644 --- a/dG3D/src/dG3DMaterialLaw.cpp +++ b/dG3D/src/dG3DMaterialLaw.cpp @@ -22,6 +22,7 @@ #include "mlawElecSMP.h" #include <iostream> #include <fstream> +#include <sstream> #include <iomanip> #include <vector> #include <string> @@ -1900,6 +1901,103 @@ double ANNBasedDG3DMaterialLaw::soundSpeed() const }; // +void ExtraVariablesAtGaussPoints::fill(const std::string fname) +{ + Msg::Info("read %s",fname.c_str()); + filled=true; + locations.clear(); + field.clear(); + + std::ifstream inputFile(fname.c_str()); + std::vector<double> numbers; + if (inputFile.is_open()) + { + std::string line; + while (std::getline(inputFile, line)) + { + std::stringstream ss(line); + std::string token; + numbers.clear(); + while (std::getline(ss, token, ';')) { + double num; + std::istringstream(token) >> num; + numbers.push_back(num); + } + int numberSize=numbers.size(); + //Msg::Info("Line of %d values",numberSize); + if(numberSize<4) + { + Msg::Error("Not enough numbers per line in %s",fname.c_str()); + Msg::Exit(0); + } + SVector3 pt(numbers[0],numbers[1],numbers[2]); + locations.push_back(pt); + std::vector<double> localfield=std::vector<double>(numbers.begin()+3,numbers.end()); + field.push_back(localfield); + } + } + else + { + Msg::Error("Unable to read %s",fname.c_str()); + Msg::Exit(0); + } + + inputFile.close(); +} + +const std::vector<double> &ExtraVariablesAtGaussPoints::getFieldAtClosestLocation(const SVector3 &gp) const +{ + std::vector<SVector3>::const_iterator it1min, it1= locations.begin(); + std::vector< std::vector<double> >::const_iterator it2min, it2 = field.begin(); + if (locations.size() == field.size()) { + double minimumdist = 1.e12; + SVector3 dist; + for (; it1 != locations.end() && it2 != field.end(); ++it1, ++it2) { + dist=gp; + dist-=*it1; + double norm=dist.norm(); + if(norm <minimumdist) + { + it1min=it1; + it2min=it2; + minimumdist=norm; + } + } + //Msg::Info("Values at Gauss point location: %f %f %f", (*it1min)[0],(*it1min)[1],(*it1min)[2]); + //std::vector<double>::const_iterator it3; + //for (it3=it2min->begin(); it3 != it2min->end(); ++it3) { + // Msg::Info("Element %f",*it3); + //} + } + else + { + Msg::Error("ExtraVariablesAtGaussPoints locations and field not of the same size"); + Msg::Exit(0); + } + return *it2min; +} + +void ExtraVariablesAtGaussPoints::print() const +{ + if (locations.size() == field.size()) { + std::vector<SVector3>::const_iterator it1 = locations.begin(); + std::vector< std::vector<double> >::const_iterator it2 = field.begin(); + for (; it1 != locations.end() && it2 != field.end(); ++it1, ++it2) { + Msg::Info("Values at Gauss point location: %f %f %f", (*it1)[0],(*it1)[1],(*it1)[2]); + std::vector<double>::const_iterator it3; + for (it3=it2->begin(); it3 != it2->end(); ++it3) { + Msg::Info("Element %f",*it3); + } + } + } + else + { + Msg::Error("ExtraVariablesAtGaussPoints locations and field not of the same size"); + Msg::Exit(0); + } +} + + torchANNBasedDG3DMaterialLaw::torchANNBasedDG3DMaterialLaw(const int num, const double rho, const int numberOfInput, const int numInternalVars, const char* nameTorch, const double EXXmean, const double EXXstd, const double EXYmean, const double EXYstd, const double EYYmean, const double EYYstd, const double EYZmean, const double EYZstd, const double EZZmean, @@ -1910,8 +2008,10 @@ torchANNBasedDG3DMaterialLaw::torchANNBasedDG3DMaterialLaw(const int num, const _EXXmean(EXXmean), _EXXstd(EXXstd), _EXYmean(EXYmean), _EXYstd(EXYstd), _EYYmean(EYYmean), _EYYstd(EYYstd), _EYZmean(EYZmean), _EYZstd(EYZstd), _EZZmean(EZZmean), _EZZstd(EZZstd), _EZXmean(EZXmean), _EZXstd(EZXstd), _SXXmean(SXXmean), _SXXstd(SXXstd), _SXYmean(SXYmean), _SXYstd(SXYstd), _SYYmean(SYYmean), _SYYstd(SYYstd), _SYZmean(SYZmean), _SYZstd(SYZstd), _SZZmean(SZZmean), _SZZstd(SZZstd), - _SZXmean(SZXmean), _SZXstd(SZXstd), _tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL), _NeedExtraNorm(false), _DoubleInput(false), - _numberOfExtraInput(0), _time_arg(-1), _ctime(0.0) + _SZXmean(SZXmean), _SZXstd(SZXstd), _tangentByPerturbation(pert), _pertTol(tol), _kinematicInput(EGL), _NeedExtraNorm(false), _DoubleInput(false), + _numberOfExtraInput(0), + _fileExtraInputsGp(""), _extraVariablesAtGaussPoints(), _time_arg(-1), _ctime(0.0) + { #if defined(HAVE_TORCH) try{ @@ -1976,9 +2076,12 @@ torchANNBasedDG3DMaterialLaw::torchANNBasedDG3DMaterialLaw(const torchANNBasedDG _EYZmean(src._EYZmean), _EYZstd(src._EYZstd), _EZZmean(src._EZZmean), _EZZstd(src._EZZstd), _EZXmean(src._EZXmean), _EZXstd(src._EZXstd), _SXXmean(src._SXXmean), _SXXstd(src._SXXstd), _SXYmean(src._SXYmean), _SXYstd(src._SXYstd), _SYYmean(src._SYYmean), _SYYstd(src._SYYstd), _SYZmean(src._SYZmean), _SYZstd(src._SYZstd), _SZZmean(src._SZZmean), _SZZstd(src._SZZstd), _SZXmean(src._SZXmean), _SZXstd(src._SZXstd), - _tangentByPerturbation(src._tangentByPerturbation), _pertTol(src._pertTol), _kinematicInput(src._kinematicInput), _NeedExtraNorm(src._NeedExtraNorm), - _DoubleInput(src._DoubleInput), _numberOfExtraInput(src._numberOfExtraInput), _initialValue_ExtraInput(src._initialValue_ExtraInput), - _normParams_ExtraInput(src._normParams_ExtraInput), _time_arg(src._time_arg), _ctime(src._ctime) + _tangentByPerturbation(src._tangentByPerturbation), _pertTol(src._pertTol), _kinematicInput(src._kinematicInput), + _NeedExtraNorm(src._NeedExtraNorm), + _DoubleInput(src._DoubleInput), + _numberOfExtraInput(src._numberOfExtraInput), _initialValue_ExtraInput(src._initialValue_ExtraInput), + _normParams_ExtraInput(src._normParams_ExtraInput), _fileExtraInputsGp(src._fileExtraInputsGp), + _extraVariablesAtGaussPoints(src._extraVariablesAtGaussPoints), _time_arg(src._time_arg), _ctime(src._ctime) { #if defined(HAVE_TORCH) module = src.module; @@ -2152,10 +2255,23 @@ void torchANNBasedDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBody // check interface element bool inter=true; const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele); + std::vector<double> value_ExtraInput=_initialValue_ExtraInput; + SPoint3 ptGlobal; + ele->pnt(GP->pt[0],GP->pt[1],GP->pt[2],ptGlobal); + SVector3 pt(ptGlobal); + for(int i=0; i<value_ExtraInput.size();i++) + { + if(_extraVariablesAtGaussPoints.isStored(i)) + { + value_ExtraInput[i]=_extraVariablesAtGaussPoints.returnFieldAtClosestLocation(i,pt); + } + } + //Msg::Info("value_ExtraInput %f %f %f", value_ExtraInput[0], value_ExtraInput[1],value_ExtraInput[2]); + if(iele==NULL) inter=false; - IPVariable* ipvi = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, _initialValue_ExtraInput); - IPVariable* ipv1 = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, _initialValue_ExtraInput); - IPVariable* ipv2 = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, _initialValue_ExtraInput); + IPVariable* ipvi = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, value_ExtraInput); + IPVariable* ipv1 = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, value_ExtraInput); + IPVariable* ipv2 = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, value_ExtraInput); if(ips != NULL) delete ips; ips = new IP3State(state_,ipvi,ipv1,ipv2); } @@ -2168,7 +2284,19 @@ void torchANNBasedDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBo if(iele == NULL){ inter=false; } - ipv = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, _initialValue_ExtraInput); + std::vector<double> value_ExtraInput=_initialValue_ExtraInput; + SPoint3 ptGlobal; + ele->pnt(GP->pt[0],GP->pt[1],GP->pt[2],ptGlobal); + SVector3 pt(ptGlobal); + for(int i=0; i<value_ExtraInput.size();i++) + { + if(_extraVariablesAtGaussPoints.isStored(i)) + { + value_ExtraInput[i]=_extraVariablesAtGaussPoints.returnFieldAtClosestLocation(i,pt); + } + } + //Msg::Info("value_ExtraInput %f %f %f", value_ExtraInput[0], value_ExtraInput[1],value_ExtraInput[2]); + ipv = new torchANNBasedDG3DIPVariable(_numberOfInternalVariables,_initialValue_h,hasBodyForce,inter, _numberOfExtraInput, value_ExtraInput); } void torchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent) @@ -2269,6 +2397,7 @@ void torchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipv // Set the current time @Mohib Extra1(0, _time_arg) = _ctime; } + //Msg::Info("value Extra1 %f %f %f", Extra1(0,0), Extra1(0,1),Extra1(0,2)); static fullMatrix<double> S(1,6), DSDE(6,6); const torch::Tensor& h0 = ipvprev->getConstRefToInternalVariables(); diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h index 2f0402cab8fade841449457847022867c3ef8daf..70a6ab8e944465ba96788be79e2ebf3baba2fcc8 100644 --- a/dG3D/src/dG3DMaterialLaw.h +++ b/dG3D/src/dG3DMaterialLaw.h @@ -467,6 +467,53 @@ class ANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ #endif //SWIG }; +class ExtraVariablesAtGaussPoints +{ + bool filled; + std::vector<SVector3> locations; + std::vector< std::vector<double> > field; + std::map <int,int> _storeExtraInputAtGp; //map between location in torch input and stored location in _extraVariablesAtGaussPoints + public: + ExtraVariablesAtGaussPoints():filled(false), field(), _storeExtraInputAtGp() {} + ExtraVariablesAtGaussPoints(const ExtraVariablesAtGaussPoints& src):filled(src.filled), locations(src.locations), field(src.field), + _storeExtraInputAtGp(src._storeExtraInputAtGp) {} + ExtraVariablesAtGaussPoints& operator = (const ExtraVariablesAtGaussPoints& src) + { + filled=src.filled; locations=src.locations; field=src.field; _storeExtraInputAtGp=src._storeExtraInputAtGp; + return *this; + } + ~ExtraVariablesAtGaussPoints(){} + void fill(const std::string fname); + const std::vector<double> &getFieldAtClosestLocation(const SVector3 &gp) const; + void print() const; + void fillExtraInputAtGp(int nb1, int nb2) + { + if(_storeExtraInputAtGp.find(nb1) == _storeExtraInputAtGp.end()) + _storeExtraInputAtGp[nb1]=nb2; + else + Msg::Error("ExtraVariablesAtGaussPoints::fillExtraInputAtGp already defined"); + } + bool isStored(int nb1) const + { + if(_storeExtraInputAtGp.find(nb1) == _storeExtraInputAtGp.end()) + return false; + else + return true; + } + double returnFieldAtClosestLocation(int nb1, const SVector3 &gp) const + { + if(!isStored(nb1)) + Msg::Error("ExtraVariablesAtGaussPoints::returnFieldAtClosestLocation field %d is not to be stored",nb1); + std::map <int,int>::const_iterator it=_storeExtraInputAtGp.find(nb1); + const std::vector<double> & field=getFieldAtClosestLocation(gp); + if(it->second>field.size()) + Msg::Error("ExtraVariablesAtGaussPoints::returnFieldAtClosestLocation field %d is not stored at location %d",nb1,it->second); + return field[it->second]; + } + +}; + + class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ public: enum KinematicInput{EGL=1,U=2,smallStrain=3,EGLInc=4,UInc=5,smallStrainInc=6}; // EGL - Green Lagraneg strain, U-Biot strain, smallStrain=small strain; Inc for increments @@ -481,6 +528,8 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ int _numberOfExtraInput; // Number of additional inputs other than 6 (or 3 for 2D) strain components @Mohib std::vector<double> _initialValue_ExtraInput; // initial values for the extra inputs @Mohib std::vector<double> _normParams_ExtraInput; // Container for additional inputs, it is needed for normalization @Mohib + std::string _fileExtraInputsGp; + ExtraVariablesAtGaussPoints _extraVariablesAtGaussPoints; double _ctime; // variable for storing current time @Mohib #if defined(HAVE_TORCH) @@ -557,6 +606,22 @@ class torchANNBasedDG3DMaterialLaw : public dG3DMaterialLaw{ // Setter for storing norm params of Extra Inputs @Mohib void setNormExtraInp(const double mean, const double std); + + + void setFileExtraInputsGp(const std::string fname) + { + _fileExtraInputsGp=fname; + _extraVariablesAtGaussPoints.fill(fname); + //_extraVariablesAtGaussPoints.print(); + //const std::vector<double> &rf=_extraVariablesAtGaussPoints.getFieldAtClosestLocation(SVector3(25,28.873,28.873)); + //double val=_extraVariablesAtGaussPoints.returnFieldAtClosestLocation(0,SVector3(25,28.873,28.873)); + //Msg::Info(" Field 0 at GP is %f",val); + } + + void fillExtraInputAtGp(int nb1, int nb2) + { + _extraVariablesAtGaussPoints.fillExtraInputAtGp(nb1,nb2); + } #ifndef SWIG torchANNBasedDG3DMaterialLaw(const torchANNBasedDG3DMaterialLaw& src); virtual ~torchANNBasedDG3DMaterialLaw();