Integrating over an interior surface

variational formulation
1

I would like to integrate a function over a surface lying within my domain. Consider the following code:

from dolfin import *
from mshr import *

omega1 = Rectangle(Point(0, 0), Point(0.5, 1))
omega2 = Rectangle(Point(0.5, 0), Point(1, 1))

domain = omega1 + omega2
domain.set_subdomain(1, omega2)

mesh = generate_mesh(domain, 32)

interior_surface = CompiledSubDomain('near(x[0], 0.5, DOLFIN_EPS)')

interior_surface_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
interior_surface_marker.set_all(0)
interior_surface.mark(interior_surface_marker,1)

g = Constant(1)
g = project(g, FunctionSpace(mesh, 'Lagrange', 1))

print(assemble(g*ds(subdomain_id = 1, subdomain_data = interior_surface_marker)))

As you can see, the mesh produced respects the interior surface \{x=0.5\}. The integral should give 1; the result, however, is 0.

On the other hand, when I replace the interior surface with a boundary surface, i.e.,

interior_surface = CompiledSubDomain('on_boundary && near(x[0], 1, DOLFIN_EPS)')

the correct result is computed.

What is the problem with my interior surface?

1 Like
2

From UFL manual ds represents exterior facet integral (that is why your code works for boundary surface). For interior surface dS is needed

assemble(avg(g)*dS(subdomain_id=1, subdomain_data=interior_surface_marker))

2 Likes
3

Thank you, Miroslav. May I expand on my previous question?

In fact, my integrand g has the form

inner(kappa*grad(T),FacetNormal(mesh))

The problem is that kappa and grad(T) are discontinuous across the interior surface and the FacetNormal is ambiguous, too.

I guess the most natural way to resolve this would be to restrict T to the subdomain omega2 and to compute the integral over the corresponding part of the boundary. How can such a restriction be done?

Explanation: I am trying to compute the energy flow across a contact surface. Mathematically and physically, the integral must be independent of the side (material 1 or material 2) that we look at. However, there will be a numerical error. Therefore, I want to compute the integrals over each side separately, and therefore I cannot make use of avg or jump.

4

Hi, consider using restriction using + and -. Here foo('+') restricts foo to the positive cell and FacetNormal('+') is the outer normal of a positive cell. Following the discussion here the +, - side can be changed by providing a form which has a dx measure term which has specified subdomain_data by some cell function g. Then a positive cell for a facet is the one for which value of g is higher. This allows for the following

from dolfin import *
from mshr import *

omega1 = Rectangle(Point(0, 0), Point(0.5, 2))
omega2 = Rectangle(Point(0.5, 0), Point(1, 2))

domain = omega1 + omega2
domain.set_subdomain(1, omega2)

mesh = generate_mesh(domain, 32)

interior_surface = CompiledSubDomain('near(x[0], 0.5, DOLFIN_EPS)')

interior_surface_marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
interior_surface_marker.set_all(0)
interior_surface.mark(interior_surface_marker,1)

V = FunctionSpace(mesh, 'CG', 1)
# Make up something with discontinuous gradient
f = interpolate(Expression('std::max(0., x[0]-0.5)', degree=1), V)

dS = Measure('dS', domain=mesh, subdomain_data=interior_surface_marker, subdomain_id=1)

n = FacetNormal(mesh)

# +, - sides are by default assigned arbitrarily so assembling below
# will not yield meaningful result
L0 = dot(grad(f)('+'), n('+'))*dS
print(assemble(L0))

L1 = dot(grad(f)('-'), n('-'))*dS
print(assemble(L1))

mf = MeshFunction('size_t', mesh, 2, mesh.domains())
# Add dummy terms that will allow for controlling restrictions
# This is (0, 0) * (1, 0) * 2
L2 = dot(grad(f)('-'), n('-'))*dS + Constant(0)*dx(domain=mesh, subdomain_data=mf)
print(assemble(L2))

