Hello everyone, i am new using fenics, i have development a short script to model a hertzian contact and stresses in a cad geometry converted to xml file using dolfin, my question is what am i doing wrong in my code because Fenics says the next error: “o Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Found different Arguments with same number and part.
Did you combine test or trial functions from different spaces?”
i provide my code here, i appreciate your helping:
from fenics import *
# Crear malla a partir de un archivo XML
# Definir la malla y el espacio de funciones
filename = '/home/jalastrag/fenics_proj/mesh.xml'
mesh = Mesh(filename)
# Definir los espacios de funciones para la deformación (u) y la presión (p)
V = VectorFunctionSpace(mesh, "P", 2)
Q = FunctionSpace(mesh, "P", 1)
# Definir las condiciones de contorno (aquí se asumen condiciones de Dirichlet homogéneas)
bc_u = DirichletBC(V, Constant((0, 0, 0)), "on_boundary")
bc_p = DirichletBC(Q, Constant(0), "on_boundary")
# Definir las funciones de prueba y de prueba de incremento
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Definir los parámetros elásticos y de Poisson
E = Constant(1e6) # Módulo de Young
nu = Constant(0.3) # Coeficiente de Poisson
# Definir la ecuación elástica lineal
sigma = 2*nu*E*sym(grad(u)) + E*(1-2*nu)*tr(sym(grad(u)))*Identity(len(u))
F_elastic = inner(sigma, grad(v))*dx
# Definir la ecuación de la presión (puede variar dependiendo de la formulación)
F_pressure = -p*q*dx
# Definir la ecuación de equilibrio variacional
F = F_elastic + F_pressure
# Resolver el problema variacional
u_solution = Function(V)
p_solution = Function(Q)
solve(F == 0, u_solution, bc_u)
solve(F_pressure == 0, p_solution, bc_p)
# Guardar resultados en archivos VTK
File("deformacion.pvd") << u_solution
File("presion.pvd") << p_solution
First, since your are solving for u and p, is better to use a MixedElement approach to avoid the error
"Found different Arguments with same number and part.
Did you combine test or trial functions from different spaces?"
Second, the error Differentiating residual from F to obtain Jacobian J = F’. comes from the fact that your are trying to use the automatic non-linear solver but using TrialFunctions.
The modified code below works and converges in 1 iteration. Note that I have used sol as a function of the mixed space V, followed by u, p = split(sol), which is typical when solving Navier-Stokes equations in Legacy Fenics.
from dolfin import *
mesh = UnitCubeMesh(5, 5, 5)
# Definir los espacios de funciones para la deformación (u) y la presión (p)
P2 = VectorElement('CG', mesh.ufl_cell(), 2)
P1 = FiniteElement('CG', mesh.ufl_cell(), 1)
element = MixedElement([P2, P1]) # Con multiplicador
V = FunctionSpace(mesh, element)
# Definir las condiciones de contorno (aquí se asumen condiciones de Dirichlet homogéneas)
bc_u = DirichletBC(V.sub(0), Constant((0, 0, 0)), "on_boundary")
bc_p = DirichletBC(V.sub(1), Constant(0), "on_boundary")
# Definir las funciones de prueba y de prueba de incremento
sol = Function(V)
u, p = split(sol)
v, q = TestFunctions(V)
# Definir los parámetros elásticos y de Poisson
E = Constant(1e6) # Módulo de Young
nu = Constant(0.3) # Coeficiente de Poisson
# Definir la ecuación elástica lineal
sigma = 2*nu*E*sym(grad(u)) + E*(1-2*nu)*tr(sym(grad(u)))*Identity(len(u))
F_elastic = inner(sigma, grad(v))*dx
# Definir la ecuación de la presión (puede variar dependiendo de la formulación)
F_pressure = -p*q*dx
# Definir la ecuación de equilibrio variacional
F = F_elastic + F_pressure
# Resolver el problema variacional
# sol = Function(V)
solve(F == 0, sol, bc_u)
uh, ph = sol.split(True)
# Guardar resultados en archivos VTK
File("deformacion.pvd") << uh
File("presion.pvd") << uh
Hello thank you for your answer sir, i apologize for don’t share my mesh because i was not be able to upload, soon i hope do it, how ever i appreciate your answer, this routine that i showed my intention is to make a nonlinear problem about wire helix spring and study behavior of them in hertzian contacts, do you have any docummentation about it where do i can find simulation using fenics of that kind of problem?.. again thank you very much for your answer
Hello Jesus thank oyu very much for your answer i used your script how ever i can’t see in paraview any answer (cube response in deformation or stresses) could you explain me please where can i change number or boundary condition to see any result? how i wrote before my intettion (i am a new user fenics) is to make a simulatin of wire helix spring studying their hertzian contact, if you have some documentation i will appreciate it, thank you for share this code, and a nother question this code repsonse means that my results are zero? : olving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.
Your system is not doing any work. You can add a body force:
fbody = Expression(('sin(x[1])', 'cos(x[0])', 'x[0]*x[1]'), degree=2)
F_body = inner(fbody, v)*dx
# Definir la ecuación de equilibrio variacional
F = F_elastic + F_pressure - F_body
On the other hand, since you are new to FEniCS, I recommend you to use the new DOLFINx, which has a really nice tutorial (or this) and a very rich documentation. There are some results here about hertzian contact in DOLFINx, and here for the Legacy FEniCS.