PlasCom2  1.0
XPACC Multi-physics simluation application
Stencil.C
Go to the documentation of this file.
1 #include "Stencil.H"
2 #include "OperatorKernels.H"
3 #include <iostream>
4 
5 double GetFunctionOrder(std::vector<double> &functionValues,std::vector<double> &independentVariable)
6 {
7  int numVals = independentVariable.size();
8  double f0 = functionValues[0];
9  double i0 = independentVariable[0];
10  double minorder = 10000.0000;
11  for(int iVar = 1;iVar < numVals;iVar++){
12  double f1 = functionValues[iVar];
13  double i1 = independentVariable[iVar];
14  double order = std::log(f0/f1)/std::log(i0/i1);
15 
16  if(order < minorder) minorder = order;
17  // std::cout << "Order[" << iVar << "] = " << order << std::endl;
18  f0 = f1;
19  i0 = i1;
20  }
21  return(minorder);
22 }
23 
24 namespace plascom2 {
25  namespace operators {
26 
27  namespace dissipation {
28 
29  int Initialize(stencilset &stencilSet1,stencilset &stencilSet2,int interiorOrder)
30  {
31 
32  if(interiorOrder == 2){
33 
34  stencilSet1.ownData = true;
35  stencilSet1.numStencils = 4;
36  stencilSet1.numValues = 7;
37  stencilSet1.boundaryDepth = 1;
38  stencilSet1.boundaryWidth = 2;
39 
40  stencilSet1.stencilSizes = new int [stencilSet1.numStencils];
41  stencilSet1.stencilStarts = new int [stencilSet1.numStencils];
42  stencilSet1.stencilWeights = new double [stencilSet1.numValues];
43  stencilSet1.stencilOffsets = new int [stencilSet1.numValues];
44 
45  stencilSet1.stencilSizes[0] = 2;
46  stencilSet1.stencilSizes[1] = 2;
47  stencilSet1.stencilSizes[2] = 2;
48  stencilSet1.stencilSizes[3] = 1;
49 
50  stencilSet1.stencilStarts[0] = 0;
51  for(int iStencil = 1;iStencil < stencilSet1.numStencils;iStencil++)
52  stencilSet1.stencilStarts[iStencil] =
53  stencilSet1.stencilStarts[iStencil-1]+ stencilSet1.stencilSizes[iStencil-1];
54 
55  // central
56  stencilSet1.stencilOffsets[0] = 0; // weight 1
57  stencilSet1.stencilOffsets[1] = -1; // weight -1
58  // left boundary
59  stencilSet1.stencilOffsets[2] = 0; // weight -1
60  stencilSet1.stencilOffsets[3] = 1; // weight 1
61  // right boundary
62  stencilSet1.stencilOffsets[4] = 0; // weight 1
63  stencilSet1.stencilOffsets[5] = -1; // weight -1
64  // blanker
65  stencilSet1.stencilOffsets[6] = 0;
66 
67  stencilSet1.stencilWeights[0] = 1.0;
68  stencilSet1.stencilWeights[1] = -1.0;
69  stencilSet1.stencilWeights[2] = -1.0;
70  stencilSet1.stencilWeights[3] = 1.0;
71  stencilSet1.stencilWeights[4] = 1.0;
72  stencilSet1.stencilWeights[5] = -1.0;
73  stencilSet1.stencilWeights[6] = 0.0;
74 
75  stencilSet1.boundaryWeight = 0.0;
76 
77  // STENCIL SET 2
78 
79  stencilSet2.ownData = true;
80  stencilSet2.numStencils = 6;
81  stencilSet2.numValues = 11;
82  stencilSet2.boundaryDepth = 2;
83  stencilSet2.boundaryWidth = 3;
84 
85  stencilSet2.stencilSizes = new int [stencilSet2.numStencils];
86  stencilSet2.stencilStarts = new int [stencilSet2.numStencils];
87  stencilSet2.stencilWeights = new double [stencilSet2.numValues];
88  stencilSet2.stencilOffsets = new int [stencilSet2.numValues];
89 
90  stencilSet2.stencilSizes[0] = 2;
91  stencilSet2.stencilSizes[1] = 2;
92  stencilSet2.stencilSizes[2] = 3;
93  stencilSet2.stencilSizes[3] = 1;
94  stencilSet2.stencilSizes[4] = 2;
95  stencilSet2.stencilSizes[5] = 1;
96 
97  stencilSet2.stencilStarts[0] = 0;
98  for(int iStencil = 1;iStencil < stencilSet2.numStencils;iStencil++)
99  stencilSet2.stencilStarts[iStencil] =
100  stencilSet2.stencilStarts[iStencil-1]+ stencilSet2.stencilSizes[iStencil-1];
101 
102  // Central
103  stencilSet2.stencilOffsets[0] = 0; // weight -1
104  stencilSet2.stencilOffsets[1] = 1; // weight 1
105  // On left boundary
106  stencilSet2.stencilOffsets[2] = 0; // weight 2
107  stencilSet2.stencilOffsets[3] = 1; // weight 2
108  // Left boundary + 1
109  stencilSet2.stencilOffsets[4] = -1; // weight -1
110  stencilSet2.stencilOffsets[5] = 0; // weight -1
111  stencilSet2.stencilOffsets[6] = 1; // weight 1
112  // Right Boundary
113  stencilSet2.stencilOffsets[7] = 0; // weight -2
114  // Right Boundary - 1
115  stencilSet2.stencilOffsets[8] = 0; // weight -1
116  stencilSet2.stencilOffsets[9] = 1; // weight 1
117  // Blanker
118  stencilSet2.stencilOffsets[10] = 0; // weight 0
119 
120  // Central
121  stencilSet2.stencilWeights[0] = -1.0;
122  stencilSet2.stencilWeights[1] = 1.0;
123  // On left boundary
124  stencilSet2.stencilWeights[2] = 2.0;
125  stencilSet2.stencilWeights[3] = 2.0;
126  // Left boundary + 1
127  stencilSet2.stencilWeights[4] = -1.0;
128  stencilSet2.stencilWeights[5] = -1.0;
129  stencilSet2.stencilWeights[6] = 1.0;
130  // Right boundary
131  stencilSet2.stencilWeights[7] = -2.0;
132  // Right boundary - 1
133  stencilSet2.stencilWeights[8] = -1.0;
134  stencilSet2.stencilWeights[9] = 1.0;
135  // Blanker
136  stencilSet2.stencilWeights[10] = 0;
137 
138  stencilSet2.boundaryWeight = 1.0;
139 
140  } else if (interiorOrder == 4) {
141 
142  stencilSet1.ownData = true;
143  stencilSet1.numStencils = 4;
144  stencilSet1.numValues = 10;
145  stencilSet1.boundaryDepth = 1;
146  stencilSet1.boundaryWidth = 3;
147 
148  stencilSet1.stencilSizes = new int [stencilSet1.numStencils];
149  stencilSet1.stencilStarts = new int [stencilSet1.numStencils];
150  stencilSet1.stencilWeights = new double [stencilSet1.numValues];
151  stencilSet1.stencilOffsets = new int [stencilSet1.numValues];
152 
153  stencilSet1.stencilSizes[0] = 3;
154  stencilSet1.stencilSizes[1] = 3;
155  stencilSet1.stencilSizes[2] = 3;
156  stencilSet1.stencilSizes[3] = 1;
157 
158  stencilSet1.stencilStarts[0] = 0;
159  for(int iStencil = 1;iStencil < stencilSet1.numStencils;iStencil++)
160  stencilSet1.stencilStarts[iStencil] =
161  stencilSet1.stencilStarts[iStencil-1]+ stencilSet1.stencilSizes[iStencil-1];
162 
163  // Central Part
164  stencilSet1.stencilOffsets[0] = -1; // weight 1
165  stencilSet1.stencilOffsets[1] = 0; // weight -2
166  stencilSet1.stencilOffsets[2] = 1; // weight 1
167 
168  // First left boundary stencil (on-boundary)
169  stencilSet1.stencilOffsets[3] = 0; // 1
170  stencilSet1.stencilOffsets[4] = 1; // -2
171  stencilSet1.stencilOffsets[5] = 2; // 1
172 
173  // First right boundary stencil (on boundary)
174  stencilSet1.stencilOffsets[6] = -2; // 1
175  stencilSet1.stencilOffsets[7] = -1; // -2
176  stencilSet1.stencilOffsets[8] = 0; // 1
177 
178  // Blanker
179  stencilSet1.stencilOffsets[9] = 0; // 0
180 
181  // Central
182  stencilSet1.stencilWeights[0] = 1.0;
183  stencilSet1.stencilWeights[1] = -2.0;
184  stencilSet1.stencilWeights[2] = 1.0;
185 
186  // Left boundary
187  stencilSet1.stencilWeights[3] = 1.0;
188  stencilSet1.stencilWeights[4] = -2.0;
189  stencilSet1.stencilWeights[5] = 1.0;
190 
191  // Right boundary
192  stencilSet1.stencilWeights[6] = 1.0;
193  stencilSet1.stencilWeights[7] = -2.0;
194  stencilSet1.stencilWeights[8] = 1.0;
195 
196  stencilSet1.stencilWeights[9] = 0.0;
197 
198  stencilSet1.boundaryWeight = 0.0;
199 
200  // STENCIL SET 2
201 
202  stencilSet2.ownData = true;
203  stencilSet2.numStencils = 10;
204  stencilSet2.numValues = 28;
205  stencilSet2.boundaryDepth = 4;
206  stencilSet2.boundaryWidth = 5;
207 
208  stencilSet2.stencilSizes = new int [stencilSet2.numStencils];
209  stencilSet2.stencilStarts = new int [stencilSet2.numStencils];
210  stencilSet2.stencilWeights = new double [stencilSet2.numValues];
211  stencilSet2.stencilOffsets = new int [stencilSet2.numValues];
212 
213  stencilSet2.stencilSizes[0] = 3;
214 
215  stencilSet2.stencilSizes[1] = 2;
216  stencilSet2.stencilSizes[2] = 3;
217  stencilSet2.stencilSizes[3] = 4;
218  stencilSet2.stencilSizes[4] = 3;
219 
220  stencilSet2.stencilSizes[5] = 2;
221  stencilSet2.stencilSizes[6] = 3;
222  stencilSet2.stencilSizes[7] = 4;
223  stencilSet2.stencilSizes[8] = 3;
224 
225  stencilSet2.stencilSizes[9] = 1;
226 
227 
228 
229  stencilSet2.stencilStarts[0] = 0;
230  for(int iStencil = 1;iStencil < stencilSet2.numStencils;iStencil++)
231  stencilSet2.stencilStarts[iStencil] =
232  stencilSet2.stencilStarts[iStencil-1]+ stencilSet2.stencilSizes[iStencil-1];
233 
234  // Central part
235  stencilSet2.stencilOffsets[0] = -1; // weight -1
236  stencilSet2.stencilOffsets[1] = 0; // weight 2
237  stencilSet2.stencilOffsets[2] = 1; // weight -1
238 
239  // Left on-boundary
240  stencilSet2.stencilOffsets[3] = 0; // weight -48/17
241  stencilSet2.stencilOffsets[4] = 1; // weight -48/17
242 
243  // Left boundary + 1
244  stencilSet2.stencilOffsets[5] = -1; // 96/59
245  stencilSet2.stencilOffsets[6] = 0; // 96/59
246  stencilSet2.stencilOffsets[7] = 1; // -48/59
247 
248  // Left boundary + 2
249  stencilSet2.stencilOffsets[8] = -2; // -48/43
250  stencilSet2.stencilOffsets[9] = -1; // -48/43
251  stencilSet2.stencilOffsets[10] = 0; // 96/43
252  stencilSet2.stencilOffsets[11] = 1; // -48/43
253 
254  // Left boundary + 3
255  stencilSet2.stencilOffsets[12] = -1; // -48/59
256  stencilSet2.stencilOffsets[13] = 0; // 96/59
257  stencilSet2.stencilOffsets[14] = 1; // -48/59
258 
259 
260 
261  // Right boundary
262  stencilSet2.stencilOffsets[15] = -1; // -48/17
263  stencilSet2.stencilOffsets[16] = 0; // -48/17
264 
265  // Right boundary - 1
266  stencilSet2.stencilOffsets[17] = -1; // -48/59
267  stencilSet2.stencilOffsets[18] = 0; // 96/59
268  stencilSet2.stencilOffsets[19] = 1; // 96/59
269 
270  // Right boundary - 2
271  stencilSet2.stencilOffsets[20] = -1; // -48/43
272  stencilSet2.stencilOffsets[21] = 0; // 96/43
273  stencilSet2.stencilOffsets[22] = 1; // -48/43
274  stencilSet2.stencilOffsets[23] = 2; // -48/43
275 
276  // Right boundary - 3
277  stencilSet2.stencilOffsets[24] = -1; // -48/49
278  stencilSet2.stencilOffsets[25] = 0; // 96/49
279  stencilSet2.stencilOffsets[26] = 1; // -48/49
280 
281  // Blanker
282  stencilSet2.stencilOffsets[27] = 0; // 0.0
283 
284 
285 
286 
287  // Central
288  stencilSet2.stencilWeights[0] = -1.0;
289  stencilSet2.stencilWeights[1] = 2.0;
290  stencilSet2.stencilWeights[2] = -1.0;
291 
292  // Left boundary
293  stencilSet2.stencilWeights[3] = -48.0/17.0;
294  stencilSet2.stencilWeights[4] = -48.0/17.0;
295 
296  // Left boundary + 1
297  stencilSet2.stencilWeights[5] = 96.0/59.0;
298  stencilSet2.stencilWeights[6] = 96.0/59.0;
299  stencilSet2.stencilWeights[7] = -48.0/59.0;
300 
301  // Left boundary + 2
302  stencilSet2.stencilWeights[8] = -48.0/43.0;
303  stencilSet2.stencilWeights[9] = -48.0/43.0;
304  stencilSet2.stencilWeights[10] = 96.0/43.0;
305  stencilSet2.stencilWeights[11] = -48.0/43.0;
306 
307  // Left boundary + 3
308  stencilSet2.stencilWeights[12] = -48.0/49.0;
309  stencilSet2.stencilWeights[13] = 96.0/49.0;
310  stencilSet2.stencilWeights[14] = -48.0/49.0;
311 
312  // Right boundary
313  stencilSet2.stencilWeights[15] = -48.0/17.0;
314  stencilSet2.stencilWeights[16] = -48.0/17.0;
315 
316  // Right boundary - 1
317  stencilSet2.stencilWeights[17] = -48.0/59.0;
318  stencilSet2.stencilWeights[18] = 96.0/59.0;
319  stencilSet2.stencilWeights[19] = 96.0/59.0;
320 
321  // Right boundary - 2
322  stencilSet2.stencilWeights[20] = -48.0/43.0;
323  stencilSet2.stencilWeights[21] = 96.0/43.0;
324  stencilSet2.stencilWeights[22] = -48.0/43.0;
325  stencilSet2.stencilWeights[23] = -48.0/43.0;
326 
327  // Right boundary - 3
328  stencilSet2.stencilWeights[24] = -48.0/49.0;
329  stencilSet2.stencilWeights[25] = 96.0/49.0;
330  stencilSet2.stencilWeights[26] = -48.0/49.0;
331 
332  // Blanker
333  stencilSet2.stencilWeights[27] = 0.0;
334 
335  stencilSet2.boundaryWeight = 1.0;
336 
337  } else if (interiorOrder == 6) {
338  return(1);
339  } else {
340  return(1);
341  }
342  return(0);
343  }
344  }
345 
346  namespace sbp {
347 
348  // Function to return whether a given
349  // overall spatial order has an SBP
350  // operator. New operators should have
351  // an entry here.
352  bool IsValidOrder(int overallOrder){
353  if(overallOrder == 2 ||
354  overallOrder == 3 ||
355  overallOrder == 4)
356  return(true);
357  return(false);
358  };
359 
360  int Initialize(base &stencilSet,int interiorOrder)
361  {
362  stencilSet.Destroy();
363 
364  int boundaryDepth;
365  int boundaryWidth;
366  double boundaryWeight;
367  std::vector<double> centralWeightsRight;
368  std::vector<std::vector<double> > leftWeights;
369  std::vector<int> leftOrders;
370 
371  if(interiorOrder == 2){
372 
373  // scalars
374  boundaryDepth = 1;
375  boundaryWidth = 2;
376  boundaryWeight = 2.0;
377 
378  // central stencil weights
379  centralWeightsRight.resize(1,.5);
380 
381  // left weights
382  std::vector<double> leftwt(2,1.0);
383  leftwt[0] = -1.0;
384  leftWeights.push_back(leftwt);
385 
386  // left boundary stencil orders
387  leftOrders.resize(1,1);
388 
389  } else if(interiorOrder == 4) {
390 
391  // scalars
392  boundaryDepth = 4;
393  boundaryWidth = 6;
394  boundaryWeight = 48.0/17.0;
395 
396  // central stencil weights
397  centralWeightsRight.resize(2,-1.0/12.0);
398  centralWeightsRight[0] *= -8.0;
399 
400  // left weights
401  leftWeights.resize(boundaryDepth);
402  double leftWeightsArr[4][6] =
403  {
404  {-24.0/17.0, 59.0/34.0 , -4.0/17.0 , -3.0/34.0, 0.0 , 0.0 },
405  {-.5 , 0.0 , .5 , 0.0 , 0.0 , 0.0 },
406  {4.0/43.0 , -59.0/86.0, 0.0 , 59.0/86.0, -4.0/43.0, 0.0 },
407  {3.0/98.0 , 0.0 , -59.0/98.0, 0.0 , 32.0/49.0, -4.0/49.0}
408  };
409 
410  leftWeights.resize(boundaryDepth);
411  for (int i=0; i < boundaryDepth; i++) {
412  leftWeights[i].resize(boundaryWidth);
413  for(int j = 0;j < boundaryWidth;j++)
414  leftWeights[i][j] = leftWeightsArr[i][j];
415  }
416 
417  // left boundary stencil orders
418  leftOrders.resize(4,2);
419 
420  } else if(interiorOrder == 6) {
421 
422  // scalars
423  boundaryDepth = 6;
424  boundaryWidth = 9;
425  boundaryWeight = 43200.0/13649.0;
426 
427  // central stencil weights
428  centralWeightsRight.resize(3);
429  centralWeightsRight[0] = 3.0/4.0;
430  centralWeightsRight[1] = -3.0/20.0;
431  centralWeightsRight[2] = 1.0/60.0;
432 
433  // left weights
434  double leftWeightsArr[6][9] =
435  {
436  {-21600.0/13649.0, 104009.0/54596.0, 30443.0/81894.0,
437  -33311.0/27298.0, 16863.0/27298.0, -15025.0/163788.0,
438  0.0, 0.0, 0.0},
439 
440  {-104009.0/240260.0, 0.0, -311.0/72078.0,
441  20229.0/24026.0, -24337.0/48052.0, 36661.0/360390.0,
442  0.0, 0.0, 0.0},
443 
444  {-30443.0/162660.0, 311.0/32532.0, 0.0,
445  -11155.0/16266.0, 41287.0/32532.0, -21999.0/54220.0,
446  0.0, 0.0, 0.0},
447 
448  {33311.0/107180.0, -20229.0/21436.0, 485.0/1398.0,
449  0.0, 4147.0/21436.0, 25427.0/321540.0,
450  72.0/5359.0, 0.0, 0.0},
451 
452  {-16863.0/78770.0, 24337.0/31508.0, -41287.0/47262.0,
453  -4147.0/15754.0, 0.0, 342523.0/472620.0,
454  -1296.0/7877.0, 144.0/7877.0, 0.0},
455 
456  {15025.0/525612.0, -36661.0/262806.0, 21999.0/87602.0,
457  -25427.0/262806.0, -342523.0/525612.0, 0.0,
458  32400.0/43801.0, -6480.0/43801.0, 720.0/43801.0}
459  };
460 
461  leftWeights.resize(boundaryDepth);
462  for (int i=0; i<boundaryDepth; i++) {
463  leftWeights[i].resize(boundaryWidth);
464  for(int j = 0;j < boundaryWidth;j++)
465  leftWeights[i][j] = leftWeightsArr[i][j];
466  }
467 
468  // left boundary stencil orders
469  leftOrders.resize(6,3);
470  leftOrders[5] = 4;
471 
472  } else if(interiorOrder == 8) {
473  return(1);
474  } else {
475  std::cout << "Invalid interiorOrder " << interiorOrder
476  << " -- no such SBP operator!" << std::endl;
477  return(1);
478  }
479 
480  OperatorSetup(stencilSet, interiorOrder, boundaryDepth, boundaryWidth,
481  boundaryWeight, centralWeightsRight, leftWeights, leftOrders);
482  return(0);
483  }
484 
485  void OperatorSetup(base &stencilSet, int interiorOrder, int boundaryDepth,
486  int boundaryWidth, double boundaryWeight, std::vector<double> centralWeightsRight,
487  std::vector<std::vector<double> > leftWeights, std::vector<int> leftOrders) {
488  // finishes setup of full operator given minimum set of coefficients
489  //
490  // Assumptions / choices:
491  // - central stencil is anti-symmetric, and total length =
492  // interiorOrder+1. Only right half coeffs are passed in.
493  // - right boundary coeffs are anti-symmetric to left boundary.
494  // - boundaryDepth = number of boundary stencils.
495  // - boundaryWeight = first element of invPdiag_block1 (PlasComCM)
496 
497  //std::cout << "--- OperatorSetup: Echo input ---" << std::endl;
498  //std::cout << "interiorOrder " << interiorOrder << std::endl;
499  //std::cout << "boundaryDepth " << boundaryDepth << std::endl;
500  //std::cout << "boundaryWidth " << boundaryWidth << std::endl;
501  //std::cout << "boundaryWeight " << boundaryWeight << std::endl;
502  //std::cout << "centralWeightsRight";
503  //for (int i=0; i<centralWeightsRight.size(); i++) {
504  // std::cout << " " << centralWeightsRight[i];
505  //}
506  //std::cout << std::endl;
507  //std::cout << "leftWeights:";
508  //for (int i=0; i<leftWeights.size(); i++) {
509  // std::cout << std::endl;
510  // for (int j=0; j<leftWeights[i].size(); j++) {
511  // std::cout << " " << leftWeights[i][j];
512  // }
513  //}
514  //std::cout << std::endl;
515  //std::cout << "leftOrders";
516  //for (int i=0; i<leftOrders.size(); i++) {
517  // std::cout << " " << leftOrders[i];
518  //}
519  //std::cout << std::endl;
520  //std::cout << "--- end input ---" << std::endl;
521 
522  stencilSet.ownData = true;
523  stencilSet.numStencils = 2 + 2*boundaryDepth;
524  stencilSet.boundaryDepth = boundaryDepth;
525  stencilSet.boundaryWidth = boundaryWidth;
526  stencilSet.overallOrder = (interiorOrder/2) + 1;
527  stencilSet.boundaryWeight = boundaryWeight;
528 
529  stencilSet.stencilSizes = new int [stencilSet.numStencils];
530  stencilSet.stencilStarts = new int [stencilSet.numStencils];
531  stencilSet.stencilOrders = new int [stencilSet.numStencils];
532 
533  // --- Sizes / Orders ---
534  // Center stencil
535  stencilSet.stencilSizes[0] = interiorOrder;
536  stencilSet.stencilOrders[0] = interiorOrder;
537 
538  // Left boundary stencils
539  std::vector<double> stencilWeightsLeft;
540  std::vector<int> stencilOffsetsLeft;
541 
542  for (int iStencil=0; iStencil<stencilSet.boundaryDepth; iStencil++) {
543  int size = 0;
544 
545  for (int i=0; i<stencilSet.boundaryWidth; i++) {
546  if (std::abs(leftWeights[iStencil][i]) > 0.0) {
547  stencilOffsetsLeft.push_back(i-iStencil);
548  stencilWeightsLeft.push_back(leftWeights[iStencil][i]);
549  size++;
550  }
551  }
552 
553  stencilSet.stencilSizes[iStencil+1] = size;
554  }
555 
556  for (int i=0; i<stencilSet.boundaryDepth; i++) {
557  stencilSet.stencilOrders[i+1] = leftOrders[i];
558  }
559 
560  // Right boundary stencils
561  for (int i=1; i<=stencilSet.boundaryDepth; i++) {
562  stencilSet.stencilSizes[i+stencilSet.boundaryDepth] = stencilSet.stencilSizes[i];
563  stencilSet.stencilOrders[i+stencilSet.boundaryDepth] = stencilSet.stencilOrders[i];
564  }
565 
566  // Hole stencil
567  stencilSet.stencilSizes[2*stencilSet.boundaryDepth+1] = 1;
568  stencilSet.stencilOrders[2*stencilSet.boundaryDepth+1] = 0;
569 
570  // --- numValues / Starts ---
571  stencilSet.numValues = stencilSet.stencilSizes[0];
572  stencilSet.stencilStarts[0] = 0;
573  for(int iStencil = 1;iStencil < stencilSet.numStencils;iStencil++){
574  stencilSet.stencilStarts[iStencil] =
575  stencilSet.stencilStarts[iStencil-1] + stencilSet.stencilSizes[iStencil-1];
576  stencilSet.numValues += stencilSet.stencilSizes[iStencil];
577  }
578 
579  // --- Offsets / Weights ---
580  stencilSet.stencilOffsets = new int [stencilSet.numValues];
581  stencilSet.stencilWeights = new double [stencilSet.numValues];
582 
583  // Central stencil
584  for (int i=1; i<=interiorOrder/2; i++) {
585  int il = interiorOrder/2 - i;
586  int ir = i-1 + interiorOrder/2;
587 
588  // right side
589  stencilSet.stencilOffsets[ir] = i;
590  stencilSet.stencilWeights[ir] = centralWeightsRight[i-1];
591 
592  // left side
593  stencilSet.stencilOffsets[il] = -i;
594  stencilSet.stencilWeights[il] = -centralWeightsRight[i-1];
595  }
596 
597  // Left boundary stencils
598  for (int i=0; i<stencilOffsetsLeft.size(); i++) {
599  stencilSet.stencilOffsets[i+interiorOrder] = stencilOffsetsLeft[i];
600  stencilSet.stencilWeights[i+interiorOrder] = stencilWeightsLeft[i];
601  }
602 
603  // Right boundary stencils
604  for (int il=1; il<=stencilSet.boundaryDepth; il++) {
605  int ir = il + stencilSet.boundaryDepth;
606 
607  int jls = stencilSet.stencilStarts[il];
608  int jle = jls + stencilSet.stencilSizes[il];
609 
610  int jrs = stencilSet.stencilStarts[ir];
611 
612  for (int jl=jls; jl<jle; jl++) {
613  int jr = jrs + (jle-1 - jl);
614  stencilSet.stencilOffsets[jr] = -stencilSet.stencilOffsets[jl];
615  stencilSet.stencilWeights[jr] = -stencilSet.stencilWeights[jl];
616  }
617  }
618 
619  // Hole
620  int iHole = stencilSet.stencilStarts[stencilSet.numStencils-1];
621  stencilSet.stencilOffsets[iHole] = 0;
622  stencilSet.stencilWeights[iHole] = 0.0;
623  }
624 
625  int BruteTest1(base &inOperator, int interiorOrder, int numDim,
626  int numComponents, int numTrials, std::vector<bool> &testResults)
627  {
628  // Tests order of accuracy on each stencil
629  // Also tests correctness of ApplyOperator
630 
631  // we want to ignore the hole stencil
632  int numStencils = inOperator.numStencils - 1;
633  int boundaryDepth = inOperator.boundaryDepth;
634 
635  // 2*numDim*numStencils:
636  // for each stencil, for each dimension, we check both accuracy and order
637  // accuracy test requires initialization to true
638  // order test will overwrite the initial value so we don't care
639  testResults.resize(2*numDim*numStencils,true);
640 
641  size_t *numX = new size_t [numDim]; // number of points in each dimension
642  double *deltaX = new double [numDim]; // spacing in each dimension
643  size_t *opInterval = new size_t [2*numDim]; // start/end indices in each dimension
644 
645  // arrays for data / stencil connectivity
646  // will be sized later based on the size of the domain for a given trial
647  double *inputData = NULL;
648  double *opData = NULL;
649  double *analyticData = NULL;
650  int *stencilID = NULL;
651 
652  // minimum value of x
653  double xmin = .5;
654  // basic length of domain
655  double baseLen = 1.0;
656  // number of points to use in domain -- decrease spacing
657  // instead of increasing number of points
658  // Using boundaryWidth instead of boundaryDepth because
659  // boundaryWidth is the max width of a boundary stencil.
660  // For several of the SBP operators, boundaryWidth > boundaryDepth
661  int baseNum = 2*inOperator.boundaryWidth+1;
662 
663  // tolerance for order of accuracy
664  double epsOrder = 0.1;
665  double epsError = 1e-14;
666 
667  for (int iStencil=0; iStencil<numStencils; iStencil++) {
668 
669  // number of points that need to be shifted to the right each
670  // trial so that the first point with stencilIndex == iStencil
671  // has the same x-value (assuming x == opDir)
672  int shift = -52;
673  if (iStencil == 0) shift = boundaryDepth; // central
674  else if (iStencil <= boundaryDepth) shift = iStencil-1; // left boundary
675  else shift = baseNum - (iStencil-boundaryDepth-1); // right boundary
676 
677  // add 1 because we multiply through by deltaX everywhere
678  int expOrder = inOperator.stencilOrders[iStencil] + 1;
679 
680  for (int iOrder=0; iOrder<=expOrder; iOrder++) {
681  for (int opDir=0; opDir<numDim; opDir++) {
682 
683  // Vector of deltaX for use in computing order
684  std::vector<double> trialSpacing(numTrials,0.0);
685  // Errors for each trial
686  std::vector<double> trialError(numTrials,0.0);
687 
688  // length of domain
689  double len = baseLen;
690  double x0 = xmin;
691 
692  for(int iTrial = 0; iTrial < numTrials; iTrial++){
693 
694  // halve length in opDir each trial
695  // baseLen stays the same
696  len *= 0.5;
697 
698  // Initialize number of grid points and spacing
699  size_t numPoints = 1;
700  for(int iDim = 0; iDim < numDim; iDim++){
701  numX[iDim] = baseNum;
702 
703  if(iDim == opDir) {
704  deltaX[iDim] = len/(numX[iDim]-1);
705  trialSpacing[iTrial] = deltaX[iDim];
706  } else {
707  deltaX[iDim] = baseLen/(numX[iDim]-1);
708  }
709 
710  numPoints *= numX[iDim];
711 
712  // set indices of interval defining domain
713  opInterval[2*iDim] = 1;
714  opInterval[2*iDim+1] = numX[iDim];
715  }
716 
717  // slide minimum x value to right so location of current
718  // stencil stays the same (specifically, 0-point of current stencil)
719  if (iTrial > 0) x0 += shift*deltaX[opDir];
720 
721  inputData = new double [numPoints]; // function values
722  opData = new double [numPoints]; // computed derivatives
723  stencilID = new int [numDim*numPoints]; // stencil connectivity
724  analyticData = new double [numPoints]; // analytical derivatives
725 
726  if(CreateStencilConnectivity(numDim,numX,opInterval,boundaryDepth,stencilID,true))
727  return(1);
728 
729  // fill arrays with function values and analytical derivatives
730  // dirIndex array is for looping over domain
731  // numIndex is set as 1 unless there exists numX for that
732  // dimension, in which case numIndex = numX
733  // basically, the intent is that for unused dimensions, the
734  // loop will only go one time with index 0, and the iPoint
735  // computation will be unaffected by the extra dimensions,
736  // but we don't have to have a 3-way if statement
737  size_t dirIndex[3] = {0,0,0};
738  size_t numIndex[3] = {1,1,1};
739  for (int iDim=0; iDim<numDim; iDim++) numIndex[iDim] = numX[iDim];
740 
741  for (dirIndex[2]=0; dirIndex[2]<numIndex[2]; dirIndex[2]++) {
742  // using numIndex, not numX, because if numDim < 2 it is still
743  // guaranteed to exist and will not mess up the computation
744  // since dirIndex[2] == 0 in that case
745  size_t zIndex = dirIndex[2]*numIndex[0]*numIndex[1];
746 
747  for (dirIndex[1]=0; dirIndex[1]<numIndex[1]; dirIndex[1]++) {
748  // using numIndex instead of numX for consistency with
749  // above, although numX[0] always exists
750  size_t yzIndex = zIndex + dirIndex[1]*numIndex[0];
751 
752  for (dirIndex[0]=0; dirIndex[0]<numIndex[0]; dirIndex[0]++) {
753  size_t iPoint = yzIndex + dirIndex[0];
754 
755  double xCurr = x0 + dirIndex[opDir]*deltaX[opDir];
756 
757  // Test function: f(x) = x^iOrder
758  inputData[iPoint] = std::pow(xCurr,static_cast<double>(iOrder));
759  //inputData[iPoint] = std::sin(xCurr);
760  // APPLYOPERATOR does not divide by deltaX, so multiply
761  // the whole equation through by deltaX
762  // Analytical derivative: f'(x) = iOrder * x^(iOrder-1)
763  analyticData[iPoint] = deltaX[opDir]*static_cast<double>(iOrder)*std::pow(xCurr,static_cast<double>(iOrder-1));
764  //analyticData[iPoint] = deltaX[opDir]*std::cos(xCurr);
765  }
766  }
767  }
768 
769  int fortOpDir = opDir+1;
770  FC_MODULE(operators,applyoperator,OPERATORS,APPLYOPERATOR)
771  (&numDim,numX,&numComponents,&numPoints,
772  &fortOpDir,opInterval,&numStencils,inOperator.stencilSizes,
773  inOperator.stencilStarts, &(inOperator.numValues),
774  inOperator.stencilWeights,inOperator.stencilOffsets,
775  &stencilID[opDir*numPoints],inputData,opData);
776 
777  // flag for whether to store error for order of accuracy
778  // calculations at a given point
779  bool storeError = true;
780 
781  for(size_t iPoint=0; iPoint<numPoints; iPoint++) {
782 
783  int stencilIndex = stencilID[iPoint + opDir*numPoints] - 1;
784  if (stencilIndex != iStencil) continue;
785 
786  double error = std::abs(opData[iPoint] - analyticData[iPoint]);
787 
788  // we store the error one time for each type of stencil
789  //
790  // note that when numDim > 1 there are many points in the
791  // domain with the same stencilIndex and same error (because
792  // the test function only varies in opDir)
793  //
794  // only check this if the test polynomial is the correct order
795  if (expOrder == iOrder && storeError) {
796  trialError[iTrial] = error;
797  storeError = false;
798  }
799 
800  // checking accuracy only for low orders -- error should be
801  // negligible since test function is low order polynomial
802  // this if statement is correct because expOrder has "+1"
803  // baked in -- i.e. if derivative operator is pth order, all
804  // polynomials up to x^p should be differentiated exactly
805  if (iOrder < expOrder) {
806  // index into testResults:
807  // stencil varies most slowly, then opDir, then accuracy/order
808  // add 0 for accuracy, 1 for order
809  testResults[stencilIndex*numDim*2 + opDir*2] = testResults[stencilIndex*numDim*2 + opDir*2] && (error < epsError);
810  }
811  }
812 
813  delete [] inputData;
814  delete [] opData;
815  delete [] stencilID;
816  delete [] analyticData;
817  } // Trial loop
818 
819  // Check order of accuracy
820  // only check this if the test polynomial is the correct order
821  if(expOrder == iOrder) {
822  double actualOrder = GetFunctionOrder(trialError, trialSpacing);
823 
824  // index into testResults:
825  // stencil varies most slowly, then opDir, then accuracy/order
826  // add 0 for accuracy, 1 for order
827  testResults[iStencil*numDim*2 + opDir*2 + 1] = (actualOrder + epsOrder > expOrder);
828  }
829  } // opDir loop
830  } // Test function order loop
831  } // Stencil loop
832 
833  delete [] numX;
834  delete [] deltaX;
835  delete [] opInterval;
836 
837  return(0);
838  }
839 
840  int CreateStencilConnectivity(int numDim,size_t *dimSizes,size_t *opInterval,int boundaryDepth,
841  int *stencilID, bool fortranInterval)
842  {
843  int yStencil = 0;
844  int zStencil = 0;
845  size_t numPoints = 1;
846  for(int iDim = 0;iDim < numDim;iDim++)
847  numPoints *= dimSizes[iDim];
848  if(numPoints == 0)
849  return(1);
850  if(numDim == 3) {
851  size_t xSize = dimSizes[0];
852  size_t ySize = dimSizes[1];
853  size_t zSize = dimSizes[2];
854  size_t plane = xSize*ySize;
855  size_t xStart = opInterval[0];
856  size_t xEnd = opInterval[1];
857  size_t yStart = opInterval[2];
858  size_t yEnd = opInterval[3];
859  size_t zStart = opInterval[4];
860  size_t zEnd = opInterval[5];
861  if(fortranInterval){
862  xStart--;
863  xEnd--;
864  yStart--;
865  yEnd--;
866  zStart--;
867  zEnd--;
868  }
869  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
870  size_t iPointZ = iZ*plane;
871  size_t rightBoundaryZone = zSize - boundaryDepth;
872  zStencil = 1;
873  if(iZ < boundaryDepth){
874  zStencil = 2+iZ;
875  }
876  if (iZ >= rightBoundaryZone) {
877  if(zStencil != 1)
878  return(1);
879  zStencil = boundaryDepth+1+zSize-iZ;
880  }
881  for(size_t iY = yStart;iY <= yEnd;iY++){
882  size_t iPointYZ = iPointZ + iY*xSize;
883  size_t rightBoundaryZone = ySize - boundaryDepth;
884  yStencil = 1;
885  if(iY < boundaryDepth){
886  yStencil = 2+iY;
887  }
888  if (iY >= rightBoundaryZone) {
889  if(yStencil != 1)
890  return(1);
891  yStencil = boundaryDepth+1+ySize-iY;
892  }
893  for(size_t iX = xStart;iX <= xEnd;iX++){
894  size_t iPoint = iPointYZ + iX;
895  size_t rightBoundaryZone = xSize - boundaryDepth;
896  stencilID[iPoint] = 1;
897  if(iX < boundaryDepth){
898  stencilID[iPoint] = 2+iX;
899  }
900  if (iX >= rightBoundaryZone){
901  stencilID[iPoint] = boundaryDepth+1+xSize-iX;
902  }
903  stencilID[numPoints+iPoint] = yStencil;
904  stencilID[2*numPoints+iPoint] = zStencil;
905  }
906  }
907  }
908  } else if (numDim == 2) {
909  size_t xSize = dimSizes[0];
910  size_t ySize = dimSizes[1];
911  size_t xStart = opInterval[0];
912  size_t xEnd = opInterval[1];
913  size_t yStart = opInterval[2];
914  size_t yEnd = opInterval[3];
915  if(fortranInterval){
916  xStart--;
917  xEnd--;
918  yStart--;
919  yEnd--;
920  }
921  for(size_t iY = yStart;iY <= yEnd;iY++){
922  size_t iPointY = iY*xSize;
923  size_t rightBoundaryZone = ySize - boundaryDepth;
924  yStencil = 1;
925  if(iY < boundaryDepth){
926  yStencil = 2 + iY;
927  }
928  if (iY >= rightBoundaryZone) {
929  yStencil = boundaryDepth+1+ySize-iY;
930  }
931  for(size_t iX = xStart;iX <= xEnd;iX++){
932  size_t iPoint = iPointY + iX;
933  size_t rightBoundaryZone = xSize - boundaryDepth;
934  stencilID[iPoint] = 1;
935  if(iX < boundaryDepth){
936  stencilID[iPoint] = 2 + iX;
937  }
938  if (iX >= rightBoundaryZone){
939  stencilID[iPoint] = boundaryDepth+1+xSize-iX;
940  }
941  stencilID[numPoints+iPoint] = yStencil;
942  }
943  }
944  } else if (numDim == 1) {
945  size_t xSize = dimSizes[0];
946  size_t xStart = opInterval[0];
947  size_t xEnd = opInterval[1];
948  if(fortranInterval){
949  xStart--;
950  xEnd--;
951  }
952  for(size_t iX = xStart;iX <= xEnd;iX++){
953  size_t iPoint = iX;
954  size_t rightBoundaryZone = xSize - boundaryDepth;
955  stencilID[iPoint] = 1;
956  if(iX < boundaryDepth){
957  stencilID[iPoint] = 2 + iX;
958  }
959  if (iX >= rightBoundaryZone){
960  stencilID[iPoint] = boundaryDepth+1+xSize-iX;
961  }
962  }
963  }
964  return(0);
965  }
966 
967  // This takes the global interval in opInterval apparently
968  int CreateStencilConnectivity(int numDim,size_t *dimSizes,size_t *opInterval,int boundaryDepth,
969  int *periodicDirs,int *stencilID, bool fortranInterval)
970  {
971 
972 
973  int yStencil = 0;
974  int zStencil = 0;
975  size_t numPoints = 1;
976  for(int iDim = 0;iDim < numDim;iDim++)
977  numPoints *= dimSizes[iDim];
978  if(numPoints == 0)
979  return(1);
980  if(numDim == 3) {
981  size_t xSize = dimSizes[0];
982  size_t ySize = dimSizes[1];
983  size_t zSize = dimSizes[2];
984  size_t plane = xSize*ySize;
985 
986  size_t xStart = opInterval[0];
987  size_t xEnd = opInterval[1];
988  size_t yStart = opInterval[2];
989  size_t yEnd = opInterval[3];
990  size_t zStart = opInterval[4];
991  size_t zEnd = opInterval[5];
992 
993  if(fortranInterval){
994  xStart--;
995  xEnd--;
996  yStart--;
997  yEnd--;
998  zStart--;
999  zEnd--;
1000  }
1001 
1002  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
1003  size_t iPointZ = iZ*plane;
1004  size_t rightBoundaryZone = zSize - boundaryDepth;
1005 
1006  zStencil = 1;
1007  if(iZ < boundaryDepth){
1008  zStencil = 2+iZ;
1009  }
1010  if (iZ >= rightBoundaryZone) {
1011  if(zStencil != 1)
1012  return(1);
1013  zStencil = boundaryDepth+1+zSize-iZ;
1014  }
1015  if(periodicDirs[2] > 0)
1016  zStencil = 1;
1017  for(size_t iY = yStart;iY <= yEnd;iY++){
1018  size_t iPointYZ = iPointZ + iY*xSize;
1019  size_t rightBoundaryZone = ySize - boundaryDepth;
1020  yStencil = 1;
1021  if(iY < boundaryDepth){
1022  yStencil = 2+iY;
1023  }
1024  if (iY >= rightBoundaryZone) {
1025  if(yStencil != 1)
1026  return(1);
1027  yStencil = boundaryDepth+1+ySize-iY;
1028  }
1029  if(periodicDirs[1] > 0)
1030  yStencil = 1;
1031  for(size_t iX = xStart;iX <= xEnd;iX++){
1032  size_t iPoint = iPointYZ + iX;
1033  size_t rightBoundaryZone = xSize - boundaryDepth;
1034  stencilID[iPoint] = 1;
1035  if(iX < boundaryDepth){
1036  stencilID[iPoint] = 2+iX;
1037  }
1038  if (iX >= rightBoundaryZone){
1039  stencilID[iPoint] = boundaryDepth+1+xSize-iX;
1040  }
1041  if(periodicDirs[0] > 0)
1042  stencilID[iPoint] = 1;
1043  stencilID[numPoints+iPoint] = yStencil;
1044  stencilID[2*numPoints+iPoint] = zStencil;
1045  }
1046  }
1047  }
1048  } else if (numDim == 2) {
1049  size_t xSize = dimSizes[0];
1050  size_t ySize = dimSizes[1];
1051  size_t xStart = opInterval[0];
1052  size_t xEnd = opInterval[1];
1053  size_t yStart = opInterval[2];
1054  size_t yEnd = opInterval[3];
1055  if(fortranInterval){
1056  xStart--;
1057  xEnd--;
1058  yStart--;
1059  yEnd--;
1060  }
1061  for(size_t iY = yStart;iY <= yEnd;iY++){
1062  size_t iPointY = iY*xSize;
1063  size_t rightBoundaryZone = ySize - boundaryDepth;
1064  yStencil = 1;
1065  if(iY < boundaryDepth){
1066  yStencil = 2 + iY;
1067  }
1068  if (iY >= rightBoundaryZone) {
1069  yStencil = boundaryDepth+1+ySize-iY;
1070  }
1071  if(periodicDirs[1] > 0)
1072  yStencil = 1;
1073  for(size_t iX = xStart;iX <= xEnd;iX++){
1074  size_t iPoint = iPointY + iX;
1075  size_t rightBoundaryZone = xSize - boundaryDepth;
1076  stencilID[iPoint] = 1;
1077  if(iX < boundaryDepth){
1078  stencilID[iPoint] = 2 + iX;
1079  }
1080  if (iX >= rightBoundaryZone){
1081  stencilID[iPoint] = boundaryDepth+1+xSize-iX;
1082  }
1083  if(periodicDirs[0] > 0)
1084  stencilID[iPoint] = 1;
1085  stencilID[numPoints+iPoint] = yStencil;
1086  }
1087  }
1088  } else if (numDim == 1) {
1089  size_t xSize = dimSizes[0];
1090  size_t xStart = opInterval[0];
1091  size_t xEnd = opInterval[1];
1092  if(fortranInterval){
1093  xStart--;
1094  xEnd--;
1095  }
1096  for(size_t iX = xStart;iX <= xEnd;iX++){
1097  size_t iPoint = iX;
1098  size_t rightBoundaryZone = xSize - boundaryDepth;
1099  stencilID[iPoint] = 1;
1100  if(iX < boundaryDepth){
1101  stencilID[iPoint] = 2 + iX;
1102  }
1103  if (iX >= rightBoundaryZone){
1104  stencilID[iPoint] = boundaryDepth+1+xSize-iX;
1105  }
1106  if(periodicDirs[0] > 0)
1107  stencilID[iPoint] = 1;
1108  }
1109  }
1110 
1111  return(0);
1112  }
1113 
1114  int BoundaryStencilConnectivity(std::vector<size_t> &bufferSizes,
1115  pcpp::IndexIntervalType &partitionInterval,
1116  pcpp::IndexIntervalType &partitionBufferInterval,
1117  pcpp::IndexIntervalType &boundaryInterval,
1118  int boundaryDepth,int boundaryDirection,
1119  int *stencilConnectivity)
1120 
1121  {
1122  int numDim = bufferSizes.size();
1123 
1124  size_t numPointsBuffer = 1;
1125  for(int iDim = 0;iDim < numDim;iDim++)
1126  numPointsBuffer *= bufferSizes[iDim];
1127  if(numPointsBuffer == 0)
1128  return(1);
1129 
1130  if(boundaryInterval.NNodes() == 0)
1131  return(0);
1132 
1133  pcpp::IndexIntervalType boundaryZoneInterval(boundaryInterval);
1134  int normalDimension = std::abs(boundaryDirection)-1;
1135 
1136  // Make sure the boundary interval is properly specified
1137  assert(boundaryInterval[normalDimension].first == boundaryInterval[normalDimension].second);
1138 
1139  // std::cout << "Boundary Interval: ";
1140  // boundaryInterval.PrettyPrint(std::cout);
1141  // std::cout << std::endl;
1142 
1143  // Create the "boundary zone" (i.e. the region in which
1144  // some boundary stencil will be needed)
1145  if(boundaryDirection > 0){
1146  boundaryZoneInterval[normalDimension].second += (boundaryDepth-1);
1147  } else {
1148  boundaryZoneInterval[normalDimension].first -= (boundaryDepth-1);
1149  }
1150 
1151  // std::cout << "Boundary Zone: ";
1152  // boundaryZoneInterval.PrettyPrint(std::cout);
1153  // std::cout << std::endl;
1154 
1155  // Collide the partitionInterval with the boundaryZone
1156  pcpp::IndexIntervalType overlapInterval;
1157  partitionInterval.Overlap(boundaryZoneInterval,overlapInterval);
1158 
1159  // No overlap, return
1160  if(overlapInterval.NNodes() == 0)
1161  return(0);
1162 
1163  // std::cout << "Partition/Boundary Overlap: ";
1164  // overlapInterval.PrettyPrint(std::cout);
1165  // std::cout << std::endl;
1166 
1167  pcpp::IndexIntervalType bufferInterval;
1168  bufferInterval.InitSimple(bufferSizes);
1169 
1170  // Page to the right beginning spot in the stencilConn
1171  size_t connectivityOffset = normalDimension*numPointsBuffer;
1172  int *connectivityPtr = &stencilConnectivity[connectivityOffset];
1173 
1174  // Determine which stencils will be applied
1175  int firstStencil = 2; // start of boundary stencils
1176  if(boundaryDirection > 0){
1177 
1178  int boundaryDistance = overlapInterval[normalDimension].first -
1179  boundaryInterval[normalDimension].first;
1180  firstStencil = 2 + boundaryDistance;
1181  } else {
1182  int boundaryDistance = boundaryInterval[normalDimension].second -
1183  overlapInterval[normalDimension].first;
1184  firstStencil = (boundaryDepth + 2) + boundaryDistance;
1185  }
1186  // std::cout << "First stencil = " << firstStencil << std::endl;
1187 
1188  // Resolve the buffer interval
1189  pcpp::IndexIntervalType overlapBufferInterval(overlapInterval);
1190  for(int iDim = 0;iDim < numDim;iDim++){
1191  overlapBufferInterval[iDim].first -= partitionInterval[iDim].first;
1192  overlapBufferInterval[iDim].second -= partitionInterval[iDim].first;
1193 
1194  overlapBufferInterval[iDim].first += partitionBufferInterval[iDim].first;
1195  overlapBufferInterval[iDim].second += partitionBufferInterval[iDim].first;
1196 
1197  }
1198 
1199  // std::cout << "Affected Buffer Interval: ";
1200  // overlapBufferInterval.PrettyPrint(std::cout);
1201  // std::cout << std::endl;
1202 
1203  // std::cout << "Updating stencil conn" << std::endl;
1204  // Update the stencil connectivity
1205  int currentStencil = firstStencil;
1206  int stencilIncrement = 1;
1207  if(boundaryDirection < 0)
1208  stencilIncrement = -1;
1209  for(size_t iLayer = overlapBufferInterval[normalDimension].first;
1210  iLayer <= overlapBufferInterval[normalDimension].second;iLayer++){
1211  pcpp::IndexIntervalType layerInterval(overlapBufferInterval);
1212  layerInterval[normalDimension].first = iLayer;
1213  layerInterval[normalDimension].second = iLayer;
1214  std::vector<size_t> layerIndices;
1215  bufferInterval.GetFlatIndices(layerInterval,layerIndices);
1216  std::vector<size_t>::iterator bufferIt = layerIndices.begin();
1217  while(bufferIt != layerIndices.end())
1218  connectivityPtr[*bufferIt++] = currentStencil;
1219  currentStencil += stencilIncrement;
1220  }
1221  // std::cout << "Done updating stencil conn" << std::endl;
1222 
1223  return(0);
1224  }
1225 
1230  int HoleStencilConnectivity(std::vector<size_t> &bufferSizes,
1231  pcpp::IndexIntervalType &partitionInterval,
1232  pcpp::IndexIntervalType &partitionBufferInterval,
1233  pcpp::IndexIntervalType &holeInterval,
1234  int holeStencil,
1235  int *stencilConnectivity)
1236 
1237  {
1238 
1239 
1240  int numDim = bufferSizes.size();
1241 
1242  size_t numPointsBuffer = 1;
1243  for(int iDim = 0;iDim < numDim;iDim++)
1244  numPointsBuffer *= bufferSizes[iDim];
1245  if(numPointsBuffer == 0)
1246  return(1);
1247 
1248  if(holeInterval.NNodes() == 0)
1249  return(0);
1250 
1251  // Collide the partitionInterval with the holeInterval
1252  pcpp::IndexIntervalType overlapInterval;
1253  partitionInterval.Overlap(holeInterval,overlapInterval);
1254 
1255  // No overlap, return
1256  if(overlapInterval.NNodes() == 0){
1257  // std::cout << "No overlap with hole detected." << std::endl;
1258  return(0);
1259  }
1260 
1261  pcpp::IndexIntervalType bufferInterval;
1262  bufferInterval.InitSimple(bufferSizes);
1263 
1264  // Resolve the buffer interval
1265  pcpp::IndexIntervalType overlapBufferInterval(overlapInterval);
1266  for(int iDim = 0;iDim < numDim;iDim++){
1267  overlapBufferInterval[iDim].first -= partitionInterval[iDim].first;
1268  overlapBufferInterval[iDim].second -= partitionInterval[iDim].first;
1269 
1270  overlapBufferInterval[iDim].first += partitionBufferInterval[iDim].first;
1271  overlapBufferInterval[iDim].second += partitionBufferInterval[iDim].first;
1272 
1273  }
1274  // This would have accomplished the same thing
1275  // overlapInterval.RelativeTranslation(partitionInterval,
1276  // partitionBufferInterval,
1277  // overlapBufferInterval);
1278 
1279  std::vector<size_t> bufferIndices;
1280  bufferInterval.GetFlatIndices(overlapBufferInterval,bufferIndices);
1281  // std::cout << "Found " << bufferInterval.NNodes() << " overlapping nodes." << std::endl;
1282  std::vector<size_t>::iterator bufferIt = bufferIndices.begin();
1283  while(bufferIt != bufferIndices.end()){
1284  size_t bufferIndex = *bufferIt++;
1285  for(int iDim = 0;iDim < numDim;iDim++){
1286  stencilConnectivity[iDim*numPointsBuffer+bufferIndex] = holeStencil;
1287  }
1288  }
1289 
1290  return(0);
1291  }
1292 
1293  size_t IntervalSize(int numDim,size_t *opInterval){
1294  size_t retVal = 1;
1295  for(int iDim = 0;iDim < numDim;iDim++)
1296  retVal *= (opInterval[2*iDim+1] - opInterval[2*iDim] + 1);
1297  return(retVal);
1298  };
1299 
1300  size_t BufferSize(int numDim,size_t *dimSize){
1301  size_t retVal = 1;
1302  for(int iDim = 0;iDim < numDim;iDim++)
1303  retVal *= dimSize[iDim];
1304  return(retVal);
1305  };
1306 
1307  int StructuredHole(int numDim,size_t *dimSizes,size_t *opInterval,size_t *holeInterval,
1308  int boundaryDepth,int boundaryStart,int holeBit,int *inMask,int *stencilID)
1309  {
1310  size_t numPoints = 1;
1311  for(int iDim = 0;iDim < numDim;iDim++)
1312  numPoints *= dimSizes[iDim];
1313  if(numPoints == 0)
1314  return(1);
1315 
1316  pcpp::IndexIntervalType bufferExtent(numDim,dimSizes);
1317  pcpp::IndexIntervalType opExtent(opInterval,numDim);
1318  pcpp::IndexIntervalType holeExtent(holeInterval,numDim);
1319  pcpp::IndexIntervalType overlapExtent(opExtent.Overlap(holeExtent));
1320 
1321  if(overlapExtent.NNodes() == 0)
1322  return(0);
1323  int rightBoundaryStart = boundaryStart + boundaryDepth;
1324 
1325  for(int sweepDir = 0;sweepDir < numDim;sweepDir++){
1326  size_t stencilOffset = sweepDir*numPoints;
1327  int *stencilConn = &stencilID[stencilOffset];
1328 
1329  pcpp::IndexIntervalType sweepInterval(overlapExtent);
1330  if(sweepInterval[sweepDir].first < boundaryDepth)
1331  sweepInterval[sweepDir].first = boundaryDepth;
1332  if(sweepInterval[sweepDir].second > (bufferExtent[sweepDir].second-boundaryDepth))
1333  sweepInterval[sweepDir].second = bufferExtent[sweepDir].second-boundaryDepth;
1334  if(sweepInterval[sweepDir].first > sweepInterval[sweepDir].second)
1335  continue;
1336 
1337  size_t sweepOffset = 1;
1338  int offDir = sweepDir;
1339  while(offDir > 0){
1340  sweepOffset *= dimSizes[offDir-1];
1341  offDir--;
1342  }
1343 
1344  // Get buffer indices of sweep region
1345  std::vector<size_t> sweepIndices;
1346  bufferExtent.GetFlatIndices(sweepInterval,sweepIndices);
1347 
1348  // Loop over all sweep indices and check for holes/boundaries
1349  std::vector<size_t>::iterator sweepIt = sweepIndices.begin();
1350  while(sweepIt != sweepIndices.end()){
1351  size_t sweepIndex = *sweepIt++;
1352 
1353  // Skip holes
1354  if(inMask[sweepIndex]&holeBit)
1355  continue;
1356 
1357  if(inMask[sweepIndex-sweepOffset]&holeBit){ // found left boundary
1358  // Set left boundary stencils
1359  for(int iBoundary = 0;iBoundary < boundaryDepth;iBoundary++){
1360  size_t connIndex = sweepIndex + iBoundary*sweepOffset;
1361  if(stencilConn[connIndex] != 1)
1362  return(1);
1363  stencilConn[connIndex] = boundaryStart+iBoundary;
1364  }
1365  }
1366 
1367  if(inMask[sweepIndex+sweepOffset]&holeBit){ // found right boundary
1368  // Set right boundary stencils
1369  for(int iBoundary = 0;iBoundary < boundaryDepth;iBoundary++){
1370  size_t connIndex = sweepIndex - iBoundary*sweepOffset;
1371  if(stencilConn[connIndex] != 1)
1372  return(1);
1373  stencilConn[connIndex] = rightBoundaryStart+iBoundary;
1374  }
1375  }
1376  }
1377  }
1378  return(0);
1379  }
1380 
1381  int MaskStencilConnectivity(int numDim,size_t *dimSizes,size_t *opInterval,
1382  int boundaryDepth,int boundaryStart,int holeBit,
1383  int *inMask,int *stencilID)
1384  {
1385 
1386  size_t numPoints = 1;
1387  for(int iDim = 0;iDim < numDim;iDim++)
1388  numPoints *= dimSizes[iDim];
1389  if(numPoints == 0)
1390  return(1);
1391 
1392  pcpp::IndexIntervalType bufferExtent(numDim,dimSizes);
1393  pcpp::IndexIntervalType opExtent(opInterval,numDim);
1394 
1395  int rightBoundaryStart = boundaryStart + boundaryDepth;
1396 
1397  for(int sweepDir = 0;sweepDir < numDim;sweepDir++){
1398  size_t stencilOffset = sweepDir*numPoints;
1399  int *stencilConn = &stencilID[stencilOffset];
1400 
1401  pcpp::IndexIntervalType sweepInterval(opExtent);
1402  if(sweepInterval[sweepDir].first < boundaryDepth)
1403  sweepInterval[sweepDir].first = boundaryDepth;
1404  if(sweepInterval[sweepDir].second > (bufferExtent[sweepDir].second-boundaryDepth))
1405  sweepInterval[sweepDir].second = bufferExtent[sweepDir].second-boundaryDepth;
1406  if(sweepInterval[sweepDir].first > sweepInterval[sweepDir].second)
1407  continue;
1408 
1409  size_t sweepOffset = 1;
1410  int offDir = sweepDir;
1411  while(offDir > 0){
1412  sweepOffset *= dimSizes[offDir-1];
1413  offDir--;
1414  }
1415 
1416  // Get buffer indices of sweep region
1417  std::vector<size_t> sweepIndices;
1418  bufferExtent.GetFlatIndices(sweepInterval,sweepIndices);
1419 
1420  // Loop over all sweep indices and check for holes/boundaries
1421  std::vector<size_t>::iterator sweepIt = sweepIndices.begin();
1422  while(sweepIt != sweepIndices.end()){
1423  size_t sweepIndex = *sweepIt++;
1424 
1425  // Skip holes
1426  if(inMask[sweepIndex]&holeBit)
1427  continue;
1428 
1429  if(inMask[sweepIndex-sweepOffset]&holeBit){ // found left boundary
1430  // Set left boundary stencils
1431  for(int iBoundary = 0;iBoundary < boundaryDepth;iBoundary++){
1432  size_t connIndex = sweepIndex + iBoundary*sweepOffset;
1433  if(stencilConn[connIndex] != 1)
1434  return(1);
1435  stencilConn[connIndex] = boundaryStart+iBoundary;
1436  }
1437  }
1438 
1439  if(inMask[sweepIndex+sweepOffset]&holeBit){ // found right boundary
1440  // Set right boundary stencils
1441  for(int iBoundary = 0;iBoundary < boundaryDepth;iBoundary++){
1442  size_t connIndex = sweepIndex - iBoundary*sweepOffset;
1443  if(stencilConn[connIndex] != 1)
1444  return(1);
1445  stencilConn[connIndex] = rightBoundaryStart+iBoundary;
1446  }
1447  }
1448  }
1449  }
1450  return(0);
1451  }
1452 
1453 
1454  int DetectHoles(int numDim,size_t *dimSizes,size_t *opExtent,
1455  int boundaryDepth,int holeBit,int *inMask,int *stencilID)
1456  {
1457 
1458  size_t numPointsBuffer = 1;
1459  for(int iDim = 0;iDim < numDim;iDim++)
1460  numPointsBuffer *= dimSizes[iDim];
1461  if(numPointsBuffer == 0)
1462  return(1);
1463 
1464  pcpp::IndexIntervalType bufferInterval(numDim,dimSizes);
1465  pcpp::IndexIntervalType opInterval(opExtent,numDim);
1466 
1467  // std::cout << "DetectHoles: BufferInterval: ";
1468  // bufferInterval.PrettyPrint(std::cout);
1469  // std::cout << std::endl << "DetectHoles: OpInterval: ";
1470  // opInterval.PrettyPrint(std::cout);
1471  // std::cout << std::endl;
1472 
1473  std::vector<size_t> nodeList;
1474  bufferInterval.GetFlatIndices(opInterval,nodeList);
1475 
1476  int numStencils = 2*(boundaryDepth+1);
1477  int holeStencil = numStencils;
1478 
1479  bool holeFound = false;
1480  if(nodeList.size() > 0){
1481  // Before the sweeps start - set the stencilIDs for all the masked-out
1482  // points.
1483 
1484  for(std::vector<size_t>::iterator nodeIt = nodeList.begin();
1485  nodeIt != nodeList.end(); nodeIt++){
1486  size_t nodeId = *nodeIt;
1487  if(inMask[nodeId] & holeBit){
1488  holeFound = true;
1489  for(int iDim = 0;iDim < numDim;iDim++)
1490  stencilID[iDim*numPointsBuffer+nodeId] = holeStencil;
1491  }
1492  }
1493 
1494  if(numDim == 1){
1495  std::vector<size_t>::iterator nodeIt = nodeList.begin();
1496  while(nodeIt != nodeList.end()){
1497  size_t nodeId = *nodeIt++;
1498  if(inMask[nodeId] & holeBit)
1499  continue;
1500 
1501  bool leftBoundary = false;
1502  bool rightBoundary = false;
1503 
1504  for(int iSearch = 1;((iSearch <= boundaryDepth));iSearch++){
1505 
1506  // search left
1507  if(iSearch <= nodeId && !leftBoundary){
1508  size_t searchNodeId = nodeId - iSearch;
1509  if((inMask[searchNodeId] & holeBit)){
1510  if(rightBoundary) return(1); // fail - not enough room for boundary stencils
1511  stencilID[nodeId] = 2 + iSearch - 1;
1512  leftBoundary = true;
1513  }
1514  }
1515 
1516  // search right
1517  if(((iSearch + nodeId) < numPointsBuffer) && !rightBoundary){
1518  size_t searchNodeId = nodeId + iSearch;
1519  if((inMask[searchNodeId] & holeBit)){
1520  if(leftBoundary) return(1); // fail - not enough room for boundary stencils
1521  stencilID[nodeId] = 2 + boundaryDepth + iSearch - 1;
1522  rightBoundary = true;
1523  }
1524  }
1525 
1526  // Other
1527  if(((nodeId < boundaryDepth) && rightBoundary) ||
1528  ((nodeId > (dimSizes[0]-boundaryDepth)) && leftBoundary))
1529  return(1);
1530 
1531  }
1532 
1533  }
1534  } else if(numDim == 2){
1535 
1536  size_t iStart = opInterval[0].first;
1537  size_t iEnd = opInterval[0].second;
1538  size_t jStart = opInterval[1].first;
1539  size_t jEnd = opInterval[1].second;
1540 
1541  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
1542  size_t jOffset = jIndex*dimSizes[0];
1543  for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
1544  size_t nodeId = jOffset + iIndex;
1545  if(inMask[nodeId] & holeBit)
1546  continue;
1547 
1548  std::vector<bool> leftBoundary(2,false);
1549  std::vector<bool> rightBoundary(2,false);
1550 
1551  for(int iSearch = 1;((iSearch <= boundaryDepth));iSearch++){
1552 
1553  // ===== SWEEP IN I =======
1554  // Scan left
1555  if(iSearch <= iIndex && !leftBoundary[0]){
1556  size_t iSearchIndex = iIndex - iSearch;
1557  size_t searchNodeId = jOffset + iSearchIndex;
1558  if(inMask[searchNodeId] & holeBit){
1559  if(rightBoundary[0]) return(1); // fail - not enough room for boundary stencils
1560  stencilID[nodeId] = 2 + iSearch - 1;
1561  leftBoundary[0] = true;
1562  }
1563  }// scan left
1564 
1565  // Scan right
1566  if(((iSearch + iIndex) < dimSizes[0]) && !rightBoundary[0]){
1567  size_t iSearchIndex = iIndex + iSearch;
1568  size_t searchNodeId = jOffset + iSearchIndex;
1569  if(inMask[searchNodeId] & holeBit){
1570  if(leftBoundary[0]) return(1); // fail - not enough room for boundary stencils
1571  stencilID[nodeId] = 2 + boundaryDepth + iSearch - 1;
1572  rightBoundary[0] = true;
1573  }
1574  } // scan right
1575 
1576  // Other
1577  if(((iIndex < boundaryDepth) && rightBoundary[0]) ||
1578  ((iIndex > (dimSizes[0]-boundaryDepth)) && leftBoundary[0]))
1579  return(1);
1580 
1581  // ======= SWEEP IN J ========
1582  // Scan left
1583  if(iSearch <= jIndex && !leftBoundary[1]){
1584  size_t jSearchIndex = jIndex - iSearch;
1585  size_t searchNodeId = jSearchIndex*dimSizes[0] + iIndex;
1586  if(inMask[searchNodeId] & holeBit){
1587  if(rightBoundary[1]) return(1); // fail - not enough room for boundary stencils
1588  stencilID[numPointsBuffer+nodeId] = 2 + iSearch - 1;
1589  leftBoundary[1] = true;
1590  }
1591  }// scan left
1592 
1593  // Scan right
1594  if(((iSearch + jIndex) < dimSizes[1]) && !rightBoundary[1]){
1595  size_t jSearchIndex = jIndex + iSearch;
1596  size_t searchNodeId = jSearchIndex*dimSizes[0] + iIndex;
1597  if(inMask[searchNodeId] & holeBit){
1598  if(leftBoundary[1]) return(1); // fail - not enough room for boundary stencils
1599  stencilID[numPointsBuffer+nodeId] = 2 + boundaryDepth + iSearch - 1;
1600  rightBoundary[1] = true;
1601  }
1602  } // scan right
1603 
1604  // Other
1605  if(((jIndex < boundaryDepth) && rightBoundary[1]) ||
1606  ((jIndex > (dimSizes[1]-boundaryDepth)) && leftBoundary[1]))
1607  return(1);
1608 
1609  } // loop over search size
1610  } // iIndex
1611  } // jIndex
1612 
1613  } else {
1614 
1615  size_t iStart = opInterval[0].first;
1616  size_t iEnd = opInterval[0].second;
1617  size_t jStart = opInterval[1].first;
1618  size_t jEnd = opInterval[1].second;
1619  size_t kStart = opInterval[2].first;
1620  size_t kEnd = opInterval[2].second;
1621  size_t nX = dimSizes[0];
1622  size_t nPlane = nX*dimSizes[1];
1623 
1624  for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
1625  size_t kOffset = kIndex*nPlane;
1626  for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
1627  size_t jOffset = jIndex*nX;
1628  size_t jkOffset = jOffset+kOffset;
1629  for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
1630 
1631  size_t nodeId = jkOffset + iIndex;
1632  if(inMask[nodeId] & holeBit)
1633  continue;
1634 
1635  std::vector<bool> leftBoundary(3,false);
1636  std::vector<bool> rightBoundary(3,false);
1637 
1638 
1639  for(int iSearch = 1;((iSearch <= boundaryDepth));iSearch++){
1640 
1641  // ===== SWEEP IN I =======
1642  // Scan left
1643  if(iSearch <= iIndex && !leftBoundary[0]){
1644  size_t iSearchIndex = iIndex - iSearch;
1645  size_t searchNodeId = jkOffset + iSearchIndex;
1646  if(inMask[searchNodeId] & holeBit){
1647  if(rightBoundary[0]) return(1); // fail - not enough room for boundary stencils
1648  stencilID[nodeId] = 2 + iSearch - 1;
1649  leftBoundary[0] = true;
1650  }
1651  }// scan left
1652 
1653  // Scan right
1654  if(((iSearch + iIndex) < dimSizes[0]) && !rightBoundary[0]){
1655  size_t iSearchIndex = iIndex + iSearch;
1656  size_t searchNodeId = jkOffset + iSearchIndex;
1657  if(inMask[searchNodeId] & holeBit){
1658  if(leftBoundary[0]) return(1); // fail - not enough room for boundary stencils
1659  stencilID[nodeId] = 2 + boundaryDepth + iSearch - 1;
1660  rightBoundary[0] = true;
1661  }
1662  } // scan right
1663 
1664  // Other
1665  if(((iIndex < boundaryDepth) && rightBoundary[0]) ||
1666  ((iIndex > (dimSizes[0]-boundaryDepth)) && leftBoundary[0]))
1667  return(1);
1668 
1669  // ======= SWEEP IN J ========
1670  // Scan left
1671  if(iSearch <= jIndex && !leftBoundary[1]){
1672  size_t jSearchIndex = jIndex - iSearch;
1673  size_t searchNodeId = jSearchIndex*nX + kOffset + iIndex;
1674  if(inMask[searchNodeId] & holeBit){
1675  if(rightBoundary[1]) return(1); // fail - not enough room for boundary stencils
1676  stencilID[numPointsBuffer+nodeId] = 2 + iSearch - 1;
1677  leftBoundary[1] = true;
1678  }
1679  }// scan left
1680 
1681  // Scan right
1682  if(((iSearch + jIndex) < dimSizes[1]) && !rightBoundary[1]){
1683  size_t jSearchIndex = jIndex + iSearch;
1684  size_t searchNodeId = jSearchIndex*nX + kOffset + iIndex;
1685  if(inMask[searchNodeId] & holeBit){
1686  if(leftBoundary[1]) return(1); // fail - not enough room for boundary stencils
1687  stencilID[numPointsBuffer+nodeId] = 2 + boundaryDepth + iSearch - 1;
1688  rightBoundary[1] = true;
1689  }
1690  } // scan right
1691 
1692  // Other
1693  if(((jIndex < boundaryDepth) && rightBoundary[1]) ||
1694  ((jIndex > (dimSizes[1]-boundaryDepth)) && leftBoundary[1]))
1695  return(1);
1696 
1697  // ======= SWEEP IN K ========
1698  // Scan left
1699  if(iSearch <= kIndex && !leftBoundary[2]){
1700  size_t kSearchIndex = kIndex - iSearch;
1701  size_t searchNodeId = kSearchIndex*nPlane + jOffset + iIndex;
1702  if(inMask[searchNodeId] & holeBit){
1703  if(rightBoundary[2]) return(1); // fail - not enough room for boundary stencils
1704  stencilID[2*numPointsBuffer+nodeId] = 2 + iSearch - 1;
1705  leftBoundary[2] = true;
1706  }
1707  }// scan left
1708 
1709  // Scan right
1710  if(((iSearch + kIndex) < dimSizes[2]) && !rightBoundary[2]){
1711  size_t kSearchIndex = kIndex + iSearch;
1712  size_t searchNodeId = kSearchIndex*nPlane + jOffset + iIndex;
1713  if(inMask[searchNodeId] & holeBit){
1714  if(leftBoundary[2]) return(1); // fail - not enough room for boundary stencils
1715  stencilID[2*numPointsBuffer+nodeId] = 2 + boundaryDepth + iSearch - 1;
1716  rightBoundary[2] = true;
1717  }
1718  } // scan right
1719 
1720  // Other
1721  if(((kIndex < boundaryDepth) && rightBoundary[2]) ||
1722  ((kIndex > (dimSizes[2]-boundaryDepth)) && leftBoundary[2]))
1723  return(1);
1724 
1725  } // loop over search size
1726  } // iIndex
1727  } // jIndex
1728  } // kIndex
1729  } // 3-dimensions
1730  } // if there are nodes in the interval
1731 
1732  // if(!holeFound)
1733  // std::cout << "DetectHoles: No holes found." << std::endl;
1734 
1735  return(0);
1736  }
1737 
1738 
1739  int InvertStencilConnectivity(int numDim,size_t *dimSizes,size_t *opInterval,int numStencils,int *stencilID,
1740  size_t *dualStencilConn,size_t *numPointsStencil,bool fortranInterval)
1741  {
1742  size_t numPointsBuffer = BufferSize(numDim,dimSizes);
1743  size_t numPointsInterval = IntervalSize(numDim,opInterval);
1744 
1745  size_t *myInterval = opInterval;
1746  if(fortranInterval){
1747  myInterval = new size_t [2*numDim];
1748  for(int iDim = 0;iDim < 2*numDim;iDim++){
1749  myInterval[iDim] = opInterval[iDim]-1;
1750  }
1751  }
1752 
1753  // Invert the connectivity
1754  // [DIM1[stencil1Points,stencil2Points....stencil(numStencils)Points]...]
1755  // size_t *dualStencilConn = new size_t [numDim*numPoints];
1756  // [DIM1[numPoints1,numPoints2....numPoints(numStencils)]...]
1757  // size_t *numPointsStencil = new size_t [numDim*numStencils];
1758 
1759 
1760  for(int iStencil = 0;iStencil < numDim*numStencils;iStencil++)
1761  numPointsStencil[iStencil] = 0;
1762  std::vector<std::list<size_t> > stencilPoints(numDim*numStencils);
1763 
1764  if(numDim == 3){
1765 
1766  size_t xStart = myInterval[0];
1767  size_t xEnd = myInterval[1];
1768  size_t yStart = myInterval[2];
1769  size_t yEnd = myInterval[3];
1770  size_t zStart = myInterval[4];
1771  size_t zEnd = myInterval[5];
1772  size_t plane = dimSizes[0]*dimSizes[1];
1773 
1774  for(size_t iZ = zStart;iZ <= zEnd;iZ++){
1775  size_t zIndex = iZ*plane;
1776  for(size_t iY = yStart;iY <= yEnd;iY++){
1777  size_t yzIndex = zIndex + iY*dimSizes[0];
1778  for(size_t iX = xStart;iX <= xEnd;iX++){
1779  size_t iPoint = yzIndex + iX;
1780 
1781  int iStencilX = stencilID[iPoint];
1782  int iStencilY = stencilID[iPoint+numPointsBuffer];
1783  int iStencilZ = stencilID[iPoint+numPointsBuffer*2];
1784 
1785  numPointsStencil[iStencilX-1]++;
1786  numPointsStencil[iStencilY+numStencils-1]++;
1787  numPointsStencil[iStencilZ+2*numStencils-1]++;
1788 
1789  stencilPoints[iStencilX-1].push_back(iPoint);
1790  stencilPoints[numStencils+iStencilY-1].push_back(iPoint);
1791  stencilPoints[2*numStencils+iStencilZ-1].push_back(iPoint);
1792  }
1793  }
1794  }
1795  } else if (numDim == 2) {
1796 
1797  size_t xStart = myInterval[0];
1798  size_t xEnd = myInterval[1];
1799  size_t yStart = myInterval[2];
1800  size_t yEnd = myInterval[3];
1801  for(size_t iY = yStart;iY <= yEnd;iY++){
1802  size_t yIndex = iY*dimSizes[0];
1803  for(size_t iX = xStart;iX <= xEnd;iX++){
1804  size_t iPoint = yIndex + iX;
1805 
1806  int iStencilX = stencilID[iPoint];
1807  int iStencilY = stencilID[iPoint+numPointsBuffer];
1808 
1809  numPointsStencil[iStencilX-1]++;
1810  numPointsStencil[iStencilY+numStencils-1]++;
1811 
1812  stencilPoints[iStencilX-1].push_back(iPoint);
1813  stencilPoints[numStencils+iStencilY-1].push_back(iPoint);
1814  }
1815  }
1816  }
1817 
1818  // Actually populate the dual
1819  for(int iDim = 0;iDim < numDim;iDim++){
1820  size_t iPoint = 0;
1821  for(int iStencil = 0;iStencil < numStencils;iStencil++){
1822  std::list<size_t> &pointList(stencilPoints[iDim*numStencils+iStencil]);
1823  std::list<size_t>::iterator stencilPointIt = pointList.begin();
1824  while(stencilPointIt != pointList.end()){
1825  dualStencilConn[iDim*numPointsInterval+iPoint++] = *stencilPointIt++;
1826  }
1827  }
1828  }
1829 
1830  if(fortranInterval){
1831  // Forticate the dual (i.e. point ids = C point index + 1)
1832  for(size_t iPoint = 0;iPoint < numDim*numPointsInterval;iPoint++)
1833  dualStencilConn[iPoint]++;
1834  delete [] myInterval;
1835  }
1836 
1837  return(0);
1838 
1839  }
1840  }
1841  }
1842 
1843 }
1844 
1845 namespace mask {
1846 
1847  void SetMask(const std::vector<size_t> &opIndices,int maskBits,int *inMask)
1848  {
1849  std::vector<size_t>::const_iterator bufferIt = opIndices.begin();
1850  while(bufferIt != opIndices.end()){
1851  size_t bufferIndex = *bufferIt++;
1852  inMask[bufferIndex] |= maskBits;
1853  }
1854  };
1855 
1856  void SetMask(const std::vector<size_t> &bufferSizes,
1857  const pcpp::IndexIntervalType &opInterval,
1858  int maskBits,int *inMask)
1859  {
1860  pcpp::IndexIntervalType bufferInterval;
1861  bufferInterval.InitSimple(bufferSizes);
1862 
1863  std::vector<size_t> opIndices;
1864  bufferInterval.GetFlatIndices(opInterval,opIndices);
1865 
1866  return(SetMask(opIndices,maskBits,inMask));
1867  };
1868 
1869  void UnSetMask(const std::vector<size_t> &opIndices,int maskBits,int *inMask)
1870  {
1871  std::vector<size_t>::const_iterator bufferIt = opIndices.begin();
1872  while(bufferIt != opIndices.end()){
1873  size_t bufferIndex = *bufferIt++;
1874  if(inMask[bufferIndex]&maskBits)
1875  inMask[bufferIndex] ^= maskBits;
1876  }
1877  };
1878 
1879  void UnSetMask(const std::vector<size_t> &bufferSizes,
1880  const pcpp::IndexIntervalType &opInterval,
1881  int maskBits,int *inMask)
1882  {
1883  pcpp::IndexIntervalType bufferInterval;
1884  bufferInterval.InitSimple(bufferSizes);
1885 
1886  std::vector<size_t> opIndices;
1887  bufferInterval.GetFlatIndices(opInterval,opIndices);
1888 
1889  return(UnSetMask(opIndices,maskBits,inMask));
1890  };
1891 }
plascom2::operators::sbp::InvertStencilConnectivity
int InvertStencilConnectivity(int numDim, size_t *dimSizes, size_t *opInterval, int numStencils, int *stencilID, size_t *dualStencilConn, size_t *numPointsStencil, bool fortranInterval=false)
Inverts stencil connectivity to populate the so-called dual stencil connectivity
Definition: Stencil.C:1739
ix::util::sizeextent::GetFlatIndices
void GetFlatIndices(const sizeextent &extent, ContainerType &indices) const
Definition: IndexUtil.H:302
plascom2::operators::stencilset::stencilSizes
int * stencilSizes
The number of weights for each stencil.
Definition: Stencil.H:102
numPoints
void const size_t * numPoints
Definition: EulerKernels.H:10
opDir
void size_t int size_t int * opDir
Definition: OperatorKernels.H:9
operators::applyoperator
subroutine applyoperator(numDim, dimSizes, numComponents, numPoints, opDir, opInterval, numStencils, stencilSizes, stencilStarts, numValues, stencilWeights, stencilOffsets, stencilID, U, dU)
applyoperator applies an operator specified as a stencil set to the provided state data ...
Definition: Operators.f90:36
plascom2::operators::sbp::BoundaryStencilConnectivity
int BoundaryStencilConnectivity(std::vector< size_t > &bufferSizes, pcpp::IndexIntervalType &partitionInterval, pcpp::IndexIntervalType &partitionBufferInterval, pcpp::IndexIntervalType &boundaryInterval, int boundaryDepth, int boundaryDirection, int *stencilConnectivity)
Update a stencil connectivity with a boundary.
Definition: Stencil.C:1114
plascom2::operators::stencilset::stencilStarts
int * stencilStarts
The starting index into the stencilWeight and stencilOffset arrays for each stencil.
Definition: Stencil.H:104
plascom2::operators::stencilset::stencilOffsets
int * stencilOffsets
The offsets wrt the grid point at which the stencil is being applied.
Definition: Stencil.H:106
plascom2::operators::sbp::IntervalSize
size_t IntervalSize(int numDim, size_t *opInterval)
Definition: Stencil.C:1293
plascom2
Definition: PC2IO.H:10
opInterval
void const size_t const size_t const size_t * opInterval
Definition: EulerKernels.H:10
plascom2::operators::sbp::BruteTest1
int BruteTest1(base &inOperator, int interiorOrder, int numDim, int numComponents, int numTrials, std::vector< bool > &testResults)
Brute-force accuracy test for SBP operators.
Definition: Stencil.C:625
plascom2::operators::sbp::BufferSize
size_t BufferSize(int numDim, size_t *dimSize)
Definition: Stencil.C:1300
ix::util::sizeextent::NNodes
size_t NNodes() const
Definition: IndexUtil.H:254
dimSizes
void double double double double double double double *void const size_t * dimSizes
Definition: SATKernels.H:48
plascom2::operators::sbp::DetectHoles
int DetectHoles(int numDim, size_t *dimSizes, size_t *opInterval, int boundaryDepth, int holeBit, int *inMask, int *stencilID)
Detect unstructured holes and set up boundary stencil id&#39;s accordingly.
Definition: Stencil.C:1454
plascom2::operators::sbp::MaskStencilConnectivity
int MaskStencilConnectivity(int numDim, size_t *dimSizes, size_t *opInterval, int boundaryDepth, int boundaryStart, int holeBit, int *inMask, int *stencilID)
Definition: Stencil.C:1381
plascom2::operators::sbp::HoleStencilConnectivity
int HoleStencilConnectivity(std::vector< size_t > &bufferSizes, pcpp::IndexIntervalType &partitionInterval, pcpp::IndexIntervalType &partitionBufferInterval, pcpp::IndexIntervalType &holeInterval, int holeStencil, int *stencilConnectivity)
Update a stencil connectivity with a hole.
Definition: Stencil.C:1230
mask::UnSetMask
void UnSetMask(const std::vector< size_t > &opIndices, int maskBits, int *inMask)
Definition: Stencil.C:1869
plascom2::operators::dissipation::Initialize
int Initialize(stencilset &stencilSet1, stencilset &stencilSet2, int interiorOrder)
Definition: Stencil.C:29
stencilID
void size_t int size_t int size_t int int int int double int int * stencilID
Definition: OperatorKernels.H:9
plascom2::operators::stencilset::stencilOrders
int * stencilOrders
Boundary weight needed by BC&#39;s.
Definition: Stencil.H:110
plascom2::operators::sbp::CreateStencilConnectivity
int CreateStencilConnectivity(int numDim, size_t *dimSizes, size_t *opInterval, int boundaryDepth, int *stencilID, bool fortranInterval=false)
Creates simple stencil connectivity assuming all domain boundaries are physical boundaries.
Definition: Stencil.C:840
plascom2::operators::stencilset
Encapsulation for a collection of operator stencils.
Definition: Stencil.H:26
plascom2::operators::stencilset::numValues
int numValues
The total number of weights for all stencils (reqd for Fortran)
Definition: Stencil.H:98
ix::util::sizeextent::Overlap
void Overlap(const sizeextent &inextent, sizeextent &outextent) const
Definition: IndexUtil.H:324
bufferSizes
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
GetFunctionOrder
double GetFunctionOrder(std::vector< double > &functionValues, std::vector< double > &independentVariable)
Definition: Stencil.C:5
bufferInterval
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t * bufferInterval
Definition: OperatorKernels.H:29
plascom2::operators::stencilset::Destroy
void Destroy()
Destroy utility destroys the stencil set memory.
Definition: Stencil.H:73
plascom2::operators::sbp::IsValidOrder
bool IsValidOrder(int overallOrder)
Definition: Stencil.C:352
plascom2::operators::stencilset::ownData
bool ownData
Indicates whether this data structure owns the stencil memory.
Definition: Stencil.H:100
error
Definition: EulerRHS.H:33
plascom2::operators::stencilset::boundaryWidth
int boundaryWidth
Boundary width is the size of the on-boundary stencil.
Definition: Stencil.H:116
plascom2::operators::stencilset::stencilWeights
double * stencilWeights
The stencil weights.
Definition: Stencil.H:108
numStencils
void size_t int size_t int size_t int * numStencils
Definition: OperatorKernels.H:9
plascom2::operators::stencilset::numStencils
int numStencils
The number of stencils (e.g. interior + boundary)
Definition: Stencil.H:93
plascom2::operators::stencilset::boundaryWeight
double boundaryWeight
The order of accuracy for each stencil.
Definition: Stencil.H:112
mask
Definition: Stencil.H:328
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
numComponents
void size_t int * numComponents
Definition: OperatorKernels.H:9
plascom2::operators::sbp::OperatorSetup
void OperatorSetup(base &stencilSet, int interiorOrder, int boundaryDepth, int boundaryWidth, double boundaryWeight, std::vector< double > centralWeightsRight, std::vector< std::vector< double > > leftWeights, std::vector< int > leftOrders)
Setup the SBP operator given the appropriate coefficients – called from Initialize.
Definition: Stencil.C:485
plascom2::operators::stencilset::boundaryDepth
int boundaryDepth
Boundary depth is the number of biased boundary stencils for one boundary.
Definition: Stencil.H:114
mask::SetMask
void SetMask(const std::vector< size_t > &opIndices, int maskBits, int *inMask)
Definition: Stencil.C:1847
plascom2::operators::stencilset::overallOrder
int overallOrder
The overall order of the scheme this stencil implements.
Definition: Stencil.H:118
ix::util::sizeextent::InitSimple
void InitSimple(const ContainerType &inSize)
Definition: IndexUtil.H:169
numPointsBuffer
void const size_t * numPointsBuffer
Definition: MetricKernels.H:19
FC_MODULE
void FC_MODULE(euler, flux1d, EULER, FLUX1D)(const int *numDim
plascom2::operators::sbp::StructuredHole
int StructuredHole(int numDim, size_t *dimSizes, size_t *opInterval, size_t *holeInterval, int boundaryDepth, int boundaryStart, int holeBit, int *inMask, int *stencilID)
Definition: Stencil.C:1307