BoxCornerDomain.cpp 12.6 KB
Newer Older
1 2 3 4
//
// Class BoxCornerDomain
//   :FIXME: add brief description
//
5
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
6
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
7
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
8
// All rights reserved
9
//
10 11 12 13 14
// 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
15 16 17 18 19 20 21 22 23 24
//
// 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/>.
25
//
kraus's avatar
kraus committed
26
#include "Solvers/BoxCornerDomain.h"
27
#include "Utilities/OpalException.h"
kraus's avatar
kraus committed
28

gsell's avatar
gsell committed
29
#include <map>
30
#include <string>
gsell's avatar
gsell committed
31 32 33 34 35
#include <cmath>
#include <iostream>

//FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)

36
BoxCornerDomain::BoxCornerDomain(double A, double B, double C,
37
                                 double L1, double L2, IntVector_t nr, Vector_t hr,
38
                                 std::string interpl)
39
    : RegularDomain(nr, hr, interpl)
40 41 42
{
    setRangeMin(Vector_t(-A, -B, L1));
    setRangeMax(Vector_t( A,  B, L1 + L2));
gsell's avatar
gsell committed
43
    C_m = C;
44 45 46

    throw OpalException("BoxCornerDomain::BoxCornerDomain()",
                        "This domain is currently not supported!");
gsell's avatar
gsell committed
47 48 49 50 51 52 53 54 55 56 57 58
}

BoxCornerDomain::~BoxCornerDomain() {
    //nothing so far
}


// for this geometry we only have to calculate the intersection on
// all x-y-planes
// for the moment we center the box corner geometry around the center of the grid
// hr holds the grid-spacings (boundary ellipse embedded in hr-grid)

59
void BoxCornerDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
gsell's avatar
gsell committed
60 61

    //there is nothing to be done if the mesh spacings have not changed
frey_m's avatar
frey_m committed
62
    //    if(hr_m[0] == getHr()[0] && hr_m[1] == getHr()[1] && hr_m[2] == getHr()[2]) {
gsell's avatar
gsell committed
63 64 65 66 67 68 69 70 71 72
    //      hasGeometryChanged_m = false;
    //      return;
    //  }

    setHr(hr);
    hasGeometryChanged_m = true;

    double bL= getB(getMinZ());
    double bH= getB(getMaxZ());

73
    actBMin_m = getYRangeMin();
74
    actBMax_m = std::max(bL,bH);
gsell's avatar
gsell committed
75 76 77 78

    //reset number of points inside domain

    // clear previous coordinate maps
79 80
    idxMap_m.clear();
    coordMap_m.clear();
gsell's avatar
gsell committed
81 82 83 84 85
    //clear previous intersection points
    IntersectYDir.clear();
    IntersectXDir.clear();

    // build a index and coordinate map
86 87
    int idx = 0;
    int x, y, z;
frey_m's avatar
frey_m committed
88 89 90
    for(x = 0; x < nr_m[0]; x++) {
        for(y = 0; y < nr_m[1]; y++) {
            for(z = 0; z < nr_m[2]; z++) {
gsell's avatar
gsell committed
91
                if(isInside(x, y, z)) {
92 93
                    idxMap_m[toCoordIdx(x, y, z)] = idx;
                    coordMap_m[idx++] = toCoordIdx(x, y, z);
gsell's avatar
gsell committed
94 95 96 97 98 99 100
                }
            }
        }
    }

    //XXX: calculate intersection on the fly
    /*
101
    switch(interpolationMethod_m) {
gsell's avatar
gsell committed
102 103 104 105 106 107 108 109

    case CONSTANT:
        break;
    case LINEAR:
    case QUADRATIC:

        // calculate intersection

frey_m's avatar
frey_m committed
110
        for(int z = 0; z < nr_m[2]; z++) {
gsell's avatar
gsell committed
111

frey_m's avatar
frey_m committed
112
            for(int x = 0; x < nr_m[0]; x++) {
gsell's avatar
gsell committed
113 114
                // the x coordinate does not change in the CornerBox geometry
                std::pair<int, int> pos(x, z);
115 116
                IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, 0.5*getXRangeMax()));
                IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, 0.5*getXRangeMin()));
gsell's avatar
gsell committed
117 118
            }

frey_m's avatar
frey_m committed
119
            for(int y = 0; y < nr_m[1]; y++) {
gsell's avatar
gsell committed
120
                std::pair<int, int> pos(y, z);
frey_m's avatar
frey_m committed
121
                double yt = getB(z*hr_m[2]);
122
                double yb = 0.5*getYRangeMin();
gsell's avatar
gsell committed
123 124 125 126 127 128 129 130
                IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, yt));
                IntersectXDir.insert(std::pair< std::pair<int, int>, double >(pos, yb));
            }
        }
    }
    */
}

