Program Listing for File tensor_product_space.cc¶
↰ Return to documentation for file (source/tensor_product_space.cc)
// ---------------------------------------------------------------------
//
// Copyright (C) 2024 by Luca Heltai
//
// This file is part of the ImmersX application, based on
// the deal.II library.
//
// The ImmersX application is free software; you can use
// it, redistribute it, and/or modify it under the terms of the Apache-2.0
// License WITH LLVM-exception as published by the Free Software Foundation;
// either version 3.0 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
// of the ImmersX distribution.
//
// ---------------------------------------------------------------------
#include "tensor_product_space.h"
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/quadrature_selector.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/particles/data_out.h>
#include <deal.II/particles/utilities.h>
#include "immersed_repartitioner.h"
#ifdef DEAL_II_WITH_VTK
# include "vtk_utils.h"
namespace
{
unsigned int
quadrature_selector_order(const std::string &quadrature_type,
const unsigned int requested_n_q_points,
const unsigned int fe_degree)
{
return quadrature_type == "gauss" ?
(requested_n_q_points == 0 ? 2 * fe_degree + 1 :
requested_n_q_points) :
0;
}
} // namespace
template <int reduced_dim, int dim, int spacedim, int n_components>
TensorProductSpaceParameters<reduced_dim, dim, spacedim, n_components>::
TensorProductSpaceParameters()
: ParameterAcceptor("Representative domain")
{
add_parameter("Finite element degree", fe_degree);
add_parameter("Quadrature type",
quadrature_type,
"1D quadrature family used by QuadratureSelector.",
this->prm,
Patterns::Selection(
QuadratureSelector<1>::get_quadrature_names()));
add_parameter("Number of quadrature points", n_q_points);
add_parameter("Number of quadrature repetitions",
n_quadrature_repetitions,
"How many times to repeat the quadrature formula.");
add_parameter("Thickness", thickness);
add_parameter("Thickness field name", thickness_field_name);
add_parameter("Reduced grid name", reduced_grid_name);
}
// Constructor for TensorProductSpace
template <int reduced_dim, int dim, int spacedim, int n_components>
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
TensorProductSpace(
const TensorProductSpaceParameters<reduced_dim, dim, spacedim, n_components>
&par,
MPI_Comm mpi_communicator)
: mpi_communicator(mpi_communicator)
, par(par)
, reference_cross_section(par.section)
, triangulation(mpi_communicator)
, fe(FE_Q<reduced_dim, spacedim>(par.fe_degree),
reference_cross_section.n_selected_basis())
, quadrature_formula(QIterated<reduced_dim>(
QuadratureSelector<1>(par.quadrature_type,
quadrature_selector_order(par.quadrature_type,
par.n_q_points,
par.fe_degree)),
par.n_quadrature_repetitions))
, dof_handler(triangulation)
, properties_dh(triangulation)
{}
// Initialize the tensor product space
template <int reduced_dim, int dim, int spacedim, int n_components>
void
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::initialize()
{
// Create the reduced grid and perform setup only if the triangulation is
// empty
if (triangulation.n_active_cells() == 0)
{
reference_cross_section.initialize();
make_reduced_grid_and_properties();
// Setup degrees of freedom
setup_dofs();
// Setup quadrature formulas
compute_points_and_weights();
}
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const ReferenceCrossSection<dim - reduced_dim, spacedim, n_components> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_reference_cross_section() const
{
return reference_cross_section;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
void
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
make_reduced_grid_and_properties()
{
// First create a serial triangulation with the VTK file
Triangulation<reduced_dim, spacedim> serial_tria;
DoFHandler<reduced_dim, spacedim> serial_properties_dh(serial_tria);
Vector<double> serial_properties;
VTKUtils::read_vtk(par.reduced_grid_name,
serial_properties_dh,
serial_properties,
properties_names);
deallog << "Read VTK file: " << par.reduced_grid_name
<< ", properties norm: " << serial_properties.l2_norm() << std::endl;
// Preprocess the serial triangulation
preprocess_serial_triangulation(serial_tria);
// Then make sure the partitioner is what the user wants
set_partitioner(triangulation);
// Once the triangulation is created, copy it to the distributed
// triangulation
triangulation.copy_triangulation(serial_tria);
if (triangulation.n_locally_owned_active_cells() == 0)
deallog << "Process " << Utilities::MPI::this_mpi_process(mpi_communicator)
<< " has no locally owned cells." << std::endl;
properties_dh.distribute_dofs(serial_properties_dh.get_fe());
properties.reinit(properties_dh.locally_owned_dofs(),
DoFTools::extract_locally_relevant_dofs(properties_dh),
mpi_communicator);
if (serial_properties_dh.n_dofs() == properties_dh.n_dofs())
{
VTKUtils::serial_vector_to_distributed_vector(serial_properties_dh,
properties_dh,
serial_properties,
properties);
}
else
{
// TODO: transfer from coarse serial grid to fine distributed grid
deallog << "PROPERTIES ARE ZERO. DO NOT REFINE INPUT GRID" << std::endl;
}
// Make sure we have ghost values
properties.update_ghost_values();
const auto &properties_fe = properties_dh.get_fe();
const auto block_indices = VTKUtils::get_block_indices(properties_fe);
for (unsigned int i = 0; i < block_indices.size(); ++i)
{
const auto &name = properties_names[i];
deallog << "Property name: " << name << ", block index: " << i
<< ", block size: " << block_indices.block_size(i)
<< ", block start: " << block_indices.block_start(i) << std::endl;
}
deallog << "Properties norm: " << properties.l2_norm() << std::endl;
deallog << "Serial properties norm: " << serial_properties.l2_norm()
<< std::endl;
AssertDimension(block_indices.total_size(), properties_fe.n_components());
AssertDimension(block_indices.size(), properties_names.size());
};
template <int reduced_dim, int dim, int spacedim, int n_components>
const DoFHandler<reduced_dim, spacedim> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::get_dof_handler()
const
{
return dof_handler;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const std::vector<Point<spacedim>> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_locally_owned_qpoints() const
{
const int n_local_qpoints = all_qpoints.size();
const int global_qpoints =
Utilities::MPI::sum(n_local_qpoints, mpi_communicator);
AssertThrow(
global_qpoints > 0,
ExcMessage(
"No quadrature points exist across all MPI ranks. You must call compute_points_and_weights() first"));
return all_qpoints;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const std::vector<std::vector<double>> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_locally_owned_weights() const
{
const int n_local_weights = all_weights.size();
const int global_weights =
Utilities::MPI::sum(n_local_weights, mpi_communicator);
AssertThrow(global_weights > 0,
ExcMessage("No weights exist across all MPI ranks. You must call"
" compute_points_and_weights() first"));
return all_weights;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const std::vector<Point<spacedim>> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_locally_owned_reduced_qpoints() const
{
const int n_local_reduced_qpoints = reduced_qpoints.size();
const int global_reduced_qpoints =
Utilities::MPI::sum(n_local_reduced_qpoints, mpi_communicator);
AssertThrow(
global_reduced_qpoints > 0,
ExcMessage(
"No reduced quadrature points exist across all MPI ranks. You must call"
" compute_points_and_weights() first"));
return reduced_qpoints;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const std::vector<std::vector<double>> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_locally_owned_reduced_weights() const
{
const int n_local_reduced_weights = reduced_weights.size();
const int global_reduced_weights =
Utilities::MPI::sum(n_local_reduced_weights, mpi_communicator);
AssertThrow(global_reduced_weights > 0,
ExcMessage(
"No reduced weights exist across all MPI ranks. You must call"
" compute_points_and_weights() first"));
return reduced_weights;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
void
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
update_local_dof_indices(
const std::map<unsigned int, IndexSet> &remote_q_point_indices)
{
auto global_cell_indices =
local_q_point_indices_to_global_cell_indices(remote_q_point_indices);
std::map<
unsigned int,
std::map<types::global_cell_index, std::vector<types::global_dof_index>>>
global_dof_indices;
for (const auto &[proc, cell_indices] : global_cell_indices)
for (const auto &id : cell_indices)
global_dof_indices[proc][id] = global_cell_to_dof_indices[id];
// Exchange the data with participating processors
auto local_dof_indices =
Utilities::MPI::some_to_some(mpi_communicator, global_dof_indices);
// update global_cell_to_dof_indices
for (const auto &[proc, cell_indices] : local_dof_indices)
{
global_cell_to_dof_indices.insert(cell_indices.begin(),
cell_indices.end());
}
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const std::vector<types::global_dof_index> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::get_dof_indices(
const types::global_cell_index cell_index) const
{
Assert(global_cell_to_dof_indices.find(cell_index) !=
global_cell_to_dof_indices.end(),
ExcMessage("Cell index not found in global cell to dof indices."));
return global_cell_to_dof_indices.at(cell_index);
}
template <int reduced_dim, int dim, int spacedim, int n_components>
std::tuple<unsigned int, unsigned int, unsigned int>
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
particle_id_to_cell_and_qpoint_indices(const unsigned int particle_id) const
{
AssertIndexRange(particle_id,
triangulation.n_global_active_cells() *
quadrature_formula.size() *
reference_cross_section.n_quadrature_points());
const unsigned int cell_index =
particle_id /
(quadrature_formula.size() * reference_cross_section.n_quadrature_points());
const unsigned int qpoint_index_in_cell =
(particle_id / reference_cross_section.n_quadrature_points()) %
quadrature_formula.size();
const unsigned int qpoint_index_in_section =
particle_id % reference_cross_section.n_quadrature_points();
AssertIndexRange(cell_index, triangulation.n_global_active_cells());
AssertIndexRange(qpoint_index_in_cell, quadrature_formula.size());
AssertIndexRange(qpoint_index_in_section,
reference_cross_section.n_quadrature_points());
return std::make_tuple(cell_index,
qpoint_index_in_cell,
qpoint_index_in_section);
}
template <int reduced_dim, int dim, int spacedim, int n_components>
std::map<unsigned int, IndexSet>
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
local_q_point_indices_to_global_cell_indices(
const std::map<unsigned int, IndexSet> &remote_q_point_indices) const
{
std::map<unsigned int, IndexSet> cell_indices;
const IndexSet &owned_cells =
triangulation.global_active_cell_index_partitioner()
.lock()
->locally_owned_range();
auto local_q_point_indices =
Utilities::MPI::some_to_some(mpi_communicator, remote_q_point_indices);
for (const auto &[proc, qpoint_indices] : local_q_point_indices)
{
IndexSet cell_indices_for_proc(triangulation.n_global_active_cells());
for (const auto &qpoint_index : qpoint_indices)
{
const auto [cell_index, q_index, i] =
particle_id_to_cell_and_qpoint_indices(qpoint_index);
cell_indices_for_proc.add_index(
owned_cells.nth_index_in_set(cell_index));
}
cell_indices_for_proc.compress();
cell_indices[proc] = cell_indices_for_proc;
}
return cell_indices;
}
// Setup degrees of freedom for the tensor product space
template <int reduced_dim, int dim, int spacedim, int n_components>
void
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::setup_dofs()
{
dof_handler.distribute_dofs(fe);
// Additional setup can be done here if needed
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
cell->get_dof_indices(dof_indices);
global_cell_to_dof_indices[cell->global_active_cell_index()] =
dof_indices;
}
}
template <int reduced_dim, int dim, int spacedim, int n_components>
IndexSet
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
locally_owned_qpoints() const
{
IndexSet locally_owned_cell_set =
triangulation.global_active_cell_index_partitioner()
.lock()
->locally_owned_range();
// Now make a tensor product of the local indices with the total number of
// quadrature points, and the number of quadrature points in the
// cross-section
const unsigned int n_qpoints_per_cell =
reference_cross_section.n_quadrature_points() * quadrature_formula.size();
IndexSet locally_owned_qpoints_set = locally_owned_cell_set.tensor_product(
complete_index_set(n_qpoints_per_cell));
return locally_owned_qpoints_set;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
IndexSet
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
locally_relevant_indices() const
{
IndexSet indices = triangulation.global_active_cell_index_partitioner()
.lock()
->locally_owned_range();
for (const auto &[cell_id, local_indices] : global_cell_to_dof_indices)
indices.add_index(cell_id);
indices.compress();
return indices;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
auto
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::get_quadrature()
const -> const Quadrature<reduced_dim> &
{
return quadrature_formula;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
void
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
compute_points_and_weights()
{
all_qpoints.reserve(triangulation.n_active_cells() *
quadrature_formula.size() *
reference_cross_section.n_quadrature_points());
all_weights.reserve(triangulation.n_active_cells() *
quadrature_formula.size() *
reference_cross_section.n_quadrature_points());
reduced_qpoints.reserve(triangulation.n_active_cells() *
quadrature_formula.size());
reduced_weights.reserve(triangulation.n_active_cells() *
quadrature_formula.size());
UpdateFlags flags =
reduced_dim == 1 ?
update_quadrature_points | update_JxW_values :
update_quadrature_points | update_normal_vectors | update_JxW_values;
FEValues<reduced_dim, spacedim> fev(fe, quadrature_formula, flags);
const auto &properties_fe = properties_dh.get_fe();
FEValues<reduced_dim, spacedim> properties_fe_values(properties_fe,
get_quadrature(),
update_values);
// Find the index of the thickness field in the properties
const unsigned int thickness_field_index =
std::distance(properties_names.begin(),
std::find(properties_names.begin(),
properties_names.end(),
par.thickness_field_name));
const auto thickness_start =
thickness_field_index >= properties_names.size() ?
numbers::invalid_unsigned_int :
VTKUtils::get_block_indices(properties_fe)
.block_start(thickness_field_index);
FEValuesExtractors::Scalar thickness(thickness_start);
// Initialize the thickness values with the constant thickness
std::vector<double> thickness_values(get_quadrature().size(), par.thickness);
for (const auto &cell : triangulation.active_cell_iterators())
if (cell->is_locally_owned())
{
fev.reinit(cell);
const auto &qpoints = fev.get_quadrature_points();
const auto &JxW = fev.get_JxW_values();
if (thickness_start != numbers::invalid_unsigned_int)
{
properties_fe_values.reinit(
cell->as_dof_handler_iterator(this->properties_dh));
properties_fe_values[thickness].get_function_values(
properties, thickness_values);
}
reduced_qpoints.insert(reduced_qpoints.end(),
qpoints.begin(),
qpoints.end());
for (const auto &w : JxW)
reduced_weights.emplace_back(std::vector<double>(1, w));
Tensor<1, spacedim> new_vertical;
if constexpr (reduced_dim == 1)
new_vertical = cell->vertex(1) - cell->vertex(0);
for (const auto &q : fev.quadrature_point_indices())
{
const auto &qpoint = qpoints[q];
if constexpr (reduced_dim == 2)
new_vertical = fev.normal_vector(q);
// [TODO] Make radius a function of the cell
auto cross_section_qpoints =
reference_cross_section.get_transformed_quadrature(
qpoint, new_vertical, thickness_values[q]);
all_qpoints.insert(all_qpoints.end(),
cross_section_qpoints.get_points().begin(),
cross_section_qpoints.get_points().end());
for (const auto &w : cross_section_qpoints.get_weights())
all_weights.emplace_back(std::vector<double>(1, w * fev.JxW(q)));
}
}
};
template <int reduced_dim, int dim, int spacedim, int n_components>
const parallel::fullydistributed::Triangulation<reduced_dim, spacedim> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_triangulation() const
{
return triangulation;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
double
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::get_scaling(
const unsigned int) const
{
return std::pow(par.thickness, -((dim - reduced_dim) / 2.0));
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const LinearAlgebra::distributed::Vector<double> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::get_properties()
const
{
return properties;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const DoFHandler<reduced_dim, spacedim> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_properties_dh() const
{
return properties_dh;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
DoFHandler<reduced_dim, spacedim> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_properties_dh()
{
return properties_dh;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
const std::vector<std::string> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_properties_names() const
{
return properties_names;
}
template <int reduced_dim, int dim, int spacedim, int n_components>
std::vector<std::string> &
TensorProductSpace<reduced_dim, dim, spacedim, n_components>::
get_properties_names()
{
return properties_names;
}
template struct TensorProductSpaceParameters<1, 2, 2, 1>;
template struct TensorProductSpaceParameters<1, 2, 3, 1>;
template struct TensorProductSpaceParameters<1, 3, 3, 1>;
template struct TensorProductSpaceParameters<2, 3, 3, 1>;
template struct TensorProductSpaceParameters<1, 2, 2, 2>;
template struct TensorProductSpaceParameters<1, 2, 3, 3>;
template struct TensorProductSpaceParameters<1, 3, 3, 3>;
template struct TensorProductSpaceParameters<2, 3, 3, 3>;
template class TensorProductSpace<1, 2, 2, 1>;
template class TensorProductSpace<1, 2, 3, 1>;
template class TensorProductSpace<1, 3, 3, 1>;
template class TensorProductSpace<2, 3, 3, 1>;
template class TensorProductSpace<1, 2, 2, 2>;
template class TensorProductSpace<1, 2, 3, 3>;
template class TensorProductSpace<1, 3, 3, 3>;
template class TensorProductSpace<2, 3, 3, 3>;
#endif // DEAL_II_WITH_VTK