Extract traction from one domain to neighbor domain

variational formulation
1

Dear all,

I’m working on a problem similar to FSI with fixed-point iteration. I want to solve the left domain as in figure and extract the traction (not nodal force) from right boundary, then apply the same traction *-1 to the right domain. The deformed structures are also shown in figure.

1
11000×1071 179 KB

I had already implemented the nodal force method after assembly, but it cause another issue at the middle corner when I have 4 domains. Thus, I hope I can directly transform the traction to avoid the corner effect. My MME is as below. I try to use \tau = \sigma(u) n as

fe.inner(sigma(u_L)*normal_L, v_R)*ds_R(0)

but I got error…

import fenics as fe
# Mesh generation
mesh_L = fe.RectangleMesh(fe.Point(0, 0), fe.Point(1, 1), 10, 10)
mesh_R = fe.RectangleMesh(fe.Point(1, 0), fe.Point(2, 1), 10, 10)
# Function Space
V_L = fe.VectorFunctionSpace(mesh_L, 'P', 1)
V_R = fe.VectorFunctionSpace(mesh_R, 'P', 1)

# Elastic parameters
lambda_, mu = fe.Constant(20.20), fe.Constant(90.91)
def epsilon(u):
    return fe.sym(fe.grad(u))
def sigma(eps_u):
    return lambda_*fe.tr(eps_u)*fe.Identity(2) + 2 * mu * eps_u

# Boundary conditions
class Top(fe.SubDomain):
    def inside(self, x, on_boundary):
        return fe.near(x[0], 1.0)
class Bottom(fe.SubDomain):
    def inside(self, x, on_boundary):
        return fe.near(x[0], 0.0)
class Interface(fe.SubDomain):
    def inside(self, x, on_boundary):
        return fe.near(x[1], 1.0)
top = Top()
bottom = Bottom()
interface = Interface()
bound_L = fe.MeshFunction('size_t', mesh_L, mesh_L.geometry().dim()-1, 10)
bound_R = fe.MeshFunction('size_t', mesh_R, mesh_R.geometry().dim()-1, 10)
top.mark(bound_L, 3)
bottom.mark(bound_L, 2)
interface.mark(bound_L, 1)
top.mark(bound_R, 3)
bottom.mark(bound_R, 2)
interface.mark(bound_R, 0)
bc_L = [fe.DirichletBC(V_L, fe.Constant((0, 0)), bound_L, 2),
        fe.DirichletBC(V_L, fe.Constant((0, 0)), bound_L, 1)]
bc_R = fe.DirichletBC(V_R, fe.Constant((0, 0)), bound_R, 2)

# Variational problem for left
u_L = fe.Function(V_L, name="displacement_L")
u_t_L, v_L = fe.TrialFunction(V_L), fe.TestFunction(V_L)
dx_L = fe.Measure("dx", domain=mesh_L)
ds_L = fe.Measure("ds", domain=mesh_L, subdomain_data=bound_L)
normal_L = fe.FacetNormal(mesh_L)
a_L = fe.inner(sigma(epsilon(u_t_L)), epsilon(v_L))*dx_L
L_Left = fe.inner(fe.Constant((0, 1.0)), v_L)*ds_L(3)

# Solve Left first
fe.solve(a_L == L_Left, u_L, bc_L)

# Variational problem for Right
u_R = fe.Function(V_R, name="displacement_R")
u_t_R, v_R = fe.TrialFunction(V_R), fe.TestFunction(V_R)
dx_R = fe.Measure("dx", domain=mesh_R)
ds_R = fe.Measure("ds", domain=mesh_R, subdomain_data=bound_R)
a_R = fe.inner(sigma(epsilon(u_t_R)), epsilon(v_R))*dx_R
n = fe.FacetNormal(mesh_R)
L_Right = fe.inner(fe.Constant((0, 1.0)), v_R)*ds_R(3) - \
    fe.inner(sigma(u_L)*normal_L, v_R)*ds_R(0)
# Solve Right
fe.solve(a_R == L_Right, u_R, bc_R)

UFLException: Trace of tensor with rank != 2 is undefined.

Any suggestions or comments are welcome.

Thanks,
Ming

2

hi, problem comes from L_Right where you should change to sigma(epsilon(u_L))

3

Hi @bleyerj,

Thanks for your help. After I corrected it, I got

Calling FFC just-in-time (JIT) compiler, this may take some time.
------------------- Start compiler output ------------------------
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp: In member function 'virtual void ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46_exterior_facet_integral_main_0::tabulate_tensor(double*, const double* const*, const double*, std::size_t, int) const':
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:104:18: error: redeclaration of 'const double J_c0'
     const double J_c0 = coordinate_dofs[0] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[2] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:93:18: note: 'const double J_c0' previously declared here
     const double J_c0 = coordinate_dofs[0] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[2] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:105:18: error: redeclaration of 'const double J_c1'
     const double J_c1 = coordinate_dofs[0] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[4] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:95:18: note: 'const double J_c1' previously declared here
     const double J_c1 = coordinate_dofs[0] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[4] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:106:18: error: redeclaration of 'const double J_c2'
     const double J_c2 = coordinate_dofs[1] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[3] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:96:18: note: 'const double J_c2' previously declared here
     const double J_c2 = coordinate_dofs[1] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[3] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:107:18: error: redeclaration of 'const double J_c3'
     const double J_c3 = coordinate_dofs[1] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[5] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~
/tmp/tmp3vtcf_xw/ffc_form_cbeed0341edf6ac3f020691fcc4ccad59675cb46.cpp:94:18: note: 'const double J_c3' previously declared here
     const double J_c3 = coordinate_dofs[1] * FE13_C0_D01_F_Q1[0][0][0] + coordinate_dofs[5] * FE13_C0_D01_F_Q1[0][0][1];
                  ^~~~

from

     64 # Solve Right
     65 fe.solve(a_R == L_Right, u_R, bc_R)
4

Hi @nate @cdaversin @kamensky, any suggestions for this question?
Thank you,
Ming