Commit 9cae2963 authored by Tilman Aleman's avatar Tilman Aleman
Browse files

find dim(K)

parent 03778be8
from ngsolve import *
from xfem.lsetcurv import *
import matplotlib.pyplot as plt
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import math
def refineMesh(mesh, phi, order, hs=10 ** 6):
with TaskManager():
mesh.SetDeformation(None)
......@@ -105,3 +106,84 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+".png")
else:
plt.show()
def computelimit(phi, order=2, iters=5, radius=1):
maxh = 0.6
cube = CSGeometry()
cubesize = 2*radius
with TaskManager():
cube.Add(OrthoBrick(Pnt(-cubesize, -cubesize, -cubesize), Pnt(cubesize, cubesize, cubesize)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True)
deformation = lsetmeshadap.CalcDeformation(phi)
lset_approx = lsetmeshadap.lset_p1
mesh.SetDeformation(deformation)
vlams = []
for k in range(3):
vlams.append([])
for k in range(max(iters,3)):
if k!=0:
lset_approx, deformation = refineMesh(mesh, phi, order)
Vh = VectorH1(mesh, order=order, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
n = Normalize(grad(lset_approx))
ntilde = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
# n=ntilde
h = specialcf.mesh_size
# epsilon=1
eta = 1 / (h) ** 2
rhou = 1.0 / h
rhob = h
# Measure on surface
ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
# Measure on the bulk around the surface
dX = dx(definedonelements=ba_IF, deformation=deformation)
u2 = VhG.TrialFunction()
v2 = VhG.TestFunction()
# Weingarten mappings
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
weing = grad(ninter)
Pmat = Pmats(n)
def E(u):
return (Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing
a2 = BilinearForm(VhG, symmetric=True)
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 += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
vlam, evecs = solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
for j in range(3):
vlams[j].append(vlam[j])
result = []
for j in range(3):
print(len(vlams[j]))
temp=0
temp += vlams[j][-1]
temp -= (vlams[j][-1]-vlams[j][-2])**2/((vlams[j][-1]-vlams[j][-2])-(vlams[j][-2]-vlams[j][-3]))
result.append(math.log(abs(temp-1),10)<math.log(abs(temp-vlams[j][-1]),10))
print("convergence for eigval ", j)
print("logdiff: ", math.log(abs(temp-vlams[j][-1]),10))
print("logex, no accel: ", math.log(abs(1-vlams[j][-1]),10))
print("logex: ", math.log(abs(temp-1),10))
print("logdiff, no accel: ",math.log(abs(vlams[j][-2]-vlams[j][-1]),10))
return result
\ No newline at end of file
......@@ -11,15 +11,16 @@ from helperfunctions import *
def getKilling(mesh, lset_approx, deformation, step, alpha=3):
uSol = (Ps * CF((-z ** 2, x, y)))
for k in range(3):
# 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
uSol2 = (Ps * CF((-z ** 2, x, y)))
for k in range(1):
kfexnorm = Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
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]
print(Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
cf=InnerProduct(uSol-uSol2, uSol-uSol2), mesh=mesh, order=5))
Vh = VectorH1(mesh, order=order, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
......@@ -130,7 +131,7 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
# gfu2.components[0].vec.data += integral*kf2.vec
kf.Set(uSol, definedonelements=ba_IF)
kf2.vec.data = gfu.vec.data - kf.vec.data
kf2.Set(kvfs[k], definedonelements=ba_IF)
for k in range(0):
kf.Set(kvfs[k], definedonelements=ba_IF)
integral = -Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
......@@ -144,7 +145,9 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
for k in range(3):
if vlams[k]<1+cuthmax**3:
if True and (vlam[k]-1)+2*(vlam[k]-1)*(cuthmax**3-2*cuthmax**(5/2))<cuthmax**(3):
print("removed ", k)
print("condition: ", (vlam[k]-1)+2*(vlam[k]-1)*(cuthmax**3-2*cuthmax**(5/2)))
sca = l2sca(gfu, evecinter[k], phi, mesh)
gfu.vec.data += -sca * evecinter[k].vec
l2err = sqrt(
......@@ -172,11 +175,12 @@ if __name__ == "__main__":
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)
parser.add_argument("--radius", dest="radius", type=float, required=False, default=0.5)
args = parser.parse_args()
interpol = True
# Mesh diameter
maxh = 0.6
maxh = 0.6*args.radius
# Geometry
geom = args.geom
# Polynomial order of FE space
......@@ -202,7 +206,7 @@ if __name__ == "__main__":
cube = CSGeometry()
if geom == "ellipsoid":
phi = Norm(CoefficientFunction((x, y, z / c))) - 1
phi = Norm(CoefficientFunction((x, y, z / c))) - args.radius
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
......@@ -280,6 +284,7 @@ if __name__ == "__main__":
print(l2errors)
print("error with no killing field: ")
print(l2errkf)
print("Limits: ", computelimit(phi, radius=args.radius))
if args.plot:
if geom=="ellipsoid" and c==1:
......
......@@ -24,7 +24,6 @@ Literature:
import argparse
from netgen.csg import CSGeometry, OrthoBrick, Pnt
from helperfunctions import *
......
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