[Getdp] Line Integrals in Y-direction result in zero

Christophe Geuzaine cgeuzaine at ulg.ac.be
Fri Mar 8 17:37:40 CET 2013


Hi Oliver,

It's not a bug, it's a "feature" ;-)

More seriously, this is related to how we perform mappings of fields in GetDP with the Jacobian tool.

Basically, "Vol" means that we map an n-dimensional entity in an n-dimensional space ; "Sur" means that we map an (n-1)-dimensional entity in an n-dimensional space, and "Lin" means that we map an (n-2)-dimensional entity in an n-dimensional space.

So in your case, for the line integrals in 2D, you should use the "Sur" Jacobian. Using "Vol" for the lines implicicitly assumes that you are solving a 1-D problem, which in GetDP is defined as mapped onto the "X-axis". 

Attached is the corrected file.

Hope this helps,

Christophe

-------------- next part --------------
A non-text attachment was scrubbed...
Name: integral.pro
Type: application/octet-stream
Size: 3356 bytes
Desc: not available
URL: <http://www.geuz.org/pipermail/getdp/attachments/20130308/1693aa5b/attachment.pro>
-------------- next part --------------



On 08 Mar 2013, at 05:56, Wackerl Oliver (ED-ENB/ENG1-NA) <Oliver.Wackerl2 at mx.bosch.com> wrote:

