Direction of boundary normals in a 2D mixed function space problem

1

Hello,

I am trying to solve L(v) term for a mixed function space problem which is supposed be:

L(v) = \int_{\Gamma} n \cdot v \nabla e^{ikd_s \cdot x} \:\mathrm{d}\Gamma

Here is how I separate it for an arbitrarly shaped boundary:

L(v) = \int_{\Gamma} \mathbf{n} \cdot (v_r + i v_i) \mathbf{\nabla} \{ \cos(k \:\mathbf{d_s} \cdot \mathbf{x}) + i \sin(k \:\mathbf{d_s} \cdot \mathbf{x}) \} \:\mathrm{d}\Gamma \quad \text{where} \quad \mathbf{d_s} = (1,0)

From here I somehow end up with some set of equations which actually is possibly wrong but I get the right answer somehow when I play around with + and - signs. I would really appreciate any help regarding this derivation. The plus and negative signs I end up with are all over the place. They don’t make sense. Here is the meaningless but the correct solution:

    d_s = Constant((1.0,0.0))
    
    L_r = \
    - dot( d_s*k, n*v_i*cos(k*x[0]))*ds(2)\
    + dot( d_s*k, n*v_r*sin(k*x[0]))*ds(2)\
    + dot( d_s*k, n*v_i*cos(k*x[0]))*ds(3)\
    - dot( d_s*k, n*v_r*sin(k*x[0]))*ds(3)
    
    L_i = \
    + dot( d_s*k, n*v_r*cos(k*x[0]))*ds(2)\
    + dot( d_s*k, n*v_i*sin(k*x[0]))*ds(2)\
    - dot( d_s*k, n*v_r*cos(k*x[0]))*ds(3)\
    - dot( d_s*k, n*v_i*sin(k*x[0]))*ds(3)

Therefore I think my real question is:

How can I know the direction of the facetnormals for an arbitrarly shaped geometry within my domain. This is not a problem with multiple subdomains. You can think of it as simple as a 2D arbitrary hole shape with two marked boundaries within some rectangular or circular single solution domain.

From some other topic I found a code that helps me to plot the boundary normals:

Here is another topic that I created sometime ago which was solved and is giving more detail about my problem and code but it was for a simple geometry. Now if I try to have a more complicated geometry, I can no longer find the correct answer:

The boundary normals for any geometry I plot always point outwards the solution domain. But the correct solution that I know and find suggests otherways. One of the boundaries should be facing another direction (or somehow should bring some -1 multiplier).

I found some other old post which is trying to find an answer to somehow a similar question:
https://fenicsproject.org/qa/13004/dirichlet-condition-with-normal-vector/

but the codes they shared are already out-dated.

2

Please format the mathematics with a double $ encapsulation.

As you have not provided a minimal example recreating the issue, it is hard to Give you any guidance.

here you are using ds which indicates that you are integrating over an exterior boundary (Where the facetnormal always will be pointing out of the domain). If you want to integrate over an interior boundary, you Need to provide a cell marker to indicate the orientation of “+” and “-”

3

Hi Dokken,

I’ve corrected my post, thank you.

Here is a geometry that I’d like to solve:

rectangle = Rectangle(Point(-0.75, -0.75), Point(0.75, 0.75))
medium = Circle(Point(0, 0), 2, 250)
domain = medium - rectangle
mesh = generate_mesh(domain,30)

I thought cell markers are used when there are cells at the both sides of a boundary? I am trying to solve for the boundaries of the rectangle in the middle. That is why I am not using dS. Should I use dS instead of ds even for this geometry?

4

If you integrate over an exterior boundary (an boundary connected to only one cell), a FacetNormal will always point out of the cell.

5

This clears it up for me. Now I am certain that the error is not coming from the normal directions but the coordinates of the quadrature points. How can I print the coordinates of the quadrature points of a marked boundary? I need to see if they are still on the global coordinate system or not. What I am trying to say is:

