[Dune] Intersection area and corners

Wed Jun 8 14:43:59 CEST 2016

Hi Marco,

Am 08.06.2016 um 14:35 schrieb Marco Cisternino:
> Hi Carsten,
> I didn't mention subgrid because the problem I have with intersection corners does not involve subgrid module.
> The toy code I sent to reproduce the error does not even depend on subgrid.
> But to test my real code with dune 2.4.1, I need subgrid module dependence. That's why I'building with it!
> I hope you can help me in building my code.
> I slightly modified subgrid module to use it in parallel. Surely, it is not a good way to do it but it works. For this reason I would like to continue to use my r424 version of subgrid module.
unfortunately I can't give support for your locally modified version of
dune-subgrid. But I doubt that a 'slightly modified' version will work
correctly in parallel. I'd expect problems especially wrt the
intersections.

Best,
Carsten

> 
> Thanks,
> Marco
> ________________________________________
> Da: Carsten Gräser <graeser at mi.fu-berlin.de>
> Inviato: lunedì 6 giugno 2016 12.47.57
> A: Marco Cisternino; Dune
> Oggetto: Re: [Dune] Intersection area and corners
> 
> 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

-- 
Prof. Dr. Carsten Gräser
Freie Universität Berlin
Institut für Mathematik
Arnimallee 6
14195 Berlin, Germany
phone: +49 30 838 72637
fax  : +49 30 838 472637
email: graeser at mi.fu-berlin.de
URL  : http://page.mi.fu-berlin.de/graeser

-------------- 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/20160608/a9e3cd84/attachment.sig>