some changes

killing.py

@@ -14,13 +14,13 @@ def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsy
     uSol2 = (Ps * CF((-z ** 2, x, y)))
 
     for k in range(maxsymex):
+        print("removing something?")
         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},
                                 cf=InnerProduct(uSol, kvfs[k]), mesh=mesh, order=5) / kfexnorm * kvfs[k]
 
-    print(Integrate(levelset_domain={"levelset": phi, "domain_type": IF},
-                                cf=InnerProduct(uSol-uSol2, uSol-uSol2), mesh=mesh, order=5))
+
     Vh = VectorH1(mesh, order=order, dirichlet=[])
 
     ci = CutInfo(mesh, lset_approx)
@@ -28,6 +28,7 @@ def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsy
     VhG = Restrict(Vh, ba_IF)
 
     X = VhG 
+
     gfu = GridFunction(X)
 
     n = Normalize(grad(lset_approx))
@@ -38,9 +39,13 @@ def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsy
 #    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)
+#    print("number of elements: ")
+#    print(len(elvol))
+#    print("hmax:")
+#    print(cuthmax)
     alpha=order+1
     epsilon = h ** alpha
-    if False:
+    if True:
         epsilon=0
     print("Epsilon:")
     print(epsilon)
@@ -65,7 +70,7 @@ def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsy
 
     #        rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
 
-    if step > 4:
+    if step > 10:
         comp = True
     else:
         comp = False
@@ -76,7 +81,7 @@ def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsy
     def E(u):
         return (Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing
 
-    a = BilinearForm(X, symmetric=True)
+    a = BilinearForm(X, symmetric=False)
     a += InnerProduct(E(u), E(v)).Compile(realcompile=comp) * ds
     a += epsilon * (Pmat * u) * (Pmat * v) * ds
     a += (eta * ((u * ntilde) * (v * ntilde))) * ds
@@ -107,7 +112,6 @@ def getKilling(mesh, lset_approx, deformation,  vlam, evecs, step,alpha=3, maxsy
     f.Assemble()
     a.Assemble()
 
-
     gfu.vec.data = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="pardiso") * f.vec
     kf = GridFunction(VhG)
     print(vlam)
@@ -193,7 +197,7 @@ if __name__ == "__main__":
 
     cube = CSGeometry()
     if geom == "ellipsoid":
-        phi = (Norm(CoefficientFunction((x, y, z / c))) - args.radius).Compile(realcompile=True)
+        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))
@@ -255,7 +259,7 @@ if __name__ == "__main__":
                                            mesh=mesh) * Integrate({"levelset": lset_approx, "domain_type": IF},
                                                                   cf=x * y ** 3 + z, mesh=mesh))
 
-    weingex = grad(Normalize(grad(phi, mesh)), mesh).Compile(realcompile=True)
+    weingex = grad(Normalize(grad(phi, mesh)), mesh).Compile()#realcompile=True)
     l2errors = []
     l2errkf = []
     vlams = []
@@ -281,6 +285,7 @@ if __name__ == "__main__":
                 if math.log(abs(1 - vlams[k][1]), 2) <-(order+1)*(i):
                     print("corresponding log: ")
                     print(math.log(abs(1 - vlams[k][1]), 2))
+                    
                     print("condition: ", -(order+1)*(i+1))
                     symmetries += 1
         elif i >= 2:
@@ -293,10 +298,12 @@ if __name__ == "__main__":
 #                if math.log(abs(1 - vlams[j][-1]), 10) < math.log(abs(temp - vlams[j][-1]), 10):
                 print("corresponding log: ")
                 print(math.log(abs(1 - temp), 2))
-                print("condition: ", -(2*order+1)*(i))
-                if False or math.log(abs(1 - temp), 2) < -(2*order+1) * (i):
+                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:
 
                     symmetries += 1
+            
         if True and i!=0:
             print("starting solve..")
             print("symmetries: ", symmetries)

test.jl

+include("vectorlaplacian.jl")
+
+Vlaplace.solvelaplace()
\ No newline at end of file

vectorlaplacian.jl

-
+module Vlaplace
 using Gridap
 using GridapEmbedded
 using LinearAlgebra
 using GridapDistributed
 using GridapGmsh: gmsh
 using GridapGmsh
+using SuiteSparse
+#using GridapPETSc
  # Define geometry
