← Back to team overview

ffc team mailing list archive

try out P2 triangle element

 

ok, so I almost have a demo that demonstrates how to roughly implement curved triangles. I checked the code against my own matlab/c++ code for the case when the triangles are straight. I tried this for a mesh of two triangles of diameter approx 0.5. The max error in the matrix entries for this straight case is 1.77E-7. I would have thought it would be smaller, but maybe my own code is not the greatest.

Anyway, I also compared the iso-param FFC implementation to my own stuff when one triangle is curved. The error was 1.11. So something is wrong. I think the problem is this. I cannot simply change

Jinv_00 ---> NEW_Jinv_00[ip]

In the formula that FFC generates. Because that quad formula is dependent on the jacobian being constant. I think. Is this correct?

- Shawn

On Mon, 22 Jun 2009, Anders Logg wrote:

I think Shawn just wants to get the quadrature points so he can
hand-modify the code for tabulate_tensor to handle higher-order mapped
elements (to make a prototype for higher-order elements).

Correct me if I'm wrong.

In that case, we won't need to make them available since they are
available to the form compiler when the code is being generated, but
it would help if they were available in a comment in the generated code.

--
Anders


On Mon, Jun 22, 2009 at 10:28:29PM +0200, Martin Sandve Alnæs wrote:
On Mon, Jun 22, 2009 at 9:10 PM, Kristian
Oelgaard<k.b.oelgaard@xxxxxxxxxx> wrote:
Quoting Shawn Walker <walker@xxxxxxxxxxxxxxx>:

Actually, I am still a little confused.  Whenever you compute an integral
by quadrature, you need to evaluate the integrand at the quadrature
points.  Those are the points I need.

If you have a quadrature element in your form

element = FiniteElement("Quadrature", triangle, 2)

you can call the
tabulate_coordinates() function of the dof_map associated with the quadrature
element. If you don't have a quadrature element there is no direct way of
getting the quadrature points.

I just need the quad rule that is used in a particular case.  I see the
weights are defined in the tabulate_tensor routine.  Are the points listed
in a comment somewhere?

No, currently not. But it might be a natural thing to have. I'll look into it.

It could be an idea to make it available through UFC,
but then again each ufc::*_integral may use different
quadrature rules for different terms...
However, an exception can be thrown if there is no unique rule in use.

Martin
_______________________________________________
FFC-dev mailing list
FFC-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/ffc-dev

References