Commit 85b9f8c4 authored by Tilman Alemán's avatar Tilman Alemán
Browse files

getting closer to killing killing

parent 8143b898
......@@ -114,7 +114,7 @@ with TaskManager():
# print("Epsilon:")
# print(epsilon)
# epsilon=1
eta = 1000.0 / (h)**2
eta = 1.0 / (h)**2
rhou = 1.0 / h
rhop = h
......@@ -141,49 +141,60 @@ with TaskManager():
def E(u):
return (Pmat * Sym(grad(u)) * Pmat)
return (Pmat * Sym(grad(u)) * Pmat)- (u * n) * weing
rhob = h
u2, v2 = VhG.TnT()
a2 = BilinearForm(VhG, symmetric=True)
a2 += InnerProduct(E(u2),E(v2)) * ds
# a2 +=(Pmat * u2) * (Pmat * v2) * ds
a2 +=(Pmat * u2) * (Pmat * v2) * ds
a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
a2 += (eta * ((u2 * n) * (v2 * n))) * ds
a2 += (eta * ((u2 * ntilde) * (v2 * ntilde))) * ds
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(u2,v2)*dX
b += InnerProduct(Pmat*u2,Pmat*v2)*ds
b += rhob*(grad(u2)*n)*(grad(v2)*n)*dX
b += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 0)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 1)
egfu = GridFunction(VhG)
egfu.vec.data = evecs.vecs[0].data
print("is it eigen?")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Pmat*Sym(ngsolve.grad(egfu))*Pmat)**2,mesh=mesh))
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Ps*Sym(ngsolve.grad(egfu))*Ps)**2,mesh=mesh))
fun = GridFunction(VhG)
kf = GridFunction(VhG)
l2errs = []
print(len(vlam))
temp = GridFunction(VhG)
temparr = []
for k in range(3):
temparr.append(GridFunction(VhG))
temparr[k].Set(kvfs[k], definedonelements=ba_IF)
for k in range(3):
lam = vlam[k]
print(k)
print("lam: "+ str(abs(lam)))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
evecnorm = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5))
kvfnorm = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kvfs[k],kvfs[k]), mesh=mesh, order=5))
scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k)/evecnorm,kvfs[k]/kvfnorm), mesh=mesh, order=5)
sgn = np.sign(scal)
l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(sgn*kvfs[k]/kvfnorm-evecs.MDComponent(k)/evecnorm)**2, mesh=mesh, order=5)/evecnorm))
temp.vec.data = evecs.vecs[k].data
for j in range(maxk):
if True:
kvfnorm = sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(kvfs[j],kvfs[j]), mesh=mesh, order=5))
scal = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),kvfs[j]), mesh=mesh, order=5)
temp.vec.data += -scal*temparr[k].vec
l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(temp)**2, mesh=mesh, order=5)))
return l2errs, evecs.MDComponent(0), evecs.MDComponent(1), evecs.MDComponent(2)
......
......@@ -184,10 +184,10 @@ with TaskManager():
print("Epsilon:")
print(epsilon)
# epsilon=1
eta = 1000.0 / (h)**2
eta = 1 / (h)**2
rhou = 1.0 / h
rhop = h
rhob = h
# Measure on surface
ds = dCut(lset_approx, IF, definedonelements=ba_IF, deformation=deformation)
# Measure on the bulk around the surface
......@@ -252,20 +252,23 @@ with TaskManager():
u2, v2 = VhG.TnT()
a2 = BilinearForm(VhG, symmetric=True)
a2 += InnerProduct(E(u2),E(v2)).Compile(realcompile=comp) * ds
# a2 +=epsilon* (Pmat * u2) * (Pmat * v2) * ds
# a2 += rhou * ((grad(u2) * n) * (grad(v2) * n)) * dX
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
a2.Assemble()
b = BilinearForm(VhG, symmetric=True)
b += InnerProduct(u2,v2)*dX
b += InnerProduct(Pmat*u2,Pmat*v2)*ds
b += rhob*(grad(u2)*n)*(grad(v2)*n)*dX
b += ((u2 * n) * (v2 * n)) * ds
b.Assemble()
evecs = GridFunction(VhG, multidim=3)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 0)
vlam = ArnoldiSolver(a2.mat,b.mat,VhG.FreeDofs(), evecs.vecs, 1)
# vlam, evecs = solvers.PINVIT(a2.mat, b.mat,pre=IdentityMatrix(), num=3)
"""
rows,cols,vals = a.mat.COO()
......@@ -304,7 +307,7 @@ with TaskManager():
print("is it eigen?")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(Pmat*Sym(ngsolve.grad(egfu))*Pmat)**2,mesh=mesh))
# print("should be smaller then: "+str(cuthmax**alpha-2*cuthmax**2))
if False and k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
if True and k<maxk:#abs(lam)<abs(cuthmax**alpha-2*cuthmax**2):
print("removing number: "+str(k))
evecnorm = Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=InnerProduct(evecs.MDComponent(k),evecs.MDComponent(k)), mesh=mesh, order=5)
......
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