add Sachs, Taylor, BetaRule and Idiart schemes

generic-behaviours/homogenization/BetaRule.mfront

@@ -0,0 +1,163 @@
+@DSL ImplicitII;
+@Behaviour BetaRule;
+@Author Antoine Martin;
+
+@Algorithm NewtonRaphson;
+@Epsilon 1.e-14;
+@Theta 1;
+
+@ModellingHypothesis Tridimensional;
+@PhysicalBounds T in [0. : *[;
+
+//parameters of beta-rule
+@Parameter real c = 50e3;
+@Parameter real DD = 1e2;
+@Parameter stress sig_0=1e4;
+
+@MaterialProperty real f;
+f.setEntryName("FirstPhaseFraction");
+
+
+@StateVariable Stensor Sig;
+Sig.setEntryName("MacroscopicStress");
+@StateVariable Stensor eto1;
+eto1.setEntryName("FirstPhaseTotalStrain");
+@StateVariable Stensor eto2;
+eto2.setEntryName("SecondPhaseTotalStrain");
+@StateVariable Stensor beta1;
+beta1.setEntryName("FirstPhaseBetaStrain");
+@StateVariable Stensor beta2;
+beta2.setEntryName("SecondPhaseBetaStrain");
+@AuxiliaryStateVariable Stensor evp1;
+evp1.setEntryName("FirstPhaseViscoPlasticStrain");
+@AuxiliaryStateVariable Stensor evp2;
+evp2.setEntryName("SecondPhaseViscoplasticStrain");
+
+
+@BehaviourVariable b1 {
+file: "MericCailletaud.mfront",
+variables_suffix: "1",
+store_gradients: false,
+external_names_prefix: "FirstPhase",
+shared_external_state_variables: {".+"}
+};
+
+@BehaviourVariable b2 {
+file: "MericCailletaud.mfront",
+variables_suffix: "2",
+store_gradients: false,
+external_names_prefix: "SecondPhase",
+shared_external_state_variables: {".+"}
+};
+
+@LocalVariable Stensor4 S1;
+@LocalVariable Stensor4 S2;
+@LocalVariable Stensor devp1;
+@LocalVariable Stensor devp2;
+
+@InitLocalVariables<Append> {
+    auto lambda1 = computeLambda(E_1,nu_1);
+    auto mu1 = computeMu(E_1,nu_1);
+    auto kap1= lambda1+real(2)/3*mu1;
+    auto J=st2tost2<3u,real>::J();
+    auto K=st2tost2<3u,real>::K();
+    S1=stress(1)/(3*kap1)*J+stress(1)/(2*mu1)*K;
+
+    auto lambda2 = computeLambda(E_2,nu_2);
+    auto mu2 = computeMu(E_2,nu_2);
+    auto kap2= lambda2+real(2)/3*mu2;
+    S2=stress(1)/(3*kap2)*J+stress(1)/(2*mu2)*K;
+}
+
+@Integrator {
+  constexpr const auto eeps = strain(1.e-14);
+  auto Id=Stensor4::Id();
+
+  const auto i1=initialize(b1);
+  if ((i1 != 0) && (i1 != 1)) {
+    return false;
+  }
+  b1.eto=eto1;
+  b1.deto=deto1;
+  constexpr auto b1_smflag = TangentOperatorTraits<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR>::STANDARDTANGENTOPERATOR;
+  const auto r1 = b1.integrate(b1_smflag,CONSISTENTTANGENTOPERATOR);
+  if ((r1 != 0) && (r1 != 1)) {
+    return false;
+  }
+  StiffnessTensor Dt1 = b1.getTangentOperator();
+
+  auto dsig1=b1.sig-sig1;
+  auto deel1=b1.eel-eel1;
+  devp1=deto1-deel1;
+  auto ndevp_1=sqrt(real(2)/3*devp1|devp1);
+  auto beta_mts_1=beta1+theta*dbeta1;
+
+  const auto i2=initialize(b2);
+  if ((i2 != 0) && (i2 != 1)) {
+    return false;
+  }
+  b2.eto=eto2;
+  b2.deto=deto2;
+  constexpr auto b2_smflag = TangentOperatorTraits<MechanicalBehaviourBase::STANDARDSTRAINBASEDBEHAVIOUR>::STANDARDTANGENTOPERATOR;
+  const auto r2 = b2.integrate(b2_smflag,CONSISTENTTANGENTOPERATOR);
+  if ((r2 != 0) && (r2 != 1)) {
+    return false;
+  }
+  StiffnessTensor Dt2 = b2.getTangentOperator();
+
+  auto dsig2=b2.sig-sig2;
+  auto deel2=b2.eel-eel2;
+  devp2=deto2-deel2;
+  real ndevp_2=sqrt(real(2)/3*devp2|devp2);
+  Stensor beta_mts_2=beta2+theta*dbeta2;
+
+  //residues
+  auto dB=(1-f)*dbeta1+f*dbeta2;
+  fSig = -deto+ (1-f)*deto1+f*deto2;
+  feto1 = (dsig1-dSig-c*(dB-dbeta1))/sig_0;
+  feto2 = (dsig2-dSig-c*(dB-dbeta2))/sig_0;
+  fbeta1 = dbeta1-devp1+DD*ndevp_1*beta_mts_1;
+  fbeta2 = dbeta2-devp2+DD*ndevp_2*beta_mts_2;
+
+
+  //jacobian matrix
+  auto Null = Stensor4{real{}};
+  dfSig_ddSig=Null;
+  dfSig_ddeto1 = (1-f)*Id;
+  dfSig_ddeto2 = f*Id;
+
+  dfeto1_ddSig = -Id/sig_0;
+  dfeto2_ddSig = -Id/sig_0;
+
+  dfeto1_ddeto1 = Dt1/sig_0;
+  dfeto2_ddeto2 = Dt2/sig_0;
+
+  dfeto1_ddbeta1 = c/sig_0*f*Id;
+  dfeto1_ddbeta2 = -c/sig_0*f*Id;
+  dfeto2_ddbeta1 = c/sig_0*(f-1)*Id;
+  dfeto2_ddbeta2 = c/sig_0*(1-f)*Id;
+
+
+  dfbeta1_ddeto1 = -(Id-S1*Dt1) + 2*DD/3*Stensor4(beta_mts_1^devp1)*Stensor4(Id-S1*Dt1)/max(ndevp_1,eeps);
+  dfbeta2_ddeto2 = -(Id-S2*Dt2) + 2*DD/3*Stensor4(beta_mts_2^devp2)*Stensor4(Id-S2*Dt2)/max(ndevp_2,eeps);
+
+  dfbeta1_ddbeta1 =Id+ theta*DD*ndevp_1*Id;
+  dfbeta2_ddbeta2 =Id+ theta*DD*ndevp_2*Id;
+}
+
+@UpdateAuxiliaryStateVariables{
+evp1+=devp1;
+evp2+=devp2;
+}
+
+@ComputeFinalStress{
+sig += dSig;
+}
+
+@TangentOperator{
+if (smt){
+  Stensor4 iJs;
+  getPartialJacobianInvert(iJs);
+  Dt=iJs;
+  }
+}

generic-behaviours/homogenization/BetaRule.mtest

@@ -0,0 +1,32 @@
+@ModellingHypothesis 'Tridimensional';
+@Behaviour<@interface@> @library@ @behaviour@;
+//@Behaviour<generic> 'src/libBehaviour.so' 'BetaRule';
+@MaterialProperty<constant> 'FirstPhaseFraction' 0.5;
+@MaterialProperty<constant> 'FirstPhaseE_'               208000;
+@MaterialProperty<constant> 'FirstPhasenu_'               0.3;
+@MaterialProperty<constant> 'FirstPhasemu_'               80000;
+@MaterialProperty<constant> 'SecondPhaseE_'               104000;
+@MaterialProperty<constant> 'SecondPhasenu_'               0.3;
+@MaterialProperty<constant> 'SecondPhasemu_'               40000;
+@MaterialProperty<constant> 'FirstPhasetau0'               36.;
+@MaterialProperty<constant> 'SecondPhasetau0'               25.;
+@MaterialProperty<constant> 'FirstPhasen'               10.;
+@MaterialProperty<constant> 'SecondPhasen'               5.;
+@ExternalStateVariable 'Temperature' 293.15;
+@ImposedStrain 'EXX' {0 : 0, 5 :1e-2};
+@ImposedStress 'SYY' {0 : 0, 1 : 0};
+@ImposedStress 'SZZ' {0 : 0, 1 : 0};
+@ImposedStress 'SXY' {0 : 0, 1 : 0};
+@ImposedStress 'SXZ' {0 : 0, 1 : 0};
+@ImposedStress 'SYZ' {0 : 0, 1 : 0};
+@Times {0, 3 in 400};
+
+//@OutputFilePrecision 14;
+
+@Test<file> @reference_file@ 'EXX' 2 1.e-14;
+@Test<file> @reference_file@ 'SXX' 8 1e-3;
+@Test<file> @reference_file@ 'SYY' 9 1e-3;
+@Test<file> @reference_file@ 'SZZ' 10 1e-3;
+@Test<file> @reference_file@ 'SXY' 11 1e-3;
+@Test<file> @reference_file@ 'SXZ' 12 1e-3;
+@Test<file> @reference_file@ 'SYZ' 13 1e-3;

generic-behaviours/homogenization/CMakeLists.txt

@@ -1,7 +1,27 @@
 mfront_behaviours_library(Homogenization
   PonteCastaneda1992
+  Taylor
+  Sachs
+  BetaRule
+  Idiart_viscoelastic
   )
 
 genericmtest(Homogenization PonteCastaneda1992
   BEHAVIOUR PonteCastaneda1992
   REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/PonteCastaneda1992.ref)
+
+genericmtest(Homogenization Taylor
+  BEHAVIOUR Taylor
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Taylor.ref)
+
+genericmtest(Homogenization Sachs
+  BEHAVIOUR Sachs
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Sachs.ref)
+
+genericmtest(Homogenization BetaRule
+  BEHAVIOUR BetaRule
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/BetaRule.ref)
+
+genericmtest(Homogenization Idiart_viscoelastic
+  BEHAVIOUR Idiart_viscoelastic
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Idiart_viscoelastic.ref)

generic-behaviours/homogenization/Idiart_viscoelastic.mfront

@@ -0,0 +1,113 @@
+@DSL Default;
+@Behaviour Idiart_viscoelastic;
+@Author Martin Antoine;
+@Date 06 / 05 / 25;
+@UseQt false;
+
+@TFELLibraries {"Material"};
+@Includes{
+#include "TFEL/Material/IsotropicEshelbyTensor.hxx"
+#include "TFEL/Material/StiffnessTensor.hxx"
+#include "TFEL/Material/IsotropicModuli.hxx"}
+
+@MaterialProperty real ka0;
+ka0.setEntryName("ka0");
+@MaterialProperty real mu0;
+mu0.setEntryName("mu0");
+@MaterialProperty real kav0;
+kav0.setEntryName("kav0");
+@MaterialProperty real muv0;
+muv0.setEntryName("muv0");
+@MaterialProperty real kar;
+kar.setEntryName("kar");
+@MaterialProperty real mur;
+mur.setEntryName("mur");
+@MaterialProperty real kavr;
+kavr.setEntryName("kavr");
+@MaterialProperty real muvr;
+muvr.setEntryName("muvr");
+@MaterialProperty real tau0;
+tau0.setEntryName("tau0");
+@MaterialProperty real fr;
+fr.setEntryName("fr");
+@MaterialProperty real e;
+e.setEntryName("e");
+
+@ModellingHypothesis Tridimensional;
+@PhysicalBounds T in [0. : *[;
+
+@StateVariable Stensor alpha_;
+alpha_.setEntryName("MacroscopicViscousStrain");
+@StateVariable Stensor alpha_s_r;
+alpha_s_r.setEntryName("InclusionSolenoidalViscousStrain");
+
+@LocalVariable Stensor4 C0;
+@LocalVariable Stensor4 M0;
+@LocalVariable Stensor4 C_;
+@LocalVariable Stensor4 M_;
+@LocalVariable Stensor4 Cs;
+@LocalVariable Stensor4 Ms;
+@LocalVariable Stensor4 P0v;
+@LocalVariable Stensor4 P0e;
+@LocalVariable Stensor4 D0;
+@LocalVariable Stensor4 I;
+@LocalVariable tfel::math::tmatrix<12,12,real> mat;
+@LocalVariable tfel::math::tmatrix<12,1,real> f_;
+
+@InitLocalVariables{
+I=tfel::math::st2tost2<3u,real>::Id();
+const auto J=tfel::math::st2tost2<3u,real>::J();
+const auto K=tfel::math::st2tost2<3u,real>::K();
+C0=3*ka0*J+2*mu0*K;
+M0=3*kav0*J+2*muv0*K;
+const auto Cr=3*kar*J+2*mur*K;
+const auto Mr=3*kavr*J+2*muvr*K;
+const auto f0=1-fr;
+const auto invdC=invert(Cr-C0);
+const auto invdM=invert(Mr-M0);
+
+using namespace tfel::material::homogenization::elasticity;
+const auto kg0=KGModuli<stress>(ka0,mu0);
+const auto kgv0=KGModuli<stress>(kav0,muv0);
+tfel::math::tvector<3u,real> n={1.,0.,0.};
+P0e=computeAxisymmetricalHillPolarisationTensor(kg0,n,e);
+P0v=computeAxisymmetricalHillPolarisationTensor(kgv0,n,e);
+
+C_=C0+fr*invert(invdC+P0e-fr*P0e);
+M_=M0+fr*invert(invdM+P0v-fr*P0v);
+Cs=fr/f0*invert(invert(C0)-P0e);
+D0=M0-tau0*C0;
+Ms=tau0*Cs+fr/f0*(D0-D0*P0v*D0);
+}
+
+@Integrator {
+//jacobian
+tfel::math::map_derivative<0,0,Stensor,Stensor>(mat) =(M_/dt+C_);
+tfel::math::map_derivative<6,0,Stensor,Stensor>(mat)=((M0+D0*P0v*(M_-M0))/dt+C0);
+tfel::math::map_derivative<0,6,Stensor,Stensor>(mat)=((M_-M0)*(I-P0v*D0)/dt+C_-C0);
+tfel::math::map_derivative<6,6,Stensor,Stensor>(mat)=-((Ms-D0*P0v*(M_-M0)*(I-P0v*D0))/dt+Cs);
+
+tfel::math::map_derivative<0,0,Stensor,real>(f_)=C_*(eto+deto)-C_*(alpha_)-(C_-C0)*(alpha_s_r);
+tfel::math::map_derivative<6,0,Stensor,real>(f_)=-(C0*alpha_-Cs*alpha_s_r-C0*(eto+deto));
+
+TinyMatrixInvert<12,real>::exe(mat);
+tfel::math::tmatrix<12,1,real> dv=mat*f_;
+
+dalpha_=tfel::math::map_derivative<0,0,Stensor,real>(dv);
+dalpha_s_r=tfel::math::map_derivative<6,0,Stensor,real>(dv);
+
+sig+=C_*(deto-dalpha_-dalpha_s_r)+C0*dalpha_s_r;
+}
+
+@TangentOperator{
+if (smt){
+Stensor4 iJ_alpha__alpha_=tfel::math::map_derivative<0,0,Stensor,Stensor>(mat);
+Stensor4 iJ_alpha__alpha_s_r=tfel::math::map_derivative<0,6,Stensor,Stensor>(mat);
+Stensor4 iJ_alpha_s_r_alpha_=tfel::math::map_derivative<6,0,Stensor,Stensor>(mat);
+Stensor4 iJ_alpha_s_r_alpha_s_r=tfel::math::map_derivative<6,6,Stensor,Stensor>(mat);
+
+const auto dalpha_deto=(iJ_alpha__alpha_*C_+iJ_alpha__alpha_s_r*C0);
+const auto dalphar_deto=(iJ_alpha_s_r_alpha_*C_+iJ_alpha_s_r_alpha_s_r*C0);
+Dt=C_-C_*(dalpha_deto+dalphar_deto)+C0*dalphar_deto;
+}
+}

generic-behaviours/homogenization/Idiart_viscoelastic.mtest

@@ -0,0 +1,30 @@
+@ModellingHypothesis 'Tridimensional';
+
+@Behaviour<@interface@> @library@ @behaviour@;
+
+@ExternalStateVariable 'Temperature' {0 : 1000,1:1000};
+@MaterialProperty<constant> 'ka0' 2.e9;
+@MaterialProperty<constant> 'mu0' 1.e9;
+@MaterialProperty<constant> 'kav0' 0.4e9;
+@MaterialProperty<constant> 'muv0' 0.2e9;
+@MaterialProperty<constant> 'kar' 20.e9;
+@MaterialProperty<constant> 'mur' 10.e9;
+@MaterialProperty<constant> 'kavr' 4e6;
+@MaterialProperty<constant> 'muvr' 2e6;
+@MaterialProperty<constant> 'tau0' 0.2;
+@MaterialProperty<constant> 'fr' 0.06557;
+@MaterialProperty<constant> 'e' 10.;
+@ImposedStrain 'EXX' {0 : 0, 6 : 6e-1};
+@ImposedStrain 'EYY' {0 : 0, 6 : 0};
+@ImposedStrain 'EZZ' {0 : 0, 6 : 0};
+@ImposedStrain 'EXY' {0 : 0, 6 : 0};
+@ImposedStrain 'EXZ' {0 : 0, 6 : 0};
+@ImposedStrain 'EYZ' {0 : 0, 6 : 0};
+@Times {0,5 in 200};
+
+@Test<file> @reference_file@ 'EXX' 2 1.e-14;
+@Test<file> @reference_file@ 'EYY' 3 1e-14;
+@Test<file> @reference_file@ 'EZZ' 4 1e-14;
+@Test<file> @reference_file@ 'EXY' 5 1e-14;
+@Test<file> @reference_file@ 'EXZ' 6 1e-14;
+@Test<file> @reference_file@ 'EYZ' 7 1e-14;