# FEM: Electric Field between two arbitrary defined shapes

I was wondering how to do the following: I would like to compute the electrostatic field between two shapes using the FEM method.

(*Define Boundaries*)
air = Rectangle[{-3, -3}, {3, 3}];
object1 = Disk[];
object2 = Rectangle[{2, 0}, {2.5, 2}];
Show[Graphics[{Blue, air}], Graphics[{Magenta, object1}],Graphics[{Green, object2}]] Calculation of the electric field at every point {x,y} in 2D space: $$r_i$$ is the vector of the point charge; $$r$$ is the vector to the point in 2D (or also 3D) space where we want to calculate the electric field.

I make a Mathematica function out of it (for the moment I omit the constant term):

eField[x_, y_] := q Sum[({x, y} - pts[[i]])/Norm[{x, y} - pts[[i]]]^3, {i, n}]


where pts[[i]] are the boundary points of the charged object and x and y are coordinates of the "air" object.

How I would proceed:

1. I calculate the electrostatic field of object 1 -> $$E_1$$

2. I calculate the electrostatic field of object 1 -> $$E_2$$

3. I use superposition to get the resultant electric field: $$E_{Total} = E_1 +E_2$$

Needs["NDSolveFEM"];

r1 = RegionDifference[air, object1];
r2 = RegionDifference[air, object2];
mesh1 = ToElementMesh[r1];
mesh2 = ToElementMesh[r2];
mesh1["Wireframe"]
mesh2["Wireframe"] I would really appreciate if someone could show me how to do it in Mathematica using finite elements (FEM).

EDIT: Basend on the excellent answer below, I would like to use the answer here Get Perimeter Region from Object to automate the finding of the region boundaries for the DirichletCondition:

Needs["NDSolveFEM"];
(*Define Boundaries*)
air = Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[];
object2 = Rectangle[{-2.5, -2.5}, {2.5, -2}];
reg12 = RegionUnion[object1, object2];
reg = RegionDifference[air, reg12]

mesh = ToElementMesh[reg, {{-5, 5}, {-5, 5}},
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]

eq = Laplacian[u[x, y], {x, y}]; V1 = 1; V2 = -2;
bc = {DirichletCondition[u[x, y] == V1,
RegionRegionProperty[RegionBoundary[object1], {x, y},
"FastDescription"][][]],
DirichletCondition[u[x, y] == V2,
RegionRegionProperty[RegionBoundary[object2], {x, y},
"FastDescription"][][]]};
U = NDSolveValue[{eq == 0, bc}, u, {x, y} \[Element] mesh];

ef = -Grad[U[x, y], {x, y}];

DensityPlot[U[x, y], {x, y} \[Element] reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotPoints -> 50,
PlotRange -> {{-4, 4}, {-4, 4}}]

StreamDensityPlot[Evaluate[ef], {x, y} \[Element] reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> {x, y}, StreamStyle -> LightGray, VectorPoints -> Fine,
PlotRange -> {{-3, 3}, {-2.5, 3}}]


EDIT 2: Just for beauty: Parallel Plate Capacitor Use this answer to make it work: FEM Simulation: Meshing two Arbitrary objects in an "air" mesh • Have you seen this or this – user21 Apr 3 at 8:59
• @james The solution depends on the electrical conductivity of the disk and the rectangle. There is one solution for a conductor, and another for a dielectric. – Alex Trounev Apr 3 at 10:33
• The comment of @AlexTrounev is very important. Are the disk and the rectangle conductor ? – andre314 Apr 3 at 10:41
• @AlexTrounev here on MSE, there is a recurrent problem of confusion between 1) dielectric versus conductor 2) potential imposed at boundary versus presence of charges. I give up, though if you want to explain the things, I could do an effort to participate. The problem is not related to Mma but some people think the Mathematica+FEM will solve all this magically. – andre314 Apr 3 at 10:54
• One point more : The formula above giving the electrical field E=..r/r^3 is wrong in the 2D case. In the 2D case it is E=..r/r^2. I have already seen this error on MSE too. – andre314 Apr 5 at 13:22

In the case of two metal objects, we can set the potential of each object as $$V_1, V_2$$. Then the code for a numerical solution in 2D is

Needs["NDSolveFEM"];
(*Define Boundaries*)air = Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[];
object2 = Rectangle[{2, 0}, {2.5, 2}]; reg12 =
RegionUnion[object1, object2];
reg = RegionDifference[air, reg12];
mesh = ToElementMesh[reg,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]
eq = Laplacian[u[x, y], {x, y}]; V1 = 1; V2 = -2;
bc = {DirichletCondition[u[x, y] == V1, x^2 + y^2 == 1],
DirichletCondition[
u[x, y] ==
V2, (x == 2 || x == 2.5 && 0 <= y <= 2) || (y == 0 ||
y == 2 && 2 <= x <= 2.5)]};
U = NDSolveValue[{eq == 0, bc}, u, {x, y} ∈ mesh];

ef = -Grad[U[x, y], {x, y}];


Visualisation of solution

{DensityPlot[U[x, y], {x, y} ∈ reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotPoints -> 50,
PlotRange -> {{-4, 4}, {-4, 4}}],
StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> {x, y}, StreamStyle -> LightGray,
VectorPoints -> Fine, PlotRange -> {{-1, 3}, {-1, 3}}]} Update 1. Next code is devoted to solve electrostatic problem for combination of dielectric and conducting objects (glass cylinder and metal strip). For dielectric we put electric charge $$q_1$$, and for metal we put potential $$V_2$$. Code:

Needs["NDSolveFEM"];
par = {eps1 -> 3.5, eps2 -> 1.0}; air =
Rectangle[{-5, -5}, {5, 5}];
object1 = Disk[]; q1 = 1; vol1 =
NIntegrate[1, {x, y} ∈ object1]; rho1 = q1/vol1;
object2 = Rectangle[{2, 0}, {2.5, 2}];
rho[x_, y_] := rho1 Boole[{x, y} ∈ object1];
eps[x_, y_] :=
eps2 + (eps1 - eps2) Boole[{x, y} ∈ object1]; reg =
RegionDifference[air, object2];
mesh = ToElementMesh[reg,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.001 (0.1 + 10 Norm[Mean[vertices]])]]
mesh["Wireframe"]
V2 = -2; eq =
Inactive[Div][
eps[x, y] Inactive[Grad][u[x, y], {x, y}], {x, y}] == -2 Pi rho[x,
y]; bc =
DirichletCondition[u[x, y] == V2, {x, y} ∈ object2];
U = NDSolveValue[{eq /. par, bc}, u, {x, y} ∈ mesh];

ef = -Grad[U[x, y], {x, y}];


Visualisation

{DensityPlot[U[x, y], {x, y} ∈ mesh,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> Automatic, PlotPoints -> 50,
PlotRange -> {{-4, 4}, {-4, 4}}],
StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
FrameLabel -> {x, y}, StreamStyle -> LightGray,
VectorPoints -> Fine, PlotRange -> {{-1, 3}, {-1, 3}}],
StreamDensityPlot[Evaluate[ef], {x, y} ∈ reg,
ColorFunction -> Hue, FrameLabel -> {x, y}, StreamStyle -> Blue,
PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, PlotLegends -> Automatic]}