| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Discretization | ||
| 10 | * \brief Contains functionality for L2-projections from one function | ||
| 11 | * space into another, which can live both on the same or different | ||
| 12 | * grids of potentially different dimensionality. | ||
| 13 | * \note In the case of the function space bases living on the same grid | ||
| 14 | * an optimized implementation could be implemented. Currently, we | ||
| 15 | * require a glue object containing the intersections between two | ||
| 16 | * grids to be provided to the free factory functions. | ||
| 17 | */ | ||
| 18 | #ifndef DUMUX_DISCRETIZATION_PROJECTOR_HH | ||
| 19 | #define DUMUX_DISCRETIZATION_PROJECTOR_HH | ||
| 20 | |||
| 21 | #include <algorithm> | ||
| 22 | #include <string> | ||
| 23 | #include <utility> | ||
| 24 | #include <type_traits> | ||
| 25 | |||
| 26 | #include <dune/common/timer.hh> | ||
| 27 | #include <dune/common/fmatrix.hh> | ||
| 28 | #include <dune/common/exceptions.hh> | ||
| 29 | #include <dune/common/promotiontraits.hh> | ||
| 30 | #include <dune/common/parametertree.hh> | ||
| 31 | |||
| 32 | #include <dune/geometry/quadraturerules.hh> | ||
| 33 | #include <dune/istl/matrixindexset.hh> | ||
| 34 | #include <dune/istl/bcrsmatrix.hh> | ||
| 35 | #include <dune/istl/bvector.hh> | ||
| 36 | |||
| 37 | #include <dumux/common/parameters.hh> | ||
| 38 | #include <dumux/linear/linearsolvertraits.hh> | ||
| 39 | #include <dumux/linear/linearalgebratraits.hh> | ||
| 40 | #include <dumux/linear/istlsolvers.hh> | ||
| 41 | #include <dumux/assembly/jacobianpattern.hh> | ||
| 42 | #include <dumux/discretization/functionspacebasis.hh> | ||
| 43 | |||
| 44 | namespace Dumux { | ||
| 45 | |||
| 46 | /*! | ||
| 47 | * \ingroup Discretization | ||
| 48 | * \brief Does an L2-projection from one discrete function space | ||
| 49 | * into another. The convenience functions makeProjectorPair | ||
| 50 | * or makeProjector can be used to create such a projection. | ||
| 51 | */ | ||
| 52 | template<class ScalarType> | ||
| 53 | class Projector | ||
| 54 | { | ||
| 55 | using MatrixBlockType = Dune::FieldMatrix<ScalarType, 1, 1>; | ||
| 56 | |||
| 57 | public: | ||
| 58 | //! Export the scalar type | ||
| 59 | using Scalar = ScalarType; | ||
| 60 | //! Export the type of the projection matrices | ||
| 61 | using Matrix = Dune::BCRSMatrix< MatrixBlockType >; | ||
| 62 | |||
| 63 | //! Parameters that can be passed to project() | ||
| 64 | struct Params | ||
| 65 | { | ||
| 66 | std::size_t maxIterations{100}; | ||
| 67 | Scalar residualReduction{1e-13}; | ||
| 68 | int verbosity{0}; | ||
| 69 | }; | ||
| 70 | |||
| 71 | //! delete default constructor | ||
| 72 | Projector() = delete; | ||
| 73 | |||
| 74 | /*! | ||
| 75 | * \brief Constructor. Receives the mass and projection | ||
| 76 | * matrix that define the linear system describing | ||
| 77 | * the L2-projection from a function space into another. | ||
| 78 | */ | ||
| 79 | 29 | Projector(Matrix&& massMatrix, Matrix&& projectionMatrix) | |
| 80 | 29 | : massMat_(std::make_shared<Matrix>(std::move(massMatrix))) | |
| 81 | 29 | , projMat_(std::make_shared<Matrix>(std::move(projectionMatrix))) | |
| 82 |
6/18✓ Branch 2 taken 29 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 29 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 29 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 29 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 29 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 29 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
87 | , numDofsTarget_(massMat_->N()) |
| 83 | { | ||
| 84 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 29 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 29 times.
|
87 | if (massMat_->N() != projMat_->N()) |
| 85 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_->N() << " vs " << projMat_->N()); | |
| 86 | |||
| 87 |
2/8✓ Branch 2 taken 29 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 29 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
29 | massMatrixSolver_.setMatrix(massMat_); |
| 88 | 29 | } | |
| 89 | |||
| 90 | /*! | ||
| 91 | * \brief Constructor for projection into a target space that occupies | ||
| 92 | * a larger geometric region than the domain space. In this case, | ||
| 93 | * the mass matrix can be chosen such that is is solved for only | ||
| 94 | * those dofs which will be populated with values from the projection. | ||
| 95 | * This requires an additional index map that maps the entries of the | ||
| 96 | * projected solution into the solution vector for the target space. | ||
| 97 | * Furthermore, the number of degrees of freedom must be specified to | ||
| 98 | * set up the coefficient vector with correct size for the target space. | ||
| 99 | */ | ||
| 100 | 1 | Projector(Matrix&& massMatrix, | |
| 101 | Matrix&& projectionMatrix, | ||
| 102 | std::vector<std::size_t>&& indexMap, | ||
| 103 | std::size_t numDofsTarget) | ||
| 104 | 1 | : massMat_(std::make_shared<Matrix>(std::move(massMatrix))) | |
| 105 | 1 | , projMat_(std::make_shared<Matrix>(std::move(projectionMatrix))) | |
| 106 | 1 | , indexMapTarget_(std::move(indexMap)) | |
| 107 |
5/16✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
|
2 | , numDofsTarget_(numDofsTarget) |
| 108 | { | ||
| 109 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
3 | if (indexMapTarget_.size() != massMat_->N()) |
| 110 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Target index map size mismatch: " << indexMapTarget_.size() << " vs " << massMat_->N()); | |
| 111 | |||
| 112 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
|
3 | if (massMat_->N() != projMat_->N()) |
| 113 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_->N() << " vs " << projMat_->N()); | |
| 114 | |||
| 115 |
4/8✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
|
4 | if (*std::max_element(indexMapTarget_.begin(), indexMapTarget_.end()) > numDofsTarget_) |
| 116 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Index map exceeds provided number of dofs in target domain!"); | |
| 117 | |||
| 118 |
2/8✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
1 | massMatrixSolver_.setMatrix(massMat_); |
| 119 | 1 | } | |
| 120 | |||
| 121 | /*! | ||
| 122 | * \brief Project a solution u into up | ||
| 123 | * \param u The solution living on the domain space | ||
| 124 | * \param params Optional parameters for mass matrix solve | ||
| 125 | * \return The projection of u into the target space | ||
| 126 | */ | ||
| 127 | template< class BlockType, std::enable_if_t<std::is_convertible<BlockType, ScalarType>::value, int> = 0 > | ||
| 128 | 18 | Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const | |
| 129 | { | ||
| 130 | // be picky about size of u | ||
| 131 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
|
36 | if ( u.size() != projMat_->M()) |
| 132 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Vector size mismatch"); | |
| 133 | |||
| 134 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
54 | Dune::BlockVector<BlockType> up(massMat_->N()); |
| 135 | |||
| 136 |
3/6✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
36 | auto rhs = up; |
| 137 | 36 | projMat_->mv(u, rhs); | |
| 138 | |||
| 139 |
2/4✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
|
36 | Dune::ParameterTree solverParams; |
| 140 |
7/18✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 18 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 18 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 18 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 18 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
18 | solverParams["maxit"] = std::to_string(params.maxIterations); |
| 141 |
7/18✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 18 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 18 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 18 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 18 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
18 | solverParams["reduction"] = std::to_string(params.residualReduction); |
| 142 |
7/18✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 18 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 18 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 18 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 18 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
|
18 | solverParams["verbose"] = std::to_string(params.verbosity); |
| 143 |
1/2✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
|
36 | auto solver = massMatrixSolver_; // copy the solver to modify the parameters |
| 144 |
3/8✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
18 | solver.setParams(solverParams); |
| 145 |
1/2✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
|
18 | solver.solve(up, rhs); |
| 146 | |||
| 147 | // if target space occupies a larger region, fill missing entries with zero | ||
| 148 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 17 times.
|
36 | if (!indexMapTarget_.empty()) |
| 149 | { | ||
| 150 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | Dune::BlockVector<BlockType> result(numDofsTarget_); |
| 151 | |||
| 152 | 1 | result = 0.0; | |
| 153 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 1 times.
|
44 | for (std::size_t i = 0; i < indexMapTarget_.size(); ++i) |
| 154 | 84 | result[indexMapTarget_[i]] = up[i]; | |
| 155 | |||
| 156 | 2 | return result; | |
| 157 | } | ||
| 158 | |||
| 159 | 17 | return up; | |
| 160 | } | ||
| 161 | |||
| 162 | /*! | ||
| 163 | * \brief Project a solution u into up | ||
| 164 | * \param u The solution living on the domain space | ||
| 165 | * \param params Optional parameters for mass matrix solve | ||
| 166 | * \return The projection of u into the target space | ||
| 167 | */ | ||
| 168 | template< class BlockType, std::enable_if_t<!std::is_convertible<BlockType, ScalarType>::value, int> = 0 > | ||
| 169 | 2 | Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const | |
| 170 | { | ||
| 171 | 2 | Dune::BlockVector<BlockType> result(numDofsTarget_); | |
| 172 | |||
| 173 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | for (int pvIdx = 0; pvIdx < BlockType::size(); ++pvIdx) |
| 174 | { | ||
| 175 |
1/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
8 | Dune::BlockVector<Dune::FieldVector<Scalar, 1>> tmp(u.size()); |
| 176 | 22596 | std::transform(u.begin(), u.end(), tmp.begin(), [pvIdx] (const auto& v) { return v[pvIdx]; }); | |
| 177 | |||
| 178 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
|
12 | const auto p = project(tmp, params); |
| 179 |
2/2✓ Branch 0 taken 1124 times.
✓ Branch 1 taken 4 times.
|
1128 | for (std::size_t i = 0; i < p.size(); ++i) |
| 180 | 4496 | result[i][pvIdx] = p[i]; | |
| 181 | } | ||
| 182 | |||
| 183 | 2 | return result; | |
| 184 | } | ||
| 185 | |||
| 186 | /*! | ||
| 187 | * \brief Returns the default parameters. | ||
| 188 | */ | ||
| 189 | static Params defaultParams() | ||
| 190 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
14 | { return {}; } |
| 191 | |||
| 192 | private: | ||
| 193 | std::shared_ptr<Matrix> massMat_; | ||
| 194 | std::shared_ptr<Matrix> projMat_; | ||
| 195 | |||
| 196 | SSORCGIstlSolver< | ||
| 197 | SeqLinearSolverTraits, LinearAlgebraTraits<Matrix, Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> | ||
| 198 | > massMatrixSolver_; | ||
| 199 | |||
| 200 | std::vector<std::size_t> indexMapTarget_; | ||
| 201 | std::size_t numDofsTarget_; | ||
| 202 | }; | ||
| 203 | |||
| 204 | /*! | ||
| 205 | * \brief Traits class stating the type of projector between to bases | ||
| 206 | */ | ||
| 207 | template<class FEBasisDomain, class FEBasisTarget> | ||
| 208 | class ProjectorTraits | ||
| 209 | { | ||
| 210 | using FiniteElementDomain = typename FEBasisDomain::LocalView::Tree::FiniteElement; | ||
| 211 | using FiniteElementTarget = typename FEBasisTarget::LocalView::Tree::FiniteElement; | ||
| 212 | using ScalarDomain = typename FiniteElementDomain::Traits::LocalBasisType::Traits::RangeFieldType; | ||
| 213 | using ScalarTarget = typename FiniteElementTarget::Traits::LocalBasisType::Traits::RangeFieldType; | ||
| 214 | using Scalar = typename Dune::PromotionTraits<ScalarDomain, ScalarTarget>::PromotedType; | ||
| 215 | public: | ||
| 216 | using Projector = Dumux::Projector< Scalar >; | ||
| 217 | }; | ||
| 218 | |||
| 219 | |||
| 220 | // Projector construction details | ||
| 221 | namespace Detail { | ||
| 222 | |||
| 223 | /*! | ||
| 224 | * \brief Reduces a mass matrix and projection matrix such that they are composed | ||
| 225 | * of only those dofs that actually take part in the projection. Simultaneously, | ||
| 226 | * a container with the index map into the complete target space is filled so that | ||
| 227 | * the entries after projection can be assigned to the corresponding dof in the | ||
| 228 | * overall target space. | ||
| 229 | */ | ||
| 230 | template<class Matrix> | ||
| 231 | 1 | void setupReducedMatrices(const Matrix& massMatrix, const Matrix& projMatrix, const std::vector<bool>& dofIsVoid, | |
| 232 | Matrix& reducedM, Matrix& reducedP, std::vector<std::size_t>& expansionMap) | ||
| 233 | { | ||
| 234 | 3 | const std::size_t numNonVoidDofs = std::count_if(dofIsVoid.begin(), dofIsVoid.end(), [] (bool v) { return !v; }); | |
| 235 | |||
| 236 | // reduce matrices to only dofs that take part and create index map | ||
| 237 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2 | std::vector<std::size_t> reductionMap(massMatrix.N()); |
| 238 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | expansionMap.resize(numNonVoidDofs); |
| 239 | |||
| 240 | std::size_t idxInReducedSpace = 0; | ||
| 241 |
4/4✓ Branch 0 taken 441 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 441 times.
✓ Branch 3 taken 1 times.
|
443 | for (std::size_t dofIdx = 0; dofIdx < dofIsVoid.size(); ++dofIdx) |
| 242 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 420 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
|
882 | if (!dofIsVoid[dofIdx]) |
| 243 | { | ||
| 244 | 21 | reductionMap[dofIdx] = idxInReducedSpace; | |
| 245 | 21 | expansionMap[idxInReducedSpace] = dofIdx; | |
| 246 | 21 | idxInReducedSpace++; | |
| 247 | } | ||
| 248 | |||
| 249 | // create reduced M/P matrix | ||
| 250 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
|
3 | Dune::MatrixIndexSet patternMReduced, patternPReduced; |
| 251 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | patternMReduced.resize(numNonVoidDofs, numNonVoidDofs); |
| 252 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | patternPReduced.resize(numNonVoidDofs, projMatrix.M()); |
| 253 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 442 times.
✓ Branch 4 taken 441 times.
✓ Branch 5 taken 1 times.
|
443 | for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt) |
| 254 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 420 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
|
882 | if (!dofIsVoid[rowIt.index()]) |
| 255 | { | ||
| 256 | 21 | const auto reducedRowIdx = reductionMap[rowIt.index()]; | |
| 257 |
5/8✗ Branch 0 not taken.
✓ Branch 1 taken 426 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 426 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 426 times.
✓ Branch 6 taken 405 times.
✓ Branch 7 taken 21 times.
|
447 | for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt) |
| 258 |
6/6✓ Branch 0 taken 81 times.
✓ Branch 1 taken 324 times.
✓ Branch 2 taken 81 times.
✓ Branch 3 taken 324 times.
✓ Branch 4 taken 81 times.
✓ Branch 5 taken 324 times.
|
1215 | if (!dofIsVoid[colIt.index()]) |
| 259 |
2/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
|
162 | patternMReduced.add(reducedRowIdx, reductionMap[colIt.index()]); |
| 260 | } | ||
| 261 | |||
| 262 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 442 times.
✓ Branch 4 taken 441 times.
✓ Branch 5 taken 1 times.
|
443 | for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt) |
| 263 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 420 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
|
882 | if (!dofIsVoid[rowIt.index()]) |
| 264 | { | ||
| 265 | 21 | const auto reducedRowIdx = reductionMap[rowIt.index()]; | |
| 266 |
5/8✗ Branch 0 not taken.
✓ Branch 1 taken 156 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 156 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 156 times.
✓ Branch 6 taken 135 times.
✓ Branch 7 taken 21 times.
|
447 | for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt) |
| 267 |
2/4✓ Branch 1 taken 135 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 135 times.
✗ Branch 5 not taken.
|
270 | patternPReduced.add(reducedRowIdx, colIt.index()); |
| 268 | } | ||
| 269 | |||
| 270 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | patternMReduced.exportIdx(reducedM); |
| 271 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | patternPReduced.exportIdx(reducedP); |
| 272 | |||
| 273 | // fill matrix entries | ||
| 274 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 442 times.
✓ Branch 4 taken 441 times.
✓ Branch 5 taken 1 times.
|
443 | for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt) |
| 275 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 420 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
|
882 | if (!dofIsVoid[rowIt.index()]) |
| 276 | { | ||
| 277 | 21 | const auto reducedRowIdx = reductionMap[rowIt.index()]; | |
| 278 |
5/8✗ Branch 0 not taken.
✓ Branch 1 taken 426 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 426 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 426 times.
✓ Branch 6 taken 405 times.
✓ Branch 7 taken 21 times.
|
447 | for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt) |
| 279 |
6/6✓ Branch 0 taken 81 times.
✓ Branch 1 taken 324 times.
✓ Branch 2 taken 81 times.
✓ Branch 3 taken 324 times.
✓ Branch 4 taken 81 times.
✓ Branch 5 taken 324 times.
|
1215 | if (!dofIsVoid[colIt.index()]) |
| 280 |
4/8✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 81 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 81 times.
✗ Branch 11 not taken.
|
324 | reducedM[reducedRowIdx][reductionMap[colIt.index()]] = *colIt; |
| 281 | } | ||
| 282 | |||
| 283 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 442 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 442 times.
✓ Branch 4 taken 441 times.
✓ Branch 5 taken 1 times.
|
443 | for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt) |
| 284 |
4/4✓ Branch 0 taken 21 times.
✓ Branch 1 taken 420 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 420 times.
|
882 | if (!dofIsVoid[rowIt.index()]) |
| 285 | { | ||
| 286 | 21 | const auto reducedRowIdx = reductionMap[rowIt.index()]; | |
| 287 |
5/8✗ Branch 0 not taken.
✓ Branch 1 taken 156 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 156 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 156 times.
✓ Branch 6 taken 135 times.
✓ Branch 7 taken 21 times.
|
447 | for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt) |
| 288 |
4/8✓ Branch 1 taken 135 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 135 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 135 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 135 times.
✗ Branch 11 not taken.
|
540 | reducedP[reducedRowIdx][colIt.index()] = *colIt; |
| 289 | } | ||
| 290 | 1 | } | |
| 291 | |||
| 292 | /*! | ||
| 293 | * \brief Creates the matrices underlying l2-projections | ||
| 294 | * \tparam doBidirectional If false, the backward projection matrix is not assembled | ||
| 295 | * \param feBasisDomain The basis to the domain finite element space | ||
| 296 | * \param feBasisTarget The basis to the target finite element space | ||
| 297 | * \param glue The glue object containing the intersections between the two grids | ||
| 298 | * \param treatDiagonalZeroes If true, zero entries on the diagonal of the matrices | ||
| 299 | * that appear if the two domains occupy different geometric regions (and some | ||
| 300 | * dofs to not take part in the projection as a result) are substituted by ones. | ||
| 301 | * This substitution will lead to those dofs being mapped to zeroes in the target space. | ||
| 302 | * \returns An std::pair of projection matrices, where the first entry stores the | ||
| 303 | * matrices of the forward projection and the second entry stores those | ||
| 304 | * of the backward projection. The entries of the returned pair are itself | ||
| 305 | * std::pairs which store the mass matrix in the first and the projection | ||
| 306 | * matrix in the second entry. | ||
| 307 | */ | ||
| 308 | template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType> | ||
| 309 | 17 | auto createProjectionMatrices(const FEBasisDomain& feBasisDomain, | |
| 310 | const FEBasisTarget& feBasisTarget, | ||
| 311 | const GlueType& glue, | ||
| 312 | bool treatDiagonalZeroes = true) | ||
| 313 | { | ||
| 314 | // we assume that target dim <= domain dimension | ||
| 315 | static constexpr int domainDim = FEBasisDomain::GridView::dimension; | ||
| 316 | static constexpr int targetDim = FEBasisTarget::GridView::dimension; | ||
| 317 | static_assert(targetDim <= domainDim, "This expects target dim < domain dim, please swap arguments"); | ||
| 318 | |||
| 319 | using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector; | ||
| 320 | using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector; | ||
| 321 | |||
| 322 | using ForwardProjectionMatrix = typename ForwardProjector::Matrix; | ||
| 323 | using BackwardProjectionMatrix = typename BackwardProjector::Matrix; | ||
| 324 | |||
| 325 |
1/2✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
|
34 | auto domainLocalView = feBasisDomain.localView(); |
| 326 |
1/2✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
|
17 | auto targetLocalView = feBasisTarget.localView(); |
| 327 | |||
| 328 | // determine mass matrix patterns (standard FE scheme pattern) | ||
| 329 |
2/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
|
51 | Dune::MatrixIndexSet backwardPatternM, forwardPatternM; |
| 330 |
1/2✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
|
17 | forwardPatternM = getFEJacobianPattern(feBasisTarget); |
| 331 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | if (doBidirectional) backwardPatternM = getFEJacobianPattern(feBasisDomain); |
| 332 | |||
| 333 | // determine projection matrix patterns | ||
| 334 |
2/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
|
51 | Dune::MatrixIndexSet backwardPatternP, forwardPatternP; |
| 335 |
4/7✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 14 times.
✗ Branch 8 not taken.
|
17 | forwardPatternP.resize(feBasisTarget.size(), feBasisDomain.size()); |
| 336 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
1 | if (doBidirectional) backwardPatternP.resize(feBasisDomain.size(), feBasisTarget.size()); |
| 337 | |||
| 338 | using std::max; | ||
| 339 | 17 | unsigned int maxBasisOrder = 0; | |
| 340 |
4/4✓ Branch 0 taken 321334 times.
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 321334 times.
✓ Branch 3 taken 15 times.
|
335480 | for (const auto& is : intersections(glue)) |
| 341 | { | ||
| 342 | // since target dim <= domain dim there is maximum one! | ||
| 343 |
3/5✓ Branch 1 taken 321312 times.
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 321312 times.
✗ Branch 5 not taken.
|
335446 | targetLocalView.bind( is.targetEntity(0) ); |
| 344 |
1/2✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
|
335446 | const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis(); |
| 345 | |||
| 346 |
4/4✓ Branch 0 taken 321356 times.
✓ Branch 1 taken 321334 times.
✓ Branch 2 taken 321356 times.
✓ Branch 3 taken 321334 times.
|
1006360 | for (unsigned int nIdx = 0; nIdx < is.numDomainNeighbors(); ++nIdx) |
| 347 | { | ||
| 348 |
2/4✓ Branch 1 taken 321356 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 321356 times.
✗ Branch 5 not taken.
|
335468 | domainLocalView.bind( is.domainEntity(nIdx) ); |
| 349 |
1/2✓ Branch 1 taken 321356 times.
✗ Branch 2 not taken.
|
335468 | const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis(); |
| 350 | |||
| 351 | // keep track of maximum basis order (used in integration) | ||
| 352 |
5/8✓ Branch 1 taken 321356 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 321356 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 321356 times.
✓ Branch 8 taken 9 times.
✓ Branch 9 taken 321347 times.
|
335468 | maxBasisOrder = max(maxBasisOrder, max(domainLocalBasis.order(), targetLocalBasis.order())); |
| 353 | |||
| 354 |
3/4✓ Branch 1 taken 1181480 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 860124 times.
✓ Branch 4 taken 321356 times.
|
1287320 | for (unsigned int i = 0; i < domainLocalBasis.size(); ++i) |
| 355 |
3/4✓ Branch 1 taken 4156944 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3296820 times.
✓ Branch 4 taken 860124 times.
|
4933104 | for (unsigned int j = 0; j < targetLocalBasis.size(); ++j) |
| 356 | { | ||
| 357 |
5/10✓ Branch 1 taken 3296820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3296820 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3296820 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3296820 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3296820 times.
✗ Branch 14 not taken.
|
19906260 | forwardPatternP.add(targetLocalView.index(j), domainLocalView.index(i)); |
| 358 |
5/10✓ Branch 1 taken 1188 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1188 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1188 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1188 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1188 times.
✗ Branch 14 not taken.
|
5940 | if (doBidirectional) backwardPatternP.add(domainLocalView.index(i), targetLocalView.index(j)); |
| 359 | } | ||
| 360 | } | ||
| 361 | } | ||
| 362 | |||
| 363 | // assemble matrices | ||
| 364 |
2/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
|
51 | ForwardProjectionMatrix forwardM, forwardP; |
| 365 |
2/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
|
17 | forwardPatternM.exportIdx(forwardM); forwardM = 0.0; |
| 366 |
2/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
|
17 | forwardPatternP.exportIdx(forwardP); forwardP = 0.0; |
| 367 | |||
| 368 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
51 | BackwardProjectionMatrix backwardM, backwardP; |
| 369 | if (doBidirectional) | ||
| 370 | { | ||
| 371 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | backwardPatternM.exportIdx(backwardM); backwardM = 0.0; |
| 372 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | backwardPatternP.exportIdx(backwardP); backwardP = 0.0; |
| 373 | } | ||
| 374 | |||
| 375 |
4/4✓ Branch 0 taken 321334 times.
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 321334 times.
✓ Branch 3 taken 15 times.
|
335480 | for (const auto& is : intersections(glue)) |
| 376 | { | ||
| 377 |
1/2✓ Branch 1 taken 321312 times.
✗ Branch 2 not taken.
|
335446 | const auto& targetElement = is.targetEntity(0); |
| 378 |
1/2✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
|
335446 | const auto& targetElementGeometry = targetElement.geometry(); |
| 379 | |||
| 380 |
1/2✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
|
335446 | targetLocalView.bind( targetElement ); |
| 381 |
1/2✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
|
335446 | const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis(); |
| 382 | |||
| 383 | // perform integration over intersection geometry | ||
| 384 | using IsGeometry = typename std::decay_t<decltype(is.geometry())>; | ||
| 385 | using ctype = typename IsGeometry::ctype; | ||
| 386 | |||
| 387 |
1/2✓ Branch 1 taken 321334 times.
✗ Branch 2 not taken.
|
335446 | const auto& isGeometry = is.geometry(); |
| 388 | 335446 | const int intOrder = maxBasisOrder + 1; | |
| 389 | 335446 | const auto& quad = Dune::QuadratureRules<ctype, IsGeometry::mydimension>::rule(isGeometry.type(), intOrder); | |
| 390 |
4/4✓ Branch 0 taken 663836 times.
✓ Branch 1 taken 321334 times.
✓ Branch 2 taken 663836 times.
✓ Branch 3 taken 321334 times.
|
2432794 | for (auto&& qp : quad) |
| 391 | { | ||
| 392 | 713228 | const auto weight = qp.weight(); | |
| 393 | 713228 | const auto ie = isGeometry.integrationElement(qp.position()); | |
| 394 | 1426456 | const auto globalPos = isGeometry.global(qp.position()); | |
| 395 | |||
| 396 |
2/3✓ Branch 0 taken 663792 times.
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
|
2139684 | std::vector< Dune::FieldVector<ctype, 1> > targetShapeVals; |
| 397 |
2/5✓ Branch 1 taken 663792 times.
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
1426412 | targetLocalBasis.evaluateFunction(targetElementGeometry.local(globalPos), targetShapeVals); |
| 398 | |||
| 399 | // mass matrix entries target domain | ||
| 400 |
3/4✓ Branch 1 taken 2999456 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2335620 times.
✓ Branch 4 taken 663836 times.
|
3387536 | for (unsigned int i = 0; i < targetLocalBasis.size(); ++i) |
| 401 | { | ||
| 402 |
1/2✓ Branch 1 taken 2335620 times.
✗ Branch 2 not taken.
|
2674308 | const auto dofIdxI = targetLocalView.index(i); |
| 403 |
7/14✓ Branch 1 taken 2335620 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2335620 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2335620 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2335620 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2335620 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2335620 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2335620 times.
✗ Branch 20 not taken.
|
18720156 | forwardM[dofIdxI][dofIdxI] += ie*weight*targetShapeVals[i]*targetShapeVals[i]; |
| 404 | |||
| 405 |
3/4✓ Branch 1 taken 6243624 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2335620 times.
✓ Branch 4 taken 3908004 times.
|
7725384 | for (unsigned int j = i+1; j < targetLocalBasis.size(); ++j) |
| 406 | { | ||
| 407 |
1/2✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
|
5051076 | const auto dofIdxJ = targetLocalView.index(j); |
| 408 |
4/8✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3908004 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3908004 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3908004 times.
✗ Branch 11 not taken.
|
20204304 | const auto value = ie*weight*targetShapeVals[i]*targetShapeVals[j]; |
| 409 |
5/10✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3908004 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3908004 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3908004 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3908004 times.
✗ Branch 14 not taken.
|
20204304 | forwardM[dofIdxI][dofIdxJ] += value; |
| 410 |
4/8✓ Branch 1 taken 3908004 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3908004 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3908004 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3908004 times.
✗ Branch 11 not taken.
|
20204304 | forwardM[dofIdxJ][dofIdxI] += value; |
| 411 | } | ||
| 412 | } | ||
| 413 | |||
| 414 | // If targetDim < domainDim, there can be several "neighbors" if | ||
| 415 | // targetElement is aligned with a facet of domainElement. In | ||
| 416 | // this case make sure the basis functions are not added | ||
| 417 | // multiple times! (division by numNeighbors) | ||
| 418 | 713228 | const auto numNeighbors = is.numDomainNeighbors(); | |
| 419 |
2/2✓ Branch 0 taken 663880 times.
✓ Branch 1 taken 663836 times.
|
1426500 | for (unsigned int nIdx = 0; nIdx < numNeighbors; ++nIdx) |
| 420 | { | ||
| 421 |
1/2✓ Branch 1 taken 663880 times.
✗ Branch 2 not taken.
|
713272 | const auto& domainElement = is.domainEntity(nIdx); |
| 422 |
1/2✓ Branch 1 taken 663880 times.
✗ Branch 2 not taken.
|
713272 | domainLocalView.bind( domainElement ); |
| 423 |
1/2✓ Branch 1 taken 663880 times.
✗ Branch 2 not taken.
|
713272 | const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis(); |
| 424 | |||
| 425 |
1/2✓ Branch 1 taken 663880 times.
✗ Branch 2 not taken.
|
2139816 | std::vector< Dune::FieldVector<ctype, 1> > domainShapeVals; |
| 426 |
1/4✓ Branch 2 taken 663880 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
1426544 | domainLocalBasis.evaluateFunction(domainElement.geometry().local(globalPos), domainShapeVals); |
| 427 | |||
| 428 | // add entries in matrices | ||
| 429 |
3/4✓ Branch 1 taken 3000160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2336280 times.
✓ Branch 4 taken 663880 times.
|
3388240 | for (unsigned int i = 0; i < domainLocalBasis.size(); ++i) |
| 430 | { | ||
| 431 |
1/2✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
|
2674968 | const auto dofIdxDomain = domainLocalView.index(i); |
| 432 |
1/2✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
|
2674968 | const auto domainShapeVal = domainShapeVals[i]; |
| 433 | if (doBidirectional) | ||
| 434 | { | ||
| 435 |
6/12✓ Branch 1 taken 792 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 792 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 792 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 792 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 792 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 792 times.
✗ Branch 17 not taken.
|
4752 | backwardM[dofIdxDomain][dofIdxDomain] += ie*weight*domainShapeVal*domainShapeVal; |
| 436 | |||
| 437 |
3/4✓ Branch 1 taken 3960 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3168 times.
✓ Branch 4 taken 792 times.
|
3960 | for (unsigned int j = i+1; j < domainLocalBasis.size(); ++j) |
| 438 | { | ||
| 439 |
1/2✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
|
3168 | const auto dofIdxDomainJ = domainLocalView.index(j); |
| 440 |
3/6✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3168 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3168 times.
✗ Branch 8 not taken.
|
9504 | const auto value = ie*weight*domainShapeVal*domainShapeVals[j]; |
| 441 |
5/10✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3168 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3168 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3168 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3168 times.
✗ Branch 14 not taken.
|
12672 | backwardM[dofIdxDomain][dofIdxDomainJ] += value; |
| 442 |
4/8✓ Branch 1 taken 3168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3168 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3168 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3168 times.
✗ Branch 11 not taken.
|
12672 | backwardM[dofIdxDomainJ][dofIdxDomain] += value; |
| 443 | } | ||
| 444 | } | ||
| 445 | |||
| 446 |
3/4✓ Branch 1 taken 12489888 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10152024 times.
✓ Branch 4 taken 2337864 times.
|
15453408 | for (unsigned int j = 0; j < targetLocalBasis.size(); ++j) |
| 447 | { | ||
| 448 |
1/2✓ Branch 1 taken 10153608 times.
✗ Branch 2 not taken.
|
12778440 | const auto dofIdxTarget = targetLocalView.index(j); |
| 449 |
3/6✓ Branch 1 taken 10153608 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10153608 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10153608 times.
✗ Branch 8 not taken.
|
38335320 | const auto entry = ie*weight*domainShapeVal*targetShapeVals[j]; |
| 450 | |||
| 451 |
5/10✓ Branch 1 taken 10153608 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10153608 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10153608 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10153608 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2376 times.
✗ Branch 14 not taken.
|
51113760 | forwardP[dofIdxTarget][dofIdxDomain] += entry/numNeighbors; |
| 452 | if (doBidirectional) | ||
| 453 |
4/8✓ Branch 1 taken 2376 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2376 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2376 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2376 times.
✗ Branch 11 not taken.
|
9504 | backwardP[dofIdxDomain][dofIdxTarget] += entry; |
| 454 | } | ||
| 455 | } | ||
| 456 | } | ||
| 457 | } | ||
| 458 | } | ||
| 459 | |||
| 460 | // maybe treat zeroes on the diagonal | ||
| 461 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
|
17 | if (treatDiagonalZeroes) |
| 462 | { | ||
| 463 | ✗ | for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget) | |
| 464 | ✗ | if (forwardM[dofIdxTarget][dofIdxTarget] == 0.0) | |
| 465 | ✗ | forwardM[dofIdxTarget][dofIdxTarget] = 1.0; | |
| 466 | |||
| 467 | if (doBidirectional) | ||
| 468 | { | ||
| 469 | ✗ | for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain) | |
| 470 | ✗ | if (backwardM[dofIdxDomain][dofIdxDomain] == 0.0) | |
| 471 | ✗ | backwardM[dofIdxDomain][dofIdxDomain] = 1.0; | |
| 472 | } | ||
| 473 | } | ||
| 474 | |||
| 475 | 17 | return std::make_pair( std::make_pair(std::move(forwardM), std::move(forwardP)), | |
| 476 |
3/6✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15 times.
✗ Branch 8 not taken.
|
68 | std::make_pair(std::move(backwardM), std::move(backwardP)) ); |
| 477 | } | ||
| 478 | |||
| 479 | /*! | ||
| 480 | * \brief Creates a projector class between two function space bases | ||
| 481 | * \tparam doBidirectional If false, the backward projection matrix is not assembled | ||
| 482 | * \returns an std::pair with the forward and backward projector | ||
| 483 | */ | ||
| 484 | template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType> | ||
| 485 | 17 | auto makeProjectorPair(const FEBasisDomain& feBasisDomain, | |
| 486 | const FEBasisTarget& feBasisTarget, | ||
| 487 | const GlueType& glue) | ||
| 488 | { | ||
| 489 | using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector; | ||
| 490 | using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector; | ||
| 491 | |||
| 492 | using ForwardProjectionMatrix = typename ForwardProjector::Matrix; | ||
| 493 | using BackwardProjectionMatrix = typename BackwardProjector::Matrix; | ||
| 494 | |||
| 495 | 17 | auto projectionMatrices = createProjectionMatrices<doBidirectional>(feBasisDomain, feBasisTarget, glue, false); | |
| 496 | 17 | auto& forwardMatrices = projectionMatrices.first; | |
| 497 | 17 | auto& backwardMatrices = projectionMatrices.second; | |
| 498 | |||
| 499 | 17 | auto& forwardM = forwardMatrices.first; | |
| 500 | 17 | auto& forwardP = forwardMatrices.second; | |
| 501 | |||
| 502 | 17 | auto& backwardM = backwardMatrices.first; | |
| 503 | 17 | auto& backwardP = backwardMatrices.second; | |
| 504 | |||
| 505 | // determine the dofs that do not take part in intersections | ||
| 506 |
2/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
|
51 | std::vector<bool> isVoidTarget(forwardM.N(), false); |
| 507 |
2/2✓ Branch 0 taken 78355 times.
✓ Branch 1 taken 15 times.
|
78934 | for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget) |
| 508 |
4/8✓ Branch 1 taken 78355 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 78355 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 78355 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 78355 times.
|
157834 | if (forwardM[dofIdxTarget][dofIdxTarget] == 0.0) |
| 509 | ✗ | isVoidTarget[dofIdxTarget] = true; | |
| 510 | |||
| 511 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
34 | std::vector<bool> isVoidDomain; |
| 512 | if (doBidirectional) | ||
| 513 | { | ||
| 514 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | isVoidDomain.resize(backwardM.N(), false); |
| 515 |
2/2✓ Branch 0 taken 441 times.
✓ Branch 1 taken 1 times.
|
442 | for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain) |
| 516 |
6/8✓ Branch 1 taken 441 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 441 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 420 times.
✓ Branch 7 taken 21 times.
✓ Branch 8 taken 420 times.
✓ Branch 9 taken 21 times.
|
882 | if (backwardM[dofIdxDomain][dofIdxDomain] == 0.0) |
| 517 | 1260 | isVoidDomain[dofIdxDomain] = true; | |
| 518 | } | ||
| 519 | |||
| 520 | 51 | const bool hasVoidTarget = std::any_of(isVoidTarget.begin(), isVoidTarget.end(), [] (bool v) { return v; }); | |
| 521 | 51 | const bool hasVoidDomain = std::any_of(isVoidDomain.begin(), isVoidDomain.end(), [] (bool v) { return v; }); | |
| 522 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1 times.
|
17 | if (!hasVoidDomain && !hasVoidTarget) |
| 523 | { | ||
| 524 | 16 | return std::make_pair(ForwardProjector(std::move(forwardM), std::move(forwardP)), | |
| 525 |
2/4✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 14 times.
✗ Branch 5 not taken.
|
32 | BackwardProjector(std::move(backwardM), std::move(backwardP))); |
| 526 | } | ||
| 527 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | else if (!hasVoidDomain && hasVoidTarget) |
| 528 | { | ||
| 529 | ✗ | std::vector<std::size_t> expansionMapTarget; | |
| 530 | ✗ | ForwardProjectionMatrix forwardMReduced, forwardPReduced; | |
| 531 | ✗ | setupReducedMatrices(forwardM, forwardP, isVoidTarget, | |
| 532 | forwardMReduced, forwardPReduced, expansionMapTarget); | ||
| 533 | |||
| 534 | ✗ | return std::make_pair( ForwardProjector(std::move(forwardMReduced), | |
| 535 | ✗ | std::move(forwardPReduced), | |
| 536 | ✗ | std::move(expansionMapTarget), | |
| 537 | forwardM.N()), | ||
| 538 | ✗ | BackwardProjector(std::move(backwardM), std::move(backwardP)) ); | |
| 539 | } | ||
| 540 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | else if (hasVoidDomain && !hasVoidTarget) |
| 541 | { | ||
| 542 | if (doBidirectional) | ||
| 543 | { | ||
| 544 |
1/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
2 | std::vector<std::size_t> expansionMapDomain; |
| 545 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
|
3 | BackwardProjectionMatrix backwardMReduced, backwardPReduced; |
| 546 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | setupReducedMatrices(backwardM, backwardP, isVoidDomain, |
| 547 | backwardMReduced, backwardPReduced, expansionMapDomain); | ||
| 548 | |||
| 549 | 1 | return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)), | |
| 550 | 1 | BackwardProjector(std::move(backwardMReduced), | |
| 551 | 1 | std::move(backwardPReduced), | |
| 552 | 1 | std::move(expansionMapDomain), | |
| 553 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | backwardM.N()) ); |
| 554 | } | ||
| 555 | else | ||
| 556 | ✗ | return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)), | |
| 557 | ✗ | BackwardProjector(std::move(backwardM), std::move(backwardP)) ); | |
| 558 | } | ||
| 559 | else | ||
| 560 | { | ||
| 561 | ✗ | std::vector<std::size_t> expansionMapTarget; | |
| 562 | ✗ | ForwardProjectionMatrix forwardMReduced, forwardPReduced; | |
| 563 | ✗ | setupReducedMatrices(forwardM, forwardP, isVoidTarget, | |
| 564 | forwardMReduced, forwardPReduced, expansionMapTarget); | ||
| 565 | |||
| 566 | if (doBidirectional) | ||
| 567 | { | ||
| 568 | ✗ | std::vector<std::size_t> expansionMapDomain; | |
| 569 | ✗ | BackwardProjectionMatrix backwardMReduced, backwardPReduced; | |
| 570 | ✗ | setupReducedMatrices(backwardM, backwardP, isVoidDomain, | |
| 571 | backwardMReduced, backwardPReduced, expansionMapDomain); | ||
| 572 | |||
| 573 | ✗ | return std::make_pair( ForwardProjector(std::move(forwardMReduced), | |
| 574 | ✗ | std::move(forwardPReduced), | |
| 575 | ✗ | std::move(expansionMapTarget), | |
| 576 | forwardM.N()), | ||
| 577 | ✗ | BackwardProjector(std::move(backwardMReduced), | |
| 578 | ✗ | std::move(backwardPReduced), | |
| 579 | ✗ | std::move(expansionMapDomain), | |
| 580 | ✗ | backwardM.N()) ); | |
| 581 | } | ||
| 582 | else | ||
| 583 | ✗ | return std::make_pair( ForwardProjector(std::move(forwardMReduced), | |
| 584 | ✗ | std::move(forwardPReduced), | |
| 585 | ✗ | std::move(expansionMapTarget), | |
| 586 | forwardM.N()), | ||
| 587 | ✗ | BackwardProjector(std::move(backwardM), std::move(backwardP)) ); | |
| 588 | } | ||
| 589 | } | ||
| 590 | |||
| 591 | } // end namespace Detail | ||
| 592 | |||
| 593 | |||
| 594 | /*! | ||
| 595 | * \ingroup Discretization | ||
| 596 | * \brief Creates a pair of projectors between the space with | ||
| 597 | * basis feBasisDomain to the space with basis feBasisTarget. | ||
| 598 | * \param feBasisDomain The domain finite element space basis | ||
| 599 | * \param feBasisTarget The target finite element space basis | ||
| 600 | * \param glue The glue object containing the intersections between the grids. | ||
| 601 | * \return An std::pair of projectors where the first is the forward | ||
| 602 | * projector from the space with basis feBasisDomain to | ||
| 603 | * the space with basis feBasisTarget and the second | ||
| 604 | * does the backward projection. | ||
| 605 | */ | ||
| 606 | template< class FEBasisDomain, class FEBasisTarget, class GlueType > | ||
| 607 | auto makeProjectorPair(const FEBasisDomain& feBasisDomain, | ||
| 608 | const FEBasisTarget& feBasisTarget, | ||
| 609 | GlueType glue) | ||
| 610 | { | ||
| 611 | // we assume that target dim <= domain dimension | ||
| 612 | static constexpr int domainDim = FEBasisDomain::GridView::dimension; | ||
| 613 | static constexpr int targetDim = FEBasisTarget::GridView::dimension; | ||
| 614 | static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments"); | ||
| 615 | |||
| 616 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return Detail::makeProjectorPair<true>(feBasisDomain, feBasisTarget, glue); |
| 617 | } | ||
| 618 | |||
| 619 | /*! | ||
| 620 | * \ingroup Discretization | ||
| 621 | * \brief Creates a forward projector from the space feBasisDomain | ||
| 622 | * to the space with basis feBasisTarget. | ||
| 623 | * \param feBasisDomain The domain finite element space basis | ||
| 624 | * \param feBasisTarget The target finite element space basis | ||
| 625 | * \param glue The glue object containing the intersections between the grids. | ||
| 626 | * \return The forward projector from the space with basis feBasisDomain | ||
| 627 | * to the space with basis feBasisTarget. | ||
| 628 | */ | ||
| 629 | template< class FEBasisDomain, class FEBasisTarget, class GlueType > | ||
| 630 | 14 | auto makeProjector(const FEBasisDomain& feBasisDomain, | |
| 631 | const FEBasisTarget& feBasisTarget, | ||
| 632 | GlueType glue) | ||
| 633 | { | ||
| 634 | // we assume that target dim <= domain dimension | ||
| 635 | static constexpr int domainDim = FEBasisDomain::GridView::dimension; | ||
| 636 | static constexpr int targetDim = FEBasisTarget::GridView::dimension; | ||
| 637 | static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments"); | ||
| 638 | |||
| 639 | 14 | return Detail::makeProjectorPair<false>(feBasisDomain, feBasisTarget, glue).first; | |
| 640 | } | ||
| 641 | |||
| 642 | /*! | ||
| 643 | * \brief Creates the matrices underlying l2-projections | ||
| 644 | * \param feBasisDomain The basis to the domain finite element space | ||
| 645 | * \param feBasisTarget The basis to the target finite element space | ||
| 646 | * \param glue The glue object containing the intersections between the two grids | ||
| 647 | * \returns An std::pair of projection matrices, where the first entry stores the | ||
| 648 | * matrices of the forward projection and the second entry stores those | ||
| 649 | * of the backward projection. The entries of the returned pair are itself | ||
| 650 | * std::pairs which store the mass matrix in the first and the projection | ||
| 651 | * matrix in the second entry. | ||
| 652 | */ | ||
| 653 | template< class FEBasisDomain, class FEBasisTarget, class GlueType > | ||
| 654 | auto makeProjectionMatricesPair(const FEBasisDomain& feBasisDomain, | ||
| 655 | const FEBasisTarget& feBasisTarget, | ||
| 656 | GlueType glue) | ||
| 657 | { | ||
| 658 | // we assume that target dim <= domain dimension | ||
| 659 | static constexpr int domainDim = FEBasisDomain::GridView::dimension; | ||
| 660 | static constexpr int targetDim = FEBasisTarget::GridView::dimension; | ||
| 661 | static_assert(targetDim <= domainDim, "makeProjectionMatrixPair() expects targetDim < domainDim, please swap arguments"); | ||
| 662 | |||
| 663 | return Detail::createProjectionMatrices<true>(feBasisDomain, feBasisTarget, glue); | ||
| 664 | } | ||
| 665 | |||
| 666 | /*! | ||
| 667 | * \brief Creates the matrices underlying l2-projections | ||
| 668 | * \param feBasisDomain The basis to the domain finite element space | ||
| 669 | * \param feBasisTarget The basis to the target finite element space | ||
| 670 | * \param glue The glue object containing the intersections between the two grids | ||
| 671 | * \returns An std::pair of matrices, which store the mass matrix in the first and the | ||
| 672 | * projection matrix in the second entry. | ||
| 673 | */ | ||
| 674 | template< class FEBasisDomain, class FEBasisTarget, class GlueType > | ||
| 675 | auto makeProjectionMatrices(const FEBasisDomain& feBasisDomain, | ||
| 676 | const FEBasisTarget& feBasisTarget, | ||
| 677 | GlueType glue) | ||
| 678 | { | ||
| 679 | // we assume that target dim <= domain dimension | ||
| 680 | static constexpr int domainDim = FEBasisDomain::GridView::dimension; | ||
| 681 | static constexpr int targetDim = FEBasisTarget::GridView::dimension; | ||
| 682 | static_assert(targetDim <= domainDim, "makeProjectionMatrixPair() expects targetDim < domainDim, please swap arguments"); | ||
| 683 | |||
| 684 | return Detail::createProjectionMatrices<false>(feBasisDomain, feBasisTarget, glue).first; | ||
| 685 | } | ||
| 686 | |||
| 687 | } // end namespace Dumux | ||
| 688 | |||
| 689 | #endif | ||
| 690 |