Notas sobre FreeFem++ 2D y 3D: traduccin del manual de F.

Hecht o
Eliseo Chacn Vera o
Departamento de Ecuaciones Diferenciales y Anlisis Numrico, a e Facultad de Matemticas, Universidad de Sevilla a

10 de julio de 2009

1.

Introduccin a FreeFem++ 2D o

Los archivos de rdenes se guardan con la terminacion .edp, el nombre o no es importante pero si lo es la terminacin. Normalmente, en un entorno o Linux, se ejecutan mediante la orden FreeFem++ [Link] y en un entorno Windows presionando dos veces con el botn derecho del o ratn sobre el archivo. La gramtica es similar a la del lenguaje C/C++. o a Por ejemplo, las variables se declaran siempre, las distintas rdenes terminan o con ; , la orden cout indica la salida por defecto, pantalla o archivo, y \n signica un salto de linea. Algunos ejemplos de aritmtica es: e Archivo [Link] real x=3.14,y; int i,j; complex c; cout<< " x = "<< x << "\n"; x=1;y=2;x=y;i=0;j=1; cout<< 1+3 << " " << 1./3 << "\n"; cout<< 10^10. << "\n"; cout<< 10^-10. << "\n"; cout<< -10^-2.+5 << "==4.99 \n"; cout<< 10^-2.+5 << "==5.01 \n"; cout<< " 8^(1/3) = "<< (8)^(1./3.) <<"\n"; cout<< "-------------complex ---------- \n"; 1

 2 GENERACION DE MALLAS cout<< 10-10i <<"\n"; cout<< " -1^(1/3) = "<< (-1+0i)^(1./3.) <<"\n"; Con [Link] y [Link] se puede ver el uso de ciclos. Observad i + + que equivale a i = i + 1. Archivo [Link]: int i; for (i=0;i<10;i=i+1) cout <<i <<\n; Archivo [Link]: int i; i=0; real eps=1; while (eps>1e-5) {eps=eps/2; if(i++>100) break; cout <<i <<\n; cout <<eps <<endl; };

2.

Generacin de mallas o

Para crear una malla de un rectngulo se puede usar una descripcin a o paramtrica de cada lado del mismo junto con una etiqueta, mediante la e orden label= que sirva para asociar condiciones de contorno. Archivo [Link]: border a(t=0,2){x=t;y=0;label=1;}; border b(t=0,1){x=2;y=t;label=2;}; border c(t=0,2){x=2-t;y=1;label=3;}; border d(t=0,1){x=0;y=1-t;label=4;}; int n=5; mesh th=buildmesh(a(2*n)+b(n)+c(2*n)+d(n)); plot(th,wait=1,ps=[Link]);

Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

Malla no uniforme generada con [Link] Con la orden border se dene un segmento del contorno de la geometr y a con las ordenes mesh y buildmesh se genera la malla. La orden square genera por defecto el rectngulo unidad [0, 1]2 . Tambin puede servir para a e generar un rectngulo cualquiera y adems se crea una malla uniforme. Se a a pueden dibujar dos mallas al mismo tiempo con la orden plot: Archivo [Link] real x0=1.2,x1=1.8; real y0=0.,y1=1.; int n=5,m=20; mesh Th=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]); mesh th=square(4,5); plot(Th,th, ps=[Link]);

Dos mallas uniformes generadas con [Link] La opcin ps = . . . genera un chero postscript. La orden square genera o una malla uniforme en el cuadrado y la numeracin por defecto de los o 3 Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

lados es 1,2,3,4 para los lados de base, derecha, superior, izquierda. La orden buildmesh genera una malla no uniforme. Se puede alterar la numeracin de los lados mediante la orden label o que dota de un nmero de referencia para poder asignar datos de contorno. u Varios lados pueden tener el mismo nmero de referencia si van a tener el u mismo dato de contorno. Con el siguiente ejemplo [Link] se crea un dominio con un hueco simplemente cambiando el sentido de orientacin de los nodos sobre la frono tera interna: Archivo [Link] real pi=4*atan(1.0); border a(t=0,2*pi){x=2*cos(t);y=sin(t);label=1;} border b(t=0,2*pi){x=.3+.3*cos(t);y=.3*sin(t);label=2;} mesh thwithouthole=buildmesh(a(50)+b(30)); mesh thwithhole=buildmesh(a(50)+b(-30)); plot(thwithouthole,wait=1,ps=[Link]); plot(thwithhole,wait=1,ps=[Link]);

Dos mallas complementarias generadas con [Link]

Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

Mallas con hueco generada con [Link] El siguiente ejemplo [Link] muestra una geometr poligonal y se usan las a rdenes savemesh y readmesh para guardar y leer una malla: o Archivo [Link] border a(t=0,1){x=t;y=0;label=1;}; border b(t=0,0.5){x=1;y=t;label=1;}; border c(t=0,0.5){x=1-t;y=0.5;label=1;}; border d(t=0.5,1){x=0.5;y=t;label=1;}; border e(t=0.5,1){x=1-t;y=1;label=1;}; border f(t=0,1){x=0;y=1-t;label=1;}; int n=5; mesh rh=buildmesh(a(n)+b(n)+c(n)+d(n)+e(n)+f(n)); plot(rh,wait=1); savemesh(rh,"[Link]"); mesh th=readmesh("[Link]"); plot(th,wait=1,ps="[Link]");

Malla para L generada con [Link] 5 Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

Resolucin de problemas variacionales o

Tres palabras clave: problem, solve y varf problem ejemplo()...; dene el problema variacional construyendo el sistema lineal y asociandole una manera de invertirlo. El problema variacional se resuelve al usar la orden ejemplo; solve ejemplo()...; dene y resuelve el problema variacional construyendo el sistema lineal y asociandole una manera de invertirlo. varf construye las partes del problema variacional y permite extraer la matriz y el termino independiente del problema. Resolutores lineales: Bsicamente tres: CG, UMFPACK y GMRES. a Por defecto siempre se usa UMFPACK pero para sistemas grandes es mejor usar GMRES. Se ja el resolutor del sistema lineal con la orden solver=... CG Gradiente conjugado, es un mtodo directo. e UMFPACK es un mtodo directo, permite manejar cualquier tipo e de matriz mediante una factorizacin. GMRES mtodo iterativo para matrices grandes y huecas. e El ejemplo [Link] resuelve el problema de Poisson en un escaln hacia adelante con datos de contorno de tipo Dirichlet: o Archivo [Link] real v1=4,v2=6,w=v1+v2; border a(t=0,v1){x=t;y=0;label=1;}; //lado y=0, x en [0,v1] border b(t=0,1){x=v1;y=-t;label=2;};//lado x=v1, y en [0,-1] border c(t=0,v2){x=v1+t;y=-1;label=3;}; //lado y=-1, x en [v1,w] border d(t=0,2){x=w;y=-1+t;label=4;};//lado x=w, y en [-1,1] border e(t=0,w){x=w-t;y=1;label=5;};//lado y=1, x en [0,w] border ff(t=0,1){x=0;y=1-t;label=6;};//lado x=0, y en [1,0] int n=10; mesh th=buildmesh(a(2*n)+b(n)+c(2*n)+d(2*n)+e(4*n)+ff(n)); plot(th,wait=1,ps="[Link]"); fespace Vh(th,P1); Vh u,v; func f=0; problem laplace(u,v) = int2d(th)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinear part 6 Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

+int2d(th)(-f*v) // right hand side +on(1,u=0) +on(2,u=0) +on(3,u=0) +on(4,u=0) +on(5,u=1) +on(6,u=0);// Dirichlet boundary condition real cpu=clock(); laplace; //solve the pde plot(u,wait=1,ps="[Link]"); Observacin 1 Las funciones int2d(th)(...) y similares, se usan o con la funcin test y la incognita o o slo con la funcin test o o pero no ambas situaciones. Es decir, es incorrecto int2d(th)(dx(u)*dx(v)+dy(u)*dy(v)-f*v) +on(1,u=0) y es correcto int2d(th)(dx(u)*dx(v)+dy(u)*dy(v)) +int2d(th)(-f*v) +on(1,u=0) // bilinear part // right hand side

Malla generada con [Link]

Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

Solucin obtenida con [Link] o El ejemplo [Link] resuelve el problema de Poisson en un escaln hacia adelante con datos de contorno de tipo Dirichlet en o todos los lados salvo en uno en el que el dato es de tipo Neumann, derivada normal cero. Esto implica que no hay especicacione sobre este trozo de frontera: Archivo [Link] real v1=4,v2=6,w=v1+v2; border a(t=0,v1){x=t;y=0;label=1;}; //lado y=0, x en [0,v1] border b(t=0,1){x=v1;y=-t;label=1;};//lado x=v1, y en [0,-1] border c(t=0,v2){x=v1+t;y=-1;label=1;}; //lado y=-1, x en [v1,w] border d(t=0,2){x=w;y=-1+t;label=2;};//lado x=w, y en [-1,1] border e(t=0,w){x=w-t;y=1;label=2;};//lado y=1, x en [0,w] border (t=0,1){x=0;y=1-t;label=1;};//lado x=0, y en [1,0] int n=10; mesh th=buildmesh(a(2*n)+b(n)+c(2*n)+d(2*n)+e(4*n)+(n)); plot(th,wait=1); fespace Vh(th,P1); Vh u,v; func f=1; problem laplace(u,v) = int2d(th)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinear part +int2d(th)(-f*v) // right hand side +on(1,u=0); // Dirichlet boundary condition laplace; //solve the pde plot(u,wait=1,value=1,ps=[Link]);

Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

