Repost Is the approach for the electrostatic Bi-linear form correct?

[phillip mobley]

Hello all,

I am unable to edit my original post. So I am creating a new post with the equations. I think that I solved issue with the equations not being displayed. The original posting will be deleted.

I have recently started to learn how to use Deal.II. I have been spending some going through some of the video tutorials for Deal.II, going through the coding tutorials, and watching some of the videos from the University of Michigan related to Finite Element Analysis. I have decided that I need to get my hands a bit dirty so that everything that I am reading and watching will start to sink.

So please, bear with me. I am still quite new to this field and I am still learning some of the basics.

With that being said, I am looking at creating my own electrostatic solver. For those unfamiliar, I am solving the PDEs:

For the electrostatic simulation, I would be solving for the Voltage scalar field and the Electric field (which is a vector).

So far, I have found that there are two paths that I can take to get my bilinear form. I have already attempted one path and I ran into a compile error with Deal.II. So, I am thinking that this bilinear form may not be appropriate. But I wanted to turn to the community here to see if my bilinear form has an issue or if my code has an issue.

For the first path, since V and E are my unknowns, I wanted both V and E to be on the left-hand side of their PDEs. With a bit of substitution, the equations that I would be solving turned into this:

Normally, there would be a - signed on the right-hand side of both equations. This will get dropped out because I am going to simplify the equation a bit further later down the road. In order to construct the bilinear form, I used the technique introduced in lecture 19. I am not going to show the entire process but rather the end result. If the entire process is needed, please post and I will respond with the complete step-by-step process that I performed.

then

This is how I created my bilinear form. Note that the p/e term dropped out because in this simplified simulation, there is no charge density. I create my simulation that is based on the tutorial step-20. There are alot of extra features in the tutorial that I did not add in to my simulation but rather focused on the FESystem portion. My solver is a CG Solver. Nothing fancy here. My first version of the solver can be found below

void SetupSolver(FESystem<2> &finiteElement, DoFHandler<2> &dofHandler, SparseMatrix<double> &masterMatrix, Vector<double> &masterRHS, Vector<double> &masterSolution)
{
    QGauss<2> quadrature(3);// Sets the form of the function
    FEValues<2> finiteElementValues(finiteElement, quadrature, update_values | update_gradients | update_JxW_values | update_quadrature_points);
    const unsigned int dofsPerCell = finiteElement.dofs_per_cell;
    const unsigned int numberQuadraturePoints = quadrature.size();

    FullMatrix<double> cell_matrix(dofsPerCell, dofsPerCell);
    Vector<double> rhsCell(dofsPerCell);

    vector<types::global_dof_index> localDofIndices(dofsPerCell);

    DoFHandler<2>::active_cell_iterator cell = dofHandler.begin_active();
    DoFHandler<2>::active_cell_iterator endc = dofHandler.end();
   
    const FEValuesExtractors::Vector electricField(2);
    const FEValuesExtractors::Scalar voltageField(0);

    cout << "Setting up Solver" << endl;

    for(; cell != endc; cell++)
    {
        // Calculates the values and gradients of the shape function for a particular cell
        finiteElementValues.reinit(cell);
        cell_matrix = 0;
        rhsCell = 0;

        cout << "working on cell: " << cell << endl;

        for(unsigned int qIndex = 0; qIndex < numberQuadraturePoints; qIndex++)
        {
            for(unsigned int i = 0; i < dofsPerCell; i++)
            {
                const double div_phi_i_e = finiteElementValues[electricField].divergence(i, qIndex);
                const Tensor<1, 2> grad_phi_i_v = finiteElementValues[voltageField].gradient(i, qIndex);
               
                for(unsigned int j = 0; j < dofsPerCell; j++)
                {
                    const Tensor<1, 2> phi_j_e = finiteElementValues[electricField].value(j, qIndex);
                    const Tensor<1, 2> grad_phi_j_v = finiteElementValues[voltageField].gradient(j, qIndex);
                   
                    cell_matrix(i, j) += ((div_phi_i_e * phi_j_e) + (grad_phi_i_v * grad_phi_j_v))  * finiteElementValues.JxW(qIndex);
                }
               
                rhsCell(i) += 0.0;
            }
        }

        cell->get_dof_indices(localDofIndices);

        for(unsigned int i = 0; i < dofsPerCell; i++)
        {
            for(unsigned int j = 0; j < dofsPerCell; j++)
            {
               masterMatrix.add(localDofIndices[i], localDofIndices[j], cell_matrix(i, j));
            }
        }

        for(unsigned int i = 0; i < dofsPerCell; i++)
        {
            masterRHS(localDofIndices[i]) += rhsCell(i);
        }

    }

    cout << "Finished setting up cells" << endl;

//    SolvingFunction test();

    map<types::global_dof_index, double> boundaryValues;

    cout << "Interperlating boundary conditions" << endl;

    VectorTools::interpolate_boundary_values(dofHandler, 1, ConstantFunction<2>(100.), boundaryValues);
    VectorTools::interpolate_boundary_values(dofHandler, 2, ZeroFunction<2>(), boundaryValues);
    VectorTools::interpolate_boundary_values(dofHandler, 0, ZeroFunction<2>(), boundaryValues);

    MatrixTools::apply_boundary_values(boundaryValues, masterMatrix, masterSolution, masterRHS);

    cout << "Finished setting up solver" << endl;
}