L = dot( d_s*k, n*v_r*cos(k*x[0]))*ds(2)

For such an equation which calls a marked boundary, ds(2) in this case, x[0] is not coming from the global coordinate system according to the results that I get but from some shifted local coordinate system. I would like to test if it actually is the case. So how can I print the x[0] values for ds(2)?

6

If x is defined using SpatialCoordinate it will be the physical coordinates at the quadrature points.

There is no nice way of printing these values, as they only exist in the global coordinate system inside the generated code for the integration kernel.

1 Like
7

Not being able to see those coordinates is bad news for me.

Yes. I am getting x (and n) from:

    x = SpatialCoordinate(mesh)
    n = FacetNormal(mesh)

But the error analysis I made with an exact solution for a shifted case and the results that I compared with another solution suggests that I need to use some shift within cos(k*x[0]) term according to the distance that I mark my boundaries - regardless of me getting x[0] from SpatialCoordinate or even from an Expression actually. I tried them both.

Do you mean something else with “physical coordinates at the quadrature points”?

8

The SpatialCoordinate and Expression are features that have been tested in DOLFIN for years, and they do what they should.

You can verify this by integrating a known expression, say

x = SpatialCoordinate(mesh)
r = x[0]**2 + x[1]**2
f = sqrt(r)*cos(x[0])
print(assemble(f*ds(marker_value))

over any known boundary (marked with the value marker_value), and check that you get the expected result.

1 Like
9

Here is the test that I am making following your suggestion:

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# generate a mesh
rectangle = Rectangle(Point(-0.5, -0.75), Point(0.75, 0.75))
cylinder  = Circle(Point(-0.5, 0), 0.5, 250)

medium = Circle(Point(0, 0), 2, 250)
domain = medium - (rectangle - cylinder)

mesh = generate_mesh(domain,30)

# mark the boundary
mf = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
ds_sub = Measure('ds', subdomain_data=mf)

boundary_1 = AutoSubDomain(lambda x: 0.5 + 0.01 >= sqrt(pow(-0.5 - x[0],2) + pow(0 - x[1],2)) >= 0.5 - 0.01 and x[1] < 0.75 + 0.01 and x[1] > -0.75 - 0.01 and x[0] > -0.75 - 0.01)
boundary_1.mark(mf, 1)

# check if the boundary is properly marked
file = File("Marked_Boundary.pvd")
file << mf

x = SpatialCoordinate(mesh)
r = x[0]**2 + x[1]**2
f = sqrt(r)*cos(x[0])

# what result do we expect from here?
assemble(f*ds_sub(1))
# result is: 0.559522861906133

What can we say about this?

10

You can compute this line integral by hand, as you know you are integrating over the half circle (x+0.5)^2+y^2=0.5^2

1 Like
11

Hi Dokken,

I’ve done it. Actually, I’ve made the problem even a little easier for hand-calculation.

For a problem like:

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

rectangle = Rectangle(Point(-1, -3.25), Point(3.25, 3.25))
cylinder  = Circle(Point(-1, 0), 3, 300)

medium = Circle(Point(0, 0), 6, 50)
domain = medium - (rectangle - cylinder)

mesh = generate_mesh(domain,50)

# mark the boundary
mf = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
ds_sub = Measure('ds', subdomain_data=mf)

boundary_1 = AutoSubDomain(lambda x:3 + 0.01 >= sqrt(pow(-1 - x[0],2) + pow(0 - x[1],2)) >= 3 - 0.01 and x[1] < 3.25 + 0.01 and x[1] > -3.25 - 0.01 and x[0] > -3.25 - 0.01)
boundary_1.mark(mf, 1)

# check if the boundary is properly marked
file = File("Marked_Boundary.pvd")
file << mf

x = SpatialCoordinate(mesh)
r = x[0]**2 + x[1]**2
f = r
assemble(f*ds_sub(1))

Answer is: 58.241830033731446

Equation for such a circle is (x-a)^2 + (y-b)^2 = r^2 where a = -1 and b = 0. So that we can write:

x = 3\cos t - 1
y = 3\sin t
x' = -3\sin t
y' = 3\cos t

For the line integral formula:

\int_C f(x, y) \:ds = \int_a^b f(x(t),y(t))\sqrt{\left(\dfrac{dx}{dt}\right)^2 + \left(\dfrac{dy}{dt}\right)^2} \:dt

Substituting the parameters in to the formula:

\int_C f(x^2 + y^2) \:ds = \int_{-\pi/2}^{\pi/2} (3 \cos t - 1)^2 + (3 \sin t)^2\sqrt{\left( -3 \sin t \right)^2 + \left( 3 \cos t \right)^2} \:dt.

Which can be solved like:

\int_C f(x^2 + y^2) \:ds = \int_{-\pi/2}^{\pi/2} \big( 9 \cos^2 t + 1 - 6 \cos t + 9 \sin^2 t \big) \sqrt{ 9 \sin^2 t + 9 \cos^2} \:dt

\int_C f(x^2 + y^2) \:ds = \int_{-\pi/2}^{\pi/2} \big( 30 - 18 \cos t \big) \:dt = \big( 30t - 18 \sin t \big) \bigg\rvert_{-\pi/2}^{\pi/2}

Finally:

\int_C f(x^2 + y^2) \:ds = 30(\pi/2) - 18 \sin (\pi/2) - (30(-\pi/2) - 18 \sin (-\pi/2))

which is:

\int_C f(x^2 + y^2) \:ds = 15\pi - 18 + 15\pi - 18 = 30\pi = 58.24

Looks like we get what we wanted. Though, this does not help me understand my main problem. Because I’ve tried another solution with shifting the geometry in FEM:

# generate a mesh
rectangle = Rectangle(Point(0, -3.25), Point(3.25, 3.25))
cylinder  = Circle(Point(0, 0), 3, 300)

medium = Circle(Point(0, 0), 6, 50)
domain = medium - (rectangle - cylinder)

mesh = generate_mesh(domain,50)

# mark the boundary
mf = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
ds_sub = Measure('ds', subdomain_data=mf)

boundary_1 = AutoSubDomain(lambda x:3 + 0.01 >= sqrt(pow(0 - x[0],2) + pow(0 - x[1],2)) >= 3 - 0.01 and x[1] < 3.25 + 0.01 and x[1] > -3.25 - 0.01 and x[0] > -3.25 - 0.01)
boundary_1.mark(mf, 1)

# check if the boundary is properly marked
file = File("Marked_Boundary.pvd")
file << mf

x = SpatialCoordinate(mesh)
r = x[0]**2 + x[1]**2
f = r
assemble(f*ds_sub(1))

Some other answer: 84.81.

As expected since we have moved our geometry. The analytical solution will require that too. I think this does not give me any information as to why my FEM requires a shift. Because in this case, solving by hand, I have to shift the function too anyways. Or am I very confused?

12

Sorry, I think you are right. I got confused and thought it does not show anything. But it actually shows that SpatialCoordinate works the way it is supposed to. So I guess I have to write my problem all the way down.

This is a 2D acoustic scattering problem. I am trying to solve the Helmholtz equation with a sound-hard boundary condition for the scattering obstacle. I am getting the correct results for this problem when the scattering obstacle is a circular cylinder, and its center is exactly on the center of the domain. But if I change the shape of the scattering obstacle or move the obstacle to some other position within the domain, the solution no longer works. I have found a way around the solution to work - shifting it. But I cannot explain how or why that mechanism works.

Here is the full code:

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

xcyl = 1
## You can also try this code for xcyl = 0 
## and produce correct results

rectangle = Rectangle(Point(-3, -3), Point(3, 3))
cylinder  = Circle(Point(xcyl, 0), 1, 50)
domain    = rectangle - cylinder

mesh = generate_mesh(domain, 50)

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

k = 5

# Define boundary subdomains

boundary_markers = MeshFunction("size_t",mesh, mesh.topology().dim()-1)

class BoundaryL(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], -3, 0.01)

