[Dune] Boundary conditions, segments and indexing

[Firmbach, Max]

Hey Carsten,


ok, then i will use the insertionIndex() to be on the save side.

Right now, Load is a vector with fixed values so I guess that

should work fine.


So for a minimal example, let's consider the domain given in the

attached picture. This could be a slender beam or plate or something

similar. On one side it is fixed at a wall (tag 1 for that B.C.), on the other

end a loadvector, which could be a force is applied (tag 3 for that B.C.).

The other boundary faces are tagged with 0, meaning that no B.C. should

be applied there.


The domain is meshed with one hexaeder element and let's consider linear

shape functions. The mesh is loaded via the gmsh-reader of DUNE taking

the boundary-segments into consideration. Printing out the vector with the

boundary tags they are given as:


B.C. tag: 0
B.C. tag: 1
B.C. tag: 3
B.C. tag: 0
B.C. tag: 0
B.C. tag: 0


so i guess they are just stored in the way appearing in the .msh file (

consecutively ordered). Now I want to check how DUNE sees the boundary
segments and which face sees which tag.


for( const auto& element : elements(gridView)) {
    for( const auto& isect : intersections(gridView, element)) {
      if(isect.boundary()) {

        auto ref = referenceElement<double, dim>(element.type());

        cout << "B.C. seg. index:  " << isect.boundarySegmentIndex() << endl;

        for(int i=0; i<ref.size(isect.indexInInside(), 1, dim); i++) {
          cout << "Coordinate of " << i << " is: "
               << isect.inside().geometry().corner(ref.subEntity(isect.indexInInside(), 1, i, dim))
               << endl;
        }
   }
}


This generates the following output listing which index

the actual boundary-segment has and lists the coordinates

of the points that span the boundary face (hopefully they are

printed in the right way):


B.C. seg. index:  2
Coordinate of 0 is: 0 0 0
Coordinate of 1 is: 5 0 0
Coordinate of 2 is: 0 0.1 0
Coordinate of 3 is: 5 0.1 0
B.C. seg. index:  0
Coordinate of 0 is: 0 0 0
Coordinate of 1 is: 5 0 0
Coordinate of 2 is: 0 0 1
Coordinate of 3 is: 5 0 1
B.C. seg. index:  4
Coordinate of 0 is: 5 0 0
Coordinate of 1 is: 5 0.1 0
Coordinate of 2 is: 5 0 1
Coordinate of 3 is: 5 0.1 1
B.C. seg. index:  5
Coordinate of 0 is: 0 0.1 0
Coordinate of 1 is: 5 0.1 0
Coordinate of 2 is: 0 0.1 0
Coordinate of 3 is: 5 0.1 1
B.C. seg. index:  1
Coordinate of 0 is: 0 0 0
Coordinate of 1 is: 0 0.1 0
Coordinate of 2 is: 0 0 1
Coordinate of 3 is: 0 0.1 1
B.C. seg. index:  3
Coordinate of 0 is: 0 0 1
Coordinate of 1 is: 5 0 1
Coordinate of 2 is: 0 0.1 1
Coordinate of 3 is: 5 0.1 1

Checking with the domain from above, the insertion index chooses

the tag for the boundary faces. So boundary segment face

with tag 1 corresponds to insertion index 1 that fits with the domain,

looking at the corresponding coordinates.

Now, tag 3 should correspond to insertion index 2, but the coordinates

behind that boundary face are not right, they don't fit with the picture

of the domain. So the insertion index, the boundary tag and the corresponding

boundary face coordinates don't match in the right way.

This results in wrong B.C. for the overall simulation, due to the mismatch.


Coming back to the "not working", this is what i mean/meant with that.

In 2d, using the same code just with the different dimension, i didn't face such a problem.


This is quite a lengthy answer, but i hope it clarifies my problem.

I'm not able to attach the msh.file somehow, if the file is needed i will post it.


Thanks for the help.

Max



________________________________
Von: Carsten Gräser <graeser at mi.fu-berlin.de>
Gesendet: Montag, 24. August 2020 10:34:50
An: Firmbach, Max; dune at lists.dune-project.org
Betreff: Re: [Dune] Boundary conditions, segments and indexing

Hi Max,

Am 24.08.20 um 09:46 schrieb Firmbach, Max:
[...]
> For the simple one element, hexaeder case the tags are read into
> the boundaryEntity vector like [0, 0, 3, 0, 1, 0]. So on four faces
> there is just a 0, which means no B.C. is applied. For front and back
> either B.C. 1 or 3 should be applied. Looping over the boundary
> intersections my insertion index for the segments is something like
> [0, 1, 2, 4, 3, 5], which is identical to the boundarySegmentIndex()
> method (somewhere it's said that they don't need to be identical
> somehow, for this case they are).
in general they are not, because a grid may renumber them
during creation. The GridFactory::insertionIndex() method
exists, to keep track of such renumberings. However, UGGrid
does not do any renumbering, as you observed. So that's not
the problem here (but you may still adjust your code to not
depend on implementation details).

There's another potential issue in your code. But since
you've only send a code fragment, it's not clear if it
really is a problem:

> for( const auto& element : elements(gridView)) {
>     for( const auto& isect : intersections(gridView, element)) {
>
>       GeometryType elementGeometry = element.type();
>       using Factory = ReferenceElements<double, dim>;
>       const auto &ref = Factory::general(elementGeometry);
[...]
>             for(int i=0; i<ref.size(isect.indexInInside(), 1, dim); i++) {
>               int row = indexSet.subIndex(element, ref.subEntity(isect.indexInInside(), 1, i, dim), dim);
>              loadVector[row] = Load;
If Load is a fixed value, that's OK. If it's computed depending
on i, then this is probably a bug, because the subentity
orientation in the reference element and grid elements
does in general not coincide. To avoid ever using i
directly better use re.subEntities(isect.indexInInside(), 1, dim)
instead.


> As written before, for 2d meshes it works pretty well,
> but for 3d it's the complete opposite and I don't really
> understand why.
Can you be a little more precise and describe what
"it does not work" means? A minimal complete example
would also be helpful.

Best,
Carsten
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20200824/141617e6/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Domain.png
Type: image/png
Size: 5871 bytes
Desc: Domain.png
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20200824/141617e6/attachment.png>