diff --git a/NonLinearSolver/CMakeLists.txt b/NonLinearSolver/CMakeLists.txt
index 250f128d614b05e7f753dd59639d353a2e5513c7..e19b96bd08b229d182e6ed05d4de50f058dfc483 100644
--- a/NonLinearSolver/CMakeLists.txt
+++ b/NonLinearSolver/CMakeLists.txt
@@ -19,6 +19,11 @@ else()
set(_DEBUG FALSE)
endif(CMAKE_BUILD_TYPE MATCHES "debug")
+option(ENABLE_CG_MPI_INTERFACE "use mpi interface between CG domains" ON)
+if(ENABLE_CG_MPI_INTERFACE)
+ set(HAVE_CG_MPI_INTERFACE TRUE)
+endif(ENABLE_CG_MPI_INTERFACE)
+
set(CMAKE_C_FLAGS " ${CMAKE_C_FLAGS} -DNONLOCALGMSH")
set(CMAKE_CXX_FLAGS " ${CMAKE_CXX_FLAGS} -DNONLOCALGMSH")
#not compatible with ABI option of torch
diff --git a/NonLinearSolver/NonLinearSolverConfig.h.in b/NonLinearSolver/NonLinearSolverConfig.h.in
index f6020412dc710f546b3388d010b3564e32f021a3..1af20012d61af6c483ce35fe8af9092d3d20bceb 100644
--- a/NonLinearSolver/NonLinearSolverConfig.h.in
+++ b/NonLinearSolver/NonLinearSolverConfig.h.in
@@ -7,7 +7,7 @@
#cmakedefine ENABLE_DG3D
#cmakedefine ENABLE_MSCH
#cmakedefine HAVE_TORCH
-
+#cmakedefine HAVE_CG_MPI_INTERFACE
#endif
diff --git a/NonLinearSolver/internalPoints/CMakeLists.txt b/NonLinearSolver/internalPoints/CMakeLists.txt
index c984cad0ab6cbb6da2921929215c6e1587d93c26..83a0969341c50a3c734b7cbb0757e9641d220707 100644
--- a/NonLinearSolver/internalPoints/CMakeLists.txt
+++ b/NonLinearSolver/internalPoints/CMakeLists.txt
@@ -48,6 +48,7 @@ set(SRC
ipUMATInterface.cpp
ipCrystalPlasticity.cpp
ipGursonUMAT.cpp
+ ipVEVPUMAT.cpp
ipAnIsotropicElecTherMech.cpp
ipHyperelastic.cpp
ipNonLocalDamageHyperelastic.cpp
diff --git a/NonLinearSolver/internalPoints/ipJ2linear.cpp b/NonLinearSolver/internalPoints/ipJ2linear.cpp
index b5d1aa6bf8567ad80215e8b64e460f5e8057b2e3..34361c9cc1ee9efbab6cb0aaa13af84e61aaa81f 100644
--- a/NonLinearSolver/internalPoints/ipJ2linear.cpp
+++ b/NonLinearSolver/internalPoints/ipJ2linear.cpp
@@ -13,7 +13,8 @@
IPJ2linear::IPJ2linear() : IPVariableMechanics(), _j2lepspbarre(0.), _j2lepsp(1.), _j2ldsy(0.), _elasticEnergy(0.),
_plasticPower(0.),_DplasticPowerDF(0.),_Ee(0.),_irreversibleEnergy(0.),_DirreversibleEnergyDF(0.),
_dissipationBlocked(false),_dissipationActive(false),_plasticEnergy(0.),
-_P(0.), _sig(0.), _svm(0.), _strain(0.), _eigenstress(0.), _eigenstress_eq(0.), _eigenstrain(0.)
+_P(0.), _sig(0.), _svm(0.), _strain(0.), _eigenstress(0.), _eigenstress_eq(0.), _eigenstrain(0.),
+_DplasticEnergyDF(0.)
{
ipvJ2IsotropicHardening=NULL;
Msg::Error("IPJ2Linear::IPJ2Linear is not initialized with a hardening IP Variable");
@@ -22,7 +23,8 @@ _P(0.), _sig(0.), _svm(0.), _strain(0.), _eigenstress(0.), _eigenstress_eq(0.),
IPJ2linear::IPJ2linear(const J2IsotropicHardening *j2IH) : IPVariableMechanics(), _j2lepspbarre(0.), _j2lepsp(1.), _j2ldsy(0.),
_elasticEnergy(0.),_plasticPower(0.),_DplasticPowerDF(0.),_Ee(0.),_irreversibleEnergy(0.),_DirreversibleEnergyDF(0.),
_dissipationBlocked(false),_dissipationActive(false),_plasticEnergy(0.),
-_P(0.), _sig(0.),_svm(0.), _strain(0.), _eigenstress(0.), _eigenstress_eq(0.), _eigenstrain(0.)
+_P(0.), _sig(0.),_svm(0.), _strain(0.), _eigenstress(0.), _eigenstress_eq(0.), _eigenstrain(0.),
+_DplasticEnergyDF(0.)
{
ipvJ2IsotropicHardening=NULL;
if(j2IH ==NULL) Msg::Error("IPJ2Linear::IPJ2Linear has no j2IH");
@@ -33,7 +35,8 @@ IPJ2linear::IPJ2linear(const IPJ2linear &source) : IPVariableMechanics(source),
_j2lepsp(source._j2lepsp), _j2ldsy(source._j2ldsy),
_elasticEnergy(source._elasticEnergy),
_plasticPower(source._plasticPower),_DplasticPowerDF(source._DplasticPowerDF), _Ee(source._Ee),_irreversibleEnergy(source._irreversibleEnergy), _DirreversibleEnergyDF(source._DirreversibleEnergyDF), _dissipationBlocked(source._dissipationBlocked),_dissipationActive(source._dissipationActive),_plasticEnergy(source._plasticEnergy),
-_P(source._P), _sig(source._sig), _svm(source._svm), _strain(source._strain), _eigenstress(source._eigenstress), _eigenstress_eq(source._eigenstress_eq), _eigenstrain(source._eigenstrain)
+_P(source._P), _sig(source._sig), _svm(source._svm), _strain(source._strain), _eigenstress(source._eigenstress), _eigenstress_eq(source._eigenstress_eq), _eigenstrain(source._eigenstrain),
+_DplasticEnergyDF(source._DplasticEnergyDF)
{
ipvJ2IsotropicHardening = NULL;
if(source.ipvJ2IsotropicHardening != NULL)
@@ -59,6 +62,7 @@ IPJ2linear& IPJ2linear::operator=(const IPVariable &source)
_dissipationBlocked = src->_dissipationBlocked;
_dissipationActive = src->_dissipationActive;
_plasticEnergy = src->_plasticEnergy;
+ _DplasticEnergyDF= src->_DplasticEnergyDF;
_P = src->_P;
_sig = src->_sig;
_svm = src->_svm;
@@ -263,6 +267,7 @@ void IPJ2linear::restart()
restartManager::restart(_dissipationBlocked);
restartManager::restart(_dissipationActive);
restartManager::restart(_plasticEnergy);
+ restartManager::restart(_DplasticEnergyDF);
restartManager::restart(_P);
restartManager::restart(_sig);
restartManager::restart(_svm);
diff --git a/NonLinearSolver/internalPoints/ipJ2linear.h b/NonLinearSolver/internalPoints/ipJ2linear.h
index 8d1120de22f12e0df3ffd670e9fa1c64af4ccab7..ede3b99b3475bbfceb08829667700ab5429d47b0 100644
--- a/NonLinearSolver/internalPoints/ipJ2linear.h
+++ b/NonLinearSolver/internalPoints/ipJ2linear.h
@@ -36,6 +36,7 @@ class IPJ2linear : public IPVariableMechanics
STensor3 _strain;
double _plasticEnergy;
+ STensor3 _DplasticEnergyDF;
double _plasticPower; // plastic power
STensor3 _DplasticPowerDF; // dplastic power DF
STensor3 _Ee; // elastic deformation
@@ -106,6 +107,9 @@ class IPJ2linear : public IPVariableMechanics
virtual double& getRefToPlasticEnergy() {return _plasticEnergy;};
virtual const double& getConstRefToPlasticEnergy() const {return _plasticEnergy;};
+
+ virtual STensor3& getRefToDPlasticEnergyDF() {return _DplasticEnergyDF;};
+ virtual const STensor3& getConstRefToDPlasticEnergyDF() const {return _DplasticEnergyDF;};
virtual double& getRefToPlasticPower() {return _plasticPower;};
virtual const double& getConstRefToPlasticPower() const {return _plasticPower;};
diff --git a/NonLinearSolver/internalPoints/ipUMATInterface.h b/NonLinearSolver/internalPoints/ipUMATInterface.h
index 148aea078607493a05e010bfb8fd30f1627c0bc5..cc4db217888149239ae64aebf6795c6cc464a87f 100644
--- a/NonLinearSolver/internalPoints/ipUMATInterface.h
+++ b/NonLinearSolver/internalPoints/ipUMATInterface.h
@@ -1,5 +1,5 @@
//
-// Description: storing class for umat interface
+// Description: storing class for umat interface
// Author: <L. NOELS>, (C) 2019
//
// Copyright: See COPYING file that comes with this distribution
@@ -32,7 +32,7 @@ class IPUMATInterface : public IPVariableMechanics
public:
IPUMATInterface(int _nsdv,double eleSize) : IPVariableMechanics(), _nsdv(_nsdv), _strain(6), _stress(6),
- _materialTensor(6,6),_elasticEne(0.),
+ _materialTensor(6,6),_elasticEne(0.),
_plasticEne(0.),_damageEne(0.), _elementSize(eleSize), _dissipationBlocked(false)
{
_statev = new double[_nsdv];
@@ -55,7 +55,7 @@ class IPUMATInterface : public IPVariableMechanics
{
IPVariableMechanics::operator=(_source);
const IPUMATInterface *source=static_cast<const IPUMATInterface *> (&_source);
- if(_nsdv !=source->_nsdv) Msg::Error("IPCrystalPlasticity ips do not have the same number of internal variables");
+ if(_nsdv !=source->_nsdv) Msg::Error("IPUMATInterface ips do not have the same number of internal variables");
for(int i = 0; i<_nsdv; i++) _statev[i]=source->_statev[i];
_strain = source->_strain;
_stress = source->_stress;
diff --git a/NonLinearSolver/internalPoints/ipVEVPUMAT.cpp b/NonLinearSolver/internalPoints/ipVEVPUMAT.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..61493e7e824af44deb2a6e4f48cf931398b8cf05
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipVEVPUMAT.cpp
@@ -0,0 +1,26 @@
+//
+// Description: storing class for vevp umat interface
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "ipVEVPUMAT.h"
+#include "ipField.h"
+#include "restartManager.h"
+
+
+double IPVEVPUMAT::get(int comp) const
+{
+ if(comp==IPField::PLASTICSTRAIN)
+ return _statev[18];
+ else
+ return IPUMATInterface::get(comp);
+}
+
+void IPVEVPUMAT::restart()
+{
+ IPUMATInterface::restart();
+ return;
+}
diff --git a/NonLinearSolver/internalPoints/ipVEVPUMAT.h b/NonLinearSolver/internalPoints/ipVEVPUMAT.h
new file mode 100644
index 0000000000000000000000000000000000000000..9e09306301e3a5434a6ee8c310a3cb4c9c410b83
--- /dev/null
+++ b/NonLinearSolver/internalPoints/ipVEVPUMAT.h
@@ -0,0 +1,46 @@
+//
+// Description: storing class for vevp umat interface
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef IPVEVPUMAT_H_
+#define IPVEVPUMAT_H_
+#include "ipUMATInterface.h"
+#include "STensor3.h"
+class IPVEVPUMAT : public IPUMATInterface
+{
+ protected:
+
+ public: // All data public to avoid the creation of function to access
+
+ public:
+ IPVEVPUMAT(int _nsdv, double eleSize) : IPUMATInterface(_nsdv, eleSize)
+ {
+
+ _statev[0]=1;
+ _statev[1]=1;
+ _statev[2]=1;
+ }
+ IPVEVPUMAT(const IPVEVPUMAT &source) : IPUMATInterface(source)
+ {
+ }
+ IPVEVPUMAT &operator = (const IPVariable &_source)
+ {
+ IPUMATInterface::operator=(_source);
+ return *this;
+ }
+ virtual ~IPVEVPUMAT()
+ {
+ }
+ virtual IPVariable* clone() const {return new IPVEVPUMAT(*this);};
+ virtual void restart();
+ virtual double get(const int i) const;
+
+ protected:
+ IPVEVPUMAT(){}
+};
+
+#endif // IPVEVPUMAT_H_
diff --git a/NonLinearSolver/materialLaw/CMakeLists.txt b/NonLinearSolver/materialLaw/CMakeLists.txt
index e2045dcc04444bf4fae5b06f1e08f30d2d0116b5..f0652edb953c17ee350e0eb36feee1c914817f8c 100644
--- a/NonLinearSolver/materialLaw/CMakeLists.txt
+++ b/NonLinearSolver/materialLaw/CMakeLists.txt
@@ -36,12 +36,14 @@ set(SRC
mlawVUMATinterface.cpp
mlawCrystalPlasticity.cpp
mlawGursonUMAT.cpp
+ mlawVEVPUMAT.cpp
additionalFortran.f
Ksub15.f
Kevpm15.f
ShearB.f
ALAMEL15.f
gurson.f
+ vevp.f
LD_utils.f
mlawAnisotropic.cpp
vumat.f
diff --git a/NonLinearSolver/materialLaw/mlaw.h b/NonLinearSolver/materialLaw/mlaw.h
index 817ff93c36706ad0c969cb12e21566b5bbe844ff..ceb1d394ed986f8defaf528b67f5cd00e679433f 100644
--- a/NonLinearSolver/materialLaw/mlaw.h
+++ b/NonLinearSolver/materialLaw/mlaw.h
@@ -35,7 +35,7 @@ class materialLaw{
ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic,nonLocalDamageGursonThermoMechanics,
localDamageJ2SmallStrain,nonLocalDamageJ2SmallStrain,cluster,tfa,ANN, DMN, torchANN, LinearElecMagTherMech, LinearElecMagInductor, hyperviscoelastic,
GenericThermoMechanics, ElecMagGenericThermoMechanics, ElecMagInductor,
- nonlineartvm,nonlinearTVE,nonlinearTVP,vevpmfh};
+ nonlineartvm,nonlinearTVE,nonlinearTVP,vevpmfh, VEVPUMAT};
protected :
diff --git a/NonLinearSolver/materialLaw/mlawJ2linear.cpp b/NonLinearSolver/materialLaw/mlawJ2linear.cpp
index 1c3888da3fb1ca67c04be6c05aff15e5e30b2bdd..f865ac2ae6d99b4058d309d392a95030c6cb48a3 100644
--- a/NonLinearSolver/materialLaw/mlawJ2linear.cpp
+++ b/NonLinearSolver/materialLaw/mlawJ2linear.cpp
@@ -891,6 +891,15 @@ void mlawJ2linear::predictorCorector(const STensor3& F0, // initial defo
else{
STensorOperation::zero(dplasticpowerDF);
}
+
+ STensor3& DplasticEnergyDF = q->getRefToDPlasticEnergyDF();
+ if (Deps > 0.){
+ DplasticEnergyDF = dpdF;
+ DplasticEnergyDF *= (Sy+H*Deps);
+ }
+ else{
+ STensorOperation::zero(DplasticEnergyDF);
+ }
// irreversible energy
STensor3& DirrEnegDF = q->getRefToDIrreversibleEnergyDF();
diff --git a/NonLinearSolver/materialLaw/mlawVEVPUMAT.cpp b/NonLinearSolver/materialLaw/mlawVEVPUMAT.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..605d2b7b0ff4f25dff648acdcb96625b99f1e4a0
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawVEVPUMAT.cpp
@@ -0,0 +1,120 @@
+//
+// C++ Interface: material law
+//
+// Description: UMAT interface for VEVP model
+//
+//
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#include <math.h>
+#include "MInterfaceElement.h"
+#include <fstream>
+#include "mlawVEVPUMAT.h"
+#include "STensorOperations.h"
+
+#if !defined(F77NAME)
+#define F77NAME(x) (x##_)
+#endif
+
+extern "C" {
+ void F77NAME(umat_vevp)(double *STRESS, double *STATEV, double *DDSDDE, double *SSE, double *SPD, double *SCD,
+ double *RPL, double *DDSDDT, double *DRPLDE, double *DRPLDT,
+ double *STRAN, double *DSTRAN, double *TIM, double *DTIME, double *TEMP, double *DTEMP, double *PREDEF, double *DPRED, const char *CMNAME,
+ int *NDI, int*NSHR, int* NTENS, int*NSTATV, double *PROPS, int *NPROPS, double *COORDS, double *DROT, double *PNEWDT,
+ double *CELENT, double *DFGRD0, double *DFGRD1, int *NOEL, int *NPT, int *LAYER, int *KSPT, int *KSTEP, int *KINC);
+}
+
+
+
+mlawVEVPUMAT::mlawVEVPUMAT(const int num, const char *propName) : mlawUMATInterface(num,0.,propName)
+{
+ //should depend on the umat: here we read the grains.inp file with 3 euler angle, 2 parameters equal to 1, and the number of internal variables (23)
+ std::ifstream in(propName);
+ if(!in) Msg::Error("Cannot open the file %s! Maybe is missing ",propName);
+ std::string line;
+ std::getline(in, line,'\n');
+ nsdv=atoi(line.c_str());
+ std::getline(in, line,'\n');
+ nprops1=atoi(line.c_str());
+ std::cout << "number of internal variables: "<< nsdv << "; number of properties: "<< nprops1 << std::endl;
+ props=new double[nprops1];
+ int i=0;
+ while (std::getline(in, line,'\n') and i<nprops1)
+ {
+ //std::cout << line << std::endl;
+ props[i]=atof(line.c_str());
+ std::cout << "property "<< i << ": "<< props[i] << std::endl;
+ i++;
+ }
+ in.close();
+
+ double Kinfty = props[1];
+ double Ginfty = props[2];
+ _rho = 0.;
+ _nu0 = (3*Kinfty - 2*Ginfty)/2./(3*Kinfty+Ginfty);
+ _mu0 = Ginfty;
+}
+mlawVEVPUMAT::mlawVEVPUMAT(const mlawVEVPUMAT &source) :
+ mlawUMATInterface(source)
+{
+
+}
+
+mlawVEVPUMAT& mlawVEVPUMAT::operator=(const materialLaw &source)
+{
+ mlawUMATInterface::operator=(source);
+ const mlawVEVPUMAT* src =static_cast<const mlawVEVPUMAT*>(&source);
+ return *this;
+};
+
+mlawVEVPUMAT::~mlawVEVPUMAT()
+{
+
+}
+void mlawVEVPUMAT::createIPState(IPVEVPUMAT *ivi, IPVEVPUMAT *iv1, IPVEVPUMAT *iv2) const
+{
+ mlawUMATInterface::createIPState(ivi, iv1, iv2);
+}
+
+void mlawVEVPUMAT::createIPVariable(IPVEVPUMAT *ipv, bool hasBodyForce, const MElement *ele,const int nbFF,const IntPt *GP, const int gpt) const
+{
+ mlawUMATInterface::createIPVariable(ipv, hasBodyForce, ele, nbFF,GP, gpt);
+}
+
+
+void mlawVEVPUMAT::callUMAT(double *stress, double *statev, double **ddsdde, double &sse, double &spd, double &scd, double &rpl,
+ double *ddsddt, double *drplde, double &drpldt, double *stran, double *dtsran,
+ double *tim, double timestep, double temperature, double deltaTemperature, double *predef, double *dpred,
+ const char *CMNAME, int &ndi, int &nshr, int tensorSize,
+ int statevSize, double *prs, int matPropSize, double *coords, double **dRot,
+ double &pnewdt, double &celent, double **F0, double **F1,
+ int &noel, int &npt, int &layer, int &kspt, int &kstep, int &kinc) const
+{
+
+
+ double dRtmp[9];
+ double F0tmp[9];
+ double F1tmp[9];
+
+ double ddsddetmp[36];
+
+ convert3x3To9((const double**)dRot, dRtmp);
+ convert3x3To9((const double**)F0, F0tmp);
+ convert3x3To9((const double**)F1, F1tmp);
+
+ convert6x6To36((const double**)ddsdde,ddsddetmp);
+
+ F77NAME(umat_vevp)(stress, statev, ddsddetmp, &sse, &spd, &scd, &rpl,
+ ddsddt, drplde, &drpldt, stran, dtsran,
+ tim, ×tep, &temperature, &deltaTemperature, predef, dpred,
+ CMNAME, &ndi, &nshr, &tensorSize, &statevSize, prs, &matPropSize, coords, dRtmp,
+ &pnewdt, &celent, F0tmp, F1tmp,
+ &noel, &npt, &layer, &kspt, &kstep, &kinc);
+
+ convert9To3x3((const double*)dRtmp,dRot);
+ convert36To6x6((const double*)ddsddetmp,ddsdde);
+}
+
diff --git a/NonLinearSolver/materialLaw/mlawVEVPUMAT.h b/NonLinearSolver/materialLaw/mlawVEVPUMAT.h
new file mode 100644
index 0000000000000000000000000000000000000000..affaabb1563ba9ffa3e5bcf51871b9290b1ecbbc
--- /dev/null
+++ b/NonLinearSolver/materialLaw/mlawVEVPUMAT.h
@@ -0,0 +1,55 @@
+//
+// C++ Interface: material law
+//
+// Description: UMAT interface for VEVP model
+//
+//
+// Author: <V.D. NGUYEN>, (C) 2023
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef MLAWVEVPUMAT_H_
+#define MLAWVEVPUMAT_H_
+#include "mlawUMATInterface.h"
+#include "ipVEVPUMAT.h"
+
+class mlawVEVPUMAT : public mlawUMATInterface
+{
+
+ public:
+ mlawVEVPUMAT(const int num, const char *propName);
+
+ #ifndef SWIG
+ mlawVEVPUMAT(const mlawVEVPUMAT &source);
+ mlawVEVPUMAT& operator=(const materialLaw &source);
+ virtual ~mlawVEVPUMAT();
+ virtual materialLaw* clone() const {return new mlawVEVPUMAT(*this);};
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{}; // do nothing
+ // function of materialLaw
+ virtual matname getType() const{return materialLaw::VEVPUMAT;}
+
+ virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const
+ {
+ Msg::Error("Cannot be called");
+ }
+
+ virtual void createIPState(IPVEVPUMAT *ivi, IPVEVPUMAT *iv1, IPVEVPUMAT *iv2) const;
+ virtual void createIPVariable(IPVEVPUMAT *ipv,bool hasBodyForce, const MElement *ele,const int nbFF, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}; // this law is initialized so nothing to do
+
+ virtual void callUMAT(double *stress, double *statev, double **ddsdde, double &sse, double &spd, double &scd, double &rpl,
+ double *ddsddt, double *drplde, double &drpldt, double *stran, double *dtsran,
+ double *tim, double timestep, double temperature, double deltaTemperature, double *predef, double *dpred,
+ const char *CMNAME, int &ndi, int &nshr, int tensorSize,
+ int statevSize, double *prs, int matPropSize, double *coords, double **dRot,
+ double &pnewdt, double &celent, double **F0, double **F1,
+ int &noel, int &npt, int &layer, int &kspt, int &kstep, int &kinc) const;
+
+ virtual const char* getCMNAME() const {return "VEVP";}
+ virtual double getDensity() const {return _rho;}
+
+ #endif // SWIG
+};
+
+#endif // MLAWVEVPUMAT_H_
diff --git a/NonLinearSolver/materialLaw/vevp.f b/NonLinearSolver/materialLaw/vevp.f
new file mode 100644
index 0000000000000000000000000000000000000000..78af446fa88ef55093758a7853bbdae22c85f48c
--- /dev/null
+++ b/NonLinearSolver/materialLaw/vevp.f
@@ -0,0 +1,1431 @@
+C UMAT - Finite strain VEVP
+C Based on Paper by V.-D.Nguyen et al. https://orbi.uliege.be/handle/2268/197898
+
+C Total no of material parameters : 51
+C Total no of internal variables (STATEV) : 108
+C
+C Note: The UMAT is hardcoded for 8 Maxwell VE Branches.
+C If you feel the need for more branches send me an email
+C I will add more.
+C
+C Note: Tollerance and max iteration for the newton raphson solver are hard coded
+C as variables TOLL_G and MAX_i respectively
+C
+C Note : Computations are done in standard matrix notation and
+C at the end changed to voit notation where needed for return to main
+C routine. I understand this isn't the most elegant way of
+C doing things but when dealing with tensors
+C with 4 or more dummy indices, I can code much
+C faster this way. If you know the algebra and have the
+C time to work in voit notation, change it and send me an
+C email.
+C
+C Note : Consistent tangents are calculated numerically.
+C Tollerance for numeric tangent is hardcoded as variable TOLL.
+
+C Mohib Mustafa - First Version - IMDEA 15 March 2022
+C Mohib Mustafa - Added more VE Branches - ULIEGE 10 Feburary 2023
+C mohibmustafa@gmail.com
+
+
+
+C !--------------------------------------------------------------
+C ! UMAT SUBROUTINE
+C !--------------------------------------------------------------
+
+ SUBROUTINE umat_vevp(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
+ 1 RPL,DDSDDT,DRPLDE,DRPLDT,
+ 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,
+ 3 DPRED,CMNAME,
+ 4 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,
+ 5 DROT,PNEWDT,
+ 6 CELENT,DFGRD0,DFGRD1,NOEL,NPT,
+ 7 LAYER,KSPT,KSTEP,KINC)
+
+C INCLUDE 'ABA_PARAM.INC'
+
+ CHARACTER*80 CMNAME
+ REAL*8 STRESS(NTENS),STATEV(NSTATV),
+ 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
+ 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
+ 3 PROPS(NPROPS),DFGRD0(3,3), DFGRD1(3,3)
+
+ REAL*8 SSE, SPD, SCD, RPL, DTIME, TEMP, CELENT, DRPLDT
+ REAL*8 PNEWDT, DTEMP, R, dR, ddR
+
+ call VEVP(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,PROPS,
+ & NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1,
+ & PNEWDT)
+
+ return
+ end
+
+
+
+C !--------------------------------------------------------------
+C ! Finite Strain VEVP Subroutine
+C !--------------------------------------------------------------
+
+ SUBROUTINE VEVP(STRESS,STATEV,DDSDDE,STRAN,NTENS,NSTATV,
+ 1 PROPS,NPROPS,DTIME,DSTRAN,KINC,KSTEP,NOEL,DFGRD0,DFGRD1,
+ 2 PNEWDT)
+
+
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: ntens,nprops,nstatv
+ INTEGER, INTENT(IN) :: kinc,kstep,noel
+ REAL(prec), INTENT(IN) :: stran(ntens), dstran(ntens)
+ REAL(prec), INTENT(IN) :: props(nprops), dtime
+ REAL(prec), INTENT(IN) :: DFGRD0(3,3), DFGRD1(3,3)
+ REAL(prec), INTENT(INOUT) :: statev(nstatv)
+ REAL(prec), INTENT(OUT) :: stress(ntens)
+ REAL(prec), INTENT(OUT) :: ddsdde(ntens,ntens), PNEWDT
+
+ !Declaration of local variables
+ INTEGER :: ii, jj, O6, order
+ REAL(prec) :: I_mat(3,3)
+ ! Internal Variables
+ REAL(prec) :: F_vp_n(3, 3), E_ve_n(3, 3)
+ REAL(prec) :: gma_0, gma_n, b_n(3, 3)
+ REAL(prec) :: AA_n(8, 3, 3), BB_n(8)
+ ! Material Params
+ REAL(prec) :: KK_inf, GG_inf, alpha, nu_p, eta, p_exp
+ REAL(prec) :: sigmac0, hc1, hc2, hcexp
+ REAL(prec) :: sigmat0, ht1, ht2, htexp
+ REAL(prec) :: hb0, hb1, hb2
+ REAL(prec) :: KK(8), k(8), GG(8), g(8), beta, k_plast
+ ! VE trial variables
+ REAL(prec) :: F_ve_tr(3, 3), C_ve_tr(3, 3), D_ve_tr(3, 3)
+ REAL(prec) :: E_ve_tr(3, 3), dlogC_ve_tr(3,3,3,3)
+ REAL(prec) :: ddlogC_ve_tr(3,3,3,3,3,3), AA_tr(8,3,3)
+ REAL(prec) :: BB_tr(8),GGe, KKe, kappa_tr(3, 3), S_tr(3, 3)
+ REAL(prec) :: temp1(3,3), tau_tr(3, 3)
+ ! Variables for trial yield
+ REAL(prec) :: phi_tr(3, 3), tr_phi_tr, dev_phi_tr(3, 3)
+ REAL(prec) :: phi_p_tr, phi_e_tr, F_tr
+ ! trial VE variables - perturbation
+ REAL(prec) :: F_hat(6, 3, 3), J, F_ve_tr_hat(6, 3, 3)
+ REAL(prec) :: C_ve_tr_hat(6, 3, 3), D_ve_tr_hat(6,3,3)
+ REAL(prec) :: E_ve_tr_hat(6,3,3), dlogC_ve_tr_hat(6,3,3,3,3)
+ REAL(prec) :: ddlogC_ve_tr_hat(6,3,3,3,3,3,3)
+ REAL(prec) :: AA_tr_hat(6,8,3,3)
+ REAL(prec) :: BB_tr_hat(6, 8), GGe_hat(6), KKe_hat(6)
+ REAL(prec) :: kappa_tr_hat(6, 3, 3), S_tr_hat(6, 3, 3)
+ REAL(prec) :: temp1_hat(6, 3, 3), tau_tr_hat(6, 3, 3)
+ REAL(prec) :: tau_tr_hat_v(6, 6)
+ ! VP variables
+ REAL(prec) :: ptilde, PhiEq, sigma_c, sigma_t, HHc, HHt, HHb
+ REAL(prec) :: m, a0, a1, a2, GAMMA, u, v, dev_phi(3, 3)
+ REAL(prec) :: dev_Q(3, 3), tr_Q, Q(3, 3), GQ(3,3)
+ REAL(prec) :: exp_GQ(3, 3), F_vp(3,3), F_vp_inv(3,3)
+ ! VE corrected variables
+ REAL(prec) :: F_ve(3, 3), C_ve(3, 3), D_ve(3, 3), E_ve(3, 3)
+ REAL(prec) :: dlogC_ve(3,3,3,3), ddlogC_ve(3,3,3,3,3,3)
+ REAL(prec) :: AA(8,3,3), BB(8), kappa(3, 3), b(3, 3)
+ REAL(prec) :: S(3, 3), temp2(3, 3), tau(3, 3)
+ ! VP variables perturbated
+ REAL(prec) :: phi_tr_hat(6,3,3), dev_phi_tr_hat(6, 3, 3)
+ REAL(prec) :: tr_phi_tr_hat(6), phi_p_tr_hat(6)
+ REAL(prec) :: phi_e_tr_hat(6), ptilde_hat(6), PhiEq_hat(6)
+ REAL(prec) :: gma_n_hat(6), sigma_c_hat(6), sigma_t_hat(6)
+ REAL(prec) :: HHc_hat(6), HHt_hat(6), HHb_hat(6), m_hat(6)
+ REAL(prec) :: a0_hat(6), a1_hat(6), a2_hat(6), F_tr_hat(6)
+ REAL(prec) :: GAMMA_hat(6), u_hat(6), v_hat(6)
+ REAL(prec) :: dev_phi_hat(6, 3, 3), dev_Q_hat(6, 3, 3)
+ REAL(prec) :: tr_Q_hat(6), Q_hat(6,3,3)
+ REAL(prec) :: GQ_hat(6,3,3), exp_GQ_hat(6,3,3)
+ REAL(prec) :: F_vp_hat(6,3,3), F_vp_inv_hat(6,3,3)
+ ! perturbated VE corrected variables
+ REAL(prec) :: F_ve_hat(6,3,3), C_ve_hat(6,3,3)
+ REAL(prec) :: D_ve_hat(6,3,3), E_ve_hat(6,3,3)
+ REAL(prec) :: dlogC_ve_hat(6,3,3,3,3)
+ REAL(prec) :: ddlogC_ve_hat(6,3,3,3,3,3,3), AA_hat(6,8,3,3)
+ REAL(prec) :: BB_hat(6,8), kappa_hat(6,3,3), S_hat(6,3,3)
+ REAL(prec) :: temp2_hat(6,3,3), tau_hat(6,3,3)
+ REAL(prec) :: tau_hat_v(6,6)
+
+ ! Tollerances
+ REAL(prec), PARAMETER :: TOLL=0.001D0, TOLL_G=0.999D-11
+ INTEGER, PARAMETER :: MAX_i=500
+
+ !Define 2nd order identity tensor
+ data I_mat(1,:) /1.D0, 0.D0, 0.D0/
+ data I_mat(2,:) /0.D0, 1.D0, 0.D0/
+ data I_mat(3,:) /0.D0, 0.D0, 1.D0/
+
+ ! This is a workaround to set the diagnol values of F_vp = 1
+ ! for the first time step (initial condition) in FFTMAD.
+ ! Can be commented out for ABAQUS, If you comment out, then
+ ! in abaqus use the inp file to set initial values for F_vp.
+C IF ((KINC .EQ. 1) .AND. (KSTEP .EQ. 1)) THEN
+C STATEV(1) = 1.D0
+C STATEV(2) = 1.D0
+C STATEV(3) = 1.D0
+C END IF
+
+ ! State variables at previous time step
+ CALL voit2mat(STATEV(1:9), F_vp_n(:,:))
+ CALL voit2mat(STATEV(10:18), E_ve_n(:,:))
+ gma_n = STATEV(19)
+ gma_0 = gma_n
+ CALL voit2mat(STATEV(20:28), b_n(:,:))
+ CALL voit2mat(STATEV(29:37), AA_n(1,:,:))
+ CALL voit2mat(STATEV(38:46), AA_n(2,:,:))
+ CALL voit2mat(STATEV(47:55), AA_n(3,:,:))
+ CALL voit2mat(STATEV(56:64), AA_n(4,:,:))
+ CALL voit2mat(STATEV(65:73), AA_n(5,:,:))
+ CALL voit2mat(STATEV(74:82), AA_n(6,:,:))
+ CALL voit2mat(STATEV(83:91), AA_n(7,:,:))
+ CALL voit2mat(STATEV(92:100), AA_n(8,:,:))
+
+ BB_n(1) = STATEV(101)
+ BB_n(2) = STATEV(102)
+ BB_n(3) = STATEV(103)
+ BB_n(4) = STATEV(104)
+ BB_n(5) = STATEV(105)
+ BB_n(6) = STATEV(106)
+ BB_n(7) = STATEV(107)
+ BB_n(8) = STATEV(108)
+
+ ! Get Material Properties
+ order =PROPS(1)
+ KK_inf =PROPS(2)
+ GG_inf =PROPS(3)
+ alpha =PROPS(4)
+ nu_p =PROPS(5)
+ eta =PROPS(6)
+ p_exp =PROPS(7)
+ sigmac0 =PROPS(8)
+
+ hc1 =PROPS(9)
+ hc2 =PROPS(10)
+ hcexp =PROPS(11)
+ sigmat0 =PROPS(12)
+ ht1 =PROPS(13)
+ ht2 =PROPS(14)
+ htexp =PROPS(15)
+ hb0 =PROPS(16)
+
+ hb1 =PROPS(17)
+ hb2 =PROPS(18)
+ KK(1) =PROPS(19)
+ KK(2) =PROPS(20)
+ KK(3) =PROPS(21)
+ KK(4) =PROPS(22)
+ KK(5) =PROPS(23)
+ KK(6) =PROPS(24)
+ KK(7) =PROPS(25)
+ KK(8) =PROPS(26)
+ k(1) =PROPS(27)
+ k(2) =PROPS(28)
+ k(3) =PROPS(29)
+ k(4) =PROPS(30)
+ k(5) =PROPS(31)
+ k(6) =PROPS(32)
+ k(7) =PROPS(33)
+ k(8) =PROPS(34)
+
+ GG(1) =PROPS(35)
+ GG(2) =PROPS(36)
+ GG(3) =PROPS(37)
+ GG(4) =PROPS(38)
+ GG(5) =PROPS(39)
+ GG(6) =PROPS(40)
+ GG(7) =PROPS(41)
+ GG(8) =PROPS(42)
+ g(1) =PROPS(43)
+ g(2) =PROPS(44)
+ g(3) =PROPS(45)
+ g(4) =PROPS(46)
+ g(5) =PROPS(47)
+ g(6) =PROPS(48)
+ g(7) =PROPS(49)
+ g(8) =PROPS(50)
+
+ ! Compute determinant of F (Deformation Gradient)
+ CALL determinant(DFGRD1(:,:), J)
+ ! Calculate VE Right Cauchy Green Tensor (C_ve) at trial State
+ CALL vevpSplit(DFGRD1(:, :), F_vp_n(:, :),
+ 1 F_ve_tr(:, :), C_ve_tr(:, :))
+
+ ! Approximate VE log Strain (power series) at trial state
+ D_ve_tr(:,:) = C_ve_tr(:, :) - I_mat(:,:)
+ CALL approx_log(D_ve_tr(:,:), order, E_ve_tr(:, :),
+ 1 dlogC_ve_tr(:,:,:,:), ddlogC_ve_tr(:,:,:,:,:,:))
+ E_ve_tr(:, :) = 0.5D0 * E_ve_tr(:, :)
+
+ ! Calculate the corotated kirchoff stress at trial state along
+ ! with VE internal variables
+ CALL vePredictor_log(E_ve_tr(:,:), E_ve_n(:,:), AA_n(:,:,:),
+ 1 BB_n(:), DTIME, PROPS, NPROPS, AA_tr(:,:,:),
+ 2 BB_tr(:), GGe, KKe, kappa_tr(:,:))
+
+ ! Check yielding at trial state => gma = gma_n, u=1, v=1, GAMMA = 0
+ ! Get phi_e_tr, phi_p_tr
+ phi_tr(:, :) = kappa_tr(:,:) - b_n(:,:)
+
+ CALL tr_dev_split(phi_tr(:,:), dev_phi_tr(:,:), tr_phi_tr)
+ phi_p_tr = tr_phi_tr / 3.D0
+
+ CALL mat2ddot(dev_phi_tr(:, :), dev_phi_tr(:, :), phi_e_tr)
+ phi_e_tr = (3.D0 / 2.D0) * phi_e_tr
+ phi_e_tr = SQRT(phi_e_tr)
+
+ ptilde = phi_p_tr
+ PhiEq = phi_e_tr
+
+ ! Update hardening variables on trial state
+ CALL hardn(PROPS, NPROPS, gma_0, gma_n, sigma_c, sigma_t,
+ 1 HHc, HHt, HHb)
+ ! Update Drucker Pragger coefficients along with ecc factor m
+ CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
+
+ F_tr = a2 * PhiEq**alpha - a1 * ptilde - a0
+
+ ! Numeric tangent
+ CALL perturb_F(DFGRD1(:, :), TOLL, F_hat(:, :, :))
+ DO O6 = 1, 6
+
+ CALL vevpSplit(F_hat(O6, :, :), F_vp_n(:, :),
+ 1 F_ve_tr_hat(O6, :, :), C_ve_tr_hat(O6, :, :))
+
+ ! Approximate VE log Strain at trial state
+ D_ve_tr_hat(O6, :,:) = C_ve_tr_hat(O6, :, :) - I_mat(:,:)
+ CALL approx_log(D_ve_tr_hat(O6, :,:), order,
+ 1 E_ve_tr_hat(O6, :, :), dlogC_ve_tr_hat(O6,:,:,:,:),
+ 2 ddlogC_ve_tr_hat(O6,:,:,:,:,:,:))
+ E_ve_tr_hat(O6,:, :) = 0.5D0 * E_ve_tr_hat(O6, :, :)
+
+ ! Calculate the corotated kirchoff stress at trial state
+ ! along with VE internal variables
+ CALL vePredictor_log(E_ve_tr_hat(O6, :,:), E_ve_n(:,:),
+ 1 AA_n(:,:,:), BB_n(:), DTIME, PROPS, NPROPS,
+ 2 AA_tr_hat(O6, :,:,:),
+ 3 BB_tr_hat(O6, :), GGe_hat(O6), KKe_hat(O6),
+ 4 kappa_tr_hat(O6, :,:))
+
+ ! Calculate 2Pk as per paper
+ CALL mat24ddot(kappa_tr_hat(O6, :, :),
+ 1 dlogC_ve_tr_hat(O6, :, :, :, :), S_tr_hat(O6, :, :))
+
+ ! Calculate kirchoff stress
+ CALL mat2dot(F_ve_tr_hat(O6, :, :), S_tr_hat(O6, :, :),
+ 1 temp1_hat(O6, :, :))
+ CALL mat2dotT(temp1_hat(O6, :, :), F_ve_tr_hat(O6, :, :),
+ 1 tau_tr_hat(O6, :, :))
+
+ END DO
+
+ ! In case the yield surface is not breached for the trial state (VE step)
+ IF(F_tr .LE. TOLL_G) THEN
+
+ ! Calculate 2Pk as per paper
+ CALL mat24ddot(kappa_tr(:, :), dlogC_ve_tr(:, :, :, :),
+ 1 S_tr(:, :))
+
+ ! Calculate kirchoff stress
+ CALL mat2dot(F_ve_tr(:, :), S_tr(:, :), temp1(:, :))
+ CALL mat2dotT(temp1(:, :), F_ve_tr(:, :), tau_tr(:, :))
+
+ ! Return cauchy stress for abaqus
+ STRESS(1) = (1.D0 / J) * tau_tr(1, 1)
+ STRESS(2) = (1.D0 / J) * tau_tr(2, 2)
+ STRESS(3) = (1.D0 / J) * tau_tr(3, 3)
+ STRESS(4) = (1.D0 / J) * tau_tr(1, 2)
+ STRESS(5) = (1.D0 / J) * tau_tr(1, 3)
+ STRESS(6) = (1.D0 / J) * tau_tr(2, 3)
+
+ ! Update internal variables for trial state. No need to update VP internal variables
+ CALL mat2voit(E_ve_tr(:,:),STATEV(10:18))
+
+ CALL mat2voit(AA_tr(1,:,:),STATEV(29:37))
+ CALL mat2voit(AA_tr(2,:,:),STATEV(38:46))
+ CALL mat2voit(AA_tr(3,:,:),STATEV(47:55))
+ CALL mat2voit(AA_tr(4,:,:),STATEV(56:64))
+ CALL mat2voit(AA_tr(5,:,:),STATEV(65:73))
+ CALL mat2voit(AA_tr(6,:,:),STATEV(74:82))
+ CALL mat2voit(AA_tr(7,:,:),STATEV(83:91))
+ CALL mat2voit(AA_tr(8,:,:),STATEV(92:100))
+
+ STATEV(101)=BB_tr(1)
+ STATEV(102)=BB_tr(2)
+ STATEV(103)=BB_tr(3)
+ STATEV(104)=BB_tr(4)
+ STATEV(105)=BB_tr(5)
+ STATEV(106)=BB_tr(6)
+ STATEV(107)=BB_tr(7)
+ STATEV(108)=BB_tr(8)
+
+ !Turn perturbated stress into voit notation
+ tau_tr_hat_v(:, :) = 0.D0
+ DO O6 = 1, 6
+ DO ii = 1, 3
+ tau_tr_hat_v(O6, ii) = tau_tr_hat(O6, ii, ii)
+ END DO
+ tau_tr_hat_v(O6, 4) = tau_tr_hat(O6, 1, 2)
+ tau_tr_hat_v(O6, 5) = tau_tr_hat(O6, 1, 3)
+ tau_tr_hat_v(O6, 6) = tau_tr_hat(O6, 2, 3)
+ END DO
+
+ !Compute Tangent for Abaqus
+ DO ii = 1, 6
+ DO jj = 1, 6
+ DDSDDE(ii, jj) = (1.D0 / (J * TOLL))
+ 1 * (tau_tr_hat_v(jj, ii) - J*STRESS(ii))
+ END DO
+ END DO
+
+
+ ! In case the yield surface is breached at trail state, corrector steps are needed
+ ELSE
+ ! Use Newton Raphson scheme to get GAMMA, gma and related variables
+ CALL nlinSolver(PROPS, NPROPS, DTIME, TOLL_G, MAX_i,
+ 1 GGe, KKe,
+ 2 F_tr, gma_0, gma_n, ptilde, PhiEq,
+ 3 GAMMA, u, v, beta, k_plast)
+
+ ! Update dev_phi (updated ptilde already returned by the subroutine)
+ dev_phi(:,:) = dev_phi_tr(:,:) / u
+
+ ! Update flow normal Q
+ dev_Q(:, :) = 3.D0 * dev_phi(:,:)
+ tr_Q = 2.D0 * beta * ptilde / 3.D0
+ Q(:,:) = dev_Q(:,:) + tr_Q * I_mat(:,:)
+
+ ! Make corrections to F_vp, F_ve and C_ve
+ GQ(:,:) = GAMMA * Q(:,:)
+ CALL approx_exp(GQ(:,:), order, exp_GQ(:,:))
+ CALL mat2dot(exp_GQ(:,:), F_vp_n(:,:), F_vp(:,:))
+ CALL matInv(F_vp(:,:), F_vp_inv(:,:))
+ CALL mat2dot(DFGRD1(:,:), F_vp_inv(:,:), F_ve(:, :))
+ CALL mat2Tdot(F_ve(:,:), F_ve(:,:), C_ve(:,:))
+
+ ! Approximate VE log Strain (power series) at the corrected state
+ D_ve(:,:) = C_ve(:, :) - I_mat(:,:)
+ CALL approx_log(D_ve(:,:), order, E_ve(:, :),
+ 1 dlogC_ve(:,:,:,:), ddlogC_ve(:,:,:,:,:,:))
+ E_ve(:, :) = 0.5D0 * E_ve(:, :)
+
+ ! Calculate the corotated kirchoff stress at corrected state along
+ ! with VE internal variables. Note : I know i am calculating the stresses again, whereas
+ ! I just need to correct the VE internal variables, it doesnt affect performance alot.
+ CALL vePredictor_log(E_ve(:,:), E_ve_n(:,:), AA_n(:,:,:),
+ 1 BB_n(:), DTIME, PROPS, NPROPS, AA(:,:,:), BB(:),
+ 2 GGe, KKe, kappa(:,:))
+
+ ! Calculate 2PK stress at corrected state
+ CALL mat24ddot(kappa(:, :), dlogC_ve(:, :, :, :), S(:, :))
+ ! Calculate kirchoff stress at corrected state
+ CALL mat2dot(F_ve(:, :), S(:, :), temp2(:, :))
+ CALL mat2dotT(temp2(:, :), F_ve(:, :), tau(:, :))
+
+ ! Return cauchy stress for abaqus
+ STRESS(1) = (1.D0 / J) * tau(1, 1)
+ STRESS(2) = (1.D0 / J) * tau(2, 2)
+ STRESS(3) = (1.D0 / J) * tau(3, 3)
+ STRESS(4) = (1.D0 / J) * tau(1, 2)
+ STRESS(5) = (1.D0 / J) * tau(1, 3)
+ STRESS(6) = (1.D0 / J) * tau(2, 3)
+
+ ! Update kin hardening variable for the corrected state
+ CALL hardn(PROPS, NPROPS, gma_0, gma_n,
+ 1 sigma_c, sigma_t, HHc, HHt, HHb)
+
+ b(:,:) = b_n(:,:) + k_plast*HHb * GQ(:,:)
+
+ ! Return updated internal variables
+ CALL mat2voit(F_vp(:,:),STATEV(1:9))
+ CALL mat2voit(E_ve(:,:),STATEV(10:18))
+ STATEV(19) = gma_n
+ CALL mat2voit(b(:,:),STATEV(20:28))
+
+ CALL mat2voit(AA(1,:,:),STATEV(29:37))
+ CALL mat2voit(AA(2,:,:),STATEV(38:46))
+ CALL mat2voit(AA(3,:,:),STATEV(47:55))
+ CALL mat2voit(AA(4,:,:),STATEV(56:64))
+ CALL mat2voit(AA(5,:,:),STATEV(65:73))
+ CALL mat2voit(AA(6,:,:),STATEV(74:82))
+ CALL mat2voit(AA(7,:,:),STATEV(83:91))
+ CALL mat2voit(AA(8,:,:),STATEV(92:100))
+
+ STATEV(101)=BB(1)
+ STATEV(102)=BB(2)
+ STATEV(103)=BB(3)
+ STATEV(104)=BB(4)
+ STATEV(105)=BB(5)
+ STATEV(106)=BB(6)
+ STATEV(107)=BB(7)
+ STATEV(108)=BB(8)
+
+ ! Numeric Tangent
+ DO O6 = 1, 6
+
+ gma_n_hat(O6) = gma_n
+
+ ! Get phi_e_tr, phi_p_tr
+ phi_tr_hat(O6, :, :) = kappa_tr_hat(O6, :,:) - b_n(:,:)
+
+ CALL tr_dev_split(phi_tr_hat(O6, :,:),
+ 1 dev_phi_tr_hat(O6, :,:), tr_phi_tr_hat(O6))
+ phi_p_tr_hat(O6) = tr_phi_tr_hat(O6) / 3.D0
+
+ CALL mat2ddot(dev_phi_tr_hat(O6, :, :),
+ 1 dev_phi_tr_hat(O6, :, :), phi_e_tr_hat(O6))
+ phi_e_tr_hat(O6) = (3.D0 / 2.D0) * phi_e_tr_hat(O6)
+ phi_e_tr_hat(O6) = SQRT(phi_e_tr_hat(O6))
+
+ ptilde_hat(O6) = phi_p_tr_hat(O6)
+ PhiEq_hat(O6) = phi_e_tr_hat(O6)
+
+ ! Update hardening variables on purturb state
+ CALL hardn(PROPS, NPROPS, gma_0, gma_n_hat(O6),
+ 1 sigma_c_hat(O6), sigma_t_hat(O6), HHc_hat(O6),
+ 2 HHt_hat(O6), HHb_hat(O6))
+ ! Update Drucker Pragger coefficients along with ecc factor m
+ CALL DPcoeff(alpha, sigma_c_hat(O6), sigma_t_hat(O6),
+ 1 m_hat(O6), a0_hat(O6), a1_hat(O6), a2_hat(O6))
+ ! Calculate yield func at purturbed state
+ F_tr_hat(O6) = a2_hat(O6) * PhiEq_hat(O6)**alpha
+ 1 - a1_hat(O6) * ptilde_hat(O6) - a0_hat(O6)
+
+ ! Solve Newton raphson to get gma and GAMMA at purturbed state
+ CALL nlinSolver(PROPS, NPROPS, DTIME, TOLL_G, MAX_i,
+ 1 GGe, KKe,
+ 2 F_tr_hat(O6), gma_0, gma_n_hat(O6), ptilde_hat(O6),
+ 3 PhiEq_hat(O6), GAMMA_hat(O6), u_hat(O6), v_hat(O6),
+ 4 beta, k_plast)
+
+ ! Update dev_phi_hat (ptilde_hat already returned by the subroutine)
+ dev_phi_hat(O6, :,:) = dev_phi_tr_hat(O6,:,:) / u_hat(O6)
+
+ ! Update flow normal Q_hat
+ dev_Q_hat(O6, :, :) = 3.D0 * dev_phi_hat(O6, :,:)
+ tr_Q_hat(O6) = 2.D0 * beta * ptilde_hat(O6) / 3.D0
+ Q_hat(O6,:,:) = dev_Q_hat(O6,:,:)
+ 1 + tr_Q_hat(O6) * I_mat(:,:)
+
+ ! Make corrections to F_vp_hat, F_ve_hat and C_ve_hat
+ GQ_hat(O6,:,:) = GAMMA_hat(O6) * Q_hat(O6,:,:)
+ CALL approx_exp(GQ_hat(O6,:,:), order, exp_GQ_hat(O6,:,:))
+ CALL mat2dot(exp_GQ_hat(O6, :,:), F_vp_n(:,:),
+ 1 F_vp_hat(O6,:,:))
+ CALL matInv(F_vp_hat(O6,:,:), F_vp_inv_hat(O6,:,:))
+ CALL mat2dot(F_hat(O6, :,:), F_vp_inv_hat(O6,:,:),
+ 1 F_ve_hat(O6, :, :))
+ CALL mat2Tdot(F_ve_hat(O6,:,:), F_ve_hat(O6,:,:),
+ 1 C_ve_hat(O6,:,:))
+
+
+ ! Approximate VE log Strain (power series) at the purturbated state
+ D_ve_hat(O6,:,:) = C_ve_hat(O6,:, :) - I_mat(:,:)
+ CALL approx_log(D_ve_hat(O6,:,:), order, E_ve_hat(O6,:, :),
+ 1 dlogC_ve_hat(O6,:,:,:,:),
+ 2 ddlogC_ve_hat(O6,:,:,:,:,:,:))
+ E_ve_hat(O6,:, :) = 0.5D0 * E_ve_hat(O6,:, :)
+
+ ! Calculate the corotated kirchoff stress at purturbed state along
+ ! with VE internal variables
+ CALL vePredictor_log(E_ve_hat(O6,:,:), E_ve_n(:,:),
+ 1 AA_n(:,:,:), BB_n(:), DTIME, PROPS, NPROPS,
+ 2 AA_hat(O6,:,:,:),
+ 3 BB_hat(O6,:), GGe, KKe, kappa_hat(O6,:,:))
+
+ ! Calculate 2PK stress at purturbed state
+ CALL mat24ddot(kappa_hat(O6,:, :),
+ 1 dlogC_ve_hat(O6,:, :, :, :), S_hat(O6,:, :))
+
+ ! Calculate kirchoff stress at purturbed state
+ CALL mat2dot(F_ve_hat(O6,:, :), S_hat(O6,:, :),
+ 1 temp2_hat(O6,:, :))
+ CALL mat2dotT(temp2_hat(O6,:, :), F_ve_hat(O6,:, :),
+ 1 tau_hat(O6,:, :))
+
+ END DO
+
+ !Turn perturbated stress into voit notation
+ tau_hat_v(:, :) = 0.D0
+ DO O6 = 1, 6
+ DO ii = 1, 3
+ tau_hat_v(O6, ii) = tau_hat(O6, ii, ii)
+ END DO
+ tau_hat_v(O6, 4) = tau_hat(O6, 1, 2)
+ tau_hat_v(O6, 5) = tau_hat(O6, 1, 3)
+ tau_hat_v(O6, 6) = tau_hat(O6, 2, 3)
+ END DO
+
+ !Tangent for Abaqus
+ DO ii = 1, 6
+ DO jj = 1, 6
+ DDSDDE(ii, jj) = (1.D0 / (J * TOLL))
+ 1 * (tau_hat_v(jj, ii) - J*STRESS(ii))
+ END DO
+ END DO
+
+ END IF
+
+ RETURN
+ END SUBROUTINE VEVP
+
+ !--------------------------------------------------------------
+ ! Helper function to convert a voit array to matrix
+ !--------------------------------------------------------------
+
+ SUBROUTINE voit2mat(voit, mat)
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: voit(9)
+ REAL(prec), INTENT(OUT):: mat(3, 3)
+
+ mat(1, 1) = voit(1)
+ mat(2, 2) = voit(2)
+ mat(3, 3) = voit(3)
+ mat(1, 2) = voit(4)
+ mat(1, 3) = voit(5)
+ mat(2, 3) = voit(6)
+ mat(2, 1) = voit(7)
+ mat(3, 1) = voit(8)
+ mat(3, 2) = voit(9)
+ RETURN
+ END SUBROUTINE
+
+ !--------------------------------------------------------------
+ ! Helper function to convert a matrix to voit array
+ !--------------------------------------------------------------
+
+ SUBROUTINE mat2voit(mat, voit)
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: mat(3, 3)
+ REAL(prec), INTENT(OUT):: voit(9)
+
+ voit(1) = mat(1, 1)
+ voit(2) = mat(2, 2)
+ voit(3) = mat(3, 3)
+ voit(4) = mat(1, 2)
+ voit(5) = mat(1, 3)
+ voit(6) = mat(2, 3)
+ voit(7) = mat(2, 1)
+ voit(8) = mat(3, 1)
+ voit(9) = mat(3, 2)
+ RETURN
+ END SUBROUTINE
+
+ !--------------------------------------------------------------
+ ! Helper function to compute Determinant of 3x3 matrix
+ !--------------------------------------------------------------
+ SUBROUTINE determinant(matrix, det)
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: matrix(3, 3)
+ REAL(prec), INTENT(OUT):: det
+
+ det = matrix(1,1) * (matrix(2,2) * matrix(3,3)
+ 1 - matrix(3,2) * matrix(2, 3)) - matrix(1,2) * (matrix(2,1)
+ 2 * matrix(3,3) - matrix(2,3) * matrix(3,1)) + matrix(1,3)
+ 3 * (matrix(2,1) * matrix(3,2) - matrix(2,2) * matrix(3,1))
+ RETURN
+ END SUBROUTINE determinant
+
+
+ !--------------------------------------------------------------
+ ! Helper function to compute inverse of 3x3 matrix
+ !--------------------------------------------------------------
+ SUBROUTINE matInv(matrix, inv)
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: matrix(3, 3)
+ REAL(prec), INTENT(OUT) :: inv(3, 3)
+ REAL(prec) :: J
+
+ inv(1,1) = matrix(2,2)*matrix(3,3) - matrix(3,2)*matrix(2,3)
+ inv(1,2) = matrix(3,2)*matrix(1,3) - matrix(1,2)*matrix(3,3)
+ inv(1,3) = matrix(1,2)*matrix(2,3) - matrix(2,2)*matrix(1,3)
+ inv(2,1) = matrix(3,1)*matrix(2,3) - matrix(2,1)*matrix(3,3)
+ inv(2,2) = matrix(1,1)*matrix(3,3) - matrix(3,1)*matrix(1,3)
+ inv(2,3) = matrix(2,1)*matrix(1,3) - matrix(1,1)*matrix(2,3)
+ inv(3,1) = matrix(2,1)*matrix(3,2) - matrix(3,1)*matrix(2,2)
+ inv(3,2) = matrix(3,1)*matrix(1,2) - matrix(1,1)*matrix(3,2)
+ inv(3,3) = matrix(1,1)*matrix(2,2) - matrix(2,1)*matrix(1,2)
+
+ CALL determinant(matrix, J)
+
+ inv(:, :) = inv(:, :)/ J
+
+ RETURN
+ END SUBROUTINE matInv
+
+
+ !------------------------------------------------------------
+ ! Helper function for c = a.b
+ ! (2nd order tensors, single contraction)
+ !------------------------------------------------------------
+
+ SUBROUTINE mat2dot(a, b, c)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: a(3, 3), b(3, 3)
+ REAL(prec), INTENT(OUT) :: c(3, 3)
+
+ INTEGER :: O1, ii, jj
+
+ c(:, :) = 0.D0
+
+ DO ii = 1, 3
+ DO jj = 1, 3
+ DO O1 = 1, 3
+ c(ii, jj) = c(ii, jj)
+ 1 + a(ii, O1) * b(O1, jj)
+ END DO
+ END DO
+ END DO
+
+ END SUBROUTINE mat2dot
+
+ !------------------------------------------------------------
+ ! Helper function for c = a^T.b
+ ! (2nd order tensors, single contraction)
+ !------------------------------------------------------------
+ SUBROUTINE mat2Tdot(a, b, c)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: a(3, 3), b(3, 3)
+ REAL(prec), INTENT(OUT) :: c(3, 3)
+
+ INTEGER :: O1, ii, jj
+
+ c(:, :) = 0.D0
+
+ DO O1 = 1, 3
+ DO ii = 1, 3
+ DO jj = 1, 3
+ c(ii, jj) = c(ii, jj)
+ 1 + a(O1, ii) * b(O1, jj)
+ END DO
+ END DO
+ END DO
+
+ END SUBROUTINE mat2Tdot
+
+ !------------------------------------------------------------
+ ! Helper function for c = a.b^T
+ ! (2nd order tensors, single contraction)
+ !------------------------------------------------------------
+ SUBROUTINE mat2dotT(a, b, c)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: a(3, 3), b(3, 3)
+ REAL(prec), INTENT(OUT) :: c(3, 3)
+
+ INTEGER :: O1, ii, jj
+
+ c(:, :) = 0.D0
+
+ DO O1 = 1, 3
+ DO ii = 1, 3
+ DO jj = 1, 3
+ c(ii, jj) = c(ii, jj)
+ 1 + a(ii, O1) * b(jj, O1)
+ END DO
+ END DO
+ END DO
+
+ END SUBROUTINE mat2dotT
+
+ !------------------------------------------------------------
+ ! Helper function for c = a:b
+ ! (2nd order tensors, double contraction)
+ !------------------------------------------------------------
+ SUBROUTINE mat2ddot(a, b, c)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: a(3, 3), b(3, 3)
+ REAL(prec), INTENT(OUT) :: c
+
+ INTEGER :: O1, ii, jj
+
+ c = 0.D0
+
+ DO ii = 1, 3
+ DO jj = 1, 3
+ c = c + a(ii, jj) * b(ii, jj)
+ END DO
+ END DO
+
+ END SUBROUTINE mat2ddot
+
+ !--------------------------------------------------------------
+ ! Helper function for 2nd order 4th order contraction
+ !--------------------------------------------------------------
+ SUBROUTINE mat24ddot(a, b, c)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: a(3, 3), b(3, 3, 3, 3)
+ REAL(prec), INTENT(OUT) :: c(3, 3)
+
+ INTEGER :: ii, jj, mm, nn
+
+ DO ii = 1, 3
+ DO jj= 1, 3
+ c(ii, jj) = 0.D0
+ DO mm = 1, 3
+ DO nn = 1, 3
+ c(ii, jj) = c(ii, jj)
+ 1 + a(mm, nn) * b(mm, nn, ii, jj)
+ END DO
+ END DO
+ END DO
+ END DO
+
+ END SUBROUTINE mat24ddot
+
+ !--------------------------------------------------------------
+ ! Helper function to get additive dev trace split
+ !--------------------------------------------------------------
+ SUBROUTINE tr_dev_split(matrix, dev, tr)
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: matrix(3, 3)
+ REAL(prec), INTENT(OUT) :: dev(3, 3), tr
+ REAL(prec) :: I_mat(3, 3)
+
+ !Define 2nd order identity tensor
+ data I_mat(1,:) /1.D0, 0.D0, 0.D0/
+ data I_mat(2,:) /0.D0, 1.D0, 0.D0/
+ data I_mat(3,:) /0.D0, 0.D0, 1.D0/
+
+ tr = matrix(1, 1) + matrix(2, 2) + matrix(3, 3)
+
+ dev(:, :) = matrix(:, :) - (1.D0 / 3.D0) * tr * I_mat(:, :)
+ RETURN
+ END SUBROUTINE tr_dev_split
+
+ !--------------------------------------------------------------
+ ! Helper function to get VEVP Multiplicative Split
+ !--------------------------------------------------------------
+ SUBROUTINE vevpSplit(F, F_vp, F_ve, C_ve)
+
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: F(3, 3), F_vp(3, 3)
+ REAL(prec), INTENT(OUT) :: F_ve(3, 3), C_ve(3, 3)
+ REAL(prec) :: F_vp_inv(3, 3)
+
+
+ ! Compute inverse of int variable F_vp
+ CALL matInv(F_vp(:, :), F_vp_inv(:, :))
+
+ ! Compute F_ve
+ CALL mat2dot(F(:, :), F_vp_inv(:, :), F_ve(:, :))
+
+ ! Compute C_ve
+ CALL mat2Tdot(F_ve(:, :), F_ve(:, :), C_ve(:, :))
+ END SUBROUTINE vevpSplit
+
+ !--------------------------------------------------------------
+ ! Function to compute perturb of matrix
+ !--------------------------------------------------------------
+
+ SUBROUTINE perturb_F(F, TOLL, F_hat)
+ IMPLICIT NONE
+
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: F(3, 3), TOLL
+ REAL(prec), INTENT(OUT) :: F_hat(6, 3, 3)
+ REAL(prec) :: DD1(3, 3, 3, 3), DD2(3, 3, 3, 3), I_mat(3, 3)
+ REAL(prec) :: delF(3, 3, 3, 3), F_hat_mat(3, 3, 3, 3)
+ INTEGER :: ii, jj, ll, mm, nn
+
+ !Decleration of constants
+ REAL(prec), PARAMETER :: ZERO=0.D0, ONE=1.D0, TWO=2.D0
+ REAL(prec), PARAMETER :: THREE=3.D0, FOUR= 4.D0, SIX=6.D0
+ REAL(prec), PARAMETER :: NINE=9.D0
+
+ !Define 4th order identity tensor
+ data I_mat(1,:) /ONE, ZERO, ZERO/
+ data I_mat(2,:) /ZERO, ONE, ZERO/
+ data I_mat(3,:) /ZERO, ZERO, ONE/
+
+ ! Calculate Perturbation of F
+ DD1(:, :, :, :) = ZERO
+ DD2(:, :, :, :) = ZERO
+ DO ii = 1, 3
+ DO jj = 1, 3
+ DO ll = 1, 3
+ DO mm = 1, 3
+ DO nn = 1, 3
+ DD1(ii, jj, ll, mm) = DD1(ii, jj, ll, mm)
+ 1 + (TOLL/TWO) * (I_mat(ll, ii) * I_mat(nn, jj)
+ 2 * F(nn, mm))
+ DD2(ii, jj, ll, mm) = DD2(ii, jj, ll, mm)
+ 1 + (TOLL/TWO) * (I_mat(ll, jj) * I_mat(nn, ii)
+ 2 * F(nn, mm))
+ END DO
+ END DO
+ END DO
+ END DO
+ END DO
+
+ delF(:, :, :, :) = ZERO
+ DO ii = 1, 3
+ DO jj = 1, 3
+ delF(ii, jj, :, :) = DD1(ii, jj, :, :)
+ 1 + DD2(ii, jj, :, :)
+ END DO
+ END DO
+
+ F_hat_mat(:, :, :, :) = ZERO
+ DO ii = 1, 3
+ DO jj = 1, 3
+ F_hat_mat(ii, jj, :, :) = F(:, :)+ delF(ii, jj, :, :)
+ END DO
+ END DO
+
+ F_hat(1, : , :) = F_hat_mat(1, 1, :, :)
+ F_hat(2, : , :) = F_hat_mat(2, 2, :, :)
+ F_hat(3, : , :) = F_hat_mat(3, 3, :, :)
+ F_hat(4, : , :) = F_hat_mat(1, 2, :, :)
+ F_hat(5, : , :) = F_hat_mat(1, 3, :, :)
+ F_hat(6, : , :) = F_hat_mat(2, 3, :, :)
+
+ RETURN
+
+ END SUBROUTINE perturb_F
+
+ !--------------------------------------------------------------
+ ! Function to compute approx log strains and derivatives
+ ! Power series approximation
+ !--------------------------------------------------------------
+
+ SUBROUTINE approx_log(AA, order, logAA, dlogAA, ddlogAA)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: order
+ REAL(prec), INTENT(IN) :: AA(3, 3)
+ REAL(prec), INTENT(OUT) :: logAA(3, 3), dlogAA(3, 3, 3, 3)
+ REAL(prec), INTENT(OUT) :: ddlogAA(3, 3, 3, 3, 3, 3)
+
+ ! List of internal variables
+ INTEGER :: ii, jj, nn, mm, ll, rr, ss, pp, qq
+ REAL(prec) :: coeffs(order + 1), I_mat(3, 3)
+ REAL(prec) :: FF(3, 3), dFF(3, 3, 3, 3)
+ REAL(prec) :: ddFF(3, 3, 3, 3, 3, 3), II_mat(3, 3, 3, 3)
+ REAL(prec) :: temp1(3, 3, 3, 3), temp2(3, 3, 3, 3)
+
+ ! Define 2nd order identity tensor
+ data I_mat(1,:) /1.D0, 0.D0, 0.D0/
+ data I_mat(2,:) /0.D0, 1.D0, 0.D0/
+ data I_mat(3,:) /0.D0, 0.D0, 1.D0/
+
+ ! Define 4th order symetric identity tensor
+ II_mat(:,:,:,:) = 0.D0
+ temp1(:,:,:,:) = 0.D0
+ temp2(:,:,:,:) = 0.D0
+ DO ii = 1, 3
+ DO jj = 1, 3
+ DO mm = 1, 3
+ DO ll = 1, 3
+ temp1(ii,jj,mm,ll) = temp1(ii,jj,mm,ll)
+ 1 + I_mat(ii,mm) * I_mat(jj,ll)
+ temp2(ii,jj,mm,ll) = temp2(ii,jj,mm,ll)
+ 1 + I_mat(ii,ll) * I_mat(jj,mm)
+ END DO
+ END DO
+ END DO
+ END DO
+ II_mat(:,:,:,:) = 0.5D0 * (temp1(:,:,:,:) + temp2(:,:,:,:))
+
+
+ ! Calculate coeffs for the power law approxiamation of log
+ DO ii = 1, order+1
+ IF (ii .EQ. 1) THEN
+ coeffs(ii) = 0.D0
+ ELSE IF (modulo(ii, 2) .EQ. 0) THEN
+ coeffs(ii) = 1.D0 / (real(ii, 8) - 1.D0)
+ ELSE
+ coeffs(ii) = -1.D0 / (real(ii, 8) - 1.D0)
+ END IF
+ END DO
+
+ ! Implementation of polynomial func as per dG3D
+ logAA(:, :) = 0.D0
+ DO ii = 1, 3
+ DO jj = 1, 3
+ logAA(ii, jj) = I_mat(ii, jj) * coeffs(order + 1)
+ END DO
+ END DO
+
+ dlogAA(:, :, :, :) = 0.D0
+ ddlogAA(:, :, :, :, :, :) = 0.D0
+
+ DO ii = 1, order
+ FF(:, :) = 0.D0
+ dFF(:, :, :, :) = 0.D0
+ ddFF(:, :, :, :, :, :) = 0.D0
+
+ nn = order + 1 - ii
+
+ DO jj = 1, 3
+ DO mm = 1, 3
+ FF(jj, mm) = coeffs(nn)*I_mat(jj, mm)
+
+ DO ll = 1, 3
+ FF(jj, mm) = FF(jj, mm) + logAA(jj, ll) * AA(ll, mm)
+
+ DO rr = 1, 3
+ DO ss = 1, 3
+ dFF(jj, mm, rr, ss) = dFF(jj, mm, rr, ss)
+ 1 + dlogAA(jj, ll, rr, ss) * AA(ll, mm)
+ 2 + logAA(jj, ll) * II_mat(ll, mm, rr, ss)
+
+ DO pp = 1, 3
+ DO qq = 1, 3
+ ddFF(jj, mm, rr, ss, pp, qq) = ddFF(jj, mm, rr, ss, pp, qq)
+ 1 + ddlogAA(jj, ll, rr, ss, pp, qq) * AA(ll, mm)
+ 2 + dlogAA(jj, ll, rr, ss) * II_mat(ll, mm, pp, qq)
+ 3 + dlogAA(jj, ll, pp, qq) * II_mat(ll, mm, rr, ss)
+ END DO
+ END DO
+
+ END DO
+ END DO
+
+ END DO
+
+ END DO
+ END DO
+ logAA(:, :) = FF(:, :)
+ dlogAA(:,:,:,:) = dFF(:,:,:,:)
+ ddlogAA(:,:,:,:,:,:) = ddFF(:, :, :, :, :, :)
+ END DO
+
+ END SUBROUTINE approx_log
+
+ !--------------------------------------------------------------
+ ! Function to compute exp of matrix
+ ! Power series approximation
+ !--------------------------------------------------------------
+
+ SUBROUTINE approx_exp(Q, order, expQ)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: order
+ REAL(prec), INTENT(IN) :: Q(3, 3)
+ REAL(prec), INTENT(OUT) :: expQ(3, 3)
+
+
+ INTEGER :: ii, jj, nn, mm, O5
+ REAL(prec) :: Q_temp(order, 3, 3), I_mat(3, 3)
+ REAL(prec) :: coeffs(order - 1)
+
+ !Define 2nd order identity tensor
+ data I_mat(1,:) /1.D0, 0.D0, 0.D0/
+ data I_mat(2,:) /0.D0, 1.D0, 0.D0/
+ data I_mat(3,:) /0.D0, 0.D0, 1.D0/
+
+
+ Q_temp(:, :, :) = 0.D0
+ Q_temp(1, :, :) = I_mat(:, :)
+ Q_temp(2, :, :) = Q(:, :)
+
+ DO O5 = 3, order
+ DO ii = 1, 3
+ DO jj = 1, 3
+ DO mm = 1, 3
+ Q_temp(O5, ii, jj) = Q_temp(O5, ii, jj)
+ 1 + Q_temp(O5 - 1, ii, mm) * Q(mm, jj)
+ END DO
+ END DO
+ END DO
+ END DO
+
+ coeffs(1) = 1.d0
+ DO O5 = 2, order
+ coeffs(O5) = coeffs(O5 - 1) * REAL(O5, 8)
+ END DO
+
+ expQ(:, :) = 0.D0
+ expQ(:, :) = Q_temp(1, :, :)
+
+ DO O5 = 2, order
+ expQ(:, :) = expQ(:, :)
+ 1 + Q_temp(O5, :, :) / coeffs(O5 - 1)
+ END DO
+ END SUBROUTINE approx_exp
+
+
+ !--------------------------------------------------------------
+ ! VE Predictor routine - 3 Maxwell Branches
+ !--------------------------------------------------------------
+
+ SUBROUTINE vePredictor_log(E_ve, E_ve_n, AA_n, BB_n,
+ 1 dt, PROPS, NPROPS, AA, BB, GGe, KKe, kappa)
+
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: NPROPS
+ REAL(prec), INTENT(IN) :: E_ve(3, 3), E_ve_n(3, 3)
+ REAL(prec), INTENT(IN) :: AA_n(8, 3, 3), BB_n(8), dt
+ REAL(prec), INTENT(IN) :: PROPS(NPROPS)
+ REAL(prec), INTENT(OUT) :: AA(8, 3, 3), BB(8), GGe, KKe
+ REAL(prec), INTENT(OUT) :: kappa(3, 3)
+
+ INTEGER :: ii
+ REAL(prec) :: dev_E_ve(3, 3), tr_E_ve, DE_ve(3, 3)
+ REAL(prec) :: dev_DE_ve(3, 3), tr_DE_ve, dtg, expmdtg, ztag
+ REAL(prec) :: dtk, expmdtk, ztak, dev_kappa(3, 3)
+ REAL(prec) :: p, I_mat(3, 3)
+
+ ! Material Params
+ REAL(prec) :: KK_inf, GG_inf
+ REAL(prec) :: KK(8), k(8), GG(8), g(8)
+
+ !Define 2nd order identity tensor
+ data I_mat(1,:) /1.D0, 0.D0, 0.D0/
+ data I_mat(2,:) /0.D0, 1.D0, 0.D0/
+ data I_mat(3,:) /0.D0, 0.D0, 1.D0/
+
+ ! Get Material Properties
+ KK_inf =PROPS(2)
+ GG_inf =PROPS(3)
+ KK(1) =PROPS(19)
+ KK(2) =PROPS(20)
+ KK(3) =PROPS(21)
+ KK(4) =PROPS(22)
+ KK(5) =PROPS(23)
+ KK(6) =PROPS(24)
+ KK(7) =PROPS(25)
+ KK(8) =PROPS(26)
+ k(1) =PROPS(27)
+ k(2) =PROPS(28)
+ k(3) =PROPS(29)
+ k(4) =PROPS(30)
+ k(5) =PROPS(31)
+ k(6) =PROPS(32)
+ k(7) =PROPS(33)
+ k(8) =PROPS(34)
+ GG(1) =PROPS(35)
+ GG(2) =PROPS(36)
+ GG(3) =PROPS(37)
+ GG(4) =PROPS(38)
+ GG(5) =PROPS(39)
+ GG(6) =PROPS(40)
+ GG(7) =PROPS(41)
+ GG(8) =PROPS(42)
+ g(1) =PROPS(43)
+ g(2) =PROPS(44)
+ g(3) =PROPS(45)
+ g(4) =PROPS(46)
+ g(5) =PROPS(47)
+ g(6) =PROPS(48)
+ g(7) =PROPS(49)
+ g(8) =PROPS(50)
+
+ CALL tr_dev_split(E_ve(:,:), dev_E_ve(:,:), tr_E_ve)
+
+ DE_ve(:, :) = E_ve(:, :) - E_ve_n(:, :)
+ CALL tr_dev_split(DE_ve(:,:), dev_DE_ve(:,:), tr_DE_ve)
+
+ GGe = GG_inf
+ AA(:,:,:)=0.D0
+
+ KKe = KK_inf
+ BB(:)=0.D0
+
+ DO ii = 1,8
+ dtg = dt / g(ii)
+ expmdtg = EXP(-dtg)
+ ztag = EXP(-dtg/2.D0)
+
+ GGe = GGe + GG(ii) * ztag
+ AA(ii,:,:) = expmdtg*AA_n(ii,:, :)
+ 1 + ztag * dev_DE_ve(:, :)
+
+
+ dtk = dt / k(ii)
+ expmdtk = EXP(-dtk)
+ ztak = EXP(-dtk/2.D0)
+
+ KKe = KKe + KK(ii) * ztak
+ BB(ii) = BB_n(ii) * expmdtk + ztak * tr_DE_ve
+ END DO
+
+
+ dev_kappa(:, :) = dev_E_ve(:, :) * 2.D0 * GG_inf
+ p = tr_E_ve * KK_inf
+
+ DO ii = 1, 8
+ dev_kappa(:, :) = dev_kappa(:, :)
+ 1 + 2.D0 * GG(ii) * AA(ii, :, :)
+
+ p = p + BB(ii) * KK(ii)
+ END DO
+
+ kappa(:, :) = dev_kappa(:, :) + p * I_mat(:, :)
+
+ RETURN
+ END SUBROUTINE vePredictor_log
+
+ !--------------------------------------------------------------
+ ! Iso Hardening Function (Linear exponential)
+ !--------------------------------------------------------------
+
+ SUBROUTINE hardn(PROPS, NPROPS, gma_0, gma,
+ 1 sigma_c, sigma_t, HHc, HHt, HHb)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: NPROPS
+ REAL(prec), INTENT(IN) :: PROPS(NPROPS), gma_0, gma
+ REAL(prec), INTENT(OUT) :: sigma_c, sigma_t, HHc, HHt,HHb
+
+ REAL(prec) :: sigmac0, hc1, hc2, hcexp, tempc
+ REAL(prec) :: sigmat0, ht1, ht2, htexp, tempt
+ REAL(prec) :: hb0, hb1, hb2
+
+ sigmac0 =PROPS(8)
+ hc1 =PROPS(9)
+ hc2 =PROPS(10)
+ hcexp =PROPS(11)
+ sigmat0 =PROPS(12)
+ ht1 =PROPS(13)
+ ht2 =PROPS(14)
+ htexp =PROPS(15)
+ hb0 =PROPS(16)
+ hb1 =PROPS(17)
+ hb2 =PROPS(18)
+
+ tempc = EXP(-hcexp * gma)
+ sigma_c = sigmac0 + hc1 * gma + hc2*(1.D0 - tempc)
+ HHc = hc1 + hc2 * hcexp * tempc
+
+ tempt = EXP(-htexp * gma)
+ sigma_t = sigmat0 + ht1 * gma + ht2*(1.D0 - tempt)
+ HHt = ht1 + ht2 * htexp * tempt
+
+ HHb = hb0 + hb1*gma_0 + hb2*gma_0*gma_0
+ END SUBROUTINE hardn
+
+
+ !--------------------------------------------------------------
+ ! Drucker - Prager coefficient update subroutine
+ !--------------------------------------------------------------
+ SUBROUTINE DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.D0)
+ INTEGER, PARAMETER :: prec=double
+
+ REAL(prec), INTENT(IN) :: alpha, sigma_c, sigma_t
+ REAL(prec), INTENT(OUT) :: m, a0, a1, a2
+
+ REAL(prec) :: mral, mp1, sigcral
+
+ m = sigma_t / sigma_c
+
+ mral = m**alpha
+ mp1 = m + 1.D0
+ sigcral = sigma_c**alpha
+
+ a2 = 1.D0 / sigcral
+
+ a1 = 3.D0 * ( (mral - 1.D0)/(mp1) ) * (1.D0 / sigma_c)
+
+ a0 = (mral + m) / mp1
+
+ END SUBROUTINE DPcoeff
+
+ !--------------------------------------------------------------
+ ! NonLinear solver (Newton Raphson)
+ !--------------------------------------------------------------
+
+ SUBROUTINE nlinSolver(PROPS, NPROPS, DTIME, TOLL, MAX_i,
+ 1 GGe, KKe, F, gma_0, gma,
+ 2 ptilde, PhiEq, GAMMA, u, v, beta, k_plast)
+
+ IMPLICIT NONE
+ INTEGER, PARAMETER :: double=kind(1.d0)
+ INTEGER, PARAMETER :: prec=double
+
+ INTEGER, INTENT(IN) :: NPROPS
+ REAL(prec), INTENT(IN) :: PROPS(NPROPS), DTIME, TOLL
+ INTEGER, INTENT(IN) :: MAX_i
+ REAL(prec), INTENT(IN) :: GGe, KKe
+ REAL(prec), INTENT(OUT) :: GAMMA, beta, k_plast, u, v
+ REAL(prec), INTENT(INOUT) :: F, gma_0, gma, ptilde, PhiEq
+
+ INTEGER :: iter_G
+ REAL(prec) :: gma_i, A, nu_p, GG_til, KK_til, sigma_c
+ REAL(prec) :: phi_e_tr, phi_p_tr, alpha, eta, p_exp
+ REAL(prec) :: sigma_t, HHc, HHt, HHb
+ REAL(prec) :: dAdGamma, dDgmadGamma, Dm, Da1Dm
+ REAL(prec) :: m, a0, a1, a2, H2, H1, H0, etaOverDt
+ REAL(prec) :: dfdDgma, DfDGamma, dGamma, Dgma, viscoterm
+
+ ! Material Props
+ alpha = PROPS(4)
+ nu_p = PROPS(5)
+ eta = PROPS(6)
+ p_exp = PROPS(7)
+
+ ! Calc derived material props
+ beta = (9.D0 / 2.D0)
+ 1 * ((1.D0 - 2.D0 * nu_p) / (nu_p + 1.D0))
+
+ k_plast = 1.D0 / ((1.D0 + 2.D0 * nu_p ** 2.D0)**(0.5D0))
+
+
+
+ phi_e_tr = PhiEq
+ phi_p_tr = ptilde
+ ! Initialize for first step
+ iter_G = 0
+
+ GAMMA = 0.0
+ u = 1.D0
+ v = 1.D0
+
+ A = (6.D0 * (PhiEq**2.D0)
+ 1 + (4.D0 / 3.D0) * (beta**2.D0)
+ 2 * (ptilde ** 2.D0)) ** 0.5D0
+
+ gma_i = gma
+
+ ! Update hardening variables on trial state
+ CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c, sigma_t,
+ 1 HHc, HHt, HHb)
+ ! Update Drucker Pragger coefficients along with ecc factor m
+ CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
+
+ GG_til = GGe + (k_plast/2.D0)*HHb
+ KK_til = KKe + (k_plast/3.D0)*HHb
+
+ DO WHILE (((ABS(F) .GE. TOLL) .OR. (iter_G .LT. 1))
+ 1 .AND. (iter_G .LE. MAX_i))
+
+ etaOverDt = eta / DTIME
+
+ dAdGamma = -(72.D0 * GG_til * PhiEq * PhiEq /u
+ 1 + 16.D0 * KK_til * beta*beta*beta*ptilde*ptilde
+ 2 /(3.D0*v)) / (2.D0 * A)
+
+ dDgmadGamma = k_plast*(A + GAMMA * dAdGamma)
+
+ Dm = (HHt*sigma_c - HHc*sigma_t) / (sigma_c * sigma_c)
+
+ Da1Dm = (3.D0 / sigma_c)
+ 1 * (alpha * (m**(alpha - 1.D0))/(m + 1.D0)
+ 2 - (((m**alpha) - 1.D0)/(m + 1.D0))/(m + 1.D0))
+
+ H2 = -alpha * (sigma_c**(-alpha - 1.D0)) * HHc
+
+ H1 = Da1Dm * Dm
+ 1 - 3.D0 * ((((m**alpha) - 1.D0) / (m + 1.D0))
+ 2 /(sigma_c * sigma_c))*HHc
+
+ H0 =((alpha * (m**(alpha - 1.D0)) + 1.D0) / (m + 1.D0)
+ 1 - (((m**alpha) + m)/(m + 1.D0))/(m + 1.D0))*Dm
+
+
+ dfdDgma = H2 * (PhiEq**alpha) -H1*ptilde - H0
+
+ DfDGamma = (dfdDgma * dDgmadGamma)
+ 1 - (alpha * a2 * 6.D0 * GG_til) * (PhiEq**alpha) / u
+ 2 + a1 * ptilde * 2.D0 * beta * KK_til / v
+
+ IF ((GAMMA .GT. 0.D0) .AND. (etaOverDt .GT. 0.D0)) THEN
+ DfDGamma = DfDGamma - (etaOverDt**p_exp) * p_exp
+ 1 * (GAMMA**(p_exp - 1.D0))
+ END IF
+
+ dGamma = -F / DfDGamma
+
+ IF ((GAMMA + dGamma) .LE. 0.D0) THEN
+ GAMMA = GAMMA / 2.D0
+
+ ELSE
+ GAMMA = GAMMA + dGamma
+ END IF
+
+ u = 1.D0 + 6.D0 * GG_til * GAMMA
+ v = 1.D0 + 2.D0 * beta * KK_til * GAMMA
+
+ PhiEq = phi_e_tr / u
+ ptilde = phi_p_tr / v
+
+ A = (6.D0 * (PhiEq**2.D0)
+ 1 + (4.D0 / 3.D0) * (beta**2.D0)
+ 2 * (ptilde ** 2.D0)) ** 0.5D0
+
+ Dgma = k_plast * GAMMA * A
+ gma_i = gma + Dgma
+
+ ! Update hardening variables on gma_i
+ CALL hardn(PROPS, NPROPS, gma_0, gma_i, sigma_c, sigma_t,
+ 1 HHc, HHt, HHb)
+ ! Update Drucker Pragger coefficients along with ecc factor m
+ CALL DPcoeff(alpha, sigma_c, sigma_t, m, a0, a1, a2)
+
+ GG_til = GGe + (k_plast/2.D0)*HHb
+ KK_til = KKe + (k_plast/3.D0)*HHb
+
+ F = a2 * PhiEq**alpha - a1 * ptilde - a0
+ viscoterm = etaOverDt * GAMMA
+
+ IF ((GAMMA .GT. 0.D0) .AND. (etaOverDt .GT. 0.D0)) THEN
+ F = F - viscoterm ** p_exp
+ END IF
+
+ iter_G = iter_G + 1
+
+ END DO
+
+
+
+ gma = gma_i
+ ! Debuging statements for N.R. Comment out for big sims
+ ! WRITE(*,*) "Iter : ", iter_G
+ ! WRITE(*,*) "GAMMA : ", GAMMA
+ ! WRITE(*,*) "RESI : ", ABS(F)
+ ! WRITE(*,*) "MAX_i : ", MAX_i
+ ! WRITE(*,*) "-------------------------------------------"
+
+ END SUBROUTINE nlinSolver
diff --git a/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp b/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp
index c2a406fe346c16d2216dd3053adf3305cbf3d409..5eed9a420edd39ad4d86e42e9e740cfbfe6e342c 100644
--- a/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp
+++ b/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp
@@ -15,14 +15,17 @@
#include <math.h>
#include "OS.h"
#include "ipFiniteStrain.h"
+#include "TrainingDeepMaterialNetworksNonLinear.h"
#if defined(HAVE_PETSC)
#include "nonLinearSystemPETSc.h"
#elif defined(HAVE_GMM)
#include "nonLinearSystemGmm.h"
#endif //HAVE_GMM
+#include <filesystem>
DeepMaterialNetwork::DeepMaterialNetwork(int tag): _isInitialized(false), _tag(tag),
-_ownInteractionData(false),_networkInteraction(NULL),_enumMinus(-1),_enumPlus(-1), _gnum(0)
+_ownInteractionData(false),_networkInteraction(NULL),_enumMinus(-1),_enumPlus(-1), _gnum(0),
+_saveDir("")
#if defined(HAVE_PETSC)
,_isAddedSparsityPattern(false)
#endif //HAVE_PETSC
@@ -33,7 +36,7 @@ _ownInteractionData(false),_networkInteraction(NULL),_enumMinus(-1),_enumPlus(-1
};
DeepMaterialNetwork::DeepMaterialNetwork(const Tree& T, int tag): _isInitialized(false), _tag(tag), _ownInteractionData(true),
-_enumMinus(-1),_enumPlus(-1), _gnum(0)
+_enumMinus(-1),_enumPlus(-1), _gnum(0), _saveDir("")
#if defined(HAVE_PETSC)
,_isAddedSparsityPattern(false)
#endif //HAVE_PETSC
@@ -46,7 +49,7 @@ _enumMinus(-1),_enumPlus(-1), _gnum(0)
};
DeepMaterialNetwork::DeepMaterialNetwork(const std::string fname, int tag ):_isInitialized(false), _tag(tag), _ownInteractionData(true),
-_enumMinus(-1),_enumPlus(-1), _gnum(0)
+_enumMinus(-1),_enumPlus(-1), _gnum(0), _saveDir("")
#if defined(HAVE_PETSC)
,_isAddedSparsityPattern(false)
#endif //HAVE_PETSC
@@ -79,6 +82,7 @@ void DeepMaterialNetwork::copyViewSetting(const DeepMaterialNetwork& src)
{
_viewIP = src._viewIP;
_averagePhase = src._averagePhase;
+ _saveDir = src._saveDir;
};
void DeepMaterialNetwork::setMicroProblemIndentification(int ele, int gpt){
@@ -93,23 +97,53 @@ void DeepMaterialNetwork::setMicroProblemIndentification(int eleMinus, int elePl
_gnum = gpt;
};
+void DeepMaterialNetwork::setSaveDir(std::string saveDir)
+{
+ Msg::Info("result will be saved in %s",saveDir.c_str());
+ _saveDir = saveDir;
+}
+
+int DeepMaterialNetwork::evaluateByPath(double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath)
+{
+ return evaluatePath(this, maximumStrainStep, strainPath, stressPath);
+};
+
std::string DeepMaterialNetwork::getFileSavingPrefix() const
{
+ if (_saveDir.size() >0)
+ {
+ bool fl = std::filesystem::is_directory(_saveDir.c_str());
+ if (!fl)
+ {
+ std::filesystem::create_directory(_saveDir.c_str());
+ Msg::Info("directory %s is created",_saveDir.c_str());
+ }
+ }
+
if ((_enumMinus != -1) && (_enumPlus != -1))
{
if (_enumMinus == _enumPlus)
{
- return "E_"+int2str(_enumMinus)+"_GP_"+int2str(_gnum)+"_";
+ if (_saveDir.size() ==0)
+ return "E_"+int2str(_enumMinus)+"_GP_"+int2str(_gnum)+"_";
+ else
+ return _saveDir+"/E_"+int2str(_enumMinus)+"_GP_"+int2str(_gnum)+"_";
}
else
{
- return "E_Interface_"+int2str(_enumMinus)+"_"+int2str(_enumPlus)+"_GP_"+int2str(_gnum)+"_";
+ if (_saveDir.size() ==0)
+ return "E_Interface_"+int2str(_enumMinus)+"_"+int2str(_enumPlus)+"_GP_"+int2str(_gnum)+"_";
+ else
+ return _saveDir+"/E_Interface_"+int2str(_enumMinus)+"_"+int2str(_enumPlus)+"_GP_"+int2str(_gnum)+"_";
}
}
else
{
- return "";
+ if (_saveDir.size() ==0)
+ return "";
+ else
+ return _saveDir+"/";
}
};
diff --git a/NonLinearSolver/modelReduction/DeepMaterialNetworks.h b/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
index 10a1203397bf0c64fa8ab0cc585548344f28c578..93f08583e1b59b562413076e92ff2d0e38582529 100644
--- a/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
+++ b/NonLinearSolver/modelReduction/DeepMaterialNetworks.h
@@ -26,6 +26,7 @@ template<class T>
class nonLinearSystemGmm;
#endif //HAVE_GMM
#include "nlsField.h"
+class PathSTensor3;
#endif //SWIG
#include "Tree.h"
#include "NetworkInteraction.h"
@@ -115,6 +116,7 @@ class DeepMaterialNetwork
std::vector<AveragePhase> _averagePhase;
// Element number and Integration point number
int _enumMinus,_enumPlus, _gnum;
+ std::string _saveDir;
#endif //SWIG
public:
@@ -132,7 +134,8 @@ class DeepMaterialNetwork
std::string getFileSavingPrefix() const;
void setMicroProblemIndentification(int ele, int gpt);
void setMicroProblemIndentification(int eleMinus, int elePlus, int gpt);
-
+ void setSaveDir(std::string saveDir);
+ int evaluateByPath(double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath);
void addInteraction(NetworkInteraction& inter);
void saveInteractionToFile(const std::string fname) const;
diff --git a/NonLinearSolver/modelReduction/NetworkInteraction.h b/NonLinearSolver/modelReduction/NetworkInteraction.h
index 3e1312070e45390a29bdd2292c616c8c0055c960..d05f611dcd84c139cdeae806e19ecb59153d6e8b 100644
--- a/NonLinearSolver/modelReduction/NetworkInteraction.h
+++ b/NonLinearSolver/modelReduction/NetworkInteraction.h
@@ -228,6 +228,24 @@ class CoefficientReduction
#endif //SWIG
};
+
+class CoefficientReductionVoid: public CoefficientReduction
+{
+ public:
+ #ifndef SWIG
+ CoefficientReductionVoid(){}
+ CoefficientReductionVoid(const CoefficientReductionVoid& src){}
+ virtual ~CoefficientReductionVoid(){}
+ virtual bool withCoefficientVars(const NetworkInteraction& interaction, const InteractionMechanism& im) const {return false;};
+ virtual int getNumberOfCoefficientVars(const NetworkInteraction& interaction, const InteractionMechanism& im) const {return 0;};
+ virtual void updateCoefficientVars(const NetworkInteraction& interaction, InteractionMechanism& im) const {};
+ virtual void updateCoefficientFromCoefficientVars(const NetworkInteraction& interaction, InteractionMechanism& im) const {};
+ virtual void getDCoefficientsDWeights(const NetworkInteraction& interaction, const InteractionMechanism& im, int nodeId, std::vector<int>&allNodes, std::vector<double>& DalphaDW) const {};
+ virtual void getDCoefficientsDCoefficientVars(const NetworkInteraction& interaction, const InteractionMechanism& im, int nodeId, std::vector<double>& DalphaVars) const {};
+ virtual CoefficientReduction* clone() const {return new CoefficientReductionVoid();};
+ #endif //SWIG
+};
+
class CoefficientReductionAll : public CoefficientReduction
{
public:
diff --git a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp
index 6f2d9b1067c9c0b2e322a8666662c3ef412a2ba5..4704a0705d63aee8357c4d3090ced332012d0e73 100644
--- a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp
+++ b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.cpp
@@ -151,6 +151,16 @@ void PathSTensor3::setVal(int step, int i, int j, double val)
F(i,j) = val;
};
+double PathSTensor3::getValue(int step, int i, int j) const
+{
+ if (step > trueSize-1)
+ {
+ Msg::Error("step %s exceeds the path length %d",step,trueSize);
+ return 0;
+ }
+ return value[step](i,j);
+};
+
void PathSTensor3::fromFile(const std::string filename)
{
clear();
@@ -240,6 +250,17 @@ void PathSTensor3::copyTo(PathSTensor3& dest, int space) const
}
};
+PathDouble PathSTensor3::determinant() const
+{
+ int Nelements = this->size();
+ PathDouble det(Nelements);
+ for (int i=0; i< Nelements; i++)
+ {
+ det.getValue(i) = STensorOperation::determinantSTensor3(this->getValue(i));
+ };
+ return det;
+};
+
//
PathDouble::PathDouble(int N): trueSize(N), value(N,0.){}
PathDouble::~PathDouble(){};
@@ -266,15 +287,7 @@ double& PathDouble::getValue(int step)
}
return value[step];
};
-const double& PathDouble::getValue(int step) const
-{
- if (step > trueSize-1)
- {
- Msg::Error("size %d exceed the maximum index %d",step,trueSize-1);
- Msg::Exit(0);
- }
- return value[step];
-};
+
std::vector<double>& PathDouble::getDValueDWeight(int step)
{
if (step > trueSize-1)
@@ -378,6 +391,24 @@ void PathDouble::setVal(int step, double val)
value[step] = val;
};
+void PathDouble::addValue(double v)
+{
+ for (int i=0; i< trueSize; i++)
+ {
+ value[i] += v;
+ }
+}
+
+double PathDouble::getValue(int step) const
+{
+ if (step > trueSize-1)
+ {
+ Msg::Error("step %s exceeds the path length %d",step,trueSize);
+ return 0;
+ }
+ return value[step];
+}
+
void PathDouble::fromFile(const std::string filename)
{
clear();
@@ -488,6 +519,16 @@ double RelativeSquareLossPathMeasure::get(const PathSTensor3& Pref, const PathST
{
DerrorDP->allocate(size);
}
+ double averageVal = 0;
+ for (int i=0; i< sizePref; i++)
+ {
+ averageVal += 0.5*STensorOperation::doubledot(Pref.getValue(i),Pref.getValue(i));
+ }
+ averageVal/=double(sizePref);
+ if (averageVal ==0)
+ {
+ averageVal = 1.; // to avoid any reference value
+ }
//
for (int i=0; i< size; i++)
{
@@ -497,12 +538,18 @@ double RelativeSquareLossPathMeasure::get(const PathSTensor3& Pref, const PathST
double valSq = 0.5*STensorOperation::doubledot(dP,dP);
double vR = 0.5*STensorOperation::doubledot(Pref.getValue(i),Pref.getValue(i));
- error += valSq/(size*vR);
+ if (vR > 0)
+ error += valSq/(size*vR);
+ else
+ error += valSq/(size*averageVal);
if (stiff)
{
STensor3& DerrorDPj = DerrorDP->getValue(i);
DerrorDPj = dP;
- DerrorDPj *= (1./(vR*size));
+ if (vR >0)
+ DerrorDPj *= (1./(vR*size));
+ else
+ DerrorDPj *= (1./(averageVal*size));
}
};
return error;
@@ -519,17 +566,34 @@ double RelativeSquareLossPathMeasure::get(const PathDouble& Pref, const PathDoub
DerrorDP->allocate(size);
}
//
+ double averageVal = 0;
+ for (int i=0; i< sizePref; i++)
+ {
+ averageVal += 0.5*(Pref.getValue(i)*Pref.getValue(i));
+ }
+ averageVal/=double(sizePref);
+ if (averageVal ==0)
+ {
+ averageVal = 1.; // to avoid any reference value
+ }
+
for (int i=0; i< size; i++)
{
double dP = P.getValue(i)-Pref.getValue(i);
double valSq = 0.5*dP*dP;
double vR = 0.5*(Pref.getValue(i)*Pref.getValue(i));
- error += valSq/(size*vR);
+ if (vR >0)
+ error += valSq/(size*vR);
+ else
+ error += valSq/(size*averageVal);
if (stiff)
{
double& DerrorDPj = DerrorDP->getValue(i);
DerrorDPj = dP;
- DerrorDPj *= (1./(vR*size));
+ if (vR >0)
+ DerrorDPj *= (1./(vR*size));
+ else
+ DerrorDPj *= (1./(averageVal*size));
}
};
return error;
@@ -598,6 +662,17 @@ double RelativeAbsoluteLossPathMeasure::get(const PathSTensor3& Pref, const Path
DerrorDP->allocate(size);
}
//
+ double averageVal = 0;
+ for (int i=0; i< sizePref; i++)
+ {
+ averageVal += sqrt(0.5*STensorOperation::doubledot(Pref.getValue(i),Pref.getValue(i)));
+ };
+ averageVal /= double(sizePref);
+ if (averageVal ==0)
+ {
+ averageVal = 1.; // to avoid any reference value
+ }
+
for (int i=0; i< size; i++)
{
static STensor3 dP;
@@ -605,12 +680,24 @@ double RelativeAbsoluteLossPathMeasure::get(const PathSTensor3& Pref, const Path
dP -= Pref.getValue(i);
double val = sqrt(0.5*STensorOperation::doubledot(dP,dP));
double vR = sqrt(0.5*STensorOperation::doubledot(Pref.getValue(i),Pref.getValue(i)));
- error += val/(vR*size);
+ if (val > 0)
+ {
+ if (vR >0)
+ error += val/(vR*size);
+ else
+ error += val/(averageVal*size);
+ }
if (stiff)
{
STensor3& DerrorDPj = DerrorDP->getValue(i);
DerrorDPj = dP;
- DerrorDPj *= (1./(2.*val*vR*size));
+ if (val > 0)
+ {
+ if (vR >0)
+ DerrorDPj *= (1./(2.*val*vR*size));
+ else
+ DerrorDPj *= (1./(2.*val*averageVal*size));
+ }
};
};
@@ -628,17 +715,40 @@ double RelativeAbsoluteLossPathMeasure::get(const PathDouble& Pref, const PathDo
DerrorDP->allocate(size);
}
//
+ double averageVal = 0;
+ for (int i=0; i< sizePref; i++)
+ {
+ averageVal += sqrt(0.5*(Pref.getValue(i)*Pref.getValue(i)));
+ };
+ averageVal /= double(sizePref);
+ if (averageVal ==0)
+ {
+ averageVal = 1.; // to avoid any reference value
+ }
+
for (int i=0; i< size; i++)
{
double dP = P.getValue(i) - Pref.getValue(i);
double val = sqrt(0.5*dP*dP);
double vR = sqrt(0.5*(Pref.getValue(i)*Pref.getValue(i)));
- error += val/(vR*size);
+ if (val > 0)
+ {
+ if (vR > 0)
+ error += val/(vR*size);
+ else
+ error += val/(averageVal*size);
+ }
if (stiff)
{
double& DerrorDPj = DerrorDP->getValue(i);
DerrorDPj = dP;
- DerrorDPj *= (1./(2.*val*vR*size));
+ if (val > 0)
+ {
+ if (vR >0)
+ DerrorDPj *= (1./(2.*val*vR*size));
+ else
+ DerrorDPj *= (1./(2.*val*averageVal*size));
+ }
};
};
@@ -667,7 +777,8 @@ double AbsoluteLossPathMeasure::get(const PathSTensor3& Pref, const PathSTensor3
{
STensor3& DerrorDPj = DerrorDP->getValue(i);
DerrorDPj = dP;
- DerrorDPj *= (1./(2.*val*size));
+ if (val >0)
+ DerrorDPj *= (1./(2.*val*size));
}
}
return error;
@@ -693,7 +804,8 @@ double AbsoluteLossPathMeasure::get(const PathDouble& Pref, const PathDouble& P,
{
double& DerrorDPj = DerrorDP->getValue(i);
DerrorDPj = dP;
- DerrorDPj *= (1./(2.*val*size));
+ if (val > 0)
+ DerrorDPj *= (1./(2.*val*size));
}
}
return error;
@@ -702,8 +814,9 @@ double AbsoluteLossPathMeasure::get(const PathDouble& Pref, const PathDouble& P,
TrainingDeepMaterialNetworkNonLinear::TrainingDeepMaterialNetworkNonLinear(DeepMaterialNetwork& dmn, const CoefficientReduction& reduction):
TrainingArbitraryDeepMaterialNetwork(dmn.getRefToNetworkInteraction(),reduction),
_trainingWithPhaseDeformation(false),
- _phaseIpComForTraining(IPField::UNDEFINED),
- _phaseIndexForTraining(0),
+ _phaseIpComForTraining(),
+ _phaseIndexForIpComForTraining(),
+ _phaseIndexForTrainingPhaseDeformation(0),
_mixedFactor(0.),
_mixedFactorIPComp(0.),
_DMN(&dmn),
@@ -727,7 +840,7 @@ void TrainingDeepMaterialNetworkNonLinear::setMaximumStrainStep(double val)
int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strainPath, PathSTensor3& stressPath, bool messageActive, bool stiff,
bool extractPhaseDefo, PathSTensor3* phaseDefoPath,
- bool extractIPComp, PathDouble* IpCompPath)
+ bool extractIPComp, std::vector<PathDouble>* IpCompPath)
{
double tstart = Cpu();
//
@@ -751,6 +864,7 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
int sizeOfWeights = 0;
int sizeOfAlphas = 0;
int sizeOfNormals = 0;
+ int sizeIpCompPath= _phaseIpComForTraining.size();
if (stiff)
{
@@ -768,7 +882,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp)> 0.)
{
- IpCompPath->allocateWithDerivatives(pathLength,sizeOfWeights, sizeOfAlphas, sizeOfNormals);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*IpCompPath)[k].allocateWithDerivatives(pathLength,sizeOfWeights, sizeOfAlphas, sizeOfNormals);
+ }
}
}
else
@@ -780,7 +897,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- IpCompPath->allocate(pathLength);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*IpCompPath)[k].allocate(pathLength);
+ }
}
}
for (int ip = 0; ip <pathLength; ip++)
@@ -823,9 +943,9 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
std::vector<STensor3>* DFPhasecurDAlpha = NULL;
std::vector<STensor3>* DFPhasecurDNormal = NULL;
- std::vector<double>* DIPCompValcurDW = NULL;
- std::vector<double>* DIPCompValcurDAlpha = NULL;
- std::vector<double>* DIPCompValcurDNormal = NULL;
+ std::vector<std::vector<double>*> DIPCompValcurDW(sizeIpCompPath,NULL);
+ std::vector<std::vector<double>*> DIPCompValcurDAlpha(sizeIpCompPath,NULL);
+ std::vector<std::vector<double>*> DIPCompValcurDNormal(sizeIpCompPath,NULL);
if (stiff)
{
DPcurDW = &(stressPath.getDValueDWeight(ip));
@@ -839,9 +959,12 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp)> 0.)
{
- DIPCompValcurDW = &(IpCompPath->getDValueDWeight(ip));
- DIPCompValcurDAlpha = &(IpCompPath->getDValueDCoefficient(ip));
- DIPCompValcurDNormal = &(IpCompPath->getDValueDNormal(ip));
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ DIPCompValcurDW[k] = &((*IpCompPath)[k].getDValueDWeight(ip));
+ DIPCompValcurDAlpha[k] = &((*IpCompPath)[k].getDValueDCoefficient(ip));
+ DIPCompValcurDNormal[k] = &((*IpCompPath)[k].getDValueDNormal(ip));
+ }
}
};
@@ -920,7 +1043,7 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
// derivatives should be here
if (extractPhaseDefo && fabs(_mixedFactor) > 0.)
{
- _DMN->getDeformationGradientPhase(_phaseIndexForTraining,phaseDefoPath->getValue(ip));
+ _DMN->getDeformationGradientPhase(_phaseIndexForTrainingPhaseDeformation,phaseDefoPath->getValue(ip));
phaseDefoPath->getValue(ip)(0,0) -= 1.;
phaseDefoPath->getValue(ip)(1,1) -= 1.;
phaseDefoPath->getValue(ip)(2,2) -= 1.;
@@ -928,7 +1051,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- IpCompPath->getValue(ip) = _DMN->getPerPhase(_phaseIndexForTraining,_phaseIpComForTraining);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*IpCompPath)[k].getValue(ip) = _ipOffsets[k]+_DMN->getPerPhase(_phaseIndexForIpComForTraining[k],_phaseIpComForTraining[k]);
+ }
};
if (stiff)
@@ -947,7 +1073,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- STensorOperation::zero((*DIPCompValcurDW)[iv]);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ STensorOperation::zero((*DIPCompValcurDW[k])[iv]);
+ }
}
}
@@ -962,7 +1091,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- STensorOperation::zero((*DIPCompValcurDAlpha)[iv]);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ STensorOperation::zero((*DIPCompValcurDAlpha[k])[iv]);
+ }
}
}
// zero first
@@ -977,7 +1109,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- STensorOperation::zero((*DIPCompValcurDNormal)[iv]);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ STensorOperation::zero((*DIPCompValcurDNormal[k])[iv]);
+ }
}
};
//
@@ -986,16 +1121,20 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
if (!_fixedWeightsOfMaterialNodes)
{
- static fullMatrix<double> DPDW, DFiDW, DIPCompDW;
+ static fullMatrix<double> DPDW, DFiDW;
+ std::vector<fullMatrix<double> > DIPCompDW(sizeIpCompPath);
static std::vector<int> allNodeId;
_DMN->computeDStressDWeights(*_coeffReduction,allNodeId,DPDW);
if (extractPhaseDefo && fabs(_mixedFactor) > 0.)
{
- _DMN->computeDStrainDWeights(_phaseIndexForTraining,*_coeffReduction,allNodeId,DFiDW);
+ _DMN->computeDStrainDWeights(_phaseIndexForTrainingPhaseDeformation,*_coeffReduction,allNodeId,DFiDW);
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- _DMN->computeDIPCompDWeights(_phaseIpComForTraining,_phaseIndexForTraining,*_coeffReduction,allNodeId,DIPCompDW);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ _DMN->computeDIPCompDWeights(_phaseIpComForTraining[k],_phaseIndexForIpComForTraining[k],*_coeffReduction,allNodeId,DIPCompDW[k]);
+ }
}
if (_fixedPhaseFraction)
{
@@ -1074,7 +1213,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- (*DIPCompValcurDW)[num] += DIPCompDW(0,j)*DfiDZj[j]/fPhase;
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*DIPCompValcurDW[k])[num] += DIPCompDW[k](0,j)*DfiDZj[j]/fPhase;
+ }
}
}
}
@@ -1109,7 +1251,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- (*DIPCompValcurDW)[num] += DIPCompDW(0,j)*DWDVar;
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*DIPCompValcurDW[k])[num] += DIPCompDW[k](0,j)*DWDVar;
+ }
}
};
}
@@ -1123,15 +1268,18 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
const InteractionMechanism* im = iinter->second;
static fullMatrix<double> DPcurDCoeffVars, DPcurDNorVars;
static fullMatrix<double> DFPhasecurDCoeffVars, DFPhasecurDNorVars;
- static fullMatrix<double> DIPCompValcurDCoeffVars, DIPCompValcurDNorVars;
+ std::vector<fullMatrix<double> >DIPCompValcurDCoeffVars(sizeIpCompPath), DIPCompValcurDNorVars(sizeIpCompPath);
_DMN->computeDStressDInteraction(*_coeffReduction, *im, DPcurDCoeffVars,DPcurDNorVars);
if (extractPhaseDefo&& fabs(_mixedFactor) > 0.)
{
- _DMN->computeDStrainDInteraction(_phaseIndexForTraining, *_coeffReduction, *im, DFPhasecurDCoeffVars,DFPhasecurDNorVars);
+ _DMN->computeDStrainDInteraction(_phaseIndexForTrainingPhaseDeformation, *_coeffReduction, *im, DFPhasecurDCoeffVars,DFPhasecurDNorVars);
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- _DMN->computeDIPCompDInteraction(_phaseIpComForTraining,_phaseIndexForTraining, *_coeffReduction, *im, DIPCompValcurDCoeffVars,DIPCompValcurDNorVars);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ _DMN->computeDIPCompDInteraction(_phaseIpComForTraining[k],_phaseIndexForIpComForTraining[k], *_coeffReduction, *im, DIPCompValcurDCoeffVars[k],DIPCompValcurDNorVars[k]);
+ }
}
std::vector<Dof> R;
if (_coeffReduction->getNumberOfCoefficientVars(*_interaction, *im) > 0)
@@ -1160,7 +1308,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- (*DIPCompValcurDAlpha)[num] += DIPCompValcurDCoeffVars(0,j);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*DIPCompValcurDAlpha[k])[num] += DIPCompValcurDCoeffVars[k](0,j);
+ }
}
}
}
@@ -1190,7 +1341,10 @@ int TrainingDeepMaterialNetworkNonLinear::evaluatePath(const PathSTensor3& strai
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- (*DIPCompValcurDNormal)[num] += DIPCompValcurDNorVars(0,j);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ (*DIPCompValcurDNormal[k])[num] += DIPCompValcurDNorVars[k](0,j);
+ }
}
};
}//
@@ -1227,17 +1381,17 @@ void TrainingDeepMaterialNetworkNonLinear::assembleOneDof(const Dof& key, double
void TrainingDeepMaterialNetworkNonLinear::computeErrorVec(std::vector<fullVector<double> >& allErrors,
const std::vector<LossPathMeasure*>& lossMeasures, const PathSTensor3& Fref, const PathSTensor3& Pref,
bool extractPhaseDefo, PathSTensor3* phaseDefoPathRef,
- bool extractIPComp, PathDouble* IpCompPathRef)
+ bool extractIPComp, std::vector<PathDouble>* IpCompPathRef)
{
if (lossMeasures.size() == 0)
{
allErrors.clear();
return;
}
-
+ int sizeIpCompPath= _phaseIpComForTraining.size();
allErrors.resize(lossMeasures.size(),fullVector<double>(3));
static PathSTensor3 P, Fi;
- static PathDouble IPCompi;
+ std::vector<PathDouble> IPCompi(sizeIpCompPath);
bool messageActive = false;
int lastConvergedIter = evaluatePath(Fref,P,messageActive,false,extractPhaseDefo,&Fi, extractIPComp, &IPCompi);
@@ -1252,7 +1406,10 @@ void TrainingDeepMaterialNetworkNonLinear::computeErrorVec(std::vector<fullVecto
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- allErrors[j](2) =lossMeasures[j]->get(*IpCompPathRef,IPCompi,lastConvergedIter+1);
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ allErrors[j](2) += (lossMeasures[j]->get((*IpCompPathRef)[k],IPCompi[k],lastConvergedIter+1))/(double)sizeIpCompPath;
+ }
}
}
};
@@ -1264,7 +1421,7 @@ double TrainingDeepMaterialNetworkNonLinear::getEquivalentError(const fullVector
{
val += _mixedFactor*allErr(1);
}
- if (_phaseIpComForTraining != IPField::UNDEFINED && fabs(_mixedFactorIPComp) > 0.)
+ if (_phaseIpComForTraining.size() >0 && fabs(_mixedFactorIPComp) > 0.)
{
val += _mixedFactorIPComp*allErr(2);
}
@@ -1274,17 +1431,18 @@ double TrainingDeepMaterialNetworkNonLinear::getEquivalentError(const fullVector
double TrainingDeepMaterialNetworkNonLinear::computeError(const LossPathMeasure& lossMeasure, const PathSTensor3& Fref, const PathSTensor3& Pref,
bool diff, fullVector<double>* grad,
bool extractPhaseDefo, PathSTensor3* phaseDefoPathRef,
- bool extractIPComp, PathDouble* IpCompPathRef)
+ bool extractIPComp, std::vector<PathDouble>* IpCompPathRef)
{
// evaluate path
+ int sizeIpCompPath= _phaseIpComForTraining.size();
static PathSTensor3 P, phaseDefo;
- static PathDouble IPcomp;
+ std::vector<PathDouble> IPcomp(sizeIpCompPath);
bool messageActive = false;
int lastConvergedIter = evaluatePath(Fref,P,messageActive,diff,extractPhaseDefo,&phaseDefo, extractIPComp,&IPcomp);
//Msg::Info("done evaluatePath lastConvergedIter = %d",lastConvergedIter);
static PathSTensor3 DerrorDP, DerrorDFi;
- static PathDouble DerrorDIPcomp;
+ std::vector<PathDouble> DerrorDIPcomp(sizeIpCompPath);
double error =lossMeasure.get(Pref,P,lastConvergedIter+1,diff,&DerrorDP);
//Msg::Info(" error first = %.16g",error);
if (extractPhaseDefo && fabs(_mixedFactor) > 0.)
@@ -1296,8 +1454,11 @@ double TrainingDeepMaterialNetworkNonLinear::computeError(const LossPathMeasure&
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- double nextError = lossMeasure.get(*IpCompPathRef,IPcomp,lastConvergedIter+1,diff,&DerrorDIPcomp);
- error += (_mixedFactorIPComp)*nextError;
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ double nextError = lossMeasure.get((*IpCompPathRef)[k],IPcomp[k],lastConvergedIter+1,diff,&DerrorDIPcomp[k]);
+ error += (_mixedFactorIPComp)*nextError/(double)sizeIpCompPath;
+ }
//Msg::Info(" error third = %.16g",nextError);
}
//Msg::Info(" error total = %.16g",error);
@@ -1333,7 +1494,10 @@ double TrainingDeepMaterialNetworkNonLinear::computeError(const LossPathMeasure&
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- val += (_mixedFactorIPComp)*DerrorDIPcomp.getValue(i)*IPcomp.getDValueDWeight(i)[itR];
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ val += (_mixedFactorIPComp)*DerrorDIPcomp[k].getValue(i)*IPcomp[k].getDValueDWeight(i)[itR]/(double)sizeIpCompPath;
+ }
}
assembleOneDof(key,val,*grad);
};
@@ -1350,7 +1514,10 @@ double TrainingDeepMaterialNetworkNonLinear::computeError(const LossPathMeasure&
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- val += (_mixedFactorIPComp)*DerrorDIPcomp.getValue(i)*IPcomp.getDValueDCoefficient(i)[itR];
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ val += (_mixedFactorIPComp)*DerrorDIPcomp[k].getValue(i)*IPcomp[k].getDValueDCoefficient(i)[itR]/(double)sizeIpCompPath;
+ }
}
assembleOneDof(key,val,*grad);
};
@@ -1366,7 +1533,10 @@ double TrainingDeepMaterialNetworkNonLinear::computeError(const LossPathMeasure&
}
if (extractIPComp && fabs(_mixedFactorIPComp) > 0.)
{
- val += (_mixedFactorIPComp)*DerrorDIPcomp.getValue(i)*IPcomp.getDValueDNormal(i)[itR];
+ for (int k=0; k< sizeIpCompPath; k++)
+ {
+ val += (_mixedFactorIPComp)*DerrorDIPcomp[k].getValue(i)*IPcomp[k].getDValueDNormal(i)[itR]/(double)sizeIpCompPath;
+ }
}
assembleOneDof(key,val,*grad);
};
@@ -1399,7 +1569,7 @@ double TrainingDeepMaterialNetworkNonLinear::computeCostFunction(const LossPathM
Msg::Error("sample %s does not exist !!!",sample);
Msg::Exit(0);
}
- if (_trainingWithPhaseDeformation && _phaseIpComForTraining != IPField::UNDEFINED)
+ if (_trainingWithPhaseDeformation && _phaseIpComForTraining.size() > 0)
{
err += computeError(loss,_XTrain[sample],_YTrain[sample],true,&gPart,true,&_YTrainPhaseDefo[sample], true, &_YTrainPhaseIPComp[sample]);
}
@@ -1407,7 +1577,7 @@ double TrainingDeepMaterialNetworkNonLinear::computeCostFunction(const LossPathM
{
err += computeError(loss,_XTrain[sample],_YTrain[sample],true,&gPart,true,&_YTrainPhaseDefo[sample]);
}
- else if (_phaseIpComForTraining != IPField::UNDEFINED)
+ else if (_phaseIpComForTraining.size() >0)
{
err += computeError(loss,_XTrain[sample],_YTrain[sample],true,&gPart,false,NULL, true, &_YTrainPhaseIPComp[sample]);
}
@@ -1429,7 +1599,7 @@ double TrainingDeepMaterialNetworkNonLinear::evaluateTrainingSet(const LossPathM
double err = 0;
for (int i=0; i< Ntrain; i++)
{
- if (_trainingWithPhaseDeformation && _phaseIpComForTraining != IPField::UNDEFINED)
+ if (_trainingWithPhaseDeformation && _phaseIpComForTraining.size() > 0)
{
err += computeError(loss,_XTrain[i],_YTrain[i],false,NULL,true,&_YTrainPhaseDefo[i],true,&_YTrainPhaseIPComp[i]);
}
@@ -1437,7 +1607,7 @@ double TrainingDeepMaterialNetworkNonLinear::evaluateTrainingSet(const LossPathM
{
err += computeError(loss,_XTrain[i],_YTrain[i],false,NULL,true,&_YTrainPhaseDefo[i]);
}
- else if (_phaseIpComForTraining != IPField::UNDEFINED)
+ else if (_phaseIpComForTraining.size() > 0)
{
err += computeError(loss,_XTrain[i],_YTrain[i],false,NULL,false,NULL,true,&_YTrainPhaseIPComp[i]);
}
@@ -1461,7 +1631,7 @@ double TrainingDeepMaterialNetworkNonLinear::evaluateTestSet(const LossPathMeasu
double err = 0;
for (int i=0; i< NTest; i++)
{
- if (_trainingWithPhaseDeformation && _phaseIpComForTraining != IPField::UNDEFINED)
+ if (_trainingWithPhaseDeformation && _phaseIpComForTraining.size() > 0)
{
err += computeError(loss,_XTest[i],_YTest[i],false,NULL,true,&_YTestPhaseDefo[i],true,&_YTestPhaseIPComp[i]);
}
@@ -1469,7 +1639,7 @@ double TrainingDeepMaterialNetworkNonLinear::evaluateTestSet(const LossPathMeasu
{
err += computeError(loss,_XTest[i],_YTest[i],false,NULL,true,&_YTestPhaseDefo[i]);
}
- else if (_phaseIpComForTraining != IPField::UNDEFINED)
+ else if (_phaseIpComForTraining.size() > 0)
{
err += computeError(loss,_XTest[i],_YTest[i],false,NULL,false,NULL,true,&_YTestPhaseIPComp[i]);
}
@@ -1497,21 +1667,48 @@ void TrainingDeepMaterialNetworkNonLinear::trainingDataSizeWithPhaseDeformation(
_YTrainPhaseDefo.resize(Ns,PathSTensor3());
//
_trainingWithPhaseDeformation = true;
- _phaseIndexForTraining = phaseIndex;
+ _phaseIndexForTrainingPhaseDeformation = phaseIndex;
_mixedFactor = mixFactor;
};
-void TrainingDeepMaterialNetworkNonLinear::trainingDataSizeWithPhaseIPComp(int Ns, int phaseIndex, int IPComp, double mixFactor)
+void TrainingDeepMaterialNetworkNonLinear::trainingDataSizeWithPhaseIPComp(int Ns, int phaseIndex, int IPComp, double mixFactor, double offset)
{
_XTrain.resize(Ns,PathSTensor3());
_YTrain.resize(Ns,PathSTensor3());
- _YTrainPhaseIPComp.resize(Ns,PathDouble());
+ _YTrainPhaseIPComp.clear();
+ _YTrainPhaseIPComp.resize(Ns, std::vector<PathDouble>(1, PathDouble()));
- _phaseIndexForTraining = phaseIndex;
- _phaseIpComForTraining = (IPField::Output)IPComp;
+ _phaseIndexForIpComForTraining.clear();
+ _phaseIndexForIpComForTraining.resize(1,phaseIndex);
+ _phaseIpComForTraining.clear();
+ _phaseIpComForTraining.resize(1,(IPField::Output)IPComp);
_mixedFactorIPComp = mixFactor;
+
+ _ipOffsets.clear();
+ _ipOffsets.resize(1,offset);
}
+void TrainingDeepMaterialNetworkNonLinear::trainingDataSizeWithTwoPhaseIPComps(int Ns, int phaseIndex, int IPComp0, int IPComp1, double mixFactor, double offset0, double offset1)
+{
+ _XTrain.resize(Ns,PathSTensor3());
+ _YTrain.resize(Ns,PathSTensor3());
+ _YTrainPhaseIPComp.clear();
+ _YTrainPhaseIPComp.resize(Ns, std::vector<PathDouble>(2,PathDouble()));
+
+
+ _phaseIndexForIpComForTraining.clear();
+ _phaseIndexForIpComForTraining.push_back(phaseIndex);
+ _phaseIndexForIpComForTraining.push_back(phaseIndex);
+ _phaseIpComForTraining.clear();
+ _phaseIpComForTraining.push_back((IPField::Output)IPComp0);
+ _phaseIpComForTraining.push_back((IPField::Output)IPComp1);
+ _mixedFactorIPComp = mixFactor;
+
+ _ipOffsets.clear();
+ _ipOffsets.push_back(offset0);
+ _ipOffsets.push_back(offset1);
+};
+
void TrainingDeepMaterialNetworkNonLinear::setTrainingSample(int row, const PathSTensor3& Fref, const PathSTensor3& Pref, int space)
{
if (row > _XTrain.size()-1)
@@ -1536,13 +1733,9 @@ void TrainingDeepMaterialNetworkNonLinear::setTrainingSamplePhaseDeformation(int
{
if (_trainingWithPhaseDeformation)
{
- if (_YTrainPhaseDefo.size() < _XTrain.size())
- {
- _YTrainPhaseDefo.resize(_XTrain.size(),PathSTensor3());
- }
- if (row > _XTrain.size()-1)
+ if (row > _YTrainPhaseDefo.size()-1)
{
- Msg::Error("row %d exceed the training data size %d",row,_XTrain.size());
+ Msg::Error("row %d exceed the training data size %d",row,_YTrainPhaseDefo.size());
Msg::Exit(0);
}
//
@@ -1557,28 +1750,32 @@ void TrainingDeepMaterialNetworkNonLinear::setTrainingSamplePhaseDeformation(int
}
};
-void TrainingDeepMaterialNetworkNonLinear::setTrainingSamplePhaseIPComp(int row, const PathDouble& IPCompRef, int space)
+void TrainingDeepMaterialNetworkNonLinear::setTrainingSamplePhaseIPComp(int row, const PathDouble& IPCompRef, int ipLoc, int space)
{
- if (_phaseIpComForTraining != IPField::UNDEFINED)
+ if (_phaseIpComForTraining.size() > 0)
{
- if (_YTrainPhaseIPComp.size() < _XTrain.size())
- {
- _YTrainPhaseIPComp.resize(_XTrain.size(),PathDouble());
- }
- if (row > _XTrain.size()-1)
+ if (row > _YTrainPhaseIPComp.size()-1)
{
- Msg::Error("row %d exceed the training data size %d",row,_XTrain.size());
+ Msg::Error("row %d exceed the training data size %d",row,_YTrainPhaseIPComp.size());
Msg::Exit(0);
}
//
if (space == 1)
{
- _YTrainPhaseIPComp[row] = IPCompRef;
+ _YTrainPhaseIPComp[row][ipLoc] = IPCompRef;
}
else
{
- IPCompRef.copyTo(_YTrainPhaseIPComp[row],space);
+ IPCompRef.copyTo(_YTrainPhaseIPComp[row][ipLoc],space);
}
+
+ Msg::Error("apply offset = %e to %d", _ipOffsets[ipLoc], IPField::ToString(_phaseIpComForTraining[ipLoc]));
+ _YTrainPhaseIPComp[row][ipLoc].addValue(_ipOffsets[ipLoc]);
+ }
+ else
+ {
+ Msg::Error("internal variables for training are empty");
+ Msg::Exit(0);
}
};
@@ -1590,34 +1787,38 @@ void TrainingDeepMaterialNetworkNonLinear::testDataSize(int Ns)
{
_YTestPhaseDefo.resize(Ns,PathSTensor3());
}
- if (_phaseIndexForTraining != IPField::UNDEFINED)
+ if (_phaseIpComForTraining.size() > 0)
{
- _YTestPhaseIPComp.resize(Ns,PathDouble());
+ _YTestPhaseIPComp.clear();
+ _YTestPhaseIPComp.resize(Ns, std::vector<PathDouble>(_phaseIpComForTraining.size(),PathDouble()));
}
};
-void TrainingDeepMaterialNetworkNonLinear::setTestSamplelePhaseIPComp(int row, const PathDouble& IPCompRef, int space)
+void TrainingDeepMaterialNetworkNonLinear::setTestSamplelePhaseIPComp(int row, const PathDouble& IPCompRef, int ipLoc, int space)
{
- if (_phaseIpComForTraining != IPField::UNDEFINED)
+ if (_phaseIpComForTraining.size() > 0)
{
- if (_YTestPhaseIPComp.size() < _XTest.size())
- {
- _YTestPhaseIPComp.resize(_XTest.size(),PathDouble());
- }
- if (row > _XTest.size()-1)
+ if (row > _YTestPhaseIPComp.size()-1)
{
- Msg::Error("row %d exceed the training data size %d",row,_XTest.size());
+ Msg::Error("row %d exceed the training data size %d",row,_YTestPhaseIPComp.size());
Msg::Exit(0);
}
//
if (space == 1)
{
- _YTestPhaseIPComp[row] = IPCompRef;
+ _YTestPhaseIPComp[row][ipLoc] = IPCompRef;
}
else
{
- IPCompRef.copyTo(_YTestPhaseIPComp[row],space);
+ IPCompRef.copyTo(_YTestPhaseIPComp[row][ipLoc],space);
}
+
+ _YTestPhaseIPComp[row][ipLoc].addValue(_ipOffsets[ipLoc]);
+ }
+ else
+ {
+ Msg::Error("internal variables for testing are empty");
+ Msg::Exit(0);
}
};
@@ -1625,13 +1826,9 @@ void TrainingDeepMaterialNetworkNonLinear::setTestSamplelePhaseDeformation(int r
{
if (_trainingWithPhaseDeformation)
{
- if (_YTestPhaseDefo.size() < _XTest.size())
+ if (row > _YTestPhaseDefo.size()-1)
{
- _YTestPhaseDefo.resize(_XTest.size(),PathSTensor3());
- }
- if (row > _XTest.size()-1)
- {
- Msg::Error("row %d exceed the testing data size %d",row,_XTest.size());
+ Msg::Error("row %d exceed the testing data size %d",row,_YTestPhaseDefo.size());
Msg::Exit(0);
}
if (space == 1)
@@ -1709,7 +1906,7 @@ void TrainingDeepMaterialNetworkNonLinear::evaluateTrainingSet(const std::vector
for (int i=0; i< Ntrain; i++)
{
static std::vector<fullVector<double> > allErrors;
- if (_trainingWithPhaseDeformation && _phaseIpComForTraining != IPField::UNDEFINED)
+ if (_trainingWithPhaseDeformation && _phaseIpComForTraining.size() >0)
{
computeErrorVec(allErrors,allMetrics,_XTrain[i],_YTrain[i],true,&_YTrainPhaseDefo[i],true,&_YTrainPhaseIPComp[i]);
}
@@ -1717,7 +1914,7 @@ void TrainingDeepMaterialNetworkNonLinear::evaluateTrainingSet(const std::vector
{
computeErrorVec(allErrors,allMetrics,_XTrain[i],_YTrain[i],true,&_YTrainPhaseDefo[i]);
}
- else if (_phaseIpComForTraining != IPField::UNDEFINED)
+ else if (_phaseIpComForTraining.size() > 0)
{
computeErrorVec(allErrors,allMetrics,_XTrain[i],_YTrain[i],false,NULL,true,&_YTrainPhaseIPComp[i]);
}
@@ -1772,7 +1969,7 @@ void TrainingDeepMaterialNetworkNonLinear::evaluateTestSet(const std::vector<std
for (int i=0; i< NTest; i++)
{
static std::vector<fullVector<double> > allErrors;
- if (_trainingWithPhaseDeformation && _phaseIpComForTraining != IPField::UNDEFINED)
+ if (_trainingWithPhaseDeformation && _phaseIpComForTraining.size() >0)
{
computeErrorVec(allErrors,allMetrics,_XTest[i],_YTest[i],true,&_YTestPhaseDefo[i], true, &_YTestPhaseIPComp[i]);
}
@@ -1780,7 +1977,7 @@ void TrainingDeepMaterialNetworkNonLinear::evaluateTestSet(const std::vector<std
{
computeErrorVec(allErrors,allMetrics,_XTest[i],_YTest[i],true,&_YTestPhaseDefo[i]);
}
- else if (_phaseIpComForTraining != IPField::UNDEFINED)
+ else if (_phaseIpComForTraining.size() >0)
{
computeErrorVec(allErrors,allMetrics,_XTest[i],_YTest[i],false,NULL, true, &_YTestPhaseIPComp[i]);
}
@@ -1827,4 +2024,13 @@ void TrainingDeepMaterialNetworkNonLinear::updateModelFromFittingParameters(cons
void TrainingDeepMaterialNetworkNonLinear::reinitializeModel()
{
_interaction->reInitialize();
-}
\ No newline at end of file
+};
+
+int evaluatePath(DeepMaterialNetwork* dmn, double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath)
+{
+ CoefficientReductionVoid reduction;
+ TrainingDeepMaterialNetworkNonLinear offTraining(*dmn,reduction);
+ offTraining.setMaximumStrainStep(maximumStrainStep);
+ int nb = offTraining.evaluatePath(strainPath,stressPath,true,false);
+ return nb;
+};
\ No newline at end of file
diff --git a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h
index 3bf011cd2019d89cda8e69f37434057c4aa7db13..5bf0de89e2d4fa74f102cff169f3b8bfeb057c25 100644
--- a/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h
+++ b/NonLinearSolver/modelReduction/TrainingDeepMaterialNetworksNonLinear.h
@@ -15,6 +15,7 @@
#include "TrainingArbitraryDeepMaterialNetworks.h"
#endif //SWIG
+class PathDouble;
class PathSTensor3
{
#ifndef SWIG
@@ -45,9 +46,11 @@ class PathSTensor3
int size() const;
void clear();
void setVal(int step, int i, int j, double val);
+ double getValue(int step, int i, int j) const;
void fromFile(const std::string filename);
void toFile(const std::string filename, int num) const;
void copyTo(PathSTensor3& dest, int space=1) const;
+ PathDouble determinant() const;
};
@@ -64,7 +67,6 @@ class PathDouble
public:
PathDouble& operator = (const PathDouble& src);
double& getValue(int step);
- const double& getValue(int step) const;
std::vector<double>& getDValueDWeight(int step);
const std::vector<double>& getDValueDWeight(int step) const;
std::vector<double>& getDValueDCoefficient(int step);
@@ -80,6 +82,8 @@ class PathDouble
void allocateWithDerivatives(int N, int sizeW, int sizeAlpha, int sizeNormal);
int size() const;
void clear();
+ void addValue(double v);
+ double getValue(int step) const;
void setVal(int step, double val);
void fromFile(const std::string filename);
void toFile(const std::string filename, int num) const;
@@ -136,61 +140,122 @@ class RelativeAbsoluteLossPathMeasure : public LossPathMeasure
};
+/*! \brief Training material networks with nonlinear paths
+*/
class TrainingDeepMaterialNetworkNonLinear : public TrainingArbitraryDeepMaterialNetwork
{
#ifndef SWIG
protected:
- DeepMaterialNetwork* _DMN; // estimator
+ DeepMaterialNetwork* _DMN; ///< deep material network as estimator
// maximal strain step
- double _maximumStrainStep;
- bool _trainingWithPhaseDeformation;
- int _phaseIndexForTraining; // phase index for training when using phase defo or internal state
- double _mixedFactor, _mixedFactorIPComp;
- IPField::Output _phaseIpComForTraining;
+ double _maximumStrainStep; ///< allowable strain step when using DMN
+ bool _trainingWithPhaseDeformation; ///< activate training with the per-phase average deformation
+ int _phaseIndexForTrainingPhaseDeformation;; ///< phase index for training when using the per-phase average deformation
+ double _mixedFactor, _mixedFactorIPComp; ///< factor to combine different norm
+ std::vector<int> _phaseIndexForIpComForTraining; ///< phase index for training when using an internal variable
+ std::vector<IPField::Output> _phaseIpComForTraining; ///< internal variable specification
+ std::vector<double> _ipOffsets; ///< offset value for internal variable specification
// offline data
- std::vector<PathSTensor3> _XTrain;
- std::vector<PathSTensor3> _YTrain;
- std::vector<PathSTensor3> _YTrainPhaseDefo;
- std::vector<PathDouble> _YTrainPhaseIPComp;
- std::vector<PathSTensor3>_XTest;
- std::vector<PathSTensor3> _YTest;
- std::vector<PathSTensor3> _YTestPhaseDefo;
- std::vector<PathDouble> _YTestPhaseIPComp;
+ std::vector<PathSTensor3> _XTrain; ///< strain path in the training dataset
+ std::vector<PathSTensor3> _YTrain; ///< stress path in the training dataset
+ std::vector<PathSTensor3> _YTrainPhaseDefo; ///< per-phase average strain path in the training dataset
+ std::vector<std::vector<PathDouble> >_YTrainPhaseIPComp; ///< internal variable path in the training dataset, row- sampleIndex
+ std::vector<PathSTensor3>_XTest; ///< strain path in the testing dataset
+ std::vector<PathSTensor3> _YTest; ///< stress path in the testing dataset
+ std::vector<PathSTensor3> _YTestPhaseDefo; ///< per-phase average strain path in testing dataset
+ std::vector<std::vector<PathDouble> >_YTestPhaseIPComp; ///< internal variable path in testing dataset , row- sampleIndex
protected:
void assembleOneDof(const Dof& key, double val, fullVector<double>& grad) const;
+
#endif //SWIG
public:
TrainingDeepMaterialNetworkNonLinear(DeepMaterialNetwork& dmn, const CoefficientReduction& reduction);
virtual ~TrainingDeepMaterialNetworkNonLinear();
+ /*! \brief set the maximum strain step for the DMN evaluation
+ @param[in] val maximum step value
+ */
void setMaximumStrainStep(double val);
- // return last sucessfuly iteration index at which the simulation converges
- int evaluatePath(const PathSTensor3& strainPath, PathSTensor3& stressPath, bool messageActive, bool stiff=false,
- bool extractPhaseDefo=false, PathSTensor3* phaseDefoPath=NULL,
- bool extractIPComp = false, PathDouble* IpCompPath=NULL);
- //
+
+ /*! \brief set number of samples in the training dataset
+ */
void trainingDataSize(int Ns);
+
+ /*! \brief activate the training with the per-phase average deformation
+ @param[in] Ns number of samples
+ @param[in] phaseIndex phase index to make the averaging operation
+ @param[in] mixFactor factor to combine the corresponding error norm
+ */
void trainingDataSizeWithPhaseDeformation(int Ns, int phaseIndex, double mixFactor=1.);
- void trainingDataSizeWithPhaseIPComp(int Ns, int phaseIndex, int IPComp, double mixFactor=1.);
- //
+
+ /*! \brief activate the training with a single internal variable
+ @param[in] Ns number of samples
+ @param[in] phaseIndex phase index to make the averaging operation
+ @param[in] IPComp an internal variable defined by IPField::Output
+ @param[in] mixFactor factor to combine the corresponding error norm
+ */
+ void trainingDataSizeWithPhaseIPComp(int Ns, int phaseIndex, int IPComp, double mixFactor=1., double offset=0);
+
+ /*! \brief activate the training with two internal variables
+ @param[in] Ns number of samples
+ @param[in] phaseIndex phase index to make the averaging operation
+ @param[in] IPComp0 an internal variable defined by IPField::Output
+ @param[in] IPComp1 an internal variable defined by IPField::Output
+ @param[in] mixFactor factor to combine the corresponding error norm
+ */
+ void trainingDataSizeWithTwoPhaseIPComps(int Ns, int phaseIndex, int IPComp0, int IPComp1, double mixFactor=1., double offset0=0, double offset1=0);
+
+ /*! \brief set training sample
+ @param[in] row sample index
+ @param[in] Fref strain path
+ @param[in] Pref stress path
+ @param[in] space slicing based on position on the path, values are taken after each "space" values
+ */
void setTrainingSample(int row, const PathSTensor3& Fref, const PathSTensor3& Pref, int space = 1);
+
+ /*! \brief set training sample for the per-phase average path
+ @param[in] row sample index
+ @param[in] FphaseRef per-phase average strain
+ @param[in] space slicing based on position on the path, values are taken after each "space" values
+ */
void setTrainingSamplePhaseDeformation(int row, const PathSTensor3& FphaseRef, int space = 1);
- void setTrainingSamplePhaseIPComp(int row, const PathDouble& IPCompRef, int space = 1);
+
+ /*! \brief set training sample for the internal variable
+ @param[in] row sample index
+ @param[in] IPCompRef internal variable path
+ @param[in] ipLoc location of IP variable in _phaseIpComForTraining
+ @param[in] space slicing based on position on the path, values are taken after each "space" values
+ */
+ void setTrainingSamplePhaseIPComp(int row, const PathDouble& IPCompRef, int ipLoc=0, int space = 1);
//
void testDataSize(int Ns);
void setTestSample(int row, const PathSTensor3& Fref, const PathSTensor3& Pref, int space = 1);
void setTestSamplelePhaseDeformation(int row, const PathSTensor3& FphaseRef, int space = 1);
- void setTestSamplelePhaseIPComp(int row, const PathDouble& IPCompRef, int space = 1);
+ void setTestSamplelePhaseIPComp(int row, const PathDouble& IPCompRef, int ipLoc=0, int space = 1);
//
#ifndef SWIG
+ /*! \brief evaluate DMN response based on a strain path
+ @param[in] strainPath strain path as input
+ @param[out] stressPath stress path as output
+ @param[in] messageActive true if the convergence history is printed in the terminal
+ @param[in] stiff true if the the tangent operator is estimated in the DMN
+ @param[in] extractPhaseDefo true if the per-phase average deformation is estimated
+ @param[out] phaseDefoPath the extracted per-phase average deformation path when extractPhaseDefo = true
+ @param[in] extractIPComp true if the internal variable is extracted
+ @param[out] IpCompPath the extracted internal variable path when extractIPComp = true
+ @return number of the last sucessfuly iteration index at which the simulation converges
+ */
+ int evaluatePath(const PathSTensor3& strainPath, PathSTensor3& stressPath, bool messageActive, bool stiff=false,
+ bool extractPhaseDefo=false, PathSTensor3* phaseDefoPath=NULL,
+ bool extractIPComp = false, std::vector<PathDouble>* IpCompPath=NULL);
void computeErrorVec(std::vector<fullVector<double> >& allErrors, const std::vector<LossPathMeasure*>& lossMeasures, const PathSTensor3& Fref, const PathSTensor3& Pref,
bool extractPhaseDefo=false, PathSTensor3* phaseDefoPathRef=NULL,
- bool extractIPComp = false, PathDouble* IpCompPathRef=NULL);
+ bool extractIPComp = false, std::vector<PathDouble>* IpCompPathRef=NULL);
double getEquivalentError(const fullVector<double>& allErr) const;
double computeError(const LossPathMeasure& lossMeasure, const PathSTensor3& Fref, const PathSTensor3& Pref,
bool diff=false, fullVector<double>* grad=NULL,
bool extractPhaseDefo=false, PathSTensor3* phaseDefoPathRef=NULL,
- bool extractIPComp = false, PathDouble* IpCompPathRef=NULL);
+ bool extractIPComp = false, std::vector<PathDouble>* IpCompPathRef=NULL);
double computeCostFunction(const LossPathMeasure& loss, const std::vector<int>& sampleIndex, fullVector<double>& grad);
double evaluateTrainingSet(const LossPathMeasure& loss);
@@ -217,4 +282,7 @@ class TrainingDeepMaterialNetworkNonLinear : public TrainingArbitraryDeepMateria
#endif //SWIG
};
+
+int evaluatePath(DeepMaterialNetwork* dmn, double maximumStrainStep, const PathSTensor3& strainPath, PathSTensor3& stressPath);
+
#endif //_TRAININGDEEPMATERIALNETWORKSNONLINEAR_H_
\ No newline at end of file
diff --git a/NonLinearSolver/nlmechsolpy.i b/NonLinearSolver/nlmechsolpy.i
index 30a05474aec3527835127ee5e25f11a20e50d45d..99fca14ad7f8e591b61bd716a4917a90fcd0bb26 100644
--- a/NonLinearSolver/nlmechsolpy.i
+++ b/NonLinearSolver/nlmechsolpy.i
@@ -74,6 +74,7 @@
#include "TrainingDeepMaterialNetworksNonLinear.h"
#include "laplaceSolver.h"
#include "JIntegralByDomainIntegration.h"
+ #include "highOrderTensor.h"
%}
%init%{
@@ -223,5 +224,6 @@
%include "TrainingDeepMaterialNetworksNonLinear.h"
%include "laplaceSolver.h"
%include "JIntegralByDomainIntegration.h"
+%include "highOrderTensor.h"
using namespace fullMatrixOperation;
diff --git a/NonLinearSolver/nlsolver/meshModification.cpp b/NonLinearSolver/nlsolver/meshModification.cpp
index 9542e437606282ad56dce29fd7ecdd400892e1f1..70f67d64b0a570a42e5e951f637382077a827fe2 100644
--- a/NonLinearSolver/nlsolver/meshModification.cpp
+++ b/NonLinearSolver/nlsolver/meshModification.cpp
@@ -173,6 +173,22 @@ bool makePhysicalByBox::box::inBox(MElement* ele) const
}
};
+void makePhysicalByBox::physicalToMesh(int dim, int phys, const std::string inputMeshFile, const std::string outputMeshFile, int newPhy)
+{
+ GModel pModel;
+ pModel.readMSH(inputMeshFile.c_str());
+ std::map<int, elementGroup> allGrps;
+ if (newPhy > 0)
+ {
+ allGrps[newPhy-1].addPhysical(dim,phys);
+ }
+ else
+ {
+ allGrps[phys-1].addPhysical(dim,phys);
+ }
+ write_MSH2(allGrps,outputMeshFile);
+};
+
void makePhysicalByBox::run(const std::string inputMeshFile, const std::string outputMeshFile,
double xmin, double xmax, int nx,
double ymin, double ymax, int ny)
diff --git a/NonLinearSolver/nlsolver/meshModification.h b/NonLinearSolver/nlsolver/meshModification.h
index d5325e3be24027128e6ec143a559f5efeaa40ad8..efda7aaefabff8a47309366dba875cc9b13e3455 100644
--- a/NonLinearSolver/nlsolver/meshModification.h
+++ b/NonLinearSolver/nlsolver/meshModification.h
@@ -111,6 +111,7 @@ class makePhysicalByBox
makePhysicalByBox(){};
~makePhysicalByBox(){};
+ void physicalToMesh(int dim, int phys, const std::string inputMeshFile, const std::string outputMeshFile, int newPhy=-1);
void run(const std::string inputMeshFile, const std::string outputMeshFile,
double xmin, double xmax, int nx,
double ymin, double ymax, int ny);
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
index 3092be1790350481ae2158442b58c13dc994de35..80e9c9ea1ff34bf1c6b1d7e246dc2a8a9c52d65e 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.cpp
@@ -190,9 +190,9 @@ void AvergageCluster::getAverageValue(const IPField* ipf, std::vector<double>& v
}
};
-TimeManager::TimeManager():
+TimeManager::TimeManager(double tol):
_timeSeries(),_numStepSeries(),
-_startTime(0.), _endTime(0.), _numStep(1), _tol(1e-6),
+_startTime(0.), _endTime(0.), _numStep(1), _tol(tol),
_lastTime(0.),
_lastIterationIndex(0.),
_timeStep(0.),
@@ -394,6 +394,7 @@ void TimeManager::computeTimeStepForNextSolving(int niteNR)
}
}
+
if (lastTimeInterVal == -1)
{
if (_timeStepAdaptation)
@@ -508,7 +509,7 @@ void TimeManager::initializeTimeSteppingPlan()
}
};
-PointwiseTimeManager::PointwiseTimeManager(double t):TimeManager(),_savePoints(),_tol(t){}
+PointwiseTimeManager::PointwiseTimeManager(double tol, double tolSave):TimeManager(tol),_savePoints(),_tolSave(tolSave){}
void PointwiseTimeManager::put(double val)
{
@@ -520,7 +521,7 @@ bool PointwiseTimeManager::willArchive(double curtime, int numstep) const
// the most stupid algorithm
for (int i=0; i< _savePoints.size(); i++)
{
- if (fabs(curtime - _savePoints[i]) < _tol)
+ if (fabs(curtime - _savePoints[i]) < _tolSave)
{
printf("saving file: \n");
return true;
@@ -529,7 +530,7 @@ bool PointwiseTimeManager::willArchive(double curtime, int numstep) const
return false;
}
-StepwiseBeginTimeManager::StepwiseBeginTimeManager(int nstepArchBegin):TimeManager(),_nstepArchBegin(nstepArchBegin){}
+StepwiseBeginTimeManager::StepwiseBeginTimeManager(double tol,int nstepArchBegin):TimeManager(tol),_nstepArchBegin(nstepArchBegin){}
bool StepwiseBeginTimeManager::willArchive(double curtime, int numstep) const
{
@@ -3150,10 +3151,7 @@ void nonLinearMechSolver::init2()
Msg::Info("Fix and number Dofs");
this->numberDofs();
- printf("total number of unknown Dofs = %d\n",pAssembler->sizeOfR());
- printf("total number of fixed Dofs = %d\n",pAssembler->sizeOfF());
- Msg::Info("Dofs are fixed and numbered");
-
+
/* MPI init */
#if defined(HAVE_MPI)
if(whatScheme!=StaticLinear and Msg::GetCommSize()>1)
@@ -3163,6 +3161,10 @@ void nonLinearMechSolver::init2()
Msg::Info("MPI initialization OK");
}
#endif // HAVE_MPI
+
+ printf("total number of unknown Dofs = %d\n",pAssembler->sizeOfR());
+ printf("total number of fixed Dofs = %d\n",pAssembler->sizeOfF());
+ Msg::Info("Dofs are fixed and numbered");
// make data as previous if exists
this->updateDataFromPreviousSolver();
@@ -3203,7 +3205,7 @@ void nonLinearMechSolver::init2()
Msg::Info("Creation of IPVariable");
_ipf = new IPField(this,vaip,ipView,ipViewInterface); // Field for GaussPoint
Msg::Info("IPVariable are created");
-
+
_ipf->compute1state(IPStateBase::initial,stiffcomputation);
_ipf->initialBroken(pModel, initbrokeninter);
_ipf->initialBrokenDomain(pModel,initbrokeninterInDomains);
@@ -5697,6 +5699,11 @@ void nonLinearMechSolver::createInterfaceElement(){
#if defined(HAVE_MPI)
+ #if !defined(HAVE_CG_MPI_INTERFACE)
+ if (isDgDomain())
+ {
+ #endif //
+
if(Msg::GetCommSize() > 1 and (_mapRanks.find(Msg::GetCommRank()) != _mapRanks.end())){
// first add interface to the interdomain defined by the user
for(std::vector<partDomain*>::iterator itdom = domainVector.begin(); itdom != domainVector.end(); ++itdom){
@@ -5774,6 +5781,9 @@ void nonLinearMechSolver::createInterfaceElement(){
}
}
}
+ #if !defined(HAVE_CG_MPI_INTERFACE)
+ }
+ #endif //
#endif // HAVE_MPI
// save to file
@@ -10218,23 +10228,60 @@ void nonLinearMechSolver::createSystem()
if(whatScheme==Explicit)
{
lsys = this->createExplicitSystem();
- pAssembler = new staticDofManager(lsys,whatScheme,false,true);
+ #if !defined(HAVE_CG_MPI_INTERFACE)
+ if (!isDgDomain())
+ {
+ pAssembler = new staticDofManagerFullCG(lsys,whatScheme,false,true);
+ }
+ else
+ #endif //
+ {
+ pAssembler = new staticDofManager(lsys,whatScheme,false,true);
+ }
}
else if (whatScheme == Implicit)
{
lsys = this->createImplicitSystem();
- pAssembler = new staticDofManager(lsys,whatScheme,true,true);
+ #if !defined(HAVE_CG_MPI_INTERFACE)
+ if (!isDgDomain())
+ {
+ pAssembler = new staticDofManagerFullCG(lsys,whatScheme,true,true);
+ }
+ else
+ #endif //
+ {
+ pAssembler = new staticDofManager(lsys,whatScheme,true,true);
+ }
}
else if (whatScheme == StaticLinear or whatScheme == StaticNonLinear)
{
lsys = this->createSNLSystem();
if (whatScheme == StaticLinear)
{
- pAssembler = new staticDofManager(lsys,whatScheme,true,true);
+ #if !defined(HAVE_CG_MPI_INTERFACE)
+ if (!isDgDomain())
+ {
+ pAssembler = new staticDofManagerFullCG(lsys,whatScheme,true,true);
+ }
+ else
+ #endif
+ {
+ pAssembler = new staticDofManager(lsys,whatScheme,true,true);
+ }
}
else
{
- pAssembler = new staticDofManager(lsys,whatScheme,true,true);
+ #if !defined(HAVE_CG_MPI_INTERFACE)
+ if (!isDgDomain())
+ {
+ pAssembler = new staticDofManagerFullCG(lsys,whatScheme,true,true);
+ }
+ else
+ #endif
+ {
+ pAssembler = new staticDofManager(lsys,whatScheme,true,true);
+
+ }
}
}
else if (whatScheme == Eigen)
diff --git a/NonLinearSolver/nlsolver/nonLinearMechSolver.h b/NonLinearSolver/nlsolver/nonLinearMechSolver.h
index f26dc08be61da4387112e99d6d2d80b3d891ff4a..2472a71678b79c870d7d8487d34b719514be2098 100644
--- a/NonLinearSolver/nlsolver/nonLinearMechSolver.h
+++ b/NonLinearSolver/nlsolver/nonLinearMechSolver.h
@@ -494,7 +494,7 @@ class TimeManager
friend class pathFollowingManager;
#endif //SWIG
public:
- TimeManager();
+ TimeManager(double tol=1e-12);
void reset();
void activateTimeStepAdaptation(bool adap, int numNROptimal, double exp, double maximalStep, double minTimeStep);
void setManageTimeStepSafeguard(const int miteNR,const int iteIncreaseTS, const double redfactor, const int maxAttemptRedFactor);
@@ -539,11 +539,11 @@ class PointwiseTimeManager : public TimeManager
{
protected:
std::vector<double> _savePoints;
- double _tol;
+ double _tolSave;
public:
- PointwiseTimeManager(double tol = 1e-6);
+ PointwiseTimeManager(double tol = 1e-12, double tolSave=1e-12);
void put(double val);
- PointwiseTimeManager(const PointwiseTimeManager& src):TimeManager(src),_savePoints(src._savePoints),_tol(src._tol){}
+ PointwiseTimeManager(const PointwiseTimeManager& src):TimeManager(src),_savePoints(src._savePoints),_tolSave(src._tolSave){}
virtual ~PointwiseTimeManager(){}
virtual bool willArchive(double curtime, int numstep) const;
virtual TimeManager* clone() const {return new PointwiseTimeManager(*this);};
@@ -554,7 +554,7 @@ class StepwiseBeginTimeManager : public TimeManager
protected:
int _nstepArchBegin;
public:
- StepwiseBeginTimeManager(int nstepArchBegin = 1);
+ StepwiseBeginTimeManager(double tol = 1e-12, int nstepArchBegin = 1);
StepwiseBeginTimeManager(const StepwiseBeginTimeManager& src):TimeManager(src),_nstepArchBegin(src._nstepArchBegin){}
virtual ~StepwiseBeginTimeManager(){}
virtual bool willArchive(double curtime, int numstep) const;
diff --git a/NonLinearSolver/nlsolver/staticDofManager.cpp b/NonLinearSolver/nlsolver/staticDofManager.cpp
index 896f5776c7ccdd9f079d0b148233a2bd9f923ec5..3dc875cbe2120f8166ccc66480c05e7a7d8e9d12 100644
--- a/NonLinearSolver/nlsolver/staticDofManager.cpp
+++ b/NonLinearSolver/nlsolver/staticDofManager.cpp
@@ -2247,5 +2247,370 @@ void staticDofManager::linearElastic_assemble(const std::vector<Dof> &R, const f
}
}
}
+};
+
+
+#if defined(HAVE_MPI)
+int staticDofManagerFullCG::sizeOfR() const
+{
+ if (Msg::GetCommSize() > 1 and _isPartitioned)
+ {
+ return _localNumberDofMPI[Msg::GetCommRank()];
+ }
+ else
+ {
+ return staticDofManager::sizeOfR();
+ }
+};
+
+
+void staticDofManagerFullCG::initMPI()
+{
+ if (!_isPartitioned)
+ {
+ return;
+ }
+
+ if(Msg::GetCommSize() > 1)
+ {
+ MPI_Status status;
+ _otherPartUnknowns.clear();
+ _localDofOtherPartUnknowns.clear();
+ if (Msg::GetCommRank() > 0)
+ {
+ // set data to other rank
+ int localUnknownSize = unknown.size();
+ std::vector<int> buffer(2*localUnknownSize);
+ int j=0;
+ for (std::map<Dof, int>::const_iterator it = unknown.begin(); it != unknown.end(); it++)
+ {
+ const Dof& D = it->first;
+ buffer[j] = D.getEntity();
+ buffer[j+1]=D.getType();
+ j += 2;
+ }
+ MPI_Send(&localUnknownSize, 1, MPI_INT, 0, Msg::GetCommSize()+ Msg::GetCommRank(), MPI_COMM_WORLD);
+ MPI_Send(&buffer[0], 2*localUnknownSize, MPI_INT, 0, Msg::GetCommRank(), MPI_COMM_WORLD);
+ }
+ else
+ {
+ std::map<int, std::map<Dof, int> > unknownRanksTrue;
+ std::map<int, std::map<Dof, int> > unknownRanks;
+ std::map<int, std::set<Dof> > allDofToSendOtherRanks;
+ std::map<Dof, int> globalDofNumbers;
+ globalDofNumbers = unknown;
+ unknownRanks[0] = unknown;
+ unknownRanksTrue[0]=unknown;
+ _localNumberDofMPI[0] = unknown.size();
+
+ for (int i=1; i < Msg::GetCommSize(); i++)
+ {
+ int localUnknownSizeOtherRank = 0;
+ MPI_Recv(&localUnknownSizeOtherRank,1,MPI_INT,i,Msg::GetCommSize()+i,MPI_COMM_WORLD,&status);
+ std::vector<int> buffer(2*localUnknownSizeOtherRank);
+ MPI_Recv(&buffer[0],2*localUnknownSizeOtherRank,MPI_INT,i,i,MPI_COMM_WORLD,&status);
+ for (int j=0; j< localUnknownSizeOtherRank; j++)
+ {
+ Dof key(buffer[2*j],buffer[2*j+1]);
+ if (globalDofNumbers.find(key) == globalDofNumbers.end())
+ {
+ int oldSize = globalDofNumbers.size();
+ globalDofNumbers[key] = oldSize;
+ unknownRanksTrue[i][key] = oldSize;
+ }
+ unknownRanks[i][key] = globalDofNumbers[key];
+ bool inOtherRank=false;
+ for (int k=0; k< i; k++)
+ {
+ // if key is present in lower rank
+ if (unknownRanksTrue[k].find(key) != unknownRanksTrue[k].end())
+ {
+ allDofToSendOtherRanks[k].insert(key);
+ inOtherRank=true;
+ }
+ }
+ if (inOtherRank)
+ {
+ if (_otherPartUnknowns.find(key) == _otherPartUnknowns.end())
+ {
+ int oldSize = _otherPartUnknowns.size();
+ _otherPartUnknowns[key]= oldSize;
+ }
+ }
+ }
+ _localNumberDofMPI[i] = globalDofNumbers.size();
+ for (int j=0; j< i; j++)
+ {
+ _localNumberDofMPI[i] -= _localNumberDofMPI[j];
+ }
+ }
+
+ // send to other rank
+ for (int i=1; i < Msg::GetCommSize(); i++)
+ {
+ int localUnknownSizeOtherRank = unknownRanks[i].size();
+ std::vector<int> buffer(4*localUnknownSizeOtherRank);
+ int j=0;
+ for (std::map<Dof, int>::const_iterator it = unknownRanks[i].begin(); it != unknownRanks[i].end(); it++)
+ {
+ buffer[j] = it->first.getEntity();
+ buffer[j+1] = it->first.getType();
+ buffer[j+2] = it->second;
+ if (allDofToSendOtherRanks[i].find(it->first) == allDofToSendOtherRanks[i].end())
+ {
+ buffer[j+3] =0;
+ }
+ else
+ {
+ buffer[j+3] =1;
+ }
+ j += 4;
+ }
+ MPI_Send(&buffer[0], 4*localUnknownSizeOtherRank, MPI_INT, i, i, MPI_COMM_WORLD);
+ }
+ _localDofOtherPartUnknowns = allDofToSendOtherRanks[0];
+ }
+
+ if (Msg::GetCommRank() > 0)
+ {
+ int localUnknownSizeOtherRank = unknown.size();
+ std::vector<int> buffer(4*localUnknownSizeOtherRank);
+ MPI_Recv(&buffer[0],4*localUnknownSizeOtherRank,MPI_INT,0, Msg::GetCommRank(),MPI_COMM_WORLD,&status);
+ for (int j=0; j< localUnknownSizeOtherRank; j++)
+ {
+ Dof key(buffer[4*j],buffer[4*j+1]);
+ unknown[key] = buffer[4*j+2];
+ if (buffer[4*j+3] == 1)
+ {
+ _localDofOtherPartUnknowns.insert(key);
+ }
+ }
+ }
+
+ MPI_Bcast(_localNumberDofMPI, Msg::GetCommSize(), MPI_INT, 0, MPI_COMM_WORLD);
+ int numDofOtherRanks = 0;
+ if (Msg::GetCommRank() ==0)
+ {
+ numDofOtherRanks = _otherPartUnknowns.size();
+ }
+ MPI_Bcast(&numDofOtherRanks, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ std::vector<int> otherPartUnknownsBuffer(3*numDofOtherRanks);
+ if (Msg::GetCommRank() == 0)
+ {
+ int j=0;
+ for (std::map<Dof, int>::const_iterator it = _otherPartUnknowns.begin(); it != _otherPartUnknowns.end(); it++)
+ {
+ otherPartUnknownsBuffer[j] = it->first.getEntity();
+ otherPartUnknownsBuffer[j+1] = it->first.getType();
+ otherPartUnknownsBuffer[j+2] = it->second;
+ j += 3;
+ }
+ }
+ MPI_Bcast(&otherPartUnknownsBuffer[0], 3*numDofOtherRanks, MPI_INT, 0, MPI_COMM_WORLD);
+ if (Msg::GetCommRank() > 0)
+ {
+ for (int j=0; j< numDofOtherRanks; j++)
+ {
+ Dof key(otherPartUnknownsBuffer[3*j],otherPartUnknownsBuffer[3*j+1]);
+ _otherPartUnknowns[key] = otherPartUnknownsBuffer[3*j+2];
+ }
+ }
+
+
+ int nbdofOtherPart = _otherPartUnknowns.size();
+ this->_otherPartPositions = new double[nbdofOtherPart];
+ this->_previousOtherPartPositions = new double[nbdofOtherPart];
+ for (int pp=0; pp < nbdofOtherPart; pp++)
+ {
+ this->_otherPartPositions[pp] = 0;
+ this->_previousOtherPartPositions[pp] = 0;
+ }
+
+ if(_scheme==nonLinearMechSolver::Explicit or _scheme == nonLinearMechSolver::Implicit)
+ {
+ this->_otherPartVelocities = new double[nbdofOtherPart];
+ this->_previousOtherPartVelocities = new double[nbdofOtherPart];
+ this->_otherPartAccelerations = new double[nbdofOtherPart];
+ this->_previousOtherPartAccelerations = new double[nbdofOtherPart];
+ for (int pp=0; pp < nbdofOtherPart; pp++)
+ {
+ _otherPartVelocities[pp] = 0.;
+ _previousOtherPartVelocities[pp] = 0.;
+ _otherPartAccelerations[pp] = 0.;
+ _previousOtherPartAccelerations[pp] = 0.;
+ }
+ }
+ }
+ _mpiInit = true;
+}
+
+void staticDofManagerFullCG::allocateSystem()
+{
+ if(!_mpiInit) this->initMPI();
+ staticDofManager::allocateSystem();
+}
+
+void staticDofManagerFullCG::setFirstRigidContactUnknowns()
+{
+ if (Msg::GetCommSize() > 1 and _isPartitioned)
+ {
+ if(!_mpiInit) this->initMPI();
+ _firstRigidContactDof = _localNumberDofMPI[Msg::GetCommRank()];
+ }
+ else
+ {
+ staticDofManager::setFirstRigidContactUnknowns();
+ }
+}
+
+void staticDofManagerFullCG::manageMPIComm(const int otherPartNum,const std::vector<Dof> &otherPartR)
+{
+ // do nothing
}
+void staticDofManagerFullCG::manageMPISparsity(const std::vector<Dof> &myR, const int otherPartNum,const std::vector<Dof> &otherPartR)
+{
+ // do nothing
+}
+
+void staticDofManagerFullCG::systemMPIComm()
+{
+ if (!_isPartitioned) {
+ return;
+ }
+ if(!_mpiInit) this->initMPI();
+
+ int rank = Msg::GetCommRank();
+ int numDofs = _localDofOtherPartUnknowns.size();
+ int nbdofPerNode =1;
+ if(_scheme == nonLinearMechSolver::Explicit or _scheme == nonLinearMechSolver::Implicit)
+ {
+ nbdofPerNode+=2;
+ }
+ int mpiBufferSize =_otherPartUnknowns.size() * nbdofPerNode;
+ std::vector<double> mpiBuffer(mpiBufferSize);
+ std::fill(mpiBuffer.begin(), mpiBuffer.end(),0.);
+ for(std::set<Dof>::const_iterator itR = _localDofOtherPartUnknowns.begin(); itR!=_localDofOtherPartUnknowns.end(); itR++)
+ {
+ int loc = _otherPartUnknowns[*itR];
+ int NR = unknown[*itR];
+ if(_scheme == nonLinearMechSolver::Explicit or _scheme == nonLinearMechSolver::Implicit)
+ {
+ _NLScurrent->getFromSolution(NR, mpiBuffer[3*loc],nonLinearBoundaryCondition::position);
+ _NLScurrent->getFromSolution(NR, mpiBuffer[3*loc+1],nonLinearBoundaryCondition::velocity);
+ _NLScurrent->getFromSolution(NR, mpiBuffer[3*loc+2],nonLinearBoundaryCondition::acceleration);
+ }
+ else
+ {
+ (this->_current)->getFromSolution(NR, mpiBuffer[loc]);
+ }
+ }
+ MPI_Allreduce(MPI_IN_PLACE,&mpiBuffer[0],mpiBufferSize,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ for (int j=0; j< _otherPartUnknowns.size(); j++)
+ {
+ _otherPartPositions[j] = mpiBuffer[nbdofPerNode*j];
+ if(_scheme == nonLinearMechSolver::Explicit or _scheme == nonLinearMechSolver::Implicit)
+ {
+ _otherPartVelocities[j]= mpiBuffer[nbdofPerNode*j+1];
+ _otherPartAccelerations[j] = mpiBuffer[nbdofPerNode*j+2];
+ }
+ }
+};
+
+void staticDofManagerFullCG::getDofValue(const Dof& key, double &val, nonLinearBoundaryCondition::whichCondition wv) const
+{
+ switch(wv)
+ {
+ case nonLinearBoundaryCondition::position:
+ {
+ // look first for fixed dof (the type can be negative)
+ std::map<Dof, double>::const_iterator it = this->fixed.find(key);
+ if (it != this->fixed.end())
+ {
+ val = it->second;
+ return ;
+ }
+ std::map<Dof,int>::const_iterator itR = this->_otherPartUnknowns.find(key);
+ if(itR != _otherPartUnknowns.end())
+ {
+ val = this->_otherPartPositions[itR->second];
+ return;
+ }
+ std::map<Dof, int>::const_iterator itK = this->unknown.find(key);
+ if (itK != this->unknown.end())
+ {
+ _current->getFromSolution(itK->second, val);
+ return;
+ }
+ std::map<Dof, DofAffineConstraint< double > >::const_iterator itC = this->constraints.find(key);
+ if (itC != this->constraints.end())
+ {
+ double tmp(val);
+ val = itC->second.shift;
+ for (unsigned i=0;i<(itC->second).linear.size();i++)
+ {
+ std::map<Dof, int>::const_iterator itu = this->unknown.find(((itC->second).linear[i]).first);
+ getDofValue(((itC->second).linear[i]).first, tmp, wv);
+ dofTraits<double>::gemm(val,((itC->second).linear[i]).second, tmp, 1, 1);
+ }
+ return ;
+ }
+ Msg::Error("position Dof not found %d %d",key.getEntity(), key.getType());
+ Msg::Exit(0);
+ }
+ break;
+ case nonLinearBoundaryCondition::velocity:
+ {
+ // look first for fixed Dof
+ std::map<Dof, FixedNodeDynamicData<double> >::const_iterator it = this->_fixedDynamic.find(key);
+ if (it != this->_fixedDynamic.end())
+ {
+ val = it->second.velocity;
+ return ;
+ }
+ std::map<Dof,int>::const_iterator itR = this->_otherPartUnknowns.find(key);
+ if (itR != _otherPartUnknowns.end())
+ {
+ val = this->_otherPartVelocities[itR->second];
+ return;
+ }
+ std::map<Dof, int>::const_iterator itK = this->unknown.find(key);
+ if (itK != this->unknown.end())
+ {
+ _NLScurrent->getFromSolution(itK->second, val,wv);
+ return;
+ }
+ Msg::Error("velocity Dof not found %d %d",key.getEntity(), key.getType());
+ Msg::Exit(0);
+ }
+ break;
+ case nonLinearBoundaryCondition::acceleration:
+ {
+ // look first for fixed Dof
+ std::map<Dof, FixedNodeDynamicData<double> >::const_iterator it = this->_fixedDynamic.find(key);
+ if (it != this->_fixedDynamic.end())
+ {
+ val = it->second.acceleration;
+ return ;
+ }
+ std::map<Dof,int>::const_iterator itR = this->_otherPartUnknowns.find(key);
+ if(itR != _otherPartUnknowns.end())
+ {
+ val = this->_otherPartAccelerations[itR->second];
+ return;
+ }
+ std::map<Dof, int>::const_iterator itK = this->unknown.find(key);
+ if (itK != this->unknown.end())
+ {
+ _NLScurrent->getFromSolution(itK->second, val,wv);
+ return;
+ }
+ Msg::Error("velocity Dof not found %d %d",key.getEntity(), key.getType());
+ Msg::Exit(0);
+ }
+ break;
+ };
+};
+
+#endif // HAVE_MPI
+
#endif // HAVE_MPI
diff --git a/NonLinearSolver/nlsolver/staticDofManager.h b/NonLinearSolver/nlsolver/staticDofManager.h
index e08c2b682f51add91b4b2b85f4569c56997a5041..c884f6579227d4312924c37a5c40f8f8b809bc2d 100644
--- a/NonLinearSolver/nlsolver/staticDofManager.h
+++ b/NonLinearSolver/nlsolver/staticDofManager.h
@@ -217,9 +217,9 @@ class staticDofManager : public nlsDofManager
// initialization of MPI communication (send the unknowns that this processor need and recieve from other rank the unknown to send
- private:
+ protected:
#if defined(HAVE_MPI)
- void initMPI();
+ virtual void initMPI();
#endif // HAVE_MPI
public:
@@ -347,5 +347,34 @@ class staticDofManager : public nlsDofManager
void linearElastic_assemble(const Dof &R, const Dof& C, const double &value);
void linearElastic_assemble(const std::vector<Dof> &R, const fullMatrix<double> &m);
};
+
+
+class staticDofManagerFullCG : public staticDofManager
+{
+ protected:
+
+ #if defined(HAVE_MPI)
+ // staticDofManager::_otherPartUnknowns all dofs at interfaces
+ std::set<Dof> _localDofOtherPartUnknowns; // local dofs at currant rank needs to be send to other ranks
+ virtual void initMPI();
+ #endif // HAVE_MPI
+
+ public:
+ staticDofManagerFullCG(linearSystem< double > *l, nonLinearMechSolver::scheme sh, bool useGlobalSystemPosition, bool part):
+ staticDofManager(l,sh,useGlobalSystemPosition,part){};
+ staticDofManagerFullCG(linearSystem<double> *l1, linearSystem<double> *l2, bool part):
+ staticDofManager(l1,l2,part){};
+ virtual ~staticDofManagerFullCG(){};
+
+ #if defined(HAVE_MPI)
+ virtual int sizeOfR() const;
+ virtual void allocateSystem();
+ virtual void setFirstRigidContactUnknowns();
+ virtual void manageMPIComm(const int otherPartNum,const std::vector<Dof> &otherPartR);
+ virtual void manageMPISparsity(const std::vector<Dof> &myR, const int otherPartNum,const std::vector<Dof> &otherPartR);
+ virtual void systemMPIComm();
+ virtual void getDofValue(const Dof& key, double &val, nonLinearBoundaryCondition::whichCondition wv) const;
+ #endif // HAVE_MPI
+};
#endif // STATICDOFMANAGERSYSTEM_H_
diff --git a/NonLinearSolver/pathFollowing/generalLocalControlBasedPathFollowingSystemPETSc.h b/NonLinearSolver/pathFollowing/generalLocalControlBasedPathFollowingSystemPETSc.h
index 9020f8a01e5b5503827a6d53a0c3cf5378f70af8..dd3136cc916bd51d99a9855057bc11a3156d2671 100644
--- a/NonLinearSolver/pathFollowing/generalLocalControlBasedPathFollowingSystemPETSc.h
+++ b/NonLinearSolver/pathFollowing/generalLocalControlBasedPathFollowingSystemPETSc.h
@@ -176,12 +176,12 @@ class generalLocalControlBasedPathFollowingSystemPETSc : public nonLinearSystem
if(!this->_flagb){
#if defined(HAVE_MPI)
if(comSize > 1){
- _try(VecAssemblyBegin(this->_q));
- _try(VecAssemblyEnd(this->_q));
- _try(VecAssemblyBegin(this->_Fint));
- _try(VecAssemblyEnd(this->_Fint));
- _try(VecAssemblyBegin(this->_Fext));
- _try(VecAssemblyEnd(this->_Fext));
+ _try(VecAssemblyBegin(this->_q));
+ _try(VecAssemblyEnd(this->_q));
+ _try(VecAssemblyBegin(this->_Fint));
+ _try(VecAssemblyEnd(this->_Fint));
+ _try(VecAssemblyBegin(this->_Fext));
+ _try(VecAssemblyEnd(this->_Fext));
}
#endif // HAVE_MPI
_try(VecCopy(this->_Fext,this->_b));
@@ -190,7 +190,9 @@ class generalLocalControlBasedPathFollowingSystemPETSc : public nonLinearSystem
}
linearSystemPETSc<scalar>::systemSolve();
_try(VecAXPY(this->_xsol,1,this->_x));
-
+
+ _try(VecAssemblyBegin(this->_g));
+ _try(VecAssemblyEnd(this->_g));
PetscScalar gTdr(0.);
_try(VecDot(_g,this->_x,&gTdr));
@@ -235,8 +237,8 @@ class generalLocalControlBasedPathFollowingSystemPETSc : public nonLinearSystem
_try(VecAssemblyEnd(this->_q));
_try(VecAssemblyBegin(this->_Fint));
_try(VecAssemblyEnd(this->_Fint));
- _try(VecAssemblyBegin(this->_Fext));
- _try(VecAssemblyEnd(this->_Fext));
+ _try(VecAssemblyBegin(this->_Fext));
+ _try(VecAssemblyEnd(this->_Fext));
}
#endif // HAVE_MPI
_try(VecCopy(this->_Fext,this->_b));
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/comparison.py b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/comparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d157fde2846147b900cd54a67bc8a0ffbc47a5f
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/comparison.py
@@ -0,0 +1,22 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+plt.figure()
+
+u = pd.read_csv("cube-fortran/NodalDisplacementPhysical43Num2comp0.csv",sep=";",header=None)
+f = pd.read_csv("cube-fortran/force30comp0.csv",sep=";",header=None)
+plt.plot(u.values[:,1],f.values[:,1],"r-",label="CUR -fortran")
+
+u = pd.read_csv("cube/NodalDisplacementPhysical43Num2comp0.csv",sep=";",header=None)
+f = pd.read_csv("cube/force30comp0.csv",sep=";",header=None)
+plt.plot(u.values[:,1],f.values[:,1],"g--",label="ref")
+
+
+plt.xlabel("disp")
+plt.ylabel("force")
+
+plt.legend()
+
+plt.show()
+
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.geo b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.geo
new file mode 100644
index 0000000000000000000000000000000000000000..9654d90c5497bf8a1ff446e39d1b9ad26a25086c
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.geo
@@ -0,0 +1,51 @@
+// Cube.geo
+mm=1.0; // Units
+n=1; // Number of layers (in thickness / along z)
+m = 1; // Number of element + 1 along y
+p = 3; // Number of element + 1 along x (min 3 to have a crack)
+L=2.*mm; // Cube lenght
+sl1=L/n;
+
+// Face z = 0 (point in anti-clockwise order)
+Point(1)={0,0,0,sl1};
+Point(2)={L,0,0,sl1};
+Point(3)={L,L,0,sl1};
+Point(4)={0,L,0,sl1};
+Line(1)={1,2};
+Line(2)={2,3};
+Line(3)={3,4};
+Line(4)={4,1};
+Line Loop(5) = {3, 4, 1, 2};
+Plane Surface(6) = {5};
+
+// Cube extrusion
+Extrude {0, 0, L} {
+ Surface{6}; Layers{n}; Recombine;
+}
+
+// Physical entities
+ // Volume
+Physical Volume(29) = {1};
+ // Surface
+Physical Surface(30) = {19}; // Left face x = 0
+Physical Surface(31) = {27}; // Right face x = L
+Physical Surface(32) = {6}; // Back face z = 0
+Physical Surface(33) = {28}; // Front face z = L
+Physical Surface(34) = {15}; // Down face y = 0
+Physical Surface(35) = {23}; // Up face y = L
+
+// Mesh
+Transfinite Line {2, 4} = m Using Progression 1;
+Transfinite Line {1, 3} = p Using Progression 1;
+Transfinite Surface {6};
+Recombine Surface{6};
+
+
+// To save cube elongation
+Physical Point(41) = {1}; // Point on left face (0,0,0)
+Physical Point(42) = {4}; // Point on left face (0,L,0)
+Physical Point(43) = {2}; // Point on right face (L,0,0)
+Physical Point(44) = {3}; // Point on right face (L,L,0)
+
+
+
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.msh b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.msh
new file mode 100644
index 0000000000000000000000000000000000000000..a7382df16647c759d6e8e08f9a57c54f24e473bd
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.msh
@@ -0,0 +1,109 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+8 12 6 1
+1 0 0 0 1 41
+2 2 0 0 1 43
+3 2 2 0 1 44
+4 0 2 0 1 42
+5 2 2 2 0
+6 0 2 2 0
+10 0 0 2 0
+14 2 0 2 0
+1 0 0 0 2 0 0 0 2 1 -2
+2 2 0 0 2 2 0 0 2 2 -3
+3 0 2 0 2 2 0 0 2 3 -4
+4 0 0 0 0 2 0 0 2 4 -1
+8 0 2 2 2 2 2 0 2 5 -6
+9 0 0 2 0 2 2 0 2 6 -10
+10 0 0 2 2 0 2 0 2 10 -14
+11 2 0 2 2 2 2 0 2 14 -5
+13 2 2 0 2 2 2 0 2 3 -5
+14 0 2 0 0 2 2 0 2 4 -6
+18 0 0 0 0 0 2 0 2 1 -10
+22 2 0 0 2 0 2 0 2 2 -14
+6 0 0 0 2 2 0 1 32 4 3 4 1 2
+15 0 2 0 2 2 2 1 34 4 3 14 -8 -13
+19 0 0 0 0 2 2 1 30 4 4 18 -9 -14
+23 0 0 0 2 0 2 1 35 4 1 22 -10 -18
+27 2 0 0 2 2 2 1 31 4 2 13 -11 -22
+28 0 0 2 2 2 2 1 33 4 8 9 10 11
+1 0 0 0 2 2 2 1 29 6 -6 28 15 19 23 27
+$EndEntities
+$Nodes
+19 12 1 12
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+2 0 0
+0 3 0 1
+3
+2 2 0
+0 4 0 1
+4
+0 2 0
+0 5 0 1
+5
+2 2 2
+0 6 0 1
+6
+0 2 2
+0 10 0 1
+7
+0 0 2
+0 14 0 1
+8
+2 0 2
+1 1 0 1
+9
+0.9999999999973842 0 0
+1 3 0 1
+10
+1.000000000004119 2 0
+1 8 0 1
+11
+1.000000000004119 2 2
+1 10 0 1
+12
+0.9999999999973842 0 2
+2 6 0 0
+2 15 0 0
+2 19 0 0
+2 23 0 0
+2 27 0 0
+2 28 0 0
+3 1 0 0
+$EndNodes
+$Elements
+11 16 1 16
+0 1 15 1
+1 1
+0 2 15 1
+2 2
+0 3 15 1
+3 3
+0 4 15 1
+4 4
+2 6 3 2
+5 3 10 9 2
+6 10 4 1 9
+2 15 3 2
+7 3 10 11 5
+8 10 4 6 11
+2 19 3 1
+9 4 1 7 6
+2 23 3 2
+10 1 9 12 7
+11 9 2 8 12
+2 27 3 1
+12 2 3 5 8
+2 28 3 2
+13 5 11 12 8
+14 11 6 7 12
+3 1 5 2
+15 3 10 9 2 5 11 12 8
+16 10 4 1 9 11 6 7 12
+$EndElements
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.py b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.py
new file mode 100644
index 0000000000000000000000000000000000000000..cc02ad36df74242ecbd9d423246b86f5f9efa950
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/cube.py
@@ -0,0 +1,88 @@
+#coding-Utf-8-*-
+
+from gmshpy import *
+from dG3Dpy import*
+from math import*
+
+#script to launch PBC problem with a python script
+
+# material law
+lawnum = 11 # unique number of law
+
+# material law
+lawnum = 11 # unique number of law
+law1 = VEVPUMATDG3DMaterialLaw(lawnum,'material.i02');
+#lawPert1 = dG3DMaterialLawWithTangentByPerturbation(law1, 1e-8)
+
+
+# geometry
+meshfile="cube.msh" # name of mesh file
+
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1 = 100
+# creation of part Domain
+nfield = 29 # number of the field (physical number of surface)
+
+myfield1 = dG3DDomain(1000,nfield,space1,lawnum,fullDg,3)
+#myfield1.matrixByPerturbation(1,1,1,1e-8)
+myfield1.stabilityParameters(beta1)
+
+# solver
+sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 100 # number of step (used only if soltype=1)
+ftime =1. # Final time (used only if soltype=1)
+tol=1.e-8 # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+
+
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+#mysolver.addMaterialLaw(lawPert1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+
+mysolver.displacementBC('Face',32,2,0.)
+mysolver.displacementBC('Face',34,1,0.)
+mysolver.displacementBC('Face',30,0,0.)
+
+#mysolver.displacementBC('Face',1,0,0)
+#mysolver.displacementBC('Face',1,1,0)
+mysolver.displacementBC('Face',31,0,1)
+
+# build view
+mysolver.internalPointBuildView("Strain_xx",IPField.STRAIN_XX, 1, 1);
+mysolver.internalPointBuildView("Strain_yy",IPField.STRAIN_YY, 1, 1);
+mysolver.internalPointBuildView("Strain_zz",IPField.STRAIN_ZZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xy",IPField.STRAIN_XY, 1, 1);
+mysolver.internalPointBuildView("Strain_yz",IPField.STRAIN_YZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xz",IPField.STRAIN_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1);
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1);
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1);
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1);
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1);
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
+mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
+mysolver.internalPointBuildView("pression",IPField.PRESSION,1,1)
+mysolver.internalPointBuildView("plastic possion ratio",IPField.PLASTIC_POISSON_RATIO,1,1)
+
+#archiving
+mysolver.archivingForceOnPhysicalGroup('Face',30,0)
+mysolver.archivingForceOnPhysicalGroup('Face',30,1)
+mysolver.archivingForceOnPhysicalGroup('Face',30,2)
+mysolver.archivingNodeDisplacement(43,0,1)
+mysolver.archivingNodeDisplacement(43,1,1)
+mysolver.archivingNodeDisplacement(43,2,1)
+# solve
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-5.878380e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-3)
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/material.i02 b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/material.i02
new file mode 100644
index 0000000000000000000000000000000000000000..9ed77540254687678d4bf96b4dea6cf504bcacfa
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube-fortran/material.i02
@@ -0,0 +1,52 @@
+108 number of internal variables
+50 number of properties to be read
+11 Power Series Strain Approximation
+833.3333333333333 Bulk Modulus
+384.6153846153846 Shear Modulus
+3.5 Yield Exponent
+0.452 Plastic Poisson Ratio
+3e4 Viscoplastic coefficient
+0.21 Viscoplasic exponent
+100.0 Initial Yield Limit - Compression
+300. Isotropic Hardening Parameter - Compression
+0.0 Isotropic Hardening Parameter - Compression
+10. Isotropic Hardening Parameter - Compression
+78.0 Initial Yield Limit - Tension
+300. Isotropic Hardening Parameter - Tension
+0. Isotropic Hardening Parameter - Tension
+10. Isotropic Hardening Parameter - Tension
+0.0 Kinematic Hardening Parameter
+0.0 Kinematic Hardening Parameter
+0.0 Kinematic Hardening Parameter
+8.333333e+01 Bulk Modulus - Maxwell branch 1
+8.333333e+01 Bulk Modulus - Maxwell branch 2
+8.333333e+01 Bulk Modulus - Maxwell branch 3
+8.333333e+01 Bulk Modulus - Maxwell branch 4
+8.333333e+01 Bulk Modulus - Maxwell branch 5
+8.333333e+01 Bulk Modulus - Maxwell branch 6
+8.333333e+01 Bulk Modulus - Maxwell branch 7
+8.333333e+01 Bulk Modulus - Maxwell branch 8
+1e-4 Retardation Time Volumetric - Maxwell Branch 1
+1e-3 Retardation Time Volumetric - Maxwell Branch 2
+1e-2 Retardation Time Volumetric - Maxwell Branch 3
+1e-1 Retardation Time Volumetric - Maxwell Branch 4
+1. Retardation Time Volumetric - Maxwell Branch 5
+1e1 Retardation Time Volumetric - Maxwell Branch 6
+1e2 Retardation Time Volumetric - Maxwell Branch 7
+1e3 Retardation Time Volumetric - Maxwell Branch 8
+3.846154e+01 Shear Modulus- Maxwell branch 1
+3.846154e+01 Shear Modulus- Maxwell branch 2
+3.846154e+01 Shear Modulus- Maxwell branch 3
+3.846154e+01 Shear Modulus- Maxwell branch 4
+3.846154e+01 Shear Modulus- Maxwell branch 5
+3.846154e+01 Shear Modulus- Maxwell branch 6
+3.846154e+01 Shear Modulus- Maxwell branch 7
+3.846154e+01 Shear Modulus- Maxwell branch 8
+1e-4 Retardation Time Deviatoric - Maxwell Branch 1
+1e-3 Retardation Time Deviatoric - Maxwell Branch 2
+1e-2 Retardation Time Deviatoric - Maxwell Branch 3
+1e-1 Retardation Time Deviatoric - Maxwell Branch 4
+1 Retardation Time Deviatoric - Maxwell Branch 5
+1e1 Retardation Time Deviatoric - Maxwell Branch 6
+1e2 Retardation Time Deviatoric - Maxwell Branch 7
+1e3 Retardation Time Deviatoric - Maxwell Branch 8
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.geo b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.geo
new file mode 100644
index 0000000000000000000000000000000000000000..9654d90c5497bf8a1ff446e39d1b9ad26a25086c
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.geo
@@ -0,0 +1,51 @@
+// Cube.geo
+mm=1.0; // Units
+n=1; // Number of layers (in thickness / along z)
+m = 1; // Number of element + 1 along y
+p = 3; // Number of element + 1 along x (min 3 to have a crack)
+L=2.*mm; // Cube lenght
+sl1=L/n;
+
+// Face z = 0 (point in anti-clockwise order)
+Point(1)={0,0,0,sl1};
+Point(2)={L,0,0,sl1};
+Point(3)={L,L,0,sl1};
+Point(4)={0,L,0,sl1};
+Line(1)={1,2};
+Line(2)={2,3};
+Line(3)={3,4};
+Line(4)={4,1};
+Line Loop(5) = {3, 4, 1, 2};
+Plane Surface(6) = {5};
+
+// Cube extrusion
+Extrude {0, 0, L} {
+ Surface{6}; Layers{n}; Recombine;
+}
+
+// Physical entities
+ // Volume
+Physical Volume(29) = {1};
+ // Surface
+Physical Surface(30) = {19}; // Left face x = 0
+Physical Surface(31) = {27}; // Right face x = L
+Physical Surface(32) = {6}; // Back face z = 0
+Physical Surface(33) = {28}; // Front face z = L
+Physical Surface(34) = {15}; // Down face y = 0
+Physical Surface(35) = {23}; // Up face y = L
+
+// Mesh
+Transfinite Line {2, 4} = m Using Progression 1;
+Transfinite Line {1, 3} = p Using Progression 1;
+Transfinite Surface {6};
+Recombine Surface{6};
+
+
+// To save cube elongation
+Physical Point(41) = {1}; // Point on left face (0,0,0)
+Physical Point(42) = {4}; // Point on left face (0,L,0)
+Physical Point(43) = {2}; // Point on right face (L,0,0)
+Physical Point(44) = {3}; // Point on right face (L,L,0)
+
+
+
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.msh b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.msh
new file mode 100644
index 0000000000000000000000000000000000000000..a7382df16647c759d6e8e08f9a57c54f24e473bd
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.msh
@@ -0,0 +1,109 @@
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$Entities
+8 12 6 1
+1 0 0 0 1 41
+2 2 0 0 1 43
+3 2 2 0 1 44
+4 0 2 0 1 42
+5 2 2 2 0
+6 0 2 2 0
+10 0 0 2 0
+14 2 0 2 0
+1 0 0 0 2 0 0 0 2 1 -2
+2 2 0 0 2 2 0 0 2 2 -3
+3 0 2 0 2 2 0 0 2 3 -4
+4 0 0 0 0 2 0 0 2 4 -1
+8 0 2 2 2 2 2 0 2 5 -6
+9 0 0 2 0 2 2 0 2 6 -10
+10 0 0 2 2 0 2 0 2 10 -14
+11 2 0 2 2 2 2 0 2 14 -5
+13 2 2 0 2 2 2 0 2 3 -5
+14 0 2 0 0 2 2 0 2 4 -6
+18 0 0 0 0 0 2 0 2 1 -10
+22 2 0 0 2 0 2 0 2 2 -14
+6 0 0 0 2 2 0 1 32 4 3 4 1 2
+15 0 2 0 2 2 2 1 34 4 3 14 -8 -13
+19 0 0 0 0 2 2 1 30 4 4 18 -9 -14
+23 0 0 0 2 0 2 1 35 4 1 22 -10 -18
+27 2 0 0 2 2 2 1 31 4 2 13 -11 -22
+28 0 0 2 2 2 2 1 33 4 8 9 10 11
+1 0 0 0 2 2 2 1 29 6 -6 28 15 19 23 27
+$EndEntities
+$Nodes
+19 12 1 12
+0 1 0 1
+1
+0 0 0
+0 2 0 1
+2
+2 0 0
+0 3 0 1
+3
+2 2 0
+0 4 0 1
+4
+0 2 0
+0 5 0 1
+5
+2 2 2
+0 6 0 1
+6
+0 2 2
+0 10 0 1
+7
+0 0 2
+0 14 0 1
+8
+2 0 2
+1 1 0 1
+9
+0.9999999999973842 0 0
+1 3 0 1
+10
+1.000000000004119 2 0
+1 8 0 1
+11
+1.000000000004119 2 2
+1 10 0 1
+12
+0.9999999999973842 0 2
+2 6 0 0
+2 15 0 0
+2 19 0 0
+2 23 0 0
+2 27 0 0
+2 28 0 0
+3 1 0 0
+$EndNodes
+$Elements
+11 16 1 16
+0 1 15 1
+1 1
+0 2 15 1
+2 2
+0 3 15 1
+3 3
+0 4 15 1
+4 4
+2 6 3 2
+5 3 10 9 2
+6 10 4 1 9
+2 15 3 2
+7 3 10 11 5
+8 10 4 6 11
+2 19 3 1
+9 4 1 7 6
+2 23 3 2
+10 1 9 12 7
+11 9 2 8 12
+2 27 3 1
+12 2 3 5 8
+2 28 3 2
+13 5 11 12 8
+14 11 6 7 12
+3 1 5 2
+15 3 10 9 2 5 11 12 8
+16 10 4 1 9 11 6 7 12
+$EndElements
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.py b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.py
new file mode 100644
index 0000000000000000000000000000000000000000..7dac5eddfa09cae08637609319e2bdb0b0b3298b
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic-fortran/cube/cube.py
@@ -0,0 +1,124 @@
+#coding-Utf-8-*-
+
+from gmshpy import *
+from dG3Dpy import*
+from math import*
+
+#script to launch PBC problem with a python script
+
+# material law
+lawnum = 11 # unique number of law
+
+E = 1000. #MPa
+nu = 0.3
+K = E/3./(1.-2.*nu) # Bulk mudulus
+mu =E/2./(1.+nu) # Shear mudulus
+print(f"K={K} mu={mu}")
+rho = 7850e-9 # Bulk mass
+
+sy0c = 100. #MPa
+hc = 300.
+
+sy0t = sy0c*0.78
+ht = 300.
+
+
+# creation of law
+hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, 0., 10.)
+hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 0., 10.)
+#hardenk = PolynomialKinematicHardening(3,1)
+#hardenk.setCoefficients(1,0.)
+
+law1 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,False,1e-8)
+law1.setCompressionHardening(hardenc)
+law1.setTractionHardening(hardent)
+#law1.setKinematicHardening(hardenk)
+
+law1.setStrainOrder(11)
+
+law1.setYieldPowerFactor(3.5)
+law1.setNonAssociatedFlow(True)
+law1.setPlasticPoissonRatio(0.452)
+
+law1.setViscoelasticMethod(0)
+
+
+N = 8;
+law1.setViscoElasticNumberOfElement(N)
+for i in range(N):
+ law1.setViscoElasticData(i+1,100.,10**(i-4))
+
+eta = constantViscosityLaw(1,3e4)
+law1.setViscosityEffect(eta,0.21)
+
+
+# geometry
+meshfile="cube.msh" # name of mesh file
+
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1 = 100
+# creation of part Domain
+nfield = 29 # number of the field (physical number of surface)
+
+myfield1 = dG3DDomain(1000,nfield,space1,lawnum,fullDg,3)
+#myfield1.matrixByPerturbation(1,1,1,1e-8)
+myfield1.stabilityParameters(beta1)
+
+# solver
+sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 100 # number of step (used only if soltype=1)
+ftime =1. # Final time (used only if soltype=1)
+tol=1.e-8 # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+
+
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+
+mysolver.displacementBC('Face',32,2,0.)
+mysolver.displacementBC('Face',34,1,0.)
+mysolver.displacementBC('Face',30,0,0.)
+
+#mysolver.displacementBC('Face',1,0,0)
+#mysolver.displacementBC('Face',1,1,0)
+mysolver.displacementBC('Face',31,0,1.)
+
+# build view
+mysolver.internalPointBuildView("Strain_xx",IPField.STRAIN_XX, 1, 1);
+mysolver.internalPointBuildView("Strain_yy",IPField.STRAIN_YY, 1, 1);
+mysolver.internalPointBuildView("Strain_zz",IPField.STRAIN_ZZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xy",IPField.STRAIN_XY, 1, 1);
+mysolver.internalPointBuildView("Strain_yz",IPField.STRAIN_YZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xz",IPField.STRAIN_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1);
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1);
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1);
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1);
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1);
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
+mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
+mysolver.internalPointBuildView("pression",IPField.PRESSION,1,1)
+mysolver.internalPointBuildView("plastic possion ratio",IPField.PLASTIC_POISSON_RATIO,1,1)
+
+#archiving
+mysolver.archivingForceOnPhysicalGroup('Face',30,0)
+mysolver.archivingForceOnPhysicalGroup('Face',30,1)
+mysolver.archivingForceOnPhysicalGroup('Face',30,2)
+mysolver.archivingNodeDisplacement(43,0,1)
+mysolver.archivingNodeDisplacement(43,1,1)
+mysolver.archivingNodeDisplacement(43,2,1)
+# solve
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-5.878380e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", 30, 0),1.e-3)
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic/cylinder.py b/dG3D/benchmarks/powerYieldViscoElastoPlastic/cylinder.py
index b7beb05958c1fa46550fd5d82255a0789bd00943..be503f1fb1a242c7d61ce258ba3ce1f025ad52d3 100644
--- a/dG3D/benchmarks/powerYieldViscoElastoPlastic/cylinder.py
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic/cylinder.py
@@ -14,6 +14,7 @@ E = 2.7E3 #MPa
nu = 0.3
K = E/3./(1.-2.*nu) # Bulk mudulus
mu =E/2./(1.+nu) # Shear mudulus
+print(f"K={K} mu={mu}")
rho = 7850e-9 # Bulk mass
sy0c = 100. #MPa
@@ -26,21 +27,26 @@ ht = 300.
# creation of law
hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, 0., 10.)
hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 0., 10.)
-hardenk = PolynomialKinematicHardening(3,1)
-hardenk.setCoefficients(1,370.)
+hardenk = PolynomialKinematicHardening(3,2)
+hardenk.setCoefficients(0,0)
+hardenk.setCoefficients(1,1.e2)
+hardenk.setCoefficients(2,1.e3)
+
law1 = HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,False,1e-8)
law1.setCompressionHardening(hardenc)
law1.setTractionHardening(hardent)
law1.setKinematicHardening(hardenk)
-law1.setStrainOrder(5)
+law1.setStrainOrder(11)
law1.setYieldPowerFactor(3.5)
law1.setNonAssociatedFlow(True)
-law1.nonAssociatedFlowRuleFactor(0.5)
+law1.setPlasticPoissonRatio(0.452)
law1.setViscoelasticMethod(0)
+
+
N = 4;
law1.setViscoElasticNumberOfElement(N)
law1.setViscoElasticData(1,1380,7202.)
@@ -121,4 +127,4 @@ mysolver.archivingNodeDisplacement(5,2,1)
mysolver.solve()
check = TestCheck()
-check.equal(-5.609926e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 1, 2),1.e-3)
+check.equal(-5.246765e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 1, 2),1.e-3)
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/comparison.py b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/comparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..d7070d78feed235392eb0a9263c304747caabe87
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/comparison.py
@@ -0,0 +1,22 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+plt.figure()
+
+u = pd.read_csv("NodalDisplacementPhysical5Num5comp2.csv",sep=";",header=None)
+f = pd.read_csv("force1comp2.csv",sep=";",header=None)
+plt.plot(u.values[:,1],f.values[:,1],"-",label="CUR -fortran")
+
+u = pd.read_csv("../NodalDisplacementPhysical5Num5comp2.csv",sep=";",header=None)
+f = pd.read_csv("../force1comp2.csv",sep=";",header=None)
+plt.plot(u.values[:,1],f.values[:,1],"--",label="ref")
+
+
+plt.xlabel("disp")
+plt.ylabel("force")
+
+plt.legend()
+
+plt.show()
+
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/cylinder.msh b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/cylinder.msh
new file mode 100644
index 0000000000000000000000000000000000000000..0f52bfc5a294857b8f304b1da6e258fb7505d374
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/cylinder.msh
@@ -0,0 +1,698 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+132
+1 0 0 0
+2 6 0 0
+3 0 6 0
+4 0 6 6
+5 0 0 6
+6 6 0 6
+7 1.552914270611066 5.795554957735497 0
+8 2.999999999993018 5.196152422710663 0
+9 4.242640687111061 4.242640687127508 0
+10 5.196152422702976 3.000000000006332 0
+11 5.795554957733712 1.552914270617726 0
+12 0 4.500000000005166 0
+13 0 3.000000000008336 0
+14 0 1.500000000004214 0
+15 1.499999999993832 0 0
+16 2.999999999988567 0 0
+17 4.499999999994237 0 0
+18 1.552914270611066 5.795554957735497 6
+19 2.999999999993018 5.196152422710663 6
+20 4.242640687111061 4.242640687127508 6
+21 5.196152422702976 3.000000000006332 6
+22 5.795554957733712 1.552914270617726 6
+23 4.499999999994237 0 6
+24 2.999999999988567 0 6
+25 1.499999999993832 0 6
+26 0 1.500000000004214 6
+27 0 3.000000000008336 6
+28 0 4.500000000005166 6
+29 0 6 1.2
+30 0 6 2.4
+31 0 6 3.6
+32 0 6 4.800000000000001
+33 6 0 1.2
+34 6 0 2.4
+35 6 0 3.6
+36 6 0 4.800000000000001
+37 0 0 1.2
+38 0 0 2.4
+39 0 0 3.6
+40 0 0 4.800000000000001
+41 2.527448697695314 2.575668272009054 0
+42 1.627841394027458 3.980250003805291 0
+43 3.914833052601809 1.621430334871464 0
+44 1.441200958338174 1.441200958338581 0
+45 3.760905731444732 3.080330269885938 0
+46 2.682328465952599 1.220004215439306 0
+47 1.220004215434672 2.682328465952725 0
+48 3.000820586356653 3.856980823767017 0
+49 1.552914270611066 5.795554957735497 1.2
+50 1.552914270611066 5.795554957735497 2.4
+51 1.552914270611066 5.795554957735497 3.6
+52 1.552914270611066 5.795554957735497 4.800000000000001
+53 2.999999999993018 5.196152422710663 1.2
+54 2.999999999993018 5.196152422710663 2.4
+55 2.999999999993018 5.196152422710663 3.6
+56 2.999999999993018 5.196152422710663 4.800000000000001
+57 4.242640687111061 4.242640687127508 1.2
+58 4.242640687111061 4.242640687127508 2.4
+59 4.242640687111061 4.242640687127508 3.6
+60 4.242640687111061 4.242640687127508 4.800000000000001
+61 5.196152422702976 3.000000000006332 1.2
+62 5.196152422702976 3.000000000006332 2.4
+63 5.196152422702976 3.000000000006332 3.6
+64 5.196152422702976 3.000000000006332 4.800000000000001
+65 5.795554957733712 1.552914270617726 1.2
+66 5.795554957733712 1.552914270617726 2.4
+67 5.795554957733712 1.552914270617726 3.6
+68 5.795554957733712 1.552914270617726 4.800000000000001
+69 1.499999999993832 0 1.2
+70 1.499999999993832 0 2.4
+71 1.499999999993832 0 3.6
+72 1.499999999993832 0 4.800000000000001
+73 2.999999999988567 0 1.2
+74 2.999999999988567 0 2.4
+75 2.999999999988567 0 3.6
+76 2.999999999988567 0 4.800000000000001
+77 4.499999999994237 0 1.2
+78 4.499999999994237 0 2.4
+79 4.499999999994237 0 3.6
+80 4.499999999994237 0 4.800000000000001
+81 0 4.500000000005166 1.2
+82 0 4.500000000005166 2.4
+83 0 4.500000000005166 3.6
+84 0 4.500000000005166 4.800000000000001
+85 0 3.000000000008336 1.2
+86 0 3.000000000008336 2.4
+87 0 3.000000000008336 3.6
+88 0 3.000000000008336 4.800000000000001
+89 0 1.500000000004214 1.2
+90 0 1.500000000004214 2.4
+91 0 1.500000000004214 3.6
+92 0 1.500000000004214 4.800000000000001
+93 2.527448697695314 2.575668272009054 6
+94 1.627841394027458 3.980250003805291 6
+95 3.914833052601809 1.621430334871464 6
+96 1.441200958338174 1.441200958338581 6
+97 3.760905731444732 3.080330269885938 6
+98 2.682328465952599 1.220004215439306 6
+99 1.220004215434672 2.682328465952725 6
+100 3.000820586356653 3.856980823767017 6
+101 2.527448697695314 2.575668272009054 1.2
+102 2.527448697695314 2.575668272009054 2.4
+103 2.527448697695314 2.575668272009054 3.6
+104 2.527448697695314 2.575668272009054 4.800000000000001
+105 1.627841394027458 3.980250003805291 1.2
+106 1.627841394027458 3.980250003805291 2.4
+107 1.627841394027458 3.980250003805291 3.6
+108 1.627841394027458 3.980250003805291 4.800000000000001
+109 3.914833052601809 1.621430334871464 1.2
+110 3.914833052601809 1.621430334871464 2.4
+111 3.914833052601809 1.621430334871464 3.6
+112 3.914833052601809 1.621430334871464 4.800000000000001
+113 1.441200958338174 1.441200958338581 1.2
+114 1.441200958338174 1.441200958338581 2.4
+115 1.441200958338174 1.441200958338581 3.6
+116 1.441200958338174 1.441200958338581 4.800000000000001
+117 3.760905731444732 3.080330269885938 1.2
+118 3.760905731444732 3.080330269885938 2.4
+119 3.760905731444732 3.080330269885938 3.6
+120 3.760905731444732 3.080330269885938 4.800000000000001
+121 2.682328465952599 1.220004215439306 1.2
+122 2.682328465952599 1.220004215439306 2.4
+123 2.682328465952599 1.220004215439306 3.6
+124 2.682328465952599 1.220004215439306 4.800000000000001
+125 1.220004215434672 2.682328465952725 1.2
+126 1.220004215434672 2.682328465952725 2.4
+127 1.220004215434672 2.682328465952725 3.6
+128 1.220004215434672 2.682328465952725 4.800000000000001
+129 3.000820586356653 3.856980823767017 1.2
+130 3.000820586356653 3.856980823767017 2.4
+131 3.000820586356653 3.856980823767017 3.6
+132 3.000820586356653 3.856980823767017 4.800000000000001
+$EndNodes
+$Elements
+557
+1 15 2 5 5 5
+2 2 2 1 5 7 42 12
+3 2 2 1 5 11 17 43
+4 2 2 1 5 16 43 17
+5 2 2 1 5 12 42 13
+6 2 2 1 5 1 44 15
+7 2 2 1 5 1 14 44
+8 2 2 1 5 3 7 12
+9 2 2 1 5 11 2 17
+10 2 2 1 5 16 46 43
+11 2 2 1 5 13 42 47
+12 2 2 1 5 7 8 42
+13 2 2 1 5 10 11 43
+14 2 2 1 5 41 42 48
+15 2 2 1 5 10 43 45
+16 2 2 1 5 15 46 16
+17 2 2 1 5 13 47 14
+18 2 2 1 5 41 47 42
+19 2 2 1 5 41 43 46
+20 2 2 1 5 41 45 43
+21 2 2 1 5 15 44 46
+22 2 2 1 5 14 47 44
+23 2 2 1 5 8 48 42
+24 2 2 1 5 9 10 45
+25 2 2 1 5 41 48 45
+26 2 2 1 5 8 9 48
+27 2 2 1 5 41 46 44
+28 2 2 1 5 41 44 47
+29 2 2 1 5 9 45 48
+30 2 2 3 17 37 69 1
+31 2 2 3 17 1 69 15
+32 2 2 3 17 38 70 37
+33 2 2 3 17 37 70 69
+34 2 2 3 17 39 71 38
+35 2 2 3 17 38 71 70
+36 2 2 3 17 40 72 39
+37 2 2 3 17 39 72 71
+38 2 2 3 17 5 25 40
+39 2 2 3 17 40 25 72
+40 2 2 3 17 69 73 15
+41 2 2 3 17 15 73 16
+42 2 2 3 17 70 74 69
+43 2 2 3 17 69 74 73
+44 2 2 3 17 71 75 70
+45 2 2 3 17 70 75 74
+46 2 2 3 17 72 76 71
+47 2 2 3 17 71 76 75
+48 2 2 3 17 25 24 72
+49 2 2 3 17 72 24 76
+50 2 2 3 17 73 77 16
+51 2 2 3 17 16 77 17
+52 2 2 3 17 74 78 73
+53 2 2 3 17 73 78 77
+54 2 2 3 17 75 79 74
+55 2 2 3 17 74 79 78
+56 2 2 3 17 76 80 75
+57 2 2 3 17 75 80 79
+58 2 2 3 17 24 23 76
+59 2 2 3 17 76 23 80
+60 2 2 3 17 77 2 17
+61 2 2 3 17 77 33 2
+62 2 2 3 17 78 33 77
+63 2 2 3 17 78 34 33
+64 2 2 3 17 79 34 78
+65 2 2 3 17 79 35 34
+66 2 2 3 17 80 35 79
+67 2 2 3 17 80 36 35
+68 2 2 3 17 23 36 80
+69 2 2 3 17 23 6 36
+70 2 2 4 21 29 81 3
+71 2 2 4 21 3 81 12
+72 2 2 4 21 30 82 29
+73 2 2 4 21 29 82 81
+74 2 2 4 21 31 83 30
+75 2 2 4 21 30 83 82
+76 2 2 4 21 32 84 31
+77 2 2 4 21 31 84 83
+78 2 2 4 21 4 28 32
+79 2 2 4 21 32 28 84
+80 2 2 4 21 81 85 12
+81 2 2 4 21 12 85 13
+82 2 2 4 21 82 86 81
+83 2 2 4 21 81 86 85
+84 2 2 4 21 83 86 82
+85 2 2 4 21 83 87 86
+86 2 2 4 21 84 88 83
+87 2 2 4 21 83 88 87
+88 2 2 4 21 28 27 84
+89 2 2 4 21 84 27 88
+90 2 2 4 21 85 89 13
+91 2 2 4 21 13 89 14
+92 2 2 4 21 86 90 85
+93 2 2 4 21 85 90 89
+94 2 2 4 21 87 91 86
+95 2 2 4 21 86 91 90
+96 2 2 4 21 88 92 87
+97 2 2 4 21 87 92 91
+98 2 2 4 21 27 26 88
+99 2 2 4 21 88 26 92
+100 2 2 4 21 89 1 14
+101 2 2 4 21 89 37 1
+102 2 2 4 21 90 37 89
+103 2 2 4 21 90 38 37
+104 2 2 4 21 91 38 90
+105 2 2 4 21 91 39 38
+106 2 2 4 21 92 39 91
+107 2 2 4 21 92 40 39
+108 2 2 4 21 26 40 92
+109 2 2 4 21 26 5 40
+110 2 2 2 22 18 94 28
+111 2 2 2 22 22 23 95
+112 2 2 2 22 24 95 23
+113 2 2 2 22 28 94 27
+114 2 2 2 22 5 96 25
+115 2 2 2 22 5 26 96
+116 2 2 2 22 4 18 28
+117 2 2 2 22 22 6 23
+118 2 2 2 22 24 98 95
+119 2 2 2 22 27 94 99
+120 2 2 2 22 18 19 94
+121 2 2 2 22 21 22 95
+122 2 2 2 22 93 94 100
+123 2 2 2 22 21 95 97
+124 2 2 2 22 25 98 24
+125 2 2 2 22 27 99 26
+126 2 2 2 22 93 99 94
+127 2 2 2 22 93 95 98
+128 2 2 2 22 93 97 95
+129 2 2 2 22 25 96 98
+130 2 2 2 22 26 99 96
+131 2 2 2 22 19 100 94
+132 2 2 2 22 20 21 97
+133 2 2 2 22 93 100 97
+134 2 2 2 22 19 20 100
+135 2 2 2 22 93 98 96
+136 2 2 2 22 93 96 99
+137 2 2 2 22 20 97 100
+138 4 2 11 1 7 12 42 105
+139 4 2 11 1 81 49 105 12
+140 4 2 11 1 12 49 105 7
+141 4 2 11 1 49 81 105 106
+142 4 2 11 1 82 50 106 49
+143 4 2 11 1 81 49 82 106
+144 4 2 11 1 50 82 106 83
+145 4 2 11 1 83 51 107 50
+146 4 2 11 1 106 50 83 107
+147 4 2 11 1 51 83 107 84
+148 4 2 11 1 84 52 108 51
+149 4 2 11 1 107 51 84 108
+150 4 2 11 1 52 84 108 28
+151 4 2 11 1 28 18 94 52
+152 4 2 11 1 108 52 28 94
+153 4 2 11 1 11 43 17 109
+154 4 2 11 1 109 65 77 11
+155 4 2 11 1 17 11 109 77
+156 4 2 11 1 65 109 77 78
+157 4 2 11 1 110 66 78 65
+158 4 2 11 1 109 65 110 78
+159 4 2 11 1 66 110 78 79
+160 4 2 11 1 111 67 79 66
+161 4 2 11 1 110 66 111 79
+162 4 2 11 1 67 111 79 80
+163 4 2 11 1 112 68 80 67
+164 4 2 11 1 111 67 112 80
+165 4 2 11 1 68 112 80 23
+166 4 2 11 1 95 22 23 68
+167 4 2 11 1 112 68 95 23
+168 4 2 11 1 16 17 43 109
+169 4 2 11 1 77 73 109 16
+170 4 2 11 1 17 16 77 109
+171 4 2 11 1 73 77 109 78
+172 4 2 11 1 78 74 110 73
+173 4 2 11 1 109 73 78 110
+174 4 2 11 1 74 78 110 79
+175 4 2 11 1 79 75 111 74
+176 4 2 11 1 110 74 79 111
+177 4 2 11 1 75 79 111 80
+178 4 2 11 1 80 76 112 75
+179 4 2 11 1 111 75 80 112
+180 4 2 11 1 76 80 112 23
+181 4 2 11 1 23 24 95 76
+182 4 2 11 1 112 76 23 95
+183 4 2 11 1 12 13 42 105
+184 4 2 11 1 85 81 105 12
+185 4 2 11 1 13 12 85 105
+186 4 2 11 1 81 85 105 106
+187 4 2 11 1 86 82 106 81
+188 4 2 11 1 85 81 86 106
+189 4 2 11 1 82 86 106 83
+190 4 2 11 1 87 83 107 106
+191 4 2 11 1 86 83 87 106
+192 4 2 11 1 83 87 107 88
+193 4 2 11 1 88 84 108 107
+194 4 2 11 1 83 84 88 107
+195 4 2 11 1 84 88 108 27
+196 4 2 11 1 27 28 94 108
+197 4 2 11 1 84 28 27 108
+198 4 2 11 1 1 15 44 113
+199 4 2 11 1 69 37 113 1
+200 4 2 11 1 15 1 69 113
+201 4 2 11 1 37 69 113 114
+202 4 2 11 1 70 38 114 37
+203 4 2 11 1 69 37 70 114
+204 4 2 11 1 38 70 114 115
+205 4 2 11 1 71 39 115 38
+206 4 2 11 1 70 38 71 115
+207 4 2 11 1 39 71 115 116
+208 4 2 11 1 72 40 116 39
+209 4 2 11 1 71 39 72 116
+210 4 2 11 1 40 72 116 96
+211 4 2 11 1 25 5 96 40
+212 4 2 11 1 72 40 25 96
+213 4 2 11 1 1 44 14 113
+214 4 2 11 1 113 37 89 1
+215 4 2 11 1 14 1 113 89
+216 4 2 11 1 37 113 89 90
+217 4 2 11 1 114 38 90 37
+218 4 2 11 1 113 37 114 90
+219 4 2 11 1 38 114 90 91
+220 4 2 11 1 115 39 91 38
+221 4 2 11 1 114 38 115 91
+222 4 2 11 1 39 115 91 92
+223 4 2 11 1 116 40 92 39
+224 4 2 11 1 115 39 116 92
+225 4 2 11 1 40 116 92 26
+226 4 2 11 1 96 5 26 40
+227 4 2 11 1 116 40 96 26
+228 4 2 11 1 3 12 7 49
+229 4 2 11 1 81 29 49 3
+230 4 2 11 1 12 3 81 49
+231 4 2 11 1 29 81 49 82
+232 4 2 11 1 82 30 50 29
+233 4 2 11 1 49 29 82 50
+234 4 2 11 1 30 82 50 83
+235 4 2 11 1 83 31 51 30
+236 4 2 11 1 50 30 83 51
+237 4 2 11 1 31 83 51 84
+238 4 2 11 1 84 32 52 31
+239 4 2 11 1 51 31 84 52
+240 4 2 11 1 32 84 52 28
+241 4 2 11 1 28 4 18 32
+242 4 2 11 1 52 32 28 18
+243 4 2 11 1 11 17 2 77
+244 4 2 11 1 77 65 33 2
+245 4 2 11 1 11 65 77 2
+246 4 2 11 1 65 77 33 78
+247 4 2 11 1 78 66 34 33
+248 4 2 11 1 65 66 78 33
+249 4 2 11 1 66 78 34 79
+250 4 2 11 1 79 67 35 66
+251 4 2 11 1 34 66 79 35
+252 4 2 11 1 67 79 35 80
+253 4 2 11 1 80 68 36 35
+254 4 2 11 1 67 68 80 35
+255 4 2 11 1 68 80 36 23
+256 4 2 11 1 23 22 6 36
+257 4 2 11 1 68 22 23 36
+258 4 2 11 1 16 43 46 121
+259 4 2 11 1 109 73 121 16
+260 4 2 11 1 43 16 109 121
+261 4 2 11 1 73 109 121 122
+262 4 2 11 1 110 74 122 73
+263 4 2 11 1 109 73 110 122
+264 4 2 11 1 74 110 122 123
+265 4 2 11 1 111 75 123 74
+266 4 2 11 1 110 74 111 123
+267 4 2 11 1 75 111 123 124
+268 4 2 11 1 112 76 124 75
+269 4 2 11 1 111 75 112 124
+270 4 2 11 1 76 112 124 98
+271 4 2 11 1 95 24 98 76
+272 4 2 11 1 112 76 95 98
+273 4 2 11 1 13 47 42 105
+274 4 2 11 1 125 85 105 13
+275 4 2 11 1 47 13 125 105
+276 4 2 11 1 85 125 105 106
+277 4 2 11 1 126 86 106 125
+278 4 2 11 1 125 86 106 85
+279 4 2 11 1 86 126 106 87
+280 4 2 11 1 127 87 107 126
+281 4 2 11 1 126 87 107 106
+282 4 2 11 1 87 127 107 128
+283 4 2 11 1 128 88 108 107
+284 4 2 11 1 87 88 128 107
+285 4 2 11 1 88 128 108 99
+286 4 2 11 1 99 27 94 108
+287 4 2 11 1 88 27 99 108
+288 4 2 11 1 7 42 8 105
+289 4 2 11 1 105 49 53 8
+290 4 2 11 1 7 49 105 8
+291 4 2 11 1 49 105 53 106
+292 4 2 11 1 106 50 54 49
+293 4 2 11 1 53 49 106 54
+294 4 2 11 1 50 106 54 107
+295 4 2 11 1 107 51 55 50
+296 4 2 11 1 54 50 107 55
+297 4 2 11 1 51 107 55 108
+298 4 2 11 1 108 52 56 51
+299 4 2 11 1 55 51 108 56
+300 4 2 11 1 52 108 56 94
+301 4 2 11 1 94 18 19 52
+302 4 2 11 1 56 52 94 19
+303 4 2 11 1 10 43 11 109
+304 4 2 11 1 109 61 65 11
+305 4 2 11 1 10 61 109 11
+306 4 2 11 1 61 109 65 110
+307 4 2 11 1 110 62 66 65
+308 4 2 11 1 61 62 110 65
+309 4 2 11 1 62 110 66 111
+310 4 2 11 1 111 63 67 66
+311 4 2 11 1 62 63 111 66
+312 4 2 11 1 63 111 67 112
+313 4 2 11 1 112 64 68 63
+314 4 2 11 1 67 63 112 68
+315 4 2 11 1 64 112 68 95
+316 4 2 11 1 95 21 22 64
+317 4 2 11 1 68 64 95 22
+318 4 2 11 1 41 48 42 105
+319 4 2 11 1 129 101 105 48
+320 4 2 11 1 48 101 105 41
+321 4 2 11 1 101 129 105 102
+322 4 2 11 1 130 102 106 105
+323 4 2 11 1 129 102 130 105
+324 4 2 11 1 102 130 106 103
+325 4 2 11 1 131 103 107 130
+326 4 2 11 1 130 103 107 106
+327 4 2 11 1 103 131 107 104
+328 4 2 11 1 132 104 108 107
+329 4 2 11 1 131 104 132 107
+330 4 2 11 1 104 132 108 93
+331 4 2 11 1 100 93 94 108
+332 4 2 11 1 132 93 100 108
+333 4 2 11 1 10 45 43 117
+334 4 2 11 1 117 61 109 10
+335 4 2 11 1 43 10 117 109
+336 4 2 11 1 61 117 109 118
+337 4 2 11 1 118 62 110 61
+338 4 2 11 1 109 61 118 110
+339 4 2 11 1 62 118 110 119
+340 4 2 11 1 119 63 111 62
+341 4 2 11 1 110 62 119 111
+342 4 2 11 1 63 119 111 120
+343 4 2 11 1 120 64 112 63
+344 4 2 11 1 111 63 120 112
+345 4 2 11 1 64 120 112 97
+346 4 2 11 1 97 21 95 64
+347 4 2 11 1 112 64 97 95
+348 4 2 11 1 15 16 46 121
+349 4 2 11 1 73 69 121 15
+350 4 2 11 1 16 15 73 121
+351 4 2 11 1 69 73 121 122
+352 4 2 11 1 74 70 122 69
+353 4 2 11 1 73 69 74 122
+354 4 2 11 1 70 74 122 123
+355 4 2 11 1 75 71 123 70
+356 4 2 11 1 74 70 75 123
+357 4 2 11 1 71 75 123 124
+358 4 2 11 1 76 72 124 71
+359 4 2 11 1 75 71 76 124
+360 4 2 11 1 72 76 124 98
+361 4 2 11 1 24 25 98 72
+362 4 2 11 1 76 72 24 98
+363 4 2 11 1 13 14 47 125
+364 4 2 11 1 89 85 125 13
+365 4 2 11 1 14 13 89 125
+366 4 2 11 1 85 89 125 90
+367 4 2 11 1 90 86 126 125
+368 4 2 11 1 85 86 90 125
+369 4 2 11 1 86 90 126 91
+370 4 2 11 1 91 87 127 126
+371 4 2 11 1 86 87 91 126
+372 4 2 11 1 87 91 127 128
+373 4 2 11 1 92 88 128 87
+374 4 2 11 1 91 87 92 128
+375 4 2 11 1 88 92 128 99
+376 4 2 11 1 26 27 99 88
+377 4 2 11 1 92 88 26 99
+378 4 2 11 1 41 42 47 105
+379 4 2 11 1 105 101 125 47
+380 4 2 11 1 41 101 105 47
+381 4 2 11 1 101 105 125 102
+382 4 2 11 1 106 102 126 125
+383 4 2 11 1 105 102 106 125
+384 4 2 11 1 102 106 126 103
+385 4 2 11 1 107 103 127 126
+386 4 2 11 1 106 103 107 126
+387 4 2 11 1 103 107 127 128
+388 4 2 11 1 108 104 128 107
+389 4 2 11 1 107 104 128 103
+390 4 2 11 1 104 108 128 99
+391 4 2 11 1 94 93 99 108
+392 4 2 11 1 108 93 99 104
+393 4 2 11 1 41 46 43 101
+394 4 2 11 1 121 101 109 43
+395 4 2 11 1 46 101 121 43
+396 4 2 11 1 101 121 109 122
+397 4 2 11 1 122 102 110 109
+398 4 2 11 1 101 102 122 109
+399 4 2 11 1 102 122 110 103
+400 4 2 11 1 123 103 111 110
+401 4 2 11 1 122 103 123 110
+402 4 2 11 1 103 123 111 124
+403 4 2 11 1 124 104 112 111
+404 4 2 11 1 103 104 124 111
+405 4 2 11 1 104 124 112 93
+406 4 2 11 1 98 93 95 112
+407 4 2 11 1 124 93 98 112
+408 4 2 11 1 41 43 45 101
+409 4 2 11 1 109 101 117 43
+410 4 2 11 1 43 101 117 45
+411 4 2 11 1 101 109 117 102
+412 4 2 11 1 110 102 118 109
+413 4 2 11 1 109 102 118 117
+414 4 2 11 1 102 110 118 103
+415 4 2 11 1 111 103 119 110
+416 4 2 11 1 110 103 119 118
+417 4 2 11 1 103 111 119 120
+418 4 2 11 1 112 104 120 111
+419 4 2 11 1 111 104 120 103
+420 4 2 11 1 104 112 120 93
+421 4 2 11 1 95 93 97 112
+422 4 2 11 1 112 93 97 120
+423 4 2 11 1 15 46 44 113
+424 4 2 11 1 121 69 113 15
+425 4 2 11 1 46 15 121 113
+426 4 2 11 1 69 121 113 122
+427 4 2 11 1 122 70 114 69
+428 4 2 11 1 113 69 122 114
+429 4 2 11 1 70 122 114 123
+430 4 2 11 1 123 71 115 70
+431 4 2 11 1 114 70 123 115
+432 4 2 11 1 71 123 115 124
+433 4 2 11 1 124 72 116 71
+434 4 2 11 1 115 71 124 116
+435 4 2 11 1 72 124 116 98
+436 4 2 11 1 98 25 96 72
+437 4 2 11 1 116 72 98 96
+438 4 2 11 1 14 44 47 113
+439 4 2 11 1 113 89 125 14
+440 4 2 11 1 47 14 113 125
+441 4 2 11 1 89 113 125 90
+442 4 2 11 1 114 90 126 113
+443 4 2 11 1 113 90 126 125
+444 4 2 11 1 90 114 126 91
+445 4 2 11 1 115 91 127 114
+446 4 2 11 1 114 91 127 126
+447 4 2 11 1 91 115 127 128
+448 4 2 11 1 116 92 128 115
+449 4 2 11 1 115 92 128 91
+450 4 2 11 1 92 116 128 99
+451 4 2 11 1 96 26 99 116
+452 4 2 11 1 116 26 99 92
+453 4 2 11 1 8 42 48 105
+454 4 2 11 1 105 53 129 8
+455 4 2 11 1 48 8 105 129
+456 4 2 11 1 53 105 129 130
+457 4 2 11 1 106 54 130 53
+458 4 2 11 1 105 53 106 130
+459 4 2 11 1 54 106 130 107
+460 4 2 11 1 107 55 131 54
+461 4 2 11 1 130 54 107 131
+462 4 2 11 1 55 107 131 132
+463 4 2 11 1 108 56 132 55
+464 4 2 11 1 107 55 108 132
+465 4 2 11 1 56 108 132 100
+466 4 2 11 1 94 19 100 56
+467 4 2 11 1 108 56 94 100
+468 4 2 11 1 9 45 10 117
+469 4 2 11 1 117 57 61 9
+470 4 2 11 1 10 9 117 61
+471 4 2 11 1 57 117 61 118
+472 4 2 11 1 118 58 62 57
+473 4 2 11 1 61 57 118 62
+474 4 2 11 1 58 118 62 119
+475 4 2 11 1 119 59 63 58
+476 4 2 11 1 62 58 119 63
+477 4 2 11 1 59 119 63 120
+478 4 2 11 1 120 60 64 63
+479 4 2 11 1 59 60 120 63
+480 4 2 11 1 60 120 64 97
+481 4 2 11 1 97 20 21 64
+482 4 2 11 1 60 20 97 64
+483 4 2 11 1 41 45 48 101
+484 4 2 11 1 117 101 129 48
+485 4 2 11 1 45 101 117 48
+486 4 2 11 1 101 117 129 102
+487 4 2 11 1 118 102 130 117
+488 4 2 11 1 117 102 130 129
+489 4 2 11 1 102 118 130 103
+490 4 2 11 1 119 103 131 130
+491 4 2 11 1 118 103 119 130
+492 4 2 11 1 103 119 131 120
+493 4 2 11 1 120 104 132 131
+494 4 2 11 1 103 104 120 131
+495 4 2 11 1 104 120 132 93
+496 4 2 11 1 97 93 100 120
+497 4 2 11 1 120 93 100 132
+498 4 2 11 1 8 48 9 129
+499 4 2 11 1 129 53 57 8
+500 4 2 11 1 9 8 129 57
+501 4 2 11 1 53 129 57 130
+502 4 2 11 1 130 54 58 53
+503 4 2 11 1 57 53 130 58
+504 4 2 11 1 54 130 58 131
+505 4 2 11 1 131 55 59 54
+506 4 2 11 1 58 54 131 59
+507 4 2 11 1 55 131 59 132
+508 4 2 11 1 132 56 60 55
+509 4 2 11 1 59 55 132 60
+510 4 2 11 1 56 132 60 100
+511 4 2 11 1 100 19 20 56
+512 4 2 11 1 60 56 100 20
+513 4 2 11 1 41 44 46 113
+514 4 2 11 1 113 101 121 46
+515 4 2 11 1 41 101 113 46
+516 4 2 11 1 101 113 121 122
+517 4 2 11 1 114 102 122 113
+518 4 2 11 1 113 102 122 101
+519 4 2 11 1 102 114 122 103
+520 4 2 11 1 115 103 123 114
+521 4 2 11 1 114 103 123 122
+522 4 2 11 1 103 115 123 124
+523 4 2 11 1 116 104 124 115
+524 4 2 11 1 115 104 124 103
+525 4 2 11 1 104 116 124 93
+526 4 2 11 1 96 93 98 116
+527 4 2 11 1 116 93 98 124
+528 4 2 11 1 41 47 44 113
+529 4 2 11 1 125 101 113 47
+530 4 2 11 1 47 101 113 41
+531 4 2 11 1 101 125 113 102
+532 4 2 11 1 126 102 114 113
+533 4 2 11 1 125 102 126 113
+534 4 2 11 1 102 126 114 103
+535 4 2 11 1 127 103 115 114
+536 4 2 11 1 126 103 127 114
+537 4 2 11 1 103 127 115 128
+538 4 2 11 1 128 104 116 115
+539 4 2 11 1 103 104 128 115
+540 4 2 11 1 104 128 116 99
+541 4 2 11 1 99 93 96 116
+542 4 2 11 1 104 93 99 116
+543 4 2 11 1 9 48 45 117
+544 4 2 11 1 129 57 117 9
+545 4 2 11 1 48 9 129 117
+546 4 2 11 1 57 129 117 130
+547 4 2 11 1 130 58 118 57
+548 4 2 11 1 117 57 130 118
+549 4 2 11 1 58 130 118 119
+550 4 2 11 1 131 59 119 58
+551 4 2 11 1 130 58 131 119
+552 4 2 11 1 59 131 119 120
+553 4 2 11 1 132 60 120 59
+554 4 2 11 1 131 59 132 120
+555 4 2 11 1 60 132 120 100
+556 4 2 11 1 100 20 97 60
+557 4 2 11 1 120 60 100 97
+$EndElements
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/cylinder.py b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/cylinder.py
new file mode 100644
index 0000000000000000000000000000000000000000..9dc9a44b5a73bdcf953909b2a7cdb280e9747902
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/cylinder.py
@@ -0,0 +1,83 @@
+#coding-Utf-8-*-
+
+from gmshpy import *
+from dG3Dpy import*
+from math import*
+
+#script to launch PBC problem with a python script
+
+# material law
+lawnum = 11 # unique number of law
+law1 = VEVPUMATDG3DMaterialLaw(lawnum,'material.i02');
+#lawPert1 = dG3DMaterialLawWithTangentByPerturbation(law1, 1e-8)
+# geometry
+meshfile="cylinder.msh" # name of mesh file
+
+fullDg = 0 #O = CG, 1 = DG
+space1 = 0 # function space (Lagrange=0)
+beta1 = 100
+# creation of part Domain
+nfield = 11 # number of the field (physical number of surface)
+
+myfield1 = dG3DDomain(1000,nfield,space1,lawnum,fullDg,3)
+#myfield1.matrixByPerturbation(1,1,1,1e-8)
+myfield1.stabilityParameters(beta1)
+
+# solver
+sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
+soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 100 # number of step (used only if soltype=1)
+ftime =1. # Final time (used only if soltype=1)
+tol=1.e-8 # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+
+
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(myfield1)
+mysolver.addMaterialLaw(law1)
+#mysolver.addMaterialLaw(lawPert1)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol)
+
+mysolver.displacementBC('Face',1,2,0.)
+mysolver.displacementBC('Face',3,1,0.)
+mysolver.displacementBC('Face',4,0,0.)
+
+mysolver.displacementBC('Face',1,0,0)
+mysolver.displacementBC('Face',1,1,0)
+mysolver.displacementBC('Face',2,2,1.)
+
+# build view
+mysolver.internalPointBuildView("Strain_xx",IPField.STRAIN_XX, 1, 1);
+mysolver.internalPointBuildView("Strain_yy",IPField.STRAIN_YY, 1, 1);
+mysolver.internalPointBuildView("Strain_zz",IPField.STRAIN_ZZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xy",IPField.STRAIN_XY, 1, 1);
+mysolver.internalPointBuildView("Strain_yz",IPField.STRAIN_YZ, 1, 1);
+mysolver.internalPointBuildView("Strain_xz",IPField.STRAIN_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1);
+mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1);
+mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1);
+mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1);
+mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1);
+mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1);
+mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
+mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
+mysolver.internalPointBuildView("pression",IPField.PRESSION,1,1)
+mysolver.internalPointBuildView("plastic possion ratio",IPField.PLASTIC_POISSON_RATIO,1,1)
+
+#archiving
+mysolver.archivingForceOnPhysicalGroup('Face',1,0)
+mysolver.archivingForceOnPhysicalGroup('Face',1,1)
+mysolver.archivingForceOnPhysicalGroup('Face',1,2)
+mysolver.archivingNodeDisplacement(5,0,1)
+mysolver.archivingNodeDisplacement(5,1,1)
+mysolver.archivingNodeDisplacement(5,2,1)
+# solve
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-5.246765e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 1, 2),1.e-3)
diff --git a/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/material.i02 b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/material.i02
new file mode 100644
index 0000000000000000000000000000000000000000..3d531a75f52b068058bab1c442bccac54c23abc2
--- /dev/null
+++ b/dG3D/benchmarks/powerYieldViscoElastoPlastic/fortran/material.i02
@@ -0,0 +1,52 @@
+108 number of internal variables
+50 number of properties to be read
+11 Power Series Strain Approximation
+2250.0 Bulk Modulus
+1038.4615384615383 Shear Modulus
+3.5 Yield Exponent
+0.452 Plastic Poisson Ratio
+30000.0 Viscoplastic coefficient
+0.21 Viscoplasic exponent
+100.0 Initial Yield Limit - Compression
+300. Isotropic Hardening Parameter - Compression
+0.0 Isotropic Hardening Parameter - Compression
+10. Isotropic Hardening Parameter - Compression
+78.0 Initial Yield Limit - Tension
+300. Isotropic Hardening Parameter - Tension
+0. Isotropic Hardening Parameter - Tension
+10. Isotropic Hardening Parameter - Tension
+1e2 Kinematic Hardening Parameter
+2e3 Kinematic Hardening Parameter
+0.0 Kinematic Hardening Parameter
+1150.0 Bulk Modulus - Maxwell branch 1
+991.67 Bulk Modulus - Maxwell branch 2
+162.5 Bulk Modulus - Maxwell branch 3
+123.33 Bulk Modulus - Maxwell branch 4
+0.0 Bulk Modulus - Maxwell branch 5
+0.0 Bulk Modulus - Maxwell branch 6
+0.0 Bulk Modulus - Maxwell branch 7
+0.0 Bulk Modulus - Maxwell branch 8
+7202.0 Retardation Time Volumetric - Maxwell Branch 1
+71.6 Retardation Time Volumetric - Maxwell Branch 2
+0.73 Retardation Time Volumetric - Maxwell Branch 3
+0.073 Retardation Time Volumetric - Maxwell Branch 4
+1.0 Retardation Time Volumetric - Maxwell Branch 5
+1.0 Retardation Time Volumetric - Maxwell Branch 6
+1.0 Retardation Time Volumetric - Maxwell Branch 7
+1.0 Retardation Time Volumetric - Maxwell Branch 8
+530.77 Shear Modulus- Maxwell branch 1
+457.7 Shear Modulus- Maxwell branch 2
+75. Shear Modulus- Maxwell branch 3
+56.92 Shear Modulus- Maxwell branch 4
+0.0 Shear Modulus- Maxwell branch 5
+0.0 Shear Modulus- Maxwell branch 6
+0.0 Shear Modulus- Maxwell branch 7
+0.0 Shear Modulus- Maxwell branch 8
+7202.0 Retardation Time Deviatoric - Maxwell Branch 1
+71.6 Retardation Time Deviatoric - Maxwell Branch 2
+0.73 Retardation Time Deviatoric - Maxwell Branch 3
+0.073 Retardation Time Deviatoric - Maxwell Branch 4
+1.0 Retardation Time Deviatoric - Maxwell Branch 5
+1.0 Retardation Time Deviatoric - Maxwell Branch 6
+1.0 Retardation Time Deviatoric - Maxwell Branch 7
+1.0 Retardation Time Deviatoric - Maxwell Branch 8
diff --git a/dG3D/src/dG3DDeepMaterialNetworks.cpp b/dG3D/src/dG3DDeepMaterialNetworks.cpp
index e01fca6ceafed8277647ba2fc96777824b568d4f..16be25cffe834d852c4dc8bffcfbcb25f89985fe 100644
--- a/dG3D/src/dG3DDeepMaterialNetworks.cpp
+++ b/dG3D/src/dG3DDeepMaterialNetworks.cpp
@@ -408,6 +408,23 @@ void dG3DDeepMaterialNetwork::getDIPCompDFLocal(const IPStateBase* ips, IPField:
val = ipv->get(IPField::VOLUMETRIC_STRAIN);
STensorOperation::unity(DvalDF);
}
+ else if (comp == IPField::PLASTIC_ENERGY)
+ {
+ const ipFiniteStrain* ipv = dynamic_cast<const ipFiniteStrain*>(ips->getState(IPStateBase::current));
+ val = ipv->plasticEnergy();
+ const IPVariable* ipvData = ipv->getInternalData();
+ const IPJ2linear* ipj2 = dynamic_cast<const IPJ2linear*>(ipvData);
+ if (ipj2 != NULL)
+ {
+ DvalDF = ipj2->getConstRefToDPlasticEnergyDF();
+ }
+ else
+ {
+ Msg::Error("IPJ2linear must be used");
+ Msg::Exit(0);
+ }
+
+ }
else
{
Msg::Error("dG3DDeepMaterialNetwork::getDIPCompDFLocal is not defined for %s",IPField::ToString(comp));
diff --git a/dG3D/src/dG3DFunctionSpace.h b/dG3D/src/dG3DFunctionSpace.h
index 564648cd11a095c031cc61f16f444e1f686ca2ed..167e0ad924da3e40f185d0ae757e605eea4c1ee0 100644
--- a/dG3D/src/dG3DFunctionSpace.h
+++ b/dG3D/src/dG3DFunctionSpace.h
@@ -57,7 +57,11 @@ class g3DLagrangeFunctionSpace : public ThreeDLagrangeFunctionSpace{
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i){
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
}
}
}
@@ -72,7 +76,11 @@ class g3DLagrangeFunctionSpace : public ThreeDLagrangeFunctionSpace{
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i){
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
}
}
}
@@ -216,8 +224,12 @@ class g3DhoDGFunctionSpace : public ThreeDLagrangeFunctionSpace{
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i)
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
- }
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
+ }
}
else
#endif // HAVE_MPI
@@ -230,7 +242,11 @@ class g3DhoDGFunctionSpace : public ThreeDLagrangeFunctionSpace{
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i)
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
}
}
}
@@ -367,10 +383,15 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpace : public g3DLagrangeFun
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i)
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
}
}
else{ // match is not ensured --> generate keys for each partitions
+ #if defined(HAVE_CG_MPI_INTERFACE)
for(int p=0;p<=Msg::GetCommSize();p++){
for (int j=0;j<_ncomp;++j)
{
@@ -382,7 +403,19 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpace : public g3DLagrangeFun
for (int i=0;i<nk;++i)
keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],p)));
}
- }
+ }
+ #else
+ for (int j=0;j<_ncomp;++j)
+ {
+ int nk=ele->getNumVertices(); // return the number of vertices
+ if (withPrimaryShapeFunction(comp[j]))
+ {
+ nk = ele->getNumPrimaryVertices();
+ }
+ for (int i=0;i<nk;++i)
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ }
+ #endif
}
}
else
@@ -398,19 +431,11 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpace : public g3DLagrangeFun
}
for (int i=0;i<nk;++i)
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
-/* if(ele->getPartition()==0) // Add the keys corresponding to this partition
- {
- for (int j=0;j<_ncomp;++j)
- for (int i=0;i<nk;++i)
- keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],Msg::GetCommRank()+1)));
- // and fix all the other partitions
- for(int p=0;p<=Msg::GetCommSize();p++)
- for (int j=0;j<_ncomp;++j)
- for (int i=0;i<nk;++i)
- keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],p)));
- }*/
-
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
}
}
}
@@ -475,7 +500,12 @@ class g3DNeumannBoundaryConditionLagrangeFunctionSpace : public g3DLagrangeFunct
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i)
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],p)));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
+
}
}
};
@@ -522,10 +552,15 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpaceGhost : public g3DDirich
nk = ele->getNumPrimaryVertices();
}
for (int i=0;i<nk;++i)
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],ele->getPartition())));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ #endif
}
}
else{ // match is not ensured --> generate keys for each partitions
+ #if defined(HAVE_CG_MPI_INTERFACE)
for(int p=0;p<=Msg::GetCommSize();p++)
{
for (int j=0;j<_ncomp;++j)
@@ -539,6 +574,19 @@ class g3DDirichletBoundaryConditionLagrangeFunctionSpaceGhost : public g3DDirich
keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j],p)));
}
}
+ #else
+ for (int j=0;j<_ncomp;++j)
+ {
+ int nk=ele->getNumVertices(); // return the number of vertices
+ if (withPrimaryShapeFunction(this->comp[j]))
+ {
+ nk = ele->getNumPrimaryVertices();
+ }
+ for (int i=0;i<nk;++i)
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),-dG3DDof3IntType::createTypeWithThreeInts(comp[j],ifield[j])));
+ }
+
+ #endif
}
}
};
diff --git a/dG3D/src/dG3DHierarchicalFunctionSpace.cpp b/dG3D/src/dG3DHierarchicalFunctionSpace.cpp
index 247175ff1c94455a4f1927a9c5bcced0e3231e40..779d5f3ee69f6d11536fd8b29d93c1effded0659 100644
--- a/dG3D/src/dG3DHierarchicalFunctionSpace.cpp
+++ b/dG3D/src/dG3DHierarchicalFunctionSpace.cpp
@@ -162,7 +162,11 @@ void g3DHierarchicalFunctionSpace::getKeysOnElement_PerProc(MElement *ele, std::
{
for (int i=0; i< nv; i++)
{
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getVertex(i)->getNum(),ghost*dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],p,0)));
+ #else
+ keys.push_back(Dof(ele->getVertex(i)->getNum(),dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],0,0)));
+ #endif
}
}
// edges
@@ -176,7 +180,12 @@ void g3DHierarchicalFunctionSpace::getKeysOnElement_PerProc(MElement *ele, std::
Msg::Error("edge %d in element %d has not been numerated in %s",i,ele->getNum(),getFunctionSpaceName().c_str());
};
int entity = allEdges[ide];
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(entity,ghost*dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],p,1+k)));
+ #else
+ keys.push_back(Dof(entity,dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],0,1+k)));
+ #endif
+
}
}
// trifaces
@@ -192,7 +201,11 @@ void g3DHierarchicalFunctionSpace::getKeysOnElement_PerProc(MElement *ele, std::
Msg::Error("face %d in element %d has not been numerated in %s",i,ele->getNum(),getFunctionSpaceName().c_str());
};
int entity = allFaces[ifa];
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(entity,ghost*dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],p,11+k)));
+ #else
+ keys.push_back(Dof(entity,dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],0,11+k)));
+ #endif
}
}
}
@@ -209,14 +222,22 @@ void g3DHierarchicalFunctionSpace::getKeysOnElement_PerProc(MElement *ele, std::
Msg::Error("face %d in element %d has not been numerated in %s",i,ele->getNum(),getFunctionSpaceName().c_str());
};
int entity = allFaces[ifa];
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(entity,ghost*dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],p,21+k)));
+ #else
+ keys.push_back(Dof(entity,dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],0,21+k)));
+ #endif
}
}
}
// volumes
for (int k=0; k< nbFunctionPerDim[j][4]; k++)
{
+ #if defined(HAVE_CG_MPI_INTERFACE)
keys.push_back(Dof(ele->getNum(),ghost*dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],p,31+k)));
+ #else
+ keys.push_back(Dof(ele->getNum(),dG3DDof4IntType::createTypeWithFourInts(comp[j],ifield[j],0,31+k)));
+ #endif
}
}
};
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 3b86f51188c6bbcfe5c72ed06beeb48e70422725..9bb8fff93510130c1ef56e3152685b1d852d9223 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -4262,6 +4262,73 @@ void crystalPlasticityDG3DIPVariable::restart()
//
+//
+VEVPUMATDG3DIPVariable::VEVPUMATDG3DIPVariable(int _nsdv, double size, const bool createBodyForceHO, const bool oninter):
+ dG3DIPVariable(createBodyForceHO, oninter,0)
+{
+ _ipv = new IPVEVPUMAT(_nsdv,size);
+}
+
+VEVPUMATDG3DIPVariable::VEVPUMATDG3DIPVariable(const VEVPUMATDG3DIPVariable &source) :
+ dG3DIPVariable(source)
+{
+ _ipv = NULL;
+ if (source.getIPVEVPUMAT() != NULL){
+ _ipv = new IPVEVPUMAT((int)source.getIPVEVPUMAT()->getNsdv(),(double)source.getIPVEVPUMAT()->elementSize());
+ _ipv->operator= (*dynamic_cast<const IPVariable*>(source.getIPVEVPUMAT()));
+ }
+
+}
+
+VEVPUMATDG3DIPVariable& VEVPUMATDG3DIPVariable::operator=(const IPVariable &source)
+{
+ dG3DIPVariable::operator=(source);
+ const VEVPUMATDG3DIPVariable* src = dynamic_cast<const VEVPUMATDG3DIPVariable*>(&source);
+ if(src != NULL)
+ {
+ if (src->getIPVEVPUMAT() != NULL){
+ if (_ipv != NULL){
+ _ipv->operator=(*dynamic_cast<const IPVariable*>(src->getIPVEVPUMAT()));
+ }
+ else{
+ _ipv = dynamic_cast<IPVEVPUMAT*>(src->getIPVEVPUMAT()->clone());
+ }
+ }
+ else{
+ if (_ipv != NULL){
+ delete _ipv;
+ _ipv = NULL;
+ }
+ }
+ }
+ return *this;
+}
+double VEVPUMATDG3DIPVariable::get(const int i) const
+{
+ double val = dG3DIPVariable::get(i);
+ if (val == 0)
+ val = getIPVEVPUMAT()->get(i);
+ return val;
+}
+double VEVPUMATDG3DIPVariable::defoEnergy() const
+{
+ return getIPVEVPUMAT()->defoEnergy();
+}
+double VEVPUMATDG3DIPVariable::plasticEnergy() const
+{
+ return getIPVEVPUMAT()->plasticEnergy();
+}
+
+
+void VEVPUMATDG3DIPVariable::restart()
+{
+ dG3DIPVariable::restart();
+ if (_ipv != NULL)
+ restartManager::restart(_ipv);
+ return;
+}
+
+
//
gursonUMATDG3DIPVariable::gursonUMATDG3DIPVariable(int _nsdv, double size, const bool createBodyForceHO, const bool oninter):
dG3DIPVariable(createBodyForceHO, oninter,0)
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 4cb8ad93dc4722cab5e1a3b90d2e30f5529d8f11..a1dc41c349b50d929b97360677f2c62694f819eb 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -34,6 +34,7 @@
#include "ipEOS.h"
#include "ipCrystalPlasticity.h"
#include "ipGursonUMAT.h"
+#include "ipVEVPUMAT.h"
#include "ipViscoelastic.h"
#include "ipAnIsotropicElecTherMech.h"
#include "ipLinearElecTherMech.h"
@@ -3484,6 +3485,51 @@ class crystalPlasticityDG3DIPVariable : public dG3DIPVariable
virtual void restart();
};
+class VEVPUMATDG3DIPVariable : public dG3DIPVariable
+{
+
+ protected:
+ IPVEVPUMAT *_ipv;
+
+ public:
+ VEVPUMATDG3DIPVariable(int _nsdv, double size, const bool createBodyForceHO, const bool oninter);
+ VEVPUMATDG3DIPVariable(const VEVPUMATDG3DIPVariable &source);
+ virtual VEVPUMATDG3DIPVariable& operator=(const IPVariable &source);
+
+ /* specific function */
+ IPVEVPUMAT* getIPVEVPUMAT(){return _ipv;}
+ const IPVEVPUMAT* getIPVEVPUMAT() const{return _ipv;}
+ virtual ~VEVPUMATDG3DIPVariable()
+ {
+ if (_ipv) delete _ipv;
+ _ipv = NULL;
+ }
+
+ virtual void setLocation(const IPVariable::LOCATION loc) {
+ dG3DIPVariable::setLocation(loc);
+ if (_ipv)
+ _ipv->setLocation(loc);
+ };
+
+ virtual IPVariable* getInternalData() {return _ipv;};
+ virtual const IPVariable* getInternalData() const {return _ipv;};
+
+ virtual bool dissipationIsActive() const {return _ipv->dissipationIsActive();};
+
+ virtual void blockDissipation(const bool fl){if (_ipv) _ipv->blockDissipation(fl);};
+ virtual bool dissipationIsBlocked() const {
+ if (_ipv) return _ipv->dissipationIsBlocked();
+ else return false;
+ }
+
+ virtual double get(const int i) const;
+ virtual double defoEnergy() const;
+ virtual double plasticEnergy() const;
+
+ virtual IPVariable* clone() const {return new VEVPUMATDG3DIPVariable(*this);};
+ virtual void restart();
+};
+
class gursonUMATDG3DIPVariable : public dG3DIPVariable
{
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 62969d4c24692f0d039b2875970bcdcadb77d4ab..f2aa0c7900b3a1e0b08a5aba9d510f08f8e0b8bf 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -6771,6 +6771,154 @@ void crystalPlasticityDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver*
};
//
+//
+VEVPUMATDG3DMaterialLaw::VEVPUMATDG3DMaterialLaw(const int num, const char *propName) :
+ dG3DMaterialLaw(num,0.,true)
+{
+ _vevplaw = new mlawVEVPUMAT(num,propName);
+ double nu = _vevplaw->poissonRatio();
+ double mu = _vevplaw->shearModulus();
+ double E = 2.*mu*(1.+nu);
+ _rho = _vevplaw->getDensity();
+ fillElasticStiffness(E, nu, elasticStiffness);
+
+
+}
+
+VEVPUMATDG3DMaterialLaw::~VEVPUMATDG3DMaterialLaw()
+{
+ if (_vevplaw !=NULL) delete _vevplaw;
+ _vevplaw = NULL;
+};
+
+
+VEVPUMATDG3DMaterialLaw::VEVPUMATDG3DMaterialLaw(const VEVPUMATDG3DMaterialLaw &source) :
+ dG3DMaterialLaw(source)
+{
+ _vevplaw = NULL;
+ if (source._vevplaw != NULL){
+ _vevplaw = dynamic_cast<mlawVEVPUMAT*>(source._vevplaw->clone());
+ }
+
+}
+
+void VEVPUMATDG3DMaterialLaw::setTime(const double t,const double dtime){
+ dG3DMaterialLaw::setTime(t,dtime);
+ _vevplaw->setTime(t,dtime);
+ return;
+}
+
+materialLaw::matname VEVPUMATDG3DMaterialLaw::getType() const {return _vevplaw->getType();}
+
+bool VEVPUMATDG3DMaterialLaw::withEnergyDissipation() const {return _vevplaw->withEnergyDissipation();};
+
+void VEVPUMATDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ // check interface element
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele==NULL) inter=false;
+
+ int nsdv = _vevplaw->getNsdv();
+ double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
+
+ IPVariable* ipvi = new VEVPUMATDG3DIPVariable(nsdv,size,hasBodyForce,inter);
+ IPVariable* ipv1 = new VEVPUMATDG3DIPVariable(nsdv,size,hasBodyForce,inter);
+ IPVariable* ipv2 = new VEVPUMATDG3DIPVariable(nsdv,size,hasBodyForce,inter);
+
+ if(ips != NULL) delete ips;
+ ips = new IP3State(state_,ipvi,ipv1,ipv2);
+ _vevplaw->createIPState((static_cast <VEVPUMATDG3DIPVariable*> (ipvi))->getIPVEVPUMAT(),
+ (static_cast <VEVPUMATDG3DIPVariable*> (ipv1))->getIPVEVPUMAT(),
+ (static_cast <VEVPUMATDG3DIPVariable*> (ipv2))->getIPVEVPUMAT());
+
+}
+
+void VEVPUMATDG3DMaterialLaw::createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
+{
+ if(ipv !=NULL) delete ipv;
+ bool inter=true;
+ const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele);
+ if(iele == NULL) inter=false;
+ double size = 2.*(( const_cast<MElement*>(ele) )->getInnerRadius());
+
+ ipv = new VEVPUMATDG3DIPVariable(_vevplaw->getNsdv(),size,hasBodyForce,inter);
+
+ IPVEVPUMAT * ipvnl = static_cast <VEVPUMATDG3DIPVariable*>(ipv)->getIPVEVPUMAT();
+ _vevplaw->createIPVariable(ipvnl,hasBodyForce, ele,nbFF_,GP,gpt);
+}
+
+
+double VEVPUMATDG3DMaterialLaw::soundSpeed() const{return _vevplaw->soundSpeed();}
+
+void VEVPUMATDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const{
+ /* get ipvariable */
+ VEVPUMATDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+ const VEVPUMATDG3DIPVariable* ipvprev; //= static_cast <const nonLocalDamageDG3DIPVariable *> (ipvp);
+
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+
+ ipvcur = static_cast<VEVPUMATDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const VEVPUMATDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<VEVPUMATDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const VEVPUMATDG3DIPVariable*>(ipvp);
+ }
+ _vevplaw->checkInternalState(ipvcur->getInternalData(),ipvprev->getInternalData());
+
+};
+
+void VEVPUMATDG3DMaterialLaw::stress(IPVariable* ipv, const IPVariable* ipvp, const bool stiff, const bool checkfrac, const bool dTangent)
+{
+ /* get ipvariable */
+ VEVPUMATDG3DIPVariable* ipvcur; //= static_cast < nonLocalDamageDG3DIPVariable *> (ipv);
+ const VEVPUMATDG3DIPVariable* ipvprev; //= static_cast <const nonLocalDamageDG3DIPVariable *> (ipvp);
+
+ FractureCohesive3DIPVariable* ipvtmp = dynamic_cast<FractureCohesive3DIPVariable*>(ipv);
+ if(ipvtmp !=NULL)
+ {
+
+ ipvcur = static_cast<VEVPUMATDG3DIPVariable*>(ipvtmp->getIPvBulk());
+ const FractureCohesive3DIPVariable *ipvtmp2 = static_cast<const FractureCohesive3DIPVariable*>(ipvp);
+ ipvprev = static_cast<const VEVPUMATDG3DIPVariable*>(ipvtmp2->getIPvBulk());
+ }
+ else
+ {
+ ipvcur = static_cast<VEVPUMATDG3DIPVariable*>(ipv);
+ ipvprev = static_cast<const VEVPUMATDG3DIPVariable*>(ipvp);
+ }
+
+ /* strain xyz */
+ const STensor3& F1 = ipvcur->getConstRefToDeformationGradient();
+ const STensor3& F0 = ipvprev->getConstRefToDeformationGradient();
+
+ /* data for J2 law */
+ IPVEVPUMAT* q1 = ipvcur->getIPVEVPUMAT();
+ const IPVEVPUMAT* q0 = ipvprev->getIPVEVPUMAT();
+
+ /* compute stress */
+ _vevplaw->constitutive(F0,F1,ipvcur->getRefToFirstPiolaKirchhoffStress(),q0,q1,ipvcur->getRefToTangentModuli(),
+ stiff, NULL,dTangent,ipvcur->getPtrTodCmdFm());
+
+ ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+
+}
+
+double VEVPUMATDG3DMaterialLaw::scaleFactor() const{return _vevplaw->scaleFactor();};
+
+void VEVPUMATDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
+ dG3DMaterialLaw::setMacroSolver(sv);
+ if (_vevplaw != NULL){
+ _vevplaw->setMacroSolver(sv);
+ }
+};
+
+
//
gursonUMATDG3DMaterialLaw::gursonUMATDG3DMaterialLaw(const int num, const double temperature, const char *propName) :
dG3DMaterialLaw(num,0.,true)
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index a0479e3ede44bb98ea2b9d51633261229cb474c7..8ba7f7a811c68766aab277f03e5bb7f2cfd7fcee 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -29,6 +29,7 @@
#include "mlawEOS.h"
#include "mlawCrystalPlasticity.h"
#include "mlawGursonUMAT.h"
+#include "mlawVEVPUMAT.h"
#include "mlawViscoelastic.h"
#include "mlawLinearElecTherMech.h"
#include "mlawAnIsotropicElecTherMech.h"
@@ -2411,6 +2412,42 @@ class crystalPlasticityDG3DMaterialLaw : public dG3DMaterialLaw{
#endif
};
+class VEVPUMATDG3DMaterialLaw : public dG3DMaterialLaw{
+ protected:
+ #ifndef SWIG
+ mlawVEVPUMAT *_vevplaw;
+ #endif //SWIG
+
+ public:
+ VEVPUMATDG3DMaterialLaw(const int num, const char *propName);
+ #ifndef SWIG
+ VEVPUMATDG3DMaterialLaw(const VEVPUMATDG3DMaterialLaw &source);
+ virtual ~VEVPUMATDG3DMaterialLaw();
+
+ // set the time of _cplaw
+ virtual void setTime(const double t,const double dtime);
+ virtual materialLaw::matname getType() const;
+ virtual void createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_=NULL,const MElement *ele=NULL, const int nbFF_=0, const IntPt *GP=NULL, const int gpt = 0) const;
+ virtual void createIPVariable(IPVariable* &ipv, bool hasBodyForce, const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const;
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw){}
+ virtual double soundSpeed() const;// or change value ??
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvp) const;
+ virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
+ virtual double scaleFactor() const;
+ virtual materialLaw* clone() const{return new VEVPUMATDG3DMaterialLaw(*this);};
+ virtual bool withEnergyDissipation() const;
+ virtual void setMacroSolver(const nonLinearMechSolver* sv);
+ virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
+ {
+ return _vevplaw;
+ }
+ virtual materialLaw *getNonLinearSolverMaterialLaw()
+ {
+ return _vevplaw;
+ }
+ #endif
+};
+
class gursonUMATDG3DMaterialLaw : public dG3DMaterialLaw{
protected:
#ifndef SWIG