BoundaryPeriodicity.hpp
Go to the documentation of this file.
1/*
2 Copyright 2011 SINTEF ICT, Applied Mathematics.
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18*/
19
20#ifndef OPM_BOUNDARYPERIODICITY_HEADER_INCLUDED
21#define OPM_BOUNDARYPERIODICITY_HEADER_INCLUDED
22
23
24#include <dune/common/fvector.hh>
25#include <array>
26#include <opm/common/ErrorMacros.hpp>
27
28#include <algorithm>
29#include <vector>
30#include <iostream>
31
32namespace Opm
33{
34
37 struct BoundaryFaceInfo
38 {
40 int face_index;
42 int bid;
49 int canon_pos;
51 int partner_face_index;
53 int partner_bid;
55 double area;
57 Dune::FieldVector<double,3> centroid;
58
62 bool operator<(const BoundaryFaceInfo& other) const
63 {
64 return cmpval() < other.cmpval();
65 }
66
67 private:
68 double cmpval() const
69 {
70 const double pi = 3.14159265358979323846264338327950288;
71 return centroid[(canon_pos/2 + 1)%3] + pi*centroid[(canon_pos/2 + 2)%3];
72 }
73 };
74
75
82 bool match(std::vector<BoundaryFaceInfo>& bfaces, int face, int lower, int upper);
83
85 template <class GridView>
86 void findPeriodicPartners(std::vector<BoundaryFaceInfo>& bfinfo,
87 std::array<double, 6>& side_areas,
88 const GridView& g,
89 const std::array<bool, 2*GridView::dimension>& is_periodic,
90 double spatial_tolerance = 1e-6)
91 {
92 // Pick out all boundary faces, simultaneously find the
93 // bounding box of their centroids, and the max id.
94 const int dim = GridView::dimension;
95 typedef typename GridView::template Codim<0>::Iterator CI;
96 typedef typename GridView::IntersectionIterator FI;
97 typedef Dune::FieldVector<double, dim> Vector;
98 std::vector<FI> bface_iters;
99 Vector low(1e100);
100 Vector hi(-1e100);
101 int max_bid = 0;
102 for (CI c = g.template begin<0>(); c != g.template end<0>(); ++c) {
103 for (FI f = g.ibegin(*c); f != g.iend(*c); ++f) {
104 if (f->boundaryId()) {
105 bface_iters.push_back(f);
106 Vector fcent = f->geometry().center();
107 for (int dd = 0; dd < dim; ++dd) {
108 low[dd] = std::min(low[dd], fcent[dd]);
109 hi[dd] = std::max(hi[dd], fcent[dd]);
110 }
111 max_bid = std::max(max_bid, f->boundaryId());
112 }
113 }
114 }
115 int num_bdy = bface_iters.size();
116 if (max_bid != num_bdy) {
117 OPM_THROW(std::runtime_error, "createPeriodic() assumes that every boundary face has a unique boundary id. That seems to be violated.");
118 }
119
120 // Store boundary face info in a suitable structure. Also find side total volumes.
121 std::fill(side_areas.begin(), side_areas.end(), 0.0);
122 bfinfo.clear();
123 bfinfo.reserve(num_bdy);
124 for (int i = 0; i < num_bdy; ++i) {
125 BoundaryFaceInfo bf;
126 bf.face_index = i;
127 bf.bid = bface_iters[i]->boundaryId();
128 bf.canon_pos = -1;
129 bf.partner_face_index = -1;
130 bf.partner_bid = 0;
131 bf.area = bface_iters[i]->geometry().volume();
132 bf.centroid = bface_iters[i]->geometry().center();
133 for (int dd = 0; dd < dim; ++dd) {
134 double coord = bf.centroid[dd];
135 if (fabs(coord - low[dd]) <= spatial_tolerance) {
136 bf.canon_pos = 2*dd;
137 break;
138 } else if (fabs(coord - hi[dd]) <= spatial_tolerance) {
139 bf.canon_pos = 2*dd + 1;
140 break;
141 }
142 }
143 if (bf.canon_pos == -1) {
144 std::cerr << "Centroid: " << bf.centroid << "\n";
145 std::cerr << "Bounding box min: " << low << "\n";
146 std::cerr << "Bounding box max: " << hi << "\n";
147 OPM_THROW(std::runtime_error, "Boundary face centroid not on bounding box. Maybe the grid is not an axis-aligned shoe-box?");
148 }
149 side_areas[bf.canon_pos] += bf.area;
150 bf.centroid[bf.canon_pos/2] = 0.0;
151 bfinfo.push_back(bf);
152 }
153 assert(bfinfo.size() == bface_iters.size());
154
155 // Sort the infos so that partners end up close.
156 std::sort(bfinfo.begin(), bfinfo.end());
157
158 // Identify partners.
159 for (int i = 0; i < num_bdy; ++i) {
160 if (bfinfo[i].partner_face_index != -1) {
161 continue;
162 }
163 if (!is_periodic[bfinfo[i].canon_pos]) {
164 continue;
165 }
166 int lower = std::max(0, i - 10);
167 int upper = std::min(num_bdy, i + 10);
168 bool ok = match(bfinfo, i, lower, upper);
169 if (!ok) {
170 // We have not found a partner.
171 ok = match(bfinfo, i, 0, num_bdy);
172 if (!ok) {
173 OPM_MESSAGE("Warning: No partner found for boundary id " << bfinfo[i].bid);
174 // OPM_THROW(std::runtime_error, "No partner found.");
175 }
176 }
177 }
178 }
179
180 /*
182 template <class GridInterface>
183 void findPeriodicPartners(std::vector<BoundaryFaceInfo>& bfinfo,
184 std::array<double, 6>& side_areas,
185 const GridInterface& g,
186 const std::array<bool, 2*GridInterface::Dimension>& is_periodic,
187 double spatial_tolerance = 1e-6)
188 {
189 // Pick out all boundary faces, simultaneously find the
190 // bounding box of their centroids, and the max id.
191 typedef typename GridInterface::CellIterator CI;
192 typedef typename CI::FaceIterator FI;
193 typedef typename GridInterface::Vector Vector;
194 std::vector<FI> bface_iters;
195 Vector low(1e100);
196 Vector hi(-1e100);
197 int max_bid = 0;
198 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
199 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
200 if (f->boundaryId()) {
201 bface_iters.push_back(f);
202 Vector fcent = f->centroid();
203 for (int dd = 0; dd < GridInterface::Dimension; ++dd) {
204 low[dd] = std::min(low[dd], fcent[dd]);
205 hi[dd] = std::max(hi[dd], fcent[dd]);
206 }
207 max_bid = std::max(max_bid, f->boundaryId());
208 }
209 }
210 }
211 int num_bdy = bface_iters.size();
212 if (max_bid != num_bdy) {
213 OPM_THROW(std::runtime_error, "createPeriodic() assumes that every boundary face has a unique boundary id. That seems to be violated.");
214 }
215
216 // Store boundary face info in a suitable structure. Also find side total volumes.
217 std::fill(side_areas.begin(), side_areas.end(), 0.0);
218 bfinfo.clear();
219 bfinfo.reserve(num_bdy);
220 for (int i = 0; i < num_bdy; ++i) {
221 BoundaryFaceInfo bf;
222 bf.face_index = i;
223 bf.bid = bface_iters[i]->boundaryId();
224 bf.canon_pos = -1;
225 bf.partner_face_index = -1;
226 bf.partner_bid = 0;
227 bf.area = bface_iters[i]->area();
228 bf.centroid = bface_iters[i]->centroid();
229 for (int dd = 0; dd < GridInterface::Dimension; ++dd) {
230 double coord = bf.centroid[dd];
231 if (fabs(coord - low[dd]) <= spatial_tolerance) {
232 bf.canon_pos = 2*dd;
233 break;
234 } else if (fabs(coord - hi[dd]) <= spatial_tolerance) {
235 bf.canon_pos = 2*dd + 1;
236 break;
237 }
238 }
239 if (bf.canon_pos == -1) {
240 std::cerr << "Centroid: " << bf.centroid << "\n";
241 std::cerr << "Bounding box min: " << low << "\n";
242 std::cerr << "Bounding box max: " << hi << "\n";
243 OPM_THROW(std::runtime_error, "Boundary face centroid not on bounding box. Maybe the grid is not an axis-aligned shoe-box?");
244 }
245 side_areas[bf.canon_pos] += bf.area;
246 bf.centroid[bf.canon_pos/2] = 0.0;
247 bfinfo.push_back(bf);
248 }
249 assert(bfinfo.size() == bface_iters.size());
250
251 // Sort the infos so that partners end up close.
252 std::sort(bfinfo.begin(), bfinfo.end());
253
254 // Identify partners.
255 for (int i = 0; i < num_bdy; ++i) {
256 if (bfinfo[i].partner_face_index != -1) {
257 continue;
258 }
259 if (!is_periodic[bfinfo[i].canon_pos]) {
260 continue;
261 }
262 int lower = std::max(0, i - 10);
263 int upper = std::min(num_bdy, i + 10);
264 bool ok = match(bfinfo, i, lower, upper);
265 if (!ok) {
266 // We have not found a partner.
267 ok = match(bfinfo, i, 0, num_bdy);
268 if (!ok) {
269 OPM_MESSAGE("Warning: No partner found for boundary id " << bfinfo[i].bid);
270 // OPM_THROW(std::runtime_error, "No partner found.");
271 }
272 }
273 }
274 }
275*/
276
277
278
279} // namespace Opm
280
281#endif // OPM_BOUNDARYPERIODICITY_HEADER_INCLUDED
Opm
Definition: BlackoilFluid.hpp:32
Opm::match
bool match(std::vector< BoundaryFaceInfo > &bfaces, int face, int lower, int upper)
Find a match (periodic partner) for the given face.
Opm::findPeriodicPartners
void findPeriodicPartners(std::vector< BoundaryFaceInfo > &bfinfo, std::array< double, 6 > &side_areas, const GridView &g, const std::array< bool, 2 *GridView::dimension > &is_periodic, double spatial_tolerance=1e-6)
Common implementation for the various createPeriodic functions.
Definition: BoundaryPeriodicity.hpp:86
Opm::BoundaryFaceInfo
Definition: BoundaryPeriodicity.hpp:38
Opm::BoundaryFaceInfo::area
double area
Face area.
Definition: BoundaryPeriodicity.hpp:55
Opm::BoundaryFaceInfo::operator<
bool operator<(const BoundaryFaceInfo &other) const
Definition: BoundaryPeriodicity.hpp:62
Opm::BoundaryFaceInfo::bid
int bid
Boundary id of this face.
Definition: BoundaryPeriodicity.hpp:42
Opm::BoundaryFaceInfo::face_index
int face_index
Face index in [0, ..., #faces - 1].
Definition: BoundaryPeriodicity.hpp:40
Opm::BoundaryFaceInfo::partner_bid
int partner_bid
Boundary id of periodic partner face, or 0 if no parner.
Definition: BoundaryPeriodicity.hpp:53
Opm::BoundaryFaceInfo::partner_face_index
int partner_face_index
Face index of periodic partner, or -1 if no partner.
Definition: BoundaryPeriodicity.hpp:51
Opm::BoundaryFaceInfo::centroid
Dune::FieldVector< double, 3 > centroid
Face centroid.
Definition: BoundaryPeriodicity.hpp:57
Opm::BoundaryFaceInfo::canon_pos
int canon_pos
Definition: BoundaryPeriodicity.hpp:49