# How to deal with the dirichlet boundary condition

If the Stokes-Darcy problem is dirichlet boundary condition, that is the speed of Stokes domain and the pressure of Darcy region are known. According to the code that previously dealt with Neumann boundary conditions we made the following modifications. However, that’s not what we want. Anybody know why?

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from fealpy.tools.show import showmultirate, show_error_table
from ngsolve import *
from netgen.geom2d import SplineGeometry
from numpy import pi
import netgen.gui

ngsglobals.msg_level=1

errorType = [‘||us - us_h||_{\Omega,0}’,
||ud - ud_h||_{\Omega, 0}
]

maxit = 1
errorMatrix = np.zeros((2, maxit), dtype=np.float)
NDof = np.zeros(maxit, dtype=np.float)

def topBottom():

geometry = SplineGeometry()

# point coordinates …

pnts = [ (0,0), (1,0), (1,0.5), (1,1), (0,1), (0, 0.5)]
pnums = [geometry.AppendPoint(*p) for p in pnts]

# start-point, end-point, boundary-condition, domain on left side, domain on right side:

lines = [ (pnums[0],pnums[1],“gammaD”,2,0), (pnums[1],pnums[2],“gammaD”,2,0),
(pnums[2],pnums[3],“gammaS”,1,0), (pnums[3],pnums[4],“gammaS”,1,0),
(pnums[4],pnums[5],“gammaS”,1,0),
(pnums[5],pnums[0],“gammaD”,2,0), (pnums[5],pnums[2],“gammaSD”,1,2)]
for p1,p2,bc,left,right in lines:
geometry.Append( [“line”, p1, p2], bc=bc, leftdomain=left, rightdomain=right)
return geometry

def GetInterfaceIndicator(mesh):
dom_indicator = GridFunction(L2(mesh,order=0))
dom_indicator.Set(CoefficientFunction([1,2]))

ffes = FacetFESpace(mesh,order=0)
dom_ind_average = GridFunction(ffes)
mf = BilinearForm(ffes,symmetric=True)
mf += SymbolicBFI(ffes.TrialFunction()ffes.TestFunction(),element_boundary=True)
mf.Assemble()
ff = LinearForm(ffes)
ff += SymbolicLFI(dom_indicator
ffes.TestFunction(),element_boundary=True)
ff.Assemble()
dom_ind_average.vec.data = mf.mat.Inverse() * ff.vec

interface_indicator = BitArray(ffes.ndof)
for i in range(len(dom_ind_average.vec)):
if (dom_ind_average.vec[i] > 1.75) or (dom_ind_average.vec[i] < 1.25):
interface_indicator[i] = False
else:
interface_indicator[i] = True
print(sum(interface_indicator))
return interface_indicator

def GetFaceIndicator(mesh):
dom_indicator = GridFunction(L2(mesh,order=0))
dom_indicator.Set(CoefficientFunction([1,2]))

ffes = FacetFESpace(mesh,order=0)
dom_ind_average = GridFunction(ffes)
mf = BilinearForm(ffes,symmetric=True)
mf += SymbolicBFI(ffes.TrialFunction()ffes.TestFunction(),element_boundary=True)
mf.Assemble()
ff = LinearForm(ffes)
ff += SymbolicLFI(dom_indicator
ffes.TestFunction(),element_boundary=True)
ff.Assemble()
dom_ind_average.vec.data = mf.mat.Inverse() * ff.vec

stokes_indicator = BitArray(ffes.ndof)
darcy_indicator = BitArray(ffes.ndof)
for i in range(len(dom_ind_average.vec)):
if (dom_ind_average.vec[i] > 1.75): # exclude darcy part
stokes_indicator[i] = False
else:
stokes_indicator[i] = True
if (dom_ind_average.vec[i] < 1.25): # exclude darcy part
darcy_indicator[i] = False
else:
darcy_indicator[i] = True
#print(“facets marked as stokes facets:”, stokes_indicator)
#print(“facets marked as darcy facets:”, darcy_indicator)
return stokes_indicator, darcy_indicator

def tang(vec):
return vec - (vec*n)*n

stokesId = [1,0]
darcyId = [0,1]

nu = 1.
kinv = 1.
gamma = (1.+4.pipi)/2.
hmax = 1.0/8

