Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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
packaged_operation.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2014 - 2020 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
31template <typename Number>
32class Vector;
33template <typename Range, typename Domain, typename Payload>
34class LinearOperator;
35template <typename Range = Vector<double>>
37#endif
38
39
104template <typename Range>
106{
107public:
114 {
115 apply = [](Range &) {
116 Assert(false,
118 "Uninitialized PackagedOperation<Range>::apply called"));
119 };
120
121 apply_add = [](Range &) {
122 Assert(false,
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
202
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 }
256
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
284
293template <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
326template <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
362template <typename Range>
364 typename Range::value_type number)
365{
366 PackagedOperation<Range> return_comp;
367
368 return_comp.reinit_vector = comp.reinit_vector;
369
370 // the trivial case: number is zero
371 if (number == 0.)
372 {
373 return_comp.apply = [](Range &v) { v = 0.; };
374
375 return_comp.apply_add = [](Range &) {};
376 }
377 else
378 {
379 return_comp.apply = [comp, number](Range &v) {
380 comp.apply(v);
381 v *= number;
382 };
383
384 return_comp.apply_add = [comp, number](Range &v) {
385 v /= number;
386 comp.apply_add(v);
387 v *= number;
388 };
389 }
390
391 return return_comp;
392}
393
402template <typename Range>
403PackagedOperation<Range> operator*(typename Range::value_type number,
404 const PackagedOperation<Range> &comp)
405{
406 return comp * number;
407}
408
417template <typename Range>
419operator+(const PackagedOperation<Range> &comp, const Range &offset)
420{
421 return comp + PackagedOperation<Range>(offset);
422}
423
432template <typename Range>
434operator+(const Range &offset, const PackagedOperation<Range> &comp)
435{
436 return PackagedOperation<Range>(offset) + comp;
437}
438
447template <typename Range>
449operator-(const PackagedOperation<Range> &comp, const Range &offset)
450{
451 return comp - PackagedOperation<Range>(offset);
452}
453
454
464template <typename Range>
466operator-(const Range &offset, const PackagedOperation<Range> &comp)
467{
468 return PackagedOperation<Range>(offset) - comp;
469}
470
472
473
478
479namespace internal
480{
481 namespace PackagedOperationImplementation
482 {
483 // Poor man's trait class that determines whether type T is a vector:
484 // FIXME: Implement this as a proper type trait - similar to
485 // isBlockVector
486
487 template <typename T>
489 {
490 template <typename C>
491 static std::false_type
492 test(...);
493
494 template <typename C>
495 static std::true_type
496 test(decltype(&C::operator+=),
497 decltype(&C::operator-=),
498 decltype(&C::l2_norm));
499
500 public:
501 // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
502 // otherwise it is std::false_type
503
504 using type = decltype(test<T>(nullptr, nullptr, nullptr));
505 }; // namespace
506 } // namespace PackagedOperationImplementation
507} // namespace internal
508
509
524template <typename Range,
525 typename = typename std::enable_if<
527 Range>::type::value>::type>
529operator+(const Range &u, const Range &v)
530{
531 PackagedOperation<Range> return_comp;
532
533 // ensure to have valid PackagedOperation objects by catching op by value
534 // u is caught by reference
535
536 return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) {
537 x.reinit(u, omit_zeroing_entries);
538 };
539
540 return_comp.apply = [&u, &v](Range &x) {
541 x = u;
542 x += v;
543 };
544
545 return_comp.apply_add = [&u, &v](Range &x) {
546 x += u;
547 x += v;
548 };
549
550 return return_comp;
551}
552
553
569template <typename Range,
570 typename = typename std::enable_if<
572 Range>::type::value>::type>
574operator-(const Range &u, const Range &v)
575{
576 PackagedOperation<Range> return_comp;
577
578 // ensure to have valid PackagedOperation objects by catching op by value
579 // u is caught by reference
580
581 return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries) {
582 x.reinit(u, omit_zeroing_entries);
583 };
584
585 return_comp.apply = [&u, &v](Range &x) {
586 x = u;
587 x -= v;
588 };
589
590 return_comp.apply_add = [&u, &v](Range &x) {
591 x += u;
592 x -= v;
593 };
594
595 return return_comp;
596}
597
598
613template <typename Range,
614 typename = typename std::enable_if<
616 Range>::type::value>::type>
618 typename Range::value_type number)
619{
620 return PackagedOperation<Range>(u) * number;
621}
622
623
638template <typename Range,
639 typename = typename std::enable_if<
641 Range>::type::value>::type>
642PackagedOperation<Range> operator*(typename Range::value_type number,
643 const Range & u)
644{
645 return number * PackagedOperation<Range>(u);
646}
647
648
665template <typename Range, typename Domain, typename Payload>
668{
669 PackagedOperation<Range> return_comp;
670
671 return_comp.reinit_vector = op.reinit_range_vector;
672
673 // ensure to have valid PackagedOperation objects by catching op by value
674 // u is caught by reference
675
676 return_comp.apply = [op, &u](Range &v) { op.vmult(v, u); };
677
678 return_comp.apply_add = [op, &u](Range &v) { op.vmult_add(v, u); };
679
680 return return_comp;
681}
682
683
700template <typename Range, typename Domain, typename Payload>
703{
704 PackagedOperation<Range> return_comp;
705
706 return_comp.reinit_vector = op.reinit_domain_vector;
707
708 // ensure to have valid PackagedOperation objects by catching op by value
709 // u is caught by reference
710
711 return_comp.apply = [op, &u](Domain &v) { op.Tvmult(v, u); };
712
713 return_comp.apply_add = [op, &u](Domain &v) { op.Tvmult_add(v, u); };
714
715 return return_comp;
716}
717
718
727template <typename Range, typename Domain, typename Payload>
730 const PackagedOperation<Domain> & comp)
731{
732 PackagedOperation<Range> return_comp;
733
734 return_comp.reinit_vector = op.reinit_range_vector;
735
736 // ensure to have valid PackagedOperation objects by catching op by value
737 // u is caught by reference
738
739 return_comp.apply = [op, comp](Domain &v) {
740 GrowingVectorMemory<Range> vector_memory;
741
742 typename VectorMemory<Range>::Pointer i(vector_memory);
743 op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
744
745 comp.apply(*i);
746 op.vmult(v, *i);
747 };
748
749 return_comp.apply_add = [op, comp](Domain &v) {
750 GrowingVectorMemory<Range> vector_memory;
751
752 typename VectorMemory<Range>::Pointer i(vector_memory);
753 op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
754
755 comp.apply(*i);
756 op.vmult_add(v, *i);
757 };
758
759 return return_comp;
760}
761
762
771template <typename Range, typename Domain, typename Payload>
775{
776 PackagedOperation<Range> return_comp;
777
778 return_comp.reinit_vector = op.reinit_domain_vector;
779
780 // ensure to have valid PackagedOperation objects by catching op by value
781 // u is caught by reference
782
783 return_comp.apply = [op, comp](Domain &v) {
784 GrowingVectorMemory<Range> vector_memory;
785
786 typename VectorMemory<Range>::Pointer i(vector_memory);
787 op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
788
789 comp.apply(*i);
790 op.Tvmult(v, *i);
791 };
792
793 return_comp.apply_add = [op, comp](Domain &v) {
794 GrowingVectorMemory<Range> vector_memory;
795
796 typename VectorMemory<Range>::Pointer i(vector_memory);
797 op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
798
799 comp.apply(*i);
800 op.Tvmult_add(v, *i);
801 };
802
803 return return_comp;
804}
805
807
809
810#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
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_vector
PackagedOperation< Range > & operator*=(typename Range::value_type number)
PackagedOperation< Range > & operator=(const Range &u)
PackagedOperation(const Range &u)
std::function< void(Range &v)> apply
PackagedOperation< Range > & operator=(const PackagedOperation< Range > &)=default
PackagedOperation< Range > & operator+=(const PackagedOperation< Range > &second_comp)
std::function< void(Range &v)> apply_add
PackagedOperation< Range > & operator-=(const Range &offset)
PackagedOperation(const PackagedOperation< Range > &)=default
PackagedOperation< Range > & operator+=(const Range &offset)
PackagedOperation< Range > & operator-=(const PackagedOperation< Range > &second_comp)
Definition: vector.h:110
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:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
#define Assert(cond, exc)
Definition: exceptions.h:1465
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)