I would like to know how to solve the following problem. I am new here, so if I do something wrong let me know. I have different Mooney-Rivlin parameters in two different domains and I have to compute the first Piola-Kirchoff tensor. Is it correct to proceed in this way? Or is it wrong because I use the command diff on a discontinuous quantity? Is there a better way to do so? If I try with close values between the domains I do not find problem, otherwise the non linear solver I use later does not converge in 50 iterations.
(Do not focus on the way it is computed the Piola-Kirchoff tensor, my problem is quite complicated and I tried to simplify it, maybe in a wrong way).
Thank you in advance.
def Pk(u, p, mu1, mu2): F = DeformationGradient(u) J = Jacobian(u) I = SecondOrderIdentity(u) C_e = F.T*F Ice, IIce, IIIce = Invariants(C_e) #Strain energy: Mooney-Rivlin psi = (mu1/2)*(Ice - 3) + (mu2/2)*(IIce - 3) gamma1 = diff(psi, Ice) + Ice*diff(psi, IIce) gamma2 = -diff(psi, IIce) P_s = 2*J*F*(gamma1*I + gamma2*C_e)*inv(F).T return variable(P_s - J*p*inv(F).T) DGe = FiniteElement("Discontinuous Lagrange", mesh.ufl_cell(), 0) mu1t = Constant(-3.59e-03) mu2t = Constant(5.52e-03) mu1h = Constant(-8.96e-03) mu2h = Constant(1e-02) mu1= Expression(" pow((x-203),2) + pow((x-287),2) + pow((x-29.5),2) - 25 >= 0 ? mu1h: mu1t", mu1h=mu1h, mu1t=mu1t, element = DGe) mu2= Expression(" pow((x-203),2) + pow((x-287),2) + pow((x-29.5),2) - 25 >= 0 ? mu2h : mu2t", mu2h=mu2h, mu2t=mu2t, element = DGe) P = Pk(u, p, mu1, mu2)