OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
RectangularDomain.cpp
Go to the documentation of this file.
1 #ifdef HAVE_SAAMG_SOLVER
2 
3 #include "Solvers/RectangularDomain.h"
4 #include "Utilities/OpalException.h"
5 
6 RectangularDomain::RectangularDomain(Vector_t nr, Vector_t hr) {
7  setNr(nr);
8  setHr(hr);
9  a_m = 0.1;
10  b_m = 0.1;
11  nxy_m = nr[0] * nr[1];
12 }
13 
14 RectangularDomain::RectangularDomain(double a, double b, Vector_t nr, Vector_t hr) {
15  a_m = a;
16  b_m = b;
17  setNr(nr);
18  setHr(hr);
19  nxy_m = nr[0] * nr[1];
20 }
21 
22 void RectangularDomain::compute(Vector_t hr){
23  setHr(hr);
24  nxy_m = nr[0] * nr[1];
25 }
26 
27 int RectangularDomain::getNumXY(int z) {
28  return nxy_m;
29 }
30 
31 void RectangularDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
32 
33  //scaleFactor = 1.0;
34 
35  //W = -1/hr[0]*1/hr[0];
36  //E = -1/hr[0]*1/hr[0];
37  //N = -1/hr[1]*1/hr[1];
38  //S = -1/hr[1]*1/hr[1];
39  //F = -1/hr[2]*1/hr[2];
40  //B = -1/hr[2]*1/hr[2];
41  //C = 2/hr[0]*1/hr[0] + 2/hr[1]*1/hr[1] + 2/hr[2]*1/hr[2];
42 
43  scaleFactor = hr[0] * hr[1] * hr[2];
44  W = -hr[1] * hr[2] / hr[0];
45  E = -hr[1] * hr[2] / hr[0];
46  N = -hr[0] * hr[2] / hr[1];
47  S = -hr[0] * hr[2] / hr[1];
48  F = -hr[0] * hr[1] / hr[2];
49  B = -hr[0] * hr[1] / hr[2];
50  C = 2 * hr[1] * hr[2] / hr[0] + 2 * hr[0] * hr[2] / hr[1] + 2 * hr[0] * hr[1] / hr[2];
51 
52  if(!isInside(x + 1, y, z))
53  E = 0.0; //we are a right boundary point
54 
55  if(!isInside(x - 1, y, z))
56  W = 0.0; //we are a left boundary point
57 
58  if(!isInside(x, y + 1, z))
59  N = 0.0; //we are a upper boundary point
60 
61  if(!isInside(x, y - 1, z))
62  S = 0.0; //we are a lower boundary point
63 
64  //dirichlet
65  if(z == 0)
66  F = 0.0;
67  if(z == nr[2] - 1)
68  B = 0.0;
69 
70  if(false /*z == 0 || z == nr[2]-1*/) {
71 
72  //case where we are on the Robin BC in Z-direction
73  //where we distinguish two cases
74  //IFF: this values should not matter because they
75  //never make it into the discretization matrix
76  if(z == 0)
77  F = 0.0;
78  else
79  B = 0.0;
80 
81  //add contribution of Robin discretization to center point
82  //d the distance between the center of the bunch and the boundary
83  //double cx = (x-(nr[0]-1)/2)*hr[0];
84  //double cy = (y-(nr[1]-1)/2)*hr[1];
85  //double cz = hr[2]*(nr[2]-1);
86  //double d = sqrt(cx*cx+cy*cy+cz*cz);
87  double d = hr[2] * (nr[2] - 1) / 2;
88  //cout << cx << " " << cy << " " << cz << " " << 2/(d*hr[2]) << endl;
89  C += 2 / (d * hr[2]);
90  //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
91 
92  //scale all stencil-points in z-plane with 0.5 (Neumann discretization)
93  W /= 2.0;
94  E /= 2.0;
95  N /= 2.0;
96  S /= 2.0;
97  C /= 2.0;
98  scaleFactor *= 0.5;
99  }
100 
101 
102  //simple check if center value of stencil is positive
103 #ifdef DEBUG
104  if(C <= 0)
105  throw OpalException("RectangularDomain::getBoundaryStencil",
106  "Stencil C is <= 0! This case should never occure!");
107 #endif
108 }
109 
110 void RectangularDomain::getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
111 
112  int x = 0, y = 0, z = 0;
113 
114  getCoord(idx, x, y, z);
115  getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
116 
117 }
118 
119 void RectangularDomain::getNeighbours(int idx, double &W, double &E, double &S, double &N, double &F, double &B) {
120 
121  int x = 0, y = 0, z = 0;
122 
123  getCoord(idx, x, y, z);
124  getNeighbours(x, y, z, W, E, S, N, F, B);
125 
126 }
127 
128 void RectangularDomain::getNeighbours(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B) {
129 
130  if(x > 0)
131  W = getIdx(x - 1, y, z);
132  else
133  W = -1;
134  if(x < nr[0] - 1)
135  E = getIdx(x + 1, y, z);
136  else
137  E = -1;
138 
139  if(y < nr[1] - 1)
140  N = getIdx(x, y + 1, z);
141  else
142  N = -1;
143  if(y > 0)
144  S = getIdx(x, y - 1, z);
145  else
146  S = -1;
147 
148  if(z > 0)
149  F = getIdx(x, y, z - 1);
150  else
151  F = -1;
152  if(z < nr[2] - 1)
153  B = getIdx(x, y, z + 1);
154  else
155  B = -1;
156 }
157 
158 /*
159 void MGPoissonSolver::getNeighbours(const int idx, int& W, int& E, int& S, int& N, int& F, int& B, int& numOutOfDomain)
160 {
161 
162  int ixy, iz, iy, ix;
163  int nx = nr_m[0];
164  int ny = nr_m[1];
165  int nz = nr_m[2];
166  int nxy = nx*ny;
167  ixy = idx % nxy;
168 
169  ix = ixy % nx;
170  iy = (ixy - ix) / nx;
171  iz = (idx - ixy) / nxy;
172  numOutOfDomain = 0;
173 
174  if (iz == 0)
175  {F = -1; numOutOfDomain++;}
176  else
177  F = idx - nxy;
178  if (iz == nz - 1)
179  {B = -1; numOutOfDomain++;}
180  else
181  B = idx + nxy;
182 
183  if (ix == 0)
184  {W = -1; numOutOfDomain++;}
185  else
186  W = ixy - 1;
187  if (ix == nx - 1)
188  {E = -1; numOutOfDomain++;}
189  else
190  E = ixy + 1;
191 
192  if (iy == 0)
193  {S = -1; numOutOfDomain++;}
194  else
195  S = ixy - nx;
196  if (iy == ny - 1)
197  {N = -1; numOutOfDomain++;}
198  else
199  N = ixy + nx;
200 
201  int cor = iz * nxy;
202  switch(W) {
203  case -1: break;
204  default: W += cor;
205  }
206  switch(E) {
207  case -1: break;
208  default: E += cor;
209  }
210  switch(S) {
211  case -1: break;
212  default: S += cor;
213  }
214  switch(N) {
215  case -1: break;
216  default: N += cor;
217  }
218 
219 }
220 
221 //at the moment this stencil has neumann everywhere except on the z=0 plane (Dirichlet)
222 //this code is experimental/not-optimized/not-final at the moment
223 inline Epetra_CrsMatrix* MGPoissonSolver::stencil3DOneSidedDirichlet(Vector_t hr)
224 {
225  Epetra_CrsMatrix* Matrix = new Epetra_CrsMatrix(Copy, *Map, 7);
226 
227  Vector_t hr2 = hr*hr;
228  double a = 2*(1/hr2[0]+1/hr2[1]+1/hr2[2]);
229  double b = -1/hr2[0];
230  double d = -1/hr2[1];
231  double f = -1/hr2[2];
232 
233  int NumMyElements = Map->NumMyElements();
234  int* MyGlobalElements = Map->MyGlobalElements();
235 
236  int W, E, S, N, F, B, numout;
237  std::vector<double> Values(6);
238  std::vector<int> Indices(6);
239 
240  for (int i = 0 ; i < NumMyElements ; ++i)
241  {
242  getNeighbours(MyGlobalElements[i], W, E, S, N, F, B, numout);
243 
244  int NumEntries = 0;
245  double diag = a;
246 
247  if(F != -1) {
248 
249  switch(numout) {
250 
251  case 3: {
252  if(W == -1) {
253  Indices[NumEntries] = E;
254  Values[NumEntries++] = 0.25*b;
255  }
256 
257  if(E == -1) {
258  Indices[NumEntries] = W;
259  Values[NumEntries++] = 0.25*b;
260  }
261 
262  if(N == -1) {
263  Indices[NumEntries] = S;
264  Values[NumEntries++] = 0.25*d;
265  }
266 
267  if(S == -1) {
268  Indices[NumEntries] = N;
269  Values[NumEntries++] = 0.25*d;
270  }
271 
272  if(B == -1) {
273  Indices[NumEntries] = F;
274  Values[NumEntries++] = 0.25*f;
275  }
276 
277  if(F == -1) {
278  Indices[NumEntries] = B;
279  Values[NumEntries++] = 0.25*f;
280  }
281 
282  diag *= 0.125;
283 
284  if(NumEntries != 3)
285  cout << "ERROR: while calc corners in stencil" << endl;
286 
287  break;
288 
289  }
290  case 2: {
291 
292  if(W != -1 && E != -1) {
293  Indices[NumEntries] = W;
294  Values[NumEntries++] = 0.25*b;
295  Indices[NumEntries] = E;
296  Values[NumEntries++] = 0.25*b;
297  }
298  if(E == -1 && W != -1) {
299  Indices[NumEntries] = W;
300  Values[NumEntries++] = 0.5*b;
301  }
302  if(W == -1 && E != -1) {
303  Indices[NumEntries] = E;
304  Values[NumEntries++] = 0.5*b;
305  }
306 
307 
308  if(N != -1 && S != -1) {
309  Indices[NumEntries] = N;
310  Values[NumEntries++] = 0.25*d;
311  Indices[NumEntries] = S;
312  Values[NumEntries++] = 0.25*d;
313  }
314  if(S == -1 && N != -1) {
315  Indices[NumEntries] = N;
316  Values[NumEntries++] = 0.5*d;
317  }
318  if(N == -1 && S != -1) {
319  Indices[NumEntries] = S;
320  Values[NumEntries++] = 0.5*d;
321  }
322 
323 
324  if(B != -1 && F != -1) {
325  Indices[NumEntries] = B;
326  Values[NumEntries++] = 0.25*f;
327  Indices[NumEntries] = F;
328  Values[NumEntries++] = 0.25*f;
329  }
330  if(B == -1 && F != -1) {
331  Indices[NumEntries] = F;
332  Values[NumEntries++] = 0.5*f;
333  }
334  if(F == -1 && B != -1) {
335  Indices[NumEntries] = B;
336  Values[NumEntries++] = 0.5*f;
337  }
338 
339  diag *= 0.25;
340 
341  if(NumEntries != 4)
342  cout << "ERROR: while calc edge stencil elements" << endl;
343 
344  break;
345  }
346  case 1: {
347 
348 
349  if(W != -1 && E != -1) {
350  Indices[NumEntries] = W;
351  Values[NumEntries++] = 0.5*b;
352  Indices[NumEntries] = E;
353  Values[NumEntries++] = 0.5*b;
354  }
355  if(E == -1 && W != -1) {
356  Indices[NumEntries] = W;
357  Values[NumEntries++] = b;
358  }
359  if(W == -1 && E != -1) {
360  Indices[NumEntries] = E;
361  Values[NumEntries++] = b;
362  }
363 
364 
365  if(N != -1 && S != -1) {
366  Indices[NumEntries] = S;
367  Values[NumEntries++] = 0.5*d;
368  Indices[NumEntries] = N;
369  Values[NumEntries++] = 0.5*d;
370  }
371  if(S == -1 && N != -1) {
372  Indices[NumEntries] = N;
373  Values[NumEntries++] = d;
374  }
375  if(N == -1 && S != -1) {
376  Indices[NumEntries] = S;
377  Values[NumEntries++] = d;
378  }
379 
380 
381  if(B != -1 && F != -1) {
382  Indices[NumEntries] = F;
383  Values[NumEntries++] = 0.5*f;
384  Indices[NumEntries] = B;
385  Values[NumEntries++] = 0.5*f;
386  }
387  if(B == -1 && F != -1) {
388  Indices[NumEntries] = F;
389  Values[NumEntries++] = f;
390  }
391  if(F == -1 && B != -1) {
392  Indices[NumEntries] = B;
393  Values[NumEntries++] = f;
394  }
395 
396  diag *= 0.5;
397 
398  if(NumEntries != 5)
399  cout << "ERROR: calc surface elements of stencil" << endl;
400 
401  break;
402  }
403 
404  case 0: {
405  if(W != -1) {
406  Indices[NumEntries] = W;
407  Values[NumEntries++] = b;
408  } if(E != -1) {
409  Indices[NumEntries] = E;
410  Values[NumEntries++] = b;
411  } if(S != -1) {
412  Indices[NumEntries] = S;
413  Values[NumEntries++] = d;
414  } if(N != -1) {
415  Indices[NumEntries] = N;
416  Values[NumEntries++] = d;
417  } if(F != -1) {
418  Indices[NumEntries] = F;
419  Values[NumEntries++] = f;
420  } if(B != -1) {
421  Indices[NumEntries] = B;
422  Values[NumEntries++] = f;
423  }
424 
425  break;
426  }
427  default: {
428  cout << "ERROR CREATING THE STENCIL: points out of domain (" << numout << ") is >3 OR <0" << endl;
429  }
430  }
431  } else {
432  //F == -1 --> this is our dirichlet boundary surface
433  switch(numout) {
434  case 3: {
435 
436  if(W == -1) {
437  Indices[NumEntries] = E;
438  Values[NumEntries++] = 0.5*b;
439  }
440 
441  if(E == -1) {
442  Indices[NumEntries] = W;
443  Values[NumEntries++] = 0.5*b;
444  }
445 
446  if(N == -1) {
447  Indices[NumEntries] = S;
448  Values[NumEntries++] = 0.5*d;
449  }
450 
451  if(S == -1) {
452  Indices[NumEntries] = N;
453  Values[NumEntries++] = 0.5*d;
454  }
455 
456  Indices[NumEntries] = B;
457  Values[NumEntries++] = 0.25*f;
458 
459  if(NumEntries != 3)
460  cout << "ERROR: calc corner (below) elements of stencil" << endl;
461 
462  diag *= 0.25;
463  break;
464 
465  }
466  case 2:{
467 
468  if(W != -1 && E != -1) {
469  Indices[NumEntries] = W;
470  Values[NumEntries++] = 0.5*b;
471  Indices[NumEntries] = E;
472  Values[NumEntries++] = 0.5*b;
473  }
474  if(W == -1 && E != -1) {
475  Indices[NumEntries] = E;
476  Values[NumEntries++] = b;
477  }
478 
479  if(E == -1 && W != -1) {
480  Indices[NumEntries] = W;
481  Values[NumEntries++] = b;
482  }
483 
484  if(N != -1 && S != -1) {
485  Indices[NumEntries] = N;
486  Values[NumEntries++] = 0.5*d;
487  Indices[NumEntries] = S;
488  Values[NumEntries++] = 0.5*d;
489  }
490  if(N == -1 && S != -1) {
491  Indices[NumEntries] = S;
492  Values[NumEntries++] = d;
493  }
494 
495  if(S == -1 && N != -1) {
496  Indices[NumEntries] = N;
497  Values[NumEntries++] = d;
498  }
499 
500  Indices[NumEntries] = B;
501  Values[NumEntries++] = 0.25*f;
502 
503  if(NumEntries != 4)
504  cout << "ERROR: calc edge (below) elements of stencil" << endl;
505 
506  diag *= 0.25;
507  break;
508 
509  }
510  case 1:{
511 
512  Indices[NumEntries] = E;
513  Values[NumEntries++] = b;
514 
515  Indices[NumEntries] = W;
516  Values[NumEntries++] = b;
517 
518  Indices[NumEntries] = S;
519  Values[NumEntries++] = d;
520 
521  Indices[NumEntries] = N;
522  Values[NumEntries++] = d;
523 
524  Indices[NumEntries] = B;
525  Values[NumEntries++] = f;
526 
527  if(NumEntries != 5)
528  cout << "ERROR: calc surface (below) elements of stencil: " << NumEntries << endl;
529 
530  break;
531 
532  }
533  default: {
534  cout << "error in dirichlet surface!" << endl;
535  }
536  }
537  }
538  // put the off-diagonal entries
539  Matrix->InsertGlobalValues(MyGlobalElements[i], NumEntries,
540  &Values[0], &Indices[0]);
541 
542  // Put in the diagonal entry
543  Matrix->InsertGlobalValues(MyGlobalElements[i], 1,
544  &diag, MyGlobalElements + i);
545 
546  }
547  Matrix->FillComplete();
548  Matrix->OptimizeStorage();
549 
550  return(Matrix);
551 }
552 
553 //stencil for longitudinal neumann and transversal dirichlet BC
554 //this code is experimental/not-optimized/not-final at the moment
555 inline Epetra_CrsMatrix* MGPoissonSolver::stencil3DLongitudinalNeumann(Vector_t hr)
556 {
557  Epetra_CrsMatrix* Matrix = new Epetra_CrsMatrix(Copy, *Map, 7);
558 
559  Vector_t hr2 = hr*hr;
560  double a = 2*(1/hr2[0]+1/hr2[1]+1/hr2[2]);
561  double b = -1/hr2[0];
562  double d = -1/hr2[1];
563  double f = -1/hr2[2];
564 
565  int NumMyElements = Map->NumMyElements();
566  int* MyGlobalElements = Map->MyGlobalElements();
567 
568  int W, E, S, N, F, B, numout;
569  std::vector<double> Values(6);
570  std::vector<int> Indices(6);
571 
572  for (int i = 0 ; i < NumMyElements ; ++i)
573  {
574  getNeighbours(MyGlobalElements[i], W, E, S, N, F, B, numout);
575 
576  int NumEntries = 0;
577  double diag = a;
578  //numout = 0;
579 
580  //transversal direction: if below == -1 || above == -1 ===> NEUMANN
581  //else dirichlet:
582  if(F != -1 && B != -1)
583  numout = 0;
584 
585  //only on left longitudinal end open bc -- else dirichlet
586  //if(F != -1)
587  // numout = 0;
588 
589  switch(numout) {
590 
591  case 3: {
592  if(W == -1) {
593  Indices[NumEntries] = E;
594  Values[NumEntries++] = 0.5*b;
595  }
596 
597  if(E == -1) {
598  Indices[NumEntries] = W;
599  Values[NumEntries++] = 0.5*b;
600  }
601 
602  if(N == -1) {
603  Indices[NumEntries] = S;
604  Values[NumEntries++] = 0.5*d;
605  }
606 
607  if(S == -1) {
608  Indices[NumEntries] = N;
609  Values[NumEntries++] = 0.5*d;
610  }
611 
612  if(B == -1) {
613  Indices[NumEntries] = F;
614  Values[NumEntries++] = f;
615  }
616 
617  if(F == -1) {
618  Indices[NumEntries] = B;
619  Values[NumEntries++] = f;
620  }
621 
622  diag *= 0.5;
623 
624  if(NumEntries != 3)
625  cout << "ERROR: while calc corners in stencil" << endl;
626 
627  break;
628 
629  }
630  case 2: { //edges
631 
632  if(E != -1) {
633  Indices[NumEntries] = E;
634  Values[NumEntries++] = 0.5*b;
635  }
636  if(W != -1) {
637  Indices[NumEntries] = W;
638  Values[NumEntries++] = 0.5*b;
639  }
640 
641  if(S != -1) {
642  Indices[NumEntries] = S;
643  Values[NumEntries++] = 0.5*d;
644  }
645  if(N != -1) {
646  Indices[NumEntries] = N;
647  Values[NumEntries++] = 0.5*d;
648  }
649 
650  if(B == -1) {
651  Indices[NumEntries] = F;
652  Values[NumEntries++] = f;
653  }
654  if(F == -1) {
655  Indices[NumEntries] = B;
656  Values[NumEntries++] = f;
657  }
658 
659  diag *= 0.5;
660 
661  if(NumEntries != 4)
662  cout << "ERROR: while calc edge stencil elements" << endl;
663 
664  break;
665  }
666  case 1: {
667 
668  if(E != -1) {
669  Indices[NumEntries] = E;
670  Values[NumEntries++] = 0.5*b;
671  }
672  if(W != -1) {
673  Indices[NumEntries] = W;
674  Values[NumEntries++] = 0.5*b;
675  }
676 
677  if(S != -1) {
678  Indices[NumEntries] = S;
679  Values[NumEntries++] = 0.5*d;
680  }
681  if(N != -1) {
682  Indices[NumEntries] = N;
683  Values[NumEntries++] = 0.5*d;
684  }
685 
686  if(B == -1) {
687  Indices[NumEntries] = F;
688  Values[NumEntries++] = f;
689  }
690  if(F == -1) {
691  Indices[NumEntries] = B;
692  Values[NumEntries++] = f;
693  }
694 
695  diag *= 0.5;
696 
697  if(NumEntries != 5)
698  cout << "ERROR: calc surface elements of stencil" << endl;
699 
700  break;
701  }
702 
703  case 0: { //interior points (& dirichlet)
704  if(W != -1) {
705  Indices[NumEntries] = W;
706  Values[NumEntries++] = b;
707  } if(E != -1) {
708  Indices[NumEntries] = E;
709  Values[NumEntries++] = b;
710  } if(S != -1) {
711  Indices[NumEntries] = S;
712  Values[NumEntries++] = d;
713  } if(N != -1) {
714  Indices[NumEntries] = N;
715  Values[NumEntries++] = d;
716  } if(F != -1) {
717  Indices[NumEntries] = F;
718  Values[NumEntries++] = f;
719  } if(B != -1) {
720  Indices[NumEntries] = B;
721  Values[NumEntries++] = f;
722  }
723 
724  break;
725  }
726  default: {
727  cout << "ERROR CREATING THE STENCIL: points out of domain (" << numout << ") is >3 OR <0" << endl;
728  }
729  }
730  // put the off-diagonal entries
731  Matrix->InsertGlobalValues(MyGlobalElements[i], NumEntries,
732  &Values[0], &Indices[0]);
733 
734  // Put in the diagonal entry
735  Matrix->InsertGlobalValues(MyGlobalElements[i], 1,
736  &diag, MyGlobalElements + i);
737 
738  }
739  Matrix->FillComplete();
740  Matrix->OptimizeStorage();
741 
742  return(Matrix);
743 }
744 */
745 
746 #endif //#ifdef HAVE_SAAMG_SOLVER
Vektor
Definition: TSVMeta.h:24
OpalException
The base class for all OPAL exceptions.
Definition: OpalException.h:28
nr
const int nr
Definition: ClassicRandom.h:24
Generated on Sat Jan 4 2020 14:48:06 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5