Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Vector tools project 1 #14101

Merged
merged 1 commit into from Aug 2, 2022
Merged

Vector tools project 1 #14101

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
132 changes: 64 additions & 68 deletions include/deal.II/numerics/vector_tools_project.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,11 +52,11 @@ namespace VectorTools
* mapping here because the function we evaluate for the DoFs is zero in
* the mapped locations as well as in the original, unmapped locations
*/
template <int dim, int spacedim, typename number>
template <int dim, int spacedim, typename Number>
void
interpolate_zero_boundary_values(
const DoFHandler<dim, spacedim> & dof_handler,
std::map<types::global_dof_index, number> &boundary_values)
std::map<types::global_dof_index, Number> &boundary_values)
{
// loop over all boundary faces
// to get all dof indices of
Expand All @@ -79,11 +79,8 @@ namespace VectorTools
// that is actually wholly on
// the boundary, not only by
// one line or one vertex
typename DoFHandler<dim, spacedim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
std::vector<types::global_dof_index> face_dof_indices;
for (; cell != endc; ++cell)
for (const auto &cell : dof_handler.active_cell_iterators())
for (auto f : GeometryInfo<dim>::face_indices())
if (cell->at_boundary(f))
{
Expand All @@ -102,6 +99,8 @@ namespace VectorTools
}
}



/**
* Compute the boundary values to be used in the project() functions.
*/
Expand All @@ -111,16 +110,16 @@ namespace VectorTools
class M_or_MC,
template <int>
class Q_or_QC,
typename number>
typename Number>
void
project_compute_b_v(
const M_or_MC<dim, spacedim> & mapping,
const DoFHandler<dim, spacedim> & dof,
const Function<spacedim, number> & function,
const Function<spacedim, Number> & function,
const bool enforce_zero_boundary,
const Q_or_QC<dim - 1> & q_boundary,
const bool project_to_boundary_first,
std::map<types::global_dof_index, number> &boundary_values)
std::map<types::global_dof_index, Number> &boundary_values)
{
if (enforce_zero_boundary == true)
// no need to project boundary
Expand All @@ -143,7 +142,7 @@ namespace VectorTools
const std::vector<types::boundary_id> used_boundary_ids =
dof.get_triangulation().get_boundary_ids();

std::map<types::boundary_id, const Function<spacedim, number> *>
std::map<types::boundary_id, const Function<spacedim, Number> *>
boundary_functions;
for (const auto used_boundary_id : used_boundary_ids)
boundary_functions[used_boundary_id] = &function;
Expand All @@ -153,6 +152,7 @@ namespace VectorTools
}



/*
* MatrixFree implementation of project() for an arbitrary number of
* components of the FiniteElement.
Expand Down Expand Up @@ -180,12 +180,8 @@ namespace VectorTools
(void)q_boundary;

AssertDimension(dof.get_fe_collection().size(), 1);

Assert(dof.get_fe(0).n_components() == function.n_components,
ExcDimensionMismatch(dof.get_fe(0).n_components(),
function.n_components));
Assert(dof.get_fe(0).n_components() == components,
ExcDimensionMismatch(components, dof.get_fe(0).n_components()));
AssertDimension(dof.get_fe(0).n_components(), function.n_components);
AssertDimension(dof.get_fe(0).n_components(), components);

Quadrature<dim> quadrature_mf;

Expand Down Expand Up @@ -329,8 +325,10 @@ namespace VectorTools



// Helper interface for the matrix-free implementation of project().
// Used to determine the number of components.
/**
* Helper interface for the matrix-free implementation of project(). Used
* to determine the number of components.
*/
template <int dim, typename Number, int spacedim>
void
project_matrix_free_component(
Expand Down Expand Up @@ -422,8 +420,7 @@ namespace VectorTools
const Quadrature<dim - 1> &q_boundary,
const bool project_to_boundary_first)
{
Assert(vec_result.size() == dof.n_dofs(),
ExcDimensionMismatch(vec_result.size(), dof.n_dofs()));
AssertDimension(vec_result.size(), dof.n_dofs());

LinearAlgebra::distributed::Vector<typename VectorType::value_type>
work_result;
Expand All @@ -437,11 +434,9 @@ namespace VectorTools
q_boundary,
project_to_boundary_first);

