automatic filtering

helperfunctions.py

@@ -4,9 +4,11 @@ import matplotlib.pyplot as plt
 from netgen.csg import CSGeometry, OrthoBrick, Pnt
 import math
 def refineMesh(mesh, phi, order, hs=10 ** 6):
+    if order>2:
+        hs=10**8
     with TaskManager():
         mesh.SetDeformation(None)
-        lsetp1 = GridFunction(H1(mesh, order=2))
+        lsetp1 = GridFunction(H1(mesh, order=order))
         InterpolateToP1(phi, lsetp1)
         RefineAtLevelSet(lsetp1)
 
@@ -106,84 +108,55 @@ 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
+def computeEV(phi, mesh, lset_approx, deformation, order=2):
+
+    vlams = []
+
+    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)#.Compile()
+
+    def E(u):
+        return ((Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing)#.Compile()
+
+    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()
+
+    return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)

killing.py

@@ -9,11 +9,11 @@ from helperfunctions import *
 
 # -------------------------------- PARAMETERS ---------------------------------
 
-def getKilling(mesh, lset_approx, deformation, step, alpha=3):
+def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsym = 0, maxsymex = 0):
     uSol = (Ps * CF((-z ** 2, x, y)))
     uSol2 = (Ps * CF((-z ** 2, x, y)))
 
-    for k in range(1):
+    for k in range(maxsymex):
         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},
@@ -35,10 +35,12 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
     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])
+#    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)
+    alpha=order+1
     epsilon = h ** alpha
+    epsilon=0
     print("Epsilon:")
     print(epsilon)
     #        epsilon=1
@@ -87,23 +89,7 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
     tmptotal = GridFunction(lpvec)
     tmptotal.Set(Ps * (-divG(tmpeps, Ps, mesh)), definedonelements=ba_IF)
     rhs1 = tmptotal
-    u2, v2 = VhG.TnT()
-    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)
+
     evecinter = []
 
     for k in range(3):
@@ -143,13 +129,13 @@ def getKilling(mesh, lset_approx, deformation, step, alpha=3):
     #                                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
 
-    for k in range(3):
+    for k in range(maxsym):
         
-        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
+        
+        print("removed ", k)
+
+        sca = l2sca(gfu, evecinter[k], phi, mesh)
+        gfu.vec.data += -sca * evecinter[k].vec
     l2err = sqrt(
         Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(gfu - uSol) ** 2, mesh=mesh))
     l2errnokf = sqrt(
@@ -175,7 +161,7 @@ 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)
+    parser.add_argument("--radius", dest="radius", type=float, required=False, default=1)
 
     args = parser.parse_args()
     interpol = True
@@ -199,14 +185,14 @@ if __name__ == "__main__":
     elif c>1 and geom=="ellipsoid":
         maxk = 1
     elif geom=="arrow":
-        maxk=1
+        maxk=0
     # Geometry
     # SetHeapSize(10**9)
     # SetNumThreads(48)
 
     cube = CSGeometry()
     if geom == "ellipsoid":
-        phi = Norm(CoefficientFunction((x, y, z / c))) - args.radius
+        phi = (Norm(CoefficientFunction((x, y, z / c))) - args.radius).Compile(realcompile=True)
 
         cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
         mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
@@ -271,20 +257,52 @@ if __name__ == "__main__":
     weingex = grad(Normalize(grad(phi, mesh)), mesh).Compile(realcompile=True)
     l2errors = []
     l2errkf = []
+    vlams = []
+    symmetries = 0
+    for k in range(3):
+        vlams.append([])
     for i in range(n_cut_ref + 1):
+        
         if i != 0:
             lset_approx, deformation = refineMesh(mesh, phi, order)
-
-        print("starting solve..")
         with TaskManager():
-            l2err, l2errnokf, uerr, uh, uSolnokf = getKilling(mesh, lset_approx, deformation, step=i, alpha=args.alpha)
-        l2errors.append(l2err)
-        l2errkf.append(l2errnokf)
+            vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
+        for k in range(3):
+            vlams[k].append(vlam[k])
+        if i < 2 and i != 0:
+            symmetries = 0
+            for k in range(3):
+#                if math.log(abs(1 - vlams[k][1]), 10) < math.log(abs(vlams[k][0] - vlams[k][1]), 10):
+                if math.log(abs(1 - vlams[k][1]), 2) <-(order+1)*(i+2):
+                    print("corresponding log: ")
+                    print(math.log(abs(1 - vlams[k][1]), 2))
+                    print("condition: ", -(order+1)*(i+1))
+                    symmetries += 1
+        elif i >= 2:
+            symmetries = 0
+            for j in range(3):
+                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]))
+#                if math.log(abs(1 - vlams[j][-1]), 10) < math.log(abs(temp - vlams[j][-1]), 10):
+                if math.log(abs(1 - temp), 2) < -(order + 1) * (i+2):
+                    print("corresponding log: ")
+                    print(math.log(abs(1 - temp), 2))
+                    print("condition: ", -(order+1)*(i+1))
+                    symmetries += 1
+        if i!=0:
+            print("starting solve..")
+            print("symmetries: ", symmetries)
+            with TaskManager():
+                l2err, l2errnokf, uerr, uh, uSolnokf = getKilling(mesh, lset_approx, deformation, vlam, ev, step=i, alpha=args.alpha, maxsym = symmetries, maxsymex = maxk)
+            l2errors.append(l2err)
+            l2errkf.append(l2errnokf)
 
     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: