ArbitraryDomain.cpp 27.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
//
6
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
7 8
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//                      2016, Daniel Winklehner, Massachusetts Institute of Technology
9
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
10
// All rights reserved
11
//
12 13 14 15 16
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
frey_m's avatar
frey_m committed
17 18 19 20 21 22 23 24 25 26
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
27 28
//

29
//#define DEBUG_INTERSECT_RAY_BOUNDARY
30

kraus's avatar
kraus committed
31
#include "Solvers/ArbitraryDomain.h"
snuverink_j's avatar
snuverink_j committed
32
#include "Structure/BoundaryGeometry.h"
kraus's avatar
kraus committed
33

34 35
#include <cmath>
#include <iostream>
snuverink_j's avatar
snuverink_j committed
36
#include <tuple>
37
#include <cassert>
38
#include "Utilities/OpalException.h"
gsell's avatar
gsell committed
39

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

49 50
    bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m);
    if (have_inside_pt == false) {
51 52 53
        throw OpalException(
            "ArbitraryDomain::ArbitraryDomain()",
            "No point inside geometry found/set!");
54
    }
55
    setNr(nr);
56
    for(int i=0; i<3; i++)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
57
        Geo_hr_m[i] = (maxCoords_m[i] - minCoords_m[i])/nr[i];
58
    setHr(Geo_hr_m);
59 60 61 62 63 64 65 66 67

    startId = 0;

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

70 71
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
72 73
}

74 75
void ArbitraryDomain::compute(Vector_t hr){

76 77
    setHr(hr);

78
    globalMeanR_m = getGlobalMeanR();
gsell's avatar
gsell committed
79

80 81 82
    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
83
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
kraus's avatar
kraus committed
84

85 86
    hasGeometryChanged_m = true;

87 88 89 90 91 92 93
    IntersectLoX.clear();
    IntersectHiX.clear();
    IntersectLoY.clear();
    IntersectHiY.clear();
    IntersectLoZ.clear();
    IntersectHiZ.clear();

kraus's avatar
kraus committed
94
    //calculate intersection
95 96 97
    Vector_t P, saveP, dir, I;
    //Reference Point
    Vector_t P0 = geomCentroid_m;
98

99 100
    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
101 102 103 104 105
        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;
106

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
107 108
                rotateWithQuaternion(P, localToGlobalQuaternion_m);
                P += geomCentroid_m;
109

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
110 111
                if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
                    P0 = P;
112

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

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

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

                    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 {
144
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
145 146
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
147
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
148
                    }
149

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
150
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
151
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
152 153
                        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
154
                    } else {
155
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
156 157
                        *gmsg << "ydir=-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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
161 162
                    rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
163
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
164 165
                        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
166
                    } else {
167
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
168 169
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
170
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
171
                    }
172

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
173
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
174
                        I -= geomCentroid_m;
kraus's avatar
kraus committed
175 176
                        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
177
                    } else {
178
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
179 180
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
181
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
182 183
                    }
                } else {
184
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
185 186
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
187
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
188 189 190 191
                }
            }
        }
    }
gsell's avatar
gsell committed
192 193
    IdxMap.clear();
    CoordMap.clear();
kraus's avatar
kraus committed
194

195 196
    int id=0;
    int idx, idy, idz;
197 198 199
    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
200
                if (isInside(idx, idy, idz)) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
201 202 203
                    IdxMap[toCoordIdx(idx, idy, idz)] = id;
                    CoordMap[id] = toCoordIdx(idx, idy, idz);
                    id++;
kraus's avatar
kraus committed
204
                }
205
            }
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
206 207
        }
    }
gsell's avatar
gsell committed
208 209
}

210
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
211

212 213
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

214
    setHr(hr);
215

216 217 218 219
    globalMeanR_m = getGlobalMeanR();

    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
220
    for (int i=1; i<4; i++)
221
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
222 223 224 225 226 227 228

    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
229

230 231 232 233 234 235 236 237 238
    hasGeometryChanged_m = true;

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

239
    // Calculate intersection
240 241 242
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
243

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
264
                //P = saveP;
265

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
266 267 268
                //rotateWithQuaternion(P, localToGlobalQuaternion_m);
                //P += geomCentroid_m; //sorry, this doesn't make sense. -DW
                //P += globalMeanR_m;
269

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
272 273 274
                    // 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;
275

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
276 277 278 279 280 281
                    // 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;
282

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
288
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
289 290
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
291
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
kraus's avatar
kraus committed
292
                        IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
293
                    } else {
294
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
295 296
                        *gmsg << "zdir=+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);
                        IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
305
                    } else {
306
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
307 308
                        *gmsg << "zdir=-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
                    }
kraus's avatar
kraus committed
311

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
327
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
328 329 330 331
                        //I -= geomCentroid_m;
                        //I -= globalMeanR_m;
                        //rotateWithQuaternion(I, globalToLocalQuaternion_m);
                        IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
332
                    } else {
333
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
334 335
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
336
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
337
                    }
338

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
339 340
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
341

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
342
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
343 344 345 346
                        //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
347
                    } else {
348
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
349 350
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
351
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
352
                    }
353

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
354
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
kraus's avatar
kraus committed
355 356 357 358
                        //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
359
                    } else {
360
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
361 362
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
363
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
364 365 366
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
367
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
368 369
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
370
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
371 372 373 374
                }
            }
        }
    }
375

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

378 379 380 381 382
    //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
383
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
384 385 386 387 388
                    numGhostNodesLeft++;
            }
        }
    }

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

