Subject: [FFC-dev] iso-parametric mappings

From: Shawn Walker <[email-masked]>
Date: Wed, 13 Aug 2008 11:56:11 -0400 (EDT)

Hello.  I know this isn't very critical on the list of things to do, but 
is there any document somewhere saying how they propose to implement 
curved meshes with the UFC/FFC code?

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Wed, 13 Aug 2008 22:10:15 +0200

> curved meshes with the UFC/FFC code?

No, I guess they have no idea how to do this. ;-)

It would require an extension of UFC and it would require someone
willing to spend the time on making it work. It's not on the top of my
TODO list, but anyone wishing to investigate is welcome.

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Wed, 13 Aug 2008 16:36:48 -0400 (EDT)

ok, but I have a couple other questions/comments?

1. Has there been a decision on what kinds of curved meshes will be 
allowed?  And what would the mesh format be?  This seems to be the first 
decision.  It would seem natural to only allow polynomial based maps 
(for now), i.e. a quadratic triangle would require a 2nd order polynomial 
mapping (instead of a 1st order map).  You could also have higher degrees 
if you wanted.  I think the 'curved' mesh feature could be implemented 
first.  And the assembly could just IGNORE the extra curvy info for each 
triangle (for now).

2. For the assembly, I think it is still possible to compute the matrices 
exactly in SOME cases.  For example, a mass matrix could be computed 
exactly, as well as \int_{\Omega} q div(u), because the det(J) would 
cancel the 1/det(J) part of the inverse of the gradient map.  Of course, 
the stiffness matrix cannot be computed exactly.  For this, I think 
quadrature is the only way.  Or one could maybe take advantage of doing an
asymptotic expansion (taylor series) of the algebraic expression and 
compute that using enough terms based on the desired number of digits of 
accuracy.  I don't know.

3. One could add another boolean property to each triangle called:

Is_Straight

This could give a switch between computing exactly, and using quadrature. 
Most meshes will have straight interior elements; you don't really need 
them inside so much.  only the triangles with a face on the boundary 
would be curved.  So this would help prevent possible performance hit of 
doing high order quadrature.

4. How hard would it be to understand UFC~FFC for someone who knows just a 
little Python and a decent amount of C++?

If there is anything written up on how FENICS wants to do this (i.e. 
desired mesh format, etc...) please tell me.

- Shawn

---

From: Anders Logg <[email-masked]>
Date: Wed, 13 Aug 2008 22:54:56 +0200

On Wed, Aug 13, 2008 at 03:45:59PM -0500, Matthew Knepley wrote:
> Before spending time implementing isoparametric stuff, I would at least look
> at the NURB elements from Tom Hughes' group. They actually capture the
> geometry (whereas polynomials rarely to a good job), although they are
> currently computationally expensive. If I wer ecoding something new, this is
> what I would do.
> 
>   Matt

There has been talk about this. I've been discussing (with our friend Trond)
to possibly extend UFC to handle the isogeometric stuff (NURB elements).
I don't know if or when this will happen, but there is interest and
also some funding (I believe).

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Wed, 13 Aug 2008 23:00:07 +0200

> compute that using enough terms based on the desired number of digits of  

I think it's be possible to do it "exactly" since FFC precomputes
integrals on the reference cell (when using its so-called tensor
representation of the form).

> doing high order quadrature.

Sounds good.

> 4. How hard would it be to understand UFC~FFC for someone who knows just 
> a little Python and a decent amount of C++?

Not that hard. Start by looking at the UFC manual. It's the only one
of our manuals that is in good shape. Then look at the generated code,
then at the form compiler.

> desired mesh format, etc...) please tell me.

No, nothing. But maybe it is better to go for isogeometric rather than
isoparametric as Matt suggests?

-- 
Anders

---

From: [email-masked]
Date: Wed, 13 Aug 2008 23:10:07 +0200 (CEST)

> ok, but I have a couple other questions/comments?

I would also very much like higher order meshes. I'll try to help out, but
I also have alot on my stack.

