1 // ---------------------------------------------------------------------
3 // Copyright (C) 2010 - 2020 by the deal.II authors
5 // This file is part of the deal.II library.
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.
14 // ---------------------------------------------------------------------
17 #ifndef dealii_mesh_worker_simple_h
18 #define dealii_mesh_worker_simple_h
20 #include <deal.II/base/config.h>
22 #include <deal.II/algorithms/any_data.h>
24 #include <deal.II/base/mg_level_object.h>
25 #include <deal.II/base/smartpointer.h>
27 #include <deal.II/lac/block_vector.h>
29 #include <deal.II/meshworker/dof_info.h>
30 #include <deal.II/meshworker/functional.h>
32 #include <deal.II/multigrid/mg_constrained_dofs.h>
35 * The header containing the classes MeshWorker::Assembler::MatrixSimple,
36 * MeshWorker::Assembler::MGMatrixSimple, MeshWorker::Assembler::ResidualSimple,
37 * and MeshWorker::Assembler::SystemSimple.
40 DEAL_II_NAMESPACE_OPEN
47 * Assemble residuals without block structure.
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()
58 template <typename VectorType>
63 * Initialize with an AnyData object holding the result of assembling.
65 * Assembling currently writes into the first vector of
69 initialize(AnyData &results);
72 * Initialize the constraints.
76 const AffineConstraints<typename VectorType::value_type> &constraints);
79 * Initialize the local data in the DoFInfo object used later for
82 * The @p info object refers to a cell if <code>!face</code>, or else to an
83 * interior or boundary face.
85 template <class DOFINFO>
87 initialize_info(DOFINFO &info, bool face) const;
90 * Assemble the local residuals into the global residuals.
92 * Values are added to the previous contents. If constraints are active,
93 * AffineConstraints::distribute_local_to_global() is used.
95 template <class DOFINFO>
97 assemble(const DOFINFO &info);
100 * Assemble both local residuals into the global residuals.
102 template <class DOFINFO>
104 assemble(const DOFINFO &info1, const DOFINFO &info2);
108 * The global residual vectors filled by assemble().
113 * A pointer to the object containing constraints.
115 SmartPointer<const AffineConstraints<typename VectorType::value_type>,
116 ResidualSimple<VectorType>>
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.
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.
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.
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,
149 * @ingroup MeshWorker
151 template <typename MatrixType>
156 * Constructor, initializing the #threshold, which limits how small
157 * numbers may be to be entered into the matrix.
159 MatrixSimple(double threshold = 1.e-12);
162 * Store the result matrix for later assembling.
165 initialize(MatrixType &m);
168 * Store several result matrices for later assembling.
171 initialize(std::vector<MatrixType> &m);
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
182 const AffineConstraints<typename MatrixType::value_type> &constraints);
185 * Initialize the local data in the DoFInfo object used later for
188 * The @p info object refers to a cell if <code>!face</code>, or else to an
189 * interior or boundary face.
191 template <class DOFINFO>
193 initialize_info(DOFINFO &info, bool face) const;
196 * Assemble the local matrices associated with a single cell into the
199 template <class DOFINFO>
201 assemble(const DOFINFO &info);
204 * Assemble all local matrices associated with an interior face in the
205 * @p info1 and @p info2 objects into the global matrix.
207 template <class DOFINFO>
209 assemble(const DOFINFO &info1, const DOFINFO &info2);
213 * The vector of global matrices being assembled.
215 std::vector<SmartPointer<MatrixType, MatrixSimple<MatrixType>>> matrix;
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
222 const double threshold;
226 * Assemble a single matrix <code>M</code> into the element at
227 * <code>index</code> in the vector #matrix.
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);
236 * A pointer to the object containing constraints.
238 SmartPointer<const AffineConstraints<typename MatrixType::value_type>,
239 MatrixSimple<MatrixType>>
245 * Assemble local matrices into level matrices without using block
248 * @todo The matrix structures needed for assembling level matrices with
249 * local refinement and continuous elements are missing.
251 * @ingroup MeshWorker
253 template <typename MatrixType>
258 * Constructor, initializing the #threshold, which limits how small
259 * numbers may be to be entered into the matrix.
261 MGMatrixSimple(double threshold = 1.e-12);
264 * Store the result matrix for later assembling.
267 initialize(MGLevelObject<MatrixType> &m);
270 * Initialize the multilevel constraints.
273 initialize(const MGConstrainedDoFs &mg_constrained_dofs);
276 * Initialize the matrices #flux_up and #flux_down used for local
277 * refinement with discontinuous Galerkin methods.
280 initialize_fluxes(MGLevelObject<MatrixType> &flux_up,
281 MGLevelObject<MatrixType> &flux_down);
284 * Initialize the matrices #interface_in and #interface_out used for
285 * local refinement with continuous Galerkin methods.
288 initialize_interfaces(MGLevelObject<MatrixType> &interface_in,
289 MGLevelObject<MatrixType> &interface_out);
291 * Initialize the local data in the DoFInfo object used later for
294 * The @p info object refers to a cell if <code>!face</code>, or else to an
295 * interior or boundary face.
297 template <class DOFINFO>
299 initialize_info(DOFINFO &info, bool face) const;
302 * Assemble the matrix DoFInfo::M1[0] into the global matrix.
304 template <class DOFINFO>
306 assemble(const DOFINFO &info);
309 * Assemble both local matrices in the @p info1 and @p info2
310 * objects into the global matrices.
312 template <class DOFINFO>
314 assemble(const DOFINFO &info1, const DOFINFO &info2);
318 * Assemble a single matrix into a global matrix.
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);
327 * Assemble a single matrix into a global matrix.
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);
337 * Assemble a single matrix into a global matrix.
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);
347 * Assemble a single matrix into a global matrix.
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);
358 * Assemble a single matrix into a global matrix.
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);
369 * Assemble a single matrix into a global matrix.
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);
380 * The global matrix being assembled.
382 SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
386 * The matrix used for face flux terms across the refinement edge,
387 * coupling coarse to fine.
389 SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
393 * The matrix used for face flux terms across the refinement edge,
394 * coupling fine to coarse.
396 SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
400 * The matrix used for face contributions for continuous elements across
401 * the refinement edge, coupling coarse to fine.
403 SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
407 * The matrix used for face contributions for continuous elements across
408 * the refinement edge, coupling fine to coarse.
410 SmartPointer<MGLevelObject<MatrixType>, MGMatrixSimple<MatrixType>>
413 * A pointer to the object containing constraints.
415 SmartPointer<const MGConstrainedDoFs, MGMatrixSimple<MatrixType>>
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
423 const double threshold;
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
434 * @ingroup MeshWorker
436 template <typename MatrixType, typename VectorType>
437 class SystemSimple : private MatrixSimple<MatrixType>,
438 private ResidualSimple<VectorType>
442 * Constructor setting the threshold value in MatrixSimple.
444 SystemSimple(double threshold = 1.e-12);
447 * Store the two objects data is assembled into.
450 initialize(MatrixType &m, VectorType &rhs);
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
461 const AffineConstraints<typename VectorType::value_type> &constraints);
464 * Initialize the local data in the DoFInfo object used later for
467 * The @p info object refers to a cell if <code>!face</code>, or else to an
468 * interior or boundary face.
470 template <class DOFINFO>
472 initialize_info(DOFINFO &info, bool face) const;
475 * Assemble the matrix DoFInfo::M1[0] into the global matrix.
477 template <class DOFINFO>
479 assemble(const DOFINFO &info);
482 * Assemble both local matrices in the @p info1 and @p info2
483 * objects into the global matrix.
485 template <class DOFINFO>
487 assemble(const DOFINFO &info1, const DOFINFO &info2);
491 * Assemble a single matrix <code>M</code> into the element at
492 * <code>index</code> in the vector #matrix.
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);
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);
509 //----------------------------------------------------------------------//
511 template <typename VectorType>
513 ResidualSimple<VectorType>::initialize(AnyData &results)
520 template <typename VectorType>
522 ResidualSimple<VectorType>::initialize(
523 const AffineConstraints<typename VectorType::value_type> &c)
530 template <typename VectorType>
531 template <class DOFINFO>
533 ResidualSimple<VectorType>::initialize_info(DOFINFO &info, bool) const
535 info.initialize_vectors(residuals.size());
540 template <typename VectorType>
541 template <class DOFINFO>
543 ResidualSimple<VectorType>::assemble(const DOFINFO &info)
545 for (unsigned int k = 0; k < residuals.size(); ++k)
547 VectorType *v = residuals.entry<VectorType *>(k);
548 for (unsigned int i = 0; i != info.vector(k).n_blocks(); ++i)
550 const std::vector<types::global_dof_index> &ldi =
551 info.vector(k).n_blocks() == 1 ? info.indices :
552 info.indices_by_block[i];
554 if (constraints != nullptr)
555 constraints->distribute_local_to_global(info.vector(k).block(i),
559 v->add(ldi, info.vector(k).block(i));
564 template <typename VectorType>
565 template <class DOFINFO>
567 ResidualSimple<VectorType>::assemble(const DOFINFO &info1,
568 const DOFINFO &info2)
575 //----------------------------------------------------------------------//
577 template <typename MatrixType>
578 inline MatrixSimple<MatrixType>::MatrixSimple(double threshold)
579 : threshold(threshold)
583 template <typename MatrixType>
585 MatrixSimple<MatrixType>::initialize(MatrixType &m)
592 template <typename MatrixType>
594 MatrixSimple<MatrixType>::initialize(std::vector<MatrixType> &m)
596 matrix.resize(m.size());
597 for (unsigned int i = 0; i < m.size(); ++i)
602 template <typename MatrixType>
604 MatrixSimple<MatrixType>::initialize(
605 const AffineConstraints<typename MatrixType::value_type> &c)
611 template <typename MatrixType>
612 template <class DOFINFO>
614 MatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
616 Assert(matrix.size() != 0, ExcNotInitialized());
618 const unsigned int n = info.indices_by_block.size();
621 info.initialize_matrices(matrix.size(), face);
624 info.initialize_matrices(matrix.size() * n * n, face);
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)
630 info.matrix(k, false).row = i;
631 info.matrix(k, false).column = j;
634 info.matrix(k, true).row = i;
635 info.matrix(k, true).column = j;
643 template <typename MatrixType>
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)
651 AssertDimension(M.m(), i1.size());
652 AssertDimension(M.n(), i2.size());
654 if (constraints == nullptr)
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));
662 constraints->distribute_local_to_global(M, i1, i2, *matrix[index]);
666 template <typename MatrixType>
667 template <class DOFINFO>
669 MatrixSimple<MatrixType>::assemble(const DOFINFO &info)
671 Assert(!info.level_cell, ExcMessage("Cell may not access level dofs"));
672 const unsigned int n = info.indices_by_block.size();
675 for (unsigned int m = 0; m < matrix.size(); ++m)
676 assemble(info.matrix(m, false).matrix, m, info.indices, info.indices);
679 for (unsigned int m = 0; m < matrix.size(); ++m)
680 for (unsigned int k = 0; k < n * n; ++k)
683 info.matrix(k + m * n * n, false).matrix,
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)
693 template <typename MatrixType>
694 template <class DOFINFO>
696 MatrixSimple<MatrixType>::assemble(const DOFINFO &info1,
697 const DOFINFO &info2)
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());
704 const unsigned int n = info1.indices_by_block.size();
708 for (unsigned int m = 0; m < matrix.size(); ++m)
710 assemble(info1.matrix(m, false).matrix,
714 assemble(info1.matrix(m, true).matrix,
718 assemble(info2.matrix(m, false).matrix,
722 assemble(info2.matrix(m, true).matrix,
730 for (unsigned int m = 0; m < matrix.size(); ++m)
731 for (unsigned int k = 0; k < n * n; ++k)
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;
737 assemble(info1.matrix(k + m * n * n, false).matrix,
739 info1.indices_by_block[row],
740 info1.indices_by_block[column]);
741 assemble(info1.matrix(k + m * n * n, true).matrix,
743 info1.indices_by_block[row],
744 info2.indices_by_block[column]);
745 assemble(info2.matrix(k + m * n * n, false).matrix,
747 info2.indices_by_block[row],
748 info2.indices_by_block[column]);
749 assemble(info2.matrix(k + m * n * n, true).matrix,
751 info2.indices_by_block[row],
752 info1.indices_by_block[column]);
758 //----------------------------------------------------------------------//
760 template <typename MatrixType>
761 inline MGMatrixSimple<MatrixType>::MGMatrixSimple(double threshold)
762 : threshold(threshold)
766 template <typename MatrixType>
768 MGMatrixSimple<MatrixType>::initialize(MGLevelObject<MatrixType> &m)
773 template <typename MatrixType>
775 MGMatrixSimple<MatrixType>::initialize(const MGConstrainedDoFs &c)
777 mg_constrained_dofs = &c;
781 template <typename MatrixType>
783 MGMatrixSimple<MatrixType>::initialize_fluxes(
784 MGLevelObject<MatrixType> &up,
785 MGLevelObject<MatrixType> &down)
792 template <typename MatrixType>
794 MGMatrixSimple<MatrixType>::initialize_interfaces(
795 MGLevelObject<MatrixType> &in,
796 MGLevelObject<MatrixType> &out)
799 interface_out = &out;
803 template <typename MatrixType>
804 template <class DOFINFO>
806 MGMatrixSimple<MatrixType>::initialize_info(DOFINFO &info, bool face) const
808 const unsigned int n = info.indices_by_block.size();
811 info.initialize_matrices(1, face);
814 info.initialize_matrices(n * n, face);
816 for (unsigned int i = 0; i < n; ++i)
817 for (unsigned int j = 0; j < n; ++j, ++k)
819 info.matrix(k, false).row = i;
820 info.matrix(k, false).column = j;
823 info.matrix(k, true).row = i;
824 info.matrix(k, true).column = j;
831 template <typename MatrixType>
833 MGMatrixSimple<MatrixType>::assemble(
835 const FullMatrix<double> & M,
836 const std::vector<types::global_dof_index> &i1,
837 const std::vector<types::global_dof_index> &i2)
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
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));
851 template <typename MatrixType>
853 MGMatrixSimple<MatrixType>::assemble(
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)
860 AssertDimension(M.m(), i1.size());
861 AssertDimension(M.n(), i2.size());
863 if (mg_constrained_dofs == nullptr)
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));
872 for (unsigned int j = 0; j < i1.size(); ++j)
873 for (unsigned int k = 0; k < i2.size(); ++k)
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)
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
887 if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
888 mg_constrained_dofs->at_refinement_edge(level, i2[k]))
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])) &&
898 G.add(i1[j], i2[k], M(j, k));
904 template <typename MatrixType>
906 MGMatrixSimple<MatrixType>::assemble_up(
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)
913 AssertDimension(M.n(), i1.size());
914 AssertDimension(M.m(), i2.size());
916 if (mg_constrained_dofs == nullptr)
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));
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));
933 template <typename MatrixType>
935 MGMatrixSimple<MatrixType>::assemble_down(
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)
942 AssertDimension(M.m(), i1.size());
943 AssertDimension(M.n(), i2.size());
945 if (mg_constrained_dofs == nullptr)
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));
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));
962 template <typename MatrixType>
964 MGMatrixSimple<MatrixType>::assemble_in(
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)
971 AssertDimension(M.m(), i1.size());
972 AssertDimension(M.n(), i2.size());
973 Assert(mg_constrained_dofs != nullptr, ExcInternalError());
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().
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]))
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]) &&
997 G.add(i1[j], i2[k], M(j, k));
1002 template <typename MatrixType>
1004 MGMatrixSimple<MatrixType>::assemble_out(
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)
1011 AssertDimension(M.n(), i1.size());
1012 AssertDimension(M.m(), i2.size());
1013 Assert(mg_constrained_dofs != nullptr, ExcInternalError());
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]))
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]) &&
1026 G.add(i1[j], i2[k], M(k, j));
1031 template <typename MatrixType>
1032 template <class DOFINFO>
1034 MGMatrixSimple<MatrixType>::assemble(const DOFINFO &info)
1036 Assert(info.level_cell, ExcMessage("Cell must access level dofs"));
1037 const unsigned int level = info.cell->level();
1039 if (info.indices_by_block.size() == 0)
1041 assemble((*matrix)[level],
1042 info.matrix(0, false).matrix,
1046 if (mg_constrained_dofs != nullptr)
1048 assemble_in((*interface_in)[level],
1049 info.matrix(0, false).matrix,
1053 assemble_out((*interface_out)[level],
1054 info.matrix(0, false).matrix,
1061 for (unsigned int k = 0; k < info.n_matrices(); ++k)
1063 const unsigned int row = info.matrix(k, false).row;
1064 const unsigned int column = info.matrix(k, false).column;
1066 assemble((*matrix)[level],
1067 info.matrix(k, false).matrix,
1068 info.indices_by_block[row],
1069 info.indices_by_block[column],
1072 if (mg_constrained_dofs != nullptr)
1074 assemble_in((*interface_in)[level],
1075 info.matrix(k, false).matrix,
1076 info.indices_by_block[row],
1077 info.indices_by_block[column],
1079 assemble_out((*interface_out)[level],
1080 info.matrix(k, false).matrix,
1081 info.indices_by_block[column],
1082 info.indices_by_block[row],
1089 template <typename MatrixType>
1090 template <class DOFINFO>
1092 MGMatrixSimple<MatrixType>::assemble(const DOFINFO &info1,
1093 const DOFINFO &info2)
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();
1100 if (info1.indices_by_block.size() == 0)
1102 if (level1 == level2)
1104 assemble((*matrix)[level1],
1105 info1.matrix(0, false).matrix,
1109 assemble((*matrix)[level1],
1110 info1.matrix(0, true).matrix,
1114 assemble((*matrix)[level1],
1115 info2.matrix(0, false).matrix,
1119 assemble((*matrix)[level1],
1120 info2.matrix(0, true).matrix,
1127 Assert(level1 > level2, ExcInternalError());
1128 // Do not add info2.M1,
1131 assemble((*matrix)[level1],
1132 info1.matrix(0, false).matrix,
1138 assemble_up((*flux_up)[level1],
1139 info1.matrix(0, true).matrix,
1143 assemble_down((*flux_down)[level1],
1144 info2.matrix(0, true).matrix,
1152 for (unsigned int k = 0; k < info1.n_matrices(); ++k)
1154 const unsigned int row = info1.matrix(k, false).row;
1155 const unsigned int column = info1.matrix(k, false).column;
1157 if (level1 == level2)
1159 assemble((*matrix)[level1],
1160 info1.matrix(k, false).matrix,
1161 info1.indices_by_block[row],
1162 info1.indices_by_block[column],
1164 assemble((*matrix)[level1],
1165 info1.matrix(k, true).matrix,
1166 info1.indices_by_block[row],
1167 info2.indices_by_block[column],
1169 assemble((*matrix)[level1],
1170 info2.matrix(k, false).matrix,
1171 info2.indices_by_block[row],
1172 info2.indices_by_block[column],
1174 assemble((*matrix)[level1],
1175 info2.matrix(k, true).matrix,
1176 info2.indices_by_block[row],
1177 info1.indices_by_block[column],
1182 Assert(level1 > level2, ExcInternalError());
1183 // Do not add info2.M1,
1186 assemble((*matrix)[level1],
1187 info1.matrix(k, false).matrix,
1188 info1.indices_by_block[row],
1189 info1.indices_by_block[column],
1193 assemble_up((*flux_up)[level1],
1194 info1.matrix(k, true).matrix,
1195 info2.indices_by_block[column],
1196 info1.indices_by_block[row],
1198 assemble_down((*flux_down)[level1],
1199 info2.matrix(k, true).matrix,
1200 info2.indices_by_block[row],
1201 info1.indices_by_block[column],
1208 //----------------------------------------------------------------------//
1210 template <typename MatrixType, typename VectorType>
1211 SystemSimple<MatrixType, VectorType>::SystemSimple(double t)
1212 : MatrixSimple<MatrixType>(t)
1216 template <typename MatrixType, typename VectorType>
1218 SystemSimple<MatrixType, VectorType>::initialize(MatrixType &m,
1222 VectorType *p = &rhs;
1223 data.add(p, "right hand side");
1225 MatrixSimple<MatrixType>::initialize(m);
1226 ResidualSimple<VectorType>::initialize(data);
1229 template <typename MatrixType, typename VectorType>
1231 SystemSimple<MatrixType, VectorType>::initialize(
1232 const AffineConstraints<typename VectorType::value_type> &c)
1234 ResidualSimple<VectorType>::initialize(c);
1238 template <typename MatrixType, typename VectorType>
1239 template <class DOFINFO>
1241 SystemSimple<MatrixType, VectorType>::initialize_info(DOFINFO &info,
1244 MatrixSimple<MatrixType>::initialize_info(info, face);
1245 ResidualSimple<VectorType>::initialize_info(info, face);
1248 template <typename MatrixType, typename VectorType>
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)
1256 AssertDimension(M.m(), indices.size());
1257 AssertDimension(M.n(), indices.size());
1259 AnyData residuals = ResidualSimple<VectorType>::residuals;
1260 VectorType *v = residuals.entry<VectorType *>(index);
1262 if (ResidualSimple<VectorType>::constraints == nullptr)
1264 for (unsigned int i = 0; i < indices.size(); ++i)
1265 (*v)(indices[i]) += vector(i);
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],
1276 ResidualSimple<VectorType>::constraints->distribute_local_to_global(
1280 *MatrixSimple<MatrixType>::matrix[index],
1286 template <typename MatrixType, typename VectorType>
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)
1295 AssertDimension(M.m(), i1.size());
1296 AssertDimension(M.n(), i2.size());
1298 AnyData residuals = ResidualSimple<VectorType>::residuals;
1299 VectorType *v = residuals.entry<VectorType *>(index);
1301 if (ResidualSimple<VectorType>::constraints == nullptr)
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],
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]);
1320 template <typename MatrixType, typename VectorType>
1321 template <class DOFINFO>
1323 SystemSimple<MatrixType, VectorType>::assemble(const DOFINFO &info)
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();
1332 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1334 assemble(info.matrix(m, false).matrix,
1335 info.vector(m).block(0),
1341 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1343 for (unsigned int k = 0; k < n * n; ++k)
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;
1350 assemble(info.matrix(k + m * n * n, false).matrix,
1351 info.vector(m).block(row),
1353 info.indices_by_block[row]);
1355 assemble(info.matrix(k + m * n * n, false).matrix,
1356 info.vector(m).block(row),
1358 info.indices_by_block[row],
1359 info.indices_by_block[column]);
1365 template <typename MatrixType, typename VectorType>
1366 template <class DOFINFO>
1368 SystemSimple<MatrixType, VectorType>::assemble(const DOFINFO &info1,
1369 const DOFINFO &info2)
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());
1376 const unsigned int n = info1.indices_by_block.size();
1380 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1383 assemble(info1.matrix(m, false).matrix,
1384 info1.vector(m).block(0),
1387 assemble(info1.matrix(m, true).matrix,
1388 info1.vector(m).block(0),
1392 assemble(info2.matrix(m, false).matrix,
1393 info2.vector(m).block(0),
1396 assemble(info2.matrix(m, true).matrix,
1397 info2.vector(m).block(0),
1405 for (unsigned int m = 0; m < MatrixSimple<MatrixType>::matrix.size();
1407 for (unsigned int k = 0; k < n * n; ++k)
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;
1415 assemble(info1.matrix(k + m * n * n, false).matrix,
1416 info1.vector(m).block(row),
1418 info1.indices_by_block[row]);
1419 assemble(info2.matrix(k + m * n * n, false).matrix,
1420 info2.vector(m).block(row),
1422 info2.indices_by_block[row]);
1426 assemble(info1.matrix(k + m * n * n, false).matrix,
1427 info1.vector(m).block(row),
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),
1434 info2.indices_by_block[row],
1435 info2.indices_by_block[column]);
1437 assemble(info1.matrix(k + m * n * n, true).matrix,
1438 info1.vector(m).block(row),
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),
1445 info2.indices_by_block[row],
1446 info1.indices_by_block[column]);
1450 } // namespace Assembler
1451 } // namespace MeshWorker
1453 DEAL_II_NAMESPACE_CLOSE