Merge branch 'main' of git.rwth-aachen.de:tilmandiego.aleman/reproduction into main

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")
@@ -163,12 +166,17 @@ 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)
         
+<<<<<<< HEAD
+        blocks = [set(d  for el in mesh[v].elements for d in VhG.GetDofNrs(el) if freedofs[d]) for v in mesh.vertices]
+        
+=======
         blocks = [set(VhG.GetDofNrs(el))  for el in VhG.Elements(VOL)]
+>>>>>>> ea4519effe74102893ea15789e6c4633f4f61bbf
 #        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):
@@ -216,26 +224,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, "bddc", inverse="pardiso")
-#        print(reduce(lambda x,y:x and y, vertexdofs and VhG.FreeDofs()))
+        pre1 = Preconditioner(a2, "direct", inverse="sparsecholesky")
+#        prejac = Preconditioner(a2, "local", inverse="pardiso")
         a2.Assemble()
 #        coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
 
 
 #        pre = Projector(VhG.FreeDofs(), True)
-#        pre2 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
-#        pre1 = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True))
-#        pre = pre1+coarsepre#+pre2+coarsepre
-        
-        #pre=SymmetricGS(pre1)
-        pre = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="pardiso")
+#        pre1 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
+#        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, 5, GramSchmidt=True,maxit = 10)
+        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)