From 1f27a85f33d293f73c286c1bd9136f133f3c36d9 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 12 Feb 2019 16:27:20 +0100
Subject: [PATCH] [temp] Integration point source box
---
.../embedded/integrationpointsource.hh | 63 +++++++++++++------
1 file changed, 44 insertions(+), 19 deletions(-)
diff --git a/dumux/multidomain/embedded/integrationpointsource.hh b/dumux/multidomain/embedded/integrationpointsource.hh
index 4569845609..b4b3889b9c 100644
--- a/dumux/multidomain/embedded/integrationpointsource.hh
+++ b/dumux/multidomain/embedded/integrationpointsource.hh
@@ -122,32 +122,57 @@ public:
// loop over all intersected elements
for (unsigned int eIdx : entities)
{
- if (FVGridGeometry::discMethod == DiscretizationMethod::box)
+ if constexpr(FVGridGeometry::discMethod == DiscretizationMethod::box)
{
// check in which subcontrolvolume(s) we are
const auto element = fvGridGeometry.boundingBoxTree().entitySet().entity(eIdx);
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
const auto globalPos = source.position();
- // loop over all sub control volumes and check if the point source is inside
- constexpr int dim = FVGridGeometry::GridView::dimension;
- Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
- for (auto&& scv : scvs(fvGeometry))
- if (intersectsPointGeometry(globalPos, scv.geometry()))
- scvIndices.push_back(scv.indexInElement());
-
- // for all scvs that where tested positiv add the point sources
- // to the element/scv to point source map
- for (auto scvIdx : scvIndices)
+
+ static const bool boxPointSourceLumping = getParam<bool>("MixedDimension.EnableBoxPointSourceLumping", true);
+ if (boxPointSourceLumping)
+ {
+ // loop over all sub control volumes and check if the point source is inside
+ constexpr int dim = FVGridGeometry::GridView::dimension;
+ Dune::ReservedVector<std::size_t, 1<<dim> scvIndices;
+ for (auto&& scv : scvs(fvGeometry))
+ if (intersectsPointGeometry(globalPos, scv.geometry()))
+ scvIndices.push_back(scv.indexInElement());
+
+ // for all scvs that where tested positiv add the point sources
+ // to the element/scv to point source map
+ for (auto scvIdx : scvIndices)
+ {
+ const auto key = std::make_pair(eIdx, scvIdx);
+ if (pointSourceMap.count(key))
+ pointSourceMap.at(key).push_back(source);
+ else
+ pointSourceMap.insert({key, {source}});
+ // split equally on the number of matched scvs
+ auto& s = pointSourceMap.at(key).back();
+ s.setEmbeddings(scvIndices.size()*s.embeddings());
+ }
+ }
+ // distribute the sources according to the basis function weights
+ else
{
- const auto key = std::make_pair(eIdx, scvIdx);
- if (pointSourceMap.count(key))
- pointSourceMap.at(key).push_back(source);
- else
- pointSourceMap.insert({key, {source}});
- // split equally on the number of matched scvs
- auto& s = pointSourceMap.at(key).back();
- s.setEmbeddings(scvIndices.size()*s.embeddings());
+ using Scalar = std::decay_t<decltype(source.values()[0])>;
+ using ShapeValues = std::vector<typename Dune::FieldVector<Scalar, 1>>;
+ ShapeValues shapeValues;
+ const auto& localBasis = fvGeometry.feLocalBasis();
+ const auto ipLocal = element.geometry().local(globalPos);
+ localBasis.evaluateFunction(ipLocal, shapeValues);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto key = std::make_pair(eIdx, scv.indexInElement());
+ if (pointSourceMap.count(key))
+ pointSourceMap.at(key).push_back(source);
+ else
+ pointSourceMap.insert({key, {source}});
+ auto& s = pointSourceMap.at(key).back();
+ s *= shapeValues[scv.indexInElement()];
+ }
}
}
else
--
GitLab