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_formulation.mfront

@@ -0,0 +1,142 @@
+@DSL ImplicitII;
+@Behaviour Affine_formulation;
+@Author Martin Antoine;
+@Date 24 / 01 / 26;
+@Description{"Affine formulation for homogenization of a viscoplastic polycrystal, Masson et al. 2001.", based on self-consistent scheme};
+@UseQt false;
+@Algorithm NewtonRaphson_NumericalJacobian;
+@PerturbationValueForNumericalJacobianComputation 1e-10;
+@Epsilon 1e-14;
+
+@TFELLibraries {"Material"};
+@Includes{
+#include "TFEL/Material/IsotropicModuli.hxx"
+#include "TFEL/Material/IsotropicEshelbyTensor.hxx"
+#include "TFEL/Material/AnisotropicEshelbyTensor.hxx"
+#include "TFEL/Material/MicrostructureDescription.hxx"
+#include "TFEL/Material/MicrostructureLinearHomogenization.hxx"
+#include "TFEL/Material/PolyCrystalsSlidingSystems.hxx"}
+
+//! number of phases
+@IntegerConstant Np = 10;
+
+//! normalization strainrate
+@Parameter real gamma0 = 1.0;
+@Parameter real k0 = 1.;
+@Parameter real 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;
+
+@StateVariable Stensor sigma[Np];
+sigma.setEntryName("PhaseReferenceStress");
+
+@AuxiliaryStateVariable real ne;
+
+@LocalVariable Stensor4 I;
+@LocalVariable real frac[Np];
+@LocalVariable real tau0[Np];
+@LocalVariable Stensor4 M[Np];
+@LocalVariable Stensor4 L[Np];
+@LocalVariable tmatrix<6*Np,6,real> A;
+@LocalVariable tmatrix<6*Np,6,real> dsigma_deto;
+
+@InitLocalVariables{
+I=tfel::math::st2tost2<3u,real>::Id();
+for (int r=0;r<Np;r++){
+if (r<5){tau0[r]=tau1;}
+else{tau0[r]=tau2;}
+}
+
+}
+
+@Integrator {
+ne=nexp;
+using PolyCrystalsSlidingSystems =
+	ExtendedPolyCrystalsSlidingSystems<Np, Affine_formulationSlipSystems<real>, real>;
+const auto& gs =
+	PolyCrystalsSlidingSystems::getPolyCrystalsSlidingSystems("polycrystal.csv");
+
+//tangent modulus and polarization/////////////////////////////////
+using namespace tfel::math;
+Stensor4 LSC;
+Stensor e[Np];
+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++){
+  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)/tau0[r], nexp-1);
+  const auto fac= nexp*gamma0/tau0[r]*puisn_1;
+  M[r]+=fac*Nkr;
+  const auto puisn = puisn_1*abs(taukr)/tau0[r];
+  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]);
+}
+
+//Self-consistent scheme///////////////////////
+using namespace tfel::material::homogenization::elasticity;
+const auto IM0=tfel::material::KGModuli<stress>(k0,mu0);
+auto micro=ParticulateMicrostructure<3u,stress>(IM0);
+const auto sphere=Sphere<length>();
+std::vector<Stensor> polarisations;
+const auto zero = Stensor::zero();
+polarisations.resize(Np+1, zero);
+
+for (int r=0 ; r<Np; r++){
+  auto distribution = SphereDistribution<stress>(sphere,frac[r],L[r]);
+  micro.addInclusionPhase(distribution);
+  polarisations[r+1]=-L[r]*e[r];
+}
+const auto hmSC=computeSelfConsistent<3u,stress>(micro,1e-6,true,0,polarisations);
+LSC = hmSC.homogenized_stiffness;
+const auto tauSC = hmSC.effective_polarisation;
+const auto A_SC=hmSC.mean_strain_localisation_tensors;
+
+auto KGSC = computeKGModuli(LSC);
+const auto PSC = computeSphereHillPolarisationTensor<stress>(KGSC);
+const auto MSC_star = invert(invert(PSC)- LSC);
+
+//residues/////////////////////////////////////
+for (int r=0;r<Np;r++){
+  const Stensor4 Ar=A_SC[r+1];
+  map_derivative<Stensor,Stensor>(A,6*r,0)=Ar;
+  fsigma[r] = M[r]*(sigma[r]+dsigma[r])-Ar*(eto+deto);
+  fsigma[r]-=M[r]*(invert(M[r]+MSC_star))*(MSC_star*tauSC - e[r]);
+}
+
+//macroscopic stress//////////////////////////////////////
+sig=LSC*(eto+deto)+tauSC;
+
+}
+
+@TangentOperator{
+if (smt){
+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_formulation.mtest

@@ -0,0 +1,24 @@
+@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' 10.;
+@MaterialProperty<constant> 'tau2' 50.;
+
+@MaterialProperty<function> 'kap' "0.001";
+
+@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 50};
+
+@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-3;
+

generic-behaviours/homogenization/MassonAffineFormulation/Affine_tensors.mfront

@@ -0,0 +1,170 @@
+@DSL ImplicitII;
+@Behaviour Affine_tensors;
+@Author Martin Antoine;
+@Date 26 / 01 / 26;
+@Description{"Affine formulation for homogenization of a viscoplastic polycrystal, (Masson et al. 2001.), based on Green interaction tensors"};
+@UseQt false;
+@Algorithm NewtonRaphson_NumericalJacobian;
+@PerturbationValueForNumericalJacobianComputation 1e-10;
+@Epsilon 1e-14;
+
+@TFELLibraries {"Material"};
+
+@Includes{
+#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 k0 = 1;
+//! 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> DELTA;
+@LocalVariable tmatrix<6*Np,6,real> dsigma_deto;
+
+@InitLocalVariables{
+I=tfel::math::st2tost2<3u,real>::Id();
+J=tfel::math::st2tost2<3u,real>::J();
+K=tfel::math::st2tost2<3u,real>::K();
+DELTA=Delta<real>::get_tensor();
+L0=r0*(3*k0*J+2*mu0*K);
+M0=invert(L0);
+
+for (int r=0;r<Np;r++){
+if (r<5){tau0[r]=tau1;}
+else{tau0[r]=tau2;}
+}
+}
+
+
+
+@Integrator {
+ne=nexp;
+using PolyCrystalsSlidingSystems =
+	ExtendedPolyCrystalsSlidingSystems<Np, Affine_tensorsSlipSystems<real>, real>;
+const auto& gs =
+	PolyCrystalsSlidingSystems::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++){
+  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)/tau0[r], nexp-1);
+  const auto fac= nexp*gamma0/tau0[r]*puisn_1;
+  M[r]+=fac*Nkr;
+  const auto puisn = puisn_1*abs(taukr)/tau0[r];
+  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)=L0;
+  const auto dM = M[r]-M0;
+  for (int s=0;s<Np;s++){
+    if (true){
+    map_derivative<Stensor,Stensor>(MAT,6*r,6*s) +=r0*DELTA[r][s]*dM;
+    map_derivative<Stensor,Stensor>(G,6*r,6*s) -=r0*DELTA[r][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]*Ar[r];
+}
+
+//residues/////////////////////////////////////
+for (int r=0;r<Np;r++){
+  fsigma[r] = M[r]*(sigma[r]+dsigma[r])-M[r]*Ar[r]*(eto+deto);
+  for (int s=0;s<Np;s++){
+    auto Brs = map_derivative<Stensor,Stensor>(B,6*r,6*s);
+    fsigma[r]-=M[r]*Brs*e[s];
+  }
+}
+
+//macroscopic stress//////////////////////////////////////
+auto tau_eff=Stensor::zero();
+for (int r=0;r<Np;r++){
+  tau_eff-=frac[r]*transpose(Ar[r])*e[r];
+}
+
+sig=Chom*(eto+deto)+tau_eff;
+}
+
+
+@TangentOperator{
+if (smt){
+for (int r=0;r<Np;r++){
+  map_derivative<Stensor,Stensor>(A,6*r,0)=M[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,24 @@
+@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' 10.;
+@MaterialProperty<constant> 'tau2' 50.;
+
+@MaterialProperty<function> 'kap' "0.001";
+@MaterialProperty<function> 'r0' "1";
+
+@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 50};
+
+@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,17 @@
+mfront_behaviours_library(MassonAffineFormulation
+  Affine_formulation
+  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_formulation
+  BEHAVIOUR Affine_formulation
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Affine_formulation.ref)
+
+genericmtest(MassonAffineFormulation Affine_tensors
+  BEHAVIOUR Affine_tensors
+  REFERENCE_FILE ${CMAKE_CURRENT_SOURCE_DIR}/references/Affine_tensors.ref)