ArbitraryDomain.cpp 20.4 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,
42
                                  Vector_t /*hr*/,
kraus's avatar
kraus committed
43
                                  std::string interpl){
44
    bgeom_m  = bgeom;
45 46 47 48 49

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

    geomCentroid_m = (min_m + max_m)/2.0;
50

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
    setNr(nr);
58
    for(int i=0; i<3; i++)
59
        Geo_hr_m[i] = (max_m[i] - min_m[i])/nr[i];
60
    setHr(Geo_hr_m);
61 62 63 64 65 66 67 68 69

    startId = 0;

    if (interpl == "CONSTANT")
        interpolationMethod = CONSTANT;
    else if(interpl == "LINEAR")
        interpolationMethod = LINEAR;
    else if(interpl == "QUADRATIC")
        interpolationMethod = QUADRATIC;
gsell's avatar
gsell committed
70 71
}

72 73
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
74 75
}

76
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
77

78 79
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

80
    setHr(hr);
81

82 83 84 85
    globalMeanR_m = getGlobalMeanR();

    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
86
    for (int i=1; i<4; i++)
87
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
88 89 90 91 92 93 94

    int zGhostOffsetLeft  = (localId[2].first()== 0) ? 0 : 1;
    int zGhostOffsetRight = (localId[2].last() == nr[2] - 1) ? 0 : 1;
    int yGhostOffsetLeft  = (localId[1].first()== 0) ? 0 : 1;
    int yGhostOffsetRight = (localId[1].last() == nr[1] - 1) ? 0 : 1;
    int xGhostOffsetLeft  = (localId[0].first()== 0) ? 0 : 1;
    int xGhostOffsetRight = (localId[0].last() == nr[0] - 1) ? 0 : 1;
kraus's avatar
kraus committed
95

96 97 98 99 100 101 102 103 104
    hasGeometryChanged_m = true;

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

105
    // Calculate intersection
106 107 108
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
109

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
115
        //saveP_old[2] = (idz - (nr[2]-1)/2.0)*hr[2];
116
        P[2] = min_m[2] + (idz + 0.5) * hr[2];
117

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
120
            //saveP_old[1] = (idy - (nr[1]-1)/2.0)*hr[1];
121
            P[1] = min_m[1] + (idy + 0.5) * hr[1];
122

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
125
                //saveP_old[0] = (idx - (nr[0]-1)/2.0)*hr[0];
126
                P[0] = min_m[0] + (idx + 0.5) * hr[0];
127

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
130
                //P = saveP;
131

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
132 133 134
                //rotateWithQuaternion(P, localToGlobalQuaternion_m);
                //P += geomCentroid_m; //sorry, this doesn't make sense. -DW
                //P += globalMeanR_m;
135

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
138 139 140
                    // 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;
141

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
142 143 144 145 146 147
                    // 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;
148

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
151 152
                    //rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 0, 1);
153

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
154
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
155 156
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
157
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
kraus's avatar
kraus committed
158
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
159
                    } else {
160
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
161 162
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
163
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
164
                    }
165

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
166
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
167 168 169 170
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
171
                    } else {
172
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
173 174
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
175
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
176
                    }
kraus's avatar
kraus committed
177

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
178 179
                    //rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 1, 0);
180

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
181 182 183 184 185 186
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
                    } else {
187
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
188 189
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
190
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
191
                    }
192

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
193
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
194 195 196 197
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
198
                    } else {
199
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
200 201
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
202
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
203
                    }
204

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
205 206
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
207

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
208
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
209 210 211 212
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
213
                    } else {
214
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
215 216
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
217
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
218
                    }
219

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
220
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
221 222 223 224
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
225
                    } else {
226
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
227 228
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
229
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
230 231 232
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
233
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
234 235
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
236
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
237 238 239 240
                }
            }
        }
    }
241

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

244 245 246 247 248
    //number of ghost nodes to the left
    int numGhostNodesLeft = 0;
    if(localId[2].first() != 0) {
        for(int idx = 0; idx < nr[0]; idx++) {
            for(int idy = 0; idy < nr[1]; idy++) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
249
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
250 251 252 253 254
                    numGhostNodesLeft++;
            }
        }
    }

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

257 258 259 260
    //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
261
        int numxy = 0;
