[Dune] tranformed geometry across periodic boundaries

[matteo.semplice at uninsubria.it]

Martin,

    thanks a lot! Your patch makes face.geometryInOutside().center() work as expected in ALUCubeGrid<2,2>  as long as the intersection is conforming. (And it doesn't break ALUSimplex...)

Thanks a lot!
     Matteo


Il 21/05/2012 17:12, Martin Nolte ha scritto:
> Hi Matteo,
>
> hopefully I found the problem in the conforming case. If the intersection is
> non-conforming, the cache for the local geometries is still lacking.
>
> Could you test the attached patch?
>
> Best,
>
> Martin
>
> On 05/21/2012 02:49 PM, Martin Nolte wrote:
>> Hi Matteo,
>>
>> you are right: the trick does not work for ALUCubeGrid in 2d. The reason is
>> that 2d ALUGrid's local intersection geometries are currently computed from
>> the intersection geometry and the outside geometry (except for conforming
>> intersections in a trianglular grid). I'll have a further look into this.
>>
>> Best,
>>
>> Martin
>>
>> On 05/21/2012 10:44 AM, Matteo Semplice wrote:
>>> Thanks Martin,
>>>
>>> I take that the trick you meant is obtaining the displacement from the
>>> boundary as
>>>
>>>
>>> GlobalCoordinate displacement =
>>>
>>> outside.geometry().center()
>>>
>>> - outside.geometry().global(face.geometryInOutside().center());
>>>
>>> and then computing the outside center as
>>>
>>> GlobalCoordinate gl_out_cx_trick = face.geometry().center() + displacement;
>>>
>>> The problem is that is does not work as expected, because also
>>> face.geometryInOutside() seems to behave
>>> in an inconsistent way across different grid implementations.
>>>
>>> I attach a test.cc code and two dgf's for periodic sqare/cube grids that I
>>> used in my tests. I also
>>> attach the output of the program, when compiled and run with different grids
>>> (some comments are in there)
>>>
>>> The code sould be compiled by defining GRIDTYPE and GRIDDIM as in the
>>> grid-howto examples.
>>> Also note that since I was surprised by the YASP grid failure to create a
>>> periodic grid at all,
>>> I tried to rule out a mistake in the dgfparser, so can also compile the code
>>> with FORCE_YASP defined
>>> and it will create a periodic 2d YASP without using the dgfparser (the results
>>> are the same, though)
>>>
>>> Here is a summary of the results
>>>
>>> yes = works as in the official interface
>>>
>>> trick = not conformant to official interface, but your trick works
>>>
>>> no = neither the offical interface method nor the trick works
>>>
>>> no b.c. = grid does not handle the periodic b.c.
>>>
>>> crash = won't compile, or crash during execution
>>>
>>> 2D 3D
>>>
>>> YASP no b.c. no b.c.
>>>
>>> SGRID no b.c. no b.c.
>>>
>>> ALBERTA trick trick
>>>
>>> ALUCUBE no crash
>>>
>>> ALUSIMPLEX trick crash
>>>
>>> Best,
>>> Matteo
>>>
>>> On 16/05/2012 13:36, Martin Nolte wrote:
>>>> Hallo Matteo,
>>>>
>>>> the problem is that discussion about periodic boundary conditions got stuck
>>>> around 3 years ago.
>>>>
>>>> The documentation of dune-grid describes the "official" interface for
>>>> implementing perioditicy and is only implemented by YaspGrid.
>>>>
>>>> AlbertaGrid, ALUGrid and SPGrid support periodic boundary conditions in a
>>>> different (unofficial) way (as you already observed): They return the original
>>>> element on the other side. Since the intersection reports itselft as boundary,
>>>> you can now choose the appropriate boundary condition (e.g., periodic). Notice
>>>> that this is interpreted as a special choice of boundary condition and you
>>>> have to do something by hand, here.
>>>>
>>>> To obtain the transformation, AlbertaGrid's intersection implementation has a
>>>> method transformation (in case that helps you). Apart from that, you can use
>>>> the geometryInOutside and the outside element geometry to obtain the world
>>>> coordinates of the intersection on the outside element.
>>>>
>>>> The big problem is that only 1 developer (me) needed periodicity for a long
>>>> time and, hence, the necessary interface discussion is still pending. I hope
>>>> this clarifies the current situation, even if it does not really help.
>>>>
>>>> So, you will have to decide which technique you like better and either use the
>>>> existing grids implementing it or implement such a grid yourself.
>>>>
>>>> Best,
>>>>
>>>> Martin
>>>>
>>>> On 05/16/2012 12:09 PM, Matteo Semplice wrote:
>>>>> Hi dune,
>>>>> I think I need some explanation on the exact meaning of the phrase
>>>>> "transformed geometry" in the documentation of Dune::Intersection
>>>>> <http://www.dune-project.org/doc/doxygen/dune-grid-html/class_dune_1_1_intersection.html#_details>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> (By the way, the Detailed Description mentions a method "outer()" which is
>>>>> not
>>>>> listed in the interface, so I assume that it means "outside()"...)
>>>>>
>>>>> Say I have a 2d grid with 9 1x1 cells organized as follows
>>>>>
>>>>> 7,8,9
>>>>> 4,5,6
>>>>> 1,2,3
>>>>>
>>>>> so cell 1 has center at [.5,.5], cell 2 at [1.5,.5] and so on.
>>>>> Suppose also that the x direction have periodic boundary conditions, i.e. my
>>>>> dgf file has the section
>>>>>
>>>>> PERIODICFACETRANSFORMATION
>>>>> 1 0 , 0 1 + 3 0
>>>>> #
>>>>>
>>>>>  From the phrase "outer() cell has a periodically transformed geometry (so
>>>>> that one does not see a jump or something like that)" in the section on
>>>>> periodic boundaries, I would expect that
>>>>>
>>>>> if the intersection "face" is the left side of cell 1, and I do
>>>>>
>>>>> EntityPointer pOut = face.outside();
>>>>> Entity out = *pOut;
>>>>>
>>>>> then mapper(out) points me to data related to cell 3, but that
>>>>> out->geometry().center()
>>>>> returned [-0.5,0.5], so that effetively I would see the data of cell 3 as if
>>>>> they were at the left of cell 1 and a PDE code does not need to know about
>>>>> the
>>>>> boundary transformation.
>>>>>
>>>>> As this is not the case (out->geometry().center() returns [2.5,0.5]) ,
>>>>> what is
>>>>> the meaning of that phrase?
>>>>> And, given the current situation, does the grid implementation provide a
>>>>> method to transform the [2.5,.5] into [-.5,.5] ?
>>>>>
>>>>> Best regards,
>>>>> Matteo
>>>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> Dune mailing list
>>> Dune at dune-project.org
>>> http://lists.dune-project.org/mailman/listinfo/dune