refactoring of code, command line params

helperfunctions.py

+from ngsolve import *
+from xfem import *
+from xfem.lsetcurv import *
+def refineMesh(mesh, phi, order, hs=10**6):
+    with TaskManager():
+        mesh.SetDeformation(None)
+        lsetp1 = GridFunction(H1(mesh, order=2))
+        InterpolateToP1(phi, lsetp1)
+        RefineAtLevelSet(lsetp1)
+
+        mesh.Refine()
+
+        lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
+
+        deformation = lsetmeshadap.CalcDeformation(phi)
+        mesh.SetDeformation(deformation)
+        lset_approx = lsetmeshadap.lset_p1
+    return lset_approx, deformation
+
+
+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)
+        return tmptotal
+    else:
+
+        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:
+        return CF(tuple([u.Diff(r) for r in dirs[mesh.dim]]))
+    else:
+        return CF(tuple([u[i].Diff(r) for i in range(u.dim) for r in dirs[mesh.dim]]), dims=(u.dim, mesh.dim))
+
+
+def grad(u, mesh=None):
+    if type(u) in [ProxyFunction, GridFunction]:
+        return ngsolve.grad(u)
+    else:
+        return coef_grad(u, mesh)
+
+
+def Pmats(n):
+    return Id(3) - OuterProduct(n, n)
+
+
+# Helper Functions for rate-of-strain-tensor
+def eps(u, Ps, mesh):
+    return Ps * Sym(grad(u, mesh)) * Ps
+
+
+def divG(u, Ps, mesh):
+    if u.dim == 3:
+        return Trace(grad(u, mesh) * Ps)
+    if u.dim == 9:
+        if u.dims[0] == 3:
+            N = 3
+            divGi = [divG(u[i, :]) for i in range(N)]
+            return CF(tuple([divGi[i] for i in range(N)]))
+        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):
+
+    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

killing.py

@@ -10,6 +10,7 @@ import scipy.sparse as sp
 from scipy.sparse.linalg import eigs
 import numpy as np
 import time
+from evsolver import *
 # -------------------------------- PARAMETERS ---------------------------------
 # Interpolation for the rate-of-strain tensor when computing the r.h.s.
 # Introduces an interpolation error, but significantly improves performance
@@ -19,30 +20,30 @@ maxh = 0.6
 # Geometry
 geom = "ellipsoid"
 # Polynomial order of FE space
-order = 2
-hs=10**8
+order = 3
+hs = 10 ** 8
 # Problem parameters
 reac_cf = 1
 diff_cf = 1
 # Ellipsoid parameter
-c = 1.5
-if c==1:
+c = 1
+if c == 1:
     maxk = 3
 else:
-    maxk=1
+    maxk = 1
 # Geometry
-#SetHeapSize(10**9)
-#SetNumThreads(48)
+# SetHeapSize(10**9)
+# SetNumThreads(48)
 with TaskManager():
     cube = CSGeometry()
     if geom == "ellipsoid":
-        phi = Norm(CoefficientFunction((x , y, z/c))) - 1
+        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 == "decocube":
         phi = (x ** 2 + y ** 2 - 4) ** 2 + (y ** 2 - 1) ** 2 + (y ** 2 + z ** 2 - 4) ** 2 + (x ** 2 - 1) ** 2 + (
-                    x ** 2 + z ** 2 - 4) ** 2 + (z ** 2 - 1) ** 2 - 13
+                x ** 2 + z ** 2 - 4) ** 2 + (z ** 2 - 1) ** 2 - 13
         cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
         mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
 
@@ -50,13 +51,11 @@ with TaskManager():
     # When interpolating, getting good results requires more refinements with complex geometries, e.g. the decocube
     n_cut_ref = int(sys.argv[1])
 
-
     # Class to compute the mesh transformation needed for higher order accuracy
     #  * order: order of the mesh deformation function
     #  * threshold: barrier for maximum deformation (to ensure shape regularity)
     if geom == "ellipsoid":
-
-        kvfs = [CoefficientFunction((y, -x, 0)),CF((z,0,-x)), CF((0,z,-y))]
+        kvfs = [CoefficientFunction((y, -x, 0)), CF((z, 0, -x)), CF((0, z, -y))]
 
     lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
     deformation = lsetmeshadap.CalcDeformation(phi)
@@ -65,8 +64,7 @@ with TaskManager():
     mesh.SetDeformation(deformation)
 
 
-# Background FESpaces
-
+    # Background FESpaces
 
     # Helper Functions for tangential projection and gradients of the exact solution
     def coef_grad(u):
