element birth using FE_Nothing and FE_Collection in time-dependent problem with adaptive mesh refinement

[luis miguel reig buades]

Hello all,

I am trying to add cells to the domain during a time-dependent calculation using mesh refinement. To do so, I made an FE collection and set some cells to FE_Nothing at the beggining and  then changed them back to FE_Q at a given time step right after executing triangulation.execute_coarsening_and_refinement();

However, I get the following error that I am not able to fix when I interpolate between meshes during the refinement:

The violated condition was:
    local_values_old[i] == number() || get_abs(local_values_old[i] - local_values[i]) <= get_abs(local_values_old[i] + local_values[i]) * 100000. * std::numeric_limits<typename numbers::NumberTraits< number>::real_type>::epsilon()
Additional information:
    Called set_dof_values_by_interpolation(), but the element to be set,
    value 0, does not match with the non-zero value 5.270338238394919e-05
    already set before.

This error does not happen if instead of activating, I deactivate cells during the simulation. Anyone has any suggestions?

I have attached my code.

Thank you for your time,

[Bruno Turcksin]

Hello,

This is not how you should activate elements. Instead **before**  you call triangulation.prepare_coarsening_and_refinement(), you call your activate_FE function but instead of using set_active_fe_index() you call set_future_fe_index().

Best,

Bruno

[luis miguel reig buades]

Thank you very much for your answer Bruno, this makes a lot of sense.

However, it does not seem to solve the interpolation problem, I am assuming the error happens because I did not assign an initial solution to the activated elements.

Is there any way to assign the activated elements a solution as I am looping through the cells and activating them? I have not been able to find this

Regards,

[Bruno Turcksin]

In case you are interested in using MPI, there is class that does that work for you: https://www.dealii.org/current/doxygen/deal.II/classparallel_1_1distributed_1_1experimental_1_1FieldTransfer.html and there are examples on how to use it (for example: https://github.com/dealii/dealii/blob/08160c07fbdc80241027691d7a0b424b95dc63d9/tests/hp/field_transfer_01.cc). If you don't want to use MPI, you can still look at how it is implemented and adapt the class.

Best,

Bruno

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/IwRWNUjnSYg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/d2cc656c-ed1d-4ba0-a28f-3401115999a0n%40googlegroups.com.

[luis miguel reig buades]

Thank you Bruno,

I will give it a try but mpi is still a bit scary for me. For the sake of learning, do you know what is causing my current error? Is it the lacking initial solution or something else?

The violated condition was:
    local_values_old[i] == number() || get_abs(local_values_old[i] - local_values[i]) <= get_abs(local_values_old[i] + local_values[i]) * 100000. * std::numeric_limits<typename numbers::NumberTraits< number>::real_type>::epsilon()
Additional information:
    Called set_dof_values_by_interpolation(), but the element to be set,
    value 0, does not match with the non-zero value 5.270338238394919e-05
    already set before.

[Wolfgang Bangerth]

On 1/18/24 03:28, luis miguel reig buades wrote:
>
> However, it does not seem to solve the interpolation problem, I am assuming
> the error happens because I did not assign an initial solution to the
> activated elements.
> Is there any way to assign the activated elements a solution as I am looping
> through the cells and activating them? I have not been able to find this

Are you saying that you continue to get the same error, even with the switched
order of operations? Can you show us the code you use, and the error you get?

Best
W.

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

[luis miguel reig buades]

Hello Wolfgang,

Yes, I am still getting the same error even using cell->set_future_fe_index(0)  before executing the refinement.

I have attached the error and all files needed to reproduce it, but anything that concerns us is in main.cc.   

the code is a mix between step-26 and step-27 to add elements during the computation.

To do so, I set some elements with FE_Nothing and the beggining and then, inside the method refine_mesh, I activate them at a given time_step using set_future_fe_index as suggested by Bruno before executing triangulation.execute_coarsening_and_refinement()

See:

    if (timestep_number == 25) activate_FEs();

    triangulation.execute_coarsening_and_refinement();

    setup_system();

    solution_trans.interpolate(previous_solution, solution);

when performing the interpolation with solution_trans.interpolate(previous_solution, solution); it gives me the error.

Thank you for the time,

[Marc Fehling]

Just a shot in the dark: Will it work if you set your FE_Nothing elements to be dominating?

https://www.dealii.org/developer/doxygen/deal.II/classFE__Nothing.html#a69940158d653b46222943f09c6b25a4c

Marc

[luis miguel reig buades]

