← Back to team overview

dolfin team mailing list archive

  1. dolfin team
  2. Mailing list archive
  3. Message #14278

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

  Thread Previous • Date Previous • • Thread Next 

On Mon, 6 Jul 2009, Anders Logg wrote:


On Mon, Jul 06, 2009 at 09:34:43AM -0400, Shawn Walker wrote:


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

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


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)
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.


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.


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
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

References

  Thread Previous • Date Previous • • Thread Next