adaptive mesh refinement doubts

[RAJAT ARORA]

Hello everyone,

I am using deal.ii to solve a 3D solid mechanics problem.

My code uses Parallel::Distributed::Triangulation and Petsc Wrappers.

Until now, the code does not support AMR, however, I want to implement it in the code. For the last couple of days, I am

looking at the example codes (esp. 31) that talk about how to do adaptive mesh refinement and interpolate the solution to the new mesh.

I have read several previous posts in this group as well.

Before actually starting to code, I wanted to ask a few questions to make sure I understand the concept and should be able to debug when the need arises. (I think this is very closely related to serialize/deserialize)

The adaptive mesh refinement works as follows:

1. Once the solution is known on a coarse mesh, and mesh is refined locally in some regions, SolutionTransfer class

is used to get the solution interpolated on the new mesh.

2. In order to have gauss point history to be interpolated to the new mesh, I can make a new vector which is numbered as per the dif handler of discontinuous Galerkin. I can then store the gauss point history in this vector. This vector can then be interpolated to the new mesh and its value at the Gauss points will give the new point history.

3. Constraints related to hanging nodes have to be taken care.

This will actually transfer the solution and point history from the old mesh to the new mesh. 

I also have 2 questions now:

1. Apart from the transfer of these things (like explained above), is there anything else one needs to keep in mind while working on AMR? 

2. After repartition, there can be 3 cases that a) a new cell (previously belonging to the different processor) becomes a part as a result of repartition.

b)  there are newly refined cells with parent cell belonging to same processors (pure refinement) c)  there are newly refined cells with parent cell belonging to same processors (refinement + repartition). 

In case of no repartition, are material id and boundary conditions flags transferred to newly generated cells? In other words, as long as no motion of cells from one processor to other takes place, the data is inherited from parent cells. Is this right?

In case of repartition,  are the material id and boundary conditions transferred along with these cells (when they move)?

The link https://www.dealii.org/8.5.0/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html says that when cells are moved from one processor to other, this information is not carried. 

Thanks.

[Wolfgang Bangerth]

Rajat,

> 1. Once the solution is known on a coarse mesh, and mesh is refined
> locally in some regions, SolutionTransfer class
> is used to get the solution interpolated on the new mesh.

Yes.

> 2. In order to have gauss point history to be interpolated to the new
> mesh, I can make a new vector which is numbered as per the dif handler
> of discontinuous Galerkin. I can then store the gauss point history in
> this vector. This vector can then be interpolated to the new mesh and
> its value at the Gauss points will give the new point history.

Yes, that's the easiest way to do this.

You could also look at the classes declared in
include/deal.II/base/quadrature_point_data.h as they do similar things.

> 3. Constraints related to hanging nodes have to be taken care.

Yes. But since you use a discontinuous field for your Gauss point data,
there are no hanging node constraints for those fields.

> I also have 2 questions now:
> 1. Apart from the transfer of these things (like explained above), is
> there anything else one needs to keep in mind while working on AMR?

I can't think of anything, but step-31/32/15 will show you what we do there.


> 2. After repartition, there can be 3 cases that *a)* a new cell

> (previously belonging to the different processor) becomes a part as a
> result of repartition.

> *b)*  there are newly refined cells with parent cell belonging to same
> processors (pure refinement) *c)*  there are newly refined cells with

> parent cell belonging to same processors (refinement + repartition).

Yes.

> In case of no repartition, are material id and boundary conditions flags
> transferred to newly generated cells? In other words, as long as no
> motion of cells from one processor to other takes place, the data is
> inherited from parent cells. Is this right?

Yes.

> In case of repartition,  are the material id and boundary conditions
> transferred along with these cells (when they move)?

No.

> The link
> https://www.dealii.org/8.5.0/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html says
> that when cells are moved from one processor to other, this information
> is not carried.

Correct.
Best
W.


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

[RAJAT ARORA]

Hello Professor,

Thanks a lot for your earlier answers. I implemented all the details and it works like a charm.

I still have a few questions that I want to ask. The questions also have solutions/suggestions.

Sorry for being too verbose. Also, I am not sure if this should be a separate post so I am adding it to this thread.

1.  My problem involves the mesh movement in every time step. But with adaptive mesh refinement and automatic repartitioning, I found that the mesh becomes non-conforming, even when the displacements that I add to the mesh are conforming (I call hanging node constraints . distribute to make displacements conforming on the new mesh). I read this post https://groups.google.com/forum/#!searchin/dealii/mesh$20refinement$20move$20%7Csort:relevance/dealii/HGnGKX8Gbo0/6BtK6pg8vKcJ

It says that the instead of mesh movement with conforming displacements, store the coordinates of the mesh, then interpolate the mesh coordinates using solution transfer like other vectors, and then set the coordinates of the new mesh to this "transferred mesh coordinate vector". This works for me. The new mesh is conforming. But I dont understand why this works? I also tried the note https://www.dealii.org/8.5.0/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#a247598f1323a9f847832e60d6c840469 but even this did not work.

2. I want to apply boundary conditions on my domain. Lets say that the domain is rectangular and I am solving for displacements on this domain. I want to fix dofs that correspond to bottom left corner (x and y both fixed), and "x dof" on  bottom right corner  is fixed. Without AMR, I get the 3 dofs by looping over cells and even when the mesh  moves, the dofs are the same. But now, with AMR, the mesh moves as well the dofs are changed. 

So, I want to ask that is there any way I can know, the map from the vertex/dof in the old configuration to the vertex/dof in the new configuration.  Configuration here refers to mesh+dofs before and after refinement.

