Reference documentation for deal.II version GIT relicensing-218-g1f873f3929 2024-03-28 15:00: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\}}\)
Loading...
Searching...
No Matches
error_estimator_1d.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2013 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
19
21
24
25#include <deal.II/fe/fe.h>
28#include <deal.II/fe/mapping.h>
29
31
35
43#include <deal.II/lac/vector.h>
44
46
47#include <algorithm>
48#include <cmath>
49#include <functional>
50#include <numeric>
51#include <vector>
52
54
55
56template <int spacedim>
57template <typename Number>
58void
60 const Mapping<1, spacedim> &mapping,
61 const DoFHandler<1, spacedim> &dof_handler,
62 const Quadrature<0> &quadrature,
63 const std::map<types::boundary_id, const Function<spacedim, Number> *>
64 &neumann_bc,
65 const ReadVector<Number> &solution,
66 Vector<float> &error,
67 const ComponentMask &component_mask,
68 const Function<spacedim> *coefficients,
69 const unsigned int n_threads,
70 const types::subdomain_id subdomain_id,
71 const types::material_id material_id,
72 const Strategy strategy)
73{
74 // just pass on to the other function
75 std::vector<const ReadVector<Number> *> solutions(1, &solution);
76 std::vector<Vector<float> *> errors(1, &error);
77 ArrayView<Vector<float> *> error_view = make_array_view(errors);
78 estimate(mapping,
79 dof_handler,
80 quadrature,
81 neumann_bc,
82 make_array_view(solutions),
83 error_view,
84 component_mask,
85 coefficients,
86 n_threads,
87 subdomain_id,
88 material_id,
89 strategy);
90}
91
92
93
94template <int spacedim>
95template <typename Number>
96void
98 const DoFHandler<1, spacedim> &dof_handler,
99 const Quadrature<0> &quadrature,
100 const std::map<types::boundary_id, const Function<spacedim, Number> *>
101 &neumann_bc,
102 const ReadVector<Number> &solution,
103 Vector<float> &error,
104 const ComponentMask &component_mask,
105 const Function<spacedim> *coefficients,
106 const unsigned int n_threads,
107 const types::subdomain_id subdomain_id,
108 const types::material_id material_id,
109 const Strategy strategy)
110{
111 const auto reference_cell = ReferenceCells::Line;
112 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
113 dof_handler,
114 quadrature,
115 neumann_bc,
116 solution,
117 error,
118 component_mask,
119 coefficients,
120 n_threads,
121 subdomain_id,
122 material_id,
123 strategy);
124}
125
126
127
128template <int spacedim>
129template <typename Number>
130void
132 const DoFHandler<1, spacedim> &dof_handler,
133 const Quadrature<0> &quadrature,
134 const std::map<types::boundary_id, const Function<spacedim, Number> *>
135 &neumann_bc,
136 const ArrayView<const ReadVector<Number> *> &solutions,
137 ArrayView<Vector<float> *> &errors,
138 const ComponentMask &component_mask,
139 const Function<spacedim> *coefficients,
140 const unsigned int n_threads,
141 const types::subdomain_id subdomain_id,
142 const types::material_id material_id,
143 const Strategy strategy)
144{
145 const auto reference_cell = ReferenceCells::Line;
146 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
147 dof_handler,
148 quadrature,
149 neumann_bc,
150 solutions,
151 errors,
152 component_mask,
153 coefficients,
154 n_threads,
155 subdomain_id,
156 material_id,
157 strategy);
158}
159
160
161
162template <int spacedim>
163template <typename Number>
164void
166 const Mapping<1, spacedim> &mapping,
167 const DoFHandler<1, spacedim> &dof_handler,
168 const hp::QCollection<0> &quadrature,
169 const std::map<types::boundary_id, const Function<spacedim, Number> *>
170 &neumann_bc,
171 const ReadVector<Number> &solution,
172 Vector<float> &error,
173 const ComponentMask &component_mask,
174 const Function<spacedim> *coefficients,
175 const unsigned int n_threads,
176 const types::subdomain_id subdomain_id,
177 const types::material_id material_id,
178 const Strategy strategy)
179{
180 // just pass on to the other function
181 std::vector<const ReadVector<Number> *> solutions(1, &solution);
182 std::vector<Vector<float> *> errors(1, &error);
183 ArrayView<Vector<float> *> error_view = make_array_view(errors);
184 estimate(mapping,
185 dof_handler,
186 quadrature,
187 neumann_bc,
188 make_array_view(solutions),
189 error_view,
190 component_mask,
191 coefficients,
192 n_threads,
193 subdomain_id,
194 material_id,
195 strategy);
196}
197
198
199
200template <int spacedim>
201template <typename Number>
202void
204 const DoFHandler<1, spacedim> &dof_handler,
205 const hp::QCollection<0> &quadrature,
206 const std::map<types::boundary_id, const Function<spacedim, Number> *>
207 &neumann_bc,
208 const ReadVector<Number> &solution,
209 Vector<float> &error,
210 const ComponentMask &component_mask,
211 const Function<spacedim> *coefficients,
212 const unsigned int n_threads,
213 const types::subdomain_id subdomain_id,
214 const types::material_id material_id,
215 const Strategy strategy)
216{
217 const auto reference_cell = ReferenceCells::Line;
218 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
219 dof_handler,
220 quadrature,
221 neumann_bc,
222 solution,
223 error,
224 component_mask,
225 coefficients,
226 n_threads,
227 subdomain_id,
228 material_id,
229 strategy);
230}
231
232
233
234template <int spacedim>
235template <typename Number>
236void
238 const DoFHandler<1, spacedim> &dof_handler,
239 const hp::QCollection<0> &quadrature,
240 const std::map<types::boundary_id, const Function<spacedim, Number> *>
241 &neumann_bc,
242 const ArrayView<const ReadVector<Number> *> &solutions,
243 ArrayView<Vector<float> *> &errors,
244 const ComponentMask &component_mask,
245 const Function<spacedim> *coefficients,
246 const unsigned int n_threads,
247 const types::subdomain_id subdomain_id,
248 const types::material_id material_id,
249 const Strategy strategy)
250{
251 const auto reference_cell = ReferenceCells::Line;
252 estimate(reference_cell.template get_default_linear_mapping<1, spacedim>(),
253 dof_handler,
254 quadrature,
255 neumann_bc,
256 solutions,
257 errors,
258 component_mask,
259 coefficients,
260 n_threads,
261 subdomain_id,
262 material_id,
263 strategy);
264}
265
266
267
268template <int spacedim>
269template <typename Number>
270void
272 const Mapping<1, spacedim> &mapping,
273 const DoFHandler<1, spacedim> &dof_handler,
274 const hp::QCollection<0> &,
275 const std::map<types::boundary_id, const Function<spacedim, Number> *>
276 &neumann_bc,
277 const ArrayView<const ReadVector<Number> *> &solutions,
278 ArrayView<Vector<float> *> &errors,
279 const ComponentMask &component_mask,
280 const Function<spacedim> *coefficient,
281 const unsigned int,
282 const types::subdomain_id subdomain_id_,
283 const types::material_id material_id,
284 const Strategy strategy)
285{
286 AssertThrow(strategy == cell_diameter_over_24, ExcNotImplemented());
287 using number = Number;
289 if (const auto *triangulation = dynamic_cast<
291 &dof_handler.get_triangulation()))
292 {
293 Assert((subdomain_id_ == numbers::invalid_subdomain_id) ||
294 (subdomain_id_ == triangulation->locally_owned_subdomain()),
296 "For distributed Triangulation objects and associated "
297 "DoFHandler objects, asking for any subdomain other than the "
298 "locally owned one does not make sense."));
299 subdomain_id = triangulation->locally_owned_subdomain();
300 }
301 else
302 {
303 subdomain_id = subdomain_id_;
304 }
305
306 const unsigned int n_components = dof_handler.get_fe(0).n_components();
307 const unsigned int n_solution_vectors = solutions.size();
308
309 // sanity checks
311 neumann_bc.end(),
312 ExcMessage("You are not allowed to list the special boundary "
313 "indicator for internal boundaries in your boundary "
314 "value map."));
315
316 for (const auto &boundary_function : neumann_bc)
317 {
318 (void)boundary_function;
319 Assert(boundary_function.second->n_components == n_components,
320 ExcInvalidBoundaryFunction(boundary_function.first,
321 boundary_function.second->n_components,
322 n_components));
323 }
324
325 Assert(component_mask.represents_n_components(n_components),
326 ExcInvalidComponentMask());
327 Assert(component_mask.n_selected_components(n_components) > 0,
328 ExcInvalidComponentMask());
329
330 Assert((coefficient == nullptr) ||
331 (coefficient->n_components == n_components) ||
332 (coefficient->n_components == 1),
333 ExcInvalidCoefficient());
334
335 Assert(solutions.size() > 0, ExcNoSolutions());
336 Assert(solutions.size() == errors.size(),
337 ExcIncompatibleNumberOfElements(solutions.size(), errors.size()));
338 for (unsigned int n = 0; n < solutions.size(); ++n)
339 Assert(solutions[n]->size() == dof_handler.n_dofs(),
340 ExcDimensionMismatch(solutions[n]->size(), dof_handler.n_dofs()));
341
342 Assert((coefficient == nullptr) ||
343 (coefficient->n_components == n_components) ||
344 (coefficient->n_components == 1),
345 ExcInvalidCoefficient());
346
347 for (const auto &boundary_function : neumann_bc)
348 {
349 (void)boundary_function;
350 Assert(boundary_function.second->n_components == n_components,
351 ExcInvalidBoundaryFunction(boundary_function.first,
352 boundary_function.second->n_components,
353 n_components));
354 }
355
356 // reserve one slot for each cell and set it to zero
357 for (unsigned int n = 0; n < n_solution_vectors; ++n)
358 (*errors[n]).reinit(dof_handler.get_triangulation().n_active_cells());
359
360 // fields to get the gradients on the present and the neighbor cell.
361 //
362 // for the neighbor gradient, we need several auxiliary fields, depending on
363 // the way we get it (see below)
364 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
365 gradients_here(n_solution_vectors,
366 std::vector<std::vector<Tensor<1, spacedim, number>>>(
367 2,
368 std::vector<Tensor<1, spacedim, number>>(n_components)));
369 std::vector<std::vector<std::vector<Tensor<1, spacedim, number>>>>
370 gradients_neighbor(gradients_here);
371 std::vector<Vector<typename ProductType<number, double>::type>>
372 grad_dot_n_neighbor(n_solution_vectors,
374 n_components));
375
376 // reserve some space for coefficient values at one point. if there is no
377 // coefficient, then we fill it by unity once and for all and don't set it
378 // any more
379 Vector<double> coefficient_values(n_components);
380 if (coefficient == nullptr)
381 for (unsigned int c = 0; c < n_components; ++c)
382 coefficient_values(c) = 1;
383
384 const QTrapezoid<1> quadrature;
385 const hp::QCollection<1> q_collection(quadrature);
386 const QGauss<0> face_quadrature(1);
387 const hp::QCollection<0> q_face_collection(face_quadrature);
388
389 const hp::FECollection<1, spacedim> &fe = dof_handler.get_fe_collection();
390
391 hp::MappingCollection<1, spacedim> mapping_collection;
392 mapping_collection.push_back(mapping);
393
394 hp::FEValues<1, spacedim> fe_values(mapping_collection,
395 fe,
396 q_collection,
398 hp::FEFaceValues<1, spacedim> fe_face_values(
399 /*mapping_collection,*/ fe, q_face_collection, update_normal_vectors);
400
401 // loop over all cells and do something on the cells which we're told to
402 // work on. note that the error indicator is only a sum over the two
403 // contributions from the two vertices of each cell.
404 for (const auto &cell : dof_handler.active_cell_iterators())
405 if (((subdomain_id == numbers::invalid_subdomain_id) ||
406 (cell->subdomain_id() == subdomain_id)) &&
407 ((material_id == numbers::invalid_material_id) ||
408 (cell->material_id() == material_id)))
409 {
410 for (unsigned int n = 0; n < n_solution_vectors; ++n)
411 (*errors[n])(cell->active_cell_index()) = 0;
412
413 fe_values.reinit(cell);
414 for (unsigned int s = 0; s < n_solution_vectors; ++s)
415 fe_values.get_present_fe_values().get_function_gradients(
416 *solutions[s], gradients_here[s]);
417
418 // loop over the two points bounding this line. n==0 is left point,
419 // n==1 is right point
420 for (unsigned int n = 0; n < 2; ++n)
421 {
422 // find left or right active neighbor
423 auto neighbor = cell->neighbor(n);
424 if (neighbor.state() == IteratorState::valid)
425 while (neighbor->has_children())
426 neighbor = neighbor->child(n == 0 ? 1 : 0);
427
428 fe_face_values.reinit(cell, n);
429 Tensor<1, spacedim> normal =
430 fe_face_values.get_present_fe_values().get_normal_vectors()[0];
431
432 if (neighbor.state() == IteratorState::valid)
433 {
434 fe_values.reinit(neighbor);
435
436 for (unsigned int s = 0; s < n_solution_vectors; ++s)
437 fe_values.get_present_fe_values().get_function_gradients(
438 *solutions[s], gradients_neighbor[s]);
439
440 fe_face_values.reinit(neighbor, n == 0 ? 1 : 0);
441 Tensor<1, spacedim> neighbor_normal =
442 fe_face_values.get_present_fe_values()
443 .get_normal_vectors()[0];
444
445 // extract the gradient in normal direction of all the
446 // components.
447 for (unsigned int s = 0; s < n_solution_vectors; ++s)
448 for (unsigned int c = 0; c < n_components; ++c)
449 grad_dot_n_neighbor[s](c) =
450 -(gradients_neighbor[s][n == 0 ? 1 : 0][c] *
451 neighbor_normal);
452 }
453 else if (neumann_bc.find(n) != neumann_bc.end())
454 // if Neumann b.c., then fill the gradients field which will be
455 // used later on.
456 {
457 if (n_components == 1)
458 {
459 const Number v =
460 neumann_bc.find(n)->second->value(cell->vertex(n));
461
462 for (unsigned int s = 0; s < n_solution_vectors; ++s)
463 grad_dot_n_neighbor[s](0) = v;
464 }
465 else
466 {
467 Vector<Number> v(n_components);
468 neumann_bc.find(n)->second->vector_value(cell->vertex(n),
469 v);
470
471 for (unsigned int s = 0; s < n_solution_vectors; ++s)
472 grad_dot_n_neighbor[s] = v;
473 }
474 }
475 else
476 // fill with zeroes.
477 for (unsigned int s = 0; s < n_solution_vectors; ++s)
478 grad_dot_n_neighbor[s] = 0;
479
480 // if there is a coefficient, then evaluate it at the present
481 // position. if there is none, reuse the preset values.
482 if (coefficient != nullptr)
483 {
484 if (coefficient->n_components == 1)
485 {
486 const double c_value = coefficient->value(cell->vertex(n));
487 for (unsigned int c = 0; c < n_components; ++c)
488 coefficient_values(c) = c_value;
489 }
490 else
491 coefficient->vector_value(cell->vertex(n),
492 coefficient_values);
493 }
494
495
496 for (unsigned int s = 0; s < n_solution_vectors; ++s)
497 for (unsigned int component = 0; component < n_components;
498 ++component)
499 if (component_mask[component] == true)
500 {
501 // get gradient here
503 grad_dot_n_here =
504 gradients_here[s][n][component] * normal;
505
506 const typename ProductType<number, double>::type jump =
507 ((grad_dot_n_here - grad_dot_n_neighbor[s](component)) *
508 coefficient_values(component));
509 (*errors[s])(cell->active_cell_index()) +=
511 typename ProductType<number,
512 double>::type>::abs_square(jump) *
513 cell->diameter();
514 }
515 }
516
517 for (unsigned int s = 0; s < n_solution_vectors; ++s)
518 (*errors[s])(cell->active_cell_index()) =
519 std::sqrt((*errors[s])(cell->active_cell_index()));
520 }
521}
522
523
524
525template <int spacedim>
526template <typename Number>
527void
529 const Mapping<1, spacedim> &mapping,
530 const DoFHandler<1, spacedim> &dof_handler,
531 const Quadrature<0> &quadrature,
532 const std::map<types::boundary_id, const Function<spacedim, Number> *>
533 &neumann_bc,
534 const ArrayView<const ReadVector<Number> *> &solutions,
535 ArrayView<Vector<float> *> &errors,
536 const ComponentMask &component_mask,
537 const Function<spacedim> *coefficients,
538 const unsigned int n_threads,
539 const types::subdomain_id subdomain_id,
540 const types::material_id material_id,
541 const Strategy strategy)
542{
543 const hp::QCollection<0> quadrature_collection(quadrature);
544 estimate(mapping,
545 dof_handler,
546 quadrature_collection,
547 neumann_bc,
548 solutions,
549 errors,
550 component_mask,
551 coefficients,
552 n_threads,
553 subdomain_id,
554 material_id,
555 strategy);
556}
557
558
559
560// explicit instantiations
561#include "error_estimator_1d.inst"
562
563
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:949
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:163
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)
Abstract base class for mapping classes.
Definition mapping.h:316
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:423
const FEValuesType & get_present_fe_values() const
Definition fe_values.h:694
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:293
void push_back(const Mapping< dim, spacedim > &new_mapping)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:277
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id invalid_subdomain_id
Definition types.h:341
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
typename internal::ProductTypeImpl< std::decay_t< T >, std::decay_t< U > >::type type