Standard way of applying a function to a finite element vector

[Lucas Myers]

Hi all,

I'd like to apply a method to a finite element function (that is, something that looks like the solution vectors) in order to get another finite element function. For example, suppose my finite element function is a symmetric rank-2 tensor in type, and I would like a finite element function which is a vector of its eigenvalues. I haven't been able to find anything in the tutorials which do this (though maybe I'm skimming too much). 

My naive approach would be to (i) construct a dealii Function which evaluates some finite element vector (function?) at a given point and then (ii) use the Interpolation or Projection functions from the VectorTools namespace. However, I am not sure whether Interpolation or Projection would be a better idea. Further, I don't know how to evaluate a finite element function at an arbitrary point -- only at quadrature points via the FEValues class. Also I would be unsure how to build that into a Function class.

Let me know if there's a better way to do this -- it seems like something that's probably already implemented.

Thanks so much,

Lucas

[Wolfgang Bangerth]

Lucas,
do I understand right that what you want is something like
v_h = f(u_h)
where f(u) is some function of the solution? If so, you're right that that
cannot hold pointwise if u_h is a finite element function and v_h should also
be a finite element function unless f has a special structure and the spaces
for u_h and v_h match. So you need
v_h = I_h f(u_h) -- interpolation
or
v_h = Pi_h f(u_h) -- projection
Which of these two you choose typically depends on the application you have.
The projection only requires you to apply f at quadrature points, and so
that's going to be easy. For I_h, you need to evaluate f(u) at the *nodal*
points of the space for v_h. That too is not complicated but you probably want
to write the interpolation function yourself, after looking at how this is
implemented in VectorTools::interpolate.

The cheap way to do the latter (but also the slow way) is to use the
FEFieldFunction class to evaluate u_h at arbitrary points. You'd then write
your own class derived from Function that gets a point as argument, calls
FEFieldFunction to evaluate u_h(x), computes f(u_h(x)), and returns that
value. This kind of function you can then provide to VectorTools::interpolate.
The problem is that FEFieldFunction is expensive. If you don't care about
speed, this is the way to go. If you do for whatever reason care about speed,
go with the approach in the previous paragraph :-)

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Lucas Myers]

Hi Wolfgang,

Thanks so much for the response, this is exactly what I was looking for! A few follow-up questions:

1. For the projection operation that you describe, are you specifically talking about this implementation of the projection() function (the [5/7] implementation)? I ask because it seems like the others require a Function object for which we need the FEFieldFunction class, rather than just something that eats a quadrature & cell index and gives back the value at the given quadrature point. If so, am I reading it correctly that this is only for a scalar-valued finite element function? How would we do this with vector-valued functions?
2. With respect to accuracy (that is, neglecting speed of computation or ease of implementation) are either of these methods better? My guess would be that projection is more (the most?) accurate in an L2 sense, while interpolation is more accurate in some other sense.
3. Is this kind of operation described in any of the tutorial steps, do you know?
4. Is this how folks typically do post-processing? Really what I would like to do is extract eigenvectors corresponding to the largest eigenvalues from my 3x3 tensor solution data and plot those as a vector field. Would this be something that is typically done with scripts in a visualization software (e.g. Paraview, Visit), or is that done in the C++ source? Just trying to get a sense for best practices...

Again, thanks so much for the help!

- Lucas

[Wolfgang Bangerth]

On 7/13/21 10:36 AM, Lucas Myers wrote:
> 1. For the projection operation that you describe, are you specifically
> talking about this implementation
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FnamespaceVectorTools.html%23a51a13e948e1296dbfdfc297167e4fd5a&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ccbbf3fb8d8c64abda0af08d9461c682f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637617910090413636%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=dqvRQzIeA4Q4JYMXQfwRjnbxzSVwpjvVJ6O79t3Txfg%3D&reserved=0> of

> the projection() function (the [5/7] implementation)? I ask because it
> seems like the others require a Function object for which we need the
> FEFieldFunction class, rather than just something that eats a quadrature &
> cell index and gives back the value at the given quadrature point. If so,
> am I reading it correctly that this is only for a scalar-valued finite
> element function? How would we do this with vector-valued functions?

If you want to use VectorTools::project(), you will have to wrap your function
into an FEFieldFunction object in the same way as I described for the
interpolation. But it's easy enough to implement a projection yourself. You
just need to build a mass matrix and solve with it.


> 2. With respect to accuracy (that is, neglecting speed of computation or ease

> of implementation) are either of these methods better? My guess would be
> that projection is more (the most?) accurate in an L2 sense, while
> interpolation is more accurate in some other sense.

We know that the L2 projection and interpolation converge at the same rate
with the mesh size h. In practice, the difference in accuracy is probably
unimportant for almost any application.


> 3. Is this kind of operation described in any of the tutorial steps, do you know?

I don't think so.


> 4. Is this how folks typically do post-processing? Really what I would like

> to do is extract eigenvectors corresponding to the largest eigenvalues
> from my 3x3 tensor solution data and plot those as a vector field. Would
> this be something that is typically done with scripts in a visualization
> software (e.g. Paraview, Visit), or is that done in the C++ source? Just
> trying to get a sense for best practices...

No, if all you care about is visualization, then it is not necessary to
project f(u) into a finite element space. What you are looking for is the
DataPostprocessor class that can be used to compute derived quantities. These
will then be output as point values that visualization programs then interpret
via bilinear/trilinear interpolation. There are numerous tutorial programs
that use DataPostprocessor. A particularly simple example is step-29.