Computing the curl of a solution vector field obtained from Nedelec elements

[David Fernández]

Hello all,

I am using Nedelec elements to solve a for a vector field. Once I have the solution I want to compute the "curl" of this solution and then plot it.
However I am having trouble computing the curl of the solution.
The code that I am using to compute this curl reads the solution defined in a vector  "Vector<double> solution" and attempts stores the computed curl in a vector of the same size and type as shown in the method that follows (this workflow is similar to the one proposed by Ross Kynch in https://groups.google.com/forum/#!searchin/dealii/nedelec/dealii/1g3YSUdPSGY/gCWwidzhUrgJ):

template <int dim>
void MaxwellProblem<dim>::compute_Curl_field()
{
    QGauss<dim> quadrature_formula(quad_order);
    const unsigned int n_q_points = quadrature_formula.size();

    FEValues<dim> fe_values(fe, quadrature_formula,
                            update_values | update_gradients | update_quadrature_points
                            | update_JxW_values);

    const unsigned int dofs_per_cell = fe.dofs_per_cell;

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

    // storage for computed H field solution:
    Tensor<1, dim> curlsol(n_q_points);

    typename DoFHandler<dim>::active_cell_iterator cell =
        dof_handler.begin_active(), endc = dof_handler.end();
    for (; cell != endc; ++cell)
    {
        fe_values.reinit(cell);

        // Determine the global indices of the cell Dof:
        cell->get_dof_indices(local_dof_indices);
        if (cell == dof_handler.begin_active())
            std::cout << "local_dofs_size: " << local_dof_indices.size() << std::endl;

        // Calculate the values of curl of the fe solution:
        // Loop over quad points to calculate solution:
        for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
        {
            curlsol[q_point] = 0.0;
            // Loop over DoFs to calculate curl of solution @ quad point
            for (unsigned int i = 0; i < dofs_per_cell; ++i)
            {
                // Construct local curl value @ quad point integrating
                curlsol += solution(local_dof_indices[i])
                           * fe_values[A_re].curl(i, q_point)
                           * fe_values.JxW(q_point);
            }
        }

        // Store the curl of the solution in H_solution
        for (unsigned int i = 0; i < n_q_points; ++i)
            H_solution(local_dof_indices[i]) = curlsol[i];
        // \-> This would seems to be included in the loop over the quadrature
        //     points, but then how would it be added into the global H_solution
        //     vector of curl's.
    }

}

The code used to output both the solution and its curl follows:

template <int dim>
void MaxwellProblem<dim>::output_results() const
{
    deallog << "Generating output... ";

    DataOut<dim> data_out;
    data_out.attach_dof_handler(dof_handler);

    std::ostringstream filename;
    filename << "solution-MagnetostaticProblem" << ".vtu";
    std::ofstream output(filename.str().c_str());

    std::vector<std::string> solution_names;
    solution_names.push_back("A_x");
    solution_names.push_back("A_y");
    solution_names.push_back("A_z");

    std::vector<std::string> curl_solution_names;
    curl_solution_names.push_back("H_x");
    curl_solution_names.push_back("H_y");
    curl_solution_names.push_back("H_z");
    data_out.add_data_vector(H_solution, curl_solution_names);

    data_out.build_patches();

    data_out.write_vtu(output);
}

Has anyone computed the curl of a solution generated with Nedelec elements and plotted them?

Regards,

David

[Ross Kynch]

Hi David,

You're close with that code but have mixed up the quadrature loop. The following will calculate the curl of the solution at each quad point within a local cell:

// Loop over quad points to calculate solution:

for (unsigned int q_point=0; q_point<n_q_points; ++q_point)

{

  // Loop over DoFs to calculate curl of solution @ quad point

  curlsol_re=0.0;

  for (unsigned int i=0; i<dofs_per_cell; ++i)

  {

    // Calculate contribute @ this quad point for this local DoF:

    curlsol_re += solution(local_dof_indices[i])*fe_values[A_re].curl(i,q_point);

  }

  // Store local curl solution for this quad point:

  local_curl(q_point) = curlsol_re;

}

If you wish to construct a global vector of points you can do so, but remember you must factor in the multiplicities of the quadrature points you're using. You may just want to divide by the multiplicity so you take equal weight from all of the cells which the point belongs to, or you may wish to weight things. 

Another point where you've written:

H_solution(local_dof_indices[i]) = curlsol[i];

This doesn't really make sense in the context of a Nedelec element as the DoFs are not at collocation points. They're associated with edges, faces and interiors on the cell, so you're best off sticking to using quadrature points (or some other set of points) when you wish to compute things in the whole domain.

For visualisation purposes, another option is to use a derived DataPostprocessor class and use the following to compute the curl within derived_postprocessor::compute_derived_quantities_vector(....):

computed_quantities[q](0) = (duh[q][2][1]-duh[q][1][2]);

computed_quantities[q](1) = (duh[q][0][2]-duh[q][2][0]);

computed_quantities[q](2) = (duh[q][1][0]-duh[q][0][1]);

I'd personally say this is the best way to do it, but there are shortcomings if you need to compute which require material information on a cell into the class - you'd have to hack your way around it.

Anyway, hope this helps

Ross

[Ross Kynch]

One addition I failed to mention:

You've declared:

// storage for computed H field solution:
    Tensor<1, dim> curlsol(n_q_points);

but for my code it would be:

    Tensor<1, dim> curlsol;

as I'm just using it as temporary storage as I run through all the DoFs.

R

[David Fernández]

Hi Ross,

Thanks for the reply!

I still have a few questions:
1. In the inner loop I am doing as you suggest but I am also adding the "fe_values.JxW(q_point)" as this is needed to compute the curl in the real cell from the referece cell values given by the "fe_values[A_re].curl(i, q_point)" function as follows:

temp_curlsol += solution(local_dof_indices[i])  * fe_values[A_re].curl(i, q_point)  * fe_values.JxW(q_point);

I understand your point in not storing these values in a global vector but only in a temporal Tensor<1,dim> so that is what temp_curlsol is.

2. This approach would also help me, but I am not sure how to go about it. How can I extract these differentials (fe_values.shape_grad_component (...)?)
per quadrature point and then compute the curl as in the code snippet you suggest?
Also, I would like to evaluate the curl in the center of the element (cell->center), so I would want to evaluate the curl from the shaper functions at this point, is there support in deal.ii for this?
Thanks!

David

[Ross Kynch]

Hi David, sorry for the slow reply - have had a bit of a break over Christmas!

1. In the inner loop I am doing as you suggest but I am also adding the "fe_values.JxW(q_point)" as this is needed to compute the curl in the real cell from the referece cell values given by the "fe_values[A_re].curl(i, q_point)" function as follows:

temp_curlsol += solution(local_dof_indices[i])  * fe_values[A_re].curl(i, q_point)  * fe_values.JxW(q_point);

JxW shouldn't be used here - the shape function/curl doesn't need to be transformed because fe_values[A_re].curl(i, q_point) returns the curl in real space already. See the discussion in the last few posts here: https://groups.google.com/forum/#!topic/dealii/ZJqmZgObysw

 

2. This approach would also help me, but I am not sure how to go about it. How can I extract these differentials (fe_values.shape_grad_component (...)?)
per quadrature point and then compute the curl as in the code snippet you suggest?

The code snippet should be within a function inside a derived DataPostprocessor class, see step-29/32/33 for examples. The duh variable is passed within this class. It contains the derivatives of the shape functions as you suspected, but you don't need to generate it.

Also, I would like to evaluate the curl in the center of the element (cell->center), so I would want to evaluate the curl from the shaper functions at this point, is there support in deal.ii for this?
Thanks!

This is possible, I think. You'd have to evaluate the derivatives of the shape functions  at the centre point, compute the components of the curl (most likely you'll be able to compute the gradient which you can use) and then sum contributions of them multiplied by the corresponding solution entry.

