Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
error_estimator_1d.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2018 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 
16 #include <deal.II/base/geometry_info.h>
17 #include <deal.II/base/quadrature.h>
18 #include <deal.II/base/quadrature_lib.h>
19 #include <deal.II/base/thread_management.h>
20 #include <deal.II/base/work_stream.h>
21 
22 #include <deal.II/distributed/tria.h>
23 
24 #include <deal.II/dofs/dof_accessor.h>
25 #include <deal.II/dofs/dof_handler.h>
26 
27 #include <deal.II/fe/fe.h>
28 #include <deal.II/fe/fe_update_flags.h>
29 #include <deal.II/fe/fe_values.h>
30 #include <deal.II/fe/mapping_q1.h>
31 
32 #include <deal.II/grid/tria_iterator.h>
33 
34 #include <deal.II/hp/fe_values.h>
35 #include <deal.II/hp/mapping_collection.h>
36 #include <deal.II/hp/q_collection.h>
37 
38 #include <deal.II/lac/block_vector.h>
39 #include <deal.II/lac/la_parallel_block_vector.h>
40 #include <deal.II/lac/la_parallel_vector.h>
41 #include <deal.II/lac/la_vector.h>
42 #include <deal.II/lac/petsc_block_vector.h>
43 #include <deal.II/lac/petsc_vector.h>
44 #include <deal.II/lac/trilinos_parallel_block_vector.h>
45 #include <deal.II/lac/trilinos_vector.h>
46 #include <deal.II/lac/vector.h>
47 
48 #include <deal.II/numerics/error_estimator.h>
49 
50 #include <algorithm>
51 #include <cmath>
52 #include <functional>
53 #include <numeric>
54 #include <vector>
55 
56 DEAL_II_NAMESPACE_OPEN
57 
58 
59 
60 template <int spacedim>
61 template <typename InputVector, typename DoFHandlerType>
62 void
64  const Mapping<1, spacedim> &mapping,
65  const DoFHandlerType & dof_handler,
66  const Quadrature<0> & quadrature,
67  const std::map<types::boundary_id,
69  & neumann_bc,
70  const InputVector & solution,
71  Vector<float> & error,
72  const ComponentMask & component_mask,
73  const Function<spacedim> *coefficients,
74  const unsigned int n_threads,
75  const types::subdomain_id subdomain_id,
76  const types::material_id material_id,
77  const Strategy strategy)
78 {
79  // just pass on to the other function
80  const std::vector<const InputVector *> solutions(1, &solution);
81  std::vector<Vector<float> *> errors(1, &error);
82  estimate(mapping,
83  dof_handler,
84  quadrature,
85  neumann_bc,
86  solutions,
87  errors,
88  component_mask,
89  coefficients,
90  n_threads,
91  subdomain_id,
92  material_id,
93  strategy);
94 }
95 
96 
97 
98 template <int spacedim>
99 template <typename InputVector, typename DoFHandlerType>
100 void
102  const DoFHandlerType &dof_handler,
103  const Quadrature<0> & quadrature,
104  const std::map<types::boundary_id,
106  & neumann_bc,
107  const InputVector & solution,
108  Vector<float> & error,
109  const ComponentMask & component_mask,
110  const Function<spacedim> *coefficients,
111  const unsigned int n_threads,
112  const types::subdomain_id subdomain_id,
113  const types::material_id material_id,
114  const Strategy strategy)
115 {
117  dof_handler,
118  quadrature,
119  neumann_bc,
120  solution,
121  error,
122  component_mask,
123  coefficients,
124  n_threads,
125  subdomain_id,
126  material_id,
127  strategy);
128 }
129 
130 
131 
132 template <int spacedim>
133 template <typename InputVector, typename DoFHandlerType>
134 void
136  const DoFHandlerType &dof_handler,
137  const Quadrature<0> & quadrature,
138  const std::map<types::boundary_id,
140  & neumann_bc,
141  const std::vector<const InputVector *> &solutions,
142  std::vector<Vector<float> *> & errors,
143  const ComponentMask & component_mask,
144  const Function<spacedim> * coefficients,
145  const unsigned int n_threads,
146  const types::subdomain_id subdomain_id,
147  const types::material_id material_id,
148  const Strategy strategy)
149 {
151  dof_handler,
152  quadrature,
153  neumann_bc,
154  solutions,
155  errors,
156  component_mask,
157  coefficients,
158  n_threads,
159  subdomain_id,
160  material_id,
161  strategy);
162 }
163 
164 
165 
166 template <int spacedim>
167 template <typename InputVector, typename DoFHandlerType>
168 void
170  const Mapping<1, spacedim> &mapping,
171  const DoFHandlerType & dof_handler,
172  const hp::QCollection<0> & quadrature,
173  const std::map<types::boundary_id,
175  & neumann_bc,
176  const InputVector & solution,
177  Vector<float> & error,
178  const ComponentMask & component_mask,
179  const Function<spacedim> *coefficients,
180  const unsigned int n_threads,
181  const types::subdomain_id subdomain_id,
182  const types::material_id material_id,
183  const Strategy strategy)
184 {
185  // just pass on to the other function
186  const std::vector<const InputVector *> solutions(1, &solution);
187  std::vector<Vector<float> *> errors(1, &error);
188  estimate(mapping,
189  dof_handler,
190  quadrature,
191  neumann_bc,
192  solutions,
193  errors,
194  component_mask,
195  coefficients,
196  n_threads,
197  subdomain_id,
198  material_id,
199  strategy);
200 }
201 
202 
203 template <int spacedim>
204 template <typename InputVector, typename DoFHandlerType>
205 void
207  const DoFHandlerType & dof_handler,
208  const hp::QCollection<0> &quadrature,
209  const std::map<types::boundary_id,
211  & neumann_bc,
212  const InputVector & solution,
213  Vector<float> & error,
214  const ComponentMask & component_mask,
215  const Function<spacedim> *coefficients,
216  const unsigned int n_threads,
217  const types::subdomain_id subdomain_id,
218  const types::material_id material_id,
219  const Strategy strategy)
220 {
222  dof_handler,
223  quadrature,
224  neumann_bc,
225  solution,
226  error,
227  component_mask,
228  coefficients,
229  n_threads,
230  subdomain_id,
231  material_id,
232  strategy);
233 }
234 
235 
236 
237 template <int spacedim>
238 template <typename InputVector, typename DoFHandlerType>
239 void
241  const DoFHandlerType & dof_handler,
242  const hp::QCollection<0> &quadrature,
243  const std::map<types::boundary_id,
245  & neumann_bc,
246  const std::vector<const InputVector *> &solutions,
247  std::vector<Vector<float> *> & errors,
248  const ComponentMask & component_mask,
249  const Function<spacedim> * coefficients,
250  const unsigned int n_threads,
251  const types::subdomain_id subdomain_id,
252  const types::material_id material_id,
253  const Strategy strategy)
254 {
256  dof_handler,
257  quadrature,
258  neumann_bc,
259  solutions,
260  errors,
261  component_mask,
262  coefficients,
263  n_threads,
264  subdomain_id,
265  material_id,
266  strategy);
267 }
268 
269 
270 
271 template <int spacedim>
272 template <typename InputVector, typename DoFHandlerType>
273 void
275  const Mapping<1, spacedim> & /*mapping*/,
276  const DoFHandlerType & /*dof_handler*/,
277  const hp::QCollection<0> &,
278  const std::map<types::boundary_id,
280  & /*neumann_bc*/,
281  const std::vector<const InputVector *> & /*solutions*/,
282  std::vector<Vector<float> *> & /*errors*/,
283  const ComponentMask & /*component_mask_*/,
284  const Function<spacedim> * /*coefficient*/,
285  const unsigned int,
286  const types::subdomain_id /*subdomain_id*/,
287  const types::material_id /*material_id*/,
288  const Strategy /*strategy*/)
289 {
290  Assert(false, ExcInternalError());
291 }
292 
293 
294 
295 template <int spacedim>
296 template <typename InputVector, typename DoFHandlerType>
297 void
299  const Mapping<1, spacedim> &mapping,
300  const DoFHandlerType & dof_handler,
301  const Quadrature<0> &,
302  const std::map<types::boundary_id,
304  & neumann_bc,
305  const std::vector<const InputVector *> &solutions,
306  std::vector<Vector<float> *> & errors,
307  const ComponentMask & component_mask,
308  const Function<spacedim> * coefficient,
309  const unsigned int,
310  const types::subdomain_id subdomain_id_,
311  const types::material_id material_id,
312  const Strategy strategy)
313 {
315  using number = typename InputVector::value_type;
316 #ifdef DEAL_II_WITH_P4EST
317  if (dynamic_cast<const parallel::distributed::Triangulation<1, spacedim> *>(
318  &dof_handler.get_triangulation()) != nullptr)
319  Assert((subdomain_id_ == numbers::invalid_subdomain_id) ||
320  (subdomain_id_ ==
321  dynamic_cast<
323  dof_handler.get_triangulation())
324  .locally_owned_subdomain()),
325  ExcMessage(
326  "For parallel distributed triangulations, the only "
327  "valid subdomain_id that can be passed here is the "
328  "one that corresponds to the locally owned subdomain id."));
329 
330  const types::subdomain_id subdomain_id =
332  &dof_handler.get_triangulation()) != nullptr) ?
334  dof_handler.get_triangulation())
335  .locally_owned_subdomain() :
336  subdomain_id_);
337 #else
338  const types::subdomain_id subdomain_id = subdomain_id_;
339 #endif
340 
341  const unsigned int n_components = dof_handler.get_fe(0).n_components();
342  const unsigned int n_solution_vectors = solutions.size();
343 
344  // sanity checks
345  Assert(neumann_bc.find(numbers::internal_face_boundary_id) ==
346  neumann_bc.end(),
347  ExcMessage("You are not allowed to list the special boundary "
348  "indicator for internal boundaries in your boundary "
349  "value map."));
350 
351  for (const auto &boundary_function : neumann_bc)
352  {
353  (void)boundary_function;
354  Assert(boundary_function.second->n_components == n_components,
355  ExcInvalidBoundaryFunction(boundary_function.first,
356  boundary_function.second->n_components,
357  n_components));
358  }
359 
360  Assert(component_mask.represents_n_components(n_components),
362  Assert(component_mask.n_selected_components(n_components) > 0,
364 
365  Assert((coefficient == nullptr) ||
366  (coefficient->n_components == n_components) ||
367  (coefficient->n_components == 1),
369 
370  Assert(solutions.size() > 0, ExcNoSolutions());
371  Assert(solutions.size() == errors.size(),
372  ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
373  for (unsigned int n = 0; n < solutions.size(); ++n)
374  Assert(solutions[n]->size() == dof_handler.n_dofs(),
375  ExcDimensionMismatch(solutions[n]->size(), dof_handler.n_dofs()));
376 
377  Assert((coefficient == nullptr) ||
378  (coefficient->n_components == n_components) ||
379  (coefficient->n_components == 1),
381 
382  for (const auto &boundary_function : neumann_bc)
383  {
384  (void)boundary_function;
385  Assert(boundary_function.second->n_components == n_components,
386  ExcInvalidBoundaryFunction(boundary_function.first,
387  boundary_function.second->n_components,
388  n_components));
389  }
390 
391  // reserve one slot for each cell and set it to zero
392  for (unsigned int n = 0; n < n_solution_vectors; ++n)
393  (*errors[n]).reinit(dof_handler.get_triangulation().n_active_cells());
394 
395  // fields to get the gradients on the present and the neighbor cell.
396  //
397  // for the neighbor gradient, we need several auxiliary fields, depending on
398  // the way we get it (see below)
399  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
400  gradients_here(n_solution_vectors,
401  std::vector<std::vector<Tensor<1, spacedim, number>>>(
402  2,
403  std::vector<Tensor<1, spacedim, number>>(n_components)));
404  std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
405  gradients_neighbor(gradients_here);
406  std::vector<Vector<typename ProductType<number, double>::type>>
407  grad_dot_n_neighbor(n_solution_vectors,
408  Vector<typename ProductType<number, double>::type>(
409  n_components));
410 
411  // reserve some space for coefficient values at one point. if there is no
412  // coefficient, then we fill it by unity once and for all and don't set it
413  // any more
414  Vector<double> coefficient_values(n_components);
415  if (coefficient == nullptr)
416  for (unsigned int c = 0; c < n_components; ++c)
417  coefficient_values(c) = 1;
418 
419  const QTrapez<1> quadrature;
420  const hp::QCollection<1> q_collection(quadrature);
421  const QGauss<0> face_quadrature(1);
422  const hp::QCollection<0> q_face_collection(face_quadrature);
423 
424  const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection();
425 
426  hp::MappingCollection<1, spacedim> mapping_collection;
427  mapping_collection.push_back(mapping);
428 
429  hp::FEValues<1, spacedim> fe_values(mapping_collection,
430  fe,
431  q_collection,
433  hp::FEFaceValues<1, spacedim> fe_face_values(
434  /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors);
435 
436  // loop over all cells and do something on the cells which we're told to
437  // work on. note that the error indicator is only a sum over the two
438  // contributions from the two vertices of each cell.
439  for (typename DoFHandlerType::active_cell_iterator cell =
440  dof_handler.begin_active();
441  cell != dof_handler.end();
442  ++cell)
443  if (((subdomain_id == numbers::invalid_subdomain_id) ||
444  (cell->subdomain_id() == subdomain_id)) &&
445  ((material_id == numbers::invalid_material_id) ||
446  (cell->material_id() == material_id)))
447  {
448  for (unsigned int n = 0; n < n_solution_vectors; ++n)
449  (*errors[n])(cell->active_cell_index()) = 0;
450 
451  fe_values.reinit(cell);
452  for (unsigned int s = 0; s < n_solution_vectors; ++s)
453  fe_values.get_present_fe_values().get_function_gradients(
454  *solutions[s], gradients_here[s]);
455 
456  // loop over the two points bounding this line. n==0 is left point,
457  // n==1 is right point
458  for (unsigned int n = 0; n < 2; ++n)
459  {
460  // find left or right active neighbor
461  typename DoFHandlerType::cell_iterator neighbor = cell->neighbor(n);
462  if (neighbor.state() == IteratorState::valid)
463  while (neighbor->has_children())
464  neighbor = neighbor->child(n == 0 ? 1 : 0);
465 
466  fe_face_values.reinit(cell, n);
467  Tensor<1, spacedim> normal = fe_face_values.get_present_fe_values()
468  .get_all_normal_vectors()[0];
469 
470  if (neighbor.state() == IteratorState::valid)
471  {
472  fe_values.reinit(neighbor);
473 
474  for (unsigned int s = 0; s < n_solution_vectors; ++s)
475  fe_values.get_present_fe_values().get_function_gradients(
476  *solutions[s], gradients_neighbor[s]);
477 
478  fe_face_values.reinit(neighbor, n == 0 ? 1 : 0);
479  Tensor<1, spacedim> neighbor_normal =
480  fe_face_values.get_present_fe_values()
481  .get_all_normal_vectors()[0];
482 
483  // extract the gradient in normal direction of all the
484  // components.
485  for (unsigned int s = 0; s < n_solution_vectors; ++s)
486  for (unsigned int c = 0; c < n_components; ++c)
487  grad_dot_n_neighbor[s](c) =
488  -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
489  neighbor_normal);
490  }
491  else if (neumann_bc.find(n) != neumann_bc.end())
492  // if Neumann b.c., then fill the gradients field which will be
493  // used later on.
494  {
495  if (n_components == 1)
496  {
497  const typename InputVector::value_type v =
498  neumann_bc.find(n)->second->value(cell->vertex(n));
499 
500  for (unsigned int s = 0; s < n_solution_vectors; ++s)
501  grad_dot_n_neighbor[s](0) = v;
502  }
503  else
504  {
505  Vector<typename InputVector::value_type> v(n_components);
506  neumann_bc.find(n)->second->vector_value(cell->vertex(n),
507  v);
508 
509  for (unsigned int s = 0; s < n_solution_vectors; ++s)
510  grad_dot_n_neighbor[s] = v;
511  }
512  }
513  else
514  // fill with zeroes.
515  for (unsigned int s = 0; s < n_solution_vectors; ++s)
516  grad_dot_n_neighbor[s] = 0;
517 
518  // if there is a coefficient, then evaluate it at the present
519  // position. if there is none, reuse the preset values.
520  if (coefficient != nullptr)
521  {
522  if (coefficient->n_components == 1)
523  {
524  const double c_value = coefficient->value(cell->vertex(n));
525  for (unsigned int c = 0; c < n_components; ++c)
526  coefficient_values(c) = c_value;
527  }
528  else
529  coefficient->vector_value(cell->vertex(n),
530  coefficient_values);
531  }
532 
533 
534  for (unsigned int s = 0; s < n_solution_vectors; ++s)
535  for (unsigned int component = 0; component < n_components;
536  ++component)
537  if (component_mask[component] == true)
538  {
539  // get gradient here
540  const typename ProductType<number, double>::type
541  grad_dot_n_here =
542  gradients_here[s][n][component] * normal;
543 
544  const typename ProductType<number, double>::type jump =
545  ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
546  coefficient_values(component));
547  (*errors[s])(cell->active_cell_index()) +=
549  typename ProductType<number,
550  double>::type>::abs_square(jump) *
551  cell->diameter();
552  }
553  }
554 
555  for (unsigned int s = 0; s < n_solution_vectors; ++s)
556  (*errors[s])(cell->active_cell_index()) =
557  std::sqrt((*errors[s])(cell->active_cell_index()));
558  }
559 }
560 
561 
562 // explicit instantiations
563 #include "error_estimator_1d.inst"
564 
565 
566 DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInvalidBoundaryFunction(types::boundary_id arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidCoefficient()
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandlerType &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)
const types::subdomain_id invalid_subdomain_id
Definition: types.h:258
const unsigned int n_components
Definition: function.h:162
static ::ExceptionBase & ExcInvalidComponentMask()
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:134
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
LinearAlgebra::distributed::Vector< Number > Vector
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: dof_tools.h:57
static ::ExceptionBase & ExcIncompatibleNumberOfElements(int arg1, int arg2)
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, 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:142
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
Definition: mpi.h:90
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:223
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:196
unsigned int boundary_id
Definition: types.h:111
static ::ExceptionBase & ExcNoSolutions()
static ::ExceptionBase & ExcInternalError()
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, 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:265