@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.
