UltraScan III
us_astfem_math.cpp
Go to the documentation of this file.
1 #include "us_astfem_math.h"
3 #include "us_math2.h"
4 #include "us_hardware.h"
5 
6 void US_AstfemMath::interpolate_C0( MfemInitial& C0, double* C1,
7  QVector< double >& xvec )
8 {
9  int ja = 0;
10  double* x = xvec .data();
11  double* radi = C0.radius .data();
12  double* conc = C0.concentration.data();
13  int rsiz = C0.radius .size();
14 
15  for ( int j = 0; j < xvec.size(); j++ )
16  {
17  int i;
18  double xs = x[ j ];
19  double xc = xs + 1.e-12;
20 
21  for ( i = ja; i < rsiz; i++ )
22  if ( radi[ i ] > xc ) break;
23 
24 
25  if ( i == 0 ) // x[j] < C0.radius[0]
26  C1[ j ] = conc[ 0 ]; // Use the first value
27 
28  else if ( i == rsiz ) // x[j] > last point in C0.radius[]
29  C1[ j ] = conc[ i - 1 ];
30 
31  else
32  {
33  double a = radi[ i - 1 ];
34  double b = radi[ i ];
35 
36  double tmp = ( xs - a ) / ( b - a );
37 
38  ja = i - 1;
39 
40  C1[ j ] = conc[ ja ] * ( 1.0 - tmp ) +
41  conc[ i ] * tmp;
42 
43  }
44  }
45 }
46 
47 // Original grid: C0, final grid: C1
48 void US_AstfemMath::interpolate_C0( MfemInitial& C0, MfemInitial& C1 )
49 {
50  int ja = 0;
51  double* C0radi = C0.radius .data();
52  double* C1radi = C1.radius .data();
53  double* C0conc = C0.concentration.data();
54  double* C1conc = C1.concentration.data();
55  int r0size = C0.radius .size();
56  int r1size = C1.radius .size();
57 
58  for ( int j = 0; j < r1size; j++ )
59  {
60  int i;
61  double xs = C1radi[ j ];
62  double xc = xs + 1.e-12;
63 
64  for ( i = ja; i < r0size; i++ )
65  if ( C0.radius[ i ] > xc ) break;
66 
67 
68  if ( i == 0 ) // x[j] < C0.radius[0]
69  C1conc[ j ] = C0conc[ 0 ]; // Use the first value
70 
71  else if ( i == r0size ) // x[j] > last point in C0.radius[]
72  C1conc[ j ] = C0conc[ i - 1 ];
73 
74  else
75  {
76  ja = i - 1;
77 
78  double a = C0radi[ ja ];
79  double b = C0radi[ i ];
80 
81  double tmp = ( xs - a ) / ( b - a );
82 
83  C1conc[ j ] = C0conc[ ja ] * ( 1.0 - tmp ) + C0conc[ i ] * tmp;
84  }
85  }
86 }
87 
88 void US_AstfemMath::zero_2d( int val1, int val2, double** matrix )
89 {
90  for ( int i = 0; i < val1; i++ )
91  for ( int j = 0; j < val2; j++ )
92  matrix[ i ][ j ] = 0.0;
93 }
94 
95 void US_AstfemMath::initialize_2d( int val1, int val2, double*** matrix )
96 {
97  *matrix = new double* [ val1 ];
98 
99  for ( int i = 0; i < val1; i++ )
100  {
101  (*matrix)[ i ] = new double [ val2 ];
102 
103  for ( int j = 0; j < val2; j++ )
104  (*matrix)[ i ][ j ] = 0.0;
105 
106  }
107 }
108 
109 void US_AstfemMath::clear_2d( int val1, double** matrix )
110 {
111  for ( int i = 0; i < val1; i++ ) delete [] matrix[i];
112 
113  delete [] matrix;
114 }
115 
116 double US_AstfemMath::minval( const QVector< double >& value )
117 {
118  const double* avalue = value.data();
119  double minimum = 1.0e300;
120 
121  for ( int i = 0; i < value.size(); i++ )
122  minimum = qMin( minimum, avalue[ i ] );
123 
124  return minimum;
125 }
126 
127 double US_AstfemMath::minval( const QVector< US_Model::SimulationComponent >& value )
128 {
129  double minimum = 1.0e300;
130 
131  for ( int i = 0; i < value.size(); i++ )
132  minimum = min( minimum, value[ i ].s );
133 
134  return minimum;
135 }
136 
137 double US_AstfemMath::maxval( const QVector< double >& value )
138 {
139  const double* avalue = value.data();
140  double maximum = -1.0e300;
141 
142  for ( int i = 0; i < value.size(); i++ )
143  maximum = qMax( maximum, avalue[ i ]);
144 
145  return maximum;
146 }
147 
148 double US_AstfemMath::maxval( const QVector< US_Model::SimulationComponent >& value )
149 {
150  double maximum = -1.0e300;
151 
152  for ( int i = 0; i < value.size(); i++ )
153  maximum = max( maximum, value[ i ].s );
154 
155  return maximum;
156 }
157 
158 void US_AstfemMath::initialize_3d(
159  int val1, int val2, int val3, double**** matrix )
160 {
161  *matrix = new double** [ val1 ];
162 
163  for ( int i = 0; i < val1; i++ )
164  {
165  (*matrix)[ i ] = new double *[ val2 ];
166 
167  for ( int j = 0; j < val2; j++ )
168  {
169  (*matrix)[ i ][ j ] = new double [ val3 ];
170 
171  for ( int k = 0; k < val3; k++ )
172  {
173  (*matrix)[ i ][ j ][ k ] = 0.0;
174  }
175  }
176  }
177 }
178 
179 void US_AstfemMath::clear_3d( int val1, int val2, double*** matrix )
180 {
181  for ( int i = 0; i < val1; i++ )
182  {
183  for ( int j = 0; j < val2; j++ )
184  {
185  delete [] matrix[ i ][ j ];
186  }
187 
188  delete [] matrix[ i ];
189  }
190 
191  delete [] matrix;
192 }
193 
194 void US_AstfemMath::tridiag( double* a, double* b, double* c,
195  double* r, double* u, int N )
196 {
197  double bet = b[ 0 ];
198 #ifdef NO_DB
199  static QVector< double > gamvec;
200  static int Nsave = 0;
201 
202  if ( N > Nsave )
203  {
204  Nsave = N;
205  gamvec.resize( N );
206  }
207 #else
208  QVector< double > gamvec( N );
209 #endif
210  double* gam = gamvec.data();
211 
212  if ( bet == 0.0 ) qDebug() << "Error 1 in tridiag";
213 
214  u[ 0 ] = r[ 0 ] / bet;
215 
216  for ( int j = 1; j < N; j++ )
217  {
218  gam[ j ] = c[ j - 1 ] / bet;
219  bet = b[ j ] - a[ j ] * gam[ j ];
220 
221  if ( bet == 0.0 ) qDebug() << "Error 2 in tridiag";
222 
223  u[ j ] = ( r[ j ] - a[ j ] * u[ j - 1 ] ) / bet;
224  }
225 
226  gam[ 0 ] = gam[ 1 ];
227 
228  for ( int j = N - 2; j >= 0; j-- )
229  u[ j ] -= gam[ j + 1 ] * u[ j + 1 ];
230 }
231 
233 //
234 // cub_root: find the positive cubic-root of a cubic polynomial
235 // p(x)=a0+a1*x+a2*x^2+x^3
236 //
237 // with: a0<=0, and a1, a2>=0;
238 //
240 double US_AstfemMath::cube_root( double a0, double a1, double a2 )
241 {
242  double x;
243 
244  double Q = ( 3.0 * a1 - sq( a2 ) ) / 9.0;
245  double S = ( 9.0 * a1 * a2 - 27.0 * a0 - 2.0 * pow( a2, 3.0 ) ) / 54.0;
246  double D = pow( Q, 3.0 ) + sq( S );
247 
248  if ( D < 0 )
249  {
250  double theta = acos( S / sqrt( pow( -Q, 3.0 ) ) );
251  x = 2.0 * sqrt( -Q ) * cos( theta / 3.0) ;
252  }
253  else
254  {
255  double B;
256  double Dc = sqrt( D );
257 
258  if ( S + Dc < 0 )
259  B = - pow( - ( S + Dc ), 1.0 / 3.0 ) - pow( Dc - S, 1.0 / 3.0 );
260 
261  else if (S - Dc < 0)
262  B = pow( S + Dc, 1.0 / 3.0) - pow( Dc - S, 1.0 / 3.0);
263 
264  else
265  B = pow( S + Dc, 1.0 / 3.0 ) + pow( S - Dc, 1.0 / 3.0 );
266 
267  double Dc2 = -3.0 * (pow(B, 2.0) + 4 * Q);
268 
269  if ( Dc2 > 0 )
270  x = max( B, 0.5 * ( -B + sqrt( Dc2 ) ) );
271  else
272  x = B;
273  }
274 
275  x = x - a2 / 3.0;
276 
277  return x;
278 }
279 
281 //
282 // find_C1_mono_Nmer: find C1 from C1 + K * C1^n = CT
283 //
285 double US_AstfemMath::find_C1_mono_Nmer( int n, double K, double CT )
286 {
287  // use Newton's method for f(x) = x + K*x^n - CT
288  // x_{i+1} = x_i - f(x_i)/f'(x_i)
289 
290  double x1;
291  double x0 = 0.0;
292  const int MaxNumIt = 1000;
293  int i;
294 
295  if ( CT <= 0.0 ) return 0.0 ;
296  if ( CT <= 1.0e-12 ) return CT;
297 
298  for ( i = 1; i < MaxNumIt; i++ )
299  {
300  x1 = ( K * ( n - 1.0 ) * pow( x0, (double)n ) + CT ) /
301  ( 1.0 + K * n * pow( x0, n - 1.0 ) );
302 
303  if ( fabs( x1 - x0 ) / ( 1.0 + fabs( x1 ) ) < 1.e-12 ) break;
304  x0 = x1;
305  }
306 
307  if ( i == MaxNumIt )
308  {
309  qDebug() << "warning: Newton's method did not coonverges "
310  "in find_C1_mono_Nmer";
311  return -1.0;
312  }
313 
314  return 0.5 * ( x0 + x1 );
315 }
316 
318 //
319 // Gass Elimination for n X n system: Ax=b
320 //
321 // return value: -1: A singular, no solution,
322 // 1: successful
323 // in return: A has been altered, and b stores the solution x
324 //
326 int US_AstfemMath::GaussElim( int n, double** a, double* b )
327 {
328  // Elimination
329  for ( int i = 0; i < n; i++ )
330  {
331  // find the pivot
332  double amax = fabs( a[ i ][ i ] );
333  int Imx = i;
334 
335  for ( int ip = i + 1; ip < n; ip++ )
336  {
337  if ( fabs( a[ ip ][ i ] ) > amax )
338  {
339  amax = fabs( a[ ip ][ i ]);
340  Imx = ip;
341  }
342  }
343 
344  if ( amax == 0 )
345  {
346  qDebug() << "Singular matrix in routine GaussElim";
347  return -1;
348  }
349 
350  double tmp;
351 
352  // interchange i-th and Imx-th row
353  if ( Imx != i )
354  {
355  double* ptmp = a[ i ];
356  a[ i ] = a[ Imx ];
357  a[ Imx ] = ptmp;
358 
359  tmp = b[ i ];
360  b[ i ] = b[ Imx ];
361  b[ Imx ] = tmp;
362  }
363 
364  // Elimination
365  tmp = a[ i ][ i ];
366 
367  for ( int j = i; j < n; j++ ) a[ i ][ j ] /= tmp;
368 
369  b[ i ] /= tmp;
370 
371  for ( int ip = i + 1; ip < n; ip++ )
372  {
373  tmp = a[ ip ][ i ];
374 
375  for ( int j = i + 1; j < n; j++ )
376  a[ ip ][ j ] -= tmp * a[ i ][ j ];
377 
378  b[ ip ] -= tmp * b[ i ];
379  }
380  }
381 
382  // Backward substitution
383  for ( int i = n - 2; i >= 0; i-- )
384  for ( int j = i + 1; j < n; j++ ) b[ i ] -= a[ i ][ j ] * b[ j ];
385 
386  return 1;
387 }
388 
389 // Interpolation routine By B. Demeler 041708
390 int US_AstfemMath::interpolate( MfemData& expdata, MfemData& simdata,
391  bool use_time )
392 {
393  // NOTE: *expdata has to be initialized to have the proper size
394  // (filled with zeros) before using this routine!
395  // The radius also has to be assigned!
396 
397  if ( expdata.scan.size() == 0 || expdata.scan[ 0 ].conc.size() == 0 ||
398  simdata.scan.size() == 0 || simdata.radius.size() == 0 )
399  return -1;
400 
401  // Iterate through all experimental data scans and find the first time point
402  // in simdata that is higher or equal to each time point in expdata:
403 
404  int fscan = -1;
405  int lscan = -1;
406  int escans = expdata.scan.size();
407  int sscans = simdata.scan.size();
408  double e_omega = 0.0;
409  double e_time = 0.0;
410 
411  if ( use_time )
412  {
413  double s_time = simdata.scan[ 0 ].time;
414  double l_time = simdata.scan[ sscans - 1 ].time;
415 
416  for ( int expscan = 0; expscan < escans; expscan++ )
417  { // First determine the scan range of current experiment data
418  e_time = expdata.scan[ expscan ].time;
419  if ( fscan < 0 && e_time >= s_time )
420  fscan = expscan;
421  if ( e_time > l_time )
422  {
423  lscan = expscan;
424  break;
425  }
426  }
427 
428 qDebug() << "MATHi: s_time e_time" << s_time << e_time
429  << "fscan lscan" << fscan << lscan;
430  fscan = ( fscan < 0 ) ? 0 : fscan;
431  lscan = ( lscan < 0 ) ? escans : qMin( lscan, escans );
432  }
433  else // Use omega^2t integral for interpolation
434  {
435  double s_omega = simdata.scan[ 0 ].omega_s_t;
436  double l_omega = simdata.scan[ sscans - 1 ].omega_s_t;
437 
438  for ( int expscan = 0; expscan < escans; expscan++ )
439  { // First determine the scan range of current experiment data
440  e_omega = expdata.scan[ expscan ].omega_s_t;
441  if ( fscan < 0 && e_omega >= s_omega )
442  fscan = expscan;
443  if ( e_omega > l_omega )
444  {
445  lscan = expscan;
446  break;
447  }
448  }
449 
450  fscan = ( fscan < 0 ) ? 0 : fscan;
451  lscan = ( lscan < 0 ) ? escans : qMin( lscan, escans );
452 //qDebug() << "MATHi: s_omega l_omega" << s_omega << e_omega
453 // << "fscan lscan" << fscan << lscan;
454  }
455 
456  if ( fscan < 0 || lscan < 0 )
457  return -1;
458 
459  // Interpolate all radial points from each scan in tmp_data onto expdata
460  interpolate( expdata, simdata, use_time, fscan, lscan );
461 
462  return 0;
463 }
464 
465 // Interpolation routine By B. Demeler 041708
466 int US_AstfemMath::interpolate( MfemData& expdata, MfemData& simdata,
467  bool use_time, int fscan, int lscan )
468 {
469 //*TIMING*
470 //static int ncall=0;
471 //static int totMs=0;
472 //QDateTime tStart = QDateTime::currentDateTime();
473 //*TIMING*
474  // NOTE: *expdata has to be initialized to have the proper size
475  // (filled with zeros) before using this routine!
476  // The radius also has to be assigned!
477  int nsscan = simdata.scan.size();
478  int nsconc = simdata.radius.size();
479  int nescan = expdata.scan.size();
480  int neconc = expdata.scan[ 0 ].conc.size();
481 
482  if ( nescan == 0 || neconc == 0 || nsscan == 0 || nsconc == 0 )
483  return -1;
484 
485  // First, create a temporary MfemData instance (tmp_data) that has the
486  // same radial grid as simdata, but the same time grid as the experimental
487  // data. The time and w2t integral values are interpolated for the tmp_data
488  // structure.
489 
490  MfemData tmp_data;
491  MfemScan tmp_scan;
492 
493  tmp_data.scan .clear();
494  tmp_data.radius.clear();
495  tmp_data.scan .reserve( nescan );
496  tmp_data.radius.reserve( nsconc );
497 
498  // Fill tmp_data.radius with radius positions from the simdata array:
499 
500  for ( int ii = 0; ii < nsconc; ii++ )
501  {
502  tmp_data.radius << simdata.radius[ ii ];
503  }
504 
505  // Iterate through all experimental data scans and find the first time point
506  // in simdata that is higher or equal to each time point in expdata:
507 
508  int simscan = 0;
509  double e_omega = 0.0;
510  double e_time = 0.0;
511  double s_omega1;
512  double s_omega2;
513  double s_time1;
514  double s_time2;
515 
516  if ( use_time )
517  {
518  int eftime = (int)expdata.scan[ 0 ].time;
519  int sltime = (int)simdata.scan[ nsscan - 1 ].time;
520 
521  for ( int expscan = fscan; expscan < lscan; expscan++ )
522  { // Interpolate where needed to a range of experiment scans
523  MfemScan* sscan1 = &simdata.scan[ simscan ];
524  MfemScan* sscan2 = sscan1;
525  MfemScan* escan = &expdata.scan[ expscan ];
526 
527  e_omega = escan->omega_s_t;
528  e_time = escan->time;
529 
530  while ( simdata.scan[ simscan ].time < e_time )
531  {
532  simscan ++;
533  // Make sure we don't overrun bounds:
534  if ( simscan == nsscan )
535  { // Sim scan count has exceeded limit
536  if ( sltime >= eftime )
537  { // Output a message if time ranges overlap
538  qDebug() << "simulation time scan[" << simscan << "]: "
539  << simdata.scan[ simscan - 1 ].time
540  << ", expdata scan time[" << expscan << "]: "
541  << expdata.scan[ expscan ].time;
542 
543  qDebug() << "The simulated data does not cover the entire "
544  "experimental time range and ends too early!\n"
545  "exiting...\n";
546  }
547 
548 #if 0
549 #if defined(USE_MPI)
550  MPI_Abort( MPI_COMM_WORLD, -1 );
551 #endif
552  exit( -1 );
553 #endif
554 #if 1
555  break;
556 #endif
557  }
558  }
559 
560  int sscm = ( simscan > 0 ) ? ( simscan - 1 ) : 1;
561  sscan1 = &simdata.scan[ sscm ];
562  sscan2 = &simdata.scan[ simscan ];
563  s_omega1 = sscan1->omega_s_t;
564  s_time1 = sscan1->time;
565  s_omega2 = sscan2->omega_s_t;
566  s_time2 = sscan2->time;
567 
568  // Check to see if the time is equal or larger:
569  if ( s_time2 == e_time )
570  { // they are the same, so take this scan
571  // and push it onto the tmp_data array.
572  tmp_data.scan << simdata.scan[ simscan ];
573  }
574  else // interpolation is needed
575  {
576  double a;
577  double b;
578  tmp_scan.conc.clear();
579  tmp_scan.conc.reserve( nsconc );
580 
581  // interpolate the concentration points:
582  for ( int ii = 0; ii < nsconc; ii++ )
583  {
584  a = ( sscan2->conc[ ii ] - sscan1->conc[ ii ] ) /
585  ( s_time2 - s_time1 );
586 
587  b = sscan2->conc[ ii ] - a * s_time2;
588 
589  tmp_scan.conc << ( a * e_time + b );
590  }
591 
592  // interpolate the omega_square_t integral data:
593 
594  a = ( s_omega2 - s_omega1 ) / ( s_time2 - s_time1 );
595 
596  b = s_omega2 - a * s_time2;
597 
598  escan->omega_s_t = a * e_time + b;
599 
600  tmp_data.scan << tmp_scan;
601  }
602  }
603  }
604  else // Use omega^2t integral for interpolation
605  {
606  for ( int expscan = fscan; expscan < lscan; expscan++ )
607  { // Interpolate where needed to a range of experiment scans
608  MfemScan* sscan1 = &simdata.scan[ simscan ];
609  MfemScan* sscan2 = sscan1;
610  MfemScan* escan = &expdata.scan[ expscan ];
611 
612  e_omega = escan->omega_s_t;
613  e_time = escan->time;
614 //qDebug() << "MATHi: somg eomg" << simdata.scan[simscan].omega_s_t << e_omega
615 // << " escn sscn" << expscan << simscan;
616 
617  while ( simdata.scan[ simscan ].omega_s_t < e_omega )
618  {
619  simscan ++;
620  // Make sure we don't overrun bounds:
621  if ( simscan == (int) simdata.scan.size() )
622  {
623  qDebug() << "simulation omega^2t scan[" << simscan << "]: "
624  << simdata.scan[ simscan - 1 ].omega_s_t
625  << ", expdata scan omega^2t[" << expscan << "]: "
626  << expdata.scan[ expscan ].omega_s_t;
627 //qDebug() << "e_omega e_time" << e_omega << e_time;
628 //int lsc=expdata.scan.size()-1;
629 //e_omega=expdata.scan[lsc].omega_s_t;
630 //e_time=expdata.scan[lsc].time;
631 //qDebug() << "e-lsc e_omega e_time" << lsc << e_omega << e_time;
632 //lsc=simdata.scan.size()-1;
633 //e_omega=simdata.scan[lsc].omega_s_t;
634 //e_time=simdata.scan[lsc].time;
635 //qDebug() << "s-lsc e_omega e_time" << lsc << e_omega << e_time;
636 #if defined(USE_MPI)
637 // MPI_Abort( MPI_COMM_WORLD, -1 );
638 #endif
639  qDebug() << "The simulated data does not cover the entire "
640  "experimental time range and ends too early!\n"
641  "exiting...";
642 // exit( -1 );
643  simscan--;
644  break;
645  }
646  }
647 
648  int sscm = ( simscan > 0 ) ? ( simscan - 1 ) : 1;
649  sscan1 = &simdata.scan[ sscm ];
650  sscan2 = &simdata.scan[ simscan ];
651  s_omega1 = sscan1->omega_s_t;
652  s_time1 = sscan1->time;
653  s_omega2 = sscan2->omega_s_t;
654  s_time2 = sscan2->time;
655 
656  // Check to see if the time is equal or larger:
657  if ( s_omega2 == e_omega )
658  { // They are the same, so take this scan and
659  // push it onto the tmp_data array.
660  tmp_data.scan << simdata.scan[ simscan ];
661  }
662  else // Interpolation is needed
663  {
664  double a;
665  double b;
666  tmp_scan.conc.clear();
667  tmp_scan.conc.reserve( nsconc );
668 
669  // Interpolate the concentration points:
670  for ( int ii = 0; ii < nsconc; ii++ )
671  {
672  a = ( sscan2->conc[ ii ] - sscan1->conc[ ii ] ) /
673  ( s_omega2 - s_omega1 );
674 
675  b = sscan2->conc[ ii ] - a * s_omega2;
676 
677  tmp_scan.conc << ( a * e_omega + b );
678  }
679 
680 #if 0
681  // Interpolate the omega_square_t integral data:
682  a = ( s_time2 - s_time1 ) / ( s_omega2 - s_omega1 );
683 
684  b = s_time2 - a * s_omega2;
685 
686  escan->time = a * e_omega + b;
687 #endif
688 
689  tmp_data.scan << tmp_scan;
690  }
691  }
692  }
693 
694  // Interpolate all radial points from each scan in tmp_data onto expdata
695 
696  if ( tmp_data.radius[ 0 ] > expdata.radius[ 0 ] )
697  {
698  qDebug() << "Radius comparison: " << tmp_data.radius[ 0 ]
699  << " (simulated), " << expdata.radius[ 0 ]
700  << " (experimental)";
701 
702  qDebug() << "jj = " << 0 << ", simdata radius: "
703  << tmp_data.radius[ 0 ] << ", expdata radius: "
704  << expdata.radius[ 0 ]; // Changed from radius[ ii ]
705  // ii out of scope
706 
707  qDebug() << "The simulated data radial range does not include the "
708  "beginning of the experimental data's radii!\n"
709  "exiting...";
710 
711 #if defined(USE_MPI)
712  MPI_Abort( MPI_COMM_WORLD, -3 );
713 #endif
714  exit( -3 );
715  }
716 
717  int jj = 0;
718 
719  for ( int ii = 0; ii < neconc; ii++ )
720  {
721  while ( tmp_data.radius[ jj ] < expdata.radius[ ii ] )
722  {
723  jj++;
724  // make sure we don't overrun bounds:
725  if ( jj == tmp_data.radius.size() )
726  {
727  qDebug() << "The simulated data does not have enough "
728  "radial points and ends too early!\n"
729  "exiting...";
730 //qDebug() << "jj ii szerad trad erad" << jj << ii << expdata.radius.size()
731 // << tmp_data.radius[jj-1] << expdata.radius[ii];
732 #if defined(USE_MPI)
733  MPI_Abort( MPI_COMM_WORLD, -2 );
734 #endif
735  exit( -2 );
736  }
737  }
738 
739  // check to see if the radius is equal or larger:
740  if ( tmp_data.radius[ jj ] == expdata.radius[ ii ] )
741  { // they are the same, so simply update the concentration value:
742  int ee = 0;
743  for ( int expscan = fscan; expscan < lscan; expscan++ )
744  expdata.scan[ expscan ].conc[ ii ] +=
745  tmp_data.scan[ ee++ ].conc[ jj ];
746  }
747  else // interpolation is needed
748  {
749  int ee = 0;
750  int mm = jj - 1;
751  double radius1 = tmp_data.radius[ mm ];
752  double radius2 = tmp_data.radius[ jj ];
753  double eradius = expdata .radius[ ii ];
754  double radrange = radius2 - radius1;
755 //qDebug() << "MATHi: esize jj" << expdata.scan.size() << jj;
756 //qDebug() << "MATHi: r1 r2 er rr" << radius1 << radius2 << eradius << radrange
757 // << "fscan lscan" << fscan << lscan;
758 
759  for ( int expscan = fscan; expscan < lscan; expscan++ )
760  {
761  MfemScan* tscan = &tmp_data.scan[ ee++ ];
762 
763  double a = ( tscan->conc[ jj ] - tscan->conc[ mm ] ) / radrange;
764 
765  double b = tscan->conc[ jj ] - a * radius2;
766 
767  expdata.scan[ expscan ].conc[ ii ] += ( a * eradius + b );
768  }
769  }
770  }
771 //*TIMING*
772 //QDateTime tEnd=QDateTime::currentDateTime();
773 //ncall++;
774 //totMs+=tStart.msecsTo(tEnd);
775 //if ((ncall%100)==1) qDebug() << "interpolate() calls time-ms" << ncall << totMs;
776 //*TIMING*
777  return 0;
778 }
779 
780 void US_AstfemMath::QuadSolver( double* ai, double* bi, double* ci,
781  double* di, double* cr, double* solu, int N )
782 {
783 // Solve Quad-diagonal system [a_i, *b_i*, c_i, d_i]*[x]=[r_i]
784 // b_i are on the main diagonal line
785 //
786 // Test
787 // n=100; a=-1+rand(100,1); b=2+rand(200,1); c=-0.7*rand(100,1); d=-0.3*rand(100,1);
788 // xs=rand(100,1);
789 // r(1)=b(1)*xs(1)+c(1)*xs(2)+d(1)*xs(3);
790 // for i=2:n-2,
791 // r(i)=a(i)*xs(i-1)+b(i)*xs(i)+c(i)*xs(i+1)+d(i)*xs(i+2);
792 // end;
793 // i=n-1; r(i)=a(i)*xs(i-1)+b(i)*xs(i)+c(i)*xs(i+1);
794 // i=n; r(i)=a(i)*xs(i-1)+b(i)*xs(i);
795 
796  QVector< double > caVec( N );
797  QVector< double > cbVec( N );
798  QVector< double > ccVec( N );
799  QVector< double > cdVec( N );
800  double* ca = caVec.data();
801  double* cb = cbVec.data();
802  double* cc = ccVec.data();
803  double* cd = cdVec.data();
804 
805  for ( int i = 0; i < N; i++ )
806  {
807  ca[ i ] = ai[ i ];
808  cb[ i ] = bi[ i ];
809  cc[ i ] = ci[ i ];
810  cd[ i ] = di[ i ];
811  }
812 
813  for ( int i = 1; i <= N - 2; i++ )
814  {
815  double tmp = ca[ i ] / cb[ i - 1 ];
816 
817  cb[ i ] = cb[ i ] -cc [ i - 1] * tmp;
818  cc[ i ] = cc[ i ] -cd [ i - 1] * tmp;
819  cr[ i ] = cr[ i ] -cr [ i - 1] * tmp;
820  }
821 
822  int i = N - 1;
823  double tmp = ca[ i ] / cb[ i - 1 ];
824 
825  cb[ i] = cb[ i ] - cc[ i - 1 ] * tmp;
826  cr[ i ] = cr[ i ] - cr[ i - 1 ] * tmp;
827 
828  solu[ N - 1 ] = cr[ N - 1 ] / cb[ N - 1 ];
829  solu[ N - 2 ] = ( cr[ N - 2 ] - cc[ N - 2 ] * solu[ N - 1] ) / cb[ N - 2 ];
830 
831  i = N - 2;
832  do
833  {
834  i--;
835  solu[ i ] = ( cr[ i ]
836  - cc[ i ] * solu[ i + 1 ]
837  - cd[ i ] * solu[ i + 2 ] ) /
838  cb[ i ];
839 
840  } while ( i != 0 );
841 }
842 
843 // old version: perform integration on supp(test function) separately
844 // on left Q and right T
845 
846 void US_AstfemMath::IntQT1( double* vx, double D, double sw2,
847  double** Stif, double dt )
848 {
849  // element to define basis functions
850 
851  int npts, i, k;
852  double x_gauss, y_gauss, dval;
853  double hh, slope, xn1, phiC, phiCx;
854  double Lx[ 3 ], Ly[ 3 ];
855  double Rx[ 4 ];
856  double Qx[ 3 ], Qy[ 3 ];
857  double Tx[ 3 ], Ty[ 3 ];
858  double** StifL;
859  double** StifR;
860  double** Lam;
861  double DJac;
862  double phiL [ 3 ];
863  double phiLx[ 3 ];
864  double phiLy[ 3 ];
865  double phiR [ 4 ];
866  double phiRx[ 4 ];
867  double phiRy[ 4 ];
868 
869  // elements for define the trial function phi
870  Lx[ 0 ] = vx[0]; // vertices of left Triangle
871  Lx[ 1 ] = vx[3];
872  Lx[ 2 ] = vx[2];
873 
874  Ly[ 0 ] = 0.0;
875  Ly[ 1 ] = dt;
876  Ly[ 2 ] = dt;
877 
878  Rx[ 0 ] = vx[0]; // vertices of Q on right quadrilateral
879  Rx[ 1 ] = vx[1];
880  Rx[ 2 ] = vx[4];
881  Rx[ 3 ] = vx[3];
882 
883  initialize_2d( 3, 2, &StifL );
884  initialize_2d( 4, 2, &StifR );
885 
886  hh = vx[ 3 ] - vx[ 2 ];
887  slope = ( vx[ 3 ] - vx[ 5 ] ) / dt;
888  npts = 28;
889  initialize_2d( npts, 4, &Lam );
890  DefineFkp( npts, (double**)Lam );
891 
892  //
893  // integration over element Q (a triangle):
894  //
895  Qx[ 0 ] = vx[0]; // vertices of Q on left
896  Qx[ 1 ] = vx[3];
897  Qx[ 2 ] = vx[2];
898 
899  Qy[ 0 ] = 0.0;
900  Qy[ 1 ] = dt;
901  Qy[ 2 ] = dt;
902 
903  for (k=0; k<npts; k++)
904  {
905  x_gauss = Lam[k][0] * Qx[0] + Lam[k][1] * Qx[1] + Lam[k][2] * Qx[2];
906  y_gauss = Lam[k][0] * Qy[0] + Lam[k][1] * Qy[1] + Lam[k][2] * Qy[2];
907  DJac = 2.0 * AreaT( Qx, Qy );
908 
909  // trace-forward point at t_n+1 from (x_g, y_g)
910  xn1 = x_gauss + slope * ( dt - y_gauss );
911 
912  //
913  // find phi, phi_x, phi_y on L and C at (x,y)
914  //
915 
916  BasisTR( Lx, Ly, x_gauss, y_gauss, phiL, phiLx, phiLy);
917  // hat function on t_n+1, =1 at vx[3]; =0 at vx[2]
918  phiC = ( xn1 - vx[ 2 ] ) / hh;
919  phiCx = 1. / hh;
920 
921  for (i=0; i<3; i++)
922  {
923  dval = Integrand(x_gauss, D, sw2, phiL[i], phiLx[i], phiLy[i],
924  1.-phiC, -phiCx );
925  StifL[i][0] += Lam[k][3] * DJac * dval;
926 
927  dval = Integrand(x_gauss, D, sw2, phiL[i], phiLx[i], phiLy[i],
928  phiC, phiCx );
929  StifL[i][1] += Lam[k][3] * DJac * dval;
930  }
931  }
932 
933  //
934  // integration over T:
935  //
936  Tx[ 0 ] = vx[0]; // vertices of T on right
937  Tx[ 1 ] = vx[5];
938  Tx[ 2 ] = vx[3];
939 
940  Ty[ 0 ] = 0.0;
941  Ty[ 1 ] = 0.0;
942  Ty[ 2 ] = dt;
943 
944  for ( k = 0; k < npts; k++ )
945  {
946  x_gauss = Lam[k][0] * Tx[0] + Lam[k][1] * Tx[1] + Lam[k][2] * Tx[2];
947  y_gauss = Lam[k][0] * Ty[0] + Lam[k][1] * Ty[1] + Lam[k][2] * Ty[2];
948  DJac = 2.0 * AreaT( Tx, Ty );
949 
950  if( DJac < 1.e-22 ) break;
951 
952  xn1 = x_gauss + slope * ( dt - y_gauss ); // trace-forward point
953  // at t_n+1 from (x_g, y_g)
954  //
955  // find phi, phi_x, phi_y on R and C at (x,y)
956  //
957 
958  BasisQR( Rx, x_gauss, y_gauss, phiR, phiRx, phiRy, dt);
959  phiC = ( xn1 - vx[ 2 ] ) / hh; // hat function on t_n+1, =1 at vx[3];
960  // =0 at vx[2]
961  phiCx = 1. / hh;
962 
963  for (i = 0; i < 4; i++ )
964  {
965  dval = Integrand( x_gauss, D, sw2, phiR[ i ], phiRx[ i ],
966  phiRy[ i ], 1. - phiC, -phiCx );
967  StifR[ i ][ 0 ] += Lam[ k ][ 3 ] * DJac * dval;
968 
969  dval = Integrand( x_gauss, D, sw2, phiR[ i ], phiRx[ i ],
970  phiRy[ i ], phiC, phiCx );
971  StifR[ i ][ 1 ] += Lam[k][3] * DJac * dval;
972  }
973  }
974 
975  clear_2d( npts, Lam );
976 
977  for ( i = 0; i < 2; i++ )
978  {
979  Stif[ 0 ][ i ] = StifL[ 0 ][ i ] + StifR[ 0 ][ i ];
980  Stif[ 1 ][ i ] = StifR[ 1 ][ i ];
981  Stif[ 2 ][ i ] = StifL[ 2 ][ i ];
982  Stif[ 3 ][ i ] = StifL[ 1 ][ i ] + StifR[ 3 ][ i ];
983  Stif[ 4 ][ i ] = StifR[ 2 ][ i ];
984  }
985 
986  clear_2d( 3, StifL );
987  clear_2d( 4, StifR );
988 }
989 
990 void US_AstfemMath::IntQTm( double* vx, double D, double sw2,
991  double** Stif, double dt )
992 {
993  // element to define basis functions
994  //
995  int npts, i, k;
996  double x_gauss, y_gauss, dval;
997  double Lx[ 4 ];
998  double Cx[ 4 ];
999  double Rx[ 4 ];
1000  double Qx[ 4 ], Qy[ 4 ];
1001  double Tx[ 3 ], Ty[ 3 ];
1002  double DJac;
1003  double** StifL;
1004  double** StifR;
1005  double** Lam;
1006  double** Gs;
1007  double phiL [ 4 ];
1008  double phiLx[ 4 ];
1009  double phiLy[ 4 ];
1010  double phiCx[ 4 ];
1011  double phiCy[ 4 ];
1012  double phiC [ 4 ];
1013  double phiR [ 4 ];
1014  double phiRx[ 4 ];
1015  double phiRy[ 4 ];
1016 
1017  Lx[ 0 ] = vx[ 0 ];
1018  Lx[ 1 ] = vx[ 1 ];
1019  Lx[ 2 ] = vx[ 4 ];
1020  Lx[ 3 ] = vx[ 3 ]; // vertices of left T
1021 
1022  Cx[ 0 ] = vx[ 6 ];
1023  Cx[ 1 ] = vx[ 7 ];
1024  Cx[ 2 ] = vx[ 4 ];
1025  Cx[ 3 ] = vx[ 3 ];
1026 
1027  Rx[ 0 ] = vx[ 1 ]; // vertices of Q on right
1028  Rx[ 1 ] = vx[ 2 ];
1029  Rx[ 2 ] = vx[ 5 ];
1030  Rx[ 3 ] = vx[ 4 ];
1031 
1032  initialize_2d( 4, 2, &StifL );
1033  initialize_2d( 4, 2, &StifR );
1034 
1035  //
1036  // integration over element Q :
1037  //
1038  Qx[ 0 ] = vx[ 6 ]; // vertices of Q on right
1039  Qx[ 1 ] = vx[ 1 ];
1040  Qx[ 2 ] = vx[ 4 ];
1041  Qx[ 3 ] = vx[ 3 ];
1042 
1043  Qy[ 0 ] = 0.0;
1044  Qy[ 1 ] = 0.0;
1045  Qy[ 2 ] = dt;
1046  Qy[ 3 ] = dt; // vertices of left T
1047 
1048  npts = 5 * 5;
1049  initialize_2d( 25, 3, &Gs );
1050  DefineGaussian( 5, (double**)Gs );
1051 
1052  double psi[ 4 ];
1053  double psi1[ 4 ];
1054  double psi2[ 4 ];
1055  double jac[ 4 ];
1056 
1057  for ( k = 0; k < npts; k++ )
1058  {
1059  BasisQS( Gs[k][0], Gs[k][1], psi, psi1, psi2 );
1060 
1061  x_gauss = 0.0;
1062  y_gauss = 0.0;
1063  for ( i = 0; i < 4; i++ )
1064  {
1065  jac[ i ] = 0.0;
1066  }
1067  for ( i = 0; i < 4; i++ )
1068  {
1069  x_gauss += psi[ i ] * Qx[ i ];
1070  y_gauss += psi[ i ] * Qy[ i ];
1071  jac[ 0 ] += Qx[ i ] * psi1[ i ];
1072  jac[ 1 ] += Qx[ i ] * psi2[ i ];
1073  jac[ 2 ] += Qy[ i ] * psi1[ i ];
1074  jac[ 3 ] += Qy[ i ] * psi2[ i ];
1075  }
1076 
1077  DJac = jac[ 0 ] * jac[ 3 ] - jac[ 1 ] * jac[ 2 ];
1078 
1079  //
1080  // find phi, phi_x, phi_y on L and C at (x,y)
1081  //
1082 
1083  BasisQR( Lx, x_gauss, y_gauss, phiL, phiLx, phiLy, dt );
1084  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
1085 
1086  for ( i = 0; i < 4; i++ )
1087  {
1088  dval = Integrand( x_gauss, D, sw2, phiL[ i ], phiLx[ i ], phiLy[ i ],
1089  phiC[ 0 ] + phiC[ 3 ], phiCx[ 0 ] + phiCx[ 3 ] );
1090  StifL[ i ][ 0 ] += Gs[ k ][ 2 ] * DJac * dval;
1091 
1092  dval = Integrand( x_gauss, D, sw2, phiL[i], phiLx[i], phiLy[i],
1093  phiC[ 1 ] + phiC[ 2 ], phiCx[ 1 ] + phiCx[ 2 ] );
1094  StifL[ i ][ 1 ] += Gs[ k ][ 2 ] * DJac * dval;
1095  }
1096  }
1097 
1098  clear_2d( npts, Gs );
1099 
1100  //
1101  // integration over T:
1102  //
1103  Tx[ 0 ] = vx[ 1 ]; // vertices of T on right
1104  Tx[ 1 ] = vx[ 7 ];
1105  Tx[ 2 ] = vx[ 4 ];
1106 
1107  Ty[ 0 ] = 0.0;
1108  Ty[ 1 ] = 0.0;
1109  Ty[ 2 ] = dt;
1110 
1111  npts = 28;
1112  initialize_2d( npts, 4, &Lam );
1113  DefineFkp( npts, (double**)Lam );
1114 
1115  for ( k = 0; k < npts; k++ )
1116  {
1117  x_gauss = Lam[ k ][ 0 ] * Tx[ 0 ] + Lam[ k ][ 1 ] * Tx[ 1 ] +
1118  Lam[ k ][ 2 ] * Tx[ 2 ];
1119  y_gauss = Lam[ k ][ 0 ] * Ty[ 0 ] + Lam[ k ][ 1 ] * Ty[ 1 ] +
1120  Lam[ k ][ 2 ] * Ty[ 2 ];
1121  DJac = 2.0 * AreaT( Tx, Ty );
1122 
1123  if ( DJac < 1.e-22 ) break;
1124 
1125  //
1126  // find phi, phi_x, phi_y on R and C at (x,y)
1127  //
1128 
1129  BasisQR( Rx, x_gauss, y_gauss, phiR, phiRx, phiRy, dt );
1130  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
1131 
1132  for (i = 0; i < 4; i++ )
1133  {
1134  dval = Integrand( x_gauss, D, sw2, phiR[ i ], phiRx[ i ], phiRy[ i ],
1135  phiC[ 0 ] + phiC[ 3 ], phiCx[ 0 ] + phiCx[ 3 ] );
1136  StifR[ i ][ 0 ] += Lam[ k ][ 3 ] * DJac * dval;
1137 
1138  dval = Integrand( x_gauss, D, sw2, phiR[ i ], phiRx[ i ], phiRy[ i ],
1139  phiC[ 1 ] + phiC[ 2 ], phiCx[ 1 ] + phiCx[ 2 ] );
1140  StifR[ i ][ 1 ] += Lam[ k ][ 3 ] * DJac * dval;
1141  }
1142  }
1143 
1144  clear_2d( npts, Lam );
1145 
1146  for ( i = 0; i < 2; i++ )
1147  {
1148  Stif[ 0 ][ i ] = StifL[ 0 ][ i ];
1149  Stif[ 1 ][ i ] = StifL[ 1 ][ i ] + StifR[ 0 ][ i ];
1150  Stif[ 2 ][ i ] = StifR[ 1 ][ i ];
1151  Stif[ 3 ][ i ] = StifL[ 3 ][ i ];
1152  Stif[ 4 ][ i ] = StifL[ 2 ][ i ] + StifR[ 3 ][ i ];
1153  Stif[ 5 ][ i ] = StifR[ 2 ][ i ];
1154  }
1155 
1156  clear_2d( 4, StifL );
1157  clear_2d( 4, StifR );
1158 }
1159 
1160 
1161 void US_AstfemMath::IntQTn2( double* vx, double D, double sw2,
1162  double** Stif, double dt )
1163 {
1164  // element to define basis functions
1165  //
1166  int npts, i, k;
1167  double x_gauss, y_gauss, dval;
1168  double Lx[ 4 ];
1169  double Cx[ 4 ];
1170  double Rx[ 3 ], Ry[ 3 ];
1171  double Qx[ 4 ], Qy[ 4 ];
1172  double Tx[ 3 ], Ty[ 3 ];
1173  double DJac;
1174  double** StifL;
1175  double** StifR;
1176  double** Gs;
1177  double** Lam;
1178  double phiL [ 4 ];
1179  double phiLx[ 4 ];
1180  double phiLy[ 4 ];
1181  double phiCx[ 4 ];
1182  double phiCy[ 4 ];
1183  double phiC [ 4 ];
1184  double phiR [ 3 ];
1185  double phiRx[ 3 ];
1186  double phiRy[ 3 ];
1187 
1188  Lx[ 0 ] = vx[ 0 ];
1189  Lx[ 1 ] = vx[ 1 ];
1190  Lx[ 2 ] = vx[ 4 ];
1191  Lx[ 3 ] = vx[ 3 ]; // vertices of left T
1192 
1193  Cx[ 0 ] = vx[ 5 ];
1194  Cx[ 1 ] = vx[ 6 ];
1195  Cx[ 2 ] = vx[ 4 ];
1196  Cx[ 3 ] = vx[ 3 ];
1197 
1198  Rx[ 0 ] = vx[ 1 ]; // vertices of Q on right
1199  Rx[ 1 ] = vx[ 2 ];
1200  Rx[ 2 ] = vx[ 4 ];
1201 
1202  Ry[ 0 ] = 0.0;
1203  Ry[ 1 ] = 0.0;
1204  Ry[ 2 ] = dt;
1205 
1206  initialize_2d( 4, 2, &StifL );
1207  initialize_2d( 4, 2, &StifR );
1208 
1209  //
1210  // integration over element Q
1211  //
1212  Qx[ 0 ] = vx[ 5 ]; // vertices of Q on right
1213  Qx[ 1 ] = vx[ 1 ];
1214  Qx[ 2 ] = vx[ 4 ];
1215  Qx[ 3 ] = vx[ 3 ];
1216 
1217  Qy[ 0 ] = 0.0;
1218  Qy[ 1 ] = 0.0;
1219  Qy[ 2 ] = dt;
1220  Qy[ 3 ] = dt;
1221 
1222  npts = 5 * 5;
1223  initialize_2d( npts, 3, &Gs );
1224  DefineGaussian( 5, (double**)Gs );
1225 
1226  double psi[ 4 ], psi1[ 4 ], psi2[ 4 ], jac[ 4 ];
1227 
1228  for ( k = 0; k < npts; k++ )
1229  {
1230  BasisQS( Gs[ k ][ 0 ], Gs[ k ][ 1 ], psi, psi1, psi2 );
1231 
1232  x_gauss = 0.0;
1233  y_gauss = 0.0;
1234 
1235  for ( i = 0; i < 4; i++ )
1236  {
1237  jac[ i ] = 0.0;
1238  }
1239 
1240  for ( i = 0; i < 4; i++ )
1241  {
1242  x_gauss += psi[ i ] * Qx[ i ];
1243  y_gauss += psi[ i ] * Qy[ i ];
1244  jac[ 0 ] += Qx[ i ] * psi1[ i ];
1245  jac[ 1 ] += Qx[ i ] * psi2[ i ];
1246  jac[ 2 ] += Qy[ i ] * psi1[ i ];
1247  jac[ 3 ] += Qy[ i ] * psi2[ i ];
1248  }
1249 
1250  DJac = jac[ 0 ] * jac[ 3 ] - jac[ 1] * jac[ 2 ];
1251 
1252  //
1253  // find phi, phi_x, phi_y on L and C at (x,y)
1254  //
1255 
1256  BasisQR( Lx, x_gauss, y_gauss, phiL, phiLx, phiLy, dt );
1257  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
1258 
1259  for ( i = 0; i < 4; i++ )
1260  {
1261  dval = Integrand( x_gauss, D, sw2, phiL[ i ], phiLx[ i ], phiLy[ i ],
1262  phiC[ 0 ] + phiC[ 3 ], phiCx[ 0 ] + phiCx[ 3 ] );
1263  StifL[ i ][ 0 ] += Gs[ k ][ 2 ] * DJac * dval;
1264 
1265  dval = Integrand( x_gauss, D, sw2, phiL[ i ], phiLx[ i ], phiLy[ i ],
1266  phiC[ 1 ] + phiC[ 2 ], phiCx[ 1 ] + phiCx[ 2 ] );
1267  StifL[ i ][ 1 ] += Gs[ k ][ 2 ] * DJac * dval;
1268  }
1269  }
1270 
1271  clear_2d( npts, Gs );
1272 
1273  //
1274  // integration over T:
1275  //
1276  Tx[ 0 ] = vx[ 1 ]; // vertices of T on right
1277  Tx[ 1 ] = vx[ 6 ];
1278  Tx[ 2 ] = vx[ 4 ];
1279 
1280  Ty[ 0 ] = 0.0;
1281  Ty[ 1 ] = 0.0;
1282  Ty[ 2 ] = dt;
1283 
1284  npts = 28;
1285  initialize_2d( npts, 4, &Lam );
1286  DefineFkp( npts, (double**)Lam );
1287 
1288  for ( k = 0; k < npts; k++ )
1289  {
1290  x_gauss = Lam[ k ][ 0 ] * Tx[ 0 ] + Lam[ k ][ 1 ] *
1291  Tx[ 1 ] + Lam[ k ][ 2 ] * Tx[ 2 ];
1292  y_gauss = Lam[ k ][ 0 ] * Ty[ 0 ] + Lam[ k ][ 1 ] *
1293  Ty[ 1 ] + Lam[ k ][ 2 ] * Ty[ 2 ];
1294  DJac = 2.0 * AreaT( Tx, Ty );
1295 
1296  if ( DJac < 1.e-22 ) break;
1297 
1298  //
1299  // find phi, phi_x, phi_y on R and C at (x,y)
1300  //
1301 
1302  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
1303  BasisTR( Rx, Ry, x_gauss, y_gauss, phiR, phiRx, phiRy );
1304 
1305  for ( i = 0; i < 3; i++ )
1306  {
1307  dval = Integrand( x_gauss, D, sw2, phiR[ i ], phiRx[ i ], phiRy[ i ],
1308  phiC[ 0 ] + phiC[ 3 ], phiCx[ 0 ] + phiCx[ 3 ] );
1309  StifR[ i ][ 0 ] += Lam[ k ][ 3 ] * DJac * dval;
1310 
1311  dval = Integrand( x_gauss, D, sw2, phiR[ i ], phiRx[ i ], phiRy[ i ],
1312  phiC[ 1 ] + phiC[ 2 ], phiCx[ 1 ] + phiCx[ 2 ] );
1313  StifR[ i ][ 1 ] += Lam[ k ][ 3 ] * DJac * dval;
1314  }
1315  }
1316 
1317  clear_2d( npts, Lam );
1318 
1319  for ( i = 0; i < 2; i++ )
1320  {
1321  Stif[ 0 ][ i ] = StifL[ 0 ][ i ];
1322  Stif[ 1 ][ i ] = StifL[ 1 ][ i ] + StifR[ 0 ][ i ] ;
1323  Stif[ 2 ][ i ] = StifR[ 1 ][ i ];
1324  Stif[ 3 ][ i ] = StifL[ 3 ][ i ];
1325  Stif[ 4 ][ i ] = StifL[ 2 ][ i ] + StifR[ 2 ][ i ];
1326  }
1327 
1328  clear_2d( 4, StifL );
1329  clear_2d( 4, StifR );
1330 }
1331 
1332 void US_AstfemMath::IntQTn1( double* vx, double D, double sw2,
1333  double** Stif, double dt )
1334 {
1335  // element to define basis functions
1336  //
1337  int npts, i, k;
1338  double x_gauss, y_gauss, dval;
1339  double Lx[ 3 ], Ly[ 3 ];
1340  double Tx[ 3 ], Ty[ 3 ];
1341  double** StifR;
1342  double** Lam;
1343  double DJac;
1344  double phiL [ 4 ];
1345  double phiLx[ 4 ];
1346  double phiLy[ 4 ];
1347 
1348  Lx[ 0 ] = vx[ 0 ];
1349  Lx[ 1 ] = vx[ 1 ];
1350  Lx[ 2 ] = vx[ 2 ];
1351 
1352  Ly[ 0 ] = 0.0;
1353  Ly[ 1 ] = 0.0;
1354  Ly[ 2 ] = dt;
1355 
1356  initialize_2d( 4, 2, &StifR );
1357 
1358  //
1359  // integration over T:
1360  //
1361  Tx[ 0 ] = vx[ 3 ]; // vertices of T on right
1362  Tx[ 1 ] = vx[ 1 ];
1363  Tx[ 2 ] = vx[ 2 ];
1364 
1365  Ty[ 0 ] = 0.0;
1366  Ty[ 1 ] = 0.0;
1367  Ty[ 2 ] = dt;
1368 
1369  npts = 28;
1370  initialize_2d( npts, 4, &Lam );
1371  DefineFkp( npts, Lam );
1372 
1373  for ( k = 0; k < npts; k++ )
1374  {
1375  x_gauss = Lam[ k ][ 0 ] * Tx[ 0 ] + Lam[ k ][ 1 ] *
1376  Tx[ 1 ] + Lam[ k ][ 2 ] * Tx[ 2 ];
1377  y_gauss = Lam[ k ][ 0 ] * Ty[ 0 ] + Lam[ k ][ 1 ] *
1378  Ty[ 1 ] + Lam[ k ][ 2 ] * Ty[ 2 ];
1379  DJac = 2.0 * AreaT( Tx, Ty );
1380 
1381  if ( DJac < 1.e-22 ) break;
1382 
1383  // find phi, phi_x, phi_y on R and C at (x,y)
1384 
1385 
1386  BasisTR( Lx, Ly, x_gauss, y_gauss, phiL, phiLx, phiLy );
1387 
1388  for ( i = 0; i < 3; i++ )
1389  {
1390  dval = Integrand( x_gauss, D, sw2, phiL[ i ], phiLx[ i ], phiLy[ i ],
1391  1.0, 0.0 );
1392  StifR[ i ][ 0 ] += Lam[ k ][ 3 ] * DJac * dval;
1393  }
1394  }
1395 
1396  for ( i = 0; i < 2; i++ )
1397  {
1398  Stif[ 0 ][ i ] = StifR[ 0 ][ i ];
1399  Stif[ 1 ][ i ] = StifR[ 1 ][ i ];
1400  Stif[ 2 ][ i ] = StifR[ 2 ][ i ];
1401  }
1402 
1403  clear_2d( npts, Lam );
1404  clear_2d( 4, StifR );
1405 }
1406 
1407 void US_AstfemMath::DefineFkp( int npts, double** Lam )
1408 {
1409  // source: http://people.scs.fsu.edu/~burkardt/datasets/
1410  // quadrature_rules_tri/quadrature_rules_tri.html
1411 
1412  switch ( npts )
1413  {
1414  case 12:
1415  // STRANG9, order 12, degree of precision 6.
1416  {
1417  Lam[ 0 ][ 0 ] = 0.873821971016996;
1418  Lam[ 1 ][ 0 ] = 0.063089014491502;
1419  Lam[ 2 ][ 0 ] = 0.063089014491502;
1420  Lam[ 3 ][ 0 ] = 0.501426509658179;
1421  Lam[ 4 ][ 0 ] = 0.249286745170910;
1422  Lam[ 5 ][ 0 ] = 0.249286745170910;
1423  Lam[ 6 ][ 0 ] = 0.636502499121399;
1424  Lam[ 7 ][ 0 ] = 0.636502499121399;
1425  Lam[ 8 ][ 0 ] = 0.310352451033785;
1426  Lam[ 9 ][ 0 ] = 0.310352451033785;
1427  Lam[ 10 ][ 0 ] = 0.053145049844816;
1428  Lam[ 11 ][ 0 ] = 0.053145049844816;
1429 
1430  Lam[ 0 ][ 1 ] = 0.063089014491502;
1431  Lam[ 1 ][ 1 ] = 0.873821971016996;
1432  Lam[ 2 ][ 1 ] = 0.063089014491502;
1433  Lam[ 3 ][ 1 ] = 0.249286745170910;
1434  Lam[ 4 ][ 1 ] = 0.501426509658179;
1435  Lam[ 5 ][ 1 ] = 0.249286745170910;
1436  Lam[ 6 ][ 1 ] = 0.310352451033785;
1437  Lam[ 7 ][ 1 ] = 0.053145049844816;
1438  Lam[ 8 ][ 1 ] = 0.636502499121399;
1439  Lam[ 9 ][ 1 ] = 0.053145049844816;
1440  Lam[ 10 ][ 1 ] = 0.636502499121399;
1441  Lam[ 11 ][ 1 ] = 0.310352451033785;
1442 
1443  Lam[ 0 ][ 3 ] = 0.050844906370207;
1444  Lam[ 1 ][ 3 ] = 0.050844906370207;
1445  Lam[ 2 ][ 3 ] = 0.050844906370207;
1446  Lam[ 3 ][ 3 ] = 0.116786275726379;
1447  Lam[ 4 ][ 3 ] = 0.116786275726379;
1448  Lam[ 5 ][ 3 ] = 0.116786275726379;
1449  Lam[ 6 ][ 3 ] = 0.082851075618374;
1450  Lam[ 7 ][ 3 ] = 0.082851075618374;
1451  Lam[ 8 ][ 3 ] = 0.082851075618374;
1452  Lam[ 9 ][ 3 ] = 0.082851075618374;
1453  Lam[ 10 ][ 3 ] = 0.082851075618374;
1454  Lam[ 11 ][ 3 ] = 0.082851075618374;
1455  break;
1456  }
1457  case 28:
1458  // TOMS612_28, order 28, degree of precision 11,
1459  // a rule from ACM TOMS algorithm #612
1460  {
1461  Lam[ 0 ][ 0 ] = 0.33333333333333333;
1462  Lam[ 1 ][ 0 ] = 0.9480217181434233;
1463  Lam[ 2 ][ 0 ] = 0.02598914092828833;
1464  Lam[ 3 ][ 0 ] = 0.02598914092828833;
1465  Lam[ 4 ][ 0 ] = 0.8114249947041546;
1466  Lam[ 5 ][ 0 ] = 0.09428750264792270;
1467  Lam[ 6 ][ 0 ] = 0.09428750264792270;
1468  Lam[ 7 ][ 0 ] = 0.01072644996557060;
1469  Lam[ 8 ][ 0 ] = 0.4946367750172147;
1470  Lam[ 9 ][ 0 ] = 0.4946367750172147;
1471  Lam[ 10 ][ 0 ] = 0.5853132347709715;
1472  Lam[ 11 ][ 0 ] = 0.2073433826145142;
1473  Lam[ 12 ][ 0 ] = 0.2073433826145142;
1474  Lam[ 13 ][ 0 ] = 0.1221843885990187;
1475  Lam[ 14 ][ 0 ] = 0.4389078057004907;
1476  Lam[ 15 ][ 0 ] = 0.4389078057004907;
1477  Lam[ 16 ][ 0 ] = 0.6779376548825902;
1478  Lam[ 17 ][ 0 ] = 0.6779376548825902;
1479  Lam[ 18 ][ 0 ] = 0.04484167758913055;
1480  Lam[ 19 ][ 0 ] = 0.04484167758913055;
1481  Lam[ 20 ][ 0 ] = 0.27722066752827925;
1482  Lam[ 21 ][ 0 ] = 0.27722066752827925;
1483  Lam[ 22 ][ 0 ] = 0.8588702812826364;
1484  Lam[ 23 ][ 0 ] = 0.8588702812826364;
1485  Lam[ 24 ][ 0 ] = 0.0000000000000000;
1486  Lam[ 25 ][ 0 ] = 0.0000000000000000;
1487  Lam[ 26 ][ 0 ] = 0.1411297187173636;
1488  Lam[ 27 ][ 0 ] = 0.1411297187173636;
1489 
1490  Lam[ 0 ][ 1 ] = 0.333333333333333333;
1491  Lam[ 1 ][ 1 ] = 0.02598914092828833;
1492  Lam[ 2 ][ 1 ] = 0.9480217181434233;
1493  Lam[ 3 ][ 1 ] = 0.02598914092828833;
1494  Lam[ 4 ][ 1 ] = 0.09428750264792270;
1495  Lam[ 5 ][ 1 ] = 0.8114249947041546;
1496  Lam[ 6 ][ 1 ] = 0.09428750264792270;
1497  Lam[ 7 ][ 1 ] = 0.4946367750172147;
1498  Lam[ 8 ][ 1 ] = 0.01072644996557060;
1499  Lam[ 9 ][ 1 ] = 0.4946367750172147;
1500  Lam[ 10 ][ 1 ] = 0.2073433826145142;
1501  Lam[ 11 ][ 1 ] = 0.5853132347709715;
1502  Lam[ 12 ][ 1 ] = 0.2073433826145142;
1503  Lam[ 13 ][ 1 ] = 0.4389078057004907;
1504  Lam[ 14 ][ 1 ] = 0.1221843885990187;
1505  Lam[ 15 ][ 1 ] = 0.4389078057004907;
1506  Lam[ 16 ][ 1 ] = 0.04484167758913055;
1507  Lam[ 17 ][ 1 ] = 0.27722066752827925;
1508  Lam[ 18 ][ 1 ] = 0.6779376548825902;
1509  Lam[ 19 ][ 1 ] = 0.27722066752827925;
1510  Lam[ 20 ][ 1 ] = 0.6779376548825902;
1511  Lam[ 21 ][ 1 ] = 0.04484167758913055;
1512  Lam[ 22 ][ 1 ] = 0.00000000000000000;
1513  Lam[ 23 ][ 1 ] = 0.1411297187173636;
1514  Lam[ 24 ][ 1 ] = 0.8588702812826364;
1515  Lam[ 25 ][ 1 ] = 0.1411297187173636;
1516  Lam[ 26 ][ 1 ] = 0.8588702812826364;
1517  Lam[ 27 ][ 1 ] = 0.0000000000000000;
1518 
1519  Lam[ 0 ][ 3 ] = 0.08797730116222190;
1520  Lam[ 1 ][ 3 ] = 0.008744311553736190;
1521  Lam[ 2 ][ 3 ] = 0.008744311553736190;
1522  Lam[ 3 ][ 3 ] = 0.008744311553736190;
1523  Lam[ 4 ][ 3 ] = 0.03808157199393533;
1524  Lam[ 5 ][ 3 ] = 0.03808157199393533;
1525  Lam[ 6 ][ 3 ] = 0.03808157199393533;
1526  Lam[ 7 ][ 3 ] = 0.01885544805613125;
1527  Lam[ 8 ][ 3 ] = 0.01885544805613125;
1528  Lam[ 9 ][ 3 ] = 0.01885544805613125;
1529  Lam[ 10 ][ 3 ] = 0.07215969754474100;
1530  Lam[ 11 ][ 3 ] = 0.07215969754474100;
1531  Lam[ 12 ][ 3 ] = 0.07215969754474100;
1532  Lam[ 13 ][ 3 ] = 0.06932913870553720;
1533  Lam[ 14 ][ 3 ] = 0.06932913870553720;
1534  Lam[ 15 ][ 3 ] = 0.06932913870553720;
1535  Lam[ 16 ][ 3 ] = 0.04105631542928860;
1536  Lam[ 17 ][ 3 ] = 0.04105631542928860;
1537  Lam[ 18 ][ 3 ] = 0.04105631542928860;
1538  Lam[ 19 ][ 3 ] = 0.04105631542928860;
1539  Lam[ 20 ][ 3 ] = 0.04105631542928860;
1540  Lam[ 21 ][ 3 ] = 0.04105631542928860;
1541  Lam[ 22 ][ 3 ] = 0.007362383783300573;
1542  Lam[ 23 ][ 3 ] = 0.007362383783300573;
1543  Lam[ 24 ][ 3 ] = 0.007362383783300573;
1544  Lam[ 25 ][ 3 ] = 0.007362383783300573;
1545  Lam[ 26 ][ 3 ] = 0.007362383783300573;
1546  Lam[ 27 ][ 3 ] = 0.007362383783300573;
1547  break;
1548  }
1549  case 37:
1550  // TOMS706_37, order 37, degree of precision 13, a rule from ACM TOMS algorithm #706.
1551  {
1552  Lam[ 0 ][ 0 ] = 0.333333333333333333333333333333;
1553  Lam[ 1 ][ 0 ] = 0.950275662924105565450352089520;
1554  Lam[ 2 ][ 0 ] = 0.024862168537947217274823955239;
1555  Lam[ 3 ][ 0 ] = 0.024862168537947217274823955239;
1556  Lam[ 4 ][ 0 ] = 0.171614914923835347556304795551;
1557  Lam[ 5 ][ 0 ] = 0.414192542538082326221847602214;
1558  Lam[ 6 ][ 0 ] = 0.414192542538082326221847602214;
1559  Lam[ 7 ][ 0 ] = 0.539412243677190440263092985511;
1560  Lam[ 8 ][ 0 ] = 0.230293878161404779868453507244;
1561  Lam[ 9 ][ 0 ] = 0.230293878161404779868453507244;
1562  Lam[ 10 ][ 0 ] = 0.772160036676532561750285570113;
1563  Lam[ 11 ][ 0 ] = 0.113919981661733719124857214943;
1564  Lam[ 12 ][ 0 ] = 0.113919981661733719124857214943;
1565  Lam[ 13 ][ 0 ] = 0.009085399949835353883572964740;
1566  Lam[ 14 ][ 0 ] = 0.495457300025082323058213517632;
1567  Lam[ 15 ][ 0 ] = 0.495457300025082323058213517632;
1568  Lam[ 16 ][ 0 ] = 0.062277290305886993497083640527;
1569  Lam[ 17 ][ 0 ] = 0.468861354847056503251458179727;
1570  Lam[ 18 ][ 0 ] = 0.468861354847056503251458179727;
1571  Lam[ 19 ][ 0 ] = 0.022076289653624405142446876931;
1572  Lam[ 20 ][ 0 ] = 0.022076289653624405142446876931;
1573  Lam[ 21 ][ 0 ] = 0.851306504174348550389457672223;
1574  Lam[ 22 ][ 0 ] = 0.851306504174348550389457672223;
1575  Lam[ 23 ][ 0 ] = 0.126617206172027096933163647918;
1576  Lam[ 24 ][ 0 ] = 0.126617206172027096933163647918;
1577  Lam[ 25 ][ 0 ] = 0.018620522802520968955913511549;
1578  Lam[ 26 ][ 0 ] = 0.018620522802520968955913511549;
1579  Lam[ 27 ][ 0 ] = 0.689441970728591295496647976487;
1580  Lam[ 28 ][ 0 ] = 0.689441970728591295496647976487;
1581  Lam[ 29 ][ 0 ] = 0.291937506468887771754472382212;
1582  Lam[ 30 ][ 0 ] = 0.291937506468887771754472382212;
1583  Lam[ 31 ][ 0 ] = 0.096506481292159228736516560903;
1584  Lam[ 32 ][ 0 ] = 0.096506481292159228736516560903;
1585  Lam[ 33 ][ 0 ] = 0.635867859433872768286976979827;
1586  Lam[ 34 ][ 0 ] = 0.635867859433872768286976979827;
1587  Lam[ 35 ][ 0 ] = 0.267625659273967961282458816185;
1588  Lam[ 36 ][ 0 ] = 0.267625659273967961282458816185;
1589 
1590  Lam[ 0 ][ 1 ] = 0.333333333333333333333333333333;
1591  Lam[ 1 ][ 1 ] = 0.024862168537947217274823955239;
1592  Lam[ 2 ][ 1 ] = 0.950275662924105565450352089520;
1593  Lam[ 3 ][ 1 ] = 0.024862168537947217274823955239;
1594  Lam[ 4 ][ 1 ] = 0.414192542538082326221847602214;
1595  Lam[ 5 ][ 1 ] = 0.171614914923835347556304795551;
1596  Lam[ 6 ][ 1 ] = 0.414192542538082326221847602214;
1597  Lam[ 7 ][ 1 ] = 0.230293878161404779868453507244;
1598  Lam[ 8 ][ 1 ] = 0.539412243677190440263092985511;
1599  Lam[ 9 ][ 1 ] = 0.230293878161404779868453507244;
1600  Lam[ 10 ][ 1 ] = 0.113919981661733719124857214943;
1601  Lam[ 11 ][ 1 ] = 0.772160036676532561750285570113;
1602  Lam[ 12 ][ 1 ] = 0.113919981661733719124857214943;
1603  Lam[ 13 ][ 1 ] = 0.495457300025082323058213517632;
1604  Lam[ 14 ][ 1 ] = 0.009085399949835353883572964740;
1605  Lam[ 15 ][ 1 ] = 0.495457300025082323058213517632;
1606  Lam[ 16 ][ 1 ] = 0.468861354847056503251458179727;
1607  Lam[ 17 ][ 1 ] = 0.062277290305886993497083640527;
1608  Lam[ 18 ][ 1 ] = 0.468861354847056503251458179727;
1609  Lam[ 19 ][ 1 ] = 0.851306504174348550389457672223;
1610  Lam[ 20 ][ 1 ] = 0.126617206172027096933163647918;
1611  Lam[ 21 ][ 1 ] = 0.022076289653624405142446876931;
1612  Lam[ 22 ][ 1 ] = 0.126617206172027096933163647918;
1613  Lam[ 23 ][ 1 ] = 0.022076289653624405142446876931;
1614  Lam[ 24 ][ 1 ] = 0.851306504174348550389457672223;
1615  Lam[ 25 ][ 1 ] = 0.689441970728591295496647976487;
1616  Lam[ 26 ][ 1 ] = 0.291937506468887771754472382212;
1617  Lam[ 27 ][ 1 ] = 0.018620522802520968955913511549;
1618  Lam[ 28 ][ 1 ] = 0.291937506468887771754472382212;
1619  Lam[ 29 ][ 1 ] = 0.018620522802520968955913511549;
1620  Lam[ 30 ][ 1 ] = 0.689441970728591295496647976487;
1621  Lam[ 31 ][ 1 ] = 0.635867859433872768286976979827;
1622  Lam[ 32 ][ 1 ] = 0.267625659273967961282458816185;
1623  Lam[ 33 ][ 1 ] = 0.096506481292159228736516560903;
1624  Lam[ 34 ][ 1 ] = 0.267625659273967961282458816185;
1625  Lam[ 35 ][ 1 ] = 0.096506481292159228736516560903;
1626  Lam[ 36 ][ 1 ] = 0.635867859433872768286976979827;
1627 
1628  Lam[ 0 ][ 3 ] = 0.051739766065744133555179145422;
1629  Lam[ 1 ][ 3 ] = 0.008007799555564801597804123460;
1630  Lam[ 2 ][ 3 ] = 0.008007799555564801597804123460;
1631  Lam[ 3 ][ 3 ] = 0.008007799555564801597804123460;
1632  Lam[ 4 ][ 3 ] = 0.046868898981821644823226732071;
1633  Lam[ 5 ][ 3 ] = 0.046868898981821644823226732071;
1634  Lam[ 6 ][ 3 ] = 0.046868898981821644823226732071;
1635  Lam[ 7 ][ 3 ] = 0.046590940183976487960361770070;
1636  Lam[ 8 ][ 3 ] = 0.046590940183976487960361770070;
1637  Lam[ 9 ][ 3 ] = 0.046590940183976487960361770070;
1638  Lam[ 10 ][ 3 ] = 0.031016943313796381407646220131;
1639  Lam[ 11 ][ 3 ] = 0.031016943313796381407646220131;
1640  Lam[ 12 ][ 3 ] = 0.031016943313796381407646220131;
1641  Lam[ 13 ][ 3 ] = 0.010791612736631273623178240136;
1642  Lam[ 14 ][ 3 ] = 0.010791612736631273623178240136;
1643  Lam[ 15 ][ 3 ] = 0.010791612736631273623178240136;
1644  Lam[ 16 ][ 3 ] = 0.032195534242431618819414482205;
1645  Lam[ 17 ][ 3 ] = 0.032195534242431618819414482205;
1646  Lam[ 18 ][ 3 ] = 0.032195534242431618819414482205;
1647  Lam[ 19 ][ 3 ] = 0.015445834210701583817692900053;
1648  Lam[ 20 ][ 3 ] = 0.015445834210701583817692900053;
1649  Lam[ 21 ][ 3 ] = 0.015445834210701583817692900053;
1650  Lam[ 22 ][ 3 ] = 0.015445834210701583817692900053;
1651  Lam[ 23 ][ 3 ] = 0.015445834210701583817692900053;
1652  Lam[ 24 ][ 3 ] = 0.015445834210701583817692900053;
1653  Lam[ 25 ][ 3 ] = 0.017822989923178661888748319485;
1654  Lam[ 26 ][ 3 ] = 0.017822989923178661888748319485;
1655  Lam[ 27 ][ 3 ] = 0.017822989923178661888748319485;
1656  Lam[ 28 ][ 3 ] = 0.017822989923178661888748319485;
1657  Lam[ 29 ][ 3 ] = 0.017822989923178661888748319485;
1658  Lam[ 30 ][ 3 ] = 0.017822989923178661888748319485;
1659  Lam[ 31 ][ 3 ] = 0.037038683681384627918546472190;
1660  Lam[ 32 ][ 3 ] = 0.037038683681384627918546472190;
1661  Lam[ 33 ][ 3 ] = 0.037038683681384627918546472190;
1662  Lam[ 34 ][ 3 ] = 0.037038683681384627918546472190;
1663  Lam[ 35 ][ 3 ] = 0.037038683681384627918546472190;
1664  Lam[ 36 ][ 3 ] = 0.037038683681384627918546472190;
1665 
1666  break;
1667  }
1668  default:
1669  {
1670  return;
1671  }
1672  }
1673 
1674  for( int i = 0; i < npts; i++ )
1675  {
1676  Lam[ i ][ 2 ] = 1. - Lam[ i ][ 0 ] - Lam[ i ][ 1 ];
1677  // To make the sum( wt ) = 0.5 = area of standard elem
1678  Lam[ i ][ 3 ] /= 2.;
1679  }
1680 }
1681 
1682 // AreaT: area of a triangle (v1, v2, v3)
1683 
1684 double US_AstfemMath::AreaT( double* xv, double* yv )
1685 {
1686  return ( 0.5 * ( ( xv[ 1 ] - xv[ 0 ] ) * ( yv[ 2 ] - yv[ 0 ] )
1687  - ( xv[ 2 ] - xv[ 0 ] ) * ( yv[ 1 ] - yv[ 0 ] ) ) );
1688 }
1689 
1690 // Computer basis on standard element
1691 
1692 void US_AstfemMath::BasisTS( double xi, double et, double* phi,
1693  double* phi1, double* phi2 )
1694 {
1695  //function [phi, phi1, phi2] = BasisTS(xi,et)
1696 
1697  phi [ 0 ] = 1.0 - xi - et;
1698  phi1[ 0 ] = -1.0;
1699  phi2[ 0 ] = -1.0;
1700 
1701  phi [ 1 ] = xi;
1702  phi1[ 1 ] = 1.0;
1703  phi2[ 1 ] = 0.0;
1704 
1705  phi [ 2 ] = et;
1706  phi1[ 2 ] = 0.0;
1707  phi2[ 2 ] = 1.0;
1708 }
1709 
1710 // Computer basis on standard element
1711 
1712 void US_AstfemMath::BasisQS( double xi, double et, double* phi,
1713  double* phi1, double* phi2 )
1714 {
1715 
1716  phi [ 0 ] = ( 1.0 - xi ) * ( 1.0 - et );
1717  phi1[ 0 ] = -1.0 * ( 1.0 - et );
1718  phi2[ 0 ] = -1.0 * ( 1.0 - xi );
1719 
1720  phi [ 1 ] = xi * ( 1.0 - et );
1721  phi1[ 1 ] = 1.0 - et;
1722  phi2[ 1 ] = -xi;
1723 
1724  phi [ 2 ] = xi * et;
1725  phi1[ 2 ] = et;
1726  phi2[ 2 ] = xi;
1727 
1728  phi [ 3 ] = ( 1.0 - xi ) * et;
1729  phi1[ 3 ] = -et;
1730  phi2[ 3 ] = 1.0 - xi;
1731 }
1732 
1733 // Function BasisTR: compute basis on real element T:
1734 // phi, phi_x, phi_t at a given (xs,ts) point
1735 // the triangular is assumed to be (x1,y1), (x2, y2), (x3, y3)
1736 
1737 void US_AstfemMath::BasisTR( double* vx, double* vy,
1738  double xs, double ys, double* phi, double* phix, double* phiy )
1739 {
1740  // find (xi,et) corresponding to (xs, ts)
1741 
1742  int i;
1743  double tempv1[ 3 ], tempv2[ 3 ];
1744 
1745  tempv1[ 0 ] = xs;
1746  tempv1[ 1 ] = vx[ 2 ];
1747  tempv1[ 2 ] = vx[ 0 ];
1748  tempv2[ 0 ] = ys;
1749  tempv2[ 1 ] = vy[ 2 ];
1750  tempv2[ 2 ] = vy[ 0 ];
1751 
1752  double AreaK = AreaT( vx, vy );
1753  double xi = AreaT( tempv1, tempv2 ) / AreaK;
1754 
1755  tempv1[ 0 ] = xs;
1756  tempv1[ 1 ] = vx[ 0 ];
1757  tempv1[ 2 ] = vx[ 1 ];
1758  tempv2[ 0 ] = ys;
1759  tempv2[ 1 ] = vy[ 0 ];
1760  tempv2[ 2 ] = vy[ 1 ];
1761 
1762  double et = AreaT( tempv1, tempv2 ) / AreaK;
1763 
1764  double phi1[ 3 ];
1765  double phi2[ 3 ];
1766 
1767  BasisTS( xi, et, phi, phi1, phi2 );
1768 
1769  double Jac[ 4 ], JacN[ 4 ], det;
1770 
1771  Jac[ 0 ] = vx[ 0 ] * phi1[ 0 ] + vx[ 1 ] * phi1[ 1 ] + vx[ 2 ] * phi1[ 2 ];
1772  Jac[ 1 ] = vx[ 0 ] * phi2[ 0 ] + vx[ 1 ] * phi2[ 1 ] + vx[ 2 ] * phi2[ 2 ];
1773  Jac[ 2 ] = vy[ 0 ] * phi1[ 0 ] + vy[ 1 ] * phi1[ 1 ] + vy[ 2 ] * phi1[ 2 ];
1774  Jac[ 3 ] = vy[ 0 ] * phi2[ 0 ] + vy[ 1 ] * phi2[ 1 ] + vy[ 2 ] * phi2[ 2 ];
1775 
1776  det = Jac[ 0 ] * Jac[ 3 ] - Jac[ 1 ] * Jac [ 2 ];
1777 
1778  JacN[ 0 ] = Jac[3]/det;
1779  JacN[ 1 ] = -Jac[1]/det;
1780  JacN[ 2 ] = -Jac[2]/det;
1781  JacN[ 3 ] = Jac[0]/det;
1782 
1783  for ( i = 0; i < 3; i++ )
1784  {
1785  phix[ i ] = JacN[ 0 ] * phi1[ i ] + JacN[ 2 ] * phi2[ i ];
1786  phiy[ i ] = JacN[ 1 ] * phi1[ i ] + JacN[ 3 ] * phi2[ i ];
1787  }
1788 }
1789 
1790 void US_AstfemMath::BasisQR( double* vx, double xs, double ts,
1791  double* phi, double* phix, double* phiy, double dt )
1792 {
1793  int i;
1794 
1795  // find (xi,et) corresponding to (xs, ts)
1796  double et = ts / dt;
1797  double A = vx[ 0 ] * ( 1.0 - et ) + vx[ 3 ] * et;
1798  double B = vx[ 1 ] * ( 1.0 - et ) + vx[ 2 ] * et;
1799  double xi = ( xs - A ) / ( B - A );
1800 
1801  double phi1[ 4 ];
1802  double phi2[ 4 ];
1803 
1804  BasisQS( xi, et, phi, phi1, phi2 );
1805 
1806  double Jac[ 4 ], JacN[ 4 ], det;
1807 
1808  Jac[ 0] = vx[ 0 ] * phi1[ 0 ] + vx[ 1 ] * phi1[ 1 ] + vx[ 2 ] *
1809  phi1[ 2 ] + vx[ 3 ] * phi1[ 3 ];
1810  Jac[ 1] = vx[ 0 ] * phi2[ 0 ] + vx[ 1 ] * phi2[ 1 ] + vx[ 2 ] *
1811  phi2[ 2 ] + vx[ 3 ] * phi2[ 3 ];
1812  Jac[ 2] = dt * phi1[ 2 ] + dt * phi1[ 3 ];
1813  Jac[ 3] = dt * phi2[ 2 ] + dt * phi2[ 3 ];
1814 
1815  det = Jac[0] * Jac[3] - Jac[1] * Jac [2];
1816 
1817  JacN[ 0 ] = Jac[ 3 ] / det;
1818  JacN[ 1 ] = -Jac[ 1 ] / det;
1819  JacN[ 2 ] = -Jac[ 2 ] / det;
1820  JacN[ 3 ] = Jac[ 0 ] / det;
1821 
1822  for ( i = 0; i < 4; i++ )
1823  {
1824  phix[ i ] = JacN[ 0 ] * phi1[ i ] + JacN[ 2 ] * phi2[ i ];
1825  phiy[ i ] = JacN[ 1 ] * phi1[ i ] + JacN[ 3 ] * phi2[ i ];
1826  }
1827 }
1828 
1829 // Integrand for Lamm equation
1830 
1831 double US_AstfemMath::Integrand( double x_gauss, double D, double sw2,
1832  double u, double ux, double ut, double v, double vx )
1833 {
1834  return ( x_gauss * ut * v + D * x_gauss * ux * vx -
1835  sw2 * sq( x_gauss ) * u * vx );
1836 }
1837 
1838 // Source: http://www.math.ntnu.no/num/nnm/Program/Numlibc/
1839 
1840 void US_AstfemMath::DefineGaussian( int nGauss, double** Gs2 )
1841 {
1842  double* Gs1 = new double [ nGauss ];
1843  double* w = new double [ nGauss ];
1844 
1845  switch ( nGauss )
1846  {
1847  case 3:
1848  {
1849  Gs1[ 0 ] = -0.774596669241483; w[ 0 ] = 5.0 / 9.0;
1850  Gs1[ 1 ] = 0.0; w[ 1 ] = 8.0 / 9.0;
1851  Gs1[ 2 ] = 0.774596669241483; w[ 2 ] = 5.0 / 9.0;
1852  break;
1853  }
1854  case 5:
1855  {
1856  Gs1[ 0 ] = 0.906179845938664 ; w[ 0 ] = 0.236926885056189;
1857  Gs1[ 1 ] = 0.538469310105683 ; w[ 1 ] = 0.478628670499366;
1858  Gs1[ 2 ] = 0.000000000000000 ; w[ 2 ] = 0.568888888888889;
1859  Gs1[ 3 ] = -Gs1[ 1 ]; w[ 3 ] = w[ 1 ];
1860  Gs1[ 4 ] = -Gs1[ 0 ]; w[ 4 ] = w[ 0 ];
1861  break;
1862  }
1863  case 7:
1864  {
1865  Gs1[ 0 ] = 0.949107912342759 ; w[ 0 ] = 0.129484966168870;
1866  Gs1[ 1 ] = 0.741531185599394 ; w[ 1 ] = 0.279705391489277;
1867  Gs1[ 2 ] = 0.405845151377397 ; w[ 2 ] = 0.381830050505119;
1868  Gs1[ 3 ] = 0.000000000000000 ; w[ 3 ] = 0.417959183673469;
1869  Gs1[ 4 ] = -Gs1[ 2 ]; w[ 4 ] = w[ 2 ];
1870  Gs1[ 5 ] = -Gs1[ 1 ]; w[ 5 ] = w[ 1 ];
1871  Gs1[ 6 ] = -Gs1[ 0 ]; w[ 6 ] = w[ 0 ];
1872  break;
1873  }
1874  case 10:
1875  {
1876  Gs1[ 0 ] = 0.973906528517172 ; w[ 0 ] = 0.066671344308688;
1877  Gs1[ 1 ] = 0.865063366688985 ; w[ 1 ] = 0.149451349150581;
1878  Gs1[ 2 ] = 0.679409568299024 ; w[ 2 ] = 0.219086362515982;
1879  Gs1[ 3 ] = 0.433395394129247 ; w[ 3 ] = 0.269266719309996;
1880  Gs1[ 4 ] = 0.148874338981631 ; w[ 4 ] = 0.295524224714753;
1881  Gs1[ 5 ] = -Gs1[ 4 ]; w[ 5 ] = w[ 4 ];
1882  Gs1[ 6 ] = -Gs1[ 3 ]; w[ 6 ] = w[ 3 ];
1883  Gs1[ 7 ] = -Gs1[ 2 ]; w[ 7 ] = w[ 2 ];
1884  Gs1[ 8 ] = -Gs1[ 1 ]; w[ 8 ] = w[ 1 ];
1885  Gs1[ 9 ] = -Gs1[ 0 ]; w[ 9 ] = w[ 0 ];
1886  break;
1887  }
1888  default:
1889  {
1890  return;
1891  }
1892  }
1893 
1894  for ( int i = 0; i < nGauss; i++ )
1895  {
1896  for ( int j = 0; j < nGauss; j++ ) // map to [0,1] x [0,1]
1897  {
1898  int k = j + i * nGauss;
1899 
1900  Gs2[ k ][ 0 ] = ( Gs1[ i ] + 1.0 ) / 2.0;
1901  Gs2[ k ][ 1 ] = ( Gs1[ j ] + 1.0 ) / 2.0;
1902 
1903  Gs2[ k ][ 2 ] = w[ i ] * w[ j ] / 4.0;
1904  }
1905  }
1906 
1907  delete [] w;
1908  delete [] Gs1;
1909 }
1910 
1911 // initialize a simulation RawData object to have sizes,ranges,controls
1912 // that mirror those of an experimental EditedData object
1913 void US_AstfemMath::initSimData( US_DataIO::RawData& simdata,
1914  US_DataIO::EditedData& editdata, double concval1=0.0 )
1915 {
1916 
1917  int nconc = editdata.pointCount();
1918  int nscan = editdata.scanCount();
1919 
1920  // copy general control variable values
1921  simdata.type[ 0 ] = ' ';
1922  simdata.type[ 1 ] = ' ';
1923 
1924  for ( int jj = 0; jj < 16; jj++ )
1925  simdata.rawGUID[ jj ] = ' ';
1926 
1927  simdata.cell = editdata.cell.toInt();
1928  simdata.channel = editdata.channel.at( 0 ).toAscii();
1929  simdata.description = editdata.description;
1930  simdata.xvalues = editdata.xvalues;
1931 
1932  simdata.scanData.resize( nscan );
1933  double rvalue = concval1; // 1st scan constant concentration value
1934 
1935  for ( int ii = 0; ii < nscan; ii++ )
1936  { // Loop to copy scans
1937  US_DataIO::Scan* escan = &editdata.scanData[ ii ];
1938  US_DataIO::Scan sscan;
1939 
1940  sscan.temperature = escan->temperature;
1941  sscan.rpm = escan->rpm;
1942  sscan.seconds = escan->seconds;
1943  sscan.omega2t = escan->omega2t;
1944  sscan.wavelength = escan->wavelength;
1945  sscan.plateau = escan->plateau;
1946  sscan.delta_r = escan->delta_r;
1947  sscan.interpolated = escan->interpolated;
1948 
1949  sscan.rvalues.fill( rvalue, nconc );
1950 
1951  simdata.scanData[ ii ] = sscan;
1952 
1953  // Set values to zero for 2nd and subsequent scans
1954  rvalue = 0.0;
1955  }
1956 
1957  return;
1958 }
1959 
1960 // Calculate variance (average difference squared) between simulated and
1961 // experimental data.
1962 double US_AstfemMath::variance( US_DataIO::RawData& simdata,
1963  US_DataIO::EditedData& editdata )
1964 {
1965  QList< int > escns;
1966 
1967  return variance( simdata, editdata, escns );
1968 }
1969 
1970 // Calculate variance (average difference squared) between simulated and
1971 // experimental data.
1972 double US_AstfemMath::variance( US_DataIO::RawData& simdata,
1973  US_DataIO::EditedData& editdata,
1974  QList< int > exclScans )
1975 {
1976  int nscan = simdata .scanCount();
1977  int kscan = editdata.scanCount();
1978  int nconc = simdata .pointCount();
1979  int kconc = editdata.pointCount();
1980  double vari = 0.0;
1981 
1982  if ( nscan != kscan )
1983  {
1984  qDebug() << "*WARNING* variance(): sim/exp scan counts differ";
1985  nscan = ( nscan < kscan ) ? nscan : kscan;
1986  }
1987 
1988  if ( nconc != kconc )
1989  {
1990  qDebug() << "*WARNING* variance(): sim/exp readings counts differ";
1991  nconc = ( nconc < kconc ) ? nconc : kconc;
1992  }
1993 
1994  kscan = 0;
1995 
1996  for ( int ii = 0; ii < nscan; ii++ )
1997  { // accumulate sum of differences squared (readings in scans)
1998 
1999  if ( exclScans.contains( ii ) ) continue;
2000 
2001  kscan++;
2002 
2003  for ( int jj = 0; jj < nconc; jj++ )
2004  {
2005  vari += sq( simdata.value( ii, jj ) - editdata.value( ii, jj ) );
2006  }
2007  }
2008 
2009  vari /= (double)( kscan * nconc ); // variance is average diff-squared
2010 
2011  return vari;
2012 }
2013 
2014 // Calculate bottom radius value using channel bottom and rotor coeff. array
2015 double US_AstfemMath::calc_bottom( double rpm, double bottom_chan,
2016  double* rotorcoefs )
2017 {
2018  return ( bottom_chan
2019  + rotorcoefs[ 0 ] * rpm
2020  + rotorcoefs[ 1 ] * sq( rpm ) );
2021 }
2022 
2023 
2024 #ifdef NEVER
2025 #if defined(DEBUG_ALLOC)
2026 
2027 struct _created_matrices {
2028  long addr;
2029  int val1;
2030  int val2;
2031 };
2032 
2033 static list <_created_matrices> created_matrices;
2034 
2035 void init_matrices_alloc()
2036 {
2037  created_matrices.clear();
2038  puts( "init_matrices_alloc()" );
2039 }
2040 
2041 void list_matrices_alloc()
2042 {
2043  if ( created_matrices.size() )
2044  {
2045  printf( "allocated matrices:\n" 0 );
2046  }
2047 
2048  list<_created_matrices>::iterator Li;
2049 
2050  for ( Li = created_matrices.begin(); Li != created_matrices.end(); Li++ )
2051  {
2052  printf( "addr %lx val1 %u val2 %u\n", Li->addr, Li->val1, Li->val2 );
2053  }
2054 }
2055 
2056 static list<_created_matrices>::iterator find_matrices_alloc( long addr )
2057 {
2058  list<_created_matrices>::iterator Li;
2059 
2060  for ( Li = created_matrices.begin(); Li != created_matrices.end(); Li++ )
2061  {
2062  if ( Li->addr == addr )
2063  {
2064  return Li;
2065  }
2066  }
2067  return Li;
2068 }
2069 
2070 #endif
2071 
2072 
2073 
2074 
2075 
2076 
2077 
2078 
2079 
2080 
2081 
2082 
2083 void IntQT1_ellam(QVector <double> vx, double D, double sw2, double **Stif, double dt)
2084 {
2085  // element to define basis functions
2086  //
2087  int npts, i, k;
2088  double x_gauss, y_gauss, dval;
2089  QVector <double> Rx, Ry, Qx, Qy;
2090  double **StifR=NULL, DJac;
2091  double hh, slope, xn1, phiC, phiCx;
2092  double phiR [4];
2093  double phiRx[4];
2094  double phiRy[4];
2095 
2096  Rx.clear();
2097  Ry.clear();
2098  Rx.push_back(vx[0]); // vertices of Q on right
2099  Rx.push_back(vx[1]);
2100  Rx.push_back(vx[3]);
2101  Rx.push_back(vx[2]);
2102 
2103  Ry.push_back(0.0);
2104  Ry.push_back(0.0);
2105  Ry.push_back(dt);
2106  Ry.push_back(dt);
2107 
2108  initialize_2d(4, 2, &StifR);
2109  hh = vx[3] - vx[2];
2110  slope = (vx[3] - vx[4])/dt;
2111 
2112  //
2113  // integration over quadrilateral element Q :
2114  //
2115  if( (vx[1]-vx[4])/(vx[1]-vx[0]) <1.e-3 ) // Q_{0,4,3,2} is almost degenerated into a triangle
2116  {
2117  // elements for integration
2118  //
2119  Qx.clear();
2120  Qy.clear();
2121  Qx.push_back(vx[0]); // vertices of Q on right
2122  Qx.push_back(vx[3]);
2123  Qx.push_back(vx[2]);
2124  Qy.push_back(0.0);
2125  Qy.push_back(dt);
2126  Qy.push_back(dt);
2127 
2128  double **Lam;
2129  npts = 28;
2130  initialize_2d(npts, 4, &Lam);
2131  DefineFkp(npts, Lam);
2132 
2133  for (k=0; k<npts; k++)
2134  {
2135  x_gauss = Lam[k][0] * Qx[0] + Lam[k][1] * Qx[1] + Lam[k][2] * Qx[2];
2136  y_gauss = Lam[k][0] * Qy[0] + Lam[k][1] * Qy[1] + Lam[k][2] * Qy[2];
2137  DJac = 2.0 * AreaT(&Qx, &Qy);
2138 
2139  //
2140  // find phi, phi_x, phi_y on R and C at (x,y)
2141  //
2142 
2143  BasisQR( Rx, x_gauss, y_gauss, phiR, phiRx, phiRy, dt );
2144 
2145  xn1 = x_gauss + slope * ( dt - y_gauss ); // trace-forward point at t_n+1 from (x_g, y_g)
2146  phiC = ( xn1 - vx[2] )/hh; // hat function on t_n+1, =1 at vx[3]; =0 at vx[2]
2147  phiCx = 1./hh;
2148 
2149  for (i=0; i<4; i++)
2150  {
2151  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i], 1.-phiC, -phiCx);
2152  StifR[i][0] += Lam[k][3] * DJac * dval;
2153 
2154  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i], phiC, phiCx);
2155  StifR[i][1] += Lam[k][3] * DJac * dval;
2156  }
2157  }
2158  clear_2d(npts, Lam);
2159 
2160  }
2161  else
2162  { // Q_{0,4,3,2} is non-degenerate
2163  // elements for integration
2164  //
2165  Qx.clear();
2166  Qy.clear();
2167  Qx.push_back(vx[0]); // vertices of Q on right
2168  Qx.push_back(vx[4]);
2169  Qx.push_back(vx[3]);
2170  Qx.push_back(vx[2]);
2171 
2172  Qy.push_back(0.0);
2173  Qy.push_back(0.0);
2174  Qy.push_back(dt);
2175  Qy.push_back(dt);
2176 
2177  double **Gs=NULL;
2178  npts = 5 * 5;
2179  initialize_2d(npts, 3, &Gs);
2180  DefineGaussian(5, Gs);
2181 
2182  double psi[4], psi1[4], psi2[4], jac[4];
2183  for (k=0; k<npts; k++)
2184  {
2185  BasisQS(Gs[k][0], Gs[k][1], psi, psi1, psi2);
2186  x_gauss = 0.0;
2187  y_gauss = 0.0;
2188  for (i=0; i<4; i++)
2189  {
2190  jac[i] = 0.0;
2191  }
2192  for (i=0; i<4; i++)
2193  {
2194  x_gauss += psi[i] * Qx[i];
2195  y_gauss += psi[i] * Qy[i];
2196  jac[0] += Qx[i] * psi1[i];
2197  jac[1] += Qx[i] * psi2[i];
2198  jac[2] += Qy[i] * psi1[i];
2199  jac[3] += Qy[i] * psi2[i];
2200  }
2201  DJac = jac[0] * jac[3] - jac[1] * jac[2];
2202  //
2203  // find phi, phi_x, phi_y on L and C at (x,y)
2204  //
2205  BasisQR( Rx, x_gauss, y_gauss, phiR, phiRx, phiRy, dt );
2206  xn1 = x_gauss + slope * ( dt - y_gauss ); // trace-forward point at t_n+1 from (x_g, y_g)
2207  phiC = ( xn1 - vx[2] )/hh; // hat function on t_n+1, =1 at vx[3]; =0 at vx[2]
2208  phiCx = 1./hh;
2209  for (i=0; i<4; i++)
2210  {
2211  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i], 1.-phiC, -phiCx);
2212  StifR[i][0] += Gs[k][2] * DJac * dval;
2213  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i], phiC, phiCx);
2214  StifR[i][1] += Gs[k][2] * DJac * dval;
2215  }
2216  }
2217  clear_2d(npts, Gs);
2218  }
2219  for (i=0; i<2; i++)
2220  {
2221  Stif[0][i] = StifR[0][i];
2222  Stif[1][i] = StifR[1][i];
2223  Stif[2][i] = StifR[3][i];
2224  Stif[3][i] = StifR[2][i];
2225  }
2226 
2227  clear_2d(4, StifR);
2228 }
2229 
2230 void IntQTm_ellam(QVector <double> vx, double D, double sw2, double **Stif, double dt)
2231 {
2232  // element to define basis functions
2233  //
2234  int npts, i, k;
2235  double x_gauss, y_gauss, dval;
2236  QVector <double> Lx, Ly, Cx, Cy, Rx, Ry, Qx, Qy, Tx, Ty;
2237  double *phiR, *phiRx, *phiRy;
2238  double **StifL=NULL, **StifR=NULL, **Lam=NULL, DJac;
2239  double *phiL, *phiLx, *phiLy, *phiCx, *phiCy, *phiC;
2240  double **Gs=NULL;
2241  double phiL [4];
2242  double phiLx[4];
2243  double phiLy[4];
2244  double phiCx[4];
2245  double phiCy[4];
2246  double phiC [4];
2247  double phiR [4];
2248  double phiRx[4];
2249  double phiRy[4];
2250  Lx.clear();
2251  Ly.clear();
2252  Cx.clear();
2253  Cy.clear();
2254  Rx.clear();
2255  Ry.clear();
2256  Qx.clear();
2257  Qy.clear();
2258  Tx.clear();
2259  Ty.clear();
2260 
2261  Lx.push_back(vx[0]);
2262  Lx.push_back(vx[1]);
2263  Lx.push_back(vx[4]);
2264  Lx.push_back(vx[3]);
2265 
2266  Ly.push_back(0.0);
2267  Ly.push_back(0.0);
2268  Ly.push_back(dt);
2269  Ly.push_back(dt); // vertices of left T
2270 
2271  Cx.push_back(vx[6]);
2272  Cx.push_back(vx[7]);
2273  Cx.push_back(vx[5]);
2274  Cx.push_back(vx[4]);
2275 
2276  Cy.push_back(0.0);
2277  Cy.push_back(0.0);
2278  Cy.push_back(dt);
2279  Cy.push_back(dt);
2280 
2281  Rx.push_back(vx[1]); // vertices of Q on right
2282  Rx.push_back(vx[2]);
2283  Rx.push_back(vx[5]);
2284  Rx.push_back(vx[4]);
2285 
2286  Ry.push_back(0.0);
2287  Ry.push_back(0.0);
2288  Ry.push_back(dt);
2289  Ry.push_back(dt);
2290 
2291 
2292  initialize_2d(4, 2, &StifL);
2293  initialize_2d(4, 2, &StifR);
2294 
2295  //
2296  // integration over triangle T:
2297  //
2298  Tx.push_back(vx[6]); // vertices of T on left
2299  Tx.push_back(vx[1]);
2300  Tx.push_back(vx[4]);
2301 
2302  Ty.push_back(0.0);
2303  Ty.push_back(0.0);
2304  Ty.push_back(dt);
2305 
2306  npts = 28;
2307  initialize_2d(npts, 4, &Lam);
2308  DefineFkp(npts, Lam);
2309 
2310  for (k=0; k<npts; k++)
2311  {
2312  x_gauss = Lam[k][0] * Tx[0] + Lam[k][1] * Tx[1] + Lam[k][2] * Tx[2];
2313  y_gauss = Lam[k][0] * Ty[0] + Lam[k][1] * Ty[1] + Lam[k][2] * Ty[2];
2314  DJac = 2.0 * AreaT(&Tx, &Ty);
2315 
2316  if(DJac<1.e-22) break;
2317 
2318  //
2319  // find phi, phi_x, phi_y on R and C at (x,y)
2320  //
2321 
2322  BasisQR( Lx, x_gauss, y_gauss, phiL, phiLx, phiLy, dt );
2323  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
2324 
2325  for (i=0; i<4; i++)
2326  {
2327  dval = Integrand(x_gauss, D, sw2, phiL[i], phiLx[i], phiLy[i],
2328  phiC[0] + phiC[3], phiCx[0] + phiCx[3]);
2329  StifL[i][0] += Lam[k][3] * DJac * dval;
2330 
2331  dval = Integrand(x_gauss, D, sw2, phiL[i], phiLx[i], phiLy[i],
2332  phiC[1] + phiC[2], phiCx[1] + phiCx[2]);
2333  StifL[i][1] += Lam[k][3] * DJac * dval;
2334  }
2335  }
2336  clear_2d(npts, Lam);
2337 
2338  //
2339  // integration over quadrilateral element Q :
2340  //
2341  if( (vx[7]-vx[1])/(vx[2]-vx[1]) <1.e-3 ) // Q_{1,7,5,4} is almost degenerated into a triangle
2342  {
2343  Qx.push_back(vx[1]); // vertices of Q on right
2344  Qx.push_back(vx[5]);
2345  Qx.push_back(vx[4]);
2346  Qy.push_back(0.0);
2347  Qy.push_back(dt);
2348  Qy.push_back(dt);
2349 
2350  npts = 28;
2351  initialize_2d(npts, 4, &Lam);
2352  DefineFkp(npts, Lam);
2353 
2354  for (k=0; k<npts; k++)
2355  {
2356  x_gauss = Lam[k][0] * Qx[0] + Lam[k][1] * Qx[1] + Lam[k][2] * Qx[2];
2357  y_gauss = Lam[k][0] * Qy[0] + Lam[k][1] * Qy[1] + Lam[k][2] * Qy[2];
2358  DJac = 2.0 * AreaT(&Qx, &Qy);
2359 
2360  //
2361  // find phi, phi_x, phi_y on R and C at (x,y)
2362  //
2363 
2364  BasisQR( Rx, x_gauss, y_gauss, phiR, phiRx, phiRy, dt );
2365  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
2366  for (i=0; i<4; i++)
2367  {
2368  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i],
2369  phiC[0] + phiC[3], phiCx[0] + phiCx[3]);
2370  StifR[i][0] += Lam[k][3] * DJac * dval;
2371 
2372  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i],
2373  phiC[1] + phiC[2], phiCx[1] + phiCx[2]);
2374  StifR[i][1] += Lam[k][3] * DJac * dval;
2375  }
2376  }
2377  clear_2d(npts, Lam);
2378  }
2379  else // Q is a non-degenerate quadrilateral
2380  {
2381  Qx.push_back(vx[1]); // vertices of Q on right
2382  Qx.push_back(vx[7]);
2383  Qx.push_back(vx[5]);
2384  Qx.push_back(vx[4]);
2385 
2386  Qy.push_back(0.0);
2387  Qy.push_back(0.0);
2388  Qy.push_back(dt);
2389  Qy.push_back(dt);
2390 
2391  npts = 5 * 5;
2392  initialize_2d(npts, 3, &Gs);
2393  DefineGaussian(5, Gs);
2394 
2395  double psi[4], psi1[4], psi2[4], jac[4];
2396  for (k=0; k<npts; k++)
2397  {
2398  BasisQS(Gs[k][0], Gs[k][1], psi, psi1, psi2);
2399 
2400  x_gauss = 0.0;
2401  y_gauss = 0.0;
2402  for (i=0; i<4; i++)
2403  {
2404  jac[i] = 0.0;
2405  }
2406  for (i=0; i<4; i++)
2407  {
2408  x_gauss += psi[i] * Qx[i];
2409  y_gauss += psi[i] * Qy[i];
2410  jac[0] += Qx[i] * psi1[i];
2411  jac[1] += Qx[i] * psi2[i];
2412  jac[2] += Qy[i] * psi1[i];
2413  jac[3] += Qy[i] * psi2[i];
2414  }
2415 
2416  DJac = jac[0] * jac[3] - jac[1] * jac[2];
2417 
2418  //
2419  // find phi, phi_x, phi_y on L and C at (x,y)
2420  //
2421 
2422  BasisQR( Rx, x_gauss, y_gauss, phiR, phiRx, phiRy, dt );
2423  BasisQR( Cx, x_gauss, y_gauss, phiC, phiCx, phiCy, dt );
2424  for (i=0; i<4; i++)
2425  {
2426  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i],
2427  phiC[0] + phiC[3], phiCx[0] + phiCx[3]);
2428  StifR[i][0] += Gs[k][2] * DJac * dval;
2429 
2430  dval = Integrand(x_gauss, D, sw2, phiR[i], phiRx[i], phiRy[i],
2431  phiC[1] + phiC[2], phiCx[1] + phiCx[2]);
2432  StifR[i][1] += Gs[k][2] * DJac * dval;
2433  }
2434  }
2435  clear_2d(npts, Gs);
2436  }
2437 
2438  for (i=0; i<2; i++)
2439  {
2440  Stif[0][i] = StifL[0][i];
2441  Stif[1][i] = StifL[1][i] + StifR[0][i];
2442  Stif[2][i] = StifR[1][i];
2443  Stif[3][i] = StifL[3][i];
2444  Stif[4][i] = StifL[2][i] + StifR[3][i];
2445  Stif[5][i] = StifR[2][i];
2446  }
2447 
2448  clear_2d(4, StifL);
2449  clear_2d(4, StifR);
2450 }
2451 
2452 
2453 void IntQTn1_ellam(QVector <double> vx, double D, double sw2, double **Stif, double dt)
2454 {
2455  // element to define basis functions
2456  //
2457  int npts, i, k;
2458  double x_gauss, y_gauss, dval;
2459  QVector <double> Lx, Ly, Tx, Ty;
2460  double **StifL=NULL, **Lam=NULL, DJac;
2461  double phiL [4];
2462  double phiLx[4];
2463  double phiLy[4];
2464  Lx.clear();
2465  Ly.clear();
2466  Tx.clear();
2467  Ty.clear();
2468 
2469  Lx.push_back(vx[0]);
2470  Lx.push_back(vx[1]);
2471  Lx.push_back(vx[3]);
2472  Lx.push_back(vx[2]);
2473 
2474  Ly.push_back(0.0);
2475  Ly.push_back(0.0);
2476  Ly.push_back(dt);
2477  Ly.push_back(dt);
2478 
2479  initialize_2d(4, 2, &StifL);
2480 
2481  //
2482  // integration over T:
2483  //
2484 
2485  Tx.push_back(vx[4]); // triangle for integration
2486  Tx.push_back(vx[1]);
2487  Tx.push_back(vx[3]);
2488 
2489  Ty.push_back(0.0);
2490  Ty.push_back(0.0);
2491  Ty.push_back(dt);
2492 
2493  npts = 28;
2494  initialize_2d(npts, 4, &Lam);
2495  DefineFkp(npts, Lam);
2496 
2497  for (k=0; k<npts; k++)
2498  {
2499  x_gauss = Lam[k][0] * Tx[0] + Lam[k][1] * Tx[1] + Lam[k][2] * Tx[2];
2500  y_gauss = Lam[k][0] * Ty[0] + Lam[k][1] * Ty[1] + Lam[k][2] * Ty[2];
2501  DJac = 2.0 * AreaT(&Tx, &Ty);
2502 
2503  if(DJac<1.e-22) break;
2504  //
2505  // find phi, phi_x, phi_y on R and C at (x,y)
2506  //
2507 
2508  BasisQR( Lx, x_gauss, y_gauss, phiL, phiLx, phiLy, dt );
2509 
2510  for (i=0; i<4; i++)
2511  {
2512  dval = Integrand(x_gauss, D, sw2, phiL[i], phiLx[i], phiLy[i], 1.0, 0.0);
2513  StifL[i][0] += Lam[k][3] * DJac * dval;
2514  }
2515  }
2516  clear_2d(npts, Lam);
2517 
2518 
2519  for (i=0; i<2; i++)
2520  {
2521  Stif[0][i] = StifL[0][i];
2522  Stif[1][i] = StifL[1][i];
2523  Stif[2][i] = StifL[3][i];
2524  Stif[3][i] = StifL[2][i];
2525  }
2526 
2527  clear_2d(4, StifL);
2528 }
2529 
2530 
2531 void QuadSolver_ellam(double *ai, double *bi, double *ci, double *di, double *cr, double *solu, int N)
2532 {
2533 //
2534 // solve Quad-diagonal system [a_i, b_i, *c_i*, d_i]*[x]=[r_i]
2535 // c_i are on the main diagonal line
2536 //
2537 
2538  int i;
2539  double tmp;
2540  QVector<double> ca, cb, cc, cd;
2541  ca.clear();
2542  cb.clear();
2543  cc.clear();
2544  cd.clear();
2545  for (i=0; i<N; i++)
2546  {
2547  ca.push_back( ai[i] );
2548  cb.push_back( bi[i] );
2549  cc.push_back( ci[i] );
2550  cd.push_back( di[i] );
2551  }
2552 
2553  for (i=N-2; i>=1; i--)
2554  {
2555  tmp = cd[i]/cc[i+1];
2556  cc[i] = cc[i]-cb[i+1]*tmp;
2557  cb[i] = cb[i]-ca[i+1]*tmp;
2558  cr[i] = cr[i]-cr[i+1]*tmp;
2559  }
2560  i=0;
2561  tmp = cd[i]/cc[i+1];
2562  cc[i] = cc[i]-cb[i+1]*tmp;
2563  cr[i] = cr[i]-cr[i+1]*tmp;
2564  solu[0] = cr[0] / cc[0];
2565  solu[1] = (cr[1] - cb[1] * solu[0]) / cc[1];
2566  i = 1;
2567  do
2568  {
2569  i++;
2570  solu[i] = (cr[i] - ca[i] * solu[i-2] - cb[i] * solu[i-1]) / cc[i];
2571  } while (i != N-1);
2572 }
2573 // ******* end of ELLAM *********************
2574 
2575 
2577 // calculate total mass
2578 // (r,u) concentration defined at r(1), ...., r(M)
2579 // M: r(1).... r(M): the interval for total mass, (M-1) subintervals
2581 
2582 double IntConcentration(QVector<double> r, double *u)
2583 {
2584 //function T=IntConcentration(r,M,u)
2585  double T = 0.0;
2586  for ( int j=0; j<r.size()-1; j++)
2587  {
2588  T += (r[j+1] - r[j]) * ((r[j+1] - r[j]) * (u[j+1]/3.0 + u[j]/6.0)
2589  + r[j] * (u[j] + u[j+1])/2.0);
2590  }
2591  return T;
2592 }
2593 
2594 
2595 
2596 void DefInitCond(double **C0, int N)
2597 {
2598  int j;
2599  for(j=0; j<N; j++)
2600  {
2601  C0[0][j] = 0.3;
2602  C0[1][j] = 0.7;
2603  }
2604 }
2605 
2606 int interpolate( MfemData *simdata, int scans, int points,
2607  float *scantimes, double *radius, double **c )
2608 {
2609 
2610 /******************************************************************************
2611  * Interpolation: *
2612  * *
2613  * First, we need to interpolate the time. Create a new array with the same *
2614  * time dimensions as the raw data and the same radius dimensions as the *
2615  * simulated data. Then find the time steps from the simulated data that *
2616  * bracket the experimental data from the left and right. Then make a linear *
2617  * interpolation for the concentration values at each radius point from the *
2618  * simulated data. Then interpolate the radius points by linear interpolation.*
2619  * *
2620  ******************************************************************************/
2621 
2622  int i, j;
2623  double slope, intercept;
2624  double **ip_array;
2625  ip_array = new double* [scans];
2626  for (i=0; i<scans; i++)
2627  {
2628  ip_array[i] = new double [(*simdata).radius.size()];
2629  }
2630  int count = 0; // counting the number of time steps of the raw data
2631  for (i=0; i<scans; i++)
2632  {
2633  while (count < (*simdata).scan.size()-1
2634  && scantimes[i] >= (*simdata).scan[count].time)
2635  {
2636  count++;
2637  }
2638  if (scantimes[i] == (*simdata).scan[count].time)
2639  {
2640  for (j=0; j<(*simdata).radius.size(); j++)
2641  {
2642  ip_array[i][j] = (*simdata).scan[count].conc[j];
2643  }
2644  }
2645  else // else, perform a linear time interpolation:
2646  {
2647  for (j=0; j<(*simdata).radius.size(); j++)
2648  {
2649  slope = ((*simdata).scan[count].conc[j] - (*simdata).scan[count-1].conc[j])
2650  /((*simdata).scan[count].time - (*simdata).scan[count-1].time);
2651  intercept = (*simdata).scan[count].conc[j] - slope * (*simdata).scan[count].time;
2652  ip_array[i][j] = slope * scantimes[i] + intercept;
2653  }
2654  }
2655  }
2656  //interpolate radius then add to concentration vector
2657  for (i=0; i<scans; i++)
2658  {
2659  c[i][0] += ip_array[i][0]; // meniscus position is identical for all scans
2660  }
2661  // all other points may need interpolation:
2662  for (i=0; i<scans; i++)
2663  {
2664  count = 1;
2665  for (j=1; j<(int) points; j++)
2666  {
2667  while (radius[j] > (*simdata).radius[count] && count < (*simdata).radius.size()-1)
2668  {
2669  count++;
2670  }
2671  if (radius[j] == (*simdata).radius[count])
2672  {
2673  c[i][j] += ip_array[i][count];
2674  }
2675  else
2676  {
2677  slope = (ip_array[i][count] - ip_array[i][count-1])/((*simdata).radius[count] - (*simdata).radius[count-1]);
2678  intercept = ip_array[i][count] - (*simdata).radius[count] * slope;
2679  c[i][j] += slope * radius[j] + intercept;
2680  }
2681  }
2682  }
2683  for (i=0; i<scans; i++)
2684  {
2685  delete [] ip_array[i];
2686  }
2687  delete [] ip_array;
2688  return 0;
2689 }
2690 
2691 
2692 void interpolate_Cfinal( MfemInitial *C0, double *cfinal, QVector <double> *x )
2693 {
2694 // This routine also considers cases where C0 starts before the meniscus or
2695 // stops after the bottom is reached. However, C0 needs to cover both meniscus
2696 // and bottom. In those cases it will fill the C0 vector with the first and/or
2697 // last value of C0, respectively, for the missing points.
2698 
2699  int i;
2700  int j;
2701  int ja;
2702  double a;
2703  double b;
2704  double xs;
2705  double tmp;
2706 
2707  ja = 0;
2708  for(j=0; j<(*C0).radius.size(); j++)
2709  {
2710  xs = (*C0).radius[j];
2711  for (i=ja; i<(*x).size(); i++)
2712  {
2713  if( (*x)[i] > xs + 1.e-12)
2714  {
2715  break;
2716  }
2717  }
2718 
2719  if ( i == 0 ) // xs < x[0]
2720  {
2721  (*C0).concentration[j] = cfinal[0]; // use the first value
2722  }
2723  else if ( i == (*x).size() ) // xs > x[N]
2724  {
2725  (*C0).concentration[j] = cfinal[ i-1 ];
2726  }
2727  else // x[i-1] < xs <x[i]
2728  {
2729  a = (*x)[i-1];
2730  b = (*x)[i];
2731  tmp = (xs-a)/(b-a);
2732  (*C0).concentration[j] = cfinal[i-1] * (1. - tmp) + cfinal[i] * tmp;
2733  ja = i-1;
2734  }
2735  }
2736 }
2737 #endif
Generated on Sun Jun 21 2015 11:40:41 for UltraScan III by   doxygen 1.8.3.1-20130324