Could different Functionspaces be defined on different subdomains in one mesh?

General
1

Hi, FEniCS family. I learned some basics of FEniCS, but need help with a complex problem. I would like to compute some models with the cohesive zone modelling, and I have learned the example of Cohesive zone modelling of debonding and bulk fracture. However, our challenge is that the CZM is only applied on the interface and the matrix of model is continuous which the phase field requires on the matrix, as the following figure. The algorithm of the above example using the functionspace of DG on the whole mesh is not allowed in our model.

image
image1242×300 48 KB

The numerical realization of the localized CZM is usually through contact or cohesion elements. We found that FEniCS does not have a suitable contact algorithm, so we aim to use the cohesion element, which requires cohesive elements with displacement discontinuities at the interface while the displacement function of the whole model is continuous of CG. My guess is that dividing the subdomains of mesh at the matrix and interface (CZM region) with different displacement functions u of CG on the matrix and DG on the interface is a viable method, but I do not know how to realize it. Thank you very much.

1 Like
2

You can create different submesh from original mesh and define the functionspace over them. Refer this

3

Thank you very much! I will try the SubMesh.

4

Now. I have set two different Functions about u_matrix and u_interface as they are DG and CG on matrix and interface SubMesh, but they are independent and need to be connected on the boundary between u_matrix and u_interface. Is there any way to establish the connection and place them in a variational formula? Thank you.

Here I provide a simple example:

from fenics import *    # version: 2019.2.0.64.dev0
import numpy as np
import matplotlib.pyplot as plt

# Generate mesh
mesh = UnitSquareMesh(10, 10)
mesh_dom = MeshFunction('size_t', mesh, mesh.topology().dim())

# SubMesh
class Interface_SubDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.4, 0.6))
mesh_dom.set_all(0)
Interface_SubDomain().mark(mesh_dom, 1)

Matrix_SubMesh = SubMesh(mesh, mesh_dom, 0)
Interface_SubMesh = SubMesh(mesh, mesh_dom, 1)
Matrix_dom = MeshFunction('size_t', Matrix_SubMesh, Matrix_SubMesh.topology().dim())
Interface_dom = MeshFunction('size_t', Interface_SubMesh, Interface_SubMesh.topology().dim())

dx_Interface = Measure("dx", subdomain_data=Interface_dom)
dx_Matrix = Measure("dx", subdomain_data=Matrix_dom)

# Functionspace that 'CG' on matrix and 'DG' on interface
U_matrix = VectorFunctionSpace(Matrix_SubMesh, 'CG', 1)
U_interface = VectorFunctionSpace(Interface_SubMesh, 'DG', 1)

# Material properties
E = 9.6e6
nu = 0.2
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

# strain and stress
def eps(u):
    return sym(grad(u))
def sigma(u):
    return lamda * tr(eps(u)) * Identity(u.geometric_dimension()) + 2 * mu * eps(u)

# Essential boundary condition
def top_boundary(x, on_boundary):
    return on_boundary and near(x[1], 1)
def bottom_boundary(x, on_boundary):
    return on_boundary and near(x[1], 0)
bc_top = DirichletBC(U_matrix, (0.0, 0.5), top_boundary)
bc_bottom = DirichletBC(U_matrix, (0.0, 0.0), bottom_boundary)
bc = [bc_top, bc_bottom]

# Function
u_matrix = Function(U_matrix, name="Displeacement on Matrix")
u_matrix_trial = TrialFunction(U_matrix)
u_matrix_test = TestFunction(U_matrix)
u_interface = Function(U_interface, name="Displeacement on Interface")
u_interface_trial = TrialFunction(U_interface)
u_interface_test = TestFunction(U_interface)

########################
# I need to connect u_matrix and u_interface and make sure they have same values on boundaries between them.
########################

# I know it is a linear problem in this example, but my model is actually a nonlinear problem on the matrix and the CZM algorithm is applied on the interface.
Eu = inner(sigma(u_matrix), eps(u_matrix_test)) * dx_Matrix + inner(sigma(u_interface), eps(u_interface_test)) * dx_Interface
Eu_J = derivative(Eu, u_matrix, u_matrix_trial) 
problem_u = NonlinearVariationalProblem(Eu, u_matrix, bc, Eu_J)
solver_u = NonlinearVariationalSolver(problem_u)
solver_u.solve()

plot(u_matrix)
5

For legacy dolfin. See MeshView and the corresponding paper and examples:

For similar abstractions in DOLFINx, see for instance

6

Thanks for your help, professor.

7

@dokken Hi, I am very grateful for your help and it is useful for me. But I need to use the MixedNonlinearVariationalSolver to solve our problems, and its explanation is very rare, and I can only draw experience from previous topics, such as Mixed-dimensional branch - Problem with Mixed Nonlinear Variational Solver and Mixed-dimensional branch - how to define NonlinearVariationalProblem?. However, I made an error, just like the following example. Did I miss something or were there some other mistakes?