const IndexSet & locally_owned_dofs = dof.locally_owned_dofs();
IndexSet::ElementIterator it = locally_owned_dofs.begin();
for (; it != locally_owned_dofs.end(); ++it)
::dealii::internal::ElementAccess<VectorType>::set(work_result(*it),
*it,
for (const auto i : dof.locally_owned_dofs())
::dealii::internal::ElementAccess<VectorType>::set(work_result(i),
i,
vec_result);
vec_result.compress(VectorOperation::insert);
}
Expand All @@ -452,11 +447,11 @@ namespace VectorTools
* Return whether the boundary values try to constrain a degree of freedom
* that is already constrained to something else
*/
template <typename number>
template <typename Number>
bool
constraints_and_b_v_are_compatible(
const AffineConstraints<number> & constraints,
std::map<types::global_dof_index, number> &boundary_values)
const AffineConstraints<Number> & constraints,
std::map<types::global_dof_index, Number> &boundary_values)
{
for (const auto &boundary_value : boundary_values)
if (constraints.is_constrained(boundary_value.first))
Expand Down Expand Up @@ -495,15 +490,12 @@ namespace VectorTools
const Q_or_QC<dim - 1> &q_boundary,
const bool project_to_boundary_first)
{
using number = typename VectorType::value_type;
Assert(dof.get_fe(0).n_components() == function.n_components,
ExcDimensionMismatch(dof.get_fe(0).n_components(),
function.n_components));
Assert(vec_result.size() == dof.n_dofs(),
ExcDimensionMismatch(vec_result.size(), dof.n_dofs()));
using Number = typename VectorType::value_type;
AssertDimension(dof.get_fe(0).n_components(), function.n_components);
AssertDimension(vec_result.size(), dof.n_dofs());

// make up boundary values
std::map<types::global_dof_index, number> boundary_values;
std::map<types::global_dof_index, Number> boundary_values;
project_compute_b_v(mapping,
dof,
function,
Expand All @@ -514,11 +506,11 @@ namespace VectorTools

// check if constraints are compatible (see below)
const bool constraints_are_compatible =
constraints_and_b_v_are_compatible<number>(constraints,
constraints_and_b_v_are_compatible<Number>(constraints,
boundary_values);

// set up mass matrix and right hand side
Vector<number> vec(dof.n_dofs());
Vector<Number> vec(dof.n_dofs());
SparsityPattern sparsity;
{
DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
Expand All @@ -529,8 +521,8 @@ namespace VectorTools

sparsity.copy_from(dsp);
}
SparseMatrix<number> mass_matrix(sparsity);
Vector<number> tmp(mass_matrix.n());
SparseMatrix<Number> mass_matrix(sparsity);
Vector<Number> tmp(mass_matrix.n());

// If the constraints object does not conflict with the given boundary
// values (i.e., it either does not contain boundary values or it contains
Expand All @@ -539,7 +531,7 @@ namespace VectorTools
// interpolate the boundary values and then condense the matrix and vector
if (constraints_are_compatible)
{
const Function<spacedim, number> *dummy = nullptr;
const Function<spacedim, Number> *dummy = nullptr;
MatrixCreator::create_mass_matrix(mapping,
dof,
quadrature,
Expand All @@ -566,10 +558,10 @@ namespace VectorTools
// steps may not be sufficient, since roundoff errors may accumulate for
// badly conditioned matrices
ReductionControl control(5 * tmp.size(), 0., 1e-12, false, false);
GrowingVectorMemory<Vector<number>> memory;
SolverCG<Vector<number>> cg(control, memory);
GrowingVectorMemory<Vector<Number>> memory;
SolverCG<Vector<Number>> cg(control, memory);

PreconditionSSOR<SparseMatrix<number>> prec;
PreconditionSSOR<SparseMatrix<Number>> prec;
prec.initialize(mass_matrix, 1.2);

cg.solve(mass_matrix, vec, tmp, prec);
Expand All @@ -584,6 +576,8 @@ namespace VectorTools
vec_result);
}



template <int dim, typename VectorType, int spacedim>
void
project_parallel(
Expand All @@ -596,11 +590,10 @@ namespace VectorTools
const unsigned int)> & func,
VectorType & vec_result)
{
using Number = typename VectorType::value_type;
Assert(dof.get_fe(0).n_components() == 1,
ExcDimensionMismatch(dof.get_fe(0).n_components(), 1));
Assert(vec_result.size() == dof.n_dofs(),
ExcDimensionMismatch(vec_result.size(), dof.n_dofs()));
using Number = typename VectorType::value_type;
using LocalVectorType = LinearAlgebra::distributed::Vector<Number>;
AssertDimension(dof.get_fe(0).n_components(), 1);
AssertDimension(vec_result.size(), dof.n_dofs());

// set up mass matrix and right hand side
typename MatrixFree<dim, Number>::AdditionalData additional_data;
Expand All @@ -615,13 +608,12 @@ namespace VectorTools
constraints,
QGauss<1>(dof.get_fe().degree + 2),
additional_data);
using MatrixType = MatrixFreeOperators::
MassOperator<dim, -1, 0, 1, LinearAlgebra::distributed::Vector<Number>>;
using MatrixType =
MatrixFreeOperators::MassOperator<dim, -1, 0, 1, LocalVectorType>;
MatrixType mass_matrix;
mass_matrix.initialize(matrix_free);
mass_matrix.compute_diagonal();

using LocalVectorType = LinearAlgebra::distributed::Vector<Number>;
LocalVectorType vec, rhs, inhomogeneities;
matrix_free->initialize_dof_vector(vec);
matrix_free->initialize_dof_vector(rhs);
Expand All @@ -640,10 +632,7 @@ namespace VectorTools
const unsigned int n_q_points = quadrature.size();
Vector<Number> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
typename DoFHandler<dim, spacedim>::active_cell_iterator
cell = dof.begin_active(),
endc = dof.end();
for (; cell != endc; ++cell)
for (const auto &cell : dof.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_rhs = 0;
Expand Down Expand Up @@ -672,7 +661,7 @@ namespace VectorTools
// badly conditioned matrices. This behavior can be observed, e.g. for
// FE_Q_Hierarchical for degree higher than three.
ReductionControl control(5 * rhs.size(), 0., 1e-12, false, false);
SolverCG<LinearAlgebra::distributed::Vector<Number>> cg(control);
SolverCG<LocalVectorType> cg(control);
typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
PreconditionJacobi<MatrixType> preconditioner;
preconditioner.initialize(mass_matrix, data);
Expand Down Expand Up @@ -707,19 +696,17 @@ namespace VectorTools
const DoFHandler<dim, spacedim> &dof =
matrix_free->get_dof_handler(fe_component);

using Number = typename VectorType::value_type;
Assert(dof.get_fe(0).n_components() == 1,
ExcDimensionMismatch(dof.get_fe(0).n_components(), 1));
Assert(vec_result.size() == dof.n_dofs(),
ExcDimensionMismatch(vec_result.size(), dof.n_dofs()));
using Number = typename VectorType::value_type;
using LocalVectorType = LinearAlgebra::distributed::Vector<Number>;
AssertDimension(dof.get_fe(0).n_components(), 1);
AssertDimension(vec_result.size(), dof.n_dofs());

using MatrixType = MatrixFreeOperators::
MassOperator<dim, -1, 0, 1, LinearAlgebra::distributed::Vector<Number>>;
using MatrixType =
MatrixFreeOperators::MassOperator<dim, -1, 0, 1, LocalVectorType>;
MatrixType mass_matrix;
mass_matrix.initialize(matrix_free, {fe_component});
mass_matrix.compute_diagonal();

using LocalVectorType = LinearAlgebra::distributed::Vector<Number>;
LocalVectorType vec, rhs, inhomogeneities;
matrix_free->initialize_dof_vector(vec, fe_component);
matrix_free->initialize_dof_vector(rhs, fe_component);
Expand Down Expand Up @@ -753,7 +740,7 @@ namespace VectorTools
// badly conditioned matrices. This behavior can be observed, e.g. for
// FE_Q_Hierarchical for degree higher than three.
ReductionControl control(5 * rhs.size(), 0., 1e-12, false, false);
SolverCG<LinearAlgebra::distributed::Vector<Number>> cg(control);
SolverCG<LocalVectorType> cg(control);
typename PreconditionJacobi<MatrixType>::AdditionalData data(0.8);
PreconditionJacobi<MatrixType> preconditioner;
preconditioner.initialize(mass_matrix, data);
Expand All @@ -771,11 +758,13 @@ namespace VectorTools
vec_result.compress(VectorOperation::insert);
}



/**
* Specialization of project() for the case dim==spacedim and with correct
* number types for MatrixFree support. Check if we actually can use the
* MatrixFree implementation or need to use the matrix based one
* nonetheless based on the number of components.
* number types for MatrixFree support. Check if we actually can use the
* MatrixFree implementation or need to use the matrix based one nonetheless
* based on the number of components.
*/
template <typename VectorType, int dim, int spacedim>
std::enable_if_t<
Expand Down Expand Up @@ -834,6 +823,8 @@ namespace VectorTools
}
}



/**
* Specialization of project() for complex numbers or `dim < spacedim`,
* for which we are sure that we cannot use the MatrixFree implementation.
Expand Down Expand Up @@ -867,8 +858,11 @@ namespace VectorTools
q_boundary,
project_to_boundary_first);
}

} // namespace internal



template <int dim, typename VectorType, int spacedim>
void
project(const Mapping<dim, spacedim> & mapping,
Expand Down Expand Up @@ -1006,6 +1000,7 @@ namespace VectorTools
}



template <int dim, typename VectorType, int spacedim>
void
project(const DoFHandler<dim, spacedim> & dof,
Expand All @@ -1027,6 +1022,7 @@ namespace VectorTools
q_boundary,
project_to_boundary_first);
}

} // namespace VectorTools

DEAL_II_NAMESPACE_CLOSE
Expand Down