Commit 3135323c authored by tilman's avatar tilman
Browse files

Merge branch 'main' of git.rwth-aachen.de:tilmandiego.aleman/reproduction

parents e89838b8 122427ff
for ORDER in 1 2 3
do
for C in 1 1.1 1.5
do
./killing.py --refinements 5 -c $C --order $ORDER --plot --saveplot
done
done
for ORDER in 1 2 3
do
for C in "1" "0.5" "0.1"
do
./killing.py --geometry arrow --refinements 5 -c $C --order $ORDER --plot --saveplot
done
done
\ No newline at end of file
......@@ -3,21 +3,27 @@ from xfem.lsetcurv import *
import matplotlib.pyplot as plt
from netgen.csg import CSGeometry, OrthoBrick, Pnt
import math
from functools import reduce
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))
print("lsetp1 interpolation..")
InterpolateToP1(phi, lsetp1)
print("lsetp1 refinement..")
RefineAtLevelSet(lsetp1)
mesh.Refine()
print("lsetmeshadap...")
lsetmeshadap = LevelSetMeshAdaptation(mesh, order=order, threshold=1000, discontinuous_qn=True, heapsize=hs)
print("calc deform...")
deformation = lsetmeshadap.CalcDeformation(phi)
mesh.SetDeformation(deformation)
print("set deform..")
if order>1:
mesh.SetDeformation(deformation)
lset_approx = lsetmeshadap.lset_p1
return lset_approx, deformation
......@@ -93,27 +99,140 @@ def plot_errors(l2errs, l2errskf, geom, saveplot, order=2, refs=2, c=1):
hvals = []
ticks = range(1,len(l2errs)+1)
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")
hvals.append((1 / 2 ** (order+1)) ** k)
hvals2 = []
for k in range(len(l2errs)):
hvals2.append((1 / 2 ** order) ** k)
plt.semilogy(ticks,l2errs, label=geom+" with KF removal")
# plt.semilogy(l2errskf, label=geom + " without KF removal")
plt.semilogy(ticks,hvals, label=r'$O(h^'+ str(order+1)+r')$')
plt.semilogy(ticks,hvals2, label=r'$O(h^'+ str(order)+r')$')
plt.ylabel(r'$L2$ error')
plt.xlabel("Refinement")
plt.title("Convergence removing KF")
plt.legend()
plt.xticks(range(len(l2errs)))
plt.xticks(ticks)
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 = []
def computeEV(phi, mesh, lset_approx, deformation, order=2):
with TaskManager():
print("starting ev comp")
vlams = []
Vh = VectorH1(mesh, order=order, dirichlet=[], low_order_space=True)
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
blocks = []
freedofs = VhG.FreeDofs()
class SymmetricGS(BaseMatrix):
def __init__ (self, smoother):
super(SymmetricGS, self).__init__()
self.smoother = smoother
def Mult (self, x, y):
y[:] = 0.0
self.smoother.Smooth(y, x)
self.smoother.SmoothBack(y,x)
def Height (self):
return self.smoother.height
def Width (self):
return self.smoother.height
vertexdofs = BitArray(Vh.ndof)
vertexdofs[:] = False
print("doing stuff for precond")
"""
for v in mesh.vertices:
# vdofs = set(d for d in VhG.GetDofNrs(el) if freedofs[d])
for d in Vh.GetDofNrs(v):
vertexdofs[d] = True
# blocks.append(vdofs)
blocks = [set(d for el in mesh[v].elements for d in VhG.GetDofNrs(el) if freedofs[d]) for v in mesh.vertices]
# blocks = [reduce(lambda x,y: x or y,[set(d for d in VhG.GetDofNrs(el) if freedofs[d])for el in mesh[v].elements]) for v in mesh.vertices]
def blockcr(FEspace):
for v in FEspace.mesh.vectices:
vdofs = set()
for el in FEspace.mesh[v].elements:
vdofs |= set(d for d in FEspace.GetDofNrs(el)
if freedofs[d])
# blocks = pool.map(partial(calcstuff, mesh=mesh), range(len(mesh.vertices)))
vertexdofs &= VhG.FreeDofs()
"""
print("done with stuff for precond")
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
pre1 = Preconditioner(a2, "direct", inverse="pardiso")
a2.Assemble()
coarsepre = a2.mat.Inverse(vertexdofs, inverse="pardiso")
# pre = Projector(VhG.FreeDofs(), True)
# pre1 = SymmetricGS(a2.mat.CreateSmoother(VhG.FreeDofs()))
# pre1 = SymmetricGS(a2.mat.CreateBlockSmoother(blocks, parallel=True))
pre = pre1#+coarsepre
# pre = a2.mat.Inverse(freedofs=VhG.FreeDofs(), inverse="pardiso")
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,maxit = 15)
def getevproblem(mesh, phi, lset_approx, deformation, order=2):
Vh = VectorH1(mesh, order=order, dirichlet=[])
ci = CutInfo(mesh, lset_approx)
ba_IF = ci.GetElementsOfType(IF)
VhG = Restrict(Vh, ba_IF)
......@@ -150,13 +269,18 @@ def computeEV(phi, mesh, lset_approx, deformation, order=2):
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()
pre = Preconditioner(a2, "local")
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
return a2, b, pre, Vh, VhG
def computeEV2(Vh, VhG, a2, b, pre):
Vh.Update()
VhG.Update()
print("updated..")
a2.Assemble()
b.Assemble()
return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True)
return solvers.PINVIT(a2.mat, b.mat, pre, 3, GramSchmidt=True,maxit = 20)
\ No newline at end of file
......@@ -40,7 +40,8 @@ def getKilling(mesh, lset_approx, deformation, vlam, evecs, step,alpha=3, maxsy
# epsilon =(cuthmax)**(alpha)
alpha=order+1
epsilon = h ** alpha
epsilon=0
if False:
epsilon=0
print("Epsilon:")
print(epsilon)
# epsilon=1
......@@ -262,18 +263,22 @@ if __name__ == "__main__":
for k in range(3):
vlams.append([])
for i in range(n_cut_ref + 1):
print("refinement is: ", i)
if i==0:
# a,b,pre, Vh, VhG = getevproblem(mesh, phi, lset_approx, deformation, order)
pass
if i != 0:
lset_approx, deformation = refineMesh(mesh, phi, order)
with TaskManager():
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
vlam, ev = computeEV(phi,mesh,lset_approx, deformation, order)
# vlam, ev = computeEV2(Vh, VhG, a, b, pre)
for k in range(3):
vlams[k].append(vlam[k])
if i < 2 and i != 0:
symmetries = 0
for k in range(3):
# if math.log(abs(1 - vlams[k][1]), 10) < math.log(abs(vlams[k][0] - vlams[k][1]), 10):
if math.log(abs(1 - vlams[k][1]), 2) <-(order+1)*(i+2):
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))
......@@ -286,12 +291,13 @@ if __name__ == "__main__":
temp -= (vlams[j][-1] - vlams[j][-2]) ** 2 / (
(vlams[j][-1] - vlams[j][-2]) - (vlams[j][-2] - vlams[j][-3]))
# if math.log(abs(1 - vlams[j][-1]), 10) < math.log(abs(temp - vlams[j][-1]), 10):
if math.log(abs(1 - temp), 2) < -(order + 1) * (i+2):
print("corresponding log: ")
print(math.log(abs(1 - temp), 2))
print("condition: ", -(order+1)*(i+1))
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):
symmetries += 1
if i!=0:
if True and i!=0:
print("starting solve..")
print("symmetries: ", symmetries)
with TaskManager():
......
using Gridap
using GridapEmbedded
using LinearAlgebra
using GridapDistributed
using GridapGmsh: gmsh
using GridapGmsh
# Define geometry
function main()
gmsh.initialize()
gmsh.clear()
box = gmsh.model.occ.addBox(
-2.,-2.,-2.,2.,2.,2.)
gmsh.model.occ.synchronize()
gmsh.model.mesh.setSize(gmsh.model.getEntities(0),0.1)
gmsh.model.mesh.generate(3)
gmsh.write("model.msh")
gmsh.finalize()
n = 10
# Background model
partition = (n,n,n)
pmin = 2*Point(-1,-1,-1)
pmax = 2*Point(1,1,1)
# bgmodel = simplexify(CartesianDiscreteModel(pmin,pmax,partition))
bgmodel =GmshDiscreteModel("model.msh")
# Select geometry
R = 0.5
ϕ(x) = norm(x)-1
# Forcing data
ud = 1
f = 10
println("doing cutter")
lscutter = GridapEmbedded.LevelSetCutters.sphere(1)
# Cut the background model
cutgeo = cut(bgmodel, lscutter)
# Setup integration meshes
Ω = Triangulation(cutgeo,PHYSICAL)
Γd = EmbeddedBoundary(cutgeo)
# Setup normal vectors
n_Γd = get_normal_vector(Γd)
#writevtk(Ω,"trian_O")
# writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(Γg,"trian_Gg",cellfields=["normal"=>n_Γg])
#writevtk(Triangulation(bgmodel),"bgtrian")
# Setup Lebesgue measures
order = 1
degree = 2*order
= Measure(Ω,degree)
dΓd = Measure(Γd,degree)
reffe = ReferenceFE(lagrangian, VectorValue{3, Float64}, order)
# Setup FESpace
Ω_act = Triangulation(cutgeo,ACTIVE)
normalspace = FESpace(Ω_act,reffe)
# ntilde = interpolate_everywhere(n_Γd, normalspace)
println("setup TnT")
V = TestFESpace(Ω_act,ReferenceFE(lagrangian,VectorValue{3,Float64},order),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)
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))))*
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
# ∫( 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)
println("solve")
ls = LinearFESolver()
uh = solve(op)
# Postprocess
# if outputfile !== nothing
writevtk(Ω,"vectorlaplace_julia",cellfields=["uh"=>uh])
end
mesh_partition=(4,4)
main()
\ No newline at end of file
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