- function main()
+ function blf(Pmat, n_Γd, nex2, dΓd, dΩ)
+  e = 1
+  h=0.3
+  η = 1.0/h^2
+  ρu = 1.0/h
+  println("setup blf")
+  #E(u) = Pmat2(Pmat((∇(u)+transpose(∇(u)))))
+  E(u) = Pmat⋅ symmetric_gradient(u)⋅Pmat
+  a1(u,v)= ∫( E(u)⊙ E(v) )*dΓd 
+  #a2(u,v) = ∫( e*((Pmat(u))⋅(Pmat(v))))*dΓd
+  a2(u,v) = ∫(e*(Pmat ⋅ u)⋅(Pmat ⋅v))*dΓd
+  a3(u,v) = ∫(η*((u ⋅ n_Γd)⋅(v⋅n_Γd)))*dΓd
+#  a4(u,v) =  ∫(ρu*((∇(u)⋅n_Γd)⋅(∇(v)⋅n_Γd)))*dΩ
+  #a4(u,v) =  ∫(ρu*((∇(u))⊙(∇(v))))*dΩ
+  a4(u,v)=∫(ρu*(∇(u)⋅ nex2)⋅(∇(v)⋅nex2))*dΩ
+  a(u,v) =(a1(u,v)+a2(u,v)+a3(u,v)+a4(u,v))
+      
+ end
+ function lf(Id)
+  println("setup lf")
+  e1 = VectorValue{3,Float64}(1,0,0)
+  e2 = VectorValue{3,Float64}(0,1,0)
+  e3 = VectorValue{3,Float64}(0,0,1)
+ #Ps = Pmats(n_Γd)
+nex(x) =1/norm(x)^2 *TensorValue{3,3,Float64}(x[1]*x[1],x[1]*x[2],x[1]*x[3],x[2]*x[1],x[2]^2,x[2]*x[3],x[3]*x[1],x[3]*x[2],x[3]^2)
+Ps(x) = Id(x)-nex(x)
+#  uex(x) = Ps(x) ⋅ VectorValue{3,Float64}(-x[3]^2,x[1],x[2])
+uex(x)=Ps(x) ⋅  VectorValue{3,Float64}(-x[3]^2,x[1],x[2])
+#  ∇uex(x) = ∇(uex(x))
+#  ∇(::typeof(uex))=∇uex
+ A(x) = (Ps(x) ⋅ symmetric_gradient(uex)(x)) ⋅ Ps(x) 
+#  A1(x) = A(x)⋅ e1(x)
+ #divgamma(x) = VectorValue{3,Float64}(tr(Ps(x) ⋅ ∇(A1)(x)),tr(Ps(x) ⋅ ∇(A⋅ e2)(x)),tr(Ps(x) ⋅ ∇(A⋅ e3)(x))) 
+A1(x) = A(x)⋅ e1
+ divgamma(x) = A1(x)
+ rhs(x) = divgamma(x) + uex(x)
+#  rhs(x)=uex(x)
+ end
+function domeshing()
+    order = 2
   gmsh.initialize()
   gmsh.clear()
   box = gmsh.model.occ.addBox(
-      -2.,-2.,-2.,2.,2.,2.)
+      -2.,-2.,-2.,4.,4.,4.)
    gmsh.model.occ.synchronize()
-   gmsh.model.mesh.setSize(gmsh.model.getEntities(0),0.1)
+   gmsh.model.mesh.setSize(gmsh.model.getEntities(0),0.5)
   gmsh.model.mesh.generate(3)
-
+gmsh.model.mesh.refine()
   gmsh.write("model.msh")
+#  bgmodel = gmsh.model.getEntities(0)
   gmsh.finalize()
      n = 10
    # Background model
@@ -24,14 +66,18 @@ using GridapGmsh
    pmax = 2*Point(1,1,1)
  #  bgmodel = simplexify(CartesianDiscreteModel(pmin,pmax,partition))
    bgmodel =GmshDiscreteModel("model.msh")
+return bgmodel
+end
   # Select geometry
+  function solvelaplace()
+    bgmodel = domeshing()
   R = 0.5
   ϕ(x) = norm(x)-1
   # Forcing data
   ud = 1
   f = 10
   println("doing cutter")
-lscutter = GridapEmbedded.LevelSetCutters.sphere(1)
+lscutter = GridapEmbedded.LevelSetCutters.sphere(1, x0=Point(0.0,0.0,0.0), name="sphere")
   # Cut the background model
   cutgeo = cut(bgmodel, lscutter)
 
@@ -41,65 +87,83 @@ lscutter = GridapEmbedded.LevelSetCutters.sphere(1)
 
   # Setup normal vectors
   n_Γd = get_normal_vector(Γd)
- 
-  #writevtk(Ω,"trian_O")
+#  reffe2 = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
+#  fe = FESpace(Ω, reffe2)
+#  n = interpolate_everywhere(n_Γd,fe)
+  writevtk(Ω,"trian_O")
 #  writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
   #writevtk(Γg,"trian_Gg",cellfields=["normal"=>n_Γg])
-  #writevtk(Triangulation(bgmodel),"bgtrian")
+  writevtk(Triangulation(bgmodel),"bgtrian")
 
   # Setup Lebesgue measures
