Is there a way to traverse all DoFs?

[Yang Liu]

Dear deal developers and users,

I would like to ask is there a way to traverse all DoFs in a FE? I would also like to save them in a particular order so that I can associate each DoF with a coefficient and later use them in the matrix assemblage. 

Best Regards,

Yang

[Wolfgang Bangerth]

On 4/28/20 5:52 AM, Yang Liu wrote:
>
> I would like to ask is there a way to traverse all DoFs in a FE?

Well, you can always do
for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
...do something...
to traverse all degrees of freedom. But I suspect that's not what you want,
but that you want to tie them to individual cells, individual components, etc.

Can you describe what you actually want to do with this traversal?

Best
W.


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

[Yang Liu]

Dear Wolfgang,

I would like to generate a covariance matrix based on these DoFs. After that I would like to produce a random vector and associate the DoFs with the random vectors, in order to compute the quadrature. 

Best Regards,

Yang

在 2020年4月28日星期二 UTC+8下午10:43:50，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 4/28/20 8:50 AM, Yang Liu wrote:
>
> I would like to generate a covariance matrix based on these DoFs. After that I
> would like to produce a random vector and associate the DoFs with the random
> vectors, in order to compute the quadrature.

I think that's still too abstract. For example, you can create a random vector
using something like this:

Vector<double> r(dof_handler.n_dofs());

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

r(i) = std::rand();

I suspect you have a more concrete application in mind, but you need to be
more specific to explain what it is for us to give you suggestions of how to
implement that.

[Yang Liu]

Hi Wolfgang,

I mean I would generate a n by n covariance matrix, based on the locations of these DoFs. 

在 2020年4月28日星期二 UTC+8下午11:45:41，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 4/28/20 9:48 AM, Yang Liu wrote:
>
> I mean I would generate a n by n covariance matrix, based on the locations of
> these DoFs.

Right, but I don't know what that means. You have a very specific application
in mind but I don't know what it is. You have to explain better what that is.
For example, write down a piece of pseudo-code that shows what you want to do,
and then we can think about how that could be implemented with deal.II. But I
have no idea *of what* you want to compute the covariance matrix, and how the
locations of DoFs enter into this.

[Yang Liu]

I mean the following:

Assume we already have a vector D of size n, which contains all the DoF points (quadrature points).

Next I would like to assemble a n by n covariance matrix A, 

for (int i = 0; i < n; ++i)

    for (int j = 0; j < n; ++j)

    {

        A[i][j] = f(D[i], D[j]); // f is some function depending on two points

     }

After that I perform Cholesky decomposition $A = LL^{*}$ and compute $y = Lz$, where $z$ is some other vectors. 

I would also like to associate the vector $y$ with the quadrature points since I would like to use them later, e.g., 

for (const auto &cell : dof_handler.active_cell_iterators())

{

for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)

{

const double current_coefficient = y[quadrature_point(q_index).global_index()]

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

{

for (unsigned int j = 0; j < dofs_per_cell; ++j)

cell_matrix(i, j) +=

(current_coefficient * // a(x_q)

fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)

fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)

fe_values.JxW(q_index)); // dx

在 2020年4月28日星期二 UTC+8下午11:52:30，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

> Next I would like to assemble a n by n covariance matrix A,
>
> for (int i = 0; i < n; ++i)
>     for (int j = 0; j < n; ++j)
>     {
>         A[i][j] = f(D[i], D[j]); // f is some function depending on two points
>      }

I see. So for this, you need the support points of all degrees of freedom
(which you denote by D[i] above). You can get those by calling
DoFTools::map_dofs_to_support_points().

> I would also like to associate the vector $y$ with the quadrature points since
> I would like to use them later, e.g.,

> [...]

The key here is that you need an index of the sort
quadrature_point(q_index).global_index()
in your code snippet. Such a function doesn't exist -- we don't usually
consider a global numbering of quadrature points, we only ever number them on
each cell.

But you can compute something like this:
const unsigned int global_q_index
= cell->active_cell_index() * n_q_points + q_index;
This just enumerates them one cell after the other. Your 'y' vector then must
have size
triangulation.n_active_cells() * quadrature.size();

[Yang Liu]

Thank you for your answer. This helps a lot. 

I noticed that the active cells are traversed following the "Z-order". And the support points are labeled in a similar way up to removing duplicate points. 

I wonder if there is a explicit "Z-order" mapping, from {cell->active_cell_index(), cell->level(), q_index} to global_q_index. 

How DEAL.II implements the ordering of the supporting points?

在 2020年4月29日星期三 UTC+8上午10:44:34，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 4/29/20 4:24 AM, Yang Liu wrote:
>
> I noticed that the active cells are traversed following the "Z-order". And the
> support points are labeled in a similar way up to removing duplicate points.
>
> I wonder if there is a explicit "Z-order" mapping, from
> {cell->active_cell_index(), cell->level(), q_index} to global_q_index.

You shouldn't assume this to be true. It depends on the mesh. What you see is
true for a square mesh, but if you read a mesh from a mesh generator, the
order of cells will be unpredictable and so will be the active_cell_index() of
cells.

> How DEAL.II implements the ordering of the supporting points?

We loop through all cells and enumerate all DoFs not yet enumerated. But that
is really an implementation detail you shouldn't rely upon.