ArbitraryDomain.cpp 26.3 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
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
122
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
123
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
124
                    }
125

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

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

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

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

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

188 189
    int id=0;
    int idx, idy, idz;
190 191 192
    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
193
                if (isInside(idx, idy, idz)) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
194 195 196
                    IdxMap[toCoordIdx(idx, idy, idz)] = id;
                    CoordMap[id] = toCoordIdx(idx, idy, idz);
                    id++;
kraus's avatar
kraus committed
197
                }
198
            }
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
199 200
        }
    }
gsell's avatar
gsell committed
201 202
}

203
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
204

205 206
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

207
    setHr(hr);
208

209 210 211 212
    globalMeanR_m = getGlobalMeanR();

    globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
    localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
213
    for (int i=1; i<4; i++)
214
        localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
215 216 217 218 219 220 221

    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
222

223 224 225 226 227 228 229 230 231
    hasGeometryChanged_m = true;

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

232
    // Calculate intersection
233 234 235
    Vector_t P, dir, I;
    // Vector_t saveP, saveP_old;
    Vector_t P0 = globalInsideP0_m;
236

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
257
                //P = saveP;
258

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
259 260 261
                //rotateWithQuaternion(P, localToGlobalQuaternion_m);
                //P += geomCentroid_m; //sorry, this doesn't make sense. -DW
                //P += globalMeanR_m;
262

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
278 279
                    //rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 0, 1);
280

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
303 304
                    //rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(0, 1, 0);
305

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
328 329
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
330

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
331
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
kraus's avatar
kraus committed
332 333 334 335
                        //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
336
                    } else {
337
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
338
                        *gmsg << "xdir=+1 " << dir << " 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

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);
                        IntersectLoX.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
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
349
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
350
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
351 352 353
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
354
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
355
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
356
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
357 358 359 360
                }
            }
        }
    }
361

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

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

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

377 378 379 380
    //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
381
        int numxy = 0;
382 383 384 385 386 387 388 389 390 391 392
        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;
393
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
394
    startIdx -= numtotal;
395

396
    // Build up index and coord map
397 398
    IdxMap.clear();
    CoordMap.clear();
399
    int index = startIdx - numGhostNodesLeft;
400

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

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

    INFOMSG(level2 << "* Done." << endl);
416
}
417

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
426 427 428 429
    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
430
}
gsell's avatar
gsell committed
431

432 433 434 435 436 437 438 439
// 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
440 441
}

442
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
443

444 445 446 447
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
451 452 453
  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];
454

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
455 456
  return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0);
  }
457 458 459
*/

/*
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500
  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;
  }
501
*/
gsell's avatar
gsell committed
502 503

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
505
    return numXY[z];
506
}
gsell's avatar
gsell committed
507

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

511 512
    getCoord(idxyz, idx, idy, idz);
    getBoundaryStencil(idx, idy, idz, W, E, S, N, F, B, C, scaleFactor);
gsell's avatar
gsell committed
513 514
}

515
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
516

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

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

535
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*/) {
536 537 538 539 540 541 542 543 544

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

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

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

    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
556
        F = 0.0;
kraus's avatar
kraus committed
557
    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
558
        B = 0.0;
559
}
560

561
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)
562
{
563
    scaleFactor = 1;
564

565 566 567
    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];
568

569
    double dx_w=hr[0];
570 571
    double dx_e=hr[0];
    double dy_n=hr[1];
572 573 574
    double dy_s=hr[1];
    double dz_f=hr[2];
    double dz_b=hr[2];
575 576
    C = 0.0;

577
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
578

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

    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;
649
}
gsell's avatar
gsell committed
650

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

653
    int idx = 0, idy = 0, idz = 0;
gsell's avatar
gsell committed
654

655 656
    getCoord(id, idx, idy, idz);
    getNeighbours(idx, idy, idz, W, E, S, N, F, B);
gsell's avatar
gsell committed
657 658
}

659
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
660

661 662 663 664 665 666
    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);
667

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

671 672 673 674
    if(!isInside(idx-1,idy,idz))
        W = -1;

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

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

kraus's avatar
kraus committed
680
    if(!isInside(idx,idy,idz-1))
kraus's avatar
kraus committed
681
        F = -1;
kraus's avatar
kraus committed
682 683

    if(!isInside(idx,idy,idz+1))
kraus's avatar
kraus committed
684
        B = -1;
gsell's avatar
gsell committed
685 686 687 688 689 690 691 692 693 694

}


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

695
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
696 697 698 699 700
    // 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
701 702 703

    v = 2.0 * dot(quaternionVectorComponent, v) * quaternionVectorComponent
        + (quaternionScalarComponent * quaternionScalarComponent
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
704
           -  dot(quaternionVectorComponent, quaternionVectorComponent)) * v
705 706 707
        + 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}

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

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

716 717 718 719
    v(1) = 2.0 * (quaternion(1) * quaternion(2) + quaternion(0) * quaternion(3));
    v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
725 726 727 728
    v(1) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            + quaternion(2) * quaternion(2)
            - quaternion(3) * quaternion(3));
729 730 731 732

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
738 739 740 741
    v(2) = (quaternion(0) * quaternion(0)
            - quaternion(1) * quaternion(1)
            - quaternion(2) * quaternion(2)
            + quaternion(3) * quaternion(3));
742 743

}
744 745 746 747 748 749 750 751

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