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>
258 *
#include <deal.II/grid/cell_data.h>
259 *
#include <deal.II/grid/tria.h>
260 *
#include <deal.II/grid/tria_iterator.h>
262 *
#include <deal.II/lac/affine_constraints.h>
263 *
#include <deal.II/lac/full_matrix.h>
264 *
#include <deal.II/lac/petsc_precondition.h>
265 *
#include <deal.II/lac/petsc_sparse_matrix.h>
266 *
#include <deal.II/lac/petsc_vector.h>
267 *
#include <deal.II/lac/slepc_solver.h>
269 *
#include <deal.II/numerics/data_out.h>
270 *
#include <deal.II/numerics/vector_tools.h>
274 * For parallelization (
using WorkStream and Intel TBB)
277 *
#include <deal.II/base/multithread_info.h>
278 *
#include <deal.II/base/work_stream.h>
280 *
#include
"petscpc.h"
284 * For Error Estimation/Indication and Smoothness Indication
287 *
#include <deal.II/fe/fe_tools.h>
289 *
#include <deal.II/numerics/error_estimator.h>
290 *
#include <deal.II/numerics/smoothness_estimator.h>
296 *
#include <deal.II/grid/grid_refinement.h>
299 *
#include <iostream>
302 *
namespace Operations
309 *
curlcurl(const ::FEValues<2> &fe_values,
310 *
const unsigned int & i,
311 *
const unsigned int & j,
312 *
const unsigned int & q_point)
314 *
auto gradu1_x1x2 = fe_values.shape_grad_component(i, q_point, 0);
315 *
auto gradu2_x1x2 = fe_values.shape_grad_component(i, q_point, 1);
317 *
auto gradv1_x1x2 = fe_values.shape_grad_component(j, q_point, 0);
318 *
auto gradv2_x1x2 = fe_values.shape_grad_component(j, q_point, 1);
319 *
return (gradu2_x1x2[0] - gradu1_x1x2[1]) *
320 *
(gradv2_x1x2[0] - gradv1_x1x2[1]);
328 *
dot_term(const ::FEValues<dim> &fe_values,
329 *
const unsigned int & i,
330 *
const unsigned int & j,
331 *
const unsigned int & q_point)
333 *
double output = 0.0;
334 *
for (
unsigned int comp = 0; comp < dim; ++comp)
336 *
output += fe_values.shape_value_component(i, q_point, comp) *
337 *
fe_values.shape_value_component(j, q_point, comp);
347 *
namespace Structures
352 *
create_L_waveguide(
Triangulation<2> &triangulation,
const double &scaling)
354 *
const unsigned int dim = 2;
356 *
const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
357 *
{scaling * 0.5, scaling * 0.0},
358 *
{scaling * 0.0, scaling * 0.5},
359 *
{scaling * 0.5, scaling * 0.5},
360 *
{scaling * 0.0, scaling * 1.0},
361 *
{scaling * 0.5, scaling * 1.0},
362 *
{scaling * 1.0, scaling * 0.5},
363 *
{scaling * 1.0, scaling * 1.0}};
365 *
const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
366 *
cell_vertices = {{{0, 1, 2, 3}}, {{2, 3, 4, 5}}, {{3, 6, 5, 7}}};
367 *
const unsigned int n_cells = cell_vertices.size();
368 *
std::vector<CellData<dim>> cells(n_cells,
CellData<dim>());
369 *
for (
unsigned int i = 0; i <
n_cells; ++i)
371 *
for (
unsigned int j = 0; j < cell_vertices[i].size(); ++j)
372 *
cells[i].vertices[j] = cell_vertices[i][j];
373 *
cells[i].material_id = 0;
375 *
triangulation.create_triangulation(vertices, cells,
SubCellData());
376 *
triangulation.refine_global(1);
382 *
const double & scaling)
384 *
const unsigned int dim = 2;
386 *
const std::vector<Point<2>> vertices = {{scaling * 0.0, scaling * 0.0},
387 *
{scaling * 0.6, scaling * 0.0},
388 *
{scaling * 0.0, scaling * 0.3},
389 *
{scaling * 0.6, scaling * 0.3}};
391 *
const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
392 *
cell_vertices = {{{0, 1, 2, 3}}};
393 *
const unsigned int n_cells = cell_vertices.size();
394 *
std::vector<CellData<dim>> cells(n_cells,
CellData<dim>());
395 *
for (
unsigned int i = 0; i <
n_cells; ++i)
397 *
for (
unsigned int j = 0; j < cell_vertices[i].size(); ++j)
398 *
cells[i].vertices[j] = cell_vertices[i][j];
399 *
cells[i].material_id = 0;
401 *
triangulation.create_triangulation(vertices, cells,
SubCellData());
402 *
triangulation.refine_global(0);
429 *
virtual unsigned int
430 *
solve_problem() = 0;
432 *
set_refinement_cycle(
const unsigned int cycle);
435 *
output_solution() = 0;
440 *
unsigned int refinement_cycle = 0;
441 *
std::unique_ptr<ParameterHandler> parameters;
442 *
unsigned int n_eigenpairs = 1;
443 *
double target = 0.0;
444 *
unsigned int eigenpair_selection_scheme;
445 *
unsigned int max_cycles = 0;
446 *
ompi_communicator_t * mpi_communicator = PETSC_COMM_SELF;
454 *
: triangulation(&coarse_grid)
455 *
, parameters(std::make_unique<ParameterHandler>())
457 *
parameters->declare_entry(
458 *
"Eigenpair selection scheme",
461 *
"The type of eigenpairs to find (0 - smallest, 1 - target)");
462 *
parameters->declare_entry(
"Number of eigenvalues/eigenfunctions",
465 *
"The number of eigenvalues/eigenfunctions "
466 *
"to be computed.");
467 *
parameters->declare_entry(
"Target eigenvalue",
470 *
"The target eigenvalue (if scheme == 1)");
472 *
parameters->declare_entry(
"Cycles number",
475 *
"The number of cycles in refinement");
476 *
parameters->parse_input(prm_file);
478 *
eigenpair_selection_scheme =
479 *
parameters->get_integer(
"Eigenpair selection scheme");
483 * The
project currently only supports selection by a target eigenvalue.
484 * Furthermore, only
one eigenpair can be computed at a time.
487 *
assert(eigenpair_selection_scheme == 1 &&
488 *
"Selection by a target is the only currently supported option!");
490 *
parameters->get_integer(
"Number of eigenvalues/eigenfunctions");
492 *
n_eigenpairs == 1 &&
493 *
"Only the computation of a single eigenpair is currently supported!");
495 *
target = parameters->get_double(
"Target eigenvalue");
496 *
max_cycles = parameters->get_integer(
"Cycles number");
497 *
if (eigenpair_selection_scheme == 1)
503 *
Base<dim>::set_refinement_cycle(
const unsigned int cycle)
505 *
refinement_cycle = cycle;
514 *
class EigenSolver :
public virtual Base<dim>
517 *
EigenSolver(
const std::string & prm_file,
519 *
const unsigned int &minimum_degree,
520 *
const unsigned int &maximum_degree,
521 *
const unsigned int &starting_degree);
523 *
virtual unsigned int
524 *
solve_problem()
override;
526 *
virtual unsigned int
529 *
template <
class SolverType>
531 *
initialize_eigensolver(SolverType &eigensolver);
540 *
const std::unique_ptr<hp::FECollection<dim>> fe_collection;
541 *
std::unique_ptr<hp::QCollection<dim>> quadrature_collection;
542 *
std::unique_ptr<
hp::QCollection<dim - 1>> face_quadrature_collection;
544 *
const unsigned int max_degree, min_degree;
547 *
for the actual solution
550 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> eigenfunctions;
551 *
std::unique_ptr<std::vector<double>>
eigenvalues;
561 *
convert_solution();
573 *
EigenSolver<dim>::EigenSolver(
const std::string & prm_file,
575 *
const unsigned int &minimum_degree,
576 *
const unsigned int &maximum_degree,
577 *
const unsigned int &starting_degree)
578 *
: Base<dim>(prm_file, triangulation)
582 *
, dof_handler(triangulation)
583 *
, max_degree(maximum_degree)
584 *
, min_degree(minimum_degree)
586 *
std::make_unique<std::vector<PETScWrappers::MPI::Vector>>())
587 *
,
eigenvalues(std::make_unique<std::vector<double>>())
589 *
for (
unsigned int degree = min_degree; degree <= max_degree; ++degree)
594 * Generate quadrature collection with
sorted quadrature weights
599 *
quadrature_collection->push_back(sorted_quadrature);
601 *
const QGauss<dim - 1> face_quadrature(degree + 1);
602 *
const QSorted<dim - 1> sorted_face_quadrature(face_quadrature);
603 *
face_quadrature_collection->push_back(sorted_face_quadrature);
607 * adjust the discretization
610 *
if (starting_degree > min_degree && starting_degree <= max_degree)
612 *
const unsigned int start_diff = starting_degree - min_degree;
614 *
cell1 = dof_handler.begin_active(),
615 *
endc1 = dof_handler.end();
616 *
for (; cell1 < endc1; ++cell1)
618 *
cell1->set_active_fe_index(start_diff);
629 *
EigenSolver<dim>::get_lambda_h()
631 *
return &(*eigenvalues)[0];
640 *
EigenSolver<dim>::get_solution()
650 *
EigenSolver<dim>::convert_solution()
652 *
solution.
reinit((*eigenfunctions)[0].
size());
653 *
for (
unsigned int i = 0; i < solution.size(); ++i)
654 *
solution[i] = (*eigenfunctions)[0][i];
664 *
template <
class SolverType>
666 *
EigenSolver<dim>::initialize_eigensolver(SolverType &eigensolver)
670 * From the parameters
class, initialize the eigensolver...
673 *
switch (this->eigenpair_selection_scheme)
676 *
eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE);
679 * eigensolver.set_target_eigenvalue(this->target);
684 *
eigensolver.set_which_eigenpairs(EPS_SMALLEST_MAGNITUDE);
688 *
eigensolver.set_problem_type(EPS_GHEP);
691 *
apply a Shift-Invert spectrum transformation
697 *
double shift_scalar = this->parameters->get_double(
"Target eigenvalue");
706 *
this->mpi_communicator, additional_data);
708 *
eigensolver.set_transformation(spectral_transformation);
709 *
eigensolver.set_target_eigenvalue(this->target);
718 *
EigenSolver<dim>::solve_problem()
728 *
this->mpi_communicator);
730 *
initialize_eigensolver(eigensolver);
737 *
eigensolver.solve(stiffness_matrix,
741 *
eigenfunctions->size());
742 *
for (
auto &entry : *eigenfunctions)
744 *
constraints.distribute(entry);
746 *
convert_solution();
748 *
return solver_control.last_step();
753 *
EigenSolver<dim>::n_dofs() const
755 *
return dof_handler.n_dofs();
765 *
EigenSolver<dim>::setup_system()
767 *
dof_handler.distribute_dofs(*fe_collection);
768 *
constraints.clear();
771 *
constraints.close();
773 *
eigenfunctions->resize(this->n_eigenpairs);
776 *
IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs();
778 *
for (
auto &entry : *eigenfunctions)
780 *
entry.
reinit(eigenfunction_index_set, MPI_COMM_WORLD);
789 *
EigenSolver<dim>::assemble_system()
792 *
*quadrature_collection,
798 * Prep the system matrices
for the solution
801 *
stiffness_matrix.reinit(dof_handler.n_dofs(),
802 *
dof_handler.n_dofs(),
803 *
dof_handler.max_couplings_between_dofs());
805 *
dof_handler.n_dofs(),
806 *
dof_handler.max_couplings_between_dofs());
809 *
std::vector<types::global_dof_index> local_dof_indices;
811 *
for (
const auto &cell : dof_handler.active_cell_iterators())
813 *
const unsigned
int dofs_per_cell = cell->get_fe().dofs_per_cell;
815 *
cell_stiffness_matrix.reinit(dofs_per_cell, dofs_per_cell);
816 *
cell_stiffness_matrix = 0;
818 *
cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
819 *
cell_mass_matrix = 0;
821 *
hp_fe_values.reinit(cell);
825 *
for (
unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
828 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
830 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
834 * Note that (in general) the Nedelec element is not
835 * primitive, namely that the shape
functions are vectorial
836 * with components in more than
one direction
842 *
cell_stiffness_matrix(i, j) +=
843 *
Operations::curlcurl(fe_values, i, j, q_point) *
844 *
fe_values.JxW(q_point);
846 *
cell_mass_matrix(i, j) +=
847 *
(Operations::dot_term(fe_values, i, j, q_point)) *
848 *
fe_values.JxW(q_point);
851 *
local_dof_indices.resize(dofs_per_cell);
852 *
cell->get_dof_indices(local_dof_indices);
855 *
constraints.distribute_local_to_global(cell_stiffness_matrix,
858 *
constraints.distribute_local_to_global(cell_mass_matrix,
865 *
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
866 *
if (constraints.is_constrained(i))
868 *
stiffness_matrix.set(i, i, 10000.0);
873 * since we have just
set individual elements, we need the following
885 *
class PrimalSolver :
public EigenSolver<dim>
888 *
PrimalSolver(
const std::string & prm_file,
890 *
const unsigned int &min_degree,
891 *
const unsigned int &max_degree,
892 *
const unsigned int &starting_degree);
897 *
virtual unsigned int
898 *
n_dofs()
const override;
902 *
PrimalSolver<dim>::PrimalSolver(
const std::string & prm_file,
904 *
const unsigned int &min_degree,
905 *
const unsigned int &max_degree,
906 *
const unsigned int &starting_degree)
907 *
: Base<dim>(prm_file, triangulation)
908 *
, EigenSolver<dim>(prm_file,
921 *
PrimalSolver<dim>::output_solution()
925 *
Vector<double> fe_degrees(this->triangulation->n_active_cells());
926 *
for (
const auto &cell : this->dof_handler.active_cell_iterators())
927 *
fe_degrees(cell->active_cell_index()) =
928 *
(*this->fe_collection)[cell->active_fe_index()].degree;
929 *
data_out.add_data_vector(fe_degrees,
"fe_degree");
930 *
data_out.add_data_vector((*this->eigenfunctions)[0],
931 *
std::string(
"eigenfunction_no_") +
934 *
std::cout <<
"Eigenvalue: " << (*this->
eigenvalues)[0]
935 * <<
" NDoFs: " << this->dof_handler.n_dofs() << std::endl;
936 *
std::ofstream eigenvalues_out(
937 *
"eigenvalues-" + std::to_string(this->refinement_cycle) +
".txt");
939 *
eigenvalues_out << std::setprecision(20) << (*this->
eigenvalues)[0] <<
" "
940 * << this->dof_handler.n_dofs() << std::endl;
942 *
eigenvalues_out.close();
945 *
data_out.build_patches();
946 *
std::ofstream output(
"eigenvectors-" +
947 *
std::to_string(this->refinement_cycle) +
".vtu");
948 *
data_out.write_vtu(output);
953 *
PrimalSolver<dim>::n_dofs() const
955 *
return EigenSolver<dim>::n_dofs();
960 * Note, that at least
for the demonstrated problem (i.e., a Hermitian problem
961 * and eigenvalue QoI), the dual problem is identical to the primal problem;
962 * however, it is convenient to separate them in
this manner (
e.g.,
for
963 * considering functionals of the eigenfunction).
967 *
class DualSolver :
public EigenSolver<dim>
970 *
DualSolver(
const std::string & prm_file,
972 *
const unsigned int &min_degree,
973 *
const unsigned int &max_degree,
974 *
const unsigned int &starting_degree);
978 *
DualSolver<dim>::DualSolver(
const std::string & prm_file,
980 *
const unsigned int &min_degree,
981 *
const unsigned int &max_degree,
982 *
const unsigned int &starting_degree)
983 *
: Base<dim>(prm_file, triangulation)
984 *
, EigenSolver<dim>(prm_file,
996 *
namespace ErrorIndicators
998 *
using namespace Maxwell;
1005 *
template <
int dim,
bool report_dual>
1006 *
class DualWeightedResidual :
public PrimalSolver<dim>,
public DualSolver<dim>
1010 *
output_eigenvalue_data(std::ofstream &os);
1012 *
output_qoi_error_estimates(std::ofstream &os);
1019 *
DualWeightedResidual(
const std::string & prm_file,
1021 *
const unsigned int &min_primal_degree,
1022 *
const unsigned int &max_primal_degree,
1023 *
const unsigned int &starting_primal_degree);
1025 *
virtual unsigned int
1026 *
solve_problem()
override;
1029 *
output_solution()
override;
1031 *
virtual unsigned int
1032 *
n_dofs()
const override;
1041 *
get_primal_DoFHandler();
1044 *
get_dual_DoFHandler();
1047 *
get_FECollection();
1050 *
get_primal_FECollection();
1052 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1053 *
get_eigenfunctions();
1055 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1056 *
get_primal_eigenfunctions();
1058 *
std::unique_ptr<std::vector<double>> &
1059 *
get_primal_eigenvalues();
1061 *
std::unique_ptr<std::vector<double>> &
1062 *
get_dual_eigenvalues();
1065 *
synchronize_discretization();
1070 *
return PrimalSolver<dim>::fe_collection->max_degree();
1072 *
double qoi_error_estimate = 0;
1093 *
std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values;
1094 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values;
1095 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor;
1096 *
std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values;
1098 *
std::unique_ptr<hp::FEValues<dim>> cell_hp_fe_values_forward;
1099 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_forward;
1100 *
std::unique_ptr<hp::FEFaceValues<dim>> face_hp_fe_values_neighbor_forward;
1101 *
std::unique_ptr<hp::FESubfaceValues<dim>> subface_hp_fe_values_forward;
1102 *
using FaceIntegrals =
1103 *
typename std::map<typename DoFHandler<dim>::face_iterator,
double>;
1106 *
solve_primal_problem();
1109 *
solve_dual_problem();
1120 *
initialize_error_estimation_data();
1123 *
estimate_on_one_cell(
1127 *
const double & lambda_h,
1129 *
FaceIntegrals & face_integrals);
1132 *
integrate_over_cell(
1136 *
const double & lambda_h,
1140 *
integrate_over_regular_face(
1142 *
const unsigned int & face_no,
1145 *
FaceIntegrals & face_integrals);
1148 *
integrate_over_irregular_face(
1150 *
const unsigned int & face_no,
1153 *
FaceIntegrals & face_integrals);
1160 *
template <
int dim,
bool report_dual>
1161 *
DualWeightedResidual<dim, report_dual>::DualWeightedResidual(
1162 *
const std::string & prm_file,
1164 *
const unsigned int &min_primal_degree,
1165 *
const unsigned int &max_primal_degree,
1166 *
const unsigned int &starting_primal_degree)
1167 *
: Base<dim>(prm_file, triangulation)
1168 *
, PrimalSolver<dim>(prm_file,
1170 *
min_primal_degree,
1171 *
max_primal_degree,
1172 *
starting_primal_degree)
1173 *
, DualSolver<dim>(prm_file,
1175 *
min_primal_degree + 1,
1176 *
max_primal_degree + 1,
1177 *
starting_primal_degree + 1)
1179 *
initialize_error_estimation_data();
1186 *
template <
int dim,
bool report_dual>
1188 *
DualWeightedResidual<dim, report_dual>::get_DoFHandler()
1191 *
return &(PrimalSolver<dim>::dof_handler);
1193 *
return &(DualSolver<dim>::dof_handler);
1198 * See above function, but to specifically output the primal
DoFHandler...
1201 *
template <
int dim,
bool report_dual>
1203 *
DualWeightedResidual<dim, report_dual>::get_primal_DoFHandler()
1205 *
return &(PrimalSolver<dim>::dof_handler);
1210 * See above function, but
for the FECollection
1213 *
template <
int dim,
bool report_dual>
1215 *
DualWeightedResidual<dim, report_dual>::get_FECollection()
1218 *
return &*(PrimalSolver<dim>::fe_collection);
1220 *
return &*(DualSolver<dim>::fe_collection);
1225 * See above function, but
for the primal FECollection
1228 *
template <
int dim,
bool report_dual>
1230 *
DualWeightedResidual<dim, report_dual>::get_primal_FECollection()
1232 *
return &*(PrimalSolver<dim>::fe_collection);
1235 *
template <
int dim,
bool report_dual>
1237 *
DualWeightedResidual<dim, report_dual>::get_dual_DoFHandler()
1239 *
return &(DualSolver<dim>::dof_handler);
1243 *
template <
int dim,
bool report_dual>
1244 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1245 *
DualWeightedResidual<dim, report_dual>::get_eigenfunctions()
1248 *
return (PrimalSolver<dim>::eigenfunctions);
1250 *
return (DualSolver<dim>::eigenfunctions);
1254 *
template <
int dim,
bool report_dual>
1255 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1256 *
DualWeightedResidual<dim, report_dual>::get_primal_eigenfunctions()
1258 *
return (PrimalSolver<dim>::eigenfunctions);
1262 *
template <
int dim,
bool report_dual>
1263 *
std::unique_ptr<std::vector<double>> &
1264 *
DualWeightedResidual<dim, report_dual>::get_primal_eigenvalues()
1266 *
return PrimalSolver<dim>::eigenvalues;
1270 *
template <
int dim,
bool report_dual>
1271 *
std::unique_ptr<std::vector<double>> &
1272 *
DualWeightedResidual<dim, report_dual>::get_dual_eigenvalues()
1274 *
return DualSolver<dim>::eigenvalues;
1277 *
template <
int dim,
bool report_dual>
1279 *
DualWeightedResidual<dim, report_dual>::output_solution()
1281 *
PrimalSolver<dim>::output_solution();
1286 * Solves the primal problem
1289 *
template <
int dim,
bool report_dual>
1291 *
DualWeightedResidual<dim, report_dual>::solve_primal_problem()
1293 *
return PrimalSolver<dim>::solve_problem();
1298 * Solves the dual problem
1301 *
template <
int dim,
bool report_dual>
1303 *
DualWeightedResidual<dim, report_dual>::solve_dual_problem()
1305 *
return DualSolver<dim>::solve_problem();
1312 *
template <
int dim,
bool report_dual>
1314 *
DualWeightedResidual<dim, report_dual>::solve_problem()
1316 *
DualWeightedResidual<dim, report_dual>::solve_primal_problem();
1317 *
return DualWeightedResidual<dim, report_dual>::solve_dual_problem();
1323 *
template <
int dim,
bool report_dual>
1325 *
DualWeightedResidual<dim, report_dual>::n_dofs() const
1327 *
return PrimalSolver<dim>::n_dofs();
1336 *
template <
int dim,
bool report_dual>
1338 *
DualWeightedResidual<dim, report_dual>::synchronize_discretization()
1352 * In
this case, we have modified the polynomial orders
for the dual;
1353 * need to update the primal
1356 *
dof1 = &(DualSolver<dim>::dof_handler);
1357 *
dof2 = &(PrimalSolver<dim>::dof_handler);
1360 *
endc1 = dof1->end();
1362 *
for (; cell1 < endc1; ++cell1, ++cell2)
1364 *
cell2->set_active_fe_index(cell1->active_fe_index());
1372 *
template <
int dim,
bool report_dual>
1374 *
DualWeightedResidual<dim, report_dual>::initialize_error_estimation_data()
1378 * initialize the cell fe_values...
1381 *
cell_hp_fe_values = std::make_unique<hp::FEValues<dim>>(
1382 *
*DualSolver<dim>::fe_collection,
1383 *
*DualSolver<dim>::quadrature_collection,
1386 *
face_hp_fe_values = std::make_unique<hp::FEFaceValues<dim>>(
1387 *
*DualSolver<dim>::fe_collection,
1388 *
*DualSolver<dim>::face_quadrature_collection,
1391 *
face_hp_fe_values_neighbor = std::make_unique<hp::FEFaceValues<dim>>(
1392 *
*DualSolver<dim>::fe_collection,
1393 *
*DualSolver<dim>::face_quadrature_collection,
1396 *
subface_hp_fe_values = std::make_unique<hp::FESubfaceValues<dim>>(
1397 *
*DualSolver<dim>::fe_collection,
1398 *
*DualSolver<dim>::face_quadrature_collection,
1407 *
template <
int dim,
bool report_dual>
1409 *
DualWeightedResidual<dim, report_dual>::normalize_solutions(
1413 *
double sum_primal = 0.0, sum_dual = 0.0;
1414 *
for (
const auto &cell :
1415 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1417 *
cell_hp_fe_values->
reinit(cell);
1421 * grab the fe_values
object
1427 *
std::vector<Vector<double>> cell_primal_values(
1429 *
cell_dual_values(fe_values.n_quadrature_points,
Vector<double>(dim));
1430 *
fe_values.get_function_values(primal_solution, cell_primal_values);
1431 *
fe_values.get_function_values(dual_weights, cell_dual_values);
1434 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1437 *
cell_primal_values[p] * cell_primal_values[p] * fe_values.JxW(p);
1439 *
cell_dual_values[p] * cell_dual_values[p] * fe_values.JxW(p);
1443 *
primal_solution /=
sqrt(sum_primal);
1444 *
dual_weights /=
sqrt(sum_dual);
1451 *
template <
int dim,
bool report_dual>
1453 *
DualWeightedResidual<dim, report_dual>::estimate_error(
1458 * The constraints could be grabbed directly, but
this is simple
1463 *
primal_hanging_node_constraints);
1464 *
primal_hanging_node_constraints.close();
1468 *
dual_hanging_node_constraints);
1469 *
dual_hanging_node_constraints.close();
1473 * First map the primal solution to the space of the dual solution
1474 * This allows us to use just
one set of
FEValues objects (rather than one
1475 *
set for the primal, one
for dual)
1481 *
Vector<double> primal_solution(DualSolver<dim>::dof_handler.n_dofs());
1483 *
embed(PrimalSolver<dim>::dof_handler,
1484 *
DualSolver<dim>::dof_handler,
1485 *
dual_hanging_node_constraints,
1486 *
*(PrimalSolver<dim>::get_solution()),
1489 *
Vector<double> &dual_solution = *(DualSolver<dim>::get_solution());
1491 *
normalize_solutions(primal_solution, dual_solution);
1493 *
Vector<double> dual_weights(DualSolver<dim>::dof_handler.n_dofs()),
1494 *
dual_weights_interm(PrimalSolver<dim>::dof_handler.n_dofs());
1498 * First
extract the dual solution to the space of the primal
1501 *
extract(DualSolver<dim>::dof_handler,
1502 *
PrimalSolver<dim>::dof_handler,
1503 *
primal_hanging_node_constraints,
1504 *
*(DualSolver<dim>::get_solution()),
1505 *
dual_weights_interm);
1509 * Now embed
this back to the space of the dual solution
1512 *
embed(PrimalSolver<dim>::dof_handler,
1513 *
DualSolver<dim>::dof_handler,
1514 *
dual_hanging_node_constraints,
1515 *
dual_weights_interm,
1521 * Subtract
this from the full dual solution
1524 *
dual_weights -= *(DualSolver<dim>::get_solution());
1525 *
dual_weights *= -1.0;
1527 *
*(DualSolver<dim>::get_solution()) -= primal_solution;
1529 *
FaceIntegrals face_integrals;
1530 *
for (
const auto &cell :
1531 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1532 *
for (
const auto &face : cell->face_iterators())
1533 *
face_integrals[face] = -1e20;
1536 *
for (
const auto &cell :
1537 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1539 *
estimate_on_one_cell(cell,
1542 *
*(PrimalSolver<dim>::get_lambda_h()),
1546 *
unsigned int present_cell = 0;
1547 *
for (
const auto &cell :
1548 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1550 *
for (
const auto &face : cell->face_iterators())
1552 *
Assert(face_integrals.find(face) != face_integrals.
end(),
1553 *
ExcInternalError());
1554 *
error_indicators(present_cell) -= 0.5 * face_integrals[face];
1561 * Now, with the error indicators computed, let us produce the
1562 * estimate of the QoI error
1565 *
this->qoi_error_estimate =
1566 *
this->get_global_QoI_error(*(DualSolver<dim>::get_solution()),
1567 *
error_indicators);
1568 *
std::cout <<
"Estimated QoI error: " << std::setprecision(20)
1569 *
<< qoi_error_estimate << std::endl;
1576 *
template <
int dim,
bool report_dual>
1578 *
DualWeightedResidual<dim, report_dual>::estimate_on_one_cell(
1582 *
const double & lambda_h,
1584 *
FaceIntegrals & face_integrals)
1586 *
integrate_over_cell(
1587 *
cell, primal_solution, dual_weights, lambda_h, error_indicators);
1588 *
for (
unsigned int face_no :
GeometryInfo<dim>::face_indices())
1590 *
if (cell->face(face_no)->at_boundary())
1592 *
face_integrals[cell->face(face_no)] = 0.0;
1595 *
if ((cell->neighbor(face_no)->has_children() ==
false) &&
1596 *
(cell->neighbor(face_no)->level() == cell->level()) &&
1597 *
(cell->neighbor(face_no)->index() < cell->index()))
1599 *
if (cell->at_boundary(face_no) ==
false)
1600 *
if (cell->neighbor(face_no)->level() < cell->level())
1602 *
if (cell->face(face_no)->has_children() ==
false)
1603 *
integrate_over_regular_face(
1604 *
cell, face_no, primal_solution, dual_weights, face_integrals);
1606 *
integrate_over_irregular_face(
1607 *
cell, face_no, primal_solution, dual_weights, face_integrals);
1614 *
template <
int dim,
bool report_dual>
1616 *
DualWeightedResidual<dim, report_dual>::integrate_over_cell(
1620 *
const double & lambda_h,
1623 *
cell_hp_fe_values->reinit(cell);
1626 * Grab the fe_values
object
1630 *
std::vector<std::vector<Tensor<2, dim, double>>> cell_hessians(
1632 *
std::vector<Vector<double>> cell_primal_values(
1634 *
cell_dual_values(fe_values.n_quadrature_points,
Vector<double>(dim));
1635 *
fe_values.get_function_values(primal_solution, cell_primal_values);
1636 *
fe_values.get_function_hessians(primal_solution, cell_hessians);
1637 *
fe_values.get_function_values(dual_weights, cell_dual_values);
1642 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1645 *
( (cell_hessians[p][1][1][0] -
1646 *
cell_hessians[p][0][1][1]) *
1647 *
(cell_dual_values[p](0)) +
1649 *
(cell_hessians[p][0][0][1] - cell_hessians[p][1][0][0]) *
1650 *
(cell_dual_values[p](1)) -
1651 *
lambda_h * (cell_primal_values[p](0) * cell_dual_values[p](0) +
1652 *
cell_primal_values[p](1) * cell_dual_values[p](1))) *
1656 *
error_indicators(cell->active_cell_index()) +=
sum;
1662 *
template <
int dim,
bool report_dual>
1664 *
DualWeightedResidual<dim, report_dual>::integrate_over_regular_face(
1666 *
const unsigned int & face_no,
1669 *
FaceIntegrals & face_integrals)
1672 *
ExcInternalError());
1673 *
const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1674 *
const auto neighbor = cell->neighbor(face_no);
1676 *
const unsigned int quadrature_index =
1677 *
std::max(cell->active_fe_index(), neighbor->active_fe_index());
1678 *
face_hp_fe_values->reinit(cell, face_no, quadrature_index);
1681 *
std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1682 *
fe_face_values_cell.n_quadrature_points,
1684 *
neighbor_primal_grads(fe_face_values_cell.n_quadrature_points,
1686 *
fe_face_values_cell.get_function_gradients(primal_solution,
1687 *
cell_primal_grads);
1689 *
face_hp_fe_values_neighbor->reinit(neighbor,
1690 *
neighbor_neighbor,
1691 *
quadrature_index);
1695 *
neighbor_primal_grads);
1696 *
const unsigned int n_q_points = fe_face_values_cell.n_quadrature_points;
1697 *
double face_integral = 0.0;
1698 *
std::vector<Vector<double>> cell_dual_values(n_q_points,
1700 *
fe_face_values_cell.get_function_values(dual_weights, cell_dual_values);
1701 *
for (
unsigned int p = 0; p < n_q_points; ++p)
1703 *
auto face_normal = fe_face_values_cell.normal_vector(p);
1706 *
(cell_primal_grads[p][1][0] - cell_primal_grads[p][0][1] -
1707 *
neighbor_primal_grads[p][1][0] + neighbor_primal_grads[p][0][1]) *
1708 *
(cell_dual_values[p][0] * face_normal[1] -
1709 *
cell_dual_values[p][1] * face_normal[0]) *
1710 *
fe_face_values_cell.JxW(p);
1712 *
Assert(face_integrals.find(cell->face(face_no)) != face_integrals.end(),
1713 *
ExcInternalError());
1714 *
Assert(face_integrals[cell->face(face_no)] == -1e20, ExcInternalError());
1715 *
face_integrals[cell->face(face_no)] = face_integral;
1721 *
template <
int dim,
bool report_dual>
1723 *
DualWeightedResidual<dim, report_dual>::integrate_over_irregular_face(
1725 *
const unsigned int & face_no,
1728 *
FaceIntegrals & face_integrals)
1732 *
cell->neighbor(face_no);
1735 *
Assert(neighbor->has_children(), ExcInternalError());
1737 *
const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor(face_no);
1738 *
for (
unsigned int subface_no = 0; subface_no < face->n_children();
1742 *
cell->neighbor_child_on_subface(face_no, subface_no);
1743 *
Assert(neighbor_child->face(neighbor_neighbor) ==
1744 *
cell->face(face_no)->child(subface_no),
1745 *
ExcInternalError());
1746 *
const unsigned int quadrature_index =
1747 *
std::max(cell->active_fe_index(), neighbor_child->active_fe_index());
1750 * initialize fe_subface values_cell
1753 *
subface_hp_fe_values->reinit(cell,
1756 *
quadrature_index);
1759 *
std::vector<std::vector<Tensor<1, dim, double>>> cell_primal_grads(
1760 *
subface_fe_values_cell.n_quadrature_points,
1762 *
neighbor_primal_grads(subface_fe_values_cell.n_quadrature_points,
1764 *
subface_fe_values_cell.get_function_gradients(primal_solution,
1765 *
cell_primal_grads);
1768 * initialize fe_face_values_neighbor
1771 *
face_hp_fe_values_neighbor->reinit(neighbor_child,
1772 *
neighbor_neighbor,
1773 *
quadrature_index);
1777 *
neighbor_primal_grads);
1778 *
const unsigned int n_q_points =
1779 *
subface_fe_values_cell.n_quadrature_points;
1780 *
std::vector<Vector<double>> cell_dual_values(n_q_points,
1782 *
face_fe_values_neighbor.get_function_values(dual_weights,
1783 *
cell_dual_values);
1785 *
double face_integral = 0.0;
1787 *
for (
unsigned int p = 0; p < n_q_points; ++p)
1789 *
auto face_normal = face_fe_values_neighbor.normal_vector(p);
1791 *
(cell_primal_grads[p][0][1] - cell_primal_grads[p][1][0] +
1792 *
neighbor_primal_grads[p][1][0] -
1793 *
neighbor_primal_grads[p][0][1]) *
1794 *
(cell_dual_values[p][0] * face_normal[1] -
1795 *
cell_dual_values[p][1] * face_normal[0]) *
1796 *
face_fe_values_neighbor.JxW(p);
1798 *
face_integrals[neighbor_child->face(neighbor_neighbor)] = face_integral;
1801 *
for (
unsigned int subface_no = 0; subface_no < face->n_children();
1804 *
Assert(face_integrals.find(face->child(subface_no)) !=
1805 *
face_integrals.end(),
1806 *
ExcInternalError());
1807 *
Assert(face_integrals[face->child(subface_no)] != -1e20,
1808 *
ExcInternalError());
1809 *
sum += face_integrals[face->child(subface_no)];
1811 *
face_integrals[face] =
sum;
1814 *
template <
int dim,
bool report_dual>
1816 *
DualWeightedResidual<dim, report_dual>::get_global_QoI_error(
1820 *
auto dual_less_primal =
1824 *
double scaling_factor = 0.0;
1825 *
for (
const auto &cell :
1826 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1828 *
cell_hp_fe_values->
reinit(cell);
1831 * grab the fe_values
object
1837 *
std::vector<Vector<double>> cell_values(fe_values.n_quadrature_points,
1839 *
fe_values.get_function_values(dual_less_primal, cell_values);
1841 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1844 *
(cell_values[p] * cell_values[p]) * fe_values.JxW(p);
1847 *
double global_QoI_error = 0.0;
1848 *
for (
const auto &indicator : error_indicators)
1850 *
global_QoI_error += indicator;
1853 *
global_QoI_error /= (1 - 0.5 * scaling_factor);
1854 *
return global_QoI_error;
1858 *
template <
int dim,
bool report_dual>
1860 *
DualWeightedResidual<dim, report_dual>::embed(
1867 *
assert(u2.size() == dof2.n_dofs() &&
"Incorrect input vector size!");
1872 *
endc1 = dof1.end();
1875 *
for (; cell1 < endc1; ++cell1, ++cell2)
1882 *
assert(fe1.degree < fe2.degree &&
"Incorrect usage of embed!");
1886 * Get the embedding_dofs
1895 *
std::vector<unsigned int> embedding_dofs =
1896 *
fe2.get_embedding_dofs(fe1.degree);
1897 *
const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1903 *
local_dof_values_1.reinit(fe1.dofs_per_cell);
1904 *
cell1->get_dof_values(solution, local_dof_values_1);
1906 *
for (
unsigned int i = 0; i < local_dof_values_1.size(); ++i)
1907 *
local_dof_values_2[embedding_dofs[i]] = local_dof_values_1[i];
1911 * Now
set this changes to the global vector
1914 *
cell2->set_dof_values(local_dof_values_2, u2);
1920 * Applies the constraints of the target finite element space
1923 *
constraints.distribute(u2);
1926 *
template <
int dim,
bool report_dual>
1928 *
DualWeightedResidual<dim, report_dual>::extract(
1937 * Maps from fe1 to fe2
1940 *
assert(u2.size() == dof2.n_dofs() &&
"Incorrect input vector size!");
1945 *
endc1 = dof1.end();
1948 *
for (; cell1 < endc1; ++cell1, ++cell2)
1955 *
assert(fe1.degree > fe2.degree &&
"Incorrect usage of extract!");
1959 * Get the embedding_dofs
1962 *
std::vector<unsigned int> embedding_dofs =
1963 *
fe1.get_embedding_dofs(fe2.degree);
1964 *
const unsigned int dofs_per_cell2 = fe2.n_dofs_per_cell();
1970 *
local_dof_values_1.reinit(fe1.dofs_per_cell);
1971 *
cell1->get_dof_values(solution, local_dof_values_1);
1973 *
for (
unsigned int i = 0; i < local_dof_values_2.size(); ++i)
1974 *
local_dof_values_2[i] = local_dof_values_1[embedding_dofs[i]];
1978 * Now
set this changes to the global vector
1981 *
cell2->set_dof_values(local_dof_values_2, u2);
1987 * Applies the constraints of the target finite element space
1990 *
constraints.distribute(u2);
1992 *
template <
int dim,
bool report_dual>
1994 *
DualWeightedResidual<dim, report_dual>::output_eigenvalue_data(
1995 *
std::ofstream &os)
1997 *
os << (*this->get_primal_eigenvalues())[0] <<
" "
1998 *
<< (this->get_primal_DoFHandler())->n_dofs() <<
" "
1999 *
<< (*this->get_dual_eigenvalues())[0] <<
" "
2000 *
<< (this->get_dual_DoFHandler())->n_dofs() << std::endl;
2002 *
template <
int dim,
bool report_dual>
2004 *
DualWeightedResidual<dim, report_dual>::output_qoi_error_estimates(
2005 *
std::ofstream &os)
2007 *
os << qoi_error_estimate << std::endl;
2014 *
template <
int dim>
2015 *
class KellyErrorIndicator :
public PrimalSolver<dim>
2024 *
output_eigenvalue_data(std::ofstream &os);
2026 *
output_qoi_error_estimates(std::ofstream &);
2027 *
KellyErrorIndicator(
const std::string & prm_file,
2029 *
const unsigned int &min_degree,
2030 *
const unsigned int &max_degree,
2031 *
const unsigned int &starting_degree);
2033 *
virtual unsigned int
2034 *
solve_problem()
override;
2037 *
output_solution()
override;
2040 *
get_FECollection();
2043 *
get_primal_FECollection();
2045 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2046 *
get_eigenfunctions();
2048 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2049 *
get_primal_eigenfunctions();
2051 *
std::unique_ptr<std::vector<double>> &
2052 *
get_primal_eigenvalues();
2056 *
synchronize_discretization();
2062 *
get_primal_DoFHandler();
2067 *
return PrimalSolver<dim>::fe_collection->max_degree();
2069 *
double qoi_error_estimate = 0;
2077 *
prune_eigenpairs(
const double &TOL);
2079 *
std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
2080 *
std::vector<const double *> eigenvalue_ptrs;
2082 *
std::vector<std::shared_ptr<Vector<float>>> errors;
2085 *
template <
int dim>
2086 *
KellyErrorIndicator<dim>::KellyErrorIndicator(
2087 *
const std::string & prm_file,
2089 *
const unsigned int &min_degree,
2090 *
const unsigned int &max_degree,
2091 *
const unsigned int &starting_degree)
2092 *
: Base<dim>(prm_file, coarse_grid)
2093 *
, PrimalSolver<dim>(prm_file,
2100 *
template <
int dim>
2102 *
KellyErrorIndicator<dim>::solve_problem()
2104 *
return PrimalSolver<dim>::solve_problem();
2107 *
template <
int dim>
2109 *
KellyErrorIndicator<dim>::get_FECollection()
2111 *
return &*(PrimalSolver<dim>::fe_collection);
2114 *
template <
int dim>
2116 *
KellyErrorIndicator<dim>::get_primal_FECollection()
2118 *
return &*(PrimalSolver<dim>::fe_collection);
2121 *
template <
int dim>
2122 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2123 *
KellyErrorIndicator<dim>::get_eigenfunctions()
2125 *
return (PrimalSolver<dim>::eigenfunctions);
2128 *
template <
int dim>
2129 *
std::unique_ptr<std::vector<double>> &
2130 *
KellyErrorIndicator<dim>::get_primal_eigenvalues()
2132 *
return PrimalSolver<dim>::eigenvalues;
2135 *
template <
int dim>
2136 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2137 *
KellyErrorIndicator<dim>::get_primal_eigenfunctions()
2139 *
return (PrimalSolver<dim>::eigenfunctions);
2142 *
template <
int dim>
2144 *
KellyErrorIndicator<dim>::get_DoFHandler()
2146 *
return &(PrimalSolver<dim>::dof_handler);
2149 *
template <
int dim>
2151 *
KellyErrorIndicator<dim>::get_primal_DoFHandler()
2153 *
return &(PrimalSolver<dim>::dof_handler);
2156 *
template <
int dim>
2158 *
KellyErrorIndicator<dim>::synchronize_discretization()
2162 * This function does
nothing for this error indicator
2168 *
template <
int dim>
2170 *
KellyErrorIndicator<dim>::output_solution()
2172 *
PrimalSolver<dim>::output_solution();
2175 *
template <
int dim>
2177 *
KellyErrorIndicator<dim>::prune_eigenpairs(
const double &TOL)
2179 *
unsigned int count = 0;
2180 *
for (
size_t eigenpair_index = 0;
2181 *
eigenpair_index < this->eigenfunctions->size();
2182 *
++eigenpair_index)
2184 *
if (count >= this->n_eigenpairs)
2186 *
if (
abs((*this->eigenvalues)[eigenpair_index]) < TOL)
2189 *
eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
2190 *
eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
2194 *
template <
int dim>
2196 *
KellyErrorIndicator<dim>::estimate_error(
Vector<double> &error_indicators)
2198 *
std::cout <<
"Marking cells via Kelly indicator..." << std::endl;
2199 *
prune_eigenpairs(1e-9);
2202 * deallocate the errors vector
2206 *
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2208 *
errors.emplace_back(
2209 *
new Vector<float>(this->triangulation->n_active_cells()));
2211 *
std::vector<Vector<float> *> estimated_error_per_cell(
2212 *
eigenfunction_ptrs.size());
2213 *
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2215 *
estimated_error_per_cell[i] = errors[i].get();
2221 *
*this->face_quadrature_collection,
2226 *
for (
auto &error_vec : errors)
2228 *
auto normalized_vec = *error_vec;
2229 *
normalized_vec /= normalized_vec.l1_norm();
2231 *
for (
unsigned int i = 0; i < error_indicators.size(); ++i)
2232 *
error_indicators(i) += double(normalized_vec(i));
2234 *
std::cout <<
"...Done!" << std::endl;
2236 *
template <
int dim>
2238 *
KellyErrorIndicator<dim>::output_eigenvalue_data(std::ofstream &os)
2240 *
os << (*this->get_primal_eigenvalues())[0] <<
" "
2241 *
<< (this->get_primal_DoFHandler())->n_dofs() << std::endl;
2243 *
template <
int dim>
2245 *
KellyErrorIndicator<dim>::output_qoi_error_estimates(std::ofstream &)
2255 *
namespace RegularityIndicators
2257 *
using namespace dealii;
2261 *
template <
int dim>
2262 *
class LegendreInfo
2266 *
class LegendreInfo<2>
2269 *
std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
2277 *
assert(fe_collection !=
nullptr && dof_handler !=
nullptr &&
2278 *
"A valid FECollection and DoFHandler must be accessible!");
2280 *
legendre_u = std::make_unique<FESeries::Legendre<2>>(
2282 *
legendre_v = std::make_unique<FESeries::Legendre<2>>(
2285 *
legendre_u->precalculate_all_transformation_matrices();
2286 *
legendre_v->precalculate_all_transformation_matrices();
2289 *
template <
class VectorType>
2291 *
compute_coefficient_decay(
const VectorType & eigenfunction,
2292 *
std::vector<double> &smoothness_indicators)
2296 * Compute the coefficients
for the u and v components of the solution
2301 *
smoothness_v(smoothness_indicators.size());
2313 *
for (
unsigned int i = 0; i < smoothness_indicators.size(); ++i)
2315 *
smoothness_indicators[i] =
std::min(smoothness_u[i], smoothness_v[i]);
2323 *
template <
int dim>
2324 *
class LegendreIndicator
2332 *
template <
class VectorType>
2334 *
estimate_smoothness(
2335 *
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2336 *
const unsigned int & index_of_goal,
2337 *
std::vector<double> & smoothness_indicators);
2343 *
template <
int dim>
2345 *
LegendreIndicator<dim>::attach_FE_info_and_initialize(
2354 *
template <
int dim>
2355 *
template <
class VectorType>
2357 *
LegendreIndicator<dim>::estimate_smoothness(
2358 *
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2359 *
const unsigned int & index_of_goal,
2360 *
std::vector<double> & smoothness_indicators)
2362 *
this->
legendre.compute_coefficient_decay((*eigenfunctions)[index_of_goal],
2363 *
smoothness_indicators);
2371 *
namespace Refinement
2373 *
using namespace dealii;
2374 *
using namespace Maxwell;
2376 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2377 *
class Refiner :
public ErrorIndicator,
public RegularityIndicator
2380 *
Refiner(
const std::string & prm_file,
2382 *
const unsigned int &min_degree,
2383 *
const unsigned int &max_degree,
2384 *
const unsigned int &starting_degree);
2387 *
execute_refinement(
const double &smoothness_threshold_fraction);
2390 *
output_solution()
override;
2394 *
std::vector<double> smoothness_indicators;
2395 *
std::ofstream eigenvalues_out;
2396 *
std::ofstream error_estimate_out;
2399 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2400 *
Refiner<dim, ErrorIndicator, RegularityIndicator>::Refiner(
2401 *
const std::string & prm_file,
2403 *
const unsigned int &min_degree,
2404 *
const unsigned int &max_degree,
2405 *
const unsigned int &starting_degree)
2406 *
: Base<dim>(prm_file, coarse_grid)
2407 *
, ErrorIndicator(prm_file,
2412 *
, RegularityIndicator()
2414 *
if (ErrorIndicator::name() ==
"DWR")
2416 *
error_estimate_out.open(
"error_estimate.txt");
2417 *
error_estimate_out << std::setprecision(20);
2420 *
eigenvalues_out.open(
"eigenvalues_" + ErrorIndicator::name() +
"_out.txt");
2421 *
eigenvalues_out << std::setprecision(20);
2426 * For generating samples of the curl of the electric field
2429 *
template <
int dim>
2433 *
CurlPostprocessor()
2440 *
std::vector<
Vector<double>> &computed_quantities)
const override
2443 *
computed_quantities.size());
2444 *
for (
unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
2446 *
computed_quantities[p](0) = input_data.solution_gradients[p][1][0] -
2447 *
input_data.solution_gradients[p][0][1];
2459 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2461 *
Refiner<dim, ErrorIndicator, RegularityIndicator>::output_solution()
2463 *
CurlPostprocessor<dim> curl_u;
2466 *
auto & output_dof = *(ErrorIndicator::get_primal_DoFHandler());
2468 *
Vector<double> fe_degrees(this->triangulation->n_active_cells());
2469 *
for (
const auto &cell : output_dof.active_cell_iterators())
2470 *
fe_degrees(cell->active_cell_index()) =
2471 *
(*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
2473 *
data_out.add_data_vector(fe_degrees,
"fe_degree");
2475 *
data_out.add_data_vector(estimated_error_per_cell,
"error");
2476 *
Vector<double> smoothness_out(this->triangulation->n_active_cells());
2477 *
for (
const auto &cell : output_dof.active_cell_iterators())
2479 *
auto i = cell->active_cell_index();
2480 *
if (!cell->refine_flag_set() && !cell->coarsen_flag_set())
2481 *
smoothness_out(i) = -1;
2483 *
smoothness_out(i) = smoothness_indicators[i];
2485 *
data_out.add_data_vector(smoothness_out,
"smoothness");
2486 *
data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2487 *
std::string(
"eigenfunction_no_") +
2489 *
data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2492 *
ErrorIndicator::output_eigenvalue_data(eigenvalues_out);
2493 *
ErrorIndicator::output_qoi_error_estimates(error_estimate_out);
2495 *
std::cout <<
"Number of DoFs: " << (this->get_primal_DoFHandler())->n_dofs()
2499 *
data_out.build_patches();
2500 *
std::ofstream output(
"eigenvectors-" + ErrorIndicator::name() +
"-" +
2501 *
std::to_string(this->refinement_cycle) + +
".vtu");
2502 *
data_out.write_vtu(output);
2511 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2513 *
Refiner<dim, ErrorIndicator, RegularityIndicator>::execute_refinement(
2514 *
const double &smoothness_threshold_fraction)
2518 * First initialize the RegularityIndicator...
2519 * Depending on the limits
set,
this may take a
while
2522 *
std::cout <<
"Initializing RegularityIndicator..." << std::endl;
2524 *
<<
"(This may take a while if the max expansion order is set too high)"
2526 *
RegularityIndicator::attach_FE_info_and_initialize(
2527 *
ErrorIndicator::get_FECollection(), ErrorIndicator::get_DoFHandler());
2528 *
std::cout <<
"Done!" << std::endl <<
"Starting Refinement..." << std::endl;
2530 *
for (
unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
2532 *
this->set_refinement_cycle(cycle);
2533 *
std::cout <<
"Cycle: " << cycle << std::endl;
2534 *
ErrorIndicator::solve_problem();
2535 *
this->estimated_error_per_cell.reinit(
2536 *
this->triangulation->n_active_cells());
2538 *
ErrorIndicator::estimate_error(estimated_error_per_cell);
2542 * Depending on the source of the error estimation/indication, these
2543 *
values might be signed, so we address that with the following
2546 *
for (
double &error_indicator : estimated_error_per_cell)
2547 *
error_indicator =
std::
abs(error_indicator);
2551 *
*this->triangulation, estimated_error_per_cell, 1. / 5., 0.000);
2555 * Now get regularity indicators
2556 * For those elements which must be refined,
swap to increasing @f$p@f$
2557 * depending on the regularity threshold...
2563 *
smoothness_indicators =
2564 *
std::vector<double>(this->triangulation->n_active_cells(),
2565 *
std::numeric_limits<double>::max());
2566 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2567 *
ErrorIndicator::PrimalSolver::max_degree)
2568 *
RegularityIndicator::estimate_smoothness(
2569 *
ErrorIndicator::get_eigenfunctions(), 0, smoothness_indicators);
2575 *
this->output_solution();
2576 *
const double threshold_smoothness = smoothness_threshold_fraction;
2577 *
unsigned int num_refined = 0, num_coarsened = 0;
2578 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2579 *
ErrorIndicator::PrimalSolver::max_degree)
2581 *
for (
const auto &cell :
2582 *
ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2584 *
if (cell->refine_flag_set())
2586 *
if (cell->coarsen_flag_set())
2588 *
if (cell->refine_flag_set() &&
2589 *
smoothness_indicators[cell->active_cell_index()] >
2590 *
threshold_smoothness &&
2591 *
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2592 *
ErrorIndicator::get_FECollection()->size())
2594 *
cell->clear_refine_flag();
2595 *
cell->set_active_fe_index(cell->active_fe_index() + 1);
2597 *
else if (cell->coarsen_flag_set() &&
2598 *
smoothness_indicators[cell->active_cell_index()] <
2599 *
threshold_smoothness &&
2600 *
cell->active_fe_index() != 0)
2602 *
cell->clear_coarsen_flag();
2604 *
cell->set_active_fe_index(cell->active_fe_index() - 1);
2608 * Here we also impose a limit on how small the cells can become
2611 *
else if (cell->refine_flag_set() && cell->diameter() < 5.0e-6)
2613 *
cell->clear_refine_flag();
2614 *
if (
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2615 *
ErrorIndicator::get_FECollection()->size())
2616 *
cell->set_active_fe_index(cell->active_fe_index() + 1);
2623 * Check what the smallest
diameter is
2626 *
double min_diameter = std::numeric_limits<double>::max();
2627 *
for (
const auto &cell :
2628 *
ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2632 *
std::cout <<
"Min diameter: " << min_diameter << std::endl;
2634 *
ErrorIndicator::synchronize_discretization();
2636 *
(this->triangulation)->execute_coarsening_and_refinement();
2642 *
main(
int argc,
char **argv)
2646 *
using namespace dealii;
2647 *
using namespace Maxwell;
2648 *
using namespace Refinement;
2649 *
using namespace ErrorIndicators;
2650 *
using namespace RegularityIndicators;
2658 *
ExcMessage(
"This program can only be run in serial, use ./maxwell-hp"));
2661 *
Structures::create_L_waveguide(triangulation_DWR, 2.0);
2662 *
Structures::create_L_waveguide(triangulation_Kelly, 2.0);
2664 *
Refiner<2, KellyErrorIndicator<2>, LegendreIndicator<2>> problem_Kelly(
2666 *
triangulation_Kelly,
2671 *
Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
2672 *
problem_DWR(
"maxwell-hp.prm",
2673 *
triangulation_DWR,
2680 * The threshold
for the
hp-decision: too small -> not enough
2681 * @f$h@f$-refinement, too large -> not enough @f$p@f$-refinement
2684 *
double smoothness_threshold = 0.75;
2686 *
std::cout <<
"Executing refinement for the Kelly strategy!" << std::endl;
2687 *
problem_Kelly.execute_refinement(smoothness_threshold);
2688 *
std::cout <<
"...Done with Kelly refinement strategy!" << std::endl;
2689 *
std::cout <<
"Executing refinement for the DWR strategy!" << std::endl;
2690 *
problem_DWR.execute_refinement(smoothness_threshold);
2691 *
std::cout <<
"...Done with DWR refinement strategy!" << std::endl;
2694 *
catch (std::exception &exc)
2696 *
std::cerr << std::endl
2698 *
<<
"----------------------------------------------------"
2700 *
std::cerr <<
"Exception on processing: " << std::endl
2701 *
<< exc.what() << std::endl
2702 *
<<
"Aborting!" << std::endl
2703 *
<<
"----------------------------------------------------"
2710 *
std::cerr << std::endl
2712 *
<<
"----------------------------------------------------"
2714 *
std::cerr <<
"Unknown exception!" << std::endl
2715 *
<<
"Aborting!" << std::endl
2716 *
<<
"----------------------------------------------------"
2721 *
std::cout << std::endl <<
" Job done." << std::endl;
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * * struct InterferenceTaperTransform *
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
active_cell_iterator begin_active(const unsigned int level=0) const
const FEFaceValues< dim, spacedim > & get_present_fe_values() const
const FESubfaceValues< dim, spacedim > & get_present_fe_values() const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) 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)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
#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.
std::vector< index_type > data
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)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
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)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
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 > &)
void swap(ObserverPointer< T, P > &t1, ObserverPointer< T, Q > &t2)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
1.9.8