On line

cell_matrix(i, j) += ((div_phi_i_e * phi_j_e) + (grad_phi_i_v * grad_phi_j_v))  * finiteElementValues.JxW(qIndex);

I am getting the error:

error: no match for 'operator+' (operand types are 'dealii::Tensor<1, 2>' and 'dealii::Tensor<0, 2, double>::tensor_type {aka double}')

Which, from what the error is stating, the + operator can not take these two arguments.

This then leads me to reason that my bilinear form has an issue. Would this be a correct assumption? Or did I code something incorrectly? (This is my main issue right now)

Path 2:

I have also been considering this approach for creating the bilinear form where the PDEs that I want to solve are as follows:

Again, following the technique introduced in Lecture 19, I would simplify my PDEs into:

then

I was wondering if this would be a better approach in constructing the bi-linear form and solving for the electric field and voltage field.

Or are both methods equally valid? I have a feeling that the first method is not valid since it is essentially 2 of the same equations. What are your thoughts? Is there a different method that I can approach in solving the PDEs for the electric and voltage fields?

Side Notes:

There is a third equation which is the curl of the electric field is 0 for electrostatics. I decided that I am not using it since these two are enough. Is this an incorrect thinking and should I also account for this equation in my bilinear form?

I am using Deal.II to create a mesh which is a parallel plate capacitor. I am then solving for the electric and voltage fields that make up the parallel plates when the voltage and electric fields are 0 on the boundary. On the boundary of one of the plates, the Voltage is set to 100. On the other plate, the Voltage is set to 0.

For reference, I am attaching the source code of my simulation. There are more improvements that I will be making to the code once I get something workable. The code that is attached is a modified version of an earlier simulation where I was only solving for the voltage field in a parallel plate. This was around 2 years ago that I created this. I am simply expanding the code to also solve for the electric field. in short, it is rough. If you do see any other mistakes or areas of improvement, feel free to let me know.

Thank you

[phillip mobley]

As another side note, in this version of the simulation, the charge density (p or rho) will always be zero. Hence, the reason it drops off in my equations.

[Wolfgang Bangerth]

> So far, I have found that there are two paths that I can take to get my
> bilinear form. I have already attempted one path and I ran into a compile
> error with Deal.II. So, I am thinking that this bilinear form may not be
> appropriate.

That's not necessarily the right conclusion. Just because you chose the wrong
syntax (= compiler error) does not mean that you chose the logically wrong
approach.

> But I wanted to turn to the community here to see if my bilinear
> form has an issue or if my code has an issue.
>
> For the first path, since V and E are my unknowns, I wanted both V and E to be
> on the left-hand side of their PDEs. With a bit of substitution, the equations
> that I would be solving turned into this:

I don't see your inlined equations in my mail program, but assuming that you
are talking about the equation

(div E, phi_E) - (grad phi_V, grad V) = 0

