Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
error_estimator_1d.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
20 
22 
25 
26 #include <deal.II/fe/fe.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping_q1.h>
30 
32 
33 #include <deal.II/hp/fe_values.h>
36 
40 #include <deal.II/lac/la_vector.h>
45 #include <deal.II/lac/vector.h>
46 
48 
49 #include <algorithm>
50 #include <cmath>
51 #include <functional>
52 #include <numeric>
53 #include <vector>
54 
56 
57 
58 
59 template <int spacedim>
60 template <typename InputVector>
61 void
63  const Mapping<1, spacedim> & mapping,
64  const DoFHandler<1, spacedim> &dof_handler,
65  const Quadrature<0> & quadrature,
66  const std::map<types::boundary_id,
68  & neumann_bc,
69  const InputVector & solution,
70  Vector<float> & error,
71  const ComponentMask & component_mask,
72  const Function<spacedim> *coefficients,
73  const unsigned int n_threads,
76  const Strategy strategy)
77 {
78  // just pass on to the other function
79  const std::vector<const InputVector *> solutions(1, &solution);
80  std::vector<Vector<float> *> errors(1, &error);
81  estimate(mapping,
82  dof_handler,
83  quadrature,
84  neumann_bc,
85  solutions,
86  errors,
87  component_mask,
88  coefficients,
89  n_threads,
90  subdomain_id,
91  material_id,
92  strategy);
93 }
94 
95 
96 
97 template <int spacedim>
98 template <typename InputVector>
99 void
101  const DoFHandler<1, spacedim> &dof_handler,
102  const Quadrature<0> & quadrature,
103  const std::map<types::boundary_id,
105  & neumann_bc,
106  const InputVector & solution,
107  Vector<float> & error,
108  const ComponentMask & component_mask,
109  const Function<spacedim> *coefficients,
110  const unsigned int n_threads,
113  const Strategy strategy)
114 {
116  dof_handler,
117  quadrature,
118  neumann_bc,
119  solution,
120  error,
121  component_mask,
122  coefficients,
123  n_threads,
124  subdomain_id,
125  material_id,
126  strategy);
127 }
128 
129 
130 
131 template <int spacedim>
132 template <typename InputVector>
133 void
135  const DoFHandler<1, spacedim> &dof_handler,
136  const Quadrature<0> & quadrature,
137  const std::map<types::boundary_id,
139  & neumann_bc,
140  const std::vector<const InputVector *> &solutions,
141  std::vector<Vector<float> *> & errors,
142  const ComponentMask & component_mask,
143  const Function<spacedim> * coefficients,
144  const unsigned int n_threads,
147  const Strategy strategy)
148 {
150  dof_handler,
151  quadrature,
152  neumann_bc,
153  solutions,
154  errors,
155  component_mask,
156  coefficients,
157  n_threads,
158  subdomain_id,
159  material_id,
160  strategy);
161 }
162 
163 
164 
165 template <int spacedim>
166 template <typename InputVector>
167 void
169  const Mapping<1, spacedim> & mapping,
170  const DoFHandler<1, spacedim> &dof_handler,
171  const hp::QCollection<0> & quadrature,
172  const std::map<types::boundary_id,
174  & neumann_bc,
175  const InputVector & solution,
176  Vector<float> & error,
177  const ComponentMask & component_mask,
178  const Function<spacedim> *coefficients,
179  const unsigned int n_threads,
182  const Strategy strategy)
183 {
184  // just pass on to the other function
185  const std::vector<const InputVector *> solutions(1, &solution);
186  std::vector<Vector<float> *> errors(1, &error);
187  estimate(mapping,
188  dof_handler,
189  quadrature,
190  neumann_bc,
191  solutions,
192  errors,
193  component_mask,
194  coefficients,
195  n_threads,
196  subdomain_id,
197  material_id,
198  strategy);
199 }
200 
201 
202 template <int spacedim>
203 template <typename InputVector>
204 void
206  const DoFHandler<1, spacedim> &dof_handler,
207  const hp::QCollection<0> & quadrature,
208  const std::map<types::boundary_id,
210  & neumann_bc,
211  const InputVector & solution,
212  Vector<float> & error,
213  const ComponentMask & component_mask,
214  const Function<spacedim> *coefficients,
215  const unsigned int n_threads,
218  const Strategy strategy)
219 {
221  dof_handler,
222  quadrature,
223  neumann_bc,
224  solution,
225  error,
226  component_mask,
227  coefficients,
228  n_threads,
229  subdomain_id,
230  material_id,
231  strategy);
232 }
233 
234 
235 
236 template <int spacedim>
237 template <typename InputVector>
238 void
240  const DoFHandler<1, spacedim> &dof_handler,
241  const hp::QCollection<0> & quadrature,
242  const std::map<types::boundary_id,
244  & neumann_bc,
245  const std::vector<const InputVector *> &solutions,
246  std::vector<Vector<float> *> & errors,
247  const ComponentMask & component_mask,
248  const Function<spacedim> * coefficients,
249  const unsigned int n_threads,
252  const Strategy strategy)
253 {
255  dof_handler,
256  quadrature,
257  neumann_bc,
258  solutions,
259  errors,
260  component_mask,
261  coefficients,
262  n_threads,
263  subdomain_id,
264  material_id,
265  strategy);
266 }
267 
268 
269 
270 template <int spacedim>
271 template <typename InputVector>
272 void
274  const Mapping<1, spacedim> & /*mapping*/,
275  const DoFHandler<1, spacedim> & /*dof_handler*/,
276  const hp::QCollection<0> &,
277  const std::map<types::boundary_id,
279  & /*neumann_bc*/,
280  const std::vector<const InputVector *> & /*solutions*/,
281  std::vector<Vector<float> *> & /*errors*/,
282  const ComponentMask & /*component_mask_*/,
283  const Function<spacedim> * /*coefficient*/,
284  const unsigned int,
285  const types::subdomain_id /*subdomain_id*/,
286  const types::material_id /*material_id*/,
287  const Strategy /*strategy*/)
288 {
289  Assert(false, ExcInternalError());
290 }
291 
292 
293 
294 template <int spacedim>
295 template <typename InputVector>
296 void
298  const Mapping<1, spacedim> & mapping,
299  const DoFHandler<1, spacedim> &dof_handler,
300  const Quadrature<0> &,
301  const std::map<types::boundary_id,
303  & neumann_bc,
304  const std::vector<const InputVector *> &solutions,
305  std::vector<Vector<float> *> & errors,
306  const ComponentMask & component_mask,
307  const Function<spacedim> * coefficient,
308  const unsigned int,
309  const types::subdomain_id subdomain_id_,
311  const Strategy strategy)
312 {
314  using number = typename InputVector::value_type;
315 #ifdef DEAL_II_WITH_P4EST
316  if (dynamic_cast<const parallel::distributed::Triangulation<1, spacedim> *>(
317  &dof_handler.get_triangulation()) != nullptr)
318  Assert((subdomain_id_ == numbers::invalid_subdomain_id) ||
319  (subdomain_id_ ==
320  dynamic_cast<
322  dof_handler.get_triangulation())
323  .locally_owned_subdomain()),
324  ExcMessage(
325  "For parallel distributed triangulations, the only "
326  "valid subdomain_id that can be passed here is the "
327  "one that corresponds to the locally owned subdomain id."));
328 
331  &dof_handler.get_triangulation()) != nullptr) ?
333  dof_handler.get_triangulation())
334  .locally_owned_subdomain() :
335  subdomain_id_);
336 #else
337  const types::subdomain_id subdomain_id = subdomain_id_;
338 #endif
339 
340  const unsigned int n_components = dof_handler.get_fe(0).n_components();
341  const unsigned int n_solution_vectors = solutions.size();
342 
343  // sanity checks
344  Assert(neumann_bc.find(numbers::internal_face_boundary_id) ==
345  neumann_bc.end(),
346  ExcMessage("You are not allowed to list the special boundary "
347  "indicator for internal boundaries in your boundary "
348  "value map."));
349 
350  for (const auto &boundary_function : neumann_bc)
351  {
352  (void)boundary_function;
353  Assert(boundary_function.second->n_components == n_components,
354  ExcInvalidBoundaryFunction(boundary_function.first,
355  boundary_function.second->n_components,
356  n_components));
357  }
358 
359  Assert(component_mask.represents_n_components(n_components),
361  Assert(component_mask.n_selected_components(n_components) > 0,
363 
364  Assert((coefficient == nullptr) ||
365  (coefficient->n_components == n_components) ||
366  (coefficient->n_components == 1),
368 
369  Assert(solutions.size() > 0, ExcNoSolutions());
370  Assert(solutions.size() == errors.size(),
371  ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
372  for (unsigned int n = 0; n < solutions.size(); ++n)
373  Assert(solutions[n]->size() == dof_handler.n_dofs(),
374  ExcDimensionMismatch(solutions[n]->size(), dof_handler.n_dofs()));
375 
376  Assert((coefficient == nullptr) ||
377  (coefficient->n_components == n_components) ||
378  (coefficient->n_components == 1),
380 
381  for (const auto &boundary_function : neumann_bc)
382  {
383  (void)boundary_function;
384  Assert(boundary_function.second->n_components == n_components,
385  ExcInvalidBoundaryFunction(boundary_function.first,
386  boundary_function.second->n_components,
387  n_components));
388  }
389 
390  // reserve one slot for each cell and set it to zero
391  for (unsigned int n = 0; n < n_solution_vectors; ++n)
392  (*errors[n]).reinit(dof_handler.get_triangulation().n_active_cells());
393 
394  // fields to get the gradients on the present and the neighbor cell.
395  //
396  // for the neighbor gradient, we need several auxiliary fields, depending on
397  // the way we get it (see below)
398  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
399  gradients_here(n_solution_vectors,
400  std::vector<std::vector<Tensor<1, spacedim, number>>>(
401  2,
402  std::vector<Tensor<1, spacedim, number>>(n_components)));
403  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
404  gradients_neighbor(gradients_here);
405  std::vector<Vector<typename ProductType<number, double>::type>>
406  grad_dot_n_neighbor(n_solution_vectors,
408  n_components));
409 
410  // reserve some space for coefficient values at one point. if there is no
411  // coefficient, then we fill it by unity once and for all and don't set it
412  // any more
413  Vector<double> coefficient_values(n_components);
414  if (coefficient == nullptr)
415  for (unsigned int c = 0; c < n_components; ++c)
416  coefficient_values(c) = 1;
417 
418  const QTrapezoid<1> quadrature;
419  const hp::QCollection<1> q_collection(quadrature);
420  const QGauss<0> face_quadrature(1);
421  const hp::QCollection<0> q_face_collection(face_quadrature);
422 
423  const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection();
424 
425  hp::MappingCollection<1, spacedim> mapping_collection;
426  mapping_collection.push_back(mapping);
427 
428  hp::FEValues<1, spacedim> fe_values(mapping_collection,
429  fe,
430  q_collection,
432  hp::FEFaceValues<1, spacedim> fe_face_values(
433  /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors);
434 
435  // loop over all cells and do something on the cells which we're told to
436  // work on. note that the error indicator is only a sum over the two
437  // contributions from the two vertices of each cell.
438  for (const auto &cell : dof_handler.active_cell_iterators())
439  if (((subdomain_id == numbers::invalid_subdomain_id) ||
440  (cell->subdomain_id() == subdomain_id)) &&
441  ((material_id == numbers::invalid_material_id) ||
442  (cell->material_id() == material_id)))
443  {
444  for (unsigned int n = 0; n < n_solution_vectors; ++n)
445  (*errors[n])(cell->active_cell_index()) = 0;
446 
447  fe_values.reinit(cell);
448  for (unsigned int s = 0; s < n_solution_vectors; ++s)
449  fe_values.get_present_fe_values().get_function_gradients(
450  *solutions[s], gradients_here[s]);
451 
452  // loop over the two points bounding this line. n==0 is left point,
453  // n==1 is right point
454  for (unsigned int n = 0; n < 2; ++n)
455  {
456  // find left or right active neighbor
457  auto neighbor = cell->neighbor(n);
458  if (neighbor.state() == IteratorState::valid)
459  while (neighbor->has_children())
460  neighbor = neighbor->child(n == 0 ? 1 : 0);
461 
462  fe_face_values.reinit(cell, n);
463  Tensor<1, spacedim> normal =
464  fe_face_values.get_present_fe_values().get_normal_vectors()[0];
465 
466  if (neighbor.state() == IteratorState::valid)
467  {
468  fe_values.reinit(neighbor);
469 
470  for (unsigned int s = 0; s < n_solution_vectors; ++s)
471  fe_values.get_present_fe_values().get_function_gradients(
472  *solutions[s], gradients_neighbor[s]);
473 
474  fe_face_values.reinit(neighbor, n == 0 ? 1 : 0);
475  Tensor<1, spacedim> neighbor_normal =
476  fe_face_values.get_present_fe_values()
477  .get_normal_vectors()[0];
478 
479  // extract the gradient in normal direction of all the
480  // components.
481  for (unsigned int s = 0; s < n_solution_vectors; ++s)
482  for (unsigned int c = 0; c < n_components; ++c)
483  grad_dot_n_neighbor[s](c) =
484  -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
485  neighbor_normal);
486  }
487  else if (neumann_bc.find(n) != neumann_bc.end())
488  // if Neumann b.c., then fill the gradients field which will be
489  // used later on.
490  {
491  if (n_components == 1)
492  {
493  const typename InputVector::value_type v =
494  neumann_bc.find(n)->second->value(cell->vertex(n));
495 
496  for (unsigned int s = 0; s < n_solution_vectors; ++s)
497  grad_dot_n_neighbor[s](0) = v;
498  }
499  else
500  {
501  Vector<typename InputVector::value_type> v(n_components);
502  neumann_bc.find(n)->second->vector_value(cell->vertex(n),
503  v);
504 
505  for (unsigned int s = 0; s < n_solution_vectors; ++s)
506  grad_dot_n_neighbor[s] = v;
507  }
508  }
509  else
510  // fill with zeroes.
511  for (unsigned int s = 0; s < n_solution_vectors; ++s)
512  grad_dot_n_neighbor[s] = 0;
513 
514  // if there is a coefficient, then evaluate it at the present
515  // position. if there is none, reuse the preset values.
516  if (coefficient != nullptr)
517  {
518  if (coefficient->n_components == 1)
519  {
520  const double c_value = coefficient->value(cell->vertex(n));
521  for (unsigned int c = 0; c < n_components; ++c)
522  coefficient_values(c) = c_value;
523  }
524  else
525  coefficient->vector_value(cell->vertex(n),
526  coefficient_values);
527  }
528 
529 
530  for (unsigned int s = 0; s < n_solution_vectors; ++s)
531  for (unsigned int component = 0; component < n_components;
532  ++component)
533  if (component_mask[component] == true)
534  {
535  // get gradient here
536  const typename ProductType<number, double>::type
537  grad_dot_n_here =
538  gradients_here[s][n][component] * normal;
539 
540  const typename ProductType<number, double>::type jump =
541  ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
542  coefficient_values(component));
543  (*errors[s])(cell->active_cell_index()) +=
545  typename ProductType<number,
546  double>::type>::abs_square(jump) *
547  cell->diameter();
548  }
549  }
550 
551  for (unsigned int s = 0; s < n_solution_vectors; ++s)
552  (*errors[s])(cell->active_cell_index()) =
553  std::sqrt((*errors[s])(cell->active_cell_index()));
554  }
555 }
556 
557 
558 // explicit instantiations
559 #include "error_estimator_1d.inst"
560 
561 
static ::ExceptionBase & ExcInvalidBoundaryFunction(types::boundary_id arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidCoefficient()
const types::subdomain_id invalid_subdomain_id
Definition: types.h:276
const unsigned int n_components
Definition: function.h:164
static ::ExceptionBase & ExcInvalidComponentMask()
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type
void push_back(const Mapping< dim, spacedim > &new_mapping)
bool represents_n_components(const unsigned int n) const
unsigned int material_id
Definition: types.h:152
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)
Definition: fe_values.cc:423
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
#define Assert(cond, exc)
Definition: exceptions.h:1465
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index n_dofs() const
Abstract base class for mapping classes.
Definition: mapping.h:303
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
static ::ExceptionBase & ExcIncompatibleNumberOfElements(int arg1, int arg2)
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &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)
Definition: fe_values.cc:293
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)
Definition: tensor.h:462
const Triangulation< dim, spacedim > & get_triangulation() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
Shape function gradients.
Normal vectors.
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
static ::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
const types::boundary_id internal_face_boundary_id
Definition: types.h:255
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
Kelly error estimator with the factor .
const types::material_id invalid_material_id
Definition: types.h:228
static ::ExceptionBase & ExcNoSolutions()
static ::ExceptionBase & ExcInternalError()
const ::FEValues< dim, spacedim > & get_present_fe_values() const
Definition: fe_values.h:693