class BoundaryR(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],  3, 0.01)

class BoundaryD(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], -3, 0.01)

class BoundaryU(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1],  3, 0.01)

class BoundaryS(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and sqrt(pow(xcyl - x[0],2) + pow(0 - x[1],2)) <= 1 + 0.01

bc_L = BoundaryL()
bc_R = BoundaryR()
bc_D = BoundaryD()
bc_U = BoundaryU()
bc_S = BoundaryS()

boundary_markers.set_all(0)
bc_L.mark(boundary_markers, 1)
bc_R.mark(boundary_markers, 2)
bc_D.mark(boundary_markers, 3)
bc_U.mark(boundary_markers, 4)
bc_S.mark(boundary_markers, 5)

# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

n = FacetNormal(mesh)
x = SpatialCoordinate(mesh)

# define variational problem
(u_r, u_i) = TrialFunction(V)
(v_r, v_i) = TestFunction(V)

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(1)\
+ k*(inner(u_r,v_i) + inner(u_i,v_r))*ds(2)\
+ k*(inner(u_r,v_i) + inner(u_i,v_r))*ds(3)\
+ k*(inner(u_r,v_i) + inner(u_i,v_r))*ds(4)

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(1)\
- k*(inner(u_r,v_r) - inner(u_i,v_i))*ds(2)\
- k*(inner(u_r,v_r) - inner(u_i,v_i))*ds(3)\
- k*(inner(u_r,v_r) - inner(u_i,v_i))*ds(4)

# This is the correct form that was supposed to work.
# But the results it produces are wrong.
L_r = \
- dot(n, v_r*nabla_grad(cos(k*x[0])))*ds(5)\
+ dot(n, v_i*nabla_grad(sin(k*x[0])))*ds(5)

L_i = \
- dot(n, v_i*nabla_grad(cos(k*(x[0]))))*ds(5)\
- dot(n, v_r*nabla_grad(sin(k*(x[0]))))*ds(5)

## If the equations below are used for L(v) instead of 
## the ones above, which shift the x[0] components,
## this code produces the correct answer - error analysis proved.
## But I have no solid explanation for why that works...
###L_r = \
###- dot(n, v_r*nabla_grad(cos(k*(x[0]-xcyl))))*ds(5)\
###+ dot(n, v_i*nabla_grad(sin(k*(x[0]-xcyl))))*ds(5)
###
###L_i = \
###- dot(n, v_i*nabla_grad(cos(k*(x[0]-xcyl))))*ds(5)\
###- dot(n, v_r*nabla_grad(sin(k*(x[0]-xcyl))))*ds(5)

a = a_r + a_i
L = L_r + L_i

# compute solution
u = Function(V)
solve(a == L, u)

My best guess was the coordinate system was changing for the marked boundaries on L(v) term. That was why shifting worked. But then the test you suggested we’ve done above shows that it is not the case. Though I did not try it for “cos(x[0])”, but I don’t think it is necessary (is it?).

But my question still remains. What causes this shift? Is facetnormal not working the way is it supposed to or am I making some other fundamental error here in my FEM?

I can open another topic which explicitely discusses this problem if necessary? I hope I am able to clarify the problem this time.

13

Looks like my understanding of the problem was wrong. The correct way to tackle this problem actually involves shifting the coordinates if the source is defined on the obstacle surface boundary as with the code I’ve written above. So, my solution was actually right and everything works as they are supposed to. Thanks to Dokken pointing out “facetnormal always points out of the cell for an exterior boundary” and showing that SpatialCoordinate working as it should with the line integral test, it has been clarified.