How to efficiently add / remove cells (maybe with hp::DoFHandler)

[Mathias Anselmann]

Hello,

I'm thinking about a good way of solving my little problem and hope someone can give me some hints:

I want to simulate how a solid moves in a fluid with a fixed mesh using Cut-FEM. I implemented this using the hp-framework and set all the DoFs where just the solid is to FE_Nothing. You can see this in (attached) two pictures: overlay_0 shows a solid circle lying in a fluid. On no_overlay_0 the overlay is removed and you see the the underlying "solution" behind the cut and the 4 complete solid cells (the FE_Nothing cells).

This works well, but the idea has always been to move now the solid in the fluid. This works until the solid gets out of the starting FE_Nothing cells. What I do in a naive way every timestep is I check where the solid is and set the active FE index accordingly but this is not enough, and the solution gets distorted. You see this in overlay_1 (distorsion top left of the solid should not be there) and in no_overlay_1 you see that see that see just activated cells are not continuous (and they should be, since I use continuous elements).

So now I'm thinking of a proper and computationally cheap way to achieve my goal and have two ideas in mind:

1.) Using the hp-framework and FE_Nothing:

After setting the right FE_Index I have call distribute_dofs(), update my constraints and reinit all my matrices and vectors. I further have to somehow transfer the old solution to the new "grid". I guess these things make this option computationally expensive.

2.) Not using FE_Nothing:

Just using the fluid FE everywhere on the mesh and set the DoF where the solid is to "0" (by using constraints?). Since I use a relatively complex block matrix (in parallel) with some non-quadratic blocks I guess this probably means fiddling around with the corresponding matrix entries to make the resulting matrix not singular?

Maybe someone can give me a hint if there is another, better way to achieve my goal or which option I should use.

Thanks so much,

Mathias

[Wolfgang Bangerth]

Mathias,

> I want to simulate how a solid moves in a fluid with a fixed mesh using
> Cut-FEM. I implemented this using the hp-framework and set all the DoFs
> where just the solid is to FE_Nothing. You can see this in (attached)
> two pictures: overlay_0 shows a solid circle lying in a fluid. On
> no_overlay_0 the overlay is removed and you see the the underlying
> "solution" behind the cut and the 4 complete solid cells (the FE_Nothing
> cells).
>
> This works well, but the idea has always been to move now the solid in
> the fluid. This works until the solid gets out of the starting
> FE_Nothing cells. What I do in a naive way every timestep is I check
> where the solid is and set the active FE index accordingly but this is
> not enough, and the solution gets distorted. You see this in overlay_1
> (distorsion top left of the solid should not be there) and in
> no_overlay_1 you see that see that see just activated cells are not
> continuous (and they should be, since I use continuous elements).

I can see that this doesn't look right. But do you understand *why*
things go wrong? I can't tell from the pictures, but I've often found
that when trying to solve a problem for which I know that there is a bug
somewhere, then it's worthwhile investing the time to *understand* why
it is wrong, rather than trying to come up with completely different
ways of doing this.

This is because I may want to use the technique that doesn't work again
at a later time, and it's useful to actually understand how to make it
work. It's also possible that the bug is not in fact where one thinks it
is, and then using a different approach isn't actually going to solve
anything.

> So now I'm thinking of a proper and computationally cheap way to achieve
> my goal and have two ideas in mind:
>
> 1.) Using the hp-framework and FE_Nothing:
> After setting the right FE_Index I have call distribute_dofs(), update
> my constraints and reinit all my matrices and vectors. I further have to
> somehow transfer the old solution to the new "grid". I guess these
> things make this option computationally expensive.

No, I would venture the guess that the cost of assembling and solving
linear systems is still the dominant factor.

> 2.) Not using FE_Nothing:
> Just using the fluid FE everywhere on the mesh and set the DoF where the
> solid is to "0" (by using constraints?). Since I use a relatively
> complex block matrix (in parallel) with some non-quadratic blocks I
> guess this probably means fiddling around with the corresponding matrix
> entries to make the resulting matrix not singular?

You should be able to put these constraints into an AffineMatrix object
for the whole DoFHandler. This will ensure that all rows and columns of
the matrix are correct.

Best
W.


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

[Mathias Anselmann]

Dear Professor Bangerth,

thanks for your reply, I hope I can answer your questions:

I can see that this doesn't look right. But do you understand *why*
things go wrong? I can't tell from the pictures, but I've often found
that when trying to solve a problem for which I know that there is a bug
somewhere, then it's worthwhile investing the time to *understand* why
it is wrong, rather than trying to come up with completely different
ways of doing this.

Oh yes, I actually did not assume my current way of doing things to work.

