Subject: [FFC-dev] get quad points

From: Shawn Walker <[email-masked]>
Date: Sat, 20 Jun 2009 18:14:12 -0400 (EDT)

This is probably a dumb question, but is it easy to get the actual 
quadrature points used in a compiled "tabulate_tensor" routine?  I need 
this to try some manual modification of the local FEM matrix.

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Sun, 21 Jun 2009 10:57:31 +0200

> this to try some manual modification of the local FEM matrix.

Not sure I understand exactly what you want, but if you want the coordinates of
the quadrature points, you can use e.g., FiniteElemet("Quadrature", 1) and
tabulate the coordinates of the dof_map.

Kristian

> - Shawn
> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: Shawn Walker <[email-masked]>
Date: Sun, 21 Jun 2009 05:34:35 -0400 (EDT)

That is what I want.

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 13:57:59 -0400 (EDT)

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.

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?

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 22 Jun 2009 21:10:15 +0200

> 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.

> weights are defined in the tabulate_tensor routine.  Are the points listed

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

Kristian

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 22 Jun 2009 22:28:29 +0200

> 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

---

From: Anders Logg <[email-masked]>
Date: Mon, 22 Jun 2009 22:41:26 +0200

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

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 16:50:40 -0400 (EDT)

> elements (to make a prototype for higher-order elements).

Yes!

> it would help if they were available in a comment in the generated code.

I agree that having them in a comment would be good.

- Shawn

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 22 Jun 2009 22:52:34 +0200

Ok. Temporarily, maybe you'll be happy with the code
generated by SFC (see the demo/*/cpp/ directories in SyFi).
There both weights and points are embedded directly in the code.

Martin

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 17:03:16 -0400 (EDT)

Well, do you know which example uses a 4-pt quad rule?

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 17:08:00 -0400 (EDT)

Actually, would it be possible for you to just give me the 4-pt quad rule 
on a triangle.  The problem is, I don't have GiNAC setup, which is needed 
to run SyFi, right?

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 17:09:26 -0400 (EDT)

- Shawn

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Mon, 22 Jun 2009 23:21:35 +0200

SFC doesn't use the same quadrature rules as
FFC, so that makes no sense in my head.
Either install SyFi properly with dependencies,
or keep using FFC and ignore my suggestion.

(If you still want them, the rules SFC use are in
site-packages/sfc/quadrature/quad_tables.py)

Martin

On Mon, Jun 22, 2009 at 11:09 PM, Shawn Walker<[email-masked]> wrote:
> Actually, would it be possible for you to just give me the 4-pt quad rule on
> a triangle.  The problem is, I don't have GiNAC setup, which is needed to
> run SyFi, right?
>
> - Shawn
>

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 17:33:23 -0400 (EDT)

Alright.  I'm looking at these quad rules, the 4-pt rule doesn't match the 
weights that show up in the Poisson.h file I have right now.  Also, the 
coordinates seem weird.  Are all of these points evaluated outside the ref 
triangle?

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 17:41:26 -0400 (EDT)

Actually, I'm a little worried about getting this initial trial to work. 
Here is the code that Poisson.h has in the inner most loop:

  // Number of operations to compute entry = 18
           A[j*6 + k] += ((Jinv_00*FE0_D10[ip][j] + 
Jinv_10*FE0_D01[ip][j])*(Jinv_00*FE0_D10[ip][k] + Jinv_10*FE0_D01[ip][k]) 
+ (Jinv_01*FE0_D10[ip][j] + 
Jinv_11*FE0_D01[ip][j])*(Jinv_01*FE0_D10[ip][k] + 
Jinv_11*FE0_D01[ip][k]))*W4[ip]*det;

My question is, can I just change Jinv_00 to Jinv_00[ip]?  That FE0_D10 
variable seems to be a special case structure.  Isn't the quadrature rule 
here built SPECIFICALLY for the case of a constant jacobian?  If I have a 
variable jacobian, won't that mess this up?

I can also just throw out the code that is there, and put in my own quad 
loop.  But that kind of defeats the purpose of this `proof-of-concept'.

- Shawn

---

From: Martin Sandve Alnæs <[email-masked]>
Date: Tue, 23 Jun 2009 00:09:09 +0200

If you're looking at SFC generated code, there
is no such thing as quad rules in the Poisson.h file.
Thus it is unclear to me what you're doing now.

There may be a translation and scaling of the
rules from the .py file to the generated code.
I think the source rules are defined on [-1, 1]
as opposed to the UFC reference cell on [0, 1].

The generated *_integral*.cpp code should have
the scaled quad rules embeddded as arrays.

Martin

On Mon, Jun 22, 2009 at 11:33 PM, Shawn Walker<[email-masked]> wrote:
>>> on
>>>>>>>> element. If you don't have a quadrature element there is no direct
>>>>>>>> way
>>>>>>>> of
>>>>>>>> No, currently not. But it might be a natural thing to have. I'll
>>>>>>>> look
>>>>>>>> into it.
>>>>>>>

---

From: Shawn Walker <[email-masked]>
Date: Mon, 22 Jun 2009 18:20:43 -0400 (EDT)

The Poisson.h file I was refering to is generated by FFC.

I guess I will wait until the FFC generated form.h files include the quad 
point locations in a comment.  Just to make sure there is no mismatch.

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 23 Jun 2009 10:31:28 +0200

> point locations in a comment.  Just to make sure there is no mismatch.

This has been added.

Kristian

> >>>>>>>>>> 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?
> >>>>>>>>>

---

From: Shawn Walker <[email-masked]>
Date: Tue, 23 Jun 2009 10:02:23 -0400 (EDT)

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