Hi, FEniCS family.
I am trying to solve problems using MixedFunctionSpace. I have tried two methods: Lagrange multiplier and penalty, but there are errors occurred in the mesh mapping. An example is as follows:
from fenics import * # version: 2019.2.0.64.dev0
import numpy as np
import matplotlib.pyplot as plt
# Material properties
E = 9.6e6
nu = 0.2
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
# Penalty stiffness
K_0 = 999
K_pen = 10 * K_0
# 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 = 0
return (1 - D) * K_0 * delta_u
# Generate mesh
mesh = UnitSquareMesh(10, 10)
mesh_dom = MeshFunction('size_t', mesh, mesh.topology().dim())
class Insterface_SubDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[1], (0.5, 1))
mesh_dom.set_all(0)
Insterface_SubDomain().mark(mesh_dom, 1)
Matrix_SubMesh = MeshView.create(mesh_dom, 0)
Interface_SubMesh = MeshView.create(mesh_dom, 1)
U_matrix = VectorFunctionSpace(Matrix_SubMesh, 'CG', 1)
U_interface = VectorFunctionSpace(Interface_SubMesh, 'CG', 1)
# Lagrange multiplier
surfaces = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
class Gamma_boundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0)
Gamma_boundary().mark(surfaces, 1)
Gamma_submesh = MeshView.create(surfaces, 1)
Gamma_space = VectorFunctionSpace(Gamma_submesh, 'CG', 1)
U_mix = MixedFunctionSpace(U_matrix, U_interface, Gamma_space)
# 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(1), (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)
dx_gamma = Measure("dx", domain=Gamma_submesh)
(u_m_test, u_i_test, Gamma_test) = TestFunctions(U_mix)
u_mix = Function(U_mix)
u_m, u_i, Gamma = u_mix.sub(0), u_mix.sub(1), u_mix.sub(2)
# (u_m_test, u_i_test) = TestFunctions(U_mix)
# u_mix = Function(U_mix)
# u_m, u_i = u_mix.sub(0), u_mix.sub(1)
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
Eu = Eu_matrix_PhaseField + Eu_interface_CZM - inner((u_i - u_m), Gamma_test) * dx_gamma - inner((u_i_test - u_m_test), Gamma) * dx_gamma
# Eu = Eu_matrix_PhaseField + Eu_interface_CZM + K_pen * inner((u_i - u_m), (u_i_test - u_m_test)) * dx_gamma
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)
prmu = solver_u.parameters
prmu["newton_solver"]["relative_tolerance"] = 1E-3
prmu["newton_solver"]["absolute_tolerance"] = 1E-3
prmu["newton_solver"]["convergence_criterion"] = "residual"
prmu["newton_solver"]["error_on_nonconvergence"] = True
prmu["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True
prmu["newton_solver"]["linear_solver"] = 'mumps'
prmu["newton_solver"]["lu_solver"]["symmetric"] = False
prmu["newton_solver"]["maximum_iterations"] = 10000
prmu["newton_solver"]["relaxation_parameter"] = 1.0
print("Completion!")
solver_u.solve()
print("Solve Completion!")
Error:
Build mapping between 25 and 8
Build mapping between 25 and 14
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Error in building the mapping (25, 14) :Index not found.
Thank you very much.