After the static case was running, it was just a one liner (calling set_fe_indices() again in each timestep) which was a really naive way of doing things.

The problem at the moment is, that I just change the fe_index from pure FE_Nothing to my actual FE_System without updating my sparsity pattern, and re-initializing my matrices / vectors.

To be honest I didn't even think, that this code runs in debug mode but throws an assertion, but it didn't...

So the two option 1 and 2 are the things, that I see, which can "fix" my problem.

No, I would venture the guess that the cost of assembling and solving
linear systems is still the dominant factor.

That's definitely right. Nevertheless I guess transferring old solutions with no degrees of freedom on an ancient FE_Nothing cell to an updated version where this cell has degrees of freedom might become tricky.

The problem with the second way is that just not using FE_Nothing but instead my normal FE_System and setting all degrees of freedom to "0" on the former FE_Nothing cells will not help either, since I use continuous elements and don't want my interior ("0"-solid-solution) to interchange with my outer (flow)-solution on the "cut-cells". I guess one way to overcome this issue is to set all DoF on a "former FE_Nothing cell" to "0" besides the DoFs that also belong to a cut-cell neighbor?

Greetings,

Mathias

[Wolfgang Bangerth]

> Oh yes, I actually did not assume my current way of doing things to work.
> After the static case was running, it was just a one liner (calling
> set_fe_indices() again in each timestep) which was a really naive way of doing
> things.
> The problem at the moment is, that I just change the fe_index from pure
> FE_Nothing to my actual FE_System without updating my sparsity pattern, and
> re-initializing my matrices / vectors.
> To be honest I didn't even think, that this code runs in debug mode but throws
> an assertion, but it didn't...

Yes, this can't work. If you change the active_fe_index but not call
distribute_dofs(), you will have inconsistent internal data structures. This
won't work.

> No, I would venture the guess that the cost of assembling and solving
> linear systems is still the dominant factor.
>
> That's definitely right. Nevertheless I guess transferring old solutions with
> no degrees of freedom on an ancient FE_Nothing cell to an updated version
> where this cell has degrees of freedom might become tricky.

Well, you have to define how you want the solution that didn't exist on the
empty cell to be interpreted if you actually have something on that cell now.
I think that more than implementation, it's a question of what you *want* to
happen.

> The problem with the second way is that just not using FE_Nothing but instead
> my normal FE_System and setting all degrees of freedom to "0" on the former
> FE_Nothing cells will not help either, since I use continuous elements and
> don't want my interior ("0"-solid-solution) to interchange with my outer
> (flow)-solution on the "cut-cells". I guess one way to overcome this issue is
> to set all DoF on a "former FE_Nothing cell" to "0" besides the DoFs that also
> belong to a cut-cell neighbor?

That's exactly the question you have to answer first. What do you *want* to
happen on these cells? If you know what you want to happen, we can talk about
implementing it :-)

[Mathias Anselmann]

Dear Professor Bangerth,

thanks again for your reply, I will start answering your questions and give you more insight into what I try to do at the moment:

That's exactly the question you have to answer first. What do you *want* to
happen on these cells? If you know what you want to happen, we can talk about
implementing it :-)

