Dear all,
I want to multiply a function by the same constant at each node and then integrate. I tried to project a constant on a scalar functionspace and use this, but the error message I obtain is “ArityMismatch: Adding expressions with non-matching form arguments?”. A slightly different equation without the scalar functionspace is solved without a problem. Can you help me figure out how to solve the arity mismatch? Thank you so much in advance.
from dolfin import *
#Mesh and functionspaces
lambda0 = 5
mesh = IntervalMesh(10,0,1)
S = FunctionSpace(mesh, "Lagrange", 1)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
MIXED = Q * Q
W = FunctionSpace(mesh, MIXED)
#Boundary conditions
x0 = CompiledSubDomain("near(x[0], side,1e-12) && on_boundary" , side = 0)
x1 = CompiledSubDomain("near(x[0], side,1e-12) && on_boundary" , side = 1)
bcbx1 = DirichletBC(W.sub(1), Expression(("mu0"),mu0=1e5, degree=1), x1)
bcbx = DirichletBC(W.sub(0), Expression(("0.0"), degree=2), x0)
bcs = [bcbx, bcbx1]
#Functions
(du, dmu) = TrialFunctions(W)
z = Function(W)
dz = TrialFunction(W)
u, mu = split(z)
(v, q) = TestFunctions(W)
#Definitions
I = Identity(1)
F = 1 + u.dx(0)
J = det(F)*lambda0**2
lambdaspace = project(lambda0, S) #This defines a constant on the whole mesh
#Equations
flux = (1-J)*F**(-1)*F**(-1)*mu.dx(0)
sigma = F-F**(-1)+((ln(1-1/J)+1/J+1/J**2)+1e5-mu)*J*F**(-1)
sigma23 = lambdaspace-1./lambdaspace+((ln(1-1/J)+1/J+1/J**2)+1e5-mu)*J*1./lambdaspace
#############################################Works##############################################
F = inner(sigma, v.dx(0))*dx+2*inner(v.dx(0),sigma23)*dx+q.dx(0)*flux*dx
#############################################Works##############################################
#############################################Does not work######################################
F = inner(sigma, v.dx(0))*dx+2*inner(lambdaspace,sigma23)*dx+q.dx(0)*flux*dx
#############################################Does not work######################################
#Solver
J = derivative(F, z, dz)
pde_system = NonlinearVariationalProblem(F, z, bcs, J=J)
solver = NonlinearVariationalSolver(pde_system)
solver.solve()
#Plot solution
F = 1 + u.dx(0)
J = det(F)*lambda0**2
plot(J)
It’s not immediately clear to me what the intended behavior is, but the root of the problem is that a nonlinear residual needs to be a linear functional of a test function. So, in your example, the first definition of F is linear in the test function (v,q), while the second definition of F is not, because of the second term, which does not contain any test function.
Terms that are linear in the test function have arity 1, terms that are bilinear in a test and trial function have arity 2, and constant terms have arity 0. In a nonlinear residual, all terms need arity 1. For example:
from dolfin import *
mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
uh = Function(V)
from ufl.algorithms.formtransformations import compute_form_arities
print(compute_form_arities(u*v*dx)) # LHS of a linear problem
print(compute_form_arities(v*dx)) # RHS of a linear problem
print(compute_form_arities(uh*v*dx)) # Residual of a nonlinear problem
# Residual of a linear problem:
linearRes = u*v*dx + v*dx
print(compute_form_arities(linearRes)) # Valid reason to have multiple arities
print(compute_form_arities(lhs(linearRes)))
print(compute_form_arities(rhs(linearRes)))
print(compute_form_arities(uh*dx)) # Integrates to a known constant
1 Like
Thank you very much for your reply. I understand all terms must be linear in the test function.The equation I’m solving is not F == 0, but F == C, where C is a constant. Is this at all possible?
Again, the interpretation of the problem is unclear to me, but, if you mean something like: Find u such that for all v\in\mathcal{V}
F(v,u) = C\neq 0\text{ ,}
this is not possible for F linear in its first argument and \mathcal{V} a linear space, since, if there existed such a u, it would also need to also satisfy
F(2v,u) = C\text{ ,}
(because v\in\mathcal{V}\Rightarrow 2v\in\mathcal{V}) but, by linearity of F in its first argument, F(v,u) = C\Rightarrow F(2v,u) = 2F(v,u) = 2C\neq C.
1 Like
Thank you very much for your reply. I’m solving a 3D problem that has a given displacement in 2 directions, therefore I only need the displacement in 1 direction. I tried to solve this with one dimension to save computational power. Solving Int (S*dv/dx) dx == 0 by hand, I obtain a scalar equation with terms linear in v and u as well as constant terms from the already known components of the stress tensor.
However, since it appears that it does not work, I will solve it in 3D with the appropriate boundary conditions instead.