load "Curvature" load "isoline" load "iovtk" verbosity=0; int[int] orderOut = [1,1,1]; string outputFileName = "output2/"; real muf=0.001,c1=1.e3,penal=1.e-12, g=0, fx=0,fy=0,fm=0,fmax=-1.e9, rhof=1000,rhos=1000; int NN=1000, saveEach=10; real t=0., T=1., dt=T/NN, mum=1.e3,lambda=1.e2; real h=1.e-5, H=3.*h, L=5.*h, l=16.*h, L1 = L-2.*h; real HH=30*h, LL=30*h; int remesh=2, m=4, C=5; // Fluid Region border a1(t=-LL,LL) {x=t; y=-HH ;label=1;}; border a2(t=-HH,HH) {x=LL; y=t ;label=1;}; border a3(t=LL,-LL) {x=t; y=HH ;label=1;}; border a4(t=HH,-HH) {x=-LL; y=t ;label=1;}; func fluidbound = a1(C*m)+a2(C*m)+a3(C*m)+a4(C*m); //plot(fluidbound, dim=2, wait=1); real xo=0, yo=5*H/8, r=H/8; // Solid Region border b1(t=-L1/2., L1/2.) {x=t; y=0. ;label=5;}; border b2(t=0,-l) {x=L1/2.; y=t ;label=5;}; border b3(t=L1/2., L1/2. + h) {x=t; y=-l ;label=5;}; border b4(t=-l, H) {x=L1/2.+h; y=t ;label=5;}; border b5(t=L/2.,-L/2.) {x=t; y=H ;label=5;}; border b6(t=H, -l) {x=-L/2; y=t;label=5;}; border b7(t=-L/2,-L1/2.) {x=t; y=-l;label=5;}; border b8(t=-l,0) {x=-L1/2.; y=t;label=5;}; border cc(t=0, 2*pi) {x=r*cos(t)+xo; y=r*sin(t)+yo; label=6;}; func solidbound = b1(L1/h*m)+b2(l/h*m)+b3(m)+b4((H+l)/h*m)+b5(2*L/h*m)+b6((H+l)/h*m)+b7(m)+b8(l/h*m)+cc(-8*pi*r/h*m); plot(solidbound, dim=2, wait=1); mesh Th = buildmesh(fluidbound + solidbound); plot(Th, wait=1); mesh Thold =Th; real xoo=-L/4.,yoo=H/2; int fluid=Th(0.9*LL, 0.9*HH).region, solid=Th(xoo, yoo).region; cout<<"region number: "<