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

[Wackerl Oliver (ED-ENB/ENG1-NA)]

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>