[Dune] jacobian construction
Martin Nolte
nolte at mathematik.uni-freiburg.de
Thu Jun 23 16:21:57 CEST 2016
Hi Marco,
I guess the correct solution is a relative error. To be most general, we could
migrate the entire comparison for zero onto the Traits class.
Now you know the exact problem, would you be so kind as to provide a really
minimal test case (based only on dune-geometry) and create an issue. We can
then discuss the best solution and put test and resolution into the master
branch. Personally, I would even consider a backport to a 2.4.2 release.
Best,
Martin
On 06/23/2016 03:54 PM, Marco Cisternino wrote:
> Hi Martin,
> that was exactly the point!!!
> Forcing the tolerance to 1.0e-20, it works.
> We also tried to change the check doing it component-by-component like this
>
> for( int i = 0; i < dim-1; ++i ){
> for( int j = 0; j < 3; ++j ){
> if(fabs(jtTop[ i ][j] - jt[ i ][j]) > Traits::tolerance()){
> return false;
> }
> }
> }
>
> Do you think is good enough?
> Thanks!
>
> Marco
> ________________________________________
> Da: Martin Nolte <nolte at mathematik.uni-freiburg.de>
> Inviato: giovedì 23 giugno 2016 08.42
> A: dune at dune-project.org
> Oggetto: Re: [Dune] jacobian construction
>
> Hi Marco,
>
> am I blind or are we talking a difference in the 7th decimal? If so, could this
> be a simple round-off error?
>
> If you suspect the geometry, you might replace the CachedMultiLinearGeometry in
> dune/grid/geometrygrid/geometry.hh by a MultiLinearGeometry.
>
> Moreover, I see that your coordinates are very close to each other. If you take
> a look at dune/geometry/multilineargeometry.hh (line 859 in master), you see the
> affine property is only accurate up to an absolute tolerance. This tolerance is
> defined by the Traits class of the MultiLinearGeometry (coming from
> GeometryGrid). You can try to change this value too.
>
> I hope these hints help a little.
>
> Best,
>
> Martin
>
> On 06/22/2016 06:58 PM, Marco Cisternino wrote:
>> Hi Martin,
>> I understand your explanation but I'm far from understanding how to debug my problem.
>> I accelerate the execution of the "small" code I sent you some days ago, but because of parametrization I cannot make it smaller.
>> I tried to better understand the problem, isolating one intersection of my conforming grid.
>> I remark that the grid is a coarse grid, it is 3D ALUCubeGrid and a GeometryGrid is built introducing a Q3 parametrization.
>> I'm still working with dune 2.3.0 and ALUGrid 1.52b. I'm porting the code to dune 2.4.1, encountering some problem.
>>
>> The problematic intersection is the one between element 42 and 43, its indices are 4 and 5 respectively.
>> Element 42 and 43 store their corner in the right way, not only they are what I expect but there is coherence between the 2 elements
>>
>> Corners element 42
>> gcorner 0 2.0000000000000000 0.2598190000000000 0.0000000000000000
>> gcorner 1 2.0000000000000000 0.2598189999999993 0.2000000000000000
>> gcorner 2 1.9806379930000000 0.2548077446999999 0.0000000000000000
>> gcorner 3 1.9806379930000000 0.2548077446999999 0.2000000000000000
>> gcorner 4 2.0000000000000000 0.2588190000000001 0.0000000000000000
>> gcorner 5 2.0000000000000000 0.2588190000000007 0.2000000000000000
>> gcorner 6 1.9806379930000000 0.2538077447000000 0.0000000000000000
>> gcorner 7 1.9806379930000000 0.2538077446999998 0.2000000000000000
>>
>> Corners element 43
>> gcorner 0 2.0000000000000000 0.2608364441000000 -0.0000000000000000
>> gcorner 1 2.0000000000000000 0.2608364440999995 0.2000000000000000
>> gcorner 2 1.9806379930000000 0.2558252375000000 0.0000000000000000
>> gcorner 3 1.9806379930000000 0.2558252374999995 0.2000000000000000
>> gcorner 4 2.0000000000000000 0.2598190000000000 0.0000000000000000
>> gcorner 5 2.0000000000000000 0.2598189999999993 0.2000000000000000
>> *gcorner 6 1.9806379930000000 0.2548077446999999 0.0000000000000000
>> *gcorner 7 1.9806379930000000 0.2548077446999999 0.2000000000000000
>>
>> Paying attention to them we can see correspondence between nodes 0,1,2,3 of element 42 and nodes 4,5,6,7 of element 43
>>
>> Therefore I asked intersection corners to both the elements, expecting exactly the same corners. But
>>
>> Corners of intersection 4 of element 42
>> gcorner 0 2.0000000000000000 0.2598190000000000 0.0000000000000000
>> gcorner 1 1.9806379930000000 0.2548077446999999 0.0000000000000000
>> gcorner 2 2.0000000000000000 0.2598189999999993 0.2000000000000000
>> gcorner 3 1.9806379930000000 0.2548077446999993 0.2000000000000000
>>
>> Corners of intersection 5 of element 43
>> gcorner 0 2.0000000000000000 0.2598190000000000 0.0000000000000000
>> *gcorner 1 1.9806379930000000 0.2548077934000000 0.0000000000000000
>> gcorner 2 2.0000000000000000 0.2598189999999995 0.2000000000000000
>> *gcorner 3 1.9806379930000000 0.2548077933999995 0.2000000000000000
>>
>> I highlight using stars the problematic corners.
>> Element 42 has no problem, while there's no coherence between the intersection 5 corners of 43 and their relative element corners.
>> If you want I can give you the faster code, but it uses dune 2.3.0.
>> But it would be better if you can drive me to point of dune or ALUGrid where the problem stems.
>> In the meanwhile, I can continue my 2.4.1 porting.
>> In any case, I use in the real code the subgrid module and I'm having some problem in building a basic project using it.
>>
>> Thank you very much for any hint or even just for your attention.
>>
>> Best regards,
>> Marco
>>
>>
>> ________________________________________
>> Da: Martin Nolte <nolte at mathematik.uni-freiburg.de>
>> Inviato: lunedì 20 giugno 2016 20.05.53
>> A: Marco Cisternino; dune at dune-project.org
>> Oggetto: Re: [Dune] jacobian construction
>>
>> Hi Marco,
>>
>> there is no call to jacobianTransposed. The jacobianTransposed is a by-product
>> of the method affine( jacobianTransposed ) on the MultiLinearGeometry.
>>
>> The affine check works as follows:
>> - a point geometry is always affine; its jacobianTransposed is trivial.
>> - a prism is affine if and only if the base geometry is affine
>> - a prism over some geometry is affine if
>> + top and bottom geometry are both affine
>> + the jacobian of top and bottom geometry equals
>>
>> You see, I potentially need the jacobianTransposed on all "subgeometries".
>> Therefore, as long as the geometry is affine, I have to compute it along the
>> way. In contrast to jacobianTransposed, I have the jacobianTransposed of the
>> base geometry or the bottom geometry already available. Therefore, the method
>> jacobianTransposed is not called at all.
>>
>> Best,
>>
>> Martin
>>
>> PS: Please keep the conversation on the mailing list.
>>
>> On 06/20/2016 06:33 PM, Marco Cisternino wrote:
>>> Hi Martin,
>>> If I have understood CachedMultiLinearGeometry in the affine case constructs the jacobianTransposed calling implicitly the constructor of the base class, i.e. MultiLinearGeometry and when it is affine the constructor of MultiLinearGeometry is called and somewhere (but I cannot see where...) this constructor is called
>>>
>>> jacobianTransposed< false >( topologyId(), integral_constant< int, mydimension >(), cit, ctype( 1 ), local, ctype( 1 ), jacobianTransposed_ );
>>>
>>> right?
>>>
>>> Thanks,
>>> Marco
>>> ________________________________________
>>> Da: Martin Nolte <nolte at mathematik.uni-freiburg.de>
>>> Inviato: lunedì 20 giugno 2016 18.09.14
>>> A: Marco Cisternino
>>> Oggetto: Re: [Dune] jacobian construction
>>>
>>> Hi Marco,
>>>
>>> I guess your question is where the CachedMultiLinearGeometry computes its
>>> jacobianTransposed. In the case of an affine mapping, this is done in the
>>> constructor. Otherwise it is not cached at all.
>>>
>>> When checking whether a MultiLinearGeometry is affine, you compute the
>>> jacobianTransposed as a by-product.
>>>
>>> Best,
>>>
>>> Martin
>>>
>>> On 06/20/2016 11:57 AM, Marco Cisternino wrote:
>>>> Good morning duners,
>>>> I'm debuging my code and I have an incoherence between the corners stored in the geometry of the element and its jacobianTransposed: the code doesn't reproduce stored corners when it re-computes them by using the jacobian.
>>>> I'm using Dune 2.3.0 and ALUGrid 1.52b, setting a GeometryGrid with a element-by-element parametrization.
>>>> I would like to know where, when and how the jacobianTransposed is computed for the first time.
>>>>
>>>> Thanks,
>>>> Marco
>>>>
>>>> _______________________________________________
>>>> Dune mailing list
>>>> Dune at dune-project.org
>>>> http://lists.dune-project.org/mailman/listinfo/dune
>>>>
>>>
>>> --
>>> Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>
>>>
>>> Universität Freiburg phone: +49-761-203-5630
>>> Abteilung für angewandte Mathematik fax: +49-761-203-5632
>>> Hermann-Herder-Straße 10
>>> 79104 Freiburg, Germany
>>>
>>
>> --
>> Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>
>>
>> Universität Freiburg phone: +49-761-203-5630
>> Abteilung für angewandte Mathematik fax: +49-761-203-5632
>> Hermann-Herder-Straße 10
>> 79104 Freiburg, Germany
>>
>
> --
> Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>
>
> Universität Freiburg phone: +49-761-203-5630
> Abteilung für angewandte Mathematik fax: +49-761-203-5632
> Hermann-Herder-Straße 10
> 79104 Freiburg, Germany
>
>
--
Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>
Universität Freiburg phone: +49-761-203-5630
Abteilung für angewandte Mathematik fax: +49-761-203-5632
Hermann-Herder-Straße 10
79104 Freiburg, Germany
More information about the Dune
mailing list