Subject: [UFL-dev] A form-based approach to isoparametry

From: Peter Brune <[email-masked]>
Date: Mon, 6 Jul 2009 04:44:17 -0500

I'm working on some problems with sub/super/isoparametric elements and have
very quickly implemented it entirely using UFL and a transformation of the
form to include a geometric coefficient.  This is done using the
transformation framework to append the geometric information to the parts of
the form that require transformation.

I had previously written calculations to include the coefficient space for
geometry by hand, but I talked with Anders at ENUMATH and he told me more
about the Transformer machinery.  Right now, a Poisson form with this looks
like:

cell = triangle

iso_element = VectorElement("Lagrange", cell, 4)
aff_element = VectorElement("Lagrange", cell, 1)
element = FiniteElement("Lagrange", cell, 4)

iso_func = Function(iso_element)
affine_func = Function(aff_element)

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

b = inner(grad(u), grad(v))
K = v*f

J = dot(inv(grad(iso_func)), grad(affine_func))
detJ = det(J)

a = apply_geometry(b, J)
L = apply_geometry(K, J)
a = a*dx
L = L*dx

Where apply_geometry applies a Transformer to the form in order to include
the geometric coefficients.

The benefits of this approach are:
1. General function spaces for the geometry -- no reliance on the somewhat
contradictory concept of "extra vertices." on a simplex.
2. Uses all the already existing mechanisms for compilation and optimization

Right now I'm stuck transforming the affinely-transformed components back to
the reference and applying the map.  The ideal would be appending something
to the measure, which is then appended to the form.  This might look like:

J = grad(iso_func)
a = inner(grad(u), grad(v))*dx(0) + inner(grad(u), grad(v))*dx(1, jacobian =
J)
L = v*f*dx(0) + v*f*dx(1, jacobian = J)

With the form compiler then omitting the generation of the affine Jacobian.
Like this we can easily have the higher-order geometry only defined on, say,
boundary cells where we have a higher-order geometry defined.  Otherwise the
affine form can be used.

Thoughts on how I should go about this?  I'm still generalizing my
transformer, but have run a demos of simple forms (Poisson, Stokes) with no
real problems; now I'm moving onto what I actually want to do like this, but
improving the interface would be nice eventually.

- Peter Brune

---

From: Shawn Walker <[email-masked]>
Date: Mon, 6 Jul 2009 09:34:43 -0400 (EDT)

On Mon, 6 Jul 2009, Peter Brune wrote:

> I'm working on some problems with sub/super/isoparametric elements and have very quickly implemented it
> entirely using UFL and a transformation of the form to include a geometric coefficient.  This is done using
> the transformation framework to append the geometric information to the parts of the form that require
> transformation.
> 
> I had previously written calculations to include the coefficient space for geometry by hand, but I talked
> with Anders at ENUMATH and he told me more about the Transformer machinery.  Right now, a Poisson form with
> this looks like:
> 
> cell = triangle
> 
> Where apply_geometry applies a Transformer to the form in order to include the geometric coefficients.
> 
> 1. General function spaces for the geometry -- no reliance on the somewhat contradictory concept of "extra
> vertices." on a simplex.

How exactly do you represent the geometry as a function space?  We had 
thought about this before and there was a problem (because you need a mesh 
to create a function space in the furst place).

> Right now I'm stuck transforming the affinely-transformed components back to the reference and applying the
> map.  The ideal would be appending something to the measure, which is then appended to the form.  This might
> look like:
> 
> J = grad(iso_func)
> a = inner(grad(u), grad(v))*dx(0) + inner(grad(u), grad(v))*dx(1, jacobian = J)
> With the form compiler then omitting the generation of the affine Jacobian.  Like this we can easily have the
> higher-order geometry only defined on, say, boundary cells where we have a higher-order geometry defined. 
> Otherwise the affine form can be used.

This is great.  There is some support for reading in higher order meshes 
in DOLFIN now.  There is even a boolean parameter for saying which 
triangles are curved and which are affine.

Please keep in mind that you may want to have triangles that are NOT on 
the boundary to also be curved.  This is necessary if the mesh is highly 
anisotropic (picture a wing with a very anisotropic curved mesh at the 
boundary).

> Thoughts on how I should go about this?  I'm still generalizing my transformer, but have run a demos of
> simple forms (Poisson, Stokes) with no real problems; now I'm moving onto what I actually want to do like
> this, but improving the interface would be nice eventually.
> 
> - Peter Brune

