ufl team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

It becomes a vector coefficient to the form.

- Peter

On Tue, Jul 7, 2009 at 12:04 AM, Shawn Walker <walker@xxxxxxxxxxxxxxx>wrote:

> 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
>
>
> 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
>>
>
> _______________________________________________
> UFL-dev mailing list
> UFL-dev@xxxxxxxxxx
> http://fenics.org/mailman/listinfo/ufl-dev
>
>