Commit 03778be8 authored by Tilman Alemán's avatar Tilman Alemán
Browse files

trying vectorlaplace

parent 7e650900
This diff is collapsed.
#!/usr/bin/python
#!/usr/bin/python3
"""
In this example we solve a surface stokes problem with a
consistent penalty method following [1]
......@@ -31,7 +31,7 @@ from helperfunctions import *
# -------------------------------- PARAMETERS ---------------------------------
def solveStokes(mesh, lset_approx, deformation, step):
def solveStokes(mesh, lset_approx, deformation, step, alpha=3):
uSol = (Ps * CF((-z ** 2, x, y)))
Vh = VectorH1(mesh, order=order, dirichlet=[])
Qh = H1(mesh, order=order - 1, dirichlet=[])
......@@ -49,9 +49,8 @@ def solveStokes(mesh, lset_approx, deformation, step):
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# n=ntilde
h = specialcf.mesh_size
# elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
# cuthmax = max([(2*vol)**(1/2) for vol in elvol])
alpha = 3
#elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
#cuthmax = max([(2*vol)**(1/2) for vol in elvol])
# epsilon =(cuthmax)**(alpha)
epsilon = h ** alpha
print("Epsilon:")
......@@ -168,6 +167,7 @@ def solveStokes(mesh, lset_approx, deformation, step):
# gfu.components[0].vec.data +=-Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
# cf=InnerProduct(gfu.components[0], kvfs[k]), mesh=mesh, order=5) / kfexnorm * kf.vec
# gfu2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso")*f2.vec
removekf(VhG, gfu, evecinter, vlam, phi, mesh, maxk=maxk)
l2err = sqrt(
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
......@@ -195,6 +195,8 @@ if __name__ == "__main__":
help="Output file name for vtk")
parser.add_argument("--plot", action='store_true')
parser.add_argument("--saveplot", action='store_true')
parser.add_argument("--alpha", dest="alpha", type=float, required=False, default=3)
args = parser.parse_args()
interpol = True
# Mesh diameter
......@@ -212,10 +214,12 @@ if __name__ == "__main__":
# Ellipsoid parameter
c = args.c
if c == 1:
if c == 1 and geom=="ellipsoid":
maxk = 3
else:
elif c>1 and geom=="ellipsoid":
maxk = 1
elif geom=="arrow":
maxk=1
# Geometry
# SetHeapSize(10**9)
# SetNumThreads(48)
......@@ -224,6 +228,11 @@ if __name__ == "__main__":
if geom == "ellipsoid":
phi = Norm(CoefficientFunction((x, y, z / c))) - 1
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
elif geom=="arrow":
phi =x**2+y**2+(z/1.5)**2+c*(x*y)**2 - 1
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
elif geom == "decocube":
......@@ -256,7 +265,7 @@ if __name__ == "__main__":
# define solution and right-hand side
if geom == "ellipsoid":
if geom == "ellipsoid" or geom=="arrow":
"""
functions = {
"extvsol1": ((-y - z) * x + y * y + z * z),
......@@ -288,7 +297,7 @@ if __name__ == "__main__":
print("starting solve..")
with TaskManager():
l2err, l2errnokf, uerr, uh, uSolnokf = solveStokes(mesh, lset_approx, deformation, step=i)
l2err, l2errnokf, uerr, uh, uSolnokf = solveStokes(mesh, lset_approx, deformation, step=i, alpha=args.alpha)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment