Skip to content

Commit

Permalink
Merge pull request #14258 from peterrum/DataOutResample_multiple_comp…
Browse files Browse the repository at this point in the history
…onents

DataOutResample for multiple components
  • Loading branch information
peterrum committed Sep 15, 2022
2 parents 11b9eb1 + d9230d5 commit d13fb96
Show file tree
Hide file tree
Showing 6 changed files with 1,003 additions and 39 deletions.
6 changes: 6 additions & 0 deletions doc/news/changes/minor/20220912Munch_1
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
Improved: The functions VectorTools::point_values()/::point_gradients() and
the class Utilities::MPI::RemotePointEvaluation now allow to specify the first
component to be selected. The feature is used in the class DataOutResample
to output multi-component solutions.
<br>
(Peter Munch, Magdalena Schreter, 2022/09/12)
52 changes: 39 additions & 13 deletions include/deal.II/numerics/vector_tools_evaluate.h
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,12 @@ namespace VectorTools
* cache, dof_handler_2, vector_2);
* @endcode
*
* The function also works with FiniteElement objects with multiple
* components. If one is interested only in a range of components, one can
* select these by the parameters @p first_selected_component and
* @p n_components. For further details on supported FiniteElement
* objects, see the documentation of FEPointEvaluation.
*
* The function can also be used to evaluate cell-data vectors. For this
* purpose, one passes in a Triangulation instead of a DoFHandler and a
* vector of size Trinagulation::n_active_cells() or a vector, which
Expand Down Expand Up @@ -142,7 +148,8 @@ namespace VectorTools
const VectorType & vector,
const std::vector<Point<spacedim>> & evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg,
const unsigned int first_selected_component = 0);

