VTK file for visualisation in paraview

from netgen import gui from ngsolve import * from ngsolve import atan2 from ngsolve import sin from ngsolve import cos from netgen.geom2d import SplineGeometry geo = SplineGeometry() geo.AddRectangle( (-3, -3), (3, 3), bc = 'wall') geo.AddCircle ( (0, 0), r=0.75, leftdomain=0, rightdomain=1, bc='cyl') mesh = Mesh( geo.GenerateMesh(maxh=0.2)) mesh.Curve(3); Draw(mesh) V = H1(mesh, order=2, dirichlet='wall|cyl') gfu = GridFunction(V) u = V.TrialFunction() w = V.TestFunction() f = LinearForm(V) f += 0 * w * dx R=0.75 uinf = 1 exact = CoefficientFunction(uinf*(y-(R**(2)*y)/(x**(2)+y**(2)))) gexact = CoefficientFunction( ( (uinf*((2*R**(2)*x*y)/(x**(2)+y**(2))**2)), (uinf*(1-(R**(2)*(x**(2)-y**(2))/(x**(2)+y**(2))**2))) ) ) gfu = GridFunction(V) gfu.Set(exact, definedon=mesh.Boundaries('wall|cyl')) a = BilinearForm(V, symmetric=True) a += grad(u)*grad(w)*dx #c = Preconditioner(a, "local") c = Preconditioner(a, "direct") a.Assemble() f.Assemble() # the solution field BVP=BVP(bf=a,lf=f,gf=gfu, pre=c).Do() ndof = V.ndof print(ndof) Draw (gfu) den1=sqrt(Integrate(InnerProduct(exact,exact), mesh)) print ("L2-error:", sqrt(Integrate((gfu-exact)**2, mesh))/den1) den2=sqrt(Integrate(InnerProduct(exact,exact)+InnerProduct(gexact,gexact), mesh)) print ("H1-error:", sqrt(Integrate((gfu-exact)**2+InnerProduct(grad(gfu)-gexact,grad(gfu)-gexact), mesh))/den2) vtk = VTKOutput(ma=mesh,coefs=[gfu], names=["sol"], filename="vtk_example3", subdivision=3) vtk.Do() vtk = VTKOutput(ma=mesh,coefs=[gfu,grad(gfu)],names=["sol","gradsol"],filename="vtk_example4",subdivision=3) vtk.Do()