ArbitraryDomain.cpp 26.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
//
frey_m's avatar
frey_m committed
6 7 8 9
// Copyright (c) 2010 - 2013, Yves Ineichen, ETH Zürich,
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//                      2016, Daniel Winklehner, Massachusetts Institute of Technology
// All rights reserved
10
//
frey_m's avatar
frey_m committed
11 12 13 14 15 16 17 18 19 20 21 22 23
// Implemented as part of the PhD thesis
// "Toward massively parallel multi-objective optimization withapplication to
// particle accelerators" (https://doi.org/10.3929/ethz-a-009792359)
//
// 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/>.
24 25
//

26
//#define DEBUG_INTERSECT_RAY_BOUNDARY
27

kraus's avatar
kraus committed
28
#include "Solvers/ArbitraryDomain.h"
snuverink_j's avatar
snuverink_j committed
29
#include "Structure/BoundaryGeometry.h"
kraus's avatar
kraus committed
30

31 32
#include <cmath>
#include <iostream>
snuverink_j's avatar
snuverink_j committed
33
#include <tuple>
34
#include <assert.h>
gsell's avatar
gsell committed
35

36 37
#include <math.h>

38
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
kraus's avatar
kraus committed
39
                                  Vector_t nr,
40
                                  Vector_t /*hr*/,
kraus's avatar
kraus committed
41
                                  std::string interpl){
42 43 44
    bgeom_m  = bgeom;
    minCoords_m = bgeom->getmincoords();
    maxCoords_m = bgeom->getmaxcoords();
45
    geomCentroid_m = (minCoords_m + maxCoords_m)/2.0;
46 47

    // TODO: THis needs to be made into OPTION of the geometry.
48 49
    // A user defined point that is INSIDE with 100% certainty. -DW
    globalInsideP0_m = Vector_t(0.0, 0.0, -0.13);
50 51

    setNr(nr);
52
    for(int i=0; i<3; i++)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
53
        Geo_hr_m[i] = (maxCoords_m[i] - minCoords_m[i])/nr[i];
54
    setHr(Geo_hr_m);
55 56 57 58 59 60 61 62 63

    startId = 0;

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

66 67
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
68 69
}

70 71
void ArbitraryDomain::compute(Vector_t hr){

72 73
    setHr(hr);

74
    globalMeanR_m = getGlobalMeanR();
gsell's avatar
gsell committed
75

76 77 78
    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
    for (int i=1; i<4; i++)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
79
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
kraus's avatar
kraus committed
80

81 82
    hasGeometryChanged_m = true;

83 84 85 86 87 88 89
    IntersectLoX.clear();
    IntersectHiX.clear();
    IntersectLoY.clear();
    IntersectHiY.clear();
    IntersectLoZ.clear();
    IntersectHiZ.clear();

kraus's avatar
kraus committed
90
    //calculate intersection
91 92 93
    Vector_t P, saveP, dir, I;
    //Reference Point
    Vector_t P0 = geomCentroid_m;
94

95 96
    for (int idz = 0; idz < nr[2] ; idz++) {
        saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
97 98 99 100 101
        for (int idy = 0; idy < nr[1] ; idy++) {
            saveP[1] = (idy - (nr[1]-1)/2.0)*hr[1];
            for (int idx = 0; idx <nr[0]; idx++) {
                saveP[0] = (idx - (nr[0]-1)/2.0)*hr[0];
                P = saveP;
102

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
103 104
                rotateWithQuaternion(P, localToGlobalQuaternion_m);
                P += geomCentroid_m;
105

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
106 107
                if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
                    P0 = P;
108

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
111 112
                    rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
113
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
114 115
                        rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
116
                    } else {
117
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
118
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
119
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
120
                    }
121

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
122
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
123
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
124 125
                        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
126
                    } else {
127
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
128
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
129
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
130 131 132 133 134 135 136 137
                    }

                    rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
                        I -= geomCentroid_m;
                        rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
                    } else {
138
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
139
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
140
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
141
                    }
142

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
143
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
144
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
145 146
                        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
147
                    } else {
148
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
149
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
150
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
151
                    }
152

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
153 154
                    rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
155
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
156 157
                        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
