I am trying to compute charge from Electric field by Integrating E over a boundary (In my problem I get E on the boundary facets as a result of the computation).
However the results show 1) Lot of Oscillation 2) Are scaled by some factor : This seems to depend on the mesh for some reason. In my actual problem even the surface charges are not computed correctly (See Paraview image)
#MWE
from mpi4py import MPI
from dolfinx import mesh
import ufl
import numpy as np
from dolfinx import fem,plot
import pyvista
from dolfinx.fem.petsc import LinearProblem
def computecharge(M,E):
V=M.mesh
dx=ufl.dx(V)
ds=ufl.ds(V)
n=ufl.FacetNormal(V)
v = ufl.TrialFunction(M)
w = ufl.TestFunction(M)
a=ufl.inner(v,w)*dx
L=ufl.inner(ufl.inner(E,n),w)*ds
problem = LinearProblem(a,L,[],petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
charge=problem.solve()
return charge
domain = mesh.create_unit_cube(MPI.COMM_WORLD, 8, 8,8, mesh.CellType.tetrahedron)
tdim=domain.geometry.dim
V1 = fem.functionspace(domain, ("Lagrange", 2))
#Higher order to avoid gibbs oscillations
E = fem.Constant(domain,(0.0,0.0,-1.0))
q=computecharge(V1,E)
pl = pyvista.Plotter()
topology, cell_types,x = plot.vtk_mesh(V1)
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid["charge"] = q.x.array
grid.set_active_scalars("charge")
pl.add_mesh(grid,show_edges=True,cmap='coolwarm')
pl.show()
Update: The Oscillations seem to reduce when projecting on to DG-1 space instead of Lagrange. Checking this in main problem now.