131
void BoxCornerDomain::linearInterpolation(int x, int y, int z, StencilValue_t& value,
132
                                          double &scaleFactor) const
133
{
gsell's avatar
gsell committed
134 135
    scaleFactor = 1.0;

frey_m's avatar
frey_m committed
136 137
    double cx = x * hr_m[0] - (nr_m[0] - 1) * hr_m[0] / 2.0;
    double cy = y * hr_m[1] - (nr_m[1] - 1) * hr_m[1] / 2.0;
gsell's avatar
gsell committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158

    //XXX: calculate intersection on the fly
    /*
    multimap< pair<int, int>, double >::iterator it;
    pair< multimap< pair<int, int>, double>::iterator, multimap< pair<int, int>, double>::iterator > ret;

    double dx = 0.0;
    std::pair<int, int> coordxz(x, z);
    ret = IntersectXDir.equal_range(coordxz);
    if(cx < 0)
        it++;
    dx = it->second;

    double dy = 0.0;
    std::pair<int, int> coordyz(y, z);
    ret = IntersectYDir.equal_range(coordyz);
    if(cy < 0)
        it++;
    dy = it->second;
    */

frey_m's avatar
frey_m committed
159 160 161 162
    double dw = hr_m[0];
    double de = hr_m[0];
    double dn = hr_m[1];
    double ds = hr_m[1];
163
    value.center = 0.0;
gsell's avatar
gsell committed
164 165 166 167

    //we are a right boundary point
    if(!isInside(x + 1, y, z)) {
        double dx = getXIntersection(cx, z);
168 169
        value.center += 1 / ((dx - cx) * de);
        value.east = 0.0;
gsell's avatar
gsell committed
170
    } else {
171 172
        value.center += 1 / (de * de);
        value.east = -1 / (de * de);
gsell's avatar
gsell committed
173 174 175 176 177
    }

    //we are a left boundary point
    if(!isInside(x - 1, y, z)) {
        double dx = getXIntersection(cx, z);
178 179
        value.center += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
        value.west = 0.0;
gsell's avatar
gsell committed
180
    } else {
181 182
        value.center += 1 / (dw * dw);
        value.west = -1 / (dw * dw);
gsell's avatar
gsell committed
183 184 185 186 187
    }

    //we are a upper boundary point
    if(!isInside(x, y + 1, z)) {
        double dy = getYIntersection(cy, z);
188 189
        value.center += 1 / ((dy - cy) * dn);
        value.north = 0.0;
gsell's avatar
gsell committed
190
    } else {
191 192
        value.center += 1 / (dn * dn);
        value.north = -1 / (dn * dn);
gsell's avatar
gsell committed
193 194 195 196 197
    }

    //we are a lower boundary point
    if(!isInside(x, y - 1, z)) {
        double dy = getYIntersection(cy, z);
198 199
        value.center += 1 / ((std::abs(dy) - std::abs(cy)) * ds);
        value.south = 0.0;
gsell's avatar
gsell committed
200
    } else {
201 202
        value.center += 1 / (ds * ds);
        value.south = -1 / (ds * ds);
gsell's avatar
gsell committed
203 204
    }

frey_m's avatar
frey_m committed
205 206 207
    value.front = -1 / (hr_m[2] * hr_m[2]);
    value.back = -1 / (hr_m[2] * hr_m[2]);
    value.center += 2 / (hr_m[2] * hr_m[2]);
gsell's avatar
gsell committed
208 209

    // handle boundary condition in z direction
frey_m's avatar
frey_m committed
210
    if(z == 0 || z == nr_m[2] - 1) {
gsell's avatar
gsell committed
211 212 213 214

        // case where we are on the NEUMAN BC in Z-direction
        // where we distinguish two cases
        if(z == 0)
215
            value.front = 0.0;
gsell's avatar
gsell committed
216
        else
217
            value.back = 0.0;
gsell's avatar
gsell committed
218

frey_m's avatar
frey_m committed
219 220 221
        //hr_m[2]*(nr_m2[2]-1)/2 = radius
        double d = hr_m[2] * (nr_m[2] - 1) / 2;
        value.center += 2 / (d * hr_m[2]);
gsell's avatar
gsell committed
222

223 224 225 226 227 228
        value.west   /= 2.0;
        value.east   /= 2.0;
        value.north  /= 2.0;
        value.south  /= 2.0;
        value.center /= 2.0;
        scaleFactor  *= 0.5;
gsell's avatar
gsell committed
229 230 231 232 233 234

    }

}

