3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
facet/box/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 3 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_BOX_FACETCOUPLING_MAPPER_HH
23#define DUMUX_BOX_FACETCOUPLING_MAPPER_HH
24
25#include <dune/common/indices.hh>
26
32
33namespace Dumux {
34
46template<class BulkFVG, class LowDimFVG, std::size_t bulkId, std::size_t lowDimId>
47class FacetCouplingMapper<BulkFVG, LowDimFVG, bulkId, lowDimId, DiscretizationMethod::box>
48: public virtual FacetCouplingMapperBase<BulkFVG, LowDimFVG, bulkId, lowDimId>
49{
51 using LowDimElement = typename LowDimFVG::GridView::template Codim<0>::Entity;
52
53 // dimensions of the two grids
54 using BulkGridView = typename BulkFVG::GridView;
55 using LowDimGridView = typename LowDimFVG::GridView;
56 static constexpr int bulkDim = BulkGridView::dimension;
57 static constexpr int lowDimDim = LowDimGridView::dimension;
58
59public:
61 using ParentType::bulkGridId;
62 using ParentType::facetGridId;
63
72 template< class Embeddings >
73 void update(const BulkFVG& bulkFvGridGeometry,
74 const LowDimFVG& lowDimFvGridGeometry,
75 std::shared_ptr<const Embeddings> embeddings)
76 {
77 // forward to update function with instantiated vertex adapter
79 update(bulkFvGridGeometry, lowDimFvGridGeometry, embeddings, GridAdapter(embeddings));
80 }
81
90 template< class Embeddings, class CodimOneGridAdapter >
91 void update(const BulkFVG& bulkFvGridGeometry,
92 const LowDimFVG& lowDimFvGridGeometry,
93 std::shared_ptr<const Embeddings> embeddings,
94 const CodimOneGridAdapter& codimOneGridAdapter)
95 {
96 // define the policy how to add map entries for given lowdim element and adjoined entity indices
97 auto addCouplingEntryPolicy = [&] (auto&& adjoinedEntityIndices,
98 const LowDimElement& lowDimElement,
99 const LowDimFVG& lowDimFvGridGeometry,
100 const BulkFVG& bulkFvGridGeometry)
101 {
102 using LowDimIndexType = typename IndexTraits<LowDimGridView>::GridIndex;
103 using BulkIndexType = typename IndexTraits<BulkGridView>::GridIndex;
104
105 const auto lowDimElemIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement);
106 auto& lowDimData = this->couplingMap_(facetGridId, bulkGridId)[lowDimElemIdx];
107
108 // determine corner indices (in bulk grid indices)
109 const auto& eg = lowDimElement.geometry();
110 const auto numElementCorners = eg.corners();
111 std::vector<BulkIndexType> elementCorners(numElementCorners);
112 for (int i = 0; i < numElementCorners; ++i)
113 elementCorners[i] = codimOneGridAdapter.bulkGridVertexIndex(lowDimElement.template subEntity<lowDimDim>(i));
114
115 // save unsorted set of corner indices and search scvfs in adjoined entities
116 const auto unsortedElemCorners = elementCorners;
117 std::sort(elementCorners.begin(), elementCorners.end());
118 for (auto bulkElemIdx : adjoinedEntityIndices)
119 {
120 const auto bulkElement = bulkFvGridGeometry.element(bulkElemIdx);
121 const auto bulkRefElem = referenceElement(bulkElement);
122
123 // find the bulk element facet that lies on this low dim element (assumes conformity!)
124 bool found = false;
125 unsigned int coupledFacetIndex;
126 std::vector<unsigned int> handledFacets;
127 for (const auto& is : intersections(bulkFvGridGeometry.gridView(), bulkElement))
128 {
129 // skip already handled facets (necessary for e.g. Dune::FoamGrid)
130 if (std::count(handledFacets.begin(), handledFacets.end(), is.indexInInside()))
131 continue;
132
133 handledFacets.push_back(is.indexInInside());
134
135 // determine if it lies on low dim element by comparing corner indices
136 const auto numCorners = is.geometry().corners();
137 std::vector<BulkIndexType> facetIndices(numCorners);
138 for (int i = 0; i < numCorners; ++i)
139 {
140 const auto vIdxLocal = bulkRefElem.subEntity(is.indexInInside(), 1, i, bulkDim);
141 facetIndices[i] = bulkFvGridGeometry.vertexMapper().vertexIndex(bulkElement.template subEntity<bulkDim>(vIdxLocal));
142 }
143
144 std::sort(facetIndices.begin(), facetIndices.end());
145 if ( std::equal(facetIndices.begin(), facetIndices.end(), elementCorners.begin(), elementCorners.end()) )
146 {
147 coupledFacetIndex = is.indexInInside();
148 found = true; break;
149 }
150 }
151
152 // ensure that we found the facet!
153 if (!found)
154 DUNE_THROW(Dune::InvalidStateException, "Could not find the bulk element coupling facet!");
155
156 // we should always find numElementCorners coupling scvfs
157 auto fvGeometry = localView(bulkFvGridGeometry);
158 fvGeometry.bindElement(bulkElement);
159
160 unsigned int foundCounter = 0;
161 std::vector<BulkIndexType> embeddedScvfIndices(numElementCorners);
162 for (const auto& scvf : scvfs(fvGeometry))
163 {
164 if (scvf.interiorBoundary() && scvf.facetIndexInElement() == coupledFacetIndex)
165 {
166 // we want to order the set of scvfs lying on the lower-dimensional element such that the i-th scvf
167 // coincides with the i-th low dim element corner. This will help later to identify which scvf fluxes
168 // enter which scv of the low dim element if the lower-dimensional domain uses the box scheme
169 const auto vIdxLocal = bulkRefElem.subEntity(coupledFacetIndex, 1, scvf.indexInElementFacet(), bulkDim);
170 const auto vIdxGlobal = bulkFvGridGeometry.vertexMapper().vertexIndex(bulkElement, vIdxLocal, bulkDim);
171 const auto it = std::find(unsortedElemCorners.begin(), unsortedElemCorners.end(), vIdxGlobal);
172 assert(it != unsortedElemCorners.end());
173 const auto lowDimElemLocalCornerIdx = std::distance(unsortedElemCorners.begin(), it);
174 embeddedScvfIndices[lowDimElemLocalCornerIdx] = scvf.index();
175 foundCounter++;
176 }
177 }
178
179 // ensure we found all scvfs
180 if (foundCounter != numElementCorners)
181 DUNE_THROW(Dune::InvalidStateException, "Found " << foundCounter << " instead of " << numElementCorners << " coupling scvfs in the bulk element");
182
183 // add each dof in the low dim element to coupling stencil of the bulk element
184 auto& bulkData = this->couplingMap_(bulkGridId, facetGridId)[bulkElemIdx];
185 const auto lowDimElementDofs = LowDimFVG::discMethod == DiscretizationMethod::box
186 ? this->extractNodalDofs_(lowDimElement, lowDimFvGridGeometry)
187 : std::vector<LowDimIndexType>({lowDimElemIdx});
188
189 for (auto lowDimDofIdx : lowDimElementDofs)
190 {
191 bulkData.couplingStencil.push_back(lowDimDofIdx);
192 auto& couplingScvfs = bulkData.dofToCouplingScvfMap[lowDimDofIdx];
193 couplingScvfs.insert(couplingScvfs.end(), embeddedScvfIndices.begin(), embeddedScvfIndices.end());
194 }
195
196 // add info on which scvfs coincide with which low dim element
197 bulkData.couplingElementStencil.push_back(lowDimElemIdx);
198 auto& elemToScvfMap = bulkData.elementToScvfMap[lowDimElemIdx];
199 elemToScvfMap.insert(elemToScvfMap.end(), embeddedScvfIndices.begin(), embeddedScvfIndices.end());
200
201 // add embedment and the bulk cell dofs to coupling stencil of low dim element
202 lowDimData.embedments.emplace_back(bulkElemIdx, std::move(embeddedScvfIndices));
203
204 const auto bulkElementDofs = this->extractNodalDofs_(bulkElement, bulkFvGridGeometry);
205 for (auto bulkDofIdx : bulkElementDofs)
206 lowDimData.couplingStencil.push_back(bulkDofIdx);
207 }
208 };
209
210 // let the parent do the update subject to the execution policy defined above
211 ParentType::update_(bulkFvGridGeometry, lowDimFvGridGeometry, embeddings, addCouplingEntryPolicy);
212
213 // coupling stencils might not be unique with the policy above
214 auto makeStencilUnique = [] (auto& data)
215 {
216 auto& cs = data.second.couplingStencil;
217 std::sort(cs.begin(), cs.end());
218 cs.erase( std::unique(cs.begin(), cs.end()), cs.end() );
219 };
220
221 auto& lowDimCouplingData = this->couplingMap_(facetGridId, bulkGridId);
222 std::for_each(lowDimCouplingData.begin(), lowDimCouplingData.end(), makeStencilUnique);
223
224 // bulk coupling stencil is only non-unique if box is used
225 if (LowDimFVG::discMethod == DiscretizationMethod::box)
226 {
227 auto& bulkCouplingData = this->couplingMap_(bulkGridId, facetGridId);
228 std::for_each(bulkCouplingData.begin(), bulkCouplingData.end(), makeStencilUnique);
229 }
230 }
231};
232
233} // end namespace Dumux
234
235#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/box/couplingmapper.hh:73
void update(const BulkFVG &bulkFvGridGeometry, const LowDimFVG &lowDimFvGridGeometry, std::shared_ptr< const Embeddings > embeddings, const CodimOneGridAdapter &codimOneGridAdapter)
Update coupling maps with a given grid adapter.
Definition: facet/box/couplingmapper.hh:91
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