So, you have done this in UFL?  Does that mean just the notation is setup? 
I had done a hand modification of a poisson demo that reads in a higher 
order mesh and computes the stiffness matrix on two triangles (only one is 
curved) and compare it to a stiffness matrix computed by other means.  I 
don't know if this will be useful for you; what I did was a little hacky. 
Have you thought about how the higher order mesh would be stored in an XML 
file?  I can resend the example I made (that shows this), but it should be 
in the archive on DOLFIN.

I am traveling right now, so I won't be able to say much on this, 
unfortunately.

- Shawn

---

From: Peter Brune <[email-masked]>
Date: Mon, 6 Jul 2009 09:20:29 -0500

On Mon, Jul 6, 2009 at 8:34 AM, Shawn Walker <[email-masked]> wrote:

>> geometry by hand, but I talked
>> with Anders at ENUMATH and he told me more about the Transformer
>> machinery.  Right now, a Poisson form with
>> this looks like:
>>
>> cell = triangle
>>

For the moment my examples involve creating a trivial coordinate function
and transforming it to match the desired geometry.  This may involve
projection onto a surface, or using the coordinates as parameters of a new
coordinate function. (think: warping a relatively smooth bluff into a square
mesh)  This is then easily interpolated into a discrete function of any
order.  Because we can analytically construct the Jacobian from this, we can
use the transformation machinery in UFL to fold the geometry into the form.

>> With the form compiler then omitting the generation of the affine
>> Jacobian.  Like this we can easily have the
> This is great.  There is some support for reading in higher order meshes in
> DOLFIN now.  There is even a boolean parameter for saying which triangles
> are curved and which are affine.
>
> Please keep in mind that you may want to have triangles that are NOT on the
> boundary to also be curved.  This is necessary if the mesh is highly

