57 template <
int spacedim>
58 template <
typename Number>
70 const unsigned int n_threads,
76 std::vector<const ReadVector<Number> *> solutions(1, &solution);
77 std::vector<Vector<float> *> errors(1, &error);
95 template <
int spacedim>
96 template <
typename Number>
107 const unsigned int n_threads,
129 template <
int spacedim>
130 template <
typename Number>
141 const unsigned int n_threads,
163 template <
int spacedim>
164 template <
typename Number>
176 const unsigned int n_threads,
182 std::vector<const ReadVector<Number> *> solutions(1, &solution);
183 std::vector<Vector<float> *> errors(1, &error);
201 template <
int spacedim>
202 template <
typename Number>
213 const unsigned int n_threads,
235 template <
int spacedim>
236 template <
typename Number>
247 const unsigned int n_threads,
269 template <
int spacedim>
270 template <
typename Number>
288 using number = Number;
297 "For distributed Triangulation objects and associated "
298 "DoFHandler objects, asking for any subdomain other than the "
299 "locally owned one does not make sense."));
308 const unsigned int n_solution_vectors = solutions.size();
313 ExcMessage(
"You are not allowed to list the special boundary "
314 "indicator for internal boundaries in your boundary "
317 for (
const auto &boundary_function : neumann_bc)
319 (void)boundary_function;
320 Assert(boundary_function.second->n_components == n_components,
322 boundary_function.second->n_components,
331 Assert((coefficient ==
nullptr) ||
337 Assert(solutions.size() == errors.size(),
339 for (
unsigned int n = 0; n < solutions.size(); ++n)
343 Assert((coefficient ==
nullptr) ||
348 for (
const auto &boundary_function : neumann_bc)
350 (void)boundary_function;
351 Assert(boundary_function.second->n_components == n_components,
353 boundary_function.second->n_components,
358 for (
unsigned int n = 0; n < n_solution_vectors; ++n)
365 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
366 gradients_here(n_solution_vectors,
370 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
371 gradients_neighbor(gradients_here);
372 std::vector<Vector<typename ProductType<number, double>::type>>
373 grad_dot_n_neighbor(n_solution_vectors,
381 if (coefficient ==
nullptr)
382 for (
unsigned int c = 0; c < n_components; ++c)
383 coefficient_values(c) = 1;
411 for (
unsigned int n = 0; n < n_solution_vectors; ++n)
412 (*errors[n])(cell->active_cell_index()) = 0;
415 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
417 *solutions[s], gradients_here[s]);
421 for (
unsigned int n = 0; n < 2; ++n)
424 auto neighbor = cell->neighbor(n);
426 while (neighbor->has_children())
427 neighbor = neighbor->child(n == 0 ? 1 : 0);
429 fe_face_values.
reinit(cell, n);
435 fe_values.
reinit(neighbor);
437 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
439 *solutions[s], gradients_neighbor[s]);
441 fe_face_values.
reinit(neighbor, n == 0 ? 1 : 0);
444 .get_normal_vectors()[0];
448 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
449 for (
unsigned int c = 0; c < n_components; ++c)
450 grad_dot_n_neighbor[s](c) =
451 -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
454 else if (neumann_bc.find(n) != neumann_bc.end())
458 if (n_components == 1)
461 neumann_bc.find(n)->second->value(cell->vertex(n));
463 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
464 grad_dot_n_neighbor[s](0) = v;
469 neumann_bc.find(n)->second->vector_value(cell->vertex(n),
472 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
473 grad_dot_n_neighbor[s] = v;
478 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
479 grad_dot_n_neighbor[s] = 0;
483 if (coefficient !=
nullptr)
487 const double c_value = coefficient->
value(cell->vertex(n));
488 for (
unsigned int c = 0; c < n_components; ++c)
489 coefficient_values(c) = c_value;
497 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
498 for (
unsigned int component = 0; component < n_components;
500 if (component_mask[component] ==
true)
505 gradients_here[s][n][component] * normal;
508 ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
509 coefficient_values(component));
510 (*errors[s])(cell->active_cell_index()) +=
513 double>::type>::abs_square(jump) *
518 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
519 (*errors[s])(cell->active_cell_index()) =
520 std::sqrt((*errors[s])(cell->active_cell_index()));
526 template <
int spacedim>
527 template <
typename Number>
539 const unsigned int n_threads,
547 quadrature_collection,
562 #include "error_estimator_1d.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
bool represents_n_components(const unsigned int n) const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number >> &gradients) const
unsigned int n_components() const
const unsigned int n_components
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &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)
@ cell_diameter_over_24
Kelly error estimator with the factor .
Abstract base class for mapping classes.
unsigned int n_active_cells() const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda >> &cell, const unsigned int face_no, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
const FEValuesType & get_present_fe_values() const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, lda >> &cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
void push_back(const Mapping< dim, spacedim > &new_mapping)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcIncompatibleNumberOfElements(int arg1, int arg2)
static ::ExceptionBase & ExcInvalidCoefficient()
static ::ExceptionBase & ExcNoSolutions()
static ::ExceptionBase & ExcInvalidBoundaryFunction(types::boundary_id arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidComponentMask()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
@ valid
Iterator points to a valid object.
constexpr const ReferenceCell Line
const types::material_id invalid_material_id
const types::boundary_id internal_face_boundary_id
const types::subdomain_id invalid_subdomain_id
unsigned int subdomain_id
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type