CG  Version 25
translationAndRotationSolution.h
Go to the documentation of this file.
1 // ==================================================================================
2 // Macro to evaluate the translation and rotation solution
3 //
4 // evalSolution : true=eval solution, false=eval error
5 // ==================================================================================
6 #beginMacro getTranslationAndRotationSolution(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  bool assignStress = s11c >=0 ;
26 
27  if( pdeVariation == SmParameters::hemp )
28  {
29  printF("\n\n **************** FIX ME: getTranslationAndRotationSolution: finish me for HEMP **********\n\n");
30  // OV_ABORT("error");
31  }
32 
33  // Here is the solution for large translation and rotation
34 
35  // hard code some numbers for now:
36  std::vector<real> & trd = parameters.dbase.get<std::vector<real> >("translationAndRotationSolutionData");
37 
38  real omega = trd[0]; // rotation rate
39  real xcenter = trd[1]; // center of rotation in the reference frame
40  real ycenter = trd[2]; // center of rotation in the reference frame
41  real zcenter = trd[3]; // center of rotation in the reference frame
42  real vcenter[3]={trd[4],trd[5],trd[6]}; // velocity of center
43  real rx[3]={cos(omega*t)-1.,-sin(omega*t), 0.};
44  real ry[3]={sin(omega*t) , cos(omega*t)-1., 0.};
45  real rxt[3]={-omega*sin(omega*t),-omega*cos(omega*t), 0.};
46  real ryt[3]={ omega*cos(omega*t),-omega*sin(omega*t), 0.};
47 
48 #define U0(x,y,z,n,t) (vcenter[n-uc]*(t) + rx[n-uc]*((x)-xcenter) + ry[n-uc]*((y)-ycenter))
49 #define U0T(x,y,z,n,t) (vcenter[n-uc] + rxt[n-uc]*((x)-xcenter) + ryt[n-uc]*((y)-ycenter))
50 #define U0X(x,y,z,n,t) ( rx[n-uc] )
51 #define U0Y(x,y,z,n,t) ( ry[n-uc] )
52 
53  if( t==0. )
54  printF("**** translationAndRotationSolution, t=%8.2e omega=%8.2e (x0,x1,x2)=(%8.2e,%8.2e,%8.2e) "
55  " (v0,v1,v2)=(%8.2e,%8.2e,%8.2e) *********\n",t,omega,xcenter,ycenter,zcenter,
56  vcenter[0],vcenter[1],vcenter[2]);
57 
58 
59  int i1,i2,i3;
60  if( mg.numberOfDimensions()==2 )
61  {
62  real z0=0.;
63  FOR_3D(i1,i2,i3,I1,I2,I3)
64  {
65  real x0 = X(i1,i2,i3,0);
66  real y0 = X(i1,i2,i3,1);
67  real u1 = U0(x0,y0,z0,uc,t);
68  real u2 = U0(x0,y0,z0,vc,t);
69 
70  if( evalSolution )
71  {
72  U(i1,i2,i3,uc) =u1;
73  U(i1,i2,i3,vc) =u2;
74  }
75  else
76  {
77  ERR(i1,i2,i3,uc) =U(i1,i2,i3,uc) - u1;
78  ERR(i1,i2,i3,vc) =U(i1,i2,i3,vc) - u2;
79  }
80 
81  }
82 
83  if( assignVelocities )
84  {
85  FOR_3D(i1,i2,i3,I1,I2,I3) // loop over all points
86  {
87  real x0 = X(i1,i2,i3,0);
88  real y0 = X(i1,i2,i3,1);
89  real v1 = U0T(x0,y0,z0,uc,t);
90  real v2 = U0T(x0,y0,z0,vc,t);
91 
92  // printF(" *** assignSpecial: v1=%e v2=%e\n",v1,v2);
93 
94  if( evalSolution )
95  {
96  U(i1,i2,i3,v1c) = v1;
97  U(i1,i2,i3,v2c) = v2;
98  }
99  else
100  {
101  ERR(i1,i2,i3,v1c) = U(i1,i2,i3,v1c) - v1;
102  ERR(i1,i2,i3,v2c) = U(i1,i2,i3,v2c) - v2;
103  }
104 
105  }
106  }
107  if( assignStress )
108  {
109  FOR_3D(i1,i2,i3,I1,I2,I3) // loop over all points
110  {
111  real x0 = X(i1,i2,i3,0);
112  real y0 = X(i1,i2,i3,1);
113  real f11 = 1. + U0X(x0,y0,z0,uc,t);
114  real f12 = U0Y(x0,y0,z0,uc,t);
115  real f21 = U0X(x0,y0,z0,vc,t);
116  real f22 = 1. + U0Y(x0,y0,z0,vc,t);
117  real e11 = .5*(f11*f11+f21*f21-1.); // this is E(i,j), symmetric
118  real e12 = .5*(f11*f12+f21*f22 );
119  real e22 = .5*(f12*f12+f22*f22-1.);
120  real trace = e11 + e22;
121  real s11 = lambda*trace + 2*mu*e11; // this is S(i,j), symmetric
122  real s12 = 2*mu*e12;
123  real s21 = s12;
124  real s22 = lambda*trace + 2*mu*e22;
125  real p11 = s11*f11 + s12*f12; // this P(i,j)
126  real p12 = s11*f21 + s12*f22;
127  real p21 = s21*f11 + s22*f12;
128  real p22 = s21*f21 + s22*f22;
129  if( evalSolution )
130  {
131  U(i1,i2,i3,s11c) = p11;
132  U(i1,i2,i3,s12c) = p12;
133  U(i1,i2,i3,s21c) = p21;
134  U(i1,i2,i3,s22c) = p22;
135  }
136  else
137  {
138  ERR(i1,i2,i3,s11c) = U(i1,i2,i3,s11c) - p11;
139  ERR(i1,i2,i3,s12c) = U(i1,i2,i3,s12c) - p12;
140  ERR(i1,i2,i3,s21c) = U(i1,i2,i3,s21c) - p21;
141  ERR(i1,i2,i3,s22c) = U(i1,i2,i3,s22c) - p22;
142  }
143  }
144  }
145  }
146  else
147  { // ***** 3D ****
148  OV_ABORT("translationAndRotationSolution: finish me for 3d");
149  FOR_3D(i1,i2,i3,I1,I2,I3)
150  {
151  real x0 = X(i1,i2,i3,0);
152  real y0 = X(i1,i2,i3,1);
153  real z0 = X(i1,i2,i3,2);
154  U(i1,i2,i3,uc) =U0(x0,y0,z0,uc,t);
155  U(i1,i2,i3,vc) =U0(x0,y0,z0,vc,t);
156  U(i1,i2,i3,wc) =U0(x0,y0,z0,wc,t);
157  }
158  }
159 
160 #undef U0
161 #undef U0T
162 #undef U0X
163 #undef U0Y
164 
165 
166 }
167 #endMacro
Generated on Fri Jan 4 2013 10:19:29 for CG by   doxygen 1.8.3