Subject: [FEniCS-users] Question of Dolfin 0.9.2

From: Jehanzeb Hameed <[email-masked]>
Date: Sun, 4 Oct 2009 15:26:30 -0500

Hello,

I am getting a different answer using ffc form, and ufl forms in
dolfin 0.9.2 . I used form2ufl to convert the .form to .ufl. Both look
fine.

In the forms I am using this element:

For FFC  = FiniteElement("Raviart-Thomas", "tetrahedron", 0)
For UFL  = FiniteElement("Raviart-Thomas", tetrahedron, 0)

On assembling a certain functional,I get different answers in both.
Any idea on why this may be? Also, if  "q" is the Raviart-Thomas
function, then inner(curl(q),curl(q)*dx yields a number for UFL, but
gives nan for FFC. Can someone please tell me what may be going wrong?

-Jehanzeb

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 05 Oct 2009 08:52:24 +0200

> gives nan for FFC. Can someone please tell me what may be going wrong?

Could you be a bit more specific and maybe post the form files?

Kristian

> -Jehanzeb
> _______________________________________________
> FEniCS-users mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/fenics-users
>

---

From: Jehanzeb Hameed <[email-masked]>
Date: Mon, 5 Oct 2009 12:36:34 -0500

I am posting the files below. Here is a bried overview of the files.
The files are A , B and C. Each file has two versions, one in ffc, one
in ufl.
A defines a bilinear form and a linear functional. We get functions
(q,u) from A. The function q belongs to Raviart-Thomas space.
B computes a functional on (q,u). The value of this functional is
widely different for .form and .ufl.
C also computes a  functional on (q,u), namely square of norm of its
curl-component. The .form files gives a negative answer in this case!

The .form files.

File A:

vecElement = FiniteElement("Raviart-Thomas", "tetrahedron", 0)
scalElement = FiniteElement("Lagrange", "tetrahedron", 1)
cellElement = FiniteElement("Discontinuous Lagrange", "tetrahedron", 0)
rhsEle = FiniteElement("Discontinuous Lagrange", "tetrahedron", 1)
rhsVecEle = VectorElement("Discontinuous Lagrange", "tetrahedron", 1)

mixed_element = vecElement + scalElement

(r, v) = TestFunctions(mixed_element)
(q, u) = TrialFunctions(mixed_element)

kappa = Function(cellElement)
epsilon = Function(cellElement)
epsilonMol = Function(cellElement)
uc = Function(rhsEle)
graduc = Function(rhsVecEle)

a =  (dot(q,r) -  epsilon*dot(q,grad(v)) - epsilon*dot(grad(u),r) +
epsilon*epsilon*dot(grad(u),grad(v)) )*dx \
	+ (div(q)*div(r) - kappa*kappa*div(q)*v - kappa*kappa*div(r)*u +
kappa*kappa*kappa*kappa*u*v)*dx

L = ( (epsilon - epsilonMol)*dot(r,graduc) - epsilon*(epsilon -
epsilonMol)*dot(grad(v),graduc) )*dx  \
	+(-kappa*kappa*uc*(-div(r)+kappa*kappa*v) )*dx

File B:

vecElement = FiniteElement("Raviart-Thomas", "tetrahedron", 0)
scalElement = FiniteElement("Lagrange", "tetrahedron", 1)
cellElement = FiniteElement("Discontinuous Lagrange", "tetrahedron", 0)
rhsEle = FiniteElement("Discontinuous Lagrange", "tetrahedron", 1)
rhsVecEle = VectorElement("Discontinuous Lagrange", "tetrahedron", 1)

mixed_element = vecElement + scalElement

(q, u) = Functions(mixed_element)

kappa = Function(cellElement)
epsilon = Function(cellElement)
epsilonMol = Function(cellElement)
uc = Function(rhsEle)
graduc = Function(rhsVecEle)

M =  ( dot(q,q) -  2*epsilon*dot(q,grad(u)) +
epsilon*epsilon*dot(grad(u),grad(u)) \
	+(epsilon-epsilonMol)*(epsilon-epsilonMol)*dot(graduc,graduc) -
2*(epsilon-epsilonMol)*dot(q,graduc)\
	+ 2*epsilon *(epsilon-epsilonMol)*dot(grad(u),graduc) )*dx \
	+ (dot(-div(q) + kappa*kappa*u + kappa*kappa*uc,-div(q) +
kappa*kappa*u + kappa*kappa*uc))*dx

File C:
vecElement = FiniteElement("Raviart-Thomas", "tetrahedron", 0)
scalElement = FiniteElement("Lagrange", "tetrahedron", 1)

mixed_element = vecElement + scalElement
(q, u) = Functions(mixed_element)

M =  dot(curl(q),curl(q))*dx

Now the .ufl files

File A:
vecElement = FiniteElement("Raviart-Thomas", tetrahedron, 0)
scalElement = FiniteElement("Lagrange", tetrahedron, 1)
cellElement = FiniteElement("Discontinuous Lagrange", tetrahedron, 0)
rhsEle = FiniteElement("Discontinuous Lagrange", tetrahedron, 1)
rhsVecEle = VectorElement("Discontinuous Lagrange", tetrahedron, 1)

mixed_element = vecElement + scalElement

(r, v) = TestFunctions(mixed_element)
(q, u) = TrialFunctions(mixed_element)

kappa = Function(cellElement)
epsilon = Function(cellElement)
epsilonMol = Function(cellElement)
uc = Function(rhsEle)
graduc = Function(rhsVecEle)

a =  (inner(q,r) -  epsilon*inner(q,grad(v)) -
epsilon*inner(grad(u),r) + epsilon*epsilon*inner(grad(u),grad(v)) )*dx
\
	+ (div(q)*div(r) - kappa*kappa*div(q)*v - kappa*kappa*div(r)*u +
kappa*kappa*kappa*kappa*u*v)*dx

L = ( (epsilon - epsilonMol)*inner(r,graduc) - epsilon*(epsilon -
epsilonMol)*inner(grad(v),graduc) )*dx  \
	+(-kappa*kappa*uc*(-div(r)+kappa*kappa*v) )*dx

File B:

vecElement = FiniteElement("Raviart-Thomas", tetrahedron, 0)
scalElement = FiniteElement("Lagrange", tetrahedron, 1)
cellElement = FiniteElement("Discontinuous Lagrange", tetrahedron, 0)
rhsEle = FiniteElement("Discontinuous Lagrange", tetrahedron, 1)
rhsVecEle = VectorElement("Discontinuous Lagrange", tetrahedron, 1)

mixed_element = vecElement + scalElement

(q, u) = Functions(mixed_element)

kappa = Function(cellElement)
epsilon = Function(cellElement)
epsilonMol = Function(cellElement)
uc = Function(rhsEle)
graduc = Function(rhsVecEle)

M =  ( inner(q,q) -  2*epsilon*inner(q,grad(u)) +
epsilon*epsilon*inner(grad(u),grad(u)) \
	+(epsilon-epsilonMol)*(epsilon-epsilonMol)*inner(graduc,graduc) -
2*(epsilon-epsilonMol)*inner(q,graduc)\
	+ 2*epsilon *(epsilon-epsilonMol)*inner(grad(u),graduc) )*dx \
	+ (inner(-div(q) + kappa*kappa*u + kappa*kappa*uc,-div(q) +
kappa*kappa*u + kappa*kappa*uc))*dx

File C:

vecElement = FiniteElement("Raviart-Thomas", tetrahedron, 0)
scalElement = FiniteElement("Lagrange", tetrahedron, 1)

mixed_element = vecElement + scalElement
(q, u) = Functions(mixed_element)

M =  inner(curl(q),curl(q))*dx

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 06 Oct 2009 10:01:33 +0200

> The files are A , B and C. Each file has two versions, one in ffc, one

When I asked for the form files, I should have added 'the simplest possible that
will reproduce the error', but that's my bad :)

