branch: master commit 7a26ce81e8f2587314bf8988daef9885bc928bc5 Author: Yves Renard <yves.ren...@insa-lyon.fr> AuthorDate: Tue Jun 28 15:49:31 2022 +0200
small improvment --- .../python/demo_fluide_structure_interaction.py | 61 +++++++++++++--------- 1 file changed, 36 insertions(+), 25 deletions(-) diff --git a/interface/tests/python/demo_fluide_structure_interaction.py b/interface/tests/python/demo_fluide_structure_interaction.py index 3a5281dc..8c10cefa 100644 --- a/interface/tests/python/demo_fluide_structure_interaction.py +++ b/interface/tests/python/demo_fluide_structure_interaction.py @@ -22,7 +22,7 @@ ############################################################################ """ Incompressible Navier-Stokes fluid in interaction with a ball in a cavity. - Middle point scheme for the fluid, Verlet scheme for the ball (not the + Middle point scheme for the fluid, Verlet's scheme for the ball (not the best but can of course be changed) This program is used to check that python-getfem is working. This is also @@ -34,12 +34,15 @@ import numpy as np import getfem as gf gf.util('trace level', 0) +version = 1 # 1-without cut-fem and Nitsche's method, + # 2-with cut-fem + # Phisical parameters (à ajuster pour être physique ...) nu = 0.002 # cinematic viscosity rho = 1. # fluid density mu = rho * nu # dynamic viscosity g = 9.81 # gravity -in_velocity = 1. # inward velocity at the center bottom +in_velocity = 4. # inward velocity at the center bottom ball_radius = 0.1 ball_mass = 0.032 ball_init_pos = np.array([0., 0.2]) @@ -47,11 +50,11 @@ ball_init_pos = np.array([0., 0.2]) # Discretization parameters dt = 0.01 # Time step T = 40. # Final time -gamma0 = 1E3 # Nitsche's method parameter +gamma0 = 1. # Nitsche's method parameter # Geometry and Mesh of the cavity -W1 = 0.1 # ____________ -W2 = 0.45 # H2 |_ _| +W1 = 0.05 # ____________ +W2 = 0.475 # H2 |_ _| W = W1+2.*W2 # | | H1 = 0.8 # | | H2 = 0.2 # H1 | | @@ -148,23 +151,30 @@ md.add_Dirichlet_condition_with_multipliers(mim, "v", mfv, OUT_RG2, "v_out2") md.add_Dirichlet_condition_with_multipliers(mim, "v", mfv, WALL_RG) # Diffusive terms -md.add_nonlinear_term(mim_out, "(1/dt)*rho*(v-v0).Test_v + mu*Grad_v:Grad_Test_v") +if (version == 2): + mim_t = mim_out; +else: + mim_t = mim +md.add_nonlinear_term(mim_t, "(1/dt)*rho*(v-v0).Test_v + mu*Grad_v:Grad_Test_v") # Nonlinear convective term -md.add_nonlinear_term(mim_out, "0.25*rho*((Grad_v+Grad_v0)*(v+v0)).Test_v") +md.add_nonlinear_term(mim_t, "0.25*rho*((Grad_v+Grad_v0)*(v+v0)).Test_v") # Pressure terms -md.add_nonlinear_term(mim_out, "-p*Div_Test_v - 0.5*Test_p*(Div_v+Div_v0)") +md.add_nonlinear_term(mim_t, "-p*Div_Test_v - 0.5*Test_p*(Div_v+Div_v0)") # Gravity term -md.add_nonlinear_term(mim_out, "-f.Test_v") +md.add_nonlinear_term(mim_t, "-f.Test_v") # Small ghost penalty term -md.add_linear_term(mim, "1E-4*(Grad_v-Interpolate(Grad_v, neighbor_element)):(Grad_Test_v-Interpolate(Grad_Test_v, neighbor_element))", INTERNAL_EDGES) -# md.add_linear_term(mim, "1E-7*(p-Interpolate(p, neighbor_element))*(Test_p-Interpolate(Test_p, neighbor_element))", INTERNAL_EDGES) -# Penalty term on internal dofs -Bv = gf.Spmat('empty',nbdofv) -ibv = md.add_explicit_matrix("v", "v", Bv) -Bp = gf.Spmat('empty',nbdofp) -ibp = md.add_explicit_matrix("p", "p", Bp) +if (version == 2): + md.add_linear_term(mim, "1E-6*(Grad_v-Interpolate(Grad_v, neighbor_element)):(Grad_Test_v-Interpolate(Grad_Test_v, neighbor_element))", INTERNAL_EDGES) + # md.add_linear_term(mim, "1E-7*(p-Interpolate(p, neighbor_element))*(Test_p-Interpolate(Test_p, neighbor_element))", INTERNAL_EDGES) + # Penalty term on internal dofs + Bv = gf.Spmat('empty',nbdofv) + ibv = md.add_explicit_matrix("v", "v", Bv) + Bp = gf.Spmat('empty',nbdofp) + ibp = md.add_explicit_matrix("p", "p", Bp) + # md.add_linear_term(mim_in, "1E-4*(v - ball_v).Test_v") # Nitsche's term on the FS interface -md.add_nonlinear_term(mim_bound, "-(mu*Grad_v-p*Id(meshdim))*Normalized(Grad_ls).Test_v + gamma * (v-ball_v).Test_v - (mu*Grad_Test_v)*Normalized(Grad_ls).(v-ball_v)") +if ((version == 1) or (version == 2)): + md.add_nonlinear_term(mim_bound, "-(mu*Grad_v-p*Id(meshdim))*Normalized(Grad_ls).Test_v + gamma * (v-ball_v).Test_v - (mu*Grad_Test_v)*Normalized(Grad_ls).(v-ball_v)") t = 0 step = 0 @@ -188,14 +198,15 @@ while t < T+1e-8: mim_out.adapt() # Penalization of ball internal dofs with no contribution on the boundary - idofv = np.setdiff1d(np.arange(nbdofv), mfv.dof_from_im(mim_out)) - idofp = np.setdiff1d(np.arange(nbdofp), mfp.dof_from_im(mim_out)) - Bv = gf.Spmat('empty', nbdofv) - for i in idofv: Bv.add(i,i,1.) - md.set_private_matrix(ibv, Bv) - Bp = gf.Spmat('empty', nbdofp) - for i in idofp: Bp.add(i,i,1.) - md.set_private_matrix(ibp, Bp) + if (version == 2): + idofv = np.setdiff1d(np.arange(nbdofv), mfv.dof_from_im(mim_out)) + idofp = np.setdiff1d(np.arange(nbdofp), mfp.dof_from_im(mim_out)) + Bv = gf.Spmat('empty', nbdofv) + for i in idofv: Bv.add(i,i,1.) + md.set_private_matrix(ibv, Bv) + Bp = gf.Spmat('empty', nbdofp) + for i in idofp: Bp.add(i,i,1.) + md.set_private_matrix(ibp, Bp) # Solve # md.solve("noisy", "lsolver", "mumps", "max_res", 1e-8)