Solving Laplace's Equation in 2D with Piecewise Permitivity

[fyrefly00]

Hi everyone,
I'm working on an electrostatics solver that calculates various properties of a given 2D microstrip geometry. My code is based heavily off of Example 1, with some modifications inspired by the Volta miniapp. I'm solving Laplace's equation on my mesh with relevant boundary conditions for conductors, ground strips, air boundaries, etc. However, I am struggling to get the correct behavior for geometries where there are two or more dielectric regions (consider for example, a microstrip sitting on a substrate with some relative permittivity and a layer of air above it, or a stripline sandwiched between two different substrates). I've created a piece wise coefficient:

Vector permitivityVector(permitivities, mesh.bdr_attributes.Max());
 PWConstCoefficient *epsCoef = new PWConstCoefficient(permitivityVector);

where permitivities is a marker array that specifically highlights the plane surfaces.
However, I've been unable to get this permittivity coeffecient to have the desired effect on the solution. I'm assuming I need to use some form of Domain Integrator for the linear form or the bilinear form, but none of my attempts with BoundaryLFIntegrator, DomainLFIntegrator, DiffusionIntegrator, MassIntegrator, etc. have yielded fruitful results.

Is there a particular way to go about approaching this problem? How do I ensure that the solution takes the permitivities of each surface into account appropriately?

Cheers!

[mlstowell]

Hi, @fyrefly00,

Permittivities would normally be assigned within each region rather than on boundaries. Consequently your permittivityVector would have a length of mesh.attributes.Max() rather than mesh.bdr_attributes.Max(). Unless I'm misunderstanding your intentions.

Best wishes,
Mark

[mlstowell]

Hi, @fyrefly00,

The choice of integrators is mainly determined by the PDE being solved as well as the particular finite element method being employed. In this case I assume you're solving the standard Poisson equation from electrostatics: $-\nabla\cdot(\epsilon\nabla\phi)=\rho$

If you're following ex1p.cpp and Volta then you're probably using the H1 basis functions as test functions. In this case the left hand side of Poisson's equation will be represented by the DiffusionIntegrator with the epsCoef exactly as you've shown above. The right hand side just needs to integrate the products of $\rho$ with the H1 basis functions. This can be accomplished with a LinearForm using the DomainLFIntegrator with a Coefficient object representing $\rho$.

As I'm sure you know, the potential $\phi$ is not unique unless you specify the value of $\phi$ in at least one location. In the discrete setting this means that you need to specify Dirichlet (a.k.a essential) boundary conditions on at least one vertex. In ex1 this is done by specifying $\phi$ on the entire boundary but any single boundary surface would be sufficient. It is also possible to choose a value for the finite element solution at just one of its degrees of freedom but I'm not sure we have a good example of how this would be done.

Any boundary surfaces which are not set using Dirichlet BCs would use the natural boundary condition associated with the DiffusionIntegrator which is $\hat{n}\cdot\nabla\phi=0$. This corresponds to an electrical insulator. It is possible to override this default using a Neumann BC by adding a BoundaryLFIntegrator to the LinearForm representing the right hand side. This is equivalent to specifying the charge density along the boundary.

Best wishes,
Mark

[fyrefly00]

Mark,
Thank you-- this is all incredibly helpful information. I am indeed solving the Poissson equation, although in the form $\nabla\phi^2 =0$. I'm ultimately using the resultant $\phi$ from the solution to calculate the total charge on a conductor. As you said, I'm using Dirichlet BCs to enforce a value of $\phi$ on the surface of my conductor element.
In light of the fact that I don't know what the value of $\rho$ is (since this is ultimately what I'm solving for), how do I go about applying a DomainLFIntegrator on the left hand size with a Coeffecient that represents $\rho$?
Thank you again for your assistance; it is tremendously helpful and much appreciated.

[fyrefly00]

Touching base on this topic-- I've yet to figure out a reliable way of setting up the right hand side of the equation to work with a DomainLFIntegrator. I tried creating a GridFunctionCoeffecient based on the solution vector, but this didn't seem to have the desired effect.

[mlstowell]

Hi, @fyrefly00,

In many cases the LinearForm on right hand side is just zero. There are two types of fields that can appear in the right hand side both are related to known charge densities. The first is a volume charge density, for example a region of material with a known distribution of charge like a sphere with a constant charge distribution. In this case you add a domain integrator using something like the following:

LinearForm b(&fespace);
b.AddDomainIntegrator(new DomainLFIntegrator(*rhoCoef));

The second case is probably even more unusual and it corresponds to a known surface charge density on some portion of the boundary of the domain. If your known surface charge is given by sigmaCoef this would be added like this:

LinearForm b(&fespace);
b.AddBoundaryIntegrator(new BoundaryLFIntegrator(*sigmaCoef), bdr_attr_marker);

Where bdr_attr_marker is an array of length mesh.bdr_attributes.Max() containing zeros and ones with a '1' indicating each attribute where the sigma coefficient should applied to the boundary.

Best wishes,
Mark

[fyrefly00]

In my case, setting the right hand side coefficient to zero did the trick perfectly. Thank you again for all your help-- I'll now mark this discussion as resolved.