I think the function shape_grad(int i, Point<dim> &p) inside the FiniteElement class will compute the gradient of the ith shape function at the point p within a particular cell, so should be all you need.

Happy new year!

Ross

[Wolfgang Bangerth]

David,
Ross has already given excellent answers but let me give a different spin on
this nonetheless:

> I am using Nedelec elements to solve a for a vector field. Once I have the
> solution I want to compute the "curl" of this solution and then plot it.
> However I am having trouble computing the curl of the solution.
> The code that I am using to compute this curl reads the solution defined in a

> vector *"Vector<double> solution"* and attempts stores the computed curl in a

> vector of the same size and type as shown in the method that follows (this
> workflow is similar to the one proposed by Ross Kynch in
> https://groups.google.com/forum/#!searchin/dealii/nedelec/dealii/1g3YSUdPSGY/gCWwidzhUrgJ):

> [...]

There are a number of misunderstandings in your code. First, you multiply by
JxW. The JxW values correspond to the dx in an integral -- you need it if you
integrate over something, but not if you want to transform one quantity into
another.

Second, you say that you want to store the curl of the solution in a vector of
the same kind. But that can't work, and it is probably best explained using
just a scalar, 1d solution: if you have a piecewise linear function, then it
is described by a vector of values at the nodal points. But take the gradient
and the function becomes a piecewise constant function that can be described
by the value of the derivative at the cell center, for example. That already
explains why you won't *naturally* be able to store the gradient in a vector
"like" that of the solution. You will also notice that the number of nodes and
cells is different, so the dimensionality of the original function and its
gradient are different -- again a reason why storing it in the same vector
won't work.

The same is true when you take the curl of a vector field. Even if the vector
field is continuous, the curl is not. And the curl has a different polynomial
degree. So this is not going to work.

As pointed out by Ross, the correct way to go about it if you want to just
compute the curl for the purpose of visualization is to use the
DataPostprocessor framework. Specifically, the curl being a vector quantity,
you want to derive a class from DataPostprocessorVector and overload the
compute_derived_quantities_vector() function. It receives, among other
arguments, the tensor of derivatives of the solution. The curl just requires
you to collect a subset of these derivatives.


Thinking about this some more, this would actually make for a nice addition to
the library. If you send me the class you get out of this, then I'll look into
putting it into the library.

Best
W.

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

[SebG]

Hi,

I am currently experimenting with curl-conforming elements and their application in electromagnetism. I am using the Nedelec element FE_Nedelec and I have difficulties to understand the output and some other features:

1. I realized that VectorTools::project only works if QGauss<dim>(p) with p >= 2 is used. AssertNoZerosOnDiagonal fails in precondition_SSOR. If I am using an FE_System<dim>(FE_Q<dim>(1), dim) as the base element, everything is fine. I don't know if this qualifies as a bug. But there is no indication in the documentation of certain requirements on the quadrature applied.
2. I tried to output the curl of the vector, A, to a vtk file and used a class derived from DataPostprocessor<dim> as suggested above. The results do not look very promising if I compare them with the output of an analytic curl of A, see screenshots attached. However the output of the potential itself works as expected.
   
3. Because of the problem related to the output of the curl I tried to compute some errors. The projection error of the function value itself is converging pretty well. I have tried to compute the L2-error of curl(A_h) - curl(A) where A_h is the discrete approximation of A and A is the exact potential. The curl error is not converging if I am refining the mesh. But I am not sure if there is a mistake in my computation of the cellwise_error.

I am aware that there is a long discussion on the github page on the implementation of the Nedelec element. It would be great to get an idea what works with the current FE_Nedelec of version 9.0.1 and what doesn't.

I think that a reliable output of the curl of the solution computed with FE_Nedelec is required because in electromagnetism you want to look at the magnetic field, i.e. curl(A), instead of the potential, A. I can't tell if the problem regarding the output is related to FE_Nedelec, the DataPostprocessor interface or the vtk file format. It would also be great to get a suggestion what goes wrong here.

Best wishes,

Sebastian

[Alexander]

You shouldn't use FE_Nedelec with non-standard meshes like hyper_ball. Try FE_NedelecSZ instead. Read their docs for more information.

Best
Alexander

[Wolfgang Bangerth]

Sebastian,

> I am currently experimenting with curl-conforming elements and their
> application in electromagnetism. I am using the Nedelec element FE_Nedelec and
> I have difficulties to understand the output and some other features:
>

> 1. I realized that VectorTools::project only works if QGauss<dim>(p) with p

> >= 2 is used. AssertNoZerosOnDiagonal fails in precondition_SSOR. If I am
> using an FE_System<dim>(FE_Q<dim>(1), dim) as the base element, everything
> is fine. I don't know if this qualifies as a bug. But there is no
> indication in the documentation of certain requirements on the quadrature
> applied.

There is a statement, but you're right that it isn't quite clear. I've
proposed this here:
https://github.com/dealii/dealii/pull/7587


> 2. I tried to output the curl of the vector, A, to a vtk file and used a

> class derived from DataPostprocessor<dim> as suggested above. The results
> do not look very promising if I compare them with the output of an
> analytic curl of A, see screenshots attached. However the output of the
> potential itself works as expected.

It's hard to see anything in these figures if you don't know what it
represents :-) I looked at the code briefly and it looks correct (though I'm
not sure about the sign of each entry of the curl.) Try this on a simple case
first where you take a simple vector-valued function (say, one whose curl is a
constant vector), project it onto the mesh, and then output both the vector
potential and the curl.


> 3. Because of the problem related to the output of the curl I tried to

> compute some errors. The projection error of the function value itself is
> converging pretty well. I have tried to compute the L2-error of curl(A_h)
> - curl(A) where A_h is the discrete approximation of A and A is the exact
> potential. The curl error is not converging if I am refining the mesh. But
> I am not sure if there is a mistake in my computation of the cellwise_error.

Does the potential A_h itself converge at the correct order? If yes, then its
curl should converge at one order less. But since it's a derived quantity,
make sure first that the quantity from which it is derived (i.e., the primary
variable A_h) converges correctly and with the correct order. If it does and
the curl doesn't, then you know that the problem lies with the way you compute
the curl, not the way you compute A_h.

Best
W.

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

[SebG]

Hi Wolfgang,

sorry about not getting back to this for while. I have addressed your second point and created an example that is much easier to understand. Now the geometry is the unit cube. The vector field is A = [0 , 0 , y] and then curl(A) = [1 , 0 , 0].

I have attached three screenshots which compare the analytic result (left-hand side) to the numerical approximation (right-hand side). The approximation of the field, A, seems fine except for a strange wiggle of the magnitude on the finest mesh, see A_magnitude.png. From looking at the graphical output in curlA_glyphs.png, the approximation of the curl of A does not work well.