158
                    } else {
159
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
160
                        *gmsg << "xdir=+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)){
165
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
166 167
                        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
168
                    } else {
169
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
170
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
171
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
172 173
                    }
                } else {
174
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
175
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
176
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
177 178 179 180
                }
            }
        }
    }
gsell's avatar
gsell committed
181 182
    IdxMap.clear();
    CoordMap.clear();
kraus's avatar
kraus committed
183

184 185
    int id=0;
    int idx, idy, idz;
186 187 188
    for (idz = 0; idz < nr[2]; idz++) {
        for (idy = 0; idy < nr[1]; idy++) {
            for (idx = 0; idx < nr[0]; idx++) {
kraus's avatar
kraus committed
189
                if (isInside(idx, idy, idz)) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
190 191 192
                    IdxMap[toCoordIdx(idx, idy, idz)] = id;
                    CoordMap[id] = toCoordIdx(idx, idy, idz);
                    id++;
kraus's avatar
kraus committed
193
                }
194
            }
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
195 196
        }
    }
gsell's avatar
gsell committed
197 198
}

199
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
200

201 202
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

203
    setHr(hr);
204

205 206 207 208
    globalMeanR_m = getGlobalMeanR();

    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
209
    for (int i=1; i<4; i++)
210
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
211 212 213 214 215 216 217

    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
218

219 220 221 222 223 224 225 226 227
    hasGeometryChanged_m = true;

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

228
    // Calculate intersection
229 230 231
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
232

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
253
                //P = saveP;
254

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
255 256 257
                //rotateWithQuaternion(P, localToGlobalQuaternion_m);
                //P += geomCentroid_m; //sorry, this doesn't make sense. -DW
                //P += globalMeanR_m;
258

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
261 262 263
                    // 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;
264

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
265 266 267 268 269 270
                    // 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;
271

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
274 275
                    //rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 0, 1);
276

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
277
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
278 279
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
280
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
kraus's avatar
kraus committed
281
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
282
                    } else {
283
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
284
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
285
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
286
                    }
287

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
288
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
289 290 291 292
                        //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
293
                    } else {
294
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
295
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
296
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
297
                    }
kraus's avatar
kraus committed
298

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
299 300
                    //rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 1, 0);
301

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
302 303 304 305 306 307
                    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 {
308
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
309
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
310
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
311
                    }
312

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
313
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
314 315 316 317
                        //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
318
                    } else {
319
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
320
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
321
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
322
                    }
323

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
324 325
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
326

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
327
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
328 329 330 331
                        //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
332
                    } else {
333
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
334
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
335
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
336
                    }
337

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
338
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
339 340 341 342
                        //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
343
                    } else {
344
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
345
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
346
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
347 348 349
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
350
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
351
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
352
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
353 354 355 356
                }
            }
        }
    }
357

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

360 361 362 363 364
    //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
365
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
366 367 368 369 370
                    numGhostNodesLeft++;
            }
        }
    }

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

373 374 375 376
    //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
377
        int numxy = 0;
378 379 380 381 382 383 384 385 386 387 388 389 390
        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;
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD);
    startIdx -= numtotal;
391

392
    // Build up index and coord map
393 394
    IdxMap.clear();
    CoordMap.clear();
395
    int index = startIdx - numGhostNodesLeft;
396

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

399 400 401
    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
402
                if(isInside(idx, idy, idz)) {
403 404 405 406
                    IdxMap[toCoordIdx(idx, idy, idz)] = index;
                    CoordMap[index] = toCoordIdx(idx, idy, idz);
                    index++;
                }
407
            }
408 409
        }
    }
410 411

    INFOMSG(level2 << "* Done." << endl);
412
}
413

414 415
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
416
    return (idz * nr[1] + idy) * nr[0]  + idx;
gsell's avatar
gsell committed
417 418
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
422 423 424 425
    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
426
}
gsell's avatar
gsell committed
427

428 429 430 431 432 433 434 435
// 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
436 437
}

438
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
439

440 441 442 443
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
447 448 449
  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];
450

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
451 452
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
453 454 455
*/

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
456 457 458 459 460 461 462 463 464 465 466 467 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
  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;
  }
497
*/
gsell's avatar
gsell committed
498 499

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
501
    return numXY[z];
502
}
gsell's avatar
gsell committed
503