> compute that using enough terms based on the desired number of digits of

I don't think it is hard to get this into the form compiler (which support
quadrature).
I think the main job will be to get higher order mesh in dolfin and make sure
that the numbering makes sense in UFC.

> little Python and a decent amount of C++?

You don't need to know much Python.

> desired mesh format, etc...) please tell me.

A list of coordinates and then an element list with its node numbers. The
nodes may
be numbered according to the UFC numbering in the same way that FFC
already does
for higher order elements.

Kent

---

From: Shawn Walker <[email-masked]>
Date: Wed, 13 Aug 2008 17:21:30 -0400 (EDT)

Actually, the reason I am even considering this, is that I have a project 
that requires 2nd order polynomial mappings for triangles with a face on 
the boundary.  It is actually critical to the stability and accuracy of 
the method because there is an interaction between the velocity field on 
the boundary and the computation of curvature (because of a semi-implicit 
time discretization).  In my method, 2nd order polynomials are used for 
the velocity field, so I NEED a 2nd order piecewise parameterization of 
the boundary.

Furthermore, the problem is a two-phase free-boundary problem.  I.e. the 
geometry is NOT known a priori.  So, I am not sure how one would use the 
iso-geometric stuff; it may be overkill or problematic.  Any 
thoughts?  Also, I have already implemented my problem in MATLAB and it 
works quite well.  Of course, I could just not worry about it, but I was 
trying to think ahead.  And I wouldn't mind having my problem implemented 
in DOLFIN.  And, no, I don't want to develop another method that would 
fit with dolfin.  I think that is the wrong attitude.  :)

I think for free boundary problems (with smooth boundaries), assuming a 
polynomial curve is not bad at all.  Anyway, this could be an option in 
the code.  I would envision FFC working like this:

1. You specify whether the incoming meshes will be straight, polynomial 
curved, or iso-geometric curved.

2. You then run the compiler.

3. When you use it in a code, there will be check 'somewhere' to check for 
compatibility of the mesh with the finite element implementation.  A 
message will be output indicating any discrepancy.

- Shawn

On Wed, 13 Aug 2008, Matthew Knepley wrote:

> On Wed, Aug 13, 2008 at 3:36 PM, Shawn Walker <[email-masked]> wrote:
>> _______________________________________________
>> FFC-dev mailing list
>> [email-masked]
>> http://www.fenics.org/mailman/listinfo/ffc-dev
>>
>
>
>
> -- 
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which
> their experiments lead.
> -- Norbert Wiener
>

---

From: Jed Brown <[email-masked]>
Date: Wed, 13 Aug 2008 21:33:48 +0000

> fit with dolfin.  I think that is the wrong attitude.  :)

So you are moving the mesh as part of the nonlinear iteration or
time-stepping scheme?  I'm curious how you ensure quality (and validity)
of the parametrically mapped mesh under deformations (and possibly
topology changes).  I've always thought this was a hard problem.

Have you considered a level set method?  It would give you a high order
representation of the interface while circumventing any mesh gymnastics.

Jed

---

From: Shawn Walker <[email-masked]>
Date: Wed, 13 Aug 2008 17:37:35 -0400 (EDT)

ok, but how would you adjust the xml mesh format?  That should be figured 
out first.  I was thinking one could just add extra data for each cell 
that specifies the location of the other control points for the triangle 
sides.  Or maybe it would be better to just supply a finite element 
function (2nd order or higher) for the x,y,z coordinates of the mesh and 
do it that way?

- Shawn

---

From: Jed Brown <[email-masked]>
Date: Wed, 13 Aug 2008 21:42:55 +0000

> function (2nd order or higher) for the x,y,z coordinates of the mesh and 

I'm a fan of storing coordinates just like any other function, but you
probably want variable order (p-adaptivity) for this since you would
like most elements to have an affine mapping (for efficiency and to
improve the accuracy of the quadrature).

Jed

---

From: Shawn Walker <[email-masked]>
Date: Wed, 13 Aug 2008 17:48:44 -0400 (EDT)

