--- LaplaceNeumann.cpp 2004-01-11 08:44:09.000000000 +0100 +++ LaplaceRobin.cpp 2004-01-27 23:30:19.000000000 +0100 @@ -27,12 +27,13 @@ /* * Labels: */ -#if 0 +#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 isRobin( x ) ((x)==2) +#define isNeumann( x ) ( false ) #define isDirichlet( x ) ((x)==1) #define MESHFILE "testfrontiere.msh" @@ -41,7 +42,8 @@ // 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 isRobin( x ) ((x)==4) +#define isNeumann( x ) ( false ) #define MESHFILE "testfrontiere_big.msh" #endif @@ -55,6 +57,7 @@ {1.,1.}, // idendite, le plus classique {1.,1.}, // idendite, le plus classique {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 @@ -90,13 +93,19 @@ //------------------------------------------------------------------------------ // // on prends l'exemple de teta(x,y) = x*x + y*y + #define rayon1 ( 4. ) #define rayon2 ( 10. ) + +#define TetaRobin ( 200. ) +#define G ( -2.*rayon2/( rayon2*rayon2-TetaRobin ) ) 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 +R UNeumann( R2 & P) { return rayon2 * 2.; } +R fOmega( R2 & P) { return -4.;} // laplacien constant = 4 +R GRobin( R2 & P ) { return G; } // coefficient d'echange pour la frontiere de Robin +R TetaExtRobin( R2 & P ) { return TetaRobin; } // temperature exterieur, frontiere de Robin //------------------------------------------------------------------------------ //Le programme principal @@ -157,30 +166,28 @@ R2 P = K[0]*w0 + K[1]*w1 + K[2]*w2; - R F = fOmega(P) * PI.a; + R F = fOmega(P) * PI.a * area; - Ug[i0] += F * w0 * area; //vecteur global partie interne - Ug[i1] += F * w1 * area; //assemble a partir du vecteur elementaire - Ug[i2] += F * w2 * area; + Ug[i0] += F * w0; //vecteur global partie interne + Ug[i1] += F * w1; //assemble a partir du vecteur elementaire + Ug[i2] += F * w2; } } - const R tgv = 1e30; - for (int i=0; i<Th.nv; i++) //Traitement des noeuds sur la frontiere Dirichlet - if ( isDirichlet( ((Label)Th(i)).lab ) ) - Ug[i] = UDirichlet( Th(i) ) * (Ag.D[i] = tgv); - - // participation dus aux conditions de type Neumann + // set the current Integration formular cf FormuleIntegration.cpp + const QuadratureFormular1d &FI(QF_GaussLegendre3); + assert(FI.n == 3); // verification 3 points d'integration + + /* + * Frontiere de 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(); + Label label = ((Label)edge).lab; - 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 + if ( !isNeumann(label) ) continue; // find out vertex indices.. didn't find a better way with present fem API int i0 = &(edge[0])-Th.vertices; @@ -195,13 +202,63 @@ 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 + R N = length * UNeumann(P) * PI.p; // contribution de ce point + + Ug[i0] += N * w0 ; //vecteur global partie interne + Ug[i1] += N * w1 ; //assemble a partir du vecteur elementaire } } + /* + * Frontiere de Robin + */ + for (int i=0; i< Th.neb; i++) //boucle sur les cotes (boundary edges) + { + const BoundaryEdge &edge(Th.bedges[i]); + R length = edge.length(); + Label label = ((Label)edge).lab; + R2 P; + + if ( !isRobin(label) ) continue; + + // 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); + + + // pas besoin de formule d'integration pour la matrice + P = ( R2 )( edge[1]+edge[0] )/2; + R asdf = GRobin(P)*length/6; + Ag.add2x2Matrice( i0, i1, asdf*2., asdf, asdf*.2 ); + + for (int j=0; j<FI.n; j++) //loop on integration's points + { + R N; + QuadratureFormular1d::Point PI = FI[ j ]; // position/poids de du point j + + R w0(PI.x), w1(1-PI.x); // valeurs des fonctions chapeaux + P = ((R2)edge[1])*w0 + ((R2)edge[0])*w1; // point courant + + // modification du vecteur Ug + N = TetaExtRobin(P) * GRobin(P) * length * PI.p; // contribution de ce point + Ug[i0] += N * w0; + Ug[i1] += N * w1; + } + } + + /* + * Frontiere de Dirichlet + */ + const R tgv = 1e30; + for (int i=0; i<Th.nv; i++) //Traitement des noeuds sur la frontiere Dirichlet + if ( isDirichlet( ((Label)Th(i)).lab ) ) + Ug[i] = UDirichlet( Th(i) ) * (Ag.D[i] = tgv); + + + Ag.cholesky(); // Ag = L tL Factorisation de la matrice globale apres C.L. Ug /= Ag; // Ug = Ag**-1 Ug Resolution du systeme lineaire @@ -218,19 +275,29 @@ closegraphique(); // test de validite - R max = 0; + R max = 0.; + R tetamax =0.; R tetamin =0.; + R tetamax_real =0.; R tetamin_real =0.; assert ( Ug.N() == Th.nv ); for ( int i = 0 ; i< Th.nv; i++ ) // on parcours l'ensemble des points { R diff; - R2 P; + R2 V; - P = (R2)(Th.vertices[i]); - diff = Ug[i]-TetaExact(P) ; + V = (R2)(Th.vertices[i]); + diff = Ug[i]-TetaExact(V) ; diff *= diff; // on prends le carre if ( diff>max ) max= diff; + if ( Ug[i]>tetamax ) tetamax= Ug[i]; + if ( Ug[i]<tetamin ) tetamin= Ug[i]; + if ( TetaExact(V)>tetamax_real) tetamax_real= TetaExact(V); + if ( TetaExact(V)<tetamin_real ) tetamin_real= TetaExact(V); } - printf("\nMax des carres obtenu : %f\n\n", max); + printf("\nMax des carres obtenu : %f\n", max); + printf("\nMin,Max des tetas obtenu : %f, %f", tetamin, tetamax); + printf("\nMin,Max des tetas reels : %f, %f", tetamin_real, tetamax_real); + printf("\n\n"); + return 0; }