ArbitraryDomain.cpp 25.4 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 7 8 9 10 11 12
// Copyright (c) 2008-2020
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
// OPAL is licensed under GNU GPL version 3.
//

13
//#define DEBUG_INTERSECT_RAY_BOUNDARY
14

kraus's avatar
kraus committed
15
#include "Solvers/ArbitraryDomain.h"
snuverink_j's avatar
snuverink_j committed
16
#include "Structure/BoundaryGeometry.h"
kraus's avatar
kraus committed
17

18 19
#include <cmath>
#include <iostream>
snuverink_j's avatar
snuverink_j committed
20
#include <tuple>
21
#include <assert.h>
gsell's avatar
gsell committed
22

23 24
#include <math.h>

25
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
kraus's avatar
kraus committed
26
                                  Vector_t nr,
27
                                  Vector_t /*hr*/,
kraus's avatar
kraus committed
28
                                  std::string interpl){
29 30 31
    bgeom_m  = bgeom;
    minCoords_m = bgeom->getmincoords();
    maxCoords_m = bgeom->getmaxcoords();
32
    geomCentroid_m = (minCoords_m + maxCoords_m)/2.0;
33 34

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

    setNr(nr);
39
    for(int i=0; i<3; i++)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
40
        Geo_hr_m[i] = (maxCoords_m[i] - minCoords_m[i])/nr[i];
41
    setHr(Geo_hr_m);
42 43 44 45 46 47 48 49 50

    startId = 0;

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

53 54
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
55 56
}

57 58
void ArbitraryDomain::compute(Vector_t hr){

59 60
    setHr(hr);

61
    globalMeanR_m = getGlobalMeanR();
gsell's avatar
gsell committed
62

63 64 65
    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
66
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
kraus's avatar
kraus committed
67

68 69
    hasGeometryChanged_m = true;

70 71 72 73 74 75 76
    IntersectLoX.clear();
    IntersectHiX.clear();
    IntersectLoY.clear();
    IntersectHiY.clear();
    IntersectLoZ.clear();
    IntersectHiZ.clear();

kraus's avatar
kraus committed
77
    //calculate intersection
78 79 80
    Vector_t P, saveP, dir, I;
    //Reference Point
    Vector_t P0 = geomCentroid_m;
81

82 83
    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
84 85 86 87 88
        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;
89

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
90 91
                rotateWithQuaternion(P, localToGlobalQuaternion_m);
                P += geomCentroid_m;
92

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
93 94
                if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
                    P0 = P;
95

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
98 99
                    rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
100
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
101 102
                        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
103
                    } else {
104
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
105
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
106
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
107
                    }
108

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

                    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 {
125
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
126
                        *gmsg << "ydir=+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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
140 141
                    rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
142
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
143 144
                        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
145
                    } else {
146
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
147
                        *gmsg << "xdir=+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
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
152
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
153 154
                        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
155
                    } else {
156
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
157
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
158
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
159 160
                    }
                } else {
161
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
162
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
163
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
164 165 166 167
                }
            }
        }
    }
gsell's avatar
gsell committed
168 169
    IdxMap.clear();
    CoordMap.clear();
kraus's avatar
kraus committed
170

171 172
    int id=0;
    int idx, idy, idz;
173 174 175
    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
176
                if (isInside(idx, idy, idz)) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
177 178 179
                    IdxMap[toCoordIdx(idx, idy, idz)] = id;
                    CoordMap[id] = toCoordIdx(idx, idy, idz);
                    id++;
kraus's avatar
kraus committed
180
                }
181
            }
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
182 183
        }
    }
gsell's avatar
gsell committed
184 185
}

186
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
187

188 189
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

190
    setHr(hr);
191

192 193 194 195
    globalMeanR_m = getGlobalMeanR();

    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
196
    for (int i=1; i<4; i++)
197
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
198 199 200 201 202 203 204

    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
205

206 207 208 209 210 211 212 213 214
    hasGeometryChanged_m = true;

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

215
    // Calculate intersection
216 217 218
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
219

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
240
                //P = saveP;
241

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
242 243 244
                //rotateWithQuaternion(P, localToGlobalQuaternion_m);
                //P += geomCentroid_m; //sorry, this doesn't make sense. -DW
                //P += globalMeanR_m;
245

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
248 249 250
                    // 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;
251

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
252 253 254 255 256 257
                    // 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;
258

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
261 262
                    //rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 0, 1);
263

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
264
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
265 266
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
267
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
kraus's avatar
kraus committed
268
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
269
                    } else {
270
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
271
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
272
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
273
                    }
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 278 279
                        //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
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
                    }
kraus's avatar
kraus committed
285

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
286 287
                    //rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 1, 0);
288

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
289 290 291 292 293 294
                    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 {
295
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
296
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
297
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
298
                    }
299

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
300
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
301 302 303 304
                        //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
305
                    } 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 312
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
313

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
314
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
315 316 317 318
                        //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
319
                    } else {
320
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
321
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
322
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
323
                    }
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);
                        IntersectLoX.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 336
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
337
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
338
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
339
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
340 341 342 343
                }
            }
        }
    }