262 263 264 265 266 267 268 269 270 271 272
        for(int idx = 0; idx < nr[0]; idx++) {
            for(int idy = 0; idy < nr[1]; idy++) {
                if(isInside(idx, idy, idz))
                    numxy++;
            }
        }
        numXY[idz-localId[2].first()] = numxy;
        numtotal += numxy;
    }

    int startIdx = 0;
273
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
274
    startIdx -= numtotal;
275

276
    // Build up index and coord map
277 278
    IdxMap.clear();
    CoordMap.clear();
279
    int index = startIdx - numGhostNodesLeft;
280

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

283 284 285
    for(int idz = localId[2].first() - zGhostOffsetLeft; idz <= localId[2].last() + zGhostOffsetRight; idz++) {
        for(int idy = 0; idy < nr[1]; idy++) {
            for(int idx = 0; idx < nr[0]; idx++) {
kraus's avatar
kraus committed
286
                if(isInside(idx, idy, idz)) {
287 288 289 290
                    IdxMap[toCoordIdx(idx, idy, idz)] = index;
                    CoordMap[index] = toCoordIdx(idx, idy, idz);
                    index++;
                }
291
            }
292 293
        }
    }
294 295

    INFOMSG(level2 << "* Done." << endl);
296
}
297

298 299
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
300
    return (idz * nr[1] + idy) * nr[0]  + idx;
gsell's avatar
gsell committed
301 302
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
306 307 308 309
    if(isInside(idx, idy, idz) && idx >= 0 && idy >= 0 && idz >= 0)
        return IdxMap[toCoordIdx(idx, idy, idz)];
    else
        return -1;
kraus's avatar
kraus committed
310
}
gsell's avatar
gsell committed
311

312 313 314 315 316 317 318 319
// Conversion from a 3D index to (x,y,z)
inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) {
    int id = CoordMap[idxyz];
    idx = id % (int)nr[0];
    id /= nr[0];
    idy = id % (int)nr[1];
    id /= nr[1];
    idz = id;
gsell's avatar
gsell committed
320 321
}

322
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
323

324 325 326 327
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

331 332 333
  P[0] = min_m[0] + (idx + 0.5) * hr[0];
  P[1] = min_m[1] + (idy + 0.5) * hr[1];
  P[2] = min_m[2] + (idz + 0.5) * hr[2];
334

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
335 336
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
337 338 339
*/

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380
  inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
  Vector_t P;

  P[0] = (idx - (nr[0]-1)/2.0) * hr[0];
  P[1] = (idy - (nr[1]-1)/2.0) * hr[1];
  P[2] = (idz - (nr[2]-1)/2.0) * hr[2];

  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;
  }
381
*/
gsell's avatar
gsell committed
382 383

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
385
    return numXY[z];
386
}
gsell's avatar
gsell committed
387

388 389
void ArbitraryDomain::getBoundaryStencil(int idx, int idy, int idz,
                                         StencilValue_t& value, double &scaleFactor)
frey_m's avatar
frey_m committed
390
{
391
    scaleFactor = 1.0;
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
392
    // determine which interpolation method we use for points near the boundary
393
    switch(interpolationMethod){
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
394
    case CONSTANT:
395
        constantInterpolation(idx,idy,idz,value,scaleFactor);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
396 397
        break;
    case LINEAR:
398
        linearInterpolation(idx,idy,idz,value,scaleFactor);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
399 400
        break;
    case QUADRATIC:
401
        //  QuadraticInterpolation(idx,idy,idz,value,scaleFactor);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
402
        break;
403 404 405
    }

    // stencil center value has to be positive!
406
    assert(value.center > 0);
407 408
}

409 410
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
                                            StencilValue_t& value, double& /*scaleFactor*/)
