Electrostatic Examples for NGSolve Beginners

I'm an electrical engineer (power systems - 20 years). I have some familiarity with commercial FEM programs. I've been interested in learning NGSolve, but am finding my background too far removed from the language and concepts of the examples and documentation. Clearly the developers are brilliant folks and have worked very hard. Even with that, I've had a hard time getting traction. That said, I'm motivated to persist as it seems NGSolve has several advantages over other platforms I've looked at for my particular applications. Perhaps there are others in the same boat.

My goal with this thread is to get some basic 2D and 3D examples going in the area of electrostatics. Hopefully over time this can serve as supplementary documentation for folks like myself who are a bit overwhelmed by trial functions, test functions, bilinear forms, Galerkin methods etc.

I'll start with the simple capacitor example in the figure below. The top and bottom of the unit square are capacitor plates.



As a first attempt I tried to adapt the Poisson example (for those unfamiliar with these concepts, the Laplace equation is a special case of Poisson equation with no source (no space charge) in the dielectric region. The faulty output result follows.




Questions:

1. What is the syntax for defining permittivity for the dielectric region?
2. Boundary conditions didn't work - how can these be implemented correctly?
3. Seems like the linearform equation should be set to zero to implement the Laplace equation, but that didn't work when I tried it (not shown here) - any suggestions?
4. Once the space potential is found can a simple grad() function be used to find the e-field?
5. Suggestions for syntax to calculate the plate surface charge density?

Hi jleman,

1) As long as you just have one material the permittivity appears as a constant in the bilinear form If you have more materials, which are named correctly you can use the syntax
2) You have to set all boundary conditions at once. A second Set would set everything back to zero from the previous Set. For different values on different boundaries you can use
3) Yes, for the Laplace equation it has to be zero. I guess you got the zero solution, when setting it to zero. The reason for this is that you have to homogenize your problem for non-homogeneous (non zero) Dirichlet data as described in the documentation , namely
4) Yes works

5) How is the plate surface charge density defined (I'm not an electrical engineer)? Is it coming from the Gauss's law described in your equations? If yes should work (neglecting the constants).

Attached you'll find the whole code.

Best,
Michael

Only a very small comment on Michael's answer:
You can visualize the E field either by clicking Visual->Draw Surfacevectors and increasing the grid on the top of the same menu or using python:
Best
Christopher

Ok, making progress. I have various questions, but will start with these:

1. Christopher's addition shows an arrow field, but the draw function for E does not show any color plot behind the arrow field. I assume this is due to E being a vector field? To see the magnitude of the electric field I tried the following (from the Navier Stokes example where velocity field magnitude is plotted). How would I correct this?:
   
   Code:

   Emag = Norm(-grad(gfu.components[0])) Draw(Emag, mesh, "Emag")


2. Per your suggestion on another thread I used the following to attempt to print numerical results for gfu (the electric potential) at all vertices. The answer's don't match up with the results on the color plot which are correct.
   Code:

   for vert in mesh.vertices: #print(vert.point) #this works when not commented out print(gfu(mesh(vert.point))) #this does not show correct results


3. I was looking in in section 1.3 and 1.8 for way to iterate over elements associated with a boundary. I think you can see the intent of this code snippet which gives an error 'ngsolve.comp.Region' object is not iterable. Is there a direct way to iterate over elements previously associated with a boundary condition?
   Code:

   for vert in mesh.Boundaries("top"): print(vert.point)

To 1)
You can plot the magnitude of the vectorfield from python using or in the Netgen GUI by clicking Visual -> Scalar Function and selecting "func (E)"

Best wishes,
Henry

For 2:
A ugly bug/feature:
the mesh expects x,y,z values as input.
If no z is given z = 0, if no y -> y =0
it allows x,y,z to be array like to evalute many points more efficiently.
vert.points is array like, so it assumes this is an array for x values and you get 2 results namely gfu evaluated at (x, 0), (y,0) thats not what you want
works.

For 3 there are 2 options.
You can either evalutate the function anywhere on the boundary:
or if you need it at the vertices collect the boundary vertices firts:
Note that el.vertices just gives a NodeId which is basically the node type (vertex) plus node number. You have to put that into the mesh again to create a node with all information (like point,...)

Best
Christopher