Skip to content
GitLab
Commit adbcd4ad authored by Tilman Aleman's avatar Tilman Aleman
Browse files

working rhs, slow lf assemble

parent 3f94a648
Hide whitespace changes
Inline Side-by-side
statstokes.py
View file @ adbcd4ad
import sys
from netgen.csg import CSGeometry, OrthoBrick, Pnt
from ngsolve import *
from ngsolve.internal import *
......@@ -7,24 +9,24 @@ from math import pi
import datetime
begin_time=datetime.datetime.now()
with TaskManager():
#if True:
# -------------------------------- PARAMETERS ---------------------------------
# Mesh diameter
maxh = 0.6
geom = "torus"
geom = "circle"
# Polynomial order of FE space
order = 2
# Problem parameters
reac_cf = 1
diff_cf = 1
c=1.3
# Geometry
cube = CSGeometry()
if geom == "circle":
phi = Norm(CoefficientFunction((x, y, z))) - 1
phi = Norm(CoefficientFunction((x/c, y, z))) - 1
cube.Add(OrthoBrick(Pnt(-1.41, -1.41, -1.41), Pnt(1.41, 1.41, 1.41)))
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
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
......@@ -36,7 +38,7 @@ with TaskManager():
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))
lsetp1 = GridFunction(H1(mesh, order=2))
InterpolateToP1(phi, lsetp1)
RefineAtLevelSet(lsetp1)
mesh.Refine()
......@@ -59,22 +61,17 @@ with TaskManager():
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))
ntilde = Interpolate(coeffscalgrad(phi),VhG)
# Tangential projection
n = Normalize(grad(lset_approx))
# n = Normalize(coeffscalgrad(lset_approx2))#/Norm(coeffscalgrad(lset_approx2))
......@@ -85,8 +82,31 @@ with TaskManager():
eta = 100.0 / (h * h)
rhou = 1.0 / h
rhop = h
def coef_grad(u):
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):
if type(u) in [ProxyFunction, GridFunction]:
return ngsolve.grad(u)
else:
return coef_grad(u)
def Pmats(n):
return Id(3) - OuterProduct(n, n)
Ps = Pmats(Normalize(grad(phi)))
# define solution and right-hand side
Pmat = Id(3) - OuterProduct(n, n)
if geom=="circle":
functions = {
"extvsol1": ((-y - z) * x + y * y + z * z),
......@@ -97,17 +117,17 @@ with TaskManager():
"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"]))
uSol = Ps*CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
pSol = CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
# extpsol = CoefficientFunction(0)
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"]))
# rhs1 = CoefficientFunction((functions["rhs1"], functions["rhs2"], functions["rhs3"]))
elif geom=="decocube":
uSol=CoefficientFunction((-z**2,y,x))
uSol=Ps*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)
# extpsol = CoefficientFunction(0)
elif geom=="torus":
rhs1 = CoefficientFunction((0,0,0))
rhs2 = CoefficientFunction(x+y)
......@@ -119,46 +139,40 @@ with TaskManager():
dx = dx(definedonelements=ba_IF, deformation=deformation)
def eps(u):
return Ps * Sym(grad(u)) * Ps
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))
def divG(u):
if u.dim == 3:
return Trace(grad(u) * Ps)
if u.dim == 9:
N = 3
divGi = [divG(u[i, :]) for i in range(N)]
return CF(tuple([divGi[i] for i in range(N)]))
ntilde = coeffscalgrad(lset_approx2) / Norm(coeffscalgrad(lset_approx2))
ntilde = n
weing =specialcf.Weingarten(3)
weing = coeffgrad(n)
weing = grad(n)
weingex = grad(Normalize(grad(phi)))
# weing.Compile()
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 = Ps*-divG(eps(uSol)-(uSol*Normalize(grad(phi)))*weingex)+uSol+grad(extpsol)
rhs1.Compile()
rhs2 = trace(coeffgrad(uSol))
rhs2.Compile()
# rhs1.Compile()
rhs2 = Trace(Ps*grad(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
Pmat=Pmats(n)
def E(u):
return Pmat*Sym(grad(u))*Pmat-(u*n)*weing
a = BilinearForm(X, symmetric=True)
# a += (InnerProduct(Pmat * Sym(grad(u)) * Pmat-(u*n)*weing, Pmat * Sym(grad(v)) * Pmat-(v*n)*weing)) * ds
# a += (InnerProduct(Pmat * Sym(grad(u)) * Pmat, Pmat * Sym(grad(v)) * Pmat)) * ds
a += InnerProduct(E(u), E(v))*ds
a += (Pmat*u) * (Pmat*v) * ds
a += (eta * ((u * ntilde) * (v * ntilde))) * ds
a += rhou * ((grad(u) * n) * (grad(v) * n)) * dx
......@@ -178,10 +192,11 @@ with TaskManager():
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:
print(sys.argv[1])
if (sys.argv[1]=='1'):
# input("Continue (press enter) to create a VTK-Output to tracefem3d.vtk")
print("writing output")
......@@ -196,6 +211,6 @@ with TaskManager():
names=["P1-levelset", "uh"],
filename="stokespenalty3d", subdivision=2)
vtk.Do()
vtk.Do()
print(datetime.datetime.now()-begin_time)
vectorlaplacian.py
View file @ adbcd4ad
......@@ -4,30 +4,33 @@ from ngsolve.internal import *
from xfem import *
from xfem.lsetcurv import *
from math import pi
import pickle
import datetime
with TaskManager():
# -------------------------------- PARAMETERS ---------------------------------
# Mesh diameter
maxh = 0.4
geom ="circle"
maxh = 0.6
geom ="decocube"
# Polynomial order of FE space
order = 2
# Problem parameters
reac_cf = 1
diff_cf = 1
c=1.2
# Geometry
cube = CSGeometry()
if geom == "circle":
phi = Norm(CoefficientFunction((x, y, z))) - 1
phi = Norm(CoefficientFunction((x, y, z/c))) - 1
cube.Add(OrthoBrick(Pnt(-1.41, -1.41, -1.41), Pnt(1.41, 1.41, 1.41)))
cube.Add(OrthoBrick(Pnt(-2, -2, -2), Pnt(2, 2, 2)))
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
print("refining..")
for i in range(n_cut_ref):
lsetp1 = GridFunction(H1(mesh, order=2))
InterpolateToP1(phi, lsetp1)
......@@ -60,12 +63,30 @@ with TaskManager():
# Tangential projection
n = Normalize(grad(lset_approx))
print(n)
h = specialcf.mesh_size
eta = 100.0 / (h * h)
rho = 1.0 / h
def coef_grad(u):
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):
if type(u) in [ProxyFunction, GridFunction]:
return ngsolve.grad(u)
else:
return coef_grad(u)
def Pmats(n):
return Id(3) - OuterProduct(n, n)
Ps = Pmats(Normalize(grad(phi)))
# define solution and right-hand side
functions = {
"extvsol1": ((-y - z) * x + y * y + z * z),
......@@ -76,10 +97,11 @@ with TaskManager():
"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"]))
uSol = Ps*CoefficientFunction((functions["extvsol1"], functions["extvsol2"], functions["extvsol3"]))
rhs = CoefficientFunction((functions["rhs1"], functions["rhs2"], functions["rhs3"]))
elif geom=="decocube":
uSol = CoefficientFunction((-(z ** 2), y, x))
uSol = Ps*CoefficientFunction((-(z ** 2), y, x))
uSol=uSol.Compile()
u, v = VhG.TnT()
......@@ -88,56 +110,68 @@ with TaskManager():
# Measure on the bulk around the surface
dx = dx(definedonelements=ba_IF, deformation=deformation)
def eps(u):
return Ps * Sym(grad(u)) * Ps
def divG(u):
if u.dim == 3:
return Trace(grad(u) * Ps)
if u.dim == 9:
N = 3
divGi = [divG(u[i, :]) for i in range(N)]
return CF(tuple([divGi[i] for i in range(N)]))
print("rhs..")
begin_time = datetime.datetime.now()
rhstemp = -Ps*divG(eps(uSol))+uSol
print("compile rhs2..")
rhs2 = rhstemp.Compile()
print("compile time:")
print(datetime.datetime.now() - begin_time)
rhs = rhstemp
# pickle.dump(rhs, open("rhs_"+geom+".p","wb"))
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
print("assemble a..")
a.Assemble()
f = LinearForm(VhG)
f += rhs * v * ds
g = LinearForm(VhG)
g += rhs2*v*ds
print("assemble f..")
begin_time = datetime.datetime.now()
print("Assemble time without compile:")
f.Assemble()
print(datetime.datetime.now() - begin_time)
begin_time = datetime.datetime.now()
print("Assemble time with compile:")
f.Assemble()
print(datetime.datetime.now() - begin_time)
gfu.vec[:]=0.0
print("invert..")
gfu.vec.data = a.mat.Inverse(VhG.FreeDofs(), inverse="pardiso") * f.vec
print("l2error : ", sqrt(Integrate((gfu - uSol)**2 * ds, mesh=mesh)))
print("l2error : ", sqrt(Integrate((gfu - Ps*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")
print("write vtk file..")
vtk = VTKOutput(ma=mesh,
coefs=[lset_approx, deformation, gfu, uSol, uerr, rhs],
names=["P1-levelset", "deform", "uh", "uex", "uerr", "rhs"],
coefs=[lset_approx, gfu, uSol, uerr],
names=["P1-levelset", "uh", "uex", "uerr"],
filename="vectorlaplacepenalty3d", subdivision=2)
vtk.Do()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment