CG  Version 25
sm/src/forcing.h
Go to the documentation of this file.
1 // ==================================================================================
2 // Macro to evaluate the traveling (shock) wave solution
3 //
4 // evalSolution : true=eval solution, false=eval error
5 // ==================================================================================
6 #beginMacro getTravelingWaveSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
7 // printF("INFO: The traveling wave solutions are combinations of p and s solutions:\n"
8 // " [u1,u2] = ap* [ k1,k2] * G( k1*(x-xa)+k2*(y-ya) - cp*t ) (p-wave)\n"
9 // " [u1,u2] = as* [-k2,k1] * G( k1*(x-xa)+k2*(y-ya) - cs*t ) (s-wave\n",
10 // " where G(xi)=0 for xi>0 and G(xi)=-1 for xi<0 \n");
11 {
12  const int v1c = parameters.dbase.get<int >("v1c");
13  const int v2c = parameters.dbase.get<int >("v2c");
14  const int v3c = parameters.dbase.get<int >("v3c");
15 
16  bool assignVelocities= v1c>=0 ;
17  const int s11c = parameters.dbase.get<int >("s11c");
18  const int s12c = parameters.dbase.get<int >("s12c");
19  const int s13c = parameters.dbase.get<int >("s13c");
20  const int s21c = parameters.dbase.get<int >("s21c");
21  const int s22c = parameters.dbase.get<int >("s22c");
22  const int s23c = parameters.dbase.get<int >("s23c");
23  const int s31c = parameters.dbase.get<int >("s31c");
24  const int s32c = parameters.dbase.get<int >("s32c");
25  const int s33c = parameters.dbase.get<int >("s33c");
26 
27  const int pc = parameters.dbase.get<int >("pc");
28 
29  bool assignStress = s11c >=0 ;
30 
31  // hard code some numbers for now:
32  // real k1=1., k2=0., k3=0.;
33  // real ap=1., xa=.5, ya=0.;
34  real cp = sqrt( (lambda+2.*mu)/rho );
35  real cs = sqrt( mu/rho );
36 
37  std::vector<real> & twd = parameters.dbase.get<std::vector<real> >("travelingWaveData");
38  const int np = int(twd[0]); // number of p wave solutions
39  const int ns = int(twd[1]); // number of s wave solutions
40 
41  if( pdeVariation == SmParameters::hemp )
42  {
43  printF("\n\n **************** FIX ME: travelingWave: finish me for HEMP **********\n\n");
44  // OV_ABORT("error");
45  }
46 
47 
48  // printF("**** travelingWave: cp=%8.2e, t=%8.2e v0=%8.2e *********\n",cp,t,v0);
49  int i1,i2,i3;
50  if( mg.numberOfDimensions()==2 )
51  {
52  real z0=0.;
53  FOR_3D(i1,i2,i3,I1,I2,I3)
54  {
55  real x0 = X(i1,i2,i3,0);
56  real y0 = X(i1,i2,i3,1);
57 
58  real u1=0., u2=0.; // back-ground field
59  real v1=0., v2=0.;
60  real s11=0., s12=0., s22=0.;
61  real div=0.;
62  // real a1=0., b1=-1., a2=0., b2=0.;
63  // u1 = a1*x0 + b1*t; // more general back ground field
64  // u2 = a2*x0 + b2*t;
65 
66  int m=2;
67  for (int n=0; n<np; n++ ) // p-wave solutions
68  {
69  real ap = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
70 
71  real xi = k1*(x0-xa) + k2*(y0-ya) - cp*t;
72 
73  if( xi <= 0. )
74  {
75  u1 += -ap*k1*xi;
76  u2 += -ap*k2*xi;
77  v1 += ap*cp*k1;
78  v2 += ap*cp*k2;
79  s11+= -ap*( lambda+2.*mu*k1*k1 );
80  s12+= -ap*( 2.*mu*k1*k2 );
81  s22+= -ap*( lambda+2.*mu*k2*k2 );
82  div+= -ap;
83  }
84  }
85  for (int n=0; n<ns; n++ ) // s-wave solutions
86  {
87  real as = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
88 
89  real xi = k1*(x0-xa) + k2*(y0-ya) - cs*t;
90 
91  if( xi <= 0. )
92  { // (-k2,k1)
93  u1 += +as*k2*xi;
94  u2 += -as*k1*xi;
95  v1 += -as*cs*k2;
96  v2 += +as*cs*k1;
97  s11+= -as*(-2.*mu*k1*k2 );
98  s12+= -as*( mu*(k1*k1-k2*k2) );
99  s22+= -as*( 2.*mu*k1*k2 );
100  }
101  }
102 
103  // printF(" (i1,i2)=(%i,%i) (x0,y0)=(%8.2e,%8.2e) xi=%8.2e (u1,u2)=(%8.2e,%8.2e)\n",i1,i2,x0,y0,xi,u1,u2);
104 
105  if( evalSolution )
106  {
107  if( pdeVariation == SmParameters::hemp )
108  {
109  U(i1,i2,i3,u1c) =u1;
110  U(i1,i2,i3,u2c) =u2;
111  U(i1,i2,i3,uc) =x0;
112  U(i1,i2,i3,vc) =y0;
113 
114  state0Local(i1,i2,i3,1) = 1.0; // density
115  /*********/
116  state0Local(i1,i2,i3,0) = state0Local(i1,i2,i3,1)*det(i1,i2,i3)*mg.gridSpacing(axis1)*mg.gridSpacing(axis2); // mass
117  if( pdeVariation == SmParameters::hemp &&
118  i1 > I1Base && i1 < I1Bound &&
119  i2 > I2Base && i2 < I2Bound )
120  {
121  real area = 0.25*(det(i1,i2,i3)+det(i1+1,i2,i3)+det(i1+1,i2+1,i3)+det(i1,i2+1,i3));
122  state0Local(i1,i2,i3,0) = area*mg.gridSpacing(axis1)*mg.gridSpacing(axis2); // mass
123  }
124 
125  }
126  else
127  {
128  U(i1,i2,i3,uc) =u1;
129  U(i1,i2,i3,vc) =u2;
130  }
131 
132  if( assignVelocities )
133  {
134  U(i1,i2,i3,v1c) = v1;
135  U(i1,i2,i3,v2c) = v2;
136  }
137  if( assignStress )
138  {
139  U(i1,i2,i3,s11c) =s11;
140  U(i1,i2,i3,s12c) =s12;
141  U(i1,i2,i3,s21c) =s12;
142  U(i1,i2,i3,s22c) =s22;
143  if( pdeVariation == SmParameters::hemp )
144  {
145  real press = -(lambda+2.0*mu/3.0)*div;
146  U(i1,i2,i3,pc) = press;
147  U(i1,i2,i3,s11c) += press;
148  U(i1,i2,i3,s22c) += press;
149  }
150 
151  }
152  }
153  else
154  {
155  if( pdeVariation == SmParameters::hemp )
156  {
157  ERR(i1,i2,i3,u1c) = U(i1,i2,i3,u1c) - u1;
158  ERR(i1,i2,i3,u2c) = U(i1,i2,i3,u2c) - u2;
159  ERR(i1,i2,i3,uc) = U(i1,i2,i3,uc) - x0;
160  ERR(i1,i2,i3,vc) = U(i1,i2,i3,vc) - y0;
161  }
162  else
163  {
164  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
165  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
166  }
167 
168  if( assignVelocities )
169  {
170  ERR(i1,i2,i3,v1c) =U(i1,i2,i3,v1c) - v1;
171  ERR(i1,i2,i3,v2c) =U(i1,i2,i3,v2c) - v2;
172  }
173  if( assignStress )
174  {
175  ERR(i1,i2,i3,s11c) =U(i1,i2,i3,s11c) -s11;
176  ERR(i1,i2,i3,s12c) =U(i1,i2,i3,s12c) -s12;
177  ERR(i1,i2,i3,s21c) =U(i1,i2,i3,s21c) -s12;
178  ERR(i1,i2,i3,s22c) =U(i1,i2,i3,s22c) -s22;
179  if( pdeVariation == SmParameters::hemp )
180  {
181  real press = -(lambda+2.0*mu/3.0)*div;
182  ERR(i1,i2,i3,s11c) -= press;
183  ERR(i1,i2,i3,s22c) -= press;
184  ERR(i1,i2,i3,pc) = U(i1,i2,i3,pc) - press;
185  }
186 
187  }
188 
189  }
190  } // end FOR_3D
191 
192  }
193  else
194  {
195  OV_ABORT("Error: finish me");
196  }
197 }
198 
199 #endMacro
200 
201 // ==================================================================================
202 // Macro to evaluate the traveling (sine) wave solution
203 //
204 // evalSolution : true=eval solution, false=eval error
205 // ==================================================================================
206 #beginMacro getPlaneTravelingWaveSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
207 // printF("INFO: The traveling wave solutions are combinations of p and s solutions:\n"
208 // " [u1,u2] = ap* [ k1,k2] * G( k1*(x-xa)+k2*(y-ya) - cp*t ) (p-wave)\n"
209 // " [u1,u2] = as* [-k2,k1] * G( k1*(x-xa)+k2*(y-ya) - cs*t ) (s-wave\n",
210 // " where G(xi)=sin(freq,xi) \n");
211 {
212  const int v1c = parameters.dbase.get<int >("v1c");
213  const int v2c = parameters.dbase.get<int >("v2c");
214  const int v3c = parameters.dbase.get<int >("v3c");
215 
216  bool assignVelocities= v1c>=0 ;
217  const int s11c = parameters.dbase.get<int >("s11c");
218  const int s12c = parameters.dbase.get<int >("s12c");
219  const int s13c = parameters.dbase.get<int >("s13c");
220  const int s21c = parameters.dbase.get<int >("s21c");
221  const int s22c = parameters.dbase.get<int >("s22c");
222  const int s23c = parameters.dbase.get<int >("s23c");
223  const int s31c = parameters.dbase.get<int >("s31c");
224  const int s32c = parameters.dbase.get<int >("s32c");
225  const int s33c = parameters.dbase.get<int >("s33c");
226 
227  const int pc = parameters.dbase.get<int >("pc");
228 
229  bool assignStress = s11c >=0 ;
230 
231  // hard code some numbers for now:
232  // real k1=1., k2=0., k3=0.;
233  // real ap=1., xa=.5, ya=0.;
234  real cp = sqrt( (lambda+2.*mu)/rho );
235  real cs = sqrt( mu/rho );
236 
237  std::vector<real> & twd = parameters.dbase.get<std::vector<real> >("travelingWaveData");
238  const int np = int(twd[0]); // number of p wave solutions
239  const int ns = int(twd[1]); // number of s wave solutions
240 
241  if( pdeVariation == SmParameters::hemp )
242  {
243  printF("\n\n **************** FIX ME: travelingWave: finish me for HEMP **********\n\n");
244  // OV_ABORT("error");
245  }
246 
247 
248  // printF("**** planeTravelingWave: cp=%8.2e, t=%8.2e *********\n",cp,t);
249 
250  int i1,i2,i3;
251  if( mg.numberOfDimensions()==2 )
252  {
253  real z0=0.;
254  FOR_3D(i1,i2,i3,I1,I2,I3)
255  {
256  real x0 = X(i1,i2,i3,0);
257  real y0 = X(i1,i2,i3,1);
258 
259  real u1=0., u2=0.; // back-ground field
260  real v1=0., v2=0.;
261  real s11=0., s12=0., s22=0.;
262  real div=0.;
263  // real a1=0., b1=-1., a2=0., b2=0.;
264  // u1 = a1*x0 + b1*t; // more general back ground field
265  // u2 = a2*x0 + b2*t;
266 
267  int m=2;
268  for (int n=0; n<np; n++ ) // p-wave solutions
269  {
270  real ap = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
271  real freq=twd[m++];
272 
273  real xi = k1*(x0-xa) + k2*(y0-ya) - cp*t;
274  real sinf = sin(freq*xi), cosf=cos(freq*xi);
275 
276  u1 += -ap*k1*sinf;
277  u2 += -ap*k2*sinf;
278  v1 += ap*cp*k1*freq*cosf;
279  v2 += ap*cp*k2*freq*cosf;
280  s11+= -ap*( lambda+2.*mu*k1*k1 )*freq*cosf;
281  s12+= -ap*( 2.*mu*k1*k2 )*freq*cosf;
282  s22+= -ap*( lambda+2.*mu*k2*k2 )*freq*cosf;
283  div+= -ap;
284  }
285  for (int n=0; n<ns; n++ ) // s-wave solutions
286  {
287  real as = twd[m++], k1=twd[m++], k2=twd[m++], k3=twd[m++], xa=twd[m++], ya=twd[m++], za=twd[m++];
288 
289  real xi = k1*(x0-xa) + k2*(y0-ya) - cs*t;
290  real freq=twd[m++];
291  real sinf = sin(freq*xi), cosf=cos(freq*xi);
292 
293  // (-k2,k1)
294  u1 += +as*k2*sinf;
295  u2 += -as*k1*sinf;
296  v1 += -as*cs*k2*freq*cosf;
297  v2 += +as*cs*k1*freq*cosf;
298  s11+= -as*(-2.*mu*k1*k2 )*freq*cosf;
299  s12+= -as*( mu*(k1*k1-k2*k2) )*freq*cosf;
300  s22+= -as*( 2.*mu*k1*k2 )*freq*cosf;
301 
302  }
303 
304  // printF(" (i1,i2)=(%i,%i) (x0,y0)=(%8.2e,%8.2e) xi=%8.2e (u1,u2)=(%8.2e,%8.2e)\n",i1,i2,x0,y0,xi,u1,u2);
305 
306  if( evalSolution )
307  {
308  if( pdeVariation == SmParameters::hemp )
309  {
310  U(i1,i2,i3,u1c) =u1;
311  U(i1,i2,i3,u2c) =u2;
312  U(i1,i2,i3,uc) =x0;
313  U(i1,i2,i3,vc) =y0;
314 
315  state0Local(i1,i2,i3,1) = 1.0; // density
316  /*********/
317  state0Local(i1,i2,i3,0) = state0Local(i1,i2,i3,1)*det(i1,i2,i3)*mg.gridSpacing(axis1)*mg.gridSpacing(axis2); // mass
318  if( pdeVariation == SmParameters::hemp &&
319  i1 > I1Base && i1 < I1Bound &&
320  i2 > I2Base && i2 < I2Bound )
321  {
322  real area = 0.25*(det(i1,i2,i3)+det(i1+1,i2,i3)+det(i1+1,i2+1,i3)+det(i1,i2+1,i3));
323  state0Local(i1,i2,i3,0) = area*mg.gridSpacing(axis1)*mg.gridSpacing(axis2); // mass
324  }
325 
326  }
327  else
328  {
329  U(i1,i2,i3,uc) =u1;
330  U(i1,i2,i3,vc) =u2;
331  }
332 
333  if( assignVelocities )
334  {
335  U(i1,i2,i3,v1c) = v1;
336  U(i1,i2,i3,v2c) = v2;
337  }
338  if( assignStress )
339  {
340  U(i1,i2,i3,s11c) =s11;
341  U(i1,i2,i3,s12c) =s12;
342  U(i1,i2,i3,s21c) =s12;
343  U(i1,i2,i3,s22c) =s22;
344  if( pdeVariation == SmParameters::hemp )
345  {
346  real press = -(lambda+2.0*mu/3.0)*div;
347  U(i1,i2,i3,pc) = press;
348  U(i1,i2,i3,s11c) += press;
349  U(i1,i2,i3,s22c) += press;
350  }
351 
352  }
353  }
354  else
355  {
356  if( pdeVariation == SmParameters::hemp )
357  {
358  ERR(i1,i2,i3,u1c) = U(i1,i2,i3,u1c) - u1;
359  ERR(i1,i2,i3,u2c) = U(i1,i2,i3,u2c) - u2;
360  ERR(i1,i2,i3,uc) = U(i1,i2,i3,uc) - x0;
361  ERR(i1,i2,i3,vc) = U(i1,i2,i3,vc) - y0;
362  }
363  else
364  {
365  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
366  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
367  }
368 
369  if( assignVelocities )
370  {
371  ERR(i1,i2,i3,v1c) =U(i1,i2,i3,v1c) - v1;
372  ERR(i1,i2,i3,v2c) =U(i1,i2,i3,v2c) - v2;
373  }
374  if( assignStress )
375  {
376  ERR(i1,i2,i3,s11c) =U(i1,i2,i3,s11c) -s11;
377  ERR(i1,i2,i3,s12c) =U(i1,i2,i3,s12c) -s12;
378  ERR(i1,i2,i3,s21c) =U(i1,i2,i3,s21c) -s12;
379  ERR(i1,i2,i3,s22c) =U(i1,i2,i3,s22c) -s22;
380  if( pdeVariation == SmParameters::hemp )
381  {
382  real press = -(lambda+2.0*mu/3.0)*div;
383  ERR(i1,i2,i3,s11c) -= press;
384  ERR(i1,i2,i3,s22c) -= press;
385  ERR(i1,i2,i3,pc) = U(i1,i2,i3,pc) - press;
386  }
387 
388  }
389 
390  }
391  } // end FOR_3D
392 
393  }
394  else
395  {
396  OV_ABORT("Error: finish me");
397  }
398 }
399 
400 #endMacro
401 
402 
403 // ==========================================================================
404 // Define a Rayleigh surface wave: (see cgDoc/sm/notes.pdf)
405 // u1 = SUM_n a_n [ exp(-b1(k_n)*y + ... ] cos( 2*pi*k_n (x-c*t) )
406 // u2 = SUM_n a_n [ ] sin( 2*pi*k_n (x-c*t) )
407 //
408 // Here we assume that the solid occupies the space y<= ySurf
409 // where ySurf is given by the user.
410 // ==========================================================================
411 #beginMacro getRayleighWaveSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
412 {
413  const int v1c = parameters.dbase.get<int >("v1c");
414  const int v2c = parameters.dbase.get<int >("v2c");
415  const int v3c = parameters.dbase.get<int >("v3c");
416 
417  bool assignVelocities= v1c>=0 ;
418  const int s11c = parameters.dbase.get<int >("s11c");
419  const int s12c = parameters.dbase.get<int >("s12c");
420  const int s13c = parameters.dbase.get<int >("s13c");
421  const int s21c = parameters.dbase.get<int >("s21c");
422  const int s22c = parameters.dbase.get<int >("s22c");
423  const int s23c = parameters.dbase.get<int >("s23c");
424  const int s31c = parameters.dbase.get<int >("s31c");
425  const int s32c = parameters.dbase.get<int >("s32c");
426  const int s33c = parameters.dbase.get<int >("s33c");
427 
428  const int pc = parameters.dbase.get<int >("pc");
429 
430  bool assignStress = s11c >=0 ;
431 
432  real cp = sqrt( (lambda+2.*mu)/rho );
433  real cs = sqrt( mu/rho );
434 
435  std::vector<real> & data = parameters.dbase.get<std::vector<real> >("RayleighWaveData");
436  const int nk = int(data[0]); // number of modes
437  const real cr = data[1]; // Rayleigh wave speed
438  const real ySurf = data[2]; // y value on surface
439  const real period= data[3]; // Rayleigh wave speed
440  const real xShift= data[4]; // shift in x for wave
441  const int mStart=5; // Fourier coeff's start at this index in the data
442 
443  if( pdeVariation == SmParameters::hemp )
444  {
445  printF("\n\n **************** FIX ME: RayelighWave: finish me for HEMP **********\n\n");
446  OV_ABORT("error");
447  }
448 
449  real cb1 = sqrt(1.-SQR(cr/cp)); // b1/k : for computing b1
450  real cb2 = sqrt(1.-SQR(cr/cs)); // b2/k : for computing b2
451 
452  real c1 = .5*SQR(cr/cs)-1.; // x/2-1 , x=cr^2/cs^2
453 
454  if( t==0. )
455  {
456  printF("**** RayleighWave: ySurf=%8.2e, cr=%8.2e, period=%8.2e, t=%8.2e *********\n",ySurf,cr,period,t);
457  int m=mStart;
458  for( int n=0; n<nk; n++ )
459  {
460  real k = data[m++]; // k=wave-number
461  real a=data[m++]; // an : amplitude
462  real b=data[m++];
463  printF(" k%i = %e, a%i=%e, b%i=%e\n",n,k,n,a,n,b);
464  }
465 
466  }
467 
468  real scale = -1./( cb1+(c1/cb2) ); // make coeff of cos() in u2 = a at y=0
469 
470 
471  int i1,i2,i3;
472  if( mg.numberOfDimensions()==2 )
473  {
474  real z0=0.;
475  FOR_3D(i1,i2,i3,I1,I2,I3)
476  {
477  real x0 = X(i1,i2,i3,0)-xShift;
478  real y0 = X(i1,i2,i3,1)-ySurf;
479 
480  real u1=0., u2=0.;
481  real v1=0., v2=0.;
482  real s11=0., s12=0., s22=0.;
483 
484  int m=mStart;
485  // --- loop over different values of k and add contributions ---
486  for( int n=0; n<nk; n++ )
487  {
488  real k = twoPi*data[m++]/period; // 2*pi*k/period , k=wave-number
489  // -- note definition of a and b so that they define the Fourier coefficients
490  // of u2 on the surface
491  real b=data[m++]*scale;
492  real a=data[m++]*scale;
493 
494  real b1= k*cb1, b2=k*cb2;
495 
496  real eb1 = exp(b1*(y0)), eb2 = exp(b2*(y0));
497 
498  real ct = cos(k*(x0-cr*t));
499  real st = sin(k*(x0-cr*t));
500 
501  u1 += ( eb1 + c1*eb2 )*( a*ct+b*st);
502  u2 += ( cb1*eb1 + (c1/cb2)*eb2 )*( a*st-b*ct);
503  if( assignVelocities )
504  {
505  v1 += ( eb1 + c1*eb2 )*( k*cr*( a*st-b*ct) );
506  v2 +=-( cb1*eb1 + (c1/cb2)*eb2 )*( k*cr*( a*ct+b*st) );
507  }
508  if( assignStress )
509  {
510  real u1x = ( eb1 + c1*eb2 )*( k*(-a*st+b*ct) );
511  real u1y = ( b1*eb1 + b2*c1*eb2 )*( ( a*ct+b*st) );
512 
513  real u2x = ( cb1*eb1 + (c1/cb2)*eb2 )*( k*(a*ct+b*st) );
514  real u2y = ( b1*cb1*eb1 + b2*(c1/cb2)*eb2 )*( (a*st-b*ct) );
515  real div=u1x+u2y;
516 
517  s11 += lambda*div+2.*mu*u1x;
518  s12 += mu*( u1y+u2x );
519  s22 += lambda*div+2.*mu*u2y;
520  }
521 
522  }
523 
524  // printF(" (i1,i2)=(%i,%i) (x0,y0)=(%8.2e,%8.2e) xi=%8.2e (u1,u2)=(%8.2e,%8.2e)\n",i1,i2,x0,y0,xi,u1,u2);
525 
526  if( evalSolution )
527  {
528 
529  U(i1,i2,i3,uc) =u1;
530  U(i1,i2,i3,vc) =u2;
531 
532  if( assignVelocities )
533  {
534  U(i1,i2,i3,v1c) = v1;
535  U(i1,i2,i3,v2c) = v2;
536  }
537  if( assignStress )
538  {
539  U(i1,i2,i3,s11c) =s11;
540  U(i1,i2,i3,s12c) =s12;
541  U(i1,i2,i3,s21c) =s12;
542  U(i1,i2,i3,s22c) =s22;
543  }
544  }
545  else
546  {
547  if( pdeVariation == SmParameters::hemp )
548  {
549  ERR(i1,i2,i3,u1c) = U(i1,i2,i3,u1c) - u1;
550  ERR(i1,i2,i3,u2c) = U(i1,i2,i3,u2c) - u2;
551  ERR(i1,i2,i3,uc) = U(i1,i2,i3,uc) - x0;
552  ERR(i1,i2,i3,vc) = U(i1,i2,i3,vc) - y0;
553  }
554  else
555  {
556  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
557  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
558  }
559 
560  if( assignVelocities )
561  {
562  ERR(i1,i2,i3,v1c) =U(i1,i2,i3,v1c) - v1;
563  ERR(i1,i2,i3,v2c) =U(i1,i2,i3,v2c) - v2;
564  }
565  if( assignStress )
566  {
567  ERR(i1,i2,i3,s11c) =U(i1,i2,i3,s11c) -s11;
568  ERR(i1,i2,i3,s12c) =U(i1,i2,i3,s12c) -s12;
569  ERR(i1,i2,i3,s21c) =U(i1,i2,i3,s21c) -s12;
570  ERR(i1,i2,i3,s22c) =U(i1,i2,i3,s22c) -s22;
571 
572  }
573 
574  }
575  } // end FOR_3D
576 
577  }
578  else
579  {
580  OV_ABORT("RayleighWave:ERROR: finish me for 3D");
581  }
582 }
583 
584 #endMacro
585 
586 // =======================================================================================
587 // The function "fg" is basically the integral appearing in the D'Alambert solution
588 // ======================================================================================
589 #defineMacro fg(x,z) ( .5*( cg1*(x)*( 1.+ (z)*( 7./p + (z)*( 21./(2.*p-1.) + (z)*( 35./(3.*p-2.) + (z)*( 35./(4.*p-3.) \
590  + (z)*( 21./(5.*p-4.) + (z)*( 7./(6.*p-5.) + z/(7.*p-6.) ))))))) ) )
591 
592 
593 // ===========================================================
594 // Evaluate the D'Alambert function "f" and it's derivative
595 // Here we assume that u(x,0)=0 and v(x,0)!=0
596 // ===========================================================
597 #beginMacro getF( x, f, fPrime )
598 {
599  real xx = -x/cp;
600  real xp1 = pow(xx,p-1);
601  real z=cg2*xp1;
602  f = cp*fg(xx,z) - .5*(a/p)*xp1*xx;
603  fPrime = .5*( -cg1*pow( 1. + z , 7.) + a*xp1/cp );
604 }
605 #endMacro
606 
607 // ===========================================================
608 // Evaluate the D'Alambert function "g"
609 // Here we assume that u(x,0)=0 and v(x,0)!=0
610 // ===========================================================
611 #beginMacro getG( x, g, gPrime )
612 {
613  if( x<=0. )
614  {
615  real xx = -x/cp;
616  real xp1 = pow(xx,p-1);
617  real z=cg2*xp1;
618  g = cp*fg(xx,z) + .5*(a/p)*xp1*xx ;
619  gPrime = +.5*( -cg1*pow( 1. + z , 7.) -a*xp1/cp );
620  }
621  else
622  {
623  // g(x) = F(x/cp) - f(-x),
624  real xx=x/cp;
625  real xp1 = pow(xx,p-1);
626  real z=cg2*xp1;
627  g = -(a/p)*xp1*xx - cp*fg(xx,z) + .5*(a/p)*xp1*xx;
628  gPrime = -a*xp1/cp + .5*( -cg1*pow( 1. + z , 7.) + a*xp1/cp );
629  }
630 }
631 #endMacro
632 
633 
634 
635 // ==========================================================================
636 // Define the pistonMotion solution (see cgDoc/mp/fluidStructure/fsm.tex)
637 // ==========================================================================
638 #beginMacro getPistonMotionSolution(evalSolution,U,ERR,X,t,I1,I2,I3)
639 {
640  const int v1c = parameters.dbase.get<int >("v1c");
641  const int v2c = parameters.dbase.get<int >("v2c");
642  const int v3c = parameters.dbase.get<int >("v3c");
643 
644  bool assignVelocities= v1c>=0 ;
645  const int s11c = parameters.dbase.get<int >("s11c");
646  const int s12c = parameters.dbase.get<int >("s12c");
647  const int s13c = parameters.dbase.get<int >("s13c");
648  const int s21c = parameters.dbase.get<int >("s21c");
649  const int s22c = parameters.dbase.get<int >("s22c");
650  const int s23c = parameters.dbase.get<int >("s23c");
651  const int s31c = parameters.dbase.get<int >("s31c");
652  const int s32c = parameters.dbase.get<int >("s32c");
653  const int s33c = parameters.dbase.get<int >("s33c");
654 
655  const int pc = parameters.dbase.get<int >("pc");
656 
657  bool assignStress = s11c >=0 ;
658 
659  const real cp = sqrt( (lambda+2.*mu)/rho );
660  const real cs = sqrt( mu/rho );
661 
662  std::vector<real> & data = parameters.dbase.get<std::vector<real> >("pistonMotionData");
663 
664  int m=0;
665  const real a =data[m++];
666  const real p =data[m++];
667  const real rhog =data[m++];
668  const real pg =data[m++];
669  const real gamma=data[m++];
670  real angle =data[m++]; // angle (in degrees) for a rotated piston
671 
672  const real a0 = sqrt( gamma*pg/rhog); // speed of sound in the gas
673 
674  if( pdeVariation == SmParameters::hemp )
675  {
676  printF("\n\n **************** FIX ME: getPistonMotionSolution: finish me for HEMP **********\n\n");
677  OV_ABORT("error");
678  }
679 
680 
681 
682 
683  if( t==0. )
684  {
685  printP("**** getPistonMotion: a=%8.2e, p=%8.2e, Gas: rho=%8.2e, p=%8.2e gamma=%8.2e angle=%5.2f(degrees)*******\n",
686  a,p,rhog,pg,gamma,angle);
687  }
688 
689  angle = angle*Pi/180.;
690  const real cosa = cos(angle), sina=sin(angle);
691 
692  const real cg1 = pg/(rho*cp*cp);
693  const real cg2 = (-a)*(gamma-1.)/(2.*a0);
694  // we assume gamma=1.4
695  assert( fabs(gamma-1.4) < REAL_EPSILON*100. );
696 
697  int i1,i2,i3;
698  if( mg.numberOfDimensions()==2 )
699  {
700  real z0=0.;
701  FOR_3D(i1,i2,i3,I1,I2,I3)
702  {
703  real x0 = X(i1,i2,i3,0);
704  real y0 = X(i1,i2,i3,1);
705 
706  real u1=0., u2=0.;
707  real v1=0., v2=0.;
708  real s11=0., s12=0., s22=0.;
709 
710  // Boundary motion:
711  // F(t) = -(a/p)*t^p
712  // F'(t) = -a*t^{p-1}
713  // Solution:
714  // u1 = f(x-cp*t) + g(x+cp*t)
715  // where
716  // f(x) = - 1/(2*cp)*( int_0^x v0(s) ds ) = -(1/2)* int_0^x G(-s/cp) ds , for x<0
717  // = (cp/2)* int_0^{-x/cp} G(u) du
718  // g(x) = + 1/(2*cp)*( int_0^x v0(s) ds ) = (1/2)* int_0^x G(-s/cp) ds , for x<0
719  // = -(cp/2)* int_0^{-x/cp} G(u) du
720  // g(x) = F(x/cp) - f(-x), for x>0
721  //
722  // v0(s) = cp*G(-t/cp) : velocity at t=0 for s<0
723  // G(t) = (pg/(rho*cp^2)) * [ 1 + (gamma-1)/(2*a0)* F'(t) ]^7 + F'(t)/cp
724  // = cg1* [ 1 + cg2*t^{p-1} ]^7 + F'(t)/cp
725  // cg1=p0/(rho*cp^2), cg2=(-a)*(gamma-1)/(2*a0)
726  //
727  // where we have assumed that gamma=1.4=7/5 so that 2*gamma/(gamma-1)= 7
728  //
729  // Int_0^t G(s) ds = cg1*[ t + (7/p)*cg2*t^{p-1}*t + (21/(2p-1))*(cg2*t^{p-1})^2*t + ) + F(t)/cp
730  // = cg1*t*[ 1 + Z*(7)/(p) + Z^2*(21)/(2p-1) + Z^3*(35)/(3p-2) + Z^4*(35)/(4p-3) + ...
731  // + Z^5*(21)/(5p-4) + Z^6*(7)/(6p-5) + Z^7*(1)/(7p-6) ]
732  // Z=cg2*t^{p-1}
733 
734  // xa = distance of point (x0,y0) from the plane (cosa,sina).(x,y)=0
735  real xa = x0*cosa + y0*sina;
736 
737  real xp = xa + cp*t;
738  real xm = xa - cp*t;
739 
740 
741  real fm, fmPrime, gp, gpPrime;
742 
743  getF( xm, fm, fmPrime );
744  getG( xp, gp, gpPrime );
745 
746  real ua = fm + gp;
747  real va = cp*( - fmPrime + gpPrime );
748  real uap=fmPrime + gpPrime;
749 
750  // Piston at an angle:
751  // n = [ cos(angle) , sin(angle) ] = [ c, s ] -- normal to the face
752  // u1 = c*ua, u2=s*ua
753  // u1.x = c*ua.x, u1.y = c*ua.y, u2.x = s*ua.x, u2.y = s*ua.y
754  //
755  // (1) ua.n = c*ua.x + s*ua.y = uap
756  // (2) ua.t =-s*ua.x + c*ua.y = 0 "tangential derivative of motion"
757  //
758  // (1) and (2) --> ua.x = c*uap, ua.y=s*uap
759  // Thus: u1.x = c*ua.x = c*c*uap, u1.y = c*ua.y = c*s*uap
760 
761  u1 = ua*cosa;
762  u2 = ua*sina;
763  v1 = va*cosa;
764  v2 = va*sina;
765 
766  real u1x = uap*cosa*cosa;
767  real u1y = uap*cosa*sina;
768 
769  real u2x = uap*sina*cosa;
770  real u2y = uap*sina*sina;
771 
772 
773  s11 = (lambda+2.*mu)*u1x + lambda*u2y;
774  s12 = mu*( u1y+u2x );
775  s22 = lambda*u1x + (lambda+2.*mu)*u2y;
776 
777  // printF("piston (i1,i2)=(%i,%i) (u1,u2)=(%8.2e,%8.2e)\n",i1,i2,u1,u2);
778 
779  if( evalSolution )
780  {
781  U(i1,i2,i3,uc) =u1;
782  U(i1,i2,i3,vc) =u2;
783 
784  if( assignVelocities )
785  {
786  U(i1,i2,i3,v1c) = v1;
787  U(i1,i2,i3,v2c) = v2;
788  }
789  if( assignStress )
790  {
791  U(i1,i2,i3,s11c) =s11;
792  U(i1,i2,i3,s12c) =s12;
793  U(i1,i2,i3,s21c) =s12;
794  U(i1,i2,i3,s22c) =s22;
795  }
796  }
797  else
798  {
799  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
800  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
801 
802  if( assignVelocities )
803  {
804  ERR(i1,i2,i3,v1c) =U(i1,i2,i3,v1c) - v1;
805  ERR(i1,i2,i3,v2c) =U(i1,i2,i3,v2c) - v2;
806  }
807  if( assignStress )
808  {
809  ERR(i1,i2,i3,s11c) =U(i1,i2,i3,s11c) -s11;
810  ERR(i1,i2,i3,s12c) =U(i1,i2,i3,s12c) -s12;
811  ERR(i1,i2,i3,s21c) =U(i1,i2,i3,s21c) -s12;
812  ERR(i1,i2,i3,s22c) =U(i1,i2,i3,s22c) -s22;
813 
814  }
815 
816  }
817  } // end FOR_3D
818 
819  }
820  else
821  {
822  OV_ABORT("getPistonMotion:ERROR: finish me for 3D");
823  }
824 }
825 
826 #endMacro
827 
828 
829 
830 
831 
832 
833 
834 
835 // (Ex).t = (1/eps)*[ (Hz).y ]
836 // (Ey).t = (1/eps)*[ -(Hz).x ]
837 // (Hz).t = (1/mu) *[ (Ex).y - (Ey).x ]
838 
839 #define exTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(-ky/(eps*cc))
840 #define eyTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*( kx/(eps*cc))
841 #define hzTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))
842 
843 #define exLaplacianTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(+ky*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
844 #define eyLaplacianTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(-kx*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
845 #define hzLaplacianTrue(x,y,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*( -(twoPi*twoPi*(kx*kx+ky*ky) ) )
846 
847 // define eyTrue(x,y,t) exp( -beta*SQR((x-x0)-c*(t)) )
848 
849 // Here is a plane wave with the shape of a Gaussian
850 // xi = kx*(x)+ky*(y)-cc*(t)
851 // cc= c*sqrt( kx*kx+ky*ky );
852 #define hzGaussianPulse(xi) exp(-betaGaussianPlaneWave*((xi)*(xi)))
853 #define exGaussianPulse(xi) hzGaussianPulse(xi)*(-ky/(eps*cc))
854 #define eyGaussianPulse(xi) hzGaussianPulse(xi)*( kx/(eps*cc))
855 
856 #define hzLaplacianGaussianPulse(xi) ((4.*betaGaussianPlaneWave*betaGaussianPlaneWave*(kx*kx+ky*ky))*xi*xi-\
857  (2.*betaGaussianPlaneWave*(kx*kx+ky*ky)))*exp(-betaGaussianPlaneWave*((xi)*(xi)))
858 #define exLaplacianGaussianPulse(xi) hzLaplacianGaussianPulse(xi,t)*(-ky/(eps*cc))
859 #define eyLaplacianGaussianPulse(xi) hzLaplacianGaussianPulse(xi,t)*( kx/(eps*cc))
860 
861 // 3D
862 // E:
863 // u.tt = (1/eps)*[ ((1/mu)*u.x).x + ((1/mu)*u.y).y + ((1/mu)*u.z).z ]
864 // div(u)=0
865 // H
866 // v.tt = (1/mu)*[ ((1/eps)*v.x).x + ((1/eps)*v.y).y + ((1/eps)*v.z).z ]
867 // Define macros for forcing functions
868 
869 
870 //
871 // (Ex).t = (1/eps)*[ (Hz).y - (Hy).z ]
872 // (Ey).t = (1/eps)*[ (Hx).z - (Hz).x ]
873 // (Ez).t = (1/eps)*[ (Hy).x - (Hx).y ]
874 // (Hx).t = (1/mu) *[ (Ey).z - (Ez).y ]
875 // (Hy).t = (1/mu) *[ (Ez).x - (Ex).z ]
876 // (Hz).t = (1/mu) *[ (Ex).y - (Ey).x ]
877 
878 // ****************** finish this -> should `rotate' the 2d solution ****************
879 
880 #define exTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)+kz*(z)-cc*(t)))*(-ky/(eps*cc))
881 #define eyTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)+kz*(z)-cc*(t)))*( kx/(eps*cc))
882 #define ezTrue3d(x,y,z,t) 0
883 
884 #define hxTrue3d(x,y,z,t) 0
885 #define hyTrue3d(x,y,z,t) 0
886 #define hzTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)+kz*(z)-cc*(t)))
887 
888 #define exLaplacianTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(+ky*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
889 #define eyLaplacianTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*(-kx*(twoPi*twoPi*(kx*kx+ky*ky))/(eps*cc))
890 #define ezLaplacianTrue3d(x,y,z,t) 0
891 
892 #define hxLaplacianTrue3d(x,y,z,t) 0
893 #define hyLaplacianTrue3d(x,y,z,t) 0
894 #define hzLaplacianTrue3d(x,y,z,t) sin(twoPi*(kx*(x)+ky*(y)-cc*(t)))*( -(twoPi*twoPi*(kx*kx+ky*ky) ) )
895 
896 
897 // ------------ macros for the plane material interface -------------------------
898 #defineMacro PMIex(x,y,t) (aa1*cos(twoPi*(kx*(x)+ky*(y)-w*(t))) + r*bb1*cos(twoPi*(-kx*(x)+ky*(y)-w*(t))))
899 #defineMacro PMIey(x,y,t) (aa2*cos(twoPi*(kx*(x)+ky*(y)-w*(t))) + r*bb2*cos(twoPi*(-kx*(x)+ky*(y)-w*(t))))
900 #defineMacro PMIhz(x,y,t) ((-twoPi/w)*( (a1*ky-a2*kx)*sin(twoPi*(kx*(x)+ky*(y)-w*(t))) \
901  + (b1*ky+b2*kx)*r*sin(twoPi*(-kx*(x)+ky*(y)-w*(t)))) )
902 
903 #defineMacro PMITex(x,y,t) (tau*dd1*cos(twoPi*(kappax*(x)+kappay*(y)-w*(t))))
904 #defineMacro PMITey(x,y,t) (tau*dd2*cos(twoPi*(kappax*(x)+kappay*(y)-w*(t))))
905 #defineMacro PMIThz(x,y,t) ((-twoPi*tau/w)*(d1*kappay-d2*kappax)*sin(twoPi*(kappax*(x)+kappay*(y)-w*(t))))
906 
907 // OPTION: initialCondition, error, boundaryCondition
908 #beginMacro setPlaneMaterialInterfaceMacro(OPTION,J1,J2,J3)
909 if( numberOfDimensions==2 )
910 {
911  int i1,i2,i3;
912 
913  // we assume there are just two grids with an interface in-between
914  const real c1=cGrid(0);
915  const real c2=cGrid(1);
916 
917  const real kNorm = sqrt(kx*kx+ky*ky);
918  const real a1=-ky/kNorm, a2= kx/kNorm;
919  const real b1=-ky/kNorm, b2=-kx/kNorm;
920  const real w = c1*kNorm;
921 
922  real kappax,kappay,kappaNorm;
923 
924  const real cr=c2/c1;
925  kappay=ky;
926  // discrimant is assumed positive (otherwise internal reflection)
927  kappax=sqrt( (kx*kx+ky*ky)/(cr*cr) - kappay*kappay );
928 
929  kappaNorm=sqrt(kappax*kappax+kappay*kappay);
930  const real d1=-kappay/kappaNorm;
931  const real d2= kappax/kappaNorm;
932 
933  const real cosTheta1=kx/kNorm;
934  const real cosTheta2=kappax/kappaNorm;
935 
936  // reflection and transmission coefficients
937  const real r = (c1*cosTheta1-c2*cosTheta2)/(c1*cosTheta1+c2*cosTheta2);
938  const real tau = (2.*c2*cosTheta1)/(c1*cosTheta1+c2*cosTheta2);
939 
940  real tm=t-dt,x,y;
941  const real x0=x0PlaneMaterialInterface[0], y0=x0PlaneMaterialInterface[1];
942  // The solution is rotated according to the normal of the plane of the material interface
943  real an0=normalPlaneMaterialInterface[0], an1=normalPlaneMaterialInterface[1];
944  real aNorm = sqrt(an0*an0+an1*an1);
945  an0/= aNorm; an1/=aNorm;
946  real aa1=an0*a1-an1*a2, aa2= an1*a1+an0*a2;
947  real bb1=an0*b1-an1*b2, bb2= an1*b1+an0*b2;
948  real dd1=an0*d1-an1*d2, dd2= an1*d1+an0*d2;
949 
950  if( grid==0 )
951  { // incident plus reflected wave.
952  FOR_3D(i1,i2,i3,J1,J2,J3)
953  {
954  real xx1 = XEP(i1,i2,i3,0) -x0;
955  real yy1 = XEP(i1,i2,i3,1) -y0;
956  // rotate the coordinate system
957  x = an0*xx1+an1*yy1;
958  y =-an1*xx1+an0*yy1;
959 
960  real u1 = PMIex(x,y,t);
961  real u2 = PMIey(x,y,t);
962  real u3 = PMIhz(x,y,t);
963 
964  #If #OPTION eq "initialCondition"
965  UEX(i1,i2,i3)= u1;
966  UEY(i1,i2,i3)= u2;
967  UHZ(i1,i2,i3)= u3;
968 
969  UMEX(i1,i2,i3)= PMIex(x,y,tm);
970  UMEY(i1,i2,i3)= PMIey(x,y,tm);
971  UMHZ(i1,i2,i3)= PMIhz(x,y,tm);
972  #Elif #OPTION eq "error"
973  ERREX(i1,i2,i3)=UEX(i1,i2,i3)-u1;
974  ERREY(i1,i2,i3)=UEY(i1,i2,i3)-u2;
975  ERRHZ(i1,i2,i3)=UHZ(i1,i2,i3)-u3;
976  #Elif #OPTION eq "boundaryCondition"
977  U(i1,i2,i3,ex)= u1;
978  U(i1,i2,i3,ey)= u2;
979  U(i1,i2,i3,hz)= u3;
980  #End
981  }
982  }
983  else
984  {
985  // transmitted wave
986  FOR_3D(i1,i2,i3,J1,J2,J3)
987  {
988  real xx1 = XEP(i1,i2,i3,0) -x0;
989  real yy1 = XEP(i1,i2,i3,1) -y0;
990  // rotate the coordinate system
991  x = an0*xx1+an1*yy1;
992  y =-an1*xx1+an0*yy1;
993 
994  real u1 = PMITex(x,y,t);
995  real u2 = PMITey(x,y,t);
996  real u3 = PMIThz(x,y,t);
997 
998  #If #OPTION eq "initialCondition"
999  UEX(i1,i2,i3)= u1;
1000  UEY(i1,i2,i3)= u2;
1001  UHZ(i1,i2,i3)= u3;
1002 
1003  UMEX(i1,i2,i3)= PMITex(x,y,tm);
1004  UMEY(i1,i2,i3)= PMITey(x,y,tm);
1005  UMHZ(i1,i2,i3)= PMIThz(x,y,tm);
1006  #Elif #OPTION eq "error"
1007  ERREX(i1,i2,i3)=UEX(i1,i2,i3)-u1;
1008  ERREY(i1,i2,i3)=UEY(i1,i2,i3)-u2;
1009  ERRHZ(i1,i2,i3)=UHZ(i1,i2,i3)-u3;
1010  #Elif #OPTION eq "boundaryCondition"
1011  U(i1,i2,i3,ex)= u1;
1012  U(i1,i2,i3,ey)= u2;
1013  U(i1,i2,i3,hz)= u3;
1014  #End
1015  }
1016  }
1017  }
1018 #endMacro
1019 
1020 //==================================================================================================
1021 // Evaluate Tom Hagstom's exact solution defined as an integral of Guassian sources
1022 //
1023 // OPTION: OPTION=solution or OPTION=error OPTION=bounary to compute the solution or the error or
1024 // the boundary condition
1025 //
1026 //==================================================================================================
1027 #beginMacro getGaussianIntegralSolution(OPTION,VEX,VEY,VHZ,t)
1028 
1029 if( initialConditionOption==gaussianIntegralInitialCondition )
1030 {
1031 
1032  double wt,wx,wy;
1033  const int nsources=1;
1034  double xs[nsources], ys[nsources], tau[nsources], var[nsources], amp[nsources];
1035  xs[0]=0.;
1036  ys[0]=1.e-8*1./3.; // should not be on a grid point
1037  tau[0]=-.95;
1038  var[0]=30.;
1039  amp[0]=1.;
1040 
1041  double period= 1.; // period in y
1042  double time=t;
1043 
1044  int i1,i2,i3;
1045 
1046  FOR_3D(i1,i2,i3,J1,J2,J3)
1047  {
1048  double x=X(i1,i2,i3,0);
1049  double y=X(i1,i2,i3,1);
1050 
1051  exmax(wt,wx,wy,nsources,xs[0],ys[0],tau[0],var[0],amp[0],period,x,y,time);
1052 
1053  #If #OPTION eq "solution"
1054  VEX(i1,i2,i3) = wy;
1055  VEY(i1,i2,i3) =-wx;
1056  VHZ(i1,i2,i3)= wt;
1057  #Elif #OPTION eq "error"
1058  ERREX(i1,i2,i3) = VEX(i1,i2,i3) - wy;
1059  ERREY(i1,i2,i3) = VEY(i1,i2,i3) + wx;
1060  ERRHZ(i1,i2,i3) = VHZ(i1,i2,i3) - wt;
1061 
1062  #Else
1063  U(i1,i2,i3,ex) = wy;
1064  U(i1,i2,i3,ey) =-wx;
1065  U(i1,i2,i3,hz) = wt;
1066  #End
1067 
1068  }
1069 }
1070 
1071 #endMacro
Generated on Fri Jan 4 2013 10:19:19 for CG by   doxygen 1.8.3