Yes, of course!  Because I intend on defining it in terms of per-subdomain
measures, we can even vary the order of the geometry from mesh-piece to
mesh-piece (by having per-order parts in the form).  My present demos don't
distinguish between affine and curved (that I'll clean up and get out) do
things like wrap a pipe around the origin, considering the "reference"
coordinates to be r and theta.  I could just as easily add the affine part
to the form.  I'll try to make this part of my cleaned-up example.

>> transformer, but have run a demos of

I'm working on some good test examples right now, but the demos I have will
involve superparametric vs. a ludicrously refined mesh.  I'm able to run
both Poisson and Stokes examples (which I'm cleaning up and error
checking).  There could be a number of ways of reading these in, from
reading in a "higher-order" mesh format directly into a function space on
the "boundary-cell" subdomain and an affine space merely defined by the
geometry as usual, to storing the coordinates separately.  Our goal is to do
some trivial projection to an "exact" geometry.

> I am traveling right now, so I won't be able to say much on this,

- Peter

---

From: Anders Logg <[email-masked]>
Date: Mon, 6 Jul 2009 17:28:31 +0200

On Mon, Jul 06, 2009 at 09:34:43AM -0400, Shawn Walker wrote:
> thought about this before and there was a problem (because you need a 
> mesh to create a function space in the furst place).
>
> So, you have done this in UFL?  Does that mean just the notation is 
> setup? I had done a hand modification of a poisson demo that reads in a 
> higher order mesh and computes the stiffness matrix on two triangles 
> (only one is curved) and compare it to a stiffness matrix computed by 
> other means.  I don't know if this will be useful for you; what I did was 
> a little hacky. Have you thought about how the higher order mesh would be 
> stored in an XML file?  I can resend the example I made (that shows 
> this), but it should be in the archive on DOLFIN.
>

I'm on vacation so I also won't have much to say about this right now,
but the form-based approach is very appealing since (1) it is a simple
layer on top of existing functionality (so we don't need to modify
the form compilers) and (2) it is more general since we can use any
Function to map the mesh from the reference geometry.

We still have the problem of reading/storing the geometry and we
concluded a while back that we didn't want to use Function for this,
since it would create a circular dependency between the Mesh and
Function classes. But if the mapping of the mesh can be stored
separately from the mesh (as a coefficient in the form) then I guess
this is no longer a problem. (?)

-- 
Anders

---

From: Shawn Walker <[email-masked]>
Date: Mon, 6 Jul 2009 17:53:23 -0400 (EDT)

> this is no longer a problem. (?)

I would like to make a comment.  It sounds like you have some specific 
geometry in mind that you (Peter) want to map to, like a circle, or 
whatever.  Maybe I have this wrong, but that is what it sounded like.

However, there is a class of problems that I have talked about a long time 
ago where this is not the case.  For a free boundary problem, the geometry 
is an unknown.  It is part of the evolution problem, and in some cases 
iso-parametric or higher order elements are necessary.

So, all I am saying is that whatever format you come up with, it should be 
able to handle these other cases.  I should be able to update the geometry 
of the mesh easily.  Maybe this is fine in what you are doing.  I just 
remember getting confused over having a function space (with a function) 
represent the higher order geometry.  Maybe it is fine, but I have to 
think about it again.

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Mon, 6 Jul 2009 23:32:14 -0400 (EDT)

> this is no longer a problem. (?)

Actually, I think Brune has the right idea (with some slight 
modification).  We were not thinking of the best way to implement general 
non-affine maps.  I will illustrate.  I hope P. Brune does not take 
offense.

Originally, we thought that in order to properly define the function 
space, the incoming mesh must be higher order to begin with.  But, if we 
wanted to use a function space for representing the higher order mesh, 
then we have a circular dependency.

Brune's idea circumvents this (with a small amount of cheating).  The 
essential point is to take your standard affine mesh and define your 
function space V (given the mesh) to be used in the bilinear form just as 
we do now.  But then also define another function space G that is 
completely unrelated to the first one and define a function g on G.  g 
will be a vector function representing the higher order mesh geometry. 
Then, when defining the actual form variable, you pass it g in addition. 
And then it is up to FFC to compute the local jacobian given g.

This way the higher order mesh is treated as a perturbation of the 
bilinear form.  However, it is a slight cheat because the function space V 
is not defined correctly; it assumes a straight mesh.  This is a problem 
if you are interpolating points that lie in the true curved mesh.

So I propose the following modification.

1. Create a member variable, in FunctionSpace, that is a pointer (or 
something) to a FunctionSpace.  I admit this sort of sounds circular, but 
there must be a way to do it.

2. Now create the geometry function space G.

3. Now, when you create the function space V (that will be used to define 
the bilinear form) you pass it G (and a function g) to the constructor.

4. Now V has access to the higher order mesh data and can use it to do 
interpolation correctly.  And when the bilinear form is defined, it will 
also have access to the higher order mesh data and can use that when 
assembling the matrix.

I think this could work. The main point is to make sure V contains the 
higher order mesh function (and space).  One way to avoid that pointer in 
(1) would be to make a class called GeometrySpace that is similar to 
FunctionSpace but has less functionality (i.e. there is no need to 
interpolate).

Does this make sense?

- Shawn

---

From: Shawn Walker <[email-masked]>
Date: Tue, 7 Jul 2009 01:04:32 -0400 (EDT)

Actually, I just realized a problem with what I was proposing.

When you assemble the matrix, you loop through each triangle.  But the 
higher order mesh function is defined on a different function space.  How 
does the tabulate_tensor routine know which triangle to access in the 
higher order mesh function?

- Shawn

---

From: Peter Brune <[email-masked]>
Date: Tue, 7 Jul 2009 01:16:51 -0500

It becomes a vector coefficient to the form.

- Peter

On Tue, Jul 7, 2009 at 12:04 AM, Shawn Walker <[email-masked]>wrote:

> does the tabulate_tensor routine know which triangle to access in the higher
> order mesh function?
>
> - Shawn
>
>
>>>> entirely using UFL and a transformation of the form to include a
>>>> geometric coefficient.  This is done using
> _______________________________________________
> UFL-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufl-dev
>
>

---

From: Shawn Walker <[email-masked]>
Date: Tue, 7 Jul 2009 09:25:30 -0400 (EDT)

ok.  So, this is why you had those two spaces defined in the .ufl file. 
That way the form knows how to handle it.

This sounds pretty good!  :)

- Shawn

On Tue, 7 Jul 2009, Peter Brune wrote:

