Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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>
29#include <deal.II/fe/mapping.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{
114 const auto reference_cell = ReferenceCells::Line;
115 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
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
131template <int spacedim>
132template <typename InputVector>
133void
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,
145 const types::subdomain_id subdomain_id,
146 const types::material_id material_id,
147 const Strategy strategy)
148{
149 const auto reference_cell = ReferenceCells::Line;
150 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
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
166template <int spacedim>
167template <typename InputVector>
168void
170 const Mapping<1, spacedim> & mapping,
171 const DoFHandler<1, spacedim> &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
204template <int spacedim>
205template <typename InputVector>
206void
208 const DoFHandler<1, spacedim> &dof_handler,
209 const hp::QCollection<0> & quadrature,
210 const std::map<types::boundary_id,
212 & neumann_bc,
213 const InputVector & solution,
214 Vector<float> & error,
215 const ComponentMask & component_mask,
216 const Function<spacedim> *coefficients,
217 const unsigned int n_threads,
218 const types::subdomain_id subdomain_id,
219 const types::material_id material_id,
220 const Strategy strategy)
221{
222 const auto reference_cell = ReferenceCells::Line;
223 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
224 dof_handler,
225 quadrature,
226 neumann_bc,
227 solution,
228 error,
229 component_mask,
230 coefficients,
231 n_threads,
232 subdomain_id,
233 material_id,
234 strategy);
235}
236
237
238
239template <int spacedim>
240template <typename InputVector>
241void
243 const DoFHandler<1, spacedim> &dof_handler,
244 const hp::QCollection<0> & quadrature,
245 const std::map<types::boundary_id,
247 & neumann_bc,
248 const std::vector<const InputVector *> &solutions,
249 std::vector<Vector<float> *> & errors,
250 const ComponentMask & component_mask,
251 const Function<spacedim> * coefficients,
252 const unsigned int n_threads,
253 const types::subdomain_id subdomain_id,
254 const types::material_id material_id,
255 const Strategy strategy)
256{
257 const auto reference_cell = ReferenceCells::Line;
258 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
259 dof_handler,
260 quadrature,
261 neumann_bc,
262 solutions,
263 errors,
264 component_mask,
265 coefficients,
266 n_threads,
267 subdomain_id,
268 material_id,
269 strategy);
270}
271
272
273
274template <int spacedim>
275template <typename InputVector>
276void
278 const Mapping<1, spacedim> & mapping,
279 const DoFHandler<1, spacedim> &dof_handler,
280 const hp::QCollection<0> &,
281 const std::map<types::boundary_id,
283 & neumann_bc,
284 const std::vector<const InputVector *> &solutions,
285 std::vector<Vector<float> *> & errors,
286 const ComponentMask & component_mask,
287 const Function<spacedim> * coefficient,
288 const unsigned int,
289 const types::subdomain_id subdomain_id_,
290 const types::material_id material_id,
291 const Strategy strategy)
292{
293 AssertThrow(strategy == cell_diameter_over_24, ExcNotImplemented());
294 using number = typename InputVector::value_type;
296 if (const auto *triangulation = dynamic_cast<
298 &dof_handler.get_triangulation()))
299 {
300 Assert((subdomain_id_ == numbers::invalid_subdomain_id) ||
301 (subdomain_id_ == triangulation->locally_owned_subdomain()),
303 "For distributed Triangulation objects and associated "
304 "DoFHandler objects, asking for any subdomain other than the "
305 "locally owned one does not make sense."));
306 subdomain_id = triangulation->locally_owned_subdomain();
307 }
308 else
309 {
310 subdomain_id = subdomain_id_;
311 }
312
313 const unsigned int n_components = dof_handler.get_fe(0).n_components();
314 const unsigned int n_solution_vectors = solutions.size();
315
316 // sanity checks
318 neumann_bc.end(),
319 ExcMessage("You are not allowed to list the special boundary "
320 "indicator for internal boundaries in your boundary "
321 "value map."));
322
323 for (const auto &boundary_function : neumann_bc)
324 {
325 (void)boundary_function;
326 Assert(boundary_function.second->n_components == n_components,
327 ExcInvalidBoundaryFunction(boundary_function.first,
328 boundary_function.second->n_components,
329 n_components));
330 }
331
332 Assert(component_mask.represents_n_components(n_components),
333 ExcInvalidComponentMask());
334 Assert(component_mask.n_selected_components(n_components) > 0,
335 ExcInvalidComponentMask());
336
337 Assert((coefficient == nullptr) ||
338 (coefficient->n_components == n_components) ||
339 (coefficient->n_components == 1),
340 ExcInvalidCoefficient());
341
342 Assert(solutions.size() > 0, ExcNoSolutions());
343 Assert(solutions.size() == errors.size(),
344 ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
345 for (unsigned int n = 0; n < solutions.size(); ++n)
346 Assert(solutions[n]->size() == dof_handler.n_dofs(),
347 ExcDimensionMismatch(solutions[n]->size(), dof_handler.n_dofs()));
348
349 Assert((coefficient == nullptr) ||
350 (coefficient->n_components == n_components) ||
351 (coefficient->n_components == 1),
352 ExcInvalidCoefficient());
353
354 for (const auto &boundary_function : neumann_bc)
355 {
356 (void)boundary_function;
357 Assert(boundary_function.second->n_components == n_components,
358 ExcInvalidBoundaryFunction(boundary_function.first,
359 boundary_function.second->n_components,
360 n_components));
361 }
362
363 // reserve one slot for each cell and set it to zero
364 for (unsigned int n = 0; n < n_solution_vectors; ++n)
365 (*errors[n]).reinit(dof_handler.get_triangulation().n_active_cells());
366
367 // fields to get the gradients on the present and the neighbor cell.
368 //
369 // for the neighbor gradient, we need several auxiliary fields, depending on
370 // the way we get it (see below)
371 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
372 gradients_here(n_solution_vectors,
373 std::vector<std::vector<Tensor<1, spacedim, number>>>(
374 2,
375 std::vector<Tensor<1, spacedim, number>>(n_components)));
376 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
377 gradients_neighbor(gradients_here);
378 std::vector<Vector<typename ProductType<number, double>::type>>
379 grad_dot_n_neighbor(n_solution_vectors,
381 n_components));
382
383 // reserve some space for coefficient values at one point. if there is no
384 // coefficient, then we fill it by unity once and for all and don't set it
385 // any more
386 Vector<double> coefficient_values(n_components);
387 if (coefficient == nullptr)
388 for (unsigned int c = 0; c < n_components; ++c)
389 coefficient_values(c) = 1;
390
391 const QTrapezoid<1> quadrature;
392 const hp::QCollection<1> q_collection(quadrature);
393 const QGauss<0> face_quadrature(1);
394 const hp::QCollection<0> q_face_collection(face_quadrature);
395
396 const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection();
397
398 hp::MappingCollection<1, spacedim> mapping_collection;
399 mapping_collection.push_back(mapping);
400
401 hp::FEValues<1, spacedim> fe_values(mapping_collection,
402 fe,
403 q_collection,
405 hp::FEFaceValues<1, spacedim> fe_face_values(
406 /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors);
407
408 // loop over all cells and do something on the cells which we're told to
409 // work on. note that the error indicator is only a sum over the two
410 // contributions from the two vertices of each cell.
411 for (const auto &cell : dof_handler.active_cell_iterators())
412 if (((subdomain_id == numbers::invalid_subdomain_id) ||
413 (cell->subdomain_id() == subdomain_id)) &&
414 ((material_id == numbers::invalid_material_id) ||
415 (cell->material_id() == material_id)))
416 {
417 for (unsigned int n = 0; n < n_solution_vectors; ++n)
418 (*errors[n])(cell->active_cell_index()) = 0;
419
420 fe_values.reinit(cell);
421 for (unsigned int s = 0; s < n_solution_vectors; ++s)
422 fe_values.get_present_fe_values().get_function_gradients(
423 *solutions[s], gradients_here[s]);
424
425 // loop over the two points bounding this line. n==0 is left point,
426 // n==1 is right point
427 for (unsigned int n = 0; n < 2; ++n)
428 {
429 // find left or right active neighbor
430 auto neighbor = cell->neighbor(n);
431 if (neighbor.state() == IteratorState::valid)
432 while (neighbor->has_children())
433 neighbor = neighbor->child(n == 0 ? 1 : 0);
434
435 fe_face_values.reinit(cell, n);
436 Tensor<1, spacedim> normal =
437 fe_face_values.get_present_fe_values().get_normal_vectors()[0];
438
439 if (neighbor.state() == IteratorState::valid)
440 {
441 fe_values.reinit(neighbor);
442
443 for (unsigned int s = 0; s < n_solution_vectors; ++s)
444 fe_values.get_present_fe_values().get_function_gradients(
445 *solutions[s], gradients_neighbor[s]);
446
447 fe_face_values.reinit(neighbor, n == 0 ? 1 : 0);
448 Tensor<1, spacedim> neighbor_normal =
449 fe_face_values.get_present_fe_values()
450 .get_normal_vectors()[0];
451
452 // extract the gradient in normal direction of all the
453 // components.
454 for (unsigned int s = 0; s < n_solution_vectors; ++s)
455 for (unsigned int c = 0; c < n_components; ++c)
456 grad_dot_n_neighbor[s](c) =
457 -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
458 neighbor_normal);
459 }
460 else if (neumann_bc.find(n) != neumann_bc.end())
461 // if Neumann b.c., then fill the gradients field which will be
462 // used later on.
463 {
464 if (n_components == 1)
465 {
466 const typename InputVector::value_type v =
467 neumann_bc.find(n)->second->value(cell->vertex(n));
468
469 for (unsigned int s = 0; s < n_solution_vectors; ++s)
470 grad_dot_n_neighbor[s](0) = v;
471 }
472 else
473 {
475 neumann_bc.find(n)->second->vector_value(cell->vertex(n),
476 v);
477
478 for (unsigned int s = 0; s < n_solution_vectors; ++s)
479 grad_dot_n_neighbor[s] = v;
480 }
481 }
482 else
483 // fill with zeroes.
484 for (unsigned int s = 0; s < n_solution_vectors; ++s)
485 grad_dot_n_neighbor[s] = 0;
486
487 // if there is a coefficient, then evaluate it at the present
488 // position. if there is none, reuse the preset values.
489 if (coefficient != nullptr)
490 {
491 if (coefficient->n_components == 1)
492 {
493 const double c_value = coefficient->value(cell->vertex(n));
494 for (unsigned int c = 0; c < n_components; ++c)
495 coefficient_values(c) = c_value;
496 }
497 else
498 coefficient->vector_value(cell->vertex(n),
499 coefficient_values);
500 }
501
502
503 for (unsigned int s = 0; s < n_solution_vectors; ++s)
504 for (unsigned int component = 0; component < n_components;
505 ++component)
506 if (component_mask[component] == true)
507 {
508 // get gradient here
510 grad_dot_n_here =
511 gradients_here[s][n][component] * normal;
512
513 const typename ProductType<number, double>::type jump =
514 ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
515 coefficient_values(component));
516 (*errors[s])(cell->active_cell_index()) +=
518 typename ProductType<number,
519 double>::type>::abs_square(jump) *
520 cell->diameter();
521 }
522 }
523
524 for (unsigned int s = 0; s < n_solution_vectors; ++s)
525 (*errors[s])(cell->active_cell_index()) =
526 std::sqrt((*errors[s])(cell->active_cell_index()));
527 }
528}
529
530
531
532template <int spacedim>
533template <typename InputVector>
534void
536 const Mapping<1, spacedim> & mapping,
537 const DoFHandler<1, spacedim> &dof_handler,
538 const Quadrature<0> & quadrature,
539 const std::map<types::boundary_id,
541 & neumann_bc,
542 const std::vector<const InputVector *> &solutions,
543 std::vector<Vector<float> *> & errors,
544 const ComponentMask & component_mask,
545 const Function<spacedim> * coefficients,
546 const unsigned int n_threads,
547 const types::subdomain_id subdomain_id,
548 const types::material_id material_id,
549 const Strategy strategy)
550{
551 const hp::QCollection<0> quadrature_collection(quadrature);
552 estimate(mapping,
553 dof_handler,
554 quadrature_collection,
555 neumann_bc,
556 solutions,
557 errors,
558 component_mask,
559 coefficients,
560 n_threads,
561 subdomain_id,
562 material_id,
563 strategy);
564}
565
566
567
568// explicit instantiations
569#include "error_estimator_1d.inst"
570
571
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 hp::FECollection< dim, spacedim > & get_fe_collection() const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) 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:317
unsigned int n_active_cells() const
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
@ valid
Iterator points to a valid object.
constexpr const ReferenceCell Line
const types::material_id invalid_material_id
Definition types.h:250
const types::boundary_id internal_face_boundary_id
Definition types.h:277
const types::subdomain_id invalid_subdomain_id
Definition types.h:298
::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