group finite element formulation

[Joshua Hansel]

Hello,

I think what I'm trying to do is referred to as a "group finite element formulation" - suppose I want to evaluate some function of my solution u--let's call it myfunc(u)--using an FE expansion. For example, to evaluate the function at quadrature points in a cell, I don't want to evaluate the solution at quadrature points and then pass those values to the function: 

FEValues<dim> fe_values_u(fe_u, quadrature, update_values);

cell = dof_handler_u.begin_active();

fe_values_u.reinit(cell);

u_q = fe_values_u.get_function_values(u_dofs);

myfunc_q = myfunc(u_q);

Instead, I want to evaluate the function directly on the degrees of freedom of u--call them u_dofs--and then use these function degree of freedom values myfunc_dofs in an FE expansion:

myfunc_dofs = myfunc(u_dofs); // need to considering dof ordering here

somehow_associate_myfunc_dofs_with_a_dof_handler_and_fe_object();

FEValues<dim> fe_values_myfunc(fe_myfunc, quadrature, update_values);

cell = dof_handler_myfunc.begin_active();

fe_values_myfunc.reinit(cell);

myfunc_q = fe_values_myfunc.get_function_values(myfunc_dofs)

How might I go about doing this? Additional details:

• My solution has multiple components
• I want to use the same FE space (same type of FE object and degree) for myfunc as for each component of u.

For a scalar problem, I might guess that the answer of my question would be:

FE_Q<dim> fe_u(fe_degree);

DoFHandler<dim> dof_handler_u(triangulation);

dof_handler_u.distribute_dofs(fe_u);

FE_Q<dim> fe_myfunc(fe_degree);

DoFHandler<dim> dof_handler_myfunc(triangulation);

dof_handler_myfunc.distribute_dofs(fe_myfunc);

myfunc_dofs = myfunc(u_dofs);

// then just compute myfunc_q as given above

This would assume that the ordering in dof_handler_u and dof_handler_myfunc were the same, but I don't know if this is a safe assumption. For vector-valued problems though, I have no idea how to do this since I don't think the ability to extract dofs directly exists:

var1_dofs = u_dofs[var1_extractor];

var2_dofs = u_dofs[var2_extractor];

myfunc_dofs = myfunc(var1_dofs, var2_dofs);

Thanks!

P.S. I couldn't find an appropriate category tag - I think a tag "finite_elements" would be good.

[Daniel Arndt]

Hello Joshua,

you can always ask a DoF for its component using FiniteElement::system_to_component_index and hence split the local dofs in those for each component.

The ordering of the DoFs in the FESystem should with respect to every component (normally) be the same as for the scalar one.

In doubt, you can check this by comparing the mapped positions of support points 

(https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element)

and match them if necessary.

Bests

Daniel

[Wolfgang Bangerth]

Joshua,

> I think what I'm trying to do is referred to as a "group finite element
> formulation" - suppose I want to evaluate some function of my solution
> u--let's call it myfunc(u)--using an FE expansion. For example, to evaluate
> the function at quadrature points in a cell, I don't want to evaluate the
> solution at quadrature points and then pass those values to the function:
>
> FEValues<dim> fe_values_u(fe_u, quadrature, update_values);
> cell = dof_handler_u.begin_active();
> fe_values_u.reinit(cell);
> u_q = fe_values_u.get_function_values(u_dofs);
> myfunc_q = myfunc(u_q);
>
> Instead, I want to evaluate the function directly on the degrees of freedom of
> u--call them u_dofs--and then use these function degree of freedom values

> myfunc_dofsin an FE expansion:

>
> myfunc_dofs = myfunc(u_dofs); // need to considering dof ordering here
> somehow_associate_myfunc_dofs_with_a_dof_handler_and_fe_object();
> FEValues<dim> fe_values_myfunc(fe_myfunc, quadrature, update_values);
> cell = dof_handler_myfunc.begin_active();
> fe_values_myfunc.reinit(cell);
> myfunc_q = fe_values_myfunc.get_function_values(myfunc_dofs)
>
> How might I go about doing this? Additional details:

I don't think I quite understand what you want to do. Can you describe what
your goal is, rather than your attempt to code a way to get there?

My best interpretation is that you'd like to find a finite element field, say
w_h(x), where at every nodal point x_i of w_h you want that the equation

w_h(x_i) = f(u_h(x_i))

is true, with a user-defined function f(...) and u_h a vector valued finite
element field. Is this correct?

Best
W.


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

[Joshua Hansel]

Yes, that's exactly what I was trying to do. I'm currently working with the shallow water equations, and apparently this group finite element formulation is an ingredient that aids the goal of achieving a "well-balanced" scheme with respect to a steady-state. For example, our test problem is a "lake at rest" problem where the bottom topography b(x) of the lake is not flat, and we'd like to keep the lake-at-rest solution, i.e., the water level h(x)+b(x) to be constant with no fluid motion: u(x)=0. This should also hold up to initial perturbations in the water height profile.

Anyway, my solution to use group finite element formulation was to write some classes that wrap around FEValues; a finite element and DoF handler are created internally for the function of interest f. My classes employ some assumptions that I've made regarding the numbering of DoFs, based on what I observed from manually checking the DoF indices via Daniel's suggestion: given that all of my solution components use the same finite elements,

