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
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;
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;
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;
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,
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;
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;
563 * convert_solution();
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)
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;
817 * cell_stiffness_matrix.reinit(dofs_per_cell, dofs_per_cell);
818 * cell_stiffness_matrix = 0;
820 * cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
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;
830 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
832 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
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);
853 * local_dof_indices.resize(dofs_per_cell);
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,
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)
910 * , EigenSolver<dim>(prm_file,
923 * PrimalSolver<dim>::output_solution()
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,
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)
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,
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;
1043 * get_primal_DoFHandler();
1046 * get_dual_DoFHandler();
1049 * get_FECollection();
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;
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();
1122 * initialize_error_estimation_data();
1125 * estimate_on_one_cell(
1129 *
const double & lambda_h,
1131 * FaceIntegrals & face_integrals);
1134 * integrate_over_cell(
1138 *
const double & lambda_h,
1142 * integrate_over_regular_face(
1144 *
const unsigned int & face_no,
1147 * FaceIntegrals & face_integrals);
1150 * integrate_over_irregular_face(
1152 *
const unsigned int & face_no,
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)
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
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);
1590 *
for (
unsigned int face_no :
GeometryInfo<dim>::face_indices())
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);
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);
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);
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);
1779 * neighbor_primal_grads);
1780 *
const unsigned int n_q_points =
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
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,
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;
2042 * get_FECollection();
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;
2079 * prune_eigenpairs(
const double &TOL);
2082 * std::vector<const ReadVector<PetscScalar> *> eigenfunction_ptrs;
2084 * std::vector<const PETScWrappers::MPI::Vector *> eigenfunction_ptrs;
2086 * std::vector<const double *> eigenvalue_ptrs;
2088 * std::vector<std::shared_ptr<Vector<float>>> errors;
2091 *
template <
int dim>
2092 * KellyErrorIndicator<dim>::KellyErrorIndicator(
2093 *
const std::string & prm_file,
2095 *
const unsigned int &min_degree,
2096 *
const unsigned int &max_degree,
2097 *
const unsigned int &starting_degree)
2098 * : Base<dim>(prm_file, coarse_grid)
2099 * , PrimalSolver<dim>(prm_file,
2106 *
template <
int dim>
2108 * KellyErrorIndicator<dim>::solve_problem()
2110 *
return PrimalSolver<dim>::solve_problem();
2113 *
template <
int dim>
2115 * KellyErrorIndicator<dim>::get_FECollection()
2117 *
return &*(PrimalSolver<dim>::fe_collection);
2120 *
template <
int dim>
2122 * KellyErrorIndicator<dim>::get_primal_FECollection()
2124 *
return &*(PrimalSolver<dim>::fe_collection);
2127 *
template <
int dim>
2128 * std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2129 * KellyErrorIndicator<dim>::get_eigenfunctions()
2131 *
return (PrimalSolver<dim>::eigenfunctions);
2134 *
template <
int dim>
2135 * std::unique_ptr<std::vector<double>> &
2136 * KellyErrorIndicator<dim>::get_primal_eigenvalues()
2138 *
return PrimalSolver<dim>::eigenvalues;
2141 *
template <
int dim>
2142 * std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2143 * KellyErrorIndicator<dim>::get_primal_eigenfunctions()
2145 *
return (PrimalSolver<dim>::eigenfunctions);
2148 *
template <
int dim>
2150 * KellyErrorIndicator<dim>::get_DoFHandler()
2152 *
return &(PrimalSolver<dim>::dof_handler);
2155 *
template <
int dim>
2157 * KellyErrorIndicator<dim>::get_primal_DoFHandler()
2159 *
return &(PrimalSolver<dim>::dof_handler);
2162 *
template <
int dim>
2164 * KellyErrorIndicator<dim>::synchronize_discretization()
2168 * This function does
nothing for this error indicator
2174 *
template <
int dim>
2176 * KellyErrorIndicator<dim>::output_solution()
2178 * PrimalSolver<dim>::output_solution();
2181 *
template <
int dim>
2183 * KellyErrorIndicator<dim>::prune_eigenpairs(
const double &TOL)
2185 *
unsigned int count = 0;
2186 *
for (
size_t eigenpair_index = 0;
2187 * eigenpair_index < this->eigenfunctions->size();
2188 * ++eigenpair_index)
2190 *
if (count >= this->n_eigenpairs)
2192 *
if (
abs((*this->eigenvalues)[eigenpair_index]) < TOL)
2195 * eigenfunction_ptrs.push_back(&(*this->eigenfunctions)[eigenpair_index]);
2196 * eigenvalue_ptrs.push_back(&(*this->eigenvalues)[eigenpair_index]);
2200 *
template <
int dim>
2202 * KellyErrorIndicator<dim>::estimate_error(
Vector<double> &error_indicators)
2204 * std::cout <<
"Marking cells via Kelly indicator..." << std::endl;
2205 * prune_eigenpairs(1e-9);
2208 * deallocate the errors vector
2212 *
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2214 * errors.emplace_back(
2217 * std::vector<Vector<float> *> estimated_error_per_cell(
2218 * eigenfunction_ptrs.size());
2219 *
for (
size_t i = 0; i < eigenfunction_ptrs.size(); ++i)
2221 * estimated_error_per_cell[i] = errors[i].get();
2228 * *this->face_quadrature_collection,
2234 * *this->face_quadrature_collection,
2236 * eigenfunction_ptrs,
2237 * estimated_error_per_cell);
2240 *
for (
auto &error_vec : errors)
2242 * auto normalized_vec = *error_vec;
2243 * normalized_vec /= normalized_vec.l1_norm();
2245 *
for (
unsigned int i = 0; i < error_indicators.size(); ++i)
2246 * error_indicators(i) += double(normalized_vec(i));
2248 * std::cout <<
"...Done!" << std::endl;
2250 *
template <
int dim>
2252 * KellyErrorIndicator<dim>::output_eigenvalue_data(std::ofstream &os)
2254 * os << (*this->get_primal_eigenvalues())[0] <<
" "
2255 * << (this->get_primal_DoFHandler())->n_dofs() << std::endl;
2257 *
template <
int dim>
2259 * KellyErrorIndicator<dim>::output_qoi_error_estimates(std::ofstream &)
2269 *
namespace RegularityIndicators
2271 *
using namespace dealii;
2275 *
template <
int dim>
2276 *
class LegendreInfo
2280 *
class LegendreInfo<2>
2283 * std::unique_ptr<FESeries::Legendre<2>> legendre_u, legendre_v;
2291 * assert(fe_collection !=
nullptr && dof_handler !=
nullptr &&
2292 *
"A valid FECollection and DoFHandler must be accessible!");
2294 * legendre_u = std::make_unique<FESeries::Legendre<2>>(
2296 * legendre_v = std::make_unique<FESeries::Legendre<2>>(
2299 * legendre_u->precalculate_all_transformation_matrices();
2300 * legendre_v->precalculate_all_transformation_matrices();
2303 *
template <
class VectorType>
2305 * compute_coefficient_decay(
const VectorType & eigenfunction,
2306 * std::vector<double> &smoothness_indicators)
2310 * Compute the coefficients
for the u and v components of the solution
2315 * smoothness_v(smoothness_indicators.size());
2327 *
for (
unsigned int i = 0; i < smoothness_indicators.size(); ++i)
2329 * smoothness_indicators[i] =
std::min(smoothness_u[i], smoothness_v[i]);
2337 *
template <
int dim>
2338 *
class LegendreIndicator
2346 *
template <
class VectorType>
2348 * estimate_smoothness(
2349 *
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2350 *
const unsigned int & index_of_goal,
2351 * std::vector<double> & smoothness_indicators);
2357 *
template <
int dim>
2359 * LegendreIndicator<dim>::attach_FE_info_and_initialize(
2368 *
template <
int dim>
2369 *
template <
class VectorType>
2371 * LegendreIndicator<dim>::estimate_smoothness(
2372 *
const std::unique_ptr<std::vector<VectorType>> &eigenfunctions,
2373 *
const unsigned int & index_of_goal,
2374 * std::vector<double> & smoothness_indicators)
2376 * this->
legendre.compute_coefficient_decay((*eigenfunctions)[index_of_goal],
2377 * smoothness_indicators);
2385 *
namespace Refinement
2387 *
using namespace dealii;
2388 *
using namespace Maxwell;
2390 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2391 *
class Refiner :
public ErrorIndicator,
public RegularityIndicator
2394 * Refiner(
const std::string & prm_file,
2396 *
const unsigned int &min_degree,
2397 *
const unsigned int &max_degree,
2398 *
const unsigned int &starting_degree);
2401 * execute_refinement(
const double &smoothness_threshold_fraction);
2404 * output_solution()
override;
2408 * std::vector<double> smoothness_indicators;
2409 * std::ofstream eigenvalues_out;
2410 * std::ofstream error_estimate_out;
2413 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2414 * Refiner<dim, ErrorIndicator, RegularityIndicator>::Refiner(
2415 *
const std::string & prm_file,
2417 *
const unsigned int &min_degree,
2418 *
const unsigned int &max_degree,
2419 *
const unsigned int &starting_degree)
2420 * : Base<dim>(prm_file, coarse_grid)
2421 * , ErrorIndicator(prm_file,
2426 * , RegularityIndicator()
2428 *
if (ErrorIndicator::name() ==
"DWR")
2430 * error_estimate_out.open(
"error_estimate.txt");
2431 * error_estimate_out << std::setprecision(20);
2434 * eigenvalues_out.open(
"eigenvalues_" + ErrorIndicator::name() +
"_out.txt");
2435 * eigenvalues_out << std::setprecision(20);
2440 * For generating samples of the curl of the electric field
2443 *
template <
int dim>
2447 * CurlPostprocessor()
2454 * std::vector<
Vector<double>> &computed_quantities)
const override
2457 * computed_quantities.size());
2458 *
for (
unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
2460 * computed_quantities[p](0) = input_data.solution_gradients[p][1][0] -
2461 * input_data.solution_gradients[p][0][1];
2473 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2475 * Refiner<dim, ErrorIndicator, RegularityIndicator>::output_solution()
2477 * CurlPostprocessor<dim> curl_u;
2480 *
auto & output_dof = *(ErrorIndicator::get_primal_DoFHandler());
2483 *
for (
const auto &cell : output_dof.active_cell_iterators())
2484 * fe_degrees(cell->active_cell_index()) =
2485 * (*ErrorIndicator::get_primal_FECollection())[cell->active_fe_index()]
2487 * data_out.add_data_vector(fe_degrees,
"fe_degree");
2489 * data_out.add_data_vector(estimated_error_per_cell,
"error");
2491 *
for (
const auto &cell : output_dof.active_cell_iterators())
2493 * auto i = cell->active_cell_index();
2494 *
if (!cell->refine_flag_set() && !cell->coarsen_flag_set())
2495 * smoothness_out(i) = -1;
2497 * smoothness_out(i) = smoothness_indicators[i];
2499 * data_out.add_data_vector(smoothness_out,
"smoothness");
2500 * data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2501 * std::string(
"eigenfunction_no_") +
2503 * data_out.add_data_vector((*ErrorIndicator::get_primal_eigenfunctions())[0],
2506 * ErrorIndicator::output_eigenvalue_data(eigenvalues_out);
2507 * ErrorIndicator::output_qoi_error_estimates(error_estimate_out);
2509 * std::cout <<
"Number of DoFs: " << (this->get_primal_DoFHandler())->n_dofs()
2513 * data_out.build_patches();
2514 * std::ofstream output(
"eigenvectors-" + ErrorIndicator::name() +
"-" +
2515 * std::to_string(this->refinement_cycle) + +
".vtu");
2516 * data_out.write_vtu(output);
2525 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2527 * Refiner<dim, ErrorIndicator, RegularityIndicator>::execute_refinement(
2528 *
const double &smoothness_threshold_fraction)
2532 * First initialize the RegularityIndicator...
2533 * Depending on the limits set,
this may take a
while
2536 * std::cout <<
"Initializing RegularityIndicator..." << std::endl;
2538 * <<
"(This may take a while if the max expansion order is set too high)"
2540 * RegularityIndicator::attach_FE_info_and_initialize(
2541 * ErrorIndicator::get_FECollection(), ErrorIndicator::get_DoFHandler());
2542 * std::cout <<
"Done!" << std::endl <<
"Starting Refinement..." << std::endl;
2544 *
for (
unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
2546 * this->set_refinement_cycle(cycle);
2547 * std::cout <<
"Cycle: " << cycle << std::endl;
2548 * ErrorIndicator::solve_problem();
2549 * this->estimated_error_per_cell.reinit(
2552 * ErrorIndicator::estimate_error(estimated_error_per_cell);
2556 * Depending on the source of the error estimation/indication, these
2557 *
values might be signed, so we address that with the following
2560 *
for (
double &error_indicator : estimated_error_per_cell)
2561 * error_indicator =
std::
abs(error_indicator);
2565 * *this->
triangulation, estimated_error_per_cell, 1. / 5., 0.000);
2569 * Now get regularity indicators
2570 * For those elements which must be refined,
swap to increasing @f$p@f$
2571 * depending on the regularity threshold...
2577 * smoothness_indicators =
2578 * std::vector<double>(this->
triangulation->n_active_cells(),
2579 * std::numeric_limits<double>::max());
2580 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2581 * ErrorIndicator::PrimalSolver::max_degree)
2582 * RegularityIndicator::estimate_smoothness(
2583 * ErrorIndicator::get_eigenfunctions(), 0, smoothness_indicators);
2589 * this->output_solution();
2590 *
const double threshold_smoothness = smoothness_threshold_fraction;
2591 *
unsigned int num_refined = 0, num_coarsened = 0;
2592 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2593 * ErrorIndicator::PrimalSolver::max_degree)
2595 *
for (
const auto &cell :
2596 * ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2598 * if (cell->refine_flag_set())
2600 *
if (cell->coarsen_flag_set())
2602 *
if (cell->refine_flag_set() &&
2603 * smoothness_indicators[cell->active_cell_index()] >
2604 * threshold_smoothness &&
2605 *
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2606 * ErrorIndicator::get_FECollection()->size())
2608 * cell->clear_refine_flag();
2609 * cell->set_active_fe_index(cell->active_fe_index() + 1);
2611 *
else if (cell->coarsen_flag_set() &&
2612 * smoothness_indicators[cell->active_cell_index()] <
2613 * threshold_smoothness &&
2614 * cell->active_fe_index() != 0)
2616 * cell->clear_coarsen_flag();
2618 * cell->set_active_fe_index(cell->active_fe_index() - 1);
2622 * Here we also impose a limit on how small the cells can become
2625 *
else if (cell->refine_flag_set() && cell->diameter() < 5.0e-6)
2627 * cell->clear_refine_flag();
2628 *
if (
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2629 * ErrorIndicator::get_FECollection()->size())
2630 * cell->set_active_fe_index(cell->active_fe_index() + 1);
2637 * Check what the smallest
diameter is
2640 *
double min_diameter = std::numeric_limits<double>::max();
2641 *
for (
const auto &cell :
2642 * ErrorIndicator::get_DoFHandler()->active_cell_iterators())
2643 * if (cell->
diameter() < min_diameter)
2646 * std::cout <<
"Min diameter: " << min_diameter << std::endl;
2648 * ErrorIndicator::synchronize_discretization();
2650 * (this->
triangulation)->execute_coarsening_and_refinement();
2656 * main(
int argc,
char **argv)
2660 *
using namespace dealii;
2661 *
using namespace Maxwell;
2662 *
using namespace Refinement;
2663 *
using namespace ErrorIndicators;
2664 *
using namespace RegularityIndicators;
2672 * ExcMessage(
"This program can only be run in serial, use ./maxwell-hp"));
2675 * Structures::create_L_waveguide(triangulation_DWR, 2.0);
2676 * Structures::create_L_waveguide(triangulation_Kelly, 2.0);
2678 * Refiner<2, KellyErrorIndicator<2>, LegendreIndicator<2>> problem_Kelly(
2680 * triangulation_Kelly,
2685 * Refiner<2, DualWeightedResidual<2, false>, LegendreIndicator<2>>
2686 * problem_DWR(
"maxwell-hp.prm",
2687 * triangulation_DWR,
2694 * The threshold
for the
hp-decision: too small -> not enough
2695 * @f$h@f$-refinement, too large -> not enough @f$p@f$-refinement
2698 *
double smoothness_threshold = 0.75;
2700 * std::cout <<
"Executing refinement for the Kelly strategy!" << std::endl;
2701 * problem_Kelly.execute_refinement(smoothness_threshold);
2702 * std::cout <<
"...Done with Kelly refinement strategy!" << std::endl;
2703 * std::cout <<
"Executing refinement for the DWR strategy!" << std::endl;
2704 * problem_DWR.execute_refinement(smoothness_threshold);
2705 * std::cout <<
"...Done with DWR refinement strategy!" << std::endl;
2708 *
catch (std::exception &exc)
2710 * std::cerr << std::endl
2712 * <<
"----------------------------------------------------"
2714 * std::cerr <<
"Exception on processing: " << std::endl
2715 * << exc.what() << std::endl
2716 * <<
"Aborting!" << std::endl
2717 * <<
"----------------------------------------------------"
2724 * std::cerr << std::endl
2726 * <<
"----------------------------------------------------"
2728 * std::cerr <<
"Unknown exception!" << std::endl
2729 * <<
"Aborting!" << std::endl
2730 * <<
"----------------------------------------------------"
2735 * std::cout << std::endl <<
" Job done." << std::endl;
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 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.
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)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > 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)
VectorType::value_type * end(VectorType &V)
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)
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)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)