[Getdp] Line Integrals in Y-direction result in zero
Wackerl Oliver (ED-ENB/ENG1-NA)
Oliver.Wackerl2 at mx.bosch.com
Fri Mar 8 05:56:32 CET 2013
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.geuz.org/pipermail/getdp/attachments/20130308/7e722435/attachment.html>