Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
packaged_operation.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2019 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 
21 #include <deal.II/base/exceptions.h>
22 
23 #include <deal.II/lac/vector_memory.h>
24 
25 #include <functional>
26 
27 DEAL_II_NAMESPACE_OPEN
28 
29 // Forward declarations:
30 template <typename Number>
31 class Vector;
32 template <typename Range, typename Domain, typename Payload>
33 class LinearOperator;
34 template <typename Range = Vector<double>>
36 
37 
103 template <typename Range>
104 class PackagedOperation
105 {
106 public:
113  {
114  apply = [](Range &) {
115  Assert(false,
116  ExcMessage(
117  "Uninitialized PackagedOperation<Range>::apply called"));
118  };
119 
120  apply_add = [](Range &) {
121  Assert(false,
122  ExcMessage(
123  "Uninitialized PackagedOperation<Range>::apply_add called"));
124  };
125 
126  reinit_vector = [](Range &, bool) {
127  Assert(false,
128  ExcMessage("Uninitialized PackagedOperation<Range>::reinit_vector "
129  "method called"));
130  };
131  }
132 
136  PackagedOperation(const PackagedOperation<Range> &) = default;
137 
147  PackagedOperation(const Range &u)
148  {
149  *this = u;
150  }
151 
156  operator=(const PackagedOperation<Range> &) = default;
157 
168  operator=(const Range &u)
169  {
170  apply = [&u](Range &v) { v = u; };
171 
172  apply_add = [&u](Range &v) { v += u; };
173 
174  reinit_vector = [&u](Range &v, bool omit_zeroing_entries) {
175  v.reinit(u, omit_zeroing_entries);
176  };
177 
178  return *this;
179  }
180 
187  operator Range() const
188  {
189  Range result_vector;
190 
191  reinit_vector(result_vector, /*bool omit_zeroing_entries=*/true);
192  apply(result_vector);
193 
194  return result_vector;
195  }
196 
201 
207  {
208  *this = *this + second_comp;
209  return *this;
210  }
211 
218  {
219  *this = *this - second_comp;
220  return *this;
221  }
222 
228  operator+=(const Range &offset)
229  {
230  *this = *this + PackagedOperation<Range>(offset);
231  return *this;
232  }
233 
239  operator-=(const Range &offset)
240  {
241  *this = *this - PackagedOperation<Range>(offset);
242  return *this;
243  }
244 
249  operator*=(typename Range::value_type number)
250  {
251  *this = *this * number;
252  return *this;
253  }
255 
260  std::function<void(Range &v)> apply;
261 
266  std::function<void(Range &v)> apply_add;
267 
275  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_vector;
276 };
277 
278 
283 
292 template <typename Range>
295  const PackagedOperation<Range> &second_comp)
296 {
297  PackagedOperation<Range> return_comp;
298 
299  return_comp.reinit_vector = first_comp.reinit_vector;
300 
301  // ensure to have valid PackagedOperation objects by catching first_comp and
302  // second_comp by value
303 
304  return_comp.apply = [first_comp, second_comp](Range &v) {
305  first_comp.apply(v);
306  second_comp.apply_add(v);
307  };
308 
309  return_comp.apply_add = [first_comp, second_comp](Range &v) {
310  first_comp.apply_add(v);
311  second_comp.apply_add(v);
312  };
313 
314  return return_comp;
315 }
316 
325 template <typename Range>
328  const PackagedOperation<Range> &second_comp)
329 {
330  PackagedOperation<Range> return_comp;
331 
332  return_comp.reinit_vector = first_comp.reinit_vector;
333 
334  // ensure to have valid PackagedOperation objects by catching first_comp and
335  // second_comp by value
336 
337  return_comp.apply = [first_comp, second_comp](Range &v) {
338  second_comp.apply(v);
339  v *= -1.;
340  first_comp.apply_add(v);
341  };
342 
343  return_comp.apply_add = [first_comp, second_comp](Range &v) {
344  first_comp.apply_add(v);
345  v *= -1.;
346  second_comp.apply_add(v);
347  v *= -1.;
348  };
349 
350  return return_comp;
351 }
352 
361 template <typename Range>
363  typename Range::value_type number)
364 {
365  PackagedOperation<Range> return_comp;
366 
367  return_comp.reinit_vector = comp.reinit_vector;
368 
369  // the trivial case: number is zero
370  if (number == 0.)
371  {
372  return_comp.apply = [](Range &v) { v = 0.; };
373 
374  return_comp.apply_add = [](Range &) {};
375  }
376  else
377  {
378  return_comp.apply = [comp, number](Range &v) {
379  comp.apply(v);
380  v *= number;
381  };
382 
383  return_comp.apply_add = [comp, number](Range &v) {
384  v /= number;
385  comp.apply_add(v);
386  v *= number;
387  };
388  }
389 
390  return return_comp;
391 }
392 
401 template <typename Range>
402 PackagedOperation<Range> operator*(typename Range::value_type number,
403  const PackagedOperation<Range> &comp)
404 {
405  return comp * number;
406 }
407 
416 template <typename Range>
418 operator+(const PackagedOperation<Range> &comp, const Range &offset)
419 {
420  return comp + PackagedOperation<Range>(offset);
421 }
422 
431 template <typename Range>
433 operator+(const Range &offset, const PackagedOperation<Range> &comp)
434 {
435  return PackagedOperation<Range>(offset) + comp;
436 }
437 
446 template <typename Range>
448 operator-(const PackagedOperation<Range> &comp, const Range &offset)
449 {
450  return comp - PackagedOperation<Range>(offset);
451 }
452 
453 
463 template <typename Range>
465 operator-(const Range &offset, const PackagedOperation<Range> &comp)
466 {
467  return PackagedOperation<Range>(offset) - comp;
468 }
469 
471 
472 
477 
478 namespace internal
479 {
480  namespace PackagedOperationImplementation
481  {
482  // Poor man's trait class that determines whether type T is a vector:
483  // FIXME: Implement this as a proper type trait - similar to
484  // isBlockVector
485 
486  template <typename T>
487  class has_vector_interface
488  {
489  template <typename C>
490  static std::false_type
491  test(...);
492 
493  template <typename C>
494  static std::true_type
495  test(decltype(&C::operator+=),
496  decltype(&C::operator-=),
497  decltype(&C::l2_norm));
498 
499  public:
500  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
501  // otherwise it is std::false_type
502 
503  using type = decltype(test<T>(nullptr, nullptr, nullptr));
504  }; // namespace
505  } // namespace PackagedOperationImplementation
506 } // namespace internal
507 
508 
523 template <typename Range,
524  typename = typename std::enable_if<
525  internal::PackagedOperationImplementation::has_vector_interface<
526  Range>::type::value>::type>
528 operator+(const Range &u, const Range &v)
529 {
530  PackagedOperation<Range> return_comp;
531 
532  // ensure to have valid PackagedOperation objects by catching op by value
533  // u is caught by reference
534 
535  return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) {
536  x.reinit(u, omit_zeroing_entries);
537  };
538 
539  return_comp.apply = [&u, &v](Range &x) {
540  x = u;
541  x += v;
542  };
543 
544  return_comp.apply_add = [&u, &v](Range &x) {
545  x += u;
546  x += v;
547  };
548 
549  return return_comp;
550 }
551 
552 
568 template <typename Range,
569  typename = typename std::enable_if<
570  internal::PackagedOperationImplementation::has_vector_interface<
571  Range>::type::value>::type>
573 operator-(const Range &u, const Range &v)
574 {
575  PackagedOperation<Range> return_comp;
576 
577  // ensure to have valid PackagedOperation objects by catching op by value
578  // u is caught by reference
579 
580  return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) {
581  x.reinit(u, omit_zeroing_entries);
582  };
583 
584  return_comp.apply = [&u, &v](Range &x) {
585  x = u;
586  x -= v;
587  };
588 
589  return_comp.apply_add = [&u, &v](Range &x) {
590  x += u;
591  x -= v;
592  };
593 
594  return return_comp;
595 }
596 
597 
612 template <typename Range,
613  typename = typename std::enable_if<
614  internal::PackagedOperationImplementation::has_vector_interface<
615  Range>::type::value>::type>
617  typename Range::value_type number)
618 {
619  return PackagedOperation<Range>(u) * number;
620 }
621 
622 
637 template <typename Range,
638  typename = typename std::enable_if<
639  internal::PackagedOperationImplementation::has_vector_interface<
640  Range>::type::value>::type>
641 PackagedOperation<Range> operator*(typename Range::value_type number,
642  const Range & u)
643 {
644  return number * PackagedOperation<Range>(u);
645 }
646 
647 
664 template <typename Range, typename Domain, typename Payload>
667 {
668  PackagedOperation<Range> return_comp;
669 
670  return_comp.reinit_vector = op.reinit_range_vector;
671 
672  // ensure to have valid PackagedOperation objects by catching op by value
673  // u is caught by reference
674 
675  return_comp.apply = [op, &u](Range &v) { op.vmult(v, u); };
676 
677  return_comp.apply_add = [op, &u](Range &v) { op.vmult_add(v, u); };
678 
679  return return_comp;
680 }
681 
682 
699 template <typename Range, typename Domain, typename Payload>
702 {
703  PackagedOperation<Range> return_comp;
704 
705  return_comp.reinit_vector = op.reinit_domain_vector;
706 
707  // ensure to have valid PackagedOperation objects by catching op by value
708  // u is caught by reference
709 
710  return_comp.apply = [op, &u](Domain &v) { op.Tvmult(v, u); };
711 
712  return_comp.apply_add = [op, &u](Domain &v) { op.Tvmult_add(v, u); };
713 
714  return return_comp;
715 }
716 
717 
726 template <typename Range, typename Domain, typename Payload>
729  const PackagedOperation<Domain> & comp)
730 {
731  PackagedOperation<Range> return_comp;
732 
733  return_comp.reinit_vector = op.reinit_range_vector;
734 
735  // ensure to have valid PackagedOperation objects by catching op by value
736  // u is caught by reference
737 
738  return_comp.apply = [op, comp](Domain &v) {
739  GrowingVectorMemory<Range> vector_memory;
740 
741  typename VectorMemory<Range>::Pointer i(vector_memory);
742  op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
743 
744  comp.apply(*i);
745  op.vmult(v, *i);
746  };
747 
748  return_comp.apply_add = [op, comp](Domain &v) {
749  GrowingVectorMemory<Range> vector_memory;
750 
751  typename VectorMemory<Range>::Pointer i(vector_memory);
752  op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
753 
754  comp.apply(*i);
755  op.vmult_add(v, *i);
756  };
757 
758  return return_comp;
759 }
760 
761 
770 template <typename Range, typename Domain, typename Payload>
774 {
775  PackagedOperation<Range> return_comp;
776 
777  return_comp.reinit_vector = op.reinit_domain_vector;
778 
779  // ensure to have valid PackagedOperation objects by catching op by value
780  // u is caught by reference
781 
782  return_comp.apply = [op, comp](Domain &v) {
783  GrowingVectorMemory<Range> vector_memory;
784 
785  typename VectorMemory<Range>::Pointer i(vector_memory);
786  op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
787 
788  comp.apply(*i);
789  op.Tvmult(v, *i);
790  };
791 
792  return_comp.apply_add = [op, comp](Domain &v) {
793  GrowingVectorMemory<Range> vector_memory;
794 
795  typename VectorMemory<Range>::Pointer i(vector_memory);
796  op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
797 
798  comp.apply(*i);
799  op.Tvmult_add(v, *i);
800  };
801 
802  return return_comp;
803 }
804 
806 
807 DEAL_II_NAMESPACE_CLOSE
808 
809 #endif
PackagedOperation< Range > & operator+=(const PackagedOperation< Range > &second_comp)
PackagedOperation(const Range &u)
std::function< void(Range &v)> apply
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_vector
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
PackagedOperation< Range > & operator-=(const PackagedOperation< Range > &second_comp)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
LinearAlgebra::distributed::Vector< Number > Vector
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
PackagedOperation< Range > & operator=(const PackagedOperation< Range > &)=default
static ::ExceptionBase & ExcMessage(std::string arg1)
PackagedOperation< Range > & operator=(const Range &u)
#define Assert(cond, exc)
Definition: exceptions.h:1407
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
PackagedOperation< Range > & operator*=(typename Range::value_type number)
std::function< void(Range &v)> apply_add
std::function< void(Range &v, const Domain &u)> vmult_add
PackagedOperation< Range > & operator+=(const Range &offset)
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
PackagedOperation< Range > & operator-=(const Range &offset)