# Electrostatic problem

**URL:** <https://fenicsproject.discourse.group/t/electrostatic-problem/4033>\
**Category:** mathematics\
**Created:** [August 31, 2020, 9:17pm UTC](https://fenicsproject.discourse.group/t/electrostatic-problem/4033 "2020-08-31T21:17:27Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![gabriel\_clemente](https://avatars.discourse-cdn.com/v4/letter/g/13edae/32.png) [@gabriel\_clemente](https://fenicsproject.discourse.group/u/gabriel_clemente)\
**Post date:** [August 31, 2020, 9:17pm UTC](https://fenicsproject.discourse.group/t/electrostatic-problem/4033/1 "2020-08-31T21:17:27Z")

</div>

Hi guys, I’m from Brazil and here no one know how to use FEniCS. I’m passing through a problem that probably you can help. I must simulate a little nanometric peace of a big metalic plane (on a potencial V) with some imperfections. Well, here is my idea: make a box that bottom surface is my irregular plane (V Dirichilet boundary condition), top surface is a simple plane (theoretically calculated V’ Dirichilet boundary condition) and side surfaces are simple planes (dV/dn = 0 vonn Neumann boundary condition). When I simulate on FEniCS I find values that no make sense (for example: my potencial energy muste be x _10^-10 but the simulation gives me y_10^-32). I think could be a scale problem. In my geometry “one unity” = 100nm. I have already changed my constants, but it wasn’t work.

Could someone help me?

Here is my idea:

 ![my_problem](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/52f797744e318eaa33a49477a86ce5658d08fe42.jpeg)

And here is my code:

```
from fenics import *
from dolfin import *
import meshio
import numpy

# MESH CONVERSION (.msh ---> .xdmf)

# read the 3D Gmsh file:
msh = meshio.read("fenics_test.msh")

# write the Gmsh file in xdmf format:
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}))

# write surface data:
meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}, cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}}))

# write volume data:
meshio.write("cf.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]}, cell_data={"tetra": {"name_to_read": msh.cell_data["tetra"]["gmsh:physical"]}}))

# create mesh:
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

# create boundary subdomain facet function:
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# create subdomain cell function:
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("cf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

# SIMULATION

# Function space:
S_1 = FunctionSpace(mesh, 'P', 1)
VS_1 = VectorFunctionSpace(mesh, 'P', 1)

# Boundary conditions:
bc_bottom = DirichletBC(S_1, Constant(1000.0), mf, 1) # 1: bottom boundary tag
bc_top = DirichletBC(S_1, Constant(904.5340), mf, 2) # 2: top boundary tag

bcs = [bc_bottom, bc_top]

# Variational form:
dx = Measure('dx', domain=mesh, subdomain_data=cf)
phi = TrialFunction(S_1)
v = TestFunction(S_1)
f = Constant(0.0)
a = inner(grad(phi), grad(v))*dx
L = inner(f, v)*dx

# Solution:
phi = Function(S_1) # electric potential [V]
solve(a == L, phi, bcs)
e_0 = 8.854*10**(-12) # permittivity of free space
pi = 3.1416
a_O2 = 1.562*10**(-30)*4*pi*e_0 # O2 atomic polarizability
a_N2 = 1.710*10**(-30)*4*pi*e_0 # N2 atomic polarizability
E = Function(VS_1)
E.assign(project((-(grad(phi))), VS_1)) # electric field [V/m]
U = Function(S_1)
U.assign(project((-(a_O2*inner(E, E))), S_1)) # potential energy [J]
F = Function(VS_1)
F.assign(project((-(grad(U))), VS_1)) # force field [N]

# Check points
print(phi(0, 0, 5.001))
print(phi(0, 0, 50))
print(E(0, 0, 5.001))
print(U(0, 0, 5.001))
print(F(0, 0, 5.001))

# Plot:
XDMFFile('potencial.xdmf').write(phi)
XDMFFile('electric_field.xdmf').write(E)
XDMFFile('energy.xdmf').write(U)
XDMFFile('force_field.xdmf').write(F)

```

I really don’t know if I’m doing something wrong on my code or another thing. I think my geometry is ok, the potencial looks good on my domain.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [August 31, 2020, 9:51pm UTC](https://fenicsproject.discourse.group/t/electrostatic-problem/4033/3 "2020-08-31T21:51:10Z")

</div>

> [@gabriel\_clemente](#):
>
> `E.assign(project((-(grad(phi))), VS_1)) `

You should not project the gradient of a CG-1 function into a CG-1 space, but a DG-0 space.

> [@gabriel\_clemente](#):
>
> `U.assign(project((-(a_O2*inner(E, E))), S_1))`

To avoid Machine precision errors, i would recommend not multiplying by a\_O2, and rather scale your expected result by this constant.