Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Draw error in 3D

More
11 months 3 weeks ago - 11 months 3 weeks 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: 11 months 3 weeks ago by melanie.
More
11 months 3 weeks 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.144 seconds