The following problem arises in Domain Decomposition algorithms which use adjacent subdomain solution results to construct a vector of values upon which a robin (or Neumann) condition is imposed in the current subdomain.
For Laplace's equation a single subdomain solution uses the weak form,
$-\int_{\Omega}( \nabla \cdot u ) \cdot ( \nabla \cdot v) - \int_{\partial \Omega} \alpha u \cdot v = - \int_{\Omega} f \cdot v - \int_{\partial \Omega} G \cdot v$
Where $G$ is the vector of values computed in the adjacent subdomain.
I have included a minimal example of my implementation of this where I create a separate vector function space with dimension 6 in order to form the integral in a 6x6 square mesh. However it appears that I can't use a different test function to make the dimension match for the $G \cdot v$ integral. Am I going about this completely wrong? I have read just about all I can read in the docs, and previously asked questions that resemble this, but I have not seen anything conclusive. This seems like a pretty common thing to have to do when the boundary can't be described by a simple Expression() statement as in all the documentation examples; I'm hopeful that someone out there has an answer for me!
Thanks,
Ben
"""
Domain Decomposition for the Poisson problem
-Laplace(u) = 0 on the unit square.
u,v = 0 on the boundary not in interface.
du/dx + pu = g on the interface
dv/dx + pv = g on the interface
"""
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, 'Lagrange', 1)
boundary_parts = \
MeshFunction("uint", mesh, mesh.topology().dim()-1)
def mark_boundaries(boundary_parts):
class RightRobinBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - 1.0)< tol
Gamma_R = RightRobinBoundary()
Gamma_R.mark(boundary_parts, 2)
return boundary_parts
boundary_parts = mark_boundaries(boundary_parts)
u = TrialFunction(V)
v = TestFunction(V)
p = 0.5
g_np = np.array([2.,2.,2.,2.,2.,2.])
Vg = VectorFunctionSpace(mesh, "R", 0, dim=6)
vg = TestFunction(Vg)
g = Function(Vg)
g.vector()[:] = g_np
Define form crashes with UFLException: An Integral without a Domain is now illegal.
a = - inner(nabla_grad(u), nabla_grad(v))dx + inner(pu,v)ds(2)
L = - inner(f,v)dx - inner(g,vg)*ds(2)
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
