coupled problems

[Pavel Kus]

Dear all,

I am working on development of Agros2D, which is a multi-purpose engineering
application. We are currently trying to replace our own hp-fem library by deal.ii.
This implies that our usage of deal is slightly different, since we are not
solving one particular problem, rather than that, we need potentially everything
together (coupled, transient problem, with hp-adaptivity, etc..)

1)
My current problem is the following. We would like to solve coupled problems, where
different physical fields are defined on different parts of the domain, usually with
some overlap. First, I would like to implement weak coupling, i.e., one field is solved
first and then used as a source for the second field. As in the figure, one of the field
may be defined on M1 only, while the other on both M1 and M2. Even though the coupling
is weak, I still need to have access to values of the first field when assembling the matrix
and rhs for the second (when looping over the mesh elements), so I guess I should follow the
suggestions of step 46 from the tutorial, use only one triangulation and use FE_Nothing on parts
of the domain, where the particular field is not considered. As I have mentioned, we are using hp-fem
approach already. So I was thinking, that I would enrich the collection of fe spaces to contain
the FE_Nothing space and select it on appropriate parts. But it somehow does not work.
This is what we had previously and what functions well:

    // Gauss quadrature and fe collection
    for (unsigned int degree = m_fieldInfo->value(FieldInfo::SpacePolynomialOrder).toInt(); degree <= DEALII_MAX_ORDER; degree++)
    {
        m_feCollection->push_back(dealii::FESystem<2>(dealii::FE_Q<2>(degree), fieldInfo->numberOfSolutions()));
        m_quadrature_formulas.push_back(dealii::QGauss<2>(degree +  QUADRATURE_ORDER_INCREASE));
        m_face_quadrature_formulas.push_back(dealii::QGauss<2-1>(degree + QUADRATURE_ORDER_INCREASE));
    }

Then, at the beginning, the lowest degree elements (when hp adaptivity is not switched on, active_fe_index is 0 through
the calculation) are selected and used for the calculation, which is than correct. But
when I do the following change:

    // first position of feCollection, quadrature_formulas and face_quadrature_formulas belongs to NONE space
    // this will be used for implementation of different meshes
    m_feCollection->push_back(dealii::FESystem<2>(dealii::FE_Nothing<2>(), fieldInfo->numberOfSolutions()));
    m_quadrature_formulas.push_back(dealii::QGauss<2>(1));
    m_face_quadrature_formulas.push_back(dealii::QGauss<2-1>(1));

    // Gauss quadrature and fe collection
    for (unsigned int degree = m_fieldInfo->value(FieldInfo::SpacePolynomialOrder).toInt(); degree <= DEALII_MAX_ORDER; degree++)
    {
        m_feCollection->push_back(dealii::FESystem<2>(dealii::FE_Q<2>(degree), fieldInfo->numberOfSolutions()));
        m_quadrature_formulas.push_back(dealii::QGauss<2>(degree +  QUADRATURE_ORDER_INCREASE));
        m_face_quadrature_formulas.push_back(dealii::QGauss<2-1>(degree + QUADRATURE_ORDER_INCREASE));
    }

    // find those elements, which are used for this field
    dealii::hp::DoFHandler<2>::active_cell_iterator cell = m_doFHandler->begin_active(), endc = m_doFHandler->end();
    for (unsigned int index = 0; cell != endc; ++cell, ++index)
    {
        cell->set_active_fe_index(1);
    }

I would guess, that those pieces of code should be equivalent, but now the solution is zero. Nothing else is changed. Do you
have any idea what I might be missing?

2)
Another question is, how can be interior faces distinguished. I need to be able to specify essential boundary conditions on lines
B1 and B2, which are on the boundary of one field, but are not boundary of the triangulation. I can do it by directly using ConstraintMatrix,
as it is suggested in step 46 tutorial, but I need to be able to distinguish them. When I try to use CellData::boundary_id directly, I get
an exception ExcInteriorLineCantBeBoundary, so I guess, that inside the triangulation, the only acceptable value is internal_face_boundary_id.
Is there another way how to specify user-defined marker? The approach followed in step 46 (looking on cells adjacent to faces and distinguishing
faces between two different materials) is not sufficient for me, because, as on the figure, both B1 and B2 are between M1 and M2, but there might
be different boundary conditions prescribed.

I am sorry for a very long question, any help would be highly appreciated.

Thank you,

Pavel Kůs

[Bruno Turcksin]

Pavel,

