midnight changes

statstokes.py

@@ -46,7 +46,7 @@ order = 2
 reac_cf = 1
 diff_cf = 1
 # Ellipsoid parameter
-c = 1.5
+c = 1.1
 # Geometry
 #SetHeapSize(10**10)
 with TaskManager():
@@ -166,8 +166,8 @@ with TaskManager():
 
         h = specialcf.mesh_size
         elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
-        cuthmax = max([(6*vol)**(1/3) for vol in elvol])
-        alpha=3/2
+        cuthmax = max([(2*vol)**(1/2) for vol in elvol])
+        alpha=2
         epsilon =(cuthmax)**(alpha)
         print("Epsilon:")
         print(epsilon)
@@ -279,19 +279,30 @@ with TaskManager():
         print("removing killing VFs")
         print(vlam)
         fun = GridFunction(VhG)
+        kf = GridFunction(VhG)
         for k, lam in enumerate(vlam):
             print("lam: "+ str(abs(lam)))
             print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
-            if abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
-                
+            if k==0:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
+                kf.vec.data =evecs.vecs[0].data
+                evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5)
                 print("removing number: "+str(k))
-                integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh)
+                integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh, order=5)/evecnorm
 
-                gfu.components[0].vec.data += integral*evecs.vecs[k]
+                print("scalar product before removing killing field: " + str(integral))
+#               integral = -Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(gfu.components[0],evecs.MDComponent(k)), mesh=mesh)
+                gfu.components[0].vec.data+=  integral*evecs.vecs[k]
+#                print(gfu.components[0].vec.data-kf.vec.data)
+                print("scalar product with killing field: ")
+
+        for k in range(3):
+            print(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF},
+                            cf = InnerProduct(gfu.components[0], evecs.MDComponent(k)), mesh=mesh, order=10))
         #gfu2.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso")*f2.vec
         l2err = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(gfu.components[0] - uSol) ** 2, mesh=mesh))
+        uerr = gfu.components[0] - uSol
         print("l2error : ", l2err)
-        return l2err
+        return l2err, uerr, gfu, kf
 
 
 
@@ -315,16 +326,17 @@ with TaskManager():
 
 
         print("starting solve..")
-        l2errors.append(solveStokes(mesh, lset_approx, deformation, step=i))
+        l2err, uerr, uh, kf = solveStokes(mesh, lset_approx, deformation, step=i)
+        l2errors.append(l2err)
     print(l2errors)
-if len(sys.argv) < 2 or not (sys.argv[1] == "testmode"):
+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")
 
 #    with TaskManager():
-    pass
-    """vtk = VTKOutput(ma=mesh,
-                    coefs=[lset_approx, gfu.components[0], uSol, uerr],
-                    names=["P1-levelset", "uh", "uex", "uerr"],
-                    filename="stokestracefem3d", subdivision=2)
 
-    vtk.Do()"""
+    vtk = VTKOutput(ma=mesh,
+                    coefs=[lset_approx, uh.components[0],kf, uerr],
+                    names=["P1-levelset", "uh","KF", "uerr"],
+                    filename="stokestracefem3d", subdivision=2, legacy = False)
+
+    vtk.Do()