ArbitraryDomain.cpp 20.2 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 51

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

    geomCentroid_m = (min_m + max_m)/2.0;
52

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

    startId = 0;
gsell's avatar
gsell committed
61 62
}

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

67
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
68

69 70
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

71
    setHr(hr);
72

73 74 75 76
    globalMeanR_m = getGlobalMeanR();

    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
77
    for (int i=1; i<4; i++)
78
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
79 80

    int zGhostOffsetLeft  = (localId[2].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
81
    int zGhostOffsetRight = (localId[2].last() == nr_m[2] - 1) ? 0 : 1;
82
    int yGhostOffsetLeft  = (localId[1].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
83
    int yGhostOffsetRight = (localId[1].last() == nr_m[1] - 1) ? 0 : 1;
84
    int xGhostOffsetLeft  = (localId[0].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
85
    int xGhostOffsetRight = (localId[0].last() == nr_m[0] - 1) ? 0 : 1;
kraus's avatar
kraus committed
86

87 88 89 90 91 92 93 94 95
    hasGeometryChanged_m = true;

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

96
    // Calculate intersection
97 98 99
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
100

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
121
                //P = saveP;
122

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
123 124 125
                //rotateWithQuaternion(P, localToGlobalQuaternion_m);
                //P += geomCentroid_m; //sorry, this doesn't make sense. -DW
                //P += globalMeanR_m;
126

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
129 130 131
                    // 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;
132

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
133 134 135 136 137 138
                    // 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;
139

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
142 143
                    //rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 0, 1);
144

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
145
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
146 147
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
148
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
kraus's avatar
kraus committed
149
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
150
                    } else {
151
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
152 153
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
154
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
155
                    }
156

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
157
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
158 159 160 161
                        //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
162
                    } else {
163
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
164 165
                        *gmsg << "zdir=-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
                    }
kraus's avatar
kraus committed
168

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
172 173 174 175 176 177
                    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 {
178
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
179 180
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
181
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
182
                    }
183

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
184
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
185 186 187 188
                        //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
189
                    } else {
190
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
191 192
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
193
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
194
                    }
195

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
196 197
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
198

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
199
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
200 201 202 203
                        //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
204
                    } else {
205
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
206 207
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
208
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
209
                    }
210

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
211
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
212 213 214 215
                        //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
216
                    } else {
217
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
218 219
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
220
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
221 222 223
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
224
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
225 226
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
227
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
228 229 230 231
                }
            }
        }
    }
232

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

235 236 237
    //number of ghost nodes to the left
    int numGhostNodesLeft = 0;
    if(localId[2].first() != 0) {
frey_m's avatar
frey_m committed
238 239
        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
240
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
241 242 243 244 245
                    numGhostNodesLeft++;
            }
        }
    }

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

248 249 250 251
    //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
252
        int numxy = 0;
frey_m's avatar
frey_m committed
253 254
        for(int idx = 0; idx < nr_m[0]; idx++) {
            for(int idy = 0; idy < nr_m[1]; idy++) {
255 256 257 258 259 260 261 262 263
                if(isInside(idx, idy, idz))
                    numxy++;
            }
        }
        numXY[idz-localId[2].first()] = numxy;
        numtotal += numxy;
    }

    int startIdx = 0;
264
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
265
    startIdx -= numtotal;
266

267
    // Build up index and coord map
268 269
    IdxMap.clear();
    CoordMap.clear();
270
    int index = startIdx - numGhostNodesLeft;
271

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

274
    for(int idz = localId[2].first() - zGhostOffsetLeft; idz <= localId[2].last() + zGhostOffsetRight; idz++) {
frey_m's avatar
frey_m committed
275 276
        for(int idy = 0; idy < nr_m[1]; idy++) {
            for(int idx = 0; idx < nr_m[0]; idx++) {
kraus's avatar
kraus committed
277
                if(isInside(idx, idy, idz)) {
278 279 280 281
                    IdxMap[toCoordIdx(idx, idy, idz)] = index;
                    CoordMap[index] = toCoordIdx(idx, idy, idz);
                    index++;
                }
282
            }
283 284
        }
    }
285 286

    INFOMSG(level2 << "* Done." << endl);
287
}
288

289 290
// 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
291
    return (idz * nr_m[1] + idy) * nr_m[0]  + idx;
gsell's avatar
gsell committed
292 293
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
297 298 299 300
    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
301
}
gsell's avatar
gsell committed
302

303 304 305
// 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];
frey_m's avatar
frey_m committed
306 307 308 309
    idx = id % (int)nr_m[0];
    id /= nr_m[0];
    idy = id % (int)nr_m[1];
    id /= nr_m[1];
