Hi All,
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 P2P1.
Now in the first term of the PDE the coefficients S(p)
and U(p)
are making the problem nonlinear.

Nevertheless, I have defined the function
S(pf)
as below:def S(pf):
return ((pf)^(1/2)+1)^1
However, I am unable to do this and the error is as follows:
UFLException: Cannot take the power of a nonscalar expression <ComponentTensor id=140377318371880>

U(pf)
is not so straight forward. It is a definite integral of the functionS(pf)
in the limitsp
and0
. I am not sure how to handle this. I have tried using sympy as below:def U(pf):
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 U(pf)
function.
NotImplementedError: Cannot take length of nonvector 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.
The functions S(pf)
and 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?
Many Thanks,
Sid