PageRenderTime 290ms CodeModel.GetById 41ms app.highlight 215ms RepoModel.GetById 16ms app.codeStats 1ms

/Stokes3D.cpp

https://bitbucket.org/cuneytsert/cfd-with-cuda
C++ | 1751 lines | 1265 code | 315 blank | 171 comment | 278 complexity | 68b7412e7727b9a392fc1c580105e83e 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;        // nonlinearIterMax is not used by the Stokes solver.
  43double density, viscosity, fx, fy, nonlinearTol, solverTol;  // nonlinearTol is not used by the Stokes solver.
  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;   // Can be float or double. K is the stiffness matrix in full
  53                    // storage used Gauss Elimination solver.
  54int bigNumber;
  55
  56int **GtoL, *rowStarts, *rowStartsSmall, *colSmall, *col, **KeKMapSmall, NNZ; 
  57real *val;          // Can be float or double
  58
  59void readInput();
  60void gaussQuad();
  61void calcShape();
  62void calcJacobian();
  63void calcGlobalSys();
  64void assemble(int e, double **Ke, double *Fe);
  65void applyBC();
  66void solve();
  67void postProcess();
  68void gaussElimination(int N, real **K, real *F, real *u, bool& err);
  69void writeTecplotFile();
  70void compressedSparseRowStorage();
  71#ifdef CUSP

  72   extern void CUSPsolver();
  73#endif

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

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