> On Tue, Jul 7, 2009 at 12:04 AM, Shawn Walker <[email-masked]> wrote:
>       When you assemble the matrix, you loop through each triangle.  But the higher order mesh function
>       is defined on a different function space.  How does the tabulate_tensor routine know which
>       triangle to access in the higher order mesh function?
>
>       - Shawn
> 
> 
>                   Right now I'm stuck transforming the affinely-transformed
>                   components back to the reference and applying the
>                   simple forms (Poisson, Stokes) with no real problems; now I'm
>                   moving onto what I actually want to do like

---

From: Peter Brune <[email-masked]>
Date: Tue, 7 Jul 2009 08:46:36 -0500

Ah, let me explain further: The reason I have the two vector spaces --
affine and iso, is so that I can reverse the affine transformation before
applying the isoparametric transformation.  Today I found a small bug in my
isoparametric transformer, and now I'm sure stokes works; the pretty example
showing the differences between refined affine and unrefined isoparametric
is coming along!  We want to get rid of having to transform back to the
reference using the affine coefficient space before this is ready to be
used.  I have been looking through UFL and FFC for a good way to drop in a
new jacobian.  Passing the new jacobian to FFC is really easy using
metadata... using it is harder.

So, anyone else have any tips on that?  Here's the MAJOR issue I can see
cropping up -- I can pepper the form itself with the new Jacobian in order
to transform function (derivatives) affinely, but other transformations
(co/contravariant Piola anyone?) will not work like this.  I NEED to give
FFC something that can be manipulated like the pullback/pushforward between
geometries to use the more interesting element types, or even more
interesting forms.   I need to be able to extract information like how a
given basis function (and its derivatives) transform at the UFL
form-transformation level, OR I need to be able to change the Jacobian
(per-measure).

Peter

On Tue, Jul 7, 2009 at 8:25 AM, Shawn Walker <[email-masked]> wrote:

>>                  I'm working on some problems with sub/super/isoparametric
>> elements
>>                  and have very quickly implemented it
>>                  the transformation framework to append the geometric
>> information to
>>                  the parts of the form that require
>>                  transformation.
>>
>>                  I had previously written calculations to include the
>> coefficient
>>                  space for geometry by hand, but I talked
>>                  Where apply_geometry applies a Transformer to the form in
>> order to
>>                  include the geometric coefficients.
>>
>>                  1. General function spaces for the geometry -- no
>> reliance on the
>>                  somewhat contradictory concept of "extra
>>                  vertices." on a simplex.
>>                  2. Uses all the already existing mechanisms for
>> compilation and
>>                  optimization
>>
>>
>>                  map.  The ideal would be appending something to the
>> measure, which
>>                  is then appended to the form.  This might
>>                  look like:
>>
>>                  J = grad(iso_func)
>>                  a = inner(grad(u), grad(v))*dx(0) + inner(grad(u),
>> grad(v))*dx(1,
>>                  jacobian = J)
>>                  higher-order geometry only defined on, say, boundary
>> cells where we
>>                  have a higher-order geometry defined.

---

From: Shawn Walker <[email-masked]>
Date: Tue, 7 Jul 2009 11:50:56 -0400 (EDT)

On Tue, 7 Jul 2009, Peter Brune wrote:

> Ah, let me explain further: The reason I have the two vector spaces -- affine and iso, is so that I can
> reverse the affine transformation before applying the isoparametric transformation.  Today I found a small
> bug in my isoparametric transformer, and now I'm sure stokes works; the pretty example showing the
> differences between refined affine and unrefined isoparametric is coming along!  We want to get rid of having
> to transform back to the reference using the affine coefficient space before this is ready to be used.  I
> have been looking through UFL and FFC for a good way to drop in a new jacobian.  Passing the new jacobian to
> FFC is really easy using metadata... using it is harder.

I'm a little unclear why you need to reverse the transformation.

> So, anyone else have any tips on that?  Here's the MAJOR issue I can see cropping up -- I can pepper the form
> itself with the new Jacobian in order to transform function (derivatives) affinely, but other transformations
> (co/contravariant Piola anyone?) will not work like this.  I NEED to give FFC something that can be
> manipulated like the pullback/pushforward between geometries to use the more interesting element types, or
> even more interesting forms.   I need to be able to extract information like how a given basis function (and
> its derivatives) transform at the UFL form-transformation level, OR I need to be able to change the Jacobian
> (per-measure).
> 
> Peter

It sounds like this is more complicated than it needs to be.  But maybe I 
am missing something.

