Subject: [FFC-dev] ordering of P2 nodes

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

Hello.  I need to know how FFC orders P2 Lagrange nodes.  I need to know 
this in order to compare a FEM matrix that I created in another program.

For example, say you have one triangle with the following nodes:

[0 1 2 3 4 5]

Where is the #3 node in relation to the #0, #1, and #2 nodes?  Is #3 
opposite the #0 node?

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 23 Jun 2009 18:01:18 +0200

> For example, say you have one triangle with the following nodes:

Yes.

I usually look at the tabulate_dofs() function, for the P2 case:

dofs[0] = c.entity_indices[0][0];
dofs[1] = c.entity_indices[0][1];
dofs[2] = c.entity_indices[0][2];
unsigned int offset = m.num_entities[0];
dofs[3] = offset + c.entity_indices[1][0];
dofs[4] = offset + c.entity_indices[1][1];
dofs[5] = offset + c.entity_indices[1][2];

So you see that the first 3 dofs are the vertices, the next 3 are the midpoints
of the edges. Then you need to know the UFC ordering of edges to figure out
that e0 is the edge which does not contain vertex 0 (opposite to v0) etc.

Kristian

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

---

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

ok.  But I have one more question.

Is there any kind of re-ordering of nodes that dolfin does internally?  I 
thought I remembered seeing this.  I have been trying to compare two 
matrices, and it looks like the entries are the same, except the ordering 
is a little different.

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 23 Jun 2009 18:15:15 +0200

> matrices, and it looks like the entries are the same, except the ordering

DOLFIN reorders the mesh such that it follows the UFC convention.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Tue, 23 Jun 2009 18:15:37 +0200

DOLFIN assumes that the mesh is ordered locally so that the local
vertices on each triangle are ordered increasingly. Take a look at the
XML files in data/meshes in DOLFIN and you'll see what I mean.

-- 
Anders

---

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

I've almost got this working.  When I compare the stiffness matrix from my 
code to the one from dolfin the max difference in the entries is machine 
precision.  This is when the mesh is straight (no curves).  If I have a 
slightly curved mesh, the error is 0.0162, but I think this is because the 
quad rule that FFC generated only has 4 points and I only have two 
triangles with diameter approx 0.5 (not small enough resolve the error).

Is there a way to tell FFC to use more quad points than it `should'?  That 
way I can verify that I am computing this properly.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Tue, 23 Jun 2009 22:36:09 +0200

Add -f quadrature_order=n when compiling for a high enough value of n.

-- 
Anders

On Tue, Jun 23, 2009 at 04:28:44PM -0400, Shawn Walker wrote:
> I've almost got this working.  When I compare the stiffness matrix from 
> my code to the one from dolfin the max difference in the entries is 
> machine precision.  This is when the mesh is straight (no curves).  If I 
> have a slightly curved mesh, the error is 0.0162, but I think this is 
> because the quad rule that FFC generated only has 4 points and I only 
> have two triangles with diameter approx 0.5 (not small enough resolve the 
> error).
>
> Is there a way to tell FFC to use more quad points than it `should'?  
> That way I can verify that I am computing this properly.
>
> - Shawn
>

---

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

Hello, I have just one minor request.  When putting the quad points in the 
comment, would it be possible to break it apart into two lists (one for 
the X and one for the Y)?  that way I can just copy it into the code I am 
testing.  Otherwise, I have to copy multiple pieces separately, which can 
be tedious if there are lots of points.  This isn't really critical, 
though it would be nice.

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 24 Jun 2009 09:57:37 +0200

> be tedious if there are lots of points.  This isn't really critical,

Yes, that is possible, but I think it is more natural to just have a list of
points.
However, I added some code to

ffc/compiler/codegeneration/quadraturegenerator.py around line 420

which will do what you want. Just uncomment the lines.

We could have a quick vote here on the list to which format we want (or we could
simply include both formats in the generated code).

I don't want to push any changes before we agree on something because I will
have to update the references if I do.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Wed, 24 Jun 2009 10:01:40 +0200

> have to update the references if I do.

I also think the first option (x, y, z), (x, y, z) is more natural.

But if it doesn't take up much space, we could have both. I have no
strong opinion.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Wed, 24 Jun 2009 09:04:02 -0400 (EDT)

> I also think the first option (x, y, z), (x, y, z) is more natural.

Well, what I meant was have something like:

x0, x1, x2, ...
y0, y1, y2, ...

Another way would be to put curly braces:

{ {x0, y0, z0}, {x1, y1, z1}, ...  }

That way I could use that define the array directly.

- Shawn

---

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 24 Jun 2009 15:20:46 +0200

> Well, what I meant was have something like:

Yes, as I wrote above I've put some code in

which does just that. It is not currently active so you need to uncomment the
lines.

The reason I did not make it active is that I would have to update the
references for the regression tests and I don't want to do that because I think
what we have now is more natural.

> That way I could use that define the array directly.

I could change the format to this if you prefer.

Kristian