Commit 7e650900 authored by Tilman Aleman's avatar Tilman Aleman
Browse files

added plotting

parent 5ef190f1
from ngsolve import *
from xfem.lsetcurv import *
import matplotlib.pyplot as plt
def refineMesh(mesh, phi, order, hs=10 ** 6):
with TaskManager():
......@@ -21,11 +21,11 @@ def refineMesh(mesh, phi, order, hs=10 ** 6):
def rhsfromsol(uSol, pSol, phi, Ps, mesh, ba_IF, weingex, interpol=True, Lpmat=None, Lpvec=None):
if interpol:
with TaskManager():
tmpeps = GridFunction(Lpmat)
tmpeps.Set(eps(uSol, Ps, mesh) - (uSol * Normalize(grad(phi, mesh))) * weingex, definedonelements=ba_IF)
tmptotal = GridFunction(Lpvec)
tmptotal.Set(Ps * (-divG(tmpeps, Ps, mesh) + grad(pSol, mesh)), definedonelements=ba_IF)
tmpeps = GridFunction(Lpmat)
tmpeps.Set(eps(uSol, Ps, mesh) - (uSol * Normalize(grad(phi, mesh))) * weingex, definedonelements=ba_IF)
tmptotal = GridFunction(Lpvec)
tmptotal.Set(Ps * (-divG(tmpeps, Ps, mesh) + grad(pSol, mesh)), definedonelements=ba_IF)
return tmptotal
else:
......@@ -86,3 +86,22 @@ def removekf(VhG, f, kvs, vlams, phi, mesh, maxk=3, tol=0, autofilter=True, inpl
sca = l2sca(f.components[0], kvs[k], phi, mesh)
f.components[0].vec.data += -sca * kvs[k].vec
return
def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
hvals = []
for k in range(len(l2errs)):
hvals.append((1 / 2 ** 3) ** k)
plt.semilogy(l2errs, label=geom+" with KF removal")
plt.semilogy(l2errskf, label=geom + " without KF removal")
plt.semilogy(hvals, label="O(h^"+ str(order+1)+")")
plt.ylabel("L2 error")
plt.xlabel("Refinement")
plt.title("Convergence removing KF")
plt.legend()
plt.xticks(range(len(l2errs)))
if saveplot:
plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+".png")
else:
plt.show()
......@@ -142,7 +142,7 @@ def solveStokes(mesh, lset_approx, deformation, step):
kf = GridFunction(VhG)
print(vlam)
removekf(VhG, gfu, evecinter, vlam, phi, mesh, maxk=maxk)
kf2 = GridFunction(VhG)
for k in range(maxk):
......@@ -158,7 +158,7 @@ def solveStokes(mesh, lset_approx, deformation, step):
# gfu2.components[0].vec.data += integral*kf2.vec
kf.Set(uSol, definedonelements=ba_IF)
kf2.vec.data = gfu.components[0].vec.data - kf.vec.data
for k in range(maxk):
for k in range(0):
kf.Set(kvfs[k], definedonelements=ba_IF)
integral = -Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kf2, kvfs[k]), mesh=mesh, order=5) / Integrate(
......@@ -168,6 +168,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))
l2errnokf = sqrt(
......@@ -184,7 +185,7 @@ if __name__ == "__main__":
# Introduces an interpolation error, but significantly improves performance
parser = argparse.ArgumentParser(description="Solve the Surface Stokes problem")
parser.add_argument("--refinements", dest="n_cut_ref", type=int, required=False, default=2)
parser.add_argument("--write-vtk", dest="vtkout", type=bool, required=False, default=False)
parser.add_argument("--write-vtk", dest="vtkout", action='store_true')
parser.add_argument("--order", dest="order", type=int, required=False, default=2)
parser.add_argument("--geometry", dest="geom", type=str, required=False, default="ellipsoid",
help="Geometry to be considered. Only ellipsoid with parameter c allowed right now")
......@@ -192,6 +193,8 @@ if __name__ == "__main__":
help="Parameter for the Level Set x^2+y^2+(z/c)^2-1=0")
parser.add_argument("-o", dest="fname", type=str, required=False, default="stokestracefem3d",
help="Output file name for vtk")
parser.add_argument("--plot", action='store_true')
parser.add_argument("--saveplot", action='store_true')
args = parser.parse_args()
interpol = True
# Mesh diameter
......@@ -284,13 +287,21 @@ if __name__ == "__main__":
lset_approx, deformation = refineMesh(mesh, phi, order)
print("starting solve..")
l2err, l2errnokf, uerr, uh, uSolnokf = solveStokes(mesh, lset_approx, deformation, step=i)
with TaskManager():
l2err, l2errnokf, uerr, uh, uSolnokf = solveStokes(mesh, lset_approx, deformation, step=i)
l2errors.append(l2err)
l2errkf.append(l2errnokf)
print(l2errors)
print("error with no killing field: ")
print(l2errkf)
if args.plot:
if geom=="ellipsoid" and c==1:
plot_errors(l2errors, l2errkf, "circle", args.saveplot, order=order, refs=n_cut_ref, c=c)
else:
plot_errors(l2errors, l2errkf, "ellipsoid", args.saveplot, order=order, refs=n_cut_ref, c=c)
if vtkout: # and len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
# input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
......
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