LifeV
FastAssemblerNS.cpp
Go to the documentation of this file.
1 #include <lifev/navier_stokes_blocks/solver/FastAssemblerNS.hpp>
2 #include <chrono>
3 #ifdef EPETRA_HAVE_OMP
4 #include <omp.h>
5 #endif
6 
7 using namespace std::chrono;
8 
9 using namespace std;
10 
11 using namespace LifeV;
12 
13 //=========================================================================
14 // Constructor //
15 FastAssemblerNS::FastAssemblerNS( const meshPtr_Type& mesh, const commPtr_Type& comm,
16  const ReferenceFE* refFE_velocity, const ReferenceFE* refFE_pressure,
17  const fespacePtr_Type& fespace_velocity, const fespacePtr_Type& fespace_pressure,
18  const qr_Type* qr ) :
19  M_mesh ( mesh ),
20  M_comm ( comm ),
21  M_referenceFE_velocity ( refFE_velocity ),
22  M_referenceFE_pressure ( refFE_pressure ),
23  M_fespace_velocity ( fespace_velocity ),
24  M_fespace_pressure ( fespace_pressure ),
25  M_qr ( qr ),
26  M_useSUPG ( false )
27 {
28 }
29 //=========================================================================
30 // Destructor //
31 FastAssemblerNS::~FastAssemblerNS()
32 {
33  //---------------------
34 
35  delete[] M_detJacobian;
36 
37  //---------------------
38 
39  for( int i = 0 ; i < M_numElements ; i++ )
40  {
41  for ( int j = 0; j < 3; j++ )
42  {
43  delete [] M_invJacobian[i][j];
44  }
45  delete [] M_invJacobian[i];
46  }
47  delete [] M_invJacobian;
48 
49  //---------------------
50 
51  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
52  {
53  delete [] M_phi_velocity[i];
54  }
55 
56  delete [] M_phi_velocity;
57 
58  //---------------------
59 
60  for( int i = 0 ; i < M_referenceFE_velocity->nbDof() ; i++ )
61  {
62  for ( int j = 0; j < M_qr->nbQuadPt(); j++ )
63  {
64  delete [] M_dphi_velocity[i][j];
65  }
66  delete [] M_dphi_velocity[i];
67  }
68  delete [] M_dphi_velocity;
69 
70  //---------------------
71 
72  for( int i = 0 ; i < M_referenceFE_velocity->nbDof() ; i++ )
73  {
74  for ( int j = 0; j < M_qr->nbQuadPt(); j++ )
75  {
76  for ( int k = 0; k < 3; k++ )
77  {
78  delete [] M_d2phi_velocity[i][j][k];
79  }
80  delete [] M_d2phi_velocity[i][j];
81  }
82  delete [] M_d2phi_velocity[i];
83  }
84  delete [] M_d2phi_velocity;
85 
86  //---------------------
87 
88  for ( int i = 0; i < M_referenceFE_pressure->nbDof(); i++ )
89  {
90  delete [] M_phi_pressure[i];
91  }
92 
93  delete [] M_phi_pressure;
94 
95  //---------------------
96 
97  for( int i = 0 ; i < M_referenceFE_pressure->nbDof() ; i++ )
98  {
99  for ( int j = 0; j < M_qr->nbQuadPt(); j++ )
100  {
101  delete [] M_dphi_pressure[i][j];
102  }
103  delete [] M_dphi_pressure[i];
104  }
105  delete [] M_dphi_pressure;
106 
107  //---------------------
108 
109  for( int i = 0 ; i < M_numElements ; i++ )
110  {
111  delete [] M_elements_velocity[i];
112  }
113  delete [] M_elements_velocity;
114 
115  //---------------------
116 
117  for( int i = 0 ; i < M_numElements ; i++ )
118  {
119  delete [] M_elements_pressure[i];
120  }
121  delete [] M_elements_pressure;
122 
123  //---------------------
124 
125  for( int i = 0 ; i < M_numElements ; i++ )
126  {
127  delete [] M_rows_velocity[i];
128  delete [] M_cols_velocity[i];
129  delete [] M_rows_pressure[i];
130  delete [] M_cols_pressure[i];
131  }
132  delete [] M_rows_velocity;
133  delete [] M_cols_velocity;
134  delete [] M_rows_pressure;
135  delete [] M_cols_pressure;
136 
137  //---------------------
138 
139  for( int i = 0 ; i < M_numElements ; i++ )
140  {
141  for ( int j = 0; j < M_referenceFE_velocity->nbDof(); j++ )
142  {
143  delete [] M_vals_00[i][j];
144  }
145  delete [] M_vals_00[i];
146  }
147  delete [] M_vals_00;
148 
149  //---------------------
150 
151  for( int i = 0 ; i < M_numElements ; i++ )
152  {
153  for ( int j = 0; j < M_referenceFE_velocity->nbDof(); j++ )
154  {
155  delete [] M_vals_01[i][j];
156  }
157  delete [] M_vals_01[i];
158  }
159  delete [] M_vals_01;
160 
161  //---------------------
162 
163  for( int i = 0 ; i < M_numElements ; i++ )
164  {
165  for ( int j = 0; j < M_referenceFE_pressure->nbDof(); j++ )
166  {
167  delete [] M_vals_10[i][j];
168  }
169  delete [] M_vals_10[i];
170  }
171  delete [] M_vals_10;
172 
173  //---------------------
174 
175  for( int i = 0 ; i < M_numElements ; i++ )
176  {
177  for ( int j = 0; j < M_referenceFE_pressure->nbDof(); j++ )
178  {
179  delete [] M_vals_11[i][j];
180  }
181  delete [] M_vals_11[i];
182  }
183  delete [] M_vals_11;
184 
185  if ( M_useSUPG )
186  {
187  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
188  {
189  for ( int k = 0; k < 3; k++ )
190  {
191  for ( int z = 0; z < 3; z++ )
192  {
193  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
194  {
195  delete [] M_vals_supg[i_elem][k][z][i];
196  }
197  delete [] M_vals_supg[i_elem][k][z];
198  }
199  delete [] M_vals_supg[i_elem][k];
200  }
201  delete [] M_vals_supg[i_elem];
202  }
203  delete [] M_vals_supg;
204 
205  for( int i = 0 ; i < M_numElements ; i++ )
206  {
207  for ( int k = 0; k < 3; k++ )
208  {
209  for ( int j = 0; j < M_referenceFE_velocity->nbDof(); j++ )
210  {
211  delete [] M_vals_supg_01[i][k][j];
212  }
213  delete [] M_vals_supg_01[i][k];
214  }
215  delete [] M_vals_supg_01[i];
216  }
217  delete [] M_vals_supg_01;
218 
219  for( int i = 0 ; i < M_numElements ; i++ )
220  {
221  for ( int k = 0; k < 3; k++ )
222  {
223  for ( int j = 0; j < M_referenceFE_pressure->nbDof(); j++ )
224  {
225  delete [] M_vals_supg_10[i][k][j];
226  }
227  delete [] M_vals_supg_10[i][k];
228  }
229  delete [] M_vals_supg_10[i];
230  }
231  delete [] M_vals_supg_10;
232 
233  for( int i = 0 ; i < M_numElements ; i++ )
234  {
235  delete [] M_rows_tmp[i];
236  delete [] M_cols_tmp[i];
237  }
238  delete [] M_rows_tmp;
239  delete [] M_cols_tmp;
240 
241 
242  for( int i = 0 ; i < M_numElements ; i++ )
243  {
244  for ( int j = 0; j < 3; j++ )
245  {
246  delete [] M_G[i][j];
247  }
248  delete [] M_G[i];
249  delete [] M_g[i];
250  delete [] M_Tau_M[i];
251  delete [] M_Tau_C[i];
252  delete [] M_Tau_M_hat[i];
253  }
254  delete [] M_G;
255  delete [] M_g;
256  delete [] M_Tau_M;
257  delete [] M_Tau_C;
258  delete [] M_Tau_M_hat;
259  }
260 
261 }
262 //=========================================================================
263 void
264 FastAssemblerNS::setConstants_NavierStokes( const Real& density, const Real& viscosity, const Real& timestep, const Real& orderBDF, const Real& C_I )
265 {
266  M_density = density;
267  M_viscosity = viscosity;
268  M_timestep = timestep;
269  M_orderBDF = orderBDF;
270  M_C_I = C_I;
271 }
272 //=========================================================================
273 void
274 FastAssemblerNS::setConstants_NavierStokes( const Real& density, const Real& viscosity, const Real& timestep, const Real& orderBDF, const Real& C_I, const Real& alpha )
275 {
276  M_density = density;
277  M_viscosity = viscosity;
278  M_timestep = timestep;
279  M_orderBDF = orderBDF;
280  M_C_I = C_I;
281  M_alpha = alpha;
282 }
283 //=========================================================================
284 void
285 FastAssemblerNS::updateGeoQuantities ( CurrentFE* current_fe )
286 {
287  for ( int i = 0; i < M_numElements; i++ )
288  {
289  current_fe->update( M_mesh->element (i), UPDATE_D2PHI );
290  M_detJacobian[i] = current_fe->detJacobian(0);
291  for ( int j = 0; j < 3; j++ )
292  {
293  for ( int k = 0; k < 3; k++ )
294  {
295  M_invJacobian[i][j][k] = current_fe->tInverseJacobian(j,k,2);
296  }
297  }
298  }
299 
300  if ( M_useSUPG )
301  {
302  for ( int i = 0; i < M_numElements; i++ )
303  {
304  M_G[i][0][0] = M_invJacobian[i][0][0]*M_invJacobian[i][0][0] + M_invJacobian[i][0][1]*M_invJacobian[i][0][1] + M_invJacobian[i][0][2]*M_invJacobian[i][0][2];
305  M_G[i][0][1] = M_invJacobian[i][0][0]*M_invJacobian[i][1][0] + M_invJacobian[i][0][1]*M_invJacobian[i][1][1] + M_invJacobian[i][0][2]*M_invJacobian[i][1][2];
306  M_G[i][0][2] = M_invJacobian[i][0][0]*M_invJacobian[i][2][0] + M_invJacobian[i][0][1]*M_invJacobian[i][2][1] + M_invJacobian[i][0][2]*M_invJacobian[i][2][2];
307 
308  M_G[i][1][0] = M_invJacobian[i][1][0]*M_invJacobian[i][0][0] + M_invJacobian[i][1][1]*M_invJacobian[i][0][1] + M_invJacobian[i][1][2]*M_invJacobian[i][0][2];
309  M_G[i][1][1] = M_invJacobian[i][1][0]*M_invJacobian[i][1][0] + M_invJacobian[i][1][1]*M_invJacobian[i][1][1] + M_invJacobian[i][1][2]*M_invJacobian[i][1][2];
310  M_G[i][1][2] = M_invJacobian[i][1][0]*M_invJacobian[i][2][0] + M_invJacobian[i][1][1]*M_invJacobian[i][2][1] + M_invJacobian[i][1][2]*M_invJacobian[i][2][2];
311 
312  M_G[i][2][0] = M_invJacobian[i][2][0]*M_invJacobian[i][0][0] + M_invJacobian[i][2][1]*M_invJacobian[i][0][1] + M_invJacobian[i][2][2]*M_invJacobian[i][0][2];
313  M_G[i][2][1] = M_invJacobian[i][2][0]*M_invJacobian[i][1][0] + M_invJacobian[i][2][1]*M_invJacobian[i][1][1] + M_invJacobian[i][2][2]*M_invJacobian[i][1][2];
314  M_G[i][2][2] = M_invJacobian[i][2][0]*M_invJacobian[i][2][0] + M_invJacobian[i][2][1]*M_invJacobian[i][2][1] + M_invJacobian[i][2][2]*M_invJacobian[i][2][2];
315  M_g[i][0] = M_invJacobian[i][0][0] + M_invJacobian[i][0][1] + M_invJacobian[i][0][2];
316  M_g[i][1] = M_invJacobian[i][1][0] + M_invJacobian[i][1][1] + M_invJacobian[i][1][2];
317  M_g[i][2] = M_invJacobian[i][2][0] + M_invJacobian[i][2][1] + M_invJacobian[i][2][2];
318  }
319  }
320 }
321 //=========================================================================
322 void
323 FastAssemblerNS::allocateSpace ( CurrentFE* current_fe_velocity, const bool& use_supg )
324 {
325  M_useSUPG = use_supg;
326 
327  M_numElements = M_mesh->numVolumes();
328 
329  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
330 
331  M_detJacobian = new double[ M_numElements ];
332 
333  M_invJacobian = new double** [ M_numElements ];
334 
335  //-------------------------------------------------------------------------------------------------
336 
337  M_invJacobian = new double** [ M_numElements];
338 
339  for ( int i = 0; i < M_numElements; i++ )
340  {
341  M_invJacobian[i] = new double* [ 3 ];
342  for ( int j = 0; j < 3 ; j++ )
343  {
344  M_invJacobian[i][j] = new double [ 3 ];
345  }
346  }
347 
348  for ( int i = 0; i < M_numElements; i++ )
349  {
350  current_fe_velocity->update( M_mesh->element (i), UPDATE_D2PHI );
351  M_detJacobian[i] = current_fe_velocity->detJacobian(0);
352  for ( int j = 0; j < 3; j++ )
353  {
354  for ( int k = 0; k < 3; k++ )
355  {
356  M_invJacobian[i][j][k] = current_fe_velocity->tInverseJacobian(j,k,2);
357  }
358  }
359  }
360 
361  //-------------------------------------------------------------------------------------------------
362  // PHI VELOCITY
363  //-------------------------------------------------------------------------------------------------
364 
365  M_phi_velocity = new double* [ M_referenceFE_velocity->nbDof() ];
366  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
367  {
368  M_phi_velocity[i] = new double [ M_qr->nbQuadPt() ];
369  }
370 
371  // PHI REF
372  for (UInt q (0); q < M_qr->nbQuadPt(); ++q)
373  {
374  for (UInt j (0); j < M_referenceFE_velocity->nbDof(); ++j)
375  {
376  M_phi_velocity[j][q] = M_referenceFE_velocity->phi (j, M_qr->quadPointCoor (q) );
377  }
378  }
379 
380  //-------------------------------------------------------------------------------------------------
381  // DPHI VELOCITY
382  //-------------------------------------------------------------------------------------------------
383 
384  M_dphi_velocity = new double** [ M_referenceFE_velocity->nbDof() ];
385 
386  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
387  {
388  M_dphi_velocity[i] = new double* [ M_qr->nbQuadPt() ];
389  for ( int j = 0; j < M_qr->nbQuadPt() ; j++ )
390  {
391  M_dphi_velocity[i][j] = new double [ 3 ];
392  }
393  }
394 
395  //DPHI REF
396  for (UInt q (0); q < M_qr->nbQuadPt(); ++q)
397  {
398  for (UInt i (0); i < M_referenceFE_velocity->nbDof(); ++i)
399  {
400  for (UInt j (0); j < 3; ++j)
401  {
402  M_dphi_velocity[i][q][j] = M_referenceFE_velocity->dPhi (i, j, M_qr->quadPointCoor (q) );
403  }
404  }
405  }
406 
407  //-------------------------------------------------------------------------------------------------
408  // D2PHI VELOCITY
409  //-------------------------------------------------------------------------------------------------
410 
411  M_d2phi_velocity = new double*** [ M_referenceFE_velocity->nbDof() ];
412 
413  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
414  {
415  M_d2phi_velocity[i] = new double** [ M_qr->nbQuadPt() ];
416  for ( int j = 0; j < M_qr->nbQuadPt() ; j++ )
417  {
418  M_d2phi_velocity[i][j] = new double* [ 3 ];
419  for ( int k = 0; k < 3 ; k++ )
420  {
421  M_d2phi_velocity[i][j][k] = new double [ 3 ];
422  }
423  }
424  }
425 
426  //D2PHI REF
427  for (UInt q (0); q < M_qr->nbQuadPt(); ++q)
428  {
429  for (UInt i (0); i < M_referenceFE_velocity->nbDof(); ++i)
430  {
431  for (UInt j (0); j < 3; ++j)
432  {
433  for (UInt k (0); k < 3; ++k)
434  {
435  M_d2phi_velocity[i][q][j][k] = M_referenceFE_velocity->d2Phi (i, j, k, M_qr->quadPointCoor (q) );
436  }
437  }
438  }
439  }
440 
441  //-------------------------------------------------------------------------------------------------
442  // PHI PRESSURE
443  //-------------------------------------------------------------------------------------------------
444 
445  M_phi_pressure = new double* [ M_referenceFE_pressure->nbDof() ];
446  for ( int i = 0; i < M_referenceFE_pressure->nbDof(); i++ )
447  {
448  M_phi_pressure[i] = new double [ M_qr->nbQuadPt() ];
449  }
450 
451  // PHI REF
452  for (UInt q (0); q < M_qr->nbQuadPt(); ++q)
453  {
454  for (UInt j (0); j < M_referenceFE_pressure->nbDof(); ++j)
455  {
456  M_phi_pressure[j][q] = M_referenceFE_pressure->phi (j, M_qr->quadPointCoor (q) );
457  }
458  }
459 
460  //-------------------------------------------------------------------------------------------------
461  // DPHI PRESSURE
462  //-------------------------------------------------------------------------------------------------
463 
464  M_dphi_pressure = new double** [ M_referenceFE_pressure->nbDof() ];
465 
466  for ( int i = 0; i < M_referenceFE_pressure->nbDof(); i++ )
467  {
468  M_dphi_pressure[i] = new double* [ M_qr->nbQuadPt() ];
469  for ( int j = 0; j < M_qr->nbQuadPt() ; j++ )
470  {
471  M_dphi_pressure[i][j] = new double [ 3 ];
472  }
473  }
474 
475  //DPHI REF
476  for (UInt q (0); q < M_qr->nbQuadPt(); ++q)
477  {
478  for (UInt i (0); i < M_referenceFE_pressure->nbDof(); ++i)
479  {
480  for (UInt j (0); j < 3; ++j)
481  {
482  M_dphi_pressure[i][q][j] = M_referenceFE_pressure->dPhi (i, j, M_qr->quadPointCoor (q) );
483  }
484  }
485  }
486 
487  //-------------------------------------------------------------------------------------------------
488  // CONNECTIVITY TEST
489  //-------------------------------------------------------------------------------------------------
490 
491  M_elements_velocity = new double* [ M_numElements ];
492 
493  for ( int i = 0; i < M_numElements; i++ )
494  {
495  M_elements_velocity[i] = new double [ M_referenceFE_velocity->nbDof() ];
496  }
497 
498  for ( int i = 0; i < M_numElements; i++ )
499  {
500  for ( int j = 0; j < M_referenceFE_velocity->nbDof(); j++ )
501  {
502  M_elements_velocity[i][j] = M_fespace_velocity->dof().localToGlobalMap (i, j);
503  }
504  }
505 
506  //-------------------------------------------------------------------------------------------------
507  // CONNECTIVITY TRIAL
508  //-------------------------------------------------------------------------------------------------
509 
510  M_elements_pressure = new double* [ M_numElements ];
511 
512  for ( int i = 0; i < M_numElements; i++ )
513  {
514  M_elements_pressure[i] = new double [ M_referenceFE_pressure->nbDof() ];
515  }
516 
517  for ( int i = 0; i < M_numElements; i++ )
518  {
519  for ( int j = 0; j < M_referenceFE_pressure->nbDof(); j++ )
520  {
521  M_elements_pressure[i][j] = M_fespace_pressure->dof().localToGlobalMap (i, j);
522  }
523  }
524 
525  //-------------------------------------------------------------------------------------------------
526  // ROWS VELOCITY
527  //-------------------------------------------------------------------------------------------------
528 
529  M_rows_velocity = new int* [M_numElements];
530 
531  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
532  {
533  M_rows_velocity[i_elem] = new int [ M_referenceFE_velocity->nbDof() ];
534  }
535 
536  //-------------------------------------------------------------------------------------------------
537  // ROWS PRESSURE
538  //-------------------------------------------------------------------------------------------------
539 
540  M_rows_pressure = new int* [M_numElements];
541 
542  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
543  {
544  M_rows_pressure[i_elem] = new int [ M_referenceFE_pressure->nbDof() ];
545  }
546 
547  //-------------------------------------------------------------------------------------------------
548  // COLS VELOCITY
549  //-------------------------------------------------------------------------------------------------
550 
551  M_cols_velocity = new int* [M_numElements];
552 
553  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
554  {
555  M_cols_velocity[i_elem] = new int [ M_referenceFE_velocity->nbDof() ];
556  }
557 
558  //-------------------------------------------------------------------------------------------------
559  // COLS PRESSURE
560  //-------------------------------------------------------------------------------------------------
561 
562  M_cols_pressure = new int* [M_numElements];
563 
564  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
565  {
566  M_cols_pressure[i_elem] = new int [ M_referenceFE_pressure->nbDof() ];
567  }
568 
569  //-------------------------------------------------------------------------------------------------
570  // VALUES BLOCK (0,0)
571  //-------------------------------------------------------------------------------------------------
572 
573  M_vals_00 = new double** [M_numElements];
574 
575  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
576  {
577  M_vals_00[i_elem] = new double* [ M_referenceFE_velocity->nbDof() ];
578  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
579  {
580  M_vals_00[i_elem][i] = new double [ M_referenceFE_velocity->nbDof() ];
581  }
582  }
583 
584  //-------------------------------------------------------------------------------------------------
585  // VALUES BLOCK (0,1)
586  //-------------------------------------------------------------------------------------------------
587 
588  M_vals_01 = new double** [M_numElements];
589 
590  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
591  {
592  M_vals_01[i_elem] = new double* [ M_referenceFE_velocity->nbDof() ];
593  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
594  {
595  M_vals_01[i_elem][i] = new double [ M_referenceFE_pressure->nbDof() ];
596  }
597  }
598 
599  //-------------------------------------------------------------------------------------------------
600  // VALUES BLOCK (1,0)
601  //-------------------------------------------------------------------------------------------------
602 
603  M_vals_10 = new double** [M_numElements];
604 
605  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
606  {
607  M_vals_10[i_elem] = new double* [ M_referenceFE_pressure->nbDof() ];
608  for ( int i = 0; i < M_referenceFE_pressure->nbDof(); i++ )
609  {
610  M_vals_10[i_elem][i] = new double [ M_referenceFE_velocity->nbDof() ];
611  }
612  }
613 
614  //-------------------------------------------------------------------------------------------------
615  // VALUES BLOCK (1,1)
616  //-------------------------------------------------------------------------------------------------
617 
618  M_vals_11 = new double** [M_numElements];
619 
620  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
621  {
622  M_vals_11[i_elem] = new double* [ M_referenceFE_pressure->nbDof() ];
623  for ( int i = 0; i < M_referenceFE_pressure->nbDof(); i++ )
624  {
625  M_vals_11[i_elem][i] = new double [ M_referenceFE_pressure->nbDof() ];
626  }
627  }
628 
629  if ( M_useSUPG )
630  {
631  //-------------------------------------------------------------------------------------------------
632  // SUPG MATRIX BLOCK (0,0) - because term with Tau_C leads to blocks on mixed components
633  //-------------------------------------------------------------------------------------------------
634 
635  M_vals_supg = new double**** [M_numElements];
636 
637  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
638  {
639  M_vals_supg[i_elem] = new double*** [ 3 ];
640  for ( int k = 0; k < 3; k++ )
641  {
642  M_vals_supg[i_elem][k] = new double** [ 3 ];
643  for ( int z = 0; z < 3; z++ )
644  {
645  M_vals_supg[i_elem][k][z] = new double* [ M_referenceFE_velocity->nbDof() ];
646  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
647  {
648  M_vals_supg[i_elem][k][z][i] = new double [ M_referenceFE_velocity->nbDof() ];
649  for ( int j = 0; j < M_referenceFE_velocity->nbDof(); j++ )
650  {
651  M_vals_supg[i_elem][k][z][i][j] = 0.0;
652  }
653  }
654  }
655  }
656  }
657 
658  M_vals_supg_01 = new double*** [M_numElements];
659  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
660  {
661  M_vals_supg_01[i_elem] = new double** [ 3 ];
662  for ( int k = 0; k < 3; k++ )
663  {
664  M_vals_supg_01[i_elem][k] = new double* [ M_referenceFE_velocity->nbDof() ];
665  for ( int i = 0; i < M_referenceFE_velocity->nbDof(); i++ )
666  {
667  M_vals_supg_01[i_elem][k][i] = new double [ M_referenceFE_pressure->nbDof() ];
668  }
669  }
670  }
671 
672  M_vals_supg_10 = new double*** [M_numElements];
673  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
674  {
675  M_vals_supg_10[i_elem] = new double** [ 3 ];
676  for ( int k = 0; k < 3; k++ )
677  {
678  M_vals_supg_10[i_elem][k] = new double* [ M_referenceFE_pressure->nbDof() ];
679  for ( int i = 0; i < M_referenceFE_pressure->nbDof(); i++ )
680  {
681  M_vals_supg_10[i_elem][k][i] = new double [ M_referenceFE_velocity->nbDof() ];
682  }
683  }
684  }
685 
686  M_rows_tmp = new int* [M_numElements];
687 
688  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
689  {
690  M_rows_tmp[i_elem] = new int [ M_referenceFE_velocity->nbDof() ];
691  }
692 
693  M_cols_tmp = new int* [M_numElements];
694 
695  for ( int i_elem = 0; i_elem < M_numElements; i_elem++ )
696  {
697  M_cols_tmp[i_elem] = new int [ M_referenceFE_velocity->nbDof() ];
698  }
699 
700  //-------------------------------------------------------------------------------------------------
701  // METRIC TENSOR AND VECTOR
702  //-------------------------------------------------------------------------------------------------
703 
704  M_G = new double** [ M_numElements];
705  M_g = new double* [ M_numElements];
706 
707  for ( int i = 0; i < M_numElements; i++ )
708  {
709  M_G[i] = new double* [ 3 ];
710  M_g[i] = new double [ 3 ];
711  for ( int j = 0; j < 3 ; j++ )
712  {
713  M_G[i][j] = new double [ 3 ];
714  }
715  }
716 
717  //-------------------------------------------------------------------------------------------------
718  // STABILIZATION PARAMETERS
719  //-------------------------------------------------------------------------------------------------
720 
721  M_Tau_M = new double* [ M_numElements];
722  M_Tau_C = new double* [ M_numElements];
723  M_Tau_M_hat = new double* [ M_numElements];
724  for ( int i = 0; i < M_numElements; i++ )
725  {
726  M_Tau_M[i] = new double [ M_qr->nbQuadPt() ];
727  M_Tau_C[i] = new double [ M_qr->nbQuadPt() ];
728  M_Tau_M_hat[i] = new double [ M_qr->nbQuadPt() ];
729  }
730 
731  //-------------------------------------------------------------------------------------------------
732 
733  for ( int i = 0; i < M_numElements; i++ )
734  {
735  M_G[i][0][0] = M_invJacobian[i][0][0]*M_invJacobian[i][0][0] + M_invJacobian[i][0][1]*M_invJacobian[i][0][1] + M_invJacobian[i][0][2]*M_invJacobian[i][0][2];
736  M_G[i][0][1] = M_invJacobian[i][0][0]*M_invJacobian[i][1][0] + M_invJacobian[i][0][1]*M_invJacobian[i][1][1] + M_invJacobian[i][0][2]*M_invJacobian[i][1][2];
737  M_G[i][0][2] = M_invJacobian[i][0][0]*M_invJacobian[i][2][0] + M_invJacobian[i][0][1]*M_invJacobian[i][2][1] + M_invJacobian[i][0][2]*M_invJacobian[i][2][2];
738 
739  M_G[i][1][0] = M_invJacobian[i][1][0]*M_invJacobian[i][0][0] + M_invJacobian[i][1][1]*M_invJacobian[i][0][1] + M_invJacobian[i][1][2]*M_invJacobian[i][0][2];
740  M_G[i][1][1] = M_invJacobian[i][1][0]*M_invJacobian[i][1][0] + M_invJacobian[i][1][1]*M_invJacobian[i][1][1] + M_invJacobian[i][1][2]*M_invJacobian[i][1][2];
741  M_G[i][1][2] = M_invJacobian[i][1][0]*M_invJacobian[i][2][0] + M_invJacobian[i][1][1]*M_invJacobian[i][2][1] + M_invJacobian[i][1][2]*M_invJacobian[i][2][2];
742 
743  M_G[i][2][0] = M_invJacobian[i][2][0]*M_invJacobian[i][0][0] + M_invJacobian[i][2][1]*M_invJacobian[i][0][1] + M_invJacobian[i][2][2]*M_invJacobian[i][0][2];
744  M_G[i][2][1] = M_invJacobian[i][2][0]*M_invJacobian[i][1][0] + M_invJacobian[i][2][1]*M_invJacobian[i][1][1] + M_invJacobian[i][2][2]*M_invJacobian[i][1][2];
745  M_G[i][2][2] = M_invJacobian[i][2][0]*M_invJacobian[i][2][0] + M_invJacobian[i][2][1]*M_invJacobian[i][2][1] + M_invJacobian[i][2][2]*M_invJacobian[i][2][2];
746  M_g[i][0] = M_invJacobian[i][0][0] + M_invJacobian[i][0][1] + M_invJacobian[i][0][2];
747  M_g[i][1] = M_invJacobian[i][1][0] + M_invJacobian[i][1][1] + M_invJacobian[i][1][2];
748  M_g[i][2] = M_invJacobian[i][2][0] + M_invJacobian[i][2][1] + M_invJacobian[i][2][2];
749  }
750  }
751 }
752 //=========================================================================
753 void
754 FastAssemblerNS::assemble_supg_terms( matrixPtr_Type& block00, matrixPtr_Type& block01, matrixPtr_Type& block10, matrixPtr_Type& block11, const vector_Type& u_h )
755 {
756  int ndof_velocity = M_referenceFE_velocity->nbDof();
757  int ndof_pressure = M_referenceFE_pressure->nbDof();
758  int NumQuadPoints = M_qr->nbQuadPt();
759  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
760  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
761 
762  double w_quad[NumQuadPoints];
763  for ( int q = 0; q < NumQuadPoints ; q++ )
764  {
765  w_quad[q] = M_qr->weight(q);
766  }
767 
768  #pragma omp parallel firstprivate( w_quad, ndof_velocity, ndof_pressure, NumQuadPoints)
769  {
770  int i_elem, i_dof, q, d1, d2, i_test, i_trial, e_idof, dim_mat, iCoor, jCoor, k1, k2;
771  double integral, integral_test, integral_trial, integral_partial, integral_lapl, partialSum;
772 
773  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3];
774  double dphi_phys_pressure[ndof_pressure][NumQuadPoints][3];
775  double d2phi_phys_velocity[ndof_velocity][NumQuadPoints][3][3];
776  double uhq[3][NumQuadPoints];
777 
778  // ELEMENTI,
779  #pragma omp for
780  for ( i_elem = 0; i_elem < M_numElements; i_elem++ )
781  {
782  // DOF
783  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
784  {
785  // QUAD
786  for ( q = 0; q < NumQuadPoints ; q++ )
787  {
788  // DIM 1
789  for ( d1 = 0; d1 < 3 ; d1++ )
790  {
791  dphi_phys_velocity[i_dof][q][d1] = 0.0;
792 
793  // DIM 2
794  for ( d2 = 0; d2 < 3 ; d2++ )
795  {
796  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
797  }
798  }
799  }
800  }
801 
802  // DOF -- nota che può essere ottimizzato rishapando il gradiente fisico. Prima q poi d1 poi d2 e poi i_dof
803  for ( i_dof = 0; i_dof < ndof_pressure; i_dof++ )
804  {
805  // QUAD
806  for ( q = 0; q < NumQuadPoints ; q++ )
807  {
808  // DIM 1
809  for ( d1 = 0; d1 < 3 ; d1++ )
810  {
811  dphi_phys_pressure[i_dof][q][d1] = 0.0;
812 
813  // DIM 2
814  for ( d2 = 0; d2 < 3 ; d2++ )
815  {
816  dphi_phys_pressure[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_pressure[i_dof][q][d2];
817  }
818  }
819  }
820  }
821 
822  // QUAD
823  for ( q = 0; q < NumQuadPoints ; q++ )
824  {
825  // DOF
826  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
827  {
828  // DIM 1
829  for ( iCoor = 0; iCoor < 3 ; iCoor++ )
830  {
831  // DIM 2
832  for ( jCoor = 0; jCoor < 3 ; jCoor++ )
833  {
834  partialSum = 0.0;
835  // DIM 1
836  for ( k1 = 0; k1 < 3 ; k1++ )
837  {
838  // DIM 2
839  for ( k2 = 0; k2 < 3 ; k2++ )
840  {
841  partialSum += M_invJacobian[i_elem][iCoor][k1]
842  * M_d2phi_velocity[i_dof][q][k1][k2]
843  * M_invJacobian[i_elem][jCoor][k2];
844  }
845  }
846  d2phi_phys_velocity[i_dof][q][iCoor][jCoor] = partialSum;
847  }
848  }
849  }
850  }
851 
852  // QUAD
853  for ( q = 0; q < NumQuadPoints ; q++ )
854  {
855  for ( d1 = 0; d1 < 3 ; d1++ )
856  {
857  uhq[d1][q] = 0.0;
858  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
859  {
860  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs ;
861  uhq[d1][q] += u_h[e_idof] * M_phi_velocity[i_dof][q];
862  }
863  }
864  }
865 
866  // STABILIZZAZIONE - coefficienti Tau_M e Tau_C
867  for ( q = 0; q < NumQuadPoints ; q++ )
868  {
869  M_Tau_M[i_elem][q] = 1.0/std::sqrt(
870  M_density*M_density*M_orderBDF*M_orderBDF/(M_timestep*M_timestep) // TAU_M_DEN_DT
871  +M_density*M_density*( uhq[0][q]*( M_G[i_elem][0][0]* uhq[0][q] + M_G[i_elem][0][1] * uhq[1][q] + M_G[i_elem][0][2] * uhq[2][q] ) +
872  uhq[1][q]*( M_G[i_elem][1][0]* uhq[0][q] + M_G[i_elem][1][1] * uhq[1][q] + M_G[i_elem][1][2] * uhq[2][q] ) +
873  uhq[2][q]*( M_G[i_elem][2][0]* uhq[0][q] + M_G[i_elem][2][1] * uhq[1][q] + M_G[i_elem][2][2] * uhq[2][q] )
874  ) // TAU_M_DEN_VEL
875  +M_C_I*M_viscosity*M_viscosity*(
876  M_G[i_elem][0][0]*M_G[i_elem][0][0] + M_G[i_elem][0][1]*M_G[i_elem][0][1] + M_G[i_elem][0][2]*M_G[i_elem][0][2] +
877  M_G[i_elem][1][0]*M_G[i_elem][1][0] + M_G[i_elem][1][1]*M_G[i_elem][1][1] + M_G[i_elem][1][2]*M_G[i_elem][1][2] +
878  M_G[i_elem][2][0]*M_G[i_elem][2][0] + M_G[i_elem][2][1]*M_G[i_elem][2][1] + M_G[i_elem][2][2]*M_G[i_elem][2][2]
879  ) // TAU_M_DEN_VISC*/
880  );
881 
882  M_Tau_C[i_elem][q] = 1.0/( M_g[i_elem][0]*M_Tau_M[i_elem][q]*M_g[i_elem][0] +
883  M_g[i_elem][1]*M_Tau_M[i_elem][q]*M_g[i_elem][1] +
884  M_g[i_elem][2]*M_Tau_M[i_elem][q]*M_g[i_elem][2]
885  );
886 
887  }
888 
889  // DOF - test
890  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
891  {
892  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
893 
894  // DOF - trial
895  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
896  {
897  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
898 
899  integral = 0.0;
900  // QUAD
901  for ( q = 0; q < NumQuadPoints ; q++ )
902  {
903  integral_test = 0.0;
904  integral_trial = 0.0;
905  integral_lapl = 0.0;
906 
907  // DIM 1
908  for ( d1 = 0; d1 < 3 ; d1++ )
909  {
910  integral_lapl += d2phi_phys_velocity[i_trial][q][d1][d1];
911  integral_test += uhq[d1][q] * dphi_phys_velocity[i_test][q][d1]; // w grad(phi_i)
912  integral_trial += uhq[d1][q] * dphi_phys_velocity[i_trial][q][d1]; // w grad(phi_j) terzo indice dphi_phys e' derivata in x,y,z
913  }
914  integral += M_Tau_M[i_elem][q] * (integral_test * integral_trial + integral_test * M_phi_velocity[i_trial][q]
915  - integral_test * integral_lapl ) * w_quad[q];
916  // above, the term "integral_test * M_phi[i_trial][q]" is the one which comes from (w grad(phi_i), phi_j)
917  }
918  M_vals_supg[i_elem][0][0][i_test][i_trial] = integral * M_detJacobian[i_elem]; // (w grad(phi_i), w grad(phi_j) )
919  M_vals_supg[i_elem][1][1][i_test][i_trial] = integral * M_detJacobian[i_elem];
920  M_vals_supg[i_elem][2][2][i_test][i_trial] = integral * M_detJacobian[i_elem];
921 
922  for ( d1 = 0; d1 < 3; d1++ )
923  {
924  for ( d2 = 0; d2 < 3; d2++ )
925  {
926  integral = 0.0;
927  // QUAD
928  for ( q = 0; q < NumQuadPoints ; q++ )
929  {
930  integral += M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][d1] * dphi_phys_velocity[i_trial][q][d2] * w_quad[q]; // div(phi_i) * div(phi_j)
931  }
932  M_vals_supg[i_elem][d1][d2][i_test][i_trial] += integral * M_detJacobian[i_elem];
933  }
934  }
935 
936  }
937 
938  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
939  {
940  for ( dim_mat = 0; dim_mat < 3 ; dim_mat++ )
941  {
942  M_cols_pressure[i_elem][i_trial] = M_elements_pressure[i_elem][i_trial];
943 
944  integral = 0.0;
945 
946  // QUAD
947  for ( q = 0; q < NumQuadPoints ; q++ )
948  {
949  integral_partial = 0.0;
950  for ( d1 = 0; d1 < 3 ; d1++ )
951  {
952  integral_partial += ( dphi_phys_velocity[i_test][q][d1] * uhq[d1][q] );
953  }
954  integral += M_Tau_M[i_elem][q] * dphi_phys_pressure[i_trial][q][dim_mat] * integral_partial * w_quad[q];
955  }
956  M_vals_supg_01[i_elem][dim_mat][i_test][i_trial] = integral * M_detJacobian[i_elem];
957  }
958  }
959 
960  }
961 
962  // DOF - test - assemblo blocchi (1,0) e (1,1)
963  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
964  {
965  M_rows_pressure[i_elem][i_test] = M_elements_pressure[i_elem][i_test];
966 
967  // DOF - trial
968  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
969  {
970  integral = 0.0;
971  // QUAD
972  for ( q = 0; q < NumQuadPoints ; q++ )
973  {
974  // DIM 1
975  for ( d1 = 0; d1 < 3 ; d1++ )
976  {
977  integral += M_Tau_M[i_elem][q] * dphi_phys_pressure[i_test][q][d1] * dphi_phys_pressure[i_trial][q][d1]*w_quad[q];
978  }
979  }
980  M_vals_11[i_elem][i_test][i_trial] = 0.0;
981  M_vals_11[i_elem][i_test][i_trial] = integral * M_detJacobian[i_elem];
982  }
983  }
984 
985  for ( dim_mat = 0; dim_mat < 3 ; dim_mat++ )
986  {
987  // DOF - test
988  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
989  {
990  // DOF - trial
991  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
992  {
993  integral = 0.0;
994 
995  // QUAD
996  for ( q = 0; q < NumQuadPoints ; q++ )
997  {
998  integral_partial = 0.0;
999  for ( d1 = 0; d1 < 3 ; d1++ )
1000  {
1001  integral_partial += ( dphi_phys_velocity[i_trial][q][d1] * uhq[d1][q] - d2phi_phys_velocity[i_trial][q][d1][d1] );
1002  }
1003  integral += M_Tau_M[i_elem][q] * ( dphi_phys_pressure[i_test][q][dim_mat] * M_phi_velocity[i_trial][q] +
1004  dphi_phys_pressure[i_test][q][dim_mat] * integral_partial ) * w_quad[q];
1005  }
1006 
1007  M_vals_supg_10[i_elem][dim_mat][i_test][i_trial] = integral * M_detJacobian[i_elem];
1008  }
1009  }
1010  }
1011  }
1012  }
1013 
1014  for ( UInt d1 = 0; d1 < 3 ; d1++ ) // row index
1015  {
1016  for ( int k = 0; k < M_numElements; ++k )
1017  {
1018  for ( UInt i = 0; i < ndof_velocity; i++ )
1019  {
1020  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
1021  }
1022  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k], ndof_velocity, M_cols_velocity[k], M_vals_supg[k][d1][0], Epetra_FECrsMatrix::ROW_MAJOR);
1023  block01->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k], ndof_pressure, M_cols_pressure[k], M_vals_supg_01[k][d1], Epetra_FECrsMatrix::ROW_MAJOR);
1024  }
1025 
1026  for ( int k = 0; k < M_numElements; ++k )
1027  {
1028  for ( UInt i = 0; i < ndof_velocity; i++ )
1029  {
1030  M_cols_tmp[k][i] = M_cols_velocity[k][i] + M_numScalarDofs;
1031  }
1032  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k], ndof_velocity, M_cols_tmp[k], M_vals_supg[k][d1][1], Epetra_FECrsMatrix::ROW_MAJOR);
1033  }
1034 
1035  for ( int k = 0; k < M_numElements; ++k )
1036  {
1037  for ( UInt i = 0; i < ndof_velocity; i++ )
1038  {
1039  M_cols_tmp[k][i] = M_cols_velocity[k][i] + 2 * M_numScalarDofs;
1040  }
1041  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k], ndof_velocity, M_cols_tmp[k], M_vals_supg[k][d1][2], Epetra_FECrsMatrix::ROW_MAJOR);
1042  }
1043  }
1044 
1045  for ( int k = 0; k < M_numElements; ++k )
1046  {
1047  block11->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k], ndof_pressure, M_cols_pressure[k], M_vals_11[k], Epetra_FECrsMatrix::ROW_MAJOR);
1048  }
1049 
1050  for ( UInt d1 = 0; d1 < 3 ; d1++ )
1051  {
1052  for ( int k = 0; k < M_numElements; ++k )
1053  {
1054  for ( UInt i = 0; i < ndof_velocity; i++ )
1055  {
1056  M_cols_tmp[k][i] = M_cols_velocity[k][i] + d1 * M_numScalarDofs;
1057  }
1058  block10->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k], ndof_velocity, M_cols_tmp[k], M_vals_supg_10[k][d1], Epetra_FECrsMatrix::ROW_MAJOR);
1059  }
1060  }
1061 }
1062 //=========================================================================
1063 void
1064 FastAssemblerNS::assemble_constant_terms( matrixPtr_Type& mass, matrixPtr_Type& stiffness, matrixPtr_Type& b01, matrixPtr_Type& b10 )
1065 {
1066  int ndof_velocity = M_referenceFE_velocity->nbDof();
1067  int ndof_pressure = M_referenceFE_pressure->nbDof();
1068  int NumQuadPoints = M_qr->nbQuadPt();
1069  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
1070  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
1071 
1072  double w_quad[NumQuadPoints];
1073  for ( int q = 0; q < NumQuadPoints ; q++ )
1074  {
1075  w_quad[q] = M_qr->weight(q);
1076  }
1077 
1078  #pragma omp parallel firstprivate( w_quad, ndof_velocity, ndof_pressure, NumQuadPoints)
1079  {
1080  int i_elem, i_dof, q, d1, d2, i_test, i_trial, dim_mat, dim_mat1, dim_mat2;
1081  double integral, integral_00, integral_01, integral_02, integral_10, integral_11, integral_12;
1082  double integral_20, integral_21, integral_22;
1083 
1084  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3]; // Gradient in the physical domain
1085 
1086  // ELEMENTI
1087  #pragma omp for
1088  for ( i_elem = 0; i_elem < M_numElements ; i_elem++ )
1089  {
1090  // DOF
1091  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1092  {
1093  // QUAD
1094  for ( q = 0; q < NumQuadPoints ; q++ )
1095  {
1096  // DIM 1
1097  for ( d1 = 0; d1 < 3 ; d1++ )
1098  {
1099  dphi_phys_velocity[i_dof][q][d1] = 0.0;
1100 
1101  // DIM 2
1102  for ( d2 = 0; d2 < 3 ; d2++ )
1103  {
1104  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
1105  }
1106  }
1107  }
1108  }
1109 
1110  // DOF - test
1111  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
1112  {
1113  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
1114 
1115  // DOF - trial
1116  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
1117  {
1118  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
1119  M_vals_00[i_elem][i_test][i_trial] = 0.0;
1120  integral = 0.0;
1121  // QUAD
1122  for ( q = 0; q < NumQuadPoints ; q++ )
1123  {
1124  integral += M_phi_velocity[i_test][q] * M_phi_velocity[i_trial][q]*w_quad[q];
1125  }
1126  M_vals_00[i_elem][i_test][i_trial] = M_density * integral * M_detJacobian[i_elem];
1127  }
1128  }
1129 
1130  // DOF - test
1131  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
1132  {
1133  // DOF - trial
1134  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
1135  {
1136  M_vals_supg[i_elem][0][0][i_test][i_trial] = 0.0;
1137  M_vals_supg[i_elem][0][1][i_test][i_trial] = 0.0;
1138  M_vals_supg[i_elem][0][2][i_test][i_trial] = 0.0;
1139  M_vals_supg[i_elem][1][0][i_test][i_trial] = 0.0;
1140  M_vals_supg[i_elem][1][1][i_test][i_trial] = 0.0;
1141  M_vals_supg[i_elem][1][2][i_test][i_trial] = 0.0;
1142  M_vals_supg[i_elem][2][0][i_test][i_trial] = 0.0;
1143  M_vals_supg[i_elem][2][1][i_test][i_trial] = 0.0;
1144  M_vals_supg[i_elem][2][2][i_test][i_trial] = 0.0;
1145  integral_00 = 0.0;
1146  integral_01 = 0.0;
1147  integral_02 = 0.0;
1148  integral_10 = 0.0;
1149  integral_11 = 0.0;
1150  integral_12 = 0.0;
1151  integral_20 = 0.0;
1152  integral_21 = 0.0;
1153  integral_22 = 0.0;
1154  // QUAD
1155  for ( q = 0; q < NumQuadPoints ; q++ )
1156  {
1157  integral_00 += ( 2*dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][0] +
1158  dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][1] +
1159  dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][2] ) * w_quad[q];
1160 
1161  integral_01 += dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][0] * w_quad[q];
1162 
1163  integral_02 += dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][0] * w_quad[q];
1164 
1165  integral_10 += dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][1] * w_quad[q];
1166 
1167  integral_11 += ( 2*dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][1] +
1168  dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][0] +
1169  dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][2] ) * w_quad[q];
1170 
1171  integral_12 += dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][1] * w_quad[q];
1172 
1173  integral_20 += dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][2] * w_quad[q];
1174 
1175  integral_21 += dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][2] * w_quad[q];
1176 
1177  integral_22 += ( 2*dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][2] +
1178  dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][0] +
1179  dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][1] ) * w_quad[q];
1180 
1181  }
1182  M_vals_supg[i_elem][0][0][i_test][i_trial] = M_viscosity * integral_00 * M_detJacobian[i_elem];
1183  M_vals_supg[i_elem][0][1][i_test][i_trial] = M_viscosity * integral_01 * M_detJacobian[i_elem];
1184  M_vals_supg[i_elem][0][2][i_test][i_trial] = M_viscosity * integral_02 * M_detJacobian[i_elem];
1185  M_vals_supg[i_elem][1][0][i_test][i_trial] = M_viscosity * integral_10 * M_detJacobian[i_elem];
1186  M_vals_supg[i_elem][1][1][i_test][i_trial] = M_viscosity * integral_11 * M_detJacobian[i_elem];
1187  M_vals_supg[i_elem][1][2][i_test][i_trial] = M_viscosity * integral_12 * M_detJacobian[i_elem];
1188  M_vals_supg[i_elem][2][0][i_test][i_trial] = M_viscosity * integral_20 * M_detJacobian[i_elem];
1189  M_vals_supg[i_elem][2][1][i_test][i_trial] = M_viscosity * integral_21 * M_detJacobian[i_elem];
1190  M_vals_supg[i_elem][2][2][i_test][i_trial] = M_viscosity * integral_22 * M_detJacobian[i_elem];
1191  }
1192  }
1193 
1194  for ( dim_mat = 0; dim_mat < 3 ; dim_mat++ )
1195  {
1196  // DOF - test
1197  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
1198  {
1199  // DOF - trial
1200  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
1201  {
1202  M_vals_supg_01[i_elem][dim_mat][i_test][i_trial] = 0.0;
1203  M_cols_pressure[i_elem][i_trial] = M_elements_pressure[i_elem][i_trial];
1204  integral = 0.0;
1205  // QUAD
1206  for ( q = 0; q < NumQuadPoints ; q++ )
1207  {
1208  integral += dphi_phys_velocity[i_test][q][dim_mat]*M_phi_pressure[i_trial][q]*w_quad[q];
1209  }
1210  M_vals_supg_01[i_elem][dim_mat][i_test][i_trial] = -1.0 * integral * M_detJacobian[i_elem];
1211  }
1212  }
1213  }
1214 
1215  for ( dim_mat = 0; dim_mat < 3 ; dim_mat++ )
1216  {
1217  // DOF - test
1218  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
1219  {
1220  M_rows_pressure[i_elem][i_test] = M_elements_pressure[i_elem][i_test];
1221 
1222  // DOF - trial
1223  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
1224  {
1225  M_vals_supg_10[i_elem][dim_mat][i_test][i_trial] = 0.0;
1226  integral = 0.0;
1227  // QUAD
1228  for ( q = 0; q < NumQuadPoints ; q++ )
1229  {
1230  integral += dphi_phys_velocity[i_trial][q][dim_mat]*M_phi_pressure[i_test][q]*w_quad[q];
1231  }
1232 
1233  M_vals_supg_10[i_elem][dim_mat][i_test][i_trial] = integral * M_detJacobian[i_elem];
1234  }
1235  }
1236  }
1237  }
1238  }
1239 
1240  for ( int k = 0; k < M_numElements; ++k )
1241  {
1242  mass->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_velocity[k],
1243  ndof_velocity, M_cols_velocity[k],
1244  M_vals_00[k],
1245  Epetra_FECrsMatrix::ROW_MAJOR);
1246 
1247  b01->matrixPtr()->InsertGlobalValues ( ndof_velocity,
1248  M_rows_velocity[k],
1249  ndof_pressure,
1250  M_cols_pressure[k],
1251  M_vals_supg_01[k][0],
1252  Epetra_FECrsMatrix::ROW_MAJOR);
1253 
1254  b10->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k],
1255  ndof_velocity, M_cols_velocity[k],
1256  M_vals_supg_10[k][0],
1257  Epetra_FECrsMatrix::ROW_MAJOR);
1258 
1259  }
1260 
1261  for ( UInt d1 = 0; d1 < 3 ; d1++ ) // row index
1262  {
1263  for ( int k = 0; k < M_numElements; ++k )
1264  {
1265  for ( UInt i = 0; i < ndof_velocity; i++ )
1266  {
1267  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
1268  }
1269  stiffness->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1270  ndof_velocity, M_cols_velocity[k],
1271  M_vals_supg[k][d1][0],
1272  Epetra_FECrsMatrix::ROW_MAJOR);
1273 
1274  }
1275 
1276  for ( int k = 0; k < M_numElements; ++k )
1277  {
1278  for ( UInt i = 0; i < ndof_velocity; i++ )
1279  {
1280  M_cols_tmp[k][i] = M_cols_velocity[k][i] + M_numScalarDofs;
1281  }
1282  stiffness->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1283  ndof_velocity, M_cols_tmp[k],
1284  M_vals_supg[k][d1][1],
1285  Epetra_FECrsMatrix::ROW_MAJOR);
1286  }
1287 
1288  for ( int k = 0; k < M_numElements; ++k )
1289  {
1290  for ( UInt i = 0; i < ndof_velocity; i++ )
1291  {
1292  M_cols_tmp[k][i] = M_cols_velocity[k][i] + 2 * M_numScalarDofs;
1293  }
1294  stiffness->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1295  ndof_velocity, M_cols_tmp[k],
1296  M_vals_supg[k][d1][2],
1297  Epetra_FECrsMatrix::ROW_MAJOR);
1298  }
1299  }
1300 
1301  for ( UInt d1 = 1; d1 < 3 ; d1++ )
1302  {
1303  for ( int k = 0; k < M_numElements; ++k )
1304  {
1305  for ( UInt i = 0; i < ndof_velocity; i++ )
1306  {
1307  M_rows_velocity[k][i] += M_numScalarDofs;
1308  M_cols_velocity[k][i] += M_numScalarDofs;
1309  }
1310  mass->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_velocity[k],
1311  ndof_velocity, M_cols_velocity[k],
1312  M_vals_00[k],
1313  Epetra_FECrsMatrix::ROW_MAJOR);
1314 
1315  b01->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_velocity[k],
1316  ndof_pressure, M_cols_pressure[k],
1317  M_vals_supg_01[k][d1],
1318  Epetra_FECrsMatrix::ROW_MAJOR);
1319 
1320  b10->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k],
1321  ndof_velocity, M_cols_velocity[k],
1322  M_vals_supg_10[k][d1],
1323  Epetra_FECrsMatrix::ROW_MAJOR);
1324 
1325  }
1326  }
1327 }
1328 //=========================================================================
1329 void
1330 FastAssemblerNS::assembleConvective( matrixPtr_Type& matrix, const vector_Type& u_h )
1331 {
1332  int ndof_velocity = M_referenceFE_velocity->nbDof();
1333  int NumQuadPoints = M_qr->nbQuadPt();
1334  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
1335  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
1336 
1337  double w_quad[NumQuadPoints];
1338  for ( int q = 0; q < NumQuadPoints ; q++ )
1339  {
1340  w_quad[q] = M_qr->weight(q);
1341  }
1342 
1343  #pragma omp parallel firstprivate( w_quad, ndof_velocity, NumQuadPoints)
1344  {
1345  int i_elem, i_dof, q, d1, d2, i_test, i_trial, e_idof;
1346  double integral;
1347 
1348  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3];
1349 
1350  double uhq[3][NumQuadPoints];
1351 
1352  // ELEMENTI,
1353  #pragma omp for
1354  for ( i_elem = 0; i_elem < M_numElements; i_elem++ )
1355  {
1356  // DOF
1357  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1358  {
1359  // QUAD
1360  for ( q = 0; q < NumQuadPoints ; q++ )
1361  {
1362  // DIM 1
1363  for ( d1 = 0; d1 < 3 ; d1++ )
1364  {
1365  dphi_phys_velocity[i_dof][q][d1] = 0.0;
1366 
1367  // DIM 2
1368  for ( d2 = 0; d2 < 3 ; d2++ )
1369  {
1370  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
1371  }
1372  }
1373  }
1374  }
1375 
1376  // QUAD
1377  for ( q = 0; q < NumQuadPoints ; q++ )
1378  {
1379  for ( d1 = 0; d1 < 3 ; d1++ )
1380  {
1381  uhq[d1][q] = 0.0;
1382  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1383  {
1384  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs ;
1385  uhq[d1][q] += u_h[e_idof] * M_phi_velocity[i_dof][q];
1386  }
1387  }
1388  }
1389 
1390  // DOF - test
1391  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
1392  {
1393  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
1394 
1395  // DOF - trial
1396  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
1397  {
1398  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
1399 
1400  M_vals_00[i_elem][i_test][i_trial] = 0.0;
1401 
1402  integral = 0.0;
1403  // QUAD
1404  for ( q = 0; q < NumQuadPoints ; q++ )
1405  {
1406  // DIM 1
1407  for ( d1 = 0; d1 < 3 ; d1++ )
1408  {
1409  integral += uhq[d1][q] * dphi_phys_velocity[i_trial][q][d1] * M_phi_velocity[i_test][q] * w_quad[q];
1410  }
1411  }
1412  M_vals_00[i_elem][i_test][i_trial] = M_density * integral * M_detJacobian[i_elem];
1413  }
1414  }
1415  }
1416  }
1417 
1418  for ( int k = 0; k < M_numElements; ++k )
1419  {
1420  matrix->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_velocity[k],
1421  ndof_velocity, M_cols_velocity[k],
1422  M_vals_00[k],
1423  Epetra_FECrsMatrix::ROW_MAJOR);
1424  }
1425 
1426  for ( UInt d1 = 1; d1 < 3 ; d1++ )
1427  {
1428  for ( int k = 0; k < M_numElements; ++k )
1429  {
1430  for ( UInt i = 0; i < ndof_velocity; i++ )
1431  {
1432  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
1433  M_cols_tmp[k][i] = M_cols_velocity[k][i] + d1 * M_numScalarDofs;
1434  }
1435  matrix->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1436  ndof_velocity, M_cols_tmp[k],
1437  M_vals_00[k],
1438  Epetra_FECrsMatrix::ROW_MAJOR);
1439  }
1440  }
1441 }
1442 //=========================================================================
1443 void
1444 FastAssemblerNS::jacobianNS( matrixPtr_Type& matrix, const vector_Type& u_h )
1445 {
1446  int ndof_velocity = M_referenceFE_velocity->nbDof();
1447  int NumQuadPoints = M_qr->nbQuadPt();
1448  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
1449  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
1450 
1451  double w_quad[NumQuadPoints];
1452  for ( int q = 0; q < NumQuadPoints ; q++ )
1453  {
1454  w_quad[q] = M_qr->weight(q);
1455  }
1456 
1457  #pragma omp parallel firstprivate( w_quad, ndof_velocity, NumQuadPoints)
1458  {
1459  int i_elem, i_dof, q, d1, d2, i_test, i_trial, e_idof;
1460  double integral, integral_00, integral_01, integral_02, integral_10, integral_11, integral_12;
1461  double integral_20, integral_21, integral_22;
1462 
1463  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3];
1464 
1465  double duhq[3][3][NumQuadPoints];
1466 
1467  // ELEMENTI,
1468  #pragma omp for
1469  for ( i_elem = 0; i_elem < M_numElements; i_elem++ )
1470  {
1471 
1472  // DOF
1473  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1474  {
1475  // QUAD
1476  for ( q = 0; q < NumQuadPoints ; q++ )
1477  {
1478  // DIM 1
1479  for ( d1 = 0; d1 < 3 ; d1++ )
1480  {
1481  dphi_phys_velocity[i_dof][q][d1] = 0.0;
1482 
1483  // DIM 2
1484  for ( d2 = 0; d2 < 3 ; d2++ )
1485  {
1486  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
1487  }
1488  }
1489  }
1490  }
1491 
1492  // QUAD
1493  for ( q = 0; q < NumQuadPoints ; q++ )
1494  {
1495  for ( d1 = 0; d1 < 3 ; d1++ )
1496  {
1497  for ( d2 = 0; d2 < 3 ; d2++ )
1498  {
1499  duhq[d1][d2][q] = 0.0;
1500  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1501  {
1502  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs ;
1503  duhq[d1][d2][q] += u_h[e_idof] * dphi_phys_velocity[i_dof][q][d2];
1504  }
1505  }
1506  }
1507  }
1508 
1509  // DOF - test
1510  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
1511  {
1512  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
1513  // DOF - trial
1514  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
1515  {
1516  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
1517 
1518  M_vals_supg[i_elem][0][0][i_test][i_trial] = 0.0;
1519  M_vals_supg[i_elem][0][1][i_test][i_trial] = 0.0;
1520  M_vals_supg[i_elem][0][2][i_test][i_trial] = 0.0;
1521  M_vals_supg[i_elem][1][0][i_test][i_trial] = 0.0;
1522  M_vals_supg[i_elem][1][1][i_test][i_trial] = 0.0;
1523  M_vals_supg[i_elem][1][2][i_test][i_trial] = 0.0;
1524  M_vals_supg[i_elem][2][0][i_test][i_trial] = 0.0;
1525  M_vals_supg[i_elem][2][1][i_test][i_trial] = 0.0;
1526  M_vals_supg[i_elem][2][2][i_test][i_trial] = 0.0;
1527  integral_00 = 0.0;
1528  integral_01 = 0.0;
1529  integral_02 = 0.0;
1530  integral_10 = 0.0;
1531  integral_11 = 0.0;
1532  integral_12 = 0.0;
1533  integral_20 = 0.0;
1534  integral_21 = 0.0;
1535  integral_22 = 0.0;
1536  // QUAD
1537  for ( q = 0; q < NumQuadPoints ; q++ )
1538  {
1539  integral_00 += M_phi_velocity[i_test][q] * duhq[0][0][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1540 
1541  integral_01 += M_phi_velocity[i_test][q] * duhq[0][1][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1542 
1543  integral_02 += M_phi_velocity[i_test][q] * duhq[0][2][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1544 
1545  integral_10 += M_phi_velocity[i_test][q] * duhq[1][0][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1546 
1547  integral_11 += M_phi_velocity[i_test][q] * duhq[1][1][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1548 
1549  integral_12 += M_phi_velocity[i_test][q] * duhq[1][2][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1550 
1551  integral_20 += M_phi_velocity[i_test][q] * duhq[2][0][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1552 
1553  integral_21 += M_phi_velocity[i_test][q] * duhq[2][1][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1554 
1555  integral_22 += M_phi_velocity[i_test][q] * duhq[2][2][q] * M_phi_velocity[i_trial][q] * w_quad[q];
1556 
1557  }
1558  M_vals_supg[i_elem][0][0][i_test][i_trial] = M_density * integral_00 * M_detJacobian[i_elem];
1559  M_vals_supg[i_elem][0][1][i_test][i_trial] = M_density * integral_01 * M_detJacobian[i_elem];
1560  M_vals_supg[i_elem][0][2][i_test][i_trial] = M_density * integral_02 * M_detJacobian[i_elem];
1561  M_vals_supg[i_elem][1][0][i_test][i_trial] = M_density * integral_10 * M_detJacobian[i_elem];
1562  M_vals_supg[i_elem][1][1][i_test][i_trial] = M_density * integral_11 * M_detJacobian[i_elem];
1563  M_vals_supg[i_elem][1][2][i_test][i_trial] = M_density * integral_12 * M_detJacobian[i_elem];
1564  M_vals_supg[i_elem][2][0][i_test][i_trial] = M_density * integral_20 * M_detJacobian[i_elem];
1565  M_vals_supg[i_elem][2][1][i_test][i_trial] = M_density * integral_21 * M_detJacobian[i_elem];
1566  M_vals_supg[i_elem][2][2][i_test][i_trial] = M_density * integral_22 * M_detJacobian[i_elem];
1567  }
1568  }
1569  }
1570  }
1571 
1572  for ( UInt d1 = 0; d1 < 3 ; d1++ ) // row index
1573  {
1574  for ( int k = 0; k < M_numElements; ++k )
1575  {
1576  for ( UInt i = 0; i < ndof_velocity; i++ )
1577  {
1578  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
1579  }
1580  matrix->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1581  ndof_velocity, M_cols_velocity[k],
1582  M_vals_supg[k][d1][0],
1583  Epetra_FECrsMatrix::ROW_MAJOR);
1584 
1585  }
1586 
1587  for ( int k = 0; k < M_numElements; ++k )
1588  {
1589  for ( UInt i = 0; i < ndof_velocity; i++ )
1590  {
1591  M_cols_tmp[k][i] = M_cols_velocity[k][i] + M_numScalarDofs;
1592  }
1593  matrix->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1594  ndof_velocity, M_cols_tmp[k],
1595  M_vals_supg[k][d1][1],
1596  Epetra_FECrsMatrix::ROW_MAJOR);
1597  }
1598 
1599  for ( int k = 0; k < M_numElements; ++k )
1600  {
1601  for ( UInt i = 0; i < ndof_velocity; i++ )
1602  {
1603  M_cols_tmp[k][i] = M_cols_velocity[k][i] + 2 * M_numScalarDofs;
1604  }
1605  matrix->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
1606  ndof_velocity, M_cols_tmp[k],
1607  M_vals_supg[k][d1][2],
1608  Epetra_FECrsMatrix::ROW_MAJOR);
1609  }
1610  }
1611 }
1612 //=========================================================================
1613 void
1614 FastAssemblerNS::supg_FI_FSI_terms ( matrixPtr_Type& block00,
1615  matrixPtr_Type& block01,
1616  matrixPtr_Type& block10,
1617  matrixPtr_Type& block11,
1618  const vector_Type& beta_km1,
1619  const vector_Type& u_km1,
1620  const vector_Type& p_km1,
1621  const vector_Type& u_bdf)
1622 {
1623  int ndof_velocity = M_referenceFE_velocity->nbDof();
1624  int ndof_pressure = M_referenceFE_pressure->nbDof();
1625  int NumQuadPoints = M_qr->nbQuadPt();
1626  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
1627  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
1628 
1629  double w_quad[NumQuadPoints];
1630  for ( int q = 0; q < NumQuadPoints ; q++ )
1631  {
1632  w_quad[q] = M_qr->weight(q);
1633  }
1634 
1635 #pragma omp parallel firstprivate( w_quad, ndof_velocity, ndof_pressure, NumQuadPoints)
1636  {
1637  int i_elem, i_dof, q, d1, d2, i_test, i_trial, e_idof, dim_mat, iCoor, jCoor, k1, k2;
1638  double integral, integral_test, integral_trial, integral_partial, integral_lapl, partialSum;
1639 
1640  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3];
1641  double dphi_phys_pressure[ndof_pressure][NumQuadPoints][3];
1642  double d2phi_phys_velocity[ndof_velocity][NumQuadPoints][3][3];
1643  double u_km1_q[3][NumQuadPoints];
1644  double beta_km1_q[3][NumQuadPoints];
1645  double u_bdf_q[3][NumQuadPoints];
1646  double p_km1_q[NumQuadPoints];
1647  double i_00, i_01, i_02, i_10, i_11, i_12, i_20, i_21, i_22;
1648  double d_ukm1_q[3][3][NumQuadPoints];
1649  double d_pkm1_q[3][NumQuadPoints];
1650 
1651  // ELEMENTI,
1652 #pragma omp for
1653  for ( i_elem = 0; i_elem < M_numElements; i_elem++ )
1654  {
1655  // DOF
1656  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1657  {
1658  // QUAD
1659  for ( q = 0; q < NumQuadPoints ; q++ )
1660  {
1661  // DIM 1
1662  for ( d1 = 0; d1 < 3 ; d1++ )
1663  {
1664  dphi_phys_velocity[i_dof][q][d1] = 0.0;
1665 
1666  // DIM 2
1667  for ( d2 = 0; d2 < 3 ; d2++ )
1668  {
1669  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
1670  }
1671  }
1672  }
1673  }
1674 
1675  // DOF -- nota che può essere ottimizzato rishapando il gradiente fisico. Prima q poi d1 poi d2 e poi i_dof
1676  for ( i_dof = 0; i_dof < ndof_pressure; i_dof++ )
1677  {
1678  // QUAD
1679  for ( q = 0; q < NumQuadPoints ; q++ )
1680  {
1681  // DIM 1
1682  for ( d1 = 0; d1 < 3 ; d1++ )
1683  {
1684  dphi_phys_pressure[i_dof][q][d1] = 0.0;
1685 
1686  // DIM 2
1687  for ( d2 = 0; d2 < 3 ; d2++ )
1688  {
1689  dphi_phys_pressure[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_pressure[i_dof][q][d2];
1690  }
1691  }
1692  }
1693  }
1694 
1695  // QUAD
1696  for ( q = 0; q < NumQuadPoints ; q++ )
1697  {
1698  for ( d1 = 0; d1 < 3 ; d1++ )
1699  {
1700  beta_km1_q[d1][q] = 0.0;
1701  u_km1_q[d1][q] = 0.0;
1702  u_bdf_q[d1][q] = 0.0;
1703  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1704  {
1705  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs;
1706  beta_km1_q[d1][q] += beta_km1[e_idof] * M_phi_velocity[i_dof][q];
1707  u_km1_q[d1][q] += u_km1[e_idof] * M_phi_velocity[i_dof][q];
1708  u_bdf_q[d1][q] += u_bdf[e_idof] * M_phi_velocity[i_dof][q];
1709  }
1710  }
1711 
1712  p_km1_q[q] = 0.0;
1713  for ( i_dof = 0; i_dof < ndof_pressure; i_dof++ )
1714  {
1715  e_idof = M_elements_pressure[i_elem][i_dof];
1716  p_km1_q[q] += p_km1[e_idof] * M_phi_pressure[i_dof][q];
1717  }
1718 
1719  }
1720 
1721  // QUAD
1722  for ( q = 0; q < NumQuadPoints ; q++ )
1723  {
1724  for ( d1 = 0; d1 < 3 ; d1++ )
1725  {
1726  for ( d2 = 0; d2 < 3 ; d2++ )
1727  {
1728  d_ukm1_q[d1][d2][q] = 0.0;
1729  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
1730  {
1731  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs ;
1732  d_ukm1_q[d1][d2][q] += u_km1[e_idof] * dphi_phys_velocity[i_dof][q][d2];
1733  }
1734  }
1735  }
1736  }
1737 
1738  // QUAD
1739  for ( q = 0; q < NumQuadPoints ; q++ )
1740  {
1741  for ( d1 = 0; d1 < 3 ; d1++ )
1742  {
1743  d_pkm1_q[d1][q] = 0.0;
1744  for ( i_dof = 0; i_dof < ndof_pressure; i_dof++ )
1745  {
1746  e_idof = M_elements_pressure[i_elem][i_dof];
1747  d_pkm1_q[d1][q] += p_km1[e_idof] * dphi_phys_pressure[i_dof][q][d1];
1748  }
1749  }
1750  }
1751 
1752  // STABILIZZAZIONE - coefficienti Tau_M e Tau_C
1753  for ( q = 0; q < NumQuadPoints ; q++ )
1754  {
1755  M_Tau_M[i_elem][q] = 1.0/std::sqrt(
1756  M_density*M_density*M_orderBDF*M_orderBDF/(M_timestep*M_timestep) // TAU_M_DEN_DT
1757  +M_density*M_density*( beta_km1_q[0][q]*( M_G[i_elem][0][0]* beta_km1_q[0][q] + M_G[i_elem][0][1] * beta_km1_q[1][q] + M_G[i_elem][0][2] * beta_km1_q[2][q] ) +
1758  beta_km1_q[1][q]*( M_G[i_elem][1][0]* beta_km1_q[0][q] + M_G[i_elem][1][1] * beta_km1_q[1][q] + M_G[i_elem][1][2] * beta_km1_q[2][q] ) +
1759  beta_km1_q[2][q]*( M_G[i_elem][2][0]* beta_km1_q[0][q] + M_G[i_elem][2][1] * beta_km1_q[1][q] + M_G[i_elem][2][2] * beta_km1_q[2][q] )
1760  ) // TAU_M_DEN_VEL
1761  +M_C_I*M_viscosity*M_viscosity*(
1762  M_G[i_elem][0][0]*M_G[i_elem][0][0] + M_G[i_elem][0][1]*M_G[i_elem][0][1] + M_G[i_elem][0][2]*M_G[i_elem][0][2] +
1763  M_G[i_elem][1][0]*M_G[i_elem][1][0] + M_G[i_elem][1][1]*M_G[i_elem][1][1] + M_G[i_elem][1][2]*M_G[i_elem][1][2] +
1764  M_G[i_elem][2][0]*M_G[i_elem][2][0] + M_G[i_elem][2][1]*M_G[i_elem][2][1] + M_G[i_elem][2][2]*M_G[i_elem][2][2]
1765  ) // TAU_M_DEN_VISC
1766  );
1767 
1768  M_Tau_C[i_elem][q] = 1.0/( M_g[i_elem][0]*M_Tau_M[i_elem][q]*M_g[i_elem][0] +
1769  M_g[i_elem][1]*M_Tau_M[i_elem][q]*M_g[i_elem][1] +
1770  M_g[i_elem][2]*M_Tau_M[i_elem][q]*M_g[i_elem][2]
1771  );
1772 
1773  }
1774 
1775  // DOF - test
1776  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
1777  {
1778  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
1779 
1780  // DOF - trial
1781  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
1782  {
1783  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
1784 
1785  integral = 0.0;
1786  i_00 = 0.0; i_01 = 0.0; i_02 = 0.0; i_10 = 0.0; i_11 = 0.0; i_12 = 0.0; i_20 = 0.0; i_21 = 0.0; i_22 = 0.0;
1787 
1788  M_vals_supg[i_elem][0][0][i_test][i_trial] = 0.0;
1789  M_vals_supg[i_elem][0][1][i_test][i_trial] = 0.0;
1790  M_vals_supg[i_elem][0][2][i_test][i_trial] = 0.0;
1791  M_vals_supg[i_elem][1][0][i_test][i_trial] = 0.0;
1792  M_vals_supg[i_elem][1][1][i_test][i_trial] = 0.0;
1793  M_vals_supg[i_elem][1][2][i_test][i_trial] = 0.0;
1794  M_vals_supg[i_elem][2][0][i_test][i_trial] = 0.0;
1795  M_vals_supg[i_elem][2][1][i_test][i_trial] = 0.0;
1796  M_vals_supg[i_elem][2][2][i_test][i_trial] = 0.0;
1797 
1798  // QUAD
1799  for ( q = 0; q < NumQuadPoints ; q++ )
1800  {
1801  integral_test = 0.0;
1802  integral_trial = 0.0;
1803  integral_partial = 0.0;
1804 
1805  // DIM 1
1806  for ( d1 = 0; d1 < 3 ; d1++ )
1807  {
1808  integral_test += beta_km1_q[d1][q] * dphi_phys_velocity[i_test][q][d1];
1809  }
1810 
1811  integral_partial = ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] +
1812  dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] +
1813  dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q]
1814  ) *
1815  ( dphi_phys_velocity[i_trial][q][0] * beta_km1_q[0][q] +
1816  dphi_phys_velocity[i_trial][q][1] * beta_km1_q[1][q] +
1817  dphi_phys_velocity[i_trial][q][2] * beta_km1_q[2][q]
1818  );
1819 
1820  i_00 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] *
1821  ( M_density * M_alpha / M_timestep * u_km1_q[0][q] +
1822  d_pkm1_q[0][q]
1823  )
1824  -dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] * u_bdf_q[0][q]
1825  +M_density * M_density *
1826  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[0][0][q] * M_phi_velocity[i_trial][q]
1827  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[0][0][q] * M_phi_velocity[i_trial][q]
1828  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[0][0][q] * M_phi_velocity[i_trial][q]
1829  )
1830  + M_density * M_density * dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] *
1831  ( d_ukm1_q[0][0][q] * beta_km1_q[0][q] + d_ukm1_q[0][1][q] * beta_km1_q[1][q] + d_ukm1_q[0][2][q] * beta_km1_q[2][q] )
1832  )
1833  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][0] )
1834  ) * w_quad[q];
1835 
1836  i_01 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] *
1837  ( M_density * M_alpha / M_timestep * u_km1_q[0][q] +
1838  d_pkm1_q[0][q]
1839  )
1840  -dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] * u_bdf_q[0][q]
1841  +M_density * M_density *
1842  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[0][1][q] * M_phi_velocity[i_trial][q]
1843  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[0][1][q] * M_phi_velocity[i_trial][q]
1844  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[0][1][q] * M_phi_velocity[i_trial][q]
1845  )
1846  + M_density * M_density * dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] *
1847  ( d_ukm1_q[0][0][q] * beta_km1_q[0][q] + d_ukm1_q[0][1][q] * beta_km1_q[1][q] + d_ukm1_q[0][2][q] * beta_km1_q[2][q] )
1848  )
1849  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][1] )
1850  ) * w_quad[q];
1851 
1852  i_02 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] *
1853  ( M_density * M_alpha / M_timestep * u_km1_q[0][q] +
1854  d_pkm1_q[0][q]
1855  )
1856  -dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] * u_bdf_q[0][q]
1857  +M_density * M_density *
1858  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[0][2][q] * M_phi_velocity[i_trial][q]
1859  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[0][2][q] * M_phi_velocity[i_trial][q]
1860  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[0][2][q] * M_phi_velocity[i_trial][q]
1861  )
1862  + M_density * M_density * dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] *
1863  ( d_ukm1_q[0][0][q] * beta_km1_q[0][q] + d_ukm1_q[0][1][q] * beta_km1_q[1][q] + d_ukm1_q[0][2][q] * beta_km1_q[2][q] )
1864  )
1865  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][2] )
1866  ) * w_quad[q];
1867 
1868  i_10 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] *
1869  ( M_density * M_alpha / M_timestep * u_km1_q[1][q] +
1870  d_pkm1_q[1][q]
1871  )
1872  -dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] * u_bdf_q[1][q]
1873  +M_density * M_density *
1874  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[1][0][q] * M_phi_velocity[i_trial][q]
1875  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[1][0][q] * M_phi_velocity[i_trial][q]
1876  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[1][0][q] * M_phi_velocity[i_trial][q]
1877  )
1878  + M_density * M_density * dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] *
1879  ( d_ukm1_q[1][0][q] * beta_km1_q[0][q] + d_ukm1_q[1][1][q] * beta_km1_q[1][q] + d_ukm1_q[1][2][q] * beta_km1_q[2][q] )
1880  )
1881  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][0] )
1882  ) * w_quad[q];
1883 
1884  i_11 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] *
1885  ( M_density * M_alpha / M_timestep * u_km1_q[1][q] +
1886  d_pkm1_q[1][q]
1887  )
1888  -dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] * u_bdf_q[1][q]
1889  +M_density * M_density *
1890  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[1][1][q] * M_phi_velocity[i_trial][q]
1891  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[1][1][q] * M_phi_velocity[i_trial][q]
1892  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[1][1][q] * M_phi_velocity[i_trial][q]
1893  )
1894  + M_density * M_density * dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] *
1895  ( d_ukm1_q[1][0][q] * beta_km1_q[0][q] + d_ukm1_q[1][1][q] * beta_km1_q[1][q] + d_ukm1_q[1][2][q] * beta_km1_q[2][q] )
1896  )
1897  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][1] )
1898  ) * w_quad[q];
1899 
1900  i_12 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] *
1901  ( M_density * M_alpha / M_timestep * u_km1_q[1][q] +
1902  d_pkm1_q[1][q]
1903  )
1904  -dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] * u_bdf_q[1][q]
1905  +M_density * M_density *
1906  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[1][2][q] * M_phi_velocity[i_trial][q]
1907  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[1][2][q] * M_phi_velocity[i_trial][q]
1908  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[1][2][q] * M_phi_velocity[i_trial][q]
1909  )
1910  + M_density * M_density * dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] *
1911  ( d_ukm1_q[1][0][q] * beta_km1_q[0][q] + d_ukm1_q[1][1][q] * beta_km1_q[1][q] + d_ukm1_q[1][2][q] * beta_km1_q[2][q] )
1912  )
1913  + M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][2] )
1914  ) * w_quad[q];
1915 
1916  i_20 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] *
1917  ( M_density * M_alpha / M_timestep * u_km1_q[2][q] +
1918  d_pkm1_q[2][q]
1919  )
1920  -dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] * u_bdf_q[2][q]
1921  +M_density * M_density *
1922  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[2][0][q] * M_phi_velocity[i_trial][q]
1923  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[2][0][q] * M_phi_velocity[i_trial][q]
1924  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[2][0][q] * M_phi_velocity[i_trial][q]
1925  )
1926  + M_density * M_density * dphi_phys_velocity[i_test][q][0] * M_phi_velocity[i_trial][q] *
1927  ( d_ukm1_q[2][0][q] * beta_km1_q[0][q] + d_ukm1_q[2][1][q] * beta_km1_q[1][q] + d_ukm1_q[2][2][q] * beta_km1_q[2][q] )
1928  )
1929  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][0] )
1930  ) * w_quad[q];
1931 
1932  i_21 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] *
1933  ( M_density * M_alpha / M_timestep * u_km1_q[2][q] +
1934  d_pkm1_q[2][q]
1935  )
1936  -dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] * u_bdf_q[2][q]
1937  +M_density * M_density *
1938  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[2][1][q] * M_phi_velocity[i_trial][q]
1939  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[2][1][q] * M_phi_velocity[i_trial][q]
1940  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[2][1][q] * M_phi_velocity[i_trial][q]
1941  )
1942  + M_density * M_density * dphi_phys_velocity[i_test][q][1] * M_phi_velocity[i_trial][q] *
1943  ( d_ukm1_q[2][0][q] * beta_km1_q[0][q] + d_ukm1_q[2][1][q] * beta_km1_q[1][q] + d_ukm1_q[2][2][q] * beta_km1_q[2][q] )
1944  )
1945  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][1] )
1946  ) * w_quad[q];
1947 
1948  i_22 += ( M_Tau_M[i_elem][q] * ( M_density * dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] *
1949  ( M_density * M_alpha / M_timestep * u_km1_q[2][q] +
1950  d_pkm1_q[2][q]
1951  )
1952  -dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] * u_bdf_q[2][q]
1953  +M_density * M_density *
1954  ( dphi_phys_velocity[i_test][q][0] * beta_km1_q[0][q] * d_ukm1_q[2][2][q] * M_phi_velocity[i_trial][q]
1955  +dphi_phys_velocity[i_test][q][1] * beta_km1_q[1][q] * d_ukm1_q[2][2][q] * M_phi_velocity[i_trial][q]
1956  +dphi_phys_velocity[i_test][q][2] * beta_km1_q[2][q] * d_ukm1_q[2][2][q] * M_phi_velocity[i_trial][q]
1957  )
1958  + M_density * M_density * dphi_phys_velocity[i_test][q][2] * M_phi_velocity[i_trial][q] *
1959  ( d_ukm1_q[2][0][q] * beta_km1_q[0][q] + d_ukm1_q[2][1][q] * beta_km1_q[1][q] + d_ukm1_q[2][2][q] * beta_km1_q[2][q] )
1960  )
1961  +M_Tau_C[i_elem][q] * ( dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][2] )
1962  ) * w_quad[q];
1963 
1964  integral += M_Tau_M[i_elem][q] * ( M_density * M_density * M_alpha / M_timestep * integral_test * M_phi_velocity[i_trial][q]
1965  + M_density * M_density * integral_partial
1966  ) * w_quad[q];
1967 
1968  }
1969 
1970  M_vals_supg[i_elem][0][0][i_test][i_trial] = (i_00 + integral) * M_detJacobian[i_elem];
1971  M_vals_supg[i_elem][0][1][i_test][i_trial] = (i_01 ) * M_detJacobian[i_elem];
1972  M_vals_supg[i_elem][0][2][i_test][i_trial] = (i_02 ) * M_detJacobian[i_elem];
1973  M_vals_supg[i_elem][1][0][i_test][i_trial] = (i_10 ) * M_detJacobian[i_elem];
1974  M_vals_supg[i_elem][1][1][i_test][i_trial] = (i_11 + integral) * M_detJacobian[i_elem];
1975  M_vals_supg[i_elem][1][2][i_test][i_trial] = (i_12 ) * M_detJacobian[i_elem];
1976  M_vals_supg[i_elem][2][0][i_test][i_trial] = (i_20 ) * M_detJacobian[i_elem];
1977  M_vals_supg[i_elem][2][1][i_test][i_trial] = (i_21 ) * M_detJacobian[i_elem];
1978  M_vals_supg[i_elem][2][2][i_test][i_trial] = (i_22 + integral) * M_detJacobian[i_elem];
1979 
1980  }
1981 
1982  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
1983  {
1984  for ( dim_mat = 0; dim_mat < 3 ; dim_mat++ )
1985  {
1986  M_cols_pressure[i_elem][i_trial] = M_elements_pressure[i_elem][i_trial];
1987  M_vals_supg_01[i_elem][dim_mat][i_test][i_trial] = 0.0;
1988 
1989  integral = 0.0;
1990 
1991  // QUAD
1992  for ( q = 0; q < NumQuadPoints ; q++ )
1993  {
1994  integral_partial = 0.0;
1995  for ( d1 = 0; d1 < 3 ; d1++ )
1996  {
1997  integral_partial += ( dphi_phys_velocity[i_test][q][d1] * beta_km1_q[d1][q] );
1998  }
1999  integral += M_Tau_M[i_elem][q] * dphi_phys_pressure[i_trial][q][dim_mat] * integral_partial * w_quad[q];
2000  }
2001  M_vals_supg_01[i_elem][dim_mat][i_test][i_trial] = integral * M_detJacobian[i_elem];
2002  }
2003  }
2004 
2005  }// end test
2006 
2007  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
2008  {
2009  M_rows_pressure[i_elem][i_test] = M_elements_pressure[i_elem][i_test];
2010 
2011  // DOF - trial
2012  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
2013  {
2014  M_vals_11[i_elem][i_test][i_trial] = 0.0;
2015  integral = 0.0;
2016  // QUAD
2017  for ( q = 0; q < NumQuadPoints ; q++ )
2018  {
2019  // DIM 1
2020  for ( d1 = 0; d1 < 3 ; d1++ )
2021  {
2022  integral += M_Tau_M[i_elem][q] * dphi_phys_pressure[i_test][q][d1] * dphi_phys_pressure[i_trial][q][d1]*w_quad[q];
2023  }
2024  }
2025  M_vals_11[i_elem][i_test][i_trial] = integral * M_detJacobian[i_elem];
2026  }
2027  }
2028 
2029  // DOF - test
2030  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
2031  {
2032  // DOF - trial
2033  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
2034  {
2035  i_00 = 0.0;
2036  i_10 = 0.0;
2037  i_20 = 0.0;
2038 
2039  M_vals_supg_10[i_elem][0][i_test][i_trial] = 0.0;
2040  M_vals_supg_10[i_elem][1][i_test][i_trial] = 0.0;
2041  M_vals_supg_10[i_elem][2][i_test][i_trial] = 0.0;
2042 
2043  // QUAD
2044  for ( q = 0; q < NumQuadPoints ; q++ )
2045  {
2046  i_00 += M_Tau_M[i_elem][q] * (
2047  M_alpha * M_density / M_timestep * dphi_phys_pressure[i_test][q][0] * M_phi_velocity[i_trial][q]
2048 
2049  +M_density * ( dphi_phys_pressure[i_test][q][0] * M_phi_velocity[i_trial][q] * d_ukm1_q[0][0][q]
2050  +dphi_phys_pressure[i_test][q][1] * M_phi_velocity[i_trial][q] * d_ukm1_q[1][0][q]
2051  +dphi_phys_pressure[i_test][q][2] * M_phi_velocity[i_trial][q] * d_ukm1_q[2][0][q]
2052  )
2053  +M_density * dphi_phys_pressure[i_test][q][0] * (
2054  dphi_phys_velocity[i_trial][q][0]*beta_km1_q[0][q]
2055  +dphi_phys_velocity[i_trial][q][1]*beta_km1_q[1][q]
2056  +dphi_phys_velocity[i_trial][q][2]*beta_km1_q[2][q]
2057  )
2058  ) * w_quad[q];
2059 
2060  i_10 += M_Tau_M[i_elem][q] * (
2061  M_alpha * M_density / M_timestep * dphi_phys_pressure[i_test][q][1] * M_phi_velocity[i_trial][q]
2062  +M_density * ( dphi_phys_pressure[i_test][q][0] * M_phi_velocity[i_trial][q] * d_ukm1_q[0][1][q]
2063  +dphi_phys_pressure[i_test][q][1] * M_phi_velocity[i_trial][q] * d_ukm1_q[1][1][q]
2064  +dphi_phys_pressure[i_test][q][2] * M_phi_velocity[i_trial][q] * d_ukm1_q[2][1][q]
2065  )
2066  +M_density * dphi_phys_pressure[i_test][q][1] * (
2067  dphi_phys_velocity[i_trial][q][0]*beta_km1_q[0][q]
2068  +dphi_phys_velocity[i_trial][q][1]*beta_km1_q[1][q]
2069  +dphi_phys_velocity[i_trial][q][2]*beta_km1_q[2][q]
2070  )
2071 
2072  ) * w_quad[q];
2073 
2074  i_20 += M_Tau_M[i_elem][q] * (
2075  M_alpha * M_density / M_timestep * dphi_phys_pressure[i_test][q][2] * M_phi_velocity[i_trial][q]
2076  +M_density * ( dphi_phys_pressure[i_test][q][0] * M_phi_velocity[i_trial][q] * d_ukm1_q[0][2][q]
2077  +dphi_phys_pressure[i_test][q][1] * M_phi_velocity[i_trial][q] * d_ukm1_q[1][2][q]
2078  +dphi_phys_pressure[i_test][q][2] * M_phi_velocity[i_trial][q] * d_ukm1_q[2][2][q]
2079  )
2080  +M_density * dphi_phys_pressure[i_test][q][2] * (
2081  dphi_phys_velocity[i_trial][q][0]*beta_km1_q[0][q]
2082  +dphi_phys_velocity[i_trial][q][1]*beta_km1_q[1][q]
2083  +dphi_phys_velocity[i_trial][q][2]*beta_km1_q[2][q]
2084  )
2085 
2086  ) * w_quad[q];
2087  }
2088 
2089  M_vals_supg_10[i_elem][0][i_test][i_trial] = i_00 * M_detJacobian[i_elem];
2090  M_vals_supg_10[i_elem][1][i_test][i_trial] = i_10 * M_detJacobian[i_elem];
2091  M_vals_supg_10[i_elem][2][i_test][i_trial] = i_20 * M_detJacobian[i_elem];
2092  }
2093  }
2094 
2095  }// end elements
2096  }// end parallel region
2097 
2098  for ( UInt d1 = 0; d1 < 3 ; d1++ ) // row index
2099  {
2100  for ( int k = 0; k < M_numElements; ++k )
2101  {
2102  for ( UInt i = 0; i < ndof_velocity; i++ )
2103  {
2104  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
2105  }
2106  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2107  ndof_velocity, M_cols_velocity[k],
2108  M_vals_supg[k][d1][0],
2109  Epetra_FECrsMatrix::ROW_MAJOR);
2110  block01->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2111  ndof_pressure, M_cols_pressure[k],
2112  M_vals_supg_01[k][d1],
2113  Epetra_FECrsMatrix::ROW_MAJOR);
2114  }
2115 
2116  for ( int k = 0; k < M_numElements; ++k )
2117  {
2118  for ( UInt i = 0; i < ndof_velocity; i++ )
2119  {
2120  M_cols_tmp[k][i] = M_cols_velocity[k][i] + M_numScalarDofs;
2121  }
2122  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2123  ndof_velocity, M_cols_tmp[k],
2124  M_vals_supg[k][d1][1],
2125  Epetra_FECrsMatrix::ROW_MAJOR);
2126  }
2127 
2128  for ( int k = 0; k < M_numElements; ++k )
2129  {
2130  for ( UInt i = 0; i < ndof_velocity; i++ )
2131  {
2132  M_cols_tmp[k][i] = M_cols_velocity[k][i] + 2 * M_numScalarDofs;
2133  }
2134  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2135  ndof_velocity, M_cols_tmp[k],
2136  M_vals_supg[k][d1][2],
2137  Epetra_FECrsMatrix::ROW_MAJOR);
2138  }
2139  }
2140 
2141  for ( int k = 0; k < M_numElements; ++k )
2142  {
2143  block11->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k],
2144  ndof_pressure, M_cols_pressure[k],
2145  M_vals_11[k],
2146  Epetra_FECrsMatrix::ROW_MAJOR);
2147  }
2148 
2149  for ( UInt d1 = 0; d1 < 3 ; d1++ )
2150  {
2151  for ( int k = 0; k < M_numElements; ++k )
2152  {
2153  for ( UInt i = 0; i < ndof_velocity; i++ )
2154  {
2155  M_cols_tmp[k][i] = M_cols_velocity[k][i] + d1 * M_numScalarDofs;
2156  }
2157  block10->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k],
2158  ndof_velocity, M_cols_tmp[k],
2159  M_vals_supg_10[k][d1],
2160  Epetra_FECrsMatrix::ROW_MAJOR);
2161  }
2162  }
2163 }
2164 
2165 void
2166 FastAssemblerNS::vmsles_semi_implicit_terms ( matrixPtr_Type& block00,
2167  matrixPtr_Type& block01,
2168  matrixPtr_Type& block10,
2169  matrixPtr_Type& block11,
2170  const std::vector<std::vector<VectorSmall<3>>>& fine_scale,
2171  const vector_Type& u_extr )
2172 {
2173  int ndof_velocity = M_referenceFE_velocity->nbDof();
2174  int ndof_pressure = M_referenceFE_pressure->nbDof();
2175  int NumQuadPoints = M_qr->nbQuadPt();
2176  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
2177  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
2178 
2179  double w_quad[NumQuadPoints];
2180  for ( int q = 0; q < NumQuadPoints ; q++ )
2181  {
2182  w_quad[q] = M_qr->weight(q);
2183  }
2184 
2185  #pragma omp parallel firstprivate( w_quad, ndof_velocity, ndof_pressure, NumQuadPoints)
2186  {
2187  int i_elem, i_dof, q, d1, d2, i_test, i_trial, e_idof, dim_mat, iCoor, jCoor, k1, k2;
2188  double integral, integral_test, integral_trial, integral_partial, integral_lapl, partialSum;
2189 
2190  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3];
2191  double dphi_phys_pressure[ndof_pressure][NumQuadPoints][3];
2192  double d2phi_phys_velocity[ndof_velocity][NumQuadPoints][3][3];
2193  double uhq[3][NumQuadPoints];
2194  double i_00, i_01, i_02, i_10, i_11, i_12, i_20, i_21, i_22;
2195 
2196 
2197  // ELEMENTI,
2198  #pragma omp for
2199  for ( i_elem = 0; i_elem < M_numElements; i_elem++ )
2200  {
2201  // DOF
2202  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
2203  {
2204  // QUAD
2205  for ( q = 0; q < NumQuadPoints ; q++ )
2206  {
2207  // DIM 1
2208  for ( d1 = 0; d1 < 3 ; d1++ )
2209  {
2210  dphi_phys_velocity[i_dof][q][d1] = 0.0;
2211 
2212  // DIM 2
2213  for ( d2 = 0; d2 < 3 ; d2++ )
2214  {
2215  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
2216  }
2217  }
2218  }
2219  }
2220 
2221  for ( i_dof = 0; i_dof < ndof_pressure; i_dof++ )
2222  {
2223  // QUAD
2224  for ( q = 0; q < NumQuadPoints ; q++ )
2225  {
2226  // DIM 1
2227  for ( d1 = 0; d1 < 3 ; d1++ )
2228  {
2229  dphi_phys_pressure[i_dof][q][d1] = 0.0;
2230 
2231  // DIM 2
2232  for ( d2 = 0; d2 < 3 ; d2++ )
2233  {
2234  dphi_phys_pressure[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_pressure[i_dof][q][d2];
2235  }
2236  }
2237  }
2238  }
2239 
2240  // QUAD
2241  for ( q = 0; q < NumQuadPoints ; q++ )
2242  {
2243  // DOF
2244  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
2245  {
2246  // DIM 1
2247  for ( iCoor = 0; iCoor < 3 ; iCoor++ )
2248  {
2249  // DIM 2
2250  for ( jCoor = 0; jCoor < 3 ; jCoor++ )
2251  {
2252  partialSum = 0.0;
2253  // DIM 1
2254  for ( k1 = 0; k1 < 3 ; k1++ )
2255  {
2256  // DIM 2
2257  for ( k2 = 0; k2 < 3 ; k2++ )
2258  {
2259  partialSum += M_invJacobian[i_elem][iCoor][k1]
2260  * M_d2phi_velocity[i_dof][q][k1][k2]
2261  * M_invJacobian[i_elem][jCoor][k2];
2262  }
2263  }
2264  d2phi_phys_velocity[i_dof][q][iCoor][jCoor] = partialSum;
2265  }
2266  }
2267  }
2268  }
2269 
2270  // QUAD
2271  for ( q = 0; q < NumQuadPoints ; q++ )
2272  {
2273  for ( d1 = 0; d1 < 3 ; d1++ )
2274  {
2275  uhq[d1][q] = 0.0;
2276  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
2277  {
2278  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs ;
2279  uhq[d1][q] += u_extr[e_idof] * M_phi_velocity[i_dof][q];
2280  }
2281  }
2282  }
2283 
2284  // STABILIZZAZIONE - coefficienti Tau_M e Tau_C
2285  for ( q = 0; q < NumQuadPoints ; q++ )
2286  {
2287  M_Tau_M[i_elem][q] = 1.0/std::sqrt(
2288  M_density*M_density*( uhq[0][q]*( M_G[i_elem][0][0]* uhq[0][q] + M_G[i_elem][0][1] * uhq[1][q] + M_G[i_elem][0][2] * uhq[2][q] ) +
2289  uhq[1][q]*( M_G[i_elem][1][0]* uhq[0][q] + M_G[i_elem][1][1] * uhq[1][q] + M_G[i_elem][1][2] * uhq[2][q] ) +
2290  uhq[2][q]*( M_G[i_elem][2][0]* uhq[0][q] + M_G[i_elem][2][1] * uhq[1][q] + M_G[i_elem][2][2] * uhq[2][q] )
2291  ) // TAU_M_DEN_VEL
2292  +M_C_I*M_viscosity*M_viscosity*(
2293  M_G[i_elem][0][0]*M_G[i_elem][0][0] + M_G[i_elem][0][1]*M_G[i_elem][0][1] + M_G[i_elem][0][2]*M_G[i_elem][0][2] +
2294  M_G[i_elem][1][0]*M_G[i_elem][1][0] + M_G[i_elem][1][1]*M_G[i_elem][1][1] + M_G[i_elem][1][2]*M_G[i_elem][1][2] +
2295  M_G[i_elem][2][0]*M_G[i_elem][2][0] + M_G[i_elem][2][1]*M_G[i_elem][2][1] + M_G[i_elem][2][2]*M_G[i_elem][2][2]
2296  ) // TAU_M_DEN_VISC*/
2297  );
2298 
2299  M_Tau_C[i_elem][q] = 1.0/( M_g[i_elem][0]*M_Tau_M[i_elem][q]*M_g[i_elem][0] +
2300  M_g[i_elem][1]*M_Tau_M[i_elem][q]*M_g[i_elem][1] +
2301  M_g[i_elem][2]*M_Tau_M[i_elem][q]*M_g[i_elem][2]
2302  );
2303 
2304  M_Tau_M_hat[i_elem][q] = 1.0/( M_alpha * M_density / M_timestep + 1.0/M_Tau_M[i_elem][q]);
2305  }
2306 
2307  // DOF - test
2308  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
2309  {
2310  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
2311 
2312  // DOF - trial
2313  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
2314  {
2315  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
2316 
2317  i_00 = 0.0; i_01 = 0.0; i_02 = 0.0; i_10 = 0.0; i_11 = 0.0; i_12 = 0.0; i_20 = 0.0; i_21 = 0.0; i_22 = 0.0;
2318 
2319  M_vals_supg[i_elem][0][0][i_test][i_trial] = 0.0;
2320  M_vals_supg[i_elem][0][1][i_test][i_trial] = 0.0;
2321  M_vals_supg[i_elem][0][2][i_test][i_trial] = 0.0;
2322  M_vals_supg[i_elem][1][0][i_test][i_trial] = 0.0;
2323  M_vals_supg[i_elem][1][1][i_test][i_trial] = 0.0;
2324  M_vals_supg[i_elem][1][2][i_test][i_trial] = 0.0;
2325  M_vals_supg[i_elem][2][0][i_test][i_trial] = 0.0;
2326  M_vals_supg[i_elem][2][1][i_test][i_trial] = 0.0;
2327  M_vals_supg[i_elem][2][2][i_test][i_trial] = 0.0;
2328 
2329  double tmp_supg = 0.0;
2330  double tmp_supg_test = 0.0;
2331  double tmp_supg_trial = 0.0;
2332 
2333  // QUAD
2334  for ( q = 0; q < NumQuadPoints ; q++ )
2335  {
2336  tmp_supg = M_alpha * M_density * M_density / M_timestep * (
2337  dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q]
2338  ) * M_phi_velocity[i_trial][q];
2339 
2340  tmp_supg_test =
2341  dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q];
2342 
2343  tmp_supg_trial =
2344  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q];
2345 
2346  i_00 += ( M_Tau_M_hat[i_elem][q] * (
2347  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] * M_phi_velocity[i_trial][q] +
2348  M_density * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] * (
2349  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2350  )
2351  -M_viscosity * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] *
2352  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2353  // VMS
2354  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][0] * uhq[0][q] * M_phi_velocity[i_trial][q] )
2355  +M_density * M_density * dphi_phys_velocity[i_test][q][0] * uhq[0][q] * (
2356  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2357  )
2358  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][0] * uhq[0][q] * (
2359  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2360  )
2361  // SUPG
2362  +tmp_supg
2363  +M_density * M_density * tmp_supg_test*tmp_supg_trial
2364  -M_density * M_viscosity * tmp_supg_test * (
2365  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2366  )
2367  )
2368  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][0]
2369 
2370  ) * w_quad[q];
2371 
2372  i_01 += ( M_Tau_M_hat[i_elem][q] * (
2373  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] * M_phi_velocity[i_trial][q] +
2374  M_density * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] * (
2375  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2376  )
2377  -M_viscosity * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] *
2378  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2379  // VMS
2380  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][1] * uhq[0][q] * M_phi_velocity[i_trial][q] )
2381  +M_density * M_density * dphi_phys_velocity[i_test][q][1] * uhq[0][q] * (
2382  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2383  )
2384  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][1] * uhq[0][q] * (
2385  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2386  )
2387  )
2388  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][1]
2389  ) * w_quad[q];
2390 
2391  i_02 += ( M_Tau_M_hat[i_elem][q] * (
2392  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] * M_phi_velocity[i_trial][q] +
2393  M_density * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] * (
2394  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2395  )
2396  -M_viscosity * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] *
2397  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2398  // VMS
2399  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][2] * uhq[0][q] * M_phi_velocity[i_trial][q] )
2400  +M_density * M_density * dphi_phys_velocity[i_test][q][2] * uhq[0][q] * (
2401  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2402  )
2403  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][2] * uhq[0][q] * (
2404  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2405  )
2406  )
2407  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][2]
2408  ) * w_quad[q];
2409 
2410  i_10 += ( M_Tau_M_hat[i_elem][q] * (
2411  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] * M_phi_velocity[i_trial][q] +
2412  M_density * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] * (
2413  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2414  )
2415  -M_viscosity * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] *
2416  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2417  // VMS
2418  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][0] * uhq[1][q] * M_phi_velocity[i_trial][q] )
2419  +M_density * M_density * dphi_phys_velocity[i_test][q][0] * uhq[1][q] * (
2420  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2421  )
2422  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][0] * uhq[1][q] * (
2423  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2424  )
2425  )
2426  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][0]
2427  ) * w_quad[q];
2428 
2429  i_11 += ( M_Tau_M_hat[i_elem][q] * (
2430  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] * M_phi_velocity[i_trial][q] +
2431  M_density * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] * (
2432  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2433  )
2434  -M_viscosity * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] *
2435  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2436  // VMS
2437  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][1] * uhq[1][q] * M_phi_velocity[i_trial][q] )
2438  +M_density * M_density * dphi_phys_velocity[i_test][q][1] * uhq[1][q] * (
2439  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2440  )
2441  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][1] * uhq[1][q] * (
2442  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2443  )
2444  // SUPG
2445  +tmp_supg
2446  +M_density * M_density * tmp_supg_test*tmp_supg_trial
2447  -M_density * M_viscosity * tmp_supg_test * (
2448  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2449  )
2450  )
2451  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][1]
2452  ) * w_quad[q];
2453 
2454  i_12 += ( M_Tau_M_hat[i_elem][q] * (
2455  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] * M_phi_velocity[i_trial][q] +
2456  M_density * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] * (
2457  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2458  )
2459  -M_viscosity * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] *
2460  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2461  // VMS
2462  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][2] * uhq[1][q] * M_phi_velocity[i_trial][q] )
2463  +M_density * M_density * dphi_phys_velocity[i_test][q][2] * uhq[1][q] * (
2464  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2465  )
2466  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][2] * uhq[1][q] * (
2467  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2468  )
2469  )
2470  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][2]
2471  ) * w_quad[q];
2472 
2473  i_20 += ( M_Tau_M_hat[i_elem][q] * (
2474  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] * M_phi_velocity[i_trial][q] +
2475  M_density * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] * (
2476  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2477  )
2478  -M_viscosity * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] *
2479  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2480  // VMS
2481  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][0] * uhq[2][q] * M_phi_velocity[i_trial][q] )
2482  +M_density * M_density * dphi_phys_velocity[i_test][q][0] * uhq[2][q] * (
2483  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2484  )
2485  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][0] * uhq[2][q] * (
2486  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2487  )
2488  )
2489  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][0]
2490  ) * w_quad[q];
2491 
2492  i_21 += ( M_Tau_M_hat[i_elem][q] * (
2493  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] * M_phi_velocity[i_trial][q] +
2494  M_density * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] * (
2495  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2496  )
2497  -M_viscosity * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] *
2498  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2499  // VMS
2500  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][1] * uhq[2][q] * M_phi_velocity[i_trial][q] )
2501  +M_density * M_density * dphi_phys_velocity[i_test][q][1] * uhq[2][q] * (
2502  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2503  )
2504  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][1] * uhq[2][q] * (
2505  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2506  )
2507  )
2508  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][1]
2509  ) * w_quad[q];
2510 
2511  i_22 += ( M_Tau_M_hat[i_elem][q] * (
2512  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] * M_phi_velocity[i_trial][q] +
2513  M_density * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] * (
2514  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2515  )
2516  -M_viscosity * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] *
2517  ( d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2] )
2518  // VMS
2519  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][2] * uhq[2][q] * M_phi_velocity[i_trial][q] )
2520  +M_density * M_density * dphi_phys_velocity[i_test][q][2] * uhq[2][q] * (
2521  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2522  )
2523  -M_density * M_viscosity * dphi_phys_velocity[i_test][q][2] * uhq[2][q] * (
2524  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2525  )
2526  // SUPG
2527  +tmp_supg
2528  +M_density * M_density * tmp_supg_test*tmp_supg_trial
2529  -M_density * M_viscosity * tmp_supg_test * (
2530  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2531  )
2532  )
2533  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][2]
2534  ) * w_quad[q];
2535 
2536 
2537  } // end quadrature
2538 
2539  M_vals_supg[i_elem][0][0][i_test][i_trial] = (i_00 ) * M_detJacobian[i_elem];
2540  M_vals_supg[i_elem][0][1][i_test][i_trial] = (i_01 ) * M_detJacobian[i_elem];
2541  M_vals_supg[i_elem][0][2][i_test][i_trial] = (i_02 ) * M_detJacobian[i_elem];
2542  M_vals_supg[i_elem][1][0][i_test][i_trial] = (i_10 ) * M_detJacobian[i_elem];
2543  M_vals_supg[i_elem][1][1][i_test][i_trial] = (i_11 ) * M_detJacobian[i_elem];
2544  M_vals_supg[i_elem][1][2][i_test][i_trial] = (i_12 ) * M_detJacobian[i_elem];
2545  M_vals_supg[i_elem][2][0][i_test][i_trial] = (i_20 ) * M_detJacobian[i_elem];
2546  M_vals_supg[i_elem][2][1][i_test][i_trial] = (i_21 ) * M_detJacobian[i_elem];
2547  M_vals_supg[i_elem][2][2][i_test][i_trial] = (i_22 ) * M_detJacobian[i_elem];
2548 
2549  } // end trial
2550 
2551  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
2552  {
2553  M_cols_pressure[i_elem][i_trial] = M_elements_pressure[i_elem][i_trial];
2554  M_vals_supg_01[i_elem][0][i_test][i_trial] = 0.0;
2555  M_vals_supg_01[i_elem][1][i_test][i_trial] = 0.0;
2556  M_vals_supg_01[i_elem][2][i_test][i_trial] = 0.0;
2557 
2558  i_00 = 0.0;
2559  i_10 = 0.0;
2560  i_20 = 0.0;
2561 
2562  // QUAD
2563  for ( q = 0; q < NumQuadPoints ; q++ )
2564  {
2565  i_00 += M_Tau_M_hat[i_elem][q] * (
2566  dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] * dphi_phys_pressure[i_trial][q][0]
2567  +dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] * dphi_phys_pressure[i_trial][q][1]
2568  +dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] * dphi_phys_pressure[i_trial][q][2]
2569  +M_density * (
2570  (dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q])
2571  * dphi_phys_pressure[i_trial][q][0]
2572  + dphi_phys_velocity[i_test][q][0] * uhq[0][q] * dphi_phys_pressure[i_trial][q][0]
2573  + dphi_phys_velocity[i_test][q][1] * uhq[0][q] * dphi_phys_pressure[i_trial][q][1]
2574  + dphi_phys_velocity[i_test][q][2] * uhq[0][q] * dphi_phys_pressure[i_trial][q][2]
2575  )
2576  ) * w_quad[q];
2577 
2578  i_10 += M_Tau_M_hat[i_elem][q] * (
2579  dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] * dphi_phys_pressure[i_trial][q][0]
2580  +dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] * dphi_phys_pressure[i_trial][q][1]
2581  +dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] * dphi_phys_pressure[i_trial][q][2]
2582  +M_density * (
2583  (dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q])
2584  * dphi_phys_pressure[i_trial][q][1]
2585  + dphi_phys_velocity[i_test][q][0] * uhq[1][q] * dphi_phys_pressure[i_trial][q][0]
2586  + dphi_phys_velocity[i_test][q][1] * uhq[1][q] * dphi_phys_pressure[i_trial][q][1]
2587  + dphi_phys_velocity[i_test][q][2] * uhq[1][q] * dphi_phys_pressure[i_trial][q][2]
2588  )
2589  ) * w_quad[q];
2590 
2591  i_20 += M_Tau_M_hat[i_elem][q] * (
2592  dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] * dphi_phys_pressure[i_trial][q][0]
2593  +dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] * dphi_phys_pressure[i_trial][q][1]
2594  +dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] * dphi_phys_pressure[i_trial][q][2]
2595  +M_density * (
2596  (dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q])
2597  * dphi_phys_pressure[i_trial][q][2]
2598  + dphi_phys_velocity[i_test][q][0] * uhq[2][q] * dphi_phys_pressure[i_trial][q][0]
2599  + dphi_phys_velocity[i_test][q][1] * uhq[2][q] * dphi_phys_pressure[i_trial][q][1]
2600  + dphi_phys_velocity[i_test][q][2] * uhq[2][q] * dphi_phys_pressure[i_trial][q][2]
2601  )
2602  ) * w_quad[q];
2603  }
2604 
2605  M_vals_supg_01[i_elem][0][i_test][i_trial] = i_00 * M_detJacobian[i_elem];
2606  M_vals_supg_01[i_elem][1][i_test][i_trial] = i_10 * M_detJacobian[i_elem];
2607  M_vals_supg_01[i_elem][2][i_test][i_trial] = i_20 * M_detJacobian[i_elem];
2608 
2609  }
2610 
2611  } // end dof test
2612 
2613  // DOF - test - assemblo blocchi (1,0) e (1,1)
2614  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
2615  {
2616  M_rows_pressure[i_elem][i_test] = M_elements_pressure[i_elem][i_test];
2617 
2618  // DOF - trial
2619  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
2620  {
2621  M_vals_11[i_elem][i_test][i_trial] = 0.0;
2622 
2623  integral = 0.0;
2624  // QUAD
2625  for ( q = 0; q < NumQuadPoints ; q++ )
2626  {
2627  // DIM 1
2628  for ( d1 = 0; d1 < 3 ; d1++ )
2629  {
2630  integral += M_Tau_M_hat[i_elem][q] * dphi_phys_pressure[i_test][q][d1] * dphi_phys_pressure[i_trial][q][d1]*w_quad[q];
2631  }
2632  }
2633  M_vals_11[i_elem][i_test][i_trial] = integral * M_detJacobian[i_elem];
2634  }
2635  }
2636 
2637  // DOF - test
2638  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
2639  {
2640  // DOF - trial
2641  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
2642  {
2643  i_00 = 0.0;
2644  i_01 = 0.0;
2645  i_02 = 0.0;
2646 
2647  // QUAD
2648  for ( q = 0; q < NumQuadPoints ; q++ )
2649  {
2650  i_00 += M_Tau_M_hat[i_elem][q] * (
2651  M_density * M_alpha / M_timestep * dphi_phys_pressure[i_test][q][0] * M_phi_velocity[i_trial][q]
2652  +M_density * dphi_phys_pressure[i_test][q][0] * (
2653  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2654  )
2655  -M_viscosity * dphi_phys_pressure[i_test][q][0] * (
2656  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2657  )
2658  ) * w_quad[q];
2659 
2660  i_01 += M_Tau_M_hat[i_elem][q] * (
2661  M_density * M_alpha / M_timestep * dphi_phys_pressure[i_test][q][1] * M_phi_velocity[i_trial][q]
2662  +M_density * dphi_phys_pressure[i_test][q][1] * (
2663  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2664  )
2665  -M_viscosity * dphi_phys_pressure[i_test][q][1] * (
2666  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2667  )
2668  ) * w_quad[q];
2669 
2670  i_02 += M_Tau_M_hat[i_elem][q] * (
2671  M_density * M_alpha / M_timestep * dphi_phys_pressure[i_test][q][2] * M_phi_velocity[i_trial][q]
2672  +M_density * dphi_phys_pressure[i_test][q][2] * (
2673  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2674  )
2675  -M_viscosity * dphi_phys_pressure[i_test][q][2] * (
2676  d2phi_phys_velocity[i_trial][q][0][0] + d2phi_phys_velocity[i_trial][q][1][1] + d2phi_phys_velocity[i_trial][q][2][2]
2677  )
2678  ) * w_quad[q];
2679  }
2680 
2681  M_vals_supg_10[i_elem][0][i_test][i_trial] = i_00 * M_detJacobian[i_elem];
2682  M_vals_supg_10[i_elem][1][i_test][i_trial] = i_01 * M_detJacobian[i_elem];
2683  M_vals_supg_10[i_elem][2][i_test][i_trial] = i_02 * M_detJacobian[i_elem];
2684  }
2685  }
2686 
2687  }
2688  }
2689 
2690  for ( UInt d1 = 0; d1 < 3 ; d1++ ) // row index
2691  {
2692  for ( int k = 0; k < M_numElements; ++k )
2693  {
2694  for ( UInt i = 0; i < ndof_velocity; i++ )
2695  {
2696  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
2697  }
2698  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2699  ndof_velocity, M_cols_velocity[k],
2700  M_vals_supg[k][d1][0],
2701  Epetra_FECrsMatrix::ROW_MAJOR);
2702 
2703  block01->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2704  ndof_pressure, M_cols_pressure[k],
2705  M_vals_supg_01[k][d1],
2706  Epetra_FECrsMatrix::ROW_MAJOR);
2707  }
2708 
2709  for ( int k = 0; k < M_numElements; ++k )
2710  {
2711  for ( UInt i = 0; i < ndof_velocity; i++ )
2712  {
2713  M_cols_tmp[k][i] = M_cols_velocity[k][i] + M_numScalarDofs;
2714  }
2715  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2716  ndof_velocity, M_cols_tmp[k],
2717  M_vals_supg[k][d1][1],
2718  Epetra_FECrsMatrix::ROW_MAJOR);
2719  }
2720 
2721  for ( int k = 0; k < M_numElements; ++k )
2722  {
2723  for ( UInt i = 0; i < ndof_velocity; i++ )
2724  {
2725  M_cols_tmp[k][i] = M_cols_velocity[k][i] + 2 * M_numScalarDofs;
2726  }
2727  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
2728  ndof_velocity, M_cols_tmp[k],
2729  M_vals_supg[k][d1][2],
2730  Epetra_FECrsMatrix::ROW_MAJOR);
2731  }
2732  }
2733 
2734  for ( int k = 0; k < M_numElements; ++k )
2735  {
2736  block11->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k], ndof_pressure, M_cols_pressure[k], M_vals_11[k], Epetra_FECrsMatrix::ROW_MAJOR);
2737  }
2738 
2739  for ( UInt d1 = 0; d1 < 3 ; d1++ )
2740  {
2741  for ( int k = 0; k < M_numElements; ++k )
2742  {
2743  for ( UInt i = 0; i < ndof_velocity; i++ )
2744  {
2745  M_cols_tmp[k][i] = M_cols_velocity[k][i] + d1 * M_numScalarDofs;
2746  }
2747  block10->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k], ndof_velocity, M_cols_tmp[k], M_vals_supg_10[k][d1], Epetra_FECrsMatrix::ROW_MAJOR);
2748  }
2749  }
2750 }
2751 
2752 void
2753 FastAssemblerNS::vmsles_semi_implicit_termsP1P1 ( matrixPtr_Type& block00,
2754  matrixPtr_Type& block01,
2755  matrixPtr_Type& block10,
2756  matrixPtr_Type& block11,
2757  const std::vector<std::vector<VectorSmall<3>>>& fine_scale,
2758  const vector_Type& u_extr )
2759 {
2760  int ndof_velocity = M_referenceFE_velocity->nbDof();
2761  int ndof_pressure = M_referenceFE_pressure->nbDof();
2762  int NumQuadPoints = M_qr->nbQuadPt();
2763  int ndof_vec = M_referenceFE_velocity->nbDof()*3;
2764  M_numScalarDofs = M_fespace_velocity->dof().numTotalDof();
2765 
2766  double w_quad[NumQuadPoints];
2767  for ( int q = 0; q < NumQuadPoints ; q++ )
2768  {
2769  w_quad[q] = M_qr->weight(q);
2770  }
2771 
2772  #pragma omp parallel firstprivate( w_quad, ndof_velocity, ndof_pressure, NumQuadPoints)
2773  {
2774  int i_elem, i_dof, q, d1, d2, i_test, i_trial, e_idof, dim_mat, iCoor, jCoor, k1, k2;
2775  double integral, integral_test, integral_trial, integral_partial, integral_lapl, partialSum;
2776 
2777  double dphi_phys_velocity[ndof_velocity][NumQuadPoints][3];
2778  double dphi_phys_pressure[ndof_pressure][NumQuadPoints][3];
2779  double uhq[3][NumQuadPoints];
2780  double i_00, i_01, i_02, i_10, i_11, i_12, i_20, i_21, i_22;
2781 
2782 
2783  // ELEMENTI,
2784  #pragma omp for
2785  for ( i_elem = 0; i_elem < M_numElements; i_elem++ )
2786  {
2787  // DOF
2788  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
2789  {
2790  // QUAD
2791  for ( q = 0; q < NumQuadPoints ; q++ )
2792  {
2793  // DIM 1
2794  for ( d1 = 0; d1 < 3 ; d1++ )
2795  {
2796  dphi_phys_velocity[i_dof][q][d1] = 0.0;
2797 
2798  // DIM 2
2799  for ( d2 = 0; d2 < 3 ; d2++ )
2800  {
2801  dphi_phys_velocity[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_velocity[i_dof][q][d2];
2802  }
2803  }
2804  }
2805  }
2806 
2807  for ( i_dof = 0; i_dof < ndof_pressure; i_dof++ )
2808  {
2809  // QUAD
2810  for ( q = 0; q < NumQuadPoints ; q++ )
2811  {
2812  // DIM 1
2813  for ( d1 = 0; d1 < 3 ; d1++ )
2814  {
2815  dphi_phys_pressure[i_dof][q][d1] = 0.0;
2816 
2817  // DIM 2
2818  for ( d2 = 0; d2 < 3 ; d2++ )
2819  {
2820  dphi_phys_pressure[i_dof][q][d1] += M_invJacobian[i_elem][d1][d2] * M_dphi_pressure[i_dof][q][d2];
2821  }
2822  }
2823  }
2824  }
2825 
2826  // QUAD
2827  for ( q = 0; q < NumQuadPoints ; q++ )
2828  {
2829  for ( d1 = 0; d1 < 3 ; d1++ )
2830  {
2831  uhq[d1][q] = 0.0;
2832  for ( i_dof = 0; i_dof < ndof_velocity; i_dof++ )
2833  {
2834  e_idof = M_elements_velocity[i_elem][i_dof] + d1*M_numScalarDofs ;
2835  uhq[d1][q] += u_extr[e_idof] * M_phi_velocity[i_dof][q];
2836  }
2837  }
2838  }
2839 
2840  // STABILIZZAZIONE - coefficienti Tau_M e Tau_C
2841  for ( q = 0; q < NumQuadPoints ; q++ )
2842  {
2843  M_Tau_M[i_elem][q] = 1.0/std::sqrt(
2844  M_density*M_density*( uhq[0][q]*( M_G[i_elem][0][0]* uhq[0][q] + M_G[i_elem][0][1] * uhq[1][q] + M_G[i_elem][0][2] * uhq[2][q] ) +
2845  uhq[1][q]*( M_G[i_elem][1][0]* uhq[0][q] + M_G[i_elem][1][1] * uhq[1][q] + M_G[i_elem][1][2] * uhq[2][q] ) +
2846  uhq[2][q]*( M_G[i_elem][2][0]* uhq[0][q] + M_G[i_elem][2][1] * uhq[1][q] + M_G[i_elem][2][2] * uhq[2][q] )
2847  ) // TAU_M_DEN_VEL
2848  +M_C_I*M_viscosity*M_viscosity*(
2849  M_G[i_elem][0][0]*M_G[i_elem][0][0] + M_G[i_elem][0][1]*M_G[i_elem][0][1] + M_G[i_elem][0][2]*M_G[i_elem][0][2] +
2850  M_G[i_elem][1][0]*M_G[i_elem][1][0] + M_G[i_elem][1][1]*M_G[i_elem][1][1] + M_G[i_elem][1][2]*M_G[i_elem][1][2] +
2851  M_G[i_elem][2][0]*M_G[i_elem][2][0] + M_G[i_elem][2][1]*M_G[i_elem][2][1] + M_G[i_elem][2][2]*M_G[i_elem][2][2]
2852  ) // TAU_M_DEN_VISC*/
2853  );
2854 
2855  M_Tau_C[i_elem][q] = 1.0/( M_g[i_elem][0]*M_Tau_M[i_elem][q]*M_g[i_elem][0] +
2856  M_g[i_elem][1]*M_Tau_M[i_elem][q]*M_g[i_elem][1] +
2857  M_g[i_elem][2]*M_Tau_M[i_elem][q]*M_g[i_elem][2]
2858  );
2859 
2860  M_Tau_M_hat[i_elem][q] = 1.0/( M_alpha * M_density / M_timestep + 1.0/M_Tau_M[i_elem][q]);
2861  }
2862 
2863  // DOF - test
2864  for ( i_test = 0; i_test < ndof_velocity; i_test++ )
2865  {
2866  M_rows_velocity[i_elem][i_test] = M_elements_velocity[i_elem][i_test];
2867 
2868  // DOF - trial
2869  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
2870  {
2871  M_cols_velocity[i_elem][i_trial] = M_elements_velocity[i_elem][i_trial];
2872 
2873  i_00 = 0.0; i_01 = 0.0; i_02 = 0.0; i_10 = 0.0; i_11 = 0.0; i_12 = 0.0; i_20 = 0.0; i_21 = 0.0; i_22 = 0.0;
2874 
2875  M_vals_supg[i_elem][0][0][i_test][i_trial] = 0.0;
2876  M_vals_supg[i_elem][0][1][i_test][i_trial] = 0.0;
2877  M_vals_supg[i_elem][0][2][i_test][i_trial] = 0.0;
2878  M_vals_supg[i_elem][1][0][i_test][i_trial] = 0.0;
2879  M_vals_supg[i_elem][1][1][i_test][i_trial] = 0.0;
2880  M_vals_supg[i_elem][1][2][i_test][i_trial] = 0.0;
2881  M_vals_supg[i_elem][2][0][i_test][i_trial] = 0.0;
2882  M_vals_supg[i_elem][2][1][i_test][i_trial] = 0.0;
2883  M_vals_supg[i_elem][2][2][i_test][i_trial] = 0.0;
2884 
2885  double tmp_supg = 0.0;
2886  double tmp_supg_test = 0.0;
2887  double tmp_supg_trial = 0.0;
2888 
2889  // QUAD
2890  for ( q = 0; q < NumQuadPoints ; q++ )
2891  {
2892  tmp_supg = M_alpha * M_density * M_density / M_timestep * (
2893  dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q]
2894  ) * M_phi_velocity[i_trial][q];
2895 
2896  tmp_supg_test =
2897  dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q];
2898 
2899  tmp_supg_trial =
2900  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q];
2901 
2902  i_00 += ( M_Tau_M_hat[i_elem][q] * (
2903  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] * M_phi_velocity[i_trial][q] +
2904  M_density * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] * (
2905  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2906  )
2907  // VMS
2908  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][0] * uhq[0][q] * M_phi_velocity[i_trial][q] )
2909  +M_density * M_density * dphi_phys_velocity[i_test][q][0] * uhq[0][q] * (
2910  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2911  )
2912  // SUPG
2913  +tmp_supg
2914  +M_density * M_density * tmp_supg_test*tmp_supg_trial
2915  )
2916  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][0]
2917 
2918  ) * w_quad[q];
2919 
2920  i_01 += ( M_Tau_M_hat[i_elem][q] * (
2921  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] * M_phi_velocity[i_trial][q] +
2922  M_density * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] * (
2923  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2924  )
2925  // VMS
2926  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][1] * uhq[0][q] * M_phi_velocity[i_trial][q] )
2927  +M_density * M_density * dphi_phys_velocity[i_test][q][1] * uhq[0][q] * (
2928  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2929  )
2930 
2931  )
2932  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][1]
2933  ) * w_quad[q];
2934 
2935  i_02 += ( M_Tau_M_hat[i_elem][q] * (
2936  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] * M_phi_velocity[i_trial][q] +
2937  M_density * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] * (
2938  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2939  )
2940 
2941  // VMS
2942  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][2] * uhq[0][q] * M_phi_velocity[i_trial][q] )
2943  +M_density * M_density * dphi_phys_velocity[i_test][q][2] * uhq[0][q] * (
2944  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2945  )
2946  )
2947  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][0] * dphi_phys_velocity[i_trial][q][2]
2948  ) * w_quad[q];
2949 
2950  i_10 += ( M_Tau_M_hat[i_elem][q] * (
2951  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] * M_phi_velocity[i_trial][q] +
2952  M_density * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] * (
2953  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2954  )
2955 
2956  // VMS
2957  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][0] * uhq[1][q] * M_phi_velocity[i_trial][q] )
2958  +M_density * M_density * dphi_phys_velocity[i_test][q][0] * uhq[1][q] * (
2959  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2960  )
2961  )
2962  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][0]
2963  ) * w_quad[q];
2964 
2965  i_11 += ( M_Tau_M_hat[i_elem][q] * (
2966  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] * M_phi_velocity[i_trial][q] +
2967  M_density * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] * (
2968  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2969  )
2970 
2971  // VMS
2972  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][1] * uhq[1][q] * M_phi_velocity[i_trial][q] )
2973  +M_density * M_density * dphi_phys_velocity[i_test][q][1] * uhq[1][q] * (
2974  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2975  )
2976  // SUPG
2977  +tmp_supg
2978  +M_density * M_density * tmp_supg_test*tmp_supg_trial
2979  )
2980  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][1]
2981  ) * w_quad[q];
2982 
2983  i_12 += ( M_Tau_M_hat[i_elem][q] * (
2984  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] * M_phi_velocity[i_trial][q] +
2985  M_density * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] * (
2986  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2987  )
2988  // VMS
2989  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][2] * uhq[1][q] * M_phi_velocity[i_trial][q] )
2990  +M_density * M_density * dphi_phys_velocity[i_test][q][2] * uhq[1][q] * (
2991  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
2992  )
2993  )
2994  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][1] * dphi_phys_velocity[i_trial][q][2]
2995  ) * w_quad[q];
2996 
2997  i_20 += ( M_Tau_M_hat[i_elem][q] * (
2998  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] * M_phi_velocity[i_trial][q] +
2999  M_density * dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] * (
3000  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3001  )
3002  // VMS
3003  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][0] * uhq[2][q] * M_phi_velocity[i_trial][q] )
3004  +M_density * M_density * dphi_phys_velocity[i_test][q][0] * uhq[2][q] * (
3005  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3006  )
3007  )
3008  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][0]
3009  ) * w_quad[q];
3010 
3011  i_21 += ( M_Tau_M_hat[i_elem][q] * (
3012  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] * M_phi_velocity[i_trial][q] +
3013  M_density * dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] * (
3014  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3015  )
3016  // VMS
3017  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][1] * uhq[2][q] * M_phi_velocity[i_trial][q] )
3018  +M_density * M_density * dphi_phys_velocity[i_test][q][1] * uhq[2][q] * (
3019  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3020  )
3021  )
3022  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][1]
3023  ) * w_quad[q];
3024 
3025  i_22 += ( M_Tau_M_hat[i_elem][q] * (
3026  M_alpha * M_density * M_density / M_timestep * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] * M_phi_velocity[i_trial][q] +
3027  M_density * dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] * (
3028  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3029  )
3030  // VMS
3031  +M_alpha * M_density * M_density / M_timestep * ( dphi_phys_velocity[i_test][q][2] * uhq[2][q] * M_phi_velocity[i_trial][q] )
3032  +M_density * M_density * dphi_phys_velocity[i_test][q][2] * uhq[2][q] * (
3033  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3034  )
3035  // SUPG
3036  +tmp_supg
3037  +M_density * M_density * tmp_supg_test*tmp_supg_trial
3038  )
3039  + M_Tau_C[i_elem][q] * dphi_phys_velocity[i_test][q][2] * dphi_phys_velocity[i_trial][q][2]
3040  ) * w_quad[q];
3041 
3042 
3043  } // end quadrature
3044 
3045  M_vals_supg[i_elem][0][0][i_test][i_trial] = (i_00 ) * M_detJacobian[i_elem];
3046  M_vals_supg[i_elem][0][1][i_test][i_trial] = (i_01 ) * M_detJacobian[i_elem];
3047  M_vals_supg[i_elem][0][2][i_test][i_trial] = (i_02 ) * M_detJacobian[i_elem];
3048  M_vals_supg[i_elem][1][0][i_test][i_trial] = (i_10 ) * M_detJacobian[i_elem];
3049  M_vals_supg[i_elem][1][1][i_test][i_trial] = (i_11 ) * M_detJacobian[i_elem];
3050  M_vals_supg[i_elem][1][2][i_test][i_trial] = (i_12 ) * M_detJacobian[i_elem];
3051  M_vals_supg[i_elem][2][0][i_test][i_trial] = (i_20 ) * M_detJacobian[i_elem];
3052  M_vals_supg[i_elem][2][1][i_test][i_trial] = (i_21 ) * M_detJacobian[i_elem];
3053  M_vals_supg[i_elem][2][2][i_test][i_trial] = (i_22 ) * M_detJacobian[i_elem];
3054 
3055  } // end trial
3056 
3057  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
3058  {
3059  M_cols_pressure[i_elem][i_trial] = M_elements_pressure[i_elem][i_trial];
3060  M_vals_supg_01[i_elem][0][i_test][i_trial] = 0.0;
3061  M_vals_supg_01[i_elem][1][i_test][i_trial] = 0.0;
3062  M_vals_supg_01[i_elem][2][i_test][i_trial] = 0.0;
3063 
3064  i_00 = 0.0;
3065  i_10 = 0.0;
3066  i_20 = 0.0;
3067 
3068  // QUAD
3069  for ( q = 0; q < NumQuadPoints ; q++ )
3070  {
3071  i_00 += M_Tau_M_hat[i_elem][q] * (
3072  dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][0] * dphi_phys_pressure[i_trial][q][0]
3073  +dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][0] * dphi_phys_pressure[i_trial][q][1]
3074  +dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][0] * dphi_phys_pressure[i_trial][q][2]
3075  +M_density * (
3076  (dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q])
3077  * dphi_phys_pressure[i_trial][q][0]
3078  + dphi_phys_velocity[i_test][q][0] * uhq[0][q] * dphi_phys_pressure[i_trial][q][0]
3079  + dphi_phys_velocity[i_test][q][1] * uhq[0][q] * dphi_phys_pressure[i_trial][q][1]
3080  + dphi_phys_velocity[i_test][q][2] * uhq[0][q] * dphi_phys_pressure[i_trial][q][2]
3081  )
3082  ) * w_quad[q];
3083 
3084  i_10 += M_Tau_M_hat[i_elem][q] * (
3085  dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][1] * dphi_phys_pressure[i_trial][q][0]
3086  +dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][1] * dphi_phys_pressure[i_trial][q][1]
3087  +dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][1] * dphi_phys_pressure[i_trial][q][2]
3088  +M_density * (
3089  (dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q])
3090  * dphi_phys_pressure[i_trial][q][1]
3091  + dphi_phys_velocity[i_test][q][0] * uhq[1][q] * dphi_phys_pressure[i_trial][q][0]
3092  + dphi_phys_velocity[i_test][q][1] * uhq[1][q] * dphi_phys_pressure[i_trial][q][1]
3093  + dphi_phys_velocity[i_test][q][2] * uhq[1][q] * dphi_phys_pressure[i_trial][q][2]
3094  )
3095  ) * w_quad[q];
3096 
3097  i_20 += M_Tau_M_hat[i_elem][q] * (
3098  dphi_phys_velocity[i_test][q][0] * fine_scale[i_elem][q][2] * dphi_phys_pressure[i_trial][q][0]
3099  +dphi_phys_velocity[i_test][q][1] * fine_scale[i_elem][q][2] * dphi_phys_pressure[i_trial][q][1]
3100  +dphi_phys_velocity[i_test][q][2] * fine_scale[i_elem][q][2] * dphi_phys_pressure[i_trial][q][2]
3101  +M_density * (
3102  (dphi_phys_velocity[i_test][q][0] * uhq[0][q] + dphi_phys_velocity[i_test][q][1] * uhq[1][q] + dphi_phys_velocity[i_test][q][2] * uhq[2][q])
3103  * dphi_phys_pressure[i_trial][q][2]
3104  + dphi_phys_velocity[i_test][q][0] * uhq[2][q] * dphi_phys_pressure[i_trial][q][0]
3105  + dphi_phys_velocity[i_test][q][1] * uhq[2][q] * dphi_phys_pressure[i_trial][q][1]
3106  + dphi_phys_velocity[i_test][q][2] * uhq[2][q] * dphi_phys_pressure[i_trial][q][2]
3107  )
3108  ) * w_quad[q];
3109  }
3110 
3111  M_vals_supg_01[i_elem][0][i_test][i_trial] = i_00 * M_detJacobian[i_elem];
3112  M_vals_supg_01[i_elem][1][i_test][i_trial] = i_10 * M_detJacobian[i_elem];
3113  M_vals_supg_01[i_elem][2][i_test][i_trial] = i_20 * M_detJacobian[i_elem];
3114 
3115  }
3116 
3117  } // end dof test
3118 
3119  // DOF - test - assemblo blocchi (1,0) e (1,1)
3120  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
3121  {
3122  M_rows_pressure[i_elem][i_test] = M_elements_pressure[i_elem][i_test];
3123 
3124  // DOF - trial
3125  for ( i_trial = 0; i_trial < ndof_pressure; i_trial++ )
3126  {
3127  M_vals_11[i_elem][i_test][i_trial] = 0.0;
3128 
3129  integral = 0.0;
3130  // QUAD
3131  for ( q = 0; q < NumQuadPoints ; q++ )
3132  {
3133  // DIM 1
3134  for ( d1 = 0; d1 < 3 ; d1++ )
3135  {
3136  integral += M_Tau_M_hat[i_elem][q] * dphi_phys_pressure[i_test][q][d1] * dphi_phys_pressure[i_trial][q][d1]*w_quad[q];
3137  }
3138  }
3139  M_vals_11[i_elem][i_test][i_trial] = integral * M_detJacobian[i_elem];
3140  }
3141  }
3142 
3143  // DOF - test
3144  for ( i_test = 0; i_test < ndof_pressure; i_test++ )
3145  {
3146  // DOF - trial
3147  for ( i_trial = 0; i_trial < ndof_velocity; i_trial++ )
3148  {
3149  i_00 = 0.0;
3150  i_01 = 0.0;
3151  i_02 = 0.0;
3152 
3153  // QUAD
3154  for ( q = 0; q < NumQuadPoints ; q++ )
3155  {
3156  i_00 += M_Tau_M_hat[i_elem][q] * (
3157  M_density * M_alpha / M_timestep * dphi_phys_pressure[i_test][q][0] * M_phi_velocity[i_trial][q]
3158  +M_density * dphi_phys_pressure[i_test][q][0] * (
3159  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3160  )
3161 
3162  ) * w_quad[q];
3163 
3164  i_01 += M_Tau_M_hat[i_elem][q] * (
3165  M_density * M_alpha / M_timestep * dphi_phys_pressure[i_test][q][1] * M_phi_velocity[i_trial][q]
3166  +M_density * dphi_phys_pressure[i_test][q][1] * (
3167  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3168  )
3169 
3170  ) * w_quad[q];
3171 
3172  i_02 += M_Tau_M_hat[i_elem][q] * (
3173  M_density * M_alpha / M_timestep * dphi_phys_pressure[i_test][q][2] * M_phi_velocity[i_trial][q]
3174  +M_density * dphi_phys_pressure[i_test][q][2] * (
3175  dphi_phys_velocity[i_trial][q][0] * uhq[0][q] + dphi_phys_velocity[i_trial][q][1] * uhq[1][q] + dphi_phys_velocity[i_trial][q][2] * uhq[2][q]
3176  )
3177  ) * w_quad[q];
3178  }
3179 
3180  M_vals_supg_10[i_elem][0][i_test][i_trial] = i_00 * M_detJacobian[i_elem];
3181  M_vals_supg_10[i_elem][1][i_test][i_trial] = i_01 * M_detJacobian[i_elem];
3182  M_vals_supg_10[i_elem][2][i_test][i_trial] = i_02 * M_detJacobian[i_elem];
3183  }
3184  }
3185 
3186  }
3187  }
3188 
3189  for ( UInt d1 = 0; d1 < 3 ; d1++ ) // row index
3190  {
3191  for ( int k = 0; k < M_numElements; ++k )
3192  {
3193  for ( UInt i = 0; i < ndof_velocity; i++ )
3194  {
3195  M_rows_tmp[k][i] = M_rows_velocity[k][i] + d1 * M_numScalarDofs;
3196  }
3197  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
3198  ndof_velocity, M_cols_velocity[k],
3199  M_vals_supg[k][d1][0],
3200  Epetra_FECrsMatrix::ROW_MAJOR);
3201 
3202  block01->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
3203  ndof_pressure, M_cols_pressure[k],
3204  M_vals_supg_01[k][d1],
3205  Epetra_FECrsMatrix::ROW_MAJOR);
3206  }
3207 
3208  for ( int k = 0; k < M_numElements; ++k )
3209  {
3210  for ( UInt i = 0; i < ndof_velocity; i++ )
3211  {
3212  M_cols_tmp[k][i] = M_cols_velocity[k][i] + M_numScalarDofs;
3213  }
3214  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
3215  ndof_velocity, M_cols_tmp[k],
3216  M_vals_supg[k][d1][1],
3217  Epetra_FECrsMatrix::ROW_MAJOR);
3218  }
3219 
3220  for ( int k = 0; k < M_numElements; ++k )
3221  {
3222  for ( UInt i = 0; i < ndof_velocity; i++ )
3223  {
3224  M_cols_tmp[k][i] = M_cols_velocity[k][i] + 2 * M_numScalarDofs;
3225  }
3226  block00->matrixPtr()->InsertGlobalValues ( ndof_velocity, M_rows_tmp[k],
3227  ndof_velocity, M_cols_tmp[k],
3228  M_vals_supg[k][d1][2],
3229  Epetra_FECrsMatrix::ROW_MAJOR);
3230  }
3231  }
3232 
3233  for ( int k = 0; k < M_numElements; ++k )
3234  {
3235  block11->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k], ndof_pressure, M_cols_pressure[k], M_vals_11[k], Epetra_FECrsMatrix::ROW_MAJOR);
3236  }
3237 
3238  for ( UInt d1 = 0; d1 < 3 ; d1++ )
3239  {
3240  for ( int k = 0; k < M_numElements; ++k )
3241  {
3242  for ( UInt i = 0; i < ndof_velocity; i++ )
3243  {
3244  M_cols_tmp[k][i] = M_cols_velocity[k][i] + d1 * M_numScalarDofs;
3245  }
3246  block10->matrixPtr()->InsertGlobalValues ( ndof_pressure, M_rows_pressure[k], ndof_velocity, M_cols_tmp[k], M_vals_supg_10[k][d1], Epetra_FECrsMatrix::ROW_MAJOR);
3247  }
3248  }
3249 }
LifeV::QuadratureRule::weight
const Real & weight(const UInt &ig) const
weight(ig) is the ig-th quadrature weight
Definition: QuadratureRule.hpp:284
LifeV::VectorEpetra::operator[]
const data_type & operator[](const UInt row) const
Access operators.
Definition: VectorEpetra.cpp:184
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::CurrentFE::tInverseJacobian
Real tInverseJacobian(UInt element1, UInt element2, UInt quadNode) const
Getter for the inverse of the jacobian.
Definition: CurrentFE.hpp:465
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::CurrentFE
CurrentFE - A primordial class for the assembly of the local matrices/vectors retaining the values on...
Definition: CurrentFE.hpp:243
LifeV::ReferenceFE
The class for a reference Lagrangian finite element.
Definition: ReferenceFE.hpp:115
LifeV::CurrentFE::detJacobian
Real detJacobian(UInt quadNode) const
Getter for the determinant of the jacobian in a given quadrature node.
Definition: CurrentFE.hpp:451
LifeV::QuadratureRule::nbQuadPt
const UInt & nbQuadPt() const
Getter for the number of quadrature points.
Definition: QuadratureRule.hpp:305
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191