Program Listing for File laplacian.h¶
↰ Return to documentation for file (include/immersx/physics/laplacian.h)
// ---------------------------------------------------------------------
//
// 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.
//
// ---------------------------------------------------------------------
/* ---------------------------------------------------------------------
*/
#ifndef dealii_distributed_lagrange_multiplier_h
#define dealii_distributed_lagrange_multiplier_h
#include <deal.II/base/function.h>
#include <deal.II/base/parsed_convergence_table.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/lac/block_linear_operator.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/linear_operator_tools.h>
#include <immersx/compatibility/reference_cell.h>
#include <immersx/core/integration_support.h>
#include <immersx/physics/detail/robin_parameter_overrides.h>
// #define MATRIX_FREE_PATH
#define FORCE_USE_OF_TRILINOS
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/parsed_function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_fe_field.h>
#include <deal.II/fe/mapping_q.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/grid/manifold_lib.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/vector.h>
#ifdef MATRIX_FREE_PATH
# include <deal.II/matrix_free/operators.h>
# include <deal.II/multigrid/mg_coarse.h>
# include <deal.II/multigrid/mg_matrix.h>
# include <deal.II/multigrid/mg_smoother.h>
# include <deal.II/multigrid/mg_tools.h>
# include <deal.II/multigrid/mg_transfer_matrix_free.h>
# include <deal.II/multigrid/multigrid.h>
#endif
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/opencascade/manifold_lib.h>
#include <deal.II/opencascade/utilities.h>
#include <immersx/algebra/matrix_free_utils.h>
#include <immersx/core/domain.h>
#include <immersx/coupling/inclusions.h>
#ifdef DEAL_II_WITH_OPENCASCADE
# include <TopoDS.hxx>
#endif
#include <cmath>
#include <fstream>
#include <iostream>
#include <memory>
namespace ImmersX
{
namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
!(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
using namespace dealii::LinearAlgebraPETSc;
# define USE_PETSC_LA
#elif defined(DEAL_II_WITH_TRILINOS)
using namespace dealii::LinearAlgebraTrilinos;
#else
# error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
#endif
} // namespace LA
template <int dim, int spacedim = dim>
class ProblemParameters : public ParameterAcceptor
{
public:
ProblemParameters();
std::string output_directory = ".";
std::string output_name = "solution";
unsigned int fe_degree = 1;
DomainParameters<spacedim> domain_parameters;
std::list<types::boundary_id> dirichlet_ids{0};
std::string refinement_strategy =
"fixed_fraction";
double coarsening_fraction = 0.0;
double refinement_fraction = 0.3;
unsigned int n_refinement_cycles = 1;
unsigned int max_cells = 20000;
mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
rhs;
mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
bc;
std::list<types::boundary_id> robin_ids;
mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
robin_coefficient;
mutable ModulatedParsedFunction<spacedim> robin_bc;
std::map<types::boundary_id,
std::shared_ptr<
ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>>>
robin_coefficients_by_id;
std::map<types::boundary_id,
std::shared_ptr<ModulatedParsedFunction<spacedim>>>
robin_bcs_by_id;
std::map<types::boundary_id, std::array<bool, 3>> robin_explicit_entries;
const Functions::ParsedFunction<spacedim> &
robin_coefficient_for(const types::boundary_id id) const
{
const auto entry = robin_coefficients_by_id.find(id);
return entry == robin_coefficients_by_id.end() ? robin_coefficient :
*entry->second;
}
const ModulatedParsedFunction<spacedim> &
robin_bc_for(const types::boundary_id id) const
{
const auto entry = robin_bcs_by_id.find(id);
return entry == robin_bcs_by_id.end() ? robin_bc : *entry->second;
}
mutable ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
exact_solution;
mutable ParameterAcceptorProxy<ReductionControl> inner_control;
mutable ParameterAcceptorProxy<ReductionControl>
outer_control;
bool output_results_before_solving = false;
bool estimate_condition_number = false;
mutable ParsedConvergenceTable convergence_table;
};
template <int dim, int spacedim>
ProblemParameters<dim, spacedim>::ProblemParameters()
: ParameterAcceptor("/Immersed Problem/")
, domain_parameters("/Immersed Problem/Domain/")
, rhs("/Immersed Problem/Right hand side")
, bc("/Immersed Problem/Dirichlet boundary conditions")
, robin_coefficient("/Immersed Problem/Functions/Robin coefficient", 1)
, robin_bc("/Immersed Problem/Functions/Robin boundary conditions", 1)
, exact_solution("/Immersed Problem/Exact solution")
, inner_control("/Immersed Problem/Solver/Inner control")
, outer_control("/Immersed Problem/Solver/Outer control")
{
domain_parameters.initial_refinement = 5;
add_parameter("FE degree", fe_degree, "", this->prm, Patterns::Integer(1));
add_parameter("Output directory", output_directory);
add_parameter("Output name", output_name);
add_parameter("Output results also before solving",
output_results_before_solving);
add_parameter("Estimate condition number", estimate_condition_number);
add_parameter("Dirichlet boundary ids", dirichlet_ids);
add_parameter("Robin boundary ids", robin_ids);
enter_subsection("Refinement and remeshing");
{
add_parameter("Strategy",
refinement_strategy,
"",
this->prm,
Patterns::Selection("fixed_fraction|fixed_number|global"));
add_parameter("Coarsening fraction", coarsening_fraction);
add_parameter("Refinement fraction", refinement_fraction);
add_parameter("Maximum number of cells", max_cells);
add_parameter("Number of refinement cycles", n_refinement_cycles);
}
leave_subsection();
this->prm.enter_subsection("Error");
convergence_table.add_parameters(this->prm);
this->prm.leave_subsection();
this->parse_parameters_call_back.connect([this]() {
robin_detail::register_parameter_overrides(*this,
robin_ids,
"/Immersed Problem/Functions/",
robin_bc,
robin_coefficients_by_id,
robin_bcs_by_id,
robin_explicit_entries);
});
}
template <int dim, int spacedim = dim>
class PoissonProblem : public EnableObserverPointer
{
public:
PoissonProblem(const ProblemParameters<dim, spacedim> &par);
double
multiplier_solution_l2_norm() const
{
return solution.block(1).l2_norm();
}
const DoFHandler<spacedim> &
bulk_dof_handler() const
{
return dh;
}
const auto &
bulk_values() const
{
return solution.block(0);
}
const auto &
bulk_rhs() const
{
return system_rhs.block(0);
}
const auto &
assembled_bulk() const
{
return stiffness_matrix;
}
void
make_grid();
void
setup_fe();
void
setup_dofs();
#ifndef MATRIX_FREE_PATH
void
assemble_poisson_system();
#else
void
assemble_rhs();
#endif
void
assemble_coupling();
void
run();
IndexSet
assemble_coupling_sparsity(DynamicSparsityPattern &dsp) const;
void
solve();
void
refine_and_transfer();
std::string
output_solution() const;
void
output_results() const;
void
print_parameters() const;
private:
const ProblemParameters<dim, spacedim> ∥
MPI_Comm mpi_communicator;
ConditionalOStream pcout;
mutable TimerOutput computing_timer;
std::unique_ptr<Domain<spacedim>> domain;
parallel::TriangulationBase<spacedim> *tria = nullptr;
std::unique_ptr<FiniteElement<spacedim>> fe;
Inclusions<spacedim> inclusions;
std::unique_ptr<Quadrature<spacedim>> quadrature;
DoFHandler<spacedim> dh;
std::vector<IndexSet> owned_dofs;
std::vector<IndexSet>
relevant_dofs;
AffineConstraints<double> constraints;
AffineConstraints<double>
inclusion_constraints;
LA::MPI::SparseMatrix multiplier_load_matrix;
LA::MPI::SparseMatrix outside_mass_matrix;
LA::MPI::SparseMatrix inclusion_matrix;
LA::MPI::SparseMatrix mass_matrix;
const Mapping<spacedim> &
geometry_mapping() const
{
AssertThrow(fe != nullptr, ExcNotInitialized());
return compatibility::linear_mapping<spacedim>(fe->reference_cell());
}
#ifdef MATRIX_FREE_PATH
using VectorType = LinearAlgebra::distributed::Vector<double>;
using BlockVectorType = LinearAlgebra::distributed::BlockVector<double>;
std::unique_ptr<MatrixFreeCouplingOperator<spacedim, double, 1>>
coupling_operator;
MatrixFreeOperators::LaplaceOperator<spacedim, -1> stiffness_matrix;
using LevelMatrixType = MatrixFreeOperators::LaplaceOperator<
spacedim,
-1,
-1,
1,
LinearAlgebra::distributed::Vector<float>>;
MGLevelObject<LevelMatrixType> mg_matrices;
MGConstrainedDoFs mg_constrained_dofs;
#else
LA::MPI::SparseMatrix stiffness_matrix;
using VectorType = LA::MPI::Vector;
using BlockVectorType = LA::MPI::BlockVector;
#endif
BlockVectorType solution;
BlockVectorType locally_relevant_solution;
BlockVectorType system_rhs;
std::vector<std::vector<BoundingBox<spacedim>>> global_bounding_boxes;
unsigned int cycle = 0;
};
} // namespace ImmersX
#endif