OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
grid.C
Go to the documentation of this file.
1 #include "grid.h"
2 #include "error.h"
3 #include <stdlib.h>
4 #include <stdio.h>
5 #include <math.h>
6
7 namespace oofem {
8 #define MIN(a, b) ( ( a ) > ( b ) ? ( b ) : ( a ) )
9 #define MAX(a, b) ( ( a ) > ( b ) ? ( a ) : ( b ) )
10
11 // NB: When using matrix indices (i,j), they run from 1 to (m,n).
12 // When using 1-dim array index 'ind', it runs from 0 to (m*n-1).
13
14 // Index offsets of all neighbors (direct and diagonal ones)
15 int iOffsets_full[] = {
16  -1, -1, -1, 0, 0, 1, 1, 1
17 };
18 int jOffsets_full[] = {
19  -1, 0, 1, -1, 1, -1, 0, 1
20 };
21
22 // Indicator which neighbors are diagonal ones
23 bool is_diag[] = {
24  true, false, true, false, false, true, false, true
25 };
26
27 // Index offsets of direct neighbours
28 int iOffsets[] = {
29  -1, 1, 0, 0
30 };
31 int jOffsets[] = {
32  0, 0, -1, 1
33 };
34
35 // Offsets used by calcTime
36 // NB: Don't change the ordering without also changing the if-else
37 // statements below that set xmin, ymin.
38 int icalcOffsets[] = {
39  -2, -1, 1, 2, 0, 0, 0, 0
40 };
41 int jcalcOffsets[] = {
42  0, 0, 0, 0, -2, -1, 1, 2
43 };
44
45
46 Grid :: Grid(int N, int M)
47 {
48  n = N;
49  m = M;
50  solutionAvailable = false;
51  T = new FloatMatrix(m, n);
52  F = new FloatMatrix(m, n);
53
54  // (calloc sets all bits of the allocated memory to zero, so the initial values are "false")
55  Frozen = ( bool * ) calloc( m * n, sizeof( bool ) );
56  narrowBand = new Heap(m *n);
57 }
58
60 {
61  delete T;
62  delete F;
63  if ( Frozen ) {
64  free(Frozen);
65  }
66  delete narrowBand;
67 }
68
69 int
71 {
72  switch ( dir ) {
73  case 1: return n;
74
75  case 2: return m;
76
77  default: return 0;
78  }
79 }
80
81 /*
82  * void
83  * Grid :: setPrescribedField(double *PF)
84  * {
85  * // For efficiency, F is passed by its pointer and from now on
86  * // will be taken care of (and deleted) by the Grid.
87  * if ( F ) {
88  * delete F;
89  * }
90  * F = PF;
91  * }
92  */
93
94 void
96 {
97  for ( int ind = 0; ind < ( m * n ); ind++ ) {
98  Frozen [ ind ] = false;
99  }
100  narrowBand->setToEmpty(m * n);
101  solutionAvailable = false;
102 }
103
104 void
106 {
107  // gridCoords contains coordinates of points at which the solution is set to zero.
108  // These coordinates must be expressed the grid coordinate system,
109  // with the first coordinate (j) from 1 to n and the second (i) from 1 to m.
110  // However, they do not need to be integers - the points do not need
111  // to coincide with grid nodes.
112
113  if ( gridCoords->giveNumberOfColumns() != 2 ) {
114  OOFEM_ERROR("Matrix gridCoords passed to Grid :: setZeroValues should have 2 columns.");
115  }
116  int nValues = gridCoords->giveNumberOfRows();
117  for ( int ival = 1; ival <= nValues; ival++ ) {
118  // get coordinates of the specified point
119  double CPx = gridCoords->at(ival, 1);
120  double CPy = gridCoords->at(ival, 2);
121  // find nearest grid node (i, j)
122  double j = ( int ) ceil(CPx - 0.5);
123  if ( j < 1 ) {
124  j = 1;
125  } else if ( j > n ) {
126  j = n;
127  }
128  double i = ( int ) ceil(CPy - 0.5);
129  if ( i < 1 ) {
130  i = 1;
131  } else if ( i > m ) {
132  i = m;
133  }
134  // determine the index of the nearest grid node
135  int CPInd = ij2ind(i, j, m);
136
137  // Calculate time-distance of nearest grid node and freeze the value
138  double Fij = F->at(i, j);
139  T->at(i, j) = sqrt( ( i - CPy ) * ( i - CPy ) + ( j - CPx ) * ( j - CPx ) ) / Fij;
140  Frozen [ CPInd ] = true;
141
142  // For four direct neighbors or all eight neighbours of the nearest grid node, do the same
143  // (depending on parameter initDiag)
144
145  int nInd, ni, nj;
146  double time;
147  if ( initDiag <= 0. ) { // initialize direct neighbors only
148  for ( int neigh = 0; neigh < 4; neigh++ ) {
149  ni = i + iOffsets [ neigh ];
150  nj = j + jOffsets [ neigh ];
151  if ( isInDomain(ni, nj, m, n) ) {
152  nInd = ij2ind(ni, nj, m);
153  if ( centDiff > 0 ) { // use the average speed
154  time = sqrt( ( ni - CPy ) * ( ni - CPy ) + ( nj - CPx ) * ( nj - CPx ) ) / ( 0.5 * ( Fij + F->at(ni, nj) ) );
155  } else { // use the speed at the arrival point
156  time = sqrt( ( ni - CPy ) * ( ni - CPy ) + ( nj - CPx ) * ( nj - CPx ) ) / F->at(ni, nj);
157  }
158  if ( Frozen [ nInd ] ) {
159  T->at(ni, nj) = MIN( time, T->at(ni, nj) );
160  } else {
161  T->at(ni, nj) = time;
162  Frozen [ nInd ] = true;
163  }
164  }
165  }
166  } else { // initialize all neighbors
167  for ( int neigh = 0; neigh < 8; neigh++ ) {
168  ni = i + iOffsets_full [ neigh ];
169  nj = j + jOffsets_full [ neigh ];
170  if ( isInDomain(ni, nj, m, n) ) {
171  nInd = ij2ind(ni, nj, m);
172  if ( centDiff > 0 ) { // use the average speed
173  time = sqrt( ( ni - CPy ) * ( ni - CPy ) + ( nj - CPx ) * ( nj - CPx ) ) / ( 0.5 * ( Fij + F->at(ni, nj) ) );
174  } else { // use the speed at the arrival point
175  time = sqrt( ( ni - CPy ) * ( ni - CPy ) + ( nj - CPx ) * ( nj - CPx ) ) / F->at(ni, nj);
176  }
177  // for diagonal neighbors, use initDiag as a scaling factor
178  if ( is_diag [ neigh ] ) {
179  time *= initDiag;
180  }
181  if ( Frozen [ nInd ] ) {
182  T->at(ni, nj) = MIN( time, T->at(ni, nj) );
183  } else {
184  T->at(ni, nj) = time;
185  Frozen [ nInd ] = true;
186  }
187  }
188  }
189  }
190  }
191  // after each external change, the solution becomes invalid
192  solutionAvailable = false;
193 }
194
195 void
196 Grid :: setSolutionValueAt(int i, int j, double value)
197 {
198  T->at(i, j) = value;
199  // after each external change, the solution becomes invalid
200  solutionAvailable = false;
201 }
202
203 double
205 {
206  if ( !solutionAvailable ) {
207  int eFlag;
208  fastMarch(eFlag);
209  // TBD: one should check eFlag and send an error message if it is > 0
210  solutionAvailable = true;
211  }
212  return T->at(i, j);
213 }
214
215 void
217 {
218  for ( int j = 1; j <= n; j++ ) {
219  printf("\n");
220  for ( int i = 1; i <= m; i++ ) {
221  printf( "%d %d %g\n", j, i, T->at(i, j) );
222  }
223  }
224 }
225
227 void
228 Grid :: fastMarch(int &eFlag)
229 {
230  double time;
231  int i, j, ni, nj, nInd;
232  int tmpFlag = 0;
233  eFlag = 0;
234
235  // Create the initial narrow band
236  if ( !narrowBand ) {
237  narrowBand = new Heap(m *n);
238  }
239
240  // Loop over all grid points (not efficient, but done only once)
241  for ( int ind = 0; ind < ( m * n ); ind++ ) {
242  if ( Frozen [ ind ] ) {
243  i = ind2i(ind, m);
244  j = ind2j(ind, m);
245
246  for ( int neigh = 0; neigh < 4; neigh++ ) {
247  ni = i + iOffsets [ neigh ];
248  nj = j + jOffsets [ neigh ];
249  nInd = ij2ind(ni, nj, m);
250
251  if ( isInDomain(ni, nj, m, n) && !Frozen [ nInd ] ) {
252  if ( !narrowBand->isInHeap(nInd) ) {
253  time = calcTime(ni, nj, F->at(ni, nj), order, tmpFlag);
254  narrowBand->insert(time, nInd);
255  if ( tmpFlag > eFlag ) {
256  eFlag = tmpFlag;
257  }
258  }
259  }
260  }
261  }
262  }
263
264  // Loop
265  int CPInd;
266  while ( narrowBand->nElems() > 0 ) {
267  // Get minimal time-distance and index of this narrow-band element
268  time = narrowBand->getSmallest(& CPInd);
269  i = ind2i(CPInd, m);
270  j = ind2j(CPInd, m);
271
272  // Freeze and set time
273  Frozen [ CPInd ] = true;
274  T->at(i, j) = time;
275
276  // For all neighbours
277  for ( int neigh = 0; neigh < 4; neigh++ ) {
278  ni = i + iOffsets [ neigh ];
279  nj = j + jOffsets [ neigh ];
280  nInd = ij2ind(ni, nj, m);
281
282  // If valid for consideration
283  if ( isInDomain(ni, nj, m, n) && !Frozen [ nInd ] ) {
284  time = calcTime(ni, nj, F->at(ni, nj), order, tmpFlag);
285  // If T(ni,nj) has not been previously calculated
286  if ( !narrowBand->isInHeap(nInd) ) {
287  narrowBand->insert(time, nInd);
288  } else {
289  narrowBand->update(time, nInd);
290  }
291
292  if ( tmpFlag > eFlag ) {
293  eFlag = tmpFlag;
294  }
295  }
296  }
297  }
298  solutionAvailable = true;
299 }
300
301 // Time-dist calculation (for a grid with unit spacing)
302 double
303 Grid :: calcTime(int i, int j, double Fij, int ord, int &eFlag) {
304  // Get infinity.
305  // NB: Is this good practice? As in the matlab code, using Inf as
306  // a flag simplifies the code, so that's why I use this. Maybe this
307  // doesn't work with some compilers?
308  //double Inf = 1.0 / 0.0;
309  double Inf = std::numeric_limits<float>::infinity();
310
311  // Temporary error flag
312  int tmpFlag = 0;
313
314  // Frozen values at neighbors
315  double CrossVals [ 8 ];
316  //mj
317  double Fx = -1.;
318  double Fy = -1.;
319
320  // time calculated
321  double time;
322
323  // Indices of neighbouring nodes
324  int ni, nj, nInd;
325
326  // Variables used in the final formula
327  double xmin1 = 0., xmin2, ymin1 = 0., ymin2;
328
329  // Get values at surrounding nodes (at those that are frozen)
330  for ( int iter = 0; iter < 8; iter++ ) {
331  ni = i + icalcOffsets [ iter ];
332  nj = j + jcalcOffsets [ iter ];
333  nInd = ij2ind(ni, nj, m);
334  if ( isInDomain(ni, nj, m, n) && Frozen [ nInd ] ) {
335  CrossVals [ iter ] = T->at(ni, nj);
336  } else {
337  CrossVals [ iter ] = Inf;
338  }
339  }
340
341  // If none of the immediate neighbors is frozen, we cannot proceed
342  if ( ( CrossVals [ 1 ] == Inf ) && ( CrossVals [ 2 ] == Inf ) &&
343  ( CrossVals [ 5 ] == Inf ) && ( CrossVals [ 6 ] == Inf ) ) {
344  eFlag = 0;
345  time = Inf;
346  return time;
347  }
348
349  // Calculate coefficients of quadratic equation
350  double a = 0.;
351  double b = 0.;
352  // mj double c = -1. / ( Fij * Fij );
353  double c = 0.;
354
355  switch ( ord ) {
356  // First-order algorithm
357  case 1:
358  // Contribution of y-direction
359  if ( !( ( CrossVals [ 1 ] == Inf ) && ( CrossVals [ 2 ] == Inf ) ) ) {
360  ymin1 = MIN(CrossVals [ 1 ], CrossVals [ 2 ]);
361  a += 1.;
362  b -= 2. * ymin1;
363  c += ymin1 * ymin1;
364  //mj
365  if ( CrossVals [ 1 ] < CrossVals [ 2 ] ) {
366  ni = i + icalcOffsets [ 1 ];
367  nj = j + jcalcOffsets [ 1 ];
368  } else {
369  ni = i + icalcOffsets [ 2 ];
370  nj = j + jcalcOffsets [ 2 ];
371  }
372  Fy = F->at(ni, nj);
373  //end mj
374  }
375
376  // Contribution of x-direction
377  if ( !( ( CrossVals [ 5 ] == Inf ) && ( CrossVals [ 6 ] == Inf ) ) ) {
378  xmin1 = MIN(CrossVals [ 5 ], CrossVals [ 6 ]);
379  a += 1.;
380  b -= 2. * xmin1;
381  c += xmin1 * xmin1;
382  //mj
383  if ( CrossVals [ 5 ] < CrossVals [ 6 ] ) {
384  ni = i + icalcOffsets [ 5 ];
385  nj = j + jcalcOffsets [ 5 ];
386  } else {
387  ni = i + icalcOffsets [ 6 ];
388  nj = j + jcalcOffsets [ 6 ];
389  }
390  Fx = F->at(ni, nj);
391  //end mj
392  }
393  break;
394
395  // Second-order algorithm
396  case 2:
397  default:
398  // Contribution of y-direction
399  if ( !( ( CrossVals [ 1 ] == Inf ) && ( CrossVals [ 2 ] == Inf ) ) ) {
400  if ( CrossVals [ 1 ] < CrossVals [ 2 ] ) {
401  ymin1 = CrossVals [ 1 ];
402  //mj
403  ni = i + icalcOffsets [ 1 ];
404  nj = j + jcalcOffsets [ 1 ];
405  Fy = F->at(ni, nj);
406  //end mj
407  if ( CrossVals [ 0 ] <= CrossVals [ 1 ] ) { //second-order formula can be applied
408  ymin2 = CrossVals [ 0 ];
409  a += 2.25;
410  b += -6. * ymin1 + 1.5 * ymin2;
411  c += 4.0 * ymin1 * ymin1 + 0.25 * ymin2 * ymin2 - 2.0 * ymin1 * ymin2;
412  } else { //first-order formula
413  a += 1.;
414  b -= 2. * ymin1;
415  c += ymin1 * ymin1;
416  }
417  } else {
418  ymin1 = CrossVals [ 2 ];
419  //mj
420  ni = i + icalcOffsets [ 2 ];
421  nj = j + jcalcOffsets [ 2 ];
422  Fy = F->at(ni, nj);
423  //end mj
424  if ( CrossVals [ 3 ] <= CrossVals [ 2 ] ) { //second-order formula can be applied
425  ymin2 = CrossVals [ 3 ];
426  a += 2.25;
427  b += -6. * ymin1 + 1.5 * ymin2;
428  c += 4.0 * ymin1 * ymin1 + 0.25 * ymin2 * ymin2 - 2.0 * ymin1 * ymin2;
429  } else { //first-order formula
430  a += 1.;
431  b -= 2. * ymin1;
432  c += ymin1 * ymin1;
433  }
434  }
435  }
436  // Contribution of x-direction
437  if ( !( ( CrossVals [ 5 ] == Inf ) && ( CrossVals [ 6 ] == Inf ) ) ) {
438  if ( CrossVals [ 5 ] < CrossVals [ 6 ] ) {
439  xmin1 = CrossVals [ 5 ];
440  //mj
441  ni = i + icalcOffsets [ 5 ];
442  nj = j + jcalcOffsets [ 5 ];
443  Fx = F->at(ni, nj);
444  //end mj
445  if ( CrossVals [ 4 ] <= CrossVals [ 5 ] ) { //second-order formula can be applied
446  xmin2 = CrossVals [ 4 ];
447  a += 2.25;
448  b += -6. * xmin1 + 1.5 * xmin2;
449  c += 4.0 * xmin1 * xmin1 + 0.25 * xmin2 * xmin2 - 2.0 * xmin1 * xmin2;
450  } else { //first-order formula
451  a += 1.;
452  b -= 2. * xmin1;
453  c += xmin1 * xmin1;
454  }
455  } else {
456  xmin1 = CrossVals [ 6 ];
457  //mj
458  ni = i + icalcOffsets [ 6 ];
459  nj = j + jcalcOffsets [ 6 ];
460  Fx = F->at(ni, nj);
461  //end mj
462  if ( CrossVals [ 7 ] <= CrossVals [ 6 ] ) { //second-order formula can be applied
463  xmin2 = CrossVals [ 7 ];
464  a += 2.25;
465  b += -6. * xmin1 + 1.5 * xmin2;
466  c += 4.0 * xmin1 * xmin1 + 0.25 * xmin2 * xmin2 - 2.0 * xmin1 * xmin2;
467  } else { //first-order formula
468  a += 1.;
469  b -= 2. * xmin1;
470  c += xmin1 * xmin1;
471  }
472  }
473  }
474  }
475
476  // for centDiff=2, the average speed is used, otherwise the speed at arrival point is used
477  double Faver;
478  if ( centDiff != 2 ) {
479  Faver = Fij;
480  } else if ( Fx >= 0. && Fy >= 0. ) {
481  Faver = 0.5 * Fij + 0.25 * ( Fx + Fy );
482  } else if ( Fx >= 0. ) {
483  Faver = 0.5 * ( Fij + Fx );
484  } else if ( Fy >= 0. ) {
485  Faver = 0.5 * ( Fij + Fy );
486  } else {
487  Faver = Fij;
488  }
489
490  c -= 1. / ( Faver * Faver );
491
493  double d = ( b * b ) - ( 4.0 * a * c );
494
495  // Error treatment and appropriate time-dist calculation
496  if ( ( d < 0.0 ) && ( ord == 2 ) ) {
497  // if second-order method did not work, try first-order formula
498  time = calcTime(i, j, Fij, 1, tmpFlag);
499  eFlag = MAX(1, tmpFlag);
500  } else if ( ( d < 0.0 ) && ( ord == 1 ) ) {
501  if ( ( CrossVals [ 1 ] == Inf ) && ( CrossVals [ 2 ] == Inf ) ) {
502  ymin1 = Inf;
503  }
504  if ( ( CrossVals [ 5 ] == Inf ) && ( CrossVals [ 6 ] == Inf ) ) {
505  xmin1 = Inf;
506  }
507
508  eFlag = 2;
509  time = MIN(xmin1, ymin1) + 1.0 / Fij;
510  } else {
511  eFlag = 0;
512  time = ( -b + sqrt(d) ) / ( 2.0 * a );
513  }
514
515  return time;
516 }
517 }
int iOffsets[]
Definition: grid.C:28
void printSolutionAsData()
Definition: grid.C:216
#define MAX(a, b)
Definition: grid.C:9
int giveNumberOfColumns() const
Returns number of columns of receiver.
Definition: floatmatrix.h:158
int iOffsets_full[]
Definition: grid.C:15
double getSmallest(int *Ind)
Definition: heap.C:176
FloatMatrix * T
Matrix storing the values of the unknown (computed) field (e.g., arrival times)
Definition: grid.h:78
double giveSolutionValueAt(int i, int j)
Output methods.
Definition: grid.C:204
int ind2j(int ind, int m)
Definition: grid.h:98
void unFreeze()
Set all flags to "unfrozen".
Definition: grid.C:95
double calcTime(int i, int j, double Fij, int ord, int &eFlag)
Auxiliary function that evaluates the tentative value at one grid point, exploited by the fast marchi...
Definition: grid.C:303
int centDiff
Definition: grid.h:72
bool solutionAvailable
Flag indicating whether the solution has been computed.
Definition: grid.h:75
int jcalcOffsets[]
Definition: grid.C:41
void insert(double time, int Ind)
Definition: heap.C:146
int jOffsets[]
Definition: grid.C:31
void fastMarch(int &eFlag)
Fast marching method, solving the eikonal equation.
Definition: grid.C:228
#define MIN(a, b)
Definition: grid.C:8
int ij2ind(int i, int j, int m)
Utility methods.
Definition: grid.h:96
FloatMatrix * F
Matrix storing the values of the prescribed field (e.g., front velocities)
Definition: grid.h:81
#define OOFEM_ERROR(...)
Definition: error.h:61
void setZeroValues(FloatMatrix *gridCoords)
Methods setting the values of the unknown field at selected points (e.g., zero times) ...
Definition: grid.C:105
Grid(int N, int M)
Constructor (N = horizontal dimension, M = vertical dimension)
Definition: grid.C:46
int n
Grid dimensions: number of grid nodes horizontally (n) and vertically (m)
Definition: grid.h:67
Class implementing a heap, which is an auxiliary data structure used for efficient sorting and exploi...
Definition: heap.h:12
#define N(p, q)
Definition: mdm.C:367
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
int m
Definition: grid.h:67
void setToEmpty(int N)
Interface with external algorithms (such as fast marching)
Definition: heap.C:37
int order
Algorithmic parameters.
Definition: grid.h:70
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
void setSolutionValueAt(int i, int j, double value)
Definition: grid.C:196
int jOffsets_full[]
Definition: grid.C:18
bool isInHeap(int Ind)
Definition: heap.C:136
bool isInDomain(int i, int j, int m, int n)
Definition: grid.h:99
int giveSize(int dir)
Definition: grid.C:70
int nElems()
Definition: heap.C:142
void update(double time, int Ind)
Definition: heap.C:167
int icalcOffsets[]
Definition: grid.C:38
double initDiag
Definition: grid.h:71
int ind2i(int ind, int m)
Definition: grid.h:97
the oofem namespace is to define a context or scope in which all oofem names are defined.
Heap * narrowBand
Heap used for efficient sorting and detection of smallest candidate.
Definition: grid.h:87
~Grid()
Destructor.
Definition: grid.C:59
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
bool * Frozen
Array storing indicators of frozen points (at which the solution is already known) ...
Definition: grid.h:84
bool is_diag[]
Definition: grid.C:23