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

15
#ifdef HAVE_SAAMG_SOLVER
16 17 18 19
#include <map>
#include <cmath>
#include <iostream>
#include <assert.h>
gsell's avatar
gsell committed
20 21

#include "ArbitraryDomain.h"
22

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
ArbitraryDomain::ArbitraryDomain(
	BoundaryGeometry * bgeom, 
	Vector_t nr, 
	Vector_t hr, 
	std::string interpl) {

    	bgeom_m  = bgeom;
    	Geo_mincoords_m = bgeom->getmincoords();
    	Geo_maxcoords_m = bgeom->getmaxcoords();
   
    	setNr(nr);
    	setHr(hr);
    	setMinMaxZ(Geo_mincoords_m[2],Geo_maxcoords_m[2]);
   	startId = 0;

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

46 47
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
48 49
}

50
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
gsell's avatar
gsell committed
51 52 53

    setHr(hr);

54 55 56 57 58 59 60
    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;
    
61 62
    hasGeometryChanged_m = true;

63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
    IntersectLoX.clear();
    IntersectHiX.clear();
    IntersectLoY.clear();
    IntersectHiY.clear();
    IntersectLoZ.clear();
    IntersectHiZ.clear();


    //calculate intersection 
    Vector_t P, dir, I;
    for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
	 P[2] = (idz - (nr[2]-1)/2.0)*hr[2];

	 for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
	     P[1] = (idy - (nr[1]-1)/2.0)*hr[1];

    	     for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
	       	  P[0] = (idx - (nr[0]-1)/2.0)*hr[0];

       			std::tuple<int, int, int> pos(idx, idy, idz);

	        	dir = Vector_t(0,0,1);
		        if (bgeom_m->intersectRayBoundary(P, dir, I))
       	      		 IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));

	        	dir = Vector_t(0,0,-1);
		        if (bgeom_m->intersectRayBoundary(P, dir, I))
       	      		 IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));

	        	dir = Vector_t(0,1,0);
		        if (bgeom_m->intersectRayBoundary(P, dir, I))
       	      		 IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));

	        	dir = Vector_t(0,-1,0);
		        if (bgeom_m->intersectRayBoundary(P, dir, I))
       	      		 IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));

	        	dir = Vector_t(1,0,0);
		        if (bgeom_m->intersectRayBoundary(P, dir, I))
       	      		 IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));

	        	dir = Vector_t(-1,0,0);
		        if (bgeom_m->intersectRayBoundary(P, dir, I))
       	      		 IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
		}
	 }
     }

    //number of ghost nodes to the right 
    int znumGhostNodesRight = 0;
    if(zGhostOffsetRight != 0) {
        for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
            for(int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
                if(isInside(idx, idy, localId[2].last() + zGhostOffsetRight))
                    znumGhostNodesRight++;
gsell's avatar
gsell committed
118 119
            }
        }
120
    }
gsell's avatar
gsell committed
121

122 123 124 125 126 127 128
    //number of ghost nodes to the left 
    int znumGhostNodesLeft = 0;
    if(zGhostOffsetLeft != 0) {
        for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
            for(int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
                    znumGhostNodesLeft++;
gsell's avatar
gsell committed
129 130
            }
        }
131
    }
gsell's avatar
gsell committed
132

133 134 135 136 137 138 139
    //number of ghost nodes to the right 
    int ynumGhostNodesRight = 0;
    if(yGhostOffsetRight != 0) {
        for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
            for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
                if(isInside(idx, localId[1].last() + yGhostOffsetRight, idz))
                    ynumGhostNodesRight++;
gsell's avatar
gsell committed
140 141 142 143
            }
        }
    }

144 145 146 147 148 149 150 151 152 153
    //number of ghost nodes to the left 
    int ynumGhostNodesLeft = 0;
    if(yGhostOffsetLeft != 0) {
        for(int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
            for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
                if(isInside(idx, localId[1].first() - yGhostOffsetLeft, idz))
                    ynumGhostNodesLeft++;
            }
        }
    }
