Laplace equation and back

General Discussion
1

Hello everyone,

I am trying to implement Laplace equation for a capacitor on which the initial voltage is applied to only a small part of the electrodes. To see if everything is right, I am doing the following: first, imposing a voltage somewhere on the electrodes, then computing the electric field and charge density, and finally going back to the potential. Comparing the potential set at the beginning with the final potential computed, I get two different things.

The equation is implemented with mixed boundary conditions: Dirichlet on the sides and Neumann on the top and bottom. At first I thought that this effect was due to the two different domains used (only the plates for the charges and potential and field on the plates and inside), but having all of the three quantities on the same domain does not seem to improve the situation. Here is the code I have by now, it uses the Cube.idp file provided in the documentation:

include "Cube.idp"

// Problem spatial data
int[int] Nxyz = [40, 10, 4];

real Dlength = 32.5e-2;
real Dheight = 20e-6;
real Dspace = 200e-6;

real[int, int] Bxyzhb = [[0., Dlength], [0., 7.5e-2], [0., Dheight]];
real[int, int] Bxyzi = [[0., Dlength], [0., 7.5e-2], [0., Dspace]];

int[int, int] Lxyzh = [[11, 12], [13, 14], [15, 16]];
int[int, int] Lxyzb = [[21, 22], [23, 24], [25, 26]];
int[int, int] Lxyzi = [[31, 32], [33, 34], [35, 36]];

// Problem parameters
real eps0 = 8.85e-12;
real epsr = 2.2;

// Create mesh
mesh3 Thh = Cube(Nxyz, Bxyzhb, Lxyzh);
Thh = movemesh(Thh, [x, y, z + (Dheight + Dspace) / 2]);
mesh3 Thb = Cube(Nxyz, Bxyzhb, Lxyzb);
Thb = movemesh(Thb, [x, y, z - (Dheight + Dspace) / 2]);
mesh3 Th = Thh + Thb;

mesh3 Thi = Cube(Nxyz, Bxyzi, Lxyzi);
mesh3 Thf = Th + Thi;

// FE basis
fespace Vhv(Thf, [P1, P1, P1]);
Vhv [Ex, Ey, Ez];

fespace Vhs(Thf, P1);
Vhs rho;
Vhs u;

fespace Vhf(Thf, P1);
Vhf phi;
Vhf v;

// Poisson (or Laplace) var. form.
problem InversePoisson(phi, v)
	= int3d(Thf)(
		(
			  dx(phi) * dx(v)
			+ dy(phi) * dy(v)
			+ dz(phi) * dz(v)
		)
	  )
	+ int3d(Thf)(
		rho * v / (eps0 * epsr)
	  )
	+ int2d(Thf, 15, 26)(
		(dx(phi) * N.x + dy(phi) * N.y + dz(phi) * N.z) * v
	  )
	+ on(11, 12, 13, 14, 21, 22, 23, 24, 31, 32, 33, 34, phi=0)
	;

func real phi0() {
	if(z >= Dspace / 2.) {
		if(x > 1e-2 && x < 2e-2) {
			return 10.;
		}
	}
}

phi = phi0();

plot(phi, wait=true);

// Step 1: Direct Poisson to get density
[Ex, Ey, Ez] = [-dx(phi), -dy(phi), -dz(phi)];
rho = eps0 * epsr * (dx(Ex) + dy(Ey) + dz(Ez));

// Step 2: Inverse Poisson to get potential back
InversePoisson;

plot(phi);

Thanks if you have a clue on how I could make it work!