> improve the accuracy of the quadrature).

Well, if DOLFIN has p-adaptive functions built into it, then this would 
not be a problem.  If it doesn't have it yet, then it could wait.  In the 
meantime, I would just have a full P2, or P3 (2nd or 3rd order) function 
for all triangles.  And, in addition, provide a boolean function on each 
cell that says whether it is straight or not.

- Shawn

---

From: "Johan Hoffman" <[email-masked]>
Date: Thu, 14 Aug 2008 12:22:45 +0200 (MEST)

> cell that says whether it is straight or not.

I am also interested in curved boundaries. If we do not want to change the
xml format, maybe it would be enough to use a boolean MeshFunction over
the boundary to achieve what you describe?

/Johan

---

From: Shawn Walker <[email-masked]>
Date: Thu, 14 Aug 2008 09:02:06 -0400 (EDT)

I'm not against changing the xml format.  Actually, it would be nice to 
have a format that is backwards compatible; that way you could ignore the 
curved information and just use the mesh as straight if you want.  That 
seems do-able.

And then add a local struct to the UFC format that is accessible by the 
tabulate_tensor routine.  This local struct would have the local nonlinear 
map for each cell and facet in place.  I need to look at the UFC manual 
some more and see how this might work...

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Thu, 14 Aug 2008 15:28:45 -0400 (EDT)

That seems reasonable.  I would propose adding a struct-variable under
ufc::cell that captures the polynomial map (without quadrature), and some 
kind of ENUM variable that indicates the order of the map.  And a boolean 
variable to ufc::cell that says whether it is straight or not.  Of course, 
a good data structure for the xml format should be chosen...

Then, I think, there would be a switch inside tabulate_tensor that chooses 
between straight and curved.  The code for straight would be the same as 
before, and for curved, it could just be quadrature (at first).

This will effect FFC.  But will it effect FIAT?  Or FErari?

- Shawn

On Thu, 14 Aug 2008, Martin Sandve Alnæs wrote:

> Without having thought this through in any way, I think this should
> rather be additions to the ufc::cell struct than adding a new struct.
>
> -- 
> Martin
>
>
>
> 2008/8/14 Shawn Walker <[email-masked]>:

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 14 Aug 2008 20:41:34 +0200

-- 
Martin

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 Aug 2008 20:46:56 +0200

> some more and see how this might work...

The mesh XML format has a data section that may contain arbitrarily
named data so it should be easy to extend (if/when we decide to
implement this).

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Sun, 17 Aug 2008 00:08:28 -0400 (EDT)

I think you are right.  I would propose the following:

1. Add a meshfunction to the mesh.xml file.  The name of it will 
correspond to a new enum variable in ufc::cell.  For example, say we have 
piecewise quadratics for the local map, so we call it "P2_map".  Maybe
later this could be made a more permanent xml data type, but for now this 
should be ok.  The mesh function would just be a vector lagrange 2nd 
order function.  And there would also be a boolean function on each cell 
to say whether the element is straight or not.  This of course rules out 
mixing curved types, i.e. having P3 and P2 in the same mesh.  But I am ok 
with that for now.

2. In ufc::cell, we have another enum variable called "map_type":

enum map_type {P1_map, P2_map, P3_map, etc...}

Of course, this will only be for lagrange type polynomial maps.  But this 
enum variable can have other types (i.e. iso-geometric).  We then create a 
variable "map_type cell_map_type;"

3. When the mesh is read in, and the element is created, I think the 
"coordinates" variable can store the vertex positions of the element, 
including the extra ones for the curved sides.  If not, then some other 
variable would need to be created to store it; actually, that may be 
better.  Also, the "cell_map_type" will be set to either P1_map, or 
P2_map, or etc...

4. The tabulate_tensor routine will have a switch to pick the correct 
routine for computing the tensor coefs.  This will require modifying FFC 
which I don't know yet.

