dolfin team mailing list archive

[Shawn Walker <walker@xxxxxxxxxxxxxxx>]

Here is an example of using P2 triangles for the Poisson problem.  This is 
VASTLY SUPERIOR to what I did before; it is more elegant and uses the 
Function class to represent the higher order geometry.  This is probably 
what Peter Brune was getting at in his emails awhile back, but the 
implementation is probably different.

You can download the demo here:

http://www.cims.nyu.edu/~walker/Higher_Order_Mesh_Demo.tar.gz

Just run scons and then ./demo.

HOW I DID IT:
-------------

1. First, I defined this form file:

Poisson.ufl
------------------------------------------
element = FiniteElement("Lagrange", "triangle", 2)
vector  = VectorElement("Lagrange", "triangle", 2)

v = TestFunction(element)
u = TrialFunction(element)
f = Function(element)
b = Function(vector)

c = 1E-15

a = inner(grad(v), grad(u))*dx + c*v*inner(b, grad(u))*dx
L = v*f*dx
------------------------------------------

Note the constant = 1E-15.  I did this to `fake' the FFC compiler so that 
it would ensure that `a' would accept an extra vector finite element 
function.  If I set the constant to 0.0, then the compiler would just 
ignore that extra piece (of course) but then I could not do this in the 
main C++ file:

// Create higher order map FunctionSpace
P2map::FunctionSpace V_map(mesh_t);
// Create vector function
Function Higher_Order_Map(V_map,"P2map_coefs.xml");

a.b = Higher_Order_Map; // <-- I had to `fake' the form compiler to get 
this to work


2. I then defined this form file:

P2map.ufl
------------------------------------------
element = VectorElement("Lagrange", "triangle", 2)
------------------------------------------

I use this to read in a P2 vector function from a .xml file.  This will be 
the P2 (higher order map) map.


3. I then manually modified the generated Poisson.h file to compute 
different Jinv coefficients (for the P2 map) based on the coefficients of 
the extra vector function `b' (see the Poisson.ufl).  This was an easy 
thing to do and should be easy to automate.  Without going into too many 
details, I basically just needed to add this loop within the quad loop:

      for (unsigned int r = 0; r < 6; r++)
      {
        NEW_J_00 += FE1_D10[ip][r]*w[0][nzc0[r]];
        NEW_J_01 += FE1_D01[ip][r]*w[0][nzc0[r]];
        NEW_J_10 += FE1_D10[ip][r]*w[0][nzc1[r]];
        NEW_J_11 += FE1_D01[ip][r]*w[0][nzc1[r]];
      }// end loop over 'r'

Where FE1_D?? are the precomputed basis function evaluations (for the 
gradient of a P2 polynomial) that FFC already generated!  This was not 
hard to implement, even by hand.


4. Next, I wrote up a main.cpp file that would read all this in, assemble 
the stiffness matrix, and compare the matrix entries to a stiffness matrix 
calculation that I did with my own separate MATLAB code.  The error is 
about 2.4E-12 for the quadrature rule that I used.


I know this seems hacky, but it can be cleaned up in obvious ways. The 
main idea is simple and takes advantage of the Function class that we 
already have.  All you have to do is tell the .ufl file that `a' should 
expect an extra vector function and pass that function to the `a' in the 
C++ code.  Thus, the higher order mesh data is stored as a simple vector 
Function.  There is no more crappy higher order mesh data to be stored in 
the mesh.xml file (this is what I did before).

Of course, to make this general, FFC would need to be modified, but would 
not be a big deal.  I think all the pieces are already there; it just 
needs some straight forward modification.  The hack I did with the small 
constant can easily be replaced by something more sensible (i.e. please 
let `a' take another function argument).

I would like to hear any comments people have about this.  I know people 
are busy with other things, but I would like to know if there is something 
wrong with what I did here.  It almost seems too easy after all the 
futzing around I did before.

- Shawn