CellDataStorage with mesh refinement

[Frederik S.]

Hello!

I've got a question on the usage of CellDataStorage for refined cells for a parallelized simulation. Below you'll find a sample of the code I use for refinement. 

I'm transferring the solution as stated in https://www.dealii.org/8.5.0/doxygen/deal.II/classparallel_1_1distributed_1_1SolutionTransfer.html and I sadly can't use ContinuousQuadratureDataTransfer, because my CellData is discontinous. Therefore I'm iterating through all cells, look at each cells parent and copy the CellData of the parent to the child cells (the blue part below; in reality the copying process is much longer, but then the code would be unnecessarily complex for this example). If I run the code with just one cpu core, everything works perfectly fine, so in this example every child cell has it's parent's CellData. But for simulations with more than one core, I just copy zeros into the child cells CellData (but I do not get an runtime-error).

Where did I make a mistake? Or is it just not possible to access parent cell's CellData for parallel meshes?

Thanks in advance!

--------------------------------------------------------

template <int dim>

void TEST<dim>::refine_grid (){

const unsigned int n_q_points    = quadrature_formula.size();

std::vector<const PETScWrappers::MPI::Vector *> solution_vectors (2);

solution_vectors[0] = &solution;

solution_vectors[1] = &old_solution;

parallel::distributed::SolutionTransfer<dim, PETScWrappers::MPI::Vector > soltrans(dof_handler);

typename parallel::distributed::Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end();

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

if (cell->subdomain_id() == this_mpi_process){

cell->set_refine_flag();

}

}

triangulation.prepare_coarsening_and_refinement();

soltrans.prepare_for_coarsening_and_refinement(solution_vectors);

triangulation.execute_coarsening_and_refinement();

setup_system (0);

typename parallel::distributed::Triangulation<dim>::active_cell_iterator cell2 = triangulation.begin_active(), endc2 = triangulation.end();

for (; cell2!=endc2; ++cell2){

point_history_accessor.initialize(cell2, n_q_points);

}

int highest_level_after_refinement=triangulation.n_global_levels()-1;

PETScWrappers::MPI::Vector interpolated_solution(system_rhs);

PETScWrappers::MPI::Vector interpolated_old_solution(system_rhs);

std::vector<PETScWrappers::MPI::Vector *> tmp (2);

tmp[0] = &(interpolated_solution);

tmp[1] = &(interpolated_old_solution);

soltrans.interpolate(tmp);

cell2 = triangulation.begin_active(), endc2 = triangulation.end();

for (; cell2!=endc2; ++cell2)

if(cell2->level()==highest_level_after_refinement){

typename parallel::distributed::Triangulation<dim>::cell_iterator cell_p=cell2->parent();

const std::vector<std::shared_ptr<PointHistory<dim> > > point_history_p = point_history_accessor.get_data(cell_p);

for (unsigned int child=0; child<8; child++){

unsigned int q_point=0;

std::vector<std::shared_ptr<PointHistory<dim> > > point_history = point_history_accessor.get_data(cell_p->child(child));

for (unsigned int q_point=0;q_point<n_q_points;q_point++){

point_history[q_point]->status=point_history_p[q_point]->status;

}

}

}

constraints.clear ();

constraints.reinit (locally_relevant_dofs);

DoFTools::make_hanging_node_constraints (dof_handler, constraints);

constraints.close ();

constraints.distribute(interpolated_solution);

constraints.distribute(interpolated_old_solution);

solution     = interpolated_solution;

old_solution = interpolated_old_solution;

dof_handler.distribute_dofs (fe);

}

[Timo Heister]

Frederik,