frey_m's avatar
frey_m committed
411
{
412 413 414 415 416 417 418
    value.west = -1/(hr[0]*hr[0]);
    value.east = -1/(hr[0]*hr[0]);
    value.north = -1/(hr[1]*hr[1]);
    value.south = -1/(hr[1]*hr[1]);
    value.front = -1/(hr[2]*hr[2]);
    value.back = -1/(hr[2]*hr[2]);
    value.center = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
419

420
    if(!isInside(idx-1,idy,idz))
421
        value.west = 0.0;
kraus's avatar
kraus committed
422
    if(!isInside(idx+1,idy,idz))
423
        value.east = 0.0;
424

425
    if(!isInside(idx,idy+1,idz))
426
        value.north = 0.0;
kraus's avatar
kraus committed
427
    if(!isInside(idx,idy-1,idz))
428
        value.south = 0.0;
kraus's avatar
kraus committed
429 430

    if(!isInside(idx,idy,idz-1))
431
        value.front = 0.0;
kraus's avatar
kraus committed
432
    if(!isInside(idx,idy,idz+1))
433
        value.back = 0.0;
434
}
435

436 437
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
                                          StencilValue_t& value, double &scaleFactor)
438
{
439
    scaleFactor = 1;
440

441 442 443
    double cx = (idx - (nr[0]-1)/2.0)*hr[0];
    double cy = (idy - (nr[1]-1)/2.0)*hr[1];
    double cz = (idz - (nr[2]-1)/2.0)*hr[2];
444

445
    double dx_w=hr[0];
446 447
    double dx_e=hr[0];
    double dy_n=hr[1];
448 449 450
    double dy_s=hr[1];
    double dz_f=hr[2];
    double dz_b=hr[2];
451
    value.center = 0.0;
452

453
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
454

455
    if (idx == nr[0]-1)
456
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
457
    if (idx == 0)
458
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
459
    if (idy == nr[1]-1)
460
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
461
    if (idy == 0)
462
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
463
    if (idz == nr[2]-1)
464
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
465
    if (idz == 0)
466
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
467 468

    if(dx_w != 0)
469
        value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
470
    else
471
        value.west = 0;
472
    if(dx_e != 0)
473
        value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
474
    else
475
        value.east = 0;
476
    if(dy_n != 0)
477
        value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
478
    else
479
        value.north = 0;
480
    if(dy_s != 0)
481
        value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
482
    else
483
        value.south = 0;
484
    if(dz_f != 0)
485
        value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
486
    else
487
        value.front = 0;
488
    if(dz_b != 0)
489
        value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
490
    else
491
        value.back = 0;
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515

    //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;

516
    value.center = 2 / hr[2];
517
    if(dx_w != 0 || dx_e != 0)
518
        value.center *= (dx_w + dx_e);
519
    if(dy_n != 0 || dy_s != 0)
520
        value.center *= (dy_n + dy_s);
521
    if(dx_w != 0 || dx_e != 0)
522
        value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
523
    if(dy_n != 0 || dy_s != 0)
524
        value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
525
}
gsell's avatar
gsell committed
526 527 528 529 530 531 532

inline void ArbitraryDomain::crossProduct(double A[], double B[], double C[]) {
    C[0] = A[1] * B[2] - A[2] * B[1];
    C[1] = A[2] * B[0] - A[0] * B[2];
    C[2] = A[0] * B[1] - A[1] * B[0];
}

533
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
534 535 536 537 538
    // rotates a Vector_t (3 elements) using a quaternion.
    // Flip direction of rotation by quaternionVectorcomponent *= -1

    Vector_t const quaternionVectorComponent = Vector_t(quaternion(1), quaternion(2), quaternion(3));
    double const quaternionScalarComponent = quaternion(0);
kraus's avatar
kraus committed
539 540 541

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
542
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
543 544 545
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

546
inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
547
    // rotates the positive xaxis using a quaternion.
kraus's avatar
kraus committed
548

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
549 550 551 552
    v(0) = (quaternion(0) * quaternion(0)
            + quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
kraus's avatar
kraus committed
553

554 555 556 557
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

558
inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
559
    // rotates the positive yaxis using a quaternion.
kraus's avatar
kraus committed
560

561
    v(0) = 2.0 * (quaternion(1) * quaternion(2) - quaternion(0) * quaternion(3));
kraus's avatar
kraus committed
562

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
563 564 565 566
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
567 568 569 570

    v(2) = 2.0 * (quaternion(2) * quaternion(3) + quaternion(0) * quaternion(1));
}

571
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
572 573
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
574
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
575

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
576 577 578 579
    v(2) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            + quaternion(3) * quaternion(3));
580 581

}
582 583 584 585 586 587 588 589

// 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: