Merge remote-tracking branch 'origin/main'

helperfunctions.py

+from ngsolve import *
+from xfem.lsetcurv import *
+import matplotlib.pyplot as plt
+from netgen.csg import CSGeometry, OrthoBrick, Pnt
+import math
+def refineMesh(mesh, phi, order, hs=10 ** 6):
+    if order>2:
+        hs=10**8
+    with TaskManager():
+        mesh.SetDeformation(None)
+        lsetp1 = GridFunction(H1(mesh, order=order))
+        InterpolateToP1(phi, lsetp1)
+        RefineAtLevelSet(lsetp1)
+
+        mesh.Refine()
+
+        lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
+
+        deformation = lsetmeshadap.CalcDeformation(phi)
+        mesh.SetDeformation(deformation)
+        lset_approx = lsetmeshadap.lset_p1
+    return lset_approx, deformation
+
+
+def rhsfromsol(uSol, pSol, phi, Ps, mesh, ba_IF, weingex, interpol=True, Lpmat=None, Lpvec=None):
+    if interpol:
+
+        tmpeps = GridFunction(Lpmat)
+        tmpeps.Set(eps(uSol, Ps, mesh) - (uSol * Normalize(grad(phi, mesh))) * weingex, definedonelements=ba_IF)
+        tmptotal = GridFunction(Lpvec)
+        tmptotal.Set(Ps * (-divG(tmpeps, Ps, mesh) + grad(pSol, mesh)), definedonelements=ba_IF)
+        return tmptotal
+    else:
+
+        return Ps * (-divG(eps(uSol, Ps).Compile())) + Ps * grad(pSol, mesh).Compile()
+
+
+def coef_grad(u, mesh):
+    dirs = {1: [x], 2: [x, y], 3: [x, y, z]}
+    if u.dim == 1:
+        return CF(tuple([u.Diff(r) for r in dirs[mesh.dim]]))
+    else:
+        return CF(tuple([u[i].Diff(r) for i in range(u.dim) for r in dirs[mesh.dim]]), dims=(u.dim, mesh.dim))
+
+
+def grad(u, mesh=None):
+    if type(u) in [ProxyFunction, GridFunction]:
+        return ngsolve.grad(u)
+    else:
+        return coef_grad(u, mesh)
+
+
+def Pmats(n):
+    return Id(3) - OuterProduct(n, n)
+
+
+# Helper Functions for rate-of-strain-tensor
+def eps(u, Ps, mesh):
+    return Ps * Sym(grad(u, mesh)) * Ps
+
+
+def divG(u, Ps, mesh):
+    if u.dim == 3:
+        return Trace(grad(u, mesh) * Ps)
+    if u.dim == 9:
+        if u.dims[0] == 3:
+            N = 3
+            divGi = [divG(u[i, :]) for i in range(N)]
+            return CF(tuple([divGi[i] for i in range(N)]))
+        else:
+            N = 3
+            r = Grad(u) * Ps
+            return CF(tuple([Trace(r[N * i:N * (i + 1), 0:N]) for i in range(N)]))
+
+
+def l2sca(f, g, phi, mesh):
+    return Integrate(levelset_domain={"levelset": phi, "domain_type": IF}, cf=InnerProduct(f, g), mesh=mesh)
+
+
+def l2norm(f, phi, mesh):
+    return sqrt(l2sca(f, f, phi, mesh))
+
+
+def removekf(VhG, f, kvs, vlams, phi, mesh, maxk=3, tol=0, autofilter=True, inplace=False):
+    print(maxk)
+    for k in range(3):
+        if (k < maxk and autofilter) or False:
+            print(l2norm(kvs[k], phi, mesh))
+            sca = l2sca(f.components[0], kvs[k], phi, mesh)
+            f.components[0].vec.data += -sca * kvs[k].vec
+    return
+def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
+
+
+    hvals = []
+    for k in range(len(l2errs)):
+        hvals.append((1 / 2 ** 3) ** k)
+    plt.semilogy(l2errs, label=geom+" with KF removal")
+    plt.semilogy(l2errskf, label=geom + " without KF removal")
+
+    plt.semilogy(hvals, label="O(h^"+ str(order+1)+")")
+    plt.ylabel("L2 error")
+    plt.xlabel("Refinement")
+    plt.title("Convergence removing KF")
+    plt.legend()
+    plt.xticks(range(len(l2errs)))
+    if saveplot:
+        plt.savefig("convergence_plot_"+geom+"_c"+str(c)+"_order"+str(order)+"_refs"+str(refs)+".png")
+    else:
+        plt.show()
+def computeEV(phi, mesh, lset_approx, deformation, order=2):
+
+    vlams = []
+
+    Vh = VectorH1(mesh, order=order, dirichlet=[])
+
+    ci = CutInfo(mesh, lset_approx)
+    ba_IF = ci.GetElementsOfType(IF)
+    VhG = Restrict(Vh, ba_IF)
+    n = Normalize(grad(lset_approx))
+    ntilde = GridFunction(VhG)
+    ntilde.Set(Normalize(grad(phi, mesh)), definedonelements=ba_IF)
+    #        n=ntilde
+    h = specialcf.mesh_size
+    #        epsilon=1
+    eta = 1 / (h) ** 2
+    rhou = 1.0 / h
+    rhob = h
+    # 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)
+
+    u2 = VhG.TrialFunction()
+    v2 = VhG.TestFunction()
+
+    # Weingarten mappings
+
+    ninter = GridFunction(VhG)
+    ninter.Set(n, definedonelements=ba_IF)
+    weing = grad(ninter)
+
+    Pmat = Pmats(n)#.Compile()
+
+    def E(u):
+        return ((Pmat * Sym(grad(u)) * Pmat) - (u * n) * weing)#.Compile()
+
+    a2 = BilinearForm(VhG, symmetric=True)
+    a2 += InnerProduct(E(u2), E(v2)) * ds
+    a2 += (Pmat * u2) * (Pmat * v2) * ds
+    a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
+    a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
+    pre = Preconditioner(a2, "direct")
+    a2.Assemble()
+
+    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, 3, GramSchmidt=True)