Hey,
I stumbled into the following error while trying to assemble a mixed variational system:
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /build/dolfin-GfTndI/dolfin-2019.2.0~git20201207.b495043/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
The system I’m trying to simulate follows the following equations:
\begin{aligned}
Eq0: &&-\nabla^2 u0 &= &1\\
Eq1: &&-\nabla^2 u0 -\nabla^2 u1 &= &1\\
\end{aligned}
with a Dirichlet Boundary at the left with u0_{left}=1, u1_{left}=1. On the other boundaries there should be no flow going outside: \frac{\delta u_0}{\delta \overrightarrow{n}}=0 and \frac{\delta u_1}{\delta \overrightarrow{n}}=0.
Hereby u1 is only defined on the yellow subdomain, whereas u0 is defined on the whole domain (purple and yellow):
Defining the TestFunctions v0 and v1 for equation 0 and 1 results in:
\begin{aligned}
Eq0: && \int_{\Omega_0} \nabla u0 \nabla v0 dx&=& &\int_{\Omega_0} v0 dx + \underbrace{\int_{\delta\Omega_0} \frac{\delta u0}{\delta \overrightarrow{n}} v0 ds}_{=0} \\
Eq1: &&\int_{\Omega_1}\nabla u0\nabla v1 dx +\int_{\Omega_1}\nabla u1\nabla v1 dx &= &&\int_{\Omega_1} v1 dx + \underbrace{\int_{\delta\Omega_1} \frac{\delta u1}{\delta \overrightarrow{n}} v1 ds}_{=0} + \underbrace{\int_{\delta\Omega_0} \frac{\delta u0}{\delta \overrightarrow{n}} v1 ds}_{\neq0}\\
\end{aligned}
Therefore my code is the following:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(64, 64)
# Create Subdomain
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return x[1] <= 0.75 + tol and x[1] >= 0.25 - tol
materials = MeshFunction("size_t", mesh, 2)
materials.set_all(0)
subdomain = Omega_0()
subdomain.mark(materials, 1)
mesh_subdomain = MeshView.create(materials, 1)
# Introduce FunctionSpaces for different variables
v0 = FunctionSpace(mesh, 'P', 1)
v1 = FunctionSpace(mesh_subdomain, 'P', 1)
V = MixedFunctionSpace(v0, v1)
## Init the functions
u0, u1 = TrialFunctions(V)
v0, v1 = TestFunctions(V)
# dx and ds for the integrals
dx_0 = Measure("dx", domain = V.sub_space(0).mesh())
dx_1 = Measure("dx", domain = V.sub_space(1).mesh())
ds_1 = Measure("ds", domain = V.sub_space(1).mesh())
# Boundary conditions
def boundary_Left(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 0., tol)
bc_0 = DirichletBC(V.sub_space(0), Constant(1.), boundary_Left)
bc_1 = DirichletBC(V.sub_space(1), Constant(1.), boundary_Left)
bcs = [bc_0, bc_1]
n_1 = FacetNormal(V.sub_space(1).mesh())
### Set up the variational problem
a_0 = dot(grad(u0), grad(v0)) * dx_0
a_1 = dot(grad(u1), grad(v1)) * dx_1 + \
dot(grad(u0), grad(v1)) * dx_1 - \
dot(n_1, grad(u0)) * v1 * ds_1
a = a_0 + a_1
L_0 = Constant(1.) * v0 * dx_0
L_1 = Constant(1.) * v1 * dx_1
L = L_0 + L_1
du_ = Function(V)
system = assemble_mixed_system(a==L, du_, bcs)
At the last line (assemble_mixed_system) I get said error.
Confusingly if I use a differnt Subdomain Omega_1, then assembling the system works.
# Create Subdomain
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return x[1] <= 0.75 + tol
Can anyone point me to what I am missing?
Thanks in advance!
Emanuel