# This is (1, 0) . (-1, 0) * 2
L3 = dot(grad(f)('+'), n('+'))*dS + Constant(0)*dx(domain=mesh, subdomain_data=mf)
print(assemble(L3))
10 Likes
5

This is working! Thank you!

1 Like
6

Hey everyone,
I’m new in the fenics community and I could use some help right now…
I’m having some troubles in the same subject, I’ve solved the variational problem already, and I wanted to do an integral on my internal boundary. so I have:

## Airfoil subdomain
dInterior = Measure("dS", subdomain_data=subdomains, subdomain_id=5)
n = FacetNormal(mesh)
tangent = as_vector([n[1], -n[0]])
L0 = ((inner(avg(vvec), tangent('+')))*dInterior)
circulation = assemble(L0)

where vvec is my velocity vector. What i’m not understanding is the circulation is always equal to zero, no matter what. When I plot the velocity vector it’s pretty clear it shoudn’t be zero. Am I missing something trivial here? I’ve tried all combinations of (’+’), (’-’) and avg() And it’s still zero.
Meddling around, I decided to try using ds, even though it’s for exterior surfaces and the answer is different than zero, but still wrong (shoudn’t be negative):

ds = Measure("ds", domain=mesh, subdomain_data=subdomains)
GammaP = ds(5)
n = FacetNormal(mesh)
tangent = as_vector([n[1], -n[0]])
circulation = assemble((dot(vvec, tangent))*GammaP)
circulation
>> -18.485904269034087

Here’s an image o my vvec just for demonstration:

image
image1418×427 176 KB

Any help will be greatly appreciated.

7

Could you be more specific of which boundary is supposed to have the marker 5 (used in the dS measure)? If this is the boundary of you object, this is an exterior boundary (interior boundaries are those between two cells in 2D). And you should therefore use ds.

1 Like
8

Hello Dokken,
Ok… seems like I got confused with the definitions… I should be using ds because the boundary with marker 5 is the airfoil profile, but still, the answer doesn’t seems to be correct… I’m still missing something.
Let me share the whole code here. I’m just solving a homogeneus laplacian with Dirichlet boundary conditions to solve for the stream function, and then getting the velocity vector from the stream. The boundary 5 corresponds to the airfoil:

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
# Create mesh
length = 4.0
diameter = 1.0
p0 = Point(np.array([0.0, 0.0]))
p1 = Point(np.array([length, diameter]))
aero_points = np.genfromtxt('naca2412.txt', delimiter = ' ')
aerofoil = [Point(i[0] + 1.5, i[1] + 0.5) for i in aero_points]
aerofoil = Polygon(aerofoil)
domain = Rectangle(p0, p1) - aerofoil
mesh = generate_mesh(domain, 60)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('230*x[1]', degree=1)

# Making a mark dictionary
mark = {"generic": 0,"lower_wall": 1,"upper_wall": 2,"left": 3,"right": 4, "airfoil": 5 }
subdomains = MeshFunction("size_t", mesh, 1)
subdomains.set_all(mark["generic"])
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], length)
class UpperWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], diameter)
class LowerWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0)
class Airfoil(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and not (near(x[1], 0) or near(x[1], diameter) or near(x[0], length) or near(x[0], 0))
# Marking the subdomains
left = Left()
left.mark(subdomains, mark["left"])
right = Right()
right.mark(subdomains, mark["right"])
upper_wall = UpperWall()
upper_wall.mark(subdomains, mark["upper_wall"])
lower_wall = LowerWall()
lower_wall.mark(subdomains, mark["lower_wall"])
airfoil = Airfoil()
airfoil.mark(subdomains, mark["airfoil"])

bc_lower_wall = DirichletBC(V, 0, subdomains, mark["lower_wall"])
bc_upper_wall = DirichletBC(V, u_D((0,diameter)), subdomains, mark["upper_wall"])
bc_left = DirichletBC(V, u_D, subdomains, mark["left"])
bc_right = DirichletBC(V, u_D, subdomains, mark["right"])
bc_airfoil = DirichletBC(V, u_D((2,diameter/2)), subdomains, mark["airfoil"])
bcs = [bc_lower_wall, bc_upper_wall, bc_left, bc_right, bc_airfoil]
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bcs)
W = VectorFunctionSpace(mesh, 'P', 1)
vx = project(u.dx(1), V)
vy = project(-u.dx(0), V)
vvec = project(as_vector([vx, vy]),W)

