Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16:50: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\}}\)
packaged_operation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 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 
16 #ifndef dealii_packaged_operation_h
17 #define dealii_packaged_operation_h
18 
19 #include <deal.II/base/config.h>
20 
22 
24 
25 #include <functional>
26 
28 
29 // Forward declarations:
30 #ifndef DOXYGEN
31 template <typename Number>
32 class Vector;
33 template <typename Range, typename Domain, typename Payload>
34 class LinearOperator;
35 template <typename Range = Vector<double>>
36 class PackagedOperation;
37 #endif
38 
39 
104 template <typename Range>
106 {
107 public:
114  {
115  apply = [](Range &) {
116  Assert(false,
117  ExcMessage(
118  "Uninitialized PackagedOperation<Range>::apply called"));
119  };
120 
121  apply_add = [](Range &) {
122  Assert(false,
123  ExcMessage(
124  "Uninitialized PackagedOperation<Range>::apply_add called"));
125  };
126 
127  reinit_vector = [](Range &, bool) {
128  Assert(false,
129  ExcMessage("Uninitialized PackagedOperation<Range>::reinit_vector "
130  "method called"));
131  };
132  }
133 
138 
148  PackagedOperation(const Range &u)
149  {
150  *this = u;
151  }
152 
158 
169  operator=(const Range &u)
170  {
171  apply = [&u](Range &v) { v = u; };
172 
173  apply_add = [&u](Range &v) { v += u; };
174 
175  reinit_vector = [&u](Range &v, bool omit_zeroing_entries) {
176  v.reinit(u, omit_zeroing_entries);
177  };
178 
179  return *this;
180  }
181 
188  operator Range() const
189  {
190  Range result_vector;
191 
192  reinit_vector(result_vector, /*bool omit_zeroing_entries=*/true);
193  apply(result_vector);
194 
195  return result_vector;
196  }
197 
208  {
209  *this = *this + second_comp;
210  return *this;
211  }
212 
219  {
220  *this = *this - second_comp;
221  return *this;
222  }
223 
229  operator+=(const Range &offset)
230  {
231  *this = *this + PackagedOperation<Range>(offset);
232  return *this;
233  }
234 
240  operator-=(const Range &offset)
241  {
242  *this = *this - PackagedOperation<Range>(offset);
243  return *this;
244  }
245 
250  operator*=(typename Range::value_type number)
251  {
252  *this = *this * number;
253  return *this;
254  }
261  std::function<void(Range &v)> apply;
262 
267  std::function<void(Range &v)> apply_add;
268 
276  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_vector;
277 };
278 
279 
293 template <typename Range>
296  const PackagedOperation<Range> &second_comp)
297 {
298  PackagedOperation<Range> return_comp;
299 
300  return_comp.reinit_vector = first_comp.reinit_vector;
301 
302  // ensure to have valid PackagedOperation objects by catching first_comp and
303  // second_comp by value
304 
305  return_comp.apply = [first_comp, second_comp](Range &v) {
306  first_comp.apply(v);
307  second_comp.apply_add(v);
308  };
309 
310  return_comp.apply_add = [first_comp, second_comp](Range &v) {
311  first_comp.apply_add(v);
312  second_comp.apply_add(v);
313  };
314 
315  return return_comp;
316 }
317 
326 template <typename Range>
329  const PackagedOperation<Range> &second_comp)
330 {
331  PackagedOperation<Range> return_comp;
332 
333  return_comp.reinit_vector = first_comp.reinit_vector;
334 
335  // ensure to have valid PackagedOperation objects by catching first_comp and
336  // second_comp by value
337 
338  return_comp.apply = [first_comp, second_comp](Range &v) {
339  second_comp.apply(v);
340  v *= -1.;
341  first_comp.apply_add(v);
342  };
343 
344  return_comp.apply_add = [first_comp, second_comp](Range &v) {
345  first_comp.apply_add(v);
346  v *= -1.;
347  second_comp.apply_add(v);
348  v *= -1.;
349  };
350 
351  return return_comp;
352 }
353 
362 template <typename Range>
365  typename Range::value_type number)
366 {
367  PackagedOperation<Range> return_comp;
368 
369  return_comp.reinit_vector = comp.reinit_vector;
370 
371  // the trivial case: number is zero
372  if (number == 0.)
373  {
374  return_comp.apply = [](Range &v) { v = 0.; };
375 
376  return_comp.apply_add = [](Range &) {};
377  }
378  else
379  {
380  return_comp.apply = [comp, number](Range &v) {
381  comp.apply(v);
382  v *= number;
383  };
384 
385  return_comp.apply_add = [comp, number](Range &v) {
386  v /= number;
387  comp.apply_add(v);
388  v *= number;
389  };
390  }
391 
392  return return_comp;
393 }
394 
403 template <typename Range>
405 operator*(typename Range::value_type number,
406  const PackagedOperation<Range> &comp)
407 {
408  return comp * number;
409 }
410 
419 template <typename Range>
421 operator+(const PackagedOperation<Range> &comp, const Range &offset)
422 {
423  return comp + PackagedOperation<Range>(offset);
424 }
425 
434 template <typename Range>
436 operator+(const Range &offset, const PackagedOperation<Range> &comp)
437 {
438  return PackagedOperation<Range>(offset) + comp;
439 }
440 
449 template <typename Range>
451 operator-(const PackagedOperation<Range> &comp, const Range &offset)
452 {
453  return comp - PackagedOperation<Range>(offset);
454 }
455 
456 
466 template <typename Range>
468 operator-(const Range &offset, const PackagedOperation<Range> &comp)
469 {
470  return PackagedOperation<Range>(offset) - comp;
471 }
472 
481 namespace internal
482 {
483  namespace PackagedOperationImplementation
484  {
485  // Poor man's trait class that determines whether type T is a vector:
486  // FIXME: Implement this as a proper type trait - similar to
487  // isBlockVector
488 
489  template <typename T>
491  {
492  template <typename C>
493  static std::false_type
494  test(...);
495 
496  template <typename C>
497  static std::true_type
498  test(decltype(&C::operator+=),
499  decltype(&C::operator-=),
500  decltype(&C::l2_norm));
501 
502  public:
503  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
504  // otherwise it is std::false_type
505 
506  using type = decltype(test<T>(nullptr, nullptr, nullptr));
507  }; // namespace
508  } // namespace PackagedOperationImplementation
509 } // namespace internal
510 
511 
526 template <
527  typename Range,
528  typename = std::enable_if_t<internal::PackagedOperationImplementation::
529  has_vector_interface<Range>::type::value>>
531 operator+(const Range &u, const Range &v)
532 {
533  PackagedOperation<Range> return_comp;
534 
535  // ensure to have valid PackagedOperation objects by catching op by value
536  // u is caught by reference
537 
538  return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) {
539  x.reinit(u, omit_zeroing_entries);
540  };
541 
542  return_comp.apply = [&u, &v](Range &x) {
543  x = u;
544  x += v;
545  };
546 
547  return_comp.apply_add = [&u, &v](Range &x) {
548  x += u;
549  x += v;
550  };
551 
552  return return_comp;
553 }
554 
555 
571 template <
572  typename Range,
573  typename = std::enable_if_t<internal::PackagedOperationImplementation::
574  has_vector_interface<Range>::type::value>>
576 operator-(const Range &u, const Range &v)
577 {
578  PackagedOperation<Range> return_comp;
579 
580  // ensure to have valid PackagedOperation objects by catching op by value
581  // u is caught by reference
582 
583  return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) {
584  x.reinit(u, omit_zeroing_entries);
585  };
586 
587  return_comp.apply = [&u, &v](Range &x) {
588  x = u;
589  x -= v;
590  };
591 
592  return_comp.apply_add = [&u, &v](Range &x) {
593  x += u;
594  x -= v;
595  };
596 
597  return return_comp;
598 }
599 
600 
615 template <
616  typename Range,
617  typename = std::enable_if_t<internal::PackagedOperationImplementation::
618  has_vector_interface<Range>::type::value>>
620 operator*(const Range &u, typename Range::value_type number)
621 {
622  return PackagedOperation<Range>(u) * number;
623 }
624 
625 
640 template <
641  typename Range,
642  typename = std::enable_if_t<internal::PackagedOperationImplementation::
643  has_vector_interface<Range>::type::value>>
645 operator*(typename Range::value_type number, const Range &u)
646 {
647  return number * PackagedOperation<Range>(u);
648 }
649 
650 
667 template <typename Range, typename Domain, typename Payload>
670 {
671  PackagedOperation<Range> return_comp;
672 
673  return_comp.reinit_vector = op.reinit_range_vector;
674 
675  // ensure to have valid PackagedOperation objects by catching op by value
676  // u is caught by reference
677 
678  return_comp.apply = [op, &u](Range &v) { op.vmult(v, u); };
679 
680  return_comp.apply_add = [op, &u](Range &v) { op.vmult_add(v, u); };
681 
682  return return_comp;
683 }
684 
685 
702 template <typename Range, typename Domain, typename Payload>
705 {
706  PackagedOperation<Range> return_comp;
707 
708  return_comp.reinit_vector = op.reinit_domain_vector;
709 
710  // ensure to have valid PackagedOperation objects by catching op by value
711  // u is caught by reference
712 
713  return_comp.apply = [op, &u](Domain &v) { op.Tvmult(v, u); };
714 
715  return_comp.apply_add = [op, &u](Domain &v) { op.Tvmult_add(v, u); };
716 
717  return return_comp;
718 }
719 
720 
729 template <typename Range, typename Domain, typename Payload>
732  const PackagedOperation<Domain> & comp)
733 {
734  PackagedOperation<Range> return_comp;
735 
736  return_comp.reinit_vector = op.reinit_range_vector;
737 
738  // ensure to have valid PackagedOperation objects by catching op by value
739  // u is caught by reference
740 
741  return_comp.apply = [op, comp](Domain &v) {
742  GrowingVectorMemory<Range> vector_memory;
743 
744  typename VectorMemory<Range>::Pointer i(vector_memory);
745  op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
746 
747  comp.apply(*i);
748  op.vmult(v, *i);
749  };
750 
751  return_comp.apply_add = [op, comp](Domain &v) {
752  GrowingVectorMemory<Range> vector_memory;
753 
754  typename VectorMemory<Range>::Pointer i(vector_memory);
755  op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
756 
757  comp.apply(*i);
758  op.vmult_add(v, *i);
759  };
760 
761  return return_comp;
762 }
763 
764 
773 template <typename Range, typename Domain, typename Payload>
777 {
778  PackagedOperation<Range> return_comp;
779 
780  return_comp.reinit_vector = op.reinit_domain_vector;
781 
782  // ensure to have valid PackagedOperation objects by catching op by value
783  // u is caught by reference
784 
785  return_comp.apply = [op, comp](Domain &v) {
786  GrowingVectorMemory<Range> vector_memory;
787 
788  typename VectorMemory<Range>::Pointer i(vector_memory);
789  op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
790 
791  comp.apply(*i);
792  op.Tvmult(v, *i);
793  };
794 
795  return_comp.apply_add = [op, comp](Domain &v) {
796  GrowingVectorMemory<Range> vector_memory;
797 
798  typename VectorMemory<Range>::Pointer i(vector_memory);
799  op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
800 
801  comp.apply(*i);
802  op.Tvmult_add(v, *i);
803  };
804 
805  return return_comp;
806 }
807 
811 
812 #endif
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
PackagedOperation< Range > & operator+=(const Range &offset)
PackagedOperation< Range > & operator=(const Range &u)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_vector
PackagedOperation(const Range &u)
std::function< void(Range &v)> apply
PackagedOperation< Range > & operator-=(const Range &offset)
std::function< void(Range &v)> apply_add
PackagedOperation(const PackagedOperation< Range > &)=default
PackagedOperation< Range > & operator-=(const PackagedOperation< Range > &second_comp)
PackagedOperation< Range > & operator=(const PackagedOperation< Range > &)=default
PackagedOperation< Range > & operator+=(const PackagedOperation< Range > &second_comp)
PackagedOperation< Range > & operator*=(typename Range::value_type number)
Definition: vector.h:109
static std::true_type test(decltype(&C::operator+=), decltype(&C::operator-=), decltype(&C::l2_norm))
decltype(test< T >(nullptr, nullptr, nullptr)) type
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcMessage(std::string arg1)
PackagedOperation< Range > operator*(const LinearOperator< Range, Domain, Payload > &op, const PackagedOperation< Domain > &comp)
PackagedOperation< Range > operator-(const PackagedOperation< Range > &comp, const Range &offset)
PackagedOperation< Domain > operator*(const PackagedOperation< Range > &comp, const LinearOperator< Range, Domain, Payload > &op)
PackagedOperation< Range > operator-(const Range &u, const Range &v)
PackagedOperation< Range > operator*(const LinearOperator< Range, Domain, Payload > &op, const Domain &u)
PackagedOperation< Domain > operator*(const Range &u, const LinearOperator< Range, Domain, Payload > &op)
PackagedOperation< Range > operator+(const Range &offset, const PackagedOperation< Range > &comp)
PackagedOperation< Range > operator*(typename Range::value_type number, const Range &u)
PackagedOperation< Range > operator-(const Range &offset, const PackagedOperation< Range > &comp)
PackagedOperation< Range > operator-(const PackagedOperation< Range > &first_comp, const PackagedOperation< Range > &second_comp)
PackagedOperation< Range > operator+(const PackagedOperation< Range > &comp, const Range &offset)
PackagedOperation< Range > operator+(const Range &u, const Range &v)
PackagedOperation< Range > operator*(const PackagedOperation< Range > &comp, typename Range::value_type number)
PackagedOperation< Range > operator*(typename Range::value_type number, const PackagedOperation< Range > &comp)
PackagedOperation< Range > operator*(const Range &u, typename Range::value_type number)
PackagedOperation< Range > operator+(const PackagedOperation< Range > &first_comp, const PackagedOperation< Range > &second_comp)
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)