Algebraic Multigrid for Compressed VectorH1 space

from netgen.geom2d import SplineGeometry from netgen.meshing import MeshingParameters from xfem import * from netgen.csg import CSGeometry, OrthoBrick, Pnt import netgen.gui from netgen.geom2d import unit_square from ngsolve import * from ngsolve.la import EigenValues_Preconditioner from xfem.lsetcurv import * def eps(u): return 0.5 * (Grad(u) + Grad(u).trans) def Setup(): cube = CSGeometry() cube.Add(OrthoBrick(Pnt(-1.5, -1.5, -1.5), Pnt(1.5, 1.5, 1.5))) mesh = Mesh(cube.GenerateMesh(maxh=1/12, quad_dominated=False)) lsetmeshadap = LevelSetMeshAdaptation(mesh, order=3, threshold=1000, discontinuous_qn=True) phi = Norm(CoefficientFunction((x, y, z))) - 1 deformation = lsetmeshadap.CalcDeformation(phi) lsetp1 = lsetmeshadap.lset_p1 mesh.SetDeformation(deformation) ci = CutInfo(mesh, lsetp1) Vhbase = VectorH1(mesh, order=2, dirichlet=[1,2,3,4,5,6]) Vhneg = Restrict(Vhbase, ci.GetElementsOfType(HASNEG)) Vhpos = Restrict(Vhbase, ci.GetElementsOfType(HASPOS)) WhGV = FESpace([Vhneg, Vhpos], flags={"dgjumps": True}) lset_neg = {"levelset": lsetp1, "domain_type": NEG} lset_pos = {"levelset": lsetp1, "domain_type": POS} lset_doms = [lset_neg, lset_pos] elements_doms = [HASNEG, HASPOS] u, v = WhGV.TnT() a = BilinearForm(WhGV) for i in [0, 1]: # loop over domains a += SymbolicBFI(lset_doms[i], form=2* InnerProduct(eps(u[i]), eps(v[i])), definedonelements=ci.GetElementsOfType(elements_doms[i])) a += SymbolicBFI(lset_doms[i], form=u[i] * v[i], definedonelements=ci.GetElementsOfType(elements_doms[i])) return mesh, WhGV, a mesh, WhGV, a = Setup() gfu = GridFunction(WhGV) gfu.components[1].Set(CoefficientFunction((x,y,z)),BND) mg = Preconditioner(a, "h1amg") # register mg to a a.Assemble() print('NDOF = ', WhGV.ndof) preI = Projector(mask=WhGV.FreeDofs(), range=True) lams = EigenValues_Preconditioner(mat=a.mat, pre=mg.mat) print("condition number", max(lams)/min(lams))