and then I try to calculate the circulation around the airfoil: (This is wrong based on your answer)

dS = Measure("dS", domain=mesh, subdomain_data=subdomains)
GammaPS = dS(5)
n = FacetNormal(mesh)
tangent = as_vector([n[1], -n[0]])
L0 = avg(dot((vvec), (tangent)))*GammaPS
circulation = assemble(L0)
circulation
>> 0.0

or: (This one’s supposed to be right)

ds = Measure("ds", domain=mesh, subdomain_data=subdomains)
GammaP = ds(5)
n = FacetNormal(mesh)
tangent = as_vector([n[1], -n[0]])
L1 = (dot(vvec, tangent))*GammaP
circulation = assemble(L1)
circulation
>> -18.485904269034172

I thought It might be me messing up the orientations and getting a negative value, but when I tried calculating the circulation for a symmetrical object, which is supposed to be zero, I got a value different than zero, which is wrong. I’m just writing it down here because I believe my reply is long enough for now, and adding more code will make it even longer…
Just to make the code 100% reproducible, here is the content from naca2412.txt:

1.0 -1.6616460838224803e-17
0.9993213820969975 0.00014499032984683925
0.9972872631225536 0.000578728511557539
0.9939028455768911 0.001297536081139998
0.9891767944945985 0.00229534808032572
0.9831212283309125 0.0035638093035125227
0.9757517054010261 0.005092404207668343
0.9670872051355998 0.006868615913445671
0.9571501032924649 0.00887810878330874
0.9459661401982883 0.011104928369901545
0.933564381088738 0.013531712116464912
0.9199771676711413 0.016139904079208404
0.9052400601453433 0.018909967133751822
0.8893917690780644 0.0218215866145346
0.8724740767219776 0.024853860094137484
0.8545317475888985 0.027985469002641854
0.8356124283115054 0.031194828968031024
0.8157665370441834 0.03446021707058515
0.7950471428460603 0.03775987558390232
0.7735098356452664 0.04107209315542856
0.7512125874929715 0.04437526569228436
0.7282156058726268 0.04764794039820733
0.7045811798320735 0.05086884739458386
0.6803735196560986 0.054016924101172765
0.6556585907012681 0.05707133800954382
0.6305039418838468 0.0600115136267799
0.6049785291587992 0.06281716918547092
0.5791525341688195 0.06546836821069141
0.5530971780935735 0.06794559022326693
0.5268845306070138 0.0702298237739368
0.5005873137694109 0.07230268369144069
0.47427870065252875 0.07414655294729396
0.4480321085293681 0.07574474795852291
0.42192098655787974 0.07708170454064288
0.3959866523958327 0.07814207063005084
0.37015695727219483 0.07885512396776277
0.3446799907251953 0.07917983840371742
0.3196303337262879 0.07911176559702549
0.295081536819036 0.07864954600904192
0.2711058222179506 0.0777950160429226
0.24777378350855334 0.07655327310524312
0.2251540869530402 0.07493269119192177
0.20331317876031155 0.07294488068029803
0.18231500284650742 0.07060458743317967
0.16222073357227657 0.0679295280389216
0.1430885276880899 0.06494015996054617
0.12497329925453382 0.061659387472332274
0.10792652064937588 0.05811220644271906
0.09199605195815633 0.05432529319262308
0.07722600011133889 0.0503265447331433
0.06365660812779954 0.046144579584377475
0.051324173805717924 0.04180821002266096
0.04026099622352848 0.03734589793152308
0.030495347529559153 0.03278520738857021
0.02205146675845987 0.028152267666707042
0.014949571856738614 0.023471260439048798
0.009205885759561888 0.018763944643771863
0.004832672255117894 0.014049231695066635
0.0018382775066629242 0.00934282254175867
0.00022717346758919707 0.004656916512940597
0.0 0.0
0.0011432917778369696 -0.0045199873743347885
0.0036398271250637857 -0.008796887680483256
0.007478987149744335 -0.012827539314411746
0.012646513506632423 -0.016608550176758007
0.019124601854193072 -0.0201364086500858
0.0268920169463866 -0.02340763583354531
0.035924225973239104 -0.026418972522308737
0.04619354613387065 -0.029167592939172374
0.05766930200591418 -0.03165133605209156
0.07031798808776175 -0.03386894448976651
0.08410343194323706 -0.03582030063062645
0.09898695366689622 -0.03750664940746651
0.11492751789365321 -0.03893079774434977
0.13188187526807194 -0.04009728131251159
0.14980469112536252 -0.04101249043088248
0.1686486600688652 -0.04168474840283598
0.18836460610365507 -0.042124337318878755
0.20890156894721532 -0.0423434682987691
0.2302068780319326 -0.04235619522042734
0.25222621649144655 -0.042178273105243136
0.27490367804250265 -0.04182696441529416
0.29818182010516364 -0.04132079848343261
0.3220017167284117 -0.04067929106222442
0.34630301489985726 -0.03992263247012663
0.37102399762528443 -0.03907135397193605
0.39610165678640813 -0.038145982808275594
0.4216445484018893 -0.03713442553609651
0.4474394282029785 -0.03599793724535724
0.4733853431045275 -0.03475223822003523
0.49941268623058876 -0.033413794802551806
0.5254514256359297 -0.031998530296521256
0.5514312851740798 -0.03052163424938244
0.5772819308714112 -0.0289974106514835
0.6029331616589602 -0.027439167727298798
0.6283151032186739 -0.025859150186468484
0.6533584036736794 -0.02426851302516897
0.6779944298892017 -0.022677334297823405
0.7021554632437267 -0.02109466278488639
0.72577489386692 -0.01952859522458461
0.7487874125070283 -0.01798637680339546
0.7711291993697603 -0.016474517946209404
0.7927381094464127 -0.014998920131944411
0.8135538540056538 -0.0135650034875589
0.8335181780473524 -0.012177829271078687
0.852575033597649 -0.01084221101582571
0.8706707487554164 -0.009562809033714818
0.8877541923789065 -0.00834420411430316
0.9037769342296041 -0.007190947548681063
0.9186934002742827 -0.00610758598798812
0.9324610226957009 -0.005098661047403128
0.9450403839900795 -0.0041686849205010874
0.9563953543501361 -0.0033220945120995644
0.9664932213616019 -0.0025631876728228263
0.9753048108941276 -0.001896045977266546
0.9828045979581557 -0.001324449092619368
0.9889708062392072 -0.0008517861112930204
0.9937854950182466 -0.0004809692585186644
0.99723463224572 -0.00021435513695233267
0.9993081526575762 -5.367815167922411e-05
1.0 1.6616460838224803e-17

