245 *
#include <deal.II/base/function_parser.h>
246 *
#include <deal.II/base/index_set.h>
247 *
#include <deal.II/base/parameter_handler.h>
248 *
#include <deal.II/base/quadrature_lib.h>
249 *
#include <deal.II/base/utilities.h>
251 *
#include <deal.II/dofs/dof_handler.h>
252 *
#include <deal.II/dofs/dof_tools.h>
254 *
#include <deal.II/fe/fe_nedelec.h>
255 *
#include <deal.II/fe/fe_series.h>
256 *
#include <deal.II/fe/fe_values.h>
259 *
#include <deal.II/grid/cell_data.h>
261 *
#include <deal.II/grid/tria.h>
262 *
#include <deal.II/grid/tria_iterator.h>
264 *
#include <deal.II/lac/affine_constraints.h>
265 *
#include <deal.II/lac/full_matrix.h>
266 *
#include <deal.II/lac/petsc_precondition.h>
267 *
#include <deal.II/lac/petsc_sparse_matrix.h>
268 *
#include <deal.II/lac/petsc_vector.h>
269 *
#include <deal.II/lac/slepc_solver.h>
271 *
#include <deal.II/numerics/data_out.h>
272 *
#include <deal.II/numerics/vector_tools.h>
276 * For parallelization (
using WorkStream and Intel TBB)
279 *
#include <deal.II/base/multithread_info.h>
280 *
#include <deal.II/base/work_stream.h>
282 *
#include
"petscpc.h"
286 * For Error Estimation/Indication and Smoothness Indication
289 *
#include <deal.II/fe/fe_tools.h>
291 *
#include <deal.II/numerics/error_estimator.h>
292 *
#include <deal.II/numerics/smoothness_estimator.h>
298 *
#include <deal.II/grid/grid_refinement.h>
301 *
#include <iostream>
304 *
namespace Operations
311 *
curlcurl(const ::FEValues<2> &fe_values,
312 *
const unsigned int & i,
313 *
const unsigned int & j,
314 *
const unsigned int & q_point)
316 *
auto gradu1_x1x2 = fe_values.shape_grad_component(i, q_point, 0);
317 *
auto gradu2_x1x2 = fe_values.shape_grad_component(i, q_point, 1);
319 *
auto gradv1_x1x2 = fe_values.shape_grad_component(j, q_point, 0);
320 *
auto gradv2_x1x2 = fe_values.shape_grad_component(j, q_point, 1);
321 *
return (gradu2_x1x2[0] - gradu1_x1x2[1]) *
322 *
(gradv2_x1x2[0] - gradv1_x1x2[1]);
330 *
dot_term(const ::FEValues<dim> &fe_values,
331 *
const unsigned int & i,
332 *
const unsigned int & j,
333 *
const unsigned int & q_point)
335 *
double output = 0.0;
336 *
for (
unsigned int comp = 0; comp < dim; ++comp)
338 *
output += fe_values.shape_value_component(i, q_point, comp) *
339 *
fe_values.shape_value_component(j, q_point, comp);
349 *
namespace Structures
354 *
create_L_waveguide(
Triangulation<2> &triangulation,
const double &scaling)
356 *
const unsigned int dim = 2;
358 *
const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
359 *
{scaling * 0.5, scaling * 0.0},
360 *
{scaling * 0.0, scaling * 0.5},
361 *
{scaling * 0.5, scaling * 0.5},
362 *
{scaling * 0.0, scaling * 1.0},
363 *
{scaling * 0.5, scaling * 1.0},
364 *
{scaling * 1.0, scaling * 0.5},
365 *
{scaling * 1.0, scaling * 1.0}};
367 *
const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
368 *
cell_vertices = {{{0, 1, 2, 3}}, {{2, 3, 4, 5}}, {{3, 6, 5, 7}}};
369 *
const unsigned int n_cells = cell_vertices.size();
370 *
std::vector<CellData<dim>> cells(n_cells,
CellData<dim>());
371 *
for (
unsigned int i = 0; i <
n_cells; ++i)
373 *
for (
unsigned int j = 0; j < cell_vertices[i].size(); ++j)
374 *
cells[i].vertices[j] = cell_vertices[i][j];
375 *
cells[i].material_id = 0;
377 *
triangulation.create_triangulation(vertices, cells,
SubCellData());
378 *
triangulation.refine_global(1);
384 *
const double & scaling)
386 *
const unsigned int dim = 2;
388 *
const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
389 *
{scaling * 0.6, scaling * 0.0},
390 *
{scaling * 0.0, scaling * 0.3},
391 *
{scaling * 0.6, scaling * 0.3}};
393 *
const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
394 *
cell_vertices = {{{0, 1, 2, 3}}};
395 *
const unsigned int n_cells = cell_vertices.size();
396 *
std::vector<CellData<dim>> cells(n_cells,
CellData<dim>());
397 *
for (
unsigned int i = 0; i <
n_cells; ++i)
399 *
for (
unsigned int j = 0; j < cell_vertices[i].size(); ++j)
400 *
cells[i].vertices[j] = cell_vertices[i][j];
401 *
cells[i].material_id = 0;
403 *
triangulation.create_triangulation(vertices, cells,
SubCellData());
404 *
triangulation.refine_global(0);
431 *
virtual unsigned int
432 *
solve_problem() = 0;
434 *
set_refinement_cycle(
const unsigned int cycle);
437 *
output_solution() = 0;
442 *
unsigned int refinement_cycle = 0;
443 *
std::unique_ptr<ParameterHandler> parameters;
444 *
unsigned int n_eigenpairs = 1;
445 *
double target = 0.0;
446 *
unsigned int eigenpair_selection_scheme;
447 *
unsigned int max_cycles = 0;
448 *
ompi_communicator_t * mpi_communicator = PETSC_COMM_SELF;
456 *
: triangulation(&coarse_grid)
457 *
, parameters(std::make_unique<ParameterHandler>())
459 *
parameters->declare_entry(
460 *
"Eigenpair selection scheme",
463 *
"The type of eigenpairs to find (0 - smallest, 1 - target)");
464 *
parameters->declare_entry(
"Number of eigenvalues/eigenfunctions",
467 *
"The number of eigenvalues/eigenfunctions "
468 *
"to be computed.");
469 *
parameters->declare_entry(
"Target eigenvalue",
472 *
"The target eigenvalue (if scheme == 1)");
474 *
parameters->declare_entry(
"Cycles number",
477 *
"The number of cycles in refinement");
478 *
parameters->parse_input(prm_file);
480 *
eigenpair_selection_scheme =
481 *
parameters->get_integer(
"Eigenpair selection scheme");
485 * The
project currently only supports selection by a target eigenvalue.
486 * Furthermore, only
one eigenpair can be computed at a time.
489 *
assert(eigenpair_selection_scheme == 1 &&
490 *
"Selection by a target is the only currently supported option!");
492 *
parameters->get_integer(
"Number of eigenvalues/eigenfunctions");
494 *
n_eigenpairs == 1 &&
495 *
"Only the computation of a single eigenpair is currently supported!");
497 *
target = parameters->get_double(
"Target eigenvalue");
498 *
max_cycles = parameters->get_integer(
"Cycles number");
499 *
if (eigenpair_selection_scheme == 1)
505 *
Base<dim>::set_refinement_cycle(
const unsigned int cycle)
507 *
refinement_cycle = cycle;
516 *
class EigenSolver :
public virtual Base<dim>
519 *
EigenSolver(
const std::string & prm_file,
520 *
Triangulation<dim> &coarse_grid,
521 *
const unsigned int &minimum_degree,
522 *
const unsigned int &maximum_degree,
523 *
const unsigned int &starting_degree);
525 *
virtual unsigned int
526 *
solve_problem()
override;
528 *
virtual unsigned int
531 *
template <
class SolverType>
533 *
initialize_eigensolver(SolverType &eigensolver);
542 *
const std::unique_ptr<hp::FECollection<dim>> fe_collection;
543 *
std::unique_ptr<hp::QCollection<dim>> quadrature_collection;
544 *
std::unique_ptr<hp::QCollection<dim - 1>> face_quadrature_collection;
545 *
DoFHandler<dim> dof_handler;
546 *
const unsigned int max_degree, min_degree;
549 *
for the actual solution
552 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> eigenfunctions;
553 *
std::unique_ptr<std::vector<double>>
eigenvalues;
554 *
Vector<double> solution;
563 *
convert_solution();
566 *
AffineConstraints<double> constraints;
567 *
PETScWrappers::SparseMatrix stiffness_matrix,
mass_matrix;
575 *
EigenSolver<dim>::EigenSolver(
const std::string & prm_file,
577 *
const unsigned int &minimum_degree,
578 *
const unsigned int &maximum_degree,
579 *
const unsigned int &starting_degree)
580 *
: Base<dim>(prm_file, triangulation)
584 *
, dof_handler(triangulation)
585 *
, max_degree(maximum_degree)
586 *
, min_degree(minimum_degree)
588 *
std::make_unique<std::vector<PETScWrappers::MPI::Vector>>())
589 *
,
eigenvalues(std::make_unique<std::vector<double>>())
591 *
for (
unsigned int degree = min_degree; degree <= max_degree; ++degree)
596 * Generate quadrature collection with
sorted quadrature weights
601 *
quadrature_collection->push_back(sorted_quadrature);
603 *
const QGauss<dim - 1> face_quadrature(degree + 1);
604 *
const QSorted<dim - 1> sorted_face_quadrature(face_quadrature);
605 *
face_quadrature_collection->push_back(sorted_face_quadrature);
609 * adjust the discretization
612 *
if (starting_degree > min_degree && starting_degree <= max_degree)
614 *
const unsigned int start_diff = starting_degree - min_degree;
616 *
cell1 = dof_handler.begin_active(),
617 *
endc1 = dof_handler.end();
618 *
for (; cell1 < endc1; ++cell1)
620 *
cell1->set_active_fe_index(start_diff);
631 *
EigenSolver<dim>::get_lambda_h()
633 *
return &(*eigenvalues)[0];
642 *
EigenSolver<dim>::get_solution()
652 *
EigenSolver<dim>::convert_solution()
654 *
solution.reinit((*eigenfunctions)[0].size());
655 *
for (
unsigned int i = 0; i < solution.size(); ++i)
656 *
solution[i] = (*eigenfunctions)[0][i];
666 *
template <
class SolverType>
668 *
EigenSolver<dim>::initialize_eigensolver(SolverType &eigensolver)
672 * From the parameters
class, initialize the eigensolver...
675 *
switch (this->eigenpair_selection_scheme)
678 *
eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
681 * eigensolver.set_target_eigenvalue(this->target);
686 *
eigensolver.set_which_eigenpairs(EPS_SMALLEST_MAGNITUDE);
690 *
eigensolver.set_problem_type(EPS_GHEP);
693 *
apply a Shift-Invert spectrum transformation
699 *
double shift_scalar = this->parameters->get_double(
"Target eigenvalue");
708 *
this->mpi_communicator, additional_data);
710 *
eigensolver.set_transformation(spectral_transformation);
711 *
eigensolver.set_target_eigenvalue(this->target);
720 *
EigenSolver<dim>::solve_problem()
730 *
this->mpi_communicator);
732 *
initialize_eigensolver(eigensolver);
739 *
eigensolver.solve(stiffness_matrix,
743 *
eigenfunctions->size());
744 *
for (
auto &entry : *eigenfunctions)
746 *
constraints.distribute(entry);
748 *
convert_solution();
750 *
return solver_control.last_step();
755 *
EigenSolver<dim>::n_dofs() const
757 *
return dof_handler.n_dofs();
767 *
EigenSolver<dim>::setup_system()
769 *
dof_handler.distribute_dofs(*fe_collection);
770 *
constraints.clear();
773 *
constraints.close();
775 *
eigenfunctions->resize(this->n_eigenpairs);
778 *
IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs();
780 *
for (
auto &entry : *eigenfunctions)
782 *
entry.reinit(eigenfunction_index_set, MPI_COMM_WORLD);
791 *
EigenSolver<dim>::assemble_system()
794 *
*quadrature_collection,
800 * Prep the system matrices
for the solution
803 *
stiffness_matrix.reinit(dof_handler.n_dofs(),
804 *
dof_handler.n_dofs(),
805 *
dof_handler.max_couplings_between_dofs());
807 *
dof_handler.n_dofs(),
808 *
dof_handler.max_couplings_between_dofs());
811 *
std::vector<types::global_dof_index> local_dof_indices;
813 *
for (
const auto &cell : dof_handler.active_cell_iterators())
815 *
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
818 *
cell_stiffness_matrix = 0;
821 *
cell_mass_matrix = 0;
823 *
hp_fe_values.reinit(cell);
827 *
for (
unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
836 * Note that (in general) the Nedelec element is not
837 * primitive, namely that the shape
functions are vectorial
838 * with components in more than
one direction
844 *
cell_stiffness_matrix(i, j) +=
845 *
Operations::curlcurl(fe_values, i, j, q_point) *
846 *
fe_values.JxW(q_point);
848 *
cell_mass_matrix(i, j) +=
849 *
(Operations::dot_term(fe_values, i, j, q_point)) *
850 *
fe_values.JxW(q_point);
854 *
cell->get_dof_indices(local_dof_indices);
857 *
constraints.distribute_local_to_global(cell_stiffness_matrix,
860 *
constraints.distribute_local_to_global(cell_mass_matrix,
867 *
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
868 *
if (constraints.is_constrained(i))
870 *
stiffness_matrix.set(i, i, 10000.0);
875 * since we have just
set individual elements, we need the following
887 *
class PrimalSolver :
public EigenSolver<dim>
890 *
PrimalSolver(
const std::string & prm_file,
891 *
Triangulation<dim> &triangulation,
892 *
const unsigned int &min_degree,
893 *
const unsigned int &max_degree,
894 *
const unsigned int &starting_degree);
899 *
virtual unsigned int
900 *
n_dofs()
const override;
904 *
PrimalSolver<dim>::PrimalSolver(
const std::string & prm_file,
906 *
const unsigned int &min_degree,
907 *
const unsigned int &max_degree,
908 *
const unsigned int &starting_degree)
909 *
: Base<dim>(prm_file, triangulation)
910 *
, EigenSolver<dim>(prm_file,
923 *
PrimalSolver<dim>::output_solution()
926 *
data_out.attach_dof_handler(this->dof_handler);
927 *
Vector<double> fe_degrees(this->triangulation->n_active_cells());
928 *
for (
const auto &cell : this->dof_handler.active_cell_iterators())
929 *
fe_degrees(cell->active_cell_index()) =
930 *
(*this->fe_collection)[cell->active_fe_index()].degree;
931 *
data_out.add_data_vector(fe_degrees,
"fe_degree");
932 *
data_out.add_data_vector((*this->eigenfunctions)[0],
933 *
std::string(
"eigenfunction_no_") +
936 *
std::cout <<
"Eigenvalue: " << (*this->eigenvalues)[0]
937 *
<<
" NDoFs: " << this->dof_handler.n_dofs() << std::endl;
938 *
std::ofstream eigenvalues_out(
939 *
"eigenvalues-" + std::to_string(this->refinement_cycle) +
".txt");
941 *
eigenvalues_out << std::setprecision(20) << (*this->eigenvalues)[0] <<
" "
942 *
<< this->dof_handler.n_dofs() << std::endl;
944 *
eigenvalues_out.close();
947 *
data_out.build_patches();
948 *
std::ofstream output(
"eigenvectors-" +
949 *
std::to_string(this->refinement_cycle) +
".vtu");
950 *
data_out.write_vtu(output);
955 *
PrimalSolver<dim>::n_dofs() const
957 *
return EigenSolver<dim>::n_dofs();
962 * Note, that at least
for the demonstrated problem (i.e., a Hermitian problem
963 * and eigenvalue QoI), the dual problem is identical to the primal problem;
964 * however, it is convenient to separate them in
this manner (
e.g.,
for
965 * considering functionals of the eigenfunction).
969 *
class DualSolver :
public EigenSolver<dim>
972 *
DualSolver(
const std::string & prm_file,
973 *
Triangulation<dim> &triangulation,
974 *
const unsigned int &min_degree,
975 *
const unsigned int &max_degree,
976 *
const unsigned int &starting_degree);
980 *
DualSolver<dim>::DualSolver(
const std::string & prm_file,
982 *
const unsigned int &min_degree,
983 *
const unsigned int &max_degree,
984 *
const unsigned int &starting_degree)
985 *
: Base<dim>(prm_file, triangulation)
986 *
, EigenSolver<dim>(prm_file,
998 *
namespace ErrorIndicators
1000 *
using namespace Maxwell;
1007 *
template <
int dim,
bool report_dual>
1008 *
class DualWeightedResidual :
public PrimalSolver<dim>,
public DualSolver<dim>
1012 *
output_eigenvalue_data(std::ofstream &os);
1014 *
output_qoi_error_estimates(std::ofstream &os);
1021 *
DualWeightedResidual(
const std::string & prm_file,
1022 *
Triangulation<dim> &triangulation,
1023 *
const unsigned int &min_primal_degree,
1024 *
const unsigned int &max_primal_degree,
1025 *
const unsigned int &starting_primal_degree);
1027 *
virtual unsigned int
1028 *
solve_problem()
override;
1031 *
output_solution()
override;
1033 *
virtual unsigned int
1034 *
n_dofs()
const override;
1037 *
estimate_error(Vector<double> &error_indicators);
1043 *
get_primal_DoFHandler();
1046 *
get_dual_DoFHandler();
1048 *
hp::FECollection<dim> *
1049 *
get_FECollection();
1051 *
hp::FECollection<dim> *
1052 *
get_primal_FECollection();
1054 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1055 *
get_eigenfunctions();
1057 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1058 *
get_primal_eigenfunctions();
1060 *
std::unique_ptr<std::vector<double>> &
1061 *
get_primal_eigenvalues();
1063 *
std::unique_ptr<std::vector<double>> &
1064 *
get_dual_eigenvalues();
1067 *
synchronize_discretization();
1072 *
return PrimalSolver<dim>::fe_collection->max_degree();
1074 *
double qoi_error_estimate = 0;
1078 *
embed(
const DoFHandler<dim> & dof1,
1079 *
const DoFHandler<dim> & dof2,
1080 *
const AffineConstraints<double> &constraints,
1081 *
const Vector<double> & solution,
1082 *
Vector<double> & u2);
1085 *
extract(
const DoFHandler<dim> & dof1,
1086 *
const DoFHandler<dim> & dof2,
1087 *
const AffineConstraints<double> &constraints,
1088 *
const Vector<double> & solution,
1089 *
Vector<double> & u2);
1095 *
std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values;
1096 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values;
1097 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor;
1098 *
std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values;
1100 *
std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values_forward;
1101 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_forward;
1102 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor_forward;
1103 *
std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values_forward;
1104 *
using FaceIntegrals =
1105 *
typename std::map<typename DoFHandler<dim>::face_iterator,
double>;
1108 *
solve_primal_problem();
1111 *
solve_dual_problem();
1114 *
normalize_solutions(Vector<double> &primal_solution,
1115 *
Vector<double> &dual_weights);
1118 *
get_global_QoI_error(Vector<double> &dual_solution,
1119 *
Vector<double> &error_indicators);
1122 *
initialize_error_estimation_data();
1125 *
estimate_on_one_cell(
1127 *
const Vector<double> & primal_solution,
1128 *
const Vector<double> & dual_weights,
1129 *
const double & lambda_h,
1130 *
Vector<double> & error_indicators,
1131 *
FaceIntegrals & face_integrals);
1134 *
integrate_over_cell(
1136 *
const Vector<double> & primal_solution,
1137 *
const Vector<double> & dual_weights,
1138 *
const double & lambda_h,
1139 *
Vector<double> & error_indicators);
1142 *
integrate_over_regular_face(
1144 *
const unsigned int & face_no,
1145 *
const Vector<double> & primal_solution,
1146 *
const Vector<double> & dual_weights,
1147 *
FaceIntegrals & face_integrals);
1150 *
integrate_over_irregular_face(
1152 *
const unsigned int & face_no,
1153 *
const Vector<double> & primal_solution,
1154 *
const Vector<double> & dual_weights,
1155 *
FaceIntegrals & face_integrals);
1162 *
template <
int dim,
bool report_dual>
1163 *
DualWeightedResidual<dim, report_dual>::DualWeightedResidual(
1164 *
const std::string & prm_file,
1166 *
const unsigned int &min_primal_degree,
1167 *
const unsigned int &max_primal_degree,
1168 *
const unsigned int &starting_primal_degree)
1169 *
: Base<dim>(prm_file, triangulation)
1170 *
, PrimalSolver<dim>(prm_file,
1172 *
min_primal_degree,
1173 *
max_primal_degree,
1174 *
starting_primal_degree)
1175 *
, DualSolver<dim>(prm_file,
1177 *
min_primal_degree + 1,
1178 *
max_primal_degree + 1,
1179 *
starting_primal_degree + 1)
1181 *
initialize_error_estimation_data();
1188 *
template <
int dim,
bool report_dual>
1190 *
DualWeightedResidual<dim, report_dual>::get_DoFHandler()
1193 *
return &(PrimalSolver<dim>::dof_handler);
1195 *
return &(DualSolver<dim>::dof_handler);
1200 * See above function, but to specifically output the primal
DoFHandler...
1203 *
template <
int dim,
bool report_dual>
1205 *
DualWeightedResidual<dim, report_dual>::get_primal_DoFHandler()
1207 *
return &(PrimalSolver<dim>::dof_handler);
1212 * See above function, but
for the FECollection
1215 *
template <
int dim,
bool report_dual>
1217 *
DualWeightedResidual<dim, report_dual>::get_FECollection()
1220 *
return &*(PrimalSolver<dim>::fe_collection);
1222 *
return &*(DualSolver<dim>::fe_collection);
1227 * See above function, but
for the primal FECollection
1230 *
template <
int dim,
bool report_dual>
1232 *
DualWeightedResidual<dim, report_dual>::get_primal_FECollection()
1234 *
return &*(PrimalSolver<dim>::fe_collection);
1237 *
template <
int dim,
bool report_dual>
1239 *
DualWeightedResidual<dim, report_dual>::get_dual_DoFHandler()
1241 *
return &(DualSolver<dim>::dof_handler);
1245 *
template <
int dim,
bool report_dual>
1246 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1247 *
DualWeightedResidual<dim, report_dual>::get_eigenfunctions()
1250 *
return (PrimalSolver<dim>::eigenfunctions);
1252 *
return (DualSolver<dim>::eigenfunctions);
1256 *
template <
int dim,
bool report_dual>
1257 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1258 *
DualWeightedResidual<dim, report_dual>::get_primal_eigenfunctions()
1260 *
return (PrimalSolver<dim>::eigenfunctions);
1264 *
template <
int dim,
bool report_dual>
1265 *
std::unique_ptr<std::vector<double>> &
1266 *
DualWeightedResidual<dim, report_dual>::get_primal_eigenvalues()
1268 *
return PrimalSolver<dim>::eigenvalues;
1272 *
template <
int dim,
bool report_dual>
1273 *
std::unique_ptr<std::vector<double>> &
1274 *
DualWeightedResidual<dim, report_dual>::get_dual_eigenvalues()
1276 *
return DualSolver<dim>::eigenvalues;
1279 *
template <
int dim,
bool report_dual>
1281 *
DualWeightedResidual<dim, report_dual>::output_solution()
1283 *
PrimalSolver<dim>::output_solution();
1288 * Solves the primal problem
1291 *
template <
int dim,
bool report_dual>
1293 *
DualWeightedResidual<dim, report_dual>::solve_primal_problem()
1295 *
return PrimalSolver<dim>::solve_problem();
1300 * Solves the dual problem
1303 *
template <
int dim,
bool report_dual>
1305 *
DualWeightedResidual<dim, report_dual>::solve_dual_problem()
1307 *
return DualSolver<dim>::solve_problem();
1314 *
template <
int dim,
bool report_dual>
1316 *
DualWeightedResidual<dim, report_dual>::solve_problem()
1318 *
DualWeightedResidual<dim, report_dual>::solve_primal_problem();
1319 *
return DualWeightedResidual<dim, report_dual>::solve_dual_problem();
1325 *
template <
int dim,
bool report_dual>
1327 *
DualWeightedResidual<dim, report_dual>::n_dofs() const
1329 *
return PrimalSolver<dim>::n_dofs();
1338 *
template <
int dim,
bool report_dual>
1340 *
DualWeightedResidual<dim, report_dual>::synchronize_discretization()
1354 * In
this case, we have modified the polynomial orders
for the dual;
1355 * need to update the primal
1358 *
dof1 = &(DualSolver<dim>::dof_handler);
1359 *
dof2 = &(PrimalSolver<dim>::dof_handler);
1362 *
endc1 = dof1->end();
1364 *
for (; cell1 < endc1; ++cell1, ++cell2)
1366 *
cell2->set_active_fe_index(cell1->active_fe_index());
1374 *
template <
int dim,
bool report_dual>
1376 *
DualWeightedResidual<dim, report_dual>::initialize_error_estimation_data()
1380 * initialize the cell fe_values...
1383 *
cell_hp_fe_values = std::make_unique<hp::FEValues<dim>>(
1384 *
*DualSolver<dim>::fe_collection,
1385 *
*DualSolver<dim>::quadrature_collection,
1388 *
face_hp_fe_values = std::make_unique<hp::FEFaceValues<dim>>(
1389 *
*DualSolver<dim>::fe_collection,
1390 *
*DualSolver<dim>::face_quadrature_collection,
1393 *
face_hp_fe_values_neighbor = std::make_unique<hp::FEFaceValues<dim>>(
1394 *
*DualSolver<dim>::fe_collection,
1395 *
*DualSolver<dim>::face_quadrature_collection,
1398 *
subface_hp_fe_values = std::make_unique<hp::FESubfaceValues<dim>>(
1399 *
*DualSolver<dim>::fe_collection,
1400 *
*DualSolver<dim>::face_quadrature_collection,
1409 *
template <
int dim,
bool report_dual>
1411 *
DualWeightedResidual<dim, report_dual>::normalize_solutions(
1415 *
double sum_primal = 0.0, sum_dual = 0.0;
1416 *
for (
const auto &cell :
1417 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1419 *
cell_hp_fe_values->reinit(cell);
1423 * grab the fe_values
object
1427 *
cell_hp_fe_values->get_present_fe_values();
1429 *
std::vector<Vector<double>> cell_primal_values(
1431 *
cell_dual_values(fe_values.n_quadrature_points,
Vector<double>(dim));
1432 *
fe_values.get_function_values(primal_solution, cell_primal_values);
1433 *
fe_values.get_function_values(dual_weights, cell_dual_values);
1436 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1439 *
cell_primal_values[p] * cell_primal_values[p] * fe_values.JxW(p);
1441 *
cell_dual_values[p] * cell_dual_values[p] * fe_values.JxW(p);
1445 *
primal_solution /=
sqrt(sum_primal);
1446 *
dual_weights /=
sqrt(sum_dual);
1453 *
template <
int dim,
bool report_dual>
1455 *
DualWeightedResidual<dim, report_dual>::estimate_error(
1460 * The constraints could be grabbed directly, but
this is simple
1465 *
primal_hanging_node_constraints);
1466 *
primal_hanging_node_constraints.close();
1470 *
dual_hanging_node_constraints);
1471 *
dual_hanging_node_constraints.close();
1475 * First map the primal solution to the space of the dual solution
1476 * This allows us to use just
one set of
FEValues objects (rather than one
1477 *
set for the primal, one
for dual)
1483 *
Vector<double> primal_solution(DualSolver<dim>::dof_handler.n_dofs());
1485 *
embed(PrimalSolver<dim>::dof_handler,
1486 *
DualSolver<dim>::dof_handler,
1487 *
dual_hanging_node_constraints,
1488 *
*(PrimalSolver<dim>::get_solution()),
1491 *
Vector<double> &dual_solution = *(DualSolver<dim>::get_solution());
1493 *
normalize_solutions(primal_solution, dual_solution);
1495 *
Vector<double> dual_weights(DualSolver<dim>::dof_handler.n_dofs()),
1496 *
dual_weights_interm(PrimalSolver<dim>::dof_handler.n_dofs());
1500 * First
extract the dual solution to the space of the primal
1503 *
extract(DualSolver<dim>::dof_handler,
1504 *
PrimalSolver<dim>::dof_handler,
1505 *
primal_hanging_node_constraints,
1506 *
*(DualSolver<dim>::get_solution()),
1507 *
dual_weights_interm);
1511 * Now embed
this back to the space of the dual solution
1514 *
embed(PrimalSolver<dim>::dof_handler,
1515 *
DualSolver<dim>::dof_handler,
1516 *
dual_hanging_node_constraints,
1517 *
dual_weights_interm,
1523 * Subtract
this from the full dual solution
1526 *
dual_weights -= *(DualSolver<dim>::get_solution());
1527 *
dual_weights *= -1.0;
1529 *
*(DualSolver<dim>::get_solution()) -= primal_solution;
1531 *
FaceIntegrals face_integrals;
1532 *
for (
const auto &cell :
1533 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1534 *
for (
const auto &face : cell->face_iterators())
1535 *
face_integrals[face] = -1e20;
1538 *
for (
const auto &cell :
1539 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1541 *
estimate_on_one_cell(cell,
1544 *
*(PrimalSolver<dim>::get_lambda_h()),
1548 *
unsigned int present_cell = 0;
1549 *
for (
const auto &cell :
1550 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1552 *
for (
const auto &face : cell->face_iterators())
1554 *
Assert(face_integrals.find(face) != face_integrals.end(),
1555 *
ExcInternalError());
1556 *
error_indicators(present_cell) -= 0.5 * face_integrals[face];
1563 * Now, with the error indicators computed, let us produce the
1564 * estimate of the QoI error
1567 *
this->qoi_error_estimate =
1568 *
this->get_global_QoI_error(*(DualSolver<dim>::get_solution()),
1569 *
error_indicators);
1570 *
std::cout <<
"Estimated QoI error: " << std::setprecision(20)
1571 *
<< qoi_error_estimate << std::endl;
1578 *
template <
int dim,
bool report_dual>
1580 *
DualWeightedResidual<dim, report_dual>::estimate_on_one_cell(
1584 *
const double & lambda_h,
1586 *
FaceIntegrals & face_integrals)
1588 *
integrate_over_cell(
1589 *
cell, primal_solution, dual_weights, lambda_h, error_indicators);
1592 *
if (cell->face(face_no)->at_boundary())
1594 *
face_integrals[cell->face(face_no)] = 0.0;
1597 *
if ((cell->neighbor(face_no)->has_children() ==
false) &&
1598 *
(cell->neighbor(face_no)->level() == cell->level()) &&
1599 *
(cell->neighbor(face_no)->index() < cell->index()))
1601 *
if (cell->at_boundary(face_no) ==
false)
1602 *
if (cell->neighbor(face_no)->level() < cell->level())
1604 *
if (cell->face(face_no)->has_children() ==
false)
1605 *
integrate_over_regular_face(
1606 *
cell, face_no, primal_solution, dual_weights, face_integrals);
1608 *
integrate_over_irregular_face(
1609 *
cell, face_no, primal_solution, dual_weights, face_integrals);
1616 *
template <
int dim,
bool report_dual>
1618 *
DualWeightedResidual<dim, report_dual>::integrate_over_cell(
1622 *
const double & lambda_h,
1625 *
cell_hp_fe_values->reinit(cell);
1628 * Grab the fe_values
object
1632 *
std::vector<std::vector<Tensor<2, dim, double>>> cell_hessians(
1634 *
std::vector<Vector<double>> cell_primal_values(
1636 *
cell_dual_values(fe_values.n_quadrature_points,
Vector<double>(dim));
1637 *
fe_values.get_function_values(primal_solution, cell_primal_values);
1638 *
fe_values.get_function_hessians(primal_solution, cell_hessians);
1639 *
fe_values.get_function_values(dual_weights, cell_dual_values);
1644 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1647 *
( (cell_hessians[p][1][1][0] -
1648 *
cell_hessians[p][0][1][1]) *
1649 *
(cell_dual_values[p](0)) +
1651 *
(cell_hessians[p][0][0][1] - cell_hessians[p][1][0][0]) *
1652 *
(cell_dual_values[p](1)) -
1653 *
lambda_h * (cell_primal_values[p](0) * cell_dual_values[p](0) +
1654 *
cell_primal_values[p](1) * cell_dual_values[p](1))) *
1658 *
error_indicators(cell->active_cell_index()) +=
sum;
1664 *
template <
int dim,
bool report_dual>
1666 *
DualWeightedResidual<dim, report_dual>::integrate_over_regular_face(
1668 *
const unsigned int & face_no,
1671 *
FaceIntegrals & face_integrals)
1674 *
ExcInternalError());
1675 *
const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1676 *
const auto neighbor = cell->neighbor(face_no);
1678 *
const unsigned int quadrature_index =
1679 *
std::max(cell->active_fe_index(), neighbor->active_fe_index());
1680 *
face_hp_fe_values->reinit(cell, face_no, quadrature_index);
1682 *
face_hp_fe_values->get_present_fe_values();
1683 *
std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1684 *
fe_face_values_cell.n_quadrature_points,
1686 *
neighbor_primal_grads(fe_face_values_cell.n_quadrature_points,
1688 *
fe_face_values_cell.get_function_gradients(primal_solution,
1689 *
cell_primal_grads);
1691 *
face_hp_fe_values_neighbor->reinit(neighbor,
1692 *
neighbor_neighbor,
1693 *
quadrature_index);
1695 *
face_hp_fe_values_neighbor->get_present_fe_values();
1696 *
fe_face_values_cell_neighbor.get_function_gradients(primal_solution,
1697 *
neighbor_primal_grads);
1698 *
const unsigned int n_q_points = fe_face_values_cell.n_quadrature_points;
1699 *
double face_integral = 0.0;
1700 *
std::vector<Vector<double>> cell_dual_values(
n_q_points,
1702 *
fe_face_values_cell.get_function_values(dual_weights, cell_dual_values);
1703 *
for (
unsigned int p = 0; p <
n_q_points; ++p)
1705 *
auto face_normal = fe_face_values_cell.normal_vector(p);
1708 *
(cell_primal_grads[p][1][0] - cell_primal_grads[p][0][1] -
1709 *
neighbor_primal_grads[p][1][0] + neighbor_primal_grads[p][0][1]) *
1710 *
(cell_dual_values[p][0] * face_normal[1] -
1711 *
cell_dual_values[p][1] * face_normal[0]) *
1712 *
fe_face_values_cell.JxW(p);
1714 *
Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
1715 *
ExcInternalError());
1716 *
Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
1717 *
face_integrals[cell->face(face_no)] = face_integral;
1723 *
template <
int dim,
bool report_dual>
1725 *
DualWeightedResidual<dim, report_dual>::integrate_over_irregular_face(
1727 *
const unsigned int & face_no,
1730 *
FaceIntegrals & face_integrals)
1734 *
cell->neighbor(face_no);
1737 *
Assert(neighbor->has_children(), ExcInternalError());
1739 *
const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1740 *
for (
unsigned int subface_no = 0; subface_no < face->n_children();
1744 *
cell->neighbor_child_on_subface(face_no, subface_no);
1745 *
Assert(neighbor_child->face(neighbor_neighbor) ==
1746 *
cell->face(face_no)->child(subface_no),
1747 *
ExcInternalError());
1748 *
const unsigned int quadrature_index =
1749 *
std::max(cell->active_fe_index(), neighbor_child->active_fe_index());
1752 * initialize fe_subface values_cell
1755 *
subface_hp_fe_values->reinit(cell,
1758 *
quadrature_index);
1760 *
subface_hp_fe_values->get_present_fe_values();
1761 *
std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1762 *
subface_fe_values_cell.n_quadrature_points,
1764 *
neighbor_primal_grads(subface_fe_values_cell.n_quadrature_points,
1766 *
subface_fe_values_cell.get_function_gradients(primal_solution,
1767 *
cell_primal_grads);
1770 * initialize fe_face_values_neighbor
1773 *
face_hp_fe_values_neighbor->reinit(neighbor_child,
1774 *
neighbor_neighbor,
1775 *
quadrature_index);
1777 *
face_hp_fe_values_neighbor->get_present_fe_values();
1778 *
face_fe_values_neighbor.get_function_gradients(primal_solution,
1779 *
neighbor_primal_grads);
1781 *
subface_fe_values_cell.n_quadrature_points;
1782 *
std::vector<Vector<double>> cell_dual_values(
n_q_points,
1784 *
face_fe_values_neighbor.get_function_values(dual_weights,
1785 *
cell_dual_values);
1787 *
double face_integral = 0.0;
1789 *
for (
unsigned int p = 0; p <
n_q_points; ++p)
1791 *
auto face_normal = face_fe_values_neighbor.normal_vector(p);
1793 *
(cell_primal_grads[p][0][1] - cell_primal_grads[p][1][0] +
1794 *
neighbor_primal_grads[p][1][0] -
1795 *
neighbor_primal_grads[p][0][1]) *
1796 *
(cell_dual_values[p][0] * face_normal[1] -
1797 *
cell_dual_values[p][1] * face_normal[0]) *
1798 *
face_fe_values_neighbor.JxW(p);
1800 *
face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
1803 *
for (
unsigned int subface_no = 0; subface_no < face->n_children();
1806 *
Assert(face_integrals.find(face->child(subface_no)) !=
1807 *
face_integrals.end(),
1808 *
ExcInternalError());
1809 *
Assert(face_integrals[face->child(subface_no)] != -1e20,
1810 *
ExcInternalError());
1811 *
sum += face_integrals[face->child(subface_no)];
1813 *
face_integrals[face] =
sum;
1816 *
template <
int dim,
bool report_dual>
1818 *
DualWeightedResidual<dim, report_dual>::get_global_QoI_error(
1822 *
auto dual_less_primal =
1826 *
double scaling_factor = 0.0;
1827 *
for (
const auto &cell :
1828 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1830 *
cell_hp_fe_values->reinit(cell);
1833 * grab the fe_values
object
1837 *
cell_hp_fe_values->get_present_fe_values();
1839 *
std::vector<Vector<double>> cell_values(fe_values.n_quadrature_points,
1841 *
fe_values.get_function_values(dual_less_primal, cell_values);
1843 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1846 *
(cell_values[p] * cell_values[p]) * fe_values.JxW(p);
1849 *
double global_QoI_error = 0.0;
1850 *
for (
const auto &indicator : error_indicators)
1852 *
global_QoI_error += indicator;
1855 *
global_QoI_error /= (1 - 0.5 * scaling_factor);
1856 *
return global_QoI_error;
1860 *
template <
int dim,
bool report_dual>
1862 *
DualWeightedResidual<dim, report_dual>::embed(
1869 *
assert(u2.size() == dof2.n_dofs() &&
"Incorrect input vector size!");
1874 *
endc1 = dof1.end();
1877 *
for (; cell1 < endc1; ++cell1, ++cell2)
1884 *
assert(fe1.degree < fe2.degree &&
"Incorrect usage of embed!");
1888 * Get the embedding_dofs
1897 *
std::vector<unsigned int> embedding_dofs =
1898 *
fe2.get_embedding_dofs(fe1.degree);
1899 *
const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1905 *
local_dof_values_1.reinit(fe1.dofs_per_cell);
1906 *
cell1->get_dof_values(solution, local_dof_values_1);
1908 *
for (
unsigned int i = 0; i < local_dof_values_1.size(); ++i)
1909 *
local_dof_values_2[embedding_dofs[i]] = local_dof_values_1[i];
1913 * Now
set this changes to the global vector
1916 *
cell2->set_dof_values(local_dof_values_2, u2);
1922 * Applies the constraints of the target finite element space
1925 *
constraints.distribute(u2);
1928 *
template <
int dim,
bool report_dual>
1930 *
DualWeightedResidual<dim, report_dual>::extract(
1939 * Maps from fe1 to fe2
1942 *
assert(u2.size() == dof2.n_dofs() &&
"Incorrect input vector size!");
1947 *
endc1 = dof1.end();
1950 *
for (; cell1 < endc1; ++cell1, ++cell2)
1957 *
assert(fe1.degree > fe2.degree &&
"Incorrect usage of extract!");
1961 * Get the embedding_dofs
1964 *
std::vector<unsigned int> embedding_dofs =
1965 *
fe1.get_embedding_dofs(fe2.degree);
1966 *
const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1972 *
local_dof_values_1.reinit(fe1.dofs_per_cell);
1973 *
cell1->get_dof_values(solution, local_dof_values_1);
1975 *
for (
unsigned int i = 0; i < local_dof_values_2.size(); ++i)
1976 *
local_dof_values_2[i] = local_dof_values_1[embedding_dofs[i]];
1980 * Now
set this changes to the global vector
1983 *
cell2->set_dof_values(local_dof_values_2, u2);
1989 * Applies the constraints of the target finite element space
1992 *
constraints.distribute(u2);
1994 *
template <
int dim,
bool report_dual>
1996 *
DualWeightedResidual<dim, report_dual>::output_eigenvalue_data(
1997 *
std::ofstream &os)
1999 *
os << (*this->get_primal_eigenvalues())[0] <<
" "
2000 *
<< (this->get_primal_DoFHandler())->n_dofs() <<
" "
2001 *
<< (*this->get_dual_eigenvalues())[0] <<
" "
2002 *
<< (this->get_dual_DoFHandler())->n_dofs() << std::endl;
2004 *
template <
int dim,
bool report_dual>
2006 *
DualWeightedResidual<dim, report_dual>::output_qoi_error_estimates(
2007 *
std::ofstream &os)
2009 *
os << qoi_error_estimate << std::endl;
2016 *
template <
int dim>
2017 *
class KellyErrorIndicator :
public PrimalSolver<dim>
2026 *
output_eigenvalue_data(std::ofstream &os);
2028 *
output_qoi_error_estimates(std::ofstream &);
2029 *
KellyErrorIndicator(
const std::string & prm_file,
2030 *
Triangulation<dim> &coarse_grid,
2031 *
const unsigned int &min_degree,
2032 *
const unsigned int &max_degree,
2033 *
const unsigned int &starting_degree);
2035 *
virtual unsigned int
2036 *
solve_problem()
override;
2039 *
output_solution()
override;
2041 *
hp::FECollection<dim> *
2042 *
get_FECollection();
2044 *
hp::FECollection<dim> *
2045 *
get_primal_FECollection();
2047 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2048 *
get_eigenfunctions();
2050 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2051 *
get_primal_eigenfunctions();
2053 *
std::unique_ptr<std::vector<double>> &
2054 *
get_primal_eigenvalues();
2058 *
synchronize_discretization();
2064 *
get_primal_DoFHandler();
2069 *
return PrimalSolver<dim>::fe_collection->max_degree();
2071 *
double qoi_error_estimate = 0;
2075 *
estimate_error(Vector<double> &error_indicators);
2079 *
prune_eigenpairs(
const double &TOL);
2081 *
std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
2082 *
std::vector<const double *> eigenvalue_ptrs;
2084 *
std::vector<std::shared_ptr<Vector<float>>> errors;
2087 *
template <
int dim>
2088 *
KellyErrorIndicator<dim>::KellyErrorIndicator(
2089 *
const std::string & prm_file,
2091 *
const unsigned int &min_degree,
2092 *
const unsigned int &max_degree,
2093 *
const unsigned int &starting_degree)
2094 *
: Base<dim>(prm_file, coarse_grid)
2095 *
, PrimalSolver<dim>(prm_file,
2102 *
template <
int dim>
2104 *
KellyErrorIndicator<dim>::solve_problem()
2106 *
return PrimalSolver<dim>::solve_problem();
2109 *
template <
int dim>
2111 *
KellyErrorIndicator<dim>::get_FECollection()
2113 *
return &*(PrimalSolver<dim>::fe_collection);
2116 *
template <
int dim>
2118 *
KellyErrorIndicator<dim>::get_primal_FECollection()
2120 *
return &*(PrimalSolver<dim>::fe_collection);
2123 *
template <
int dim>
2124 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2125 *
KellyErrorIndicator<dim>::get_eigenfunctions()
2127 *
return (PrimalSolver<dim>::eigenfunctions);
2130 *
template <
int dim>
2131 *
std::unique_ptr<std::vector<double>> &
2132 *
KellyErrorIndicator<dim>::get_primal_eigenvalues()
2134 *
return PrimalSolver<dim>::eigenvalues;
2137 *
template <
int dim>
2138 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2139 *
KellyErrorIndicator<dim>::get_primal_eigenfunctions()
2141 *
return (PrimalSolver<dim>::eigenfunctions);
2144 *
template <
int dim>
2146 *
KellyErrorIndicator<dim>::get_DoFHandler()
2148 *
return &(PrimalSolver<dim>::dof_handler);
2151 *
template <
int dim>
2153 *
KellyErrorIndicator<dim>::get_primal_DoFHandler()
2155 *
return &(PrimalSolver<dim>::dof_handler);
2158 *
template <
int dim>
2160 *
KellyErrorIndicator<dim>::synchronize_discretization()
2164 * This function does
nothing for this error indicator
2170 *
template <
int dim>
2172 *
KellyErrorIndicator<dim>::output_solution()
2174 *
PrimalSolver<dim>::output_solution();
2177 *
template <
int dim>
2179 *
KellyErrorIndicator<dim>::prune_eigenpairs(
const double &TOL)
2181 *
unsigned int count = 0;
2182 *
for (
size_t eigenpair_index = 0;
2183 *
eigenpair_index < this->eigenfunctions->size();
2184 *
++eigenpair_index)
2186 *
if (count >= this->n_eigenpairs)
2188 *
if (
abs((*this->eigenvalues)[eigenpair_index]) < TOL)
2191 *
eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
2192 *
eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
2196 *
template <
int dim>
2198 *
KellyErrorIndicator<dim>::estimate_error(
Vector<double> &error_indicators)
2200 *
std::cout <<
"Marking cells via Kelly indicator..." << std::endl;
2201 *
prune_eigenpairs(1e-9);
2204 * deallocate the errors vector
2208 *
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2210 *
errors.emplace_back(
2211 *
new Vector<float>(this->triangulation->n_active_cells()));
2213 *
std::vector<Vector<float> *> estimated_error_per_cell(
2214 *
eigenfunction_ptrs.size());
2215 *
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2217 *
estimated_error_per_cell[i] = errors[i].get();
2223 *
*this->face_quadrature_collection,
2228 *
for (
auto &error_vec : errors)
2230 *
auto normalized_vec = *error_vec;
2231 *
normalized_vec /= normalized_vec.l1_norm();
2233 *
for (
unsigned int i = 0; i < error_indicators.size(); ++i)
2234 *
error_indicators(i) += double(normalized_vec(i));
2236 *
std::cout <<
"...Done!" << std::endl;
2238 *
template <
int dim>
2240 *
KellyErrorIndicator<dim>::output_eigenvalue_data(std::ofstream &os)
2242 *
os << (*this->get_primal_eigenvalues())[0] <<
" "
2243 *
<< (this->get_primal_DoFHandler())->n_dofs() << std::endl;
2245 *
template <
int dim>
2247 *
KellyErrorIndicator<dim>::output_qoi_error_estimates(std::ofstream &)
2257 *
namespace RegularityIndicators
2259 *
using namespace dealii;
2263 *
template <
int dim>
2264 *
class LegendreInfo
2268 *
class LegendreInfo<2>
2271 *
std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
2279 *
assert(fe_collection !=
nullptr && dof_handler !=
nullptr &&
2280 *
"A valid FECollection and DoFHandler must be accessible!");
2282 *
legendre_u = std::make_unique<FESeries::Legendre<2>>(
2284 *
legendre_v = std::make_unique<FESeries::Legendre<2>>(
2287 *
legendre_u->precalculate_all_transformation_matrices();
2288 *
legendre_v->precalculate_all_transformation_matrices();
2291 *
template <
class VectorType>
2293 *
compute_coefficient_decay(
const VectorType & eigenfunction,
2294 *
std::vector<double> &smoothness_indicators)
2298 * Compute the coefficients
for the u and v components of the solution
2303 *
smoothness_v(smoothness_indicators.size());
2315 *
for (
unsigned int i = 0; i < smoothness_indicators.size(); ++i)
2317 *
smoothness_indicators[i] =
std::min(smoothness_u[i], smoothness_v[i]);
2325 *
template <
int dim>
2326 *
class LegendreIndicator
2334 *
template <
class VectorType>
2336 *
estimate_smoothness(
2337 *
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2338 *
const unsigned int & index_of_goal,
2339 *
std::vector<double> & smoothness_indicators);
2345 *
template <
int dim>
2347 *
LegendreIndicator<dim>::attach_FE_info_and_initialize(
2356 *
template <
int dim>
2357 *
template <
class VectorType>
2359 *
LegendreIndicator<dim>::estimate_smoothness(
2360 *
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2361 *
const unsigned int & index_of_goal,
2362 *
std::vector<double> & smoothness_indicators)
2364 *
this->
legendre.compute_coefficient_decay((*eigenfunctions)[index_of_goal],
2365 *
smoothness_indicators);
2373 *
namespace Refinement
2375 *
using namespace dealii;
2376 *
using namespace Maxwell;
2378 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2379 *
class Refiner :
public ErrorIndicator,
public RegularityIndicator
2382 *
Refiner(
const std::string & prm_file,
2383 *
Triangulation<dim> &coarse_grid,
2384 *
const unsigned int &min_degree,
2385 *
const unsigned int &max_degree,
2386 *
const unsigned int &starting_degree);
2389 *
execute_refinement(
const double &smoothness_threshold_fraction);
2392 *
output_solution()
override;
2395 *
Vector<double> estimated_error_per_cell;
2396 *
std::vector<double> smoothness_indicators;
2397 *
std::ofstream eigenvalues_out;
2398 *
std::ofstream error_estimate_out;
2401 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2402 *
Refiner<dim, ErrorIndicator, RegularityIndicator>::Refiner(
2403 *
const std::string & prm_file,
2405 *
const unsigned int &min_degree,
2406 *
const unsigned int &max_degree,
2407 *
const unsigned int &starting_degree)
2408 *
: Base<dim>(prm_file, coarse_grid)
2409 *
, ErrorIndicator(prm_file,
2414 *
, RegularityIndicator()
2416 *
if (ErrorIndicator::name() ==
"DWR")
2418 *
error_estimate_out.open(
"error_estimate.txt");
2419 *
error_estimate_out << std::setprecision(20);
2422 *
eigenvalues_out.open(
"eigenvalues_" + ErrorIndicator::name() +
"_out.txt");
2423 *
eigenvalues_out << std::setprecision(20);
2428 * For generating samples of the curl of the electric field
2431 *
template <
int dim>
2435 *
CurlPostprocessor()
2441 *
const DataPostprocessorInputs::Vector<dim> &input_data,
2442 *
std::vector<Vector<double>> &computed_quantities)
const override
2445 *
computed_quantities.size());
2446 *
for (
unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
2448 *
computed_quantities[p](0) = input_data.solution_gradients[p][1][0] -
2449 *
input_data.solution_gradients[p][0][1];
2461 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2463 *
Refiner<dim, ErrorIndicator, RegularityIndicator>::output_solution()
2465 *
CurlPostprocessor<dim> curl_u;
2468 *
auto & output_dof = *(ErrorIndicator::get_primal_DoFHandler());
2469 *
data_out.attach_dof_handler(output_dof);
2470 *
Vector<double> fe_degrees(this->triangulation->n_active_cells());
2471 *
for (
const auto &cell : output_dof.active_cell_iterators())
2472 *
fe_degrees(cell->active_cell_index()) =
2473 *
(*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
2475 *
data_out.add_data_vector(fe_degrees,
"fe_degree");
2477 *
data_out.add_data_vector(estimated_error_per_cell,
"error");
2478 *
Vector<double> smoothness_out(this->triangulation->n_active_cells());
2479 *
for (
const auto &cell : output_dof.active_cell_iterators())
2481 *
auto i = cell->active_cell_index();
2482 *
if (!cell->refine_flag_set() && !cell->coarsen_flag_set())
2483 *
smoothness_out(i) = -1;
2485 *
smoothness_out(i) = smoothness_indicators[i];
2487 *
data_out.add_data_vector(smoothness_out,
"smoothness");
2488 *
data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2489 *
std::string(
"eigenfunction_no_") +
2491 *
data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2494 *
ErrorIndicator::output_eigenvalue_data(eigenvalues_out);
2495 *
ErrorIndicator::output_qoi_error_estimates(error_estimate_out);
2497 *
std::cout <<
"Number of DoFs: " << (this->get_primal_DoFHandler())->n_dofs()
2501 *
data_out.build_patches();
2502 *
std::ofstream output(
"eigenvectors-" + ErrorIndicator::name() +
"-" +
2503 *
std::to_string(this->refinement_cycle) + +
".vtu");
2504 *
data_out.write_vtu(output);
2513 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2515 *
Refiner<dim, ErrorIndicator, RegularityIndicator>::execute_refinement(
2516 *
const double &smoothness_threshold_fraction)
2520 * First initialize the RegularityIndicator...
2521 * Depending on the limits
set,
this may take a
while
2524 *
std::cout <<
"Initializing RegularityIndicator..." << std::endl;
2526 *
<<
"(This may take a while if the max expansion order is set too high)"
2528 *
RegularityIndicator::attach_FE_info_and_initialize(
2529 *
ErrorIndicator::get_FECollection(), ErrorIndicator::get_DoFHandler());
2530 *
std::cout <<
"Done!" << std::endl <<
"Starting Refinement..." << std::endl;
2532 *
for (
unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
2534 *
this->set_refinement_cycle(cycle);
2535 *
std::cout <<
"Cycle: " << cycle << std::endl;
2536 *
ErrorIndicator::solve_problem();
2537 *
this->estimated_error_per_cell.reinit(
2538 *
this->triangulation->n_active_cells());
2540 *
ErrorIndicator::estimate_error(estimated_error_per_cell);
2544 * Depending on the source of the error estimation/indication, these
2545 * values might be signed, so we address that with the following
2548 *
for (
double &error_indicator : estimated_error_per_cell)
2549 *
error_indicator =
std::abs(error_indicator);
2553 *
*this->triangulation, estimated_error_per_cell, 1. / 5., 0.000);
2557 * Now get regularity indicators
2558 * For those elements which must be refined,
swap to increasing @f$p@f$
2559 * depending on the regularity threshold...
2565 *
smoothness_indicators =
2566 *
std::vector<double>(this->triangulation->n_active_cells(),
2567 *
std::numeric_limits<double>::max());
2568 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2569 *
ErrorIndicator::PrimalSolver::max_degree)
2570 *
RegularityIndicator::estimate_smoothness(
2571 *
ErrorIndicator::get_eigenfunctions(), 0, smoothness_indicators);
2577 *
this->output_solution();
2578 *
const double threshold_smoothness = smoothness_threshold_fraction;
2579 *
unsigned int num_refined = 0, num_coarsened = 0;
2580 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2581 *
ErrorIndicator::PrimalSolver::max_degree)
2583 *
for (
const auto &cell :
2584 *
ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2586 *
if (cell->refine_flag_set())
2588 *
if (cell->coarsen_flag_set())
2590 *
if (cell->refine_flag_set() &&
2591 *
smoothness_indicators[cell->active_cell_index()] >
2592 *
threshold_smoothness &&
2593 *
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2594 *
ErrorIndicator::get_FECollection()->size())
2596 *
cell->clear_refine_flag();
2597 *
cell->set_active_fe_index(cell->active_fe_index() + 1);
2599 *
else if (cell->coarsen_flag_set() &&
2600 *
smoothness_indicators[cell->active_cell_index()] <
2601 *
threshold_smoothness &&
2602 *
cell->active_fe_index() != 0)
2604 *
cell->clear_coarsen_flag();
2606 *
cell->set_active_fe_index(cell->active_fe_index() - 1);
2610 * Here we also impose a limit on how small the cells can become
2613 *
else if (cell->refine_flag_set() && cell->diameter() < 5.0e-6)
2615 *
cell->clear_refine_flag();
2616 *
if (
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2617 *
ErrorIndicator::get_FECollection()->size())
2618 *
cell->set_active_fe_index(cell->active_fe_index() + 1);
2625 * Check what the smallest
diameter is
2628 *
double min_diameter = std::numeric_limits<double>::max();
2629 *
for (
const auto &cell :
2630 *
ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2631 *
if (cell->diameter() < min_diameter)
2632 *
min_diameter = cell->diameter();
2634 *
std::cout <<
"Min diameter: " << min_diameter << std::endl;
2636 *
ErrorIndicator::synchronize_discretization();
2638 *
(this->triangulation)->execute_coarsening_and_refinement();
2644 *
main(
int argc,
char **argv)
2648 *
using namespace dealii;
2649 *
using namespace Maxwell;
2650 *
using namespace Refinement;
2651 *
using namespace ErrorIndicators;
2652 *
using namespace RegularityIndicators;
2660 *
ExcMessage(
"This program can only be run in serial, use ./maxwell-hp"));
2663 *
Structures::create_L_waveguide(triangulation_DWR, 2.0);
2664 *
Structures::create_L_waveguide(triangulation_Kelly, 2.0);
2666 *
Refiner<2, KellyErrorIndicator<2>, LegendreIndicator<2>> problem_Kelly(
2668 *
triangulation_Kelly,
2673 *
Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
2674 *
problem_DWR(
"maxwell-hp.prm",
2675 *
triangulation_DWR,
2682 * The threshold
for the
hp-decision: too small -> not enough
2683 * @f$h@f$-refinement, too large -> not enough @f$p@f$-refinement
2686 *
double smoothness_threshold = 0.75;
2688 *
std::cout <<
"Executing refinement for the Kelly strategy!" << std::endl;
2689 *
problem_Kelly.execute_refinement(smoothness_threshold);
2690 *
std::cout <<
"...Done with Kelly refinement strategy!" << std::endl;
2691 *
std::cout <<
"Executing refinement for the DWR strategy!" << std::endl;
2692 *
problem_DWR.execute_refinement(smoothness_threshold);
2693 *
std::cout <<
"...Done with DWR refinement strategy!" << std::endl;
2696 *
catch (std::exception &exc)
2698 *
std::cerr << std::endl
2700 *
<<
"----------------------------------------------------"
2702 *
std::cerr <<
"Exception on processing: " << std::endl
2703 *
<< exc.what() << std::endl
2704 *
<<
"Aborting!" << std::endl
2705 *
<<
"----------------------------------------------------"
2712 *
std::cerr << std::endl
2714 *
<<
"----------------------------------------------------"
2716 *
std::cerr <<
"Unknown exception!" << std::endl
2717 *
<<
"Aborting!" << std::endl
2718 *
<<
"----------------------------------------------------"
2723 *
std::cout << std::endl <<
" Job done." << std::endl;
* * * struct InterferenceTaperTransform *
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
DataPostprocessorScalar(const std::string &name, const UpdateFlags update_flags)
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
const FEValues< dim, spacedim > & get_present_fe_values() const
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
#define DEAL_II_VERSION_GTE(major, minor, subminor)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
typename ActiveSelector::cell_iterator cell_iterator
typename ActiveSelector::face_iterator face_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_zero_boundary_constraints(const DoFHandler< dim, spacedim > &dof, const types::boundary_id boundary_id, AffineConstraints< number > &zero_boundary_constraints, const ComponentMask &component_mask={})
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ valid
Iterator points to a valid object.
constexpr types::blas_int one
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, ShapeDataMemorySpace > shape_data, const ViewTypeIn in, ViewTypeOut out)
FESeries::Legendre< dim, spacedim > default_fe_series(const hp::FECollection< dim, spacedim > &fe_collection, const unsigned int component=numbers::invalid_unsigned_int)
void coefficient_decay(FESeries::Legendre< dim, spacedim > &fe_legendre, const DoFHandler< dim, spacedim > &dof_handler, const VectorType &solution, Vector< float > &smoothness_indicators, const VectorTools::NormType regression_strategy=VectorTools::Linfty_norm, const double smallest_abs_coefficient=1e-10, const bool only_flagged_cells=false)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void save(Archive &ar, const ::std_cxx26::inplace_vector< T, N > &vec, const unsigned int)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(& functions)(const void *v1, const void *v2)
double legendre(unsigned int l, double x)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
* const unsigned int dofs_per_cell
* const unsigned int n_q_points
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
1.16.1