3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
facet/cellcentered/mpfa/couplingmapper.hh
Go to the documentation of this file.
1/*****************************************************************************
2 * See the file COPYING for full copying permissions. *
3 * *
4 * This program is free software: you can redistribute it and/or modify *
5 * it under the terms of the GNU General Public License as published by *
6 * the Free Software Foundation, either version 2 of the License, or *
7 * (at your option) any later version. *
8 * *
9 * This program is distributed in the hope that it will be useful, *
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12 * GNU General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
16 *****************************************************************************/
22#ifndef DUMUX_CCMPFA_FACETCOUPLING_MAPPER_HH
23#define DUMUX_CCMPFA_FACETCOUPLING_MAPPER_HH
24
25#include <dune/common/indices.hh>
26#include <dune/common/float_cmp.hh>
27
33
34namespace Dumux {
35
48template<class BulkFVG, class LowDimFVG, std::size_t bulkId, std::size_t lowDimId>
49class FacetCouplingMapper<BulkFVG, LowDimFVG, bulkId, lowDimId, DiscretizationMethod::ccmpfa>
50: public virtual FacetCouplingMapperBase<BulkFVG, LowDimFVG, bulkId, lowDimId>
51{
53
54 using BulkFVElementGeometry = typename BulkFVG::LocalView;
55 using BulkSubControlVolumeFace = typename BulkFVG::SubControlVolumeFace;
56
57 using BulkIndexType = typename IndexTraits<typename BulkFVG::GridView>::GridIndex;
58 using LowDimIndexType = typename IndexTraits<typename LowDimFVG::GridView>::GridIndex;
59
60 using LowDimElement = typename LowDimFVG::GridView::template Codim<0>::Entity;
61 using GlobalPosition = typename LowDimElement::Geometry::GlobalCoordinate;
62 static constexpr int lowDimDim = LowDimFVG::GridView::dimension;
63
64public:
66 using ParentType::bulkGridId;
67 using ParentType::facetGridId;
68
77 template< class Embeddings >
78 void update(const BulkFVG& bulkFvGridGeometry,
79 const LowDimFVG& lowDimFvGridGeometry,
80 std::shared_ptr<const Embeddings> embeddings)
81 {
82 // forward to update function with instantiated vertex adapter
84 update(bulkFvGridGeometry, lowDimFvGridGeometry, embeddings, GridAdapter(embeddings));
85 }
86
96 template< class Embeddings, class CodimOneGridAdapter >
97 void update(const BulkFVG& bulkFvGridGeometry,
98 const LowDimFVG& lowDimFvGridGeometry,
99 std::shared_ptr<const Embeddings> embeddings,
100 const CodimOneGridAdapter& codimOneGridAdapter)
101 {
102 // define the policy how to add map entries for given lowdim element and adjoined entity indices
103 auto addCouplingEntryPolicy = [&] (auto&& adjoinedEntityIndices,
104 const LowDimElement& lowDimElement,
105 const LowDimFVG& lowDimFvGridGeometry,
106 const BulkFVG& bulkFvGridGeometry)
107 {
108 const auto lowDimElemIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement);
109 auto& lowDimData = this->couplingMap_(facetGridId, bulkGridId)[lowDimElemIdx];
110
111 // determine corner indices (in bulk grid indices)
112 const auto& lowDimGeometry = lowDimElement.geometry();
113 const auto numElementCorners = lowDimElement.subEntities(lowDimDim);
114 std::vector<BulkIndexType> elemCornerIndices(numElementCorners);
115 for (int i = 0; i < numElementCorners; ++i)
116 elemCornerIndices[i] = codimOneGridAdapter.bulkGridVertexIndex(lowDimElement.template subEntity<lowDimDim>(i));
117
118 // Find the scvfs in the adjoined entities coinciding with the low dim element
119 for (auto bulkElemIdx : adjoinedEntityIndices)
120 {
121 const auto bulkElement = bulkFvGridGeometry.element(bulkElemIdx);
122
123 auto fvGeometry = localView(bulkFvGridGeometry);
124 fvGeometry.bind(bulkElement);
125
126 std::vector<BulkIndexType> embeddedScvfIndices;
127 for (const auto& scvf : scvfs(fvGeometry))
128 {
129 // for non-boundary faces, it suffices to check if one
130 // of the outside scv indices is in the set of embedments
131 if (!scvf.boundary())
132 {
133 if ( std::count(adjoinedEntityIndices.begin(), adjoinedEntityIndices.end(), scvf.outsideScvIdx()) )
134 embeddedScvfIndices.push_back(scvf.index());
135 }
136
137 // otherwise, do float comparison of element and scvf facet corner
138 else
139 {
140 const auto eps = lowDimGeometry.volume()*1e-8;
141 const auto diffVec = lowDimGeometry.center()-scvf.facetCorner();
142
143 if ( Dune::FloatCmp::eq<GlobalPosition, Dune::FloatCmp::CmpStyle::absolute>(diffVec, GlobalPosition(0.0), eps) )
144 embeddedScvfIndices.push_back(scvf.index());
145 }
146 }
147
148 // Error tracking. The boundary scvf detection might has to be improved for very fine grids!?
149 if ( embeddedScvfIndices.size() != numElementCorners )
150 DUNE_THROW(Dune::InvalidStateException, "Could not find all coupling scvfs in embedment");
151
152 // add each dof in the low dim element to coupling stencil of the bulk element and vice versa
153 auto& bulkData = this->couplingMap_(bulkGridId, facetGridId)[bulkElemIdx];
154 const auto lowDimElementDofs = LowDimFVG::discMethod == DiscretizationMethod::box
155 ? this->extractNodalDofs_(lowDimElement, lowDimFvGridGeometry)
156 : std::vector<LowDimIndexType>( {lowDimElemIdx} );
157
158 // add coupling info to all elements/scvfs in interaction volumes
159 for (auto dofIdx : lowDimElementDofs)
160 for (auto scvfIdx : embeddedScvfIndices)
161 this->addCouplingsFromIV_(bulkFvGridGeometry, fvGeometry.scvf(scvfIdx), fvGeometry, lowDimElemIdx, dofIdx);
162
163 // sort the scvfs according to the corners of the low dim element if box is used
164 // that allows identifying which scvf flux enters which low dim scv later
165 if (LowDimFVG::discMethod == DiscretizationMethod::box)
166 {
167 const auto copy = embeddedScvfIndices;
168
169 for (unsigned int i = 0; i < numElementCorners; ++i)
170 {
171 const auto& scvf = fvGeometry.scvf(copy[i]);
172 auto it = std::find(elemCornerIndices.begin(), elemCornerIndices.end(), scvf.vertexIndex());
173 assert(it != elemCornerIndices.end());
174 embeddedScvfIndices[ std::distance(elemCornerIndices.begin(), it) ] = copy[i];
175 }
176 }
177
178 // add info on which scvfs coincide with which low dim element
179 auto& elemToScvfMap = bulkData.elementToScvfMap[lowDimElemIdx];
180 elemToScvfMap.insert(elemToScvfMap.end(), embeddedScvfIndices.begin(), embeddedScvfIndices.end());
181
182 // add embedment
183 lowDimData.embedments.emplace_back( bulkElemIdx, std::move(embeddedScvfIndices) );
184 }
185 };
186
187 // let the parent do the update subject to the execution policy defined above
188 ParentType::update_(bulkFvGridGeometry, lowDimFvGridGeometry, embeddings, addCouplingEntryPolicy);
189
190 // lambda to make a container unique
191 auto makeUnique = [] (auto& c)
192 {
193 std::sort(c.begin(), c.end());
194 c.erase( std::unique(c.begin(), c.end()), c.end() );
195 };
196
197 // lambda to make bulk coupling map entry unique
198 auto makeBulkMapEntryUnique = [&makeUnique] (auto& data)
199 {
200 makeUnique(data.second.couplingStencil);
201 makeUnique(data.second.couplingElementStencil);
202 std::for_each(data.second.dofToCouplingScvfMap.begin(),
203 data.second.dofToCouplingScvfMap.end(),
204 [&makeUnique] (auto& pair) { makeUnique(pair.second); });
205 };
206
207 // make bulk map unique
208 auto& bulkCouplingData = this->couplingMap_(bulkGridId, facetGridId);
209 std::for_each(bulkCouplingData.begin(), bulkCouplingData.end(), makeBulkMapEntryUnique);
210
211 // coupling stencils might not be unique in lowdim domain
212 auto& lowDimCouplingData = this->couplingMap_(facetGridId, bulkGridId);
213 std::for_each(lowDimCouplingData.begin(),
214 lowDimCouplingData.end(),
215 [&makeUnique] (auto& pair) { makeUnique(pair.second.couplingStencil); });
216 }
217
218private:
219 void addCouplingsFromIV_(const BulkFVG& bulkFvGridGeometry,
220 const BulkSubControlVolumeFace& scvf,
221 const BulkFVElementGeometry& fvGeometry,
222 LowDimIndexType lowDimElemIdx,
223 LowDimIndexType lowDimDofIdx)
224 {
225 const auto& gridIvIndexSets = bulkFvGridGeometry.gridInteractionVolumeIndexSets();
226 if (bulkFvGridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
227 addCouplingsFromIVIndexSet_(gridIvIndexSets.secondaryIndexSet(scvf), fvGeometry, lowDimElemIdx, lowDimDofIdx);
228 else
229 addCouplingsFromIVIndexSet_(gridIvIndexSets.primaryIndexSet(scvf), fvGeometry, lowDimElemIdx, lowDimDofIdx);
230 }
231
232 template< class IVIndexSet >
233 void addCouplingsFromIVIndexSet_(const IVIndexSet& indexSet,
234 const BulkFVElementGeometry& fvGeometry,
235 LowDimIndexType lowDimElemIdx,
236 LowDimIndexType lowDimDofIdx)
237 {
238 auto& lowDimData = this->couplingMap_(facetGridId, bulkGridId)[lowDimElemIdx];
239 for (auto bulkElemIdx : indexSet.gridScvIndices())
240 {
241 auto& bulkData = this->couplingMap_(bulkGridId, facetGridId)[bulkElemIdx];
242
243 lowDimData.couplingStencil.push_back( bulkElemIdx );
244 bulkData.couplingStencil.push_back( lowDimDofIdx );
245 bulkData.couplingElementStencil.push_back( lowDimElemIdx );
246 }
247
248 // insert scvf couplings stemming from this interaction volume
249 for (auto scvfIdx : indexSet.gridScvfIndices())
250 {
251 const auto& scvf = fvGeometry.scvf(scvfIdx);
252 auto& bulkData = this->couplingMap_(bulkGridId, facetGridId)[scvf.insideScvIdx()];
253 auto& couplingScvfs = bulkData.dofToCouplingScvfMap[lowDimDofIdx];
254 couplingScvfs.insert(couplingScvfs.end(), scvfIdx);
255 }
256 }
257};
258
259} // end namespace Dumux
260
261#endif
Defines the index types used for grid and local indices.
The available discretization methods in Dumux.
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
Definition: adapt.hh:29
typename GridView::IndexSet::IndexType GridIndex
Definition: indextraits.hh:39
void update(const BulkFVG &bulkFvGridGeometry, const LowDimFVG &lowDimFvGridGeometry, std::shared_ptr< const Embeddings > embeddings)
Update coupling maps. This is the standard interface required by any mapper implementation.
Definition: facet/cellcentered/mpfa/couplingmapper.hh:78
void update(const BulkFVG &bulkFvGridGeometry, const LowDimFVG &lowDimFvGridGeometry, std::shared_ptr< const Embeddings > embeddings, const CodimOneGridAdapter &codimOneGridAdapter)
Update coupling maps. This is the standard interface required by any mapper implementation.
Definition: facet/cellcentered/mpfa/couplingmapper.hh:97
Adapter that allows retrieving information on a d-dimensional grid for entities of a (d-1)-dimensiona...
Definition: codimonegridadapter.hh:52
BulkIndexType bulkGridVertexIndex(const FacetGridVertex &v) const
Returns the index within the d-dimensional grid of a vertex of the (d-1)-dimensional grid.
Definition: codimonegridadapter.hh:145
Implementation for the coupling mapper that sets up and stores the coupling maps between two domains ...
Definition: facet/couplingmapper.hh:52
Base class for the coupling mapper that sets up and stores the coupling maps between two domains of d...
Definition: couplingmapperbase.hh:51