diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 299bcd77e153502a1352b8948963c05ae7e5eaec..922e0fdd50d563c98d63c07ee68fd0e32a3ea7f4 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -691,7 +691,10 @@ std::string IPField::ToString(const int i){
else if (i == TEMPERATURE) return "TEMPERATURE";
else if (i == THERMALFLUX_X) return "THERMALFLUX_X";
else if (i == THERMALFLUX_Y) return "THERMALFLUX_Y";
- else if (i == THERMALFLUX_Z) return "THERMALFLUX_Z";
+ else if (i == THERMALFLUX_Z) return "THERMALFLUX_Z";
+ else if (i == TEMP_GRAD_X) return "TEMP_GRAD_X";
+ else if (i == TEMP_GRAD_Y) return "TEMP_GRAD_Y";
+ else if (i == TEMP_GRAD_Z) return "TEMP_GRAD_Z";
else if (i == THERMALSOURCE) return "THERMALSOURCE";
else if (i == VOLTAGE) return "VOLTAGE";
else if (i == ELECTRICALFLUX_X) return "ELECTRICALFLUX_X";
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index 2065f0f4a20c132d9d133decffd1f297fea0bbb4..c19e5143813beb01aab955ba95b145b295ea7bb4 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -260,7 +260,7 @@ class IPField : public elementsField {
ISO_YIELD,ISO_YIELD_TENSILE, ISO_YIELD_COMPRESSION, ISO_YIELD_SHEAR,KIN_YIELD,
MTX_SVM,MTX_SIG_XX,MTX_SIG_YY,MTX_SIG_ZZ,MTX_SIG_XY,MTX_SIG_YZ,MTX_SIG_XZ,
INC_SVM,INC_SIG_XX,INC_SIG_YY,INC_SIG_ZZ,INC_SIG_XY,INC_SIG_YZ,INC_SIG_XZ,
- TEMPERATURE,THERMALFLUX_X,THERMALFLUX_Y,THERMALFLUX_Z,THERMALSOURCE,
+ TEMPERATURE,THERMALFLUX_X,THERMALFLUX_Y,THERMALFLUX_Z,THERMALSOURCE,TEMP_GRAD_X,TEMP_GRAD_Y,TEMP_GRAD_Z,
VOLTAGE,ELECTRICALFLUX_X,ELECTRICALFLUX_Y,ELECTRICALFLUX_Z,ELECTRICALSOURCE,
LOCAL_POROSITY, NONLOCAL_POROSITY, CORRECTED_POROSITY, YIELD_POROSITY,
NUCLEATED_POROSITY_TOT, NUCLEATED_POROSITY_0, NUCLEATED_POROSITY_1, NUCLEATED_POROSITY_2,
diff --git a/NonLinearSolver/internalPoints/ipFiniteStrain.h b/NonLinearSolver/internalPoints/ipFiniteStrain.h
index b5dc338a7f95c3b9a5f6a507ab411a514f91db6e..a09d0cc1781e4ca43f3bd5cfc75819a37536b110 100644
--- a/NonLinearSolver/internalPoints/ipFiniteStrain.h
+++ b/NonLinearSolver/internalPoints/ipFiniteStrain.h
@@ -76,6 +76,59 @@ class ipFiniteStrain : public IPVariableMechanics{
virtual const SVector3 &getConstRefToJump() const=0;
virtual SVector3 &getRefToJump()=0;
+ // FLE
+ //extraDof
+ virtual int getNumConstitutiveExtraDofDiffusionVariable() const=0;
+ virtual double getConstRefToField(const int idex) const=0;
+ virtual double &getRefToField(const int idex)=0;
+ virtual const std::vector<SVector3> &getConstRefToGradField() const=0;
+ virtual std::vector<SVector3> &getRefToGradField()=0;
+ virtual const std::vector<SVector3> &getConstRefToFlux() const=0;
+ virtual std::vector<SVector3> &getRefToFlux()=0;
+
+ virtual const std::vector<STensor3> &getConstRefTodPdField() const=0;
+ virtual std::vector<STensor3> &getRefTodPdField()=0;
+ virtual const std::vector<STensor33> &getConstRefTodPdGradField() const=0;
+ virtual std::vector<STensor33> &getRefTodPdGradField()=0;
+ virtual const std::vector<std::vector<STensor3> > &getConstRefTodFluxdGradField() const=0;
+ virtual std::vector<std::vector<STensor3> > &getRefTodFluxdGradField()=0;
+ virtual const std::vector<std::vector<SVector3> > &getConstRefTodFluxdField() const=0;
+ virtual std::vector<std::vector<SVector3> > &getRefTodFluxdField()=0;
+ virtual const std::vector<STensor33> &getConstRefTodFluxdF() const=0;
+ virtual std::vector<STensor33> &getRefTodFluxdF()=0;
+
+ virtual const std::vector<std::vector<const STensor3*> > &getConstRefToLinearK() const =0;
+ virtual void setConstRefToLinearK(const STensor3 &eT, int i, int j) =0;
+ virtual const std::vector<std::vector<std::vector<STensor3> > > &getConstRefTodLinearKdField() const=0;
+ virtual std::vector<std::vector<std::vector<STensor3> > > &getRefTodLinearKdField()=0;
+ virtual const std::vector< std::vector<STensor43> > &getConstRefTodLinearKdF() const=0;
+ virtual std::vector< std::vector<STensor43> > &getRefTodLinearKdF()=0;
+
+ virtual const std::vector< std::vector<STensor43> > &getConstRefTodFluxdGradFielddF() const=0;
+ virtual std::vector<std::vector<STensor43> > &getRefTodFluxdGradFielddF()=0;
+ virtual const std::vector<std::vector<std::vector<STensor3> > > &getConstRefTodFluxdGradFielddField() const=0;
+ virtual std::vector<std::vector<std::vector<STensor3> > > &getRefTodFluxdGradFielddField()=0;
+ virtual const fullVector<double> &getConstRefToFieldSource() const=0;
+ virtual fullVector<double> &getRefToFieldSource()=0;
+ virtual const fullMatrix<double> &getConstRefTodFieldSourcedField() const=0;
+ virtual fullMatrix<double> &getRefTodFieldSourcedField()=0;
+ virtual const std::vector<std::vector<SVector3> > &getConstRefTodFieldSourcedGradField() const=0;
+ virtual std::vector<std::vector<SVector3> > &getRefTodFieldSourcedGradField()=0;
+ virtual const std::vector<STensor3> &getConstRefTodFieldSourcedF() const=0;
+ virtual std::vector<STensor3> &getRefTodFieldSourcedF()=0;
+ virtual const fullVector<double> &getConstRefToMechanicalSource() const=0;
+ virtual fullVector<double> &getRefToMechanicalSource()=0;
+ virtual const fullMatrix<double> &getConstRefTodMechanicalSourcedField() const=0;
+ virtual fullMatrix<double> &getRefTodMechanicalSourcedField()=0;
+ virtual const std::vector<std::vector<SVector3> > &getConstRefTodMechanicalSourcedGradField() const=0;
+ virtual std::vector<std::vector<SVector3> > &getRefTodMechanicalSourcedGradField()=0;
+ virtual const std::vector<STensor3> &getConstRefTodMechanicalSourcedF() const=0;
+ virtual std::vector<STensor3> &getRefTodMechanicalSourcedF()=0;
+
+ virtual const fullVector<double> &getConstRefToExtraDofFieldCapacityPerUnitField() const=0;
+ virtual fullVector<double> &getRefToExtraDofFieldCapacityPerUnitField()=0;
+ // FLE
+
virtual bool hasBodyForceForHO() const =0;
virtual const STensor33 &getConstRefToGM()const=0;
virtual STensor33 &getRefToGM()=0;
diff --git a/NonLinearSolver/materialLaw/mlaw.h b/NonLinearSolver/materialLaw/mlaw.h
index 4329b2659e603d7caf858e2f615033f59751a692..c9fa30748ee0daca3cd9b3e2ef7aab0afa9e36f3 100644
--- a/NonLinearSolver/materialLaw/mlaw.h
+++ b/NonLinearSolver/materialLaw/mlaw.h
@@ -33,7 +33,7 @@ class materialLaw{
hyperelastic, powerYieldLaw, powerYieldLawWithFailure,nonLocalDamagePowerYieldHyper,
localDamagePowerYieldHyperWithFailure,nonLocalDamagePowerYieldHyperWithFailure,ElecSMP,
ThermalConducter,AnIsotropicTherMech, localDamageJ2Hyper,linearElastic,nonLocalDamageGursonThermoMechanics,
- localDamageJ2SmallStrain,nonLocalDamageJ2SmallStrain,cluster,tfa,ANN, DMN, torchANN, NonlocalDamageTorchANN, StochDMN, LinearElecMagTherMech, LinearElecMagInductor, hyperviscoelastic, GenericResidualStrain,
+ localDamageJ2SmallStrain,nonLocalDamageJ2SmallStrain,cluster,tfa,ANN, DMN, torchANN, NonlocalDamageTorchANN, StochDMN, StochTMDMN, LinearElecMagTherMech, LinearElecMagInductor, hyperviscoelastic, GenericResidualStrain,
GenericThermoMechanics, ElecMagGenericThermoMechanics, ElecMagInductor,
nonlineartvm,nonlinearTVE,nonlinearTVP,vevpmfh, VEVPUMAT, IMDEACPUMAT, Hill48,nonlinearTVEnonlinearTVP,nonlinearTVEnonlinearTVP2,
TVEVPwithFailure,nonLocalDamageTVEVP,localDamageTVEVPWithFailure,nonLocalDamageTVEVPWithFailure,
diff --git a/NonLinearSolver/modelReduction/NetworkInteraction.cpp b/NonLinearSolver/modelReduction/NetworkInteraction.cpp
index dcfc63fc2f0076f0ea42cc60d3a954f4263ee1a5..755d286b09f1fd4a5f6044e566782791d596aad9 100644
--- a/NonLinearSolver/modelReduction/NetworkInteraction.cpp
+++ b/NonLinearSolver/modelReduction/NetworkInteraction.cpp
@@ -408,6 +408,148 @@ void InteractionMechanism::getDInducedStrainDCoefficient(int nodeId, STensor3& D
};
+// FLE
+const double& InteractionMechanismExtraDOF::getConstRefToIncompatibleExtraDOF(const IPStateBase::whichState ws) const
+{
+ if (ws == IPStateBase::current)
+ {
+ return tcur;
+ }
+ else if (ws == IPStateBase::previous)
+ {
+ return tprev;
+ }
+ else if (ws == IPStateBase::initial)
+ {
+ return tinit;
+ }
+ else
+ {
+ Msg::Error("state %s has not been defined in InteractionMechanism::getConstRefToIncompatibleExtraDOF",ws);
+ Msg::Exit(0);
+ }
+};
+
+double& InteractionMechanismExtraDOF::getRefToIncompatibleExtraDOF(const IPStateBase::whichState ws)
+{
+ if (ws == IPStateBase::current)
+ {
+ return tcur;
+ }
+ else if (ws == IPStateBase::previous)
+ {
+ return tprev;
+ }
+ else if (ws == IPStateBase::initial)
+ {
+ return tinit;
+ }
+ else
+ {
+ Msg::Error("state %s has not been defined in InteractionMechanism::getRefToIncompatibleExtraDOF",ws);
+ Msg::Exit(0);
+ }
+};
+
+void InteractionMechanismExtraDOF::resetIncompatibleExtraDOF()
+{
+ tcur = 0.; tprev = 0.; tinit = 0.;
+};
+
+void InteractionMechanismExtraDOF::nextStepExtraDOF()
+{
+ tprev = tcur;
+};
+void InteractionMechanismExtraDOF::resetToPreviousStepExtraDOF()
+{
+ tcur = tprev;
+};
+
+void InteractionMechanismExtraDOF::copyExtraDOF(const IPStateBase::whichState ws1, const IPStateBase::whichState ws2)
+{
+ getRefToIncompatibleExtraDOF(ws2) = getRefToIncompatibleExtraDOF(ws1);
+};
+
+void InteractionMechanismExtraDOF::getInducedExtraDOFGrad(int nodeId, SVector3& Hi) const
+{
+ std::map<int,double>::const_iterator itF = coefficients.find(nodeId);
+ if (itF != coefficients.end())
+ {
+ static SVector3 direction;
+ getDirection(direction);
+ double alpha = itF->second;
+ for (int i=0; i<3; i++)
+ {
+ Hi(i) = alpha*tcur*direction(i);
+ }
+ }
+ else
+ {
+ Hi = 0.;
+ }
+};
+
+void InteractionMechanismExtraDOF::getDInducedExtraDOFGradDa(int nodeId, SVector3& DHiDt) const
+{
+ std::map<int,double>::const_iterator itF = coefficients.find(nodeId);
+ if (itF != coefficients.end())
+ {
+ static SVector3 direction;
+ getDirection(direction);
+ double alpha = itF->second;
+ for (int r=0; r<3; r++)
+ {
+ DHiDt(r) = direction(r)*alpha;
+ }
+ }
+ else
+ {
+ DHiDt = 0.;
+ }
+};
+
+void InteractionMechanismExtraDOF::getDInducedExtraDOFGradDDirection(int nodeId, STensor3& DHiDdir) const
+{
+ std::map<int,double>::const_iterator itF = coefficients.find(nodeId);
+ if (itF != coefficients.end())
+ {
+ static STensor3 I(1.);
+ double alpha = itF->second;
+ for (int r=0; r<3; r++)
+ {
+ for (int c=0; c<3; c++)
+ {
+ DHiDdir(r,c) = alpha*tcur*I(r,c);
+ }
+ };
+ }
+ else
+ {
+ STensorOperation::zero(DHiDdir);
+ }
+};
+
+void InteractionMechanismExtraDOF::getDInducedExtraDOFGradDCoefficient(int nodeId, SVector3& DHiDalpha) const
+{
+ std::map<int,double>::const_iterator itF = coefficients.find(nodeId);
+ if (itF != coefficients.end())
+ {
+ static SVector3 direction;
+ getDirection(direction);
+ double alpha = itF->second;
+ for (int r=0; r<3; r++)
+ {
+ DHiDalpha(r) = tcur*direction(r);
+ };
+ }
+ else
+ {
+ DHiDalpha = 0.;
+ }
+};
+
+// FLE
+
NetworkInteraction::NetworkInteraction(int tag): _tag(tag)
{
diff --git a/NonLinearSolver/modelReduction/NetworkInteraction.h b/NonLinearSolver/modelReduction/NetworkInteraction.h
index a737379823fb847bcbd6243c98577667afe05fe9..309890b95c4231b1bdcb4ec847e4b1d18468f25e 100644
--- a/NonLinearSolver/modelReduction/NetworkInteraction.h
+++ b/NonLinearSolver/modelReduction/NetworkInteraction.h
@@ -78,6 +78,37 @@ class InteractionMechanism
#endif //SWIG
};
+// FLE
+class InteractionMechanismExtraDOF: public InteractionMechanism // for tempGrad
+{
+ protected:
+ #ifndef SWIG
+ double tcur, tprev, tinit; // incompatible temperature scalar
+ #endif //SWIG
+ public:
+ InteractionMechanismExtraDOF():InteractionMechanism(),tcur(),tprev(),tinit(){}
+
+ #ifndef SWIG
+ InteractionMechanismExtraDOF(const InteractionMechanismExtraDOF& src): InteractionMechanism(src),tcur(src.tcur),tprev(src.tprev),tinit(src.tinit){}
+ InteractionMechanismExtraDOF* clone() const {return new InteractionMechanismExtraDOF(*this);}
+ ~InteractionMechanismExtraDOF(){}
+
+ const double& getConstRefToIncompatibleExtraDOF(const IPStateBase::whichState ws) const;
+ double& getRefToIncompatibleExtraDOF(const IPStateBase::whichState ws);
+
+ void nextStepExtraDOF();
+ void resetIncompatibleExtraDOF();
+ void resetToPreviousStepExtraDOF();
+ void copyExtraDOF(const IPStateBase::whichState ws1, const IPStateBase::whichState ws2);
+ // induce tempGrad from interaction
+ void getInducedExtraDOFGrad(int nodeIndex, SVector3& Hi) const;
+ void getDInducedExtraDOFGradDa(int nodeIndex, SVector3& DHiDt) const;
+ void getDInducedExtraDOFGradDDirection(int nodeIndex, STensor3& DHiDdir) const;
+ void getDInducedExtraDOFGradDCoefficient(int nodeIndex, SVector3& DHiDalpha) const;
+ #endif //SWIG
+};
+// FLE
+
class NetworkInteraction
{
#ifndef SWIG
@@ -175,7 +206,7 @@ class NetworkInteraction
//
const STensor3& getStress(const MaterialNode* node) const;
const STensor3& getDeformationGradient(const MaterialNode* node) const;
- const STensor43& getTangent(const MaterialNode* node ) const;
+ const STensor43& getTangent(const MaterialNode* node ) const; // dPdF
void getStress(STensor3& P) const;
void getDeformationGradientPhase(int phaseIndex, STensor3& F) const;
//
@@ -209,6 +240,132 @@ class NetworkInteraction
#endif //SWIG
};
+/*
+// FLE
+class NetworkInteractionExtraDOF : public NetworkInteraction
+{
+ public:
+ NetworkInteractionExtraDOF();
+ NetworkInteractionExtraDOF(const NetworkInteraction& src);
+ NetworkInteractionExtraDOF* clone() const {return new NetworkInteractionExtraDOF(*this);};
+ ~NetworkInteractionExtraDOF();
+
+ int getTag() const;
+ void clearData();
+ void getInteractionFromTree(const Tree& T, bool standardize=true);
+
+ void resetIncompatibleExtraDOF();
+
+ void copyFrom(const NetworkInteraction& src);
+ void mergeNodesInteraction(InteractionMechanism& ie, const std::vector<int>& allNodes) const;
+ void mergeNodesAllInteractions(const std::vector<int>& allNodes);
+ void mergeNodes(int phaseIndex, // phase index
+ int numberNodes, // only interactions with numberNodes nodes
+ NetworkInteraction& newInteraction);
+
+ int getNumberOfMaterialNodes() const;
+ int getDimension() const;
+ double getTotalWeight() const;
+ double getTotalWeightPhase(int iphase) const;
+ void getTotalWeightAllPhases(std::map<int, double>& allPhase) const;
+ void getPhaseFractions(std::map<int, double>& allPhase) const;
+ double getPhaseFraction(int phase) const;
+ int getPhaseIndex(int nodeId) const;
+ double getWeight(const MaterialNode* node) const;
+ double getWeight(int nodeId) const;
+ void getWeights(const InteractionMechanism* im, std::vector<double>& allWeights) const;
+ double get(int comp, int eleVal) const;
+ double getPerPhase(int phaseIndex, int comp, int eleVal) const;
+
+ void saveDataToFile(std::string filename) const;
+ void loadDataFromFile(std::string filename);
+ void loadCohesiveDataFromFile(std::string filename);
+ void convertToLinearActivation(std::string outputFileName);
+ void createInteractionFromPairsData(std::string filename, int dim, std::string affunc="relu");
+ void createInteractionFullByDivison(int numPhases, int nodePerPhase, int numInteraction, int dim, std::string affunc="relu");
+ void initializeRandomFullInteraction(int numberMaterialNodes, int numberInteractions, int dim,
+ const fullVector<double>& phaseFraction, const std::string affunc="relu",
+ bool rand=true);
+
+ void createTreeBasedFullInteraction(int Nlayers, int nodePerPhaseFullInteraction,
+ int numFullInteraction, int dim,
+ const fullVector<double>& phaseFraction,
+ const std::string affunc="relu",
+ bool rand=true, bool triple=false);
+ void initializeRandomPointwiseInteraction(int numberMaterialNodes, int numPhases, int dim, bool rand=true, std::string affunc="relu");
+ void normalizeWeights(std::string newaffunc="linear");
+ void standardizeInteractionDirections();
+ void normalizeAllInteractionCoefficients();
+ bool reInitialize();
+ bool checkValidity(bool message=true) const;
+
+ #ifndef SWIG
+ void getMaterialNodesByMaterialIndex(int matIndex, std::vector<const MaterialNode*>& allNode) const;
+ void getMaterialNodesByMaterialIndex(int matIndex, std::vector<MaterialNode*>& allNode);
+ MaterialNodeContainer& getMaterialNodeContainer(){return _allMaterialNodes;}
+ const MaterialNodeContainer& getMaterialNodeContainer() const {return _allMaterialNodes;}
+ InteractionContainer& getInteractionContainer() {return _allInteractions;};
+ const InteractionContainer& getInteractionContainer() const {return _allInteractions;};
+ NodeMapContainer& getNodeMapContainer() {return _nodeInteractionMap;}
+ const NodeMapContainer& getNodeMapContainer() const {return _nodeInteractionMap;}
+
+ void getAllMaterialNodes(std::vector<const MaterialNode*>& allNode) const;
+ void getAllMaterialNodesForMat(int matNum, std::vector<const MaterialNode*>& allNode) const;
+ void getAllMaterialNodes(std::vector<MaterialNode*>& allNode);
+
+ MaterialNode* getMaterialNode(const int number);
+ const MaterialNode* getMaterialNode(const int number) const;
+ void getAllMaterialNodeIds(std::vector<int>& allIds) const;
+ void getAllMaterialNodeIdsForMat(int matNum, std::vector<int>& allIds) const;
+ void getNumberMaterialNodesPerMat(std::map<int, int>& nodePerMat) const;
+
+ const std::vector<const InteractionMechanism*>* getInteractionsForNode(const MaterialNode* node) const;
+
+ void getDHDExtraDOF(const MaterialNode* node, std::vector<SVector3>& DHiDt) const;
+ void getDHDNormal(const MaterialNode* node, std::vector<STensor3>& DHiDn) const;
+ void getDHDCoefficient(const MaterialNode* node, std::vector<SVector3>& DHiDalpha) const;
+
+ void getDfluxDIncompatibleExtraDOF(const MaterialNode* node, std::vector<SVector3>& DqiDt) const;
+ void getDfluxDNormal(const MaterialNode* node, std::vector<STensor3>& DqiDn) const;
+ void getDfluxDCoefficient(const MaterialNode* node, std::vector<SVector3>& DqiDalpha) const;
+ //
+ const STensor3& getFlux(const MaterialNode* node) const;
+ const STensor3& getExtraDOFGrad(const MaterialNode* node) const;
+ const STensor43& getTangent(const MaterialNode* node ) const; // dqDGradT // CHECK
+ void getFlux(SVector3& q) const;
+ void getExtraDOFGradientPhase(int phaseIndex, SVector3& H) const;
+ //
+ void nextStep();
+ void resetToPreviousStep();
+ void copy(const IPStateBase::whichState ws1, const IPStateBase::whichState ws2);
+ //
+ // linear term for im with respect to jm
+ void getDLinearTermDCoefficient(const InteractionMechanism* im, const InteractionMechanism* jm, fullMatrix<double>& mat) const;
+ // linear term for im with respect to jm
+ void getDLinearTermDNormalVars(const InteractionMechanism* im, const InteractionMechanism* jm, fullMatrix<double>& mat) const;
+ // with repsecto to material weights
+ void getDLinearTermDWeights(const InteractionMechanism* im, std::vector<int>& materialIds, fullMatrix<double>& DvecDWeights) const;
+ //
+ void computeLocalDeformationGradient(const STensor3& Fmacro);
+ // get force
+ void getLinearTerm(const InteractionMechanism* im, fullVector<double>& vec) const;
+ // get DforceDa
+ void getBilinearTerm(const InteractionMechanism* im, std::vector<const InteractionMechanism*>& others, fullMatrix<double>& mat) const;
+ // get DforceDFbar
+ void getTangentTerm(const InteractionMechanism* im, fullMatrix<double>& mat) const;
+ //
+ void getLinearInterfaceTerm(const MaterialNode* node, fullVector<double>& vec) const;
+ //
+ void getBiLinearInterfaceTerm(const MaterialNode* node, fullMatrix<double>& mat) const;
+ //
+ void getBilinearCrossTermBulkInterface(const InteractionMechanism* im, std::vector<const MaterialNode*>& materialNodes, fullMatrix<double>& mat) const;
+ //
+ void getBilinearCrossTermInterfaceBulk(const MaterialNode* node, std::vector<const InteractionMechanism*>& others,
+ fullMatrix<double>& mat) const;
+ #endif //SWIG
+};
+*/
+// FLE
class CoefficientReduction
{
diff --git a/dG3D/src/computeWithMaterialLaw.cpp b/dG3D/src/computeWithMaterialLaw.cpp
index e5cb775c14d567fda9a4af1cdeea5f1c0f4eb48f..df5aeaa61872f21954664183b73e6816d3ce6a92 100644
--- a/dG3D/src/computeWithMaterialLaw.cpp
+++ b/dG3D/src/computeWithMaterialLaw.cpp
@@ -13,8 +13,8 @@
#include "nonLinearMechSolver.h"
#include "STensorOperations.h"
-computeMaterialLaw::computeMaterialLaw(dG3DMaterialLaw& mlaw): _materialLaw(&mlaw), _ips(NULL), _step(0)
-{
+computeMaterialLaw::computeMaterialLaw(dG3DMaterialLaw& mlaw, bool flag_isothermal, bool flag_microTempFixed): _materialLaw(&mlaw),
+_ips(NULL), _step(0), _flag_isothermal(flag_isothermal), _flag_microTempFixed(flag_microTempFixed){
static nonLinearMechSolver solver(1000);
_materialLaw->setMacroSolver(&solver);
static bool state = true;
@@ -78,6 +78,42 @@ void computeMaterialLaw::setDeformationGradient(int i, int j, double val)
F(i,j) = val;
}
+void computeMaterialLaw::setTemperatureGradient(int i, double val)
+{
+ if (!_flag_isothermal){
+ // store current for next step
+ ThermoMechanicsDG3DIPVariableBase* dgip = dynamic_cast<ThermoMechanicsDG3DIPVariableBase*>(_ips->getState(IPStateBase::current));
+ if (!dgip)
+ {
+ Msg::Error("ThermoMechanicsDG3DIPVariableBase must be created from material law computeMaterialLaw::setTemperatureGradient");
+ Msg::Exit(0);
+ }
+ SVector3& H = dgip->getRefToGradT();
+ H(i) = val;
+ }
+ else{
+ Msg::Error("computeMaterialLaw::setTemperatureGradient can be activated only if flag_isothermal is set to false in the constructor.");
+ }
+}
+
+void computeMaterialLaw::setTemperature(double val)
+{
+ if (!_flag_isothermal){
+ // store current for next step
+ ThermoMechanicsDG3DIPVariableBase* dgip = dynamic_cast<ThermoMechanicsDG3DIPVariableBase*>(_ips->getState(IPStateBase::current));
+ if (!dgip)
+ {
+ Msg::Error("ThermoMechanicsDG3DIPVariableBase must be created from material law computeMaterialLaw::setTemperature");
+ Msg::Exit(0);
+ }
+ double& T = dgip->getRefToTemperature();
+ T = val;
+ }
+ else{
+ Msg::Error("computeMaterialLaw::setTemperature can be activated only if flag_isothermal is set to false in the constructor.");
+ }
+}
+
void computeMaterialLaw::setTime(double t, double dt)
{
_materialLaw->setTime(t,dt);
@@ -108,7 +144,6 @@ void computeMaterialLaw::computeStressState_uniaxialStress(bool stiff)
}
};
-
void computeMaterialLaw::computeStressState_PlaneStrainUniaxialStress( char* NoZeroPlane, bool stiff )
{
IPVariable* ipvcur = _ips->getState(IPStateBase::current);
@@ -139,10 +174,6 @@ void computeMaterialLaw::computeStressState_PlaneStrainUniaxialStress( char* NoZ
};
-
-
-
-
double computeMaterialLaw::getTime() const
{
return _materialLaw->getTime();
@@ -397,6 +428,32 @@ void computeMaterialLaw::getStress(std::vector<double>& Pstress){
}
}
+void computeMaterialLaw::getHeatFlux(std::vector<double>& Qstress){
+ if(!_flag_isothermal){
+ IPVariable* ipvcur = _ips->getState(IPStateBase::current);
+ static SVector3 Q;
+ Q(0) = ipvcur->get(IPField::THERMALFLUX_X);
+ Q(1) = ipvcur->get(IPField::THERMALFLUX_Y);
+ Q(2) = ipvcur->get(IPField::THERMALFLUX_Z);
+ for(int i=0; i<3; i++){
+ Qstress.push_back(Q(i));
+ }
+ }
+ else{
+ Msg::Error("computeMaterialLaw::getHeatFlux can be used only if flag_isothermal is set to false in the constructor.");
+ }
+}
+
+void computeMaterialLaw::get3rdDependentVariable(double& Cstress){
+ if(!_flag_microTempFixed){
+ IPVariable* ipvcur = _ips->getState(IPStateBase::current);
+ Cstress = ipvcur->get(IPField::FIELD_CAPACITY); // specific heat
+ }
+ else{
+ Msg::Error("computeMaterialLaw::get3rdDependentVariable can be used only if _flag_microTempFixed is set to false in the constructor.");
+ }
+};
+
void computeMaterialLaw::getDeformationGradient(std::vector<double>& Fstrain){
IPVariable* ipvcur = _ips->getState(IPStateBase::current);
static STensor3 F;
@@ -416,6 +473,32 @@ void computeMaterialLaw::getDeformationGradient(std::vector<double>& Fstrain){
}
}
+void computeMaterialLaw::getTemperatureGradient(std::vector<double>& Hstrain){
+ if(!_flag_isothermal){
+ IPVariable* ipvcur = _ips->getState(IPStateBase::current);
+ static SVector3 H;
+ H(0) = ipvcur->get(IPField::TEMP_GRAD_X);
+ H(1) = ipvcur->get(IPField::TEMP_GRAD_Y);
+ H(2) = ipvcur->get(IPField::TEMP_GRAD_Z);
+ for(int i=0; i<3; i++){
+ Hstrain.push_back(H(i));
+ }
+ }
+ else{
+ Msg::Error("computeMaterialLaw::getTemperatureGradient can be used only if flag_isothermal is set to false in the constructor.");
+ }
+}
+
+void computeMaterialLaw::getTemperature(double& T){
+ if(!_flag_isothermal){
+ IPVariable* ipvcur = _ips->getState(IPStateBase::current);
+ T = ipvcur->get(IPField::TEMPERATURE);
+ }
+ else{
+ Msg::Error("computeMaterialLaw::getTemperatureGradient can be used only if flag_isothermal is set to false in the constructor.");
+ }
+}
+
void computeMaterialLaw::getParameterDerivative(std::vector< fullMatrix<double> >& dPdPn, std::vector< fullMatrix<double> >&dPdPv) const{
const IPVariable* ipvcur = _ips->getState(IPStateBase::current);
_materialLaw->getParameterDerivative(ipvcur, dPdPn, dPdPv);
diff --git a/dG3D/src/computeWithMaterialLaw.h b/dG3D/src/computeWithMaterialLaw.h
index 1d1ece95bd25d54a8311397cb50623cd5e9bcfa8..de66420e1893fefc1091da55184145a6031d2837 100644
--- a/dG3D/src/computeWithMaterialLaw.h
+++ b/dG3D/src/computeWithMaterialLaw.h
@@ -24,12 +24,16 @@ class computeMaterialLaw
dG3DMaterialLaw* _materialLaw;
IPStateBase* _ips;
int _step;
+ bool _flag_isothermal;
+ bool _flag_microTempFixed;
#endif //SWIG
public:
- computeMaterialLaw(dG3DMaterialLaw& mlaw);
+ computeMaterialLaw(dG3DMaterialLaw& mlaw, bool flag_isothermal = true, bool flag_microTempFixed = true);
virtual ~computeMaterialLaw();
void resetState(dG3DMaterialLaw& mlaw);
void setDeformationGradient(int i, int, double val);
+ void setTemperatureGradient(int i, double val);
+ void setTemperature(double val);
void setTime(double t, double dt);
void nextStep();
void computeStressState(bool stiff=false);
@@ -40,7 +44,11 @@ class computeMaterialLaw
double getTime() const;
double getTimeStep() const;
void getStress(std::vector<double>& Pstress);
+ void getHeatFlux(std::vector<double>& Qstress);
+ void get3rdDependentVariable(double& Cstress);
void getDeformationGradient(std::vector<double>& Fstrain);
+ void getTemperatureGradient(std::vector<double>& Hstrain);
+ void getTemperature(double& T);
void getParameterDerivative(std::vector< fullMatrix<double> >& dPdPn, std::vector< fullMatrix<double> >&dPdPv) const;
void reset_Parameter(dG3DMaterialLaw& mlaw, std::vector<std::vector<double>>& Para_Norm, std::vector<std::vector<double>>& Para_Wt, const double Vf=0.0);
void reset_Parameter(dG3DMaterialLaw& mlaw, const char* Para);
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index d29ab2917974a60884fe98ee77b294444df88f43..06d0c8844edbdbaff9e46650d729358e4b16a032 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -2165,134 +2165,6 @@ double ANNBasedDG3DIPVariable::get(const int comp) const
}
}
-
-StochDMNDG3DIPVariable::StochDMNDG3DIPVariable(const int NRoot,const int NInterface, const bool createBodyForceHO, const bool oninter):
- dG3DIPVariable(createBodyForceHO, oninter), _NRoot(NRoot), _NInterface(NInterface), _IPVector(NRoot,NULL){
- bool resizeFlag;
- int entryFP = 9*_NRoot;
- int entry_a = 3*_NInterface;
- resizeFlag = Fvec.resize(entryFP, true);
- resizeFlag = Pvec.resize(entryFP, true);
- resizeFlag = a.resize(entry_a, true);
- }
-
-StochDMNDG3DIPVariable::StochDMNDG3DIPVariable(const StochDMNDG3DIPVariable& src):
- dG3DIPVariable(src), _NRoot(src._NRoot), Fvec(src.Fvec), Pvec(src.Pvec), a(src.a),_IPVector(src._NRoot,NULL){
- for (int j=0; j< src._NRoot; j++){
- if (src._IPVector[j] != NULL){
- _IPVector[j] = src._IPVector[j]->clone();
- }
- }
- }
-
-
-StochDMNDG3DIPVariable& StochDMNDG3DIPVariable::operator =(const IPVariable& src)
-{
- dG3DIPVariable::operator=(src);
- const StochDMNDG3DIPVariable* psrc = dynamic_cast<const StochDMNDG3DIPVariable*>(&src);
- if (psrc!=NULL)
- {
- _NRoot = psrc->_NRoot;
- Fvec = psrc->Fvec;
- Pvec = psrc->Pvec;
- a = psrc->a;
- for (int j=0; j< _NRoot; j++)
- {
- if ((_IPVector[j]!=NULL) and (psrc->_IPVector[j] != NULL))
- {
- _IPVector[j]->operator =(*(psrc->_IPVector[j]));
- }
- else
- {
- Msg::Error("cannot perform equal operator in StochDMNDG3DIPVariable::operator =");
- }
- }
- }
- return *this;
-}
-
-StochDMNDG3DIPVariable::~StochDMNDG3DIPVariable()
-{
- for (int j=0; j< _NRoot; j++)
- {
- if (_IPVector[j] != NULL)
- {
- delete _IPVector[j];
- }
- }
-}
-
-void StochDMNDG3DIPVariable::addIPv(int loc, IPVariable* ipv)
-{
- _IPVector[loc] = ipv;
-}
-
-double StochDMNDG3DIPVariable::get(const int comp) const
-{
- if (comp == IPField::DEFO_ENERGY ||
- comp == IPField::PLASTIC_ENERGY ||
- comp == IPField::DAMAGE_ENERGY ||
- comp == IPField::DAMAGE )
- {
- double v = 0;
- for (int j=0; j< _NRoot; j++)
- {
- v += _IPVector[j]->get(comp);
- }
- return v;
- }
- else
- {
- return dG3DIPVariable::get(comp);
- }
-};
-
-
-void StochDMNDG3DIPVariable::restart()
-{
- dG3DIPVariable::restart();
- restartManager::restart(_NRoot);
- restartManager::restart(Fvec.getDataPtr(),9*_NRoot);
- restartManager::restart(Pvec.getDataPtr(),9*_NRoot);
- restartManager::restart(a.getDataPtr(),3*_NInterface);
- restartManager::restart(_IPVector);
- for (int j=0; j< _NRoot; j++){
- _IPVector[j]->restart();
- }
-};
-
-
-
-double StochDMNDG3DIPVariable::defoEnergy() const
-{
- double v = 0;
- for (int j=0; j< _NRoot; j++)
- {
- v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->defoEnergy();
- }
- return v;
-}
-double StochDMNDG3DIPVariable::plasticEnergy() const
-{
- double v = 0;
- for (int j=0; j< _NRoot; j++)
- {
- v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->plasticEnergy();
- }
- return v;
-}
-double StochDMNDG3DIPVariable::damageEnergy() const
-{
- double v = 0;
- for (int j=0; j< _NRoot; j++)
- {
- v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->damageEnergy();
- }
- return v;
-}
-
-
-
//torchANNBasedDG3DIPVariable::torchANNBasedDG3DIPVariable(const int n, const double initial_h, const bool createBodyForceHO, const bool oninter):
// dG3DIPVariable(createBodyForceHO, oninter), _n(n), _initial_h(initial_h), restart_internalVars(1,n), _kinematicVariables(1,6)
//{
@@ -4156,6 +4028,18 @@ double ThermoMechanicsDG3DIPVariableBase::get(const int i) const
else if (i == IPField::MECHANICAL_SOURCE){
return getConstRefToMechanicalSource()(0);
}
+ else if (i == IPField::TEMP_GRAD_X){
+ return getConstRefToGradField()[0](0);
+ }
+ else if (i == IPField::TEMP_GRAD_Y){
+ return getConstRefToGradField()[0](1);
+ }
+ else if (i == IPField::TEMP_GRAD_Z){
+ return getConstRefToGradField()[0](2);
+ }
+ else if (i == IPField::FIELD_CAPACITY){
+ return getConstitutiveExtraDofFieldCapacityPerUnitField(0);
+ }
else
{
return dG3DIPVariable::get(i);
@@ -4168,6 +4052,250 @@ void ThermoMechanicsDG3DIPVariableBase::restart()
return;
}
+StochDMNDG3DIPVariable::StochDMNDG3DIPVariable(const int NRoot,const int NInterface, const bool createBodyForceHO, const bool oninter):
+ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter), _NRoot(NRoot), _NInterface(NInterface), _IPVector(NRoot,NULL){
+ static STensor3 _I(1.); _R = _I;
+ bool resizeFlag;
+ int entryFP = 9*_NRoot;
+ int entry_a = 3*_NInterface;
+ resizeFlag = Fvec.resize(entryFP, true);
+ resizeFlag = Pvec.resize(entryFP, true);
+ resizeFlag = a.resize(entry_a, true);
+ }
+
+StochDMNDG3DIPVariable::StochDMNDG3DIPVariable(const StochDMNDG3DIPVariable& src):
+ThermoMechanicsDG3DIPVariableBase(src), _NRoot(src._NRoot), _NInterface(src._NInterface), Fvec(src.Fvec), Pvec(src.Pvec), a(src.a),_IPVector(src._NRoot,NULL), _R(src._R){
+ for (int j=0; j< src._NRoot; j++){
+ if (src._IPVector[j] != NULL){
+ _IPVector[j] = src._IPVector[j]->clone();
+ }
+ }
+ }
+
+
+StochDMNDG3DIPVariable& StochDMNDG3DIPVariable::operator =(const IPVariable& src)
+{
+ ThermoMechanicsDG3DIPVariableBase::operator=(src);
+ const StochDMNDG3DIPVariable* psrc = dynamic_cast<const StochDMNDG3DIPVariable*>(&src);
+ if (psrc!=NULL)
+ {
+ _NRoot = psrc->_NRoot;
+ _NInterface = psrc->_NInterface;
+ Fvec = psrc->Fvec;
+ Pvec = psrc->Pvec;
+ a = psrc->a;
+ _R = psrc->_R;
+ for (int j=0; j< _NRoot; j++)
+ {
+ if ((_IPVector[j]!=NULL) and (psrc->_IPVector[j] != NULL))
+ {
+ _IPVector[j]->operator =(*(psrc->_IPVector[j]));
+ }
+ else
+ {
+ Msg::Error("cannot perform equal operator in StochDMNDG3DIPVariable::operator =");
+ }
+ }
+ }
+ return *this;
+}
+
+StochDMNDG3DIPVariable::~StochDMNDG3DIPVariable()
+{
+ for (int j=0; j< _NRoot; j++)
+ {
+ if (_IPVector[j] != NULL)
+ {
+ delete _IPVector[j];
+ }
+ }
+}
+
+void StochDMNDG3DIPVariable::addIPv(int loc, IPVariable* ipv)
+{
+ _IPVector[loc] = ipv;
+}
+
+double StochDMNDG3DIPVariable::get(const int comp) const
+{
+ if (comp == IPField::DEFO_ENERGY ||
+ comp == IPField::PLASTIC_ENERGY ||
+ comp == IPField::DAMAGE_ENERGY ||
+ comp == IPField::DAMAGE )
+ {
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += _IPVector[j]->get(comp);
+ }
+ return v;
+ }
+ else
+ {
+ return ThermoMechanicsDG3DIPVariableBase::get(comp);
+ }
+};
+
+
+void StochDMNDG3DIPVariable::restart()
+{
+ ThermoMechanicsDG3DIPVariableBase::restart();
+ restartManager::restart(_NRoot);
+ restartManager::restart(_NInterface);
+ restartManager::restart(Fvec.getDataPtr(),9*_NRoot);
+ restartManager::restart(Pvec.getDataPtr(),9*_NRoot);
+ restartManager::restart(a.getDataPtr(),3*_NInterface);
+ restartManager::restart(_R);
+ restartManager::restart(_IPVector);
+ for (int j=0; j< _NRoot; j++){
+ _IPVector[j]->restart();
+ }
+};
+
+
+
+double StochDMNDG3DIPVariable::defoEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->defoEnergy();
+ }
+ return v;
+}
+double StochDMNDG3DIPVariable::plasticEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->plasticEnergy();
+ }
+ return v;
+}
+double StochDMNDG3DIPVariable::damageEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const dG3DIPVariableBase*>(_IPVector[j])->damageEnergy();
+ }
+ return v;
+}
+
+// FLE
+StochTMDMNDG3DIPVariable::StochTMDMNDG3DIPVariable(const int NRoot,const int NInterface, const int col_Di, const int row_Di,
+ const bool createBodyForceHO, const bool oninter):
+ StochDMNDG3DIPVariable(NRoot,NInterface,createBodyForceHO,oninter),_col_Di(col_Di), _row_Di(row_Di){
+ bool resizeFlag;
+ int entryFHPQ = _col_Di*_NRoot;
+ int entry_ab = _row_Di*_NInterface;
+ resizeFlag = FHvec.resize(entryFHPQ, true);
+ resizeFlag = PQvec.resize(entryFHPQ, true);
+ resizeFlag = ab.resize(entry_ab, true);
+ resizeFlag = Cpvec.resize(_NRoot, true);
+ resizeFlag = thermSrcvec.resize(_NRoot, true);
+ resizeFlag = mechSrcvec.resize(_NRoot, true);
+ }
+
+StochTMDMNDG3DIPVariable::StochTMDMNDG3DIPVariable(const StochTMDMNDG3DIPVariable& src): StochDMNDG3DIPVariable(src),
+ _col_Di(src._col_Di), _row_Di(src._row_Di),
+ FHvec(src.FHvec), PQvec(src.PQvec), ab(src.ab), Cpvec(src.Cpvec),
+ thermSrcvec(src.thermSrcvec), mechSrcvec(src.mechSrcvec){
+ }
+
+StochTMDMNDG3DIPVariable& StochTMDMNDG3DIPVariable::operator =(const IPVariable& src)
+{
+ StochDMNDG3DIPVariable::operator=(src);
+ const StochTMDMNDG3DIPVariable* psrc = dynamic_cast<const StochTMDMNDG3DIPVariable*>(&src);
+ if (psrc!=NULL)
+ {
+ _col_Di = psrc->_col_Di;
+ _row_Di = psrc->_row_Di;
+ FHvec = psrc->FHvec;
+ PQvec = psrc->PQvec;
+ ab = psrc->ab;
+ Cpvec = psrc->Cpvec;
+ thermSrcvec = psrc->thermSrcvec;
+ mechSrcvec = psrc->mechSrcvec;
+ }
+ return *this;
+}
+
+StochTMDMNDG3DIPVariable::~StochTMDMNDG3DIPVariable()
+{
+}
+
+void StochTMDMNDG3DIPVariable::addIPv(int loc, IPVariable* ipv)
+{
+ _IPVector[loc] = ipv;
+}
+
+double StochTMDMNDG3DIPVariable::get(const int comp) const
+{
+ if (comp == IPField::_thermalEnergy)
+ {
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += _IPVector[j]->get(comp);
+ }
+ return v;
+ }
+ else{
+ return StochDMNDG3DIPVariable::get(comp);
+ }
+};
+
+void StochTMDMNDG3DIPVariable::restart()
+{
+ StochDMNDG3DIPVariable::restart();
+ restartManager::restart(_col_Di);
+ restartManager::restart(_row_Di);
+ restartManager::restart(FHvec.getDataPtr(),_col_Di*_NRoot);
+ restartManager::restart(PQvec.getDataPtr(),_col_Di*_NRoot);
+ restartManager::restart(ab.getDataPtr(),_row_Di*_NInterface);
+};
+
+double StochTMDMNDG3DIPVariable::defoEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(_IPVector[j])->defoEnergy();
+ }
+ return v;
+}
+double StochTMDMNDG3DIPVariable::plasticEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(_IPVector[j])->plasticEnergy();
+ }
+ return v;
+}
+double StochTMDMNDG3DIPVariable::damageEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(_IPVector[j])->damageEnergy();
+ }
+ return v;
+}
+
+double StochTMDMNDG3DIPVariable::thermalEnergy() const
+{
+ double v = 0;
+ for (int j=0; j< _NRoot; j++)
+ {
+ v += dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(_IPVector[j])->getInternalEnergyExtraDofDiffusion();
+ }
+ return v;
+}
+
+// FLE
+
LinearThermoMechanicsDG3DIPVariable::LinearThermoMechanicsDG3DIPVariable(double tinitial, const bool createBodyForceHO, const bool oninter) :
ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter)
{
@@ -6094,6 +6222,7 @@ double NonLinearTVMDG3DIPVariable::get(const int i) const{
else if (i == IPField::PRESSURE){return _ipViscoElastic->_pressure;}
else if (i == IPField::viscousDissipatedEnergy){return _ipViscoElastic->_viscousDissipatedEnergy;}
else if (i == IPField::mullinsDissipatedEnergy){return _ipViscoElastic->_mullinsDissipatedEnergy;}
+ else if (i == IPField::_thermalEnergy){return _ipViscoElastic->_thermalEnergy;}
else if (i == IPField::Eve1_XX){
if (_ipViscoElastic->_A.size() > 0)
return _ipViscoElastic->_A[0](0,0) + _ipViscoElastic->_B[0]/3.;
@@ -6286,6 +6415,7 @@ double NonLinearTVPDG3DIPVariable::get(const int i) const{
else if (i == IPField::Dev_corKir_Inf_ZZ){return _ipViscoElastoPlastic->_devCorKirinf_TVE(2,2);}
else if (i == IPField::viscousDissipatedEnergy){return _ipViscoElastoPlastic->_viscousDissipatedEnergy;}
else if (i == IPField::mullinsDissipatedEnergy){return _ipViscoElastoPlastic->_mullinsDissipatedEnergy;}
+ else if (i == IPField::_thermalEnergy){return _ipViscoElastoPlastic->_thermalEnergy;}
else if (i == IPField::Eve1_XX){
if (_ipViscoElastoPlastic->_A.size() > 0)
return _ipViscoElastoPlastic->_A[0](0,0) + _ipViscoElastoPlastic->_B[0]/3.;
diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h
index 492538adbaccafc523676ca0d47f9612fb1115ba..6bfee1b0890d80b3f275ae3dfba8a543631724a7 100644
--- a/dG3D/src/dG3DIPVariable.h
+++ b/dG3D/src/dG3DIPVariable.h
@@ -2395,48 +2395,6 @@ class ANNBasedDG3DIPVariable :public dG3DIPVariable
virtual void restart();
};
-
-class StochDMNDG3DIPVariable : public dG3DIPVariable
-{
- protected:
- int _NRoot, _NInterface;
- std::vector<IPVariable*> _IPVector;
- fullVector<double> Fvec;
- fullVector<double> Pvec;
- fullVector<double> a;
-
- public:
- StochDMNDG3DIPVariable(const int NRoot,const int NInterface, const bool createBodyForceHO, const bool oninter);
- StochDMNDG3DIPVariable(const StochDMNDG3DIPVariable& src);
- virtual StochDMNDG3DIPVariable& operator =(const IPVariable& src);
- virtual ~StochDMNDG3DIPVariable();
-
- void InitializeFluctuationVector() {a.setAll(0.0);};
- virtual double defoEnergy() const;
- virtual double plasticEnergy() const;
- virtual double damageEnergy() const;
-
- void addIPv(int i, IPVariable* ipv);
- IPVariable* getIPv(int i) {return _IPVector[i];};
- const IPVariable* getIPv(int i) const {return _IPVector[i];};
-
- virtual double get(const int comp) const;
- virtual fullVector<double>& getRefToDeformationGradientVect() {return Fvec;};
- virtual const fullVector<double>& getConstRefToDeformationGradientVect() const {return Fvec;};
- virtual fullVector<double>& getRefToFirstPiolaKirchhoffStressVect() {return Pvec;};
- virtual const fullVector<double>& getConstRefToFirstPiolaKirchhoffStressVect() const {return Pvec;};
- virtual fullVector<double>& getRefToInterfaceFluctuationVect() {return a;};
- virtual const fullVector<double>& getConstRefToInterfaceFluctuationVect() const {return a;};
-
-
-
- virtual IPVariable* getInternalData() {return NULL;};
- virtual const IPVariable* getInternalData() const {return NULL;};
- virtual IPVariable* clone() const {return new StochDMNDG3DIPVariable(*this);};
- virtual void restart();
-};
-
-
class torchANNBasedDG3DIPVariable :public dG3DIPVariable
{
protected:
@@ -3476,6 +3434,99 @@ class ThermoMechanicsDG3DIPVariableBase : public dG3DIPVariable, public extraDof
virtual void restart();
};
+class StochDMNDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase
+{
+ protected:
+ int _NRoot, _NInterface;
+ std::vector<IPVariable*> _IPVector;
+ fullVector<double> Fvec;
+ fullVector<double> Pvec;
+ fullVector<double> a;
+ STensor3 _R;
+
+ public:
+ StochDMNDG3DIPVariable(const int NRoot,const int NInterface, const bool createBodyForceHO, const bool oninter);
+ StochDMNDG3DIPVariable(const StochDMNDG3DIPVariable& src);
+ virtual StochDMNDG3DIPVariable& operator =(const IPVariable& src);
+ virtual ~StochDMNDG3DIPVariable();
+
+ void InitializeFluctuationVector() {a.setAll(0.0);};
+ virtual double defoEnergy() const;
+ virtual double plasticEnergy() const;
+ virtual double damageEnergy() const;
+ virtual double getInternalEnergyExtraDofDiffusion() const {return 0.;};
+
+ void addIPv(int i, IPVariable* ipv);
+ IPVariable* getIPv(int i) {return _IPVector[i];};
+ const IPVariable* getIPv(int i) const {return _IPVector[i];};
+
+ virtual double get(const int comp) const;
+ virtual fullVector<double>& getRefToDeformationGradientVect() {return Fvec;};
+ virtual const fullVector<double>& getConstRefToDeformationGradientVect() const {return Fvec;};
+ virtual fullVector<double>& getRefToFirstPiolaKirchhoffStressVect() {return Pvec;};
+ virtual const fullVector<double>& getConstRefToFirstPiolaKirchhoffStressVect() const {return Pvec;};
+ virtual fullVector<double>& getRefToInterfaceFluctuationVect() {return a;};
+ virtual const fullVector<double>& getConstRefToInterfaceFluctuationVect() const {return a;};
+
+ virtual STensor3& getReftoRotationMatrix() {return _R;};
+ virtual const STensor3& getConstReftoRotationMatrix() const {return _R;};
+
+ virtual IPVariable* getInternalData() {return NULL;};
+ virtual const IPVariable* getInternalData() const {return NULL;};
+ virtual IPVariable* clone() const {return new StochDMNDG3DIPVariable(*this);};
+ virtual void restart();
+};
+
+// FLE
+class StochTMDMNDG3DIPVariable : public StochDMNDG3DIPVariable
+{
+ protected:
+ int _col_Di, _row_Di;
+
+ // combined vectors
+ fullVector<double> FHvec; // F- defo, H -tempGrad
+ fullVector<double> PQvec; // P- 1st PK, Q-heatFlux
+ fullVector<double> ab; // a- disp, b-temp
+ fullVector<double> Cpvec, thermSrcvec, mechSrcvec;
+
+ public:
+ StochTMDMNDG3DIPVariable(const int NRoot, const int NInterface, const int col_Di, const int row_Di, const bool createBodyForceHO, const bool oninter);
+ StochTMDMNDG3DIPVariable(const StochTMDMNDG3DIPVariable& src);
+ virtual StochTMDMNDG3DIPVariable& operator =(const IPVariable& src);
+ virtual ~StochTMDMNDG3DIPVariable();
+
+ void InitializeFluctuationVector() {ab.setAll(0.0);};
+ virtual double defoEnergy() const;
+ virtual double plasticEnergy() const;
+ virtual double damageEnergy() const;
+ virtual double thermalEnergy() const;
+ virtual double getInternalEnergyExtraDofDiffusion() const {return this->thermalEnergy();};
+
+ void addIPv(int i, IPVariable* ipv);
+ IPVariable* getIPv(int i) {return _IPVector[i];};
+ const IPVariable* getIPv(int i) const {return _IPVector[i];};
+
+ virtual double get(const int comp) const;
+ virtual fullVector<double>& getRefToCombinedGradientVect() {return FHvec;};
+ virtual const fullVector<double>& getConstRefToCombinedGradientVect() const {return FHvec;};
+ virtual fullVector<double>& getRefToCombinedStressVect() {return PQvec;};
+ virtual const fullVector<double>& getConstRefToCombinedStressVect() const {return PQvec;};
+ virtual fullVector<double>& getRefToCombinedInterfaceFluctuationVect() {return ab;};
+ virtual const fullVector<double>& getConstRefToCombinedInterfaceFluctuationVect() const {return ab;};
+ virtual fullVector<double>& getRefToCpVect() {return Cpvec;};
+ virtual const fullVector<double>& getConstRefToCpVect() const {return Cpvec;};
+ virtual fullVector<double>& getRefToThermalSourceVect() {return thermSrcvec;};
+ virtual const fullVector<double>& getConstRefToThermalSourceVect() const {return thermSrcvec;};
+ virtual fullVector<double>& getRefToMechSourceVect() {return mechSrcvec;};
+ virtual const fullVector<double>& getConstRefToMechSourceVect() const {return mechSrcvec;};
+
+ virtual IPVariable* getInternalData() {return NULL;};
+ virtual const IPVariable* getInternalData() const {return NULL;};
+ virtual IPVariable* clone() const {return new StochTMDMNDG3DIPVariable(*this);};
+ virtual void restart();
+};
+// FLE
+
class LinearThermoMechanicsDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase
{
@@ -5688,7 +5739,6 @@ public:
// Take Care here, nonlocal variables cannot be directly defined with the ThermoMechanicsDG3DIPVariable Base Class (This is a pure virtual base class).
-// class NonLinearTVMDG3DIPVariable : public dG3DIPVariable{
class NonLinearTVMDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase{ // Changed
protected:
@@ -5696,8 +5746,6 @@ class NonLinearTVMDG3DIPVariable : public ThermoMechanicsDG3DIPVariableBase{ /
NonLinearTVMDG3DIPVariable(const bool createBodyForceHO, const bool oninter):ThermoMechanicsDG3DIPVariableBase(createBodyForceHO, oninter) {Msg::Error("NonLinearTVMDG3DIPVariable: Cannot use the default constructor of ThermoMechanicsDG3DIPVariableBase");}
public:
-// Constructor Changed
- // NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, const bool oninter=false, int nbExtraDof=0); // ?? - WHAT to add/change?
NonLinearTVMDG3DIPVariable(const mlawNonLinearTVM& viscoLaw, double Tinitial, const bool createBodyForceHO, const bool oninter); // Changed
NonLinearTVMDG3DIPVariable(const NonLinearTVMDG3DIPVariable& src);
virtual NonLinearTVMDG3DIPVariable& operator = (const IPVariable& src);
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 61e231436d125eab80eb4a5fe0a318438e75fd54..a8a50547267dee52bfbebec6e067c7a75a5ceb6e 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -3633,21 +3633,27 @@ double torchANNBasedDG3DMaterialLaw::soundSpeed() const
-StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree, const bool porous,const double tol):
- dG3DMaterialLaw(num,rho,false), _ReadTree(ReadTree), _porous(porous),_tol(tol){
- if(_ReadTree){
- read_TreeData(ParaFile);
- }
- else{
- read_parameter(ParaFile);
- }
- fill_Matrices();
+StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree, const bool porous,const double tol, const bool flag_isothermal):
+ dG3DMaterialLaw(num,rho,false), _ReadTree(ReadTree), _porous(porous),_tol(tol), _flag_isothermal(flag_isothermal){
+
+ _alpha = 0.; _beta = 0.; _gamma = 0.;
+
+ if(_flag_isothermal){
+ if(_ReadTree){
+ read_TreeData(ParaFile);
+ }
+ else{
+ read_parameter(ParaFile);
+ }
+ fill_Matrices();
+ }
}
StochDMNDG3DMaterialLaw::StochDMNDG3DMaterialLaw(const StochDMNDG3DMaterialLaw &src): dG3DMaterialLaw(src), _ReadTree(src._ReadTree), _porous(src._porous),
_LevelNode(src._LevelNode), _NodeLevel(src._NodeLevel), _Mat_Id(src._Mat_Id), _Parent(src._Parent), _Child(src._Child), _TotalLevel(src._TotalLevel), _NLeaf(src._NLeaf),
_NInterface(src._NInterface), _N_Node(src._N_Node), _NPara_Wt(src._NPara_Wt), _NPara_Norm(src._NPara_Norm), _Dim(src._Dim), _Vf(src._Vf), _tol(src._tol),
- _Nv(src._Nv), _NTV(src._NTV), _V(src._V), _Para_Wt(src._Para_Wt), _Para_Norm(src._Para_Norm), _VA(src._VA), _VI(src._VI), _mapLaw(src._mapLaw){ }
+ _Nv(src._Nv), _NTV(src._NTV), _V(src._V), _Para_Wt(src._Para_Wt), _Para_Norm(src._Para_Norm), _VA(src._VA), _VI(src._VI), _mapLaw(src._mapLaw),
+ _flag_isothermal(src._flag_isothermal), _alpha(src._alpha), _beta(src._beta), _gamma(src._gamma){ }
StochDMNDG3DMaterialLaw::~StochDMNDG3DMaterialLaw(){};
@@ -3872,6 +3878,9 @@ void StochDMNDG3DMaterialLaw::getLeafOfInterface(const int pos, int &node_start,
}
void StochDMNDG3DMaterialLaw::fill_NLocal(const int pos, double N[][3])const{
+
+ // code to fill a matrix related to local normals
+
double n0 = _Para_Norm[pos][0];
double n1 = _Para_Norm[pos][1];
double n2 = _Para_Norm[pos][2];
@@ -3897,7 +3906,7 @@ void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos, double W[][2]){
double VI = _VI[pos];
double Wtmp0 = _Para_Wt[pos][0];
double Wtmp1 = _Para_Wt[pos][1];
- double WA = Wtmp0*(1.0-VI)+ Wtmp1*VI;
+ double WA = Wtmp0*(1.0-VI)+ Wtmp1*VI; // Eq. 39 in Ling's draft
for(int i=0; i<2; i++){
for(int j=0; j<2; j++){
@@ -3906,13 +3915,14 @@ void StochDMNDG3DMaterialLaw::fill_W_WI(const int pos, double W[][2]){
}
W[0][0]= WA;
W[0][1]= 1.0-WA;
- W[1][0]= (Wtmp1*VI)/WA;
- W[1][1]= ((1.0-Wtmp1)*VI)/(1.0-WA);
+ W[1][0]= (Wtmp1*VI)/WA; // Eq. 40 in Ling's draft, v_1A
+ W[1][1]= ((1.0-Wtmp1)*VI)/(1.0-WA); // Eq. 40 in Ling's draft, v_1B
}
// fill the matrices which keep the topological information
-void StochDMNDG3DMaterialLaw::fill_Matrices(){
+void StochDMNDG3DMaterialLaw::fill_Matrices(){
+
double Norm[9][3];
double W_WI[2][2];
int pos0, pos1;
@@ -4116,7 +4126,7 @@ void StochDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
}
}
- Msg::Info("Finish reading parameter filefor Complete Binary Tree, level = %d, Dim = %d, NLeaf = %d", _TotalLevel, _Dim, _NLeaf);
+ Msg::Info("Finish reading parameter file for Complete Binary Tree, level = %d, Dim = %d, NLeaf = %d", _TotalLevel, _Dim, _NLeaf);
}
// end of read_parameter
@@ -4385,7 +4395,7 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
}
}
- _NTV.mult(Pvec, Res_node);
+ _NTV.mult(Pvec, Res_node); // Eq 60 in Ling's draft, Res_node = residual at the node
r = Res_node.norm()/_NInterface;
_NTV.mult(_C, NTVC);
@@ -4418,7 +4428,7 @@ void StochDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, c
if(stiff){
NTVC.mult(_I, dRdF);
- Jacobian.mult(dRdF, dadF);
+ Jacobian.mult(dRdF, dadF); // Eq. 68 in Ling's draft
_Nv.mult(dadF, tmp);
tmp.add(_I);
@@ -4994,7 +5004,956 @@ void StochDMNDG3DMaterialLaw::setTime(const double t,const double dtime){
}
}
-//
+void StochDMNDG3DMaterialLaw::eulerZXZToRotationMatrix(const double phi, const double theta, const double psi, STensor3& rotationMatrix) {
+ double cosPhi = cos(phi);
+ double sinPhi = sin(phi);
+ double cosTheta = cos(theta);
+ double sinTheta = sin(theta);
+ double cosPsi = cos(psi);
+ double sinPsi = sin(psi);
+
+ rotationMatrix(0,0) = cosPsi * cosPhi - cosTheta * sinPhi * sinPsi;
+ rotationMatrix(0,1) = cosPsi * sinPhi + cosTheta * cosPhi * sinPsi;
+ rotationMatrix(0,2) = sinPsi * sinTheta;
+
+ rotationMatrix(1,0) = -sinPsi * cosPhi - cosTheta * sinPhi * cosPsi;
+ rotationMatrix(1,1) = -sinPsi * sinPhi + cosTheta * cosPhi * cosPsi;
+ rotationMatrix(1,2) = cosPsi * sinTheta;
+
+ rotationMatrix(2,0) = sinTheta * sinPhi;
+ rotationMatrix(2,1) = -sinTheta * cosPhi;
+ rotationMatrix(2,2) = cosTheta;
+}
+
+void StochDMNDG3DMaterialLaw::setZXZRotationMatrix_in_IP(StochDMNDG3DIPVariable* ipv){
+ STensor3& R = ipv->getReftoRotationMatrix();
+ eulerZXZToRotationMatrix(_alpha,_beta,_gamma,R);
+};
+
+// FLE
+StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree,
+ const bool porous, const bool flag_isothermal, const bool flag_microTempFixed, const double tol):
+ StochDMNDG3DMaterialLaw(num,rho,E,nu,ParaFile,ReadTree,porous,tol,flag_isothermal),_flag_microTempFixed(flag_microTempFixed){
+
+ // Meanings of _col_Di and _row_Di are interchanged in this class.
+ if (_flag_isothermal && _flag_microTempFixed){
+ _col_Di = 9; _row_Di = 3;
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ _col_Di = 12; _row_Di = 4;
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){ // Deprecated option - Impossible to implement without length scale.
+ _col_Di = 13; _row_Di = 5;
+ }
+ else{
+ Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
+ }
+ if(_ReadTree){
+ StochTMDMNDG3DMaterialLaw::read_TreeData(ParaFile);
+ }
+ else{
+ StochTMDMNDG3DMaterialLaw::read_parameter(ParaFile);
+ }
+
+ StochTMDMNDG3DMaterialLaw::fill_Matrices();
+}
+
+StochTMDMNDG3DMaterialLaw::StochTMDMNDG3DMaterialLaw(const StochTMDMNDG3DMaterialLaw &source): StochDMNDG3DMaterialLaw(source),
+ _col_Di(source._col_Di), _row_Di(source._row_Di), _flag_microTempFixed(source._flag_microTempFixed), _Para_Alpha(source._Para_Alpha){};
+
+StochTMDMNDG3DMaterialLaw::~StochTMDMNDG3DMaterialLaw(){};
+
+void StochTMDMNDG3DMaterialLaw::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;
+ IPVariable* ipvi = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
+ IPVariable* ipv1 = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
+ IPVariable* ipv2 = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
+
+ StochTMDMNDG3DIPVariable* ipvi_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipvi);
+ StochTMDMNDG3DIPVariable* ipv1_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipv1);
+ StochTMDMNDG3DIPVariable* ipv2_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipv2);
+ int mat_IP;
+ for (int i=0; i< _NLeaf; i++){
+ mat_IP = _Mat_Id[_NInterface + i];
+ IPStateBase* ips_i = NULL;
+ _mapLaw[mat_IP]->createIPState(ips_i, hasBodyForce, state_, ele, nbFF_, GP, gpt);
+ std::vector<IPVariable*> ip_all;
+ ips_i->getAllIPVariable(ip_all);
+ ipvi_Leaf->addIPv(i, ip_all[0]);
+ ipv1_Leaf->addIPv(i, ip_all[1]);
+ ipv2_Leaf->addIPv(i, ip_all[2]);
+ }
+ if(ips != NULL) delete ips;
+ ips = new IP3State(state_,ipvi,ipv1,ipv2);
+}
+
+void StochTMDMNDG3DMaterialLaw::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;
+ }
+ ipv = new StochTMDMNDG3DIPVariable(_NLeaf, _NInterface, _col_Di, _row_Di, hasBodyForce,inter);
+ StochTMDMNDG3DIPVariable* ipv_Leaf = static_cast<StochTMDMNDG3DIPVariable*>(ipv);
+ int mat_IP;
+ for (int i=0; i< _NLeaf; i++){
+ IPVariable* ipv_i = NULL;
+ mat_IP = _Mat_Id[_NInterface + i];
+ _mapLaw[mat_IP]->createIPVariable(ipv_i, hasBodyForce, ele, nbFF_, GP, gpt);
+ ipv_Leaf->addIPv(i, ipv_i);
+ }
+};
+
+void StochTMDMNDG3DMaterialLaw::initialIPVariable(IPVariable* ipv, bool stiff){
+ StochTMDMNDG3DIPVariable* ipv_all = dynamic_cast<StochTMDMNDG3DIPVariable*>(ipv);
+ ipv_all->InitializeFluctuationVector();
+ if (ipv_all == NULL){
+ Msg::Error("StochTMDMNDG3DIPVariable must be used in StochTMDMNDG3DMaterialLaw::initialIPVariable");
+ }
+ int mat_IP;
+ for (int i=0; i< _NLeaf; i++){
+ IPVariable* ipv_i = ipv_all->getIPv(i);
+ mat_IP = _Mat_Id[_NInterface + i];
+ _mapLaw[mat_IP]->initialIPVariable(ipv_i, stiff);
+ }
+}
+
+void StochTMDMNDG3DMaterialLaw::checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{
+ StochTMDMNDG3DIPVariable* ipv_all = dynamic_cast<StochTMDMNDG3DIPVariable*>(ipv);
+ const StochTMDMNDG3DIPVariable* ipvprev_all = dynamic_cast<const StochTMDMNDG3DIPVariable*>(ipvprev);
+ if (ipv_all == NULL){
+ Msg::Error("StochTMDMNDG3DIPVariable must be used in StochTMDMNDG3DMaterialLaw::checkInternalState");
+ }
+ int mat_IP;
+ for (int i=0; i< _NLeaf; i++){
+ IPVariable* ipv_i = ipv_all->getIPv(i);
+ const IPVariable* ipvprev_i = ipvprev_all->getIPv(i);
+ mat_IP = _Mat_Id[_NInterface + i];
+ _mapLaw[mat_IP]->checkInternalState(ipv_i, ipvprev_i);
+ }
+}
+
+void StochTMDMNDG3DMaterialLaw::setMacroSolver(const nonLinearMechSolver* sv){
+ dG3DMaterialLaw::setMacroSolver(sv);
+ for (std::vector<dG3DMaterialLaw*>::iterator ite = _mapLaw.begin(); ite!= _mapLaw.end(); ite++){
+ (*ite)->setMacroSolver(sv);
+ }
+}
+
+void StochTMDMNDG3DMaterialLaw::fill_NLocal(const int pos, fullMatrix<double> N)const{
+
+ // code to fill a matrix related to local normals
+
+ double n0 = _Para_Norm[pos][0];
+ double n1 = _Para_Norm[pos][1];
+ double n2 = _Para_Norm[pos][2];
+
+ if (_flag_isothermal && _flag_microTempFixed){
+ N.set(0,0,n0); N.set(1,1,n0); N.set(2,2,n0);
+ N.set(3,0,n1); N.set(4,1,n1); N.set(5,2,n1);
+ N.set(6,0,n2); N.set(7,1,n2); N.set(8,2,n2);
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ N.set(0,0,n0); N.set(1,1,n0); N.set(2,2,n0);
+ N.set(3,0,n1); N.set(4,1,n1); N.set(5,2,n1);
+ N.set(6,0,n2); N.set(7,1,n2); N.set(8,2,n2);
+ N.set(9,3,n0); N.set(10,3,n1); N.set(11,3,n2);
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){
+ N.set(0,0,n0); N.set(1,1,n0); N.set(2,2,n0);
+ N.set(3,0,n1); N.set(4,1,n1); N.set(5,2,n1);
+ N.set(6,0,n2); N.set(7,1,n2); N.set(8,2,n2);
+ N.set(9,3,n0); N.set(10,3,n1); N.set(11,3,n2);
+ N.set(12,4,1.);
+ }
+ else{
+ Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
+ }
+}
+
+void StochTMDMNDG3DMaterialLaw::initLaws(const std::map<int,materialLaw*> &maplaw){
+ // To define this function, StochDMNDG3DMaterialLaw::TwoPhasesInteract has to be redefined in this class.
+}
+
+// fill the matrices which keep the topological information
+void StochTMDMNDG3DMaterialLaw::fill_Matrices(){
+
+ fullMatrix<double> Norm; bool resizeFlag = Norm.resize(_col_Di, _row_Di, true);
+ double W_WI[2][2];
+ int pos0, pos1;
+ int L_ch0, L_ch1;
+ int node_start0, node_end0, node_start1, node_end1;
+ int Col_start, Col_start0, Col_start1, Row_start, Row_start0, Row_start1;
+ double tmp;
+ int col = _col_Di; // = 9, 12; // set the size of rows and cols in Di matrix
+
+ _VA.resize(_N_Node,0.0); // Vector for position of nodes for Stress
+ _VA[0] = 1.0;
+ _VI.resize(_NInterface, 0.0); // fill matrix for equilibrium of interface
+ _VI[0] = _Vf;
+ int l;
+
+ for(int pos=0; pos <_NInterface; pos++){
+ l = _NodeLevel[pos];
+ StochTMDMNDG3DMaterialLaw::fill_NLocal(pos,Norm); // size of Norm depends on the flags set inside this function
+ pos0 = _Child[pos][0];
+ pos1 = _Child[pos][1];
+ L_ch0 = _NodeLevel[pos0];
+ L_ch1 = _NodeLevel[pos1];
+ Row_start = pos*_row_Di; // 3, 4;
+ Col_start0 = (pos0-1)*col;
+ Col_start1 = (pos1-1)*col;
+
+ double S_ba = 1.0;
+ if(l ==_TotalLevel-1 and _porous) S_ba = _Para_Wt[pos][1];
+
+ // information of interface norm
+ for(int nc=0; nc<col;nc++){
+ for(int nr=0; nr<_row_Di; nr++){
+ _NT.set(Row_start+nr, Col_start0+nc, Norm(nc,nr));
+ _NT.set(Row_start+nr, Col_start1+nc,-S_ba*Norm(nc,nr));
+ }
+ }
+
+ if(l ==_TotalLevel-1){ // Last level
+ if(_porous){
+ _VA[pos0] = (1.0-_VI[pos])*_Para_Wt[pos][0];
+ _VA[pos1] = (1.0-_VI[pos])*(1.0-_Para_Wt[pos][0]);
+ }
+ else{
+ _VA[pos0] = 1.0-_VI[pos];
+ _VA[pos1] = _VI[pos];
+ }
+ }
+ else{ // Composite nodes - all other levels
+ StochDMNDG3DMaterialLaw::fill_W_WI(pos,W_WI);
+ _VA[pos0] = W_WI[0][0];
+ _VA[pos1] = W_WI[0][1];
+ if(L_ch0 < _TotalLevel) _VI[pos0]=W_WI[1][0];
+ if(L_ch1 < _TotalLevel) _VI[pos1]=W_WI[1][1];
+ // Msg::Info("Node= %d,%d,%d, VI = %f,%f,%f",pos,pos0,pos1,_VI[pos],W_WI[1][0],W_WI[1][1]);
+ }
+ //information of homogenous strain to local strain
+ Col_start = pos*_row_Di; // *3, 4 // Note that the meanings of _row_Di and _col_Di have been swapped
+ StochDMNDG3DMaterialLaw::getLeafOfInterface(pos0, node_start0, node_end0);
+ StochDMNDG3DMaterialLaw::getLeafOfInterface(pos1, node_start1, node_end1);
+ for(int r=node_start0; r<=node_end0; r++){
+ Row_start0 = (r-_NInterface)*col;
+ for(int nr=0; nr<col;nr++){
+ for(int nc=0; nc<_row_Di;nc++){
+ _Nv.set(Row_start0+nr, Col_start+nc, Norm(nr,nc)/_VA[pos0]);
+ }
+ }
+ }
+ for(int r=node_start1; r<=node_end1; r++){
+ Row_start1 = (r-_NInterface)*col;
+ for(int nr=0; nr<col;nr++){
+ for(int nc=0; nc<_row_Di;nc++){
+ _Nv.set(Row_start1+nr, Col_start+nc,-S_ba*Norm(nr,nc)/_VA[pos1]);
+ }
+ }
+ }
+ }
+
+ // information of volume fraction to stress on node
+ Row_start = col*(_N_Node-1-_NLeaf);
+ for(int nc=0; nc<col*_NLeaf; nc++){
+ _V.set(Row_start+nc,nc,1.0);
+ }
+
+ for(int pos = _NInterface-1; pos > 0; pos--){
+ pos0 = _Child[pos][0];
+ pos1 = _Child[pos][1];
+ Row_start = (pos-1)*col;
+ Row_start0 = (pos0-1)*col;
+ Row_start1 = (pos1-1)*col;
+ for(int nr=0; nr<col;nr++){
+ for(int nc=0; nc<col*_NLeaf ; nc++){
+ tmp = _VA[pos0]*_V(Row_start0+nr,nc) + _VA[pos1]*_V(Row_start1+nr,nc);
+ _V.set(Row_start+nr,nc, tmp);
+ }
+ }
+ }
+ _NT.mult(_V, _NTV);
+
+ for(int i=0; i<col*_NLeaf;i++){
+ _I(i, i%col) = 1.0;
+ }
+}
+ // end fill_matrices
+
+// read parameter from input file
+void StochTMDMNDG3DMaterialLaw::read_parameter(const char *ParaFile){
+
+ double sin0, cos0, sin1, cos1;
+
+ FILE *Para = fopen(ParaFile, "r");
+ if ( Para != NULL ){
+ int okf = fscanf(Para, "%d %d\n", &_Dim, &_TotalLevel);
+
+ bool resizeFlag;
+ _NLeaf = int(pow(2,_TotalLevel));
+ _NInterface = _NLeaf-1;
+ _N_Node = _NLeaf + _NInterface;
+
+ int row = _col_Di*_NLeaf;
+ int col = _row_Di*_NInterface;
+ resizeFlag = _Nv.resize(row, col, true);
+
+ resizeFlag = _I.resize(row, _col_Di, true);
+
+ row = _col_Di*(_N_Node-1);
+ col = _col_Di*_NLeaf;
+ resizeFlag = _V.resize(row, col, true);
+
+ row = _row_Di*_NInterface;
+ resizeFlag = _NTV.resize(row, col, true);
+
+ resizeFlag = _NT.resize(row, _col_Di*(_N_Node-1), true);
+
+ okf = fscanf(Para, "%lf\n", &_Vf);
+
+ _NPara_Norm = _NInterface;
+
+ resizeFlag = _C.resize(col,col,true);
+ resizeFlag = Jacobian.resize(row,row,true);
+
+ resizeFlag = _ParaAngle.resize(_NInterface,2,true);
+
+ double ParaN[2]; // W_alpha and w_beta -> Eq. 37 in Ling's draft
+ double Pi(3.14159265359);
+ SPoint3 Para_Norm;
+ SPoint2 Para_Wt;
+
+ for(int i=0;i<_NInterface;i++){
+ if(_Dim==2){
+ okf = fscanf(Para, "%lf\n", &ParaN[0]);
+ ParaN[1] = 0.0;
+ }
+ else if(_Dim==3){
+ okf = fscanf(Para, "%lf %lf\n", &ParaN[0], &ParaN[1]);
+ }
+ sin0 = sin(2.0*Pi*ParaN[0]);
+ cos0 = cos(2.0*Pi*ParaN[0]);
+ sin1 = sin(Pi*ParaN[1]);
+ cos1 = cos(Pi*ParaN[1]);
+ Para_Norm[0] = cos1*sin0;
+ Para_Norm[1] = cos1*cos0;
+ Para_Norm[2] = sin1;
+ _Para_Norm.push_back(Para_Norm);
+ _ParaAngle(i,0) = ParaN[0];
+ _ParaAngle(i,1) = ParaN[1];
+ }
+
+ if(_porous){
+ _NPara_Wt = _NInterface;
+ }
+ else{
+ _NPara_Wt = int(pow(2,(_TotalLevel-1)))-1;
+ }
+ for(int i=0;i<_NPara_Wt;i++){
+ okf = fscanf(Para, "%lf %lf\n", &Para_Wt[0], &Para_Wt[1]);
+ _Para_Wt.push_back(Para_Wt);
+ }
+ }
+ // creat data as for tree data
+ int n = 0;
+ int m = -10;
+ for(int l=0; l<=_TotalLevel; l++){
+ std::vector<int> ln;
+ _LevelNode.push_back(ln);
+ while(n < pow(2,l+1)-1){
+ _LevelNode[l].push_back(n);
+ _NodeLevel.push_back(l);
+ if(n==0){
+ _Parent.push_back(0);
+ }
+ else{
+ _Parent.push_back(int((n-1)/2));
+ }
+ if(l==_TotalLevel){
+ _Mat_Id.push_back((n+1)%2);
+ }
+ else{
+ std::vector<int> ch;
+ ch.push_back(2*n+1);
+ ch.push_back(2*n+2);
+ _Child.push_back(ch);
+ _Mat_Id.push_back(m);
+ }
+ n +=1;
+ }
+ }
+
+ Msg::Info("Finished reading parameter file for Complete Binary Tree, level = %d, Dim = %d, NLeaf = %d", _TotalLevel, _Dim, _NLeaf);
+}
+ // end of read_parameter
+
+void StochTMDMNDG3DMaterialLaw::read_TreeData(const char *TreeFile){
+ std::ifstream file(TreeFile);
+ if(!file.is_open()) {
+ Msg::Error( "Error: Could not open the tree file." );
+ return ;
+ }
+ int level_value, node_value, parent_value;
+ double sin0, cos0, sin1,cos1;
+ double ParaN[2];
+ int pn = 0;
+ double Pi(3.14159265359);
+ SPoint3 Para_Norm;
+ SPoint2 Para_Wt;
+ bool resizeFlag;
+
+ std::string line;
+ int DataId = 0;
+ while(std::getline(file, line)) {
+ if(line.empty()) continue;
+ std::istringstream stream(line);
+ std::string value;
+ std::vector<std::string> row;
+ while (std::getline(stream, value, ',')) {
+ row.push_back(value);
+ }
+ if (std::isalpha(row[0][0])){
+ DataId +=1;
+ if(DataId == 4){
+ _NLeaf = static_cast<int>(_LevelNode[_TotalLevel].size());
+ _N_Node = static_cast<int>(_Parent.size());
+ _NInterface = _N_Node - _NLeaf;
+ resizeFlag = _ParaAngle.resize(_NInterface,2,true);
+ }
+ }
+ else{
+ if(DataId == 1){ // read the tree structure------------------
+ _TotalLevel = std::stoi(row[0]);
+ for(int i=0; i<=_TotalLevel; i++){
+ std::vector<int> LN;
+ _LevelNode.push_back(LN);
+ }
+ }
+ else if(DataId == 2){
+ level_value = std::stoi(row[0]);
+ node_value = std::stoi(row[1]);
+ parent_value = std::stoi(row[2]);
+ _Mat_Id.push_back(std::stoi(row[3]));
+
+ _LevelNode[level_value].push_back(node_value);
+ _NodeLevel.push_back(level_value);
+ _Parent.push_back(parent_value);
+
+ if(level_value < _TotalLevel){
+ std::vector<int> ch;
+ _Child.push_back(ch);
+ }
+ if(parent_value != node_value) _Child[parent_value].push_back(node_value);
+ }
+ else if(DataId == 3){ // read the DMN Parameter -----------------------
+ _Dim = std::stoi(row[0]);
+ _Vf = std::stod(row[1]);
+ }
+ else if(DataId == 4){
+ if(_Dim==2){
+ ParaN[0] = std::stod(row[0]);
+ ParaN[1] = 0.0;
+ }
+ else if(_Dim==3){
+ ParaN[0] = std::stod(row[0]);
+ ParaN[1] = std::stod(row[1]);
+ }
+ sin0 = sin(2.0*Pi*ParaN[0]);
+ cos0 = cos(2.0*Pi*ParaN[0]);
+ sin1 = sin(Pi*ParaN[1]);
+ cos1 = cos(Pi*ParaN[1]);
+ Para_Norm[0] = cos1*sin0;
+ Para_Norm[1] = cos1*cos0;
+ Para_Norm[2] = sin1;
+ _Para_Norm.push_back(Para_Norm);
+ _ParaAngle(pn,0) = ParaN[0];
+ _ParaAngle(pn,1) = ParaN[1];
+ pn +=1;
+ }
+ else if(DataId == 5){
+ Para_Wt[0] = std::stod(row[0]);
+ Para_Wt[1] = std::stod(row[1]);
+ _Para_Wt.push_back(Para_Wt);
+ }
+ }
+ }
+
+ int row = _col_Di*_NLeaf;
+ int col = _row_Di*_NInterface;
+ resizeFlag = _Nv.resize(row, col, true);
+ resizeFlag = _I.resize(row, _col_Di, true);
+
+ row = _col_Di*(_N_Node-1);
+ col = _col_Di*_NLeaf;
+ resizeFlag = _V.resize(row, col, true);
+
+ row = _row_Di*_NInterface;
+ resizeFlag = _NT.resize(row, _col_Di*(_N_Node-1), true);
+ resizeFlag = _NTV.resize(row, col, true);
+
+ resizeFlag = _C.resize(col, col,true);
+ resizeFlag = Jacobian.resize(row, row,true);
+ _NPara_Norm = _NInterface;
+ _NPara_Wt = static_cast<int>(_Para_Wt.size());
+
+ Msg::Info("Finish reading Tree Data and parameter file, level = %d, Dim = %d, NLeaf = %d", _TotalLevel, _Dim, _NLeaf);
+}
+
+void StochTMDMNDG3DMaterialLaw::reset_Parameter(std::vector<std::vector<double>>& Para_Norm, std::vector<std::vector<double>>& Para_Wt, std::vector<double>& Para_Alpha, const double Vf){
+
+ double sin0, cos0, sin1,cos1;
+ double Pi(3.14159265359);
+ for(int i=0;i<_NInterface;i++){
+ if(_Dim==2){
+ sin1 = 0.0;
+ cos1 = 1.0;
+ _ParaAngle(i,0) = Para_Norm[i][0];
+ }
+ else if(_Dim==3){
+ sin1 = sin(Pi*Para_Norm[i][1]);
+ cos1 = cos(Pi*Para_Norm[i][1]);
+ _ParaAngle(i,0) = Para_Norm[i][0];
+ _ParaAngle(i,1) = Para_Norm[i][1];
+ }
+ sin0 = sin(2.0*Pi*Para_Norm[i][0]);
+ cos0 = cos(2.0*Pi*Para_Norm[i][0]);
+ _Para_Norm[i][0] = cos1*sin0;
+ _Para_Norm[i][1] = cos1*cos0;
+ _Para_Norm[i][2] = sin1;
+ }
+
+ for(int i=0;i<_NPara_Wt;i++){
+ _Para_Wt[i][0] = Para_Wt[i][0];
+ _Para_Wt[i][1] = Para_Wt[i][1];
+ if(_NodeLevel[_Child[i][1]]==_TotalLevel) _Para_Wt[i][1] = 1.0;
+ }
+ for(int i=0;i<Para_Alpha.size();i++){
+ _Para_Alpha[i] = _Para_Alpha[i];
+ _Para_Alpha[i] = _Para_Alpha[i];
+ }
+ _VA.clear();
+ _VI.clear();
+ _Nv.setAll(0.0);
+ _NTV.setAll(0.0);
+ _NT.setAll(0.0);
+ _V.setAll(0.0);
+ _I.setAll(0.0);
+ _C.setAll(0.0);
+ Jacobian.setAll(0.0);
+ if(Vf!= 0.0){
+ _Vf = Vf;
+ }
+ StochTMDMNDG3DMaterialLaw::fill_Matrices();
+ //Msg::Info("Finish update parameters of DMN ");
+}
+
+void StochTMDMNDG3DMaterialLaw::reset_Parameter(const char* Para){
+ _VA.clear();
+ _VI.clear();
+ _Para_Norm.clear();
+ _Para_Wt.clear();
+ _Para_Alpha.clear();
+ _LevelNode.clear();
+ _NodeLevel.clear();
+ _Mat_Id.clear();
+ _Parent.clear();
+ _Child.clear();
+ if(_ReadTree){
+ StochTMDMNDG3DMaterialLaw::read_TreeData(Para);
+ }
+ else{
+ StochTMDMNDG3DMaterialLaw::read_parameter(Para);
+ }
+ StochTMDMNDG3DMaterialLaw::fill_Matrices();
+}
+
+void StochTMDMNDG3DMaterialLaw::stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff, const bool checkfrac, const bool dTangent){
+
+ // get ipvariable
+ StochTMDMNDG3DIPVariable* ipvcur = static_cast<StochTMDMNDG3DIPVariable*>(ipv);
+ const StochTMDMNDG3DIPVariable* ipv_prev = static_cast<const StochTMDMNDG3DIPVariable*>(ipvprev);
+ const STensor3& F = ipvcur->getConstRefToDeformationGradient();
+ const SVector3& H = ipvcur->getConstRefToGradT();
+ const double T = ipvcur->getConstRefToTemperature();
+
+ STensor3& P = ipvcur->getRefToFirstPiolaKirchhoffStress(); STensorOperation::zero(P);
+ SVector3& Q = ipvcur->getRefToThermalFlux(); STensorOperation::zero(Q);
+ double& w = ipvcur->getRefToThermalSource(); w = 0.;
+ double& Cp = ipvcur->getRefToExtraDofFieldCapacityPerUnitField()(0);
+ double& mechSource = ipvcur->getRefToMechanicalSource()(0); mechSource = 0.;
+
+ STensor43& L = ipvcur->getRefToTangentModuli(); STensorOperation::zero(L);
+ STensor3& dPdT = ipvcur->getRefTodPdT(); STensorOperation::zero(dPdT);
+ STensor3& dqdgradT = ipvcur->getRefTodThermalFluxdGradT(); STensorOperation::zero(dqdgradT);
+ SVector3& dqdT = ipvcur->getRefTodThermalFluxdT(); STensorOperation::zero(dqdT);
+ STensor33& dqdF = ipvcur->getRefTodThermalFluxdF(); STensorOperation::zero(dqdF);
+ double& dwdt = ipvcur->getRefTodThermalSourcedField();
+ STensor3& dwdf = ipvcur->getRefTodThermalSourcedF(); STensorOperation::zero(dwdf);
+ double& dmechSourcedt = ipvcur->getRefTodMechanicalSourcedField()(0,0); dmechSourcedt = 0.;
+ STensor3& dmechSourcedf = ipvcur->getRefTodMechanicalSourcedF()[0]; STensorOperation::zero(dmechSourcedf);
+ double dCpdT(0.);
+ STensor3 dCpdF; STensorOperation::zero(dCpdF);
+ STensor33 dPdgradT; STensorOperation::zero(dPdgradT);
+
+ // The tangent C_hom is combined with thermal tangents.
+ fullMatrix<double> C_hom(_col_Di,_col_Di);
+ // if(C_hom.size1()!=_col_Di || C_hom.size2()!=_col_Di){C_hom.resize(_col_Di,_col_Di);};
+ C_hom.setAll(0.);
+ fullVector<double> PQ_hom(_col_Di);
+ // if(PQ_hom.size()!=_col_Di){PQ_hom.resize(_col_Di);};
+ PQ_hom.setAll(0.);
+ fullVector<double> dPQdT_hom(_col_Di);
+ // if(dPQdT_hom.size()!=_col_Di){dPQdT_hom.resize(_col_Di);};
+ dPQdT_hom.setAll(0.);
+ double Cp_hom(0.), thermSrc_hom(0.), mechSource_hom(0.), dwdt_hom(0.), dmechSourcedt_hom(0.);
+ fullVector<double> dwdf_hom(9), dmechSourcedf_hom(9); dwdf_hom.setAll(0.); dmechSourcedf_hom.setAll(0.);
+
+ fullVector<double>& FHvec = ipvcur->getRefToCombinedGradientVect();
+ fullVector<double>& PQvec = ipvcur->getRefToCombinedStressVect(); PQvec.setAll(0.0);
+ fullVector<double>& ab_cur = ipvcur->getRefToCombinedInterfaceFluctuationVect();
+ fullVector<double>& Cpvec = ipvcur->getRefToCpVect(); Cpvec.setAll(0.0);
+ fullVector<double>& thermSrcvec = ipvcur->getRefToThermalSourceVect(); thermSrcvec.setAll(0.0);
+ fullVector<double>& mechSrcvec = ipvcur->getRefToMechSourceVect(); mechSrcvec.setAll(0.0);
+ fullVector<double> dwdtvec, dmechSrcdTvec, dwdFvec, dmechSrcdFvec;
+ dwdtvec.resize(_NLeaf); dwdtvec.setAll(0.); dmechSrcdTvec.resize(_NLeaf); dmechSrcdTvec.setAll(0.);
+ dwdFvec.resize(9*_NLeaf); dwdFvec.setAll(0.); dmechSrcdFvec.resize(9*_NLeaf); dmechSrcdFvec.setAll(0.);
+
+ fullVector<double> Res_node(_row_Di*_NInterface);
+ // if(Res_node.size()!=_row_Di*_NInterface){Res_node.resize(_row_Di*_NInterface);};
+ Res_node.setAll(0.);
+ fullVector<double> ab_step(_row_Di*_NInterface);
+ // if(ab_step.size()!=_row_Di*_NInterface){ab_step.resize(_row_Di*_NInterface);};
+ ab_step.setAll(0.0);
+ fullMatrix<double> NTVC(_row_Di*_NInterface, _col_Di*_NLeaf);
+ // if(NTVC.size1()!=_row_Di*_NInterface || NTVC.size2()!=_col_Di*_NLeaf){NTVC.resize(_row_Di*_NInterface, _col_Di*_NLeaf);};
+ NTVC.setAll(0.0);
+ fullMatrix<double> Jacobian_inv(_row_Di*_NInterface, _row_Di*_NInterface);
+ // if(Jacobian_inv.size1()!=_row_Di*_NInterface || Jacobian_inv.size2()!=_row_Di*_NInterface){Jacobian_inv.resize(_row_Di*_NInterface, _row_Di*_NInterface);};
+ Jacobian_inv.setAll(0.0);
+ fullMatrix<double> tmp(_col_Di*_NLeaf, _col_Di);
+ // if(tmp.size1()!=_col_Di*_NLeaf || tmp.size2()!=_col_Di){tmp.resize(_col_Di*_NLeaf, _col_Di);};
+ tmp.setAll(0.0);
+ fullMatrix<double> dRdF(_row_Di*_NInterface, _col_Di);
+ // if(dRdF.size1()!=_row_Di*_NInterface || dRdF.size2()!=_col_Di){dRdF.resize(_row_Di*_NInterface, _col_Di);};
+ dRdF.setAll(0.0);
+ fullMatrix<double> dadF(_row_Di*_NInterface, _col_Di);
+ // if(dadF.size1()!=_row_Di*_NInterface || dadF.size2()!=_col_Di){dadF.resize(_row_Di*_NInterface, _col_Di);};
+ dadF.setAll(0.0);
+
+ _C.setAll(0.0);
+ int ite = 0;
+ double r = 1.0;
+ int pos_start = _NInterface;
+ int max_iter= 15;
+ bool stiff_loc = true;
+
+
+ while(r>_tol and ite < max_iter){
+ ite +=1;
+ _Nv.mult(ab_cur, FHvec);
+ for(int i=0; i<_NLeaf; i++){
+
+ // Msg::Error("i is %d", i);
+
+ if (_VA[pos_start+i] > 1.e-6){
+ ThermoMechanicsDG3DIPVariableBase* ipv_i = dynamic_cast<ThermoMechanicsDG3DIPVariableBase*>(ipvcur->getIPv(i));
+ const ThermoMechanicsDG3DIPVariableBase* ipvprev_i = dynamic_cast<const ThermoMechanicsDG3DIPVariableBase*>(ipv_prev->getIPv(i));
+ // std::cout << "Object type: " << typeid(*ipv_i).name() << std::endl;
+
+ STensor3& Floc = ipv_i->getRefToDeformationGradient();
+ SVector3& Hloc = ipv_i->getRefToGradT();
+ double& Tloc = ipv_i->getRefToTemperature();
+
+ if (_flag_isothermal && _flag_microTempFixed){
+ Tloc = T;
+ Hloc = H;
+ for(int m=0; m<3; m++){
+ for(int n=0; n<3; n++){
+ Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*_row_Di);
+ }
+ }
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ Tloc = T;
+ for(int m=0; m<3; m++){
+ Hloc(m) = H(m) + FHvec(i*_col_Di+ 9+m);
+ for(int n=0; n<3; n++){
+ Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*3);
+ }
+ }
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){
+ Tloc = T + FHvec((i+1)*_col_Di-1);
+ for(int m=0; m<3; m++){
+ Hloc(m) = H(m) + FHvec(i*_col_Di+ 9+m);
+ for(int n=0; n<3; n++){
+ Floc(m,n) = F(m,n) + FHvec(i*_col_Di + m + n*3);
+ }
+ }
+ }
+ else{
+ Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
+ }
+
+ int mat_IP = _Mat_Id[_NInterface + i];
+ _mapLaw[mat_IP]->stress(ipv_i, ipvprev_i, stiff_loc, checkfrac, dTangent);
+
+ const STensor3& P_i = ipv_i->getConstRefToFirstPiolaKirchhoffStress();
+ const SVector3& Q_i = ipv_i->getConstRefToThermalFlux();
+ const double& w_i = ipv_i->getConstRefToThermalSource();
+ const double& Cp_i = ipv_i->getConstRefToExtraDofFieldCapacityPerUnitField()(0);
+ const double& mechSource_i = ipv_i->getConstRefToMechanicalSource()(0);
+
+ if (_flag_isothermal && _flag_microTempFixed){
+ Cpvec(i) = Cp_i; thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i;
+ for(int m=0; m<3; m++){
+ for(int n=0; n<3; n++){
+ PQvec(i*_col_Di + m + n*_row_Di) = P_i(m,n);
+ }
+ }
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ Cpvec(i) = Cp_i; thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i;
+ for(int m=0; m<3; m++){
+ PQvec(i*_col_Di+ 9+m) = Q_i(m);
+ for(int n=0; n<3; n++){
+ PQvec(i*_col_Di + m + n*3) = P_i(m,n);
+ }
+ }
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){
+ thermSrcvec(i) = w_i; mechSrcvec(i) = mechSource_i;
+ PQvec((i+1)*_col_Di-1) = Cp_i;
+ for(int m=0; m<3; m++){
+ PQvec(i*_col_Di+ 9+m) = Q_i(m);
+ for(int n=0; n<3; n++){
+ PQvec(i*_col_Di + m + n*3) = P_i(m,n);
+ }
+ }
+ }
+ else{
+ Msg::Error("This combination of flags doesn't exist. Redefine in the constructor.");
+ }
+
+
+ if(stiff_loc){
+ const STensor43& L_i = ipv_i->getConstRefToTangentModuli();
+ const STensor3& dPdT_i = ipv_i->getRefTodPdT();
+ const STensor3& dqdgradT_i = ipv_i->getRefTodThermalFluxdGradT();
+ const SVector3& dqdT_i = ipv_i->getRefTodThermalFluxdT();
+ const STensor33& dqdF_i = ipv_i->getRefTodThermalFluxdF();
+ const double& dwdt_i = ipv_i->getRefTodThermalSourcedField();
+ const STensor3& dwdf_i = ipv_i->getRefTodThermalSourcedF(); STensorOperation::zero(dwdf);
+ const double& dmechSourcedt_i = ipv_i->getRefTodMechanicalSourcedField()(0,0);
+ const STensor3& dmechSourcedf_i = ipv_i->getRefTodMechanicalSourcedF()[0];
+
+ if (_flag_isothermal && _flag_microTempFixed){
+ dwdtvec(i) = 0.; dmechSrcdTvec(i) = 0.;
+ for(int m=0; m<3; m++)
+ for(int n=0; n<3; n++){
+ dwdFvec(i*_col_Di+m+_row_Di*n) = 0.;
+ dmechSrcdFvec(i*_col_Di+m+_row_Di*n) = 0.;
+ for(int k=0; k<3; k++)
+ for(int l=0; l<3; l++)
+ _C(i*_col_Di+m+_row_Di*n, i*_col_Di+k+_row_Di*l) = L_i(m,n,k,l);
+ }
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ dwdtvec(i) = dwdt_i; dmechSrcdTvec(i) = dmechSourcedt_i;
+ for(int m=0; m<3; m++)
+ for(int n=0; n<3; n++){
+ _C(i*_col_Di + 9+m, i*_col_Di + 9+n) = dqdgradT_i(m,n);
+ dwdFvec(i*9+m+_row_Di*n) = dwdf_i(m,n);
+ dmechSrcdFvec(i*9+m+_row_Di*n) = dmechSourcedf_i(m,n);
+ for(int k=0; k<3; k++){
+ // _C(i*_col_Di+m+_row_Di*n, i*_col_Di + (k+1)*_row_Di-1) = 0.; // This is for dPdgradT_i
+ _C(i*_col_Di + 9+m, i*_col_Di+n+3*k) = dqdF_i(m,n,k);
+ for(int l=0; l<3; l++)
+ _C(i*_col_Di+m+3*n, i*_col_Di+k+3*l) = L_i(m,n,k,l);
+ }
+ }
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){
+ dwdtvec(i) = dwdt_i; dmechSrcdTvec(i) = dmechSourcedt_i;
+ for(int m=0; m<3; m++){
+ _C(i*_col_Di + (m+1)*(_row_Di-1)-1, (i+1)*_col_Di-1) = dqdT_i(m);
+ for(int n=0; n<3; n++){
+ _C(i*_col_Di + 9+m, i*_col_Di + 9+m) = dqdgradT_i(m,n);
+ _C(i*_col_Di + m + 3*n, (i+1)*_col_Di-1) = dPdT_i(m,n);
+ dwdFvec(i*9+m+_row_Di*n) = dwdf_i(m,n);
+ dmechSrcdFvec(i*9+m+_row_Di*n) = dmechSourcedf_i(m,n);
+ for(int k=0; k<3; k++){
+ // _C(i*_col_Di+m+(_row_Di-1)*n, i*_col_Di + (k+1)*(_row_Di-1)-1) = 0.; // This is for dPdgradT_i
+ _C(i*_col_Di + 9+m, i*_col_Di+n+3*k) = dqdF_i(m,n,k);
+ for(int l=0; l<3; l++)
+ _C(i*_col_Di+m+3*n, i*_col_Di+k+3*l) = L_i(m,n,k,l);
+ }
+ }
+ }
+ }
+
+ }
+ }
+ // Msg::Error("Sizes of PQVec are %d", PQvec.size());
+ }
+ // Msg::Error("Sizes of PQVec are %d", PQvec.size());
+ Msg::Error("Sizes of _NTV are %d and %d", _NTV.size1(),_NTV.size2());
+ _NTV.mult(PQvec, Res_node);
+ r = Res_node.norm()/_NInterface;
+
+ _NTV.mult(_C, NTVC);
+ NTVC.mult(_Nv, Jacobian);
+ Jacobian.invertInPlace();
+
+ Jacobian.mult(Res_node, ab_step);
+ ab_step.scale(-1.0);
+ ab_cur.axpy(ab_step, 1.0);
+ }
+
+ Jacobian.scale(-1.0);
+ int pos0 = _Child[0][0];
+ int pos1 = _Child[0][1];
+ if(r>_tol){
+ Msg::Error("TMDMN residual (= %lf) didn't converge to tolenerce after %d iteration",r,ite);
+ return;
+ }
+ else{
+ for(int m=0;m<_col_Di;m++){
+ PQ_hom(m) = 0.0;
+ for(int i=0;i<_NLeaf;i++){
+ PQ_hom(m) += (_VA[pos0]*_V((pos0-1)*_col_Di+m, _col_Di*i+m) + _VA[pos1]*_V((pos1-1)*_col_Di+m, _col_Di*i+m))*PQvec(_col_Di*i+m); // Stress average
+ }
+ }
+
+ double temp = 0.;
+ for(int i=0;i<_NLeaf;i++){
+ temp = (_VA[pos0]*_V(0, _col_Di*i) + _VA[pos1]*_V(_col_Di, _col_Di*i));
+ Cp_hom += temp*Cpvec(i); // Cp average
+ thermSrc_hom += temp*thermSrcvec(i); // w average
+ mechSource_hom += temp*mechSrcvec(i); // mechSrc average
+ }
+
+ if (_flag_isothermal && _flag_microTempFixed){
+ Cp = Cp_hom; w = 0.; mechSource = 0.;
+ for(int i=0;i<3;i++){
+ Q(i) = 0.;
+ for(int j=0;j<3;j++){
+ P(i,j) = PQ_hom(i+_row_Di*j);
+ }
+ }
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ Cp = Cp_hom; w = thermSrc_hom; mechSource = mechSource_hom;
+ for(int i=0;i<3;i++){
+ Q(i) = PQ_hom(9+i);
+ for(int j=0;j<3;j++){
+ P(i,j) = PQ_hom(i+3*j);
+ }
+ }
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){
+ Cp = PQ_hom(_col_Di-1); w = thermSrc_hom; mechSource = mechSource_hom;
+ for(int i=0;i<3;i++){
+ Q(i) = PQ_hom(9+i);
+ for(int j=0;j<3;j++){
+ P(i,j) = PQ_hom(i+3*j);
+ }
+ }
+ }
+
+ if(stiff){
+
+ NTVC.mult(_I, dRdF);
+ Jacobian.mult(dRdF, dadF);
+ _Nv.mult(dadF, tmp);
+ tmp.add(_I);
+
+ for(int m=0;m<_col_Di;m++){
+ for(int n=0;n<_col_Di;n++){
+ C_hom(m,n) = 0.0;
+ for(int i=0;i<_NLeaf;i++){
+ for(int j=0;j< _col_Di*_NLeaf;j++){
+ C_hom(m,n) += (_VA[pos0]*_V((pos0-1)*_col_Di+m,m+_col_Di*i) + _VA[pos1]*_V((pos1-1)*_col_Di+m,m+_col_Di*i))*_C(m+_col_Di*i,j)*tmp(j,n);
+ }
+ }
+ }
+ }
+
+ double temp(0.), dTidT(1.); // = Cp_hom/(temp*Cpvec(i))
+ for(int i=0;i<_NLeaf;i++){
+ temp = (_VA[pos0]*_V(0, _col_Di*i) + _VA[pos1]*_V(_col_Di, _col_Di*i));
+ dwdt_hom += temp*dwdtvec(i)*dTidT; // dwdT average
+ dmechSourcedt_hom += temp*dmechSrcdTvec(i)*dTidT; // dmechSourcedt average
+ for(int m=0;m<9;m++)
+ for(int n=0;n<9;n++){
+ dwdf_hom(n) += temp*dwdFvec(m+9*i)*tmp(9*i,n);
+ dmechSourcedf_hom(n) += temp*dmechSrcdFvec(m+9*i)*tmp(9*i,n);
+ }
+ }
+
+ if (_flag_isothermal && _flag_microTempFixed){
+ for(int m=0; m<3; m++)
+ for(int n=0; n<3; n++)
+ for(int k=0; k<3; k++)
+ for(int l=0; l<3; l++)
+ L(m,n,k,l) = C_hom(m+_row_Di*n,k+_row_Di*l);
+ }
+ else if(!_flag_isothermal && _flag_microTempFixed){
+ dwdt = dwdt_hom; dmechSourcedt = dmechSourcedt_hom;
+ for(int m=0; m<3; m++)
+ for(int n=0; n<3; n++){
+ dqdgradT(m,n) = C_hom(9+m,9+n);
+ dwdf(m,n) = dwdf_hom(m+3*n);
+ dmechSourcedf(m,n) = dmechSourcedf_hom(m+3*n);
+ for(int k=0; k<3; k++){
+ dqdF(m,n,k) = C_hom(9+m,n+3*k);
+ // dPdgradT(m,n,l) = C_hom[m+_row_Di*n][(k+1)*_row_Di-1]; // 0. This is for dPdgradT
+ for(int l=0; l<3; l++)
+ L(m,n,k,l) = C_hom(m+3*n,k+3*l);
+ }
+ }
+ }
+ else if(!_flag_isothermal && !_flag_microTempFixed){
+ for(int m=0; m<3; m++){
+ dqdT(m) = C_hom( (m+1)*(_row_Di-1)-1, _col_Di-1 );
+ for(int n=0; n<3; n++){
+ dPdT(m,n) = C_hom( m + n*(_row_Di-1), _col_Di-1 );
+ dqdgradT(m,n) = C_hom( 9+m, 9+n);
+ dwdf(m,n) = dwdf_hom(m+3*n);
+ dmechSourcedf(m,n) = dmechSourcedf_hom(m+3*n);
+ for(int k=0; k<3; k++){
+ dqdF(m,n,k) = C_hom(9+m, n+3*k);
+ // C_hom[m+(_row_Di-1)*n][(k+1)*(_row_Di-1)-1] = 0.; This is for dPdgradT
+ for(int l=0; l<3; l++)
+ L(m,n,k,l) = C_hom(m+3*n, k+3*l);
+ }
+ }
+ }
+ }
+ }
+ }
+ Msg::Error("Sizes of PQVec are %d", PQvec.size());
+ ipvcur->setRefToDGElasticTangentModuli(this->elasticStiffness);
+ Msg::Error("Sizes");
+ Msg::Error("Sizes of PQVec are %d", PQvec.size());
+}
+
+void StochTMDMNDG3DMaterialLaw::setTime(const double t,const double dtime){
+ dG3DMaterialLaw::setTime(t,dtime);
+ for (std::vector<dG3DMaterialLaw*>::iterator it = _mapLaw.begin(); it != _mapLaw.end(); it++) {
+ (*it)->setTime(t,dtime);
+ }
+}
+
+// end of stress
+
+//
J2SmallStrainDG3DMaterialLaw::J2SmallStrainDG3DMaterialLaw(const int num,const double rho,
double E,const double nu,const double sy0,const double h,
const double tol,const bool tangentByPert, const double pert) : dG3DMaterialLaw(num,rho,true),
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index a1ef00808fe2fdcd0021e214799e34bce3621928..a059ea83f427fefe42e9533bfb256401ea48810e 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -697,20 +697,21 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
int _NLeaf;
int _NInterface;
int _N_Node;
- int _NPara_Wt;
- int _NPara_Norm;
+ int _NPara_Wt;
+ int _NPara_Norm;
int _Dim;
double _Vf;
double _tol;
+ double _alpha, _beta, _gamma; // Euler Angles
std::vector<std::vector<int>> _LevelNode; // keep the Nodes at each level
std::vector<int> _NodeLevel; // _NodeLevel[Node] = Level of Node
std::vector<int> _Mat_Id;
std::vector<int> _Parent; //Parent[Node] = Parent of Node
std::vector<std::vector<int>> _Child; // Child[Node][0] = left child ; Child[Node][1] = right child of Node
-
+
- std::vector<SPoint2> _Para_Wt;
- std::vector<SPoint3> _Para_Norm;
+ std::vector<SPoint2> _Para_Wt; // vector of weights
+ std::vector<SPoint3> _Para_Norm; // vector of normals
fullMatrix<double> _ParaAngle;
std::vector<dG3DMaterialLaw*> _mapLaw;
fullMatrix<double> _Nv;
@@ -722,12 +723,14 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
fullMatrix<double> _C;
fullMatrix<double> Jacobian;
-
- void fill_NLocal(const int pos, double N[][3])const;
+
+ bool _flag_isothermal;
+
+ virtual void fill_NLocal(const int pos, double N[][3])const;
void fill_W_WI(const int pos, double W[][2]);
- void fill_Matrices();
- void read_parameter(const char *ParaFile);
- void read_TreeData(const char *TreeFile);
+ virtual void fill_Matrices();
+ virtual void read_parameter(const char *ParaFile);
+ virtual void read_TreeData(const char *TreeFile);
virtual void getLeafOfInterface(const int pos, int &node_start, int &node_end) const;
virtual void TwoPhasesInteract(const std::vector<STensor43>& C, fullMatrix<double> &mat, const int pos);
virtual void getdVAdP(std::vector< std::vector <SPoint2> >& dVAdPv) const;
@@ -735,9 +738,10 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
virtual void get_matrixdVdPv(std::vector< std::vector <SPoint2> >& dVAdPv, std::vector<std::vector<std::vector<std::vector<double>>>>& d_VdPv) const;
virtual void get_matrixdNdP(std::vector< std::vector <SPoint2> >& dVAdPv, std::vector<std::vector<std::vector<std::vector<double>>>>& d_NvdPv, std::vector<std::vector<std::vector<std::vector<double>>>>& d_NvdPn, std::vector<std::vector<std::vector<std::vector<double>>>>& d_NTdPn) const;
virtual void get_matrixdNTdPv(std::vector<std::vector<std::vector<std::vector<double>>>>& d_NTdPv) const;
+ virtual void eulerZXZToRotationMatrix(const double alpha, const double beta, const double gamma, STensor3& rotationMatrix);
#endif //SWIG
public:
- StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree=false, const bool porous=false,const double tol=1e-6);
+ StochDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree=false, const bool porous=false,const double tol=1e-6,const bool flag_isothermal=true);
void addLaw(dG3DMaterialLaw* law);
#ifndef SWIG
StochDMNDG3DMaterialLaw(const StochDMNDG3DMaterialLaw &source);
@@ -770,10 +774,63 @@ class StochDMNDG3DMaterialLaw : public dG3DMaterialLaw{
virtual void setMacroSolver(const nonLinearMechSolver* sv);
virtual void getParameterDerivative(const IPVariable*ipv, std::vector< fullMatrix<double> >& dPdPn, std::vector< fullMatrix<double> >&dPdPv) const;
virtual void reset_Parameter(std::vector<std::vector<double>>& Para_Norm, std::vector<std::vector<double>>& Para_Wt, const double Vf = 0.0);
- virtual void reset_Parameter(const char* Para);
+ virtual void reset_Parameter(const char* Para);
+ virtual void setZXZRotationMatrix_in_IP(StochDMNDG3DIPVariable* ipv);
#endif //SWIG
+ virtual void setEulerAngles_rotation_DMN_response(const double alpha, const double beta, const double gamma){_alpha = alpha; _beta = beta; _gamma = gamma;};
};
+// FLE
+class StochTMDMNDG3DMaterialLaw : public StochDMNDG3DMaterialLaw{
+ protected:
+ #ifndef SWIG
+ int _col_Di, _row_Di;
+ bool _flag_microTempFixed;
+ std::vector<double> _Para_Alpha; // vector of coefficients
+ void fill_Matrices();
+ void read_parameter(const char *ParaFile);
+ void read_TreeData(const char *TreeFile);
+ void fill_NLocal(const int pos, fullMatrix<double> N) const;
+ #endif //SWIG
+ public:
+ StochTMDMNDG3DMaterialLaw(const int num, const double rho, const double E,const double nu, const char *ParaFile, const bool ReadTree=false,
+ const bool porous = false, const bool flag_isothermal = false, const bool flag_microTempFixed = true, const double tol=1e-6);
+ #ifndef SWIG
+ StochTMDMNDG3DMaterialLaw(const StochTMDMNDG3DMaterialLaw &source);
+ virtual ~StochTMDMNDG3DMaterialLaw();
+
+ virtual void setTime(const double t,const double dtime);
+ virtual materialLaw::matname getType() const {return materialLaw::StochTMDMN;}
+
+ 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 initialIPVariable(IPVariable* ipv, bool stiff);
+ virtual void initLaws(const std::map<int,materialLaw*> &maplaw);
+ virtual double scaleFactor() const {return 1.;};
+ virtual double soundSpeed() const {return 0.;};
+ virtual materialLaw* clone() const {return new StochTMDMNDG3DMaterialLaw(*this);};
+ virtual void stress(IPVariable*ipv, const IPVariable*ipvprev, const bool stiff=true, const bool checkfrac=true, const bool dTangent=false);
+
+ virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const;
+ virtual bool withEnergyDissipation() const {return false;};
+ virtual const materialLaw *getConstNonLinearSolverMaterialLaw() const
+ {
+ Msg::Error("getConstNonLinearSolverMaterialLaw in StochTMDMNDG3DMaterialLaw does not exist");
+ return NULL;
+ }
+ virtual materialLaw *getNonLinearSolverMaterialLaw()
+ {
+ Msg::Error("getNonLinearSolverMaterialLaw in StochTMDMNDG3DMaterialLaw does not exist");
+ return NULL;
+ }
+ virtual void setMacroSolver(const nonLinearMechSolver* sv);
+ virtual void reset_Parameter(const char* Para);
+ virtual void reset_Parameter(std::vector<std::vector<double>>& Para_Norm, std::vector<std::vector<double>>& Para_Wt, std::vector<double>& Para_Alpha, const double Vf);
+ #endif //SWIG
+ virtual void setEulerAngles_rotation_DMN_response(const double alpha, const double beta, const double gamma){_alpha = alpha; _beta = beta; _gamma = gamma;};
+ };
+// FLE
+
class J2SmallStrainDG3DMaterialLaw : public dG3DMaterialLaw
{
protected: