I am working on PDE of the below weak form: (not the exact problem I am dealing with, but I have simplified it without worrying about its physical sense)
a = (p*S(p)-U(p))*inner(epsilon(v),Identity(ndim))*dx + q*tr(epsilon(u))*dx
u,p are trial functions
v,q are test functions
I am handling this problem using a mixed finite element. (
u is vector and
p is scalar solutions) Correspondingly the mixed element is P2-P1.
Now in the first term of the PDE the coefficients
U(p) are making the problem nonlinear.
Nevertheless, I have defined the function
However, I am unable to do this and the error is as follows:
UFLException: Cannot take the power of a non-scalar expression <ComponentTensor id=140377318371880>
U(pf)is not so straight forward. It is a definite integral of the function
S(pf)in the limits
0. I am not sure how to handle this. I have tried using sympy as below:
z = sympy.Symbol("z")
return sympy.integrate(S(z),(z,0,pf)) - pf*S(pf)
With this I get an error as below when evaluating the return of the
NotImplementedError: Cannot take length of non-vector expression.
I have tried another way using the assemble function in FEniCS as below:
def U(pf): mesh = IntervalMesh(10,pf,0.0) z = SpatialCoordinate(mesh) return assemble(S(z)*dz(degree=2)) - pf*S(pf)
With the second method I get an error on the evaluation of the IntervalMesh definition as below:
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported: 1. dolfin.cpp.generation.IntervalMesh(arg0: int, arg1: float, arg2: float) 2. dolfin.cpp.generation.IntervalMesh(arg0: MPICommWrapper, arg1: int, arg2: float, arg3: float) Invoked with: 10, Indexed(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 3), MixedElement(VectorElement(FiniteElement('Lagrange', triangle, 2), dim=2), FiniteElement('Lagrange', triangle, 1), FiniteElement('Lagrange', triangle, 1), FiniteElement('Lagrange', triangle, 1))), 1, None), MultiIndex((FixedIndex(2),))), 0.0
I believe these errors are related to the solution
p being passed directly for function evaluation.
U(pf) need to be evaluated at every node i.e. for every value of
p in the domain and should be returned to the weak form.
I believe that I am making a basic mistake and cannot get my head around it. Can anyone of you help with this?