[Getdp] Bug in Plugin Integrate ?

[Olivier Castany]

Hello,

I think there is a bug in the plugin "Integrate" :

if I want to measure a surface (with ComputeLevelsetPositive = 1), the 
result is sometimes wrong.

It seems that it is wrong when the surface has a hole in it.

I provide in the attached files a minimal example for which 
the problem happens.

-- 
O.C.


-------------- next part --------------
l_ext = 20;
lc_ext = 2; 

L_cond = 10;
l_cond = 1;
lc_cond = 1;

xc1 = 10; 
yc1 = 12; 

Point(1) = {0,0,0,lc_ext};
Point(2) = {0,l_ext,0,lc_ext};
Point(3) = {l_ext,l_ext,0,lc_ext};
Point(4) = {l_ext,0,0,lc_ext};

Point(5) = {xc1-L_cond/2,yc1-l_cond/2,0,lc_cond};
Point(6) = {xc1+L_cond/2,yc1-l_cond/2,0,lc_cond};
Point(7) = {xc1+L_cond/2,yc1+l_cond/2,0,lc_cond};
Point(8) = {xc1-L_cond/2,yc1+l_cond/2,0,lc_cond};


Line(1) = {1,4};
Line(2) = {4,3};
Line(3) = {3,2};
Line(4) = {2,1};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,5};

Line Loop(13) = {1,2,3,4};
Line Loop(15) = {5,6,7,8};
Plane Surface(16) = {13,15};

Physical Line(18) = {5,6,7,8};		// ?lectrode 
Physical Line(19) = {1,2,3,4};		// bord ext?rieur

Physical Surface(20) = {16};		// surface o? vit le potentiel

-------------- next part --------------
/**********************
Bug dans le plugin "Integrate" ?

J'ai fait un exemple minimal d'une surface avec un trou.

-> faire la r??solution 
-> lancer le plugin "Integrate" sur le champ "phi"
-> mettre : ComputeLevelsetPositive = 1 
-> Le r??sultat obtenu est 385 alors que la surface o?? vit "phi" est de 390 (=20*20-10*1). Ceci semble ??tre un bug.

Rq : 
- sans le trou, il y a bien ??galit?? entre le r??sultat du plugin et la surface.

************************/

Group{
	D = Region[20];    // surface o?? vit le potentiel
	C1 = Region[18];   // ??lectrode ?? l'int??rieur de la boite
	S = Region[19];    // bord ext??rieur de la boite
}

Function {
	epsilon[D] = 1.;
}

Constraint {
	{ Name potentiels_imposes ; Type Assign;
	  Case {
	    { Region S ; Value 0. ; }
	    { Region C1 ; Value 1. ; }
	  }
  	}
}

FunctionSpace {
	{ Name potentiel_dans_boite ; Type Form0 ;
	  BasisFunction {
	  	{ Name w_n ; NameOfCoef v_n ; Function BF_Node ;
		  Support D ; Entity NodesOf[ All ] ; }
	  }
	  Constraint {
	  	{ NameOfCoef v_n ; EntityType NodesOf ; 
		  NameOfConstraint potentiels_imposes ; }
	  }
	}
}

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

Integration {
  { Name I ;
    Case { 
      { Type Gauss ;
        Case { 
          { GeoElement Triangle    ; NumberOfPoints  4 ; }
          { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
        }
      }
    }
  }
}

Formulation {
  { Name formulation_du_probleme ; Type FemEquation ;
    Quantity { 
      { Name phi ; Type Local ; NameOfSpace potentiel_dans_boite ; }
    }
    Equation {
      Galerkin { [ epsilon[] * Dof{d phi} , {d phi} ] ;  
                 In D ; Jacobian J ; Integration I ; }
    }
  }
}

Resolution {
  { Name resolution_du_probleme ;
    System {
      { Name A ; NameOfFormulation formulation_du_probleme ; }
    }
    Operation { 
      Generate[A] ; Solve[A] ; SaveSolution[A] ; 
    }
  }
}

PostProcessing {
  { Name post_processing ; NameOfFormulation formulation_du_probleme ;
    Quantity {
      { Name phi; Value { Local { [ {phi} ] ; In D ; Jacobian J ; }}}
      { Name E; Value { Local { [ -{d phi} ] ; In D ; Jacobian J ; }}}

// Remarque : en calculant la surface comme ce qui suit, on trouve bien 390 :
// (?? ce propos, y a-t-il un moyen d'envoyer ce r??sultat sur la console sans
//  l'afficher dans gmsh ???)

      { Name surf; Value { Integral { [ 1 ] ; In D ;
                   Integration I ; Jacobian J ; }}}

    }
  }
}

PostOperation {
  { Name post_operations ; NameOfPostProcessing post_processing;
    Operation {
      Print[ phi, OnElementsOf D, File "phi.pos"] ;
      Print[ E,   OnElementsOf D, File "E.pos" ] ;
      Print[ surf [D], OnGlobal , File "surf.result", Format Table];
    }
  }
}