OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
EllipticDomain.cpp
Go to the documentation of this file.
1 #ifdef HAVE_SAAMG_SOLVER
2 
3 #include "Solvers/EllipticDomain.h"
4 
5 #include <map>
6 #include <cmath>
7 #include <iostream>
8 #include <assert.h>
9 
10 //FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)
11 
12 
13 EllipticDomain::EllipticDomain(Vector_t nr, Vector_t hr) {
14  setNr(nr);
15  setHr(hr);
16 }
17 
18 EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl) {
19  SemiMajor = semimajor;
20  SemiMinor = semiminor;
21  setNr(nr);
22  setHr(hr);
23 
24  if(interpl == "CONSTANT")
25  interpolationMethod = CONSTANT;
26  else if(interpl == "LINEAR")
27  interpolationMethod = LINEAR;
28  else if(interpl == "QUADRATIC")
29  interpolationMethod = QUADRATIC;
30 }
31 
32 EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl) {
33  SemiMajor = bgeom->getA();
34  SemiMinor = bgeom->getB();
35  setMinMaxZ(bgeom->getS(), bgeom->getLength());
36  setNr(nr);
37  setHr(hr);
38 
39  if(interpl == "CONSTANT")
40  interpolationMethod = CONSTANT;
41  else if(interpl == "LINEAR")
42  interpolationMethod = LINEAR;
43  else if(interpl == "QUADRATIC")
44  interpolationMethod = QUADRATIC;
45 }
46 
47 EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl) {
48  SemiMajor = bgeom->getA();
49  SemiMinor = bgeom->getB();
50  setMinMaxZ(bgeom->getS(), bgeom->getLength());
51  Vector_t hr_m;
52  hr_m[0] = (getXRangeMax()-getXRangeMin())/nr[0];
53  hr_m[1] = (getYRangeMax()-getYRangeMin())/nr[1];
54  hr_m[2] = (getZRangeMax()-getZRangeMin())/nr[2];
55  setHr(hr_m);
56 
57  if(interpl == "CONSTANT")
58  interpolationMethod = CONSTANT;
59  else if(interpl == "LINEAR")
60  interpolationMethod = LINEAR;
61  else if(interpl == "QUADRATIC")
62  interpolationMethod = QUADRATIC;
63 }
64 
65 EllipticDomain::~EllipticDomain() {
66  //nothing so far
67 }
68 
69 
70 // for this geometry we only have to calculate the intersection on
71 // one x-y-plane
72 // for the moment we center the ellipse around the center of the grid
73 void EllipticDomain::compute(Vector_t hr){
74  //there is nothing to be done if the mesh spacings have not changed
75  if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
76  hasGeometryChanged_m = false;
77  return;
78  }
79  setHr(hr);
80  hasGeometryChanged_m = true;
81  //reset number of points inside domain
82  nxy_m = 0;
83 
84  // clear previous coordinate maps
85  IdxMap.clear();
86  CoordMap.clear();
87  //clear previous intersection points
88  IntersectYDir.clear();
89  IntersectXDir.clear();
90 
91  // build a index and coordinate map
92  int idx = 0;
93  int x, y;
94 
95  for(x = 0; x < nr[0]; x++) {
96  for(y = 0; y < nr[1]; y++) {
97  if(isInside(x, y, 1)) {
98  IdxMap[toCoordIdx(x, y)] = idx;
99  CoordMap[idx++] = toCoordIdx(x, y);
100  nxy_m++;
101  }
102 
103  }
104  }
105 
106  switch(interpolationMethod) {
107  case CONSTANT:
108  break;
109  case LINEAR:
110  case QUADRATIC:
111 
112  double smajsq = SemiMajor * SemiMajor;
113  double sminsq = SemiMinor * SemiMinor;
114  double yd = 0.0;
115  double xd = 0.0;
116  double pos = 0.0;
117  double mx = (nr[0] - 1) * hr[0] / 2.0;
118  double my = (nr[1] - 1) * hr[1] / 2.0;
119 
120  //calculate intersection with the ellipse
121  for(x = 0; x < nr[0]; x++) {
122  pos = x * hr[0] - mx;
123  if (pos <= -SemiMajor || pos >= SemiMajor)
124  {
125  IntersectYDir.insert(std::pair<int, double>(x, 0));
126  IntersectYDir.insert(std::pair<int, double>(x, 0));
127  }else{
128  yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
129  IntersectYDir.insert(std::pair<int, double>(x, yd));
130  IntersectYDir.insert(std::pair<int, double>(x, -yd));
131  }
132 
133  }
134 
135  for(y = 0; y < nr[1]; y++) {
136  pos = y * hr[1] - my;
137  if (pos <= -SemiMinor || pos >= SemiMinor)
138  {
139  IntersectXDir.insert(std::pair<int, double>(y, 0));
140  IntersectXDir.insert(std::pair<int, double>(y, 0));
141  }else{
142  xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
143  IntersectXDir.insert(std::pair<int, double>(y, xd));
144  IntersectXDir.insert(std::pair<int, double>(y, -xd));
145  }
146  }
147  }
148 }
149 
150 void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
151  //there is nothing to be done if the mesh spacings have not changed
152  if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
153  hasGeometryChanged_m = false;
154  return;
155  }
156  setHr(hr);
157  hasGeometryChanged_m = true;
158  //reset number of points inside domain
159  nxy_m = 0;
160 
161  // clear previous coordinate maps
162  IdxMap.clear();
163  CoordMap.clear();
164  //clear previous intersection points
165  IntersectYDir.clear();
166  IntersectXDir.clear();
167 
168  // build a index and coordinate map
169  int idx = 0;
170  int x, y;
171 
172  for(x = localId[0].first(); x<= localId[0].last(); x++) {
173  for(y = localId[1].first(); y <= localId[1].last(); y++) {
174  if(isInside(x, y, 1)) {
175  IdxMap[toCoordIdx(x, y)] = idx;
176  CoordMap[idx++] = toCoordIdx(x, y);
177  nxy_m++;
178  }
179  }
180  }
181 
182  switch(interpolationMethod) {
183  case CONSTANT:
184  break;
185  case LINEAR:
186  case QUADRATIC:
187 
188  double smajsq = SemiMajor * SemiMajor;
189  double sminsq = SemiMinor * SemiMinor;
190  double yd = 0.0;
191  double xd = 0.0;
192  double pos = 0.0;
193  double mx = (nr[0] - 1) * hr[0] / 2.0;
194  double my = (nr[1] - 1) * hr[1] / 2.0;
195 
196  //calculate intersection with the ellipse
197  for(x = localId[0].first(); x <= localId[0].last(); x++) {
198  pos = x * hr[0] - mx;
199  if (pos <= -SemiMajor || pos >= SemiMajor)
200  {
201  IntersectYDir.insert(std::pair<int, double>(x, 0));
202  IntersectYDir.insert(std::pair<int, double>(x, 0));
203  }else{
204  yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
205  IntersectYDir.insert(std::pair<int, double>(x, yd));
206  IntersectYDir.insert(std::pair<int, double>(x, -yd));
207  }
208 
209  }
210 
211  for(y = localId[0].first(); y < localId[1].last(); y++) {
212  pos = y * hr[1] - my;
213  if (pos <= -SemiMinor || pos >= SemiMinor)
214  {
215  IntersectXDir.insert(std::pair<int, double>(y, 0));
216  IntersectXDir.insert(std::pair<int, double>(y, 0));
217  }else{
218  xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
219  IntersectXDir.insert(std::pair<int, double>(y, xd));
220  IntersectXDir.insert(std::pair<int, double>(y, -xd));
221  }
222  }
223  }
224 }
225 
226 void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
227  scaleFactor = 1.0;
228 
229  // determine which interpolation method we use for points near the boundary
230  switch(interpolationMethod) {
231  case CONSTANT:
232  constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
233  break;
234  case LINEAR:
235  linearInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
236  break;
237  case QUADRATIC:
238  quadraticInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
239  break;
240  }
241 
242  // stencil center value has to be positive!
243  assert(C > 0);
244 }
245 
246 void EllipticDomain::getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
247 
248 
249  int x = 0, y = 0, z = 0;
250  getCoord(idx, x, y, z);
251  getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
252 }
253 
254 
255 void EllipticDomain::getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B) {
256 
257  int x = 0, y = 0, z = 0;
258  getCoord(idx, x, y, z);
259  getNeighbours(x, y, z, W, E, S, N, F, B);
260 
261 }
262 
263 void EllipticDomain::getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) {
264 
265  if(x > 0)
266  W = getIdx(x - 1, y, z);
267  else
268  W = -1;
269  if(x < nr[0] - 1)
270  E = getIdx(x + 1, y, z);
271  else
272  E = -1;
273 
274  if(y < nr[1] - 1)
275  N = getIdx(x, y + 1, z);
276  else
277  N = -1;
278  if(y > 0)
279  S = getIdx(x, y - 1, z);
280  else
281  S = -1;
282 
283  if(z > 0)
284  F = getIdx(x, y, z - 1);
285  else
286  F = -1;
287  if(z < nr[2] - 1)
288  B = getIdx(x, y, z + 1);
289  else
290  B = -1;
291 
292 }
293 
294 void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV, double &EV, double &SV, double &NV, double &FV, double &BV, double &CV, double &scaleFactor) {
295 
296  scaleFactor = 1.0;
297 
298  WV = -1/(hr[0]*hr[0]);
299  EV = -1/(hr[0]*hr[0]);
300  NV = -1/(hr[1]*hr[1]);
301  SV = -1/(hr[1]*hr[1]);
302  FV = -1/(hr[2]*hr[2]);
303  BV = -1/(hr[2]*hr[2]);
304  CV = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
305 
306  // we are a right boundary point
307  if(!isInside(x + 1, y, z))
308  EV= 0.0;
309 
310  // we are a left boundary point
311  if(!isInside(x - 1, y, z))
312  WV = 0.0;
313 
314  // we are a upper boundary point
315  if(!isInside(x, y + 1, z))
316  NV = 0.0;
317 
318  // we are a lower boundary point
319  if(!isInside(x, y - 1, z))
320  SV = 0.0;
321 
322  if(z == 1 || z == nr[2] - 2) {
323 
324  // case where we are on the Robin BC in Z-direction
325  // where we distinguish two cases
326  // IFF: this values should not matter because they
327  // never make it into the discretization matrix
328  if(z == 1)
329  FV = 0.0;
330  else
331  BV = 0.0;
332 
333  // add contribution of Robin discretization to center point
334  // d the distance between the center of the bunch and the boundary
335  //double cx = (x-(nr[0]-1)/2)*hr[0];
336  //double cy = (y-(nr[1]-1)/2)*hr[1];
337  //double cz = hr[2]*(nr[2]-1);
338  //double d = sqrt(cx*cx+cy*cy+cz*cz);
339  double d = hr[2] * (nr[2] - 1) / 2;
340  CV += 2 / (d * hr[2]);
341  //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
342 
343  // scale all stencil-points in z-plane with 0.5 (Robin discretization)
344  WV /= 2.0;
345  EV /= 2.0;
346  NV /= 2.0;
347  SV /= 2.0;
348  CV /= 2.0;
349  scaleFactor *= 0.5;
350  }
351 
352 }
353 
354 void EllipticDomain::linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
355 
356  scaleFactor = 1.0;
357 
358  double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
359  double cy = y * hr[1] - (nr[1] - 1) * hr[1] / 2.0;
360 
361  double dx = 0.0;
362  std::multimap<int, double>::iterator it = IntersectXDir.find(y);
363  if(cx < 0)
364  ++it;
365  dx = it->second;
366 
367  double dy = 0.0;
368  it = IntersectYDir.find(x);
369  if(cy < 0)
370  ++it;
371  dy = it->second;
372 
373 
374  double dw = hr[0];
375  double de = hr[0];
376  double dn = hr[1];
377  double ds = hr[1];
378  C = 0.0;
379 
380  //we are a right boundary point
381  if(!isInside(x + 1, y, z)) {
382  C += 1 / ((dx - cx) * de);
383  E = 0.0;
384  } else {
385  C += 1 / (de * de);
386  E = -1 / (de * de);
387  }
388 
389  //we are a left boundary point
390  if(!isInside(x - 1, y, z)) {
391  C += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
392  W = 0.0;
393  } else {
394  C += 1 / (dw * dw);
395  W = -1 / (dw * dw);
396  }
397 
398  //we are a upper boundary point
399  if(!isInside(x, y + 1, z)) {
400  C += 1 / ((dy - cy) * dn);
401  N = 0.0;
402  } else {
403  C += 1 / (dn * dn);
404  N = -1 / (dn * dn);
405  }
406 
407  //we are a lower boundary point
408  if(!isInside(x, y - 1, z)) {
409  C += 1 / ((std::abs(dy) - std::abs(cy)) * ds);
410  S = 0.0;
411  } else {
412  C += 1 / (ds * ds);
413  S = -1 / (ds * ds);
414  }
415 
416  F = -1 / (hr[2] * hr[2]);
417  B = -1 / (hr[2] * hr[2]);
418  C += 2 / (hr[2] * hr[2]);
419 
420  // handle boundary condition in z direction
421 /*
422  if(z == 0 || z == nr[2] - 1) {
423 
424  // case where we are on the NEUMAN BC in Z-direction
425  // where we distinguish two cases
426  if(z == 0)
427  F = 0.0;
428  else
429  B = 0.0;
430 
431  //hr[2]*(nr2[2]-1)/2 = radius
432  double d = hr[2] * (nr[2] - 1) / 2;
433  C += 2 / (d * hr[2]);
434 
435  W /= 2.0;
436  E /= 2.0;
437  N /= 2.0;
438  S /= 2.0;
439  C /= 2.0;
440  scaleFactor *= 0.5;
441 
442  }
443 */
444 }
445 
446 void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
447 
448  double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
449  double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
450 
451  // since every vector for elliptic domains has ALWAYS size 2 we
452  // can catch all cases manually
453  double dx = 0.0;
454  std::multimap<int, double>::iterator it = IntersectXDir.find(y);
455  if(cx < 0)
456  ++it;
457  dx = it->second;
458 
459  double dy = 0.0;
460  it = IntersectYDir.find(x);
461  if(cy < 0)
462  ++it;
463  dy = it->second;
464 
465  double dw = hr[0];
466  double de = hr[0];
467  double dn = hr[1];
468  double ds = hr[1];
469  W = 1.0;
470  E = 1.0;
471  N = 1.0;
472  S = 1.0;
473  F = 1.0;
474  B = 1.0;
475  C = 0.0;
476 
477  //TODO: = cx+hr[0] > dx && cx > 0
478  //if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
481  //de = dx-cx;
482  //}
483 
484  //if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
487  //dw = std::abs(dx)-std::abs(cx);
488  //}
489 
490  //if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
493  //dn = dy-cy;
494  //}
495 
496  //if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
499  //ds = std::abs(dy)-std::abs(cy);
500  //}
501 
502  //TODO: = cx+hr[0] > dx && cx > 0
503  //if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
504  //we are a right boundary point
505  if(!isInside(x + 1, y, z)) {
506  de = dx - cx;
507  E = 0.0;
508  }
509 
510  //if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
511  //we are a left boundary point
512  if(!isInside(x - 1, y, z)) {
513  dw = std::abs(dx) - std::abs(cx);
514  W = 0.0;
515  }
516 
517  //if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
518  //we are a upper boundary point
519  if(!isInside(x, y + 1, z)) {
520  dn = dy - cy;
521  N = 0.0;
522  }
523 
524  //if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
525  //we are a lower boundary point
526  if(!isInside(x, y - 1, z)) {
527  ds = std::abs(dy) - std::abs(cy);
528  S = 0.0;
529  }
530 
531  //2/dw*(dw_de)
532  W *= -1.0 / (dw * (dw + de));
533  E *= -1.0 / (de * (dw + de));
534  N *= -1.0 / (dn * (dn + ds));
535  S *= -1.0 / (ds * (dn + ds));
536  F = -1 / (hr[2] * (hr[2] + hr[2]));
537  B = -1 / (hr[2] * (hr[2] + hr[2]));
538 
539  //TODO: problem when de,dw,dn,ds == 0
540  //is NOT a regular BOUND PT
541  C += 1 / de * 1 / dw;
542  C += 1 / dn * 1 / ds;
543  C += 1 / hr[2] * 1 / hr[2];
544  scaleFactor = 0.5;
545 
546 
547  //for regular gridpoints no problem with symmetry, just boundary
548  //z direction is right
549  //implement isLastInside(dir)
550  //we have LOCAL x,y coordinates!
551 
552  /*
553  if(dw != 0 && !wIsB)
554  W = -1/dw * (dn+ds) * 2*hr[2];
555  else
556  W = 0;
557  if(de != 0 && !eIsB)
558  E = -1/de * (dn+ds) * 2*hr[2];
559  else
560  E = 0;
561  if(dn != 0 && !nIsB)
562  N = -1/dn * (dw+de) * 2*hr[2];
563  else
564  N = 0;
565  if(ds != 0 && !sIsB)
566  S = -1/ds * (dw+de) * 2*hr[2];
567  else
568  S = 0;
569  F = -(dw+de)*(dn+ds)/hr[2];
570  B = -(dw+de)*(dn+ds)/hr[2];
571  */
572 
573  //if(dw != 0)
574  //W = -2*hr[2]*(dn+ds)/dw;
575  //else
576  //W = 0;
577  //if(de != 0)
578  //E = -2*hr[2]*(dn+ds)/de;
579  //else
580  //E = 0;
581  //if(dn != 0)
582  //N = -2*hr[2]*(dw+de)/dn;
583  //else
584  //N = 0;
585  //if(ds != 0)
586  //S = -2*hr[2]*(dw+de)/ds;
587  //else
588  //S = 0;
589  //F = -(dw+de)*(dn+ds)/hr[2];
590  //B = -(dw+de)*(dn+ds)/hr[2];
591 
594  //scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr[2]);
595 
596  // catch the case where a point lies on the boundary
597  //FIXME: do this more elegant!
598  //double m1 = dw*de;
599  //double m2 = dn*ds;
600  //if(de == 0)
601  //m1 = dw;
602  //if(dw == 0)
603  //m1 = de;
604  //if(dn == 0)
605  //m2 = ds;
606  //if(ds == 0)
607  //m2 = dn;
610  //C = 2/hr[2];
611  //if(dw != 0 || de != 0)
612  //C *= (dw+de);
613  //if(dn != 0 || ds != 0)
614  //C *= (dn+ds);
615  //if(dw != 0 || de != 0)
616  //C += (2*hr[2])*(dn+ds)*(dw+de)/m1;
617  //if(dn != 0 || ds != 0)
618  //C += (2*hr[2])*(dw+de)*(dn+ds)/m2;
619 
620  //handle Neumann case
621  //if(z == 0 || z == nr[2]-1) {
622 
623  //if(z == 0)
624  //F = 0.0;
625  //else
626  //B = 0.0;
627 
629  //W = W/2.0;
630  //E = E/2.0;
631  //N = N/2.0;
632  //S = S/2.0;
633  //C /= 2.0;
634 
635  //scaleFactor /= 2.0;
636  //}
637 
638  // handle boundary condition in z direction
639  if(z == 0 || z == nr[2] - 1) {
640 
641  // case where we are on the NEUMAN BC in Z-direction
642  // where we distinguish two cases
643  if(z == 0)
644  F = 0.0;
645  else
646  B = 0.0;
647 
648  //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
649  //hr[2]*(nr2[2]-1)/2 = radius
650  double d = hr[2] * (nr[2] - 1) / 2;
651  C += 2 / (d * hr[2]);
652 
653  W /= 2.0;
654  E /= 2.0;
655  N /= 2.0;
656  S /= 2.0;
657  C /= 2.0;
658  scaleFactor /= 2.0;
659 
660  }
661 }
662 
663 
664 #endif //#ifdef HAVE_SAAMG_SOLVER
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
Vektor
Definition: TSVMeta.h:24
nr
const int nr
Definition: ClassicRandom.h:24
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
mslang::iterator
std::string::iterator iterator
Definition: MSLang.h:16
Generated on Sat Jan 4 2020 14:48:05 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5