How to use CellId class

[Zhidong Brian Zhang]

Hello, everyone,

In short, my question is how to get the cell_id properly by using CellId class ?

I want to assign each cell with a different value (its pesudo density in topology optimization). All the values are stored in a cell-based vector x, which is initialized as follows,

opt->x.reinit(opt->mpi_communicator,
                  opt->triangulation.n_global_active_cells(),
                  opt->triangulation.n_locally_owned_active_cells(),
                  false);

Now I want to access to x value in each cell. My basic idea is to get the cell_id and then use x(cell_id). In order to get cell_id, I found I should use CellId class, specifically, cell->id().to_string().

std::string cell_id_string;
cell_id_string = cell->id().to_string();

When I printed cell_id_string to check its value

std::cout << cell_id_string << std::endl;

I found the results are in this form,

0_3:000
0_3:200
0_3:003
0_3:006
0_3:406
0_3:606
0_3:206
0_3:007
0_3:407
0_3:607
0_3:207

I am not sure how to deal with these values.

I attempted to solve it by trying to convert it to unsigned int. I achieve this not elegantly as follows,

unsigned int cell_id = std::stoi(cell_id_string.substr(4, 1)) * 8 * 8
                                   + std::stoi(cell_id_string.substr(4, 1)) * 8
                                   + std::stoi(cell_id_string.substr(6, 1));

Incidentally, I can get the cell_id.  But the number of digits in  cell_id_string can be changing with different mesh, for example, when I use finer mesh, the output becomes,

0_5:00370
0_5:00371
0_5:00372
0_5:00373
0_5:00374
0_5:00375
0_5:00376
0_5:00377
0_5:00400
0_5:00401
0_5:00402

Thus, the previous method is not working anymore.

So may I ask for your help on this problem? How can I intepret the cell_id with decimal format properly in order to use it as indices?

Do I misunderstand the CellId class? To be honest, I am not sure how this class is working?

I am pretty new to deal.II. Please excuse me if I ask improper questions.

Looking forward to your suggestions!

Thank you very much for your help in advance!

Best regards,

Z. B. Zhang

[Wolfgang Bangerth]

Zhidong,

> In short, my question is how to get the cell_id properly by using CellId class ?
>
> I want to assign each cell with a different value (its pesudo density in
> topology optimization). All the values are stored in a cell-based vector x,
> which is initialized as follows,
>
> opt->x.reinit(opt->mpi_communicator,
>                   opt->triangulation.n_global_active_cells(),
>                   opt->triangulation.n_locally_owned_active_cells(),
>                   false);
>
> Now I want to access to x value in each cell. My basic idea is to get the
> cell_id and then use x(cell_id). In order to get cell_id, I found I should use
> CellId class, specifically, cell->id().to_string().

I think this is the wrong approach. The cell ids are used to globally identify
cells, but there is no easy way to translate between cell ids and an index
into vectors. What you want to do is create another finite element field that
is piecewise constant, i.e., using FE_DGQ(0). Such fields have exactly one
value per cell, and they can be accessed like any other field-based quantity.

Best
W.

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

[Zhidong Brian Zhang]

Hi, Dr. W. Bangerth,

Thank you so much for your timely help!

I will try it and come back if I have further questions.

Best regards,

Zhidong Brian Zhang

[Zhidong Brian Zhang]

Hi, Dr. W. Bangerth,

I have tried the your suggestion that is to use FE_DGQ(0) to build another finite element field.

It works even though there is till some small problems. However, I prefer to double check the implementation with you:

I built two fem field objects.

FESystem<dim> fe;   // for describing the vector-valued field (eg. displacement)
FE_DGQ<dim> fe_dg; // for describing the cell field (eg. pesudo density in topology optimization)

And initialize them in the initialization list of the class,

fe(FE_Q<dim>(1), dim),
fe_dg(0),

And I initialize the cell-based vector (the pesudo density in topology optimization) as follows,

opt->locally_owned_cells= opt->dof_handler_dg.locally_owned_dofs();

DoFTools::extract_locally_relevant_dofs(opt->dof_handler_dg,
                                            opt->locally_relevant_cells);

opt->x.reinit(opt->locally_owned_cells,
                  opt->locally_relevant_cells,
                  opt->mpi_communicator);

Is this process correct?

Moreover, I wonder how to access opt->x by the cell ("cell" represents the cell iterator)?

My attempting solution is to use the syntax of PETSc directly as follows,

PetscScalar *xp;
VecGetArray(opt->x, &xp);  // getting a local array of the global vector

But when I access it by

xp[cell->active_cell_index()]

I got error:

