CG  Version 25
sphereEigenmode.h
Go to the documentation of this file.
1 // ==================================================================================
2 // Macro to evaluate eigen modes of the sphere
3 //
4 // evalSolution : true=eval solution, false=eval error
5 // ==================================================================================
6 #beginMacro getSphereEigenmode(evalSolution,U,ERR,X,t,I1,I2,I3)
7 {
8  const int v1c = parameters.dbase.get<int >("v1c");
9  const int v2c = parameters.dbase.get<int >("v2c");
10  const int v3c = parameters.dbase.get<int >("v3c");
11 
12  bool assignVelocities= v1c>=0 ;
13  const int s11c = parameters.dbase.get<int >("s11c");
14  const int s12c = parameters.dbase.get<int >("s12c");
15  const int s13c = parameters.dbase.get<int >("s13c");
16  const int s21c = parameters.dbase.get<int >("s21c");
17  const int s22c = parameters.dbase.get<int >("s22c");
18  const int s23c = parameters.dbase.get<int >("s23c");
19  const int s31c = parameters.dbase.get<int >("s31c");
20  const int s32c = parameters.dbase.get<int >("s32c");
21  const int s33c = parameters.dbase.get<int >("s33c");
22 
23  const int pc = parameters.dbase.get<int >("pc");
24 
25  const real & rho=parameters.dbase.get<real>("rho");
26  const real & mu = parameters.dbase.get<real>("mu");
27  const real & lambda = parameters.dbase.get<real>("lambda");
28  const RealArray & muGrid = parameters.dbase.get<RealArray>("muGrid");
29  const RealArray & lambdaGrid = parameters.dbase.get<RealArray>("lambdaGrid");
30 
31 
32  bool assignStress = s11c >=0 ;
33 
34  if( pdeVariation == SmParameters::hemp )
35  {
36  printF("\n\n **************** FIX ME: getSphereEigenmode: finish me for HEMP **********\n\n");
37  // OV_ABORT("error");
38  }
39 
40  assert( mg.numberOfDimensions()==3 );
41 
42  // parameters are stored here:
43  std::vector<real> & data = parameters.dbase.get<std::vector<real> >("sphereEigenmodeData");
44 
45  const int vibrationClass = (int)data[0];
46  const int n = (int)data[1];
47  const int m = (int)data[2];
48  const real rad = data[3];
49 
50  const real cp =sqrt( (lambda+2.*mu)/rho );
51  const real cs = sqrt( mu/rho );
52 
53  // The eigenvalues z = kappa*a satisfy (n=mode)
54  // (n-1)*psi_n(z) + z*psi_n'(z) = 0
55 
56  // cgDoc/sm/sphere.maple:
57  // n=1 : m=1 x/Pi = 1.83456604098843e+00
58  // n=1 : m=2 x/Pi = 2.89503202144422e+00
59 
60  // n=2 : m=1 x/Pi = 7.96135239733083e-01
61  // n=2 : m=2 x/Pi = 2.27146214644857e+00
62 
63  // kappa = omega /cs
64  // h = omega/cp
65 
66  real omega;
67  real cPhi=1., cOmega=1.;
68  if( vibrationClass==1 )
69  {
70  // eigenvalues are kappa*a
71  real radialEigenvalue;
72  if( n==1 )
73  {
74  if( m==1 )
75  radialEigenvalue = 1.83456604098843*Pi;
76  else if( m==2 )
77  radialEigenvalue = 2.89503202144422*Pi;
78  else
79  {
80  OV_ABORT("getSphereEigenmode: invalid value for m");
81  }
82  }
83  else if( n==2 )
84  {
85  assert( m==1 );
86 
87  if( m==1 )
88  radialEigenvalue = 7.96135239733083e-01*Pi;
89  else if( m==2 )
90  radialEigenvalue = 2.27146214644857*Pi;
91  else
92  {
93  OV_ABORT("getSphereEigenmode: invalid value for m");
94  }
95  }
96  omega = (radialEigenvalue/rad)*cs;
97 
98  }
99  else if( vibrationClass==2 )
100  {
101  // See cgDoc/sm/sphere2.maple
102  // % sphere1.maple:
103  // % Class II : n=0 : root m=1 x=4.43999821235653e+00, x/Pi = 1.41329532563144e+00, h*a/Pi=8.15966436697752e-01
104  // % Class II : n=0 : root m=2 x=1.04939244112140e+01, x/Pi = 3.34031988495483e+00, h*a/Pi=1.92853458475813e+00
105  // % Class II : n=0 : root m=3 x=1.60731909374045e+01, x/Pi = 5.11625557789555e+00, h*a/Pi=2.95387153514092e+00
106  // % Class II : n=0 : root m=4 x=2.15793450854558e+01, x/Pi = 6.86891887807218e+00, h*a/Pi=3.96577216329668e+00
107  real radialEigenvalue;
108  if( n==0 )
109  { // radial vibrations
110  // h*a
111  const real eig[4] ={ //
112  8.15966436697752e-01,
113  1.92853458475813e+00,
114  2.95387153514092e+00,
115  3.96577216329668e+00
116  };
117 
118  if( m>=1 && m<=4 )
119  {
120  radialEigenvalue = eig[m-1]*Pi; // Love: m=0 : .8160
121  }
122  else
123  {
124  printF("sphereEigenmode: ERROR: n=%i m=%i not available\n",n,m);
125  OV_ABORT("ERROR");
126  }
127 
128  omega = (radialEigenvalue/rad)*cp;
129  }
130  else if( n==2 )
131  {
132  // % n=2 : root m=1 x=2.63986927790186e+00, x/Pi = 8.40296489389027e-01, kappa*a/Pi=8.40296489389027e-01
133  // % an =-1.13478082789981e-02, cn=-3.02305786589451e-02, -an/cn=-3.75375159272393e-01
134  // % bn =-2.04041203329361e-02, dn=-5.43566078599509e-02, -bn/dn=-3.75375159272393e-01
135  // %
136  // % n=2 : root m=2 x=4.86527284993742e+00, x/Pi = 1.54866444711667e+00, kappa*a/Pi=1.54866444711667e+00
137  // % an =1.50294520720980e-02, cn=1.88125428532884e-01, -an/cn=-7.98905931500425e-02
138  // % bn =-1.22064204309416e-02, dn=-1.52789207710809e-01, -bn/dn=-7.98905931500425e-02
139  // %
140  // % n=2 : root m=3 x=8.32919545905501e+00, x/Pi = 2.65126525857435e+00, kappa*a/Pi=2.65126525857435e+00
141  // % an =4.57627801659678e-03, cn=-9.86226192484144e-02, -an/cn=4.64019111586348e-02
142  // % bn =-9.34189056911829e-04, dn=2.01325556121623e-02, -bn/dn=4.64019111586348e-02
143  // %
144  // % n=2 : root m=4 x=9.78016346034290e+00, x/Pi = 3.11312271792062e+00, kappa*a/Pi=3.11312271792062e+00
145  // % an =7.28789419129781e-04, cn=4.50066777991950e-02, -an/cn=-1.61929174684121e-02
146  // % bn =1.21547211256669e-03, dn=7.50619593373295e-02, -bn/dn=-1.61929174684121e-02
147 
148  // kappa*a
149  const real eig[4] ={ //
150  8.40296489389027e-01,
151  1.54866444711667e+00,
152  2.65126525857435e+00,
153  3.11312271792062e+00
154  };
155  const real cwp[4] ={ //
156  -3.75375159272393e-01,
157  -7.98905931500425e-02,
158  4.64019111586348e-02,
159  -1.61929174684121e-02
160  };
161  // radialEigenvalue= 2.63986927790039;
162  // radialEigenvalue = .840*Pi; // Love p286
163  if( m>=1 && m<=4 )
164  {
165  radialEigenvalue=eig[m-1]*Pi;
166  cPhi=cwp[m-1];
167  }
168  else
169  {
170  printF("sphereEigenmode: ERROR: n=%i m=%i not available\n",n,m);
171  OV_ABORT("ERROR");
172  }
173  omega = (radialEigenvalue/rad)*cs;
174 
175  }
176  else
177  {
178 
179  OV_ABORT("finish me for vibrationClass==2");
180  }
181 
182  }
183  else
184  {
185  OV_ABORT("ERROR: vibrationClass");
186  }
187 
188 
189  real kappa,h;
190  // kappa = omega/cs
191  kappa = omega/cs;
192  // h = omega/cp
193  h = omega/cp;
194 
195  const real h2= h*h, h3=h*h*h, h4=h*h*h*h, kappa2=kappa*kappa, kappa3=kappa*kappa*kappa, kappa4=kappa*kappa*kappa*kappa;
196 
197 
198 
199  if( t==0. )
200  printF("**** getSphereEigenmode: t=%8.2e class=%i n=%i m=%i radius=%9.3e omega=%9.3e cp=%8.2e cs=%8.2e"
201  " h*a/pi=%9.3e kappa*a/pi=%9.3e\n",
202  t,vibrationClass,n,m,rad,omega,cp,cs,h*rad/Pi,kappa*rad/Pi);
203 
204  // printF("**** getSphereEigenmode: t=%8.2e class=%i evalSolution=%i \n",t,vibrationClass,evalSolution);
205 
206  real cost = cos(omega*t);
207  real sint = sin(omega*t);
208 
209  // vibrationClass==1 solution:
210  //
211 
212  if( vibrationClass==1 )
213  {
214  real amp=10./rad; // amplitude
215  if( n==2 )
216  amp=100./rad;
217 
218  const real a0xy =1.; // shear in the x-y plane
219  const real a0yz =.8; // shear in the y-z plane
220  const real a0zx =.6; // shear in the z-x plane
221 #define U1(x,y,z) (amp*cost*psi*( a0xy*(y) -a0zx*(z) ))
222 #define V1(x,y,z) (amp*cost*psi*(-a0xy*(x)+a0yz*(z) ))
223 #define W1(x,y,z) (amp*cost*psi*( -a0yz*(y)+a0zx*(x) ))
224 
225  // psi has an asymptotic series for small argument:
226  real psi0, psi1= 1./(2.*(2.*n+3.));
227  if( n==1 )
228  {
229  psi0 = -1./(1.*3.);
230  }
231  else if( n==2 )
232  {
233  psi0 = 1./(1.*3.*5.);
234  }
235  else
236  {
237  OV_ABORT("finish me for n>2");
238  }
239 
240 
241  const real rEps = pow(REAL_EPSILON,.25); // use small r approximation when kappa*r < rEps
242  // const real rEps = pow(REAL_EPSILON,.5); // use small r approximation when kappa*r < rEps
243 
244  int i1,i2,i3;
245  real x0,y0,z0,r,kr,psi,u1,u2,u3;
246 
247  FOR_3D(i1,i2,i3,I1,I2,I3)
248  {
249  x0 = X(i1,i2,i3,0);
250  y0 = X(i1,i2,i3,1);
251  z0 = X(i1,i2,i3,2);
252 
253  r = sqrt( x0*x0 + y0*y0 + z0*z0 );
254  kr = kappa*r;
255 
256  if( n==1 )
257  {
258  // Here is psi_1
259  if( kr < rEps ) // use taylor series: (Love p279)
260  psi = psi0*( 1. - psi1*kr*kr ); // + O( kr^4 )
261  else
262  psi = (kr*cos(kr)-sin(kr))/(kr*kr*kr);
263 
264  u1 = U1(x0,y0,z0);
265  u2 = V1(x0,y0,z0);
266  u3 = W1(x0,y0,z0);
267  }
268  else if( n==2 )
269  {
270 #define U2(x,y,z) (amp*cost*psi*( a0xy*(y)*(z) -a0zx*(z)*(y) ))
271 #define V2(x,y,z) (amp*cost*psi*(-a0xy*(x)*(z)+a0yz*(z)*(x) ))
272 #define W2(x,y,z) (amp*cost*psi*( -a0yz*(y)*(x) +a0zx*(x)*(y) ))
273 
274  // Here is psi_2
275  if( kr < rEps ) // use taylor series: (Love p279)
276  psi = psi0*( 1. - psi1*kr*kr ); // + O( kr^4 )
277  else
278  psi = ( 3.*kr*cos(kr)+ (kr*kr-3.)*sin(kr) )/(kr*kr*kr*kr*kr);
279 
280  u1 = U2(x0,y0,z0);
281  u2 = V2(x0,y0,z0);
282  u3 = W2(x0,y0,z0);
283  }
284  else
285  {
286  OV_ABORT("un-implemented value for n");
287  }
288 
289 
290  if( evalSolution )
291  {
292  U(i1,i2,i3,uc) =u1;
293  U(i1,i2,i3,vc) =u2;
294  U(i1,i2,i3,wc) =u3;
295  }
296  else
297  {
298  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
299  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
300  ERR(i1,i2,i3,wc) =U(i1,i2,i3,wc) - u3;
301  }
302 
303  }
304 
305  if( assignVelocities )
306  {
307 
308  OV_ABORT("getSphereEigenmode: finish me");
309 
310  FOR_3D(i1,i2,i3,I1,I2,I3) // loop over all points
311  {
312  if( evalSolution )
313  {
314  }
315  else
316  {
317  }
318 
319  }
320  }
321  if( assignStress )
322  {
323 
324  OV_ABORT("getSphereEigenmode: finish me");
325 
326  FOR_3D(i1,i2,i3,I1,I2,I3) // loop over all points
327  {
328  if( evalSolution )
329  {
330  }
331  else
332  {
333  }
334  }
335  }
336  }
337  else if( vibrationClass==2 )
338  {
339  const real amp=50./rad; // amplitude -- this makes |u| about 1
340 
341  const real rEps = pow(REAL_EPSILON,.25); // use small r approximation when kappa*r < rEps
342  // const real rEps = pow(REAL_EPSILON,.5); // use small r approximation when kappa*r < rEps
343 
344  // psi has an asymptotic series for small argument:
345  // psi_n(x) = psin0*( 1. - psin1*x^2 + ... )
346  const real psi10 = -1./(1.*3.);
347  const real psi20 = 1./(1.*3.*5.);
348  const real psi30 = -1./(1.*3.*5.*7.);
349  const real psi40 = 1./(1.*3.*5.*7.*9.);
350 
351  const real psi11 = 1./(2.*(2.*1.+3.));
352  const real psi21 = 1./(2.*(2.*2.+3.));
353  const real psi31 = 1./(2.*(2.*3.+3.));
354  const real psi41 = 1./(2.*(2.*4.+3.));
355 
356  int m1=0; // m value for first solid harmonic "wn"
357  int m2=0; // m value for second solid harmonic "phi"
358 
359  const real n2p1 = 2.*n+1.;
360 
361  int i1,i2,i3;
362  real x0,y0,z0,r;
363  real kr,kr2,kr3,kr4,kr5,kr7,kr9;
364  real hr,hr2,hr3,hr4,hr5,hr7,hr9;
365  real u1,u2,u3;
366  real v1,v2,v3;
367  real u1x,u1y,u1z, u2x,u2y,u2z, u3x,u3y,u3z,div;
368  real s11,s12,s13,s21,s22,s23,s31,s32,s33;
369  real psi1hr, psi2hr, psi3hr, psi4hr;
370  real psi1kr, psi2kr, psi3kr, psi4kr;
371 
372  FOR_3D(i1,i2,i3,I1,I2,I3)
373  {
374  x0 = X(i1,i2,i3,0);
375  y0 = X(i1,i2,i3,1);
376  z0 = X(i1,i2,i3,2);
377 
378  r = sqrt( x0*x0 + y0*y0 + z0*z0 );
379 
380  kr =kappa*r;
381  kr2=kr*kr;
382  kr3=kr2*kr;
383  kr5=kr3*kr2;
384  kr7=kr5*kr2;
385  real coskr = cos(kr);
386  real sinkr = sin(kr);
387 
388  hr = h*r;
389  hr2=hr*hr;
390  hr3=hr2*hr;
391  hr5=hr3*hr2;
392  hr7=hr5*hr2;
393  real coshr = cos(hr);
394  real sinhr = sin(hr);
395 
396  if( n==0 )
397  {
398  // radial vibrations
399 
400  if( hr < rEps ) // use taylor series: (Love p279)
401  {
402  psi1hr = psi10*( 1. - psi11*hr2 ); // + O( hr^4 )
403  }
404  else
405  {
406  psi1hr = ( hr *coshr - sinhr )/hr3;
407  }
408 
409  u1 = amp*cost*( x0*psi1hr );
410  u2 = amp*cost*( y0*psi1hr );
411  u3 = amp*cost*( z0*psi1hr );
412  if( assignVelocities )
413  {
414  v1 = -amp*omega*sint*( x0*psi1hr );
415  v2 = -amp*omega*sint*( y0*psi1hr );
416  v3 = -amp*omega*sint*( z0*psi1hr );
417  }
418  if( assignStress )
419  {
420  if( hr < rEps ) // use taylor series: (Love p279)
421  psi2hr = psi20*( 1. - psi21*hr2 ); // + O( hr^4 )
422  else
423  psi2hr = - ( hr2*sinhr + 3.*hr*coshr - 3.*sinhr)/hr5;
424 
425  // u1 = x*psi1(hr)
426  // u1x = psi1(hr) + x*h*x/r*psi1'(hr) = psi1r + (h*x)^2 psi2(hr)
427  u1x = amp*cost*( psi1hr + h*h*x0*x0*psi2hr );
428  u1y = amp*cost*( h*h*x0*y0*psi2hr );
429  u1z = amp*cost*( h*h*x0*z0*psi2hr );
430 
431  u2x = amp*cost*( h*h*y0*x0*psi2hr );
432  u2y = amp*cost*( psi1hr + h*h*y0*y0*psi2hr );
433  u2z = amp*cost*( h*h*y0*z0*psi2hr );
434 
435  u3x = amp*cost*( h*h*z0*x0*psi2hr );
436  u3y = amp*cost*( h*h*z0*y0*psi2hr );
437  u3z = amp*cost*( psi1hr + h*h*z0*z0*psi2hr );
438  }
439 
440  }
441  else if( n==2 )
442  {
443  // cgDoc/sm/sphericalHarmonics.maple
444  // > lprint(psi1);
445  //1/x^3*(cos(x)*x-sin(x))
446  //> lprint(psi2);
447  //-(3*cos(x)*x-3*sin(x)+sin(x)*x^2)/x^5
448  //> lprint(psi3);
449  //-(x^3*cos(x)-6*sin(x)*x^2-15*cos(x)*x+15*sin(x))/x^7
450 
451  // We need
452  // psi2(hr), psi3(hr)
453  // psi1(kr), psi3(kr)
454 
455  if( hr < rEps ) // use taylor series: (Love p279)
456  {
457  psi2hr = psi20*( 1. - psi21*hr2 ); // + O( hr^4 )
458  psi3hr = psi30*( 1. - psi31*hr2 ); // + O( hr^4 )
459  }
460  else
461  {
462  psi2hr = - ( hr2*sinhr + 3.*hr*coshr - 3.*sinhr)/hr5;
463  psi3hr = - ( hr3*coshr - 6.*hr2*sinhr - 15.*hr*coshr +15.*sinhr)/hr7;
464  }
465 
466  if( kr < rEps ) // use taylor series: (Love p279)
467  {
468  psi1kr = psi10*( 1. - psi11*kr2 ); // + O( kr^4 )
469  psi3kr = psi30*( 1. - psi31*kr2 ); // + O( kr^4 )
470  }
471  else
472  {
473  psi1kr = ( kr *coskr - sinkr )/kr3;
474  psi3kr = - ( kr3*coskr - 6.*kr2*sinkr - 15.*kr*coskr +15.*sinkr )/kr7;
475  }
476 
477  // wn = r^2*Y_2^0 = ( 2*z^2 - x^2 - y^2 )
478  // wn = r^2*Y_2^1 = ( z*(a*x+b*y) )
479  // wn = r^2*Y_2^2 = ( a*(x^2-y^2) + b*x*y )
480  real w,wx,wy,wz;
481  real wxx,wxy,wxz, wyy,wyz,wzz;
482 
483  if( m1==0 )
484  {
485  w = 2.*z0*z0 - x0*x0 - y0*y0;
486  wx = -2.*x0;
487  wy = -2.*y0;
488  wz = 4.*z0;
489  wxx = -2.; wxy=0.; wxz=0.;
490  wyy = -2.; wyz=0.;
491  wzz = 4.;
492  }
493  else
494  {
495  OV_ABORT("ERROR: unexpected m1");
496  }
497 
498  // Question: do we need to have the correct leading coeff for the solid harmonics?
499 
500  // NOTE:
501  // an*wn+cn*pn = 0
502  // bn*wn+dn*pn = 0
503  //
504  // p/w = -an/cn =-bn/dn
505  // sphere2.maple:
506  // an =-1.13478082790558e-02, cn=-3.02305786593820e-02, -an/cn=-3.75375159268876e-01
507  // bn =-2.04041203329458e-02, dn=-5.43566078598750e-02, -bn/dn=-3.75375159273096e-01
508 
509 //% % ---Digits 25:
510 //% n=2 : root m=1 x=2.63986927790186e+00, x/Pi = 8.40296489389027e-01, kappa*a/Pi=8.40296489389027e-01
511 //% an =-1.13478082789981e-02, cn=-3.02305786589451e-02, -an/cn=-3.75375159272393e-01
512 //% bn =-2.04041203329361e-02, dn=-5.43566078599509e-02, -bn/dn=-3.75375159272393e-01
513 //%
514 //% n=2 : root m=2 x=4.86527284993742e+00, x/Pi = 1.54866444711667e+00, kappa*a/Pi=1.54866444711667e+00
515 //% an =1.50294520720980e-02, cn=1.88125428532884e-01, -an/cn=-7.98905931500425e-02
516 //% bn =-1.22064204309416e-02, dn=-1.52789207710809e-01, -bn/dn=-7.98905931500425e-02
517 //%
518 //% n=2 : root m=3 x=8.32919545905501e+00, x/Pi = 2.65126525857435e+00, kappa*a/Pi=2.65126525857435e+00
519 //% an =4.57627801659678e-03, cn=-9.86226192484144e-02, -an/cn=4.64019111586348e-02
520 //% bn =-9.34189056911829e-04, dn=2.01325556121623e-02, -bn/dn=4.64019111586348e-02
521 //%
522 //% n=2 : root m=4 x=9.78016346034290e+00, x/Pi = 3.11312271792062e+00, kappa*a/Pi=3.11312271792062e+00
523 //% an =7.28789419129781e-04, cn=4.50066777991950e-02, -an/cn=-1.61929174684121e-02
524 //% bn =1.21547211256669e-03, dn=7.50619593373295e-02, -bn/dn=-1.61929174684121e-02
525 
526 // const real cw=1.;
527 // const real cPhi=-3.75375159268876e-01; // *********
528 
529 
530  real p,px,py,pz;
531  real pxx,pxy,pxz, pyy,pyz,pzz;
532  if( m2==0 )
533  {
534  p = cPhi*(2.*z0*z0 - x0*x0 - y0*y0);
535  px = cPhi*( -2.*x0 );
536  py = cPhi*( -2.*y0 );
537  pz = cPhi*( 4.*z0 );
538 
539  pxx = cPhi*(-2. ); pxy=0.; pxz=0.;
540  pyy = cPhi*(-2. ); pyz=0.;
541  pzz = cPhi*( 4. );
542  }
543  else if( m2==1 )
544  {
545  // here is a linear combination of the real and im parts
546  // p = cPhi*( z0*( x0+ y0) );
547  // px = cPhi*( z0 );
548  // py = cPhi*( z0 );
549  // pz= cPhi*( x0+y0 );
550  }
551  else
552  {
553  OV_ABORT("ERROR: unexpected m2");
554  }
555 
556 #define VB2A(xa,wxa,pxa) \
557  ( (-1./(h*h))*( psi2hr + (hr2/n2p1)*psi3hr )*(wxa) \
558  + (1./n2p1)*psi3hr*( r*r*(wxa) - n2p1*(xa)*w ) \
559  + psi1kr*(pxa) - (n/(n+1.))*psi3kr*kappa*kappa*( r*r*(pxa) - n2p1*(xa)*p ) )
560 
561 // simpler:
562 #define VB2(xa,wxa,pxa) \
563  ( (-1./(h*h))*( psi2hr*(wxa) + h*h*(xa)*psi3hr*w ) \
564  + psi1kr*(pxa) - (n/(n+1.))*psi3kr*kappa*kappa*( r*r*(pxa) - n2p1*(xa)*p ) )
565 
566  real vb2x = VB2(x0,wx,px);
567  real vb2y = VB2(y0,wy,py);
568  real vb2z = VB2(z0,wz,pz);
569 
570  u1 = amp*cost*vb2x;
571  u2 = amp*cost*vb2y;
572  u3 = amp*cost*vb2z;
573 
574  if( assignVelocities )
575  {
576  v1 = -amp*omega*sint*vb2x;
577  v2 = -amp*omega*sint*vb2y;
578  v3 = -amp*omega*sint*vb2z;
579  }
580  if( assignStress )
581  {
582  // we need psi4(hr) and psi4(kr)
583  // > lprint(psi4);
584  // 1/x^9*(10*x^3*cos(x) -45*sin(x)*x^2+ x^4*sin(x) -105*cos(x)*x+105*sin(x))
585  kr4=kr2*kr2;
586  kr9=kr7*kr2;
587 
588  hr4=hr2*hr2;
589  hr9=hr7*hr2;
590 
591  if( hr < rEps ) // use taylor series: (Love p279)
592  psi4hr = psi40*( 1. - psi41*hr2 ); // + O( hr^4 )
593  else
594  psi4hr = ( (hr4-45.*hr2+105.)*sinhr + (10.*hr3-105.*hr)*coshr )/hr9;
595 
596  if( kr < rEps ) // use taylor series: (Love p279)
597  {
598  psi2kr = psi20*( 1. - psi21*kr2 ); // + O( kr^4 )
599  psi4kr = psi40*( 1. - psi41*kr2 ); // + O( kr^4 )
600  }
601  else
602  {
603  psi2kr = - ( kr2*sinkr + 3.*kr*coskr - 3.*sinkr )/kr5;
604  psi4kr = ( (kr4-45.*kr2+105.)*sinkr + (10.*kr3-105.*kr)*coskr )/kr9;
605  }
606 
607  // D( u_j )/D( x_a ) **check me **
608 #define VB2X(xj,wj,pj, xa,wa,pa, wja,pja, deltaja )\
609  ( (-1./(h2))*( h2*xa*psi3hr*wj + psi2hr*wja + h2*deltaja*psi3hr*w + h4*xj*xa*psi4hr*w + h2*xj*psi3hr*wa )\
610  + kappa2*xa*psi2kr*pj + psi1kr*pja \
611  - (n/(n+1.))*kappa4*xa*psi4kr*( r*r*pj - n2p1*xj*p )\
612  - (n/(n+1.))*kappa2*psi3kr*( 2.*xa*pj - n2p1*(deltaja)*p + r*r*pja - n2p1*xj*pa ) )
613 
614  u1x = amp*cost*( VB2X(x0,wx,px, x0,wx,px, wxx,pxx, 1. ) );
615  u1y = amp*cost*( VB2X(x0,wx,px, y0,wy,py, wxy,pxy, 0. ) );
616  u1z = amp*cost*( VB2X(x0,wx,px, z0,wz,pz, wxz,pxz, 0. ) );
617 
618  u2x = amp*cost*( VB2X(y0,wy,py, x0,wx,px, wxy,pxy, 0. ) );
619  u2y = amp*cost*( VB2X(y0,wy,py, y0,wy,py, wyy,pyy, 1. ) );
620  u2z = amp*cost*( VB2X(y0,wy,py, z0,wz,pz, wyz,pyz, 0. ) );
621 
622  u3x = amp*cost*( VB2X(z0,wz,pz, x0,wx,px, wxz,pxz, 0. ) );
623  u3y = amp*cost*( VB2X(z0,wz,pz, y0,wy,py, wyz,pyz, 0. ) );
624  u3z = amp*cost*( VB2X(z0,wz,pz, z0,wz,pz, wzz,pzz, 1. ) );
625 
626  }
627 
628 #undef VB2
629 
630  }
631  else
632  {
633  OV_ABORT("ERROR: unexpected n");
634  }
635 
636 
637  if( evalSolution )
638  {
639  U(i1,i2,i3,uc) =u1;
640  U(i1,i2,i3,vc) =u2;
641  U(i1,i2,i3,wc) =u3;
642  }
643  else
644  {
645  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
646  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
647  ERR(i1,i2,i3,wc) =U(i1,i2,i3,wc) - u3;
648  }
649 
650  if( assignVelocities )
651  {
652  if( evalSolution )
653  {
654  U(i1,i2,i3,v1c) =v1;
655  U(i1,i2,i3,v2c) =v2;
656  U(i1,i2,i3,v3c) =v3;
657  }
658  else
659  {
660  ERR(i1,i2,i3,v1c) =U(i1,i2,i3,v1c) - v1;
661  ERR(i1,i2,i3,v2c) =U(i1,i2,i3,v2c) - v2;
662  ERR(i1,i2,i3,v3c) =U(i1,i2,i3,v3c) - v3;
663  }
664 
665 
666  }
667  if( assignStress )
668  {
669  div = u1x+u2y+u3z;
670  s11 = lambda*div + 2.*mu*u1x;
671  s12 = mu*( u1y+u2x );
672  s13 = mu*( u1z+u3x );
673  s21 = s12;
674  s22 = lambda*div + 2.*mu*u2y;
675  s23 = mu*( u2z + u3y );
676  s31 = s13;
677  s32 = s23;
678  s33 = lambda*div + 2.*mu*u3z;
679 
680  if( evalSolution )
681  {
682 
683  U(i1,i2,i3,s11c) =s11;
684  U(i1,i2,i3,s12c) =s12;
685  U(i1,i2,i3,s13c) =s13;
686  U(i1,i2,i3,s21c) =s21;
687  U(i1,i2,i3,s22c) =s22;
688  U(i1,i2,i3,s23c) =s23;
689  U(i1,i2,i3,s31c) =s31;
690  U(i1,i2,i3,s32c) =s32;
691  U(i1,i2,i3,s33c) =s33;
692  }
693  else
694  {
695  ERR(i1,i2,i3,s11c) =U(i1,i2,i3,s11c) - s11;
696  ERR(i1,i2,i3,s12c) =U(i1,i2,i3,s12c) - s12;
697  ERR(i1,i2,i3,s13c) =U(i1,i2,i3,s13c) - s13;
698  ERR(i1,i2,i3,s21c) =U(i1,i2,i3,s21c) - s21;
699  ERR(i1,i2,i3,s22c) =U(i1,i2,i3,s22c) - s22;
700  ERR(i1,i2,i3,s23c) =U(i1,i2,i3,s23c) - s23;
701  ERR(i1,i2,i3,s31c) =U(i1,i2,i3,s31c) - s31;
702  ERR(i1,i2,i3,s32c) =U(i1,i2,i3,s32c) - s32;
703  ERR(i1,i2,i3,s33c) =U(i1,i2,i3,s33c) - s33;
704  }
705  }
706 
707 
708  } // end for3d
709 
710 
711  }
712  else
713  {
714  OV_ABORT("ERROR: unexpected vibrationClass");
715  }
716 
717 
718 #undef U2
719 #undef V2
720 #undef W2
721 
722 
723 }
724 
725 #endMacro
Generated on Fri Jan 4 2013 10:19:29 for CG by   doxygen 1.8.3