344

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

347 348 349 350 351
    //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
352
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
353 354 355 356 357
                    numGhostNodesLeft++;
            }
        }
    }

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

360 361 362 363
    //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
364
        int numxy = 0;
365 366 367 368 369 370 371 372 373 374 375 376 377
        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;
378

379
    // Build up index and coord map
380 381
    IdxMap.clear();
    CoordMap.clear();
382
    int index = startIdx - numGhostNodesLeft;
383

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

386 387 388
    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
389
                if(isInside(idx, idy, idz)) {
390 391 392 393
                    IdxMap[toCoordIdx(idx, idy, idz)] = index;
                    CoordMap[index] = toCoordIdx(idx, idy, idz);
                    index++;
                }
394
            }
395 396
        }
    }
397 398

    INFOMSG(level2 << "* Done." << endl);
399
}
400

401 402
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
403
    return (idz * nr[1] + idy) * nr[0]  + idx;
gsell's avatar
gsell committed
404 405
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
409 410 411 412
    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
413
}
gsell's avatar
gsell committed
414

415 416 417 418 419 420 421 422
// 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
423 424
}

425
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
426

427 428 429 430
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
434 435 436
  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];
437

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
438 439
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
440 441 442
*/

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
443 444 445 446 447 448 449 450 451 452 453 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
  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;
  }
484
*/
gsell's avatar
gsell committed
485 486

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
488
    return numXY[z];
489
}
gsell's avatar
gsell committed
490

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

494 495
    getCoord(idxyz, idx, idy, idz);
    getBoundaryStencil(idx, idy, idz, W, E, S, N, F, B, C, scaleFactor);
gsell's avatar
gsell committed
496 497
}

498
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
499

500
    scaleFactor = 1.0;
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
501
    // determine which interpolation method we use for points near the boundary
502
    switch(interpolationMethod){
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
503 504 505 506 507 508 509 510 511
    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;
512 513 514 515 516 517
    }

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

518
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*/) {
519 520 521 522 523 524 525 526 527

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

528
    if(!isInside(idx-1,idy,idz))
529
        W = 0.0;
kraus's avatar
kraus committed
530
    if(!isInside(idx+1,idy,idz))
531
        E = 0.0;
532

533
    if(!isInside(idx,idy+1,idz))
534
        N = 0.0;
kraus's avatar
kraus committed
535
    if(!isInside(idx,idy-1,idz))
536
        S = 0.0;
kraus's avatar
kraus committed
537 538

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

544
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)
545
{
546
    scaleFactor = 1;
547

548 549 550
    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];
551

552
    double dx_w=hr[0];
553 554
    double dx_e=hr[0];
    double dy_n=hr[1];
555 556 557
    double dy_s=hr[1];
    double dz_f=hr[2];
    double dz_b=hr[2];
558 559
    C = 0.0;

560
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
561

562
    if (idx == nr[0]-1)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
563
        dx_e = fabs(IntersectHiX.find(coordxyz)->second - cx);
564
    if (idx == 0)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
565
        dx_w = fabs(IntersectLoX.find(coordxyz)->second - cx);
566
    if (idy == nr[1]-1)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
567
        dy_n = fabs(IntersectHiY.find(coordxyz)->second - cy);
568
    if (idy == 0)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
569
        dy_s = fabs(IntersectLoY.find(coordxyz)->second - cy);
570
    if (idz == nr[2]-1)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
571
        dz_b = fabs(IntersectHiZ.find(coordxyz)->second - cz);
572
    if (idz == 0)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
573
        dz_f = fabs(IntersectLoZ.find(coordxyz)->second - cz);
574 575 576 577 578 579 580 581 582 583 584 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

    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;
632
}
gsell's avatar
gsell committed
633

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

636
    int idx = 0, idy = 0, idz = 0;
gsell's avatar
gsell committed
637

638 639
    getCoord(id, idx, idy, idz);
    getNeighbours(idx, idy, idz, W, E, S, N, F, B);
gsell's avatar
gsell committed
640 641
}

642
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
643

644 645 646 647 648 649
    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);
650

kraus's avatar
kraus committed
651
    if(!isInside(idx+1,idy,idz))
gsell's avatar
gsell committed
652 653
        E = -1;

654 655 656 657
    if(!isInside(idx-1,idy,idz))
        W = -1;

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

kraus's avatar
kraus committed
660
    if(!isInside(idx,idy-1,idz))
gsell's avatar
gsell committed
661 662
        S = -1;

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

    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
667
        B = -1;
gsell's avatar
gsell committed
668 669 670 671 672 673 674 675 676 677

}


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

678
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
679 680 681 682 683
    // 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
684 685 686

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
687
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
688 689 690
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
694 695 696 697
    v(0) = (quaternion(0) * quaternion(0)
            + quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
kraus's avatar
kraus committed
698

699 700 701 702
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

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

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

716
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
717 718
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
719
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
720

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

}
727 728 729 730 731 732 733 734

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