Thank you very much for your attention, Dokken,

9

Consider the following problem for a circular object:

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

def compute_circulation(N):
    # Create mesh
    length = 4.0
    diameter = 1.0
    p0 = Point(np.array([0.0, 0.0]))
    p1 = Point(np.array([length, diameter]))
    r = 0.2
    obj = Circle(Point(length/2, diameter/2), r)
    domain = Rectangle(p0, p1) - obj
    mesh = generate_mesh(domain, N)
    V = FunctionSpace(mesh, 'P', 1)
    # Define boundary condition
    u_D = Expression('230*x[1]', degree=1)

    # Making a mark dictionary
    mark = {"generic": 0,"lower_wall": 1,"upper_wall": 2,"left": 3,"right": 4, "airfoil": 5 }
    subdomains = MeshFunction("size_t", mesh, 1)
    subdomains.set_all(mark["generic"])
    class Left(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], 0)
    class Right(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0], length)
    class UpperWall(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[1], diameter)
    class LowerWall(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[1], 0)
    class Airfoil(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (x[0]-length/2)**2+(x[1]-diameter/2)**2<1.1*r**2
    File("mesh.pvd") << mesh
    # Marking the subdomains
    left = Left()
    left.mark(subdomains, mark["left"])
    right = Right()
    right.mark(subdomains, mark["right"])
    upper_wall = UpperWall()
    upper_wall.mark(subdomains, mark["upper_wall"])
    lower_wall = LowerWall()
    lower_wall.mark(subdomains, mark["lower_wall"])
    airfoil = Airfoil()
    airfoil.mark(subdomains, mark["airfoil"])

    bc_lower_wall = DirichletBC(V, 0, subdomains, mark["lower_wall"])
    bc_upper_wall = DirichletBC(V, u_D((0,diameter)), subdomains, mark["upper_wall"])
    bc_left = DirichletBC(V, u_D, subdomains, mark["left"])
    bc_right = DirichletBC(V, u_D, subdomains, mark["right"])
    bc_airfoil = DirichletBC(V, u_D((2,diameter/2)), subdomains, mark["airfoil"])
    bcs = [bc_lower_wall, bc_upper_wall, bc_left, bc_right, bc_airfoil]
    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(0)
    a = dot(grad(u), grad(v))*dx
    L = f*v*dx
    # Compute solution
    u = Function(V)
    solve(a == L, u, bcs)
    W = VectorFunctionSpace(mesh, 'P', 1)
    vx = project(u.dx(1), V)
    vy = project(-u.dx(0), V)
    vvec = project(as_vector([vx, vy]),W)

    ds = Measure("ds", domain=mesh, subdomain_data=subdomains)
    GammaP = ds(5)
    n = FacetNormal(mesh)
    tangent = as_vector([n[1], -n[0]])
    L1 = (dot(vvec, tangent))*GammaP

    # Project tangent to CG 1 space for visualization
    VV = VectorFunctionSpace(mesh, "CG", 1)
    q,r = TrialFunction(VV), TestFunction(VV)
    aV = inner(q, r)*ds(5)
    lV = inner(tangent, r)*ds(5)
    tang_proj = Function(VV)
    AV = assemble(aV, keep_diagonal=True)
    AV.ident_zeros()
    LV = assemble(lV)
    solve(AV, tang_proj.vector(), LV)
    File("tang.pvd") << tang_proj
    File("circ.pvd") << vvec

    circulation = assemble(L1)
    print(N, circulation, assemble(inner(vvec, as_vector((1,0)))*ds(5)))


compute_circulation(50)
compute_circulation(100)
compute_circulation(200)
compute_circulation(400)

yielding:

fenics@3a01410489e7:~/shared$ python3 circulation.py 
Solving linear variational problem.
50 -2.814902850116198 300.78413974646304
100 -0.35969509592035376 314.1475947113785
200 -0.3435394737202848 323.7011340127153
400 0.03041307517778638 328.270430678979

Here I have made the circulation a function of the mesh resolution (for a circular object). As you can observe, on a coarse mesh, the mesh circulation is far from 0 on the coarsest meshes. Also note that last value is simply the integral of the x-component of the velocity field, is of quite a large magnitude (and the flow field around the obstacle does not satisfy a no penetration constraint.

Therefore I would suggest using a finer mesh.

Note that I have also added a projection of the tangent vector to a CG 1 space, such that it can be visualized with either matplotlib or Paraview.

The accuracy of the circulation integral should also be greatly improved by using a second order mesh. There is limited support for higher order geometries in dolfin, while this has been one of the main focuses of dolfinx

2 Likes
10

Dokken,
Thank you so much for answering! Your answer was perfectly clear and made everything much clearer for me now. Thank you again for your time. Fenics seems to be an awesome platform!
Oh… and I apologize for the off-topic question, since it’s in a interior boundary question. Perhaps i should create a new question about circulation and add your answer?

11

I think we should leave the posts as they are, to avoid duplication;)