without looking at your code (sorry!), I think you can use
Triangulation::register_data_attach and notify_ready_to_unpack to work
with any kind of data that you want to transfer around (this is the
functionality that is used by SolutionTransfer and the
ContinuousQuadratureDataTransfer.

See https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_blob_master_tests_mpi_attach-5Fdata-5F01.cc&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4zqq8LzalJka1USeMsXKUwuzA9ZYVYDocY6Qd-sHk98&s=nC5wgqjDzcurV-G_a89DWCjEZbxP9MHaiIWP1yIKd4w&e=
for an example.

On Thu, Oct 26, 2017 at 12:50 PM, 'Frederik S.' via deal.II User Group
<dea...@googlegroups.com> wrote:
> Hello!
>
> I've got a question on the usage of CellDataStorage for refined cells for a
> parallelized simulation. Below you'll find a sample of the code I use for
> refinement.
>
> I'm transferring the solution as stated in

> https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.5.0_doxygen_deal.II_classparallel-5F1-5F1distributed-5F1-5F1SolutionTransfer.html&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4zqq8LzalJka1USeMsXKUwuzA9ZYVYDocY6Qd-sHk98&s=Nc9Mb-wVDF_hg0oMht1A43zN_F7kpOv6vDewaKyyvzI&e=

> --
> The deal.II project is located at https://urldefense.proofpoint.com/v2/url?u=http-3A__www.dealii.org_&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4zqq8LzalJka1USeMsXKUwuzA9ZYVYDocY6Qd-sHk98&s=vf2x5H8OxETNRY5Xkg_2UlqLPYffjI7fsx3-Ng2CfHo&e=
> For mailing list/forum options, see
> https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_forum_dealii-3Fhl-3Den&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4zqq8LzalJka1USeMsXKUwuzA9ZYVYDocY6Qd-sHk98&s=na9F3U-UB7JStYIGEje_0_pNnfNGJoooknwYVOSoG6c&e=
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
> For more options, visit https://urldefense.proofpoint.com/v2/url?u=https-3A__groups.google.com_d_optout&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4zqq8LzalJka1USeMsXKUwuzA9ZYVYDocY6Qd-sHk98&s=uYqlSNAOmeNTEXlF1E5WamZiNmM6rO8HIgUxzSTtkn4&e= .



--
Timo Heister
http://www.math.clemson.edu/~heister/

[Jean-Paul Pelteret]

Hi Frederik,

I’m being wildly presumptuous here (so feel free to say if I’ve presumed incorrectly), but unless you have a special set of data that needs to be considered then the ContinuousQuadratureDataTransfer class may still be relevant to your use case. There is no requirement that your global FE space be continuous to use this class. Whats relevant is this snippet from the documentation:

"the data located at data_quadrature of each cell is L2-projected to the continuous space defined by a single FiniteElement projection_fe “

What happens is that, for each cell, the specified QP data is projected to the finite element space specified by the input projection_fe, and then interpolated onto the new cell (or children cells in the case of refinement) of the refined triangulation. So I think that it might be possible that you've misinterpreted is what exactly the “Continuous” part of the class name refers to: It is assumed that for each individual cell the data is continuous within the cell (i.e. can be directly interpolated using the input FE), while the data between cells may not be. An example case where data is truly discontinuous within a cell (and therefore where this class is not applicable) would be the internal variables for yield or damage in plasticity. Would you care to elaborate more on what you’re using this for?

Best regards,

Jean-Paul

P.S. You really want to use the GeometryInfo class to get this value:

for (unsigned int child=0; child<8; child++){

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en

---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Frederik S.]

Hi Jean-Paul,

Thank you for your answer!

I have a truly discontinuous case, the case you mentioned with internal yield variables. Do you think it will be easier to use the method Timo mentioned "from scratch" or will it be easier to copy the source files of ContinuousQuadratureDataTransfer class and alter them to DiscontinuousQuadratureDataTransfer with some kind of nearest neighbor Interpolation?

Probably  it will still be reasonable to use ContinuousQuadratureDataTransfer for the continuous data I have. I yet don't absolutely understand which order to use for the mass_quadrature element? data_quadrature is the quadrature rule I'm using in my model and projection_fe has to be of an order high enough to interpolate this data correctly, am I at least right there?

Best regards,

Frederik

PS: Thanks, I forgot GeometryInfo at this point

[Frederik S.]

Hi Timo,

Thank you for your answer, this sounds really good!

I tried to use the class, but I still have some Problems. In the example you are just transferring an integer variable. Is it also possible to transfer a whole CellDataStorage element (or at least some packed vector with ist elements) from a parent cell to each child? Do you also have an example how this works?

Thanks in advance,

Frederik

[Timo Heister]

> I tried to use the class, but I still have some Problems. In the example you
> are just transferring an integer variable. Is it also possible to transfer a
> whole CellDataStorage element (or at least some packed vector with ist
> elements) from a parent cell to each child? Do you also have an example how
> this works?

You can transfer whatever you want and the "status" argument will tell
you if you have to convert your data to children or parents.

Examples? It is used in SolutionTransfer and ContinuousQuadratureDataTransfer:
https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_blob_60bfd3170a465002ac1404dc36a5ba9671254b71_source_distributed_solution-5Ftransfer.cc-23L84&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4DhDywRdwgpQcqLQhc9e3vXnNh9twLFXh6cAv-rYOtU&s=HDE5UAeI6pVJZbKAVdMR4mt-Hb-RJxrDcUsUjCX4myo&e=
https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_dealii_dealii_blob_a0073ac1c0e88a30702a52f6522cbfffe1834ed1_include_deal.II_base_quadrature-5Fpoint-5Fdata.h-23L735&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4DhDywRdwgpQcqLQhc9e3vXnNh9twLFXh6cAv-rYOtU&s=j5xv1X624c05TmUCkWaj5qBdYefoSr6JBbEZ1-PIhcU&e=

We also use it in ASPECT to transfer particles:
https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_geodynamics_aspect_blob_0eec247edd99cd782a8fc66b36b140ce7a4105ae_source_particle_particle-5Fhandler.cc-23L858&d=DwIBaQ&c=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4&r=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw&m=4DhDywRdwgpQcqLQhc9e3vXnNh9twLFXh6cAv-rYOtU&s=OmRYoVuKEIL1dlaaCOgQLqNoCmHaHzJulO4dsMuAq2o&e=

[Jean-Paul Pelteret]

Dear Frederik,

Sorry for taking a little time to get back to you. What you’ve mentioned in terms of a DiscontinuousQuadratureDataTransfer is exactly the concept that we laid out in this issue when we first implemented this class (see the last point in the check list).

https://github.com/dealii/dealii/issues/2555

In my opinion, inheriting a value from the closest QP (as if there are patches on a cell where the internal variable values are constant) is good first attempt to deal with this difficult situation. If you’re willing to try to implement it then I’d be happy to provide assistance if you need it.

Best regards,

Jean-Paul

[Frederik S.]

Hello Timo!

Sorry it took so long, I was on a conference this week and didn't come to answer.

Many thanks for the example! I read through all the documentation pages, but it didn't come to my mind to just check the source code of this classes... Thanks again, this will help me a lot!

PS: see also my answer to Jean-Pauls message.

[Frederik S.]

Hey Jean-Paul!

Sorry it took so long to answer, I was on a conference this week and didn't come to write you.

Thank you for offering your help! I will definitely try to implement this and I'll try my best, to have it as close to the current layout of ContinuousQuadratureDataTransfer as possible. I'm thinking of giving the prepare_for_coarsening_and_refinement-function some kind of switch (only continuous data - only discontinuous data - both types of data), implementing the corresponding discontinuous interpolation and adding a pack_discontiuous_values/unpack_discontiuous_values-function to TransferableQuadraturePointData, to see which values are continuous (-> pack_values) and which are discontinuous (-> pack_discontiuous_values). I think two seperate pack functions would be easier to handle and more straightforward than one pack-function with some masking-variable. 

What do you think, does this sound like a good concept?