In my current Cut FEM approach I'm interested in the solution on the "physical domain", in my case that's the solution on the (moving) flow domain including all the boundaries (also on the solid moving ball). Since I don't adjust the grid to the problem I have also a kind of "flow domain" inside the solid, which is nonphysically. I don't want to see any effect of these "nonphysical cells" on my "physical cells", so the natural approach to me was using FE_Nothing - and this is perfect for the static case and does the job very well (I don't have a fluid domain inside the solid and additionally I save some DoFs inside the solid, which is a nice side effect, but not crucial to me). Therefore I would like to also use this possibility in the dynamic case, so I would like to go for the 1.) approach in my first post.

Well, you have to define how you want the solution that didn't exist on the
empty cell to be interpreted if you actually have something on that cell now.
I think that more than implementation, it's a question of what you *want* to
happen.

 

That's indeed the crucial point that I have to answer. What I definitely need, are more DoFs, when the solid moves outside a FE_Nothing cell. The question is what the value of my old solution on these cells shall be. I guess there are several possibilities e.g. a "fancy one" like extrapolating the solution from the physical domain from the step before to the now activated cell. But this is not what I want to do in this first step now. I just want to set the solution on these cell to "0" (knowing that probably an extrapolation would be better as a first guess, but since it is just a starting value for my solver, I guess the "0"-approach is good enough).

This behavior is of course not implemented in deal.ii at the moment. What I did to reflect the above answer is that I improved my naive approach, which I described in my first post here, a little bit and do the following now:

1. initialize a parallel::distributed::SolutionTransfer object and add my old (ghosted) solution vectors via prepare_for_coarsening_and_refinement
   
2. check if we have a cell that get's "activated"(change from FE_Nothing  to my FE_System for my flow problem) or "deactivated" (other way around). If this is the case, I do:
1. update the fe_indices
2. call execute_coarsening_and_refinement() on the triangulation 
3. distribute the new DoFs
   
4. create a new sparsity pattern and reinitialize my system matrix with it
5. create a new solution vector with the new (locally owned) DoFs
   
6. interpolate the solution to this new solution vector using SolutionTransfer::interpolate()
   

This works again in the static case, also if I omit the check in step 2 and do the whole process even if it's not necessary. In the dynamic case step 6 fails:

SolutionTransfer::interpolate() leads to the SolutionTransfer::unpack_callback() function called by the tria, which calls set_dof_values_by_interpolation() on the cell.

This calls set_dof_values() and the following assertion fails:

  Assert(values.size() == this->get_dof_handler().n_dofs(),
         typename DoFCellAccessor::ExcVectorDoesNotMatch());

which is of course understandable. At the moment I'm trying to understand the internals of the unpack process and look for a nice way to implement what I described above - so any help or advice is greatly welcomed :)

Greetings,

Mathias

[Wolfgang Bangerth]

Mathias,

> thanks again for your reply, I will start answering your questions and
> give you more insight into what I try to do at the moment:
>
> That's exactly the question you have to answer first. What do you
> *want* to
> happen on these cells? If you know what you want to happen, we can
> talk about
> implementing it :-)
>
>
> In my current Cut FEM approach I'm interested in the solution on the
> "physical domain", in my case that's the solution on the (moving) flow
> domain including all the boundaries (also on the solid moving ball).
> Since I don't adjust the grid to the problem I have also a kind of "flow
> domain" inside the solid, which is nonphysically. I don't want to see
> any effect of these "nonphysical cells" on my "physical cells", so the
> natural approach to me was using FE_Nothing - and this is perfect for
> the static case and does the job very well (I don't have a fluid domain
> inside the solid and additionally I save some DoFs inside the solid,
> which is a nice side effect, but not crucial to me). Therefore I would
> like to also use this possibility in the dynamic case, so I would like
> to go for the 1.) approach in my first post.

What I meant to say is this: Say your domain moves and a cell that used
to be in the solid (=nonphysical) now becomes part of the fluid
(=physical). If you have a time-dependent problem (or a nonlinear
problem), your current solution will have degrees of freedom on this
cell, but the previous solution does not. But you need to compute a time
derivative, or a Newton update on this cell -- you will somehow have to
define what the *previous* solution should be on the cell in question
when there were no previous degrees of freedom there.

You answered that question here:

> That's indeed the crucial point that I have to answer. What I definitely
> need, are more DoFs, when the solid moves outside a FE_Nothing cell. The
> question is what the value of my old solution on these cells shall be. I
> guess there are several possibilities e.g. a "fancy one" like
> extrapolating the solution from the physical domain from the step before
> to the now activated cell. But this is not what I want to do in this
> first step now. I just want to set the solution on these cell to "0"
> (knowing that probably an extrapolation would be better as a first
> guess, but since it is just a starting value for my solver, I guess the
> "0"-approach is good enough).
>
> This behavior is of course not implemented in deal.ii at the moment.
> What I did to reflect the above answer is that I improved my naive
> approach, which I described in my first post here, a little bit and do
> the following now:
>

> 1. initialize a parallel::distributed::SolutionTransfer object and add

> my old (ghosted) solution vectors via
> prepare_for_coarsening_and_refinement

> 2. check if we have a cell that get's "activated"(change from

> FE_Nothing  to my FE_System for my flow problem) or "deactivated"
> (other way around). If this is the case, I do:

> 1. update the fe_indices
> 2. call execute_coarsening_and_refinement() on the triangulation
> 3. distribute the new DoFs
> 4. create a new sparsity pattern and reinitialize my system matrix
> with it
> 5. create a new solution vector with the new (locally owned) DoFs
> 6. interpolate the solution to this new solution vector using
> SolutionTransfer::interpolate()

That seems like a reasonable workflow.

> This works again in the static case, also if I omit the check in step 2
> and do the whole process even if it's not necessary. In the dynamic case
> step 6 fails:
>
> SolutionTransfer::interpolate() leads to the
> SolutionTransfer::unpack_callback() function called by the tria, which
> calls set_dof_values_by_interpolation() on the cell.
> This calls set_dof_values() and the following assertion fails:
>
>   Assert(values.size() == this->get_dof_handler().n_dofs(),
>          typename DoFCellAccessor::ExcVectorDoesNotMatch());
>
>
> which is of course understandable. At the moment I'm trying to
> understand the internals of the unpack process and look for a nice way
> to implement what I described above - so any help or advice is greatly
> welcomed :)

