ArbitraryDomain.cpp 20.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,
42
                                  Vector_t /*hr*/,
kraus's avatar
kraus committed
43
                                  std::string interpl){
44 45 46
    bgeom_m  = bgeom;
    minCoords_m = bgeom->getmincoords();
    maxCoords_m = bgeom->getmaxcoords();
47
    geomCentroid_m = (minCoords_m + maxCoords_m)/2.0;
48

49 50
    bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m);
    if (have_inside_pt == false) {
51 52 53
        throw OpalException(
            "ArbitraryDomain::ArbitraryDomain()",
            "No point inside geometry found/set!");
54
    }
55
    setNr(nr);
56
    for(int i=0; i<3; i++)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
57
        Geo_hr_m[i] = (maxCoords_m[i] - minCoords_m[i])/nr[i];
58
    setHr(Geo_hr_m);
59 60 61 62 63 64 65 66 67

    startId = 0;

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

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

74
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
75

76 77
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

78
    setHr(hr);
79

80 81 82 83
    globalMeanR_m = getGlobalMeanR();

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

    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
93

94 95 96 97 98 99 100 101 102
    hasGeometryChanged_m = true;

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

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

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
128
                //P = saveP;
129

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
164
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
165 166 167 168
                        //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
169
                    } else {
170
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
171 172
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
173
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
174
                    }
kraus's avatar
kraus committed
175

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
179 180 181 182 183 184
                    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 {
185
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
186 187
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
188
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
189
                    }
190

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
191
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
192 193 194 195
                        //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
196
                    } else {
197
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
198 199
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
200
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
201
                    }
202

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
206
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
207 208 209 210
                        //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
211
                    } else {
212
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
213 214
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
215
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
216
                    }
217

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
218
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
219 220 221 222
                        //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
223
                    } else {
224
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
225 226
                        *gmsg << "xdir=-1 " << -dir << " 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
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
231
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
232 233
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
234
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
235 236 237 238
                }
            }
        }
    }
239

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

242 243 244 245 246
    //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
247
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
248 249 250 251 252
                    numGhostNodesLeft++;
            }
        }
    }

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

255 256 257 258
    //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
259
        int numxy = 0;
260 261 262 263 264 265 266 267 268 269 270
        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;
271
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
272
    startIdx -= numtotal;
273

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

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

281 282 283
    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
284
                if(isInside(idx, idy, idz)) {
285 286 287 288
                    IdxMap[toCoordIdx(idx, idy, idz)] = index;
                    CoordMap[index] = toCoordIdx(idx, idy, idz);
                    index++;
                }
289
            }
290 291
        }
    }
292 293

    INFOMSG(level2 << "* Done." << endl);
294
}
295

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
304 305 306 307
    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
308
}
gsell's avatar
gsell committed
309

310 311 312 313 314 315 316 317
// 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
318 319
}

320
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
321

322 323 324 325
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
329 330 331
  P[0] = minCoords_m[0] + (idx + 0.5) * hr[0];
  P[1] = minCoords_m[1] + (idy + 0.5) * hr[1];
  P[2] = minCoords_m[2] + (idz + 0.5) * hr[2];
332

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

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
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 372 373 374 375 376 377 378
  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;
  }
379
*/
gsell's avatar
gsell committed
380 381

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

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

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

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

407 408
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
                                            StencilValue_t& value, double& /*scaleFactor*/)
frey_m's avatar
frey_m committed
409
{
410 411 412 413 414 415 416
    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]);
417

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

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

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

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

439 440 441
    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];
442

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

451
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
452

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

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

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

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

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

531
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
532 533 534 535 536
    // 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
537 538 539

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

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

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

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

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

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

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

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

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

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

}
580 581 582 583 584 585 586 587

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