IsoValue -0.0263321 0.0131661 0.0394982 0.0658303 0.0921624 0.118494 0.144827 0.171159 0.197491 0.223823 0.250155 0.276487 0.302819 0.329151 0.355483 0.381816 0.408148 0.43448 0.460812 0.526642

Soluc obtenida con [Link] on El siguiente cdigo [Link] resuelve mediante elementos nitos o P2 el problem variacional ( u, v) = (f, v) con dato de contorno u = 0 en una elipse: Archivo [Link] border a(t=0,2*pi){x=2*cos(t);y=sin(t);label=5;}; mesh Th= buildmesh (a(150)); plot(Th,wait=1,ps=[Link]); fespace Vh(Th,P2); Vh u, v; func f=sin(x*y); problem laplace(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinear part -int2d(Th)(f*v) // right hand side +on(5,u=0); // Dirichlet boundary condition real d=cpu(); Get current CPU clock time laplace; //solve the pde cout<< CPU = << clock() - d << endl; Gives CPU time for the process plot(u,wait=1,ps=[Link]); a Observacin 2 El tiempo usado en el clculo se puede obtener con o real d=cpu();\\Get current CPU clock time ... cout$<<$ "\, CPU = " $<<$ clock() - d $<<$ endl; \\Gives CPU time \\for the process 9 Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

se obtiene el tiempo de CPU empleado por todas las instrucciones entre ambas l neas.

Malla para elipse generada con [Link]

Solucin obtenida con [Link] o Se puede realizar una adaptacin de mallado. El cdigo [Link] o o es un ejemplo: Se resuelve en un rectngulo y mediante elementos nitos P1 a el problem variacional ( u, v) = (f, v) con dato de contorno u = g. Archivo [Link] verbosity=0; int n=5; mesh Th=square(n,n,[10*x,5*y]);// Malla inicial uniforme en [0,10]x[0,5] plot(Th,wait=1,ps="[Link]"); fespace Vh(Th,P1); Vh u=0,v=0,zero=0; 10 Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

func f=-sin(x*y); func g=x+y; int i=0; real error=0.1, coef=0.1^(1./5.); real cpu=clock(); problem Problem1(u,v,solver=CG,init=i,eps=-1.0e-6)= int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v)) +int2d(Th)( v*f) +on(1,2,3,4,u=g); Problem1; plot(u,zero,wait=1,ps="[Link]"); Th= adaptmesh(Th,u,inquire=1,err=error); real d=clock(); cout << " CPU = "<<cpu -d <<endl; for (i=1;i<5;i++) { Problem1; plot(Th,u,zero,wait=1); Th= adaptmesh(Th,u,inquire=1,err=error); cout << " CPU = "<<clock() -d <<endl; d=clock(); error=error*coef; }; cout << " CPU = "<<clock() -cpu <<endl; plot(Th,u,zero,wait=1,ps="[Link]"); plot(Th,wait=1,ps="[Link]"); La opcin init=i con i = 0 fuerza a reconstruir el sistema lineal cada vez. o Esto es lgico puesto que se ha remallado. o El parmetro eps=... indica un criterio de parada. a Si eps < 0 el criterio es el error absoluto menor que eps Ax b |eps|. Si eps > 0 el criterio es el error relativo Ax b eps Ax0 b .

11

Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

Malla inicial generada con [Link]

Solucion inicial obtenida con [Link]

Malla nal generada con [Link]

12

Variaciones sobre manual de Frdric Hecht e e

 2 GENERACION DE MALLAS

Solucion nal obtenida con [Link] Observad como los comentarios se escriben despus de // e La orden fspace Vh(Th,P1) construye un espacio de elementos nitos P1 de acuerdo a la malla especicada Th y de nombre Vh. Con la orden Vh u,v,. . . espacio vectorial. se declaran distintos elementos de este

Con la orden func g=. . . se puede declarar una funcion de acuerdo a una expresin anal o tica. La orden problem dene un problema tipo y se le asigna un nombre, en este caso Problem1. Los sistemas lineales se resuelven mediante CG que equivale a Conjugate Gradient o mediante LU que equivale a la factorizacin LU. o La malla se renueva mediante la orden adaptmesh Con la orden adaptmesh(Th,f ) se reconstruye T h con respecto a la funcin o f usando su Hessiano. Por ejemplo, la funcin o f (x, y) = 10 x3 + y 3 + atan1 [ /(sin(5 y) 2 x)], = 0,0001

cambia muy rapido y una malla uniforme no puede verla. Observemos como se adapta la malla: Cdigo [Link] o real eps = 0.0001; real h=1; real hmin=0.05; 13 Variaciones sobre manual de Frdric Hecht e e

 3 GRAFICAS

func f = 10.0*x^3+y^3+h*atan2(eps,sin(5.0*y)-2.0*x); mesh Th=square(5,5,[-1+2*x,-1+2*y]); fespace Vh(Th,P1); Vh fh=f; plot(Th,fh, ,ps="[Link]"); for (int i=0;i<5;i++) { Th=adaptmesh(Th,fh); fh=f; // old mesh is deleted plot(Th,fh,wait=1,ps="[Link]"); }

Malla y funcin inicial generada con [Link] o

Malla nal y funcin generada con [Link] o

3.

Grcas a

Existen varias opciones una vez que se tiene una grca en pantalla. a Con la tecla ? se obtiene un menu de informacion con las distintas opciones. 14 Variaciones sobre manual de Frdric Hecht e e

 3 GRAFICAS

Algunas de ellas son: p se muestra la grca anterior. a f click del ratn avanza el proceso o o +,- hace un zoom + o - en torno a la posicin del ratn o o = restaura el zoom g cambia de gr a color s m muestra la malla 3 muestra una versin 3D de la grca o a b cambia de color a blanco y negro c y C aumenta o disminuye el coeciente de las echas f rellena color entre l neas de corriente v muestra los valores numricos en las l e neas de corriente El siguiente cdigo [Link] hace lo mismo pero con o un problema mas complicado. Se resuelve en un rectngulo y mediante elea mentos nitos P2 el problem variacional ( u, v) + (g u, v) = (f, v) con dato de contorno u = 0. Archivo [Link] real a=1,b=1; int n=4,m=4; mesh Th=square(n,m,[a*x,b*y]);// Malla uniforme en [0,a]x[0,b] plot(Th,wait=1,ps="[Link]"); fespace Vh(Th,P2); Vh u=0,v=0; func f=1; func g=10*x*y; int i=0; real error=0.1, coef=0.1^(1./5.); problem Problem1(u,v,solver=CG,init=i,eps=-1.0e-6)= int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v)+g*u*v) -int2d(Th)( v*f)+on(1,u=0) 15 Variaciones sobre manual de Frdric Hecht e e

 3 GRAFICAS

+on(2,u=1)+on(3,4,u=0); real cpu=clock(); Problem1; plot(u,wait=1,ps="[Link]"); for (i=0;i<5;i++){ real d=clock(); Problem1; cout << " CPU = " << clock() -d << endl; plot(u,wait=1); Th= adaptmesh(Th,u,inquire=1,err=error); error=error*coef; }; plot(u,wait=1,ps="[Link]"); plot(Th,wait=1,ps="[Link]");

Malla inicial generada con [Link]

Solucion inicial obtenida con [Link]

16

Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS VARIACIONALES

Malla nal generada con [Link]

Solucion nal obtenida con [Link]

4.

Problemas variacionales
FreeFem provee de varios tipos de elementos nitos: P0 P1 P1dc ( P1 discontinuo ) P1b ( P1 + burbuja ) P2 P2dc RT0 ( Raviart-Thomas )

17

Variaciones sobre manual de Frdric Hecht e e

4.1

Orden de convergencia

4 PROBLEMAS VARIACIONALES

Ejemplo de como se dene el espacio de elementos nitos una vez que la malla se tiene: real a=1.0,b=1.0; int n=10,m=10; mesh Th=square(n,m,[a*x,b*y]);// Malla uniforme en [0,a]x[0,b] plot(Th,wait=1); fespace Vh(Th,P2); Vh u,v,zero; Vh[int] uh(5); ... La ultima orden permite denir una funcin vectorial de 5 componentes o donde cada componente es un elemento de Vh.

4.1.

Orden de convergencia

Se puede comprobar el orden de convergencia de los elementos nitos programados con el siguiente cdigo [Link] o verbosity =0; real a=1.5,b=1.; //border Gamma(t=0,2*pi){x =a*cos(t);y=b*sin(t);} real ratio; int npasos=10; func phiexact=sin(x)*sin(y); func fexact=2*phiexact; real[int] ErrorH1(npasos),ErrorL2(npasos),ta(npasos); for(int n=0;n<npasos;n++) { // mesh Th=buildmesh(Gamma(10*(n+1))); mesh Th=square(10*(n+1),10*(n+1)); int nbverts=[Link]; fespace Mh(Th,P0); Mh hhh = hTriangle; cout << " mesh size = "<< hhh[].max <<" nb of vertices = " << nbverts<<endl; ta[n]=hhh[].max; fespace Vh(Th,P2); Vh phi,w,errorh,phih,fh; phih=phiexact; 18 Variaciones sobre manual de Frdric Hecht e e

4.1

Orden de convergencia

4 PROBLEMAS VARIACIONALES

//

fh=fexact; solve laplace(phi,w)=int2d(Th)(dx(phi)*dx(w) + dy(phi)*dy(w)) - int2d(Th)(fh*w) +on(Gamma,phi=phih); +on(1,2,3,4,phi=phih); errorh=phi-phih; ErrorH1[n]= int2d(Th)(dx(errorh)^2+dy(errorh)^2); ErrorL2[n]= int2d(Th)(errorh^2);

}; cout << "++++++++++++++++++++++++++++++++++++++ " <<endl; for(int n=0;n<npasos;n++) { cout << " Error seminorm H1 " << n << " = "<<sqrt( ErrorH1[n]) <<endl; }; for(int n=0;n<npasos-1;n++) { ratio=log(sqrt(ErrorH1[n])/sqrt(ErrorH1[n+1])); ratio=ratio/log(ta[n]/ta[n+1]); cout <<" convergence rate error semi norm H1= "<< ratio <<endl; }; cout << "++++++++++++++++++++++++++++++++++++++ " <<endl; for(int n=0;n<npasos;n++) { cout << " ErrorL2 " << n << " = "<< sqrt(ErrorL2[n]) <<endl; }; for(int n=0;n<npasos-1;n++) { ratio=log(sqrt(ErrorL2[n])/sqrt(ErrorL2[n+1])); ratio=ratio/log(ta[n]/ta[n+1]); cout <<" convergence rate error L2= "<< ratio <<endl; }; for(int n=0;n<npasos;n++) { cout << " Error norm H1 " << n << " = "<< sqrt(ErrorL2[n]+ErrorH1[n]) <<endl; }; for(int n=0;n<npasos-1;n++) { ratio=log(sqrt(ErrorL2[n]+ErrorH1[n])/sqrt(ErrorL2[n+1]+ErrorH1[n+1])); ratio=ratio/log(ta[n]/ta[n+1]); cout <<" convergence rate error norm H1= "<< ratio <<endl; 19 Variaciones sobre manual de Frdric Hecht e e

4.1

Orden de convergencia

4 PROBLEMAS VARIACIONALES

}; los resultados para elementos nitos P2 sobre una triangulacin uniforme o son size of mesh = 0.141421 nb of vertices = 121 size of mesh = 0.0707107 nb of vertices = 441 size of mesh = 0.0471405 nb of vertices = 961 size of mesh = 0.0353553 nb of vertices = 1681 size of mesh = 0.0282843 nb of vertices = 2601 size of mesh = 0.0235702 nb of vertices = 3721 size of mesh = 0.0202031 nb of vertices = 5041 size of mesh = 0.0176777 nb of vertices = 6561 size of mesh = 0.0157135 nb of vertices = 8281 size of mesh = 0.0141421 nb of vertices = 10201 ++++++++++++++++++++++++++++++++++++++ Error seminorm H1 0 = 1.29502e-005 Error seminorm H1 1 = 1.69264e-006 Error seminorm H1 2 = 5.08725e-007 Error seminorm H1 3 = 2.1613e-007 Error seminorm H1 4 = 1.11122e-007 Error seminorm H1 5 = 6.44848e-008 Error seminorm H1 6 = 4.06886e-008 Error seminorm H1 7 = 2.72984e-008 Error seminorm H1 8 = 1.91945e-008 Error seminorm H1 9 = 1.40056e-008 convergence rate error semi norm H1= 2.93563 convergence rate error semi norm H1= 2.96483 convergence rate error semi norm H1= 2.9756 convergence rate error semi norm H1= 2.98129 convergence rate error semi norm H1= 2.98481 convergence rate error semi norm H1= 2.98722 convergence rate error semi norm H1= 2.98896 convergence rate error semi norm H1= 2.99028 convergence rate error semi norm H1= 2.99133 ++++++++++++++++++++++++++++++++++++++ ErrorL2 0 = 2.89877e-007 ErrorL2 1 = 1.85985e-008 ErrorL2 2 = 3.69895e-009 ErrorL2 3 = 1.17389e-009 20 Variaciones sobre manual de Frdric Hecht e e

