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

[Anna Avdeeva]

Hi All,

I have successfully computed the secondary electric field for 3D case as described in the paper by Grayver&Buerg, 2014. These field is declared as BlockVector<double> solution, where the block correspond to real and imaginary parts.
Now I would like to compute total electric E (secondary plus primary) and magnetic field (coef* curlE)
1) I can use PostProcessor class to do this
template <int dim>
class Postprocessor : public DataPostprocessor<dim>
{
const SmartPointer<EPrimary<dim> > eprimary;
public:
Postprocessor (const SmartPointer<EPrimary<dim> > eprimary) :
eprimary (eprimary)
{
}
virtual
void
evaluate_vector_field( ...
...
}
But then I need to extract values of the total fields at set of receiver locations.
If I had total fields in the form similar to solution, I could use
for (unsigned int i = 0; i < receiver_locations.size (); ++i)
{
VectorTools::point_value (dof_handler, solution,
receiver_locations[i],
solution_at_receiver); ...
}
But since total fields are computed with help of postprocessor, I do not know how to extract values at receivers.

2) The other option is to follow the advice by Ross Kynch, and have loop over cells, quadrature points and dofs to get local contributions to total electric and magnetic field.
In this case I am not sure how then I get a global vector.

Please let me know what is the best way.

I would very much appreciate you reply,
Anna

[Wolfgang Bangerth]

Anna,

> I have successfully computed the secondary electric field for 3D case as described in the paper by Grayver&Buerg, 2014. These field is declared as BlockVector<double> solution, where the block correspond to real and imaginary parts.
> Now I would like to compute total electric E (secondary plus primary) and magnetic field (coef* curlE)
> 1) I can use PostProcessor class to do this
> template <int dim>
> class Postprocessor : public DataPostprocessor<dim>
> {
> const SmartPointer<EPrimary<dim> > eprimary;
> public:
> Postprocessor (const SmartPointer<EPrimary<dim> > eprimary) :
> eprimary (eprimary)
> {
> }
> virtual
> void
> evaluate_vector_field( ...
> ...
> }

Yes, this is definitely the way to create the total electric field or
magnetic field for purposes of visualization.

> But then I need to extract values of the total fields at set of receiver locations.

This is a separate issue because your receiver locations are typically
located at random locations, not necessarily at the vertices of the mesh
(where we generate visualization data). Evaluating something at random
locations is generally more complicated.

> If I had total fields in the form similar to solution, I could use
> for (unsigned int i = 0; i < receiver_locations.size (); ++i)
> {
> VectorTools::point_value (dof_handler, solution,
> receiver_locations[i],
> solution_at_receiver); ...
> }
> But since total fields are computed with help of postprocessor, I do not know how to extract values at receivers.

Correct. But you can use VectorTools::point_value/gradient() to evaluate
the values or gradient of your solution, and then repeat essentially the
same piece of code that you already have in your DataPostprocessor class.

Another option is to use the PointValueHistory class that combines
evaluating at arbitrary locations with using the DataPostprocessor
facilities.

Does this help?

Best
W.


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

[Anna Avdeeva]

Thank you for the prompt reply.

I will implement the first option with VectorTools::point_value/gradient() .

Could you please also point out to me where I can find clarifications to the following:

1) There are many versions of VectorTools::point_value function. At the moment I use

void VectorTools::point_value (

const DoFHandler< dim, spacedim > & dof, const VectorType & fe_function,
const Point< spacedim > & point, Vector< typename VectorType::value_type > & value )

Will the result be better if I use a version with mapping argument? 

2) I am interested to know how exactly the VectorTools::point_value/gradient() interpolate vector field values from quadrature points to arbitrary location. I suspect that it uses only one cell?

3) For some of the applications I am interested in the receivers are located at the interface between the air(or sea) and the ground. In other words receivers are at the boundary of the cell. Does it mean I should shift the receiver location slightly downwards, to insure that it lies inside the cell and to get more predictable result?

Finally, do you think using PointValueHistory class is a better way?

Thank you very much in advance,

Anna

[Wolfgang Bangerth]

> Could you please also point out to me where I can find clarifications to the
> following:
> 1) There are many versions of VectorTools::point_value function. At the moment
> I use
>
> void VectorTools::point_value (
> const DoFHandler< dim, spacedim > & dof, const VectorType & fe_function,
> const Point< spacedim > & point, Vector< typename VectorType::value_type > &
> value )
>
> Will the result be better if I use a version with mapping argument?

If your cells have curved boundaries, then yes. Otherwise, all mappings are
equivalent to not specifying a mapping at all.

> 2) I am interested to know how exactly the VectorTools::point_value/gradient()
> interpolate vector field values from quadrature points to arbitrary location.
> I suspect that it uses only one cell?

Yes, it finds the cell a point is in, and then interpolates the shape
functions on that cell.

> 3) For some of the applications I am interested in the receivers are located
> at the interface between the air(or sea) and the ground. In other words
> receivers are at the boundary of the cell. Does it mean I should shift the
> receiver location slightly downwards, to insure that it lies inside the cell
> and to get more predictable result?

Maybe. But if your solution is continuous, then it shouldn't matter what cell
you evaluate the point from. Of course, if your solution has a steep gradient
on one side of the interface, then the matter may be different.

> Finally, do you think using PointValueHistory class is a better way?

No. The only advantage is that it can make use of DataPostprocessor facilities.

[Anna Avdeeva]

Thank you very much for all your answers. 

[Guido Kanschat]

Dear Anna,

> I have successfully computed the secondary electric field for 3D case as described in the paper by Grayver&Buerg, 2014. These field is declared as BlockVector<double> solution, where the block correspond to real and imaginary parts.
> Now I would like to compute total electric E (secondary plus primary) and magnetic field (coef* curlE)

I am not sure what the terms "primary" and "secondary" fields mean. If
they are electric field and electric displacement, be aware that they
are in different function spaces!

> 2) The other option is to follow the advice by Ross Kynch, and have loop over cells, quadrature points and dofs to get local contributions to total electric and magnetic field.
> In this case I am not sure how then I get a global vector.

When you do this, please do not forget the de Rham complex: the curl
of a Nedelec element is a Raviart-Thomas element. So, if you do the
postprocessing, you should either choose such an element to store your
values, or a DG element containing the Raviart-Thomas space.

Best,
Guido


--
Prof. Dr. Guido Kanschat
Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
Universität Heidelberg
Im Neuenheimer Feld 368, 69120 Heidelberg