ArbitraryDomain.cpp 21.1 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

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

    // stencil center value has to be positive!
    assert(C > 0);
}

frey_m's avatar
frey_m committed
408 409 410 411
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz, double& W,
                                            double& E, double& S, double& N, double& F,
                                            double& B, double& C, double& /*scaleFactor*/)
{
412 413 414 415 416 417 418 419
    W = -1/(hr[0]*hr[0]);
    E = -1/(hr[0]*hr[0]);
    N = -1/(hr[1]*hr[1]);
    S = -1/(hr[1]*hr[1]);
    F = -1/(hr[2]*hr[2]);
    B = -1/(hr[2]*hr[2]);
    C = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);

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

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

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

frey_m's avatar
frey_m committed
436 437 438
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, double& W,
                                          double& E, double& S, double& N, double& F,
                                          double& B, double& C, double &scaleFactor)
439
{
440
    scaleFactor = 1;
441

442 443 444
    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];
445

446
    double dx_w=hr[0];
447 448
    double dx_e=hr[0];
    double dy_n=hr[1];
449 450 451
    double dy_s=hr[1];
    double dz_f=hr[2];
    double dz_b=hr[2];
452 453
    C = 0.0;

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

456
    if (idx == nr[0]-1)
457
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
458
    if (idx == 0)
459
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
460
    if (idy == nr[1]-1)
461
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
462
    if (idy == 0)
463
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
464
    if (idz == nr[2]-1)
465
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
466
    if (idz == 0)
467
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525

    if(dx_w != 0)
        W = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
    else
        W = 0;
    if(dx_e != 0)
        E = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
    else
        E = 0;
    if(dy_n != 0)
        N = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
    else
        N = 0;
    if(dy_s != 0)
        S = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
    else
        S = 0;
    if(dz_f != 0)
        F = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
    else
        F = 0;
    if(dz_b != 0)
        B = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
    else
        B = 0;

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

    C = 2 / hr[2];
    if(dx_w != 0 || dx_e != 0)
        C *= (dx_w + dx_e);
    if(dy_n != 0 || dy_s != 0)
        C *= (dy_n + dy_s);
    if(dx_w != 0 || dx_e != 0)
        C += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
    if(dy_n != 0 || dy_s != 0)
        C += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
526
}
gsell's avatar
gsell committed
527

frey_m's avatar
frey_m committed
528 529 530
void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W,
                                    int &E, int &S, int &N, int &F, int &B)
{
531 532 533 534 535 536
    W = getIdx(idx - 1, idy, idz);
    E = getIdx(idx + 1, idy, idz);
    N = getIdx(idx, idy + 1, idz);
    S = getIdx(idx, idy - 1, idz);
    F = getIdx(idx, idy, idz - 1);
    B = getIdx(idx, idy, idz + 1);
537

kraus's avatar
kraus committed
538
    if(!isInside(idx+1,idy,idz))
gsell's avatar
gsell committed
539 540
        E = -1;

541 542 543 544
    if(!isInside(idx-1,idy,idz))
        W = -1;

    if(!isInside(idx,idy+1,idz))
gsell's avatar
gsell committed
545
        N = -1;
546

kraus's avatar
kraus committed
547
    if(!isInside(idx,idy-1,idz))
gsell's avatar
gsell committed
548 549
        S = -1;

kraus's avatar
kraus committed
550
    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
551
        F = -1;
kraus's avatar
kraus committed
552 553

    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
554
        B = -1;
gsell's avatar
gsell committed
555 556 557 558 559 560 561 562 563 564

}


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

565
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
566 567 568 569 570
    // 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
571 572 573

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
574
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
575 576 577
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
581 582 583 584
    v(0) = (quaternion(0) * quaternion(0)
            + quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
kraus's avatar
kraus committed
585

586 587 588 589
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
595 596 597 598
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
599 600 601 602

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

603
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
604 605
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
606
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
607

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
608 609 610 611
    v(2) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            + quaternion(3) * quaternion(3));
612 613

}
614 615 616 617 618 619 620 621

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