4.1

Orden de convergencia

4 PROBLEMAS VARIACIONALES

ErrorL2 4 = 4.81636e-010 ErrorL2 5 = 2.3251e-010 ErrorL2 6 = 1.25613e-010 ErrorL2 7 = 7.36696e-011 ErrorL2 8 = 4.60503e-011 ErrorL2 9 = 3.02595e-011 convergence rate error L2= 3.96218 convergence rate error L2= 3.98316 convergence rate error L2= 3.98955 convergence rate error L2= 3.99246 convergence rate error L2= 3.99434 convergence rate error L2= 3.99433 convergence rate error L2= 3.99618 convergence rate error L2= 3.98916 convergence rate error L2= 3.98559 Error norm H1 0 = 1.29534e-005 Error norm H1 1 = 1.69274e-006 Error norm H1 2 = 5.08739e-007 Error norm H1 3 = 2.16133e-007 Error norm H1 4 = 1.11123e-007 Error norm H1 5 = 6.44852e-008 Error norm H1 6 = 4.06888e-008 Error norm H1 7 = 2.72985e-008 Error norm H1 8 = 1.91946e-008 Error norm H1 9 = 1.40056e-008 convergence rate error norm H1= 2.9359 convergence rate error norm H1= 2.96491 convergence rate error norm H1= 2.97564 convergence rate error norm H1= 2.98131 convergence rate error norm H1= 2.98483 convergence rate error norm H1= 2.98723 convergence rate error norm H1= 2.98897 convergence rate error norm H1= 2.99029 convergence rate error norm H1= 2.99133 Observacin 3 Se ve perfectamente el orden de convergencia en el caso o de una triangulacin uniforme. Pero no se puede estimar de la misma o forma cuando la triangulacin deja de serlo. o La informacin sobre la triangulacin se obtiene con este conjunto de ordenes o o 21 Variaciones sobre manual de Frdric Hecht e e

4.2

Salida y lectura a archivos de datosPROBLEMAS VARIACIONALES 4

mesh Th=square(10*(n+1),10*(n+1)); int nbvertices=[Link]; fespace Mh(Th,P0); Mh hhh = hTriangle; cout << "size of mesh = "<< hhh[].max <<" nb of vertices = " << nbvertices<<endl;; Se usa que hTriangle es el tamao de la arista ms larga de cada tringulo n a a y que [Link] es el nmero de vrtices de la triangulacin. u e o

4.2.

Salida y lectura a archivos de datos

Para guardar un resultado usamos { ofstream u1data("[Link]"); u1data << uu1[]; ofstream u2data("[Link]"); u2data << uu2[]; } y para leer un dato de un chero y guardarlo en una funcin de elementos o nitos usamos { ifstream u1data("[Link]"); u1data >> uu1[]; ifstream u2data("[Link]"); u2data >> uu2[]; } donde uu1 o uu2 son las funciones de elementos nitos y uu1[] el vector de valores asociado a cada funcion de elementos nitos. Observacin 4 Al igual que en C/C++ el entorno o { .... } plantea un uso local de memoria, es decir, las variables denidas dentro son locales.

22

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

5.

Ecuaciones de Stokes y Navier-Stokes

El siguiente cdigo [Link] resuelve el problema de Stokes o en un cuadrado [0, 1] [0, 1] con una malla rectangular y con los datos de contorno del problema de la cavidad. Se usa el par de elementos nitos (P2,P1) para las incognitas velocidad-presin. Tambin se obtienen las lio e neas de corriente dadas por tal que rot() = u, o mejor = rot(u) = y u1 x u2 .

Observacin 5 La compatibilidad de la condicin inf-sup a nivel continuo o o 1 ())d y L2 (). Por lo tanto, el par (P2,P1) no es la tienen los espacios (H0 0 compatible. La compatibilidad correcta se encuentra con el par (P2,P1/R). Para evitar esta situacin, que tambin producir adems una matriz llena, o e a a se le a ade el trmino p q que hace que el problema continuo se plantee n e 1 en los espacios (H0 ())d y L2 () siendo ya compatible el sistema y distando la solucin de la correcta un orden de . o Observacin 6 UMFPACK elimina los autovalores cero de la matriz del o sistema lineal que resuelve. Es por ello que en el caso = 0 tambin da la e solucin correcta. o Observacin 7 UMFPACK utiliza la tcnica de bloqueo mediante un valor o e grande para jar los datos de tipo Dirichlet. Archivo [Link]: mesh Th=square(10,10); fespace Xh(Th,P2); fespace Mh(Th,P1); Xh u1,v1,u2,v2; Mh p,q; solve Stokes ([u1,u2,p],[v1,v2,q]) = int2d(Th)( ( dx(u1)*dx(v1)+dy(u1)*dy(v1) + dx(u2)*dx(v2)+dy(u2)*dy(v2) ) +p*q*0.000001 -p*dx(v1)-p*dy(v2) 23 Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

+dx(u1)*q+dy(u2)*q ) +on(1,2,4,u1=0,u2=0)+on(3,u1=1,u2=0); plot(coef=.5,cmm=" Velocidad [u1,u2] ", value=true,[u1,u2],wait=1,ps="[Link]"); plot(cmm=" Presion ",p,value=true,wait=1,ps="[Link]"); plot([u1,u2],p,wait=1,value=true,coef=0.5,ps="[Link]"); Xh psi,phi; solve streamlines(psi,phi) = int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi) ) +int2d(Th)( -phi*(dy(u1)-dx(u2)) ) +on(1,2,3,4,psi=0); plot(cmm=" Streamlines "psi,wait=1, ps="[Link]");

