# Problem calculating charges from Electric field

**URL:** <https://fenicsproject.discourse.group/t/problem-calculating-charges-from-electric-field/15786>\
**Category:** General\
**Created:** [September 12, 2024, 5:22am UTC](https://fenicsproject.discourse.group/t/problem-calculating-charges-from-electric-field/15786 "2024-09-12T05:22:02Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![bmitra265](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bmitra265/32/7430_2.png) [@bmitra265](https://fenicsproject.discourse.group/u/bmitra265)\
**Post date:** [September 12, 2024, 5:22am UTC](https://fenicsproject.discourse.group/t/problem-calculating-charges-from-electric-field/15786/1 "2024-09-12T05:22:02Z")

</div>

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)

 ![Pyvista Output](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4ead3e3d06e7bfa0e851b70de469afe5d993ea30.png)  
 ![Screenshot 2024-09-12 at 9.24.03 AM](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/4e4fbf9434f453f681e9e6298d9c60c01b6f108d.jpeg)

```auto
#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.