Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10: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 - 2022 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.h>
30 
32 
33 #include <deal.II/hp/fe_values.h>
36 
44 #include <deal.II/lac/vector.h>
45 
47 
48 #include <algorithm>
49 #include <cmath>
50 #include <functional>
51 #include <numeric>
52 #include <vector>
53 
55 
56 
57 template <int spacedim>
58 template <typename Number>
59 void
61  const Mapping<1, spacedim> &mapping,
62  const DoFHandler<1, spacedim> &dof_handler,
63  const Quadrature<0> &quadrature,
64  const std::map<types::boundary_id, const Function<spacedim, Number> *>
65  &neumann_bc,
66  const ReadVector<Number> &solution,
67  Vector<float> &error,
68  const ComponentMask &component_mask,
69  const Function<spacedim> *coefficients,
70  const unsigned int n_threads,
73  const Strategy strategy)
74 {
75  // just pass on to the other function
76  std::vector<const ReadVector<Number> *> solutions(1, &solution);
77  std::vector<Vector<float> *> errors(1, &error);
78  ArrayView<Vector<float> *> error_view = make_array_view(errors);
79  estimate(mapping,
80  dof_handler,
81  quadrature,
82  neumann_bc,
83  make_array_view(solutions),
84  error_view,
85  component_mask,
86  coefficients,
87  n_threads,
90  strategy);
91 }
92 
93 
94 
95 template <int spacedim>
96 template <typename Number>
97 void
99  const DoFHandler<1, spacedim> &dof_handler,
100  const Quadrature<0> &quadrature,
101  const std::map<types::boundary_id, const Function<spacedim, Number> *>
102  &neumann_bc,
103  const ReadVector<Number> &solution,
104  Vector<float> &error,
105  const ComponentMask &component_mask,
106  const Function<spacedim> *coefficients,
107  const unsigned int n_threads,
110  const Strategy strategy)
111 {
113  estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
114  dof_handler,
115  quadrature,
116  neumann_bc,
117  solution,
118  error,
119  component_mask,
120  coefficients,
121  n_threads,
122  subdomain_id,
123  material_id,
124  strategy);
125 }
126 
127 
128 
129 template <int spacedim>
130 template <typename Number>
131 void
133  const DoFHandler<1, spacedim> &dof_handler,
134  const Quadrature<0> &quadrature,
135  const std::map<types::boundary_id, const Function<spacedim, Number> *>
136  &neumann_bc,
137  const ArrayView<const ReadVector<Number> *> &solutions,
138  ArrayView<Vector<float> *> &errors,
139  const ComponentMask &component_mask,
140  const Function<spacedim> *coefficients,
141  const unsigned int n_threads,
144  const Strategy strategy)
145 {
147  estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
148  dof_handler,
149  quadrature,
150  neumann_bc,
151  solutions,
152  errors,
153  component_mask,
154  coefficients,
155  n_threads,
156  subdomain_id,
157  material_id,
158  strategy);
159 }
160 
161 
162 
163 template <int spacedim>
164 template <typename Number>
165 void
167  const Mapping<1, spacedim> &mapping,
168  const DoFHandler<1, spacedim> &dof_handler,
169  const hp::QCollection<0> &quadrature,
170  const std::map<types::boundary_id, const Function<spacedim, Number> *>
171  &neumann_bc,
172  const ReadVector<Number> &solution,
173  Vector<float> &error,
174  const ComponentMask &component_mask,
175  const Function<spacedim> *coefficients,
176  const unsigned int n_threads,
179  const Strategy strategy)
180 {
181  // just pass on to the other function
182  std::vector<const ReadVector<Number> *> solutions(1, &solution);
183  std::vector<Vector<float> *> errors(1, &error);
184  ArrayView<Vector<float> *> error_view = make_array_view(errors);
185  estimate(mapping,
186  dof_handler,
187  quadrature,
188  neumann_bc,
189  make_array_view(solutions),
190  error_view,
191  component_mask,
192  coefficients,
193  n_threads,
194  subdomain_id,
195  material_id,
196  strategy);
197 }
198 
199 
200 
201 template <int spacedim>
202 template <typename Number>
203 void
205  const DoFHandler<1, spacedim> &dof_handler,
206  const hp::QCollection<0> &quadrature,
207  const std::map<types::boundary_id, const Function<spacedim, Number> *>
208  &neumann_bc,
209  const ReadVector<Number> &solution,
210  Vector<float> &error,
211  const ComponentMask &component_mask,
212  const Function<spacedim> *coefficients,
213  const unsigned int n_threads,
216  const Strategy strategy)
217 {
219  estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
220  dof_handler,
221  quadrature,
222  neumann_bc,
223  solution,
224  error,
225  component_mask,
226  coefficients,
227  n_threads,
228  subdomain_id,
229  material_id,
230  strategy);
231 }
232 
233 
234 
235 template <int spacedim>
236 template <typename Number>
237 void
239  const DoFHandler<1, spacedim> &dof_handler,
240  const hp::QCollection<0> &quadrature,
241  const std::map<types::boundary_id, const Function<spacedim, Number> *>
242  &neumann_bc,
243  const ArrayView<const ReadVector<Number> *> &solutions,
244  ArrayView<Vector<float> *> &errors,
245  const ComponentMask &component_mask,
246  const Function<spacedim> *coefficients,
247  const unsigned int n_threads,
250  const Strategy strategy)
251 {
253  estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
254  dof_handler,
255  quadrature,
256  neumann_bc,
257  solutions,
258  errors,
259  component_mask,
260  coefficients,
261  n_threads,
262  subdomain_id,
263  material_id,
264  strategy);
265 }
266 
267 
268 
269 template <int spacedim>
270 template <typename Number>
271 void
273  const Mapping<1, spacedim> &mapping,
274  const DoFHandler<1, spacedim> &dof_handler,
275  const hp::QCollection<0> &,
276  const std::map<types::boundary_id, const Function<spacedim, Number> *>
277  &neumann_bc,
278  const ArrayView<const ReadVector<Number> *> &solutions,
279  ArrayView<Vector<float> *> &errors,
280  const ComponentMask &component_mask,
281  const Function<spacedim> *coefficient,
282  const unsigned int,
283  const types::subdomain_id subdomain_id_,
285  const Strategy strategy)
286 {
288  using number = Number;
290  if (const auto *triangulation = dynamic_cast<
292  &dof_handler.get_triangulation()))
293  {
294  Assert((subdomain_id_ == numbers::invalid_subdomain_id) ||
295  (subdomain_id_ == triangulation->locally_owned_subdomain()),
296  ExcMessage(
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."));
300  subdomain_id = triangulation->locally_owned_subdomain();
301  }
302  else
303  {
304  subdomain_id = subdomain_id_;
305  }
306 
307  const unsigned int n_components = dof_handler.get_fe(0).n_components();
308  const unsigned int n_solution_vectors = solutions.size();
309 
310  // sanity checks
311  Assert(neumann_bc.find(numbers::internal_face_boundary_id) ==
312  neumann_bc.end(),
313  ExcMessage("You are not allowed to list the special boundary "
314  "indicator for internal boundaries in your boundary "
315  "value map."));
316 
317  for (const auto &boundary_function : neumann_bc)
318  {
319  (void)boundary_function;
320  Assert(boundary_function.second->n_components == n_components,
321  ExcInvalidBoundaryFunction(boundary_function.first,
322  boundary_function.second->n_components,
323  n_components));
324  }
325 
326  Assert(component_mask.represents_n_components(n_components),
328  Assert(component_mask.n_selected_components(n_components) > 0,
330 
331  Assert((coefficient == nullptr) ||
332  (coefficient->n_components == n_components) ||
333  (coefficient->n_components == 1),
335 
336  Assert(solutions.size() > 0, ExcNoSolutions());
337  Assert(solutions.size() == errors.size(),
338  ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
339  for (unsigned int n = 0; n < solutions.size(); ++n)
340  Assert(solutions[n]->size() == dof_handler.n_dofs(),
341  ExcDimensionMismatch(solutions[n]->size(), dof_handler.n_dofs()));
342 
343  Assert((coefficient == nullptr) ||
344  (coefficient->n_components == n_components) ||
345  (coefficient->n_components == 1),
347 
348  for (const auto &boundary_function : neumann_bc)
349  {
350  (void)boundary_function;
351  Assert(boundary_function.second->n_components == n_components,
352  ExcInvalidBoundaryFunction(boundary_function.first,
353  boundary_function.second->n_components,
354  n_components));
355  }
356 
357  // reserve one slot for each cell and set it to zero
358  for (unsigned int n = 0; n < n_solution_vectors; ++n)
359  (*errors[n]).reinit(dof_handler.get_triangulation().n_active_cells());
360 
361  // fields to get the gradients on the present and the neighbor cell.
362  //
363  // for the neighbor gradient, we need several auxiliary fields, depending on
364  // the way we get it (see below)
365  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
366  gradients_here(n_solution_vectors,
367  std::vector<std::vector<Tensor<1, spacedim, number>>>(
368  2,
369  std::vector<Tensor<1, spacedim, number>>(n_components)));
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,
375  n_components));
376 
377  // reserve some space for coefficient values at one point. if there is no
378  // coefficient, then we fill it by unity once and for all and don't set it
379  // any more
380  Vector<double> coefficient_values(n_components);
381  if (coefficient == nullptr)
382  for (unsigned int c = 0; c < n_components; ++c)
383  coefficient_values(c) = 1;
384 
385  const QTrapezoid<1> quadrature;
386  const hp::QCollection<1> q_collection(quadrature);
387  const QGauss<0> face_quadrature(1);
388  const hp::QCollection<0> q_face_collection(face_quadrature);
389 
390  const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection();
391 
392  hp::MappingCollection<1, spacedim> mapping_collection;
393  mapping_collection.push_back(mapping);
394 
395  hp::FEValues<1, spacedim> fe_values(mapping_collection,
396  fe,
397  q_collection,
399  hp::FEFaceValues<1, spacedim> fe_face_values(
400  /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors);
401 
402  // loop over all cells and do something on the cells which we're told to
403  // work on. note that the error indicator is only a sum over the two
404  // contributions from the two vertices of each cell.
405  for (const auto &cell : dof_handler.active_cell_iterators())
407  (cell->subdomain_id() == subdomain_id)) &&
409  (cell->material_id() == material_id)))
410  {
411  for (unsigned int n = 0; n < n_solution_vectors; ++n)
412  (*errors[n])(cell->active_cell_index()) = 0;
413 
414  fe_values.reinit(cell);
415  for (unsigned int s = 0; s < n_solution_vectors; ++s)
417  *solutions[s], gradients_here[s]);
418 
419  // loop over the two points bounding this line. n==0 is left point,
420  // n==1 is right point
421  for (unsigned int n = 0; n < 2; ++n)
422  {
423  // find left or right active neighbor
424  auto neighbor = cell->neighbor(n);
425  if (neighbor.state() == IteratorState::valid)
426  while (neighbor->has_children())
427  neighbor = neighbor->child(n == 0 ? 1 : 0);
428 
429  fe_face_values.reinit(cell, n);
430  Tensor<1, spacedim> normal =
431  fe_face_values.get_present_fe_values().get_normal_vectors()[0];
432 
433  if (neighbor.state() == IteratorState::valid)
434  {
435  fe_values.reinit(neighbor);
436 
437  for (unsigned int s = 0; s < n_solution_vectors; ++s)
439  *solutions[s], gradients_neighbor[s]);
440 
441  fe_face_values.reinit(neighbor, n == 0 ? 1 : 0);
442  Tensor<1, spacedim> neighbor_normal =
443  fe_face_values.get_present_fe_values()
444  .get_normal_vectors()[0];
445 
446  // extract the gradient in normal direction of all the
447  // components.
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] *
452  neighbor_normal);
453  }
454  else if (neumann_bc.find(n) != neumann_bc.end())
455  // if Neumann b.c., then fill the gradients field which will be
456  // used later on.
457  {
458  if (n_components == 1)
459  {
460  const Number v =
461  neumann_bc.find(n)->second->value(cell->vertex(n));
462 
463  for (unsigned int s = 0; s < n_solution_vectors; ++s)
464  grad_dot_n_neighbor[s](0) = v;
465  }
466  else
467  {
468  Vector<Number> v(n_components);
469  neumann_bc.find(n)->second->vector_value(cell->vertex(n),
470  v);
471 
472  for (unsigned int s = 0; s < n_solution_vectors; ++s)
473  grad_dot_n_neighbor[s] = v;
474  }
475  }
476  else
477  // fill with zeroes.
478  for (unsigned int s = 0; s < n_solution_vectors; ++s)
479  grad_dot_n_neighbor[s] = 0;
480 
481  // if there is a coefficient, then evaluate it at the present
482  // position. if there is none, reuse the preset values.
483  if (coefficient != nullptr)
484  {
485  if (coefficient->n_components == 1)
486  {
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;
490  }
491  else
492  coefficient->vector_value(cell->vertex(n),
493  coefficient_values);
494  }
495 
496 
497  for (unsigned int s = 0; s < n_solution_vectors; ++s)
498  for (unsigned int component = 0; component < n_components;
499  ++component)
500  if (component_mask[component] == true)
501  {
502  // get gradient here
503  const typename ProductType<number, double>::type
504  grad_dot_n_here =
505  gradients_here[s][n][component] * normal;
506 
507  const typename ProductType<number, double>::type jump =
508  ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
509  coefficient_values(component));
510  (*errors[s])(cell->active_cell_index()) +=
512  typename ProductType<number,
513  double>::type>::abs_square(jump) *
514  cell->diameter();
515  }
516  }
517 
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()));
521  }
522 }
523 
524 
525 
526 template <int spacedim>
527 template <typename Number>
528 void
530  const Mapping<1, spacedim> &mapping,
531  const DoFHandler<1, spacedim> &dof_handler,
532  const Quadrature<0> &quadrature,
533  const std::map<types::boundary_id, const Function<spacedim, Number> *>
534  &neumann_bc,
535  const ArrayView<const ReadVector<Number> *> &solutions,
536  ArrayView<Vector<float> *> &errors,
537  const ComponentMask &component_mask,
538  const Function<spacedim> *coefficients,
539  const unsigned int n_threads,
542  const Strategy strategy)
543 {
544  const hp::QCollection<0> quadrature_collection(quadrature);
545  estimate(mapping,
546  dof_handler,
547  quadrature_collection,
548  neumann_bc,
549  solutions,
550  errors,
551  component_mask,
552  coefficients,
553  n_threads,
554  subdomain_id,
555  material_id,
556  strategy);
557 }
558 
559 
560 
561 // explicit instantiations
562 #include "error_estimator_1d.inst"
563 
564 
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition: array_view.h:950
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
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, 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.
Definition: mapping.h:317
Definition: tensor.h:516
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:695
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
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:1631
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:1732
@ 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
Definition: types.h:278
const types::boundary_id internal_face_boundary_id
Definition: types.h:313
const types::subdomain_id invalid_subdomain_id
Definition: types.h:342
unsigned int subdomain_id
Definition: types.h:44
unsigned int material_id
Definition: types.h:168
unsigned int boundary_id
Definition: types.h:145
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type