# Difficulties and errors when trying to apply Robin Boundary Conditions with data from another domain (Eigenvalue Problem)

Hi. I’m currently trying to implement the Schwarz Alternating Algorithm for an EigenValue problem. The algorithm works in a domain that has an overlap, such as:

We start by solving the eigenvalue problem on each domain separately, which has the bilinear form:

a_n(u,v) = \int_{\Omega_n} \nabla u_n \cdot \nabla v d\Omega_n -\omega^2 \int_{\Omega_n} u_n v d\Omega

With homogeneous boundary conditions in the non-overlapping boundary, and free conditions on the overlapping boundaries:

u_n = 0 \text{ on } \partial \Omega \setminus\Gamma_n

After solving for both Subdomains, we start an iterative algorithm that applies a Robin Boundary Condition in \Gamma _n , which is a function of the other subdomain. The bilinear form of this step is given by:

a_n (u_n, v) = \int_{\Omega_n} \nabla u_n \cdot \nabla v d\Omega_n -\omega^2 \int_{\Omega_n} u_n v d\Omega - \int_ {\Gamma _n} (1 - g_n) (\nabla u_n \cdot \vec{n}) v d \Gamma_n + \int_ {\Gamma_n} f_n u_n v d {\Gamma _n}

The Schwarz Alternating Algorithm is defined by:

1. Do f_1 = \nabla u_2 \cdot n_2 and g_1 = u_2 |_{\Gamma _1}
2. Find u_1 ^{i+1} that is solution to a_1(u_1^{i+1},v) = 0
3. Do f_2 = \nabla u_1 \cdot n_1 and g_2 = u_1 |_{\Gamma _2}
4. Find u_2 ^{i+1} that is solution to a_1(u_2^{i+1},v) = 0
5. Until || u_1 ^{i+1} - u_2 ^{i+1} ||_{\Omega_{12}} \leq \epsilon

I’m currently having trouble in applying this Robin Boundary Condition, as I have to access information in the other Subdomain. Currently, I have this MWE to show:

from fenics import *
from dolfin import *
from petsc4py import PETSc
## Geometries and Meshes Definition
nx = 24
ny = 24

o = 0.7 # Total Overlapping Area
b = 1.5*(1-o)/2 # Total non-Overlapping Area

omega1 = RectangleMesh(Point(0.0,0.0), Point(b+o*1.5,1), nx, ny, "right/left") # First SubDomain
omega2 = RectangleMesh(Point(b,0.0), Point(1.5,1), nx, ny, "right/left") # Second SubDomain

# Function Spaces
V1 = FunctionSpace(omega1,"Lagrange",1)
V2 = FunctionSpace(omega2,"Lagrange",1)

