Subject: [FEniCS-users] trying to find gradients

From: phil marinier <[email-masked]>
Date: Mon, 17 Aug 2009 23:33:50 -0400

I am trying to find the gradients of a function over my domain. As suggested I used this .ufl file:

element = VectorElement("Lagrange", triangle, 1)
 
 v = TestFunction(element)
 u = TrialFunction(element)
 f = Function(element)
 
 a = (v, du)
 L = (v, grad(f))
I compiled this using ffc and ran the variational problem with no boundary condition.
When I did:

GradX = Gradients[0];
GradX.vector().size();

I noticed that the solution has exactly the same amount of values as I have nodes in my domain. I assumed this meant that the gradients were deinfed at the nodes. I thought that the gradients should have been defined over the elements, not at the nodes, so I decided to double check the results using a simple test case from a previous grad student working on a similar problem. The results did not compare very well at all. They were mostly in the same general ball park, but i was averaging maybe 20 or 30% difference in the areas with steep gradients. 

The way I  did this comparision was to assume the values I was getting were nodal values, and then averaging these values out over the elements. This way, if the dolfin values really were nodal values, the results should have been close.

When I compared the values the other grad student had gotten at the elements with the values I got using dolfin, they were very close, definately within acceptable error, which leads me to the conclusion that the values I am getting from dolfin must the the elemental, as opposed to the nodal, gradients. But if that is the case, I should be getting a size that is equal to my number of elements, or a whole number multiple thereof.

When I extract the vector from the gradient in a specific driection, am I getting nodal values or elemental values? How do you explain the problem associated with each scenario? I need to end up with the nodal values, but i don't mind writing code to average it out. How do you suggest I do this?

Thank you for your help

Phil

---

From: Anders Logg <[email-masked]>
Date: Tue, 18 Aug 2009 10:10:57 +0200

On Mon, Aug 17, 2009 at 11:33:50PM -0400, phil marinier wrote:
> I noticed that the solution has exactly the same amount of values as I have
> nodes in my domain. I assumed this meant that the gradients were deinfed at the
> nodes. I thought that the gradients should have been defined over the elements,
> not at the nodes, so I decided to double check the results using a simple test
> case from a previous grad student working on a similar problem. The results did
> not compare very well at all. They were mostly in the same general ball park,
> but i was averaging maybe 20 or 30% difference in the areas with steep
> gradients.

Yes, the above projects the gradient to piecewise linears, so you get
2 values per vertex in 2D. If you want element values, then just
project to piecewise constants:

P1 = FiniteElement("CG", triangle, 1)
P0 = VectorElement("DG", triangle, 0)

v = TestFunction(P0)
u = TrialFunction(P0)
f = Function(P1)

a = (v, du)
L = (v, grad(f))

--
Anders

> The way I  did this comparision was to assume the values I was getting were
> nodal values, and then averaging these values out over the elements. This way,
> if the dolfin values really were nodal values, the results should have been
> close.
>
> When I compared the values the other grad student had gotten at the elements
> with the values I got using dolfin, they were very close, definately within
> acceptable error, which leads me to the conclusion that the values I am getting
> from dolfin must the the elemental, as opposed to the nodal, gradients. But if
> that is the case, I should be getting a size that is equal to my number of
> elements, or a whole number multiple thereof.
>
> When I extract the vector from the gradient in a specific driection, am I
> getting nodal values or elemental values? How do you explain the problem
> associated with each scenario? I need to end up with the nodal values, but i
> don't mind writing code to average it out. How do you suggest I do this?
>
> Thank you for your help
>
> Phil
>
>

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

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 18 Aug 2009 12:27:30 +0200

Quoting Anders Logg <[email-masked]>:

> > I am trying to find the gradients of a function over my domain. As
> suggested I
> > used this .ufl file:
> >
> >

Or you can get the gradients at the quadrature points (if you have 2nd order
displacements).

P1 = FiniteElement("CG", triangle, 2)
P0 = VectorElement("Quadrature", triangle, 2)

v = TestFunction(P0)
u = TrialFunction(P0)
f = Function(P1)

a = (v, du)
L = (v, grad(f))

Where the '2' in the quadrature element means that it will use enough points to
integrate af 2nd order polynomial (a 2x2 gauss scheme).
Note that changing '2' -> '3' will result in the same 2x2 scheme, you need to
change '2' -> '4' to get a 3x3 scheme.

Kristian