Regarding your third point, I have created the following convergence table:

cycle cells dofs  field               curl
                  error          rate error    rate
    0     1    54 0.1866            - 1.4410      -
    1     8   300 0.1967         0.95 3.0428   0.47
    2    64  1944 0.1576         1.25 4.8742   0.62
    3   512 13872 0.1429         1.10 8.8342   0.55

Regarding the field everything looks good. For the curl, the order seems to be one half although the error increases. This is a bit surprising for me because the overall error increases with the number of dofs. Apart from this, I am asking myself if the errors should not be much smaller because the field is linear and its curl is constant. Hence, these functions should be an element of the approximating finite element space. If I am not wrong, a single element should be sufficient for an approximation in this case. With FE_NedelecSZ the results are similar.

I would be grateful to get an advice if there is something wrong with my interpretation of the output or with the code.

Best,

Sebastian

[Wolfgang Bangerth]

On 2/12/19 6:15 PM, SebG wrote:
>
> sorry about not getting back to this for while. I have addressed your second
> point and created an example that is much easier to understand. Now the
> geometry is the unit cube. The vector field is A = [0 , 0 , y] and then
> curl(A) = [1 , 0 , 0].
>
> I have attached three screenshots which compare the analytic result (left-hand
> side) to the numerical approximation (right-hand side). The approximation of
> the field, A, seems fine except for a strange wiggle of the magnitude on the
> finest mesh, see A_magnitude.png. From looking at the graphical output

> in**curlA_glyphs.png, the approximation of the curl of A does not work well.

>
> Regarding your third point, I have created the following convergence table:
>
> cycle cells dofs  field               curl
>                   error          rate error    rate
>     0     1    54 0.1866            - 1.4410      -
>     1     8   300 0.1967         0.95 3.0428   0.47
>     2    64  1944 0.1576         1.25 4.8742   0.62
>     3   512 13872 0.1429         1.10 8.8342   0.55
>
> Regarding the field everything looks good. For the curl, the order seems to be
> one half although the error increases. This is a bit surprising for me because
> the overall error increases with the number of dofs. Apart from this, I am
> asking myself if the errors should not be much smaller because the field is
> linear and its curl is constant. Hence, these functions should be an element
> of the approximating finite element space. If I am not wrong, a single element
> should be sufficient for an approximation in this case. With FE_NedelecSZ the
> results are similar.

Seb -- as you may have noticed, none of the main developers work on curl
problems with the Nedelec elements, and so it is difficult to get definitive
answers about questions.

If you're still chasing answers, let's first focus on the potential (which I
think the table above labels as "field error". I don't know how the error rate
is computed, but the numbers for the rates can't be correct: I think they show
a reduction *factor*, not a rate. The rate is the exponent in
error = C h^r
where from one line to the next, h is reduced by a factor of 2. So for the
last two rows, for example, you need to have

0.1429 = (1/2)^r 0.1576

or r = -log(0.1429/0.1576)/log(2) = 0.14.

So what I see here is that the convergence of the potential (field) is already
not correct unless you happen to have a very non-smooth right hand side of the
equation you are trying to solve.

When you later compute the error of the curl of the potential, one typically
expects that the convergence rate is one lower for the derivative than the one
for the field itself. This would suggest that you would get a rate of -0.86 --
which turns out to be almost exactly what you have in your table above when
one computes the rate with the formula above.