I don't understand how we arrive at this exact place, but I can think of
many scenarios where the class is simply not prepared to interpolate
from an FE_Q to an FE_Nothing or the other way around. Do you think you
can create an artificial testcase that shows exactly the kind of problem
you have? I.e., a situation where you create a mesh with just a few
cells (as few as possible), distribute DoFs on it, create a vector, and
then change the active_fe_index pattern before you call
SolutionTransfer? The test doesn't need to compute anything -- just
leave the vectors you have at zero, or initialize them with random data.
The point is simply to come up with as simple a testcase that shows the
assertion you mention above.

[Mathias Anselmann]

Dear Professor Bangerth,

thank you very much for your reply and please excuse my late answer. Unfortunately I didn't have the time the last couple of days to work on this problem or further investigate it.

I have edited the pack_callback() function in solution_transfer.cc a little bit myself before your reply and made further progress on this problem (and in return created some other problems...). I will do some more research the next couple of days on this topic and report back here.

Thanks for your help,

Mathias

[Mathias Anselmann]

Hello,

so finally I have made some progress on this topic and implemented a first very naive approach in this branch of deal.II.

I just want to recap what my goal is:

- I have a solid ball in a fluid, that moves to the right:

I use a non-fitted fluid background mesh which stays the same, independently of where the solid ball is.

On cells, where just the pure solid is, I define a pure FE_Nothing system - I don't want and don't need a solution there:

On the cells that are partly fluid and partly solid I do a cut-cell integration just over the fluid part - the "solution" in the solid part is never seen by my solver.

The ball moves to the right and in the next step the solution looks like this:

As you see a cell, which was FE_Nothing on the picture before got now activated and is part of the solution.

Here you see the full solution, without overlay:

So here are just 2 cells pure "FE_Nothing" cells.

In order to get to this solution, I use Newton's method and need a good starting value for my solver.

So I somehow want from my old solution:

a meaningful starting value for the solution above - and this is what I wanted to achieve as a starting value:

So I basically used the old solution, activated two "FE_Nothing cells" and set the values on the boundary to the old values of the previous solution.

In order to get there, I do the following, as mentioned above:

1. initialize a parallel::distributed::SolutionTransfer object and add my old (ghosted) solution vectors via prepare_for_coarsening_and_refinement
   

1. check if we have a cell that get's "activated"(change from FE_Nothing  to my FE_System for my flow problem) or "deactivated" (other way around). If this is the case, I do:

1. update the fe_indices and flag the cells that change from FE_Nothing to a real FE_System with the user_flag
   

1. call execute_coarsening_and_refinement() on the triangulation 

1. distribute the new DoFs
   

1. create a new sparsity pattern and reinitialize my system matrix with it

1. create a new solution vector with the new (locally owned) DoFs
   

1. interpolate the solution to this new solution vector using SolutionTransfer::interpolate()

If I do this with the current master branch of deal.II the assertion in "error_1.txt" is thrown, because the FE_Nothing cells don't create any values that are packed.

Therefore I changed some parts of solution_transfer.cc in order to overcome some issues: I just create 0-vectors of the right size, but just for the FE_Nothing cells, that are activated (user_flag is set) --> line 372 in solution_transfer.cc

In order to overcome the issue, that the library also tries to pack values on pure FE_Nothing cells I added the check in line 74, so cells with no DoFs pack actually nothing.

In order to get the unpack process running, I added the check in line 462, so that we just unpack data, if we are not on a cell that changes from FE_Nothing to a real FE_System (user flag is set) or ir we are on a pure FE_Nothing cell.

The second check here (cell->user_flag_set() == false) also omits unpacking anything on the cells that change, so these DoFs get initialized with 0 and afterwards the DoFs on the surrounding of these cells overwrite the values on the faces to the ones which I need. If I would omit this check and would start with this solution:

The resulting "interpolated" version would be:

which is not what I would like. You see here that the unpacking of the old FE_Nothing cells (where I packed just "0" vectors) results in a wrong interpolation. Just omitting the unpack process completely results in the needed behavior.

So that's where I'm at the moment.

The problem that I see is, that these changes to solution_transfer.cc are all nasty and I don't  like them. I ask for user flags and dofs_per_cell and that's all very specific to my problem - it's just awkward.

I don't know if anybody else needs such a functionality in deal.II and if it's worth thinking of how to implement it in the library directly but can anybody here give me some advice on how to achieve the things that I describe above in a nicer way? I don't know where to do such a thing in the library or if I should overwrite a deal.II class - so any hint or help here would really be very nice.