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 > 1.75) or (dom_ind_average.vec < 1.25):
interface_indicator = False
else:
interface_indicator = 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 > 1.75): # exclude darcy part
stokes_indicator = False
else:
stokes_indicator = True
if (dom_ind_average.vec < 1.25): # exclude darcy part
darcy_indicator = False
else:
darcy_indicator = 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.*pi*pi)/2.
hmax = 1.0/8

u1 = sin(pi*x)*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./pi*u2))
pexS = CoefficientFunction(-1./pi*u2)
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,cst2*u2))
darcySource = CoefficientFunction(cst3*u2)

dcFlag = False


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

order = 1
orderp = order-1
penalty = 10*order*(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))

## DEBUGGING
# 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()