initial commit

statstokes.py

+from netgen.csg import CSGeometry, OrthoBrick, Pnt
+from ngsolve import *
+from ngsolve.internal import *
+from xfem import *
+from xfem.lsetcurv import *
+from math import pi
+import datetime
+begin_time=datetime.datetime.now()
+with TaskManager():
+
+    # -------------------------------- PARAMETERS ---------------------------------
+    # Mesh diameter
+    maxh = 0.6
+    geom = "torus"
+    # Polynomial order of FE space
+    order = 2
+
+    # Problem parameters
+    reac_cf = 1
+    diff_cf = 1
+
+    # Geometry
+    cube = CSGeometry()
+    if geom == "circle":
+        phi = Norm(CoefficientFunction((x, y, z))) - 1
+
+        cube.Add(OrthoBrick(Pnt(-1.41, -1.41, -1.41), Pnt(1.41, 1.41, 1.41)))
+        mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
+    elif geom=="decocube":
+        phi=(x**2+y**2-4)**2+(y**2-1)**2+(y**2+z**2-4)**2+(x**2-1)**2+(x**2+z**2-4)**2+(z**2-1)**2-13
+        cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
+        mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
+    elif geom=="torus":
+        phi = sqrt(z**2+(sqrt(x**2+y**2)-1)**2)-0.5
+        cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
+        mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
+    n_cut_ref = 3
+    for i in range(n_cut_ref):
+        lsetp1 = GridFunction(H1(mesh, order=1))
+        InterpolateToP1(phi, lsetp1)
+        RefineAtLevelSet(lsetp1)
+        mesh.Refine()
+    lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order)
+    deformation = lsetmeshadap.CalcDeformation(phi)
+    lset_approx = lsetmeshadap.lset_p1
+
+    mesh.SetDeformation(deformation)
+    # Class to compute the mesh transformation needed for higher order accuracy
+    #  * order: order of the mesh deformation function
+    #  * threshold: barrier for maximum deformation (to ensure shape regularity)
+
+    # Background FESpace
+    Vh = VectorH1(mesh, order=order, dirichlet=[])
+    Qh = H1(mesh, order=order-1, dirichlet = [])
+    ci = CutInfo(mesh, lset_approx)
+    ba_IF = ci.GetElementsOfType(IF)
+    VhG = Restrict(Vh, ba_IF)
+    L2G = Restrict(Qh, ba_IF)
+    print("pre product space")
+    X=VhG*L2G
+    gfu = GridFunction(X)
+    lset_approx2 = Interpolate(phi,VhG)
+
+    # Coefficients / parameters:
+    def coeffgrad(u):
+        var = [x,y,z]
+        temp = []
+        for v in var:
+            for k in range(3):
+                temp.append(u[k].Diff(v))
+        return Pmat*CoefficientFunction(tuple(temp), dims=(3,3))
+    def coeffscalgrad(n):
+        temp = []
+        var = [x,y,z]
+        for v in var:
+            temp.append(n.Diff(v))
+        return CoefficientFunction(tuple(temp))
+    # Tangential projection
+    n = Normalize(grad(lset_approx))
+#    n = Normalize(coeffscalgrad(lset_approx2))#/Norm(coeffscalgrad(lset_approx2))
+    n.Compile()
+    print(n.dim)
+    h = specialcf.mesh_size
+
+    eta = 100.0 / (h * h)
+    rhou = 1.0 / h
+    rhop = h
+    # define solution and right-hand side
+    Pmat = Id(3) - OuterProduct(n, n)
+    if geom=="circle":
+        functions = {
+            "extvsol1": ((-y - z) * x + y * y + z * z),
+            "extvsol2": ((-x - z) * y + x * x + z * z),
+            "extvsol3": ((-x - y) * z + x * x + y * y),
+            "rhs1": -((y + z) * x - y * y - z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
+            "rhs2": ((-x - z) * y + x * x + z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
+            "rhs3": ((-x - y) * z + x * x + y * y) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
+        }
+        extpsol = (x*y**3+z*(x**2+y**2+z**2)**(3/2))/((x**2+y**2+z**2)**2)
+        uSol = CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
+        pSol = CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
+#        extpsol = CoefficientFunction(0)
+#        extpsol = CoefficientFunction((x*y**3+z*(x**2+y**2+z**2)**(3/2))/(x**2+y**2+z**2)**2)
+        rhs1 = CoefficientFunction((functions["rhs1"], functions["rhs2"], functions["rhs3"]))
+    elif geom=="decocube":
+        uSol=CoefficientFunction((-z**2,y,x))
+        pSol = CoefficientFunction(x*y**3+z-1/Integrate({"levelset": lset_approx, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh)*Integrate({"levelset": lset_approx, "domain_type":IF},cf=x*y**3+z,mesh=mesh))
+        extpsol = pSol
+        extpsol.Compile()
+#        extpsol = CoefficientFunction(1)
+    elif geom=="torus":
+        rhs1 = CoefficientFunction((0,0,0))
+        rhs2 = CoefficientFunction(x+y)
+    u,p = X.TrialFunction()
+    v, q = X.TestFunction()
+    # Measure on surface
+    ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
+    # Measure on the bulk around the surface
+    dx = dx(definedonelements=ba_IF, deformation=deformation)
+
+
+
+
+    def coeffgradtrans(u):
+        return coeffgrad(u).trans
+    def tangsymgrad(u):
+        return 0.5*Pmat*(coeffgrad(u)+coeffgradtrans(u))
+    def tangdiv(u):
+        b = Pmat*tangsymgrad(u)
+        return b[0,0]+b[1,1]+b[2,2]
+    def trace(A):
+        return A[0,0]+A[1,1]+A[2,2]
+    def divtens(A):
+        result = []
+        for k in range(3):
+            vec = A[k,0:3]
+            gradvec = coeffgrad(vec)
+            result.append(gradvec[0,0]+gradvec[1,1]+gradvec[2,2])
+
+        return CoefficientFunction(tuple(result))
+
+
+    ntilde = coeffscalgrad(lset_approx2) / Norm(coeffscalgrad(lset_approx2))
+    ntilde = n
+    weing =specialcf.Weingarten(3)
+    weing = coeffgrad(n)
+    print("start rhs")
+    if geom!="torus":
+    #    rhs1 = -Pmat*divtens(tangsymgrad(uSol)-(uSol*n)*weing)+Pmat*uSol+Pmat*CoefficientFunction((extpsol.Diff(x),extpsol.Diff(y),extpsol.Diff(z)))
+        rhs1 = -Pmat*divtens(tangsymgrad(uSol))+Pmat*uSol+Pmat*CoefficientFunction((extpsol.Diff(x),extpsol.Diff(y),extpsol.Diff(z)))
+
+        rhs1.Compile()
+        rhs2 = trace(coeffgrad(uSol))
+        rhs2.Compile()
+        print("end rhs")
+#    rhs2 = 0
+    # bilinear forms:
+    a = BilinearForm(X, symmetric=True)
+#    a += (InnerProduct(0.5*Pmat * Sym(grad(u)) * Pmat-(u*n)*weing, 0.5*Pmat * Sym(grad(v)) * Pmat-(u*n)*weing)) * ds
+    a += (InnerProduct(0.5*Pmat * Sym(grad(u)) * Pmat, 0.5*Pmat * Sym(grad(v)) * Pmat)) * ds
+
+    a += (Pmat*u) * (Pmat*v) * ds
+    a += (eta * ((u * ntilde) * (v * ntilde))) * ds
+    a += rhou * ((grad(u) * n) * (grad(v) * n)) * dx
+    a += InnerProduct(u,Pmat*grad(q))*ds
+    a += InnerProduct(v, Pmat*grad(p)) * ds
+    a += -rhop*InnerProduct(n*grad(p),n*grad(q))*dx
+    print("preassembling")
+    a.Assemble()
+    print("assembled a")
+    f = LinearForm(X)
+    f += rhs1 * v * ds
+    f += -rhs2*q*ds
+    f.Assemble()
+    print("assembled f, inverting")
+    gfu.vec.data = a.mat.Inverse(freedofs = X.FreeDofs(), inverse="pardiso") * f.vec
+    if geom != "torus":
+        print("computing error")
+        print("l2error : ", sqrt(Integrate((gfu.components[0] - uSol)**2 * ds, mesh=mesh)))
+        uerr = (gfu.components[0]-uSol)
+    Draw(deformation, mesh, "deformation")
+    Draw(gfu, mesh, "u")
+
+    if (len(sys.argv) < 2 or not (sys.argv[1] == "testmode")) and True:
+#        input("Continue (press enter) to create a VTK-Output to tracefem3d.vtk")
+
+        print("writing output")
+        if geom != "torus":
+            vtk = VTKOutput(ma=mesh,
+                        coefs=[lset_approx, gfu.components[0], uSol, uerr],
+                        names=["P1-levelset", "uh", "uex", "uerr"],
+                        filename="stokespenalty3d", subdivision=2)
+        else:
+            vtk = VTKOutput(ma=mesh,
+                            coefs=[lset_approx, gfu.components[0]],
+                            names=["P1-levelset", "uh"],
+                            filename="stokespenalty3d", subdivision=2)
+
+    vtk.Do()
+print(datetime.datetime.now()-begin_time)
+

vectorlaplacian.py

+from netgen.csg import CSGeometry, OrthoBrick, Pnt
+from ngsolve import *
+from ngsolve.internal import *
+from xfem import *
+from xfem.lsetcurv import *
+from math import pi
+with TaskManager():
+    # -------------------------------- PARAMETERS ---------------------------------
+    # Mesh diameter
+    maxh = 0.4
+    geom ="circle"
+    # Polynomial order of FE space
+    order = 2
+
+    # Problem parameters
+    reac_cf = 1
+    diff_cf = 1
+
+    # Geometry
+    cube = CSGeometry()
+    if geom == "circle":
+        phi = Norm(CoefficientFunction((x, y, z))) - 1
+
+        cube.Add(OrthoBrick(Pnt(-1.41, -1.41, -1.41), Pnt(1.41, 1.41, 1.41)))
+        mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
+    elif geom=="decocube":
+        phi=(x**2+y**2-4)**2+(y**2-1)**2+(y**2+z**2-4)**2+(x**2-1)**2+(x**2+z**2-4)**2+(z**2-1)**2-13
+        cube.Add(OrthoBrick(Pnt(-3, -3, -3), Pnt(3, 3, 3)))
+        mesh = Mesh(cube.GenerateMesh(maxh=maxh, quad_dominated=False))
+    n_cut_ref = 2
+    for i in range(n_cut_ref):
+        lsetp1 = GridFunction(H1(mesh, order=2))
+        InterpolateToP1(phi, lsetp1)
+        RefineAtLevelSet(lsetp1)
+        mesh.Refine()
+    lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order)
+    deformation = lsetmeshadap.CalcDeformation(phi)
+#    deformation = None
+    lset_approx = lsetmeshadap.lset_p1
+#    lsetp1 = GridFunction(H1(mesh, order=1))
+#    InterpolateToP1(phi, lsetp1)
+
+#    lset_approx = lsetp1
+    mesh.SetDeformation(deformation)
+    # Class to compute the mesh transformation needed for higher order accuracy
+    #  * order: order of the mesh deformation function
+    #  * threshold: barrier for maximum deformation (to ensure shape regularity)
+
+    # Background FESpace
+    Vh = VectorH1(mesh, order=order, dirichlet=[])
+
+    ci = CutInfo(mesh, lset_approx)
+    ba_IF = ci.GetElementsOfType(IF)
+    VhG = Restrict(Vh, ba_IF)
+
+    gfu = GridFunction(VhG)
+
+    # Coefficients / parameters:
+
+
+    # Tangential projection
+    n = Normalize(grad(lset_approx))
+    print(n)
+    h = specialcf.mesh_size
+
+    eta = 100.0 / (h * h)
+    rho = 1.0 / h
+
+    # define solution and right-hand side
+    functions = {
+        "extvsol1": ((-y - z) * x + y * y + z * z),
+        "extvsol2": ((-x - z) * y + x * x + z * z),
+        "extvsol3": ((-x - y) * z + x * x + y * y),
+        "rhs1": -((y + z) * x - y * y - z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
+        "rhs2": ((-x - z) * y + x * x + z * z) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
+        "rhs3": ((-x - y) * z + x * x + y * y) * (x * x + y * y + z * z + 1) / (x * x + y * y + z * z),
+    }
+    if geom=="circle":
+        uSol = CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
+        rhs = CoefficientFunction((functions["rhs1"], functions["rhs2"], functions["rhs3"]))
+    elif geom=="decocube":
+        uSol = CoefficientFunction((-(z ** 2), y, x))
+
+    u, v = VhG.TnT()
+
+    # Measure on surface
+    ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
+    # Measure on the bulk around the surface
+    dx = dx(definedonelements=ba_IF, deformation=deformation)
+
+    Pmat = Id(3) - OuterProduct(n, n)
+    def coeffgrad(u):
+        var = [x,y,z]
+        temp = []
+        for v in var:
+            for k in range(3):
+                temp.append(u[k].Diff(v))
+        return Pmat*CoefficientFunction(tuple(temp), dims=(3,3))*Pmat
+    def coeffgradtrans(u):
+        return coeffgrad(u).trans
+    def tangsymgrad(u):
+        return (coeffgrad(u)+coeffgradtrans(u))
+    def tangdiv(u):
+        b = Pmat*tangsymgrad(u)
+        return b[0,0]+b[1,1]+b[2,2]
+    def divtens(A):
+        result = []
+        for k in range(3):
+            vec = A[k,0:3]
+            gradvec = coeffgrad(vec)
+            result.append(gradvec[0,0]+gradvec[1,1]+gradvec[2,2])
+
+        return CoefficientFunction(tuple(result))
+#    rhs = -Pmat*divtens(tangsymgrad(uSol))+Pmat*uSol
+    rhs.Compile()
+    # bilinear forms:
+    a = BilinearForm(VhG, symmetric=True)
+    a += (InnerProduct(Pmat * Sym(grad(u)) * Pmat, Pmat * Sym(grad(v)) * Pmat)) * ds
+    a += InnerProduct(u,  v) * ds
+    a += (eta * ((u * n) * (v * n))) * ds
+    a += rho * ((grad(u) * n) * (grad(v) * n)) * dx
+    a.Assemble()
+
+    f = LinearForm(VhG)
+    f += rhs * v * ds
+    f.Assemble()
+    gfu.vec[:]=0.0
+    gfu.vec.data = a.mat.Inverse(VhG.FreeDofs(), inverse="pardiso") * f.vec
+
+    print("l2error : ", sqrt(Integrate((gfu - uSol)**2 * ds, mesh=mesh)))
+#    Draw(deformation, mesh, "deformation")
+    Draw(gfu, mesh, "u")
+    uerr = gfu-uSol
+    if (len(sys.argv) < 2 or not (sys.argv[1] == "testmode")) and False:
+        input("Continue (press enter) to create a VTK-Output to tracefem3d.vtk")
+
+
+        vtk = VTKOutput(ma=mesh,
+                        coefs=[lset_approx, deformation, gfu, uSol, uerr, rhs],
+                        names=["P1-levelset", "deform", "uh", "uex", "uerr", "rhs"],
+                        filename="vectorlaplacepenalty3d", subdivision=2)
+        vtk.Do()
+
+