I am using the mixed-dimensional branch and I’m building up on the problem I had earlier.
My the code appears to work well when I use MeshView to create a mesh that is on the edge of the source mesh. But as soon as I try to move away from it sharing 3 edges, which I need to do for my code, then it seems fail.
I have attached my code and the output error. Any help would be greatly appreciated.
from dolfin import *
import matplotlib.pyplot as plt
from typing import List, Tuple
set_log_level(LogLevel.TRACE)
# Width of total domain
width = 1.0
# Height of total domain
height = 1.0
# Total resolution
resolution = 20
f = Constant(10)
a_top = 1
b_top = 1
def newton_solver_parameters():
return{"nonlinear_solver": "newton",
"newton_solver": {"linear_solver": "gmres"}}
def plott(x, title):
plt.figure()
im = plot(x)
ax = plt.gca()
cbar = ax.figure.colorbar(im, ax=ax)
cbar.ax.set_ylabel("", rotation=-90, va="bottom")
plt.title(title)
plt.show()
# Define SubDomain equations
class BottomSubDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[1], (0, 0.5))
class MiddleSubDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[1], (0.25, 0.75))
# Define Boundary equations
class ThreeQuartersBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.75)
class MiddleBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.5)
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)
def run_problem_with(mesh_full: Mesh, mesh_mid: Mesh, domains: MeshFunction,
boundaries_full: MeshFunction, boundaries_sub: MeshFunction):
plott(domains, 'domains')
# Create the function spaces
W_a = FunctionSpace(mesh_full, "Lagrange", 1)
W_b = FunctionSpace(mesh_mid, "Lagrange", 1)
W = MixedFunctionSpace(W_a, W_b)
# Define boundary conditions
bcs = [DirichletBC(W.sub_space(0), a_top, boundaries_full, 1),
DirichletBC(W.sub_space(1), b_top, boundaries_sub, 1),
]
u = Function(W)
ua, ub = u.sub(0), u.sub(1)
va, vb = TestFunctions(W)
ua_0 = 1.0
ub_0 = 1.0
ua.vector()[:] = ua_0
ub.vector()[:] = ub_0
dxa = Measure("dx", domain=mesh_full, subdomain_data=domains)
dxb = Measure("dx", domain=mesh_mid)
F_a = inner(nabla_grad(ua), nabla_grad(va)) * dxa
F_b = inner(nabla_grad(ub), nabla_grad(vb)) * dxb - f * vb * dxb
# This one works for both
# F_a_couple = - f * va * dxa
# This one does not work when the b domain is in the middle
F_a_couple = - f * va * dxa(0) - ub * va * dxa(1)
F = F_a + F_b + F_a_couple
solve(F == 0, u, bcs, solver_parameters=newton_solver_parameters())
plott(ua, "ua")
plott(ub, "ub")
def create_subdomain(mesh_full: Mesh, sub_domain: SubDomain, top_boundary: SubDomain) -> Tuple[Mesh, MeshFunction, MeshFunction]:
# Define domains
domains = MeshFunction("size_t", mesh_full, 2)
domains.set_all(0)
# Mark area for sub domain
sub_domain.mark(domains, 1)
# Create mesh for the sub domain
mesh_sub = MeshView.create(domains, 1)
# Create boundary function for sub domain
boundaries_sub = MeshFunction("size_t", mesh_sub, 1)
boundaries_sub.set_all(0)
top_boundary.mark(boundaries_sub, 1)
return mesh_sub, boundaries_sub, domains
# Define full mesh
mesh_full = RectangleMesh(Point(0, 0),
Point(width, height),
resolution,
resolution)
# Define boundaries on full mesh
boundaries_full = MeshFunction("size_t", mesh_full, 1)
boundaries_full.set_all(0)
top_boundary = TopBoundary()
top_boundary.mark(boundaries_full, 1)
# Run the problem with the coupled domain in the middle of the full mesh
mesh_bot, boundaries_bot, domains_bot = create_subdomain(mesh_full, BottomSubDomain(), MiddleBoundary())
run_problem_with(mesh_full, mesh_bot, domains_bot, boundaries_full, boundaries_bot)
# Run the problem with the coupled domain in the middle of the full mesh
mesh_mid, boundaries_mid, domains_mid = create_subdomain(mesh_full, MiddleSubDomain(), ThreeQuartersBoundary())
run_problem_with(mesh_full, mesh_mid, domains_mid, boundaries_full, boundaries_mid)
Based on the error message it doesn’t seem to like the local index, though I have no idea why.
Error Message
PetscDolfinErrorHandler: line ‘723’, function ‘ISLocalToGlobalMappingApply’, file ‘/build/petsc-Y5jkE2/petsc-3.10.5+dfsg1/src/vec/is/utils/isltog.c’,
: error code ‘63’ (Argument out of range), message follows:
Local index 48437696 too large 230 (max) at 0
PetscDolfinErrorHandler: line ‘2151’, function ‘MatSetValuesLocal’, file ‘/build/petsc-Y5jkE2/petsc-3.10.5+dfsg1/src/mat/interface/matrix.c’,
: error code ‘63’ (Argument out of range), message follows:
PETSc error in ‘/home/alex/fenics-md/dolfin/dolfin/la/PETScMatrix.cpp’, ‘MatSetValuesLocal’
PETSc error code ‘63’ (Argument out of range), message follows:
Local index 48437696 too large 230 (max) at 0
Elapsed wall, usr, sys time: 0.000238802, 0, 0 (Assemble cells)
Traceback (most recent call last):
File “/home/alex/.local/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3325, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File “”, line 1, in
runfile(’/home/alex/Dropbox/Programming/Python/Multi Domain/stackTest.py’, wdir=’/home/alex/Dropbox/Programming/Python/Multi Domain’)
File “/home/alex/Pycharm/helpers/pydev/_pydev_bundle/pydev_umd.py”, line 197, in runfile
pydev_imports.execfile(filename, global_vars, local_vars) # execute the script
File “/home/alex/Pycharm/helpers/pydev/_pydev_imps/_pydev_execfile.py”, line 18, in execfile
exec(compile(contents+"\n", file, ‘exec’), glob, loc)
File “/home/alex/Dropbox/Programming/Python/Multi Domain/stackTest.py”, line 143, in
run_problem_with(mesh_full, mesh_mid, domains_mid, boundaries_full, boundaries_mid)
File “/home/alex/Dropbox/Programming/Python/Multi Domain/stackTest.py”, line 99, in run_problem_with
solve(F == 0, u, bcs, solver_parameters=newton_solver_parameters())
File “/home/alex/.local/lib/python3.7/site-packages/dolfin/fem/solving.py”, line 233, in solve
_solve_varproblem(*args, **kwargs)
File “/home/alex/.local/lib/python3.7/site-packages/dolfin/fem/solving.py”, line 298, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /home/alex/fenics-md/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: cca9a4d76bbd15801aaf965ac110d13d1fafa3c2
*** -------------------------------------------------------------------------