ArbitraryDomain.cpp 15.5 KB
Newer Older
1
//
2
// Class ArbitraryDomain
kraus's avatar
kraus committed
3
//   Interface to iterative solver and boundary geometry
4 5
//   for space charge calculation
//
6
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
7 8
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//                      2016, Daniel Winklehner, Massachusetts Institute of Technology
9
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
10
// All rights reserved
11
//
12 13 14 15 16
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
frey_m's avatar
frey_m committed
17 18 19 20 21 22 23 24 25 26
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
27 28
//

29
//#define DEBUG_INTERSECT_RAY_BOUNDARY
30

kraus's avatar
kraus committed
31
#include "Solvers/ArbitraryDomain.h"
snuverink_j's avatar
snuverink_j committed
32
#include "Structure/BoundaryGeometry.h"
kraus's avatar
kraus committed
33

34 35
#include <cmath>
#include <iostream>
snuverink_j's avatar
snuverink_j committed
36
#include <tuple>
37
#include <cassert>
38
#include "Utilities/OpalException.h"
gsell's avatar
gsell committed
39

40
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
kraus's avatar
kraus committed
41
                                  Vector_t nr,
frey_m's avatar
frey_m committed
42 43
                                  Vector_t hr,
                                  std::string interpl)
44
    : IrregularDomain(nr, hr, interpl)
frey_m's avatar
frey_m committed
45
{
46
    bgeom_m  = bgeom;
47 48 49 50

    setRangeMin(bgeom->getmincoords());
    setRangeMax(bgeom->getmaxcoords());

51 52
    bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m);
    if (have_inside_pt == false) {
53 54 55
        throw OpalException(
            "ArbitraryDomain::ArbitraryDomain()",
            "No point inside geometry found/set!");
56
    }
57 58

    startId = 0;
gsell's avatar
gsell committed
59 60
}

61 62
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
63 64
}

65
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
66

67 68
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

69
    setHr(hr);
70

71
    int zGhostOffsetLeft  = (localId[2].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
72
    int zGhostOffsetRight = (localId[2].last() == nr_m[2] - 1) ? 0 : 1;
73
    int yGhostOffsetLeft  = (localId[1].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
74
    int yGhostOffsetRight = (localId[1].last() == nr_m[1] - 1) ? 0 : 1;
75
    int xGhostOffsetLeft  = (localId[0].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
76
    int xGhostOffsetRight = (localId[0].last() == nr_m[0] - 1) ? 0 : 1;
kraus's avatar
kraus committed
77

78 79 80 81 82 83 84 85 86
    hasGeometryChanged_m = true;

    IntersectLoX.clear();
    IntersectHiX.clear();
    IntersectLoY.clear();
    IntersectHiY.clear();
    IntersectLoZ.clear();
    IntersectHiZ.clear();

87
    // Calculate intersection
88 89 90
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
91

92
    // We cannot assume that the geometry is symmetric about the xy, xz, and yz planes!
93 94
    // In my spiral inflector simulation, this is not the case for z direction for
    // example (-0.13 to +0.025). -DW
95
    for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
96

frey_m's avatar
frey_m committed
97
        //saveP_old[2] = (idz - (nr_m[2]-1)/2.0)*hr[2];
98
        P[2] = min_m[2] + (idz + 0.5) * hr[2];
99

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
100
        for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
101

frey_m's avatar
frey_m committed
102
            //saveP_old[1] = (idy - (nr_m[1]-1)/2.0)*hr[1];
103
            P[1] = min_m[1] + (idy + 0.5) * hr[1];
104

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
105
            for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
106

frey_m's avatar
frey_m committed
107
                //saveP_old[0] = (idx - (nr_m[0]-1)/2.0)*hr[0];
108
                P[0] = min_m[0] + (idx + 0.5) * hr[0];
109

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
110
                // *gmsg << "Now working on point " << saveP << " (original was " << saveP_old << ")" << endl;
111

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
112
                //P = saveP;
113

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
114
                if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
kraus's avatar
kraus committed
115

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
116 117 118
                    // Fill the map with true or false values for very fast isInside tests
                    // during the rest of the fieldsolve.
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = true;
119

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
120 121 122 123 124 125
                    // Replace the old reference point with the new point (which we know is
                    // inside because we just tested for it. This makes the algorithm faster
                    // because fastIsInside() creates a number of segments that depends on the
                    // distance between P and P0. Using the previous P as the new P0
                    // assures the smallest possible distance in most cases. -DW
                    P0 = P;
126

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
127
                    std::tuple<int, int, int> pos(idx, idy, idz);
128

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
129
                    dir = Vector_t(0, 0, 1);
130

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
131
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
132
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
133
                    } else {
134
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
135 136
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
137
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
138
                    }
139

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
140
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
141
                        IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
142
                    } else {
143
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
144 145
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
146
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
147
                    }
kraus's avatar
kraus committed
148

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
149
                    dir = Vector_t(0, 1, 0);
150

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
151 152 153
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
                        IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
                    } else {
154
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
155 156
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
157
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
158
                    }
