Dot product of boundary normal on a 2D mixed function-space problem

1

Hi all,

I am trying to implement the boundary integral,

\begin{equation}
\int_{\Gamma} n \cdot v\nabla u :\mathrm{d}\Gamma = \int_{\Gamma} n \cdot (v_r \nabla u_r - v_i \nabla u_i) :\mathrm{d}\Gamma + i\int_{\Gamma} n \cdot (v_r \nabla u_i + v_i \nabla u_r) :\mathrm{d}\Gamma
\end{equation}

to my problem. But I have to solve my problem using a mixed function space since it involves complex numbers.

I am doing this to calculate that integral only on the relevant boundary (on the surface of the inner circle). I actually don’t want other boundaries to be involved but I hope it is handling it with the boundary markings.

Here is my attempt:

from fenics import *
from mshr import *
import numpy as np
from dolfin import *

domain = Circle(Point(0, 0), 2, 120) - Circle(Point(0, 0), 1, 60)
mesh   = generate_mesh(domain,20)

Er = FiniteElement('P', triangle, 2)
Ei = FiniteElement('P', triangle, 2)
Ec = Er * Ei
V = FunctionSpace(mesh,Ec)

k = 6

boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
    
class In(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and pow(x[0],2)+pow(x[1],2)<=pow(1,2)+.01

class Out(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and pow(x[0],2)+pow(x[1],2)>=pow(2,2)-.01

inc = In()
out = Out()
        
boundaries.set_all(0)
inc.mark(boundaries, 1)
out.mark(boundaries, 2)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
   
(u_r, u_i) = TrialFunction(V)
(v_r, v_i) = TestFunction(V)
    
n = FacetNormal(mesh)

a_r = \
+ inner(nabla_grad(u_r), nabla_grad(v_r))*dx\
- inner(nabla_grad(u_i), nabla_grad(v_i))*dx\
- pow(k,2)*( inner(u_r,v_r) - inner(u_i,v_i))*dx\
+ k*(inner(u_r,v_i) + inner(u_i,v_r))*ds(2)

a_i = \
+ inner(nabla_grad(u_r), nabla_grad(v_i))*dx\
+ inner(nabla_grad(u_i), nabla_grad(v_r))*dx\
- pow(k,2)*( inner(u_r,v_i) + inner(u_i,v_r))*dx\
- k*(inner(u_r,v_r) - inner(u_i,v_i))*ds(2)

a = a_r + a_i

# This is where I am having trouble. 

L_r = dot(n, v_r*nabla_grad(u_r))*ds(1) - dot(n, v_i*nabla_grad(u_i))*ds(1)
    
L_i = dot(n, v_r*nabla_grad(u_i))*ds(1) + dot(n, v_i*nabla_grad(u_r))*ds(1)
    
L = L_r + L_i

u = Function(V)
solve(a == L, u )

But it fails. The error is:
“Unable to define linear variational problem a(u, v) = L(v) for all v”

Why does that happen?

2

This is because your L is bilinear (contains u_r and u_i, and should therefore be part of a

1 Like
3

Thank you Dokken. You are right, I changed it. But now, I have no terms in L(v) and it still gives me an errror. Here is the modified code:

Er = FiniteElement('P', triangle, 2)
Ei = FiniteElement('P', triangle, 2)
Ec = Er * Ei
V = FunctionSpace(mesh,Ec)

k = 10

boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)

class In(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and pow(x[0],2)+pow(x[1],2)<=pow(1,2)+.01

class Out(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and pow(x[0],2)+pow(x[1],2)>=pow(2,2)-.01

inc = In()
out = Out()
        
boundaries.set_all(0)
inc.mark(boundaries, 1)
out.mark(boundaries, 2)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
   
(u_r, u_i) = TrialFunction(V)
(v_r, v_i) = TestFunction(V)
  
n = FacetNormal(mesh)

a_r = \
+ inner(nabla_grad(u_r), nabla_grad(v_r))*dx\
- inner(nabla_grad(u_i), nabla_grad(v_i))*dx\
- pow(k,2)*( inner(u_r,v_r) - inner(u_i,v_i))*dx\
+ k*(inner(u_r,v_i) + inner(u_i,v_r))*ds(2)\
+ dot(n, v_r*nabla_grad(u_r))*ds(1)\
- dot(n, v_i*nabla_grad(u_i))*ds(1)

a_i = \
+ inner(nabla_grad(u_r), nabla_grad(v_i))*dx\
+ inner(nabla_grad(u_i), nabla_grad(v_r))*dx\
- pow(k,2)*( inner(u_r,v_i) + inner(u_i,v_r))*dx\
- k*(inner(u_r,v_r) - inner(u_i,v_i))*ds(2)\
+ dot(n, v_r*nabla_grad(u_i))*ds(1)\
+ dot(n, v_i*nabla_grad(u_r))*ds(1)

a = a_r + a_i

L_r = Constant(0)
L_i = Constant(0)

L = L_r + L_i

u = Function(V)
solve(a == L, u )

The error is still the same: Unable to define nonlinear variational problem F(u; v) = 0 for all v.

4

You have now set the RHS to be zero (and not a form). You should either

  1. Define u= Function(V), (u_r, u_i) = split(u) and solve solve(a==0, u)
  2. Define L=inner(Constant(0), v_r)*dx(99) and solve(a==L, u) (where u is still the trial function.
1 Like