Draw error in 3D

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
 
Last edit: 1 year 4 months ago by melanie.
More
1 year 4 months ago #4829 by joachim
Replied by joachim on topic Draw error in 3D
If you use the Draw from the Netgen-gui you have to give also a name for the function:

Draw (u-uex, mesh, "error")

Joachim
 
Time to create page: 0.116 seconds