20#ifndef OPM_BOUNDARYPERIODICITY_HEADER_INCLUDED
21#define OPM_BOUNDARYPERIODICITY_HEADER_INCLUDED
24#include <dune/common/fvector.hh>
26#include <opm/common/ErrorMacros.hpp>
64 return cmpval() < other.cmpval();
70 const double pi = 3.14159265358979323846264338327950288;
82 bool match(std::vector<BoundaryFaceInfo>& bfaces,
int face,
int lower,
int upper);
85 template <
class Gr
idView>
87 std::array<double, 6>& side_areas,
89 const std::array<bool, 2*GridView::dimension>& is_periodic,
90 double spatial_tolerance = 1e-6)
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;
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]);
111 max_bid = std::max(max_bid, f->boundaryId());
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.");
121 std::fill(side_areas.begin(), side_areas.end(), 0.0);
123 bfinfo.reserve(num_bdy);
124 for (
int i = 0; i < num_bdy; ++i) {
127 bf.
bid = bface_iters[i]->boundaryId();
131 bf.
area = bface_iters[i]->geometry().volume();
132 bf.
centroid = bface_iters[i]->geometry().center();
133 for (
int dd = 0; dd < dim; ++dd) {
135 if (fabs(coord - low[dd]) <= spatial_tolerance) {
138 }
else if (fabs(coord - hi[dd]) <= spatial_tolerance) {
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?");
151 bfinfo.push_back(bf);
153 assert(bfinfo.size() == bface_iters.size());
156 std::sort(bfinfo.begin(), bfinfo.end());
159 for (
int i = 0; i < num_bdy; ++i) {
160 if (bfinfo[i].partner_face_index != -1) {
163 if (!is_periodic[bfinfo[i].canon_pos]) {
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);
171 ok =
match(bfinfo, i, 0, num_bdy);
173 OPM_MESSAGE(
"Warning: No partner found for boundary id " << bfinfo[i].bid);
Definition: BlackoilFluid.hpp:32
bool match(std::vector< BoundaryFaceInfo > &bfaces, int face, int lower, int upper)
Find a match (periodic partner) for the given face.
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
Definition: BoundaryPeriodicity.hpp:38
double area
Face area.
Definition: BoundaryPeriodicity.hpp:55
bool operator<(const BoundaryFaceInfo &other) const
Definition: BoundaryPeriodicity.hpp:62
int bid
Boundary id of this face.
Definition: BoundaryPeriodicity.hpp:42
int face_index
Face index in [0, ..., #faces - 1].
Definition: BoundaryPeriodicity.hpp:40
int partner_bid
Boundary id of periodic partner face, or 0 if no parner.
Definition: BoundaryPeriodicity.hpp:53
int partner_face_index
Face index of periodic partner, or -1 if no partner.
Definition: BoundaryPeriodicity.hpp:51
Dune::FieldVector< double, 3 > centroid
Face centroid.
Definition: BoundaryPeriodicity.hpp:57
int canon_pos
Definition: BoundaryPeriodicity.hpp:49