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
// Implemented as part of the PhD thesis
frey_m's avatar
frey_m committed
12
// "Toward massively parallel multi-objective optimization with application to
frey_m's avatar
frey_m committed
13 14 15 16 17 18 19 20 21 22 23
// 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
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
kraus's avatar
kraus committed
37
                                  Vector_t nr,
38
                                  Vector_t /*hr*/,
kraus's avatar
kraus committed
39
                                  std::string interpl){
40 41 42
    bgeom_m  = bgeom;
    minCoords_m = bgeom->getmincoords();
    maxCoords_m = bgeom->getmaxcoords();
43
    geomCentroid_m = (minCoords_m + maxCoords_m)/2.0;
44 45

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

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

    startId = 0;

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

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

68 69
void ArbitraryDomain::compute(Vector_t hr){

70 71
    setHr(hr);

72
    globalMeanR_m = getGlobalMeanR();
gsell's avatar
gsell committed
73

74 75 76
    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
77
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
kraus's avatar
kraus committed
78

79 80
    hasGeometryChanged_m = true;

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

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

93 94
    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
95 96 97 98 99
        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;
100

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
109 110
                    rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
111
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
112 113
                        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
114
                    } else {
115
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
116
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
117
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
118
                    }
119

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

                    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 {
136
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
137
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
138
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
139
                    }
140

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
151 152
                    rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
153
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
154 155
                        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
156
                    } else {
157
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
158
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
159
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
160
                    }
161

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

182 183
    int id=0;
    int idx, idy, idz;
184 185 186
    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
187
                if (isInside(idx, idy, idz)) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
188 189 190
                    IdxMap[toCoordIdx(idx, idy, idz)] = id;
                    CoordMap[id] = toCoordIdx(idx, idy, idz);
                    id++;
kraus's avatar
kraus committed
191
                }
192
            }
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
193 194
        }
    }
gsell's avatar
gsell committed
195 196
}

197
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
198

199 200
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

201
    setHr(hr);
202

203 204 205 206
    globalMeanR_m = getGlobalMeanR();

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

    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
216

217 218 219 220 221 222 223 224 225
    hasGeometryChanged_m = true;

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

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

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
251
                //P = saveP;
252

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    INFOMSG(level2 << "* Done." << endl);
410
}
411

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

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

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

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

436
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
437

438 439 440 441
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

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

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

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
454 455 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
  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;
  }
495
*/
gsell's avatar
gsell committed
496 497

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

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

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

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

509
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
510

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

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

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

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

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

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

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

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

559 560 561
    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];
562

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

571
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
572

573
    if (idx == nr[0]-1)
574
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
575
    if (idx == 0)
576
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
577
    if (idy == nr[1]-1)
578
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
579
    if (idy == 0)
580
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
581
    if (idz == nr[2]-1)
582
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
583
    if (idz == 0)
584
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
585 586 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

    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;
643
}
gsell's avatar
gsell committed
644

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

647
    int idx = 0, idy = 0, idz = 0;
gsell's avatar
gsell committed
648

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

653
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
654

655 656 657 658 659 660
    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);
661

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

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

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

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

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

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

}


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

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

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

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

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

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

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

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

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

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

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

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

}
738 739 740 741 742 743 744

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