71 add massonaffineformulation behaviours

README.md

@@ -19,6 +19,7 @@ It is mainly decomposed in two parts:
     - `plasticity`
     - `viscoelasticity`
     - `viscoplasticity`
+    - `homogenization`
 - `materials`: this part gathers to specific materials.
 
 The last purpose of this project is to show how to build a compilation

generic-behaviours/homogenization/CMakeLists.txt

@@ -25,3 +25,5 @@ genericmtest(Homogenization BetaRule
 genericmtest(Homogenization Idiart_viscoelastic
   BEHAVIOUR Idiart_viscoelastic
   REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Idiart_viscoelastic.ref)
+
+add_subdirectory(MassonAffineFormulation)

generic-behaviours/homogenization/MassonAffineFormulation/Affine_tensors.mfront

@@ -0,0 +1,169 @@
+@DSL ImplicitII;
+@Behaviour Affine_tensors;
+@Author Martin Antoine;
+@Date 06 / 03 / 26;
+@Description{"Affine formulation for homogenization of a viscoplastic polycrystal, (Masson et al. 2001.), based on Green interaction tensors"};
+@UseQt false;
+@Algorithm NewtonRaphson;
+@Epsilon 1e-14;
+
+@TFELLibraries {"Material"};
+
+@Includes{
+#include "TFEL/Material/IsotropicModuli.hxx"
+#include "extra-headers/TFEL/Material/tensors.hxx"
+#include "TFEL/Material/PolyCrystalsSlidingSystems.hxx"}
+
+//! number of phases
+@IntegerConstant Np =10;
+
+//! normalization strainrate
+@Parameter real gamma0 = 1.0;
+//! bulk modulus of the reference medium used for the interaction tensors
+@Parameter stress mu0 = 1;
+
+@ModellingHypothesis Tridimensional;
+@OrthotropicBehaviour;
+@CrystalStructure HCP;
+@SlidingSystems{<1, 1, -2, 0>{1, -1, 0, 0},
+                <-2, 1, 1, 3>{1, -1, 0, 1},
+                <-2, 1, 1, 0>{0, 0, 0, 1},
+                <1, 1, -2, 0>{1, -1, 0, 1}};
+
+@MaterialProperty real nexp;
+@MaterialProperty real tau1;
+@MaterialProperty real tau2;
+@MaterialProperty real kap;
+//@MaterialProperty real r0;
+
+@StateVariable Stensor sigma[Np];
+sigma.setEntryName("PhaseReferenceStress");
+
+@AuxiliaryStateVariable real ne;
+
+@LocalVariable Stensor4 I;
+@LocalVariable Stensor4 J;
+@LocalVariable Stensor4 K;
+@LocalVariable Stensor4 L0;
+@LocalVariable Stensor4 M0;
+@LocalVariable real frac[Np];
+@LocalVariable real tau0[Np];
+@LocalVariable Stensor4 M[Np];
+@LocalVariable Stensor4 L[Np];
+@LocalVariable Stensor4 Ar[Np];
+@LocalVariable tmatrix<6*Np,6,real> A;
+@LocalVariable std::array<std::array<tfel::math::st2tost2<3u,real>,Np>,Np> GAMMA;
+@LocalVariable tmatrix<6*Np,6,real> dsigma_deto;
+@LocalVariable real r0;
+@LocalVariable real muhom;
+
+@InitLocalVariables{
+I=tfel::math::st2tost2<3u,real>::Id();
+J=tfel::math::st2tost2<3u,real>::J();
+K=tfel::math::st2tost2<3u,real>::K();
+GAMMA=Gamma<real>::get_tensor();
+r0=1;
+L0=2*mu0*K;
+M0=1/(2*mu0)*K;
+r0=max(real(1),r0);
+}
+
+
+
+@Integrator {
+ne=nexp;
+using ExtendedPolyCrystalsSlidingSystems =
+	ExtendedPolyCrystalsSlidingSystems<Np, Affine_tensorsSlipSystems<real>, real>;
+const auto& gs =
+	ExtendedPolyCrystalsSlidingSystems::getPolyCrystalsSlidingSystems("polycrystal.csv");
+
+//tangent modulus and polarization/////////////////////////////////
+using namespace tfel::math;
+Stensor e[Np];
+Stensor4 Chom;
+Stensor dpsi_dsigma[Np];
+
+for (int r=0;r<Np;r++){
+ M[r]=Stensor4::zero();
+ dpsi_dsigma[r]=Stensor::zero();
+ frac[r]=gs.volume_fractions[r];
+ for (int k=0;k<Nss;k++){
+ auto tau0rk = tau1;
+ if (k>int(Nss/12)){tau0rk=tau1;}
+ else{tau0rk=tau2;}
+  const auto Nkr = gs.mus[r][k]^gs.mus[r][k];
+  const auto taukr = gs.mus[r][k] | (sigma[r]+dsigma[r]);
+  const auto puisn_1 = pow(abs(taukr)/tau0rk, nexp-1);
+  const auto fac= nexp*gamma0/tau0rk*puisn_1;
+  M[r]+=fac*Nkr;
+  const auto puisn = puisn_1*abs(taukr)/tau0rk;
+  const auto sgn= taukr< 0 ? -real(1) : real(1);
+  dpsi_dsigma[r]+=sgn*gamma0*puisn*gs.mus[r][k]; 
+ }
+ e[r]=dpsi_dsigma[r]-M[r]*(sigma[r]+dsigma[r]);
+ M[r]=M[r]+kap*I;
+ L[r]=invert(M[r]);
+}
+
+//Operators A and B/////////////////////////////////
+tmatrix<6*Np,6*Np,real> MAT;
+tmatrix<6*Np,6*Np,real> G;
+tmatrix<6*Np,6,real> E;
+for (int r=0;r<Np;r++){
+  map_derivative<Stensor,Stensor>(E,6*r,0)=I;
+  for (int s=0;s<Np;s++){
+   const auto dL = L[s]-r0*L0;
+    map_derivative<Stensor,Stensor>(MAT,6*r,6*s) =1/r0*GAMMA[r][s]*dL;
+    map_derivative<Stensor,Stensor>(G,6*r,6*s) =1/r0*GAMMA[r][s]*L[s];
+  }
+  map_derivative<Stensor,Stensor>(MAT,6*r,6*r) +=I;
+}
+
+TinyMatrixInvert<6*Np,real>::exe(MAT);
+A = MAT*E;
+tmatrix<6*Np,6*Np,real> B = MAT*G;
+
+Chom=Stensor4::zero();
+for (int r=0;r<Np;r++){
+  Ar[r]=map_derivative<Stensor,Stensor>(A,6*r,0);
+  Chom+=frac[r]*L[r]*Ar[r];
+}
+
+
+auto KGhom=tfel::material::computeKGModuli(Chom);
+muhom=KGhom.mu;
+r0=muhom/mu0;
+
+auto tau_eff=Stensor::zero();
+//residues/////////////////////////////////////
+for (int r=0;r<Np;r++){
+  tau_eff-=frac[r]*L[r]*e[r];
+  fsigma[r] = sigma[r]+dsigma[r]-L[r]*Ar[r]*(eto+deto)+L[r]*e[r];
+  for (int s=0;s<Np;s++){
+    auto Brs = map_derivative<Stensor,Stensor>(B,6*r,6*s);
+    tau_eff+=frac[r]*L[r]*Brs*e[s];
+    fsigma[r]-=L[r]*Brs*e[s];
+  }
+  dfsigma_ddsigma(r,r)=I;
+}
+
+//macroscopic stress//////////////////////////////////////
+sig=Chom*(eto+deto)+tau_eff;
+
+}
+
+
+@TangentOperator{
+for (int r=0;r<Np;r++){
+map_derivative<Stensor,Stensor>(A,6*r,0)=L[r]*Ar[r];
+}
+tmatrix<6*Np,6*Np,real> iJ = jacobian;
+TinyMatrixInvert<6*Np,real>::exe(iJ);
+dsigma_deto=iJ*A;
+Dt=Stensor4::zero();
+for (int r=0;r<Np;r++){
+   const auto dsigmar_deto = map_derivative<Stensor,Stensor>(dsigma_deto,6*r,0);
+   Dt+=frac[r]*dsigmar_deto;
+}
+
+}

generic-behaviours/homogenization/MassonAffineFormulation/Affine_tensors.mtest

@@ -0,0 +1,23 @@
+@ModellingHypothesis 'Tridimensional';
+@Behaviour<@interface@> @library@ @behaviour@;
+@ExternalStateVariable 'Temperature' {0 : 1000,1:1000};
+
+@MaterialProperty<function> 'nexp' "1+0.9*(exp(2.4*t)-1)";
+@MaterialProperty<constant> 'tau1' 5.;
+@MaterialProperty<constant> 'tau2' 1;
+
+@MaterialProperty<function> 'kap' "0.05";
+
+@ImposedStrain 'EXX' {0 : 1, 1 : 1};
+@ImposedStrain 'EYY' {0 : -1, 1 : -1};
+@ImposedStrain 'EZZ' {0 : 0, 1 : 0};
+@ImposedStrain 'EXY' {0 : 0, 1 : 0};
+@ImposedStrain 'EXZ' {0 : 0, 1 : 0};
+@ImposedStrain 'EYZ' {0 : 0, 1 : 0};
+@Times {0.,1 in 100};
+
+@OutputFilePrecision 14;
+
+@Test<file> @reference_file@ 'SXX' 8 1.e-3;
+@Test<file> @reference_file@ 'SYY' 9 1.e-3;
+@Test<file> @reference_file@ 'ne' 74 1e-5;

generic-behaviours/homogenization/MassonAffineFormulation/CMakeLists.txt

@@ -0,0 +1,12 @@
+mfront_behaviours_library(MassonAffineFormulation
+  Affine_tensors
+  )
+
+STRING(COMPARE EQUAL "${CMAKE_SOURCE_DIR}" "${CMAKE_BINARY_DIR}" _insource)
+if(NOT _insource)
+  configure_file(polycrystal.csv polycrystal.csv)
+endif()
+
+genericmtest(MassonAffineFormulation Affine_tensors
+  BEHAVIOUR Affine_tensors
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Affine_tensors.ref)