159

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
160
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
161
                        IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
162
                    } else {
163
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
164 165
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
166
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
167
                    }
168

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
169
                    dir = Vector_t(1, 0, 0);
170

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
171
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
172
                        IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
173
                    } else {
174
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
175 176
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
177
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
178
                    }
179

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
180
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
181
                        IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
182
                    } else {
183
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
184 185
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
186
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
187 188 189
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
190
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
191 192
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
193
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
194 195 196 197
                }
            }
        }
    }
198

199 200
    INFOMSG(level2 << "* Finding number of ghost nodes to the left..." << endl);

201 202 203
    //number of ghost nodes to the left
    int numGhostNodesLeft = 0;
    if(localId[2].first() != 0) {
frey_m's avatar
frey_m committed
204 205
        for(int idx = 0; idx < nr_m[0]; idx++) {
            for(int idy = 0; idy < nr_m[1]; idy++) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
206
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
207 208 209 210 211
                    numGhostNodesLeft++;
            }
        }
    }

212 213
    INFOMSG(level2 << "* Finding number of xy points in each plane along z..." << endl);

214 215 216 217
    //xy points in z plane
    int numtotal = 0;
    numXY.clear();
    for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
snuverink_j's avatar
snuverink_j committed
218
        int numxy = 0;
frey_m's avatar
frey_m committed
219 220
        for(int idx = 0; idx < nr_m[0]; idx++) {
            for(int idy = 0; idy < nr_m[1]; idy++) {
221 222 223 224 225 226 227 228 229
                if(isInside(idx, idy, idz))
                    numxy++;
            }
        }
        numXY[idz-localId[2].first()] = numxy;
        numtotal += numxy;
    }

    int startIdx = 0;
230
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
231
    startIdx -= numtotal;
232

233
    // Build up index and coord map
234 235
    idxMap_m.clear();
    coordMap_m.clear();
236
    int index = startIdx - numGhostNodesLeft;
237

238 239
    INFOMSG(level2 << "* Building up index and coordinate map..." << endl);

240
    for(int idz = localId[2].first() - zGhostOffsetLeft; idz <= localId[2].last() + zGhostOffsetRight; idz++) {
frey_m's avatar
frey_m committed
241 242
        for(int idy = 0; idy < nr_m[1]; idy++) {
            for(int idx = 0; idx < nr_m[0]; idx++) {
kraus's avatar
kraus committed
243
                if(isInside(idx, idy, idz)) {
244 245
                    idxMap_m[toCoordIdx(idx, idy, idz)] = index;
                    coordMap_m[index] = toCoordIdx(idx, idy, idz);
246 247
                    index++;
                }
248
            }
249 250
        }
    }
251 252

    INFOMSG(level2 << "* Done." << endl);
253
}
254

255 256
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
frey_m's avatar
frey_m committed
257
    return (idz * nr_m[1] + idy) * nr_m[0]  + idx;
gsell's avatar
gsell committed
258 259
}