//FIXME: this probably needs some cleanup/rewriting
235
void BoxCornerDomain::quadraticInterpolation(int x, int y, int z, StencilValue_t& value,
236
                                             double &scaleFactor) const
237
{
frey_m's avatar
frey_m committed
238 239
    double cx = (x - (nr_m[0] - 1) / 2.0) * hr_m[0];
    double cy = (y - (nr_m[1] - 1) / 2.0) * hr_m[1];
gsell's avatar
gsell committed
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263

    double dx = getXIntersection(cx, z);
    double dy = getYIntersection(cy, z);

    //XXX: calculate intersection on the fly
    /*
    multimap< pair<int, int>, double >::iterator it;
    pair< multimap< pair<int, int>, double>::iterator, multimap< pair<int, int>, double>::iterator > ret;

    double dx = 0.0;
    std::pair<int, int> coordxz(x, z);
    ret = IntersectXDir.equal_range(coordxz);
    if(cx < 0)
        it++;
    dx = it->second;

    double dy = 0.0;
    std::pair<int, int> coordyz(y, z);
    ret = IntersectYDir.equal_range(coordyz);
    if(cy < 0)
        it++;
    dy = it->second;
    */

frey_m's avatar
frey_m committed
264 265 266 267
    double dw = hr_m[0];
    double de = hr_m[0];
    double dn = hr_m[1];
    double ds = hr_m[1];
268 269 270 271 272 273 274
    value.west   = 1.0;
    value.east   = 1.0;
    value.north  = 1.0;
    value.south  = 1.0;
    value.front  = 1.0;
    value.back   = 1.0;
    value.center = 0.0;
gsell's avatar
gsell committed
275

frey_m's avatar
frey_m committed
276 277
    //TODO: = cx+hr_m[0] > dx && cx > 0
    //if((x-nr_m[0]/2.0+1)*hr_m[0] > dx && cx > 0) {
gsell's avatar
gsell committed
278 279 280 281 282
    ////we are a right boundary point
    ////if(!isInside(x+1,y,z)) {
    //de = dx-cx;
    //}

frey_m's avatar
frey_m committed
283
    //if((x-nr_m[0]/2.0-1)*hr_m[0] < dx && cx < 0) {
gsell's avatar
gsell committed
284 285 286 287 288
    ////we are a left boundary point
    ////if(!isInside(x-1,y,z)) {
    //dw = std::abs(dx)-std::abs(cx);
    //}

frey_m's avatar
frey_m committed
289
    //if((y-nr_m[1]/2.0+1)*hr_m[1] > dy && cy > 0) {
gsell's avatar
gsell committed
290 291 292 293 294
    ////we are a upper boundary point
    ////if(!isInside(x,y+1,z)) {
    //dn = dy-cy;
    //}

frey_m's avatar
frey_m committed
295
    //if((y-nr_m[1]/2.0-1)*hr_m[1] < dy && cy < 0) {
gsell's avatar
gsell committed
296 297 298 299 300
    ////we are a lower boundary point
    ////if(!isInside(x,y-1,z)) {
    //ds = std::abs(dy)-std::abs(cy);
    //}

frey_m's avatar
frey_m committed
301 302
    //TODO: = cx+hr_m[0] > dx && cx > 0
    //if((x-nr_m[0]/2.0+1)*hr_m[0] > dx && cx > 0) {
gsell's avatar
gsell committed
303 304 305
    //we are a right boundary point
    if(!isInside(x + 1, y, z)) {
        de = dx - cx;
306
        value.east = 0.0;
gsell's avatar
gsell committed
307 308
    }

frey_m's avatar
frey_m committed
309
    //if((x-nr_m[0]/2.0-1)*hr_m[0] < dx && cx < 0) {
gsell's avatar
gsell committed
310 311 312
    //we are a left boundary point
    if(!isInside(x - 1, y, z)) {
        dw = std::abs(dx) - std::abs(cx);
313
        value.west = 0.0;
gsell's avatar
gsell committed
314 315
    }

frey_m's avatar
frey_m committed
316
    //if((y-nr_m[1]/2.0+1)*hr_m[1] > dy && cy > 0) {
gsell's avatar
gsell committed
317 318 319
    //we are a upper boundary point
    if(!isInside(x, y + 1, z)) {
        dn = dy - cy;
320
        value.north = 0.0;
gsell's avatar
gsell committed
321 322
    }

frey_m's avatar
frey_m committed
323
    //if((y-nr_m[1]/2.0-1)*hr_m[1] < dy && cy < 0) {
gsell's avatar
gsell committed
324 325 326
    //we are a lower boundary point
    if(!isInside(x, y - 1, z)) {
        ds = std::abs(dy) - std::abs(cy);
327
        value.south = 0.0;
gsell's avatar
gsell committed
328 329 330
    }

    //2/dw*(dw_de)
331 332 333 334
    value.west  *= -1.0 / (dw * (dw + de));
    value.east  *= -1.0 / (de * (dw + de));
    value.north *= -1.0 / (dn * (dn + ds));
    value.south *= -1.0 / (ds * (dn + ds));
frey_m's avatar
frey_m committed
335 336
    value.front = -1 / (hr_m[2] * (hr_m[2] + hr_m[2]));
    value.back  = -1 / (hr_m[2] * (hr_m[2] + hr_m[2]));
gsell's avatar
gsell committed
337 338 339

    //TODO: problem when de,dw,dn,ds == 0
    //is NOT a regular BOUND PT
340 341
    value.center += 1 / de * 1 / dw;
    value.center += 1 / dn * 1 / ds;
frey_m's avatar
frey_m committed
342
    value.center += 1 / hr_m[2] * 1 / hr_m[2];
gsell's avatar
gsell committed
343 344 345 346 347 348 349 350 351 352
    scaleFactor = 0.5;


    //for regular gridpoints no problem with symmetry, just boundary
    //z direction is right
    //implement isLastInside(dir)
    //we have LOCAL x,y coordinates!

    /*
       if(dw != 0 && !wIsB)
frey_m's avatar
frey_m committed
353
       W = -1/dw * (dn+ds) * 2*hr_m[2];
gsell's avatar
gsell committed
354 355 356
       else
       W = 0;
       if(de != 0 && !eIsB)
frey_m's avatar
frey_m committed
357
       E = -1/de * (dn+ds) * 2*hr_m[2];
gsell's avatar
gsell committed
358 359 360
       else
       E = 0;
       if(dn != 0 && !nIsB)
frey_m's avatar
frey_m committed
361
       N = -1/dn * (dw+de) * 2*hr_m[2];
gsell's avatar
gsell committed
362 363 364
       else
       N = 0;
       if(ds != 0 && !sIsB)
frey_m's avatar
frey_m committed
365
       S = -1/ds * (dw+de) * 2*hr_m[2];
gsell's avatar
gsell committed
366 367
       else
       S = 0;
frey_m's avatar
frey_m committed
368 369
       F = -(dw+de)*(dn+ds)/hr_m[2];
       B = -(dw+de)*(dn+ds)/hr_m[2];
gsell's avatar
gsell committed
370 371 372
       */

    //if(dw != 0)
frey_m's avatar
frey_m committed
373
    //W = -2*hr_m[2]*(dn+ds)/dw;
gsell's avatar
gsell committed
374 375 376
    //else
    //W = 0;
    //if(de != 0)
frey_m's avatar
frey_m committed
377
    //E = -2*hr_m[2]*(dn+ds)/de;
gsell's avatar
gsell committed
378 379 380
    //else
    //E = 0;
    //if(dn != 0)
frey_m's avatar
frey_m committed
381
    //N = -2*hr_m[2]*(dw+de)/dn;
gsell's avatar
gsell committed
382 383 384
    //else
    //N = 0;
    //if(ds != 0)
frey_m's avatar
frey_m committed
385
    //S = -2*hr_m[2]*(dw+de)/ds;
gsell's avatar
gsell committed
386 387
    //else
    //S = 0;
frey_m's avatar
frey_m committed
388 389
    //F = -(dw+de)*(dn+ds)/hr_m[2];
    //B = -(dw+de)*(dn+ds)/hr_m[2];
gsell's avatar
gsell committed
390 391 392

    //// RHS scaleFactor for current 3D index
    //// Factor 0.5 results from the SW/quadratic extrapolation
frey_m's avatar
frey_m committed
393
    //scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr_m[2]);
gsell's avatar
gsell committed
394 395 396 397 398 399 400 401 402 403 404 405 406 407

    // catch the case where a point lies on the boundary
    //FIXME: do this more elegant!
    //double m1 = dw*de;
    //double m2 = dn*ds;
    //if(de == 0)
    //m1 = dw;
    //if(dw == 0)
    //m1 = de;
    //if(dn == 0)
    //m2 = ds;
    //if(ds == 0)
    //m2 = dn;
    ////XXX: dn+ds || dw+de can be 0
frey_m's avatar
frey_m committed
408 409
    ////C = 2*(dn+ds)*(dw+de)/hr_m[2];
    //C = 2/hr_m[2];
gsell's avatar
gsell committed
410 411 412 413 414
    //if(dw != 0 || de != 0)
    //C *= (dw+de);
    //if(dn != 0 || ds != 0)
    //C *= (dn+ds);
    //if(dw != 0 || de != 0)
frey_m's avatar
frey_m committed
415
    //C += (2*hr_m[2])*(dn+ds)*(dw+de)/m1;
gsell's avatar
gsell committed
416
    //if(dn != 0 || ds != 0)
frey_m's avatar
frey_m committed
417
    //C += (2*hr_m[2])*(dw+de)*(dn+ds)/m2;
gsell's avatar
gsell committed
418 419

    //handle Neumann case
frey_m's avatar
frey_m committed
420
    //if(z == 0 || z == nr_m[2]-1) {
gsell's avatar
gsell committed
421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437

    //if(z == 0)
    //F = 0.0;
    //else
    //B = 0.0;

    ////neumann stuff
    //W = W/2.0;
    //E = E/2.0;
    //N = N/2.0;
    //S = S/2.0;
    //C /= 2.0;

    //scaleFactor /= 2.0;
    //}

    // handle boundary condition in z direction
frey_m's avatar
frey_m committed
438
    if(z == 0 || z == nr_m[2] - 1) {
gsell's avatar
gsell committed
439 440 441 442

        // case where we are on the NEUMAN BC in Z-direction
        // where we distinguish two cases
        if(z == 0)
443
            value.front = 0.0;
gsell's avatar
gsell committed
444
        else
445
            value.back = 0.0;
gsell's avatar
gsell committed
446

frey_m's avatar
frey_m committed
447 448 449 450
        //value.center += 2/((hr_m[2]*(nr_m[2]-1)/2.0) * hr_m[2]);
        //hr_m[2]*(nr_m2[2]-1)/2 = radius
        double d = hr_m[2] * (nr_m[2] - 1) / 2;
        value.center += 2 / (d * hr_m[2]);
451 452 453 454 455 456 457

        value.west   /= 2.0;
        value.east   /= 2.0;
        value.north  /= 2.0;
        value.south  /= 2.0;
        value.center /= 2.0;
        scaleFactor  /= 2.0;
gsell's avatar
gsell committed
458 459 460

    }
}