Hi, I want to integrate cos(6\theta) on an inner surface which is a quarter of a cylinder, that the angle is from 0 to \pi/2. The length of the cylinder is 6 and its radius is 0.665929. The integration result should be zero, but it gives 0.0188. Part of my code is:

V = FunctionSpace(domain, ("CG", 4))
u = Function(V)
u.interpolate(lambda x: (32*(x[0]/a)**6-48*(x[0]/a)**4+18*(x[0]/a)**2-1))
area_form_A = form(u*ds(35))
A = assemble_scalar(area_form_A)
print(A)

Here cos(6\theta) = 32cos(\theta)^6-48cos(\theta)^4+18cos(\theta)^2-1 and a=0.665929. What’s the problem?

Assuming you’re integrating with exact quadrature (to machine precision), then the remaining error is likely coming from imprecise resolution of the curve. You could check that your integral converges to zero at a rate of \mathcal{O}(h^{p+1}) where p is the order of the polynomials used on each cell to approximate the boundary. See for example, gmsh’s set_order function for meshing higher order boundaries.

Dear nate, thanks for your help. Now I have to set the order of the polynomials to 2 to make the surface integration accurately. Is it also necessary to set the order to 2 if I don’t care about the surface integration for example for a Laplace problem? Another question it’s the derivative to the potential solution of a static electric problem is not quite smooth. Do you know how this could be resolved?