260 261
// Conversion from (x,y,z) to index on the 3D grid
int ArbitraryDomain::getIdx(int idx, int idy, int idz) {
262

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
263
    if(isInside(idx, idy, idz) && idx >= 0 && idy >= 0 && idz >= 0)
264
        return idxMap_m[toCoordIdx(idx, idy, idz)];
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
265 266
    else
        return -1;
kraus's avatar
kraus committed
267
}
gsell's avatar
gsell committed
268

269 270
// Conversion from a 3D index to (x,y,z)
inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) {
271
    int id = coordMap_m[idxyz];
frey_m's avatar
frey_m committed
272 273 274 275
    idx = id % (int)nr_m[0];
    id /= nr_m[0];
    idy = id % (int)nr_m[1];
    id /= nr_m[1];
276
    idz = id;
gsell's avatar
gsell committed
277 278
}

279
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
280

281 282 283 284
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
285 286
  inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
  Vector_t P;
287

frey_m's avatar
frey_m committed
288 289 290
  P[0] = min_m[0] + (idx + 0.5) * hr_m[0];
  P[1] = min_m[1] + (idy + 0.5) * hr_m[1];
  P[2] = min_m[2] + (idz + 0.5) * hr_m[2];
291

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
292 293
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
294 295 296
*/

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
297 298 299
  inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
  Vector_t P;

frey_m's avatar
frey_m committed
300 301 302
  P[0] = (idx - (nr_m[0]-1)/2.0) * hr_m[0];
  P[1] = (idy - (nr_m[1]-1)/2.0) * hr_m[1];
  P[2] = (idz - (nr_m[2]-1)/2.0) * hr_m[2];
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337

  bool ret = false;
  int  countH, countL;
  std::multimap < std::tuple<int, int, int>, double >::iterator itrH, itrL;
  std::tuple<int, int, int> coordxyz(idx, idy, idz);

  //check if z is inside with x,y coords
  itrH = IntersectHiZ.find(coordxyz);
  itrL = IntersectLoZ.find(coordxyz);

  countH = IntersectHiZ.count(coordxyz);
  countL = IntersectLoZ.count(coordxyz);
  if(countH == 1 && countL == 1)
  ret = (P[2] <= itrH->second) && (P[2] >= itrL->second);

  //check if y is inside with x,z coords
  itrH = IntersectHiY.find(coordxyz);
  itrL = IntersectLoY.find(coordxyz);

  countH = IntersectHiY.count(coordxyz);
  countL = IntersectLoY.count(coordxyz);
  if(countH == 1 && countL == 1)
  ret = ret && (P[1] <= itrH->second) && (P[1] >= itrL->second);

  //check if x is inside with y,z coords
  itrH = IntersectHiX.find(coordxyz);
  itrL = IntersectLoX.find(coordxyz);

  countH = IntersectHiX.count(coordxyz);
  countL = IntersectLoX.count(coordxyz);
  if(countH == 1 && countL == 1)
  ret = ret && (P[0] <= itrH->second) && (P[0] >= itrL->second);

  return ret;
  }
338
*/
gsell's avatar
gsell committed
339 340

int ArbitraryDomain::getNumXY(int z) {
kraus's avatar
kraus committed
341

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
342
    return numXY[z];
343
}
gsell's avatar
gsell committed
344

345 346
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
                                            StencilValue_t& value, double& /*scaleFactor*/)
