Hello,
I am trying to use the mixed-dimensional branch to solve a non-linear problem where, at the surface of a cylinder, I have a different constitutive than at the inside. I am defining my problem as a MixedNonlinearVariationalProblem and using MixedNonlinearVariationalSolver to solve it. The problem seems to be defined fine, but at the first iteration, the solver fails:
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 56 (No support for this operation for this object type).
*** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: d7350469ca8a7a2d3ba3b5efd3cb34a5bea91853
*** -------------------------------------------------------------------------
I did a minimal example, where I just use the same constitutive for the surface and the volume and where I do not add the Lagrange multiplier to constrained that both fields should be the same at the surface.
import dolfin as d
import numpy as np
import mshr
# parameters
ri = 3.
h = 27
E = 250.
nu = 0.3
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
#%% Mesh definition
cylinder = mshr.Cylinder(d.Point(0, 0, 0), d.Point(0, 0, h), ri, ri)
geometry = cylinder
mesh = mshr.generate_mesh(geometry, 40)
marker = d.MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
class Surface(d.SubDomain):
def inside(self,x,on_boundary):
return d.near(np.linalg.norm(x[0:2]), ri, eps = 0.1) and on_boundary
surf = Surface()
surf.mark(marker, 10)
submesh = d.MeshView.create(marker, 10)
#%% Function Spaces
V = d.VectorFunctionSpace(mesh, 'CG', 1)
SV = d.VectorFunctionSpace(submesh, 'CG', 1)
W = d.MixedFunctionSpace(V, SV)
w = d.Function(W)
u, us = w.sub(0), w.sub(1)
(du, dus) = d.TrialFunctions(W)
(delta_u, delta_us) = d.TestFunctions(W)
dV = d.Measure("dx", domain = W.sub_space(0).mesh())
dS = d.Measure("dx", domain = W.sub_space(1).mesh())
#%% Boundary conditions
class BoundaryUP(d.SubDomain):
def inside(self,x,on_boundary):
return d.near(x[2], 27.) and on_boundary
class BoundaryDOWN(d.SubDomain):
def inside(self,x,on_boundary):
return d.near(x[2], 0.) and on_boundary
boundary_markers = d.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
boundary_markers.set_all(0)
bu = BoundaryUP()
bd = BoundaryDOWN()
bu.mark(boundary_markers, 1)
bd.mark(boundary_markers, 2)
# define value for BC
uD = d.Constant((0.0, 0.0, 0.1))
uU = d.Constant((0.0, 0.0, 0.0))
# define BC
bcu = d.DirichletBC(W.sub_space(0), uU, boundary_markers, 1)
bcd = d.DirichletBC(W.sub_space(0), uD, boundary_markers, 2)
bcs = [bcu, bcd]
#%% Form non linear problem
# Volume
dim = len(u)
I = d.Identity(dim) # Identity tensor
F = I + d.grad(u) # Deformation gradient
J = d.det(F)
FiT = d.inv(F.T)
PK1 = mu*F + (lmbda*d.ln(J) - mu)*FiT
gradu = d.grad(delta_u)
G = d.inner(PK1, gradu)*dV
# Surface
dim = len(us)
Is = d.Identity(dim) # Identity tensor
Fs = I + d.grad(us) # Deformation gradient
Js = d.det(Fs)
FiTs = d.inv(Fs.T)
PK1s = mu*Fs + (lmbda*d.ln(Js) - mu)*FiTs
gradus = d.grad(delta_us)
G += d.inner(PK1s, gradus)*dS
#%% Extract blocks and find Jacobian
Gi = d.extract_blocks(G)
dGi = []
for J in Gi:
for wi in w.split():
dg = d.derivative(J, wi)
dGi.append(dg)
u_sol = w.split()
problem = d.MixedNonlinearVariationalProblem(Gi, u_sol, bcs, dGi)
solver = d.MixedNonlinearVariationalSolver(problem)
solver.solve()
I am using docker to run the program. Does anyone has any idea why this fails?
Thank you in advance!