Abort: Memory balance (p4est)
Abort: ../../zzd_Software_Installation_Packages/p4est-2.2/sc/src/sc.c:705
Abort
Abort: Memory balance (p4est)
Abort: ../../zzd_Software_Installation_Packages/p4est-2.2/sc/src/sc.c:705
Abort
Abort: Memory balance (p4est)
Abort: ../../zzd_Software_Installation_Packages/p4est-2.2/sc/src/sc.c:705
Abort
Abort: Memory balance (p4est)
Abort: ../../zzd_Software_Installation_Packages/p4est-2.2/sc/src/sc.c:705
Abort

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 7914 RUNNING AT zzd-tnk-4755
=   EXIT CODE: 6
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Aborted (signal 6)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions

May I ask you to give me some hints on this problem?

Thank you very much in advance!

Best regards,

Zhidong Brian Zhang

On Saturday, October 19, 2019 at 10:56:47 AM UTC-4, Wolfgang Bangerth wrote:

[Konrad Simon]

Hi Zhidong,

I don't know what exactly you want to do but I found it quite useful that the CellId class has an order relation. That makes it usable in a std::map. Use the CellId (you can retrieve it through cell-->id() if cell is a cell accessor like the one you use to loop over all your cells) as a key then and connect to it an object that, for example, contains specific information about your cell(s).

Best,
Konrad  

[Wolfgang Bangerth]

On 10/19/19 3:01 PM, Zhidong Brian Zhang wrote:
>
> And I initialize the cell-based vector (the pesudo density in topology
> optimization) as follows,
>
> opt->locally_owned_cells= opt->dof_handler_dg.locally_owned_dofs();
> DoFTools::extract_locally_relevant_dofs(opt->dof_handler_dg,
>                                             opt->locally_relevant_cells);
> opt->x.reinit(opt->locally_owned_cells,
>                   opt->locally_relevant_cells,
>                   opt->mpi_communicator);
>
> Is this process correct?

Yes. Though it seems confusing to me that you call these variables
locally_*_cells, because they really are index sets for degrees of freedom,
not cells.

> Moreover, I wonder how to access opt->x by the cell ("cell" represents the
> cell iterator)?

You need to find out which DoF lives on the cell. So you can do
cell_dg->get_dof_indices(...)
and obtain an array of the local DoF indices (which for a DGQ(0) element has a
single element) that contains the global index of the DoF that lives on this cell.

> My attempting solution is to use the syntax of PETSc directly as follows,
>
> PetscScalar *xp;
> VecGetArray(opt->x, &xp);  // getting a local array of the global vector
>
> But when I access it by
>
> xp[cell->active_cell_index()]

Yes, this does not work. That's because you confuse the index of the cell and
the index of the DoF that lives on it. You want the latter.

[Zhidong Brian Zhang]

Thank you very much for your prompt reply, Konrad!

My confusion is the output of cell->id(), for example,

0_3:000
0_3:200
0_3:003
0_3:006
0_3:406
0_3:606
0_3:206
0_3:007
0_3:407
0_3:607

0_3:207.

What type I need to use as the key in std::map? And is the map parallel-distributed? Because my model requires that.

I have some preliminary results based on Dr. W. Bangerth's latest suggest.

Best regards,

Z. B. Zhang

[Zhidong Brian Zhang]

Hi, Dr. Bangerth,

Thank you very much for prompt reply!

> Yes. Though it seems confusing to me that you call these variables
> locally_*_cells, because they really are index sets for degrees of freedom,
> not cells.

OK, I misunderstood it. I thought the numbering of the dof is same as that of the cell since the degrees of freedom is 1 in FE_DGQ(0).

They are still different, right?

> You need to find out which DoF lives on the cell. So you can do
>    cell_dg->get_dof_indices(...)
> and obtain an array of the local DoF indices (which for a DGQ(0) element has a
> single element) that contains the global index of the DoF that lives on this cell.

OK, got it. I should get the index of dofs instead of that of cells.

> Yes, this does not work. That's because you confuse the index of the cell and
> the index of the DoF that lives on it. You want the latter.

Again now I understand they are different.

However, for my project, I need to access both

DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_dg;

at the same time.

If I use the following loop,

for (const auto &cell : opt->dof_handler.active_cell_iterators())

{

   ...

}

how can I access to dof_handler_dg?

My attempting solution is to set up multiple iterator, just as in the tutorial:

https://www.dealii.org/current/doxygen/deal.II/group__Iterators.html,

typename Triangulation<dim>::active_cell_iterator ti =
                    opt->triangulation.begin_active();
typename DoFHandler<dim>::active_cell_iterator di1 =
                    opt->dof_handler.begin_active();