310
    idz = id;
gsell's avatar
gsell committed
311 312
}

313
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
314

315 316 317 318
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

frey_m's avatar
frey_m committed
322 323 324
  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];
325

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
326 327
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
328 329 330
*/

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

frey_m's avatar
frey_m committed
334 335 336
  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
337 338 339 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

  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;
  }
372
*/
gsell's avatar
gsell committed
373 374

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
376
    return numXY[z];
377
}
gsell's avatar
gsell committed
378

379 380
void ArbitraryDomain::getBoundaryStencil(int idx, int idy, int idz,
                                         StencilValue_t& value, double &scaleFactor)
frey_m's avatar
frey_m committed
381
{
382
    scaleFactor = 1.0;
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
383
    // determine which interpolation method we use for points near the boundary
384
    switch(interpolationMethod_m){
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
385
    case CONSTANT:
386
        constantInterpolation(idx,idy,idz,value,scaleFactor);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
387 388
        break;
    case LINEAR:
389
        linearInterpolation(idx,idy,idz,value,scaleFactor);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
390 391
        break;
    case QUADRATIC:
392
        //  QuadraticInterpolation(idx,idy,idz,value,scaleFactor);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
393
        break;
394 395 396
    }

    // stencil center value has to be positive!
397
    assert(value.center > 0);
398 399
}

400 401
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
                                            StencilValue_t& value, double& /*scaleFactor*/)
frey_m's avatar
frey_m committed
402
{
frey_m's avatar
frey_m committed
403 404 405 406 407 408 409
    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]);
410

411
    if(!isInside(idx-1,idy,idz))
412
        value.west = 0.0;
kraus's avatar
kraus committed
413
    if(!isInside(idx+1,idy,idz))
414
        value.east = 0.0;
415

416
    if(!isInside(idx,idy+1,idz))
417
        value.north = 0.0;
kraus's avatar
kraus committed
418
    if(!isInside(idx,idy-1,idz))
419
        value.south = 0.0;
kraus's avatar
kraus committed
420 421

    if(!isInside(idx,idy,idz-1))
422
        value.front = 0.0;
kraus's avatar
kraus committed
423
    if(!isInside(idx,idy,idz+1))
424
        value.back = 0.0;
425
}
426

427 428
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
                                          StencilValue_t& value, double &scaleFactor)
429
{
430
    scaleFactor = 1;
431

frey_m's avatar
frey_m committed
432 433 434
    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];
435

frey_m's avatar
frey_m committed
436 437 438 439 440 441
    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];
442
    value.center = 0.0;
443

444
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
445

frey_m's avatar
frey_m committed
446
    if (idx == nr_m[0]-1)
447
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
448
    if (idx == 0)
449
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
frey_m's avatar
frey_m committed
450
    if (idy == nr_m[1]-1)
451
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
452
    if (idy == 0)
453
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
frey_m's avatar
frey_m committed
454
    if (idz == nr_m[2]-1)
455
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
456
    if (idz == 0)
457
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
458 459

    if(dx_w != 0)
460
        value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
461
    else
462
        value.west = 0;
463
    if(dx_e != 0)
464
        value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
465
    else
466
        value.east = 0;
467
    if(dy_n != 0)
468
        value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
469
    else
470
        value.north = 0;
471
    if(dy_s != 0)
472
        value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
473
    else
474
        value.south = 0;
475
    if(dz_f != 0)
476
        value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
477
    else
478
        value.front = 0;
479
    if(dz_b != 0)
480
        value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
481
    else
482
        value.back = 0;
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506

    //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
507
    value.center = 2 / hr_m[2];
508
    if(dx_w != 0 || dx_e != 0)
509
        value.center *= (dx_w + dx_e);
510
    if(dy_n != 0 || dy_s != 0)
511
        value.center *= (dy_n + dy_s);
512
    if(dx_w != 0 || dx_e != 0)
513
        value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
514
    if(dy_n != 0 || dy_s != 0)
515
        value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
516
}
gsell's avatar
gsell committed
517 518 519 520 521 522 523

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];
}

524
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
525 526 527 528 529
    // 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
530 531 532

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
533
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
534 535 536
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
540 541 542 543
    v(0) = (quaternion(0) * quaternion(0)
            + quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
kraus's avatar
kraus committed
544

545 546 547 548
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
554 555 556 557
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
558 559 560 561

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

562
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
563 564
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
565
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
566

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

}
573 574 575 576 577 578 579 580

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