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

@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 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);
  ∂σ₁∕∂Δe₁ = λ₁ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁ ⋅ I₄;
  ∂σ₁∕∂Δe₂ = λ₁₂ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ₁₂ ⋅ I₄;
  ∂σ₂∕∂Δe₁ = ∂σ₁∕∂Δe₂;
  ∂σ₂∕∂Δe₂ = λ₂ ⋅ (I₂⊗I₂) + 2 ⋅ μ₂ ⋅ I₄;
  σ₁ = ∂σ₁∕∂Δe₁ ⋅ e₁ + ∂σ₁∕∂Δe₂ ⋅ e₂;
  σ₂ = ∂σ₂∕∂Δe₁ ⋅ e₁ + ∂σ₂∕∂Δe₂ ⋅ e₂;
  I = kappa*V;
  ∂I∕∂ΔV = kappa ⋅ I₄;
}