# Defining Boundary Domains
# Omega 1
class DirichletBoundaryOmega1(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] < b+1.5*o -DOLFIN_EPS)
class RobinBoundaryOmega1(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] > b+1.5*o -DOLFIN_EPS)
# Boundary Markers
omega1_markers = MeshFunction('size_t', omega1, omega1.topology().dim()-1,0)
omega1_markers.set_all(9999)
bc_D1 = DirichletBoundaryOmega1()
bc_R1 = RobinBoundaryOmega1()
bc_D1.mark(omega1_markers, 0)
bc_R1.mark(omega1_markers, 1)
# Omega 2
class DirichletBoundaryOmega2(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] > b+DOLFIN_EPS)
class RobinBoundaryOmega2(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (x[0] < b+DOLFIN_EPS)
# Boundary Markers
omega2_markers = MeshFunction('size_t', omega2, omega2.topology().dim()-1,0)
omega2_markers.set_all(9999)
bc_D2 = DirichletBoundaryOmega2()
bc_D2.mark(omega2_markers, 0)
bc_R2 = RobinBoundaryOmega2()
bc_R2.mark(omega2_markers, 1)

## Assemblies

# Trial/Test Functions
u1 = TrialFunction(V1)
u2 = TrialFunction(V2)

v1 = TestFunction(V1)
v2 = TestFunction(V2)

# Bilinear Forms
m1 = u1*v1*dx

m2 = u2*v2*dx

# Matrices
# Stiffness
K1 = PETScMatrix()
K2 = PETScMatrix()

assemble(k1, tensor=K1)
assemble(k2, tensor=K2)

# Applying Boundary Conditions
bcs1 = DirichletBC(V1, 0.0, bc_D1)
bcs1.apply(K1)

bcs2 = DirichletBC(V2, 0.0, bc_D2)
bcs2.apply(K2)

# Mass
M1 = PETScMatrix()
M2 = PETScMatrix()

assemble(m1, tensor=M1)
assemble(m2, tensor=M2)

# EigenSolver
# Create the eigensolvers
omega1solver = SLEPcEigenSolver(K1,M1)
omega1solver.parameters['spectrum'] = "smallest magnitude"
omega1solver.parameters['tolerance'] = 1e-6
omega1solver.parameters['problem_type'] = "pos_gen_non_hermitian"

omega2solver = SLEPcEigenSolver(K2,M2)
omega2solver.parameters['spectrum'] = "smallest magnitude"
omega2solver.parameters['tolerance'] = 1e-6
omega2solver.parameters['problem_type'] = "pos_gen_non_hermitian"

# Solve the initial problem
omega1solver.solve(1)
omega2solver.solve(1)

# Eigenvalue Extraction
r1, c1, rx1, cx1 = omega1solver.get_eigenpair(0)
r2, c2, rx2, cx2 = omega2solver.get_eigenpair(0)

uu = Function(V)
uu.vector()[:] = rx
uu2 = Function(V2)
uu2.vector()[:] = rx2

## Schwarz Alternating Algorithm
# First Step: Do F1 = du2/dn2 G1 = u2
W2 = VectorFunctionSpace(omega2, 'P', V2.ufl_element().degree())
print(uu2(0.8,0.8))

u1 = TrialFunction(V1)
v1 = TestFunction(V1)

dx = Measure('dx', domain = omega1)
ds = Measure('ds', subdomain_data = omega1_markers)

class OmegaF1(UserExpression):
def eval(self, value, x):
value[0] = f[0]
value[1] = f[1]
class OmegaG1(UserExpression):
def eval(self, value, x):
g = uu2(x[0], x[1])
value[0] = g[0]

f1 = OmegaF1()
g1 = OmegaG1()
n = FacetNormal(omega1)

am1 = u1*v1*dx

K1 = PETScMatrix()
bcs1 = DirichletBC(V1, 0.0, bc_D1)
bcs1.apply(K1)
assemble(ak1, tensor=K1)
M1 = PETScMatrix()
assemble(am1, tensor=M1)

omega1solver = SLEPcEigenSolver(K1,M1)
omega1solver.parameters['spectrum'] = "smallest magnitude"
omega1solver.parameters['tolerance'] = 1e-6
omega1solver.parameters['problem_type'] = "pos_gen_non_hermitian"
omega1solver.solve(1)
r1, c1, rx1, cx1 = omega1solver.get_eigenpair(0)
print("First Iteration Omega 1 Frequency: ", sqrt(r1))


With this code, I run into the following problems:

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
Traceback (most recent call last):
File "SCHWARZ_OverlappingMembraneEigenValueProblem.py", line 196, in <module>
File "/home/careca/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/exproperators.py", line 449, in _getitem
all_indices, slice_indices, repeated_indices = create_slice_indices(component, shape, self.ufl_free_indices)
File "/home/careca/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/ufl/index_combination_utils.py", line 163, in create_slice_indices
if int(ind) >= shape[len(all_indices)]:
IndexError: tuple index out of range


I assume it’s something with the W2 = VectorFunctionSpace function, but I can’t find the solution. Can anyone help me?

I would suggest using the MixedDimensional Finite element method implemented by @cdaversin, whole paper can be found at: https://arxiv.org/pdf/1911.01166.pdf

Using a MeshView for each domain should make it possible to transfer data between the two domains.

You could Also consider interpolating between the non-matching meshes, see: Data projection from one mesh to another mesh

I’m going for the LagrangeInterpolator with this part of the code changed:

## Schwarz Alternating Algorithm
# First Step: Do F1 = du2/dn2 G1 = u2
g1 = Function(V1)
LagrangeInterpolator.interpolate(uu2,g1)
f1 = Function(V1)


I then received the following error:

Building point search tree to accelerate distance queries.
Computed bounding box tree with 2303 nodes for 1152 points.
Traceback (most recent call last):
File "SCHWARZ_OverlappingMembraneEigenValueProblem.py", line 174, in <module>
AttributeError: 'Grad' object has no attribute '_cpp_object'


You should interpolate the function uu2 to the new mesh, and the take the gradient of the projected function, rather than trying to interpolate the gradient on the original mesh to the new mesh.
If you want to take the gradient of the function on the original mesh, i suggest projecting that to a suitable space (usually DG, degree-1), then interpolating the resulting function.

## Schwarz Alternating Algorithm
# First Step: Do F1 = du2/dn2 G1 = u2
V1 = FunctionSpace(omega1,"Lagrange", 1)
V2 = FunctionSpace(omega2, "Lagrange",1)

g1 = Function(V1)
LagrangeInterpolator.interpolate(uu2,g1)
f1 = Function(V1)

ff2 = project(uu2.dx(0) +uu2.dx(1),V2)
LagrangeInterpolator.interpolate(ff2,f1)


Now works! Thanks!

Could you put the correct code all over here please, or send me in private!

I tried to insert those last lines but I couldn’t make it work.