> Hello
> I tried to implement line integrals in the solutions and I didn't get correct answers compared to other FE-Programs. I reduced the problem and tracked it down to the line integral behavior in Y-direction.
> Attached a modified to the Post-Integral example:
> I try to integrate 1 over the line length s/2 in two cases: A line in X and Y direction.
> The expected result should be in both cases s/2. Only in X direction I can observe the the result ;in Y direction the result is zero.
> Integration of areas works as expected
> Post_integral.geo:
> s=1;
> Point(1) = {0,0,0,s/10};
> Point(2) = {s,0,0,s/10};
> Point(3) = {2*s,0,0,s/10};
> Point(4) = {0,s,0,s/10};
> Point(5) = {s,s,0,s/10};
> Point(6) = {2*s,s,0,s/10};
> Point(7) = {s/4,s/2,0,s/10};
> Point(8) = {s*3/4,s/2,0,s/10};
> Point(9) = {s/2,s/4,0,s/10};
> Point(10)= {s/2,s*3/4,0,s/10};
> Point(11)= {s/2,s/2,0,s/10};
> Line(1) = {1,2};
> Line(2) = {2,3};
> Line(3) = {3,6};
> Line(4) = {6,5};
> Line(5) = {5,4};
> Line(6) = {4,1};
> Line(7) = {2,5};
> Line Loop(8) = {1,7,5,6};
> Plane Surface(9) = {8};
> Line Loop(10) = {2,3,4,-7};
> Plane Surface(11) = {10};
> Line(12) = {7,8};
> Line(13) = {9,10};
> Circle(14) = {7,11,10};
> Circle(15) = {10,11,8};
> Circle(16) = {8,11,9};
> Circle(17) = {9,11,7};
>  
> Physical Surface(1) = 9;
> Physical Surface(2) = 11;
> Physical Line(10) = 6; // left
> Physical Line(11) = 3; // right
> Physical Line(12) = 12; // Line x-direction
> Physical Line(13) = 13; // Line y-direction
> Physical Line(14) = {14,15,16,17}; // Circle
>  
>  
> Post_integral.pro:
> Group {
>   // Domain [0,2]x[0,1]
>   // Left half of Omega
>   Omega1 = Region[ {1} ];
>   // Right half of Omega
>   Omega2 = Region[ {2} ];
>  
>   Omega = Region[ {Omega1,Omega2} ];
>  
>   // Left and right boundaries
>   Gamma1 = Region[ {10} ];
>   Gamma2 = Region[ {11} ];
>   Gamma3 = Region[ {12} ];
>   Gamma4 = Region[ {13} ];
>   Gamma5 = Region[ {14} ];
>  
>   Gamma_All = Region [ {Gamma3,Gamma4,Gamma5} ];
> }
>  
> Jacobian {
>   { Name Vol ; Case { { Region All ; Jacobian Vol ; } } }
> }
>  
> Integration {
>   { Name Int ; Case { { Type Gauss ;
>     Case { { GeoElement Triangle ; NumberOfPoints 3 ; } } } } }
>   { Name IntLine ; Case { { Type Gauss ;
>     Case { { GeoElement Line ; NumberOfPoints 2 ; } } } } }
> }
>  
> Constraint {
>   { Name Dirichlet ;
>     Case {
>       { Region Gamma1 ; Value -1. ; }
>       { Region Gamma2 ; Value 1. ; }
>     }
>   }
> }
>  
> FunctionSpace {
>   { Name H1 ; Type Form0 ;
>     BasisFunction {
>       { Name sn ; NameOfCoef vn ; Function BF_Node ;
>         Support Omega ; Entity NodesOf[ All ] ; }
>     }
>     Constraint {
>       { NameOfCoef vn; EntityType NodesOf ; NameOfConstraint Dirichlet; }
>     }
>   }
> }
>  
> Formulation {
>   { Name Laplace ; Type FemEquation ;
>     Quantity {
>       { Name v ; Type Local ; NameOfSpace H1 ; }
>     }
>     Equation {
>       Galerkin { [ Dof{d v} , {d v} ] ; In Omega ; Jacobian Vol ; Integration Int ; }
>     }
>   }
> }
>  
> Resolution {
>   { Name Laplace ;
>     System {
>       { Name A ; NameOfFormulation Laplace ; }
>     }
>     Operation {
>       Generate[A] ; Solve[A] ; SaveSolution[A] ;
>     }
>   }
> }
>  
> PostProcessing {
>   { Name Laplace ; NameOfFormulation Laplace  ;
>     Quantity {
>       { Name v ; Value { Term { [ {v} ] ; In Omega ; Jacobian Vol; } } }
>       // defines the integral of v in a generic way, over the whole domain Omega
>       { Name intv ; Value { Integral { [ {v} ] ; In Omega ; Jacobian Vol; Integration Int; } } }
>          
>           // Area
>           { Name Area ; Value { Integral { [ 1 ] ; In Omega ; Jacobian Vol; Integration Int; } } }
>           // LineLength
>           { Name LineLength; Value { Integral { [ 1 ]; In Gamma_All; Jacobian Vol; Integration IntLine; } } }
>     }
>   }
> }
>  
> PostOperation{
>   { Name v ; NameOfPostProcessing Laplace;
>     Operation {
>       Print[ v , OnElementsOf Omega , File "v.pos" ] ;
>       // Integral of v over Omega (should be 0)
>       Print[ intv[Omega] , OnGlobal, Format Table, File   "Result.txt" ] ;
>       // Integral of v over Omega1 (should be -0.5)
>       Print[ intv[Omega1] , OnGlobal, Format Table, File  > "Result.txt" ] ;
>       // Integral of v over Omega2 (should be 0.5)
>       Print[ intv[Omega2] , OnGlobal, Format Table, File   >"Result.txt" ] ;
>           // Area Integral of 1 over Omega1 (should be s^2=1)
>       Print[ Area[Omega1] , OnGlobal, Format Table, File   >"Result.txt" ] ;
>           // Area Integral of 1 over Omega2 (should be s^2=1)
>       Print[ Area[Omega2] , OnGlobal, Format Table, File   >"Result.txt" ] ;
>           // Line Integral of 1 over Gamma3 (should be s/2=0.5)
>       Print[ LineLength[Gamma3] , OnGlobal, Format Table, File   >"Result.txt" ] ;
>           // Line Integral of 1 over Gamma4 (should be s/2=0.5)
>       Print[ LineLength[Gamma4] , OnGlobal, Format Table, File   >"Result.txt" ] ;
>           // Line Integral of 1 over Gamma5 (should be 2*pi*(s/4) = pi/2)
>       Print[ LineLength[Gamma5] , OnGlobal, Format Table, File   >"Result.txt" ] ;       
>     }
>   }
> }
>  
> Best
>  
> Oliver
>  
> _______________________________________________
> getdp mailing list
> getdp at geuz.org
> http://www.geuz.org/mailman/listinfo/getdp

-- 
Prof. Christophe Geuzaine
University of Liege, Electrical Engineering and Computer Science 
http://www.montefiore.ulg.ac.be/~geuzaine