[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Yves Renard |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Tue, 28 Jun 2022 09:51:00 -0400 (EDT) |
branch: master
commit 7a26ce81e8f2587314bf8988daef9885bc928bc5
Author: Yves Renard <Yves.Renard@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)