gsell's avatar
gsell committed
154 155


156 157 158 159 160 161 162
    //number of ghost nodes to the right 
    int xnumGhostNodesRight = 0;
    if(xGhostOffsetRight != 0) {
	for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
            for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
                if(isInside(localId[0].last() + xGhostOffsetRight, idy, idz))
                    xnumGhostNodesRight++;
gsell's avatar
gsell committed
163 164 165 166
            }
        }
    }

167 168 169 170 171 172 173 174 175 176
    //number of ghost nodes to the left 
    int xnumGhostNodesLeft = 0;
    if(xGhostOffsetLeft != 0) {
       	for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
            for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
                if(isInside(localId[0].first() - xGhostOffsetLeft, idy, idz))
                    xnumGhostNodesLeft++;
            }
        }
    }
gsell's avatar
gsell committed
177
    //xy points in z plane
178 179
    int numxy; 
    int numtotalxy = 0;
gsell's avatar
gsell committed
180

181
    numXY.clear();
gsell's avatar
gsell committed
182

183 184 185 186 187 188
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
	numxy =0;
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
            for (int idx =localId[0].first(); idx <= localId[0].last(); idx++) {
                if (isInside(idx, idy, idz))
                   numxy++;
gsell's avatar
gsell committed
189 190
            }
        }
191
        numtotalxy += numxy;
gsell's avatar
gsell committed
192 193
    }

194 195
    startId = 0;
    MPI_Scan(&numtotalxy, &startId, 1, MPI_INTEGER, MPI_SUM, Ippl::getComm());
gsell's avatar
gsell committed
196

197
    startId -= numtotalxy;
gsell's avatar
gsell committed
198 199 200 201

    //build up index and coord map
    IdxMap.clear();
    CoordMap.clear();
202 203 204 205 206 207 208 209 210 211 212 213 214
    
    register int id = startId - xnumGhostNodesLeft - ynumGhostNodesLeft - znumGhostNodesLeft;
     for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
    	 for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
    	     for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
	            if (isInside(idx, idy, idz)) {
                    IdxMap[toCoordIdx(idx, idy, idz)] = id;
                    CoordMap[id] = toCoordIdx(idx, idy, idz);
                    id++;
                 }
             }
         }
     }
gsell's avatar
gsell committed
215 216
}

217 218 219
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
	return (idz * nr[1] + idy) * nr[0]  + idx;
gsell's avatar
gsell committed
220 221
}

222 223 224 225 226 227 228
// Conversion from (x,y,z) to index on the 3D grid
int ArbitraryDomain::getIdx(int idx, int idy, int idz) {
	if(isInside(idx, idy, idz) && idx>=0 && idy >=0 && idz >=0 )
       	       	 	return IdxMap[toCoordIdx(idx, idy, idz)];
    	else
        	return -1;
} 
gsell's avatar
gsell committed
229

230 231
// Conversion from a 3D index to (x,y,z)
inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) {
gsell's avatar
gsell committed
232

233
    int id = CoordMap[idxyz];
gsell's avatar
gsell committed
234

235 236 237 238 239
    idx = id % (int)nr[0];
    id /= nr[0];
    idy = id % (int)nr[1];
    id /= nr[1];
    idz = id;
gsell's avatar
gsell committed
240 241
}

242 243 244 245 246
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
 /* Expensive computation to check */
 /*
	Vector_t P0, P;
    	P0 = Vector_t(0,0,0);
gsell's avatar
gsell committed
247

248 249 250 251 252
	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[z];
	return (bgeom_m->fastIsInside(P0, P) % 2 == 0); 
 */
gsell's avatar
gsell committed
253
    bool ret = false;
254 255 256
    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];
gsell's avatar
gsell committed
257

258 259 260 261
    int    countHz, countLz, countHy, countLy, countHx, countLx;
    std::multimap < std::tuple<int, int, int>, double >::iterator itrHz, itrLz;
    std::multimap < std::tuple<int, int, int>, double >::iterator itrHy, itrLy;
    std::multimap < std::tuple<int, int, int>, double >::iterator itrHx, itrLx;
gsell's avatar
gsell committed
262

263 264
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
             
gsell's avatar
gsell committed
265
    //check if z is inside with x,y coords
266 267 268 269 270 271 272 273 274
    itrHz = IntersectHiZ.find(coordxyz);
    itrLz = IntersectLoZ.find(coordxyz);

    countHz = IntersectHiZ.count(coordxyz);
    countLz = IntersectLoZ.count(coordxyz);

     //check if y is inside with x,z coords
    itrHy = IntersectHiY.find(coordxyz);
    itrLy = IntersectLoY.find(coordxyz);
gsell's avatar
gsell committed
275

276 277
    countHy = IntersectHiY.count(coordxyz);
    countLy = IntersectLoY.count(coordxyz);
gsell's avatar
gsell committed
278

279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
    //check if x is inside with y,z coords
    itrHx = IntersectHiX.find(coordxyz);
    itrLx = IntersectLoX.find(coordxyz);

    countHx = IntersectHiX.count(coordxyz);
    countLx = IntersectLoX.count(coordxyz);

    if(countHz == 1 && countLz == 1)
        ret = (cz <= itrHz->second) && (cz >= itrLz->second);
    else if(countHy == 1 && countLy == 1)
        ret = ret && (cy <= itrHy->second) && (cy >= itrLy->second);
    else if(countHx == 1 && countLx == 1)
        ret = ret && (cx <= itrHx->second) && (cx >= itrLx->second);
       
    return ret; 
gsell's avatar
gsell committed
294 295 296
}

int ArbitraryDomain::getNumXY(int z) {
297 298 299
    
	return numXY[z];
}
gsell's avatar
gsell committed
300 301


302 303 304 305 306
void ArbitraryDomain::getBoundaryStencil(int idxyz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
    int x = 0, y = 0, z = 0;

    getCoord(idxyz, x, y, z);
    getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
gsell's avatar
gsell committed
307 308 309 310
}

void ArbitraryDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {

311
    scaleFactor = 1.0;
312
   // determine which interpolation method we use for points near the boundary
313 314 315 316 317 318 319 320 321 322
    switch(interpolationMethod){
    	case CONSTANT:
        	ConstantInterpolation(x,y,z,W,E,S,N,F,B,C,scaleFactor);
        	break;
    	case LINEAR:
        //	LinearInterpolation(x,y,z,W,E,S,N,F,B,C,scaleFactor);
        	break;
    	case QUADRATIC:
	    //  QuadraticInterpolation(x,y,z,W,E,S,N,F,B,C,scaleFactor);
        	break;
323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338
    }

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

void ArbitraryDomain::ConstantInterpolation(int x, int y, int z, double& W, double& E, double& S, double& N, double& F, double& B, double& C, double &scaleFactor) {

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

339
    if(!isInside(x+1,y,z)) 
340 341
        E = 0.0;

342
    if(!isInside(x-1,y,z))
343 344
        W = 0.0;

345
    if(!isInside(x,y+1,z))
346 347
        N = 0.0;

348
    if(!isInside(x,y-1,z)) 
349
        S = 0.0;
350 351 352
    
    if(!isInside(x,y,z-1)) 
	F = 0.0;	
353

354 355
    if(!isInside(x,y,z+1)) 
	B = 0.0;	
356 357 358 359 360 361 362 363 364 365

}
/*
void ArbitraryDomain::LinearInterpolation(int x, int y, int z, double& W, double& E, double& S, double& N, double& F, double& B, double& C, double &scaleFactor) 
{

    scaleFactor = 1.0;

    double dx=-1, dy=-1;

366 367
    double cx = x * hr[0];
    double cy = y * hr[1];
368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456

    std::multimap< std::pair<int, int>, double >::iterator it;
    std::pair< std::multimap< std::pair<int, int>, double>::iterator, std::multimap< std::pair<int, int>, double>::iterator > ret;

    std::pair<int, int> coordyz(y, z);
    ret = IntersectXDir.equal_range(coordyz);
    for(it = ret.first; it != ret.second; ++it) {
        if(fabs(it->second - cx) < hr[0]) {
            dx = it->second;
            break;
        }
    }

    std::pair<int, int> coordxz(x, z);
    ret = IntersectYDir.equal_range(coordxz);
    for(it = ret.first; it != ret.second; ++it) {
        if(fabs(it->second - cy) < hr[1]) {
            dy = it->second;
            break;
        }
    }


    double dw=hr[0];
    double de=hr[0];
    double dn=hr[1];
    double ds=hr[1];
    C = 0.0;

    // we are a right boundary point
    if(dx >= 0 && dx > cx)
    {	
        C += 1/((dx-cx)*de);
        E = 0.0;
    } else {
        C += 1/(de*de);
        E = -1/(de*de);
    }

    // we are a left boundary point
    if(dx <= 0 && dx < cx)
    {
        C += 1/((cx-dx)*dw);
        W = 0.0;
    } else {
        C += 1/(dw*dw);
        W = -1/(dw*dw);
    }

    // we are a upper boundary point
    if(dy >= 0 && dy > cy)
    {
        C += 1/((dy-cy)*dn);
        N = 0.0;
    } else {
        C += 1/(dn*dn);
        N = -1/(dn*dn);
    }

    // we are a lower boundary point
    if(dy <= 0 && dy < cy)
    {
        C += 1/((cy-dy)*ds);
        S = 0.0;
    } else {
        C += 1/(ds*ds);
        S = -1/(ds*ds);
    }

    F = -1/(hr[2]*hr[2]); 
    B = -1/(hr[2]*hr[2]);

    //XXX: In stand-alone only Dirichlet for validation purposes
    if(z == 0 || z == nr[2]-1) {

        // Dirichlet
        C += 2/hr[2]*1/hr[2];

        //C += 1/hr[2]*1/hr[2];

        // case where we are on the Neumann BC in Z-direction
        // where we distinguish two cases  
        if(z == 0)
            F = 0.0;
        else 
            B = 0.0;

        //for test no neumann 
        //C += 2/((hr[2]*nr[2]/2.0) * hr[2]);
457 458
        //
        //   double d = hr[2]*(nr[2])/2;
459 460 461
        //   C += 2/(d * hr[2]);


462 463 464 465 466
        ////neumann stuff
        //W /= 2.0;
        //E /= 2.0;
        //N /= 2.0;
        //S /= 2.0;
467 468
        //C /= 2.0;

469
        scaleFactor *= 0.5;
470 471 472 473

    } else 
        C += 2*1/hr[2]*1/hr[2];
}
474
        */
gsell's avatar
gsell committed
475

476
void ArbitraryDomain::getNeighbours(int id, int &W, int &E, int &S, int &N, int &F, int &B) {
gsell's avatar
gsell committed
477 478 479

    int x = 0, y = 0, z = 0;

480
    getCoord(id, x, y, z);
gsell's avatar
gsell committed
481 482 483
    getNeighbours(x, y, z, W, E, S, N, F, B);
}

484
void ArbitraryDomain::getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) {
gsell's avatar
gsell committed
485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521

    if(x > 0)
        W = getIdx(x - 1, y, z);
    else
        W = -1;
    if(x < nr[0] - 1)
        E = getIdx(x + 1, y, z);
    else
        E = -1;

    if(y < nr[1] - 1)
        N = getIdx(x, y + 1, z);
    else
        N = -1;
    if(y > 0)
        S = getIdx(x, y - 1, z);
    else
        S = -1;

    if(z > 0)
        F = getIdx(x, y, z - 1);
    else
        F = -1;
    if(z < nr[2] - 1)
        B = getIdx(x, y, z + 1);
    else
        B = -1;

}


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

522
#endif //#ifdef HAVE_SAAMG_SOLVER