/**
* Given a (distributed) solution vector @p vector, evaluate the values at
Expand Down Expand Up @@ -171,7 +178,8 @@ namespace VectorTools
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg,
const unsigned int first_selected_component = 0);

/**
* Given a (distributed) solution vector @p vector, evaluate the gradients at
Expand Down Expand Up @@ -200,7 +208,8 @@ namespace VectorTools
const VectorType & vector,
const std::vector<Point<spacedim>> & evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg,
const unsigned int first_selected_component = 0);

/**
* Given a (distributed) solution vector @p vector, evaluate the gradients at
Expand Down Expand Up @@ -228,7 +237,8 @@ namespace VectorTools
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg,
const unsigned int first_selected_component = 0);



Expand All @@ -252,11 +262,13 @@ namespace VectorTools
const VectorType & vector,
const std::vector<Point<spacedim>> &evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const EvaluationFlags::EvaluationFlags flags)
const EvaluationFlags::EvaluationFlags flags,
const unsigned int first_selected_component)
{
cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);

return point_values<n_components>(cache, mesh, vector, flags);
return point_values<n_components>(
cache, mesh, vector, flags, first_selected_component);
}


Expand All @@ -277,11 +289,13 @@ namespace VectorTools
const VectorType & vector,
const std::vector<Point<spacedim>> &evaluation_points,
Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const EvaluationFlags::EvaluationFlags flags)
const EvaluationFlags::EvaluationFlags flags,
const unsigned int first_selected_component)
{
cache.reinit(evaluation_points, mesh.get_triangulation(), mapping);

return point_gradients<n_components>(cache, mesh, vector, flags);
return point_gradients<n_components>(
cache, mesh, vector, flags, first_selected_component);
}


Expand Down Expand Up @@ -396,6 +410,7 @@ namespace VectorTools
const VectorType & vector,
const UpdateFlags update_flags,
const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
const unsigned int first_selected_component,
const std::function<
value_type(const FEPointEvaluation<n_components,
dim,
Expand Down Expand Up @@ -438,7 +453,10 @@ namespace VectorTools
dim,
spacedim,
typename VectorType::value_type>>(
cache.get_mapping(), cell->get_fe(), update_flags);
cache.get_mapping(),
cell->get_fe(),
update_flags,
first_selected_component);
auto &evaluator = *evaluators[cell->active_fe_index()];

evaluator.reinit(cell, unit_points);
Expand Down Expand Up @@ -532,6 +550,7 @@ namespace VectorTools
const VectorType & vector,
const UpdateFlags,
const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
const unsigned int first_selected_component,
const std::function<
value_type(const FEPointEvaluation<n_components,
dim,
Expand All @@ -547,8 +566,9 @@ namespace VectorTools
typename VectorType::value_type>>> &)
{
(void)evaluation_flags;
(void)first_selected_component;

Assert(n_components == 1,
Assert(n_components == 1 && first_selected_component == 0,
ExcMessage(
"A cell-data vector can only have a single component."));

Expand Down Expand Up @@ -581,7 +601,8 @@ namespace VectorTools
const MeshType & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags,
const UpdateFlags update_flags,
const unsigned int first_selected_component,
const UpdateFlags update_flags,
const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
const std::function<
value_type(const FEPointEvaluation<n_components,
Expand Down Expand Up @@ -627,6 +648,7 @@ namespace VectorTools
vector,
update_flags,
evaluation_flags,
first_selected_component,
process_quadrature_point,
values,
solution_values,
Expand Down Expand Up @@ -688,7 +710,8 @@ namespace VectorTools
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags)
const EvaluationFlags::EvaluationFlags flags,
const unsigned int first_selected_component)
{
return internal::evaluate_at_points<
n_components,
Expand All @@ -704,6 +727,7 @@ namespace VectorTools
mesh,
vector,
flags,
first_selected_component,
update_values,
dealii::EvaluationFlags::values,
[](const auto &evaluator, const auto &q) {
Expand All @@ -726,7 +750,8 @@ namespace VectorTools
const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
const MeshType<dim, spacedim> & mesh,
const VectorType & vector,
const EvaluationFlags::EvaluationFlags flags)
const EvaluationFlags::EvaluationFlags flags,
const unsigned int first_selected_component)
{
return internal::evaluate_at_points<
n_components,
Expand All @@ -743,6 +768,7 @@ namespace VectorTools
mesh,
vector,
flags,
first_selected_component,
update_gradients,
dealii::EvaluationFlags::gradients,
[](const auto &evaluator, const unsigned &q) {
Expand Down
62 changes: 36 additions & 26 deletions source/numerics/data_out_resample.cc
Original file line number Diff line number Diff line change
Expand Up @@ -151,32 +151,42 @@ DataOutResample<dim, patch_dim, spacedim>::build_patches(

const auto &dh = *data_ptr->dof_handler;

// TODO: enable more components
AssertDimension(dh.get_fe_collection().n_components(), 1);

const auto values =
VectorTools::point_values<1>(rpe, dh, data_ptr->vector);

vectors.emplace_back(
std::make_shared<LinearAlgebra::distributed::Vector<double>>(
partitioner));

for (unsigned int j = 0; j < values.size(); ++j)
vectors.back()->local_element(point_to_local_vector_indices[j]) =
values[j];

vectors.back()->set_ghost_state(true);

// we can give the vectors arbitrary names ("temp_*") here, since these
// are only used internally (by patch_data_out) but not later on during
// the actual output to file
patch_data_out.add_data_vector(
*vectors.back(),
std::string("temp_" + std::to_string(counter)),
DataOut_DoFData<patch_dim, patch_dim, spacedim, spacedim>::
DataVectorType::type_dof_data);

counter++;
#ifdef DEBUG
for (const auto &fe : dh.get_fe_collection())
Assert(
fe.n_base_elements() == 1,
ExcMessage(
"This class currently only supports scalar elements and elements "
"with a single base element."));
#endif

for (unsigned int comp = 0; comp < dh.get_fe_collection().n_components();
++comp)
{
const auto values = VectorTools::point_values<1>(
rpe, dh, data_ptr->vector, VectorTools::EvaluationFlags::avg, comp);

vectors.emplace_back(
std::make_shared<LinearAlgebra::distributed::Vector<double>>(
partitioner));

for (unsigned int j = 0; j < values.size(); ++j)
vectors.back()->local_element(point_to_local_vector_indices[j]) =
values[j];

vectors.back()->set_ghost_state(true);

// we can give the vectors arbitrary names ("temp_*") here, since
// these are only used internally (by patch_data_out) but not later on
// during the actual output to file
patch_data_out.add_data_vector(
*vectors.back(),
std::string("temp_" + std::to_string(counter)),
DataOut_DoFData<patch_dim, patch_dim, spacedim, spacedim>::
DataVectorType::type_dof_data);

counter++;
}
}

patch_data_out.build_patches(*patch_mapping,
Expand Down
124 changes: 124 additions & 0 deletions tests/remote_point_evaluation/data_out_resample_02.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2021 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

// Test DataOutResample to create a slice for FESystem with multiple
// components.

#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi_remote_point_evaluation.h>

#include <deal.II/distributed/tria.h>

#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q1.h>

#include <deal.II/grid/grid_generator.h>

#include <deal.II/numerics/data_out_resample.h>
#include <deal.II/numerics/vector_tools.h>

#include "../tests.h"

using namespace dealii;

template <int dim>
class AnalyticalFunction : public Function<dim>
{
public:
AnalyticalFunction(const unsigned n_components)
: Function<dim>(n_components)
{}

virtual double
value(const Point<dim> &p, const unsigned int component = 0) const
{
return p[component] * (component + 1.0) + component;
}
};


template <int dim, int spacedim>
std::shared_ptr<const Utilities::MPI::Partitioner>
create_partitioner(const DoFHandler<dim, spacedim> &dof_handler)
{
IndexSet locally_relevant_dofs;

DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);

return std::make_shared<const Utilities::MPI::Partitioner>(
dof_handler.locally_owned_dofs(),
locally_relevant_dofs,
dof_handler.get_communicator());
}



int
main(int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
MPILogInitAll all;

const int dim = 3;
const int patch_dim = 2;
const int spacedim = 3;
const int n_components = dim;

const unsigned int n_refinements_1 = 3;
const unsigned int n_refinements_2 = 3;
const MPI_Comm comm = MPI_COMM_WORLD;

parallel::distributed::Triangulation<patch_dim, spacedim> tria_slice(comm);
GridGenerator::hyper_cube(tria_slice, -1.0, +1.0);
tria_slice.refine_global(n_refinements_2);

MappingQ1<patch_dim, spacedim> mapping_slice;

parallel::distributed::Triangulation<dim, spacedim> tria_backround(comm);
GridGenerator::hyper_cube(tria_backround, -1.0, +1.0);
tria_backround.refine_global(n_refinements_1);

DoFHandler<dim, spacedim> dof_handler(tria_backround);
dof_handler.distribute_dofs(
FESystem<dim, spacedim>(FE_Q<dim, spacedim>{1}, n_components));

MappingQ1<dim, spacedim> mapping;

LinearAlgebra::distributed::Vector<double> vector(
create_partitioner(dof_handler));

VectorTools::interpolate(mapping,
dof_handler,
AnalyticalFunction<dim>(n_components),
vector);

vector.update_ghost_values();

DataOutResample<dim, patch_dim, spacedim> data_out(tria_slice, mapping_slice);
data_out.add_data_vector(dof_handler, vector, "solution_0");
data_out.add_data_vector(dof_handler, vector, "solution_1");
data_out.update_mapping(mapping);
data_out.build_patches();

#if 1
data_out.write_vtk(deallog.get_file_stream());
#else
data_out.write_vtu_with_pvtu_record("./", "data_out_01", 0, comm, 1, 1);
#endif
}

0 comments on commit d13fb96

Please sign in to comment.