Here is a test case about fluid-structure interaction from FreeFEM. After solving the elasticity displacement of the solid, I have to remap the mesh or DOFs according to the displacement solution on the interaction surface. FreeFEM has a simple function movemesh to do this, but I can’t find a similar function in getfem. Could anyone please tell me how to achieve this in getfem?
Currently because of a bug in the update of mesh dependencies you will need to rebuild your model from scratch with the deformed mesh.
You can create the new deformed mesh with a command like
m_new = gf.Mesh("clone", m_old)
m_new.set_pts(md.interpolation("u", m_old.pts(), m_old).reshape((-1,2)).T)
Alternatively you can setup your problem using an updated Lagrangian description, where the reference (stress free) configuration is defined by a field variable u_0 containing the mesh deformation.
Thanks for your reply! But I am using C++, could you please tell me how to achieve this procedure in C++?
let’s say you have a variable “u” in a model with some values in it:
getfem::model md;
md.add_fem_variable("u", mfu);
gmm::copy(U, md.set_real_variable("u"));
To deform the points of a mesh m_old, put all the points of the mesh in a mesh_trans_inv structure:
getfem::mesh_trans_inv mti(m_old);
for (dal::bv_visitor j(m_old.points_index()); !j.finished(); ++j)
mti.add_point(m_old.points()[j]);
Then interpolate “u” on this points collection and store the result in a vector pts_def with
getfem::base_vector pts_def(2*mti.nb_points());
getfem::ga_interpolation_mti(md, "u", mti, pts_def);
Then create a new copy of the old mesh and deform its nodes with the values from pts_def
getfem::mesh m_new(m_old);
for (dal::bv_visitor j(m_new.points_index()); !j.finished(); ++j) {
(m_new.points()[j])[0] += pts_def[2*j];
(m_new.points()[j])[1] += pts_def[2*j+1];
}
m_new.optimize_structure();
Thank you so much! These codes will help me a lot 
Dear Konstantinos,
I wang to solve the problem shown in the above figure. There is a interface [7] between the beam (solid) and the fluid inside a cavity. The beam and fluid are two regions in the mesh. Firstly I solve the deformation of the beam under its own weight and have a solution “u” of the beam. Now I want to solve the Stokes equation of the fluid with the boundary [1,2,3,7], it seems that I have to regenerate the mesh of the fluid region. Can I regenerate the fluid mesh in getfem directly by Delaunay triangulation or export the coordinates on the regions [1,2,3,7] and regenerate the mesh with gmsh and reimport gmsh-format mesh file as a new mesh objective?
this problem is quite easy to solve in GetFEM:
this code is without applying the fluid pressure as a load back to the beam. But if necessary it is very simple to add the fluid structure interaction load to the beam.
#!/usr/bin/env python
# -*- coding: UTF8 -*-
#
############################################################################
import numpy as np
import getfem as gf
gf.util_trace_level(1)
gf.util_warning_level(1)
π = np.pi
E = 1000.
ν = 0.4
L = 10.
H_liquid = 5.
H_beam = 0.5
beam_load = 10.
μ_liquid = 0.3
v_in = 1.
p_out =0.
NX_liquid = 20
NX_beam = 25
NY_liquid = 10
NY_beam = 3
# regions
XP_RG = 101
XM_RG = 102
YP_RG = 103
YM_RG = 104
#DIR_RG = 105
m1 = gf.Mesh("import", "structured",
f"GT='GT_QK(2,2)';ORG=[0,0];SIZES=[{L},{H_liquid}];NSUBDIV=[{NX_liquid},{NY_liquid}]")
m2 = gf.Mesh("import", "structured",
f"GT='GT_QK(2,2)';ORG=[0,{H_liquid}];SIZES=[{L},{H_beam}];NSUBDIV=[{NX_beam},{NY_beam}]")
for m in (m1,m2):
m.set_region(XP_RG, m.outer_faces_with_direction([1,0], π/6))
m.set_region(XM_RG, m.outer_faces_with_direction([-1,0], π/6))
m.set_region(YP_RG, m.outer_faces_with_direction([0,1], π/6))
m.set_region(YM_RG, m.outer_faces_with_direction([0,-1], π/6))
#m2.region_merge(DIR_RG, XP_RG)
#m2.region_merge(DIR_RG, XM_RG)
mim1 = gf.MeshIm(m1, 5)
mim2 = gf.MeshIm(m2, 5)
mfu1_ = gf.MeshFem(m1, 2)
mfu1_.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))
mfv1_ = gf.MeshFem(m1, 2)
mfv1_.set_classical_fem(2)
mfp1 = gf.MeshFem(m1, 1)
mfp1.set_classical_fem(1)
mfu2_ = gf.MeshFem(m2, 2)
mfu2_.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))
# remove fixed dofs
kept_dofs = list(set(range(mfu1_.nbdof()))
-set(mfu1_.basic_dof_on_region(XP_RG))
-set(mfu1_.basic_dof_on_region(XM_RG))
-set(mfu1_.basic_dof_on_region(YM_RG)))
mfu1 = gf.MeshFem("partial", mfu1_, kept_dofs)
kept_dofs = list(set(range(mfv1_.nbdof()))
-set(mfv1_.basic_dof_on_region(YP_RG))
-set(mfv1_.basic_dof_on_region(YM_RG)))
mfv1 = gf.MeshFem("partial", mfv1_, kept_dofs)
kept_dofs = list(set(range(mfu2_.nbdof()))
-set(mfu2_.basic_dof_on_region(XP_RG))
-set(mfu2_.basic_dof_on_region(XM_RG)))
mfu2 = gf.MeshFem("partial", mfu2_, kept_dofs)
md = gf.Model("real")
md.add_fem_variable("v", mfv1)
md.add_fem_variable("p", mfp1)
md.add_fem_variable("u1", mfu1)
md.add_fem_variable("u2", mfu2)
md.add_filtered_fem_variable("mult", mfu1, YP_RG)
md.add_interpolate_transformation_from_expression("beam", m1, m2, "X")
md.add_initialized_data("K", E/(3*(1-2*ν)))
md.add_initialized_data("G", E/(2*(1+ν)))
md.add_macro("Dev33(A)", "A-Id(2)/3*(Trace(A)+1)")
md.add_macro("F(u)", "Id(2)+Grad(u)")
md.add_macro("J(u)", "Det(F(u))")
md.add_macro("tauD(u)", "G*pow(J(u),-2/3)*Dev33(F(u)*F(u)')")
md.add_macro("tauH(u)", "K*log(J(u))")
md.add_nonlinear_term(mim1, "((tauH(u1)*Id(2)+tauD(u1))*Inv(F(u1)')):Grad(Test_u1)")
md.add_nonlinear_term(mim2, "((tauH(u2)*Id(2)+tauD(u2))*Inv(F(u2)')):Grad(Test_u2)")
md.add_nonlinear_term(mim1, "mult.Test_u1", YP_RG) # "load" on the top of the liquid domain
md.add_nonlinear_term(mim2, f"{beam_load}*Test_u2(2)", YP_RG) # uniform load on the top of the beam
md.add_nonlinear_term(mim1, "(u1-Interpolate(u2,beam)).Test_mult", YP_RG) # bond liquid top surface to beam bottom surface
md.add_initialized_data("mu", μ_liquid)
md.add_initialized_data("v_in", v_in)
md.add_initialized_data("p_out", p_out)
if False: # Flow in undeformed domain
md.add_nonlinear_term(mim1, "mu*Grad(v):Grad(Test_v)-p*Div(Test_v)" # stokes problem
"+Div(v)*Test_p")
md.add_nonlinear_term(mim1, f"100*(v+v_in*Normal*(1-sqr(X(2)/{H_liquid/2}-1))).Test_v", XM_RG)
else: # Flow in deformed domain
md.add_nonlinear_term(mim1, "mu*(Grad(v)*Inv(F(u1))):(Grad(Test_v)*Inv(F(u1)))-p*Trace(Grad(Test_v)*Inv(F(u1)))" # stokes problem
"+Trace(Grad(v)*Inv(F(u1)))*Test_p")
md.add_nonlinear_term(mim1, f"100*(v+v_in*Normalized(Inv(F(u1))*Normal)*(1-sqr(X(2)/{H_liquid/2}-1))).Test_v", XM_RG)
md.add_nonlinear_term(mim1, "100*(p-p_out)*Test_p", XP_RG)
md.solve("noisy", "max_iter", 20, "max_res", 1e-8,
"lsearch", "simplest", "alpha max ratio", 3., "alpha min", 0.02, "alpha mult", 0.6)
mfu1.export_to_vtu("liquid.vtu", mfu1, md.variable("u1"), "Displacements",
mfv1, md.variable("v"), "Velocities",
mfp1, md.variable("p"), "Pressure")
mfu2.export_to_vtu("beam.vtu", mfu2, md.variable("u2"), "Displacements")
1 Like
Thanks for your reply and codes, they are very helpful to me. It seems that the macro “F(u)” takes into account the deformation of the mesh and introduces this effect into the physical term in the model and J(u) is the Jacobian. If my guess is correct, the function “add_interpolate_transformation_from_expression” will keep the displacement constraint and pressure load at the fluid and beam interface. I will try this in C++, but there are still some python codes I don’t know how to achieve in C++. I will be very happy if you write these in C++, but only these python codes are quite good for me. Thanks again!
Moreover, I also have a question: when an unknown quantity is a vector u=[ux, uy], I want to set Dirichlet condition of [uy] on a boundary, how can I do this? I used to set the quantity as two scalar fields as [ux] and [uy], but it is inconvenient to express divergence or related expressions in the model. So far, I still prefer to use scalar fields because I don’t know how to set boundary conditions on a component of a vector field.
this “solution” sounds like starting chemotherapy against a flu. Don’t do that. If “u” is a vector field with 2 components, these can be accessed in GWFL with u(1) and u(2).
Invest some time to study the available examples that use GWFL and learn some good practices.
Look for “modern” post-2014 scripts like the ones in this folder:
and try to understand them in depth. Python or C++ does not matter, the modelling techniques and the math are important, not the programming.
Thank you very much.