@@ -91,12 +89,11 @@ with TaskManager():
     # Coefficients / parameters:
 
     # Exact tangential projection
-    #.Compile(realcompile=True)
+    # .Compile(realcompile=True)
     def doKilling(mesh, lset_approx, deformation, step):
         Ps = Pmats(Normalize(grad(phi)))
         Vh = VectorH1(mesh, order=order, dirichlet=[])
 
-
         ci = CutInfo(mesh, lset_approx)
         ba_IF = ci.GetElementsOfType(IF)
         VhG = Restrict(Vh, ba_IF)
@@ -107,13 +104,13 @@ with TaskManager():
 
         h = specialcf.mesh_size
 
-        alpha=3
-#        epsilon =(cuthmax)**(alpha)
-        epsilon=h**alpha
-#        print("Epsilon:")
-#        print(epsilon)
-#        epsilon=1
-        eta = 1000.0 / (h)**2
+        alpha = 3
+        #        epsilon =(cuthmax)**(alpha)
+        epsilon = h ** alpha
+        #        print("Epsilon:")
+        #        print(epsilon)
+        #        epsilon=1
+        eta = 1000.0 / (h) ** 2
         rhou = 1.0 / h
         rhop = h
 
@@ -124,117 +121,126 @@ with TaskManager():
 
         # Weingarten mappings
 
-
         ninter = GridFunction(VhG)
         ninter.Set(n, definedonelements=ba_IF)
         weing = ngsolve.grad(ninter)
 
-#        rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
+        #        rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
 
-        if step>4:
-            comp=True
+        if step > 4:
+            comp = True
         else:
-            comp=False
+            comp = False
 
         Pmat = Pmats(n)
 
         def E(u):
-            return (Pmat * Sym(grad(u)) * Pmat)- (u * n) * weing
+            return (Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing
 
-        
         rhob = h
         u2, v2 = VhG.TnT()
-        a2 = BilinearForm(VhG, symmetric=True)
+        a2 = BilinearForm(VhG)
 
-        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()
-#        a2inv = a2.Inverse(dofs=VhG.FreeDofs(), inverse="pardiso")
-        b = BilinearForm(VhG, symmetric=True)
-        b += InnerProduct(Pmat*u2,Pmat*v2)*ds 
-        b += rhob*(grad(u2)*n)*(grad(v2)*n)*dX
+        #        a2inv = a2.Inverse(dofs=VhG.FreeDofs(), inverse="pardiso")
+        b = BilinearForm(VhG)
+        b += InnerProduct(Pmat * u2, Pmat * v2) * ds
+        b += rhob * (grad(u2) * n) * (grad(v2) * n) * dX
         b += ((u2 * n) * (v2 * n)) * ds
+
+
         b.Assemble()
-        def bh(u,v):
-            ret = InnerProduct(Pmat*u, Pmat*v)
-            ret += rhob*(Grad(u)*n)*(Grad(v)*n)
+
+  #      print(preJpoint.COO())
+        def bh(u, v):
+            ret = InnerProduct(Pmat * u, Pmat * v)
+            ret += rhob * (Grad(u) * n) * (Grad(v) * n)
             ret += ((u * n) * (v * n))
             return ret
-        def l2sca(u,v):
-            return Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(u,v), mesh=mesh)
+
+        def l2sca(u, v):
+            return Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=bh(u, v), mesh=mesh)
+
         def l2norm(u):
-            return l2sca(u,u)
+            return l2sca(u, u)
+
         def orthproj(f, kfs):
             res = GridFunction(VhG)
             res.vec.data = f.vec.data
             for k in range(maxk):
-                sca = l2sca(f,kfs[k])#/l2norm(kfs[k])
-                res.vec.data += -sca*kfs[k].vec
+                sca = l2sca(f, kfs[k]) /l2norm(kfs[k])
+                print(sca)
+                res.vec.data += -sca * kfs[k].vec
             return res
-        evecs = GridFunction(VhG, multidim=3)
 
-        evecs2 = GridFunction(VhG, multidim=3)
-        evecs3 = GridFunction(VhG, multidim=3)
-        vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs2.vecs, 1)
+
+
+        vlam, evecs = solvers.PINVIT(a2.mat,b.mat,pre,6,GramSchmidt=True)
         evecinter = []
+        """
         for k in range(3):
             evecs3.vecs[k].data = evecs2.vecs[k]
+            """
         for k in range(3):
