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
<https://dealii.org/developer/doxygen/deal.II/classFEValuesBase.html#ad1f4e0deb5d982e8172d82141c634a67>(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
> On 23. Apr 2021, at 12:03, Sylvain Mathonnière
> <[email protected]> wrote:
>
> 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.
>
> El viernes, 23 de abril de 2021 a la(s) 10:07:23 UTC+2, Sylvain Mathonnière
> escribió:
> 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
>
> <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,
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <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 [email protected]
> <mailto:[email protected]>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/e29faa8d-fe5a-4cd0-a6f3-561bf4686bf9n%40googlegroups.com
>
> <https://groups.google.com/d/msgid/dealii/e29faa8d-fe5a-4cd0-a6f3-561bf4686bf9n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/D9AF46BA-7E42-4E1A-9CB6-E11CD12EB062%40gmail.com.