The earlier way of looping over cells and comparing coordinates to fix the correponding dofs will not work as the mesh moves with time. One way could be that I store the coordinates to be fixed and update that with time as well but I want to know if there is a map that can tell me that dof no. "45" in the old mesh is dof no "145" in the new mesh where the mapping is based on physical coordinates of the mesh.

3. This is most important question of all (for me :), pardon me if I am unable to explain this well). Uptil now I have seen and understood how to transfer smooth/continuous solutions from old mesh to new refined mesh. 

These solution vector just have the properties that either they are continuous through out or inside the elemen and the dofs value at the hanging nodes depends on these constraints only.

However, think of situation where the constraint is that the sum of all nodal dofs have to be constant plus there is only one value at the nodes for elements sharing the vertex (unlike discontinuous Galerkin). 

For example, take a single 1d element with continuous linear shape functions. This single element is now refined to prduce a total of 2 elements.

In such a case, we cannot say that the value at the hanging node (middle one) is just the average of the value at the end nodes of the single element to begin with. 

Now one may ask, where does such a situation arises? Think of that one element (2 after refining), as the boundary of the 2d domain/elements and a loading function is applied on that "1d surface". Now, when the equations are discretized, you get some discrete nodal traction dofs vector at the nodes -- the sum which is equal to the total force applied. Now, think that you donot know the loading function and this nodal vector with this property (sum = total force applied = constant) is given and the mesh is refined. In this scenario we cannot say that the value at the center node is the average of the two end nodes. 

Example:  in 1d with constant loading function such that total force = F, the nodal value of traction dof vector is F/2, F/2 on both nodes. Now, with the new mesh with 2 elements of equal length, this becomes, F/4 (left node), 2F/4 (center node), F/4 (right node).  However, if we simply interpolate the solution, we get F/2, F/2, F/2 on the 3 nodes, and this increases the total force on the system. With each refinement this gets worse.

(When I said that everything works like a charm, I am doing traction free case, so such a problem is not there for this case :) )

For my case, the loading is not given but comes from the solution of another differential equation. Now, I think that the interpolation of such a vector based on continuity is not well suited. (Please correct me if I am wrong).

I would like to know how such a vector can be interpolated after refinement. (Of course, 1 trivial solution is to never refine the cells on the boundary).

Thanks a lot for your attention.

[Wolfgang Bangerth]

Rajat,

> 1.  My problem involves the mesh movement in every time step. But with
> adaptive mesh refinement and automatic repartitioning, I found that the
> mesh becomes non-conforming, even when the displacements that I add to
> the mesh are conforming (I call hanging node constraints . distribute to
> make displacements conforming on the new mesh). I read this post
> https://groups.google.com/forum/#!searchin/dealii/mesh$20refinement$20move$20%7Csort:relevance/dealii/HGnGKX8Gbo0/6BtK6pg8vKcJ
> <ttps://groups.google.com/forum/#!searchin/dealii/mesh$20refinement$20move$20%7Csort:relevance/dealii/HGnGKX8Gbo0/6BtK6pg8vKcJ>
>
> It says that the instead of mesh movement with conforming displacements,
> store the coordinates of the mesh, then interpolate the mesh coordinates
> using solution transfer like other vectors, and then set the coordinates
> of the new mesh to this "transferred mesh coordinate vector". This works
> for me. The new mesh is conforming. But I dont understand why this
> works? I also tried the note
> https://www.dealii.org/8.5.0/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#a247598f1323a9f847832e60d6c840469 but
> even this did not work.

I don't know the answer to this question. You will need to construct a
small testcase with just a few cells that demonstrates the issue and
that people can look at.

> 2. I want to apply boundary conditions on my domain. Lets say that the
> domain is rectangular and I am solving for displacements on this domain.
> I want to fix dofs that correspond to bottom left corner (x and y both
> fixed), and "x dof" on  bottom right corner  is fixed. Without AMR, I
> get the 3 dofs by looping over cells and even when the mesh  moves, the
> dofs are the same. But now, with AMR, the mesh moves as well the dofs
> are changed.
> So, I want to ask that is there any way I can know, the map from the
> vertex/dof in the old configuration to the vertex/dof in the new
> configuration.  Configuration here refers to mesh+dofs before and after
> refinement.

You could store the cell->vertex_index() that corresponds to these
corner points at the beginning of the computation. Later, you would then
go over all cells, over all vertices of the cell, and compare the
vertex_index. The vertex_index is a piece of geometric information and
has nothing to do with the DoF index.

> 3. This is most important question of all (for me :), pardon me if I am

> unable to explain this well). [...]

> In such a case, we cannot say that the value at the hanging node (middle
> one) is just the average of the value at the end nodes of the
> single element to begin with.

The argument you are making above is the wrong one. You are thinking in
terms of nodal values and degrees of freedom, but you should be thinking
of functions of x,y that happen to be parameterized by nodal values but
at the end of the day are still just functions.

So, a constraint that says that the sum of the nodal values should be
zero does not make sense. What you want is that the *average* (or
integral) of the function is zero; this happens to be equivalent to the
sum of nodal values if you have a uniform mesh, but on adaptive meshes
the nodes are weighted differently (as are, by the way, nodes at the
boundary).

Similarly, if you want to think of the finite element function as
*continuous*, then the *only way of achieving this* is by requiring that
the hanging node in the middle is the average of the two adjacent ones.
If it does not satisfy this property, then the function is not continuous.

In your example, you are representing the *displacement* as a continuous
function, but the *traction* is related to the derivative of the
displacement and consequently is not, in general, continuous -- so it
should not have hanging nodes, and the 1/2 rule does not apply.