from fenics import *    # version: 2019.2.0.64.dev0
import numpy as np
import matplotlib.pyplot as plt

# Generate mesh
mesh = UnitSquareMesh(10, 10)
mesh_dom = MeshFunction('size_t', mesh, mesh.topology().dim())

# SubMesh
class Interface_SubDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.4, 0.6))
mesh_dom.set_all(0)
Interface_SubDomain().mark(mesh_dom, 1)

Matrix_SubMesh = MeshView.create(mesh_dom, 0)
Interface_SubMesh = MeshView.create(mesh_dom, 1)

# Functionspace that 'CG' on matrix and 'DG' on interface
U_matrix = VectorFunctionSpace(Matrix_SubMesh, 'CG', 1)
U_interface = VectorFunctionSpace(Interface_SubMesh, 'DG', 1)
facet_interface = MeshFunction('size_t', Interface_SubMesh, Interface_SubMesh.topology().dim() - 1)
U_Mix = MixedFunctionSpace(U_matrix, U_interface)

# Boundary conditions
def top_boundary(x, on_boundary):
    return on_boundary and near(x[1], 1)
def bottom_boundary(x, on_boundary):
    return on_boundary and near(x[1], 0)
bc_top = DirichletBC(U_Mix.sub_space(0), (0.0, 0.5), top_boundary)
bc_bottom = DirichletBC(U_Mix.sub_space(0), (0.0, 0.0), bottom_boundary)
bcs = [bc_top, bc_bottom]

dx_m = Measure("dx", subdomain_data=Matrix_SubMesh)
dx_i = Measure("dx", subdomain_data=Interface_SubMesh)
dS_i = Measure("dS", subdomain_data=facet_interface, metadata={"quadrature_degree": 2})

# functions
(u_m_test, u_i_test) = TestFunctions(U_Mix)
(u_i_trial, u_i_trial) = TrialFunctions(U_Mix)
u_mix = Function(U_Mix)
u_m, u_i = u_mix.sub(0), u_mix.sub(1)

# Material properties
E = 9.6e6
nu = 0.2
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

# strain and stress
def eps(u):
    return sym(grad(u))
def sigma(u):
    return lamda * tr(eps(u)) * Identity(u.geometric_dimension()) + 2 * mu * eps(u)

# Phase-field for the matrix
def Omega_func(d=0):
    # Here is the algorithm for the phase-field
    return Constant(1.0)

# Cohesive zone modelling for the interface
def T_czm(delta_u):
    # Here is the algorithm for the CZM
    D, K = 0, 99999
    return (1 - D) * K * delta_u

Eu_matrix_PhaseField = inner(Omega_func() * sigma(u_m), eps(u_m_test)) * dx_m
Eu_interface_CZM = inner(sigma(u_i), eps(u_i_test)) * dx_i + inner(T_czm(jump(u_i)), jump(u_i_test)) * dS_i

Eu = Eu_matrix_PhaseField + Eu_interface_CZM

Jacobian = []
Eu_blocks = extract_blocks(Eu)
for Li in Eu_blocks:
    for ui in u_mix.split():
        block_d = derivative(Li, ui)
        Jacobian.append(block_d)

problem_u = MixedNonlinearVariationalSolver(Eu_blocks, u_mix.split(), bcs, J=Jacobian)
solver_u = MixedNonlinearVariationalSolver(problem_u)

solver_u.solve()

This is its error. I think the error is on Jacobian.

image
image2162×796 55.3 KB

8

Hi,
I was able to reproduce your error. First of all, there is a mistake in the definition of problem_u, I think you meant MixedNonlinearVariationalProblem and not MixedNonlinearVariationalSolver:

problem_u = MixedNonlinearVariationalProblem(Eu_blocks, u_mix.split(), bcs, J=Jacobian)

Also, your Measure objects are not correctly defined. When defining a dx-type Measure, you can give the Mesh corresponding to you integration domain as domain=<Mesh object>, and subdomain_data takes a MeshFunction defined from this mesh, where you can mark specific entities and integrate over a part of your mesh. For your measures, I think that what you want is :

dx_m = Measure("dx", domain=Matrix_SubMesh)
dx_i = Measure("dx", domain=Interface_SubMesh)

This should fix the error you had, but there is still an issue with the jump term in the variational formulation that I wasn’t able to understand. The corresponding Measure should be :

dS_i = Measure("dS", domain=Interface_SubMesh, subdomain_data=facet_interface)

but I don’t quite understand the definition of facet_interface, which is of dimension 0 ?

9

