ArbitraryDomain.cpp 25.7 KB
Newer Older
1 2 3 4 5 6
// ------------------------------------------------------------------------
// $Version: 1.2.1 $
// ------------------------------------------------------------------------
// Copyright & License: See Copyright.readme in src directory
// ------------------------------------------------------------------------
// Class ArbitraryDomain
kraus's avatar
kraus committed
7
//   Interface to iterative solver and boundary geometry
8 9 10 11 12 13
//   for space charge calculation
//
// ------------------------------------------------------------------------
// $Author: kaman $
// $Date: 2014 $
// ------------------------------------------------------------------------
14
//#define DEBUG_INTERSECT_RAY_BOUNDARY
15

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

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

24 25
#include <math.h>

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

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

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

    startId = 0;

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

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

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

60 61
    setHr(hr);

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

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

69 70
    hasGeometryChanged_m = true;

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
110
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
111
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
112 113
                        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
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 120 121 122 123 124 125
                    }

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
152
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
153
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
154 155
                        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
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
                    }
                } else {
162
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
163
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
164
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
165 166 167 168
                }
            }
        }
    }
gsell's avatar
gsell committed
169 170
    IdxMap.clear();
    CoordMap.clear();
kraus's avatar
kraus committed
171

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

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

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

191
    setHr(hr);
192

193 194 195 196
    globalMeanR_m = getGlobalMeanR();

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

    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
206

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
276
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
277 278 279 280
                        //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
281
                    } else {
282
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
283
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
284
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
285
                    }
kraus's avatar
kraus committed
286

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
312 313
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
314

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
326
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
327 328 329 330
                        //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
331
                    } else {
332
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
333
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
334
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
335 336 337
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
338
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
339
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
340
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
341 342 343 344
                }
            }
        }
    }
345

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
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 484
  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;
  }
485
*/
gsell's avatar
gsell committed
486 487

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

}


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

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

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

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

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

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

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

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

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

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

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

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

}
728 729 730 731 732 733 734 735

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