284 const std::vector<const InputVector *> &solutions,
294 using number =
typename InputVector::value_type;
303 "For distributed Triangulation objects and associated "
304 "DoFHandler objects, asking for any subdomain other than the "
305 "locally owned one does not make sense."));
310 subdomain_id = subdomain_id_;
314 const unsigned int n_solution_vectors = solutions.size();
319 ExcMessage(
"You are not allowed to list the special boundary "
320 "indicator for internal boundaries in your boundary "
323 for (
const auto &boundary_function : neumann_bc)
325 (void)boundary_function;
326 Assert(boundary_function.second->n_components == n_components,
327 ExcInvalidBoundaryFunction(boundary_function.first,
328 boundary_function.second->n_components,
333 ExcInvalidComponentMask());
335 ExcInvalidComponentMask());
337 Assert((coefficient ==
nullptr) ||
340 ExcInvalidCoefficient());
342 Assert(solutions.size() > 0, ExcNoSolutions());
343 Assert(solutions.size() == errors.size(),
344 ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
345 for (
unsigned int n = 0; n < solutions.size(); ++n)
349 Assert((coefficient ==
nullptr) ||
352 ExcInvalidCoefficient());
354 for (
const auto &boundary_function : neumann_bc)
356 (void)boundary_function;
357 Assert(boundary_function.second->n_components == n_components,
358 ExcInvalidBoundaryFunction(boundary_function.first,
359 boundary_function.second->n_components,
364 for (
unsigned int n = 0; n < n_solution_vectors; ++n)
371 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
372 gradients_here(n_solution_vectors,
376 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
377 gradients_neighbor(gradients_here);
378 std::vector<Vector<typename ProductType<number, double>::type>>
379 grad_dot_n_neighbor(n_solution_vectors,
387 if (coefficient ==
nullptr)
388 for (
unsigned int c = 0; c < n_components; ++c)
389 coefficient_values(c) = 1;
413 (cell->subdomain_id() == subdomain_id)) &&
415 (cell->material_id() == material_id)))
417 for (
unsigned int n = 0; n < n_solution_vectors; ++n)
418 (*errors[n])(cell->active_cell_index()) = 0;
421 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
423 *solutions[s], gradients_here[s]);
427 for (
unsigned int n = 0; n < 2; ++n)
430 auto neighbor = cell->neighbor(n);
432 while (neighbor->has_children())
433 neighbor = neighbor->child(n == 0 ? 1 : 0);
435 fe_face_values.
reinit(cell, n);
441 fe_values.
reinit(neighbor);
443 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
445 *solutions[s], gradients_neighbor[s]);
447 fe_face_values.
reinit(neighbor, n == 0 ? 1 : 0);
450 .get_normal_vectors()[0];
454 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
455 for (
unsigned int c = 0; c < n_components; ++c)
456 grad_dot_n_neighbor[s](c) =
457 -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
460 else if (neumann_bc.find(n) != neumann_bc.end())
464 if (n_components == 1)
466 const typename InputVector::value_type v =
467 neumann_bc.find(n)->second->value(cell->vertex(n));
469 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
470 grad_dot_n_neighbor[s](0) = v;
475 neumann_bc.find(n)->second->vector_value(cell->vertex(n),
478 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
479 grad_dot_n_neighbor[s] = v;
484 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
485 grad_dot_n_neighbor[s] = 0;
489 if (coefficient !=
nullptr)
493 const double c_value = coefficient->
value(cell->vertex(n));
494 for (
unsigned int c = 0; c < n_components; ++c)
495 coefficient_values(c) = c_value;
503 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
504 for (
unsigned int component = 0; component < n_components;
506 if (component_mask[component] ==
true)
511 gradients_here[s][n][component] * normal;
514 ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
515 coefficient_values(component));
516 (*errors[s])(cell->active_cell_index()) +=
519 double>::type>::abs_square(jump) *
524 for (
unsigned int s = 0; s < n_solution_vectors; ++s)
525 (*errors[s])(cell->active_cell_index()) =
526 std::sqrt((*errors[s])(cell->active_cell_index()));
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, typename InputVector::value_type > * > &neumann_bc, const InputVector &solution, Vector< float > &error, const ComponentMask &component_mask=ComponentMask(), 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)