@DSL DefaultGenericBehaviour;
@Behaviour MultiphaseModel;
@ModellingHypotheses {PlaneStrain, Tridimensional};

@Gradient StrainStensor e₁;
e₁.setEntryName("MatrixStrain");
@Flux StressStensor σ₁;
σ₁.setEntryName("MatrixStress");

@Gradient StrainStensor e₂;
e₂.setEntryName("InclusionStrain");
@Flux StressStensor σ₂;
σ₂.setEntryName("InclusionStress");

@Gradient StrainStensor V;
V.setEntryName("RelativeDisplacement");
@Flux StressStensor I;
I.setEntryName("InteractionForce");

@TangentOperatorBlock ∂σ₁∕∂Δe₁;
@AdditionalTangentOperatorBlock ∂σ₁∕∂Δe₂;
@AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₁;
@AdditionalTangentOperatorBlock ∂σ₂∕∂Δe₂;
@AdditionalTangentOperatorBlock ∂I∕∂ΔV;

@MaterialProperty stress Y1;
Y1.setEntryName("MatrixYoungModulus");
@MaterialProperty real nu1;
nu1.setEntryName("MatrixPoissonRatio");
@MaterialProperty stress Y2;
Y2.setEntryName("FiberYoungModulus");
@MaterialProperty real nu2;
nu2.setEntryName("FiberPoissonRatio");
@MaterialProperty real rho;
rho.setEntryName("FiberVolumeFraction");
@MaterialProperty real s;
s.setEntryName("Size");

@ProvidesTangentOperator;
@Integrator {
  // remove useless warnings, as we always compute the tangent operator
  static_cast<void>(computeTangentOperator_);
const Stensor4 Cₕ = {
      1, 2, 3, 4, 5, 6, // 1ère ligne
      5, 6, 7, 8, 9, 10, // 2ème ligne
      1, 2, 3, 4, 11, 12, // 3ème ligne
      1, 2, 3, 4, 13, 14, // 4ème ligne
      1, 2, 3, 4, 13, 14, // 4ème ligne
      1, 2, 3, 4, 13, 14 // 4ème ligne
  };
  const auto λ₁ = computeLambda(Y1,nu1);
  const auto μ₁ = computeMu(Y1,nu1);
  const auto λ₂ = rho*computeLambda(Y2,nu2);
  const auto μ₂ = rho*computeMu(Y2,nu2);
  constexpr const auto λ₁₂ = 0.;
  constexpr const auto μ₁₂ = 0.;
  const auto kappa = μ₁/12/(1-rho)/pow(s,2);
  const auto C1 = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄;
  const auto C2 = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄;
  const auto idC = invert(C2-C1);
  const auto avgC = rho*C1+(1-rho)*C2;
  const auto D11 = (1-rho)*C1-(C1*idC*(Cₕ-avgC)*idC*C1);
  ∂σ₁∕∂Δe₁ = C1;
  ∂σ₁∕∂Δe₂ = λ₁₂ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁₂ ⋅ I₄;
  ∂σ₂∕∂Δe₁ = ∂σ₁∕∂Δe₂;
  ∂σ₂∕∂Δe₂ = C2;
  σ₁ = ∂σ₁∕∂Δe₁ ⋅ e₁ + ∂σ₁∕∂Δe₂ ⋅ e₂;
  σ₂ = ∂σ₂∕∂Δe₁ ⋅ e₁ + ∂σ₂∕∂Δe₂ ⋅ e₂;
  I = kappa*V;
  ∂I∕∂ΔV = kappa ⋅ I₄;
}