GitLab

Here is a little self-contained example on how to integrate quantities in post-processing. Note that GetDP can only integrate quantities on their domain of definition (specified in the function space). To integrate on arbitrary regions you can use, e.g., Plugin(Integrate) in Gmsh.

post_integral.pro

Group {
  // Domain [0,2]x[0,1]
  Omega = Region[ {1,2} ];
  // Left half of Omega
  Omega1 = Region[ {1} ];
  // Right half of Omega
  Omega2 = Region[ {2} ];
  // Left and right boundaries
  Gamma1 = Region[ {10} ];
  Gamma2 = Region[ {11} ];
}

Jacobian {
  { Name Vol ; Case { { Region All ; Jacobian Vol ; } } }
}

Integration {
  { Name Int ; Case { { Type Gauss ; 
    Case { { GeoElement Triangle ; NumberOfPoints 3 ; } } } } }
}

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 {
      Integral { [ 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; } } }
    }
  }
}

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 ] ;
      // Integral of v over Omega1 (should be -0.5)
      Print[ intv[Omega1] , OnGlobal ] ; 
      // Integral of v over Omega2 (should be 0.5)
      Print[ intv[Omega2] , OnGlobal ] ;
    }
  }
}

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};
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};

Physical Surface(1) = 9; 
Physical Surface(2) = 11;
Physical Line(10) = 6; // left
Physical Line(11) = 3; // right