> widely different for .form and .ufl.

In the previous version of FFC, UFL support was only experimental, and there was
a bug when determining the correct quadrature scheme for forms involving RT and
Nedelec elements. This should have been fixed in the new FFC release (this was
actually a FIAT issue). This might explain the discrepancy?

> curl-component. The .form files gives a negative answer in this case!

Strange, if it is possible can you reduce the complexity of the form files and
resend them? Also include the main.cpp file which you use to compute the
result. Then it is more likely that someone will have a look. I can't promise
that I have time to look at it just now though.

Kristian

---

From: Jehanzeb Hameed <[email-masked]>
Date: Tue, 6 Oct 2009 15:30:54 -0500

> resend them?

I will try and do it sometime. While working on this issues, I have
come across another related issue which I cant resolve. This time its
the norm of curl vector lagrange elements. I have written two
equivalent forms. But they give different answers. The two form files
begin by:

vecElement = VectorElement("Lagrange", "tetrahedron", 1)
scalElement = FiniteElement("Lagrange", "tetrahedron", 1)

then I compute square of curl norm either as

M =  dot(curl(q),curl(q))*dx
OR
M =  (  D(q[1],0)  - D(q[0],1) )*( D(q[1],0)  - D(q[0],1)  )*dx + \
	(  D(q[2],1)  - D(q[1],2) )*(  D(q[2],1)  - D(q[1],2) )*dx + \
	(  D(q[0],2)  - D(q[2],0) )*(  D(q[0],2)  - D(q[2],0) )*dx

