Reference documentation for deal.II version 9.4.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\}}\)
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
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>
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
403template <typename Range>
405operator*(typename Range::value_type number,
406 const PackagedOperation<Range> &comp)
407{
408 return comp * number;
409}
410
419template <typename Range>
421operator+(const PackagedOperation<Range> &comp, const Range &offset)
422{
423 return comp + PackagedOperation<Range>(offset);
424}
425
434template <typename Range>
436operator+(const Range &offset, const PackagedOperation<Range> &comp)
437{
438 return PackagedOperation<Range>(offset) + comp;
439}
440
449template <typename Range>
451operator-(const PackagedOperation<Range> &comp, const Range &offset)
452{
453 return comp - PackagedOperation<Range>(offset);
454}
455
456
466template <typename Range>
468operator-(const Range &offset, const PackagedOperation<Range> &comp)
469{
470 return PackagedOperation<Range>(offset) - comp;
471}
472
474
475
480
481namespace 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
526template <typename Range,
527 typename = typename std::enable_if<
529 Range>::type::value>::type>
531operator+(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
571template <typename Range,
572 typename = typename std::enable_if<
574 Range>::type::value>::type>
576operator-(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
615template <typename Range,
616 typename = typename std::enable_if<
618 Range>::type::value>::type>
620operator*(const Range &u, typename Range::value_type number)
621{
622 return PackagedOperation<Range>(u) * number;
623}
624
625
640template <typename Range,
641 typename = typename std::enable_if<
643 Range>::type::value>::type>
645operator*(typename Range::value_type number, const Range &u)
646{
647 return number * PackagedOperation<Range>(u);
648}
649
650
667template <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
702template <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
729template <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
773template <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
809
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
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: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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
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)