-            evn =   sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(evecs3.MDComponent(k))**2, mesh=mesh) )
-            evecs.vecs[k].data =  evecs3.vecs[k].data
+
+#            evn = sqrt(Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(evecs3.MDComponent(k)) ** 2,
+#                                 mesh=mesh))
+#            evecs.vecs[k].data = evecs3.vecs[k].data
             evecinter.append(GridFunction(VhG))
-            evecinter[k].Set(evecs.MDComponent(k), definedonelements=ba_IF)
+            evecinter[k].vec.data = evecs[k].data
+#            evecinter[k].Set(evecs[k], definedonelements=ba_IF)
             print("evecsgrad: ")
-#            print(l2norm(Grad(evecinter[k])))  
-        print("IP: ", InnerProduct(evecs.vecs[1],evecs2.vecs[0]))
+        #            print(l2norm(Grad(evecinter[k])))
+#        print("IP: ", InnerProduct(evecs.vecs[1], evecs2.vecs[0]))
 
-        
         if False:
-                    
+
             sizeA = a2.mat.height
             print("sizeA: ", sizeA)
-            
-            
-            atemp = a2.mat.CreateMatrix()    
-            atemp.AsVector().data = a2.mat.AsVector() -b.mat.AsVector()
-        #    atemp.AsVector().data = m.mat.AsVector()
+
+            atemp = a2.mat.CreateMatrix()
+            atemp.AsVector().data = a2.mat.AsVector() - b.mat.AsVector()
+            #    atemp.AsVector().data = m.mat.AsVector()
             ainv = atemp.Inverse(VhG.FreeDofs(), inverse="pardiso")
-            
+
             x1 = ninter.vec.CreateVector()
             y1 = ninter.vec.CreateVector()
-            
+
             def A_inv(v):
                 for j in range(sizeA):
-                    y1[j] = v[j]         
-                x1.data = ainv*y1
+                    y1[j] = v[j]
+                x1.data = ainv * y1
                 for j in range(sizeA):
-                    v[j] = x1[j]         
+                    v[j] = x1[j]
                 return v
 
-            Ainv = sp.linalg.LinearOperator( (sizeA,sizeA), A_inv)
-            
-            rowsA,colsA,valsA = a2.mat.COO()
-            Amattemp = sp.csr_matrix((valsA,(rowsA,colsA))) 
-            Amat = Amattemp + Amattemp.transpose() - sp.diags(Amattemp.diagonal(),0)
-            rowsM,colsM,valsM = b.mat.COO()
-            Mmattemp = sp.csr_matrix((valsM,(rowsM,colsM))) 
-            Mmat = Mmattemp + Mmattemp.transpose() - sp.diags(Mmattemp.diagonal(),0)
-            
-            vals, vecs = sp.linalg.eigsh(Amat, k=3, M=Mmat, sigma=1, which='LM', OPinv=Ainv)
+            Ainv = sp.linalg.LinearOperator((sizeA, sizeA), A_inv)
+
+            rowsA, colsA, valsA = a2.mat.COO()
+            Amattemp = sp.csr_matrix((valsA, (rowsA, colsA)))
+            Amat = Amattemp + Amattemp.transpose() - sp.diags(Amattemp.diagonal(), 0)
+            rowsM, colsM, valsM = b.mat.COO()
+            Mmattemp = sp.csr_matrix((valsM, (rowsM, colsM)))
+            Mmat = Mmattemp + Mmattemp.transpose() - sp.diags(Mmattemp.diagonal(), 0)
+
+            vals, vecs = sp.linalg.eigsh(Amat, k=30, M=Mmat, sigma=1, which='LM', OPinv=Ainv)
             ui = []
             for k in range(maxk):
                 ui.append(GridFunction(VhG))
-                ui[k].vec.FV().NumPy()[:] = vecs[:,k]
-                evecs.vecs[k].data = ui[0].vec.data
+                ui[k].vec.FV().NumPy()[:] = vecs[:, k]
+                evecinter[k].vec.data = ui[k].vec.data
+                evecs.vecs[k].data = ui[k].vec.data
                 print("evecsnorms: ")
-#                print(l2norm(evecs.MDComponent(k)))
+        #                print(l2norm(evecs.MDComponent(k)))
 
         egfu = GridFunction(VhG)
-        egfu.vec.data = evecs.vecs[0].data
-            
+        egfu.vec.data = evecs[0].data
+
         def divG(u):
             if u.dim == 3:
                 return Trace(grad(u) * Ps)
@@ -246,8 +252,9 @@ with TaskManager():
                 else:
                     N = 3
                     r = Grad(u) * Ps