However, I get different answers.

Any hints on this?

-Jehanzeb

---

From: Marie Rognes <[email-masked]>
Date: Wed, 07 Oct 2009 08:14:24 +0200

> However, I get different answers.

What do you mean by different answers?

I'll take a look if you send a simple, runnable test case.

--
Marie

---

From: Jehanzeb Hameed <[email-masked]>
Date: Wed, 7 Oct 2009 00:42:11 -0500

> What do you mean by different answers?

The norm of curl is a scalar value. That value is different for both
functionals.  In my original forms, one value was around 180, the
other 200.

> I'll take a look if you send a simple, runnable test case.

I have created a simple runnable case. For the test case the values of
curl norm is 235.379 and  231.563. The case consists of a main
program, 3 form files, and a mesh. Shall I email these to you
directly?

Thanks,
-Jehanzeb

---

From: Marie Rognes <[email-masked]>
Date: Wed, 07 Oct 2009 18:21:44 +0200

> program, 3 form files, and a mesh. Shall I email these to you

Send it to the list and I'll take a look.

--
Marie

---

From: Marie Rognes <[email-masked]>
Date: Fri, 09 Oct 2009 06:16:59 +0200

Jehanzeb Hameed wrote:
> Attaching the files.

I get the same result for the two form files with a  DOLFIN > 0.9.3:

    Curl error is 223.93
    Alt  Curl error is 223.93

In other words, try upgrading :)

For future reference: Since the issue at hand was comparing the two
curl implementations, the test could have been further simplified
by for instance saving the solution u and just sending the solution,
the two ResCurl*.form files and a simplified main.cpp

--
Marie

> On Wed, Oct 7, 2009 at 11:21 AM, Marie Rognes <[email-masked]> wrote:
>   
>> Jehanzeb Hameed wrote:
>>     
>>>>>> a bug when determining the correct quadrature scheme for forms
>>>>>> involving
>>>>>> RT and

---

From: Johan Hake <[email-masked]>
Date: Thu, 8 Oct 2009 21:27:01 +0200

> the two ResCurl*.form files and a simplified main.cpp

While at it, could you take a look at the curl-curl demo too? On my computer 
it produces highly different results for cpp and python. Python takes forever 
and gives weird results.

Johan

> >>>>>>> I am posting the files below. Here is a bried overview of the
> >>>>>>> files. The files are A , B and C. Each file has two versions, one
> >>>>>>> in ffc, one in ufl.
> >>>>>>
> >>>>>>> C also computes a  functional on (q,u), namely square of norm of
> >>>>>>> its curl-component. The .form files gives a negative answer in this
> >>>>>>> case!
> >>>>>>
> >>>>>>>     + (div(q)*div(r) - kappa*kappa*div(q)*v - kappa*kappa*div(r)*u
> >>>>>>> + kappa*kappa*kappa*kappa*u*v)*dx
> >>>>>>>
> >>>>>>> cellElement = FiniteElement("Discontinuous Lagrange", tetrahedron,
> >>>>>>> 0) rhsEle = FiniteElement("Discontinuous Lagrange", tetrahedron, 1)
> >>>>>>> epsilon*inner(grad(u),r) + epsilon*epsilon*inner(grad(u),grad(v))
> >>>>>>> )*dx \
> >>>>>>>     + (div(q)*div(r) - kappa*kappa*div(q)*v - kappa*kappa*div(r)*u
> >>>>>>> + kappa*kappa*kappa*kappa*u*v)*dx
> >>>>>>>
> >>>>>>> cellElement = FiniteElement("Discontinuous Lagrange", tetrahedron,
> >>>>>>> 0) rhsEle = FiniteElement("Discontinuous Lagrange", tetrahedron, 1)
> >>>>>>>     +(epsilon-epsilonMol)*(epsilon-epsilonMol)*inner(graduc,graduc)
> >>>>>>> - 2*(epsilon-epsilonMol)*inner(q,graduc)\
> >>>>>>>>> On assembling a certain functional,I get different answers in
> >>>>>>>>> both. Any idea on why this may be? Also, if  "q" is the
> >>>>>>>>> Raviart-Thomas function, then inner(curl(q),curl(q)*dx yields a
> >>>>>>>>> number for UFL, but gives nan for FFC. Can someone please tell me
> >>>>>>>>> what may be going wrong?
> >>>>>>>>