24

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Velocidad [u1,u2]

Vec Value 0 0.0526343 0.105269 0.157903 0.210537 0.263172 0.315806 0.36844 0.421075 0.473709 0.526343 0.578978 0.631612 0.684247 0.736881 0.789515 0.84215 0.894784 0.947418 1.00005

Velocidad obtenida con [Link]

25

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Presion

IsoValue -49.71 -44.697 -39.6841 -34.6711 -29.6581 -24.6452 -19.6322 -14.6192 -9.60627 -4.5933 0.41967 5.43264 10.4456 15.4586 20.4715 25.4845 30.4975 35.5104 40.5234 45.5364

Presin obtenida con [Link] o

26

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

IsoValue -49.71 -44.697 -39.6841 -34.6711 -29.6581 -24.6452 -19.6322 -14.6192 -9.60627 -4.5933 0.41967 5.43264 10.4456 15.4586 20.4715 25.4845 30.4975 35.5104 40.5234 45.5364 Vec Value 0 0.0526343 0.105269 0.157903 0.210537 0.263172 0.315806 0.36844 0.421075 0.473709 0.526343 0.578978 0.631612 0.684247 0.736881 0.789515 0.84215 0.894784 0.947418 1.00005

Velocidad y Presin obtenida con [Link] o

27

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Lineas de corriente

IsoValue 0.00242788 0.00729546 0.012163 0.0170306 0.0218982 0.0267658 0.0316333 0.0365009 0.0413685 0.0462361 0.0511036 0.0559712 0.0608388 0.0657064 0.0705739 0.0754415 0.0803091 0.0851767 0.0900442 0.0949118

Funcin de corriente obtenida con [Link] o Observar como la orden square genera el cuadrado [0, 1] [0, 1] y adems etiqueta los a contornos por defecto por 1,2,3,4 empezando por la base y en sentido contrario al del giro de un reloj. coef= . . . da la razn de aspecto de los mdulos de los vectores con o o respecto a la dimensin del dominio. o cmm= . . . pone texto en la grca a La funcin de corriente es solucin del problema o o = u

Usando el resultado anterior calculamos ahora el ujo en una cavidad obteniendo primero el ujo de Stokes como dato inicial. 28 Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