5. Another thing that should change will be the basis function 
evaluations.  This may be a pain.  Computing the inverse map to get the 
location of the coordinates in the reference element may be inconvenient. 
This can be ignored for now; this isn't really critical for what I would 
like to do.

I am a little confused on where the ufc::cell data gets set by the data in 
the mesh.xml file.  I'm still very new to this code.  Any help would be 
appreciated.  If I am totally off on this, please say so!

- Shawn

---

From: Kent-Andre Mardal <[email-masked]>
Date: Mon, 18 Aug 2008 11:14:01 +0200

> variable "map_type cell_map_type;"

Yes, we probably would need one extra variable for the map. 

The way I see it, if one assume that all elements in the mesh has the
same order, the main work would be to create a dolfin mesh containing
a higher order geometry. Once this is done, you would need to ensure
that the coordinates in ufc::cell follows the UFC convention locally. 
Finally, the computation of the Jacobian, the inverse of the Jacobian
and the determinant of the Jacobian has to be moved inside the
quadrature loop. This is the only modification to the form compiler, I
believe. The finite elements, its derivatives, dofs etc
stay the same. 

Hence, once a higher order mesh is up an running, there should only 
be minimal modifications to the form compilers and UFC. 

> appreciated.  If I am totally off on this, please say so!

Have a look in UFCCell. 

Kent

---

From: Shawn Walker <[email-masked]>
Date: Mon, 18 Aug 2008 09:01:22 -0400 (EDT)

Sure, but I do want to have the option of some elements being curved and 
others not.  This is actually practical since only the elements near the 
boundary really need to be curved.

Also, it is possible to evaluate some forms exactly, even if the local map 
is a polynomial.  Of course, this can be ignored for now...

- Shawn

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 18 Aug 2008 15:24:51 +0200

Kents point is that there are some assumptions in the current code about
having the same element throughout the whole domain, which will require
some more changes. It will be easier to fix all the other issues first.

--
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Mon, 18 Aug 2008 15:25:15 +0200

> boundary really need to be curved.

Sure, I agree totally on this. But this requires more. Neither UFC, the
form compilers (FFC or SFC) nor the FEM basis function generators,  FIAT
and SyFi, support different polynomial degree at different places. 

Of course, in your case I guess it is only the geometrical mapping that
changes degree throughout the mesh and this may simplify alot eg. the 
finite element basis functions will be the same throughout the mesh. The
dofs can also be reused. 

I only suggest starting with the implementation of a higher-order mesh
because it seems as a feasible first step. 

Kent

---

From: Shawn Walker <[email-masked]>
Date: Mon, 18 Aug 2008 09:38:31 -0400 (EDT)

ok.  Having all the elements have 6 vertices (for P2) is not so bad.  I 
will just do that first.  I assume this only means modifying (in the 
obvious way) the mesh.xml file so that each triangle has 6 vertices 
instead of 3.  The coordinates should be automatically read in correctly. 
Of course, I will test this.

I can then (hopefully) setup FFC to just generate the right quadrature for 
this.  I'll make sure it works by comparing with my own MATLAB/C++ code. 
THEN I can put in the 'cell_map_type' and boolean function, and then 
modify FFC to deal with that.

- Shawn.

---

From: Kent-Andre Mardal <[email-masked]>
Date: Mon, 18 Aug 2008 15:43:44 +0200

> Of course, I will test this.

You might need to add something to the C++ class Mesh in Dolfin to
generate a higher order mesh. I guess Anders knows. 

Kent

---

From: Anders Logg <[email-masked]>
Date: Tue, 19 Aug 2008 18:52:41 +0200

> generate a higher order mesh. I guess Anders knows. 

Yes, a good place to start would be to think about how to extend the
Mesh class to do higher order meshes.

It has been the plan all the way to support higher order meshes. The
Mesh consists basically or two things: MeshTopology and MeshGeometry.
The MeshGeometry class is very simple and just contains a list of
vertex coordinates. One could think of either extending this class
or having separate implementations handle different type meshes.

-- 
Anders