frey_m's avatar
frey_m committed
347
{
frey_m's avatar
frey_m committed
348 349 350 351 352 353 354
    value.west = -1/(hr_m[0]*hr_m[0]);
    value.east = -1/(hr_m[0]*hr_m[0]);
    value.north = -1/(hr_m[1]*hr_m[1]);
    value.south = -1/(hr_m[1]*hr_m[1]);
    value.front = -1/(hr_m[2]*hr_m[2]);
    value.back = -1/(hr_m[2]*hr_m[2]);
    value.center = 2/(hr_m[0]*hr_m[0]) + 2/(hr_m[1]*hr_m[1]) + 2/(hr_m[2]*hr_m[2]);
355

356
    if(!isInside(idx-1,idy,idz))
357
        value.west = 0.0;
kraus's avatar
kraus committed
358
    if(!isInside(idx+1,idy,idz))
359
        value.east = 0.0;
360

361
    if(!isInside(idx,idy+1,idz))
362
        value.north = 0.0;
kraus's avatar
kraus committed
363
    if(!isInside(idx,idy-1,idz))
364
        value.south = 0.0;
kraus's avatar
kraus committed
365 366

    if(!isInside(idx,idy,idz-1))
367
        value.front = 0.0;
kraus's avatar
kraus committed
368
    if(!isInside(idx,idy,idz+1))
369
        value.back = 0.0;
370
}
371

372 373
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
                                          StencilValue_t& value, double &scaleFactor)
374
{
375
    scaleFactor = 1;
376

frey_m's avatar
frey_m committed
377 378 379
    double cx = (idx - (nr_m[0]-1)/2.0)*hr_m[0];
    double cy = (idy - (nr_m[1]-1)/2.0)*hr_m[1];
    double cz = (idz - (nr_m[2]-1)/2.0)*hr_m[2];
380

frey_m's avatar
frey_m committed
381 382 383 384 385 386
    double dx_w=hr_m[0];
    double dx_e=hr_m[0];
    double dy_n=hr_m[1];
    double dy_s=hr_m[1];
    double dz_f=hr_m[2];
    double dz_b=hr_m[2];
387
    value.center = 0.0;
388

389
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
390

frey_m's avatar
frey_m committed
391
    if (idx == nr_m[0]-1)
392
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
393
    if (idx == 0)
394
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
frey_m's avatar
frey_m committed
395
    if (idy == nr_m[1]-1)
396
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
397
    if (idy == 0)
398
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
frey_m's avatar
frey_m committed
399
    if (idz == nr_m[2]-1)
400
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
401
    if (idz == 0)
402
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
403 404

    if(dx_w != 0)
405
        value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
406
    else
407
        value.west = 0;
408
    if(dx_e != 0)
409
        value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
410
    else
411
        value.east = 0;
412
    if(dy_n != 0)
413
        value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
414
    else
415
        value.north = 0;
416
    if(dy_s != 0)
417
        value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
418
    else
419
        value.south = 0;
420
    if(dz_f != 0)
421
        value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
422
    else
423
        value.front = 0;
424
    if(dz_b != 0)
425
        value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
426
    else
427
        value.back = 0;
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451

    //RHS scaleFactor for current 3D index
    //0.5* comes from discretiztaion
    //scaleFactor = 0.5*(dw+de)*(dn+ds)*(df+db);
    scaleFactor = 0.5;
    if(dx_w + dx_e != 0)
        scaleFactor *= (dx_w + dx_e);
    if(dy_n + dy_s != 0)
        scaleFactor *= (dy_n + dy_s);
    if(dz_f + dz_b != 0)
        scaleFactor *= (dz_f + dz_b);

    //catch the case where a point lies on the boundary
    double m1 = dx_w * dx_e;
    double m2 = dy_n * dy_s;
    if(dx_e == 0)
        m1 = dx_w;
    if(dx_w == 0)
        m1 = dx_e;
    if(dy_n == 0)
        m2 = dy_s;
    if(dy_s == 0)
        m2 = dy_n;

frey_m's avatar
frey_m committed
452
    value.center = 2 / hr_m[2];
453
    if(dx_w != 0 || dx_e != 0)
454
        value.center *= (dx_w + dx_e);
455
    if(dy_n != 0 || dy_s != 0)
456
        value.center *= (dy_n + dy_s);
457
    if(dx_w != 0 || dx_e != 0)
458
        value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
459
    if(dy_n != 0 || dy_s != 0)
460
        value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
461
}
gsell's avatar
gsell committed
462

463 464 465 466 467 468 469
// vi: set et ts=4 sw=4 sts=4:
// Local Variables:
// mode:c
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End: