fem.dvi

Size: px

Start display at page:

Download "fem.dvi"

ゆみかきせんばる
5 years ago
Views:

1 2.3 FreeFEM FreeFEM++ 4] FreeFEM 6 O. (Pironneau) F. (Hecht, FreeFEM++ FreeFEM++ FreeFEM++-cs lehyaric/ffcs/index.htm FreeFEM++-cs W w-triangulation.edp w-triangulation.msh w-triangulation.edp border y x : W 5 e koyama/w.html 32

2 int n=5; 2.16: W w-triangulation.edp border Gamma1(t = 0, 1) {x = -t -3; y = 4;} border Gamma2(t = 0, 1) {x = -t +4; y = 4;} border Gamma3(t = 0, 1) {x = 3*t -4; y = -6*t + 4;} border Gamma4(t = 0, 1) {x = 3*t +1; y = 6*t - 2;} border Gamma5(t = 0, 1) {x = -2*t -1; y = 4*t;} border Gamma6(t = 0, 1) {x = -2*t +3; y = -4*t + 4;} border Gamma7(t = 0, 1) {x = -t; y = -2*t + 2;} border Gamma8(t = 0, 1) {x = -t+1; y = 2*t;} border Gamma9(t = 0, 1) {x = t-1; y = 2*t-2;} border Gamma10(t = 0, 1) {x = t; y = -2*t;} mesh Th = buildmesh(gamma1(n)+gamma2(n)+gamma3(6*n)+gamma4(6*n) +Gamma5(4*n)+Gamma6(4*n)+Gamma7(2*n)+Gamma8(2*n) +Gamma9(2*n)+Gamma10(2*n)); plot(th, wait=true, ps="w-triangulation.eps"); savemesh(th, "w-triangulation.msh"); 33

3 Th savemesh w-triangulation.msh w-triangulation.msh : x y

4 2.3.5 Poisson (Laplace ) 2.15 W Laplace Δu = 0 in u = 1 on 1, (P ) u = 0 on 2, u n = 0 on 1, 2. Find u V (g) such that ( (Π) u v x x + u ) v dxdy =0 v V. y y V (g) := { } w H 1 () w = g on 1 2, { 1 on 1, (2.50) g := 0 on 2, V := { } v H 1 () v =0on w-laplace.edp n=10 35

5 int n=10; w-laplace.edp border Gamma1(t = 0, 1) {x = -t -3; y = 4;} border Gamma2(t = 0, 1) {x = -t +4; y = 4;} border Gamma3(t = 0, 1) {x = 3*t -4; y = -6*t + 4;} border Gamma4(t = 0, 1) {x = 3*t +1; y = 6*t - 2;} border Gamma5(t = 0, 1) {x = -2*t -1; y = 4*t;} border Gamma6(t = 0, 1) {x = -2*t +3; y = -4*t + 4;} border Gamma7(t = 0, 1) {x = -t; y = -2*t + 2;} border Gamma8(t = 0, 1) {x = -t+1; y = 2*t;} border Gamma9(t = 0, 1) {x = t-1; y = 2*t-2;} border Gamma10(t = 0, 1) {x = t; y = -2*t;} mesh Th = buildmesh(gamma1(n)+gamma2(n)+gamma3(6*n)+gamma4(6*n) +Gamma5(4*n)+Gamma6(4*n)+Gamma7(2*n)+Gamma8(2*n) +Gamma9(2*n)+Gamma10(2*n)); fespace Vh(Th, P1); Vh u, v; solve laplace(u, v) = int2d(th)(dx(u)*dx(v) + dy(u)*dy(v)) + on(gamma1, u=1) + on(gamma2, u=0); plot(u, wait=true, value=true, fill=true, ps="w-laplace.eps"); u (x, t) Δu(x, t) = f(x, t) in (0, T], t u(x, t) = g(x, t) on D (0, T], (P ) u n (x, t) = 0 on N (0, T], u(x, 0) = u 0 (x) in. D N Find u :0,T] H 1 () such that d (u(t), v)+a(u(t), v) = (f(t), v) v V, (Π) dt u(t) = g(t) on D, u(0) = u 0 in. 36

6 IsoValue : w-laplace.edp (u(t), v) := u(x, t)v(x) dx, a(u(t), v) := u(x, t) v(x) dx, V := { v H 1 () v =0on D }. q 1, q 2,..., q N q 1,..., q N 0 q N 0+1,..., q N 0+N 2 N q N +1,..., q N +N 1 D N := N 0 + N 2 N = N + N 1 q i ϕ i ϕ i ϕ i (q j )=δ ij (1 i, j N) 1 V h := span{ϕ i 1 i N }, W h := span{ϕ i 1 i N} 6 6 V h = {v h W h v h =0on D } 37

7 (Π) Find u h :0,T] W h such that d (Π h ) dt (u h(t), v h )+a(u h (t), v h ) = (f(t), v h ) v h V h u h (t) = g h (t) on D, u h (0) = u 0 h in g h u 0 h g u0 (Π h ) g h g h (x, t)= N +N 1 j=n +1 g(q j,t)ϕ j (x) (Π h ) u h N u h (x, t)= c j (t)ϕ j (x)+ j=1 N +N 1 j=n +1 g(q j,t)ϕ j (x) c j (t) (1 j N ) (Π h ) (Π h ) v h ϕ i (1 i N ) (Π h ) N j=1 ] N dc j b ij dt (t)+a +N 1 ijc j (t) = f i (t) a ij := a(ϕ j,ϕ i ) (1 i, j N), b ij := (ϕ j,ϕ i ) (1 i, j N), f i (t) := (f(t), ϕ i ) (1 i N) A := (a ij ) 1 i, j N, B := (b ij ) 1 i, j N, ( f(t) := f i (t) c(t) := (c i (t)) 1 i N N +N 1 j=n +1 j=n +1 ] g b ij t (qj,t)+a ij g(q j,t) ] ) g b ij t (qj,t)+a ij g(q j,t) (Π h ) Find c :0,T] R N such that (S h ) B dc (t)+ac(t) = f(t), dt c(0) = c 0. 1 i N, (1 i N ). 38

8 u 0 h(x) = N j=1 c 0 := ( ) c 0 i c 0 jϕ j (x)+ 1 i N N +N 1 j=n +1 g(q j, 0)ϕ j (x) d/dt Euler τ t n := nτ (n =0, 1, 2...) c(t n ) c n, f n := f(t n ) (S h ) For each n =1, 2..., find c n R N (Sh) τ ( ) such that c n+1 c n B + Ac n+1 = f n+1. τ (S τ h ) (2.51) (B + τa)c n+1 = Bc n + τf n+1 (Π h ) u h (t n ) u n h (Π h) For each n =1, 2...,find u n h W h such that ( u n+1 (Π τ h h) (Π τ h ) u n+1 h v h dxdy + = 0 v h V h ) u n h,v h + a(u n+1 h,v h ) τ = (f(t n+1 ),v h ) v h V h u n+1 h = g h (t n+1 ) on D, u h (0) = u 0 h in τ ( u n+1 h x v h x + un+1 h y ) v h dxdy u n y hv h dxdy + τ f(t n+1 )v h dxdy FreeFEM++ (P ) 2.15 W D N f 0 g (2.50) u 0 0 (Π τ h ) w-heat.edp 2.12 ] c Find :0,T] R N such that c D ] ] ] B O d c(t) A O (S h ) + O O dt c D (t) O I 39 c(t) c D (t) c(0) c D (0) ] ] = = f(t) g(t) c 0 g(0) ], ].

9 g(t) :=(g(q j,t)) N +1 j N +N N N Drichlet (2.51) ] ] ] ] ] B + τa O c n+1 B O c n f n+1 O τi c n+1 = + τ D O O c n D g n+1 g n := g(t n ) int n=10; real T = 40, tau = 0.1; w-heat.edp border Gamma1(t = 0, 1) {x = -t -3; y = 4;} border Gamma2(t = 0, 1) {x = -t +4; y = 4;} border Gamma3(t = 0, 1) {x = 3*t -4; y = -6*t + 4;} border Gamma4(t = 0, 1) {x = 3*t +1; y = 6*t - 2;} border Gamma5(t = 0, 1) {x = -2*t -1; y = 4*t;} border Gamma6(t = 0, 1) {x = -2*t +3; y = -4*t + 4;} border Gamma7(t = 0, 1) {x = -t; y = -2*t + 2;} border Gamma8(t = 0, 1) {x = -t+1; y = 2*t;} border Gamma9(t = 0, 1) {x = t-1; y = 2*t-2;} border Gamma10(t = 0, 1) {x = t; y = -2*t;} mesh Th = buildmesh(gamma1(n)+gamma2(n)+gamma3(6*n)+gamma4(6*n) +Gamma5(4*n)+Gamma6(4*n)+Gamma7(2*n)+Gamma8(2*n) +Gamma9(2*n)+Gamma10(2*n)); fespace Vh(Th, P1); Vh u=0, v, uold; problem heat(u, v) = int2d(th)(u*v + tau*(dx(u)*dx(v) + dy(u)*dy(v))) - int2d(th)(uold*v) + on(gamma1, u=1) + on(gamma2, u=0); for(real t=0; t<t; t+=tau){ uold = u; heat; plot(u, fill=true); } 40

