[Dune] Intersection area and corners

[Carsten Gräser]

Hi Marco,

Am 01.06.2016 um 16:26 schrieb Marco Cisternino:
> Hi Martin,
> 
> I turned on the dune assertion but I have no interruption using dune 2.3
> 
> I tried to build my code using dune 2.4.1 and dune-subgrid r424 but
> dunecontrol tells me that it cannot find a configuration file for
> sune-subgrid.
you did not mention that you're using dune-subgrid before.
There's three possibilities:

a)You're only marking complete levels for the subgrid.
  Then it will inherit the conformity of the host grid.
b)You're marking only parts of each level. Then you will
  in general have _non-conforming_ grids.
c)You're very carefull in marking parts of each grid
  level, such that the subgrid is still conforming.

Since c) is very unlikely, I assume that either a) or b)
is true. In case of a), please try to give a test case
without using dune-subgrid. In case of b) please double
check if your grid is conforming, having in mind that
adaptively refined subgrids are in general non-conforming
in contrast to the hostgrid.

Best,
Carsten


> 
> Is there a way to avoid this build error?
> 
> Thanks
> 
> 
> Bests,
> Marco
> ------------------------------------------------------------------------
> *Da:* Marco Cisternino <marco.cisternino at optimad.it>
> *Inviato:* mercoledì 1 giugno 2016 16.01
> *A:* Martin Nolte; Dune
> *Oggetto:* Re: [Dune] Intersection area and corners
>  
> 
> Hi Martin,
> 
> I have a restart procedure, but I try to give the simplest version of
> the code reproducing the error.
> 
> However, I'll try to rebuild dune with assertion in order to see if even
> with dune 2.3 I have the same assertion interruption.
> 
> I'll try to move to 2.4, but I need to say if the dune-subgrid di
> compatible with this version, because I need it.
> 
> Did you need to change a lot of things in my code to run with dune 2.4?
> 
> What do you mean with bugs in initialization order??
> 
> Thanks a lot for your help.
> 
> 
> Marco
> ------------------------------------------------------------------------
> *Da:* Martin Nolte <nolte at mathematik.uni-freiburg.de>
> *Inviato:* mercoledì 1 giugno 2016 12.29
> *A:* Dune
> *Oggetto:* Re: [Dune] Intersection area and corners
>  
> Hi Macro,
> 
> I ported you code to the master branch (which was not too painful). For me,
> the program results in an assertion (see output below), so I cannot
> reproduce
> your bug.
> 
> Unfortunately, your "small example" is really huge (> 4000 lines of code).
> This is quite a problem, because I have to understand (and trust) your
> example
> code, before I can start hunting bugs in GeometryGrid. At the moment, I do
> neither understand the code, nor do I trust it (my compiler spat out a
> ton of
> warnings, including bugs in initialization order). It also seems to me that
> the code can be further reduced.
> 
> Another problem is that the code runs quite some time in debug mode (to
> compute some parametrizations). Couldn't they be precompiled offline (on
> your
> machine) and be shipped as data?
> 
> To be honest, I am a bit at a loss how to help you. We're not supporting 2.3
> anymore (since 2.4 is released) and I cannot reproduce your bug in the
> master
> branch. Instead I get a curious assertion, which indicates you are
> dereferencing an end intersection iterator.
> 
> Best,
> 
> Martin
> 
> PS: This is the output I get:
> 
> Hello World! This is intersectionBug.
> I am rank 0 of 1 processes!
> Starting folder selection...
>  The folder exists. With write permissions.
>  The folder exists. With write permissions.
> End of folder selection...
> The number of element parameters is: 0
> The number of vertex parameters is: 0
> Has boundary parameters? 1
> 4018 vertices read.
> 0 is the size of vtxParams_
> 1920 elements read: 0 simplices, 0 pyramids, 0 prisms, 1920 cubes.
> 0 is the size of elParams_
> WARNING (ignored): Could not open file 'alugrid.cfg', using default
> values 0 <
> [balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.
> 
> You are using DUNE-ALUGrid, please don't forget to cite the paper:
> Alkaemper, Dedner, Kloefkorn, Nolte. The DUNE-ALUGrid Module, 2016.
> 
> Created parallel ALUGrid<3,3,cube,nonconforming> from macro grid file
> 'testgrid2'.
> 
> Computing parametrizations...
> End of parametrizations computation...
> intersectionBug:
> /home/nolte/numerics/master/dune-alugrid/dune/alugrid/3d/topology.hh:201:
> static int Dune::ElementTopologyMapping<type>::dune2aluFace(int) [with
> Dune::ALU3dGridElementType type = (Dune::ALU3dGridElementType)7u]: Assertion
> `index >= 0 && index < numFaces' failed.
> Aborted
> 
> 
> On 05/30/2016 07:14 PM, Marco Cisternino wrote:
>> Hi Martin,
>> 
>> I found where the error is produced:
>> 
>> 
>> it is exaclty in dune/grid/geometrygrid/cornerstorage.hh in the
>> IntersectionCoordVector method 
>> 
>> template< std::size_t size >
>>       void calculate ( array< Coordinate, size > (&corners) ) const
>> 
>> Here corners are computed by using
>> 
>>  corners[ i ] = elementGeometry_.global( hostLocalGeometry_.corner( i ) )
>> 
>> 
>> hostLocalGeometry_.corner( i ) gives the right coordinates in reference
>> element of corner i of the intersection, but the global() method gives the
>> wrong result
>> 
>> In global() method this statement is used to compute the return 
>> mapping_->global( local )
>> 
>> if the mapping is affine (and this is the case) the method starts from
>> physical position of 0-corner of the element and compute the position of the
>> current corner by
>> jacobianTransposed_.umtv( local, global )
>> 
>> this is where the error is made: local is right and input global too, but
>> for only one element, namely element 43 and for only 2 corner, namely
>> {0,1,1}, i.e. corner 6 ,and {1,1,1}, i.e. corner 7 the results is wrong
>> I would say wrong jacobian but it works fine for all the other 6 corner of the
>> element (and the algorithm to compute it works for all the other intersections
>> in the grid!!)
>> 
>> Let me remark that elementGeometry_ member in IntersectionCoordVector has the
>> right physical corners relatives to {0,1,1}, i.e. corner 6 ,and {1,1,1}, i.e.
>> corner 7. So, it is a little bit weird to recompute intersection corners
>> starting from corner 0, instead of using corners stored in elementGeometry_.
>> Could you explain me why you chose this way?
>> 
>> Thanks for any help!
>> Bests,
>> 
>> Marco
>> 
>> 
>> ------------------------------------------------------------------------------
>> *Da:* Marco Cisternino <marco.cisternino at optimad.it>
>> *Inviato:* venerdì 27 maggio 2016 19.13
>> *A:* Martin Nolte; Dune
>> *Oggetto:* Re: [Dune] Intersection area and corners
>>  
>> 
>> Hi Martin,
>> 
>> I prepared the smallest case I could do.
>> 
>> It is a bit slow because of the computation of the element parametrizations.
>> 
>> It will print at stdout the corners of the element (index 43) having the sick
>> intersection from the geometry grid and from the host one and the corners of
>> the sick intersection from both grids.
>> 
>> The intersection has index 3 in this element.
>> 
>> As you can see in the geometry grid case corners 0 and 1 of the intersection
>> are not at the same position of corners 6 and 7 of the element.
>> 
>> In the case of the host grid everything works fine.
>> 
>> I would no say the bug is in dune, but I cannot understand why the
>> intersection is not coherent with its element.
>> 
>> It is not evident from this code but the neighbour through this intersection
>> is not coherent too. Its index is 1025.
>> 
>> 
>> I remark that using the leaf mapper of the geometry grid crashes the run
>> giving segmantation fault, see line 142 of intersectionBug.cc.
>> All the test I made are serial.
>> 
>> If you need more feel free to ask.
>> Thanks for your help,
>> 
>> Marco
>> ------------------------------------------------------------------------------
>> *Da:* Marco Cisternino <marco.cisternino at optimad.it>
>> *Inviato:* giovedì 26 maggio 2016 14.12
>> *A:* Martin Nolte; Dune
>> *Oggetto:* Re: [Dune] Intersection area and corners
>>  
>> 
>> Hi Martin,
>> 
>> I'll try to set up a small example.
>> 
>> The weird thing is that I have 2k elements in my coarse grid and this problem
>> occurs for only 1 intersection. That's why I would like to see how the
>> intersection is built, in order to understand if my code wrongly changes that
>> intersection object.
>> 
>> Thanks.
>> 
>> 
>> Marco
>> 
>> 
>> 
>> ------------------------------------------------------------------------------
>> *Da:* Martin Nolte <nolte at mathematik.uni-freiburg.de>
>> *Inviato:* mercoledì 25 maggio 2016 17.47
>> *A:* Dune
>> *Oggetto:* Re: [Dune] Intersection area and corners
>>  
>> Hi Marco,
>> 
>> the intersection of GeometryGrid is based on the intersection of the host
>> grid. Then, you use the geometryInInside and the element geometry (of the
>> inside element) to build the global geometry of the intersection. So, in a
>> conforming situation, the corners of the intersection should be a subset of
>> the corners of the inside element.
>> 
>> Therefore, I consider your question a bug. Can you provide a small example?
>> 
>> Best,
>> 
>> Martin
>> 
>> On 05/25/2016 03:59 PM, Marco Cisternino wrote:
>>> Thank you Martin for your quick reply.
>>> 
>>> No, for the configuration I'm debugging I'm not dealing with non-conforming
>>> elements. That's why I cannot really understand why intersection and element
>>> are not coherent.
>>> 
>>> How the intersection is built? I could look in the constructor to see why its
>>> corners are not among the element corners.
>>> 
>>> Thanks again
>>> 
>>> 
>>> Marco
>>> 
>>> 
>>> 
>>> 
>>> ------------------------------------------------------------------------------
>>> *Da:* Martin Nolte <nolte at mathematik.uni-freiburg.de>
>>> *Inviato:* mercoledì 25 maggio 2016 11.00
>>> *A:* Dune
>>> *Oggetto:* Re: [Dune] Intersection area and corners
>>>  
>>> Hi Marco,
>>> 
>>> I am not sure whether I get your question correctly. You seem to be comparing
>>> the geometry of the inside element to the geometry of the intersection.
>>> 
>>> This is only reasonable, if the intersection is conforming. Otherwise I expect
>>> the corners to be different. Is it possible, that your intersection is
>>> non-conforming?
>>> 
>>> Please note the using GeometryGrid with non-conforming grid views is a bit
>>> problematic. Basically you are creating a first order Lagrange function over a
>>> non-conforming grid. If you don't constrain hanging nodes, the position of a
>>> hanging node might not be on the edge of the bigger element but slightly moved.
>>> 
>>> You might also want to have a look at the following merge request:
>>> 
>>> https://gitlab.dune-project.org/core/dune-grid/merge_requests/60
> <https://gitlab.dune-project.org/core/dune-grid/merge_requests/60>
> 	
> Bugfix/geometrygrid various (!60) · Merge Requests · Core Modules /
> dune-grid <https://gitlab.dune-project.org/core/dune-grid/merge_requests/60>
> gitlab.dune-project.org
> Grid Interface and Implementations
> 
> 
> 
>> <https://gitlab.dune-project.org/core/dune-grid/merge_requests/60>
>>        
>> Bugfix/geometrygrid various (!60) · Merge Requests · Core Modules / dune-grid
>> <https://gitlab.dune-project.org/core/dune-grid/merge_requests/60>
>> gitlab.dune-project.org
>> Grid Interface and Implementations
>> 
>> 
>> 
>>> <https://gitlab.dune-project.org/core/dune-grid/merge_requests/60>
>>>        
>>> Bugfix/geometrygrid various (!60) · Merge Requests · Core Modules / dune-grid
>>> <https://gitlab.dune-project.org/core/dune-grid/merge_requests/60>
>>> gitlab.dune-project.org
>>> Grid Interface and Implementations
>>> 
>>> 
>>> 
>>> 
>>> It deals with normals in non-conforming grids, which might also be buggy in
>>> DUNE 2.3.
>>> 
>>> Best,
>>> 
>>> Martin
>>> 
>>> On 05/25/2016 10:44 AM, Marco Cisternino wrote:
>>>> Good morning,
>>>> I have a problem with one intersection in my dune application.
>>>> The area I get is not what I expect.
>>>> So I tried to look inside the Intersection object using gdb
>>>> It is important to say that my grid is a GeometryGrid with a DiscreteFunction to parametrize the elements.
>>>> 
>>>> If I look at these two things:
>>>> 
>>>> - *(in.real.insideGeo_.mapping_)
>>>> - *(iGeo.realGeometry.mapping_)
>>>> 
>>>> where the in variable type is 
>>>> GeometryGrid::LeafGridView::IntersectionIterator::Intersection
>>>> and iGeo is 
>>>> GeometryGrid::LeafGridView::IntersectionIterator::Intersection::Geometry
>>>> 
>>>> I see different corners stored in. I mean two of the 4 corners given by *(iGeo.realGeometry.mapping_) are different from the two correspondant corners given by *(in.real.insideGeo_.mapping_). All the other corners are the same.
>>>> 
>>>> If I compute the intersection area using corners given by *(in.real.insideGeo_.mapping_), I get the right area. On the other hand, if I compute the intersection area using the corners given by *(iGeo.realGeometry.mapping_), the area is wrong.
>>>> 
>>>> I want to specify that the corners given by *(in.real.insideGeo_.mapping_) are the same I get using e.geometry().corners(i) where e is of type
>>>> GeometryGrid::LeafGridView::template Codim<0>::template Partition<Dune::InteriorBorder_Partition>::Iterator::Entity
>>>> 
>>>> Could anyone help me please? Am I doing something wrong?
>>>> Ask for more details if you need.
>>>> 
>>>> I'm using DUNE 2.3 with ALUGrid 1.52. I know it is old, but at the moment I cannot migrate to the new Dune releases.
>>>> 
>>>> Thanks a lot!
>>>> 
>>>> Marco

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune/attachments/20160606/f03cfe0b/attachment.sig>