integrating over boundaries

[Sylvain Mathonnière]

Hello everyone !

I have a very basic question but somehow I did not manage to find the answer in tutorials, neither looking at the documentation or in the this user group. Or more probably I did not understand it when I saw it since I am fairly new to deal.II and to C++.

My main goal : I am trying to integrate over a boundary (a line since I am in 2D) the heat flux to get the total heat flux through my boundary.

What I did :

1 - I looked at this documentation https://www.dealii.org/current/doxygen/deal.II/classDataPostprocessorVector.html to calculate first the heat flux on my boundary which worked fine.

2 - I tried to integrate the norm of the heat flux over my boundary in the function evaluate_scalar_field() function but I could not. My best effort is this :

template <int dim>

class HeatFluxPostprocessor : public DataPostprocessorVector<dim>

{

public:

  HeatFluxPostprocessor ()

    :

    DataPostprocessorVector<dim> ("_k_grad_T",

                                  update_gradients | update_quadrature_points)

  {}

  virtual void evaluate_scalar_field(

    const DataPostprocessorInputs::Scalar<dim> &input_data,

    std::vector<Vector<double> >               &computed_quantities) const

    {

      // I recover the value of conductivity

      MaterialData<dim> material_data;

      const typename DoFHandler<dim>::cell_iterator current_cell = input_data.template get_cell<DoFHandler<dim>>();

      FullMatrix<double>  conductivity = material_data.get_th_conductivity(current_cell->material_id());    

      double heat_flux_cell_norm = 0.;

      for (auto face_index : GeometryInfo<dim>::face_indices())

      {

          if (current_cell->face(face_index)->boundary_id() == 2)

          {

            double long_cell = 0;

            long_cell = (current_cell->face(face_index)->vertex(0) - current_cell->face(face_index)->vertex(1)).norm();

            for (unsigned int d=0; d<dim; ++d)

            {

              heat_flux_cell_norm += -conductivity[d][d]*input_data.solution_gradients[face_index][d] * 

                                      conductivity[d][d]*input_data.solution_gradients[face_index][d];

            }

            std::cout << "total heat flux = " << heat_flux_cell_norm << std::endl;

            heat_flux_cell_norm *= long_cell;

          }

      }

Why I could not go further:

The problem is that since the evaluate_scalar_field() loops over the cells internally (I believe. Or is it add_data_vector() ?), I cannot sum my heat_flux_cell_norm double to obtain my integral.

What I tried :

So then I tried creating a separate function in my class HeatFluxPostProcessor and do it myself. This is my best attempt :

  void output_heat_flux(const DataPostprocessorInputs::Scalar<dim> &input_data)

    {

    MaterialData<dim> material_data;

    const typename DoFHandler<dim>::cell_iterator current_cell = input_data.template get_cell<DoFHandler<dim>>();

    FullMatrix<double>  conductivity = material_data.get_th_conductivity(current_cell->material_id());

    double heat_flux_tot = 0;

    for (current_cell.begin(); current_cell.end(); ++current_cell)

      {

        double long_cell = 0;

        double heat_flux_cell_norm = 0;

            for (auto face_index : GeometryInfo<dim>::face_indices())

            {

                if (current_cell->face(face_index)->boundary_id() == 2)

                {

                  long_cell = (current_cell->face(face_index)->vertex(0) - current_cell->face(face_index)->vertex(1)).norm();

                  for (unsigned int d=0; d<dim; ++d)

                  {

                    heat_flux_cell_norm += - conductivity[d][d]*input_data.solution_gradients[face_index][d] * 

                                             conductivity[d][d]*input_data.solution_gradients[face_index][d];

                  }

                }

            }

          heat_flux_tot += heat_flux_cell_norm * long_cell;

          std::cout << "total heat flux = " << heat_flux_tot << std::endl;

      }

    }

Why I could no go further:

However I actually do not know how to get the argument

const DataPostprocessorInputs::Scalar<dim> &input_data from my ouptut data to launch this function... I tried to check how evaluate_scalar_field() does but could not work my way through the documentation to identify it (again I am fairly new to C++). What should be the generic argument for this function based on my solution vector (Temperature) ?

I am probably missing something pretty obvious to do what I want since it is rather easy, any help would be appreciated. I believe I could integrate more properly using quadrature point but I am not too sure how to do it. Again any help/directions about it would be greatly appreciated.

Best regards,

[Sylvain Mathonnière]

By the way, I am aware than I forgot to take the square root of heat_flux_cell_norm. I am also aware that using input_data.solution_gradients[face_index][d]  is weird in this context since the first indices should refer to a vertex indices of the cell and not the face index and it only works because quadrangle have 4 faces and 4 vertices but so far this was not my main issue but it is definitively on my list of fixing.

[Jean-Paul Pelteret]

Hi Sylvain,

Thanks for detailing your attempts so thoroughly. With this

      // I recover the value of conductivity

      MaterialData<dim> material_data;

      const typename DoFHandler<dim>::cell_iterator current_cell = input_data.template get_cell<DoFHandler<dim>>();

      FullMatrix<double>  conductivity = material_data.get_th_conductivity(current_cell->material_id());    

      double heat_flux_cell_norm = 0.;

      for (auto face_index : GeometryInfo<dim>::face_indices())

      {

          if (current_cell->face(face_index)->boundary_id() == 2)

          {

            double long_cell = 0;

            long_cell = (current_cell->face(face_index)->vertex(0) - current_cell->face(face_index)->vertex(1)).norm();

            for (unsigned int d=0; d<dim; ++d)

            {

              heat_flux_cell_norm += -conductivity[d][d]*input_data.solution_gradients[face_index][d] * 

                                      conductivity[d][d]*input_data.solution_gradients[face_index][d];

            }

            std::cout << "total heat flux = " << heat_flux_cell_norm << std::endl;

            heat_flux_cell_norm *= long_cell;

          }

      }

you’ve got most of the structure for your calculation completely correct, but embedding it within a DataPostprocessor is not what you want to do. These classes are there to help visualise quantities via DataOut (i.e. do some per-cell calculations), rather than computing some boundary integral like you want. Instead, you should encapsulate the above within your own function and the standard DoFHandler active cell loop. You would then replace “input_data” with some FEFaceValues object (reinitialised on the given cell and face_index) and then call "fe_face_values.get_function_gradients(solution_vector, gradients);” to store the solution gradients at each face quadrature point and from that compute the heat flux at each face quadrature point. But I think that you’re also just wanting to consider the component of the heat flux that is normal to the boundary, right? 

I think that your calculation should have this sort of structure:

double total_heat_flux = 0.0;

FEFaceValues<dim> fe_face_values(fe, face_quadrature, update_gradients | update_normal_vectors | update_JxW_values);

std::vector<Tensor<1,dim>> temperature_gradient(face_quadrature.size());

  for (const auto cell : dof_handler.active_cell_iterators())

    for (const auto face_index : GeometryInfo<dim>::face_indices())

      If (… at boundary of interest …)

      {

       fe_face_values.reinit(cell, face_index);

       fe_face_values.reinit.get_solution_gradents(solution_vector, temperature_gradient);

        for (const auto q_point : fe_face_values.quadrature_point_indices())

        {

          const Tensor<1,dim> heat_flux_q = -conductivity*temperature_gradient[q_point]; // Q = -K.Grad(T)

          total_heat_flux += heat_flux_q * fe_face_values.normal_vector(q_point) * fe_face_values.JxW(q_point); // q = Q.N. ; q_total = \int q dA

      }

    }

If I’ve interpreted things correctly then that (or something close to it) should achieve what you’re looking to do.

I hope that this help!

Best,

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.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/e29faa8d-fe5a-4cd0-a6f3-561bf4686bf9n%40googlegroups.com.

[Sylvain Mathonnière]

Thank you very much for your help ! Indeed it is simpler outside of the DataPostprocessor class.

I compared the integration using quadrature points and using constant value along cell face (as I initially intended) and got similar results which comfort me that I have the correct implementation.

The key enabler was really to get the gradient of the solution vector using "fe_face_values.get_function_gradients(solution_vector, gradients)" for me. Thanks again for the help.

Best regards,

Sylvain Mathonnière