then this is not correct. If you start out from the equation
div E - Delta V = 0
then you will want to multiply this (scalar) equation with a test function and
integrate by parts where necessary. But it's *one* test function, so you need
to end up multiplying *both* terms with the same test function phi_V, not
different test functions for the two terms.

That said, here are two questions:
* In essence, you want to solve the Laplace equation for V, but derive a mixed
form, which with your variables is usually stated as
E + grad V = 0
whereas you choose
div E + Delta V = 0
This raises the question of why you want to do it this way? You are
requiring more differentiability than necessary, negating the advantages of
the mixed form.

* If you chose to go with the formulation you have, I suspect that you want to
integrate by parts:
- (E, grad phi_V) - (grad phi_V, grad V) = 0
so that all derivatives are on the V test functions, not on E. Any good
reasons not to do so? I'm not sure that your formulation is uniquely
solvable at all.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[phillip mobley]

Hello Wolfgang,

Thank you for your reply. What you described above is the first approach that I took. There are some differences but overall, it is the same thing. I am not familiar with the form:

div E - Delta V = 0

but, I was able to figure out how you came this. Instead, I kept the equations separate and applied the technique described in lecture 19. For reference, the equations that I am referring to are:

-grad E = 0

delta V = 0

where E is the electric field vector.

I saw a new approach that I could take which I will describe here. Coming back to the original equations:

delta V = 0 

-grad V = E

And then applying lecture 19 technique:

(delta phi_V, V) - (grad phi_V, V) = (E, phi_V)

However, I run into the problem which I believe you briefly touched upon on your answer which is I am multiplying by a single test function; phi_V. So I am also thinking that this is the wrong approach.

There is a third equation that I have not touched yet which is the curl E = 0. If I were to account for this equation, then this would become the third approach to the problem. Thus the "complete" set of equations would be

delta V = 0 

-grad V = E

curl E = 0

If curl E = 0 is taken into account, then I would be working with the form:

(delta phi_V, V) - (grad phi_V, V) + (curl phi_E, E) = (phi_V, E)

Now, it would seem that i am working with two test functions, phi_V and phi_E. So it would appear that this form is more appropriate to work with. Is this the direction that I should be taking or is there another approach that I should consider?

[Jean-Paul Pelteret]

Dear Philip,

It looks to me like you’re mixing things up a little bit. In all of the equations that you presented there is a mixture of governing equations, assumed relationships and the consequence of combining them. From what I can tell you’re actually trying to use the same equation more than once which is problematic. (I hope that I have clearly understood what you’re trying to do, otherwise what follows is a bit of a waste of space and time…)

I think that if you work from the irrefutable, then there should be little confusion. So let me see if I can help by sketching out electrostatics from the basics, as I know it:

The starting point is Maxwell’s equations. Let’s immediately simplify things by saying that we already assume a quasi-static state — this means that the electric and magnetic fields are now decoupled. We then get

(1) div d = \rho^{f}  [Gauss’ law] 

(2) curl e = 0          [Faraday’s law]

on B. Here d is the electric displacement vector, e is the electric field vector, and \rho^{f} is the free charge density. 

There is a further fundamental constitutive law that links d and e, namely

(3) d = \epsilon_{0} \epsilon_{r} e

where \epsilon_{0} is the electric permittivity of free space and \epsilon_{r} is the relative permittivity of the material. Here I have assumed that this material is linearly polarisable, so that is that the polarisation is in the same direction as the electric field (that’s where the scaling factor \epsilon_{r} comes from).

There are also the continuity conditions for Maxwells equations, which for this simplified scenario would be

(4) n x [[e]] = 0

(5) n . [[d]] = 0    (assuming no surface charges)

on dB. Here n is a surface normal and [[ ]] denotes the jump of a quantity, “.” the dot product and “x” the cross product.

So, going back to Maxwell’s equations, you now need to choose which variable is going to be the primary variable, and which one follows from it (so use (3) in (1) or (2)). Since in general you wish to include source terms \rho^{f}, it is appropriate to choose the electric field to be the primary variable. what comes next is that you substitute (3) into (1) thereby eliminating the electric displacement from the equation:

(6) div [\epsilon_{0} \epsilon_{r} e] = \rho^{f} 