typename DoFHandler<dim>::active_cell_iterator di2 =
                    opt->dof_handler_dg.begin_active();

Therefore, I set up the loop as follows,

while (cell != opt->triangulation.end())

{

    // do someting

   

    ++ti;

    ++di1;

    ++di2;

}

(Note: now I can access to the cell-based vector x successfully as follows,

di2->get_dof_indices(local_dof_indices_dg);
unsigned int local_dof_indices_dg_int = local_dof_indices_dg.at(0);

opt->x(local_dof_indices_dg_int);

)

However, I found this method is very time-consuming compared to "for (const auto &cell : opt->dof_handler.active_cell_iterators())", like 11s vs. 1 s.

My question now is: is there other method to access multiple DoFHandler or/and Triangulation?

Thank you a lot in advance!

Best regards,

Z. B. Zhang

[Wolfgang Bangerth]

On 10/20/19 5:42 PM, Zhidong Brian Zhang wrote:
> Thank you very much for your prompt reply, Konrad!
>
> My confusion is the output of cell->id(), for example,
> 0_3:000
> 0_3:200
> 0_3:003
> 0_3:006
> 0_3:406
> 0_3:606
> 0_3:206
> 0_3:007
> 0_3:407
> 0_3:607
> 0_3:207.

These strings are simply a representation of the id of a cell -- they're only
meant to be written to screen, but not interpreted by humans:
https://github.com/dealii/dealii/pull/8934/files

> What type I need to use as the key in std::map? And is the map
> parallel-distributed? Because my model requires that.

CellId can be used in a std::map without further ado because it has an
operator<. But if you use std::map, then it only stores what you put in -- it
has no concept of anything being distributed across processors.

[Wolfgang Bangerth]

> > Yes. Though it seems confusing to me that you call these variables
> > locally_*_cells, because they really are index sets for degrees of freedom,
> > not cells.
>
> OK, I misunderstood it. I thought the numbering of the dof is same as that of
> the cell since the degrees of freedom is 1 in FE_DGQ(0).
> They are still different, right?

They are *conceptually* different. Whether or not they are different *in
practice* is a different question, but you should never *assume* that they are
the same.

> > Yes, this does not work. That's because you confuse the index of the cell and
> > the index of the DoF that lives on it. You want the latter.
>
> Again now I understand they are different.
>
> However, for my project, I need to access both
>
> DoFHandler<dim> dof_handler;
> DoFHandler<dim> dof_handler_dg;
>
> at the same time.
>
> If I use the following loop,
>
> for (const auto &cell : opt->dof_handler.active_cell_iterators())
> {
>    ...
> }
>
> how can I access to dof_handler_dg?

You can have two iterators into the two DoFHandlers that run in lock step just
like you do:

> My attempting solution is to set up multiple iterator, just as in the tutorial:
> https://www.dealii.org/current/doxygen/deal.II/group__Iterators.html,
>
> typename Triangulation<dim>::active_cell_iterator ti =
>                     opt->triangulation.begin_active();
> typename DoFHandler<dim>::active_cell_iterator di1 =
>                     opt->dof_handler.begin_active();
> typename DoFHandler<dim>::active_cell_iterator di2 =
>                     opt->dof_handler_dg.begin_active();
>
> Therefore, I set up the loop as follows,
>
> while (cell != opt->triangulation.end())
> {
>     // do someting
>
>     ++ti;
>     ++di1;
>     ++di2;
> }
>
> (Note: now I can access to the cell-based vector x successfully as follows,
> di2->get_dof_indices(local_dof_indices_dg);
> unsigned int local_dof_indices_dg_int = local_dof_indices_dg.at(0);
> opt->x(local_dof_indices_dg_int);
> )
>
>
> However, I found this method is very time-consuming compared to "for (const
> auto &cell : opt->dof_handler.active_cell_iterators())", like 11s vs. 1 s.

You mean time consuming to program, or to execute? The loop itself is cheap;
if it's slow, you need to figure out which part of the body is slow.

> My question now is: is there other method to access multiple DoFHandler or/and
> Triangulation?

You could just do everything in one DoFHandler. That's what I usually do: Put
all solution variables into the same DoFHandler.

[Zhidong Brian Zhang]

Hi, Dr. W. Bangerth and everyone,

Thank you for your help, Dr. W. Bangerth!

I have figured the solution out. I summarize it as follows in case that it can help some one else as well.

Problem summary: In topology optimization, each cell has a scalar variable (pesudo density). So the density will form a vector, but the vector dimensions are different from those of the vector, say LA::MPI::Vector locally_relevant_solution (refer to Step-40). Because it only has one related value per cell. How to build up such "cell-based" vector?

