Draw error in 3D
- melanie
- Topic Author
- New Member
Less
More
1 year 4 months ago - 1 year 4 months ago #4828
by melanie
Draw error in 3D was created by melanie
Hi,
I would like to draw the error of my eigenfunction in 3D (as I've already done in 2D), but when I try to do that with "Draw(u-exact,mesh)" I get an error (see attachments).
What should I do in order to get the picture of the error also in 3D?
My code:
import sys
import numpy
from ngsolve import *
from ngsolve.webgui import Draw
import netgen.geom2d as geom2d
from netgen.geom2d import unit_square
import math
import scipy.linalg
from numpy import random
import netgen.gui
from ngsolve import Draw, Redraw, Mesh # just for visualization
from netgen.csg import *
left = Plane(Pnt(0,0,0), Vec(-1,0,0))
right = Plane(Pnt(1,1,1), Vec( 1,0,0))
front = Plane(Pnt(0,0,0), Vec(0,-1,0))
back = Plane(Pnt(1,1,1), Vec(0, 1,0))
bot = Plane(Pnt(0,0,0), Vec(0,0,-1))
top = Plane(Pnt(1,1,1), Vec(0,0, 1))
cube = left * right * front * back * bot * top
geo = CSGeometry()
geo.Add (cube)
exact = (-4/math.sqrt(2)) * sin(math.pi * x) * sin(math.pi * y) * sin(math.pi * z)
for meshsize in range(1,4):
maxh = 2 ** (-meshsize/2)
mesh = Mesh(geo.GenerateMesh(maxh = maxh))
print("Number of Tetrahedra", mesh.ne)
#mesh = Mesh(geo.unit_cube(maxh = 2**(-4)))
fes = FESpace("nonconforming",mesh, dirichlet=".*")
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
pre = Preconditioner(a, "bddc")
m = BilinearForm(fes)
m += u*v*dx
a.Assemble()
m.Assemble()
u = GridFunction(fes)
r = u.vec.CreateVector()
w = u.vec.CreateVector()
Mu = u.vec.CreateVector()
Au = u.vec.CreateVector()
Mw = u.vec.CreateVector()
Aw = u.vec.CreateVector()
r.FV().NumPy()[:] = random.rand(fes.ndof)
u.vec.data = Projector(fes.FreeDofs(), True) * r
for i in range(20):
Au.data = a.mat * u.vec
Mu.data = m.mat * u.vec
auu = InnerProduct(Au, u.vec)
muu = InnerProduct(Mu, u.vec)
# Rayleigh quotient
lam = auu/muu
#print (lam / (math.pi**2))
# residual
r.data = Au - lam * Mu
w.data = pre.mat * r.data
w.data = 1/Norm(w) * w
Aw.data = a.mat * w
Mw.data = m.mat * w
# setup and solve 2x2 small eigenvalue problem
asmall = Matrix(2,2)
asmall[0,0] = auu
asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)
asmall[1,1] = InnerProduct(Aw, w)
msmall = Matrix(2,2)
msmall[0,0] = muu
msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)
msmall[1,1] = InnerProduct(Mw, w)
# print ("asmall =", asmall, ", msmall = ", msmall)
eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)
# print (eval, evec)
u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w
z = GridFunction (fes, name="z")
z.Update()
solver = m.mat.Inverse(fes.FreeDofs(), inverse = 'sparsecholesky')
normB_r = sqrt(InnerProduct(solver * r, r))
#print('norm Br', normB_r)
max=0
for face in mesh.facets:
point1 = numpy.array(mesh.__getitem__(face.vertices[0]).point)
point2 = numpy.array(mesh.__getitem__(face.vertices[1]).point)
#print(mesh.__getitem__(face.vertices[0]).point,mesh.__getitem__(face.vertices[1]).point)
dist = numpy.linalg.norm(point1 - point2)
if dist > max:
max = dist
#Draw (u)
Draw(u-exact, mesh)
Bests,
Melanie
I would like to draw the error of my eigenfunction in 3D (as I've already done in 2D), but when I try to do that with "Draw(u-exact,mesh)" I get an error (see attachments).
What should I do in order to get the picture of the error also in 3D?
My code:
import sys
import numpy
from ngsolve import *
from ngsolve.webgui import Draw
import netgen.geom2d as geom2d
from netgen.geom2d import unit_square
import math
import scipy.linalg
from numpy import random
import netgen.gui
from ngsolve import Draw, Redraw, Mesh # just for visualization
from netgen.csg import *
left = Plane(Pnt(0,0,0), Vec(-1,0,0))
right = Plane(Pnt(1,1,1), Vec( 1,0,0))
front = Plane(Pnt(0,0,0), Vec(0,-1,0))
back = Plane(Pnt(1,1,1), Vec(0, 1,0))
bot = Plane(Pnt(0,0,0), Vec(0,0,-1))
top = Plane(Pnt(1,1,1), Vec(0,0, 1))
cube = left * right * front * back * bot * top
geo = CSGeometry()
geo.Add (cube)
exact = (-4/math.sqrt(2)) * sin(math.pi * x) * sin(math.pi * y) * sin(math.pi * z)
for meshsize in range(1,4):
maxh = 2 ** (-meshsize/2)
mesh = Mesh(geo.GenerateMesh(maxh = maxh))
print("Number of Tetrahedra", mesh.ne)
#mesh = Mesh(geo.unit_cube(maxh = 2**(-4)))
fes = FESpace("nonconforming",mesh, dirichlet=".*")
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
pre = Preconditioner(a, "bddc")
m = BilinearForm(fes)
m += u*v*dx
a.Assemble()
m.Assemble()
u = GridFunction(fes)
r = u.vec.CreateVector()
w = u.vec.CreateVector()
Mu = u.vec.CreateVector()
Au = u.vec.CreateVector()
Mw = u.vec.CreateVector()
Aw = u.vec.CreateVector()
r.FV().NumPy()[:] = random.rand(fes.ndof)
u.vec.data = Projector(fes.FreeDofs(), True) * r
for i in range(20):
Au.data = a.mat * u.vec
Mu.data = m.mat * u.vec
auu = InnerProduct(Au, u.vec)
muu = InnerProduct(Mu, u.vec)
# Rayleigh quotient
lam = auu/muu
#print (lam / (math.pi**2))
# residual
r.data = Au - lam * Mu
w.data = pre.mat * r.data
w.data = 1/Norm(w) * w
Aw.data = a.mat * w
Mw.data = m.mat * w
# setup and solve 2x2 small eigenvalue problem
asmall = Matrix(2,2)
asmall[0,0] = auu
asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)
asmall[1,1] = InnerProduct(Aw, w)
msmall = Matrix(2,2)
msmall[0,0] = muu
msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)
msmall[1,1] = InnerProduct(Mw, w)
# print ("asmall =", asmall, ", msmall = ", msmall)
eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)
# print (eval, evec)
u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w
z = GridFunction (fes, name="z")
z.Update()
solver = m.mat.Inverse(fes.FreeDofs(), inverse = 'sparsecholesky')
normB_r = sqrt(InnerProduct(solver * r, r))
#print('norm Br', normB_r)
max=0
for face in mesh.facets:
point1 = numpy.array(mesh.__getitem__(face.vertices[0]).point)
point2 = numpy.array(mesh.__getitem__(face.vertices[1]).point)
#print(mesh.__getitem__(face.vertices[0]).point,mesh.__getitem__(face.vertices[1]).point)
dist = numpy.linalg.norm(point1 - point2)
if dist > max:
max = dist
#Draw (u)
Draw(u-exact, mesh)
Bests,
Melanie
Attachments:
Last edit: 1 year 4 months ago by melanie.
Time to create page: 0.116 seconds