El codigo [Link] resuelve el problema de la cavidad para las ecuaciones de Navier-Stokes tambin con el par P2-P1 para velocidade presin y tomando la solucin del problema de Stokes como dato inicial o o Observacin 8 La expresin (a %b) en C++ equivale a la divisin en o o o mdulo de a entre b y ! es el operador lgico de negacin. La expresin o o o o !(i %10) permite mostrar grca cada vez que i es mltiplo de 10. a u if (!(i %10)) { plot(. . . ) }; La orden nbiso=. . . indica el nmero de isolvalores que se van a represenu tar. Por defecto slo son 20. o mesh Th=square(60,60); fespace Xh(Th,P2); fespace Mh(Th,P1); Xh u2,v2, u1,v1,up1,up2; Mh p,q; func f1=0; func f2=0; real nu=0.001; solve Stokes (u1,u2,p,v1,v2,q) = int2d(Th)( nu*( dx(u1)*dx(v1)+dy(u1)*dy(v1) + dx(u2)*dx(v2)+dy(u2)*dy(v2) ) + p*q*(0.00000001) -p*dx(v1)-p*dy(v2) +dx(u1)*q+dy(u2)*q) - int2d(Th)(f1*v1-f2*v2) +on(1,2,4,u1=0,u2=0)+on(3,u1=1,u2=0); int i=0; real dt=0.1; real alpha=1/dt; problem NS(u1,u2,p,v1,v2,q,solver=UMFPACK,init=i) = int2d(Th)( 29 Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

alpha*(u1*v1+u2*v2) +nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1) +dx(u2)*dx(v2)+dy(u2)*dy(v2) ) + p*q*(0.00000001) -p*dx(v1)- p*dy(v2) +dx(u1)*q+ dy(u2)*q ) + int2d(Th)( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 -f1*v1-f2*v2 ) +on(1,2,4,u1=0,u2=0) +on(3,u1=1,u2=0); for (i=0;i<=100;i++) { up1=u1; up2=u2; NS; //if ( !(i % 20)) plot(coef=0.7,cmm=" [u1,u2] y p ",p,[u1,u2],fill=0,value=1); //plot(coef=.5, fill=1,cmm=" Velocidad [u1,u2] ",[u1,u2]); }; plot(coef=.7,cmm=" Velocidad [u1,u2] ",[u1,u2],wait=1,ps="[Link]"); plot(cmm=" Presion ",p,wait=1,nbiso=50,ps="[Link]");

30

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Velocidad [u1,u2]

Velocidad obtenida con [Link]

Presion

Presin obtenida con [Link] o

31

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Con [Link] obtenemos el ujo en un escalon hacia adelante y adems usamos dos mallas distintas para ver mas cerca la solua cin o real w1=4,w2=26,w=w1+w2,wcut=6,altura=1.5,step=0.75; // // e // y=1.5-> |----------------------------------| // ff| | // y=0.75-> |___________ | d // a | | // b| | // y=0--> |______________________| // x=0 x=4 c x=30 // // border ff(t=0,step){x=0;y=altura-t;label=1;};// x=0, y en [step,altura] border a(t=0,w1){x=t;y=step;label=2;}; // y=step,x en [0,w1] border b(t=0,step){x=w1;y=step-t;label=3;} // x=w1,y en [0,step] border c(t=0,w2){x=w1+t;y=0;label=4;}; //y=0,x en [w1,w1+w2] border d(t=0,altura){x=w;y=t;label=5;}; //x=w, y en [0,altura] border e(t=0,w){x=w-t;y=altura;label=6;}; // y=1.5,x en [0,w] border ccut(t=0,wcut){x=w1+t;y=0;label=4;}; // y=0, x en [w1,w1+w2] border dcut(t=0,altura){x=w1+wcut;y=t;label=5;}; //x=w, y en [0,altura] border ecut(t=0,w1+wcut){x=w1+wcut-t;y=altura;label=6;}; // y=1.5, // x en [0,w] int n=25; mesh th=buildmesh(ff(n)+a(2*n)+b(n)+c(6*n)+d(2*n)+e(12*n)); mesh thcut=buildmesh(ff(n)+a(2*n)+b(n)+ccut(6*n)+dcut(2*n)+ecut(12*n)); plot(th,wait=0,ps="[Link]"); fespace fespace fespace fespace Xh(th,P2); Mh(th,P1); Xhcut(thcut,P2); Mhcut(thcut,P1);

Xh u1,v1,u2,v2,up1,up2; Mh p,q; Xhcut u1cut,v1cut; 32 Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Mhcut pcut; func f1=0; func f2=0; func g=(y-0.75)*(1.5-y)/((2.25/2-0.75)*(1.5-2.25/2));// flujo de entrada real nu=0.0025; solve Stokes (u1,u2,p,v1,v2,q) = int2d(th)(nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1)+ dx(u2)*dx(v2)+dy(u2)*dy(v2)) -p*dx(v1)-p*dy(v2) +dx(u1)*q+dy(u2)*q ) +int2d(th)(-f1*v1-f2*v2) +on(2,u1=0,u2=0) +on(3,u1=0,u2=0) +on(4,u1=0,u2=0) +on(6,u1=0,u2=0) +on(1,u1=g,u2=0); // Dirichlet boundary condition int i=0,nmax=100; real dt=0.1,alpha=1/dt; problem ns(u1,u2,p,v1,v2,q,init=i) = int2d(th)( alpha*(u1*v1+u2*v2)+ nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1) // bilinear part +dx(u2)*dx(v2)+dy(u2)*dy(v2) ) +p*dx(v1)+p*dy(v2)+q*dx(u1)+q*dy(u2) ) // bilinear part +int2d(th)( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 -f1*v1-f2*v2 ) // right hand side +on(2,u1=0,u2=0) +on(3,u1=0,u2=0) +on(4,u1=0,u2=0) +on(6,u1=0,u2=0) +on(1,u1=g,u2=0); // Dirichlet boundary condition

33

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

for (i=0;i<=nmax;i++) { up1=u1; up2=u2; ns; cout << " Valor de i = " << i << " \n"; if( !(i%3) ) { u1cut=u1; plot( u1cut); }; }; plot(cmm= "u1 ", u1,value=0,nbiso=10,wait=1,ps="[Link]"); plot(cmm= "presion ", p,nbiso=10,value=0,ps="[Link]"); u1cut=u1; pcut=p; plot(cmm= "u1cut ", u1cut,value=0,wait=1,ps="[Link]"); plot(cmm= "presioncut ", pcut,value=0,ps="[Link]");

Malla usada con [Link]

34

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

u1

Velocidad obtenida con [Link]

presion

Presin obtenida con [Link] o

35

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

u1cut

Ampliacion de la velocidad obtenida con [Link]

presioncut

Ampliacion de la presin obtenida con [Link] o Se usa la orden convect para calcular la derivada material por el mtodo e de las caracter sticas: t + u = 0, (x, t = 0) = 0 (x).

Ya hemos visto su para Navier-Stokes, otro ejemplo es [Link]

36

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

real pi=4*atan(1.0),dt=0.01; int i; border a(t=0,2*pi){x=cos(t);y=sin(t);}; mesh th=buildmesh(a(300)); plot(th,wait=1); fespace vh(th,P2); vh u,u1,u2,phi; func phi0=exp(-20*((x-0.1)^2+(y-0.3)^2)); u=phi0; plot(u,wait=1,ps="[Link]"); u1=-y; u2=x; for (i=0; i<100;i++) { phi=convect([u1,u2],-dt,u); plot(phi,value=1,ps="[Link]"); u=phi; };

Solucin inicial para [Link] o

37

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

Solucin nal con [Link] o En [Link] se resuelve Navier-Stokes en un canal con un obstculo a real pi=4*atan(1.0),w=20; // // c (3) // y=4 --------------------------------------// | __ | // d (4) | |__| obstaculo | b (2) // | | // y=0 --------------------------------------// x=0 a (1) x=w // // border a(t=0,w){x=t;y=0;label=1;}; border b(t=0,4){x=w;y=t;label=2;}; border c(t=w,0){x=t;y=4;label=3;}; border d(t=4,0){x=0;y=t;label=4;}; int n=20; border circulo(t=0,2*pi){x=3+.5*cos(t);y=2+.2*sin(t);label=5;} mesh th=buildmesh(a(4*n)+b(n)+c(4*n)+d(n)+circulo(-5*n)); plot(th,ps="[Link]"); fespace Xh(th,P2); 38 Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

fespace Mh(th,P1); Xh u1,v1,u2,v2,up1,up2; Mh p,q; int i=0,nmax=500; real nu=0.0025,dt=0.1,alpha=1/dt; func f1=0; func f2=0; func g=-y*(y-4.0)/4.0; solve stokes(u1,u2,p,v1,v2,q) = int2d(th)( nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1) // +dx(u2)*dx(v2)+dy(u2)*dy(v2)) +p*q*0.000000001 -p*dx(v1)-p*dy(v2)+ q*dx(u1)+q*dy(u2) ) // bilinear part +int2d(th)(-f1*v1-f2*v2) //right hand side +on(1,u1=0,u2=0) +on(3,u1=0,u2=0) +on(4,u1=g,u2=0) +on(5,u1=0,u2=0) ; // Dirichlet boundary condition problem ns(u1,u2,p,v1,v2,q,init=i) = int2d(th)( alpha*(u1*v1+u2*v2)+ nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1) // +dx(u2)*dx(v2)+dy(u2)*dy(v2)) +p*q*0.000000001 -p*dx(v1)-p*dy(v2)+ q*dx(u1)+q*dy(u2) ) // bilinear part +int2d(th)(-alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 -f1*v1-f2*v2) //right hand side +on(1,u1=0,u2=0) +on(3,u1=0,u2=0) +on(4,u1=g,u2=0) +on(5,u1=0,u2=0) ; // Dirichlet boundary condition

bilinear part

bilinear part

39

Variaciones sobre manual de Frdric Hecht e e

5 ECUACIONES DE STOKES Y NAVIER-STOKES

for (i=0;i<=nmax;i++) { up1=u1; up2=u2; ns; cout<< "Valor de i = "<< i <<"\n"; plot(u1); //plot(coef=0.5,value=0,cmm= "Velocidad [u1, u2] ", [u1,u2]); }; plot(u1,wait=1,ps="[Link]"); plot(cmm= "presion ", p,value=0,wait=1,ps="[Link]");

Malla usada con [Link]

Velocidad obtenida con [Link]

Velocidad horizontal con [Link] 40 Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

Velocidad vertical con [Link]

presion

Presin obtenida con [Link] o

6.

Descomposicin de dominios o

Con [Link] podemos generar una malla en un c rculo grande menos un c rculo pequeo. Las triangulaciones se generan a la izquierda del sentido n de recorrido de la curva. Por lo tanto, un cambio de signo invierte el sentido de recorrido.

Dos mallas complementarias generadas con [Link]

41

Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

Mallas con hueco generada con [Link] Solapamiento con [Link]: real pi=4*atan(1.0); int n=5; border a(t=0,1){x=t;y=0;}; border a1(t=1,2){x=t;y=0;}; border b(t=0,1){x=2;y=t;}; border c(t=2,0){x=t;y=1;}; border d(t=1,0){x=0;y=t;}; border e(t=0,pi/2){x=cos(t);y=sin(t);}; border e1(t=pi/2,2*pi){x=cos(t);y=sin(t);}; mesh sh=buildmesh(a(n)+a1(n)+b(n)+c(n)+d(n)); mesh SH=buildmesh(e(5*n)+e1(5*n)); plot(sh,SH,wait=1,ps=[Link]);

Dos mallas superpuestas generadas con [Link] Aqu las mallas se superponen pero no son coincidentes (el origen de esta geometr se remonta a los tiempos en los que resolver en un c a rculo o en un cuadrado era mas fcil) En el siguiente ejemplo [Link] se a consiguen dos mallas coincidentes sobre una interfaz. int n=5; // // Denicion del rectangulo [0,1]x[0,1] // 42 Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

border a(t=0,1){x=t;y=0;}; border b(t=0,1){x=1;y=t;}; border c(t=0,1){x=1-t;y=1;};//Interfaz border d(t=0,1){x=0;y=1-t;}; // // Denicion de un rectangulo sobre [0,1]x[0,1] // border c1(t=0,1){x=t;y=1;};//Interfaz border e(t=0,0.2){x=1;y=1+t;}; border f(t=0,1){x=1-t;y=1.2;}; border g(t=0,0.2){x=0;y=1.2-t;}; mesh th=buildmesh(a(10*n)+b(10*n)+c(10*n)+d(10*n)); mesh TH=buildmesh(e(5*n)+f(10*n)+g(5*n)+c1(10*n)); plot(th,TH,wait=1,ps=[Link]);

Dos mallas coincidentes en la interfaz generadas con [Link] lo mismo, pero mas general, se tiene en [Link] y usando mallas uniformes // // Rectangulo [x0,x1]x[y0,y1] junto con el [x0,x1]x[y1,y2] // real x0=0,x1=1; real y0=0.,y1=0.5,y2=1; int n=10,m=10; mesh Th=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]); mesh th=square(n,m,[x0+(x1-x0)*x,y1+(y2-y1)*y]); plot(Th,th,wait=1);

43

Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

Dos mallas coincidentes en la interfaz generadas con [Link] En el ejemplo [Link] se resuelve un problema de Poisson mediante el mtodo de descomposicin de dominios de Schwarz con solapamiento: e o int inside = 2; int outside = 1; real pi=4*atan(1.0),velH1,velH1old=0,uH1,error,tol=1e-7; int n=5; border a(t=1,2){x=t;y=0;label=outside;}; border b(t=0,1){x=2;y=t;label=outside;}; border c(t=2,0){x=t;y=1;label=outside;}; border d(t=1,0){x=1-t;y=t;label=inside;}; border e(t=0,pi/2){x=cos(t);y=sin(t);label=inside;}; border e1(t=pi/2,2*pi){x=cos(t);y=sin(t);label=outside;}; mesh th=buildmesh(a(5*n)+b(5*n)+c(10*n)+d(5*n)); mesh TH=buildmesh(e(5*n)+e1(25*n)); plot(th,TH,wait=1); fespace vh(th,P1); fespace VH(TH,P1); vh u=0,v,uold=0; VH U,V,Uold=0; func f=x-y; int i=0; // // Definicion del problema en el circulo // problem PB(U,V,init=i,solver=Cholesky)= int2d(TH)(dx(U)*dx(V)+dy(U)*dy(V)) +int2d(TH)(-f*V) 44 Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

+on(inside,U=u) +on(outside,U=0); // // Definicion del problema en el rectangulo // problem pb(u,v,init=i,solver=Cholesky)= int2d(th)(dx(u)*dx(v)+dy(u)*dy(v)) +int2d(th)(-f*v) +on(inside,u=U) +on(outside,u=0); // //Ciclo de calculo // for (int j=0; j<10;j++) { PB; pb; // // Calculo del error // vh uerr=u-uold; VH Uerr=U-Uold; uH1=int2d(th)(dx(uerr)^2+dy(uerr)^2); uH1=uH1+int2d(TH)(dx(Uerr)^2+dy(Uerr)^2); velH1=sqrt(uH1); cout<< "-----------------------------------------------"<<"\n"; cout<< "---- Iteracion j= "<<j <<"\n"; cout<< "Error H1 velocidad = "<< velH1 <<"\n"; cout<< "Errores consecutivos= "<< velH1old << " "<<velH1<<"\n"; cout<<" diferencia = "<<abs(velH1-velH1old) <<"\n"; cout<< "-----------------------------------------------"<<"\n"; error=abs(velH1-velH1old); velH1old=velH1; if (error<tol) { cout<< "Velocidades H1 = "<< velH1 <<"\n"; break;}; uold=u; Uold=U; plot(U,u,wait=0); 45 Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

cout<< "Iteracion = "<< j <<"\n"; };

Dos mallas solapadas generadas con [Link]

Solucin obtenida con [Link] o En el ejemplo [Link] se resuelve un problema de Poisson mediante un multiplicador para la forzar el salto a cero (aqui la inf sup uniforme no se cumple): int inside = 2; int outside = 1; real pi=4*atan(1.0),tol=1e-6,velH1,velH1old=0,error,uH1; int n=5; border a(t=1,2){x=t;y=0;label=outside;}; border b(t=0,1){x=2;y=t;label=outside;}; border c(t=2,0){x=t;y=1;label=outside;}; border d(t=1,0){x=1-t;y=t;label=inside;}; 46 Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

border e(t=0,1){x=1-t;y=t;label=inside;}; border e1(t=pi/2,2*pi){x=cos(t);y=sin(t);label=outside;}; mesh th=buildmesh(a(5*n)+b(5*n)+c(10*n)+d(5*n)); mesh TH=buildmesh(e(5*n)+e1(25*n)); plot(TH,th,wait=1); fespace vh(th,P1); fespace VH(TH,P1); vh u=0,v,uold=0; VH U,V,Uold=0; vh lambda=0; func f=-x+y; int i=0,nmax=100; // // Definicion del problema en el circulo // problem PB(U,V,init=i,solver=Cholesky)= int2d(TH)(dx(U)*dx(V)+dy(U)*dy(V)) +int2d(TH)(-f*V) +int1d(TH,inside)(-lambda*V) +on(outside,U=0); // // Definicion del problema en el rectangulo // problem pb(u,v,init=i,solver=Cholesky)= int2d(th)(dx(u)*dx(v)+dy(u)*dy(v)) +int2d(th)(-f*v) +int1d(th,inside)(+lambda*v) +on(outside,u=0); // //Ciclo de calculo // i=1; for (int j=0; j<nmax;j++) { PB; pb; // // Calculo del error // vh uerr=u-uold; 47 Variaciones sobre manual de Frdric Hecht e e

DESCOMPOSICION DE DOMINIOS

VH Uerr=U-Uold; uH1=int2d(th,qft=qf2pT)(dx(uerr)^2+dy(uerr)^2); uH1=uH1+int2d(TH,qft=qf2pT)(dx(Uerr)^2+dy(Uerr)^2); velH1=sqrt(uH1); cout<< "-----------------------------------------------"<<"\n"; cout<< "---- Iteracion j= "<<j <<"\n"; cout<< "Error H1 velocidad = "<< velH1 <<"\n"; cout<< "Errores consecutivos= "<< velH1old << " "<<velH1<<"\n"; cout<<" diferencia = "<<abs(velH1-velH1old) <<"\n"; cout<< "-----------------------------------------------"<<"\n"; error=abs(velH1-velH1old); velH1old=velH1; if (error<tol) { cout<< "Velocidades H1 = "<< velH1 <<"\n"; break;}; uold=u; Uold=U; lambda=lambda+(u-U)/2; plot(U,u,wait=0); };

