| 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 EmbeddedCoupling | ||
| 10 | * \brief Coupling manager for low-dimensional domains embedded in the bulk | ||
| 11 | * domain. Point sources on each integration point are computed by an AABB tree. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGERBASE_HH | ||
| 15 | #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGERBASE_HH | ||
| 16 | |||
| 17 | #include <iostream> | ||
| 18 | #include <fstream> | ||
| 19 | #include <string> | ||
| 20 | #include <utility> | ||
| 21 | #include <unordered_map> | ||
| 22 | |||
| 23 | #include <dune/common/timer.hh> | ||
| 24 | #include <dune/geometry/quadraturerules.hh> | ||
| 25 | |||
| 26 | #include <dumux/io/format.hh> | ||
| 27 | #include <dumux/parallel/parallel_for.hh> | ||
| 28 | #include <dumux/common/properties.hh> | ||
| 29 | #include <dumux/common/numeqvector.hh> | ||
| 30 | #include <dumux/assembly/coloring.hh> | ||
| 31 | #include <dumux/geometry/distance.hh> | ||
| 32 | #include <dumux/geometry/intersectingentities.hh> | ||
| 33 | #include <dumux/discretization/method.hh> | ||
| 34 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 35 | #include <dumux/multidomain/glue.hh> | ||
| 36 | #include <dumux/multidomain/embedded/pointsourcedata.hh> | ||
| 37 | #include <dumux/multidomain/embedded/integrationpointsource.hh> | ||
| 38 | #include <dumux/multidomain/fvassembler.hh> | ||
| 39 | |||
| 40 | namespace Dumux { | ||
| 41 | |||
| 42 | //! the default point source traits | ||
| 43 | template<class MDTraits> | ||
| 44 | struct DefaultPointSourceTraits | ||
| 45 | { | ||
| 46 | private: | ||
| 47 | template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag; | ||
| 48 | template<std::size_t i> using GridGeometry = GetPropType<SubDomainTypeTag<i>, Properties::GridGeometry>; | ||
| 49 | template<std::size_t i> using NumEqVector = Dumux::NumEqVector<GetPropType<SubDomainTypeTag<i>, Properties::PrimaryVariables>>; | ||
| 50 | public: | ||
| 51 | //! export the point source type for domain i | ||
| 52 | template<std::size_t i> | ||
| 53 | using PointSource = IntegrationPointSource<typename GridGeometry<i>::GlobalCoordinate, NumEqVector<i>>; | ||
| 54 | |||
| 55 | //! export the point source helper type for domain i | ||
| 56 | template<std::size_t i> | ||
| 57 | using PointSourceHelper = IntegrationPointSourceHelper; | ||
| 58 | |||
| 59 | //! export the point source data type | ||
| 60 | using PointSourceData = Dumux::PointSourceData<MDTraits>; | ||
| 61 | }; | ||
| 62 | |||
| 63 | /*! | ||
| 64 | * \ingroup EmbeddedCoupling | ||
| 65 | * \brief Manages the coupling between bulk elements and lower dimensional elements | ||
| 66 | * Point sources on each integration point are computed by an AABB tree. | ||
| 67 | */ | ||
| 68 | template<class MDTraits, class Implementation, class PSTraits = DefaultPointSourceTraits<MDTraits>> | ||
| 69 | class EmbeddedCouplingManagerBase | ||
| 70 | : public CouplingManager<MDTraits> | ||
| 71 | { | ||
| 72 | using ParentType = CouplingManager<MDTraits>; | ||
| 73 | using Scalar = typename MDTraits::Scalar; | ||
| 74 | static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index(); | ||
| 75 | static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index(); | ||
| 76 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 77 | using PointSourceData = typename PSTraits::PointSourceData; | ||
| 78 | |||
| 79 | // the sub domain type tags | ||
| 80 | template<std::size_t id> using PointSource = typename PSTraits::template PointSource<id>; | ||
| 81 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 82 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 83 | template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>; | ||
| 84 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 85 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 86 | template<std::size_t id> using ElementMapper = typename GridGeometry<id>::ElementMapper; | ||
| 87 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 88 | template<std::size_t id> using ElementSeed = typename GridView<id>::Grid::template Codim<0>::EntitySeed; | ||
| 89 | template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex; | ||
| 90 | template<std::size_t id> using CouplingStencil = std::vector<GridIndex<id>>; | ||
| 91 | |||
| 92 | static constexpr int bulkDim = GridView<bulkIdx>::dimension; | ||
| 93 | static constexpr int lowDimDim = GridView<lowDimIdx>::dimension; | ||
| 94 | static constexpr int dimWorld = GridView<bulkIdx>::dimensionworld; | ||
| 95 | |||
| 96 | template<std::size_t id> | ||
| 97 | static constexpr bool isBox() | ||
| 98 | { return GridGeometry<id>::discMethod == DiscretizationMethods::box; } | ||
| 99 | |||
| 100 | using GlobalPosition = typename Element<bulkIdx>::Geometry::GlobalCoordinate; | ||
| 101 | using GlueType = MultiDomainGlue<GridView<bulkIdx>, GridView<lowDimIdx>, ElementMapper<bulkIdx>, ElementMapper<lowDimIdx>>; | ||
| 102 | |||
| 103 | public: | ||
| 104 | //! export traits | ||
| 105 | using MultiDomainTraits = MDTraits; | ||
| 106 | //! export the point source traits | ||
| 107 | using PointSourceTraits = PSTraits; | ||
| 108 | //! export stencil types | ||
| 109 | template<std::size_t id> using CouplingStencils = std::unordered_map<GridIndex<id>, CouplingStencil<id>>; | ||
| 110 | |||
| 111 | /*! | ||
| 112 | * \brief call this after grid adaption | ||
| 113 | */ | ||
| 114 | ✗ | void updateAfterGridAdaption(std::shared_ptr<const GridGeometry<bulkIdx>> bulkGridGeometry, | |
| 115 | std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimGridGeometry) | ||
| 116 | { | ||
| 117 | ✗ | glue_ = std::make_shared<GlueType>(); | |
| 118 | ✗ | } | |
| 119 | |||
| 120 | /*! | ||
| 121 | * \brief Constructor | ||
| 122 | */ | ||
| 123 | 28 | EmbeddedCouplingManagerBase(std::shared_ptr<const GridGeometry<bulkIdx>> bulkGridGeometry, | |
| 124 | std::shared_ptr<const GridGeometry<lowDimIdx>> lowDimGridGeometry) | ||
| 125 | 28 | { | |
| 126 |
3/10✓ Branch 3 taken 28 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 28 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
56 | updateAfterGridAdaption(bulkGridGeometry, lowDimGridGeometry); |
| 127 | 28 | } | |
| 128 | |||
| 129 | /*! | ||
| 130 | * \brief Methods to be accessed by main | ||
| 131 | */ | ||
| 132 | // \{ | ||
| 133 | |||
| 134 | 28 | void init(std::shared_ptr<Problem<bulkIdx>> bulkProblem, | |
| 135 | std::shared_ptr<Problem<lowDimIdx>> lowDimProblem, | ||
| 136 | const SolutionVector& curSol) | ||
| 137 | { | ||
| 138 | 28 | this->updateSolution(curSol); | |
| 139 | 56 | this->setSubProblems(std::make_tuple(bulkProblem, lowDimProblem)); | |
| 140 | |||
| 141 | 28 | integrationOrder_ = getParam<int>("MixedDimension.IntegrationOrder", 1); | |
| 142 | 28 | asImp_().computePointSourceData(integrationOrder_); | |
| 143 | 28 | } | |
| 144 | |||
| 145 | // \} | ||
| 146 | |||
| 147 | /*! | ||
| 148 | * \brief Methods to be accessed by the assembly | ||
| 149 | */ | ||
| 150 | // \{ | ||
| 151 | |||
| 152 | /*! | ||
| 153 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 154 | * that couple with / influence the element residual of the given element of domain i | ||
| 155 | * | ||
| 156 | * \param domainI the domain index of domain i | ||
| 157 | * \param element the coupled element of domain à | ||
| 158 | * \param domainJ the domain index of domain j | ||
| 159 | * | ||
| 160 | * \note The element residual definition depends on the discretization scheme of domain i | ||
| 161 | * box: a container of the residuals of all sub control volumes | ||
| 162 | * cc : the residual of the (sub) control volume | ||
| 163 | * fem: the residual of the element | ||
| 164 | * \note This function has to be implemented by all coupling managers for all combinations of i and j | ||
| 165 | */ | ||
| 166 | template<std::size_t i, std::size_t j> | ||
| 167 | ✗ | const CouplingStencil<j>& couplingStencil(Dune::index_constant<i> domainI, | |
| 168 | const Element<i>& element, | ||
| 169 | Dune::index_constant<j> domainJ) const | ||
| 170 | { | ||
| 171 | static_assert(i != j, "A domain cannot be coupled to itself!"); | ||
| 172 | |||
| 173 | ✗ | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); | |
| 174 | ✗ | if (couplingStencils(domainI).count(eIdx)) | |
| 175 | ✗ | return couplingStencils(domainI).at(eIdx); | |
| 176 | else | ||
| 177 | ✗ | return emptyStencil(domainI); | |
| 178 | } | ||
| 179 | |||
| 180 | /*! | ||
| 181 | * \brief evaluates the element residual of a coupled element of domain i which depends on the variables | ||
| 182 | * at the degree of freedom with index dofIdxGlobalJ of domain j | ||
| 183 | * | ||
| 184 | * \param domainI the domain index of domain i | ||
| 185 | * \param localAssemblerI the local assembler assembling the element residual of an element of domain i | ||
| 186 | * \param domainJ the domain index of domain j | ||
| 187 | * \param dofIdxGlobalJ the index of the degree of freedom of domain j which has an influence on the element residual of domain i | ||
| 188 | * | ||
| 189 | * \note we only need to evaluate the source contribution to the residual here as the coupling term is the source | ||
| 190 | * \return the element residual | ||
| 191 | */ | ||
| 192 | template<std::size_t i, std::size_t j, class LocalAssemblerI> | ||
| 193 | 29820292 | decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI, | |
| 194 | const LocalAssemblerI& localAssemblerI, | ||
| 195 | Dune::index_constant<j> domainJ, | ||
| 196 | std::size_t dofIdxGlobalJ) | ||
| 197 | { | ||
| 198 | static_assert(i != j, "A domain cannot be coupled to itself!"); | ||
| 199 | |||
| 200 | 29820292 | typename LocalAssemblerI::LocalResidual::ElementResidualVector residual; | |
| 201 | |||
| 202 | 29820292 | const auto& element = localAssemblerI.element(); | |
| 203 | 29820292 | const auto& fvGeometry = localAssemblerI.fvGeometry(); | |
| 204 | 29820292 | const auto& curElemVolVars = localAssemblerI.curElemVolVars(); | |
| 205 | |||
| 206 | 59640584 | residual.resize(fvGeometry.numScv()); | |
| 207 |
5/6✓ Branch 0 taken 16257444 times.
✓ Branch 1 taken 14910146 times.
✓ Branch 2 taken 16257444 times.
✓ Branch 3 taken 14910146 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 14613108 times.
|
124670360 | for (const auto& scv : scvs(fvGeometry)) |
| 208 | { | ||
| 209 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 16257444 times.
✓ Branch 3 taken 9487872 times.
✗ Branch 4 not taken.
|
32514888 | auto couplingSource = this->problem(domainI).scvPointSources(element, fvGeometry, curElemVolVars, scv); |
| 210 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 16257444 times.
✓ Branch 3 taken 9487872 times.
✗ Branch 4 not taken.
|
33583644 | couplingSource += this->problem(domainI).source(element, fvGeometry, curElemVolVars, scv); |
| 211 | 68318448 | couplingSource *= -GridGeometry<i>::Extrusion::volume(fvGeometry, scv)*curElemVolVars[scv].extrusionFactor(); | |
| 212 | 65029776 | residual[scv.indexInElement()] = couplingSource; | |
| 213 | } | ||
| 214 | 29820292 | return residual; | |
| 215 | } | ||
| 216 | |||
| 217 | // \} | ||
| 218 | |||
| 219 | /* \brief Compute integration point point sources and associated data | ||
| 220 | * | ||
| 221 | * This method uses grid glue to intersect the given grids. Over each intersection | ||
| 222 | * we later need to integrate a source term. This method places point sources | ||
| 223 | * at each quadrature point and provides the point source with the necessary | ||
| 224 | * information to compute integrals (quadrature weight and integration element) | ||
| 225 | * \param order The order of the quadrature rule for integration of sources over an intersection | ||
| 226 | * \param verbose If the point source computation is verbose | ||
| 227 | */ | ||
| 228 | 3 | void computePointSourceData(std::size_t order = 1, bool verbose = false) | |
| 229 | { | ||
| 230 | // initialize the maps | ||
| 231 | // do some logging and profiling | ||
| 232 | 3 | Dune::Timer watch; | |
| 233 | 6 | std::cout << "Initializing the point sources..." << std::endl; | |
| 234 | |||
| 235 | // clear all internal members like pointsource vectors and stencils | ||
| 236 | // initializes the point source id counter | ||
| 237 | 3 | clear(); | |
| 238 | |||
| 239 | // precompute the vertex indices for efficiency for the box method | ||
| 240 | 3 | this->precomputeVertexIndices(bulkIdx); | |
| 241 | 3 | this->precomputeVertexIndices(lowDimIdx); | |
| 242 | |||
| 243 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
|
3 | const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); |
| 244 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); |
| 245 | |||
| 246 | // intersect the bounding box trees | ||
| 247 | 3 | glueGrids(); | |
| 248 | |||
| 249 | 3 | pointSourceData_.reserve(this->glue().size()); | |
| 250 | 3 | averageDistanceToBulkCell_.reserve(this->glue().size()); | |
| 251 |
4/4✓ Branch 0 taken 22608 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 22608 times.
✓ Branch 3 taken 3 times.
|
22611 | for (const auto& is : intersections(this->glue())) |
| 252 | { | ||
| 253 | // all inside elements are identical... | ||
| 254 | 22608 | const auto& inside = is.targetEntity(0); | |
| 255 | // get the intersection geometry for integrating over it | ||
| 256 | 22608 | const auto intersectionGeometry = is.geometry(); | |
| 257 | |||
| 258 | // get the Gaussian quadrature rule for the local intersection | ||
| 259 | 22608 | const auto& quad = Dune::QuadratureRules<Scalar, lowDimDim>::rule(intersectionGeometry.type(), order); | |
| 260 | 45216 | const std::size_t lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside); | |
| 261 | |||
| 262 | // iterate over all quadrature points | ||
| 263 |
4/4✓ Branch 0 taken 22628 times.
✓ Branch 1 taken 22608 times.
✓ Branch 2 taken 22628 times.
✓ Branch 3 taken 22608 times.
|
90452 | for (auto&& qp : quad) |
| 264 | { | ||
| 265 | // compute the coupling stencils | ||
| 266 |
4/4✓ Branch 0 taken 22748 times.
✓ Branch 1 taken 22628 times.
✓ Branch 2 taken 22748 times.
✓ Branch 3 taken 22628 times.
|
68004 | for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) |
| 267 | { | ||
| 268 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | const auto& outside = is.domainEntity(outsideIdx); |
| 269 | 45496 | const std::size_t bulkElementIdx = bulkGridGeometry.elementMapper().index(outside); | |
| 270 | |||
| 271 | // each quadrature point will be a point source for the sub problem | ||
| 272 | 45496 | const auto globalPos = intersectionGeometry.global(qp.position()); | |
| 273 | 22748 | const auto id = idCounter_++; | |
| 274 | 22748 | const auto qpweight = qp.weight(); | |
| 275 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | const auto ie = intersectionGeometry.integrationElement(qp.position()); |
| 276 | 45496 | pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, ie, bulkElementIdx); | |
| 277 |
2/4✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22748 times.
✗ Branch 5 not taken.
|
45496 | pointSources(bulkIdx).back().setEmbeddings(is.numDomainNeighbors()); |
| 278 | 45496 | pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, lowDimElementIdx); | |
| 279 |
2/4✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22748 times.
✗ Branch 5 not taken.
|
45496 | pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); |
| 280 | |||
| 281 | // pre compute additional data used for the evaluation of | ||
| 282 | // the actual solution dependent source term | ||
| 283 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | PointSourceData psData; |
| 284 | |||
| 285 | if constexpr (isBox<lowDimIdx>()) | ||
| 286 | { | ||
| 287 | using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >; | ||
| 288 | const auto lowDimGeometry = this->problem(lowDimIdx).gridGeometry().element(lowDimElementIdx).geometry(); | ||
| 289 | ShapeValues shapeValues; | ||
| 290 | this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).gridGeometry(), lowDimGeometry, globalPos, shapeValues); | ||
| 291 | psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); | ||
| 292 | } | ||
| 293 | else | ||
| 294 | { | ||
| 295 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | psData.addLowDimInterpolation(lowDimElementIdx); |
| 296 | } | ||
| 297 | |||
| 298 | // add data needed to compute integral over the circle | ||
| 299 | if constexpr (isBox<bulkIdx>()) | ||
| 300 | { | ||
| 301 | using ShapeValues = std::vector<Dune::FieldVector<Scalar, 1> >; | ||
| 302 | const auto bulkGeometry = this->problem(bulkIdx).gridGeometry().element(bulkElementIdx).geometry(); | ||
| 303 | ShapeValues shapeValues; | ||
| 304 | this->getShapeValues(bulkIdx, this->problem(bulkIdx).gridGeometry(), bulkGeometry, globalPos, shapeValues); | ||
| 305 | psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); | ||
| 306 | } | ||
| 307 | else | ||
| 308 | { | ||
| 309 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | psData.addBulkInterpolation(bulkElementIdx); |
| 310 | } | ||
| 311 | |||
| 312 | // publish point source data in the global vector | ||
| 313 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | this->pointSourceData().emplace_back(std::move(psData)); |
| 314 | |||
| 315 | // compute average distance to bulk cell | ||
| 316 |
2/4✓ Branch 3 taken 22748 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22748 times.
✗ Branch 7 not taken.
|
45496 | averageDistanceToBulkCell_.push_back(averageDistancePointGeometry(globalPos, outside.geometry())); |
| 317 | |||
| 318 | // export the lowdim coupling stencil | ||
| 319 | // we insert all vertices / elements and make it unique later | ||
| 320 | if (isBox<bulkIdx>()) | ||
| 321 | { | ||
| 322 | const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); | ||
| 323 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), | ||
| 324 | vertices.begin(), vertices.end()); | ||
| 325 | } | ||
| 326 | else | ||
| 327 | { | ||
| 328 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); |
| 329 | } | ||
| 330 | |||
| 331 | // export the bulk coupling stencil | ||
| 332 | // we insert all vertices / elements and make it unique later | ||
| 333 | if (isBox<lowDimIdx>()) | ||
| 334 | { | ||
| 335 | const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx); | ||
| 336 | this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), | ||
| 337 | vertices.begin(), vertices.end()); | ||
| 338 | |||
| 339 | } | ||
| 340 | else | ||
| 341 | { | ||
| 342 |
1/2✓ Branch 1 taken 22748 times.
✗ Branch 2 not taken.
|
22748 | this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); |
| 343 | } | ||
| 344 | } | ||
| 345 | } | ||
| 346 | } | ||
| 347 | |||
| 348 | // make stencils unique | ||
| 349 | using namespace Dune::Hybrid; | ||
| 350 | 27 | forEach(integralRange(std::integral_constant<std::size_t, MDTraits::numSubDomains>{}), [&](const auto domainIdx) | |
| 351 | { | ||
| 352 |
4/4✓ Branch 0 taken 616 times.
✓ Branch 1 taken 3 times.
✓ Branch 3 taken 3016 times.
✓ Branch 4 taken 3 times.
|
3644 | for (auto&& stencil : this->couplingStencils(domainIdx)) |
| 353 | { | ||
| 354 | 10896 | std::sort(stencil.second.begin(), stencil.second.end()); | |
| 355 | 18160 | stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); | |
| 356 | } | ||
| 357 | }); | ||
| 358 | |||
| 359 | 9 | std::cout << "took " << watch.elapsed() << " seconds." << std::endl; | |
| 360 | 3 | } | |
| 361 | |||
| 362 | /*! | ||
| 363 | * \brief Methods to be accessed by the subproblems | ||
| 364 | */ | ||
| 365 | // \{ | ||
| 366 | |||
| 367 | //! Return a reference to the pointSource data | ||
| 368 | const PointSourceData& pointSourceData(std::size_t id) const | ||
| 369 |
3/4✓ Branch 0 taken 670548 times.
✓ Branch 1 taken 322387 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3175871 times.
|
620421217 | { return pointSourceData_[id]; } |
| 370 | |||
| 371 | //! Return a reference to the bulk problem | ||
| 372 | template<std::size_t id> | ||
| 373 | ✗ | const GridView<id>& gridView(Dune::index_constant<id> domainIdx) const | |
| 374 |
9/17✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 20 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✓ Branch 16 taken 7 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 7 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 7 times.
|
118 | { return this->problem(domainIdx).gridGeometry().gridView(); } |
| 375 | |||
| 376 | //! Return data for a bulk point source with the identifier id | ||
| 377 | PrimaryVariables<bulkIdx> bulkPriVars(std::size_t id) const | ||
| 378 |
8/8✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1018734 times.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 1152069 times.
✓ Branch 4 taken 7 times.
✓ Branch 5 taken 464201 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 324926 times.
|
627015231 | { return pointSourceData_[id].interpolateBulk(this->curSol(bulkIdx)); } |
| 379 | |||
| 380 | //! Return data for a low dim point source with the identifier id | ||
| 381 | PrimaryVariables<lowDimIdx> lowDimPriVars(std::size_t id) const | ||
| 382 |
6/6✓ Branch 0 taken 670552 times.
✓ Branch 1 taken 7375777 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 10769177 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 1367059 times.
|
327887574 | { return pointSourceData_[id].interpolateLowDim(this->curSol(lowDimIdx)); } |
| 383 | |||
| 384 | //! return the average distance to the coupled bulk cell center | ||
| 385 | Scalar averageDistance(std::size_t id) const | ||
| 386 |
8/8✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1012794 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 1012794 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 324926 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 324926 times.
|
2675452 | { return averageDistanceToBulkCell_[id]; } |
| 387 | |||
| 388 | //! Return reference to bulk point sources | ||
| 389 | const std::vector<PointSource<bulkIdx>>& bulkPointSources() const | ||
| 390 |
2/4✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
|
52 | { return std::get<bulkIdx>(pointSources_); } |
| 391 | |||
| 392 | //! Return reference to low dim point sources | ||
| 393 | const std::vector<PointSource<lowDimIdx>>& lowDimPointSources() const | ||
| 394 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
56 | { return std::get<lowDimIdx>(pointSources_); } |
| 395 | |||
| 396 | //! Return the point source if domain i | ||
| 397 | template<std::size_t i> | ||
| 398 | const std::vector<PointSource<i>>& pointSources(Dune::index_constant<i> dom) const | ||
| 399 | { return std::get<i>(pointSources_); } | ||
| 400 | |||
| 401 | //! Return reference to bulk coupling stencil member of domain i | ||
| 402 | template<std::size_t i> | ||
| 403 | ✗ | const CouplingStencils<i>& couplingStencils(Dune::index_constant<i> dom) const | |
| 404 |
0/8✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
4680 | { return std::get<i>(couplingStencils_); } |
| 405 | |||
| 406 | //! Return reference to point source data vector member | ||
| 407 | const std::vector<PointSourceData>& pointSourceData() const | ||
| 408 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 295099408 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 300816742 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1405896 times.
|
597322046 | { return pointSourceData_; } |
| 409 | |||
| 410 | //! Return a reference to an empty stencil | ||
| 411 | template<std::size_t i> | ||
| 412 | ✗ | const CouplingStencil<i>& emptyStencil(Dune::index_constant<i> dom) const | |
| 413 |
4/8✓ Branch 0 taken 1506340 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1506340 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 9594 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 9594 times.
|
8445344 | { return std::get<i>(emptyStencil_); } |
| 414 | |||
| 415 | /*! | ||
| 416 | * \brief the solution vector of the subproblem | ||
| 417 | * \param domainIdx The domain index | ||
| 418 | * \note in case of numeric differentiation the solution vector always carries the deflected solution | ||
| 419 | */ | ||
| 420 | template<std::size_t i> | ||
| 421 | ✗ | auto& curSol(Dune::index_constant<i> domainIdx) | |
| 422 |
3/6✓ Branch 0 taken 1457656 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1506340 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 48684 times.
✗ Branch 5 not taken.
|
22980836 | { return ParentType::curSol(domainIdx); } |
| 423 | |||
| 424 | /*! | ||
| 425 | * \brief the solution vector of the subproblem | ||
| 426 | * \param domainIdx The domain index | ||
| 427 | * \note in case of numeric differentiation the solution vector always carries the deflected solution | ||
| 428 | */ | ||
| 429 | template<std::size_t i> | ||
| 430 | ✗ | const auto& curSol(Dune::index_constant<i> domainIdx) const | |
| 431 |
16/16✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1018734 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 1018734 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 7375777 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 8603561 times.
✓ Branch 8 taken 7 times.
✓ Branch 9 taken 1691985 times.
✓ Branch 10 taken 7 times.
✓ Branch 11 taken 464201 times.
✓ Branch 12 taken 7 times.
✓ Branch 13 taken 464201 times.
✓ Branch 14 taken 7 times.
✓ Branch 15 taken 464201 times.
|
1262293270 | { return ParentType::curSol(domainIdx); } |
| 432 | |||
| 433 | /*! | ||
| 434 | * \brief Compute colors for multithreaded assembly | ||
| 435 | */ | ||
| 436 | 26 | void computeColorsForAssembly() | |
| 437 | { | ||
| 438 | 26 | Dune::Timer timer; | |
| 439 | |||
| 440 | // start with elements connected within the subdomains | ||
| 441 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
|
26 | const auto& bulkGG = asImp_().problem(bulkIdx).gridGeometry(); |
| 442 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
|
26 | const auto& lowDimGG = asImp_().problem(lowDimIdx).gridGeometry(); |
| 443 | |||
| 444 | 52 | auto connectedElementsBulk = Detail::computeConnectedElements(bulkGG); | |
| 445 |
1/2✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
|
26 | auto connectedElementsLowDim = Detail::computeConnectedElements(lowDimGG); |
| 446 | |||
| 447 | // to the elements from the own domain, we have to add the extended source stencil (non-local couplings) | ||
| 448 | // we consider this only in the bulk domain | ||
| 449 |
14/18✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 953896 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 3198 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 3198 times.
✓ Branch 12 taken 8224 times.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 11422 times.
✓ Branch 15 taken 1 times.
✓ Branch 17 taken 8224 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 8224 times.
✗ Branch 21 not taken.
|
1912822 | for (const auto& element : elements(bulkGG.gridView())) |
| 450 | { | ||
| 451 |
2/4✓ Branch 1 taken 11422 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11422 times.
✗ Branch 5 not taken.
|
1924192 | const auto eIdx = bulkGG.elementMapper().index(element); |
| 452 | 986320 | const auto& extendedSourceStencil = asImp_().extendedSourceStencil(eIdx); | |
| 453 |
4/4✓ Branch 0 taken 72633 times.
✓ Branch 1 taken 962096 times.
✓ Branch 2 taken 72633 times.
✓ Branch 3 taken 962096 times.
|
2934697 | for (auto dofIdx : extendedSourceStencil) |
| 454 | { | ||
| 455 | 72633 | const auto& elems = connectedElementsBulk[dofIdx]; | |
| 456 |
4/4✓ Branch 4 taken 64364 times.
✓ Branch 5 taken 8269 times.
✓ Branch 6 taken 64364 times.
✓ Branch 7 taken 8269 times.
|
290532 | if (std::find(elems.begin(), elems.end(), eIdx) == elems.end()) |
| 457 |
2/4✓ Branch 1 taken 64364 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64364 times.
✗ Branch 5 not taken.
|
128728 | connectedElementsBulk[dofIdx].push_back(eIdx); |
| 458 | } | ||
| 459 | } | ||
| 460 | |||
| 461 | // then we also need a similar container that tells us | ||
| 462 | // which elements _from the other domain_ modify our dofs. | ||
| 463 | // That, we can obtain via the coupling stencils | ||
| 464 | 52 | std::array<std::vector<std::vector<std::size_t>>, 2> connectedElementsCoupling; | |
| 465 | |||
| 466 | using namespace Dune::Hybrid; | ||
| 467 |
1/2✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
|
546 | forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto i) |
| 468 | { | ||
| 469 |
4/16✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
104 | connectedElementsCoupling[i()].resize(asImp_().problem(i).gridGeometry().numDofs()); |
| 470 |
0/4✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4168792 | forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto j) |
| 471 | { | ||
| 472 | if constexpr (i != j) | ||
| 473 | { | ||
| 474 |
24/28✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✓ Branch 5 taken 953896 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 3198 times.
✓ Branch 8 taken 950675 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 3222 times.
✓ Branch 11 taken 8224 times.
✓ Branch 12 taken 1 times.
✓ Branch 13 taken 8225 times.
✓ Branch 14 taken 6173 times.
✓ Branch 15 taken 24 times.
✓ Branch 16 taken 14396 times.
✓ Branch 17 taken 1762 times.
✓ Branch 18 taken 1 times.
✓ Branch 19 taken 16134 times.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 1738 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
✓ Branch 31 taken 44 times.
✓ Branch 32 taken 1 times.
✓ Branch 33 taken 44 times.
✓ Branch 34 taken 1 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 44 times.
|
1937251 | for (const auto& element : elements(asImp_().problem(j).gridGeometry().gridView())) |
| 475 | { | ||
| 476 |
5/11✗ Branch 0 not taken.
✓ Branch 1 taken 962096 times.
✓ Branch 3 taken 8224 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 16134 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8224 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 44 times.
|
970050 | const auto eIdx = asImp_().problem(j).gridGeometry().elementMapper().index(element); |
| 477 |
11/11✓ Branch 1 taken 30559 times.
✓ Branch 2 taken 953872 times.
✓ Branch 3 taken 24503 times.
✓ Branch 4 taken 962096 times.
✓ Branch 5 taken 2168 times.
✓ Branch 6 taken 53343 times.
✓ Branch 7 taken 7910 times.
✓ Branch 8 taken 47580 times.
✓ Branch 9 taken 7954 times.
✓ Branch 10 taken 2461 times.
✓ Branch 11 taken 44 times.
|
1042133 | for (auto dofIdx : asImp_().couplingStencil(j, element, i)) |
| 478 | { | ||
| 479 | 144166 | const auto& elems = connectedElementsCoupling[i()][dofIdx]; | |
| 480 |
4/8✓ Branch 4 taken 24503 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 24503 times.
✗ Branch 7 not taken.
✓ Branch 12 taken 47580 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 47580 times.
✗ Branch 15 not taken.
|
288332 | if (std::find(elems.begin(), elems.end(), eIdx) == elems.end()) |
| 481 |
3/6✓ Branch 1 taken 2168 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2168 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2168 times.
✗ Branch 8 not taken.
|
216249 | connectedElementsCoupling[i()][dofIdx].push_back(eIdx); |
| 482 | } | ||
| 483 | } | ||
| 484 | } | ||
| 485 | }); | ||
| 486 | }); | ||
| 487 | |||
| 488 |
1/2✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
|
5804730 | forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto i) |
| 489 | { | ||
| 490 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 26 times.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
|
52 | const auto& gg = asImp_().problem(i).gridGeometry(); |
| 491 | |||
| 492 |
7/17✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 10 taken 25 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 25 times.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 25 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 25 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 25 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
|
208 | std::vector<int> colors(gg.gridView().size(0), -1); |
| 493 | |||
| 494 | // pre-reserve some memory for helper arrays to avoid reallocation | ||
| 495 |
6/16✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 11 taken 26 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 26 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 26 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
156 | std::vector<int> neighborColors; neighborColors.reserve(200); |
| 496 |
7/13✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 25 times.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
|
156 | std::vector<bool> colorUsed; colorUsed.reserve(200); |
| 497 | |||
| 498 |
1/2✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
|
52 | auto& elementSets = std::get<i>(elementSets_); |
| 499 |
2/4✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
|
104 | const auto& connectedElements = std::get<i>(std::tie(connectedElementsBulk, connectedElementsLowDim)); |
| 500 | |||
| 501 |
20/24✓ Branch 2 taken 7954 times.
✓ Branch 3 taken 26 times.
✓ Branch 4 taken 7954 times.
✓ Branch 5 taken 26 times.
✓ Branch 6 taken 7928 times.
✓ Branch 7 taken 26 times.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 953896 times.
✓ Branch 16 taken 2 times.
✓ Branch 17 taken 3198 times.
✓ Branch 18 taken 950651 times.
✓ Branch 19 taken 3222 times.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 8224 times.
✓ Branch 22 taken 3199 times.
✓ Branch 23 taken 8224 times.
✓ Branch 24 taken 1 times.
✓ Branch 26 taken 8224 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 8224 times.
✗ Branch 30 not taken.
|
1937303 | for (const auto& element : elements(gg.gridView())) |
| 502 | { | ||
| 503 | // compute neighbor colors based on discretization-dependent stencil | ||
| 504 |
4/4✓ Branch 0 taken 7928 times.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 962070 times.
✓ Branch 3 taken 26 times.
|
970050 | neighborColors.clear(); |
| 505 |
2/4✓ Branch 1 taken 7954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 962096 times.
✗ Branch 5 not taken.
|
970050 | Detail::addNeighborColors(gg, element, colors, connectedElements, neighborColors); |
| 506 | |||
| 507 | 962096 | if constexpr (i == bulkIdx) | |
| 508 | { | ||
| 509 | // make sure we also add all elements in the extended source stencil | ||
| 510 | // Detail::addNeighborColors only adds direct neighbors | ||
| 511 | // add everyone that also modifies the extended source stencil dofs | ||
| 512 |
2/4✓ Branch 1 taken 11422 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11422 times.
✗ Branch 5 not taken.
|
1924192 | const auto eIdx = bulkGG.elementMapper().index(element); |
| 513 | 986320 | const auto& extendedSourceStencil = asImp_().extendedSourceStencil(eIdx); | |
| 514 |
4/4✓ Branch 0 taken 72633 times.
✓ Branch 1 taken 962096 times.
✓ Branch 2 taken 72633 times.
✓ Branch 3 taken 962096 times.
|
2934697 | for (auto dofIdx : extendedSourceStencil) |
| 515 |
4/4✓ Branch 0 taken 2023554 times.
✓ Branch 1 taken 72633 times.
✓ Branch 2 taken 2023554 times.
✓ Branch 3 taken 72633 times.
|
2314086 | for (const auto nIdx : connectedElementsBulk[dofIdx]) |
| 516 |
2/4✓ Branch 1 taken 2023554 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2023554 times.
✗ Branch 5 not taken.
|
4047108 | neighborColors.push_back(colors[nIdx]); |
| 517 | } | ||
| 518 | |||
| 519 | // add neighbor colors based on coupling stencil | ||
| 520 |
2/4✓ Branch 1 taken 7954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 962096 times.
✗ Branch 5 not taken.
|
8059694 | forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto j) |
| 521 | { | ||
| 522 | 970050 | if constexpr (i != j) | |
| 523 | { | ||
| 524 | // add all elements that manipulate the same coupled dofs | ||
| 525 |
8/8✓ Branch 1 taken 47580 times.
✓ Branch 2 taken 7954 times.
✓ Branch 3 taken 47580 times.
✓ Branch 4 taken 7954 times.
✓ Branch 6 taken 24503 times.
✓ Branch 7 taken 962096 times.
✓ Branch 8 taken 24503 times.
✓ Branch 9 taken 962096 times.
|
1042133 | for (auto dofIdx : asImp_().couplingStencil(i, element, j)) |
| 526 |
8/8✓ Branch 0 taken 1062304 times.
✓ Branch 1 taken 47580 times.
✓ Branch 2 taken 1062304 times.
✓ Branch 3 taken 47580 times.
✓ Branch 4 taken 470335 times.
✓ Branch 5 taken 24503 times.
✓ Branch 6 taken 470335 times.
✓ Branch 7 taken 24503 times.
|
1893054 | for (auto eIdx : connectedElementsCoupling[j][dofIdx]) |
| 527 | 3065278 | neighborColors.push_back(colors[eIdx]); | |
| 528 | } | ||
| 529 | }); | ||
| 530 | |||
| 531 | // find smallest color (positive integer) not in neighborColors | ||
| 532 |
2/4✓ Branch 1 taken 7954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 962096 times.
✗ Branch 5 not taken.
|
970050 | const auto color = Detail::smallestAvailableColor(neighborColors, colorUsed); |
| 533 | |||
| 534 | // assign color to element | ||
| 535 |
12/12✓ Branch 0 taken 1720 times.
✓ Branch 1 taken 15 times.
✓ Branch 2 taken 7560 times.
✓ Branch 3 taken 394 times.
✓ Branch 4 taken 1720 times.
✓ Branch 5 taken 11437 times.
✓ Branch 6 taken 942116 times.
✓ Branch 7 taken 558 times.
✓ Branch 8 taken 19409 times.
✓ Branch 9 taken 13 times.
✓ Branch 10 taken 11258 times.
✓ Branch 11 taken 164 times.
|
1941835 | colors[gg.elementMapper().index(element)] = color; |
| 536 | |||
| 537 | // add element to the set of elements with the same color | ||
| 538 |
8/8✓ Branch 0 taken 7560 times.
✓ Branch 1 taken 394 times.
✓ Branch 2 taken 7560 times.
✓ Branch 3 taken 394 times.
✓ Branch 4 taken 961361 times.
✓ Branch 5 taken 735 times.
✓ Branch 6 taken 961361 times.
✓ Branch 7 taken 735 times.
|
1940100 | if (color < elementSets.size()) |
| 539 |
5/10✓ Branch 2 taken 7560 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 7560 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 961361 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 8082 times.
✓ Branch 12 taken 953279 times.
✗ Branch 13 not taken.
|
977003 | elementSets[color].push_back(element.seed()); |
| 540 | else | ||
| 541 |
14/28✓ Branch 1 taken 394 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 394 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 394 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 394 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 394 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 735 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 142 times.
✓ Branch 21 taken 593 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 142 times.
✓ Branch 24 taken 593 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 142 times.
✓ Branch 27 taken 593 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 142 times.
✓ Branch 30 taken 593 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
|
3387 | elementSets.push_back(std::vector<ElementSeed<i>>{ element.seed() }); |
| 542 | } | ||
| 543 | }); | ||
| 544 | |||
| 545 |
3/8✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
26 | std::cout << Fmt::format("Colored in {} seconds:\n", timer.elapsed()); |
| 546 |
1/2✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
|
364 | forEach(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto i) |
| 547 | { | ||
| 548 |
12/26✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 26 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 15 taken 25 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 25 times.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 25 times.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 1 times.
✓ Branch 25 taken 25 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 25 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 25 times.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
|
208 | std::cout << Fmt::format( |
| 549 | "-- {} elements in subdomain {} with {} colors\n", | ||
| 550 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
|
52 | asImp_().problem(i).gridGeometry().gridView().size(0), |
| 551 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 26 times.
|
104 | i(), std::get<i>(elementSets_).size() |
| 552 | ); | ||
| 553 | }); | ||
| 554 | 26 | } | |
| 555 | |||
| 556 | /*! | ||
| 557 | * \brief Execute assembly kernel in parallel | ||
| 558 | * | ||
| 559 | * \param domainId the domain index of domain i | ||
| 560 | * \param assembleElement kernel function to execute for one element | ||
| 561 | */ | ||
| 562 | template<std::size_t i, class AssembleElementFunc> | ||
| 563 | 1168 | void assembleMultithreaded(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const | |
| 564 | { | ||
| 565 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 584 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 584 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 584 times.
|
3504 | if (std::get<i>(elementSets_).empty()) |
| 566 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Call computeColorsForAssembly before assembling in parallel!"); | |
| 567 | |||
| 568 | // make this element loop run in parallel | ||
| 569 | // for this we have to color the elements so that we don't get | ||
| 570 | // race conditions when writing into the global matrix or modifying grid variable caches | ||
| 571 | // each color can be assembled using multiple threads | ||
| 572 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 584 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 584 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 584 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 584 times.
|
1168 | const auto& grid = this->problem(domainId).gridGeometry().gridView().grid(); |
| 573 |
4/4✓ Branch 0 taken 25556 times.
✓ Branch 1 taken 584 times.
✓ Branch 2 taken 25556 times.
✓ Branch 3 taken 584 times.
|
108064 | for (const auto& elements : std::get<i>(elementSets_)) |
| 574 | { | ||
| 575 |
2/2✓ Branch 2 taken 24657 times.
✓ Branch 3 taken 899 times.
|
5315500 | Dumux::parallelFor(elements.size(), [&](const std::size_t n) |
| 576 | { | ||
| 577 |
0/20✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
|
5871196 | const auto element = grid.entity(elements[n]); |
| 578 |
2/13✓ Branch 1 taken 328960 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 328960 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
2606638 | assembleElement(element); |
| 579 | }); | ||
| 580 | } | ||
| 581 | 1168 | } | |
| 582 | |||
| 583 | /*! | ||
| 584 | * \brief Extended source stencil (for the bulk domain) | ||
| 585 | * \note The default returns an empty stencil but some of the embedded | ||
| 586 | * coupling managers use this | ||
| 587 | */ | ||
| 588 | ✗ | const CouplingStencil<bulkIdx>& extendedSourceStencil(std::size_t eIdx) const | |
| 589 | 96896 | { return asImp_().emptyStencil(bulkIdx); } | |
| 590 | |||
| 591 | protected: | ||
| 592 | using ParentType::curSol; | ||
| 593 | |||
| 594 | //! computes the vertex indices per element for the box method | ||
| 595 | template<std::size_t id> | ||
| 596 | 12 | void precomputeVertexIndices(Dune::index_constant<id> domainIdx) | |
| 597 | { | ||
| 598 | // fill helper structure for box discretization | ||
| 599 | if constexpr (isBox<domainIdx>()) | ||
| 600 | { | ||
| 601 | 12 | this->vertexIndices(domainIdx).resize(gridView(domainIdx).size(0)); | |
| 602 |
13/16✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✓ Branch 3 taken 15489 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 40 times.
✓ Branch 6 taken 15448 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 40 times.
✓ Branch 9 taken 8224 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 8224 times.
✓ Branch 12 taken 1 times.
✓ Branch 14 taken 8224 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 8224 times.
✗ Branch 18 not taken.
|
111224 | for (const auto& element : elements(gridView(domainIdx))) |
| 603 | { | ||
| 604 | 47420 | constexpr int dim = GridView<domainIdx>::dimension; | |
| 605 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 23710 times.
✓ Branch 3 taken 8224 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 8224 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8224 times.
✗ Branch 10 not taken.
|
47420 | const auto eIdx = this->problem(domainIdx).gridGeometry().elementMapper().index(element); |
| 606 |
1/2✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
|
47420 | this->vertexIndices(domainIdx, eIdx).resize(element.subEntities(dim)); |
| 607 |
5/5✓ Branch 0 taken 123648 times.
✓ Branch 1 taken 56606 times.
✓ Branch 2 taken 123648 times.
✓ Branch 3 taken 48382 times.
✓ Branch 4 taken 8224 times.
|
638776 | for (int i = 0; i < element.subEntities(dim); ++i) |
| 608 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 156544 times.
✓ Branch 3 taken 32896 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 32896 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 32896 times.
✗ Branch 10 not taken.
|
313088 | this->vertexIndices(domainIdx, eIdx)[i] = this->problem(domainIdx).gridGeometry().vertexMapper().subIndex(element, i, dim); |
| 609 | } | ||
| 610 | } | ||
| 611 | 12 | } | |
| 612 | |||
| 613 | //! compute the shape function for a given point and geometry | ||
| 614 | template<std::size_t i, class FVGG, class Geometry, class ShapeValues> | ||
| 615 | 81883 | void getShapeValues(Dune::index_constant<i> domainI, const FVGG& gridGeometry, const Geometry& geo, const GlobalPosition& globalPos, ShapeValues& shapeValues) | |
| 616 | { | ||
| 617 | if constexpr (FVGG::discMethod == DiscretizationMethods::box) | ||
| 618 | { | ||
| 619 | 81883 | const auto ipLocal = geo.local(globalPos); | |
| 620 |
1/12✗ Branch 0 not taken.
✓ Branch 1 taken 62275 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
81883 | const auto& localBasis = this->problem(domainI).gridGeometry().feCache().get(geo.type()).localBasis(); |
| 621 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
81883 | localBasis.evaluateFunction(ipLocal, shapeValues); |
| 622 | } | ||
| 623 | else | ||
| 624 | DUNE_THROW(Dune::InvalidStateException, "Shape values requested for other discretization than box!"); | ||
| 625 | 81883 | } | |
| 626 | |||
| 627 | //! Clear all internal data members | ||
| 628 | 27 | void clear() | |
| 629 | { | ||
| 630 | 81 | pointSources(bulkIdx).clear(); | |
| 631 | 81 | pointSources(lowDimIdx).clear(); | |
| 632 | 81 | couplingStencils(bulkIdx).clear(); | |
| 633 | 81 | couplingStencils(lowDimIdx).clear(); | |
| 634 | 81 | vertexIndices(bulkIdx).clear(); | |
| 635 | 81 | vertexIndices(lowDimIdx).clear(); | |
| 636 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
|
27 | pointSourceData_.clear(); |
| 637 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
|
27 | averageDistanceToBulkCell_.clear(); |
| 638 | |||
| 639 | 27 | idCounter_ = 0; | |
| 640 | 27 | } | |
| 641 | |||
| 642 | //! compute the intersections between the two grids | ||
| 643 | 27 | void glueGrids() | |
| 644 | { | ||
| 645 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 27 times.
|
27 | const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); |
| 646 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
|
27 | const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); |
| 647 | |||
| 648 | // intersect the bounding box trees | ||
| 649 | 108 | glue_->build(bulkGridGeometry.boundingBoxTree(), lowDimGridGeometry.boundingBoxTree()); | |
| 650 | 27 | } | |
| 651 | |||
| 652 | //! Return reference to point source data vector member | ||
| 653 | std::vector<PointSourceData>& pointSourceData() | ||
| 654 |
4/7✓ Branch 1 taken 61202 times.
✓ Branch 2 taken 3038 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
64249 | { return pointSourceData_; } |
| 655 | |||
| 656 | //! Return reference to average distances to bulk cell | ||
| 657 | std::vector<Scalar>& averageDistanceToBulkCell() | ||
| 658 |
2/3✓ Branch 1 taken 10310 times.
✓ Branch 2 taken 380 times.
✗ Branch 3 not taken.
|
10695 | { return averageDistanceToBulkCell_; } |
| 659 | |||
| 660 | //! Return the point source if domain i | ||
| 661 | template<std::size_t i> | ||
| 662 | ✗ | std::vector<PointSource<i>>& pointSources(Dune::index_constant<i> dom) | |
| 663 |
21/32✗ Branch 0 not taken.
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 27 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 28 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 27 times.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 84736 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✓ Branch 12 taken 84736 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✓ Branch 15 taken 84736 times.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 84737 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 81705 times.
✓ Branch 22 taken 1 times.
✓ Branch 23 taken 7 times.
✓ Branch 24 taken 81698 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 48640 times.
✓ Branch 27 taken 33058 times.
✓ Branch 28 taken 48640 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 33058 times.
✗ Branch 31 not taken.
|
174050 | { return std::get<i>(pointSources_); } |
| 664 | |||
| 665 | //! Return reference to bulk coupling stencil member of domain i | ||
| 666 | template<std::size_t i> | ||
| 667 | ✗ | CouplingStencils<i>& couplingStencils(Dune::index_constant<i> dom) | |
| 668 |
10/18✓ Branch 5 taken 90113 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 90113 times.
✓ Branch 9 taken 3038 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 90113 times.
✓ Branch 12 taken 3038 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 90113 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 58950 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 58950 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 236 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 236 times.
✗ Branch 27 not taken.
|
1037784 | { return std::get<i>(couplingStencils_); } |
| 669 | |||
| 670 | //! Return a reference to the vertex indices | ||
| 671 | template<std::size_t i> | ||
| 672 | ✗ | std::vector<GridIndex<i>>& vertexIndices(Dune::index_constant<i> dom, GridIndex<i> eIdx) | |
| 673 |
15/30✓ Branch 4 taken 42753 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 62209 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 62209 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 30006 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 10626 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 10626 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 236 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 236 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 76 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 76 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 76 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 76 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 76 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 76 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 76 times.
✗ Branch 47 not taken.
|
743128 | { return std::get<i>(vertexIndices_)[eIdx]; } |
| 674 | |||
| 675 | //! Return a reference to the vertex indices container | ||
| 676 | template<std::size_t i> | ||
| 677 | ✗ | std::vector<std::vector<GridIndex<i>>>& vertexIndices(Dune::index_constant<i> dom) | |
| 678 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
|
66 | { return std::get<i>(vertexIndices_); } |
| 679 | |||
| 680 | const GlueType& glue() const | ||
| 681 | 121 | { return *glue_; } | |
| 682 | |||
| 683 | //! Returns the implementation of the problem (i.e. static polymorphism) | ||
| 684 | Implementation &asImp_() | ||
| 685 | { return *static_cast<Implementation *>(this); } | ||
| 686 | |||
| 687 | //! \copydoc asImp_() | ||
| 688 | const Implementation &asImp_() const | ||
| 689 | { return *static_cast<const Implementation *>(this); } | ||
| 690 | |||
| 691 | //! id generator for point sources | ||
| 692 | std::size_t idCounter_ = 0; | ||
| 693 | |||
| 694 | private: | ||
| 695 | |||
| 696 | //! the point source in both domains | ||
| 697 | std::tuple<std::vector<PointSource<bulkIdx>>, std::vector<PointSource<lowDimIdx>>> pointSources_; | ||
| 698 | std::vector<PointSourceData> pointSourceData_; | ||
| 699 | std::vector<Scalar> averageDistanceToBulkCell_; | ||
| 700 | |||
| 701 | //! Stencil data | ||
| 702 | std::tuple<std::vector<std::vector<GridIndex<bulkIdx>>>, | ||
| 703 | std::vector<std::vector<GridIndex<lowDimIdx>>>> vertexIndices_; | ||
| 704 | std::tuple<CouplingStencils<bulkIdx>, CouplingStencils<lowDimIdx>> couplingStencils_; | ||
| 705 | std::tuple<CouplingStencil<bulkIdx>, CouplingStencil<lowDimIdx>> emptyStencil_; | ||
| 706 | |||
| 707 | //! The glue object | ||
| 708 | std::shared_ptr<GlueType> glue_; | ||
| 709 | |||
| 710 | //! integration order for coupling source | ||
| 711 | int integrationOrder_ = 1; | ||
| 712 | |||
| 713 | //! coloring for multithreaded assembly | ||
| 714 | std::tuple< | ||
| 715 | std::deque<std::vector<ElementSeed<bulkIdx>>>, | ||
| 716 | std::deque<std::vector<ElementSeed<lowDimIdx>>> | ||
| 717 | > elementSets_; | ||
| 718 | }; | ||
| 719 | |||
| 720 | } // end namespace Dumux | ||
| 721 | |||
| 722 | #endif | ||
| 723 |