-                    return CF(tuple([Trace(r[N*i:N*(i+1),0:N]) for i in range(N)]))       
-#        print("is killing? ", l2norm(divG(evecs2.MDComponent(0))-evecs2.MDComponent(0)))
+                    return CF(tuple([Trace(r[N * i:N * (i + 1), 0:N]) for i in range(N)]))
+                #        print("is killing? ", l2norm(divG(evecs2.MDComponent(0))-evecs2.MDComponent(0)))
+
         fun = GridFunction(VhG)
         kf = GridFunction(VhG)
         l2errs = []
@@ -260,37 +267,43 @@ with TaskManager():
         for k in range(3):
             lam = vlam[k]
             print(k)
-            print("lam: "+ str(abs(lam)))
-#                print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
-            if k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
-                temp.vec.data = evecs.vecs[k].data    
+            print("lam: " + str(abs(lam)))
+            #                print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
+            if k < maxk:  # abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
+                temp.vec.data = evecinter[k].vec.data
                 for j in range(maxk):
-                    
+
                     if True:
-                        kvfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kvfs[j],kvfs[j]), mesh=mesh, order=5)
-#                        kvfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(kvfs[j],kvfs[j]), mesh=mesh, order=5)
-
-                        scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs2.MDComponent(k),kvfs[j]), mesh=mesh, order=5)/kvfnorm
-#                        scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(evecs.MDComponent(k),kvfs[j]), mesh=mesh, order=5)/kvfnorm
-                        temp.vec.data += -scal*temparr[j].vec
-                        
-                l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(temp)**2, mesh=mesh, order=5)))
-                        
+                        kvfnorm = Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
+                                            cf=InnerProduct(kvfs[j], kvfs[j]), mesh=mesh, order=5)
+                        #                        kvfnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(kvfs[j],kvfs[j]), mesh=mesh, order=5)
+
+                        scal = Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
+                                         cf=kvfs[j]*evecinter[k] , mesh=mesh, order=5) / kvfnorm
+                        #                        scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=bh(evecs.MDComponent(k),kvfs[j]), mesh=mesh, order=5)/kvfnorm
+                        temp.vec.data += -scal * temparr[j].vec
+
+                l2errs.append(sqrt(
+                    Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=(temp) ** 2, mesh=mesh,
+                              order=5)))
+
         for k in range(3):
             lam = vlam[k]
             print(k)
-            print("lam: "+ str(abs(lam)))
-#                print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
-            if k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
-                    
-                
-                print("reverse filtering: ",sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(orthproj(temparr[k],evecinter))**2, mesh=mesh, order=5)))
-        return l2errs, evecs.MDComponent(0), evecs.MDComponent(1), evecs.MDComponent(2)
+            print("lam: " + str(abs(lam)))
+            #                print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
+            if k < maxk:  # abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
+                tempfun = orthproj( temparr[k],evecinter)
+                print("reverse filtering: ", sqrt(Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
+                                                            cf=(tempfun) ** 2, mesh=mesh,
+                                                            order=5)))
+        return l2errs, tempfun, evecinter[k], evecinter[k]#, evecs.MDComponent(2)
+
 
     l2errors = []
 
-    for i in range(n_cut_ref+1):
-        if i!=0:
+    for i in range(n_cut_ref + 1):
+        if i != 0:
             mesh.SetDeformation(None)
             lsetp1 = GridFunction(H1(mesh, order=2))
             InterpolateToP1(phi, lsetp1)
@@ -306,28 +319,27 @@ with TaskManager():
 
             mesh.SetDeformation(deformation)
 
-
         print("starting solve..")
         l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step=i)
         l2errors.append(l2err)
 
     print("L2 errors:")
     print(l2errors)
-if False:# and len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
-#    input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
+if True:  # and len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
+    #    input("Continue (press enter) to create a VTK-Output to stokestracefem3d.vtk")
 
-#    with TaskManager():
-    if maxk==1:
+    #    with TaskManager():
+    if maxk == 1:
         vtk = VTKOutput(ma=mesh,
                         coefs=[lset_approx, ev1, kvfs[0]],
                         names=["P1-levelset", "ev", "kf"],
-                        filename="killing", subdivision=0, legacy = False)
+                        filename="killing", subdivision=0, legacy=False)
 
         vtk.Do()
     else:
         vtk = VTKOutput(ma=mesh,
                         coefs=[lset_approx, ev1, ev2, ev3, kvfs[0], kvfs[1], kvfs[2]],
                         names=["P1-levelset", "ev1", "ev2", "ev3", "kf1", "kf2", "kf3"],
-                        filename="killing", subdivision=0, legacy = False)
+                        filename="killing", subdivision=0, legacy=False)
 
         vtk.Do()