Code:
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()
Thank you for any help.
. Where the problem?