—>

(7) div [e] = \rho^{f} / [\epsilon_{0} \epsilon_{r}]

So now (7) represents an amended form of (1) that needs to be solved, but we still need to satisfy (2). So what we do is to exploit the identity

(a) curl (grad (s)) = 0 

for all arbitrary scalars s. So we postulate the existence of an electric scalar potential field V that is linked to the electric field by

(8) e = -grad V

Using the identity (a), you can see that this satisfies (2):

(b) curl e = curl (- grad V) == 0

Now we can put (8) into (7), knowing that if we solve 

(9) div [-grad V] = \rho^{f} / [\epsilon_{0} \epsilon_{r}]

—>

(10) - delta V = \rho^{f} / [\epsilon_{0} \epsilon_{r}]

for V then we automatically satisfy (2). Remember that (10) is nothing other that (1) in disguise. So this actually satisfies both of Maxwell’s electrostatic governing equations at the same time (since we solve for the electric scalar potential).

Now, about those continuity conditions… So when you discretise the field V using continuous finite elements, then (4) is automatically satisfied because the continuity of the solution ensures that there is no jump in the solution V, and therefore no tangential jump in grad V. You can see this from

(11) 0 = n x [[e]] =(8)= n x [[-grad V]] = n x grad [[V]]

Finally, (5) is satisfied automatically through this same assumption

(12) 0 = n . [[d]] =(3)= n . [[\epsilon_{0} \epsilon_{r} e]] =(8)= \epsilon_{0} \epsilon_{r} n . [[-grad V]] = -\epsilon_{0} \epsilon_{r} n . grad[[V]]

So, in summary, the strong form of the governing equation to be implemented is (10), with the solution for V leading to the electric field by (8) and subsequently the electric displacement by (3). The continuity conditions (4,5) are satisfied if the solution for V is continuous. 

Does this make sense, and does it help clear up things? 

Best regards,

Jean-Paul

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[phillip mobley]

Hello Jean,

Thank you very much for your derivation. It does help clear up many questions that I had. I greatly appreciate it!

The strong form that I will be using are simplified versions of the form that you derived in the description.

- delta V = \rho^{f} / [\epsilon_{0} \epsilon_{r}]

becomes

- delta V = \rho / [\epsilon_{0}]

when f = 1 and when the material is free space (eplison_r = 1 in air).

I am starting off with a more simpler simulation to learn the ropes better. For now, I am assuming that the charge density will be zero (\rho). The strong form then becomes:

- delta V = 0

Later, I plan on expanding this for the full delta V equation (- delta V = \rho^{f} / [\epsilon_{0} \epsilon_{r}]). Since I am just starting out, I would like to keep things "simple".

I was originally thinking to obtain the electric field (e = -grad V) separately. First, i would solve the simulation for V using the equation - delta V = 0. Then solve for the electric field. But, I do not know how I can do this in deal.ii. In fact, I wasn't even sure if this was possible. Am I able to solve for the electric field using the results I obtained when I run the solver for the voltage? If so, would you be able to point me to some resources on how I can implement this using deal.ii?

Again, at the time, I didn't think this was possible. So, I decided to treat the equations as 2 separate equations and solve for them at the same time. Which lead me to where I am at currently.

Right now, I am not concerned about solving the electric displacement. Unless, it is necessary for me to solve this.

You have also confirmed my idea that solving for

- delta V = 0

div E = 0

will not work since they are the same equation.

Basically, from the description that you provided, I should solve for the voltage first and then solve for the electric field using the results I obtained from the first simulation. If this is the case, would you be able to point to me some resources that I can use in order to implement this in deal.ii? Concerning the continuity equations, should I check my results against them to make sure the results are valid?

Thank you

[phillip mobley]

Hello Jean, I was going through the deail.ii documentation and I came across this class:

https://dealii.org/9.0.0/doxygen/deal.II/classDataPostprocessorVector.html

Is this what you we referring to?

On Wednesday, September 5, 2018 at 3:17:14 AM UTC-4, Jean-Paul Pelteret wrote:

[Jean-Paul Pelteret]

Yes, either that data post-processor or via fe_value.get_gradients().