changed iteration

helperfunctions.py

@@ -136,19 +136,22 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
         blocks = []
         freedofs = VhG.FreeDofs()
         class SymmetricGS(BaseMatrix):
-            def __init__ (self, smoother):
+            def __init__ (self, smoother, inv):
                 super(SymmetricGS, self).__init__()
+                self.k=0
                 self.smoother = smoother
+
             def Mult (self, x, y):
                 y[:] = 0.0
                 self.smoother.Smooth(y, x)
                 self.smoother.SmoothBack(y,x)
+
             def Height (self):
                 return self.smoother.height
             def Width (self):
                 return self.smoother.height
         
-        vertexdofs = BitArray(Vh.ndof)
+        vertexdofs = BitArray(VhG.ndof)
         vertexdofs[:] = False
         print("doing stuff for precond")
         """
@@ -162,11 +165,13 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
                             
             
 #            vdofs = set(d for d in VhG.GetDofNrs(el) if freedofs[d])
-            for d in Vh.GetDofNrs(v):
+            for d in VhG.GetDofNrs(v):
                 vertexdofs[d] = True
     
 #            blocks.append(vdofs)
+        
         blocks = [set(d  for el in mesh[v].elements for d in VhG.GetDofNrs(el) if freedofs[d]) for v in mesh.vertices]
+        
 #        blocks = [reduce(lambda x,y: x or y,[set(d for d in VhG.GetDofNrs(el) if freedofs[d])for el in mesh[v].elements])   for v in mesh.vertices]
 
         def blockcr(FEspace):
@@ -213,24 +218,25 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
         a2 += (Pmat * u2) * (Pmat * v2) * ds
         a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
         a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
-        pre1 = Preconditioner(a2, "direct", inverse="pardiso")
-
+        pre1 = Preconditioner(a2, "direct", inverse="sparsecholesky")
+#        prejac = Preconditioner(a2, "local", inverse="pardiso")
         a2.Assemble()
-        coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
+#        coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
 
 
 #        pre = Projector(VhG.FreeDofs(), True)
 #        pre1 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
-#        pre1 = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True))
-        pre = pre1#+coarsepre
-#        pre = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="pardiso")
+#        inv = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="sparsecholesky")
+#        prejac = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True), inv)
+        pre = pre1#+prejac
+
         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,maxit = 15)
+        return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True,maxit = 20)
 def getevproblem(mesh, phi, lset_approx, deformation, order=2):
     Vh = VectorH1(mesh, order=order, dirichlet=[])
     ci = CutInfo(mesh, lset_approx)

killing.py

@@ -300,7 +300,7 @@ if __name__ == "__main__":
                 print(math.log(abs(1 - temp), 2))
                 print("condition: ", -(2*order+1)*(i+1))
                 print("log10: ", math.log(abs(1 - temp), 10))
-                if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i+1)+1:
+                if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i+1)+3:
 
                     symmetries += 1