]> https://gitweb.dealii.org/ - dealii.git/blob - include/deal.II/meshworker/simple.h
Update copyright years
[dealii.git] / include / deal.II / meshworker / simple.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 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
17 #ifndef dealii_mesh_worker_simple_h
18 #define dealii_mesh_worker_simple_h
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/algorithms/any_data.h>
23
24 #include <deal.II/base/mg_level_object.h>
25 #include <deal.II/base/smartpointer.h>
26
27 #include <deal.II/lac/block_vector.h>
28
29 #include <deal.II/meshworker/dof_info.h>
30 #include <deal.II/meshworker/functional.h>
31
32 #include <deal.II/multigrid/mg_constrained_dofs.h>
33
34 /*
35  * The header containing the classes MeshWorker::Assembler::MatrixSimple,
36  * MeshWorker::Assembler::MGMatrixSimple, MeshWorker::Assembler::ResidualSimple,
37  * and MeshWorker::Assembler::SystemSimple.
38  */
39
40 DEAL_II_NAMESPACE_OPEN
41
42 namespace MeshWorker
43 {
44   namespace Assembler
45   {
46     /**
47      * Assemble residuals without block structure.
48      *
49      * The data structure for this Assembler class is a simple vector on each
50      * cell with entries from zero to FiniteElementData::dofs_per_cell and a
51      * simple global vector with entries numbered from zero to
52      * DoFHandler::n_dofs(). No BlockInfo is required and the global vector
53      * may be any type of vector having element access through <tt>operator()
54      * (unsigned int)</tt>
55      *
56      * @ingroup MeshWorker
57      */
58     template <typename VectorType>
59     class ResidualSimple
60     {
61     public:
62       /**
63        * Initialize with an AnyData object holding the result of assembling.
64        *
65        * Assembling currently writes into the first vector of
66        * <tt>results</tt>.
67        */
68       void
69       initialize(AnyData &results);
70
71       /**
72        * Initialize the constraints.
73        */
74       void
75       initialize(
76         const AffineConstraints<typename VectorType::value_type> &constraints);
77
78       /**
79        * Initialize the local data in the DoFInfo object used later for
80        * assembling.
81        *
82        * The @p info object refers to a cell if <code>!face</code>, or else to an
83        * interior or boundary face.
84        */
85       template <class DOFINFO>
86       void
87       initialize_info(DOFINFO &info, bool face) const;
88
89       /**
90        * Assemble the local residuals into the global residuals.
91        *
92        * Values are added to the previous contents. If constraints are active,
93        * AffineConstraints::distribute_local_to_global() is used.
94        */
95       template <class DOFINFO>
96       void
97       assemble(const DOFINFO &info);
98
99       /**
100        * Assemble both local residuals into the global residuals.
101        */
102       template <class DOFINFO>
103       void
104       assemble(const DOFINFO &info1, const DOFINFO &info2);
105
106     protected:
107       /**
108        * The global residual vectors filled by assemble().
109        */
110       AnyData residuals;
111
112       /**
113        * A pointer to the object containing constraints.
114        */
115       SmartPointer<const AffineConstraints<typename VectorType::value_type>,
116                    ResidualSimple<VectorType>>
117         constraints;
118     };
119
120
121     /**
122      * Assemble local matrices into a single global matrix or several global
123      * matrices associated with the same DoFHandler. If these global matrix
124      * have a block structure, this structure is not used, but rather the
125      * global numbering of degrees of freedom.
126      *
127      * After being initialized with a SparseMatrix object (or another matrix
128      * offering the same functionality as SparseMatrix::add()) or a vector of
129      * such, this class can be used in a MeshWorker::loop() to assemble the
130      * cell and face matrices into the global matrix.
131      *
132      * If a AffineConstraints has been provided during initialization, this
133      * matrix will be used (AffineConstraints::distribute_local_to_global(), to
134      * be precise) to enter the local matrix into the global sparse matrix.
135      *
136      * The assembler can handle two different types of local data. First, by
137      * default, the obvious choice of taking a single local matrix with
138      * dimensions equal to the number of degrees of freedom of the cell.
139      * Alternatively, a local block structure can be initialized in DoFInfo.
140      * After this, the local data will be arranged as an array of n by n
141      * FullMatrix blocks (n being the number of blocks in the FESystem used by
142      * the DoFHandler in DoFInfo), which are ordered lexicographically with
143      * column index fastest in DoFInfo. If the matrix was initialized with a
144      * vector of several matrices and local block structure is used, then the
145      * first n<sup>2</sup> matrices in LocalResults will be used for the first
146      * matrix in this vector, the second set of n<sup>2</sup> for the second,
147      * and so on.
148      *
149      * @ingroup MeshWorker
150      */
151     template <typename MatrixType>
152     class MatrixSimple
153     {
154     public:
155       /**
156        * Constructor, initializing the #threshold, which limits how small
157        * numbers may be to be entered into the matrix.
158        */
159       MatrixSimple(double threshold = 1.e-12);
160
161       /**
162        * Store the result matrix for later assembling.
163        */
164       void
165       initialize(MatrixType &m);
166
167       /**
168        * Store several result matrices for later assembling.
169        */
170       void
171       initialize(std::vector<MatrixType> &m);
172
173       /**
174        * Initialize the constraints. After this function has been called with
175        * a valid AffineConstraints object, the function
176        * AffineConstraints::distribute_local_to_global() will be used by
177        * assemble() to distribute the cell and face matrices into a global
178        * sparse matrix.
179        */
180       void
181       initialize(
182         const AffineConstraints<typename MatrixType::value_type> &constraints);
183
184       /**
185        * Initialize the local data in the DoFInfo object used later for
186        * assembling.
187        *
188        * The @p info object refers to a cell if <code>!face</code>, or else to an
189        * interior or boundary face.
190        */
191       template <class DOFINFO>
192       void
193       initialize_info(DOFINFO &info, bool face) const;
194
195       /**
196        * Assemble the local matrices associated with a single cell into the
197        * global matrix.
198        */
199       template <class DOFINFO>
200       void
201       assemble(const DOFINFO &info);
202
203       /**
204        * Assemble all local matrices associated with an interior face in the
205        * @p info1 and @p info2 objects into the global matrix.
206        */
207       template <class DOFINFO>
208       void
209       assemble(const DOFINFO &info1, const DOFINFO &info2);
210
211     protected:
212       /**
213        * The vector of global matrices being assembled.
214        */
215       std::vector<SmartPointer<MatrixType, MatrixSimple<MatrixType>>> matrix;
216
217       /**
218        * The smallest positive number that will be entered into the global
219        * matrix. All smaller absolute values will be treated as zero and will
220        * not be assembled.
221        */
222       const double threshold;
223
224     private:
225       /**
226        * Assemble a single matrix <code>M</code> into the element at
227        * <code>index</code> in the vector #matrix.
228        */
229       void
230       assemble(const FullMatrix<double> &                  M,
231                const unsigned int                          index,
232                const std::vector<types::global_dof_index> &i1,
233                const std::vector<types::global_dof_index> &i2);
234
235       /**
236        * A pointer to the object containing constraints.
237        */
238       SmartPointer<const AffineConstraints<typename MatrixType::value_type>,
239                    MatrixSimple<MatrixType>>
240         constraints;
241     };
242
243
244     /**
245      * Assemble local matrices into level matrices without using block
246      * structure.
247      *
248      * @todo The matrix structures needed for assembling level matrices with
249      * local refinement and continuous elements are missing.
250      *
251      * @ingroup MeshWorker
252      */
253     template <typename MatrixType>
254     class MGMatrixSimple
255     {
256     public:
257       /**
258        * Constructor, initializing the #threshold, which limits how small
259        * numbers may be to be entered into the matrix.
260        */
261       MGMatrixSimple(double threshold = 1.e-12);
262
263       /**
264        * Store the result matrix for later assembling.
265        */
266       void
267       initialize(MGLevelObject<MatrixType> &m);
268
269       /**
270        * Initialize the multilevel constraints.
271        */
272       void
273       initialize(const MGConstrainedDoFs &mg_constrained_dofs);
274
275       /**
276        * Initialize the matrices #flux_up and #flux_down used for local
277        * refinement with discontinuous Galerkin methods.
278        */
279       void
280       initialize_fluxes(MGLevelObject<MatrixType> &flux_up,
281                         MGLevelObject<MatrixType> &flux_down);
282
283       /**
284        * Initialize the matrices #interface_in and #interface_out used for
285        * local refinement with continuous Galerkin methods.
286        */
287       void
288       initialize_interfaces(MGLevelObject<MatrixType> &interface_in,
289                             MGLevelObject<MatrixType> &interface_out);
290       /**
291        * Initialize the local data in the DoFInfo object used later for
292        * assembling.
293        *
294        * The @p info object refers to a cell if <code>!face</code>, or else to an
295        * interior or boundary face.
296        */
297       template <class DOFINFO>
298       void
299       initialize_info(DOFINFO &info, bool face) const;
300
301       /**
302        * Assemble the matrix DoFInfo::M1[0] into the global matrix.
303        */
304       template <class DOFINFO>
305       void
306       assemble(const DOFINFO &info);
307
308       /**
309        * Assemble both local matrices in the @p info1 and @p info2
310        * objects into the global matrices.
311        */
312       template <class DOFINFO>
313       void
314       assemble(const DOFINFO &info1, const DOFINFO &info2);
315
316     private:
317       /**
318        * Assemble a single matrix into a global matrix.
319        */
320       void
321       assemble(MatrixType &                                G,
322                const FullMatrix<double> &                  M,
323                const std::vector<types::global_dof_index> &i1,
324                const std::vector<types::global_dof_index> &i2);
325
326       /**
327        * Assemble a single matrix into a global matrix.
328        */
329       void
330       assemble(MatrixType &                                G,
331                const FullMatrix<double> &                  M,
332                const std::vector<types::global_dof_index> &i1,
333                const std::vector<types::global_dof_index> &i2,
334                const unsigned int                          level);
335
336       /**
337        * Assemble a single matrix into a global matrix.
338        */
339
340       void
341       assemble_up(MatrixType &                                G,
342                   const FullMatrix<double> &                  M,
343                   const std::vector<types::global_dof_index> &i1,
344                   const std::vector<types::global_dof_index> &i2,
345                   const unsigned int level = numbers::invalid_unsigned_int);
346       /**
347        * Assemble a single matrix into a global matrix.
348        */
349
350       void
351       assemble_down(MatrixType &                                G,
352                     const FullMatrix<double> &                  M,
353                     const std::vector<types::global_dof_index> &i1,
354                     const std::vector<types::global_dof_index> &i2,
355                     const unsigned int level = numbers::invalid_unsigned_int);
356
357       /**
358        * Assemble a single matrix into a global matrix.
359        */
360
361       void
362       assemble_in(MatrixType &                                G,
363                   const FullMatrix<double> &                  M,
364                   const std::vector<types::global_dof_index> &i1,
365                   const std::vector<types::global_dof_index> &i2,
366                   const unsigned int level = numbers::invalid_unsigned_int);
367
368       /**
369        * Assemble a single matrix into a global matrix.
370        */
371
372       void
373       assemble_out(MatrixType &                                G,
374                    const FullMatrix<double> &                  M,
375                    const std::vector<types::global_dof_index> &i1,
376                    const std::vector<types::global_dof_index> &i2,
377                    const unsigned int level = numbers::invalid_unsigned_int);
378
379       /**
380        * The global matrix being assembled.
381        */
382       SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
383         matrix;
384
385       /**
386        * The matrix used for face flux terms across the refinement edge,
387        * coupling coarse to fine.
388        */
389       SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
390         flux_up;
391
392       /**
393        * The matrix used for face flux terms across the refinement edge,
394        * coupling fine to coarse.
395        */
396       SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
397         flux_down;
398
399       /**
400        * The matrix used for face contributions for continuous elements across
401        * the refinement edge, coupling coarse to fine.
402        */
403       SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
404         interface_in;
405
406       /**
407        * The matrix used for face contributions for continuous elements across
408        * the refinement edge, coupling fine to coarse.
409        */
410       SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
411         interface_out;
412       /**
413        * A pointer to the object containing constraints.
414        */
415       SmartPointer<const MGConstrainedDoFs, MGMatrixSimple<MatrixType>>
416         mg_constrained_dofs;
417
418       /**
419        * The smallest positive number that will be entered into the global
420        * matrix. All smaller absolute values will be treated as zero and will
421        * not be assembled.
422        */
423       const double threshold;
424     };
425
426
427     /**
428      * Assemble a simple matrix and a simple right hand side at once. We use a
429      * combination of MatrixSimple and ResidualSimple to achieve this. Cell
430      * and face operators should fill the matrix and vector objects in
431      * LocalResults and this class will assemble them into matrix and vector
432      * objects.
433      *
434      * @ingroup MeshWorker
435      */
436     template <typename MatrixType, typename VectorType>
437     class SystemSimple : private MatrixSimple<MatrixType>,
438                          private ResidualSimple<VectorType>
439     {
440     public:
441       /**
442        * Constructor setting the threshold value in MatrixSimple.
443        */
444       SystemSimple(double threshold = 1.e-12);
445
446       /**
447        * Store the two objects data is assembled into.
448        */
449       void
450       initialize(MatrixType &m, VectorType &rhs);
451
452       /**
453        * Initialize the constraints. After this function has been called with
454        * a valid AffineConstraints object, the function
455        * AffineConstraints::distribute_local_to_global() will be used by
456        * assemble() to distribute the cell and face matrices into a global
457        * sparse matrix.
458        */
459       void
460       initialize(
461         const AffineConstraints<typename VectorType::value_type> &constraints);
462
463       /**
464        * Initialize the local data in the DoFInfo object used later for
465        * assembling.
466        *
467        * The @p info object refers to a cell if <code>!face</code>, or else to an
468        * interior or boundary face.
469        */
470       template <class DOFINFO>
471       void
472       initialize_info(DOFINFO &info, bool face) const;
473
474       /**
475        * Assemble the matrix DoFInfo::M1[0] into the global matrix.
476        */
477       template <class DOFINFO>
478       void
479       assemble(const DOFINFO &info);
480
481       /**
482        * Assemble both local matrices in the @p info1 and @p info2
483        * objects into the global matrix.
484        */
485       template <class DOFINFO>
486       void
487       assemble(const DOFINFO &info1, const DOFINFO &info2);
488
489     private:
490       /**
491        * Assemble a single matrix <code>M</code> into the element at
492        * <code>index</code> in the vector #matrix.
493        */
494       void
495       assemble(const FullMatrix<double> &                  M,
496                const Vector<double> &                      vector,
497                const unsigned int                          index,
498                const std::vector<types::global_dof_index> &indices);
499
500       void
501       assemble(const FullMatrix<double> &                  M,
502                const Vector<double> &                      vector,
503                const unsigned int                          index,
504                const std::vector<types::global_dof_index> &i1,
505                const std::vector<types::global_dof_index> &i2);
506     };
507
508
509     //----------------------------------------------------------------------//
510
511     template <typename VectorType>
512     inline void
513     ResidualSimple<VectorType>::initialize(AnyData &results)
514     {
515       residuals = results;
516     }
517
518
519
520     template <typename VectorType>
521     inline void
522     ResidualSimple<VectorType>::initialize(
523       const AffineConstraints<typename VectorType::value_type> &c)
524     {
525       constraints = &c;
526     }
527
528
529
530     template <typename VectorType>
531     template <class DOFINFO>
532     inline void
533     ResidualSimple<VectorType>::initialize_info(DOFINFO &info, bool) const
534     {
535       info.initialize_vectors(residuals.size());
536     }
537
538
539
540     template <typename VectorType>
541     template <class DOFINFO>
542     inline void
543     ResidualSimple<VectorType>::assemble(const DOFINFO &info)
544     {
545       for (unsigned int k = 0; k < residuals.size(); ++k)
546         {
547           VectorType *v = residuals.entry<VectorType *>(k);
548           for (unsigned int i = 0; i != info.vector(k).n_blocks(); ++i)
549             {
550               const std::vector<types::global_dof_index> &ldi =
551                 info.vector(k).n_blocks() == 1 ? info.indices :
552                                                  info.indices_by_block[i];
553
554               if (constraints != nullptr)
555                 constraints->distribute_local_to_global(info.vector(k).block(i),
556                                                         ldi,
557                                                         *v);
558               else
559                 v->add(ldi, info.vector(k).block(i));
560             }
561         }
562     }
563
564     template <typename VectorType>
565     template <class DOFINFO>
566     inline void
567     ResidualSimple<VectorType>::assemble(const DOFINFO &info1,
568                                          const DOFINFO &info2)
569     {
570       assemble(info1);
571       assemble(info2);
572     }
573
574
575     //----------------------------------------------------------------------//
576
577     template <typename MatrixType>
578     inline MatrixSimple<MatrixType>::MatrixSimple(double threshold)
579       : threshold(threshold)
580     {}
581
582
583     template <typename MatrixType>
584     inline void
585     MatrixSimple<MatrixType>::initialize(MatrixType &m)
586     {
587       matrix.resize(1);
588       matrix[0] = &m;
589     }
590
591
592     template <typename MatrixType>
593     inline void
594     MatrixSimple<MatrixType>::initialize(std::vector<MatrixType> &m)
595     {
596       matrix.resize(m.size());
597       for (unsigned int i = 0; i < m.size(); ++i)
598         matrix[i] = &m[i];
599     }
600
601
602     template <typename MatrixType>
603     inline void
604     MatrixSimple<MatrixType>::initialize(
605       const AffineConstraints<typename MatrixType::value_type> &c)
606     {
607       constraints = &c;
608     }
609
610
611     template <typename MatrixType>
612     template <class DOFINFO>
613     inline void
614     MatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
615     {
616       Assert(matrix.size() != 0, ExcNotInitialized());
617
618       const unsigned int n = info.indices_by_block.size();
619
620       if (n == 0)
621         info.initialize_matrices(matrix.size(), face);
622       else
623         {
624           info.initialize_matrices(matrix.size() * n * n, face);
625           unsigned int k = 0;
626           for (unsigned int m = 0; m < matrix.size(); ++m)
627             for (unsigned int i = 0; i < n; ++i)
628               for (unsigned int j = 0; j < n; ++j, ++k)
629                 {
630                   info.matrix(k, false).row    = i;
631                   info.matrix(k, false).column = j;
632                   if (face)
633                     {
634                       info.matrix(k, true).row    = i;
635                       info.matrix(k, true).column = j;
636                     }
637                 }
638         }
639     }
640
641
642
643     template <typename MatrixType>
644     inline void
645     MatrixSimple<MatrixType>::assemble(
646       const FullMatrix<double> &                  M,
647       const unsigned int                          index,
648       const std::vector<types::global_dof_index> &i1,
649       const std::vector<types::global_dof_index> &i2)
650     {
651       AssertDimension(M.m(), i1.size());
652       AssertDimension(M.n(), i2.size());
653
654       if (constraints == nullptr)
655         {
656           for (unsigned int j = 0; j < i1.size(); ++j)
657             for (unsigned int k = 0; k < i2.size(); ++k)
658               if (std::fabs(M(j, k)) >= threshold)
659                 matrix[index]->add(i1[j], i2[k], M(j, k));
660         }
661       else
662         constraints->distribute_local_to_global(M, i1, i2, *matrix[index]);
663     }
664
665
666     template <typename MatrixType>
667     template <class DOFINFO>
668     inline void
669     MatrixSimple<MatrixType>::assemble(const DOFINFO &info)
670     {
671       Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
672       const unsigned int n = info.indices_by_block.size();
673
674       if (n == 0)
675         for (unsigned int m = 0; m < matrix.size(); ++m)
676           assemble(info.matrix(m, false).matrix, m, info.indices, info.indices);
677       else
678         {
679           for (unsigned int m = 0; m < matrix.size(); ++m)
680             for (unsigned int k = 0; k < n * n; ++k)
681               {
682                 assemble(
683                   info.matrix(k + m * n * n, false).matrix,
684                   m,
685                   info.indices_by_block[info.matrix(k + m * n * n, false).row],
686                   info.indices_by_block[info.matrix(k + m * n * n, false)
687                                           .column]);
688               }
689         }
690     }
691
692
693     template <typename MatrixType>
694     template <class DOFINFO>
695     inline void
696     MatrixSimple<MatrixType>::assemble(const DOFINFO &info1,
697                                        const DOFINFO &info2)
698     {
699       Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
700       Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
701       AssertDimension(info1.indices_by_block.size(),
702                       info2.indices_by_block.size());
703
704       const unsigned int n = info1.indices_by_block.size();
705
706       if (n == 0)
707         {
708           for (unsigned int m = 0; m < matrix.size(); ++m)
709             {
710               assemble(info1.matrix(m, false).matrix,
711                        m,
712                        info1.indices,
713                        info1.indices);
714               assemble(info1.matrix(m, true).matrix,
715                        m,
716                        info1.indices,
717                        info2.indices);
718               assemble(info2.matrix(m, false).matrix,
719                        m,
720                        info2.indices,
721                        info2.indices);
722               assemble(info2.matrix(m, true).matrix,
723                        m,
724                        info2.indices,
725                        info1.indices);
726             }
727         }
728       else
729         {
730           for (unsigned int m = 0; m < matrix.size(); ++m)
731             for (unsigned int k = 0; k < n * n; ++k)
732               {
733                 const unsigned int row = info1.matrix(k + m * n * n, false).row;
734                 const unsigned int column =
735                   info1.matrix(k + m * n * n, false).column;
736
737                 assemble(info1.matrix(k + m * n * n, false).matrix,
738                          m,
739                          info1.indices_by_block[row],
740                          info1.indices_by_block[column]);
741                 assemble(info1.matrix(k + m * n * n, true).matrix,
742                          m,
743                          info1.indices_by_block[row],
744                          info2.indices_by_block[column]);
745                 assemble(info2.matrix(k + m * n * n, false).matrix,
746                          m,
747                          info2.indices_by_block[row],
748                          info2.indices_by_block[column]);
749                 assemble(info2.matrix(k + m * n * n, true).matrix,
750                          m,
751                          info2.indices_by_block[row],
752                          info1.indices_by_block[column]);
753               }
754         }
755     }
756
757
758     //----------------------------------------------------------------------//
759
760     template <typename MatrixType>
761     inline MGMatrixSimple<MatrixType>::MGMatrixSimple(double threshold)
762       : threshold(threshold)
763     {}
764
765
766     template <typename MatrixType>
767     inline void
768     MGMatrixSimple<MatrixType>::initialize(MGLevelObject<MatrixType> &m)
769     {
770       matrix = &m;
771     }
772
773     template <typename MatrixType>
774     inline void
775     MGMatrixSimple<MatrixType>::initialize(const MGConstrainedDoFs &c)
776     {
777       mg_constrained_dofs = &c;
778     }
779
780
781     template <typename MatrixType>
782     inline void
783     MGMatrixSimple<MatrixType>::initialize_fluxes(
784       MGLevelObject<MatrixType> &up,
785       MGLevelObject<MatrixType> &down)
786     {
787       flux_up   = &up;
788       flux_down = &down;
789     }
790
791
792     template <typename MatrixType>
793     inline void
794     MGMatrixSimple<MatrixType>::initialize_interfaces(
795       MGLevelObject<MatrixType> &in,
796       MGLevelObject<MatrixType> &out)
797     {
798       interface_in  = &in;
799       interface_out = &out;
800     }
801
802
803     template <typename MatrixType>
804     template <class DOFINFO>
805     inline void
806     MGMatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
807     {
808       const unsigned int n = info.indices_by_block.size();
809
810       if (n == 0)
811         info.initialize_matrices(1, face);
812       else
813         {
814           info.initialize_matrices(n * n, face);
815           unsigned int k = 0;
816           for (unsigned int i = 0; i < n; ++i)
817             for (unsigned int j = 0; j < n; ++j, ++k)
818               {
819                 info.matrix(k, false).row    = i;
820                 info.matrix(k, false).column = j;
821                 if (face)
822                   {
823                     info.matrix(k, true).row    = i;
824                     info.matrix(k, true).column = j;
825                   }
826               }
827         }
828     }
829
830
831     template <typename MatrixType>
832     inline void
833     MGMatrixSimple<MatrixType>::assemble(
834       MatrixType &                                G,
835       const FullMatrix<double> &                  M,
836       const std::vector<types::global_dof_index> &i1,
837       const std::vector<types::global_dof_index> &i2)
838     {
839       AssertDimension(M.m(), i1.size());
840       AssertDimension(M.n(), i2.size());
841       Assert(mg_constrained_dofs == 0, ExcInternalError());
842       // TODO: Possibly remove this function all together
843
844       for (unsigned int j = 0; j < i1.size(); ++j)
845         for (unsigned int k = 0; k < i2.size(); ++k)
846           if (std::fabs(M(j, k)) >= threshold)
847             G.add(i1[j], i2[k], M(j, k));
848     }
849
850
851     template <typename MatrixType>
852     inline void
853     MGMatrixSimple<MatrixType>::assemble(
854       MatrixType &                                G,
855       const FullMatrix<double> &                  M,
856       const std::vector<types::global_dof_index> &i1,
857       const std::vector<types::global_dof_index> &i2,
858       const unsigned int                          level)
859     {
860       AssertDimension(M.m(), i1.size());
861       AssertDimension(M.n(), i2.size());
862
863       if (mg_constrained_dofs == nullptr)
864         {
865           for (unsigned int j = 0; j < i1.size(); ++j)
866             for (unsigned int k = 0; k < i2.size(); ++k)
867               if (std::fabs(M(j, k)) >= threshold)
868                 G.add(i1[j], i2[k], M(j, k));
869         }
870       else
871         {
872           for (unsigned int j = 0; j < i1.size(); ++j)
873             for (unsigned int k = 0; k < i2.size(); ++k)
874               {
875                 // Only enter the local values into the global matrix,
876                 //  if the value is larger than the threshold
877                 if (std::fabs(M(j, k)) < threshold)
878                   continue;
879
880                 // Do not enter, if either the row or the column
881                 // corresponds to an index on the refinement edge. The
882                 // level problems are solved with homogeneous
883                 // Dirichlet boundary conditions, therefore we
884                 // eliminate these rows and columns. The corresponding
885                 // matrix entries are entered by assemble_in() and
886                 // assemble_out().
887                 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
888                     mg_constrained_dofs->at_refinement_edge(level, i2[k]))
889                   continue;
890
891                 // At the boundary, only enter the term on the
892                 // diagonal, but not the coupling terms
893                 if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
894                      mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
895                     (i1[j] != i2[k]))
896                   continue;
897
898                 G.add(i1[j], i2[k], M(j, k));
899               }
900         }
901     }
902
903
904     template <typename MatrixType>
905     inline void
906     MGMatrixSimple<MatrixType>::assemble_up(
907       MatrixType &                                G,
908       const FullMatrix<double> &                  M,
909       const std::vector<types::global_dof_index> &i1,
910       const std::vector<types::global_dof_index> &i2,
911       const unsigned int                          level)
912     {
913       AssertDimension(M.n(), i1.size());
914       AssertDimension(M.m(), i2.size());
915
916       if (mg_constrained_dofs == nullptr)
917         {
918           for (unsigned int j = 0; j < i1.size(); ++j)
919             for (unsigned int k = 0; k < i2.size(); ++k)
920               if (std::fabs(M(k, j)) >= threshold)
921                 G.add(i1[j], i2[k], M(k, j));
922         }
923       else
924         {
925           for (unsigned int j = 0; j < i1.size(); ++j)
926             for (unsigned int k = 0; k < i2.size(); ++k)
927               if (std::fabs(M(k, j)) >= threshold)
928                 if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
929                   G.add(i1[j], i2[k], M(k, j));
930         }
931     }
932
933     template <typename MatrixType>
934     inline void
935     MGMatrixSimple<MatrixType>::assemble_down(
936       MatrixType &                                G,
937       const FullMatrix<double> &                  M,
938       const std::vector<types::global_dof_index> &i1,
939       const std::vector<types::global_dof_index> &i2,
940       const unsigned int                          level)
941     {
942       AssertDimension(M.m(), i1.size());
943       AssertDimension(M.n(), i2.size());
944
945       if (mg_constrained_dofs == nullptr)
946         {
947           for (unsigned int j = 0; j < i1.size(); ++j)
948             for (unsigned int k = 0; k < i2.size(); ++k)
949               if (std::fabs(M(j, k)) >= threshold)
950                 G.add(i1[j], i2[k], M(j, k));
951         }
952       else
953         {
954           for (unsigned int j = 0; j < i1.size(); ++j)
955             for (unsigned int k = 0; k < i2.size(); ++k)
956               if (std::fabs(M(j, k)) >= threshold)
957                 if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
958                   G.add(i1[j], i2[k], M(j, k));
959         }
960     }
961
962     template <typename MatrixType>
963     inline void
964     MGMatrixSimple<MatrixType>::assemble_in(
965       MatrixType &                                G,
966       const FullMatrix<double> &                  M,
967       const std::vector<types::global_dof_index> &i1,
968       const std::vector<types::global_dof_index> &i2,
969       const unsigned int                          level)
970     {
971       AssertDimension(M.m(), i1.size());
972       AssertDimension(M.n(), i2.size());
973       Assert(mg_constrained_dofs != nullptr, ExcInternalError());
974
975       for (unsigned int j = 0; j < i1.size(); ++j)
976         for (unsigned int k = 0; k < i2.size(); ++k)
977           if (std::fabs(M(j, k)) >= threshold)
978             // Enter values into matrix only if j corresponds to a
979             // degree of freedom on the refinement edge, k does
980             // not, and both are not on the boundary. This is part
981             // the difference between the complete matrix with no
982             // boundary condition at the refinement edge and
983             // the matrix assembled above by assemble().
984
985             // Thus the logic is: enter the row if it is
986             // constrained by hanging node constraints (actually,
987             // the whole refinement edge), but not if it is
988             // constrained by a boundary constraint.
989             if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
990                 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
991               {
992                 if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
993                      !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
994                     (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
995                      mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
996                      i1[j] == i2[k]))
997                   G.add(i1[j], i2[k], M(j, k));
998               }
999     }
1000
1001
1002     template <typename MatrixType>
1003     inline void
1004     MGMatrixSimple<MatrixType>::assemble_out(
1005       MatrixType &                                G,
1006       const FullMatrix<double> &                  M,
1007       const std::vector<types::global_dof_index> &i1,
1008       const std::vector<types::global_dof_index> &i2,
1009       const unsigned int                          level)
1010     {
1011       AssertDimension(M.n(), i1.size());
1012       AssertDimension(M.m(), i2.size());
1013       Assert(mg_constrained_dofs != nullptr, ExcInternalError());
1014
1015       for (unsigned int j = 0; j < i1.size(); ++j)
1016         for (unsigned int k = 0; k < i2.size(); ++k)
1017           if (std::fabs(M(k, j)) >= threshold)
1018             if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
1019                 !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
1020               {
1021                 if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1022                      !mg_constrained_dofs->is_boundary_index(level, i2[k])) ||
1023                     (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
1024                      mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
1025                      i1[j] == i2[k]))
1026                   G.add(i1[j], i2[k], M(k, j));
1027               }
1028     }
1029
1030
1031     template <typename MatrixType>
1032     template <class DOFINFO>
1033     inline void
1034     MGMatrixSimple<MatrixType>::assemble(const DOFINFO &info)
1035     {
1036       Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
1037       const unsigned int level = info.cell->level();
1038
1039       if (info.indices_by_block.size() == 0)
1040         {
1041           assemble((*matrix)[level],
1042                    info.matrix(0, false).matrix,
1043                    info.indices,
1044                    info.indices,
1045                    level);
1046           if (mg_constrained_dofs != nullptr)
1047             {
1048               assemble_in((*interface_in)[level],
1049                           info.matrix(0, false).matrix,
1050                           info.indices,
1051                           info.indices,
1052                           level);
1053               assemble_out((*interface_out)[level],
1054                            info.matrix(0, false).matrix,
1055                            info.indices,
1056                            info.indices,
1057                            level);
1058             }
1059         }
1060       else
1061         for (unsigned int k = 0; k < info.n_matrices(); ++k)
1062           {
1063             const unsigned int row    = info.matrix(k, false).row;
1064             const unsigned int column = info.matrix(k, false).column;
1065
1066             assemble((*matrix)[level],
1067                      info.matrix(k, false).matrix,
1068                      info.indices_by_block[row],
1069                      info.indices_by_block[column],
1070                      level);
1071
1072             if (mg_constrained_dofs != nullptr)
1073               {
1074                 assemble_in((*interface_in)[level],
1075                             info.matrix(k, false).matrix,
1076                             info.indices_by_block[row],
1077                             info.indices_by_block[column],
1078                             level);
1079                 assemble_out((*interface_out)[level],
1080                              info.matrix(k, false).matrix,
1081                              info.indices_by_block[column],
1082                              info.indices_by_block[row],
1083                              level);
1084               }
1085           }
1086     }
1087
1088
1089     template <typename MatrixType>
1090     template <class DOFINFO>
1091     inline void
1092     MGMatrixSimple<MatrixType>::assemble(const DOFINFO &info1,
1093                                          const DOFINFO &info2)
1094     {
1095       Assert(info1.level_cell, ExcMessage("Cell must access level dofs"));
1096       Assert(info2.level_cell, ExcMessage("Cell must access level dofs"));
1097       const unsigned int level1 = info1.cell->level();
1098       const unsigned int level2 = info2.cell->level();
1099
1100       if (info1.indices_by_block.size() == 0)
1101         {
1102           if (level1 == level2)
1103             {
1104               assemble((*matrix)[level1],
1105                        info1.matrix(0, false).matrix,
1106                        info1.indices,
1107                        info1.indices,
1108                        level1);
1109               assemble((*matrix)[level1],
1110                        info1.matrix(0, true).matrix,
1111                        info1.indices,
1112                        info2.indices,
1113                        level1);
1114               assemble((*matrix)[level1],
1115                        info2.matrix(0, false).matrix,
1116                        info2.indices,
1117                        info2.indices,
1118                        level1);
1119               assemble((*matrix)[level1],
1120                        info2.matrix(0, true).matrix,
1121                        info2.indices,
1122                        info1.indices,
1123                        level1);
1124             }
1125           else
1126             {
1127               Assert(level1 > level2, ExcInternalError());
1128               // Do not add info2.M1,
1129               // which is done by
1130               // the coarser cell
1131               assemble((*matrix)[level1],
1132                        info1.matrix(0, false).matrix,
1133                        info1.indices,
1134                        info1.indices,
1135                        level1);
1136               if (level1 > 0)
1137                 {
1138                   assemble_up((*flux_up)[level1],
1139                               info1.matrix(0, true).matrix,
1140                               info2.indices,
1141                               info1.indices,
1142                               level1);
1143                   assemble_down((*flux_down)[level1],
1144                                 info2.matrix(0, true).matrix,
1145                                 info2.indices,
1146                                 info1.indices,
1147                                 level1);
1148                 }
1149             }
1150         }
1151       else
1152         for (unsigned int k = 0; k < info1.n_matrices(); ++k)
1153           {
1154             const unsigned int row    = info1.matrix(k, false).row;
1155             const unsigned int column = info1.matrix(k, false).column;
1156
1157             if (level1 == level2)
1158               {
1159                 assemble((*matrix)[level1],
1160                          info1.matrix(k, false).matrix,
1161                          info1.indices_by_block[row],
1162                          info1.indices_by_block[column],
1163                          level1);
1164                 assemble((*matrix)[level1],
1165                          info1.matrix(k, true).matrix,
1166                          info1.indices_by_block[row],
1167                          info2.indices_by_block[column],
1168                          level1);
1169                 assemble((*matrix)[level1],
1170                          info2.matrix(k, false).matrix,
1171                          info2.indices_by_block[row],
1172                          info2.indices_by_block[column],
1173                          level1);
1174                 assemble((*matrix)[level1],
1175                          info2.matrix(k, true).matrix,
1176                          info2.indices_by_block[row],
1177                          info1.indices_by_block[column],
1178                          level1);
1179               }
1180             else
1181               {
1182                 Assert(level1 > level2, ExcInternalError());
1183                 // Do not add info2.M1,
1184                 // which is done by
1185                 // the coarser cell
1186                 assemble((*matrix)[level1],
1187                          info1.matrix(k, false).matrix,
1188                          info1.indices_by_block[row],
1189                          info1.indices_by_block[column],
1190                          level1);
1191                 if (level1 > 0)
1192                   {
1193                     assemble_up((*flux_up)[level1],
1194                                 info1.matrix(k, true).matrix,
1195                                 info2.indices_by_block[column],
1196                                 info1.indices_by_block[row],
1197                                 level1);
1198                     assemble_down((*flux_down)[level1],
1199                                   info2.matrix(k, true).matrix,
1200                                   info2.indices_by_block[row],
1201                                   info1.indices_by_block[column],
1202                                   level1);
1203                   }
1204               }
1205           }
1206     }
1207
1208     //----------------------------------------------------------------------//
1209
1210     template <typename MatrixType, typename VectorType>
1211     SystemSimple<MatrixType, VectorType>::SystemSimple(double t)
1212       : MatrixSimple<MatrixType>(t)
1213     {}
1214
1215
1216     template <typename MatrixType, typename VectorType>
1217     inline void
1218     SystemSimple<MatrixType, VectorType>::initialize(MatrixType &m,
1219                                                      VectorType &rhs)
1220     {
1221       AnyData     data;
1222       VectorType *p = &rhs;
1223       data.add(p, "right hand side");
1224
1225       MatrixSimple<MatrixType>::initialize(m);
1226       ResidualSimple<VectorType>::initialize(data);
1227     }
1228
1229     template <typename MatrixType, typename VectorType>
1230     inline void
1231     SystemSimple<MatrixType, VectorType>::initialize(
1232       const AffineConstraints<typename VectorType::value_type> &c)
1233     {
1234       ResidualSimple<VectorType>::initialize(c);
1235     }
1236
1237
1238     template <typename MatrixType, typename VectorType>
1239     template <class DOFINFO>
1240     inline void
1241     SystemSimple<MatrixType, VectorType>::initialize_info(DOFINFO &info,
1242                                                           bool     face) const
1243     {
1244       MatrixSimple<MatrixType>::initialize_info(info, face);
1245       ResidualSimple<VectorType>::initialize_info(info, face);
1246     }
1247
1248     template <typename MatrixType, typename VectorType>
1249     inline void
1250     SystemSimple<MatrixType, VectorType>::assemble(
1251       const FullMatrix<double> &                  M,
1252       const Vector<double> &                      vector,
1253       const unsigned int                          index,
1254       const std::vector<types::global_dof_index> &indices)
1255     {
1256       AssertDimension(M.m(), indices.size());
1257       AssertDimension(M.n(), indices.size());
1258
1259       AnyData     residuals = ResidualSimple<VectorType>::residuals;
1260       VectorType *v         = residuals.entry<VectorType *>(index);
1261
1262       if (ResidualSimple<VectorType>::constraints == nullptr)
1263         {
1264           for (unsigned int i = 0; i < indices.size(); ++i)
1265             (*v)(indices[i]) += vector(i);
1266
1267           for (unsigned int j = 0; j < indices.size(); ++j)
1268             for (unsigned int k = 0; k < indices.size(); ++k)
1269               if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
1270                 MatrixSimple<MatrixType>::matrix[index]->add(indices[j],
1271                                                              indices[k],
1272                                                              M(j, k));
1273         }
1274       else
1275         {
1276           ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1277             M,
1278             vector,
1279             indices,
1280             *MatrixSimple<MatrixType>::matrix[index],
1281             *v,
1282             true);
1283         }
1284     }
1285
1286     template <typename MatrixType, typename VectorType>
1287     inline void
1288     SystemSimple<MatrixType, VectorType>::assemble(
1289       const FullMatrix<double> &                  M,
1290       const Vector<double> &                      vector,
1291       const unsigned int                          index,
1292       const std::vector<types::global_dof_index> &i1,
1293       const std::vector<types::global_dof_index> &i2)
1294     {
1295       AssertDimension(M.m(), i1.size());
1296       AssertDimension(M.n(), i2.size());
1297
1298       AnyData     residuals = ResidualSimple<VectorType>::residuals;
1299       VectorType *v         = residuals.entry<VectorType *>(index);
1300
1301       if (ResidualSimple<VectorType>::constraints == nullptr)
1302         {
1303           for (unsigned int j = 0; j < i1.size(); ++j)
1304             for (unsigned int k = 0; k < i2.size(); ++k)
1305               if (std::fabs(M(j, k)) >= MatrixSimple<MatrixType>::threshold)
1306                 MatrixSimple<MatrixType>::matrix[index]->add(i1[j],
1307                                                              i2[k],
1308                                                              M(j, k));
1309         }
1310       else
1311         {
1312           ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1313             vector, i1, i2, *v, M, false);
1314           ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1315             M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
1316         }
1317     }
1318
1319
1320     template <typename MatrixType, typename VectorType>
1321     template <class DOFINFO>
1322     inline void
1323     SystemSimple<MatrixType, VectorType>::assemble(const DOFINFO &info)
1324     {
1325       AssertDimension(MatrixSimple<MatrixType>::matrix.size(),
1326                       ResidualSimple<VectorType>::residuals.size());
1327       Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
1328       const unsigned int n = info.indices_by_block.size();
1329
1330       if (n == 0)
1331         {
1332           for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1333                ++m)
1334             assemble(info.matrix(m, false).matrix,
1335                      info.vector(m).block(0),
1336                      m,
1337                      info.indices);
1338         }
1339       else
1340         {
1341           for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1342                ++m)
1343             for (unsigned int k = 0; k < n * n; ++k)
1344               {
1345                 const unsigned int row = info.matrix(k + m * n * n, false).row;
1346                 const unsigned int column =
1347                   info.matrix(k + m * n * n, false).column;
1348
1349                 if (row == column)
1350                   assemble(info.matrix(k + m * n * n, false).matrix,
1351                            info.vector(m).block(row),
1352                            m,
1353                            info.indices_by_block[row]);
1354                 else
1355                   assemble(info.matrix(k + m * n * n, false).matrix,
1356                            info.vector(m).block(row),
1357                            m,
1358                            info.indices_by_block[row],
1359                            info.indices_by_block[column]);
1360               }
1361         }
1362     }
1363
1364
1365     template <typename MatrixType, typename VectorType>
1366     template <class DOFINFO>
1367     inline void
1368     SystemSimple<MatrixType, VectorType>::assemble(const DOFINFO &info1,
1369                                                    const DOFINFO &info2)
1370     {
1371       Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs"));
1372       Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs"));
1373       AssertDimension(info1.indices_by_block.size(),
1374                       info2.indices_by_block.size());
1375
1376       const unsigned int n = info1.indices_by_block.size();
1377
1378       if (n == 0)
1379         {
1380           for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1381                ++m)
1382             {
1383               assemble(info1.matrix(m, false).matrix,
1384                        info1.vector(m).block(0),
1385                        m,
1386                        info1.indices);
1387               assemble(info1.matrix(m, true).matrix,
1388                        info1.vector(m).block(0),
1389                        m,
1390                        info1.indices,
1391                        info2.indices);
1392               assemble(info2.matrix(m, false).matrix,
1393                        info2.vector(m).block(0),
1394                        m,
1395                        info2.indices);
1396               assemble(info2.matrix(m, true).matrix,
1397                        info2.vector(m).block(0),
1398                        m,
1399                        info2.indices,
1400                        info1.indices);
1401             }
1402         }
1403       else
1404         {
1405           for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1406                ++m)
1407             for (unsigned int k = 0; k < n * n; ++k)
1408               {
1409                 const unsigned int row = info1.matrix(k + m * n * n, false).row;
1410                 const unsigned int column =
1411                   info1.matrix(k + m * n * n, false).column;
1412
1413                 if (row == column)
1414                   {
1415                     assemble(info1.matrix(k + m * n * n, false).matrix,
1416                              info1.vector(m).block(row),
1417                              m,
1418                              info1.indices_by_block[row]);
1419                     assemble(info2.matrix(k + m * n * n, false).matrix,
1420                              info2.vector(m).block(row),
1421                              m,
1422                              info2.indices_by_block[row]);
1423                   }
1424                 else
1425                   {
1426                     assemble(info1.matrix(k + m * n * n, false).matrix,
1427                              info1.vector(m).block(row),
1428                              m,
1429                              info1.indices_by_block[row],
1430                              info1.indices_by_block[column]);
1431                     assemble(info2.matrix(k + m * n * n, false).matrix,
1432                              info2.vector(m).block(row),
1433                              m,
1434                              info2.indices_by_block[row],
1435                              info2.indices_by_block[column]);
1436                   }
1437                 assemble(info1.matrix(k + m * n * n, true).matrix,
1438                          info1.vector(m).block(row),
1439                          m,
1440                          info1.indices_by_block[row],
1441                          info2.indices_by_block[column]);
1442                 assemble(info2.matrix(k + m * n * n, true).matrix,
1443                          info2.vector(m).block(row),
1444                          m,
1445                          info2.indices_by_block[row],
1446                          info1.indices_by_block[column]);
1447               }
1448         }
1449     }
1450   } // namespace Assembler
1451 } // namespace MeshWorker
1452
1453 DEAL_II_NAMESPACE_CLOSE
1454
1455 #endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.