• numbering in this scenario is in order of DoF support points, i.e., all DoFs at a support point are numbered, then all DoFs at the next, etc.
• within a support point, numbering is by ordering in the finite element system FESystem,
• creating another DoFHandler on the same triangulation (but different finite element object) uses the same ordering of DoF support points.

So an example of numbering for the vector-valued solution (u,v) for Q1 elements on 3 support points, ordered (x_b,x_a,x_c) would be [u_b,v_b,u_a,v_a,u_c,v_c].

Then if I made an FE object for my scalar function f, I could assume that its DoFs are numbered like [f_b,f_a,f_c]. Similarly, I allow for vector-valued functions f.

Are these assumptions safe? Currently I just made a unit test for this to make sure that these assumptions do not become invalidated.

Anyway, if these assumptions are safe, and my classes would be useful to others, I'd be happy to share them with others, i.e., apply to push them into the deal.II repository, although they may need some work in the way of generalization to other finite element systems, since I currently only need systems of Q1 elements.

-Josh

[Wolfgang Bangerth]

On 10/13/2015 10:30 AM, Joshua Hansel wrote:
> Yes, that's exactly what I was trying to do.

OK. Then let me give you a bit of a different perspective. If you want a
finite element field w_h(x) so that

w_h(x_i) = f(u_h(x_i))

for every nodal point, then what you are saying is that you want w_h(x)
to be the interpolant of f(u_h(x)) (i.e., it *interpolated* f(u_h(x)) at
the nodal points:

w_h = I_h f(u_h)

Here, the interpolant interpolates onto the space in which w_h lives,
which may or may not be the same space the component(s) of u_h live in.

My next question would be: where do you use w_h in your mathematical
formulation?

> So an example of numbering for the vector-valued solution /(u,v)/ for
> Q1 elements on 3 support points, ordered /(x_b,x_a,x_c)/ would be
> [u_b,v_b,u_a,v_a,u_c,v_c].
>
> Then if I made an FE object for my scalar function /f/, I could

> assume that its DoFs are numbered like [f_b,f_a,f_c]. Similarly, I

> allow for vector-valued functions /f/.

Unless you specifically choose a different numbering scheme, this
assumption is correct.

[Joshua Hansel]

My next question would be: where do you use w_h in your mathematical
formulation?

I add a viscous regularization term to the continuity and momentum equations:

RHS_continuity := RHS_continuity - div(viscosity*gradient(h))

RHS_momentum := RHS_momentum - div(viscosity*gradient(h*u))

The viscosity here is entropy viscosity, which requires that I compute the entropy residual at all the cell quadrature points, as well as the entropy flux gradient jumps on the face quadrature points. Right now, I only use w_h(x) for evaluating the entropy flux gradient jumps, simply because evaluating the gradient analytically can be error-prone since the entropy flux is a relatively complicated function of the solution. Later, I expect to need w_h(x) in the context of an FCT scheme for these shallow water equations, perhaps applied to the conservation law fluxes themselves, i.e., the f(u) in the conservation law du/dt + div(f(u)) = 0. However, I'm not there yet and thus I don't know much about this yet. 

[Wolfgang Bangerth]

On 10/13/2015 11:40 AM, Joshua Hansel wrote:
>
> I add a viscous regularization term to the continuity and momentum
> equations:
>

> /RHS_continuity := RHS_continuity - div(viscosity*gradient(h))/
> /RHS_momentum := RHS_momentum - div(viscosity*gradient(h*u))//
> /
> /
> /

> The viscosity here is entropy viscosity, which requires that I compute
> the entropy residual at all the cell quadrature points, as well as the
> entropy flux gradient jumps on the face quadrature points. Right now, I

> only use /w_h(x)/ for evaluating the entropy flux gradient jumps, simply

> because evaluating the gradient analytically can be error-prone since
> the entropy flux is a relatively complicated function of the solution.

So w_h is the viscosity computed from the entropy residual of the
solution? In that case I don't understand why you need it to be a finite
element field to begin with (i.e., something you can expand in terms of
shape functions). All you need to be able to do is evaluate its values
at quadrature points. This you can do by simply computing

w(x) = f(u_h(x))

at quadrature points x=x_q. It is not necessary to interpolate f(u_h) to
a finite element space.

But maybe I miss something?

[Joshua Hansel]

I currently use w_h(x) to evaluate the entropy flux, not for the entropy viscosity itself - as far as I understand, there is no reason why I can't compute w(x) as you suggested for this purpose. Rather, I use it here as a matter of convenience: what I'm after are the gradients/divergences of the entropy flux, and evaluating this the way you suggest would require a lot of chain rule application, which would leave more room for error on my part. I believe it may also be more efficient because, assuming that f(u) takes a bit of work, I end up calling that function less if I give the entropy flux its own finite element field.

However, I do believe I will actually need w_h(x) when I make the extension to FCT for the system; I believe I will need to apply it to the conservation law fluxes at that point, though I have not read enough yet to know why.

Also, for the "lake at rest" problem where we have a nonzero bottom topography b(x), we'd like to have the water level h(x,t)+b(x) stay constant. To do this, it is necessary to interpolate b(x) in the same finite element space as h(x,t). b(x) of course is not a function of the solution u(x,t), but this is the same idea here.

-Josh