Dos mallas nosolapadas generadas con [Link]

48

Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

Solucin obtenida con [Link] o Observacin 9 Se puede observar como la funcin sobre la frontera lambda o o se dene sobre todo V h. Esto es porque no se pueden denir trazas y se usan los valores de U y u sobre la frontera para actualizar lambda.

7.

Problemas de evolucin en tiempo o

Podemos resolver la ecuacin del calor o t u u = f, u(x, 0) = u0 , u(x, t) = g, n u(x, t) = 0, (x, t) (0, T ) x (x, t) D (0, T ) (x, t) N (0, T )

siendo = N D y usando elementos nitos en espacio y diferencias nitas en tiempo. Para la funcin o w(x, y, t) = sin(x) cos(y)et el siguiente cdigo [Link] calcula w resolviendo o t u u = f u(x, y, 0) = sin(x) cos(y) u(x, y, t) = sin(x) cos(y)e
t

(x, y, t) (0, T ) (x, y) (x, y, t) (0, T )

siendo f = (1 + 2) u y = (0, 1) (0, 1) con T = 3. Archivo [Link]:

49

Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

mesh Th=square(32,32); fespace Vh(Th,P2); Vh u,v,uu,f,g; real dt=0.1, nu=0.001,error; problem dcalor(u,v) = int2d(Th)( u*v+dt*nu*(dx(u)*dx(v) + dy(u)*dy(v))) +int2d(Th)(-uu*v-dt*f*v) + on(1,2,3,4,u=g); real t=0; uu=sin(x)*cos(y)*exp(t); for (int m=0;m<=5/dt;m++) { t=t+dt; //f=(1+2*nu)*sin(x)*cos(y)*exp(t); //g=sin(x)*cos(y)*exp(t); f=1; g=0; dcalor; plot(u,wait=0,value=1); error=sqrt(int2d(Th)((u-sin(x)*cos(y)*exp(t))^2)); uu=u; cout<< "t= "<< t << " L2-Error = " << error << endl; }; plot(u,wait=1,value=1,ps="[Link]");
IsoValue 0.142811 0.428434 0.714057 0.99968 1.2853 1.57093 1.85655 2.14217 2.42779 2.71342 2.99904 3.28466 3.57029 3.85591 4.14153 4.42716 4.71278 4.9984 5.28402 5.56965

Solucion nal obtenida con [Link] Otro ejemplo lo tenemos en la siguiente aplicacin prctica: Buscamos la o a distribucin de temperatura en una placa 3D de seccin rectangular = o o (0, 6) (0, 1). La placa recibe y mantiene una fuente de calor u = u0 desde 50 Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

sus extremos laterales y est rodeada de aire a temperatura ambiente uc . a La temperatura cambia poco con la coordenada z y por lo tanto podemos considerar el problema como 2D. Si ponemos 1 = {y = 0, y = 1} y 2 = {x = 0, y = 6}, entonces se plantea como resolver la ecuacin del calor o t u ( u) = 0, u(x, y, 0) = u0 , con condiciones de contorno n u + (u uc ) = 0, u = u0 (x, t) 1 (0, T ) (x, t) 2 (0, T ) (x, t) (0, T ) x

y en donde el coeciente de difusin toma dos valores, = 0.2 sobre o y = 0.5 y = 2 bajo y = 0.5. Aplicamos el mtodo de Euler impl e cito y tendremos ((un un1 )/dt, v) + (k un , v) + ((un uc ), v) = 0 Archivo [Link]: func u0 =10+90*x/6;//temp. inicial func k = 1.8*(y<0.5)+0.2; real ue = 25;//temp. exterior real alpha=0.25, T=10, dt=0.1; mesh Th=square(30,5,[6*x,y]);//[0,6]x[0,1] plot(Th,wait=0); fespace Vh(Th,P2); Vh u=u0,v,uold; problem thermic(u,v)= int2d(Th)(u*v/dt+k*(dx(u)*dx(v) +dy(u) * dy(v))) + int1d(Th,1,3)(alpha*u*v) - int1d(Th,1,3)(alpha*ue*v) - int2d(Th)(uold*v/dt) + on(2,4,u=u0) ; for(real t=0;t<T;t+=dt){ uold=u; thermic; plot(u,wait=0); }; plot(u,wait=1,value=1,ps="[Link]"); 51 Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

Solucion nal obtenida con [Link] Finalmente, otro ejemplo lo constituye el calor en torno a un perl de acero: Archivo [Link] real S=99; border C(t=0,2*pi) { x=3*cos(t); y=3*sin(t);} border Splus(t=0,1){ x = t; y = 0.17735*sqrt(t)-0.075597*t - 0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4); label=S;} border Sminus(t=1,0){ x =t; y= -(0.17735*sqrt(t)-0.075597*t -0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4)); label=S;} mesh Th= buildmesh(C(50)+Splus(70)+Sminus(70)); plot(Th,wait=1,ps="[Link]"); fespace Vh(Th,P2); Vh psi,w; solve potential(psi,w)= int2d(Th)(dx(psi)*dx(w)+dy(psi)*dy(w))+ on(C,psi = y) + on(S,psi=0); plot(psi,fill=1,wait=1,ps="[Link]"); border D(t=0,2){x=1+t;y=0;} // added to have a fine mesh at trail mesh Sh = buildmesh(C(25)+Splus(-90)+Sminus(-90)+D(200)); plot(Sh,wait=1,ps="[Link]"); fespace Wh(Sh,P1); Wh v,vv; 52 Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

int steel=Sh(0.5,0).region, air=Sh(-1,0).region; fespace W0(Sh,P0); W0 k=0.01*(region==air)+0.1*(region==steel); W0 u1=dy(psi)*(region==air), u2=-dx(psi)*(region==air); Wh vold = 120*(region==steel); real dt=0.1, nbT=20; int i; problem thermic(v,vv,init=i,solver=LU)= int2d(Sh)(v*vv/dt+ k*(dx(v) * dx(vv) + dy(v) * dy(vv)) + 10*(u1*dx(v)+u2*dy(v))*vv ) -int2d(Sh)(vold*vv/dt); for(i=0;i<nbT;i++){ v=vold; thermic; plot(v, fill=1); } plot(v,wait=1,fill=1,value=1,ps="[Link]");

Red inicial obtenida con [Link]

53

Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

Potencial en torno obtenido con [Link]

Malla total para [Link]

54

Variaciones sobre manual de Frdric Hecht e e

PROBLEMAS DE EVOLUCION EN TIEMPO

IsoValue -7.65279 -1.43095 2.71694 6.86483 11.0127 15.1606 19.3085 23.4564 27.6043 31.7522 35.9001 40.048 44.1959 48.3437 52.4916 56.6395 60.7874 64.9353 69.0832 79.4529

