Commit a04dcd64 authored by tilman's avatar tilman
Browse files

scipy weirdness

parent 85b9f8c4
......@@ -101,7 +101,8 @@ with TaskManager():
VhG = Restrict(Vh, ba_IF)
n = Normalize(grad(lset_approx))
ntilde = Interpolate(Normalize(grad(phi)), VhG)
ntilde = GridFunction(VhG)
ntilde.Set(Normalize(grad(phi)), definedonelements=ba_IF)
h = specialcf.mesh_size
elvol = Integrate(levelset_domain={"levelset":phi, "domain_type":IF},cf=CoefficientFunction(1),mesh=mesh,element_wise=True)
......@@ -127,7 +128,7 @@ with TaskManager():
ninter = GridFunction(VhG)
ninter.Set(n, definedonelements=ba_IF)
ninter.Set(Normalize(grad(phi)), definedonelements=ba_IF)
weing = ngsolve.grad(ninter)
# rhs1ex = rhsfromsol(uSol, pSol, interpol=False)
......@@ -164,12 +165,26 @@ with TaskManager():
evecs = GridFunction(VhG, multidim=3)
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=(Ps*Sym(ngsolve.grad(egfu))*Ps)**2,mesh=mesh))
fun = GridFunction(VhG)
rows, cols, vals = a2.mat.COO()
A = sp.csr_matrix((vals, (rows, cols)))
rows2, cols2, vals2 = b.mat.COO()
B = sp.csr_matrix((vals2, (rows2, cols2)))
spevs ,spevecs = eigs(A,3,M=B,sigma=1,which='LM')
print("sp evals: ", spevs)
funs = []
print(len(spevecs[:,0]))
for k in range(3):
fun = GridFunction(VhG)
evc = spevecs[:,k]
fun.vec.FV().NumPy()[:] = evc
funs.append(fun)
print(len(fun.vec.FV().NumPy()[:]))
print("what does scipy say:")
print(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(evecs.MDComponent(0)-fun)**2, mesh=mesh))
kf = GridFunction(VhG)
# kf.vec.data = spevecs
l2errs = []
print(len(vlam))
temp = GridFunction(VhG)
......@@ -183,18 +198,25 @@ with TaskManager():
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):
temp.vec.data = evecs.vecs[k].data
# temp.Set( Pmat*evecs.MDComponent(k), definedonelements=ba_IF)
# temp.vec.data = evecs.vecs[k].data
temp.vec.data = funs[k].vec.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)
scal = Integrate(levelset_domain = {"levelset": phi, "domain_type": IF},
cf = InnerProduct(funs[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)))
# l2errs.append(sqrt(Integrate(levelset_domain={"levelset":phi, "domain_type":IF}, cf=(temp)**2, mesh=mesh, order=5)))
l2errs.append(
sqrt(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF}, cf = (temp) ** 2, mesh = mesh, order = 5)))
# print("scipy ev problem solution")
# print(sqrt(Integrate(levelset_domain = {"levelset": phi, "domain_type": IF}, cf = (temp) ** 2, mesh = mesh)))
return l2errs, evecs.MDComponent(0), evecs.MDComponent(1), evecs.MDComponent(2)
......@@ -219,7 +241,9 @@ with TaskManager():
print("starting solve..")
l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step=i)
# l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step=i)
l2err, ev1, ev2, ev3 = doKilling(mesh, lset_approx, deformation, step = i)
l2errors.append(l2err)
print("L2 errors:")
......
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