Accessing nodal values of a FEM solution

[Xuefeng Li]

Hi, there!

I have some general questions about deal.ii, and here is the background info concerning my question.

There are two functions u and v, defined over domain \Omega in 1D/2D/3D. We use deal.ii to solve for the numerical approximation of function u, using 1st degree polynomial for finite elements. That means that DOF on each cell is the same as values of function u at nodal points on each cell. In the mathematical model, we update the nodal values of function v using the nodal values of function u, pointwise, i.e., a nodal value of function v strictly depends on a nodal value of function u at the same nodal point. And that is the end of the process of solving for numerical approximations 

of functions u and v.

After reading through many tutorial examples at deal.ii, I have not come across with answers to my following questions. They are,

Is it possible to access nodal point values of a solution in deal.ii? And is it possible to revise nodal values of a solution pointwise?

The tutorial examples show only how to access values of the solution at the quadrature points within each cell.

Thanks for your insight on this question.

[Daniel Arndt]

Is it possible to access nodal point values of a solution in deal.ii? And is it possible to revise nodal values of a solution pointwise?

Yes, of course, the elements of your solution vector correspond to degrees of freedoms which are nodal values in your case.

You simply access them (both for read and write) via operator().

DoFTools::map_dofs_to_support_points (https://www.dealii.org/current/doxygen/deal.II/namespaceDoFTools.html#a5514e4f59ea659f63953d62ca429eaff)

might be useful if you need to know the support point for a given degree of freedom.

The tutorial examples show only how to access values of the solution at the quadrature points within each cell.

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

Best,

Daniel

[Wolfgang Bangerth]

On 7/19/20 6:28 PM, Daniel Arndt wrote:
>
> The tutorial examples show only how to access values of the solution at
> the quadrature points within each cell.
>
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fwiki%2FFrequently-Asked-Questions%23how-to-get-the-mapped-position-of-support-points-of-my-element&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc61491e1703e46c5be0f08d82c43defe%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637308017305040965&sdata=TM3AYA4IqcrckUMqRTWfCJWUz3MGrUx%2FRFGvRrsg40Y%3D&reserved=0>
> might also be helpful.

You might also be interested in looking at the do_half_phase_step() function
of step-58:

https://dealii.org/developer/doxygen/deal.II/step_58.html#ImplementingtheStrangsplittingsteps

Best
W.


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

[Xuefeng Li]

On Sun, Jul 19, 2020 at 7:36 PM Wolfgang Bangerth <bang...@colostate.edu> wrote:

On 7/19/20 6:28 PM, Daniel Arndt wrote:
>
>     The tutorial examples show only how to access values of the solution at
>     the quadrature points within each cell.
>
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fwiki%2FFrequently-Asked-Questions%23how-to-get-the-mapped-position-of-support-points-of-my-element&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc61491e1703e46c5be0f08d82c43defe%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637308017305040965&sdata=TM3AYA4IqcrckUMqRTWfCJWUz3MGrUx%2FRFGvRrsg40Y%3D&reserved=0>
> might also be helpful.

You might also be interested in looking at the do_half_phase_step() function
of step-58:

https://dealii.org/developer/doxygen/deal.II/step_58.html#ImplementingtheStrangsplittingsteps

Best
  W.

Thank both of you for your replies. The Step-58 function is very helpful in showing how to access nodal point values of a solution vector object. I assume I can access the nodal point coordinates in the similar way, where a vector of nodal points is created using the 

   DoFTools::map_dofs_to_support_points()

function.

However, I have two more related questions. BTW, I am a newbie in C++ programming. So my questions may seem absurd.

1. Now that we have a vector holding all nodal point coordinates, and another vector holding all nodal point values of the solution. How do we access every nodal point coordinates, and at the same time, the associated nodal point value within the same loop?
2. In addition to nodal point values of a solution, are the partial derivatives (or gradient) of the solution at each nodal point available in deal.ii?

Thanks again for any assistance!

--

Stay put, practice social distancing, and be safe!

Best,

--Xuefeng Li, (504)865-3340(phone)
   Like floating clouds, the heart rests easy
   Like flowing water, the spirit stays free

      Loyola University New Orleans

   New Orleans, Louisiana (504)865-2051(fax)

[Daniel Arndt]

However, I have two more related questions. BTW, I am a newbie in C++ programming. So my questions may seem absurd.

1. Now that we have a vector holding all nodal point coordinates, and another vector holding all nodal point values of the solution. How do we access every nodal point coordinates, and at the same time, the associated nodal point value within the same loop?

Can you clarify in a short pseudocode example what you are trying to achieve?

1. In addition to nodal point values of a solution, are the partial derivatives (or gradient) of the solution at each nodal point available in deal.ii?

That's basically what DataPostprocessor is doing. In general, you can loop over all cells and calculate the derivatives locally. Have a look at FEValuesBase::get_function_gradients (https://www.dealii.org/current/doxygen/deal.II/classFEValuesBase.html#ad1f4e0deb5d982e8172d82141c634a67).

Best,

Daniel

[Xuefeng Li]

On Thu, Jul 23, 2020 at 12:43 PM Daniel Arndt <d.arnd...@gmail.com> wrote:

However, I have two more related questions. BTW, I am a newbie in C++ programming. So my questions may seem absurd.

1. Now that we have a vector holding all nodal point coordinates, and another vector holding all nodal point values of the solution. How do we access every nodal point coordinates, and at the same time, the associated nodal point value within the same loop?

Can you clarify in a short pseudocode example what you are trying to achieve?

We need to update function v based on solution u in the following loop.

for ( x in vector_of_nodal_points )

  v(x) = f(x, u(x))  

where v(x) and u(x) are the nodal point value for functions v and u, respectively, and f() is some function depending on function u as well as the location of the nodal point. So in this case, we need to access nodal point coordinates and nodal point values of the solution in the same loop. 

 

1. In addition to nodal point values of a solution, are the partial derivatives (or gradient) of the solution at each nodal point available in deal.ii?

That's basically what DataPostprocessor is doing. In general, you can loop over all cells and calculate the derivatives locally. Have a look at FEValuesBase::get_function_gradients (https://www.dealii.org/current/doxygen/deal.II/classFEValuesBase.html#ad1f4e0deb5d982e8172d82141c634a67).

Well, the above function calculates the gradients of a finite element at the quadrature points of a cell, not at the nodal points of a cell.

Such a need arises in the following situation.

for ( x in vector_of_nodal_points )

  v(x) = g(x, u(x), grad u(x))  

[Daniel Arndt]

We need to update function v based on solution u in the following loop.

for ( x in vector_of_nodal_points )

  v(x) = f(x, u(x))  

where v(x) and u(x) are the nodal point value for functions v and u, respectively, and f() is some function depending on function u as well as the location of the nodal point. So in this case, we need to access nodal point coordinates and nodal point values of the solution in the same loop. 

You can just loop over all degrees of freedom and do

for (unsigned int i=0; i<n_dofs; ++i)

  v(i) = f(map_dof_to_support_points[i], u(i))

 

or loop over the cells yourself and do

Quadrature q(fe.get_unit_support_points());
FEValues fe_values (..., q, update_q_points);

for (const auto& cell)

  ...

 points = fe_values.get_quadrature_points();

 for (unsigned int i=0; i<n_dofs_per_cell;  ++i)

   v(local_dof_indices[i]) = f(points[i],  u(local_dof_indices));

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

 

Well, the above function calculates the gradients of a finite element at the quadrature points of a cell, not at the nodal points of a cell.

Such a need arises in the following situation.

for ( x in vector_of_nodal_points )

  v(x) = g(x, u(x), grad u(x))  

You can do similarly,

Quadrature q(fe.get_unit_support_points());
FEValues fe_values (..., q, update_q_points);

for (const auto& cell)

  ...

  points = fe_values.get_quadrature_points();

  fe_values.get_function_values(values);

  fe_values.get_function_gradients(gradients);

  for (unsigned int i=0; i<n_dofs_per_cell;  ++i)

    v(local_dof_indices[i]) = f(points[i],  values(i), gradients(i));

Best,

Daniel

[Xuefeng Li]

On Thu, Jul 23, 2020 at 5:14 PM Daniel Arndt <d.arnd...@gmail.com> wrote:

You can do similarly,

Quadrature q(fe.get_unit_support_points());
FEValues fe_values (..., q, update_q_points);

for (const auto& cell)

  ...

  points = fe_values.get_quadrature_points();

  fe_values.get_function_values(values);

  fe_values.get_function_gradients(gradients);

  for (unsigned int i=0; i<n_dofs_per_cell;  ++i)

    v(local_dof_indices[i]) = f(points[i],  values(i), gradients(i));

Yes! Declaring a quadrature on support points is a great solution. Thanks a lot!

[Wolfgang Bangerth]

On 7/23/20 12:07 PM, Xuefeng Li wrote:
>
> Well, the above function calculates the gradients of a finite element at the
> quadrature points of a cell, not at the nodal points of a cell.
> Such a need arises in the following situation.
>
> for ( x in vector_of_nodal_points )
>   v(x) = g(x, u(x), grad u(x))

It's worth pointing out, however, that for the common FE_Q elements, the
function values u(x) are continuous and so it doesn't matter how exactly you
compute u(x) at node points. On the other hand, grad u(x) is in general
discontinuous and so trying to evaluate it at node points is not actually
possible: You will either get the values from one adjacent cell or the value
from another.

In other words, if you want to compute a function that depends on 'grad u',
you need to think about what exactly you mean by that. In the formulation
above, v(x) will in general be a discontinuous function, and you need to think
about whether using FE_Q (a continuous finite element space) is really what
you want to do.

[Xuefeng Li]

On Fri, Jul 24, 2020 at 9:58 PM Wolfgang Bangerth <bang...@colostate.edu> wrote:

On 7/23/20 12:07 PM, Xuefeng Li wrote:
>
> Well, the above function calculates the gradients of a finite element at the
> quadrature points of a cell, not at the nodal points of a cell.
> Such a need arises in the following situation.
>
> for ( x in vector_of_nodal_points )
>    v(x) = g(x, u(x), grad u(x))

It's worth pointing out, however, that for the common FE_Q elements, the
function values u(x) are continuous and so it doesn't matter how exactly you
compute u(x) at node points. On the other hand, grad u(x) is in general
discontinuous and so trying to evaluate it at node points is not actually
possible: You will either get the values from one adjacent cell or the value
from another.

In other words, if you want to compute a function that depends on 'grad u',
you need to think about what exactly you mean by that. In the formulation
above, v(x) will in general be a discontinuous function, and you need to think
about whether using FE_Q (a continuous finite element space) is really what
you want to do.

Best
  W.

Indeed, grad u would be discontinuous under normal conditions when u is approximated by FE_Q. I remember vaguely that such an issue was discussed in one or more of the tutorial Step examples. Thanks for your follow-up.