Calor obtenido con [Link] Otro problema es la ecuacion del calor con una turbina border border border border border a(t=pi/2,2*pi){x=1.5*cos(t);y=1.5*sin(t);label=1;};//exterior b(t=1.5,4){x=t;y=0;label=1;};//exterior c(t=0,1.5){x=4;y=t;label=2;};//salida turbina d(t=4,0){x=t;y=1.5;label=1;};//exterior e(t=2*pi,0){x=0.25*cos(t);y=0.25*sin(t);label=3;};//Interior //de turbina mesh Th= buildmesh (a(20)+b(10)+c(6)+d(10)+e(20)); plot(Th,wait=1,ps="[Link]"); fespace Vh(Th,P1); Vh u, v,uold,f,g,s,ue=1; func u0 =0; func k=1; real N=100; real T=5*pi, dt=T/N; uold=u0; problem turbina(u,v)= int2d(Th)(u*v/dt +k*(dx(u)*dx(v)+dy(u)*dy(v)) ) - int2d(Th)(uold*v/dt) +on(1,u=ue)+ on(3,u=g); for(real tt=0;tt<T;tt+=dt) { g=abs(sin(tt+pi/2)); Variaciones sobre manual de Frdric Hecht e e

55

8 FREEFEM 3D

turbina; uold=u; plot(u,wait=0,value=0); }; plot(u,wait=1,fill=0,value=1,ps="[Link]");

Red inicial obtenida con [Link]

IsoValue 0.960811 0.962821 0.964831 0.96684 0.96885 0.97086 0.972869 0.974879 0.976889 0.978898 0.980908 0.982918 0.984927 0.986937 0.988947 0.990956 0.992966 0.994976 0.996985 0.998995

Solucion obtenida con [Link]

8.

FreeFem 3d

Con este primer ejemplo vemos como generar una malla 3D a partir de una 2D sobre un triangulo con una sola capa: Algoritmo [Link]: 56 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

load "msh3" load "medit" border a(t=0,1){x=t;y=0;}; border b(t=0,1){x=1-t;y=t;}; border c(t=0,1){x=0;y=1-t;}; int nn=1; mesh Th2=buildmesh(a(nn)+b(nn)+c(nn)); int ncapas=1; real zmin=0,zmax=2; mesh3 Th=buildlayers(Th2,ncapas,zbound=[zmin,zmax]); medit(" Triangulo1capa ", Th); savemesh(Th,"[Link]");

Triangulacin bsica de un tringulo con [Link] y una sla capa o a a o Se puede obervar que la capa est constituida por a 6 vrtices e 8 tringulos a 3 tetraedros los parmetros zmin=0 y zmax=2 sirven para indicar la altura de la capa. a Adems, el archivo [Link] contiene toda la descripcin de la malla a o 57 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

MeshVersionFormatted 1 Dimension 3 Vertices 6 0 1 0 3 0 1 1 3 0 0 0 3 0 0 1 3 1 0 0 2 1 0 1 2 Tetrahedra 3 1 6 3 5 0 1 4 3 6 0 1 6 2 4 0 Triangles 8 2 4 6 0 5 3 1 1 1 3 4 2 4 2 1 2 5 1 6 3 2 6 1 3 3 5 6 4 6 4 3 4 End

58

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Otra vista donde se observan los tres tetraedros con [Link] Ahora se pueden ver las numeraciones de las caras

59

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Numeracion de caras

60

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Triangulacin bsica de un tringulo con [Link] y dos capas o a a Cubo con slo una capa y con la particin en triangulos de la base lo mas o o simple posible: Algoritmo [Link]: // build the mesh // --------------load "msh3" load "medit" int nn=1; mesh Th2=square(nn,nn); int ncapas=1; real zmin=0,zmax=2; mesh3 Th=buildlayers(Th2,ncapas,zbound=[zmin,zmax]); medit(" Sol ", Th,order=1);// to see the sol in P1

61

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Triangulacin bsica de un cubo con [Link] o a Se puede obervar que la capa est constituida por a 8 vrtices e 12 tringulos a 6 tetraedros

Triangulacin de un cubo con [Link] y una capa o

62

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Triangulacin de un cubo con [Link] y dos capas o En total se tienen ahora 12 vrtices e 20 tringulos a 12 tetraedros Un cubo ms grande se muestra en

63

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Triangulacin de un cubo con [Link] y siete capas. o

Uso de Tetgen
Otra forma de generar mallas 3D es mediante la funcion tetgen. Para eso hay que construir primero las supercies laterales del volumen 3D que se quiere partir en tetraedros y triangularlas de manera conforme con respecto a las interfaces. A continuacin se crea el volumen y la particin. o o Cdigo [Link]: o load "msh3" load "tetgen" load "medit" // // Rectangulo [1,2]x[0,2*pi]: tapas laterales en la direccion z // en z=0 y z=1.5 // real x0,x1,y0,y1; x0=1.; x1=2.; y0=0.; y1=2*pi; mesh Thsq1 = square(5,15,[x0+(x1-x0)*x,y0+(y1-y0)*y]); func ZZ1min = 0; func ZZ1max = 1.5; func XX1 = x; func YY1 = y; 64 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

mesh3 Th31h = movemesh2D3Dsurf(Thsq1,transfo=[XX1,YY1,ZZ1max]); mesh3 Th31b = movemesh2D3Dsurf(Thsq1,transfo=[XX1,YY1,ZZ1min]); // // Rectangulo [1,2]x[0,1.5]: tapas laterales en la direccion y // en y=0 e y=2*pi // x0=1.; x1=2.; y0=0.; y1=1.5; mesh Thsq2 = square(5,8,[x0+(x1-x0)*x,y0+(y1-y0)*y]); func ZZ2 = y; func XX2 = x; func YY2min = 0.; func YY2max = 2*pi; mesh3 Th32h = movemesh2D3Dsurf(Thsq2,transfo=[XX2,YY2max,ZZ2]); mesh3 Th32b = movemesh2D3Dsurf(Thsq2,transfo=[XX2,YY2min,ZZ2]); // // Rectangulo [0,2*pi]x[0,1.5]: tapas laterales en la direccion x // en x=1 x=2 // x0=0.; x1=2*pi; y0=0.; y1=1.5; mesh Thsq3 = square(15,8,[x0+(x1-x0)*x,y0+(y1-y0)*y]); func XX3min = 1.; func XX3max = 2.; func YY3 = x; func ZZ3 = y; mesh3 Th33h = movemesh2D3Dsurf(Thsq3,transfo=[XX3max,YY3,ZZ3]); mesh3 Th33b = movemesh2D3Dsurf(Thsq3,transfo=[XX3min,YY3,ZZ3]); // // Se construye la malla de la superficie del cubo sumando las mallas // mesh3 Th33 = Th31h+Th31b+Th32h+Th32b+Th33h+Th33b; // // Se graba la malla // savemesh(Th33,"[Link]"); // // Se construye finalmente la malla con tetraedros de un cubo paralelo // al eje y usando TetGen // real[int] domain =[1.5,pi,0.75,145,0.0025]; mesh3 Thfinal =tetg(Th33,switch="paAAQY",nbofregions=1,regionlist=domain); 65 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

savemesh(Thfinal,"[Link]"); medit("mallafinal",Thfinal); // // Construimos ahora una malla de la mitad de una placa cilindrica // de radio interior 1., radio exterior 2 y altura 1.5 func mv2x = x*cos(y); func mv2y = x*sin(y); func mv2z = z; mesh3 Thmv2 = movemesh3D(Thfinal, transfo=[mv2x,mv2y,mv2z]); savemesh(Thmv2,"[Link]"); medit("mallafinalmodificada",Thmv2);

Particin de un cubo con [Link] usando tetgen o

Placa cilindrica usando tetgen 66 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Uso de buildlayermesh
Se obtiene una malla 3D por capas a partir de una 2D. Tenemos 3D = 2D [zmin, zmax] donde zmin y zmax son funciones denidas sobre 2D que determinan las capas superior e inferior del dominio 3D . En general, los elementos que se usan son prismas. Los argumentos principales son una malla 2D y el nmero de capas M. u Observacin 10 Cada supercie queda numerada por defecto con el o mismo nmero que la arista sobre la que descansa. Siendo 0 la supercie de u base que se extiende por capas y 1 la ultima capa. Siendo la regin de base o el 0, entonces rup=[0,2] pone la cara superior con la etiqueta 2 rdlow=[0,1] pone la cara inferior con la etiqueta 1 rmid=[1,1,2,1,3,1,4,1] pone las caras laterales, con etiquetas 1,2,3,4, todas con con la etiqueta 1. Cdigo [Link]: o load "msh3" load "medit" border C0(t=0,2*pi){ x=7*cos(t); y=10*sin(t);}; border C1(t=0,2*pi){ x=3+2*cos(t); y=0.5*sin(t);}; mesh Th=buildmesh(C0(100)+C1(-20)); func zmin=-1+sin(x*y)/3; func zmax=2; int MaxLayer=10; mesh3 Th3=buildlayers(Th, MaxLayer,zbound=[zmin,zmax]); savemesh(Th3,"[Link]"); medit("regionconbujero",Th3);

67

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link]

68

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link]

69

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link] Mismo ejemplo con distintos mximos y m a nimos para la malla 3D load "msh3" load "medit" border C0(t=0,2*pi){ x=7*cos(t); y=10*sin(t);}; border C1(t=0,2*pi){ x=3+2*cos(t); y=0.5*sin(t);}; mesh Th=buildmesh(C0(100)+C1(-20)); func zmin=10*((x/10)^2+(y/5)^2-1.2); func zmax=1; int MaxLayer=10; mesh3 Th3=buildlayers(Th, MaxLayer,zbound=[zmin,zmax]); savemesh(Th3,"[Link]"); medit("regionconbujero",Th3); 70 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link]