504
void ArbitraryDomain::getBoundaryStencil(int idxyz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
505
    int idx = 0, idy = 0, idz = 0;
506

507 508
    getCoord(idxyz, idx, idy, idz);
    getBoundaryStencil(idx, idy, idz, W, E, S, N, F, B, C, scaleFactor);
gsell's avatar
gsell committed
509 510
}

511
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) {
gsell's avatar
gsell committed
512

513
    scaleFactor = 1.0;
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
514
    // determine which interpolation method we use for points near the boundary
515
    switch(interpolationMethod){
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
516 517 518 519 520 521 522 523 524
    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;
525 526 527 528 529 530
    }

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

531
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*/) {
532 533 534 535 536 537 538 539 540

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

541
    if(!isInside(idx-1,idy,idz))
542
        W = 0.0;
kraus's avatar
kraus committed
543
    if(!isInside(idx+1,idy,idz))
544
        E = 0.0;
545

546
    if(!isInside(idx,idy+1,idz))
547
        N = 0.0;
kraus's avatar
kraus committed
548
    if(!isInside(idx,idy-1,idz))
549
        S = 0.0;
kraus's avatar
kraus committed
550 551

    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
552
        F = 0.0;
kraus's avatar
kraus committed
553
    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
554
        B = 0.0;
555
}
556

557
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)
558
{
559
    scaleFactor = 1;
560

561 562 563
    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];
564

565
    double dx_w=hr[0];
566 567
    double dx_e=hr[0];
    double dy_n=hr[1];
568 569 570
    double dy_s=hr[1];
    double dz_f=hr[2];
    double dz_b=hr[2];
571 572
    C = 0.0;

573
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
574

575
    if (idx == nr[0]-1)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
576
        dx_e = fabs(IntersectHiX.find(coordxyz)->second - cx);
577
    if (idx == 0)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
578
        dx_w = fabs(IntersectLoX.find(coordxyz)->second - cx);
579
    if (idy == nr[1]-1)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
580
        dy_n = fabs(IntersectHiY.find(coordxyz)->second - cy);
581
    if (idy == 0)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
582
        dy_s = fabs(IntersectLoY.find(coordxyz)->second - cy);
583
    if (idz == nr[2]-1)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
584
        dz_b = fabs(IntersectHiZ.find(coordxyz)->second - cz);
585
    if (idz == 0)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
586
        dz_f = fabs(IntersectLoZ.find(coordxyz)->second - cz);
587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644

    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;
645
}
gsell's avatar
gsell committed
646

647
void ArbitraryDomain::getNeighbours(int id, int &W, int &E, int &S, int &N, int &F, int &B) {
gsell's avatar
gsell committed
648

649
    int idx = 0, idy = 0, idz = 0;
gsell's avatar
gsell committed
650

651 652
    getCoord(id, idx, idy, idz);
    getNeighbours(idx, idy, idz, W, E, S, N, F, B);
gsell's avatar
gsell committed
653 654
}

655
void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W, int &E, int &S, int &N, int &F, int &B) {
gsell's avatar
gsell committed
656

657 658 659 660 661 662
    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);
663

kraus's avatar
kraus committed
664
    if(!isInside(idx+1,idy,idz))
gsell's avatar
gsell committed
665 666
        E = -1;

667 668 669 670
    if(!isInside(idx-1,idy,idz))
        W = -1;

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

kraus's avatar
kraus committed
673
    if(!isInside(idx,idy-1,idz))
gsell's avatar
gsell committed
674 675
        S = -1;

kraus's avatar
kraus committed
676
    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
677
        F = -1;
kraus's avatar
kraus committed
678 679

    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
680
        B = -1;
gsell's avatar
gsell committed
681 682 683 684 685 686 687 688 689 690

}


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

691
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
692 693 694 695 696
    // 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
697 698 699

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
700
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
701 702 703
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
707 708 709 710
    v(0) = (quaternion(0) * quaternion(0)
            + quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
kraus's avatar
kraus committed
711

712 713 714 715
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
721 722 723 724
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
725 726 727 728

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

729
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
730 731
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
732
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
733

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
734 735 736 737
    v(2) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            + quaternion(3) * quaternion(3));
738 739

}
740 741 742 743 744 745 746 747

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