Currently, in the FFC generated form file, the jacobian inverse terms 
enter into the quadrature loop as constants.  All one has to do is replace 
those constants by the variable jacobian terms that are determined by your 
local map.   I did a manual modification of this and it seemed to work.

My view of this is that FFC would need to be extended to generate this 
code, which shouldn't be too bad given that you can pass a 
Function/FunctionSpace to the form file that represents the higher order 
mesh geometry.

The main point here is that you can pass a Function and FunctionSpace as a 
separate object to the form file, which we missed before.

- Shawn

>       ok.  So, this is why you had those two spaces defined in the .ufl file. That way the form knows
>       how to handle it.
>

---

From: Peter Brune <[email-masked]>
Date: Tue, 7 Jul 2009 11:43:17 -0500

Reversing the affine transformation is a temporary measure; it's because I
haven't modified FFC to stop generating the affine jacobian.

It all depends on what level things are done at as to what gets complicated;
It seems that generating this code gets complicated either way, so here's
how it looks:

My present musings work entirely within UFL, and I can easily add geometric
terms at the form level like this, but will only work with
affinely-transforming elements -- a large class, but we're going for the
general case.  There are two ways we could go from there -- the one I was
considering with this question was keeping this in UFL and extracting the
transformation type from the basis functions (for non-affine basis
functions, at least).  This appears to be possible, but really quite a lot
of code, as done in the likes of the quadraturetransformer in FFC.

The other option is to generate code for the jacobian in FFC/Syfi/etc.  I
examined the generation code and determined that the second is probably the
best option for generality's sake.  However, it's still not pretty.  There's
a lot of code wrapped up in the affine transformation from FIAT onward, as
you are probably quite aware. :)

- Peter

On Tue, Jul 7, 2009 at 10:50 AM, Shawn Walker <[email-masked]>wrote:

>> affine and iso, is so that I can
>> cropping up -- I can pepper the form
>> FFC something that can be
> Currently, in the FFC generated form file, the jacobian inverse terms enter
> into the quadrature loop as constants.  All one has to do is replace those
> constants by the variable jacobian terms that are determined by your local
> map.   I did a manual modification of this and it seemed to work.
>
>>           When you assemble the matrix, you loop through each triangle.
>>  But the higher order
>>      mesh function
>>           is defined on a different function space.  How does the
>> tabulate_tensor routine know
>>      which
>>                       I'm working on some problems with
>> sub/super/isoparametric elements
>>                       Where apply_geometry applies a Transformer to the
>> form in order to

---

From: Shawn Walker <[email-masked]>
Date: Tue, 7 Jul 2009 14:12:51 -0400 (EDT)

On Tue, 7 Jul 2009, Peter Brune wrote:

> Reversing the affine transformation is a temporary measure; it's because I haven't modified FFC to stop
> generating the affine jacobian.
> 
> It all depends on what level things are done at as to what gets complicated; It seems that generating this
> code gets complicated either way, so here's how it looks: 
> 
> My present musings work entirely within UFL, and I can easily add geometric terms at the form level like
> this, but will only work with affinely-transforming elements -- a large class, but we're going for the
> general case.  There are two ways we could go from there -- the one I was considering with this question was
> keeping this in UFL and extracting the transformation type from the basis functions (for non-affine basis
> functions, at least).  This appears to be possible, but really quite a lot of code, as done in the likes of
> the quadraturetransformer in FFC.
> 
> The other option is to generate code for the jacobian in FFC/Syfi/etc.  I examined the generation code and
> determined that the second is probably the best option for generality's sake.  However, it's still not
> pretty.  There's a lot of code wrapped up in the affine transformation from FIAT onward, as you are probably
> quite aware. :)
> 
> - Peter

Actually, I am not that aware of the innards of FFC.  I am a humble man! 
:)    But I thought that since the application of the jacobian is already 
structurally there, having a variable jacobian is not such a big deal.  We 
need an FFC person to come to our aid!  I was also hoping Anders would 
comment on this.

- Shawn