391 392 393 394
    //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
395
        int numxy = 0;
396 397 398 399 400 401 402 403 404 405 406
        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;
407
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
408
    startIdx -= numtotal;
409

410
    // Build up index and coord map
411 412
    IdxMap.clear();
    CoordMap.clear();
413
    int index = startIdx - numGhostNodesLeft;
414

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

417 418 419
    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
420
                if(isInside(idx, idy, idz)) {
421 422 423 424
                    IdxMap[toCoordIdx(idx, idy, idz)] = index;
                    CoordMap[index] = toCoordIdx(idx, idy, idz);
                    index++;
                }
425
            }
426 427
        }
    }
428 429

    INFOMSG(level2 << "* Done." << endl);
430
}
431

432 433
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
434
    return (idz * nr[1] + idy) * nr[0]  + idx;
gsell's avatar
gsell committed
435 436
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
440 441 442 443
    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
444
}
gsell's avatar
gsell committed
445

446 447 448 449 450 451 452 453
// 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
454 455
}

456
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
457

458 459 460 461
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
465 466 467
  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];
468

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
469 470
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
471 472 473
*/

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514
  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;
  }
515
*/
gsell's avatar
gsell committed
516 517

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
519
    return numXY[z];
520
}
gsell's avatar
gsell committed
521

frey_m's avatar
frey_m committed
522 523 524 525
void ArbitraryDomain::getBoundaryStencil(int idxyz, double &W, double &E, double &S,
                                         double &N, double &F, double &B, double &C,
                                         double &scaleFactor)
{
526
    int idx = 0, idy = 0, idz = 0;
527

528 529
    getCoord(idxyz, idx, idy, idz);
    getBoundaryStencil(idx, idy, idz, W, E, S, N, F, B, C, scaleFactor);
gsell's avatar
gsell committed
530 531
}

frey_m's avatar
frey_m committed
532 533 534 535
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)
{
536
    scaleFactor = 1.0;
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
537
    // determine which interpolation method we use for points near the boundary
538
    switch(interpolationMethod){
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
539 540 541 542 543 544 545 546 547
    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;
548 549 550 551 552 553
    }

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

frey_m's avatar
frey_m committed
554 555 556 557
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*/)
{
558 559 560 561 562 563 564 565
    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]);

566
    if(!isInside(idx-1,idy,idz))
567
        W = 0.0;
kraus's avatar
kraus committed
568
    if(!isInside(idx+1,idy,idz))
569
        E = 0.0;
570

571
    if(!isInside(idx,idy+1,idz))
572
        N = 0.0;
kraus's avatar
kraus committed
573
    if(!isInside(idx,idy-1,idz))
574
        S = 0.0;
kraus's avatar
kraus committed
575 576

    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
577
        F = 0.0;
kraus's avatar
kraus committed
578
    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
579
        B = 0.0;
580
}
581

frey_m's avatar
frey_m committed
582 583 584
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)
585
{
586
    scaleFactor = 1;
587

588 589 590
    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];
591

592
    double dx_w=hr[0];
593 594
    double dx_e=hr[0];
    double dy_n=hr[1];
595 596 597
    double dy_s=hr[1];
    double dz_f=hr[2];
    double dz_b=hr[2];
598 599
    C = 0.0;

600
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
601

602
    if (idx == nr[0]-1)
603
        dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx);
604
    if (idx == 0)
605
        dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx);
606
    if (idy == nr[1]-1)
607
        dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy);
608
    if (idy == 0)
609
        dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy);
610
    if (idz == nr[2]-1)
611
        dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz);
612
    if (idz == 0)
613
        dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz);
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671

    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;
672
}
gsell's avatar
gsell committed
673

frey_m's avatar
frey_m committed
674 675 676
void ArbitraryDomain::getNeighbours(int id, int &W, int &E, int &S,
                                    int &N, int &F, int &B)
{
677
    int idx = 0, idy = 0, idz = 0;
gsell's avatar
gsell committed
678

679 680
    getCoord(id, idx, idy, idz);
    getNeighbours(idx, idy, idz, W, E, S, N, F, B);
gsell's avatar
gsell committed
681 682
}

frey_m's avatar
frey_m committed
683 684 685
void ArbitraryDomain::getNeighbours(int idx, int idy, int idz, int &W,
                                    int &E, int &S, int &N, int &F, int &B)
{
686 687 688 689 690 691
    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);
692

kraus's avatar
kraus committed
693
    if(!isInside(idx+1,idy,idz))
gsell's avatar
gsell committed
694 695
        E = -1;

696 697 698 699
    if(!isInside(idx-1,idy,idz))
        W = -1;

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

kraus's avatar
kraus committed
702
    if(!isInside(idx,idy-1,idz))
gsell's avatar
gsell committed
703 704
        S = -1;

kraus's avatar
kraus committed
705
    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
706
        F = -1;
kraus's avatar
kraus committed
707 708

    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
709
        B = -1;
gsell's avatar
gsell committed
710 711 712 713 714 715 716 717 718 719

}


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

720
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
721 722 723 724 725
    // 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
726 727 728

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
729
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
730 731 732
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

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

741 742 743 744
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
750 751 752 753
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
754 755 756 757

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

758
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
759 760
    // rotates the positive zaxis using a quaternion.
    v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
kraus's avatar
kraus committed
761
    v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
762

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
763 764 765 766
    v(2) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            + quaternion(3) * quaternion(3));
767 768

}
769 770 771 772 773 774 775 776

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