Commit 5ef190f1 authored by Tilman Aleman's avatar Tilman Aleman
Browse files

reformat code

parent 6198aee0
from ngsolve import *
from xfem import *
from xfem.lsetcurv import *
def refineMesh(mesh, phi, order, hs=10**6):
def refineMesh(mesh, phi, order, hs=10 ** 6):
with TaskManager():
mesh.SetDeformation(None)
lsetp1 = GridFunction(H1(mesh, order=2))
......@@ -18,8 +19,7 @@ def refineMesh(mesh, phi, order, hs=10**6):
return lset_approx, deformation
def rhsfromsol(uSol, pSol, phi,Ps, mesh, ba_IF, weingex, interpol=True, Lpmat=None, Lpvec=None):
def rhsfromsol(uSol, pSol, phi, Ps, mesh, ba_IF, weingex, interpol=True, Lpmat=None, Lpvec=None):
if interpol:
with TaskManager():
tmpeps = GridFunction(Lpmat)
......@@ -31,6 +31,7 @@ def rhsfromsol(uSol, pSol, phi,Ps, mesh, ba_IF, weingex, interpol=True, Lpmat=No
return Ps * (-divG(eps(uSol, Ps).Compile())) + Ps * grad(pSol, mesh).Compile()
def coef_grad(u, mesh):
dirs = {1: [x], 2: [x, y], 3: [x, y, z]}
if u.dim == 1:
......@@ -66,17 +67,22 @@ def divG(u, Ps, mesh):
else:
N = 3
r = Grad(u) * Ps
return CF(tuple([Trace(r[N*i:N*(i+1),0:N]) for i in range(N)]))
def l2sca(f,g, phi, mesh):
return Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(f,g), mesh=mesh)
def l2norm(f,phi,mesh):
return sqrt(l2sca(f,f,phi,mesh))
def removekf(VhG,f, kvs, vlams, phi, mesh, maxk=3, tol=0, autofilter=True, inplace=False):
return CF(tuple([Trace(r[N * i:N * (i + 1), 0:N]) for i in range(N)]))
def l2sca(f, g, phi, mesh):
return Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=InnerProduct(f, g), mesh=mesh)
def l2norm(f, phi, mesh):
return sqrt(l2sca(f, f, phi, mesh))
def removekf(VhG, f, kvs, vlams, phi, mesh, maxk=3, tol=0, autofilter=True, inplace=False):
print(maxk)
for k in range(3):
if (k<maxk and autofilter) or False:
print(l2norm(kvs[k],phi,mesh))
sca = l2sca(f.components[0],kvs[k], phi, mesh)
f.components[0].vec.data += -sca*kvs[k].vec
return
\ No newline at end of file
if (k < maxk and autofilter) or False:
print(l2norm(kvs[k], phi, mesh))
sca = l2sca(f.components[0], kvs[k], phi, mesh)
f.components[0].vec.data += -sca * kvs[k].vec
return
......@@ -21,27 +21,21 @@ Literature:
"""
# ------------------------------ LOAD LIBRARIES -------------------------------
import sys
from netgen.csg import CSGeometry, OrthoBrick, Pnt
from ngsolve import *
import argparse
from ngsolve.internal import *
from xfem import *
from xfem.lsetcurv import *
from math import pi
import scipy.sparse as sp
from scipy.sparse.linalg import eigs
import time
from netgen.csg import CSGeometry, OrthoBrick, Pnt
from helperfunctions import *
# -------------------------------- PARAMETERS ---------------------------------
def solveStokes(mesh, lset_approx, deformation, step):
uSol = (Ps* CF((-z**2,x,y)))
uSol = (Ps * CF((-z ** 2, x, y)))
Vh = VectorH1(mesh, order=order, dirichlet=[])
Qh = H1(mesh, order=order - 1, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
......@@ -53,17 +47,17 @@ def solveStokes(mesh, lset_approx, deformation, step):
n = Normalize(grad(lset_approx))
ntilde = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# n=ntilde
# 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
# epsilon =(cuthmax)**(alpha)
epsilon=h**alpha
# 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
# epsilon =(cuthmax)**(alpha)
epsilon = h ** alpha
print("Epsilon:")
print(epsilon)
# epsilon=1
eta = 1 / (h)**2
# epsilon=1
eta = 1 / (h) ** 2
rhou = 1.0 / h
rhop = h
rhob = h
......@@ -77,59 +71,53 @@ def solveStokes(mesh, lset_approx, deformation, step):
# Weingarten mappings
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = grad(ninter)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
L21 = L2(mesh, order=order+1, dim=9)
L22 = L2(mesh, order=order+1, dim=3)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
L21 = L2(mesh, order=order + 1, dim=9)
L22 = L2(mesh, order=order + 1, dim=3)
lpmat = Restrict(L21, ba_IF)
lpvec = Restrict(L22, ba_IF)
if step>4:
comp=True
if step > 4:
comp = True
else:
comp=False
comp = False
print("came here ")
rhs1 = rhsfromsol(uSol, pSol,phi,Ps, mesh,ba_IF,weingex, interpol, lpmat, lpvec).Compile(realcompile=comp)
rhs1 = rhsfromsol(uSol, pSol, phi, Ps, mesh, ba_IF, weingex, interpol, lpmat, lpvec).Compile(realcompile=comp)
rhs2 = Trace(Ps * grad(uSol, mesh)).Compile(realcompile=comp)
Pmat = Pmats(n)
def E(u):
return (Pmat * Sym(grad(u)) * Pmat)- (u * n) * weing
return (Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing
a = BilinearForm(X, symmetric=True)
a += InnerProduct(E(u), E(v)).Compile(realcompile=comp) * ds
a +=epsilon* (Pmat * u) * (Pmat * v) * ds
a += epsilon * (Pmat * u) * (Pmat * v) * ds
a += (eta * ((u * ntilde) * (v * ntilde))) * ds
a += rhou * ((grad(u) * n) * (grad(v) * n)) * dX
a += InnerProduct(u, Pmat * grad(q)) * ds
a += InnerProduct(v, Pmat * grad(p)) * ds
a += -rhop * InnerProduct(n * grad(p), n * grad(q)) * dX
u2, v2 = VhG.TnT()
a2 = BilinearForm(VhG, symmetric=True)
a2 += InnerProduct(E(u2),E(v2)) * ds
a2 +=(Pmat * u2) * (Pmat * v2) * ds
a2 += InnerProduct(E(u2), E(v2)) * ds
a2 += (Pmat * u2) * (Pmat * v2) * ds
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
pre = Preconditioner(a2, "direct")
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(Pmat*u2,Pmat*v2)*ds
b += rhob*(grad(u2)*n)*(grad(v2)*n)*dX
b += InnerProduct(Pmat * u2, Pmat * v2) * ds
b += rhob * (grad(u2) * n) * (grad(v2) * n) * dX
b += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
vlam, evecs = solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
evecinter = []
......@@ -154,55 +142,56 @@ def solveStokes(mesh, lset_approx, deformation, step):
kf = GridFunction(VhG)
print(vlam)
removekf(VhG,gfu, evecinter,vlam, phi, mesh, maxk=maxk)
removekf(VhG, gfu, evecinter, vlam, phi, mesh, maxk=maxk)
kf2 = GridFunction(VhG)
for k in range(maxk):
# kf2.Set(kvfs[k],definedonelements=ba_IF)
# kfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kf2,kf2), mesh=mesh, order=5)
# integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu2.components[0],kf2), mesh=mesh, order=5)/kfnorm
# kf2.Set(kvfs[k],definedonelements=ba_IF)
# kfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kf2,kf2), mesh=mesh, order=5)
# integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu2.components[0],kf2), mesh=mesh, order=5)/kfnorm
kfexnorm = Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kvfs[k], kvfs[k]), mesh=mesh, order=5)
cf=InnerProduct(kvfs[k], kvfs[k]), mesh=mesh, order=5)
uSol = uSol - Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(uSol, kvfs[k]), mesh=mesh, order=5) / kfexnorm * kvfs[k]
cf=InnerProduct(uSol, kvfs[k]), mesh=mesh, order=5) / kfexnorm * kvfs[k]
# 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
# 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):
kf.Set(kvfs[k],definedonelements=ba_IF)
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(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kvfs[k], kvfs[k]), mesh=mesh, order=5)
kf2.vec.data += integral*kf.vec
# 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
l2err = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
l2errnokf = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kf2,kf2) , mesh=mesh))
cf=InnerProduct(kf2, kvfs[k]), mesh=mesh, order=5) / Integrate(
levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(kvfs[k], kvfs[k]), mesh=mesh, order=5)
kf2.vec.data += integral * kf.vec
# 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
l2err = sqrt(
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
l2errnokf = sqrt(
Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=InnerProduct(kf2, kf2), mesh=mesh))
uerr = gfu.components[0] - uSol
print("l2error : ", l2err)
print("error without killing error: ", l2errnokf)
return l2err, l2errnokf, kf2, gfu, uSol
if __name__=="__main__":
if __name__ == "__main__":
# Interpolation for the rate-of-strain tensor when computing the r.h.s.
# 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("--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("--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")
parser.add_argument("-c", dest="c", type=float, required=False, default=1, 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("--geometry", dest="geom", type=str, required=False, default="ellipsoid",
help="Geometry to be considered. Only ellipsoid with parameter c allowed right now")
parser.add_argument("-c", dest="c", type=float, required=False, default=1,
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")
args = parser.parse_args()
interpol = True
# Mesh diameter
......@@ -210,7 +199,7 @@ if __name__=="__main__":
# Geometry
geom = args.geom
# Polynomial order of FE space
order=args.order#=args["order"]
order = args.order # =args["order"]
vtkout = args.vtkout
n_cut_ref = args.n_cut_ref
if order == 2:
......
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