PageRenderTime 1670ms CodeModel.GetById 222ms app.highlight 1136ms RepoModel.GetById 133ms app.codeStats 1ms

/navierStokes3D.cpp

https://bitbucket.org/cuneytsert/cfd-with-cuda
C++ | 1837 lines | 1330 code | 343 blank | 164 comment | 279 complexity | d2fd647d1ad19eccc9b5834cd0590691 MD5 | raw file

Large files files are truncated, but you can click here to view the full file

   1/*****************************************************************
   2*        This code is a part of the CFD-with-CUDA project        *
   3*             http://code.google.com/p/cfd-with-cuda             *
   4*                                                                *
   5*              Dr. Cuneyt Sert and Mahmut M. Gocmen              *
   6*                                                                *
   7*              Department of Mechanical Engineering              *
   8*                Middle East Technical University                *
   9*                         Ankara, Turkey                         *
  10*            http://www.me.metu.edu.tr/people/cuneyt             *
  11*****************************************************************/
  12
  13#include <stdio.h>

  14#include <string>

  15#include <iostream>

  16#include <iomanip>

  17#include <fstream>

  18#include <ctime>

  19#include <cmath>

  20#include <algorithm>

  21
  22using namespace std;
  23
  24#ifdef SINGLE

  25  typedef float real;
  26#else

  27  typedef double real;
  28#endif

  29
  30ifstream meshfile;
  31ifstream problemFile;
  32ofstream outputFile;
  33ofstream outputControl;
  34
  35string problemNameFile, whichProblem;
  36string problemName = "ProblemName.txt";
  37string controlFile = "Control_Output.txt";
  38string inputExtension  = ".inp";              // Change this line to specify the extension of input file (with dot).
  39string outputExtension  = ".dat";             // Change this line to specify the extension of output file (with dot).
  40
  41int name, eType, NE, NN, NGP, NEU, Ndof;
  42int NCN, NENv, NENp, nonlinearIterMax, solverIterMax;
  43double density, viscosity, fx, fy, nonlinearTol, solverTol;
  44int **LtoG, **velNodes, **pressureNodes;
  45double **coord;
  46int nBC, nVelNodes, nPressureNodes;
  47double *BCtype, **BCstrings;
  48double axyFunc, fxyFunc;
  49double **GQpoint, *GQweight;
  50double **Sp, ***DSp, **Sv, ***DSv;
  51double **detJacob, ****gDSp, ****gDSv;
  52real **K, *F, *u, *uOld;   // Can be float or double. K is the stiffness matrix
  53                           // in full storage used Gauss Elimination solver.
  54int bigNumber;
  55
  56double *Fe, **Ke;
  57double *Fe_1, *Fe_2, *Fe_3, *Fe_4;
  58double **Ke_11, **Ke_12, **Ke_13, **Ke_14, **Ke_22, **Ke_23, **Ke_24, **Ke_33, **Ke_34, **Ke_41, **Ke_42, **Ke_43, **Ke_44;
  59double **Ke_11_add, **Ke_12_add, **Ke_13_add, **Ke_14_add, **Ke_22_add, **Ke_23_add, **Ke_24_add;
  60double **Ke_33_add, **Ke_34_add, **Ke_41_add, **Ke_42_add, **Ke_43_add, **Ke_44_add;
  61double *uNodal, *vNodal, *wNodal;
  62double u0, v0, w0;
  63
  64int **GtoL, *rowStarts, *rowStartsSmall, *colSmall, *col, **KeKMapSmall, NNZ; 
  65real *val;
  66
  67void readInput();
  68void gaussQuad();
  69void calcShape();
  70void calcJacobian();
  71void initGlobalSysVariables();
  72void calcGlobalSys();
  73void assemble(int e, double **Ke, double *Fe);
  74void applyBC();
  75void solve();
  76void postProcess();
  77void gaussElimination(int N, real **K, real *F, real *u, bool& err);
  78void writeTecplotFile();
  79void compressedSparseRowStorage();
  80#ifdef CUSP

  81   extern void CUSPsolver();
  82#endif

  83
  84
  85
  86//------------------------------------------------------------------------------
  87int main()
  88//------------------------------------------------------------------------------
  89{
  90   cout << "\n********************************************************";
  91   cout << "\n*   Navier-Stokes 3D - Part of CFD with CUDA project   *";
  92   cout << "\n*       http://code.google.com/p/cfd-with-cuda         *";
  93   cout << "\n********************************************************\n\n";
  94
  95   cout << "The program is started." << endl ;
  96
  97   time_t start, end;
  98   time (&start);     // Start measuring execution time.
  99
 100   readInput();                   cout << "Input file is read." << endl ;
 101   compressedSparseRowStorage();  cout << "CSR vectors are created." << endl ;
 102   gaussQuad();
 103   calcShape();
 104   calcJacobian();
 105   initGlobalSysVariables();
 106   solve();
 107   //postProcess();
 108   writeTecplotFile();            cout << "A DAT file is created for Tecplot." << endl ;
 109
 110   time (&end);      // Stop measuring execution time.
 111   cout << endl << "Elapsed wall clock time is " << difftime (end,start) << " seconds." << endl;
 112   
 113   cout << endl << "The program is terminated successfully.\nPress a key to close this window...";
 114
 115   cin.get();
 116   return 0;
 117
 118} // End of function main()
 119
 120
 121
 122
 123//------------------------------------------------------------------------------
 124void readInput()
 125//------------------------------------------------------------------------------
 126{
 127
 128   string dummy, dummy2, dummy4, dummy5;
 129   int dummy3, i;
 130
 131   problemFile.open(problemName.c_str(), ios::in);   // Read problem name
 132   problemFile >> whichProblem;                                  
 133   problemFile.close();
 134   
 135   meshfile.open((whichProblem + inputExtension).c_str(), ios::in);
 136     
 137   meshfile.ignore(256, '\n'); // Read and ignore the line
 138   meshfile.ignore(256, '\n'); // Read and ignore the line
 139   meshfile >> dummy >> dummy2 >> eType;
 140   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 141   meshfile >> dummy >> dummy2 >> NE;
 142   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 143   meshfile >> dummy >> dummy2 >> NCN;
 144   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 145   meshfile >> dummy >> dummy2 >> NN;
 146   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 147   meshfile >> dummy >> dummy2 >> NENv;
 148   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 149   meshfile >> dummy >> dummy2 >> NENp;   
 150   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 151   meshfile >> dummy >> dummy2 >> NGP;
 152   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 153   meshfile >> dummy >> dummy2 >> nonlinearIterMax;
 154   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 155   meshfile >> dummy >> dummy2 >> nonlinearTol;
 156   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 157   meshfile >> dummy >> dummy2 >> solverIterMax;
 158   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 159   meshfile >> dummy >> dummy2 >> solverTol;
 160   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 161   meshfile >> dummy >> dummy2 >> density;
 162   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 163   meshfile >> dummy >> dummy2 >> viscosity;
 164   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 165   meshfile >> dummy >> dummy2 >> fx;
 166   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 167   meshfile >> dummy >> dummy2 >> fy;
 168   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 169   meshfile.ignore(256, '\n'); // Read and ignore the line
 170   meshfile.ignore(256, '\n'); // Read and ignore the line  
 171   
 172   NEU = 3*NENv + NENp;
 173   Ndof = 3*NN + NCN;
 174   
 175   // Read node coordinates
 176   coord = new double*[NN];
 177   
 178   if (eType == 2 || eType == 1) {
 179      for (i = 0; i < NN; i++) {
 180         coord[i] = new double[2];
 181      }
 182   	
 183      for (i=0; i<NN; i++){
 184         meshfile >> dummy3 >> coord[i][0] >> coord[i][1] ;
 185         meshfile.ignore(256, '\n'); // Ignore the rest of the line    
 186      }  
 187   } else {
 188      for (i = 0; i < NN; i++) {
 189         coord[i] = new double[3];
 190      }
 191   	
 192      for (i=0; i<NN; i++){
 193         meshfile >> dummy3 >> coord[i][0] >> coord[i][1] >> coord[i][2] ;
 194         meshfile.ignore(256, '\n'); // Ignore the rest of the line    
 195      }  
 196   }
 197   
 198   meshfile.ignore(256, '\n'); // Read and ignore the line
 199   meshfile.ignore(256, '\n'); // Read and ignore the line 
 200    
 201   // Read element connectivity, i.e. LtoG
 202   LtoG = new int*[NE];
 203
 204   for (i=0; i<NE; i++) {
 205      LtoG[i] = new int[NENv];
 206   }
 207
 208   for (int e = 0; e<NE; e++){
 209      meshfile >> dummy3;
 210      for (i = 0; i<NENv; i++){
 211         meshfile >> LtoG[e][i];
 212      }
 213      meshfile.ignore(256, '\n'); // Read and ignore the line 
 214   } 
 215
 216   meshfile.ignore(256, '\n'); // Read and ignore the line
 217   meshfile.ignore(256, '\n'); // Read and ignore the line 
 218
 219   // Read number of different BC types and details of each BC
 220   meshfile >> dummy >> dummy2 >> nBC;
 221   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 222      
 223   BCtype = new double[nBC];
 224   BCstrings = new double*[nBC];
 225   for (i=0; i<nBC; i++) {
 226      BCstrings[i] = new double[3];
 227   }   
 228    
 229   for (i = 0; i<nBC-1; i++){
 230      meshfile >> dummy >> BCtype[i] >> dummy2 >> dummy3 >> BCstrings[i][0] >> dummy4 >> BCstrings[i][1] >> dummy5 >> BCstrings[i][2];
 231      meshfile.ignore(256, '\n'); // Ignore the rest of the line
 232   }
 233   meshfile >> dummy >> BCtype[2]  >> dummy >> BCstrings[i][0];
 234   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 235   // Read EBC data 
 236    
 237   meshfile.ignore(256, '\n'); // Read and ignore the line 
 238   meshfile >> dummy >> dummy2 >> nVelNodes;
 239   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 240   meshfile >> dummy >> dummy2 >> nPressureNodes;
 241   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 242   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 243   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 244      
 245   if (nVelNodes!=0){
 246      velNodes = new int*[nVelNodes];
 247      for (i = 0; i < nVelNodes; i++){
 248         velNodes[i] = new int[2];
 249      }
 250      for (i = 0; i < nVelNodes; i++){
 251         meshfile >> velNodes[i][0] >> velNodes[i][1];   
 252         meshfile.ignore(256, '\n'); // Ignore the rest of the line
 253      }
 254   }
 255
 256   meshfile.ignore(256, '\n'); // Ignore the rest of the line  
 257   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 258   meshfile.ignore(256, '\n'); // Ignore the rest of the line
 259   
 260   if (nVelNodes!=0){
 261      pressureNodes = new int*[nPressureNodes];
 262      for (i = 0; i < nPressureNodes; i++){
 263         pressureNodes[i] = new int[2];
 264      }
 265      for (i = 0; i < nPressureNodes; i++){
 266         meshfile >> pressureNodes[i][0] >> pressureNodes[i][1];   
 267         meshfile.ignore(256, '\n'); // Ignore the rest of the line
 268      }
 269   }
 270  
 271   meshfile.close();
 272
 273} // End of function readInput()
 274
 275
 276
 277
 278//------------------------------------------------------------------------------
 279void compressedSparseRowStorage()
 280//------------------------------------------------------------------------------
 281{
 282//GtoL creation
 283int i, j, k, m, x, y, valGtoL, check, temp, *checkCol, noOfColGtoL, *GtoLCounter;
 284
 285   GtoL = new int*[NN];          // stores which elements connected to the node
 286   noOfColGtoL = 8;              // for 3D meshes created by our mesh generator max 8 elements connect to one node,
 287   GtoLCounter = new int[NN];    // but for real 3D problems this must be a big enough value!
 288   
 289   for (i=0; i<NN; i++) {     
 290      GtoL[i] = new int[noOfColGtoL];   
 291   }
 292
 293   for(i=0; i<NN; i++) {
 294      for(j=0; j<noOfColGtoL; j++) {
 295		   GtoL[i][j] = -1;      // for the nodes that didn't connect 4 nodes
 296	   }
 297      GtoLCounter[i] = 0;
 298   }
 299   
 300   for (i=0; i<NE; i++) {
 301      for(j=0; j<NENv; j++) {
 302         GtoL[LtoG[i][j]][GtoLCounter[LtoG[i][j]]] = i;
 303         GtoLCounter[LtoG[i][j]] += 1;
 304      }
 305   }         
 306   
 307   delete [] GtoLCounter;
 308   // Su anda gereksiz 0'lar oluþturuluyor. Ýleride düzeltilebilir.
 309   // for(i=0; i<nVelNodes; i++) {   // extracting EBC values from GtoL with making them "-1"
 310      // for(j=0; j<noOfColGtoL; j++) {
 311         // GtoL[velNodes[i][0]][j] = -1;
 312      // }   
 313   // } 
 314
 315
 316// Finding size of col vector, creation of rowStarts & rowStartsSmall
 317
 318   rowStarts = new int[Ndof+1];   // how many non zeros at rows of [K]
 319   rowStartsSmall = new int[NN];  // rowStarts for a piece of K(only for "u" velocity in another words 1/16 of the K(if NENv==NENp)) 
 320   checkCol = new int[1000];      // for checking the non zero column overlaps 
 321                                  // stores non-zero column number for rows (must be a large enough value)
 322   rowStarts[0] = 0;
 323
 324   for(i=0; i<NN; i++) {
 325      NNZ = 0;
 326      
 327      if(GtoL[i][0] == -1) {
 328         NNZ = 1;
 329      } else {
 330      
 331         for(k=0; k<1000; k++) {      // prepare checkCol for new row
 332            checkCol[k] = -1;
 333         }
 334      
 335         for(j=0; j<noOfColGtoL; j++) {
 336            valGtoL = GtoL[i][j];
 337            if(valGtoL != -1) {
 338               for(x=0; x<NENp; x++) {
 339                  check = 1;         // for checking if column overlap occurs or not
 340                  for(y=0; y<NNZ; y++) {
 341                     if(checkCol[y] == (LtoG[valGtoL][x])) {   // this column was created
 342                        check = 0;
 343                     }
 344                  }
 345                  if (check) {
 346                     checkCol[NNZ]=(LtoG[valGtoL][x]);         // adding new non zero number to checkCol
 347                     NNZ++;
 348                  }
 349               }
 350            }   
 351         }
 352         
 353      }
 354      
 355      rowStarts[i+1] = NNZ + rowStarts[i];
 356   }
 357
 358   // creation of rowStarts from rowStarsSmall
 359   for (i=0; i<=NN; i++) {
 360      rowStartsSmall[i] = rowStarts[i];
 361   }  
 362 
 363   for (i=1; i<=NN; i++) {
 364      rowStarts[i] = rowStarts[i]*4;
 365   } 
 366   
 367   for (i=1; i<=NN*3; i++) {
 368      rowStarts[NN+i] = rowStarts[NN] + rowStarts[i];
 369   }
 370   
 371   delete [] checkCol;
 372
 373
 374   // col & colSmall creation
 375
 376   col = new int[rowStarts[Ndof]];   // stores which non zero columns at which row data
 377   colSmall = new int[rowStarts[NN]/4]; // col for a piece of K(only for "u" velocity in another words 1/16 of the K(if NENv==NENp)) 
 378
 379   for(i=0; i<NN; i++) {
 380      NNZ = 0;
 381
 382      if(GtoL[i][0] == -1) {
 383         colSmall[rowStartsSmall[i]] = i;
 384      } else {
 385      
 386         for(j=0; j<noOfColGtoL; j++) {
 387            valGtoL = GtoL[i][j];
 388            if(valGtoL != -1) {
 389               for(x=0; x<NENp; x++) {
 390                  check = 1;
 391                  for(y=0; y<NNZ; y++) {
 392                     if(colSmall[rowStartsSmall[i]+y] == (LtoG[valGtoL][x])) {   // for checking if column overlap occurs or not
 393                        check = 0;
 394                     }
 395                  }
 396                  if (check) {
 397                     colSmall[rowStartsSmall[i]+NNZ] = (LtoG[valGtoL][x]);
 398                     NNZ++;
 399                  }
 400               }
 401            }   
 402         }
 403         
 404         for(k=1; k<NNZ; k++) {           // sorting the column vector values
 405            for(m=1; m<NNZ; m++) {        // for each row from smaller to bigger
 406               if(colSmall[rowStartsSmall[i]+m] < colSmall[rowStartsSmall[i]+m-1]) {
 407                  temp = colSmall[rowStartsSmall[i]+m];
 408                  colSmall[rowStartsSmall[i]+m] = colSmall[rowStartsSmall[i]+m-1];
 409                  colSmall[rowStartsSmall[i]+m-1] = temp;
 410               }
 411            }   
 412         }
 413
 414      }      
 415   }
 416   
 417   // creation of col from colSmall
 418   m=0;
 419   for (i=0; i<NN; i++) {
 420      for (k=0; k<(4*NN); k=k+NN) {
 421         for (j=rowStartsSmall[i]; j<rowStartsSmall[i+1]; j++) {
 422            col[m]= colSmall[j] + k;
 423            m++;
 424         }
 425      }  
 426   }   
 427   
 428   for (i=0; i<rowStarts[NN]*3; i++) {
 429      col[i+rowStarts[NN]]=col[i];
 430   }
 431
 432
 433   ////----------------------CONTROL---------------------------------------
 434
 435   //cout << endl;
 436   //for (i=0; i<5; i++) { 
 437   //   for (j=rowStartsSmall[i]; j<rowStartsSmall[i+1]; j++) {
 438   //      cout << colSmall[j] << " " ;   
 439   //   }
 440   //   cout << endl;
 441   //}
 442   //cout << endl;
 443   //cout << endl;
 444   //for (i=0; i<5; i++) { 
 445   //   for (j=rowStarts[i]; j<rowStarts[i+1]; j++) {
 446   //      cout << col[j] << " " ;   
 447   //   }
 448   //   cout << endl;
 449   //}
 450
 451   ////----------------------CONTROL---------------------------------------
 452
 453
 454   NNZ = rowStarts[Ndof];
 455
 456   //initializing val vector
 457   val = new real[rowStarts[Ndof]];
 458
 459   for(i=0; i<rowStarts[Ndof]; i++) {
 460      val[i] = 0;
 461   }
 462
 463   for (i = 0; i<NN; i++) {   //deleting the unnecessary arrays for future
 464      delete GtoL[i];
 465   }
 466   delete [] GtoL;
 467
 468} // End of function compressedSparseRowStorage()
 469
 470
 471
 472
 473//------------------------------------------------------------------------------
 474void gaussQuad()
 475//------------------------------------------------------------------------------
 476{
 477   // Generates the NGP-point Gauss quadrature points and weights.
 478
 479   GQpoint = new double*[NGP];
 480   for (int i=0; i<NGP; i++) {
 481      GQpoint[i] = new double[3];     
 482   }
 483   GQweight = new double[NGP];
 484
 485   if (eType == 3) {    // 3D Hexahedron element
 486   
 487      GQpoint[0][0] = -sqrt(1.0/3);   GQpoint[0][1] = -sqrt(1.0/3);  GQpoint[0][2] = -sqrt(1.0/3); 
 488      GQpoint[1][0] = sqrt(1.0/3);    GQpoint[1][1] = -sqrt(1.0/3);  GQpoint[1][2] = -sqrt(1.0/3);
 489      GQpoint[2][0] = sqrt(1.0/3);    GQpoint[2][1] = sqrt(1.0/3);   GQpoint[2][2] = -sqrt(1.0/3); 
 490      GQpoint[3][0] = -sqrt(1.0/3);   GQpoint[3][1] = sqrt(1.0/3);   GQpoint[3][2] = -sqrt(1.0/3); 
 491      GQpoint[4][0] = -sqrt(1.0/3);   GQpoint[4][1] = -sqrt(1.0/3);  GQpoint[4][2] = sqrt(1.0/3); 
 492      GQpoint[5][0] = sqrt(1.0/3);    GQpoint[5][1] = -sqrt(1.0/3);  GQpoint[5][2] = sqrt(1.0/3);
 493      GQpoint[6][0] = sqrt(1.0/3);    GQpoint[6][1] = sqrt(1.0/3);   GQpoint[6][2] = sqrt(1.0/3); 
 494      GQpoint[7][0] = -sqrt(1.0/3);   GQpoint[7][1] = sqrt(1.0/3);   GQpoint[7][2] = sqrt(1.0/3);      
 495      GQweight[0] = 1.0;
 496      GQweight[1] = 1.0;
 497      GQweight[2] = 1.0;
 498      GQweight[3] = 1.0;
 499      GQweight[4] = 1.0;
 500      GQweight[5] = 1.0;
 501      GQweight[6] = 1.0;
 502      GQweight[7] = 1.0;
 503       
 504   }
 505   else if (eType == 1) {       // Quadrilateral element  3 ile 4 ün yeri deðiþmeyecek mi sor!
 506      if (NGP == 1) {      // One-point quadrature
 507         GQpoint[0][0] = 0.0;            GQpoint[0][1] = 0.0;
 508         GQweight[0] = 4.0;
 509      } else if (NGP == 4) {  // Four-point quadrature
 510         GQpoint[0][0] = -sqrt(1.0/3);   GQpoint[0][1] = -sqrt(1.0/3);
 511         GQpoint[1][0] = sqrt(1.0/3);    GQpoint[1][1] = -sqrt(1.0/3);
 512         GQpoint[2][0] = -sqrt(1.0/3);   GQpoint[2][1] = sqrt(1.0/3);
 513         GQpoint[3][0] = sqrt(1.0/3);    GQpoint[3][1] = sqrt(1.0/3);
 514         GQweight[0] = 1.0;
 515         GQweight[1] = 1.0;
 516         GQweight[2] = 1.0;
 517         GQweight[3] = 1.0;
 518      } else if (NGP == 9) {  // Nine-point quadrature
 519         GQpoint[0][0] = -sqrt(3.0/5);  GQpoint[0][1] = -sqrt(3.0/5);
 520         GQpoint[1][0] = 0.0;           GQpoint[1][1] = -sqrt(3.0/5);
 521         GQpoint[2][0] = sqrt(3.0/5);   GQpoint[2][1] = -sqrt(3.0/5);
 522         GQpoint[3][0] = -sqrt(3.0/5);  GQpoint[3][1] = 0.0;
 523         GQpoint[4][0] = 0.0;           GQpoint[4][1] = 0.0;
 524         GQpoint[5][0] = sqrt(3.0/5);   GQpoint[5][1] = 0.0;
 525         GQpoint[6][0] = -sqrt(3.0/5);  GQpoint[6][1] = sqrt(3.0/5);
 526         GQpoint[7][0] = 0.0;           GQpoint[7][1] = sqrt(3.0/5);
 527         GQpoint[8][0] = sqrt(3.0/5);   GQpoint[8][1] = sqrt(3.0/5);
 528
 529         GQweight[0] = 5.0/9 * 5.0/9;
 530         GQweight[1] = 8.0/9 * 5.0/9;
 531         GQweight[2] = 5.0/9 * 5.0/9;
 532         GQweight[3] = 5.0/9 * 8.0/9;
 533         GQweight[4] = 8.0/9 * 8.0/9;
 534         GQweight[5] = 5.0/9 * 8.0/9;
 535         GQweight[6] = 5.0/9 * 5.0/9;
 536         GQweight[7] = 8.0/9 * 5.0/9;
 537         GQweight[8] = 5.0/9 * 5.0/9;
 538      } else if (NGP == 16) { // Sixteen-point quadrature
 539         GQpoint[0][0] = -0.8611363116;   GQpoint[0][1] = -0.8611363116;
 540         GQpoint[1][0] = -0.3399810435;   GQpoint[1][1] = -0.8611363116;
 541         GQpoint[2][0] =  0.3399810435;   GQpoint[2][1] = -0.8611363116;
 542         GQpoint[3][0] =  0.8611363116;   GQpoint[3][1] = -0.8611363116;
 543         GQpoint[4][0] = -0.8611363116;   GQpoint[4][1] = -0.3399810435;
 544         GQpoint[5][0] = -0.3399810435;   GQpoint[5][1] = -0.3399810435;
 545         GQpoint[6][0] =  0.3399810435;   GQpoint[6][1] = -0.3399810435;
 546         GQpoint[7][0] =  0.8611363116;   GQpoint[7][1] = -0.3399810435;
 547         GQpoint[8][0] = -0.8611363116;   GQpoint[8][1] =  0.3399810435;
 548         GQpoint[9][0] = -0.3399810435;   GQpoint[9][1] =  0.3399810435;
 549         GQpoint[10][0]=  0.3399810435;   GQpoint[10][1]=  0.3399810435;
 550         GQpoint[11][0]=  0.8611363116;   GQpoint[11][1]=  0.3399810435;
 551         GQpoint[12][0]= -0.8611363116;   GQpoint[12][1]=  0.8611363116;
 552         GQpoint[13][0]= -0.3399810435;   GQpoint[13][1]=  0.8611363116;
 553         GQpoint[14][0]=  0.3399810435;   GQpoint[14][1]=  0.8611363116;
 554         GQpoint[15][0]=  0.8611363116;   GQpoint[15][1]=  0.8611363116;
 555
 556         GQweight[0] = 0.3478548451 * 0.3478548451;
 557         GQweight[1] = 0.3478548451 * 0.6521451548;
 558         GQweight[2] = 0.3478548451 * 0.6521451548;
 559         GQweight[3] = 0.3478548451 * 0.3478548451;
 560         GQweight[4] = 0.6521451548 * 0.3478548451;
 561         GQweight[5] = 0.6521451548 * 0.6521451548;
 562         GQweight[6] = 0.6521451548 * 0.6521451548;
 563         GQweight[7] = 0.6521451548 * 0.3478548451;
 564         GQweight[8] = 0.6521451548 * 0.3478548451;
 565         GQweight[9] = 0.6521451548 * 0.6521451548;
 566         GQweight[10] = 0.6521451548 * 0.6521451548;
 567         GQweight[11] = 0.6521451548 * 0.3478548451;
 568         GQweight[12] = 0.3478548451 * 0.3478548451;
 569         GQweight[13] = 0.3478548451 * 0.6521451548;
 570         GQweight[14] = 0.3478548451 * 0.6521451548;
 571         GQweight[15] = 0.3478548451 * 0.3478548451;
 572      }
 573   } else if (eType == 2) {  // Triangular element
 574      if (NGP == 1) {          // One-point quadrature
 575         GQpoint[0][0] = 1.0/3;  GQpoint[0][1] = 1.0/3;
 576         GQweight[0] = 0.5;
 577      } else if (NGP == 3) {   // Two-point quadrature
 578         GQpoint[0][0] = 0.5;   GQpoint[0][1] = 0.0;
 579         GQpoint[1][0] = 0.0;   GQpoint[1][1] = 0.5;
 580         GQpoint[2][0] = 0.5;   GQpoint[2][1] = 0.5;
 581
 582         GQweight[0] = 1.0/6;
 583         GQweight[1] = 1.0/6;
 584         GQweight[2] = 1.0/6;
 585      } else if (NGP == 4) {   // Four-point quadrature
 586         GQpoint[0][0] = 1.0/3;   GQpoint[0][1] = 1.0/3;
 587         GQpoint[1][0] = 0.6;     GQpoint[1][1] = 0.2;
 588         GQpoint[2][0] = 0.2;     GQpoint[2][1] = 0.6;
 589         GQpoint[3][0] = 0.2;     GQpoint[3][1] = 0.2;
 590
 591         GQweight[0] = -27.0/96;
 592         GQweight[1] = 25.0/96;
 593         GQweight[2] = 25.0/96;
 594         GQweight[3] = 25.0/96;
 595      } else if (NGP == 7) {  // Seven-point quadrature
 596         GQpoint[0][0] = 1.0/3;               GQpoint[0][1] = 1.0/3;
 597         GQpoint[1][0] = 0.059715871789770;   GQpoint[1][1] = 0.470142064105115;
 598         GQpoint[2][0] = 0.470142064105115;   GQpoint[2][1] = 0.059715871789770;
 599         GQpoint[3][0] = 0.470142064105115;   GQpoint[3][1] = 0.470142064105115;
 600         GQpoint[4][0] = 0.101286507323456;   GQpoint[4][1] = 0.797426985353087;
 601         GQpoint[5][0] = 0.101286507323456;   GQpoint[5][1] = 0.101286507323456;
 602         GQpoint[6][0] = 0.797426985353087;   GQpoint[6][1] = 0.101286507323456;
 603
 604         GQweight[0] = 0.225 / 2;
 605         GQweight[1] = 0.132394152788 / 2;
 606         GQweight[2] = 0.132394152788 / 2;
 607         GQweight[3] = 0.132394152788 / 2;
 608         GQweight[4] = 0.125939180544 / 2;
 609         GQweight[5] = 0.125939180544 / 2;
 610         GQweight[6] = 0.125939180544 / 2;
 611      }
 612   }
 613
 614} // End of function gaussQuad()
 615
 616
 617
 618
 619//------------------------------------------------------------------------------
 620void calcShape()
 621//------------------------------------------------------------------------------
 622{
 623   // Calculates the values of the shape functions and their derivatives with
 624   // respect to ksi and eta at GQ points.
 625
 626   double ksi, eta, zeta;
 627   
 628   if (eType == 3) { // 3D Hexahedron Element
 629      if (NENp == 8) {
 630      
 631         Sp = new double*[8];
 632         for (int i=0; i<NGP; i++) {
 633            Sp[i] = new double[NGP];
 634         }
 635         
 636         DSp = new double **[3];	
 637         for (int i=0; i<3; i++) {
 638            DSp[i] = new double *[8];
 639            for (int j=0; j<8; j++) {
 640               DSp[i][j] = new double[NGP];
 641            }
 642         }     
 643         
 644         for (int k = 0; k<NGP; k++) {
 645            ksi = GQpoint[k][0];
 646            eta = GQpoint[k][1];
 647            zeta = GQpoint[k][2];
 648            
 649            Sp[0][k] = 0.25*(1-ksi)*(1-eta)*(1-zeta);
 650            Sp[1][k] = 0.25*(1+ksi)*(1-eta)*(1-zeta);
 651            Sp[2][k] = 0.25*(1+ksi)*(1+eta)*(1-zeta);
 652            Sp[3][k] = 0.25*(1-ksi)*(1+eta)*(1-zeta);   
 653            Sp[4][k] = 0.25*(1-ksi)*(1-eta)*(1+zeta);
 654            Sp[5][k] = 0.25*(1+ksi)*(1-eta)*(1+zeta);
 655            Sp[6][k] = 0.25*(1+ksi)*(1+eta)*(1+zeta);
 656            Sp[7][k] = 0.25*(1-ksi)*(1+eta)*(1+zeta); 
 657
 658            DSp[0][0][k] = -0.25*(1-eta)*(1-zeta);  // ksi derivative of the 1st shape funct. at k-th GQ point.
 659            DSp[1][0][k] = -0.25*(1-ksi)*(1-zeta);  // eta derivative of the 1st shape funct. at k-th GQ point.  
 660            DSp[2][0][k] = -0.25*(1-ksi)*(1-eta);   // zeta derivative of the 1st shape funct. at k-th GQ point.
 661            DSp[0][1][k] =  0.25*(1-eta)*(1-zeta);   
 662            DSp[1][1][k] = -0.25*(1+ksi)*(1-zeta);
 663            DSp[2][1][k] = -0.25*(1+ksi)*(1-eta);         
 664            DSp[0][2][k] =  0.25*(1+eta)*(1-zeta);   
 665            DSp[1][2][k] =  0.25*(1+ksi)*(1-zeta);
 666            DSp[2][2][k] = -0.25*(1+ksi)*(1+eta);         
 667            DSp[0][3][k] = -0.25*(1+eta)*(1-zeta);   
 668            DSp[1][3][k] =  0.25*(1-ksi)*(1-zeta);
 669            DSp[2][3][k] = -0.25*(1-ksi)*(1+eta); 
 670            DSp[0][4][k] = -0.25*(1-eta)*(1+zeta);  // ksi derivative of the 1st shape funct. at k-th GQ point.
 671            DSp[1][4][k] = -0.25*(1-ksi)*(1+zeta);  // eta derivative of the 1st shape funct. at k-th GQ point.  
 672            DSp[2][4][k] = 0.25*(1-ksi)*(1-eta);   // zeta derivative of the 1st shape funct. at k-th GQ point.
 673            DSp[0][5][k] =  0.25*(1-eta)*(1+zeta);   
 674            DSp[1][5][k] = -0.25*(1+ksi)*(1+zeta);
 675            DSp[2][5][k] = 0.25*(1+ksi)*(1-eta);         
 676            DSp[0][6][k] =  0.25*(1+eta)*(1+zeta);   
 677            DSp[1][6][k] =  0.25*(1+ksi)*(1+zeta);
 678            DSp[2][6][k] = 0.25*(1+ksi)*(1+eta);         
 679            DSp[0][7][k] = -0.25*(1+eta)*(1+zeta);   
 680            DSp[1][7][k] =  0.25*(1-ksi)*(1+zeta);
 681            DSp[2][7][k] = 0.25*(1-ksi)*(1+eta); 
 682         }
 683      }
 684      
 685      if (NENv == 8) {
 686      
 687         Sv = new double*[8];
 688         for (int i=0; i<NGP; i++) {
 689            Sv[i] = new double[NGP];
 690         }
 691      
 692         DSv = new double **[3];	
 693         for (int i=0; i<3; i++) {
 694            DSv[i] = new double *[8];
 695            for (int j=0; j<8; j++) {
 696               DSv[i][j] = new double[NGP];
 697            }
 698         }
 699         
 700         for (int i=0; i<NGP; i++) {
 701            Sv[i] = Sp[i];
 702         }
 703         
 704         for (int i=0; i<3; i++) {
 705            for (int j=0; j<8; j++) {
 706               DSv[i][j] = DSp[i][j];
 707            }
 708         }         
 709               
 710      }
 711   }   
 712      
 713   for (int i = 0; i<8; i++) {   //deleting the unnecessary arrays for future
 714      delete GQpoint[i];
 715   }
 716   delete [] GQpoint;
 717      
 718} // End of function calcShape()
 719
 720
 721
 722
 723//------------------------------------------------------------------------------
 724void calcJacobian()
 725//------------------------------------------------------------------------------
 726{
 727   // Calculates the Jacobian matrix, its inverse and determinant for all
 728   // elements at all GQ points. Also evaluates and stores derivatives of shape
 729   // functions wrt global coordinates x and y.
 730
 731   int e, i, j, k, x, iG; 
 732   double **e_coord;
 733   double **Jacob, **invJacob;
 734   double temp;
 735   
 736   if (eType == 3) {
 737   
 738      e_coord = new double*[NENp];
 739   
 740      for (i=0; i<NENp; i++) {
 741         e_coord[i] = new double[3];
 742      }
 743   
 744      Jacob = new double*[3];
 745      invJacob = new double*[3];
 746      for (i=0; i<3; i++) {
 747         Jacob[i] = new double[3];
 748         invJacob[i] = new double[3];
 749      }
 750      
 751      detJacob = new double*[NE];
 752      for (i=0; i<NE; i++) {
 753         detJacob[i] = new double[NGP];
 754      }
 755      
 756      gDSp = new double***[NE];
 757
 758      for (i=0; i<NE; i++) {
 759         gDSp[i] = new double**[3];
 760         for(j=0; j<3; j++) {
 761            gDSp[i][j] = new double*[NENv];
 762            for(k=0; k<NENp; k++) {
 763               gDSp[i][j][k] = new double[NGP];
 764            }
 765         }	
 766      }
 767      
 768      gDSv = new double***[NE];
 769
 770      for (i=0; i<NE; i++) {
 771         gDSv[i] = new double**[3];
 772         for(j=0; j<3; j++) {
 773            gDSv[i][j] = new double*[NENv];
 774            for(k=0; k<NENv; k++) {
 775               gDSv[i][j][k] = new double[NGP];
 776            }
 777         }	
 778      }      
 779   
 780      for (e = 0; e<NE; e++){
 781         // To calculate Jacobian matrix of an element we need e_coord matrix of
 782         // size NEN*2. Each row of it stores x, y & z coordinates of the nodes of
 783         // an element.
 784         for (i = 0; i<NENp; i++){
 785            iG = LtoG[e][i];
 786            e_coord[i][0] = coord[iG][0]; 
 787            e_coord[i][1] = coord[iG][1];
 788            e_coord[i][2] = coord[iG][2];
 789         }
 790      
 791         // For each GQ point calculate 3*3 Jacobian matrix, its inverse and its
 792         // determinant. Also calculate derivatives of shape functions wrt global
 793         // coordinates x and y & z. These are the derivatives that we'll use in
 794         // evaluating K and F integrals. 
 795   	
 796         for (k = 0; k<NGP; k++) {
 797            for (i = 0; i<3; i++) {
 798               for (j = 0; j<3; j++) {
 799                  temp = 0;
 800                  for (x = 0; x<NENp; x++) {
 801                     temp = DSp[i][x][k] * e_coord[x][j]+temp;
 802                  }
 803                  Jacob[i][j] = temp;
 804               }
 805            }
 806            
 807            invJacob[0][0] = Jacob[1][1]*Jacob[2][2]-Jacob[2][1]*Jacob[1][2];
 808            invJacob[0][1] = -(Jacob[0][1]*Jacob[2][2]-Jacob[0][2]*Jacob[2][1]);
 809            invJacob[0][2] = Jacob[0][1]*Jacob[1][2]-Jacob[1][1]*Jacob[0][2];
 810            invJacob[1][0] =  -(Jacob[1][0]*Jacob[2][2]-Jacob[1][2]*Jacob[2][0]);
 811            invJacob[1][1] = Jacob[2][2]*Jacob[0][0]-Jacob[2][0]*Jacob[0][2];
 812            invJacob[1][2] = -(Jacob[1][2]*Jacob[0][0]-Jacob[1][0]*Jacob[0][2]);
 813            invJacob[2][0] = Jacob[1][0]*Jacob[2][1]-Jacob[2][0]*Jacob[1][1];
 814            invJacob[2][1] = -(Jacob[2][1]*Jacob[0][0]-Jacob[2][0]*Jacob[0][1]);
 815            invJacob[2][2] = Jacob[1][1]*Jacob[0][0]-Jacob[1][0]*Jacob[0][1];
 816
 817            detJacob[e][k] = Jacob[0][0]*(Jacob[1][1]*Jacob[2][2]-Jacob[2][1]*Jacob[1][2]) +
 818                            Jacob[0][1]*(Jacob[1][2]*Jacob[2][0]-Jacob[1][0]*Jacob[2][2]) +
 819                            Jacob[0][2]*(Jacob[1][0]*Jacob[2][1]-Jacob[1][1]*Jacob[2][0]);
 820            
 821            for (i = 0; i<3; i++){
 822               for (j = 0; j<3; j++){
 823                  invJacob[i][j] = invJacob[i][j]/detJacob[e][k];
 824               }    
 825            }
 826         
 827            for (i = 0; i<3; i++){
 828               for (j = 0; j<8; j++){
 829                  temp = 0;
 830                  for (x = 0; x<3; x++){
 831                     temp = invJacob[i][x] * DSp[x][j][k]+temp;
 832                  }
 833                  gDSp[e][i][j][k] = temp;
 834                  gDSv[e][i][j][k] = gDSp[e][i][j][k];
 835               }
 836            }
 837         }
 838      }
 839   } 
 840   
 841}  // End of function calcJacobian()
 842
 843
 844
 845
 846//------------------------------------------------------------------------------
 847void initGlobalSysVariables()
 848//------------------------------------------------------------------------------
 849{
 850   // Allocating the memory for stiffness matrices and force vectors
 851   
 852   int i;
 853   
 854   F = new real[Ndof];
 855   K = new real*[Ndof];	
 856   for (i=0; i<Ndof; i++) {
 857      K[i] = new real[Ndof];
 858   }
 859   
 860   Fe = new double[NEU];
 861   Ke = new double*[NEU];	
 862   for (i=0; i<NEU; i++) {
 863      Ke[i] = new double[NEU];
 864   }
 865   
 866   Fe_1 = new double[NENv];
 867   Fe_2 = new double[NENv];
 868   Fe_3 = new double[NENv];
 869   Fe_4 = new double[NENp];
 870   
 871   Ke_11 = new double*[NENv];
 872   Ke_12 = new double*[NENv];
 873	Ke_13 = new double*[NENv];
 874   Ke_14 = new double*[NENv];
 875   // Ke_21 is the transpose of Ke_12
 876   Ke_22 = new double*[NENv];
 877   Ke_23 = new double*[NENv];
 878	Ke_24 = new double*[NENv]; 
 879   // Ke_31 is the transpose of Ke_13
 880   // Ke_32 is the transpose of Ke_23
 881	Ke_33 = new double*[NENv];
 882	Ke_34 = new double*[NENv];
 883   Ke_41 = new double*[NENv];
 884   Ke_42 = new double*[NENv];
 885	Ke_43 = new double*[NENv];
 886	Ke_44 = new double*[NENp];
 887   
 888   Ke_11_add = new double*[NENv];
 889   Ke_12_add = new double*[NENv];
 890   Ke_13_add = new double*[NENv];
 891   Ke_14_add = new double*[NENv];
 892   Ke_22_add = new double*[NENv];
 893   Ke_23_add = new double*[NENv];
 894   Ke_24_add = new double*[NENv];
 895   Ke_33_add = new double*[NENv];
 896   Ke_34_add = new double*[NENv];
 897   Ke_41_add = new double*[NENp];
 898   Ke_42_add = new double*[NENp];
 899   Ke_43_add = new double*[NENp];
 900   Ke_44_add = new double*[NENp];
 901   for (i=0; i<NENv; i++) {
 902      Ke_11[i] = new double[NENv];
 903      Ke_12[i] = new double[NENv];
 904      Ke_13[i] = new double[NENv];
 905      Ke_14[i] = new double[NENp];
 906      Ke_22[i] = new double[NENv];
 907      Ke_23[i] = new double[NENv];
 908      Ke_24[i] = new double[NENp];
 909      Ke_33[i] = new double[NENv];
 910      Ke_34[i] = new double[NENp];
 911      
 912      Ke_11_add[i] = new double[NENv];
 913      Ke_12_add[i] = new double[NENv];
 914      Ke_13_add[i] = new double[NENv];
 915      Ke_14_add[i] = new double[NENp];
 916      Ke_22_add[i] = new double[NENv];
 917      Ke_23_add[i] = new double[NENv];
 918      Ke_24_add[i] = new double[NENp];
 919      Ke_33_add[i] = new double[NENv];
 920      Ke_34_add[i] = new double[NENp];
 921   }
 922   for (i=0; i<NENp; i++) {
 923      Ke_41[i] = new double[NENv];
 924      Ke_42[i] = new double[NENv];
 925	   Ke_43[i] = new double[NENv];
 926      Ke_44[i] = new double[NENp]; 
 927      Ke_41_add[i] = new double[NENp];
 928      Ke_42_add[i] = new double[NENp];
 929      Ke_43_add[i] = new double[NENp];
 930      Ke_44_add[i] = new double[NENp];
 931   }   
 932   
 933   uNodal = new double[NENv];
 934   vNodal = new double[NENv];
 935   wNodal = new double[NENv];
 936   
 937}   
 938
 939
 940
 941
 942//------------------------------------------------------------------------------
 943void calcGlobalSys()
 944//------------------------------------------------------------------------------
 945{
 946   // Calculates Ke and Fe one by one for each element and assembles them into
 947   // the global K and F.
 948
 949   int e, i, j, k, m, n, node;
 950   double Tau;
 951   //   double x, y, axy, fxy;
 952
 953   // Initialize the arrays
 954
 955   for (i=0; i<Ndof; i++) {
 956      F[i] = 0;
 957      for (j=0; j<Ndof; j++) {
 958         K[i][j] = 0;
 959	   }
 960   }
 961   
 962   // Calculating the elemental stiffness matrix(Ke) and force vector(Fe)
 963
 964   for (e = 0; e<NE; e++) {
 965      // Intitialize Ke and Fe to zero.
 966      for (i=0; i<NEU; i++) {
 967         Fe[i] = 0;
 968         for (j=0; j<NEU; j++) {
 969            Ke[i][j] = 0;
 970         }
 971      }
 972
 973      for (i=0; i<NENv; i++) {
 974         Fe_1[i] = 0;
 975         Fe_2[i] = 0;
 976         Fe_3[i] = 0;
 977         Fe_4[i] = 0;
 978         for (j=0; j<NENv; j++) {
 979            Ke_11[i][j] = 0;
 980            Ke_12[i][j] = 0;
 981            Ke_13[i][j] = 0;
 982            Ke_22[i][j] = 0;
 983            Ke_23[i][j] = 0;
 984            Ke_33[i][j] = 0;
 985         }
 986         for (j=0; j<NENp; j++) {
 987            Ke_14[i][j] = 0;
 988            Ke_24[i][j] = 0;
 989            Ke_34[i][j] = 0;
 990         }         
 991      }   
 992      for (i=0; i<NENp; i++) {
 993         for (j=0; j<NENp; j++) {
 994            Ke_41[i][j] = 0;
 995            Ke_42[i][j] = 0;
 996            Ke_43[i][j] = 0;
 997            Ke_44[i][j] = 0;
 998         }
 999      }
1000      
1001      for (i=0; i<NENv; i++) {
1002         uNodal[i] = u[LtoG[e][i]]; 
1003         vNodal[i] = u[LtoG[e][i]+NN];
1004         wNodal[i] = u[LtoG[e][i]+2*NN];
1005      }
1006      
1007      for (k = 0; k<NGP; k++) {   // Gauss quadrature loop
1008            
1009         u0 = 0;
1010         v0 = 0;
1011         w0 = 0;
1012         
1013         for (i=0; i<NENv; i++) {
1014            u0 = u0 + Sp[i][k] * uNodal[i];
1015            v0 = v0 + Sp[i][k] * vNodal[i];
1016            w0 = w0 + Sp[i][k] * wNodal[i];            
1017         }
1018         
1019         for (i=0; i<NENv; i++) {
1020            for (j=0; j<NENv; j++) {
1021               Ke_11_add[i][j] = 0;
1022               Ke_12_add[i][j] = 0;
1023               Ke_13_add[i][j] = 0;
1024               Ke_22_add[i][j] = 0;
1025               Ke_23_add[i][j] = 0;
1026               Ke_33_add[i][j] = 0;
1027            }
1028            for (j=0; j<NENp; j++) {
1029               Ke_14_add[i][j] = 0;
1030               Ke_24_add[i][j] = 0;
1031               Ke_34_add[i][j] = 0;
1032            }         
1033         }
1034         for (i=0; i<NENp; i++) {
1035            for (j=0; j<NENp; j++) {
1036               Ke_41_add[i][j] = 0;
1037               Ke_42_add[i][j] = 0;
1038               Ke_43_add[i][j] = 0;
1039               Ke_44_add[i][j] = 0;
1040            }
1041         }
1042
1043         for (i=0; i<NENv; i++) {
1044            for (j=0; j<NENv; j++) {
1045               Ke_11_add[i][j] = Ke_11_add[i][j] + viscosity * (
1046                     2 * gDSv[e][0][i][k] * gDSv[e][0][j][k] + 
1047                     gDSv[e][1][i][k] * gDSv[e][1][j][k] +
1048                     gDSv[e][2][i][k] * gDSv[e][2][j][k]) +
1049                     density * Sv[i][k] * (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);
1050                     //density * Sv[i][k];
1051                     
1052               Ke_12_add[i][j] = Ke_12_add[i][j] + viscosity * gDSv[e][1][i][k] * gDSv[e][0][j][k];
1053               
1054               Ke_13_add[i][j] = Ke_13_add[i][j] + viscosity * gDSv[e][2][i][k] * gDSv[e][0][j][k];
1055               
1056               Ke_22_add[i][j] = Ke_22_add[i][j] + viscosity * (
1057                     gDSv[e][0][i][k] * gDSv[e][0][j][k] + 
1058                     2 * gDSv[e][1][i][k] * gDSv[e][1][j][k] +
1059                     gDSv[e][2][i][k] * gDSv[e][2][j][k]) +
1060                     density * Sv[i][k] * (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);
1061                     
1062               Ke_23_add[i][j] = Ke_23_add[i][j] + viscosity * gDSv[e][2][i][k] * gDSv[e][1][j][k];
1063               
1064               Ke_33_add[i][j] = Ke_33_add[i][j] + viscosity * (
1065                     gDSv[e][0][i][k] * gDSv[e][0][j][k] + 
1066                     gDSv[e][1][i][k] * gDSv[e][1][j][k] +
1067                     2 * gDSv[e][2][i][k] * gDSv[e][2][j][k]) +
1068                     density * Sv[i][k] * (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);
1069                     
1070            }
1071            for (j=0; j<NENp; j++) {
1072               Ke_14_add[i][j] = Ke_14_add[i][j] - gDSv[e][0][i][k] * Sp[j][k];
1073               Ke_24_add[i][j] = Ke_24_add[i][j] - gDSv[e][1][i][k] * Sp[j][k];
1074               Ke_34_add[i][j] = Ke_34_add[i][j] - gDSv[e][2][i][k] * Sp[j][k];
1075            }         
1076         }
1077         
1078         for (i=0; i<NENv; i++) {
1079            for (j=0; j<NENp; j++) {
1080               Ke_41_add[j][i] = Ke_14_add[i][j] ;
1081               Ke_42_add[j][i] = Ke_24_add[i][j] ;
1082               Ke_43_add[j][i] = Ke_34_add[i][j] ;
1083            }
1084         }
1085
1086  // Apply GLS stabilization for linear elements with NENv = NENp
1087    
1088         Tau = ((1.0/12.0)*2) / viscosity;  // GLS parameter
1089
1090         for (i=0; i<NENv; i++) {
1091            for (j=0; j<NENv; j++) {
1092              Ke_11_add[i][j] = Ke_11_add[i][j] + Tau * density * density *
1093                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) *
1094                             (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);                             
1095             
1096              Ke_22_add[i][j] = Ke_22_add[i][j] + Tau * density * density * 
1097                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) *
1098                             (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);   
1099                             
1100              Ke_33_add[i][j] = Ke_33_add[i][j] + Tau * density * density * 
1101                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) *
1102                             (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);                                
1103            } 
1104         } 
1105
1106         for (i=0; i<NENv; i++) {
1107            for (j=0; j<NENp; j++) {
1108              Ke_14_add[i][j] = Ke_14_add[i][j] + Tau * density * 
1109                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) * gDSp[e][0][j][k]; 
1110             
1111              Ke_24_add[i][j] = Ke_24_add[i][j] + Tau * density * 
1112                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) * gDSp[e][1][j][k]; 
1113                             
1114              Ke_34_add[i][j] = Ke_34_add[i][j] + Tau * density * 
1115                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) * gDSp[e][2][j][k];                              
1116            } 
1117         } 
1118
1119         for (i=0; i<NENp; i++) {
1120            for (j=0; j<NENv; j++) {
1121              Ke_41_add[i][j] = Ke_41_add[i][j] - Tau * density * 
1122                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) * gDSp[e][0][j][k]; 
1123             
1124              Ke_42_add[i][j] = Ke_42_add[i][j] - Tau * density * 
1125                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) * gDSp[e][1][j][k]; 
1126                             
1127              Ke_43_add[i][j] = Ke_43_add[i][j] - Tau * density * 
1128                             (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) * gDSp[e][2][j][k];        
1129            } 
1130         } 
1131         
1132         for (i=0; i<NENp; i++) {
1133            for (j=0; j<NENp; j++) {
1134               Ke_44_add[i][j] = Ke_44_add[i][j] - Tau *
1135                              ( gDSv[e][0][i][k] * gDSv[e][0][j][k] + gDSv[e][1][i][k] * gDSv[e][1][j][k] + gDSv[e][2][i][k] * gDSv[e][2][j][k] );
1136            }
1137         }
1138
1139  // Apply GLS stabilization for linear elements with NENv = NENp
1140            
1141         for (i=0; i<NENv; i++) {
1142            for (j=0; j<NENv; j++) {            
1143               Ke_11[i][j] += Ke_11_add[i][j] * detJacob[e][k] * GQweight[k];
1144               Ke_12[i][j] += Ke_12_add[i][j] * detJacob[e][k] * GQweight[k];
1145               Ke_13[i][j] += Ke_13_add[i][j] * detJacob[e][k] * GQweight[k];
1146               Ke_22[i][j] += Ke_22_add[i][j] * detJacob[e][k] * GQweight[k];
1147               Ke_23[i][j] += Ke_23_add[i][j] * detJacob[e][k] * GQweight[k];
1148               Ke_33[i][j] += Ke_33_add[i][j] * detJacob[e][k] * GQweight[k];
1149            }
1150            for (j=0; j<NENp; j++) {               
1151               Ke_14[i][j] += Ke_14_add[i][j] * detJacob[e][k] * GQweight[k];
1152               Ke_24[i][j] += Ke_24_add[i][j] * detJacob[e][k] * GQweight[k];
1153               Ke_34[i][j] += Ke_34_add[i][j] * detJacob[e][k] * GQweight[k];
1154            }
1155         }
1156         for (i=0; i<NENv; i++) {
1157            for (j=0; j<NENv; j++) {
1158               Ke_41[i][j] += Ke_41_add[i][j] * detJacob[e][k] * GQweight[k];
1159               Ke_42[i][j] += Ke_42_add[i][j] * detJacob[e][k] * GQweight[k];
1160               Ke_43[i][j] += Ke_43_add[i][j] * detJacob[e][k] * GQweight[k];
1161               Ke_44[i][j] += Ke_44_add[i][j] * detJacob[e][k] * GQweight[k];
1162            }
1163         }
1164
1165      }   // End GQ loop  
1166
1167      // Assembly of Fe
1168
1169      i=0;
1170      for (j=0; j<NENv; j++) {
1171         Fe[i]=Fe_1[j];
1172         i++;
1173      }
1174      for (j=0; j<NENv; j++) {
1175         Fe[i]=Fe_2[j];
1176         i++;
1177      }
1178      for (j=0; j<NENv; j++) {
1179         Fe[i]=Fe_3[j];
1180         i++;
1181      }
1182      for (j=0; j<NENp; j++) {
1183         Fe[i]=Fe_4[j];
1184         i++;
1185      }         
1186
1187      // Assembly of Ke
1188
1189      i=0;
1190      j=0;    
1191      for (m=0; m<NENv; m++) {
1192         j=0;
1193         for (n=0; n<NENv; n++) {       
1194            Ke[i][j] =  Ke_11[m][n];
1195            j++;
1196         }
1197         i++;
1198      }
1199      
1200      i=i-NENv;
1201      for (m=0; m<NENv; m++) {
1202         j=0;
1203         for (n=0; n<NENv; n++) {       
1204            Ke[i][j+NENv] =  Ke_12[m][n];
1205            j++;
1206         }
1207         i++;
1208      }
1209
1210      i=i-NENv;
1211      for (m=0; m<NENv; m++) {
1212         j=0;
1213         for (n=0; n<NENv; n++) {       
1214            Ke[i][j+2*NENv] =  Ke_13[m][n];
1215            j++;
1216         }
1217         i++;
1218      }
1219      
1220      i=i-NENv;
1221      for (m=0; m<NENv; m++) {
1222         j=0;
1223         for (n=0; n<NENp; n++) {       
1224            Ke[i][j+3*NENv] =  Ke_14[m][n];
1225            j++;
1226         }
1227         i++;
1228      }
1229      
1230      for (m=0; m<NENv; m++) {
1231         j=0;
1232         for (n=0; n<NENv; n++) {  
1233            Ke[i][j] =  Ke_12[n][m];
1234            j++;
1235         }
1236         i++;
1237      }
1238
1239      i=i-NENv;
1240      for (m=0; m<NENv; m++) {
1241         j=0;
1242         for (n=0; n<NENv; n++) {       
1243            Ke[i][j+NENv] =  Ke_22[m][n];
1244            j++;
1245         }
1246         i++;
1247      }
1248
1249      i=i-NENv;
1250      for (m=0; m<NENv; m++) {
1251         j=0;
1252         for (n=0; n<NENv; n++) {       
1253            Ke[i][j+2*NENv] =  Ke_23[m][n];
1254            j++;
1255         }
1256         i++;
1257      }
1258
1259      i=i-NENv;
1260      for (m=0; m<NENv; m++) {
1261         j=0;
1262         for (n=0; n<NENp; n++) {       
1263            Ke[i][j+3*NENv] =  Ke_24[m][n];
1264            j++;
1265         }
1266         i++;
1267      }
1268
1269      for (m=0; m<NENv; m++) {
1270         j=0;
1271         for (n=0; n<NENv; n++) {       
1272            Ke[i][j] =  Ke_13[n][m];
1273            j++;
1274         }
1275         i++;
1276      }      
1277
1278      i=i-NENv;
1279      for (m=0; m<NENv; m++) {
1280         j=0;
1281         for (n=0; n<NENv; n++) {       
1282            Ke[i][j+NENv] =  Ke_23[n][m];
1283            j++;
1284         }
1285         i++;
1286      } 
1287
1288      i=i-NENv;
1289      for (m=0; m<NENv; m++) {
1290         j=0;
1291         for (n=0; n<NENv; n++) {       
1292            Ke[i][j+2*NENv] =  Ke_33[m][n];
1293            j++;
1294         }
1295         i++;
1296      }
1297
1298      i=i-NENv;
1299      for (m=0; m<NENv; m++) {
1300         j=0;
1301         for (n=0; n<NENp; n++) {       
1302            Ke[i][j+3*NENv] =  Ke_34[m][n];
1303            j++;
1304         }
1305         i++;
1306      }
1307
1308      j=0;
1309      for (m=0; m<NENv; m++) {
1310         j=0;
1311         for (n=0; n<NENp; n++) {       
1312            Ke[i][j] =  Ke_41[m][n];
1313            j++;
1314         }
1315         i++;
1316      }
1317      
1318      i=i-NENv;
1319      for (m=0; m<NENv; m++) {
1320         j=0;
1321         for (n=0; n<NENp; n++) {       
1322            Ke[i][j+NENv] =  Ke_42[m][n];
1323            j++;
1324         }
1325         i++;
1326      } 
1327
1328      i=i-NENv;
1329      for (m=0; m<NENv; m++) {
1330         j=0;
1331         for (n=0; n<NENp; n++) {       
1332            Ke[i][j+2*NENv] =  Ke_43[m][n];
1333            j++;
1334         }
1335         i++;
1336      } 
1337
1338      i=i-NENv;
1339      for (m=0; m<NENp; m++) {
1340         j=0;
1341         for (n=0; n<NENp; n++) {       
1342            Ke[i][j+3*NENv] = Ke_44[m][n];
1343            j++;
1344         }
1345         i++;
1346      }
1347
1348      assemble(e, Ke, Fe);  // sending Ke & Fe for assembly
1349
1350
1351   }   // End element loop
1352 
1353} // End of function calcGlobalSys()
1354
1355
1356
1357
1358//------------------------------------------------------------------------------
1359void assemble(int e, double **Ke, double *Fe)
1360//------------------------------------------------------------------------------
1361{
1362   // Inserts Ke and Fe into proper locations of K and F.
1363
1364   int i, j, m, n, iG, jG; 
1365   
1366   // Inserts [Ke] and {Fe} into proper locations of [K] and {F} by the help of LtoG array. 
1367   
1368   for (i = 0; i<NEU; i++) {
1369         
1370      iG = LtoG[e][i%NENv] + NN*(i/NENv);
1371      F[iG] = F[iG] + Fe[i];
1372      
1373      for (j = 0; j<NEU; j++) {   
1374         jG = LtoG[e][j%NENv] + NN*(j/NENv);
1375         K[iG][jG] = K[iG][jG] + Ke[i][j];
1376      }
1377   }   
1378
1379   // Assembly process for compressed sparse storage 
1380   // KeKMapSmall creation
1381
1382   int shiftRow, shiftCol;
1383   int *nodeData, p, q, k;
1384   // int *nodeDataEBC;
1385    
1386   // nodeDataEBC = new int[NENv];     // Elemental node data(LtoG data) (modified with EBC)    // BC implementationu sonraya býrak ineff ama bir anda zor
1387   nodeData = new int[NENv];           // Stores sorted LtoG data
1388   for(k=0; k<NENv; k++) {
1389      nodeData[k] = (LtoG[e][k]);      // Takes node data from LtoG
1390   } 
1391
1392   // for(i=0; i<NEN; i++) {
1393      // nodeData[i]=nodeDataEBC[i];
1394      // if (GtoL[nodeDataEBC[i]][0] == -1) {
1395         // val[rowStarts[nodeDataEBC[i]]] = 1*bigNumber;      // Must be fixed, val[x] = 1 repeating for every element that contains the node!
1396         // nodeDataEBC[i] = -1;
1397      // }
1398   // }
1399
1400   KeKMapSmall = new int*[NENv];        // NENv X NENv                                          
1401   for(j=0; j<NENv; j++) {              // Stores map data between K elemental and value vector
1402      KeKMapSmall[j] = new int[NENv];
1403   }
1404
1405   for(i=0; i<NENv; i++) {
1406      for(j=0; j<NENv; j++) {
1407         q=0;
1408         for(p=rowStartsSmall[nodeData[i]]; p<rowStartsSmall[nodeData[i]+1]; p++) {  // p is the location of the col vector(col[x], p=x) 
1409            if(colSmall[p] == nodeData[j]) {                                         // Selection process of the KeKMapSmall data from the col vector
1410               KeKMapSmall[i][j] = q; 
1411               break;
1412            }
1413            q++;
1414         }
1415      }
1416   }
1417
1418   // Creating val vector
1419   for(shiftRow=1; shiftRow<5; shiftRow++) {
1420      for(shiftCol=1; shiftCol<5; shift

Large files files are truncated, but you can click here to view the full file