Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20:02+00:00
\(\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 - 2021 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 template <int spacedim>
59 template <typename InputVector>
60 void
62  const Mapping<1, spacedim> & mapping,
63  const DoFHandler<1, spacedim> &dof_handler,
64  const Quadrature<0> & quadrature,
65  const std::map<types::boundary_id,
67  & neumann_bc,
68  const InputVector & solution,
69  Vector<float> & error,
70  const ComponentMask & component_mask,
71  const Function<spacedim> *coefficients,
72  const unsigned int n_threads,
75  const Strategy strategy)
76 {
77  // just pass on to the other function
78  const std::vector<const InputVector *> solutions(1, &solution);
79  std::vector<Vector<float> *> errors(1, &error);
80  estimate(mapping,
81  dof_handler,
82  quadrature,
83  neumann_bc,
84  solutions,
85  errors,
86  component_mask,
87  coefficients,
88  n_threads,
91  strategy);
92 }
93 
94 
95 
96 template <int spacedim>
97 template <typename InputVector>
98 void
100  const DoFHandler<1, spacedim> &dof_handler,
101  const Quadrature<0> & quadrature,
102  const std::map<types::boundary_id,
104  & neumann_bc,
105  const InputVector & solution,
106  Vector<float> & error,
107  const ComponentMask & component_mask,
108  const Function<spacedim> *coefficients,
109  const unsigned int n_threads,
112  const Strategy strategy)
113 {
115  dof_handler,
116  quadrature,
117  neumann_bc,
118  solution,
119  error,
120  component_mask,
121  coefficients,
122  n_threads,
123  subdomain_id,
124  material_id,
125  strategy);
126 }
127 
128 
129 
130 template <int spacedim>
131 template <typename InputVector>
132 void
134  const DoFHandler<1, spacedim> &dof_handler,
135  const Quadrature<0> & quadrature,
136  const std::map<types::boundary_id,
138  & neumann_bc,
139  const std::vector<const InputVector *> &solutions,
140  std::vector<Vector<float> *> & errors,
141  const ComponentMask & component_mask,
142  const Function<spacedim> * coefficients,
143  const unsigned int n_threads,
146  const Strategy strategy)
147 {
149  dof_handler,
150  quadrature,
151  neumann_bc,
152  solutions,
153  errors,
154  component_mask,
155  coefficients,
156  n_threads,
157  subdomain_id,
158  material_id,
159  strategy);
160 }
161 
162 
163 
164 template <int spacedim>
165 template <typename InputVector>
166 void
168  const Mapping<1, spacedim> & mapping,
169  const DoFHandler<1, spacedim> &dof_handler,
170  const hp::QCollection<0> & quadrature,
171  const std::map<types::boundary_id,
173  & neumann_bc,
174  const InputVector & solution,
175  Vector<float> & error,
176  const ComponentMask & component_mask,
177  const Function<spacedim> *coefficients,
178  const unsigned int n_threads,
181  const Strategy strategy)
182 {
183  // just pass on to the other function
184  const std::vector<const InputVector *> solutions(1, &solution);
185  std::vector<Vector<float> *> errors(1, &error);
186  estimate(mapping,
187  dof_handler,
188  quadrature,
189  neumann_bc,
190  solutions,
191  errors,
192  component_mask,
193  coefficients,
194  n_threads,
195  subdomain_id,
196  material_id,
197  strategy);
198 }
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_,
287  const Strategy strategy)
288 {
290  using number = typename InputVector::value_type;
292  if (const auto *triangulation = dynamic_cast<
294  &dof_handler.get_triangulation()))
295  {
296  Assert((subdomain_id_ == numbers::invalid_subdomain_id) ||
297  (subdomain_id_ == triangulation->locally_owned_subdomain()),
298  ExcMessage(
299  "For distributed Triangulation objects and associated "
300  "DoFHandler objects, asking for any subdomain other than the "
301  "locally owned one does not make sense."));
302  subdomain_id = triangulation->locally_owned_subdomain();
303  }
304  else
305  {
306  subdomain_id = subdomain_id_;
307  }
308 
309  const unsigned int n_components = dof_handler.get_fe(0).n_components();
310  const unsigned int n_solution_vectors = solutions.size();
311 
312  // sanity checks
313  Assert(neumann_bc.find(numbers::internal_face_boundary_id) ==
314  neumann_bc.end(),
315  ExcMessage("You are not allowed to list the special boundary "
316  "indicator for internal boundaries in your boundary "
317  "value map."));
318 
319  for (const auto &boundary_function : neumann_bc)
320  {
321  (void)boundary_function;
322  Assert(boundary_function.second->n_components == n_components,
323  ExcInvalidBoundaryFunction(boundary_function.first,
324  boundary_function.second->n_components,
325  n_components));
326  }
327 
328  Assert(component_mask.represents_n_components(n_components),
330  Assert(component_mask.n_selected_components(n_components) > 0,
332 
333  Assert((coefficient == nullptr) ||
334  (coefficient->n_components == n_components) ||
335  (coefficient->n_components == 1),
337 
338  Assert(solutions.size() > 0, ExcNoSolutions());
339  Assert(solutions.size() == errors.size(),
340  ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
341  for (unsigned int n = 0; n < solutions.size(); ++n)
342  Assert(solutions[n]->size() == dof_handler.n_dofs(),
343  ExcDimensionMismatch(solutions[n]->size(), dof_handler.n_dofs()));
344 
345  Assert((coefficient == nullptr) ||
346  (coefficient->n_components == n_components) ||
347  (coefficient->n_components == 1),
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  // reserve one slot for each cell and set it to zero
360  for (unsigned int n = 0; n < n_solution_vectors; ++n)
361  (*errors[n]).reinit(dof_handler.get_triangulation().n_active_cells());
362 
363  // fields to get the gradients on the present and the neighbor cell.
364  //
365  // for the neighbor gradient, we need several auxiliary fields, depending on
366  // the way we get it (see below)
367  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
368  gradients_here(n_solution_vectors,
369  std::vector<std::vector<Tensor<1, spacedim, number>>>(
370  2,
371  std::vector<Tensor<1, spacedim, number>>(n_components)));
372  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
373  gradients_neighbor(gradients_here);
374  std::vector<Vector<typename ProductType<number, double>::type>>
375  grad_dot_n_neighbor(n_solution_vectors,
377  n_components));
378 
379  // reserve some space for coefficient values at one point. if there is no
380  // coefficient, then we fill it by unity once and for all and don't set it
381  // any more
382  Vector<double> coefficient_values(n_components);
383  if (coefficient == nullptr)
384  for (unsigned int c = 0; c < n_components; ++c)
385  coefficient_values(c) = 1;
386 
387  const QTrapezoid<1> quadrature;
388  const hp::QCollection<1> q_collection(quadrature);
389  const QGauss<0> face_quadrature(1);
390  const hp::QCollection<0> q_face_collection(face_quadrature);
391 
392  const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection();
393 
394  hp::MappingCollection<1, spacedim> mapping_collection;
395  mapping_collection.push_back(mapping);
396 
397  hp::FEValues<1, spacedim> fe_values(mapping_collection,
398  fe,
399  q_collection,
401  hp::FEFaceValues<1, spacedim> fe_face_values(
402  /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors);
403 
404  // loop over all cells and do something on the cells which we're told to
405  // work on. note that the error indicator is only a sum over the two
406  // contributions from the two vertices of each cell.
407  for (const auto &cell : dof_handler.active_cell_iterators())
409  (cell->subdomain_id() == subdomain_id)) &&
411  (cell->material_id() == material_id)))
412  {
413  for (unsigned int n = 0; n < n_solution_vectors; ++n)
414  (*errors[n])(cell->active_cell_index()) = 0;
415 
416  fe_values.reinit(cell);
417  for (unsigned int s = 0; s < n_solution_vectors; ++s)
419  *solutions[s], gradients_here[s]);
420 
421  // loop over the two points bounding this line. n==0 is left point,
422  // n==1 is right point
423  for (unsigned int n = 0; n < 2; ++n)
424  {
425  // find left or right active neighbor
426  auto neighbor = cell->neighbor(n);
427  if (neighbor.state() == IteratorState::valid)
428  while (neighbor->has_children())
429  neighbor = neighbor->child(n == 0 ? 1 : 0);
430 
431  fe_face_values.reinit(cell, n);
432  Tensor<1, spacedim> normal =
433  fe_face_values.get_present_fe_values().get_normal_vectors()[0];
434 
435  if (neighbor.state() == IteratorState::valid)
436  {
437  fe_values.reinit(neighbor);
438 
439  for (unsigned int s = 0; s < n_solution_vectors; ++s)
441  *solutions[s], gradients_neighbor[s]);
442 
443  fe_face_values.reinit(neighbor, n == 0 ? 1 : 0);
444  Tensor<1, spacedim> neighbor_normal =
445  fe_face_values.get_present_fe_values()
446  .get_normal_vectors()[0];
447 
448  // extract the gradient in normal direction of all the
449  // components.
450  for (unsigned int s = 0; s < n_solution_vectors; ++s)
451  for (unsigned int c = 0; c < n_components; ++c)
452  grad_dot_n_neighbor[s](c) =
453  -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
454  neighbor_normal);
455  }
456  else if (neumann_bc.find(n) != neumann_bc.end())
457  // if Neumann b.c., then fill the gradients field which will be
458  // used later on.
459  {
460  if (n_components == 1)
461  {
462  const typename InputVector::value_type v =
463  neumann_bc.find(n)->second->value(cell->vertex(n));
464 
465  for (unsigned int s = 0; s < n_solution_vectors; ++s)
466  grad_dot_n_neighbor[s](0) = v;
467  }
468  else
469  {
471  neumann_bc.find(n)->second->vector_value(cell->vertex(n),
472  v);
473 
474  for (unsigned int s = 0; s < n_solution_vectors; ++s)
475  grad_dot_n_neighbor[s] = v;
476  }
477  }
478  else
479  // fill with zeroes.
480  for (unsigned int s = 0; s < n_solution_vectors; ++s)
481  grad_dot_n_neighbor[s] = 0;
482 
483  // if there is a coefficient, then evaluate it at the present
484  // position. if there is none, reuse the preset values.
485  if (coefficient != nullptr)
486  {
487  if (coefficient->n_components == 1)
488  {
489  const double c_value = coefficient->value(cell->vertex(n));
490  for (unsigned int c = 0; c < n_components; ++c)
491  coefficient_values(c) = c_value;
492  }
493  else
494  coefficient->vector_value(cell->vertex(n),
495  coefficient_values);
496  }
497 
498 
499  for (unsigned int s = 0; s < n_solution_vectors; ++s)
500  for (unsigned int component = 0; component < n_components;
501  ++component)
502  if (component_mask[component] == true)
503  {
504  // get gradient here
505  const typename ProductType<number, double>::type
506  grad_dot_n_here =
507  gradients_here[s][n][component] * normal;
508 
509  const typename ProductType<number, double>::type jump =
510  ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
511  coefficient_values(component));
512  (*errors[s])(cell->active_cell_index()) +=
514  typename ProductType<number,
515  double>::type>::abs_square(jump) *
516  cell->diameter();
517  }
518  }
519 
520  for (unsigned int s = 0; s < n_solution_vectors; ++s)
521  (*errors[s])(cell->active_cell_index()) =
522  std::sqrt((*errors[s])(cell->active_cell_index()));
523  }
524 }
525 
526 
527 
528 template <int spacedim>
529 template <typename InputVector>
530 void
532  const Mapping<1, spacedim> & mapping,
533  const DoFHandler<1, spacedim> &dof_handler,
534  const Quadrature<0> & quadrature,
535  const std::map<types::boundary_id,
537  & neumann_bc,
538  const std::vector<const InputVector *> &solutions,
539  std::vector<Vector<float> *> & errors,
540  const ComponentMask & component_mask,
541  const Function<spacedim> * coefficients,
542  const unsigned int n_threads,
545  const Strategy strategy)
546 {
547  const hp::QCollection<0> quadrature_collection(quadrature);
548  estimate(mapping,
549  dof_handler,
550  quadrature_collection,
551  neumann_bc,
552  solutions,
553  errors,
554  component_mask,
555  coefficients,
556  n_threads,
557  subdomain_id,
558  material_id,
559  strategy);
560 }
561 
562 
563 
564 // explicit instantiations
565 #include "error_estimator_1d.inst"
566 
567 
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 unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Definition: fe_values.cc:3973
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type >> &gradients) const
Definition: fe_values.cc:3495
unsigned int n_components() const
const unsigned int n_components
Definition: function.h:164
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, 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)
@ cell_diameter_over_24
Kelly error estimator with the factor .
Abstract base class for mapping classes.
Definition: mapping.h:311
Definition: tensor.h:503
unsigned int n_active_cells() const
Definition: vector.h:110
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:424
const FEValuesType & get_present_fe_values() const
Definition: fe_values.h:693
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:294
void push_back(const Mapping< dim, spacedim > &new_mapping)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
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)
Definition: exceptions.h:1583
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
@ valid
Iterator points to a valid object.
const types::material_id invalid_material_id
Definition: types.h:233
const types::boundary_id internal_face_boundary_id
Definition: types.h:260
const types::subdomain_id invalid_subdomain_id
Definition: types.h:281
unsigned int subdomain_id
Definition: types.h:43
unsigned int material_id
Definition: types.h:152
unsigned int boundary_id
Definition: types.h:129
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type