UltraScan III
us_stiffbase.cpp
Go to the documentation of this file.
1 
3 #include <QtCore>
4 #include <math.h>
5 #include "us_stiffbase.h"
6 
7 US_StiffBase::US_StiffBase()
8 {
9  n_gaussT = 28;
10  n_basisT = 3;
11 
12  xgT = new Gauss2D [ n_gaussT ];
13  phiT = new double* [ n_gaussT ];
14  phiT1 = new double* [ n_gaussT ];
15  phiT2 = new double* [ n_gaussT ];
16 
17  for ( int i = 0; i < n_gaussT; i++ )
18  {
19  phiT [ i ] = new double [ n_basisT ];
20  phiT1[ i ] = new double [ n_basisT ];
21  phiT2[ i ] = new double [ n_basisT ];
22  }
23 
24  n_gaussQ = 25;
25  n_basisQ = 4;
26 
27  xgQ = new Gauss2D [ n_gaussQ ];
28  phiQ = new double* [ n_gaussQ ];
29  phiQ1 = new double* [ n_gaussQ ];
30  phiQ2 = new double* [ n_gaussQ ];
31 
32  for ( int i = 0; i < n_gaussQ; i++ )
33  {
34  phiQ [ i ] = new double [ n_basisQ ];
35  phiQ1[ i ] = new double [ n_basisQ ];
36  phiQ2[ i ] = new double [ n_basisQ ];
37  }
38 
39  SetGauss();
40  LinearBasis();
41 
42 }
43 
44 US_StiffBase::~US_StiffBase()
45 {
46  delete [] xgT;
47 
48  for ( int i = 0; i < n_gaussT; i++ )
49  {
50  delete [] phiT [ i ];
51  delete [] phiT1[ i ];
52  delete [] phiT2[ i ];
53  }
54 
55  delete [] phiT;
56  delete [] phiT1;
57  delete [] phiT2;
58 
59 
60  delete [] xgQ;
61 
62  for ( int i = 0; i < n_gaussQ; i++ )
63  {
64  delete [] phiQ [ i ];
65  delete [] phiQ1[ i ];
66  delete [] phiQ2[ i ];
67  }
68 
69  delete [] phiQ;
70  delete [] phiQ1;
71  delete [] phiQ2;
72 }
73 
74 void US_StiffBase::CompLocalStif( int NK, double xd[4][2],
75  double D, double sw2, double** Stif )
76 {
77  double xg [ 2 ];
78  double jcbv[ 5 ];
79  double phi [ 4 ];
80  double phix[ 4 ];
81  double phiy[ 4 ];
82 
83  int n_gauss = ( NK == 3 ) ? n_gaussT : n_gaussQ;
84 
85  for ( int i = 0; i < 4; i++ )
86  for ( int j = 0; j < 4; j++ )
87  Stif[ i ][ j ] = 0.0;
88 
89  for ( int k = 0; k < n_gauss; k++ )
90  {
91  AffineMapping( NK, xd, k, xg );
92  Jacobian ( NK, xd, k, jcbv );
93 
94  double wt = ( NK == 3 ) ? xgT[ k ].w : xgQ[ k ].w;
95 
96  if ( NK == 3 )
97  {
98  for ( int i = 0; i < NK; i++ )
99  {
100  phi [ i ] = phiT [ k ][ i ];
101 
102  phix[ i ] = phiT1[ k ][ i ] * jcbv[ 1 ] +
103  phiT2[ k ][ i ] * jcbv[ 3 ];
104 
105  phiy[ i ] = phiT1[ k ][ i ] * jcbv[ 2 ] +
106  phiT2[ k ][ i ] * jcbv[ 4 ];
107  }
108  }
109  else
110  {
111  for ( int i = 0; i < NK; i++ )
112  {
113  phi [ i ] = phiQ [ k ][ i ];
114 
115  phix[ i ] = phiQ1[ k ][ i ] * jcbv[ 1 ] +
116  phiQ2[ k ][ i ] * jcbv[ 3 ];
117 
118  phiy[ i ] = phiQ1[ k ][ i ] * jcbv[ 2 ] +
119  phiQ2[ k ][ i ] * jcbv[ 4 ];
120  }
121  }
122 
123  for ( int i = 0; i < NK; i++ )
124  {
125  for ( int j = 0; j < NK; j++ )
126  {
127  double tmp = xg[ 0 ] * phiy[ j ] * phi[ i ] +
128  ( D * xg[ 0 ] * phix[ j ] -
129  sw2 * xg[ 0 ] * xg[ 0 ] * phi[ j ] ) * phix[ i ];
130 
131  Stif[ j ][ i ] += tmp * jcbv[ 0 ] * wt;
132  }
133  }
134  }
135 }
136 
137 void US_StiffBase::LambdaG( uint kk, double lam1, double lam2, double w,
138  Gauss2D* Lm )
139 {
140  double lam3 = 1.0 - lam1 - lam2;
141 
142  switch ( kk )
143  {
144  case 1: // one point
145  Lm[ 0 ].x = lam1;
146  Lm[ 0 ].y = lam1;
147  Lm[ 0 ].w = w;
148  break;
149 
150  case 3: // 3 points
151  Lm[ 0 ].x = lam1;
152  Lm[ 0 ].y = lam2;
153  Lm[ 0 ].w = w;
154 
155  Lm[ 1 ].x = lam3;
156  Lm[ 1 ].y = lam1;
157  Lm[ 1 ].w = w;
158 
159  Lm[ 2 ].x = lam2;
160  Lm[ 2 ].y = lam3;
161  Lm[ 2 ].w = w;
162  break;
163 
164  case 6: // 6 points
165  Lm[ 0 ].x = lam1;
166  Lm[ 0 ].y = lam2;
167  Lm[ 0 ].w = w;
168 
169  Lm[ 1 ].x = lam1;
170  Lm[ 1 ].y = lam3;
171  Lm[ 1 ].w = w;
172 
173  Lm[ 2 ].x = lam2;
174  Lm[ 2 ].y = lam1;
175  Lm[ 2 ].w = w;
176 
177  Lm[ 3 ].x = lam2;
178  Lm[ 3 ].y = lam3;
179  Lm[ 3 ].w = w;
180 
181  Lm[ 4 ].x = lam3;
182  Lm[ 4 ].y = lam1;
183  Lm[ 4 ].w = w;
184 
185  Lm[ 5 ].x = lam3;
186  Lm[ 5 ].y = lam2;
187  Lm[ 5 ].w = w;
188  break;
189  }
190 }
191 
192 void US_StiffBase::SetGauss( void )
193 {
194  double Gs1D[ 5 ];
195  double Gs1w[ 5 ]; // Gauss pts and weights in 1D
196 
197  // Set Gass Quadrature points on Quad elem [0,1]x[0,1]
198  switch( n_gaussQ )
199  {
200  case 25:
201 
202  Gs1D[ 0 ] = -sqrt ( 5.0 + 2.0 * sqrt( 10.0 / 7.0 ) ) / 3.0;
203  Gs1w[ 0 ] = 0.3 * ( 0.7 + 5.0 * sqrt( 0.7 ) ) /
204  ( 2.0 + 5.0 * sqrt( 0.7 ) );
205 
206  Gs1D[ 1 ] = -sqrt ( 5.0 - 2.0 * sqrt( 10.0 / 7.0 ) ) / 3.0;
207  Gs1w[ 1 ] = 0.3 * ( 0.7 - 5.0 * sqrt( 0.7 ) ) /
208  ( 2.0 - 5.0 * sqrt( 0.7 ) );
209 
210  Gs1D[ 2 ] = 0.0;
211  Gs1w[ 2 ] = 128.0 / 225.0;
212 
213  Gs1D[ 3 ] = -Gs1D[ 1 ];
214  Gs1w[ 3 ] = Gs1w[ 1 ];
215 
216  Gs1D[ 4 ] = -Gs1D[ 0 ];
217  Gs1w[ 4 ] = Gs1w[ 0 ];
218 
219  for ( int i = 0; i < 5; i++ )
220  {
221  for ( int j = 0; j < 5; j++ )
222  {
223  int ind = j + 5 * i;
224 
225  xgQ[ ind ].x = ( Gs1D[ i ] + 1.0 ) / 2.0;
226  xgQ[ ind ].y = ( Gs1D[ j ] + 1.0 ) / 2.0;
227  xgQ[ ind ].w = Gs1w[ i ] * Gs1w[ j ] / 4.0;
228  }
229  }
230  break;
231 
232  default:
233  qDebug() << "Gauss Quadrature on Quad elem with "
234  << QString::number( n_gaussQ )
235  << " points not implemented" ;
236  return;
237  }
238 
239  // Set Gauss Quadrature points on Triangular elem [0,1]x[0,1]
240  switch ( n_gaussT )
241  {
242  double lam1;
243  double lam2;
244  double w;
245 
246  case 28:
247  lam1 = 0.3333333333;
248  lam2 = 0.3333333333;
249  w = 0.2178563571;
250  LambdaG( 1, lam1, lam2, w / 4, xgT );
251 
252  lam1 = 0.1063354684;
253  lam2 = 0.1063354684;
254  w = 0.1104193374;
255  LambdaG( 3, lam1, lam2, w / 4, xgT + 1 );
256 
257  lam1 = 0.0000000000;
258  lam2 = 0.5000000002;
259  w = 0.0358939762;
260  LambdaG( 3, lam1, lam2, w / 4, xgT + 4 );
261 
262  lam1 = 0.0000000000;
263  lam2 = 0.0000000000;
264  w = 0.0004021278;
265  LambdaG( 3, lam1, lam2, w / 4, xgT + 7 );
266 
267  lam1 = 0.1171809171;
268  lam2 = 0.3162697959;
269  w = 0.1771348660;
270  LambdaG( 6, lam1, lam2, w / 4, xgT + 10 );
271 
272  lam1 = 0.0000000000;
273  lam2 = 0.2655651402;
274  w = 0.0272344079;
275  LambdaG( 6, lam1, lam2, w / 4, xgT + 16 );
276 
277  lam1 = 0.0000000000;
278  lam2 = 0.0848854223;
279  w = 0.0192969460;
280  LambdaG( 6, lam1, lam2, w / 4, xgT + 22 );
281  break;
282 
283  default:
284  qDebug() << "Gauss Quadrature on Triangular elem with "
285  << QString::number( n_gaussT )
286  << " not implemented";
287  }
288 }
289 
290 
291 void US_StiffBase::AffineMapping( int NK, double xd[4][2], int gauss_ind,
292  double x[5] )
293 {
294  x[ 0 ] = 0.0;
295  x[ 1 ] = 0.0;
296 
297  if ( NK == 3 )
298  {
299  // triangular element
300  for ( int i = 0; i < 3; i++ )
301  {
302  x[ 0 ] += xd[ i ][ 0 ] * phiT[ gauss_ind ][ i ];
303  x[ 1 ] += xd[ i ][ 1 ] * phiT[ gauss_ind ][ i ];
304  }
305  }
306  else
307  { // quadrilateral element
308  for ( int i = 0; i < 4; i++ )
309  {
310  x[ 0 ] += xd[ i ][ 0 ] * phiQ[ gauss_ind ][ i ];
311  x[ 1 ] += xd[ i ][ 1 ] * phiQ[ gauss_ind ][ i ];
312  }
313  }
314 }
315 
316 void US_StiffBase::Jacobian( int NK, double xd[4][2], int gauss_ind,
317  double jcbv[5] )
318 {
319  double J11 = 0.0;
320  double J12 = 0.0;
321  double J21 = 0.0;
322  double J22 = 0.0; // jcb = d_x/d_xi
323 
324  if ( NK == 3 )
325  {
326  // triangular element
327  for ( int i = 0; i < 3; i++ )
328  {
329  J11 += xd[ i ][ 0 ] * phiT1[ gauss_ind ][ i ];
330  J12 += xd[ i ][ 0 ] * phiT2[ gauss_ind ][ i ];
331  J21 += xd[ i ][ 1 ] * phiT1[ gauss_ind ][ i ];
332  J22 += xd[ i ][ 1 ] * phiT2[ gauss_ind ][ i ];
333  }
334  }
335  else
336  {
337  // quadrilateral element
338  for( int i = 0; i < 4; i++ )
339  {
340  J11 += xd[ i ][ 0 ] * phiQ1[ gauss_ind ][ i ];
341  J12 += xd[ i ][ 0 ] * phiQ2[ gauss_ind ][ i ];
342  J21 += xd[ i ][ 1 ] * phiQ1[ gauss_ind ][ i ];
343  J22 += xd[ i ][ 1 ] * phiQ2[ gauss_ind ][ i ];
344  }
345  }
346 
347  jcbv[ 0 ] = J11 * J22 - J12 * J21; // = det(jcb)
348  jcbv[ 1 ] = J22 / jcbv[ 0 ]; // = d_xi / d_x
349  jcbv[ 2 ] = -J12 / jcbv[ 0 ]; // = d_xi / d_y
350  jcbv[ 3 ] = -J21 / jcbv[ 0 ]; // = d_et / d_x
351  jcbv[ 4 ] = J11 / jcbv[ 0 ]; // = d_et / d_y
352 }
353 
354 void US_StiffBase::LinearBasis( void )
355 {
356  // Linear basis for triangular element: (0,0), (1,0), (0,1)
357 
358  for ( int i = 0; i < n_gaussT; i++ )
359  {
360  double x = xgT[ i ].x;
361  double y = xgT[ i ].y;
362 
363  phiT [ i ][ 0 ] = 1.0 - x - y;
364  phiT1[ i ][ 0 ] = -1.0;
365  phiT2[ i ][ 0 ] = -1.0;
366 
367  phiT [ i ][ 1 ] = x;
368  phiT1[ i ][ 1 ] = 1.0;
369  phiT2[ i ][ 1 ] = 0.0;
370 
371  phiT [ i ][ 2 ] = y;
372  phiT1[ i ][ 2 ] = 0.0;
373  phiT2[ i ][ 2 ] = 1.0;
374  }
375 
376  // Linear basis for Quadrilateral element: (0,0), (1,0), (1,1), (0,1)
377 
378  for ( int i = 0; i < n_gaussQ; i++ )
379  {
380  double x = xgQ[ i ].x;
381  double y = xgQ[ i ].y;
382 
383  phiQ [ i ][ 0 ] = ( 1.0 - x ) * ( 1.0 - y );
384  phiQ1[ i ][ 0 ] = -( 1.0 - y );
385  phiQ2[ i ][ 0 ] = -( 1.0 - x );
386 
387  phiQ [ i ][ 1 ] = ( 1.0 - y ) * x;
388  phiQ1[ i ][ 1 ] = 1.0 - y;
389  phiQ2[ i ][ 1 ] = -x;
390 
391  phiQ [ i ][ 2 ] = x * y;
392  phiQ1[ i ][ 2 ] = y;
393  phiQ2[ i ][ 2 ] = x;
394 
395  phiQ [ i ][ 3 ] = ( 1.0 - x ) * y;
396  phiQ1[ i ][ 3 ] = -y;
397  phiQ2[ i ][ 3 ] = 1.0 - x;
398  }
399 }
Generated on Sun Jun 21 2015 11:40:43 for UltraScan III by   doxygen 1.8.3.1-20130324