u1 = sin(pix)exp(y/2.)
u2 = cos(pi
x)exp(y/2.)
uexS = CoefficientFunction((-1./2./pi/pi
u1,1./pi
u2))
uexD = CoefficientFunction((-2.u1,1./piu2))
pexS = CoefficientFunction(-1./piu2)
pexD = CoefficientFunction(-2./pi
u2)

cst1 = 1.+nu*(-.5+1./8./pi/pi)
cst2 = -.5/pi+nu*(pi-.25/pi)
cst3 = 0.5/pi-2.pi
stokesSource = CoefficientFunction((cst1
u1,cst2u2))
darcySource = CoefficientFunction(cst3
u2)

dcFlag = False

n = specialcf.normal(2)
h = specialcf.mesh_size

order = 1
orderp = order-1
penalty = 10order(order+1)/h
InitialMeshSize = 1

for i in range(maxit):
hmax = InitialMeshSize/2**i
mesh = Mesh(topBottom().GenerateMesh (maxh=hmax))

``````dom1 = GridFunction(L2(mesh,order=0))
dom2 = GridFunction(L2(mesh,order=0))
dom1.Set(CoefficientFunction(stokesId))
dom2.Set(CoefficientFunction(darcyId))

V1 = HDiv(mesh, order=order, dirichlet="gammaS")
``````

# FIXME define trace space only on stokes domain

``````V2 = FESpace ("vectorfacet", mesh, order=order, dirichlet="gammaS",definedon = [1], flags
={"highest_order_dc": dcFlag})

V3 = L2(mesh, order=orderp)

V = FESpace([V1,V2, V3], flags={"dgjumps": True})

u,uhat, p = V.TrialFunction()
v,vhat, q = V.TestFunction()

crossTerm = -div(u)*q-div(v)*p

du = grad(u)
dv = grad(v)
gradu = CoefficientFunction ( (du,), dims=(2,2) )
gradv = CoefficientFunction ( (dv,), dims=(2,2) )
epsu = 0.5*(gradu+gradu.trans)
epsv = 0.5*(gradv+gradv.trans)

interface_indicator = GetInterfaceIndicator(mesh)

a = BilinearForm(V)
``````

# dg stokes

``````stokesCell = 2*nu*InnerProduct(epsu,epsv)
stokesFacet1 = 2*nu*InnerProduct ( epsu*n,  tang(vhat-v) )
stokesFacet2 = 2*nu*InnerProduct ( epsv*n,  tang(uhat-u) )
stokesFacet3 = penalty*nu*InnerProduct (tang(vhat-v),tang(uhat-u) )
``````

# Hdiv darcy

``````darcyCell = kinv*u*v
``````

# interface (only stokes part is needed)

``````interfaceTerm = gamma*sqrt(kinv)*(tang(uhat)*tang(vhat))
``````

# Stokes is correct

``````a += SymbolicBFI ( stokesCell+dom2*darcyCell + crossTerm, definedon=[1])

a += SymbolicBFI ( darcyCell + crossTerm, definedon = [2])
a += SymbolicBFI ( (stokesFacet1+stokesFacet2+stokesFacet3), element_boundary=True, definedon=[1])

interfacebfi = SymbolicBFI (interfaceTerm, VOL, skeleton=True)
interfacebfi.SetDefinedOnElements(interface_indicator)
a += interfacebfi

a.Assemble()

c = Preconditioner(type="direct", bf=a, flags = { "inverse" : "pardiso" } )

f = LinearForm(V)
f += SymbolicLFI(dom1*stokesSource*v-dom2*darcySource*q)
f += -pexD*(v.Trace()*n)*ds(definedon=mesh.Boundaries("gammaD"))
f += -pexS*(v.Trace()*n)*ds(definedon=mesh.Boundaries("gammaS"))

f.Assemble()

u = GridFunction(V)

u.components[0].Set(dom1*uexS,definedon=mesh.Boundaries("gammaS"))
u.components[1].Set(uexS,definedon=mesh.Boundaries("gammaS"))

bvp =  BVP(bf=a,lf=f,gf=u,pre=c, maxsteps=1000, prec=1e-11)

errS = dom1*(u.components[0]-uexS)*(u.components[0]-uexS)
errD = dom2*(u.components[0]-uexD)*(u.components[0]-uexD)
erru1 = sqrt( Integrate(errS,mesh,order = 2*order) )
erru2 = sqrt( Integrate(errD,mesh,order = 2*order) )
errorMatrix[0, i] = erru1
errorMatrix[1, i] = erru2
``````

show_error_table(NDof, errorType, errorMatrix)

plt.show()