312 *
const unsigned int & i,
313 *
const unsigned int &
j,
314 *
const unsigned int &
q_point)
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;
338 *
output += fe_values.shape_value_component(i,
q_point,
comp) *
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
443 *
std::unique_ptr<ParameterHandler> parameters;
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);
481 *
parameters->get_integer(
"Eigenpair selection scheme");
490 *
"Selection by a target is the only currently supported option!");
492 *
parameters->get_integer(
"Number of eigenvalues/eigenfunctions");
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");
505 *
Base<dim>::set_refinement_cycle(
const unsigned int cycle)
521 *
const unsigned int &minimum_degree,
525 *
virtual unsigned int
528 *
virtual unsigned int
531 *
template <
class SolverType>
542 *
const std::unique_ptr<hp::FECollection<dim>> fe_collection;
544 *
std::unique_ptr<
hp::QCollection<dim - 1>> face_quadrature_collection;
552 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>>
eigenfunctions;
553 *
std::unique_ptr<std::vector<double>>
eigenvalues;
577 *
const unsigned int &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)
603 *
const QGauss<dim - 1> face_quadrature(degree + 1);
616 *
cell1 = dof_handler.begin_active(),
617 *
endc1 = dof_handler.end();
633 *
return &(*eigenvalues)[0];
655 *
for (
unsigned int i = 0; i < solution.size(); ++i)
666 *
template <
class SolverType>
675 *
switch (this->eigenpair_selection_scheme)
699 *
double shift_scalar = this->parameters->get_double(
"Target eigenvalue");
708 *
this->mpi_communicator, additional_data);
730 *
this->mpi_communicator);
746 *
constraints.distribute(entry);
750 *
return solver_control.last_step();
757 *
return dof_handler.n_dofs();
769 *
dof_handler.distribute_dofs(*fe_collection);
770 *
constraints.
clear();
773 *
constraints.close();
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;
823 *
hp_fe_values.reinit(cell);
825 *
const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
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)
853 *
local_dof_indices.resize(dofs_per_cell);
854 *
cell->get_dof_indices(local_dof_indices);
867 *
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
868 *
if (constraints.is_constrained(i))
893 *
const unsigned int &max_degree,
899 *
virtual unsigned int
900 *
n_dofs()
const override;
907 *
const unsigned int &max_degree,
926 *
data_out.attach_dof_handler(this->dof_handler);
928 *
for (
const auto &cell :
this->dof_handler.active_cell_iterators())
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;
942 * << this->dof_handler.n_dofs() << std::endl;
947 *
data_out.build_patches();
948 *
std::ofstream output(
"eigenvectors-" +
950 *
data_out.write_vtu(output);
975 *
const unsigned int &max_degree,
983 *
const unsigned int &max_degree,
1000 *
using namespace Maxwell;
1007 *
template <
int dim,
bool report_dual>
1027 *
virtual unsigned int
1033 *
virtual unsigned int
1034 *
n_dofs()
const override;
1054 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1057 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1060 *
std::unique_ptr<std::vector<double>> &
1063 *
std::unique_ptr<std::vector<double>> &
1105 *
typename std::map<typename DoFHandler<dim>::face_iterator,
double>;
1108 *
solve_primal_problem();
1111 *
solve_dual_problem();
1144 *
const unsigned int &
face_no,
1152 *
const unsigned int &
face_no,
1162 *
template <
int dim,
bool report_dual>
1188 *
template <
int dim,
bool report_dual>
1203 *
template <
int dim,
bool report_dual>
1215 *
template <
int dim,
bool report_dual>
1230 *
template <
int dim,
bool report_dual>
1237 *
template <
int dim,
bool report_dual>
1245 *
template <
int dim,
bool report_dual>
1246 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1256 *
template <
int dim,
bool report_dual>
1257 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
1264 *
template <
int dim,
bool report_dual>
1265 *
std::unique_ptr<std::vector<double>> &
1272 *
template <
int dim,
bool report_dual>
1273 *
std::unique_ptr<std::vector<double>> &
1279 *
template <
int dim,
bool report_dual>
1291 *
template <
int dim,
bool report_dual>
1303 *
template <
int dim,
bool report_dual>
1314 *
template <
int dim,
bool report_dual>
1325 *
template <
int dim,
bool report_dual>
1338 *
template <
int dim,
bool report_dual>
1366 *
cell2->set_active_fe_index(
cell1->active_fe_index());
1374 *
template <
int dim,
bool report_dual>
1380 * initialize
the cell fe_values...
1409 *
template <
int dim,
bool report_dual>
1416 *
for (
const auto &cell :
1417 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1436 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1453 *
template <
int dim,
bool report_dual>
1532 *
for (
const auto &cell :
1533 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1538 *
for (
const auto &cell :
1539 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1548 *
unsigned int present_cell = 0;
1549 *
for (
const auto &cell :
1550 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1555 *
ExcInternalError());
1570 *
std::cout <<
"Estimated QoI error: " << std::setprecision(20)
1578 *
template <
int dim,
bool report_dual>
1592 *
if (cell->face(
face_no)->at_boundary())
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)
1616 *
template <
int dim,
bool report_dual>
1632 *
std::vector<std::vector<Tensor<2, dim, double>>>
cell_hessians(
1644 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1664 *
template <
int dim,
bool report_dual>
1668 *
const unsigned int &
face_no,
1674 *
ExcInternalError());
1676 *
const auto neighbor = cell->neighbor(
face_no);
1679 *
std::max(cell->active_fe_index(), neighbor->active_fe_index());
1703 *
for (
unsigned int p = 0; p < n_q_points; ++p)
1715 *
ExcInternalError());
1723 *
template <
int dim,
bool report_dual>
1727 *
const unsigned int &
face_no,
1737 *
Assert(neighbor->has_children(), ExcInternalError());
1747 *
ExcInternalError());
1770 * initialize fe_face_values_neighbor
1780 *
const unsigned int n_q_points =
1789 *
for (
unsigned int p = 0; p < n_q_points; ++p)
1808 *
ExcInternalError());
1810 *
ExcInternalError());
1816 *
template <
int dim,
bool report_dual>
1827 *
for (
const auto &cell :
1828 *
DualSolver<dim>::dof_handler.active_cell_iterators())
1839 *
std::vector<Vector<double>>
cell_values(fe_values.n_quadrature_points,
1843 *
for (
unsigned int p = 0; p < fe_values.n_quadrature_points; ++p)
1860 *
template <
int dim,
bool report_dual>
1869 *
assert(
u2.size() ==
dof2.n_dofs() &&
"Incorrect input vector size!");
1884 *
assert(
fe1.degree <
fe2.degree &&
"Incorrect usage of embed!");
1898 *
fe2.get_embedding_dofs(
fe1.degree);
1925 *
constraints.distribute(
u2);
1928 *
template <
int dim,
bool report_dual>
1942 *
assert(
u2.size() ==
dof2.n_dofs() &&
"Incorrect input vector size!");
1957 *
assert(
fe1.degree >
fe2.degree &&
"Incorrect usage of extract!");
1965 *
fe1.get_embedding_dofs(
fe2.degree);
1992 *
constraints.distribute(
u2);
1994 *
template <
int dim,
bool report_dual>
1997 *
std::ofstream &os)
2004 *
template <
int dim,
bool report_dual>
2007 *
std::ofstream &os)
2016 *
template <
int dim>
2032 *
const unsigned int &max_degree,
2035 *
virtual unsigned int
2047 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2050 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2053 *
std::unique_ptr<std::vector<double>> &
2088 *
std::vector<std::shared_ptr<Vector<float>>>
errors;
2091 *
template <
int dim>
2096 *
const unsigned int &max_degree,
2106 *
template <
int dim>
2113 *
template <
int dim>
2120 *
template <
int dim>
2127 *
template <
int dim>
2128 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2134 *
template <
int dim>
2135 *
std::unique_ptr<std::vector<double>> &
2141 *
template <
int dim>
2142 *
std::unique_ptr<std::vector<PETScWrappers::MPI::Vector>> &
2148 *
template <
int dim>
2155 *
template <
int dim>
2162 *
template <
int dim>
2174 *
template <
int dim>
2181 *
template <
int dim>
2185 *
unsigned int count = 0;
2190 *
if (
count >= this->n_eigenpairs)
2200 *
template <
int dim>
2204 *
std::cout <<
"Marking cells via Kelly indicator..." << std::endl;
2228 *
*this->face_quadrature_collection,
2234 *
*this->face_quadrature_collection,
2248 *
std::cout <<
"...Done!" << std::endl;
2250 *
template <
int dim>
2257 *
template <
int dim>
2271 *
using namespace dealii;
2275 *
template <
int dim>
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>
2337 *
template <
int dim>
2346 *
template <
class VectorType>
2349 *
const std::unique_ptr<std::vector<VectorType>> &
eigenfunctions,
2357 *
template <
int dim>
2368 *
template <
int dim>
2369 *
template <
class VectorType>
2372 *
const std::unique_ptr<std::vector<VectorType>> &
eigenfunctions,
2385 *
namespace Refinement
2387 *
using namespace dealii;
2388 *
using namespace Maxwell;
2390 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2397 *
const unsigned int &max_degree,
2413 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2418 *
const unsigned int &max_degree,
2428 *
if (ErrorIndicator::name() ==
"DWR")
2434 *
eigenvalues_out.open(
"eigenvalues_" + ErrorIndicator::name() +
"_out.txt");
2443 *
template <
int dim>
2452 *
evaluate_vector_field(
2458 *
for (
unsigned int p = 0; p <
input_data.solution_gradients.size(); ++p)
2473 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2483 *
for (
const auto &cell :
output_dof.active_cell_iterators())
2487 *
data_out.add_data_vector(
fe_degrees,
"fe_degree");
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())
2501 *
std::string(
"eigenfunction_no_") +
2513 *
data_out.build_patches();
2514 *
std::ofstream output(
"eigenvectors-" + ErrorIndicator::name() +
"-" +
2516 *
data_out.write_vtu(output);
2525 *
template <
int dim,
class ErrorIndicator,
class RegularityIndicator>
2536 *
std::cout <<
"Initializing RegularityIndicator..." << std::endl;
2538 *
<<
"(This may take a while if the max expansion order is set too high)"
2542 *
std::cout <<
"Done!" << std::endl <<
"Starting Refinement..." << std::endl;
2544 *
for (
unsigned int cycle = 0; cycle <= this->max_cycles; ++cycle)
2547 *
std::cout <<
"Cycle: " << cycle << std::endl;
2549 *
this->estimated_error_per_cell.reinit(
2565 *
*this->
triangulation, estimated_error_per_cell, 1. / 5., 0.000);
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)
2592 *
if (ErrorIndicator::PrimalSolver::min_degree !=
2593 *
ErrorIndicator::PrimalSolver::max_degree)
2595 *
for (
const auto &cell :
2598 *
if (cell->refine_flag_set())
2600 *
if (cell->coarsen_flag_set())
2602 *
if (cell->refine_flag_set() &&
2605 *
static_cast<unsigned int>(cell->active_fe_index() + 1) <
2608 *
cell->clear_refine_flag();
2609 *
cell->set_active_fe_index(cell->active_fe_index() + 1);
2611 *
else if (cell->coarsen_flag_set() &&
2614 *
cell->active_fe_index() != 0)
2616 *
cell->clear_coarsen_flag();
2618 *
cell->set_active_fe_index(cell->active_fe_index() - 1);
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) <
2630 *
cell->set_active_fe_index(cell->active_fe_index() + 1);
2640 *
double min_diameter = std::numeric_limits<double>::max();
2641 *
for (
const auto &cell :
2646 *
std::cout <<
"Min diameter: " <<
min_diameter << std::endl;
2650 *
(this->
triangulation)->execute_coarsening_and_refinement();
2660 *
using namespace dealii;
2661 *
using namespace Maxwell;
2662 *
using namespace Refinement;
2672 *
ExcMessage(
"This program can only be run in serial, use ./maxwell-hp"));
2700 *
std::cout <<
"Executing refinement for the Kelly strategy!" << std::endl;
2702 *
std::cout <<
"...Done with Kelly refinement strategy!" << std::endl;
2703 *
std::cout <<
"Executing refinement for the DWR strategy!" << std::endl;
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)
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)
Number l1_norm(const Tensor< 2, dim, Number > &t)
#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)