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>
gsell's avatar
gsell committed
38

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

    // TODO: THis needs to be made into OPTION of the geometry.
49
    // A user defined point that is INSIDE with 100% certainty. -DW
50 51 52 53
    bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m);
    if (have_inside_pt == false) {
        globalInsideP0_m = Vector_t(0.0, 0.0, -0.13);
    }
54
    setNr(nr);
55
    for(int i=0; i<3; i++)
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
56
        Geo_hr_m[i] = (maxCoords_m[i] - minCoords_m[i])/nr[i];
57
    setHr(Geo_hr_m);
58 59 60 61 62 63 64 65 66

    startId = 0;

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

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

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

75 76
    setHr(hr);

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

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

84 85
    hasGeometryChanged_m = true;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

206
    setHr(hr);
207

208 209 210 211
    globalMeanR_m = getGlobalMeanR();

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

    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
221

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
316
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
kraus's avatar
kraus committed
317 318 319 320
                        //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
321
                    } else {
322
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
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 328
                    //rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
                    dir = Vector_t(1, 0, 0);
329

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

}


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

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

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

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

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

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

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

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

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

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

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

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

}
743 744 745 746 747 748 749 750

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