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