On Thursday, January 29, 2015 at 4:11:46 AM UTC-6, Pavel Kus wrote:

    // first position of feCollection, quadrature_formulas and face_quadrature_formulas belongs to NONE space
    // this will be used for implementation of different meshes
    m_feCollection->push_back(dealii::FESystem<2>(dealii::FE_Nothing<2>(), fieldInfo->numberOfSolutions()));
    m_quadrature_formulas.push_back(dealii::QGauss<2>(1));
    m_face_quadrature_formulas.push_back(dealii::QGauss<2-1>(1));

    // Gauss quadrature and fe collection
    for (unsigned int degree = m_fieldInfo->value(FieldInfo::SpacePolynomialOrder).toInt(); degree <= DEALII_MAX_ORDER; degree++)
    {
        m_feCollection->push_back(dealii::FESystem<2>(dealii::FE_Q<2>(degree), fieldInfo->numberOfSolutions()));
        m_quadrature_formulas.push_back(dealii::QGauss<2>(degree +  QUADRATURE_ORDER_INCREASE));
        m_face_quadrature_formulas.push_back(dealii::QGauss<2-1>(degree + QUADRATURE_ORDER_INCREASE));
    }

    // find those elements, which are used for this field
    dealii::hp::DoFHandler<2>::active_cell_iterator cell = m_doFHandler->begin_active(), endc = m_doFHandler->end();
    for (unsigned int index = 0; cell != endc; ++cell, ++index)
    {
        cell->set_active_fe_index(1);
    }

I would guess, that those pieces of code should be equivalent, but now the solution is zero. Nothing else is changed. Do you
have any idea what I might be missing?

This looks alright. Are you sure that you call m_doFHandler->distribute_dofs(*m_feCollection) after set_active_fe_index ? Can you send a small code where you have this problem so I can look at it ?

For your second question, take a look at https://dealii.org/8.2.1/doxygen/deal.II/DEALGlossary.html#GlossUserData You should be able to add an index to any face you like.

Best,

Bruno

 

[Wolfgang Bangerth]

Pavel,

> 1) [...]

> I would guess, that those pieces of code should be equivalent, but now the
> solution is zero. Nothing else is changed. Do you
> have any idea what I might be missing?

As Bruno already said, I think this looks alright to the degree this can be
said without looking at the code as a whole. What does dof_handler.n_dofs()
return?

> 2)
> Another question is, how can be interior faces distinguished. I need to be
> able to specify essential boundary conditions on lines
> B1 and B2, which are on the boundary of one field, but are not boundary of the
> triangulation. I can do it by directly using ConstraintMatrix,
> as it is suggested in step 46 tutorial, but I need to be able to distinguish
> them. When I try to use CellData::boundary_id directly, I get
> an exception ExcInteriorLineCantBeBoundary, so I guess, that inside the
> triangulation, the only acceptable value is internal_face_boundary_id.

Correct, at least for boundary_id.

> Is there another way how to specify user-defined marker? The approach followed
> in step 46 (looking on cells adjacent to faces and distinguishing
> faces between two different materials) is not sufficient for me, because, as
> on the figure, both B1 and B2 are between M1 and M2, but there might
> be different boundary conditions prescribed.

Bruno's suggested approach with the user flags or user indices is the way to
go. It has the drawback that it isn't inherited from the mother face to the
children upon refinement, but it shouldn't be very hard to rectify this by an
additional loop of the kind

for (cell=tria.begin(); ...) // loop over all cells, not just active ones
for (face=0; ...)
if (cell->face(f)->user_index() == something_special) // B1 and B2
if (cell->face(f)->has_children())
for (unsigned int c=0; c<cell->face(f)->n_children(); ++c)
cell->face(f)->child(c)->set_user_index(something_special);

right after mesh refinement.

Best
W.

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

[Pavel Kus]

Dear all,

thank you very much for your replies. The first issue seems to be solved, I did call m_doFHandler->distribute_dofs(

*m_feCollection), but than I used  cell->face(face)->get_dof_indices (local_face_dof_indices, INDEX) with INDEX == 0, which was wrong.

As for the boundary id issue, using set_user_index() would be great, but I am afraid that I cannot do that. I need to store information prior to triangulation creation, in the process of loading the information provided by a mesh generator (from our structures to deal.ii structures). Please see the following (simplified) piece of code:

......
.... (cell data creation)
.....
        // boundary markers
        dealii::SubCellData subcelldata;
        for (int edge_i = 0; edge_i < edgeList.count(); edge_i++)
        {
            if (edgeList[edge_i].marker == -1)
                continue;

            dealii::CellData<1> cell_data;
            cell_data.vertices[0] = edgeList[edge_i].node[0];
            cell_data.vertices[1] = edgeList[edge_i].node[1];
            cell_data.boundary_id = edgeList[edge_i].marker + 1;  // this cannot be used for internal edges

            subcelldata.boundary_lines.push_back(cell_data);
        }

        dealii::GridTools::delete_unused_vertices(vertices, cells, subcelldata);
        dealii::GridReordering<2>::invert_all_cells_of_negative_grid(vertices, cells);
        dealii::GridReordering<2>::reorder_cells(cells);

        triangulation->create_triangulation_compatibility(vertices, cells, subcelldata);

Before this call, I think I cannot iterate over the cells. But After it, cell data are reordered and the information stored in my edgeList[edge_i].marker cannot be used due to inconsistent numbering. Is there a way out? Again, I need to be able to identify some of the internal edges for integral evaluation or because they are boundaries of one of the fields that are defined only on a part of the mesh.

Thank you for any advise.

Best,

Pavel



Dne čtvrtek 29. ledna 2015 11:11:46 UTC+1 Pavel Kus napsal(a):

