dolfin team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

Yep! It's (almost) that simple.  There are still issues of what the proper
quadrature order for this kind of thing should be, but that's another thing
to talk about.  I've been getting the form-based approach to work with
corner-cases for a while here, but what you're doing would work fine for
now.  However, the way the code is generated for facet integrals will make
this a lot tougher than handling it in-form like I am.  I also prefer to be
able to edit my forms and then avoid touching the outputted c code, but this
works.

I very recently got a version of the primal solver for Unicorn working like
this.  I should go back to a much easier example and give it to you guys, I
suppose. heh.

This bring us to what we really should have here; the option to insert an
analytical coordinate function through the toolchain UFL -> FFC -> DOLFIN.
The one part that kept me from doing this a month ago after I talked to
Anders at ENUMATH was my inability to read the compiler part of FFC, which
appears to be exploded over several dozen files.  This would of course allow
for arbitrary order cell and facet jacobians.

- Peter

On Tue, Aug 18, 2009 at 4:42 PM, Shawn Walker <walker@xxxxxxxxxxxxxxx>wrote:

>
> 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<http://www.cims.nyu.edu/%7Ewalker/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
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>