ArbitraryDomain.cpp 19.6 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::constantInterpolation(int idx, int idy, int idz,
                                            StencilValue_t& value, double& /*scaleFactor*/)
frey_m's avatar
frey_m committed
381
{
frey_m's avatar
frey_m committed
382 383 384 385 386 387 388
    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]);
389

390
    if(!isInside(idx-1,idy,idz))
391
        value.west = 0.0;
kraus's avatar
kraus committed
392
    if(!isInside(idx+1,idy,idz))
393
        value.east = 0.0;
394

395
    if(!isInside(idx,idy+1,idz))
396
        value.north = 0.0;
kraus's avatar
kraus committed
397
    if(!isInside(idx,idy-1,idz))
398
        value.south = 0.0;
kraus's avatar
kraus committed
399 400

    if(!isInside(idx,idy,idz-1))
401
        value.front = 0.0;
kraus's avatar
kraus committed
402
    if(!isInside(idx,idy,idz+1))
403
        value.back = 0.0;
404
}
405

406 407
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
                                          StencilValue_t& value, double &scaleFactor)
408
{
409
    scaleFactor = 1;
410

frey_m's avatar
frey_m committed
411 412 413
    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];
414

frey_m's avatar
frey_m committed
415 416 417 418 419 420
    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];
421
    value.center = 0.0;
422

423
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
424

frey_m's avatar
frey_m committed
425
    if (idx == nr_m[0]-1)
426
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
427
    if (idx == 0)
428
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
frey_m's avatar
frey_m committed
429
    if (idy == nr_m[1]-1)
430
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
431
    if (idy == 0)
432
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
frey_m's avatar
frey_m committed
433
    if (idz == nr_m[2]-1)
434
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
435
    if (idz == 0)
436
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
437 438

    if(dx_w != 0)
439
        value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
440
    else
441
        value.west = 0;
442
    if(dx_e != 0)
443
        value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
444
    else
445
        value.east = 0;
446
    if(dy_n != 0)
447
        value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
448
    else
449
        value.north = 0;
450
    if(dy_s != 0)
451
        value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
452
    else
453
        value.south = 0;
454
    if(dz_f != 0)
455
        value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
456
    else
457
        value.front = 0;
458
    if(dz_b != 0)
459
        value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
460
    else
461
        value.back = 0;
462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485

    //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
486
    value.center = 2 / hr_m[2];
487
    if(dx_w != 0 || dx_e != 0)
488
        value.center *= (dx_w + dx_e);
489
    if(dy_n != 0 || dy_s != 0)
490
        value.center *= (dy_n + dy_s);
491
    if(dx_w != 0 || dx_e != 0)
492
        value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
493
    if(dy_n != 0 || dy_s != 0)
494
        value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
495
}
gsell's avatar
gsell committed
496 497 498 499 500 501 502

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

503
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
504 505 506 507 508
    // 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
509 510 511

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
512
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
513 514 515
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
519 520 521 522
    v(0) = (quaternion(0) * quaternion(0)
            + quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
kraus's avatar
kraus committed
523

524 525 526 527
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
533 534 535 536
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
537 538 539 540

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

541
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
542 543
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
544
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
545

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

}
552 553 554 555 556 557 558 559

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