Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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>
30
32
36
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
58template <int spacedim>
59template <typename InputVector>
60void
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,
73 const types::subdomain_id subdomain_id,
74 const types::material_id material_id,
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,
89 subdomain_id,
90 material_id,
91 strategy);
92}
93
94
95
96template <int spacedim>
97template <typename InputVector>
98void
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,
110 const types::subdomain_id subdomain_id,
111 const types::material_id material_id,
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
130template <int spacedim>
131template <typename InputVector>
132void
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,
144 const types::subdomain_id subdomain_id,
145 const types::material_id material_id,
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
164template <int spacedim>
165template <typename InputVector>
166void
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,
179 const types::subdomain_id subdomain_id,
180 const types::material_id material_id,
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
202template <int spacedim>
203template <typename InputVector>
204void
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,
216 const types::subdomain_id subdomain_id,
217 const types::material_id material_id,
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
236template <int spacedim>
237template <typename InputVector>
238void
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,
250 const types::subdomain_id subdomain_id,
251 const types::material_id material_id,
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
270template <int spacedim>
271template <typename InputVector>
272void
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 AssertThrow(strategy == cell_diameter_over_24, ExcNotImplemented());
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()),
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
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),
329 ExcInvalidComponentMask());
330 Assert(component_mask.n_selected_components(n_components) > 0,
331 ExcInvalidComponentMask());
332
333 Assert((coefficient == nullptr) ||
334 (coefficient->n_components == n_components) ||
335 (coefficient->n_components == 1),
336 ExcInvalidCoefficient());
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),
348 ExcInvalidCoefficient());
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())
408 if (((subdomain_id == numbers::invalid_subdomain_id) ||
409 (cell->subdomain_id() == subdomain_id)) &&
410 ((material_id == numbers::invalid_material_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)
418 fe_values.get_present_fe_values().get_function_gradients(
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)
440 fe_values.get_present_fe_values().get_function_gradients(
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
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
528template <int spacedim>
529template <typename InputVector>
530void
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,
543 const types::subdomain_id subdomain_id,
544 const types::material_id material_id,
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 FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
const hp::FECollection< dim, spacedim > & get_fe_collection() const
const Triangulation< dim, spacedim > & get_triangulation() const
types::global_dof_index n_dofs() 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, 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)
Abstract base class for mapping classes.
Definition: mapping.h:311
Definition: tensor.h:503
unsigned int n_active_cells() const
Definition: vector.h:109
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
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
@ 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
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
typename internal::ProductTypeImpl< typename std::decay< T >::type, typename std::decay< U >::type >::type type