-  order = 1
+order=2
   degree = 2*order
   dΩ = Measure(Ω,degree)
   dΓd = Measure(Γd,degree)
   reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, order)
-
+  println("temp?")
+  function temp(x)
+    return TensorValue{3,3, Float64}(1,0,0,0,1,0,0,0,1)
+  end
+  println("id?")
+  Id = temp
   # Setup FESpace
   Ω_act = Triangulation(cutgeo,ACTIVE)
-  normalspace = FESpace(Ω_act,reffe)
+  writevtk(Ω_act, "trian_O_act")
+  #normalspace = FESpace(Γd,reffe,conformity=:H1)
+  
 #  ntilde = interpolate_everywhere(n_Γd, normalspace) 
   println("setup TnT")
-  V = TestFESpace(Ω_act,ReferenceFE(lagrangian,VectorValue{3,Float64},order),conformity=:H1)
+  V = TestFESpace(Ω_act,ReferenceFE(lagrangian,VectorValue{3,Float64},order, space=:P),conformity=:H1)
   U = TrialFESpace(V)
 
   # Weak form
-    nex = ∇(ϕ)
-  h = (pmax - pmin)[1] / partition[1]
-  Pmats(n) = (x)->identity(x)-(n⊗ n) ⋅ x
-  Pmats2(n) = (x)-> identity(x)-x⋅ (n⊗n)
+ #   nex = ∇(ϕ)
+ nex2(x) = VectorValue{3,Float64}(x[1]/norm(x),x[2]/norm(x),x[3]/norm(x))
+  h = 0.3
+#  Pmats(n) = (x)->identity(x)-(n⊗ n) ⋅ x
+#  Pmats2(n) = (x)-> identity(x)-x⋅ (n⊗n)
+  Pmats(n) = Id-(n⊗ n)
   Pmat = Pmats(n_Γd)
-  Pmat2 = Pmats2(n_Γd)
-  
-  ϵ = 1
-  η = 1.0/h^2
-  ρu = 1.0/h
-  println("setup blf")
-  E(u) = Pmat2(Pmat((∇(u)+transpose(∇(u)))))
-  a1(u,v)= ∫( E(u)⊙ E(v) )*dΓd 
-  a2(u,v) = ∫( ϵ*((Pmat(u))⋅(Pmat(v))))*dΓd
-  a3(u,v) = ∫(η*((u ⋅ n_Γd)⋅(v⋅n_Γd)))*dΓd
- # a4(u,v) =  ∫(ρu*((∇(u)⋅n_Γd)⋅(∇(v)⋅n_Γd)))*dΩ
-  a4(u,v) =  ∫(ρu*((∇(u))⊙(∇(v))))*dΩ
-  
-  a(u,v) =(a1(u,v)+a2(u,v)+a3(u,v)+a4(u,v))
-      
-   
-println("setup lf")
-  l(v) =∫(VectorValue{3,Float64}(1,1,1)⋅ v)*dΓd
+ # Pmat2 = Pmats2(n_Γd)
+#  Pmat = (x)->x
+#  Pmat2 = (x)->x
+
+  # A = Gridap.Algebra.sparsecsr(a)
+
+  #l(v) =∫(Pmat⋅(VectorValue{3,Float64}(1,2,3))⋅ v)*dΓd
+  rhs = lf(Id)
+  a = blf(Pmat,n_Γd,nex2, dΓd,dΩ)
+  l(v)=∫(rhs ⋅ v)*dΓd
 #    ∫( v*f ) * dΩ +
 #    ∫( (γd/h)*v*ud - (n_Γd⋅∇(v))*ud ) * dΓd
 
   # FE problem
   println("setup FEop")
-   op = AffineFEOperator(a,l,U,V)
-   @time AffineFEOperator(a,l,U,V)
+    @time op = AffineFEOperator(a,l,U,V)
+    println("done setting up")
+   A3=get_matrix(op)
+   
+   f=get_vector(op)
+   print("factorize")
+   @time qra=factorize(A3)
+   uh = qra\f
+#   uh=A\f
+   println(length(f))
+
+    gh = FEFunction(V,uh,get_dirichlet_dof_values(V))
+#   @time AffineFEOperator(a,l,U,V)
+
   println("solve")
-  ls = LinearFESolver()
-  uh = solve(op)
+ # println(sqrt(length(A)))
+  #ls = LinearFESolver()#PETScLinearSolver())
+#  uh = solve(op)#, ls)
 
   # Postprocess
  # if outputfile !== nothing
-    writevtk(Ω,"vectorlaplace_julia",cellfields=["uh"=>uh])
-  end
-  mesh_partition=(4,4)
-main()
\ No newline at end of file
+    writevtk(Ω,"vectorlaplace_julia",order=2,cellfields=["uh"=>gh])
+    
+
+end
+
+end
\ No newline at end of file