@cdaversin Hi, I am grateful that you pointed out my mistake. Referring to Cohesive zone modelling of debonding and bulk fracture, dS is used to compute the cohesive force through the jumping displacement of elements.

However, the algorithm in the above example is suitable for using cohesive force models globally, but in practical applications, it usually only needs to be added at some places, otherwise it will damage the overall continuity (so that the Phase Field or other algorithms cannot work) and cause errors. My goal is to compress the model of the above example into the interface area of my model, and I think that is one of the relatively convenient methods to achieve local addition of cohesive forces using FENICS.

Returning to the topic, it still cannot work after being fixed as you said. So if I don’t consider the impact of dS as following, the program, which is reduced to a basic elastic stress-strain problem, can work but cannot iterate correctly with lots of nan. I not sure if matrix and interface displacements of u_m and u_i are connected at their boundaries, If not, the ‘u’ value is not transmitted at the boundary between matrix and interface.

Eu_interface_CZM = inner(sigma(u_i), eps(u_i_test)) * dx_i

If I consider ‘dS’, the program will crash, possibly due to the u_i don’t have the value. Here is the complete example code.

from fenics import *    # version: 2019.2.0.64.dev0
import numpy as np
import matplotlib.pyplot as plt

# Generate mesh
mesh = UnitSquareMesh(10, 10)
mesh_dom = MeshFunction('size_t', mesh, mesh.topology().dim())

# SubMesh
class Interface_SubDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (0.4, 0.6))
mesh_dom.set_all(0)
Interface_SubDomain().mark(mesh_dom, 1)

Matrix_SubMesh = MeshView.create(mesh_dom, 0)
Interface_SubMesh = MeshView.create(mesh_dom, 1)

# Functionspace that 'CG' on matrix and 'DG' on interface
U_matrix = VectorFunctionSpace(Matrix_SubMesh, 'CG', 1)
U_interface = VectorFunctionSpace(Interface_SubMesh, 'DG', 1)
facet_interface = MeshFunction('size_t', Interface_SubMesh, Interface_SubMesh.topology().dim() - 1)

U_Mix = MixedFunctionSpace(U_matrix, U_interface)

# Boundary conditions
def top_boundary(x, on_boundary):
    return on_boundary and near(x[1], 1)
def bottom_boundary(x, on_boundary):
    return on_boundary and near(x[1], 0)
bc_top = DirichletBC(U_Mix.sub_space(0), (0.0, 0.1), top_boundary)
bc_bottom = DirichletBC(U_Mix.sub_space(0), (0.0, 0.0), bottom_boundary)
bcs = [bc_top, bc_bottom]

dx_m = Measure("dx", domain=Matrix_SubMesh)
dx_i = Measure("dx", domain=Interface_SubMesh)
dS_i = Measure("dS", domain=Interface_SubMesh, subdomain_data=facet_interface, metadata={"quadrature_degree": 2})

# functions
(u_m_test, u_i_test) = TestFunctions(U_Mix)
(u_i_trial, u_i_trial) = TrialFunctions(U_Mix)
u_mix = Function(U_Mix)
u_m, u_i = u_mix.sub(0), u_mix.sub(1)

# Material properties
E = 9.6e6
nu = 0.2
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

# strain and stress
def eps(u):
    return sym(grad(u))
def sigma(u):
    return lamda * tr(eps(u)) * Identity(u.geometric_dimension()) + 2 * mu * eps(u)

# Phase-field for the matrix
def Omega_func(d=0):
    # Here is the algorithm for the phase-field
    return Constant(1.0)

# Cohesive zone modelling for the interface
def T_czm(delta_u):
    # Here is the algorithm for the CZM
    D, K = 0, 99999
    return (1 - D) * K * delta_u

Eu_matrix_PhaseField = inner(Omega_func() * sigma(u_m), eps(u_m_test)) * dx_m
# Eu_interface_CZM = inner(sigma(u_i), eps(u_i_test)) * dx_i + inner(T_czm(jump(u_i)), jump(u_i_test)) * dS_i
Eu_interface_CZM = inner(sigma(u_i), eps(u_i_test)) * dx_i

Eu = Eu_matrix_PhaseField + Eu_interface_CZM

Jacobian = []
Eu_blocks = extract_blocks(Eu)
for Li in Eu_blocks:
    for ui in u_mix.split():
        block_d = derivative(Li, ui)
        Jacobian.append(block_d)

problem_u = MixedNonlinearVariationalProblem(Eu_blocks, u_mix.split(), bcs, J=Jacobian)
solver_u = MixedNonlinearVariationalSolver(problem_u)

solver_u.solve()

Thanks for your help.
Regards.

10

See here for the FEniCSx version
https://bleyerj.github.io/comet-fenicsx/tours/interfaces/czm_interface_only/czm_interface_only.html

1 Like