Thank you very much for the suggestion Marc, but no luck :(

I activated domination like this:  fe_collection.push_back(FE_Nothing<dim>(1, true));

But got the same interpolation error:

An error occurred in line <155> of file </build/deal.ii-o0wOgt/deal.ii-9.5.1/source/dofs/dof_accessor_set.cc> in function
    void dealii::internal::set_dof_values(const dealii::DoFCellAccessor<dim, spacedim, level_dof_access>&, const dealii::Vector<number>&, OutputVector&, bool) [with int dim = 2; int spacedim = 2; bool lda = false; OutputVector = dealii::Vector<double>; number = double]

The violated condition was:
    local_values_old[i] == number() || get_abs(local_values_old[i] - local_values[i]) <= get_abs(local_values_old[i] + local_values[i]) * 100000. * std::numeric_limits<typename numbers::NumberTraits< number>::real_type>::epsilon()
Additional information:
    Called set_dof_values_by_interpolation(), but the element to be set,

    value 0, does not match with the non-zero value 0.001508036671432333
    already set before.

[luis miguel reig buades]

I have even tried to do add an initial solution like in this test: https://github.com/dealii/dealii/blob/e741ba5c864b5b273dead339480d9145badcadf5/tests/hp/solution_transfer_16.cc#L4

    triangulation.execute_coarsening_and_refinement();

    if (timestep_number == 10){

      activate_FEs();

    }

    setup_system();

    if (timestep_number == 10){

      solution.reinit(dof_handler.n_dofs());

      solution = 1.;

    }

    solution_trans.interpolate(previous_solution, solution);

but still the same error

[Magdalena]

Hi, 

we recently faced a similar problem and do not have a clean solution yet. What you could try is to use a parallel::distributed::Triangulation and the corresponding parallel::distributed::SolutionTransfer class (see https://www.dealii.org/developer/doxygen/deal.II/classparallel_1_1distributed_1_1SolutionTransfer.html#aebd38d82f946a96118f73e4e26defd61). Compared to the serial SolutionTransfer, it offers an additional constructor argument 

`const bool average_values=false` which you can set to true. Then the assert will hopefully not be triggered, because the conflicting values are averaged. Of course you have to decide yourself if the obtained values are still feasible for your computations.

Best,

Magdalena

[luis miguel reig buades]

That worked! And it was fairly easy to implement. Thank you very much Magdalena.

However, it does not seem to apply the initial solution I specify to the recently activated elements ( solution = 50.  ), it initializes them to solution = 0.

This is how I am doing it:

    parallel::distributed::SolutionTransfer< dim, Vector<double>> solution_trans(dof_handler, true);

    Vector<double> previous_solution;

    previous_solution = solution;

    triangulation.prepare_coarsening_and_refinement();

    solution_trans.prepare_for_coarsening_and_refinement(previous_solution);

    if (timestep_number == 25){

      activate_FEs(); // using set_future_FE_index()

    }

    triangulation.execute_coarsening_and_refinement();

    setup_system();

    if (timestep_number == 25){

      solution = 50.;

    }

solution_trans.interpolate(solution);

I guess I will have to set the initial solution in the activated elements myself with a function.

Regards,

[Magdalena]

In general, it would be nice to have a complete code example to better understand what you are doing. Based on what you have shown, I try my best to answer your question: Your code transfers `previous_solution` -- for which I am missing appropriate resizing according to DoFHandler layout (e.g. `Vector<double> previous_solution(dof_handler.n_dofs())`) and subsequent assignment of initial values -- to the new DoFHandler layout. The new values are stored in `solution` (described in https://dealii.org/developer/doxygen/deal.II/classparallel_1_1distributed_1_1SolutionTransfer.html#ae5abfde68fb2dc24e5c4690c8e5e247f). So in your case, setting `solution=50` has no effect because it will be overwritten within `solution_trans.interpolate(solution)`. In summary, I think you mixed up output and input vectors within `parallel::distributed::SolutionTransfer`. I hope this helps.

--
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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/IwRWNUjnSYg/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/0d4c1a42-050f-4614-af5a-c12b87cc9832n%40googlegroups.com.

--

Magdalena Schreter-Fleischhacker

[Wolfgang Bangerth]

On 1/24/24 03:41, luis miguel reig buades wrote:
>
> However, it does not seem to apply the initial solution I specify to the
> recently activated elements ( solution = 50.  ), it initializes them to
> solution = 0.

Are you running in debug mode? The statement
vector = value;
is only allowed if value==0, and vectors should be checking this in
debug mode.

Best
W.