10 2.3.7 Stokes R 2 Stokes 7 Δu + p = f in, div u = 0 in, (P ) u(x, t) = g on D, u n pn = t on N. D N n := (n 1,n 2 ) (P ) f := (f 1,f 2 ), g := (g 1,g 2 ), t := (t 1,t 2 ) u := (u 1,u 2 ) p Find {u, p} V (g) Q such that u : v dx p div v dx = f v dx + t v d (Π) N v V, N q div u dx = 0 q Q. u : v := V (g) := V := Q := 2 u i v i, i=1 {w H 1 () ]2 } w = g on D, {v H 1 () ]2 } v =0on D, { } q L 2 () qdx = W 2.18 N := 2, D f = o, t = o, { ( a ) g = 4 (y 4)(y 5), 0 on 11, o on D (Π) ( u1 v 1 x x + u 1 v 1 y y + u 2 v 2 x x + u ) 2 v 2 dxdy y y ( v1 p x + v ) ( 2 u1 dxdy + q y x + u ) 2 dxdy = 0 (v 1,v 2,q) V Q. y FreeFEM Stokes V Q V h Q h 7 Stokes 41

11 y x : W V h :P2, Q h :P1 V h :P1, Q h :P1 Inf-Sup β >0 such that h inf sup q h div v h dx β. q h Q h \{0} v h V q h Q v h V h \{0} β h 2.14 FreeFEM++ w-stokes.edpqh q qdxdy=0 1 1 solver=lu 42

12 int n=5; real a=10; func ud = -a*0.25*(y-4)*(y-5); w-stokes.edp border Gamma1(t = 0, 1) {x = -1.5*t -3.5; y = 5;} border Gamma2(t = 0, 1) {x = -t +4; y = 4;} border Gamma3(t = 0, 1) {x = 3*t -4; y = -6*t + 4;} border Gamma4(t = 0, 1) {x = 3*t +1; y = 6*t - 2;} border Gamma5(t = 0, 1) {x = -2.5*t -1; y = 5*t;} border Gamma6(t = 0, 1) {x = -2*t +3; y = -4*t + 4;} border Gamma7(t = 0, 1) {x = -t; y = -2*t + 2;} border Gamma8(t = 0, 1) {x = -t+1; y = 2*t;} border Gamma9(t = 0, 1) {x = t-1; y = 2*t-2;} border Gamma10(t = 0, 1) {x = t; y = -2*t;} border Gamma11(t = 0, 1) {x = -5; y = -t+5;} border Gamma12(t = 0, 1) {x = t-5; y = 4;} mesh Th = buildmesh(gamma1(1.5*n)+gamma2(n)+gamma3(6*n)+gamma4(6*n) +Gamma5(5*n)+Gamma6(4*n)+Gamma7(2*n)+Gamma8(2*n) +Gamma9(2*n)+Gamma10(2*n)+Gamma11(n)+Gamma12(n)); fespace Vh(Th,P2); Vh u1,u2,v1,v2; fespace Qh(Th,P1); Qh p,q; solve stokes(u1,u2,p],v1,v2,q],solver=umfpack) = int2d(th)(dx(u1)*dx(v1)+dy(u1)*dy(v1) + dx(u2)*dx(v2)+ dy(u2)*dy(v2) - p*(dx(v1)+dy(v2)) + q*(dx(u1)+dy(u2))) + on(1,3,4,5,6,7,8,9,10,12,u1=0,u2=0) + on(11,u1=ud,u2=0); plot(u1,u2],p,wait=1,ps="w-stokes.eps"); //plot(th, wait=true, ps="w-stokes-th.eps"); //savemesh(th, "w-stokes-th.msh"); 43

13 2.19: w-stokes.edp 44

C ( 2 1.x ) FreeFem++ (2014) 4 FreeFem++ 2 FreeFem++ 5 (2016/6/4,5) index.html F

C ( 2 1.x ) FreeFem++ (2014) 4 FreeFem++ 2 FreeFem++ 5 (2016/6/4,5) index.html F FreeFem++ 2007 9 26, 2016 2 14, 2017 7 16 http://nalab.mind.meiji.ac.jp/~mk/labo/text/welcome-to-freefem/ http://nalab.mind.meiji.ac.jp/~mk/labo/text/welcome-to-freefem.pdf FreeFem++ ( ) FreeFem++ 1 1