Solution steps: (inspired by the suggestion from Dr. W. Bangerth)

1. Build two fem field objects.

FESystem<dim> fe;   // for describing the vector-valued field (eg. displacement)
FE_DGQ<dim> fe_dg; // for describing the cell field (eg. pesudo density in topology optimization)

And initialize them in the initialization list of the class,

fe(FE_Q<dim>(1), dim),
fe_dg(0),

And initialize the cell-based vector (the pesudo density in topology optimization) as follows,

opt->locally_owned_cells= opt->dof_handler_dg.locally_owned_dofs();  

// it is named with "cells" because I want to distinguish it from "opt->locally_owned_dofs" which is for the first FE field (fe(FE_Q<dim>(1), dim)).

DoFTools::extract_locally_relevant_dofs(opt->dof_handler_dg,
                                            opt->locally_relevant_cells);

opt->x.reinit(opt->locally_owned_cells,
                  opt->locally_relevant_cells,
                  opt->mpi_communicator);

2. I need to access both FE fields at the same time.

DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_dg;

If I use the default loop,

for (const auto &cell : opt->dof_handler.active_cell_iterators())

{

   ...

}

how can I access to the second FE field, which is dof_handler_dg?

My attempting solution is to set up multiple iterators, just as in the tutorial:

https://www.dealii.org/current/doxygen/deal.II/group__Iterators.html,

typename Triangulation<dim>::active_cell_iterator ti =
                    opt->triangulation.begin_active();
typename DoFHandler<dim>::active_cell_iterator di1 =
                    opt->dof_handler.begin_active();

typename DoFHandler<dim>::active_cell_iterator di2 =
                    opt->dof_handler_dg.begin_active();

Therefore, I set up the loop as follows,

while (cell != opt->triangulation.end())

{

    // do someting

   

    ++ti;

    ++di1;

    ++di2;

}

3. Now I can access to the cell-based vector x successfully as follows,

di2->get_dof_indices(local_dof_indices_dg);
unsigned int local_dof_indices_dg_int = local_dof_indices_dg.at(0);

opt->x(local_dof_indices_dg_int);

It was very time-consuming to execute because I put the vector access, say opt->x(local_dof_indices_dg_int), in the most inner loop. I don't need to access it for each quadrature point nor each degree of freedom. Because for each cell, the opt->x is fixed. So I put opt->x in the outside of those loops, and the computing time is normal now.

In sum, the parallel distributed density vector has been set up.

Thank you very much for Dr. W. Bangerth's help!

Best regards,

Z. B. Zhang

 

[Wolfgang Bangerth]

On 10/21/19 9:11 AM, Zhidong Brian Zhang wrote:
>
> It was very time-consuming to execute because I put the vector access,
> say opt->x(local_dof_indices_dg_int), in the most inner loop. I don't
> need to access it for each quadrature point nor each degree of freedom.
> Because for each cell, the opt->x is fixed. So I put opt->x in the
> outside of those loops, and the computing time is normal now.

Glad to hear. And yes, moving things out of inner loops is a good strategy!

[Konrad Simon]

Hi Zhidong,

On Monday, October 21, 2019 at 1:42:30 AM UTC+2, Zhidong Brian Zhang wrote:

Thank you very much for your prompt reply, Konrad!

My confusion is the output of cell->id(), for example,

0_3:000
0_3:200
0_3:003
0_3:006
0_3:406
0_3:606
0_3:206
0_3:007
0_3:407
0_3:607
0_3:207.

What type I need to use as the key in std::map? And is the map parallel-distributed? Because my model requires that.

You can use the CellId object as a key itself: 

    std::map<CellId, ClassWithCellInfo> cell_info_map;

If you need the map distributed across nodes according to the distribution of cells then you can initialize the map (on each node separately) through

        for (; cell!=endc; ++cell)

{

if (cell->is_locally_owned())

{

                        CellId current_cell_id(cell->id());

std::pair<typename std::map<CellId, ClassWithCellInfo>::iterator, bool > result;

result = cell_basis_map.insert(std::make_pair(cell->id(),

cell_info_map));

                         // ClassWithCellInfo must have a copy constructor

Assert(result.second,

ExcMessage ("Insertion of cell_info_map into std::map failed. "

"Problem with copy constructor?"));

                 } // end if

        } // end ++cell

Best,
Konrad

[Zhidong Brian Zhang]

Thanks a lot, Konrad!

I will try to compare your "map" method to the "FE_DGQ" method in terms of efficiency. And may come back to report later.

Best regards,

Z. B. Zhang