[Timo Heister]

You can always identify lines by the two cells around it. The cells
can be tracked using CellID or level/index. You can also reconstruct a
map for all the lines from the cell information after every refinement
cycle.

> --
> 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.



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

[Pavel Kus]

I see, so is it true that reorder_cells() reorders numbering of vertices and edges,
but it always keeps the order of cells as provided by the user in
std::vector<dealii::CellData<2> > cells ? If so, I could do everything I need...

Thank you,

Pavel


Dne sobota 14. února 2015 18:00:47 UTC+1 Timo Heister napsal(a):

[Timo Heister]

> I see, so is it true that reorder_cells() reorders numbering of vertices and
> edges,
> but it always keeps the order of cells as provided by the user in
> std::vector<dealii::CellData<2> > cells ? If so, I could do everything I
> need...

I think that is true. But what stops you from identifying the cells
using user indices or material ids?

[Pavel Kus]

I admit that reading the templated code related to triangulation is hard for me.
But if I understand it well, I am not able to set_user_indices before creation (and
possible reordering) of the mesh. If both cells and their faces were reordered, my prior
information would be useless. The only things that I know how to use is boundary_id
and material_id. However, boundary_id cannot be assigned freely to interior faces,
Identifying faces through materials of neighbouring elements is not enough and I cannot
abuse material_id to contain element index before reorder, since it is unsigned char only.

If, however, cells are not reordered, I can reconstruct all the information I need.

Thanks,

Pavel

Dne neděle 15. února 2015 17:35:30 UTC+1 Timo Heister napsal(a):

[Pavel Kus]

Hi,
I have some problems with the following code, which you suggested for propagating
the user index to refined elements.

   for (cell=tria.begin(); ...)   // loop over all cells, not just active ones
     for (face=0; ...)
       if (cell->face(f)->user_index() == something_special)  // B1 and B2
         if (cell->face(f)->has_children())
           for (unsigned int c=0; c<cell->face(f)->n_children(); ++c)
             cell->face(f)->child(c)->set_user_index(something_special);

right after mesh refinement.

But according to my experience, after calling tria.refine_global(), even the
local_indices of faces of original elements are set to 0. Thus line 3 of the
previous code is never true, even though just before the call to refine_global()
it has been. Is that possible?

Thanks,

Pavel

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

[Timo Heister]

It looks like user_index of faces doesn't persist through refinement,
yes. This seems to be an oversight. You can work around it by using
the cell's user_index (it does persist!).

> --
> 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.

[Wolfgang Bangerth]

On 02/16/2015 01:55 AM, Pavel Kus wrote:
> I admit that reading the templated code related to triangulation is hard
> for me. But if I understand it well, I am not able to set_user_indices
> before creation (and possible reordering) of the mesh. If both cells and
> their faces were reordered, my prior information would be useless. The only
> things that I know how to use is boundary_id and material_id. However,
> boundary_id cannot be assigned freely to interior faces, Identifying faces
> through materials of neighbouring elements is not enough and I cannot abuse
> material_id to contain element index before reorder, since it is unsigned
> char only.
>
> If, however, cells are not reordered, I can reconstruct all the information
> I need.

I believe that it is true that they come out in the same order as they go into
reordering. If not (or in any case), you can abuse the material_id of each
cell to not actually store a material id, but an index that is different for
each cell. You can then do something like

for (cell=...) {
int index = cell->material_id();
cell->set_subdomain_id (previously_filled_array[index]);
...get more information using this index...

cell->set_material_id (numbers::invalid_material_id); // reset to default
}

Best
W.


--

[Pavel Kus]

Dne sobota 21. února 2015 5:08:19 UTC+1 bangerth napsal(a):

On 02/16/2015 01:55 AM, Pavel Kus wrote:
> I admit that reading the templated code related to triangulation is hard
> for me. But if I understand it well, I am not able to set_user_indices
> before creation (and possible reordering) of the mesh. If both cells and
> their faces were reordered, my prior information would be useless. The only
> things that I know how to use is boundary_id and material_id. However,
> boundary_id cannot be assigned freely to interior faces, Identifying faces
> through materials of neighbouring elements is not enough and I cannot abuse
> material_id to contain element index before reorder, since it is unsigned
> char only.
>
> If, however, cells are not reordered, I can reconstruct all the information
> I need.

I believe that it is true that they come out in the same order as they go into
reordering. If not (or in any case), you can abuse the material_id of each
cell to not actually store a material id, but an index that is different for
each cell. You can then do something like

   for (cell=...) {
     int index = cell->material_id();
     cell->set_subdomain_id (previously_filled_array[index]);
     ...get more information using this index...

     cell->set_material_id (numbers::invalid_material_id); // reset to default
   }

Best
  W.

I was thinking about that, but it turned out that material_id is of type unsigned char, so (ab)using it as cell id is hardly an option. It seems, however, that the order of cells is not changed in the process of triangulation creation and thus this issue seems to be solved for me. Thank you and all others very much for all advise.

Best,

Pavel