verbosity=0; real L=2,H=2; // half length and height of the domain real r=1,gamma=0.5; // diameter of the ellipse, anisotropy real theta=pi/8.,s1=0,s2=0; // theta=angle of rotation of the solid. s=[s1,s2]=translation of the solid real mu=1,rho=0.5; // viscosity and density of the fluid real rho0=1; // density of the solid macro g [0,-1] // gravitational potential real dt=1,T=500; // time step, Final time // definition of the mesh // int Dirichlet=1,Gamma=2; border bottom(t=-L,L){x=t;y=-H;label=Dirichlet;} border right(t=-H,H){x=L;y=t;label=Dirichlet;} border top(t=L,-L){x=t;y=H;label=Dirichlet;} border left(t=H,-H){x=-L;y=t;label=Dirichlet;} border ball(t=0,2.*pi){x=r*(cos(theta)*cos(t)-gamma*sin(theta)*sin(t))+s1;y=r*(gamma*cos(theta)*sin(t)+sin(theta)*cos(t))+s2;label=Gamma;} real dn=10,h=L/dn; // scale of the mesh mesh Th=buildmesh(bottom(dn*2*L)+right(dn*2*H)+top(dn*2*L)+left(2*dn*H)+ball(2*r*pi*dn)); int Ball=Th(s1,s2).region; // label of the region that contains the ball int Fluid=Th(L,H).region; // label of the fluid region real I=int2d(Th,Ball)(x'*x); // Inertia real V=int2d(Th,Ball)(1.); // Volume fespace Vh(Th,[P2,P2]),Qh(Th,P1); // FE spaces for the velocity/pressure real delta=1.e-5; // Penalization of e(u)=0 on the Ball and of the pression // Some useful macros macro s [s1,s2]// macro u [u1,u2]// macro v [v1,v2]// real sqrt2=sqrt(2.); macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/sqrt2,dy(u[1])]// macro div(u) (dx(u[0])+dy(u[1]))// macro wedge(u,v) ((u)[0]*(v)[1]-(u)[1]*(v)[0])// Vh u,v; Qh p,q; problem Stokes(u1,u2,p,v1,v2,q)= int2d(Th)(mu*e(u)'*e(v) + div(u)*q+div(v)*p) +int2d(Th,Ball)(e(u)'*e(v)/delta) -int2d(Th,Fluid)(rho*g'*v) -int2d(Th,Ball)(rho0*g'*v) +on(Dirichlet,u1=0,u2=0); for(real t=0;t