> On Tue, Jul 7, 2009 at 10:50 AM, Shawn Walker <[email-masked]> wrote:
>
> It sounds like this is more complicated than it needs to be.  But maybe I am missing something.
> 
> Currently, in the FFC generated form file, the jacobian inverse terms enter into the quadrature loop as
> constants.  All one has to do is replace those constants by the variable jacobian terms that are
> determined by your local map.   I did a manual modification of this and it seemed to work.
> 
> My view of this is that FFC would need to be extended to generate this code, which shouldn't be too bad
> given that you can pass a Function/FunctionSpace to the form file that represents the higher order mesh
> geometry.
> 
> The main point here is that you can pass a Function and FunctionSpace as a separate object to the form
> file, which we missed before.
> 
> - Shawn
> 
>
>                 When you assemble the matrix, you loop through each triangle.  But the higher
>       order
>            mesh function
>                 is defined on a different function space.  How does the tabulate_tensor routine
>       know
>            which

---

From: Anders Logg <[email-masked]>
Date: Tue, 7 Jul 2009 22:20:06 +0200

On Tue, Jul 07, 2009 at 02:12:51PM -0400, Shawn Walker wrote:
> structurally there, having a variable jacobian is not such a big deal.  
> We need an FFC person to come to our aid!  I was also hoping Anders would 
> comment on this.

I'm on vacation so I only have time to reply when something is either
(1) extremely urgent or (2) simple to respond to. This is neither so
it will have to wait... :-)

-- 
Anders

---

From: Peter Brune <[email-masked]>
Date: Tue, 7 Jul 2009 15:32:36 -0500

It's cool!  I'm making baby steps towards being able to do this.  In the
mean time it seems that I can compute a lot of things with what I have now
(I'll send a patch once my demo is nice enough)  You pretty much pointed me
towards where to look at ENUMATH.  Have a good vacation.

- Peter

On Tue, Jul 7, 2009 at 3:20 PM, Anders Logg <[email-masked]> wrote:

> >> It all depends on what level things are done at as to what gets
> complicated; It seems that generating this
> >>       Ah, let me explain further: The reason I have the two vector
> spaces -- affine and iso, is
> >>       so that I can
> >>       reverse the affine transformation before applying the
> isoparametric transformation.  Today
> >>       I found a small
> >>       bug in my isoparametric transformer, and now I'm sure stokes
> works; the pretty example
> >>       showing the
> >>       differences between refined affine and unrefined isoparametric is
> coming along!  We want to
> >>       get rid of having
> >>       to transform back to the reference using the affine coefficient
> space before this is ready
> >>       to be used.  I
> >>       have been looking through UFL and FFC for a good way to drop in a
> new jacobian.  Passing
> >>       the new jacobian to
> >>       So, anyone else have any tips on that?  Here's the MAJOR issue I
> can see cropping up -- I
> >>       can pepper the form
> >>       itself with the new Jacobian in order to transform function
> (derivatives) affinely, but
> >>       other transformations
> >>       (co/contravariant Piola anyone?) will not work like this.  I NEED
> to give FFC something
> >>       that can be
> >>       manipulated like the pullback/pushforward between geometries to
> use the more interesting
> >>       element types, or
> >>       even more interesting forms.   I need to be able to extract
> information like how a given
> >>       basis function (and
> >>       its derivatives) transform at the UFL form-transformation level,
> OR I need to be able to
> >>       change the Jacobian
> >>       (per-measure).
> >>
> >>       Peter
> >>
> >>
> enter into the quadrature loop as
> code, which shouldn't be too bad
> >>            ok.  So, this is why you had those two spaces defined in the
> .ufl file. That way the
> >>       form knows
> >>            how to handle it.
> >>
> >>                             the transformation framework to append the
> geometric information to
> >>                             I had previously written calculations to
> include the coefficient
> >>                             Where apply_geometry applies a Transformer
> to the form in order to
> >>                             1. General function spaces for the geometry
> -- no reliance on the
> >>                             2. Uses all the already existing mechanisms
> for compilation and
> >>                             optimization
> >>
> >>
> >>                             map.  The ideal would be appending something
> to the measure, which
> >>                             a = inner(grad(u), grad(v))*dx(0) +
> inner(grad(u), grad(v))*dx(1,
> >>                             jacobian = J)
> >>                             higher-order geometry only defined on, say,
> boundary cells where we
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkpTrfYACgkQTuwUCDsYZdGOvQCdEjt/J8LsaGmtrgQyRzQvH4Wq
> fMQAn29qPJEi4c0Ns7rouOBwAZ6CFnoX
> =iCb3
> -----END PGP SIGNATURE-----
>
>