Hi,
I am working on a two-phase stationary Stokes problem. I am trying to modify the geometry in stokesxfem.py in the py_tutorial folder to a unit sphere in the 333 cube scenario. I modified the setting for geometry
cube = CSGeometry()
cube.Add (OrthoBrick(Pnt(-1.5,-1.5,-1.5), Pnt(1.5,1.5,1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=0.5, quad_dominated=False))
names = [ "bottom", "top", "left", "right", "front", "back" ]
for i, name in enumerate(names):
mesh.ngmesh.SetBCName(i, name)
phi = Norm(CoefficientFunction((x,y,z)))**2-1
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2, threshold=10.5, discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
ci = CutInfo(mesh, lsetp1)
d = mesh.dim
and added a part to perform [tex]p^+{new}=p^+{old}-\frac{\int_{\Omega^+}p_{old}\ d\mathscr{H}^{3}}{|\Omega^+|} [/tex], where [tex]\Omega^+[/tex] is the cube minus (unit ball).
const_eliminate=GridFunction(WhG)
p_eliminator=const_eliminate.components[1]
ones_pressure = [1, 1]
area=[Integrate(lset_doms[i], gfp.components[i], mesh=mesh)/Integrate(lset_doms[i],1, mesh=mesh) for i in [0, 1]]
p_eliminator.components[0].Set(area[0])
p_eliminator.components[1].Set(area[1])
L0_gfp=[gfp.components[i]-p_eliminator.components[i] for i in [0, 1]]
When I print out the result, L2 error of velocity and L2 error of pressure within the unit sphere seems fine. But L2 error of pressure outside unit sphere seems blows up. Here is the result I am getting,
For maxh=0.5
L2 Error of velocity: 0.013408092886879695
L2 Error of pressure within unit sphere: 0.018069442530407473
L2 Error of pressure outside unit sphere: 9289563.362909812
For maxh=0.25
L2 Error of velocity: 0.0014547447652420093
L2 Error of pressure within unit sphere: 0.003051588344220756
L2 Error of pressure outside unit sphere: 148769.50294472638
I plot the picture of pressure I got, it seems p^+ blows up near the corner of the cube.
[attachment=undefined]Screen Shot 2020-11-30 at 8.53.29 PM.png[/attachment]
Am I missing some settings for geometry or FEM space? Thank you very much.
Here I attached my code:
[code]from netgen.geom2d import SplineGeometry
from netgen.meshing import MeshingParameters
from ngsolve import *
from xfem import *
visualization stuff
from ngsolve.internal import *
basic xfem functionality
basic geometry features (for the background mesh)
from netgen.csg import CSGeometry, OrthoBrick, Pnt
pi
from math import pi
For LevelSetAdaptationMachinery
from xfem.lsetcurv import *
3D: unit sphere configuration
cube = CSGeometry()
cube.Add (OrthoBrick(Pnt(-1.5,-1.5,-1.5), Pnt(1.5,1.5,1.5)))
mesh = Mesh(cube.GenerateMesh(maxh=0.25, quad_dominated=False))
names = [ “bottom”, “top”, “left”, “right”, “front”, “back” ]
for i, name in enumerate(names):
mesh.ngmesh.SetBCName(i, name)
phi = Norm(CoefficientFunction((x,y,z)))**2-1
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=2, threshold=10.5, discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(phi)
lsetp1 = lsetmeshadap.lset_p1
ci = CutInfo(mesh, lsetp1)
d = mesh.dim
mu=[1.0,1.0]
set analytic problem
aneg = 1 / (2mu[0])(x2+y2)
apos = 1/ (2mu[1])(x2+y2)
gammaf=0
vel_exact = [a*CoefficientFunction((-y, x ,0)) for a in [aneg, apos]]
pres_exact =[x,x]
define RHS
coef_g = [CoefficientFunction((1+4y, -4x,0))for i in [0, 1]]
stabilization parameter for ghost-penalty
gamma_stab_vel = 0.05 # if set to zero: no GP-stabilization for velocity
gamma_stab_pres = 0.05
order=2
stabilization parameter for Nitsche
lambda_nitsche = 0.5 * (mu[0] + mu[1]) * 20 * order * order
Discretization
Vhbase = VectorH1(mesh, order=order, dirichlet=“bottom|top|left|right|front|back”)
Qhbase = H1(mesh, order=order - 1)
Vhneg = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASNEG)))
Vhpos = Compress(Vhbase, GetDofsOfElements(Vhbase, ci.GetElementsOfType(HASPOS)))
Qhneg = Compress(Qhbase, GetDofsOfElements(Qhbase, ci.GetElementsOfType(HASNEG)))
Qhpos = Compress(Qhbase, GetDofsOfElements(Qhbase, ci.GetElementsOfType(HASPOS)))
WhG = FESpace([Vhneg * Vhpos, Qhneg * Qhpos, NumberSpace(mesh)], flags={“dgjumps”: True})
gfup = GridFunction(WhG)
gfu, gfp, gfn = gfup.components
gfu_neg, gfu_pos = gfu.components
gfp_neg, gfp_pos = gfp.components
h = specialcf.mesh_size
lset_neg = {“levelset”: lsetp1, “domain_type”: NEG}
lset_pos = {“levelset”: lsetp1, “domain_type”: POS}
lset_doms = [lset_neg, lset_pos]
lset_if = {“levelset”: lsetp1, “domain_type”: IF}
element, facet and dof marking w.r.t. boundary approximation with lsetp1:
ba_facets = [GetFacetsWithNeighborTypes(mesh, a=ci.GetElementsOfType(HASNEG), b=ci.GetElementsOfType(IF)),
GetFacetsWithNeighborTypes(mesh, a=ci.GetElementsOfType(HASPOS), b=ci.GetElementsOfType(IF))]
n_lset = 1.0 / Norm(grad(lsetp1)) * grad(lsetp1)
kappa = [CutRatioGF(ci), 1.0 - CutRatioGF(ci)]
a = BilinearForm(WhG, symmetric=False)
f = LinearForm(WhG)
u, p, n = WhG.TrialFunction()
v, q, m = WhG.TestFunction()
some helper expressions:
def eps(u):
return 0.5 * (Grad(u) + Grad(u).trans)
def sigma(i, u, p):
return - 2 * mu[i] * eps(u[i]) + p[i] * Id(mesh.dim)
def average_flux(u, p):
return sum([kappa[i] * sigma(i, u, p) * n_lset for i in range(2)])
def average_inv(u):
return sum([kappa[1 - i] * u[i] for i in range(2)])
def jump(u):
return u[0] - u[1]
Stokes variational formulation:
for i in [0, 1]: # loop over domains
a += SymbolicBFI(lset_doms[i], form=2 * mu[i] * InnerProduct(eps(u[i]), eps(v[i])))
a += SymbolicBFI(lset_doms[i], form=- div(u[i]) * q[i] - div(v[i]) * p[i])
f += SymbolicLFI(lset_doms[i], form=coef_g[i] * v[i])
constraining the pressure integral
a += SymbolicBFI(lset_neg, form=n * q[0] + m * p[0])
Nitsche parts:
a += SymbolicBFI(lset_if, form=average_flux(u, p) * jump(v))
a += SymbolicBFI(lset_if, form=average_flux(v, p) * jump(u))
a += SymbolicBFI(lset_if, form=lambda_nitsche / h * jump(u) * jump(v))
surface tension:
f += SymbolicLFI(lset_if, form=-gammaf * average_inv(v) * n_lset)
ghost penalty terms:
for i in range(2):
if gamma_stab_vel > 0:
a += SymbolicFacetPatchBFI(form=gamma_stab_vel / h ** 2 * (u[i] - u[i].Other()) * (v[i] - v[i].Other()),
skeleton=False, definedonelements=ba_facets[i])
if gamma_stab_pres > 0:
a += SymbolicFacetPatchBFI(form=- gamma_stab_pres * (p[i] - p[i].Other()) * (q[i] - q[i].Other()),
skeleton=False, definedonelements=ba_facets[i])
apply mesh adaptation
mesh.SetDeformation(deformation)
a.Assemble()
f.Assemble()
gfu.components[0].Set(vel_exact[0])
gfu.components[1].Set(vel_exact[1],BND)
Solve linear system
f.vec.data -= a.mat * gfup.vec
gfup.vec.data += a.mat.Inverse(WhG.FreeDofs()) * f.vec
enforce p\in L_0^2, by setting L0_gfp^+=gfp^±\integrate(gfp^+)/(\integrate(\Omega^+))
const_eliminate=GridFunction(WhG)
p_eliminator=const_eliminate.components[1]
ones_pressure = [1, 1]
area=[Integrate(lset_doms[i], gfp.components[i], mesh=mesh)/Integrate(lset_doms[i],1, mesh=mesh) for i in [0, 1]]
p_eliminator.components[0].Set(area[0])
p_eliminator.components[1].Set(area[1])
L0_gfp=[gfp.components[i]-p_eliminator.components[i] for i in [0, 1]]
#measure the error
vl2error = sqrt(sum([Integrate(lset_doms[i], Norm(gfu.components[i] - vel_exact[i]) ** 2, mesh=mesh) for i in [0, 1]]))
pl2error_inner =sqrt(sum([Integrate(lset_doms[i], (L0_gfp[i] - pres_exact[i]) ** 2, mesh=mesh) for i in [0]]))
pl2error_outer =sqrt(sum([Integrate(lset_doms[i], (L0_gfp[i] - pres_exact[i]) ** 2, mesh=mesh) for i in [1]]))
print(“L2 Error of velocity: {0}”.format(vl2error))
print(“L2 Error of pressure within unit sphere: {0}”.format(pl2error_inner))
print(“L2 Error of pressure outside unit sphere: {0}”.format(pl2error_outer))
[/code]