--- LaplaceConductivite.cpp	2004-01-27 23:04:11.000000000 +0100
+++ LaplaceNeumann.cpp	2004-01-11 08:44:09.000000000 +0100
@@ -24,6 +24,26 @@
 
 using namespace Fem2D;
 
+/*
+ * Labels:
+ */
+#if 0 
+//  testfrontiere.msh : 
+//      1   {est le NUMERO ou LABEL dans LIGNE      de } C1                      
+//      2   {est le NUMERO ou LABEL dans LIGNE      de } C2                      
+// que du robin
+#define isNeumann( x ) ((x)==2)
+#define isDirichlet( x ) ((x)==1)
+#define MESHFILE "testfrontiere.msh"
+
+#else
+//  testfrontiere_big.msh : 
+//    3   {est le NUMERO ou LABEL dans LIGNE      de } C1B
+//    4   {est le NUMERO ou LABEL dans LIGNE      de } C2B
+#define isDirichlet( x ) ((x)==3)
+#define isNeumann(  x ) ((x)==4)
+#define MESHFILE "testfrontiere_big.msh"
+#endif
 
 /*
  * Donnees du tenseur de conductivite
@@ -33,9 +53,11 @@
 #define MaxLabel 4
 R tenseurConductivite[][ 2 ] = {
 	{1.,1.}, // idendite, le plus classique
-	{.5,2.}, // surtout sur l'axe des y (anisotrope)
-	{2.,.5}, // surtout sur l'axe des x (anisotrope)
-	{0.,1.}, // rien selon l'axe des x
+	{1.,1.}, // idendite, le plus classique
+	{1.,1.}, // idendite, le plus classique
+	{.5,2.}, // surtout sur l'axe des y
+	{2.,.5}, // surtout sur l'axe des x
+	{0.,1.}, // rien selon l'ace des x
 };
 
 //------------------------------------------------------------------------------
@@ -66,10 +88,15 @@
 }
 
 //------------------------------------------------------------------------------
-//Donnees de U Dirichlet (Solution exacte) et F Omega
-
-R UDirichlet( R2 & P) { return  3*P.x+5*P.y;}//Condition de Dirichlet = solution
-R fOmega( R2 & P) { return  1;}              //Source sur Omega
+//
+// on prends l'exemple de teta(x,y) = x*x + y*y
+#define rayon1 ( 4. )
+#define rayon2 ( 10. )
+R TetaExact ( R2 & P ) { return P.x*P.x + P.y*P.y; }
+
+R UDirichlet( R2 & P) { return rayon1*rayon1; }
+R UNeumann(  R2 & P) { return rayon2 * 2; }
+R fOmega( R2 & P) { return  -4.;} // laplacien constant = -4
 
 //------------------------------------------------------------------------------
 //Le programme principal
@@ -82,7 +109,7 @@
                   //trace avec X11-Windows System du MIT
 
   //Le nom du fichier de maillage
-  const char * fn= (argc<2) ? "carrecercle.msh" : argv[1]; 
+  const char * fn= (argc<2) ? MESHFILE : argv[1]; 
   //peut etre le 2-eme mot de la commande Exemple:  Laplace MonMaillage.msh
 
   Mesh Th(fn);  //grep -in 'class Mesh' *.hpp => fem.hpp pour trouver definition
@@ -140,9 +167,41 @@
 
   const R tgv = 1e30;
   for (int i=0; i<Th.nv; i++) //Traitement des noeuds sur la frontiere Dirichlet
-    if ( !!(Label) Th(i) )    //just for boundary Dirichlet vertex 
+    if ( isDirichlet( ((Label)Th(i)).lab ) )
       Ug[i] =  UDirichlet( Th(i) ) * (Ag.D[i] = tgv);
 
+  // participation dus aux conditions de type Neumann
+  for (int i=0; i< Th.neb; i++) //boucle sur les cotes (boundary edges)
+  {
+    const BoundaryEdge &edge(Th.bedges[i]);
+    R length = edge.length();
+
+    if ( !isNeumann( ((Label)edge).lab ) ) continue; // not neumann -> continue the loop
+
+    // set the current Integration formular cf FormuleIntegration.cpp
+    const QuadratureFormular1d &FI(QF_GaussLegendre3);
+    assert(FI.n == 3); // verification 3 points d'integration
+
+    // find out vertex indices.. didn't find a better way with present fem API
+    int i0 =  &(edge[0])-Th.vertices;
+    assert(i0>=0); assert(i0<Th.nv);
+    int i1 =  &(edge[1])-Th.vertices;
+    assert(i1>=0); assert(i1<Th.nv);
+
+
+    for (int j=0; j<FI.n; j++) //loop on integration's points
+    {
+      QuadratureFormular1d::Point PI = FI[ j ];  // position/poids de du point j
+
+      R  w0(PI.x), w1(1-PI.x); // valeurs des fonctions chapeaux
+      R2 P = ((R2)edge[0])*w0 + ((R2)edge[1])*w1; // point courant
+      R  N = UNeumann(P) * PI.p; // contribution de ce point
+
+      Ug[i0] += N * w0 * length;  //vecteur global partie interne
+      Ug[i1] += N * w1 * length;  //assemble a partir du vecteur elementaire
+    }
+  }
+
   Ag.cholesky();  // Ag = L tL   Factorisation de la matrice globale apres C.L.
   Ug /= Ag;       // Ug = Ag**-1 Ug Resolution du systeme lineaire
 
@@ -157,6 +216,21 @@
   closePS();
   rattente(1);
   closegraphique();
+
+  // test de validite
+  R max = 0;
+  assert (  Ug.N() == Th.nv );
+  for ( int i = 0 ; i< Th.nv; i++ ) // on parcours l'ensemble des points
+  {
+	  R diff;
+	  R2 P;
+
+	  P = (R2)(Th.vertices[i]);
+	  diff = Ug[i]-TetaExact(P) ;
+	  diff *= diff; // on prends le carre
+	  if ( diff>max ) max= diff;
+  }
+  printf("\nMax des carres obtenu : %f\n\n", max);
   
   return 0;
 }