71

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link]

72

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link]

73

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link] Mismo ejemplo con distintos mximos y m a nimos para la malla 3D load "msh3" load "medit" real L=6; verbosity=0; border aa(t=0,1){x=t; y=0 ;}; border bb(t=0,14){x=1+t; y= - 0.1*t ;}; border cc(t=-1.4,L){x=15; y=t ;}; border dd(t=15,0){x= t ; y = L;}; border ee(t=L,0.5){ x=0; y=t ;}; border ff(t=0.5,0){ x=0; y=t ;}; border C1(t=0,2*pi){ x=3+0.5*cos(t); y=3+0.5*sin(t);}; 74 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

int n=6; mesh Th2d=buildmesh(aa(n)+bb(9*n) + cc(4*n) + dd(10*n)+ee(6*n) + ff(n)+C1(-6*n)); //plot(Th2d,wait=1); func zmin=-1.5*(x<3)+(+1.5-x)*(x>3)*(x<7)-5.5*(x>7); func zmax=1; int MaxLayer=10; mesh3 Th3=buildlayers(Th, MaxLayer,zbound=[zmin,zmax]); savemesh(Th3,"[Link]"); medit("regionconbujero",Th3);

Distintas vistas de la region generada con [Link]

75

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link]

76

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link] Otra forma de obtener profundidad y orilla variable es posible: Cdigo [Link] o load "msh3" load "medit" // // Orilla con grosor eps // real eps=0.1; // int nn=10; border cc(t=0,2*pi){x=cos(t);y=sin(t);}; mesh Th2= buildmesh(cc(100)); 77 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

// // zmin define el fondo // func zmin= 2-sqrt(4-(x*x+y*y)); // // zmax define la superficie // func zmax= eps+ 2-sqrt(3.); mesh3 Th=buildlayers(Th2,nn,coef= max((zmax-zmin)/zmax,1./(nn*2)), zbound=[zmin,zmax]); medit("orilla",Th);

78

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Distintas vistas de la region generada con [Link] para un grosor de 0.01 y de 0.1 de la orilla. Se resuelve ahora el Laplaciano en un cubo. Se debe de observar que por defecto se etiquetan las caras laterales como sigue: Cara superior... 0 Cara inferior... 1 Cara laterales... 2,3,4,5 Cdigo [Link]: o // //Construccion de la malla // load "msh3" load "medit" int nn=8; mesh Th2=square(nn,nn); // real zmin=0,zmax=2; mesh3 Th=buildlayers(Th2,nn, zbound=[zmin,zmax]); // func f=1. ; fespace Vh(Th,P23d); Vh u,v; // macro Grad3(u) [dx(u),dy(u),dz(u)] // EOM 79 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

// problem Lap3d(u,v,solver=CG) = int3d(Th)(Grad3(v) *Grad3(u)) - int3d(Th)(f*v) + on(1,u=0); Lap3d; medit(" Laplaciano ", Th, u ,order=1); // para ver la solucion en P1

Distintas vistas de la solucion u = 1 con u = 0 en partes de la frontera. Se resuelve ahora el problema de Stokes en un cubo unidad con los datos del problema de la cavidad. Se muestran cortes transversales de la solucin y para nn = 15 el mtodo directo UMFPACK consume toda la memoria o e y se debe usar GMRES. Cdigo [Link]: o 80 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

load "msh3"; load "medit"; // // Malla base // int nn=15; mesh Th2=square(nn,nn); // // Espacio de elementos finitos para los cortes transversales // fespace Vh2(Th2,P2); Vh2 ux,uz,p2; // // Malla 3D por capas // real zmin=0,zmax=1; mesh3 Th=buildlayers(Th2,nn, zbound=[zmin,zmax]); // // Espacio de elementos finitos para la solucion de Stokes // fespace VVh(Th,[P23d,P23d,P23d,P13d]); VVh [u1,u2,u3,p]; VVh [v1,v2,v3,q]; macro Grad(u) [dx(u),dy(u),dz(u)]// EOM macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM real mu=1,yy; varf vStokes([u1,u2,u3,p],[v1,v2,v3,q]) = int3d(Th)( mu*(Grad(u1)*Grad(v1)+Grad(u2)*Grad(v2)+Grad(u3)*Grad(v3)) - div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p ) +on(1,2,3,4,5,u1=0,u2=0,u3=0) + on(0,u1=1.,u2=0,u3=0) ; cout << "Construccion de la matriz de rigidez " << endl; matrix A=vStokes(VVh,VVh); cout << "le asociamos un resolutor lineal " << endl; set(A,solver=GMRES); 81 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

cout << "Construimos el termino independiente " << endl; real[int] b= vStokes(0,VVh); cout << "Invertimos el sistema lineal ....." << endl; p[]= A^-1 * b; cout << "Hecho" << endl; cout <<" Cortes seccionales para distintos valores de y ....." << endl; for(int i=1;i<10;i++) { yy=i/10.; ux= u1(x,yy,y); uz= u3(x,yy,y); p2= p(x,yy,y); plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1); if (i==5) {plot([ux,uz],p2,cmm=" cut y = "+yy,wait= 1,ps="[Link]");}; }; medit("cavidad3d",Th,p,order=1); medit("cavidad3d",Th,u3,order=1); medit("cavidad3d",Th,u2,order=1); medit("cavidad3d",Th,u1,order=1);
cut y = 0.5

Corte del problema de la cavidad para Stokes con = 1.

82

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Corte del problema de la cavidad para Stokes con = 1.

83

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Corte del problema de la cavidad para Stokes con = 1.

84

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Corte del problema de la cavidad para Stokes con = 1. Otro ejemplo con Navier-Stokes, una geometr extraa y viscosidad anistropa. a n o Cdigo [Link]: o load "msh3" load "medit" verbosity=0; real nux=0.25,nuy=0.25,nuz=0.0025,dt=0.1; real alpha=1./dt; real L=2,w=5; real a=1; verbosity=0; border aa(t=0,1){x=t; y=0 ;}; border bb(t=0,w){x=1+t; y= - 0.1*t ;}; border cc(t=-0.1*w,L){x=w+1; y=t ;}; border dd(t=w+1,0){x= t ; y = L;}; border ee(t=L,0){ x=0; y=t ;}; // 85 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

//border C1(t=0,2*pi){ x=2+0.5*cos(t); y=2+0.5*sin(t);}; // int n=2,i; mesh Th2d=buildmesh(aa(n)+bb(9*n) + cc(4*n) + dd(10*n)+ee(7*n));// +C1(-6*n)); //plot(Th2d,wait=1); func zmin=-1.5*(x<3)+(+1.5-x)*(x>3)*(x<7)-5.5*(x>7); func zmax=1; int MaxLayer=4; mesh3 Th3d=buildlayers(Th2d, MaxLayer,zbound=[zmin,zmax]); fespace VVh(Th3d,[P23d,P23d,P23d,P13d]); fespace MMh(Th3d,P13d); VVh [up1,up2,up3,pp]; VVh [u1,u2,u3,p]; VVh [v1,v2,v3,q]; macro Grad(u) [dx(u),dy(u),dz(u)]// EOM macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM solve Stokes([up1,up2,up3,pp],[v1,v2,v3,q],solver=UMFPACK) = int3d(Th3d)(nux*(Grad(up1)*Grad(v1))+nuy*(Grad(up2)*Grad(v2)) + nuz*(Grad(up3)*Grad(v3)) +a*( -up2*v1+up1*v2) -div(up1,up2,up3)*q - div(v1,v2,v3)*pp + 0.000000001*pp*q ) + on(1,2,4,5,6,7,up1=0,up2=0,up3=0) + on(0,up1=1,up2=0,up3=0) ; varf vNS([u1,u2,u3,p],[v1,v2,v3,q],init=i) = int3d(Th3d)( alpha*(u1*v1+u2*v2+u3*v3) + nux*(Grad(u1)*Grad(v1))+nuy*(Grad(u2)*Grad(v2)) +nuz*(Grad(u3)*Grad(v3)) +a*( -u2*v1+u1*v2) -div(u1,u2,u3)*q - div(v1,v2,v3)*p + 0.000000001*p*q ) + int3d(Th3d)(-alpha*(convect([up1,up2,up3],-dt,up1)*v1 + convect([up1,up2,up3],-dt,up2)*v2 + convect([up1,up2,up3],-dt,up3)*v3 ) ) + on(1,2,4,5,6,7,u1=0,u2=0,u3=0) 86 Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

+ on(0,u1=1,u2=0,u3=0) ; cout << " Construimos A" << endl; matrix A = vNS(VVh,VVh); cout << " Asociamos un solver a A" << endl; set(A,solver=UMFPACK); real t=0; for(i=0;i<21;++i) {t += dt; cout << " Calculamos termino independiente " << endl; real[int] b=vNS(0,VVh); cout << " iteration " << i << " t = " << t << endl; cout << "Invertimos sistema lineal .... " << endl; u1[]= A^-1*b; up1[]=u1[]; }; { medit("tank3d",Th3d,p,order=1); medit("tank3d",Th3d,u3,order=1); medit("tank3d",Th3d,u2,order=1); medit("tank3d",Th3d,u1,order=1); };

87

Variaciones sobre manual de Frdric Hecht e e

8 FREEFEM 3D

Corte Navier-Stokes.

88

Variaciones sobre manual de Frdric Hecht e e