PageRenderTime 78ms CodeModel.GetById 39ms RepoModel.GetById 1ms app.codeStats 0ms

/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. 

  2. *        This code is a part of the CFD-with-CUDA project        *

    3. 

  3. *             http://code.google.com/p/cfd-with-cuda             *

    4. 

  4. *                                                                *

    5. 

  5. *              Dr. Cuneyt Sert and Mahmut M. Gocmen              *

    6. 

  6. *                                                                *

    7. 

  7. *              Department of Mechanical Engineering              *

    8. 

  8. *                Middle East Technical University                *

    9. 

  9. *                         Ankara, Turkey                         *

    10. 

  10. *            http://www.me.metu.edu.tr/people/cuneyt             *

    11. 

  11. *****************************************************************/

    12. 

  12. 

    13. 

  13. #include <stdio.h>

    14. 

  14. #include <string>

    15. 

  15. #include <iostream>

    16. 

  16. #include <iomanip>

    17. 

  17. #include <fstream>

    18. 

  18. #include <ctime>

    19. 

  19. #include <cmath>

    20. 

  20. #include <algorithm>

    21. 

  21. 

    22. 

  22. using namespace std;

    23. 

  23. 

    24. 

  24. #ifdef SINGLE

    25. 

  25.   typedef float real;

    26. 

  26. #else

    27. 

  27.   typedef double real;

    28. 

  28. #endif

    29. 

  29. 

    30. 

  30. ifstream meshfile;

    31. 

  31. ifstream problemFile;

    32. 

  32. ofstream outputFile;

    33. 

  33. ofstream outputControl;

    34. 

  34. 

    35. 

  35. string problemNameFile, whichProblem;

    36. 

  36. string problemName = "ProblemName.txt";

    37. 

  37. string controlFile = "Control_Output.txt";

    38. 

  38. string inputExtension  = ".inp";              // Change this line to specify the extension of input file (with dot).

    39. 

  39. string outputExtension  = ".dat";             // Change this line to specify the extension of output file (with dot).

    40. 

  40. 

    41. 

  41. int name, eType, NE, NN, NGP, NEU, Ndof;

    42. 

  42. int NCN, NENv, NENp, nonlinearIterMax, solverIterMax;

    43. 

  43. double density, viscosity, fx, fy, nonlinearTol, solverTol;

    44. 

  44. int **LtoG, **velNodes, **pressureNodes;

    45. 

  45. double **coord;

    46. 

  46. int nBC, nVelNodes, nPressureNodes;

    47. 

  47. double *BCtype, **BCstrings;

    48. 

  48. double axyFunc, fxyFunc;

    49. 

  49. double **GQpoint, *GQweight;

    50. 

  50. double **Sp, ***DSp, **Sv, ***DSv;

    51. 

  51. double **detJacob, ****gDSp, ****gDSv;

    52. 

  52. real **K, *F, *u, *uOld;   // Can be float or double. K is the stiffness matrix

    53. 

  53.                            // in full storage used Gauss Elimination solver.

    54. 

  54. int bigNumber;

    55. 

  55. 

    56. 

  56. double *Fe, **Ke;

    57. 

  57. double *Fe_1, *Fe_2, *Fe_3, *Fe_4;

    58. 

  58. double **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;

    59. 

  59. double **Ke_11_add, **Ke_12_add, **Ke_13_add, **Ke_14_add, **Ke_22_add, **Ke_23_add, **Ke_24_add;

    60. 

  60. double **Ke_33_add, **Ke_34_add, **Ke_41_add, **Ke_42_add, **Ke_43_add, **Ke_44_add;

    61. 

  61. double *uNodal, *vNodal, *wNodal;

    62. 

  62. double u0, v0, w0;

    63. 

  63. 

    64. 

  64. int **GtoL, *rowStarts, *rowStartsSmall, *colSmall, *col, **KeKMapSmall, NNZ; 

    65. 

  65. real *val;

    66. 

  66. 

    67. 

  67. void readInput();

    68. 

  68. void gaussQuad();

    69. 

  69. void calcShape();

    70. 

  70. void calcJacobian();

    71. 

  71. void initGlobalSysVariables();

    72. 

  72. void calcGlobalSys();

    73. 

  73. void assemble(int e, double **Ke, double *Fe);

    74. 

  74. void applyBC();

    75. 

  75. void solve();

    76. 

  76. void postProcess();

    77. 

  77. void gaussElimination(int N, real **K, real *F, real *u, bool& err);

    78. 

  78. void writeTecplotFile();

    79. 

  79. void compressedSparseRowStorage();

    80. 

  80. #ifdef CUSP

    81. 

  81.    extern void CUSPsolver();

    82. 

  82. #endif

    83. 

  83. 

    84. 

  84. 

    85. 

  85. 

    86. 

  86. //------------------------------------------------------------------------------

    87. 

  87. int main()

    88. 

  88. //------------------------------------------------------------------------------

    89. 

  89. {

    90. 

  90.    cout << "\n********************************************************";

    91. 

  91.    cout << "\n*   Navier-Stokes 3D - Part of CFD with CUDA project   *";

    92. 

  92.    cout << "\n*       http://code.google.com/p/cfd-with-cuda         *";

    93. 

  93.    cout << "\n********************************************************\n\n";

    94. 

  94. 

    95. 

  95.    cout << "The program is started." << endl ;

    96. 

  96. 

    97. 

  97.    time_t start, end;

    98. 

  98.    time (&start);     // Start measuring execution time.

    99. 

  99. 

    100. 

  100.    readInput();                   cout << "Input file is read." << endl ;

    101. 

  101.    compressedSparseRowStorage();  cout << "CSR vectors are created." << endl ;

    102. 

  102.    gaussQuad();

    103. 

  103.    calcShape();

    104. 

  104.    calcJacobian();

    105. 

  105.    initGlobalSysVariables();

    106. 

  106.    solve();

    107. 

  107.    //postProcess();

    108. 

  108.    writeTecplotFile();            cout << "A DAT file is created for Tecplot." << endl ;

    109. 

  109. 

    110. 

  110.    time (&end);      // Stop measuring execution time.

    111. 

  111.    cout << endl << "Elapsed wall clock time is " << difftime (end,start) << " seconds." << endl;

    112. 

  112.    

    113. 

  113.    cout << endl << "The program is terminated successfully.\nPress a key to close this window...";

    114. 

  114. 

    115. 

  115.    cin.get();

    116. 

  116.    return 0;

    117. 

  117. 

    118. 

  118. } // End of function main()

    119. 

  119. 

    120. 

  120. 

    121. 

  121. 

    122. 

  122. 

    123. 

  123. //------------------------------------------------------------------------------

    124. 

  124. void readInput()

    125. 

  125. //------------------------------------------------------------------------------

    126. 

  126. {

    127. 

  127. 

    128. 

  128.    string dummy, dummy2, dummy4, dummy5;

    129. 

  129.    int dummy3, i;

    130. 

  130. 

    131. 

  131.    problemFile.open(problemName.c_str(), ios::in);   // Read problem name

    132. 

  132.    problemFile >> whichProblem;                                  

    133. 

  133.    problemFile.close();

    134. 

  134.    

    135. 

  135.    meshfile.open((whichProblem + inputExtension).c_str(), ios::in);

    136. 

  136.      

    137. 

  137.    meshfile.ignore(256, '\n'); // Read and ignore the line

    138. 

  138.    meshfile.ignore(256, '\n'); // Read and ignore the line

    139. 

  139.    meshfile >> dummy >> dummy2 >> eType;

    140. 

  140.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    141. 

  141.    meshfile >> dummy >> dummy2 >> NE;

    142. 

  142.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    143. 

  143.    meshfile >> dummy >> dummy2 >> NCN;

    144. 

  144.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    145. 

  145.    meshfile >> dummy >> dummy2 >> NN;

    146. 

  146.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    147. 

  147.    meshfile >> dummy >> dummy2 >> NENv;

    148. 

  148.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    149. 

  149.    meshfile >> dummy >> dummy2 >> NENp;   

    150. 

  150.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    151. 

  151.    meshfile >> dummy >> dummy2 >> NGP;

    152. 

  152.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    153. 

  153.    meshfile >> dummy >> dummy2 >> nonlinearIterMax;

    154. 

  154.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    155. 

  155.    meshfile >> dummy >> dummy2 >> nonlinearTol;

    156. 

  156.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    157. 

  157.    meshfile >> dummy >> dummy2 >> solverIterMax;

    158. 

  158.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    159. 

  159.    meshfile >> dummy >> dummy2 >> solverTol;

    160. 

  160.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    161. 

  161.    meshfile >> dummy >> dummy2 >> density;

    162. 

  162.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    163. 

  163.    meshfile >> dummy >> dummy2 >> viscosity;

    164. 

  164.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    165. 

  165.    meshfile >> dummy >> dummy2 >> fx;

    166. 

  166.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    167. 

  167.    meshfile >> dummy >> dummy2 >> fy;

    168. 

  168.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    169. 

  169.    meshfile.ignore(256, '\n'); // Read and ignore the line

    170. 

  170.    meshfile.ignore(256, '\n'); // Read and ignore the line  

    171. 

  171.    

    172. 

  172.    NEU = 3*NENv + NENp;

    173. 

  173.    Ndof = 3*NN + NCN;

    174. 

  174.    

    175. 

  175.    // Read node coordinates

    176. 

  176.    coord = new double*[NN];

    177. 

  177.    

    178. 

  178.    if (eType == 2 || eType == 1) {

    179. 

  179.       for (i = 0; i < NN; i++) {

    180. 

  180.          coord[i] = new double[2];

    181. 

  181.       }

    182. 

  182.    	

    183. 

  183.       for (i=0; i<NN; i++){

    184. 

  184.          meshfile >> dummy3 >> coord[i][0] >> coord[i][1] ;

    185. 

  185.          meshfile.ignore(256, '\n'); // Ignore the rest of the line    

    186. 

  186.       }  

    187. 

  187.    } else {

    188. 

  188.       for (i = 0; i < NN; i++) {

    189. 

  189.          coord[i] = new double[3];

    190. 

  190.       }

    191. 

  191.    	

    192. 

  192.       for (i=0; i<NN; i++){

    193. 

  193.          meshfile >> dummy3 >> coord[i][0] >> coord[i][1] >> coord[i][2] ;

    194. 

  194.          meshfile.ignore(256, '\n'); // Ignore the rest of the line    

    195. 

  195.       }  

    196. 

  196.    }

    197. 

  197.    

    198. 

  198.    meshfile.ignore(256, '\n'); // Read and ignore the line

    199. 

  199.    meshfile.ignore(256, '\n'); // Read and ignore the line 

    200. 

  200.     

    201. 

  201.    // Read element connectivity, i.e. LtoG

    202. 

  202.    LtoG = new int*[NE];

    203. 

  203. 

    204. 

  204.    for (i=0; i<NE; i++) {

    205. 

  205.       LtoG[i] = new int[NENv];

    206. 

  206.    }

    207. 

  207. 

    208. 

  208.    for (int e = 0; e<NE; e++){

    209. 

  209.       meshfile >> dummy3;

    210. 

  210.       for (i = 0; i<NENv; i++){

    211. 

  211.          meshfile >> LtoG[e][i];

    212. 

  212.       }

    213. 

  213.       meshfile.ignore(256, '\n'); // Read and ignore the line 

    214. 

  214.    } 

    215. 

  215. 

    216. 

  216.    meshfile.ignore(256, '\n'); // Read and ignore the line

    217. 

  217.    meshfile.ignore(256, '\n'); // Read and ignore the line 

    218. 

  218. 

    219. 

  219.    // Read number of different BC types and details of each BC

    220. 

  220.    meshfile >> dummy >> dummy2 >> nBC;

    221. 

  221.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    222. 

  222.       

    223. 

  223.    BCtype = new double[nBC];

    224. 

  224.    BCstrings = new double*[nBC];

    225. 

  225.    for (i=0; i<nBC; i++) {

    226. 

  226.       BCstrings[i] = new double[3];

    227. 

  227.    }   

    228. 

  228.     

    229. 

  229.    for (i = 0; i<nBC-1; i++){

    230. 

  230.       meshfile >> dummy >> BCtype[i] >> dummy2 >> dummy3 >> BCstrings[i][0] >> dummy4 >> BCstrings[i][1] >> dummy5 >> BCstrings[i][2];

    231. 

  231.       meshfile.ignore(256, '\n'); // Ignore the rest of the line

    232. 

  232.    }

    233. 

  233.    meshfile >> dummy >> BCtype[2]  >> dummy >> BCstrings[i][0];

    234. 

  234.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    235. 

  235.    // Read EBC data 

    236. 

  236.     

    237. 

  237.    meshfile.ignore(256, '\n'); // Read and ignore the line 

    238. 

  238.    meshfile >> dummy >> dummy2 >> nVelNodes;

    239. 

  239.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    240. 

  240.    meshfile >> dummy >> dummy2 >> nPressureNodes;

    241. 

  241.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    242. 

  242.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    243. 

  243.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    244. 

  244.       

    245. 

  245.    if (nVelNodes!=0){

    246. 

  246.       velNodes = new int*[nVelNodes];

    247. 

  247.       for (i = 0; i < nVelNodes; i++){

    248. 

  248.          velNodes[i] = new int[2];

    249. 

  249.       }

    250. 

  250.       for (i = 0; i < nVelNodes; i++){

    251. 

  251.          meshfile >> velNodes[i][0] >> velNodes[i][1];   

    252. 

  252.          meshfile.ignore(256, '\n'); // Ignore the rest of the line

    253. 

  253.       }

    254. 

  254.    }

    255. 

  255. 

    256. 

  256.    meshfile.ignore(256, '\n'); // Ignore the rest of the line  

    257. 

  257.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    258. 

  258.    meshfile.ignore(256, '\n'); // Ignore the rest of the line

    259. 

  259.    

    260. 

  260.    if (nVelNodes!=0){

    261. 

  261.       pressureNodes = new int*[nPressureNodes];

    262. 

  262.       for (i = 0; i < nPressureNodes; i++){

    263. 

  263.          pressureNodes[i] = new int[2];

    264. 

  264.       }

    265. 

  265.       for (i = 0; i < nPressureNodes; i++){

    266. 

  266.          meshfile >> pressureNodes[i][0] >> pressureNodes[i][1];   

    267. 

  267.          meshfile.ignore(256, '\n'); // Ignore the rest of the line

    268. 

  268.       }

    269. 

  269.    }

    270. 

  270.   

    271. 

  271.    meshfile.close();

    272. 

  272. 

    273. 

  273. } // End of function readInput()

    274. 

  274. 

    275. 

  275. 

    276. 

  276. 

    277. 

  277. 

    278. 

  278. //------------------------------------------------------------------------------

    279. 

  279. void compressedSparseRowStorage()

    280. 

  280. //------------------------------------------------------------------------------

    281. 

  281. {

    282. 

  282. //GtoL creation

    283. 

  283. int i, j, k, m, x, y, valGtoL, check, temp, *checkCol, noOfColGtoL, *GtoLCounter;

    284. 

  284. 

    285. 

  285.    GtoL = new int*[NN];          // stores which elements connected to the node

    286. 

  286.    noOfColGtoL = 8;              // for 3D meshes created by our mesh generator max 8 elements connect to one node,

    287. 

  287.    GtoLCounter = new int[NN];    // but for real 3D problems this must be a big enough value!

    288. 

  288.    

    289. 

  289.    for (i=0; i<NN; i++) {     

    290. 

  290.       GtoL[i] = new int[noOfColGtoL];   

    291. 

  291.    }

    292. 

  292. 

    293. 

  293.    for(i=0; i<NN; i++) {

    294. 

  294.       for(j=0; j<noOfColGtoL; j++) {

    295. 

  295. 		   GtoL[i][j] = -1;      // for the nodes that didn't connect 4 nodes

    296. 

  296. 	   }

    297. 

  297.       GtoLCounter[i] = 0;

    298. 

  298.    }

    299. 

  299.    

    300. 

  300.    for (i=0; i<NE; i++) {

    301. 

  301.       for(j=0; j<NENv; j++) {

    302. 

  302.          GtoL[LtoG[i][j]][GtoLCounter[LtoG[i][j]]] = i;

    303. 

  303.          GtoLCounter[LtoG[i][j]] += 1;

    304. 

  304.       }

    305. 

  305.    }         

    306. 

  306.    

    307. 

  307.    delete [] GtoLCounter;

    308. 

  308.    // Su anda gereksiz 0'lar oluþturuluyor. Ýleride düzeltilebilir.

    309. 

  309.    // for(i=0; i<nVelNodes; i++) {   // extracting EBC values from GtoL with making them "-1"

    310. 

  310.       // for(j=0; j<noOfColGtoL; j++) {

    311. 

  311.          // GtoL[velNodes[i][0]][j] = -1;

    312. 

  312.       // }   

    313. 

  313.    // } 

    314. 

  314. 

    315. 

  315. 

    316. 

  316. // Finding size of col vector, creation of rowStarts & rowStartsSmall

    317. 

  317. 

    318. 

  318.    rowStarts = new int[Ndof+1];   // how many non zeros at rows of [K]

    319. 

  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. 

  320.    checkCol = new int[1000];      // for checking the non zero column overlaps 

    321. 

  321.                                   // stores non-zero column number for rows (must be a large enough value)

    322. 

  322.    rowStarts[0] = 0;

    323. 

  323. 

    324. 

  324.    for(i=0; i<NN; i++) {

    325. 

  325.       NNZ = 0;

    326. 

  326.       

    327. 

  327.       if(GtoL[i][0] == -1) {

    328. 

  328.          NNZ = 1;

    329. 

  329.       } else {

    330. 

  330.       

    331. 

  331.          for(k=0; k<1000; k++) {      // prepare checkCol for new row

    332. 

  332.             checkCol[k] = -1;

    333. 

  333.          }

    334. 

  334.       

    335. 

  335.          for(j=0; j<noOfColGtoL; j++) {

    336. 

  336.             valGtoL = GtoL[i][j];

    337. 

  337.             if(valGtoL != -1) {

    338. 

  338.                for(x=0; x<NENp; x++) {

    339. 

  339.                   check = 1;         // for checking if column overlap occurs or not

    340. 

  340.                   for(y=0; y<NNZ; y++) {

    341. 

  341.                      if(checkCol[y] == (LtoG[valGtoL][x])) {   // this column was created

    342. 

  342.                         check = 0;

    343. 

  343.                      }

    344. 

  344.                   }

    345. 

  345.                   if (check) {

    346. 

  346.                      checkCol[NNZ]=(LtoG[valGtoL][x]);         // adding new non zero number to checkCol

    347. 

  347.                      NNZ++;

    348. 

  348.                   }

    349. 

  349.                }

    350. 

  350.             }   

    351. 

  351.          }

    352. 

  352.          

    353. 

  353.       }

    354. 

  354.       

    355. 

  355.       rowStarts[i+1] = NNZ + rowStarts[i];

    356. 

  356.    }

    357. 

  357. 

    358. 

  358.    // creation of rowStarts from rowStarsSmall

    359. 

  359.    for (i=0; i<=NN; i++) {

    360. 

  360.       rowStartsSmall[i] = rowStarts[i];

    361. 

  361.    }  

    362. 

  362.  

    363. 

  363.    for (i=1; i<=NN; i++) {

    364. 

  364.       rowStarts[i] = rowStarts[i]*4;

    365. 

  365.    } 

    366. 

  366.    

    367. 

  367.    for (i=1; i<=NN*3; i++) {

    368. 

  368.       rowStarts[NN+i] = rowStarts[NN] + rowStarts[i];

    369. 

  369.    }

    370. 

  370.    

    371. 

  371.    delete [] checkCol;

    372. 

  372. 

    373. 

  373. 

    374. 

  374.    // col & colSmall creation

    375. 

  375. 

    376. 

  376.    col = new int[rowStarts[Ndof]];   // stores which non zero columns at which row data

    377. 

  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. 

  378. 

    379. 

  379.    for(i=0; i<NN; i++) {

    380. 

  380.       NNZ = 0;

    381. 

  381. 

    382. 

  382.       if(GtoL[i][0] == -1) {

    383. 

  383.          colSmall[rowStartsSmall[i]] = i;

    384. 

  384.       } else {

    385. 

  385.       

    386. 

  386.          for(j=0; j<noOfColGtoL; j++) {

    387. 

  387.             valGtoL = GtoL[i][j];

    388. 

  388.             if(valGtoL != -1) {

    389. 

  389.                for(x=0; x<NENp; x++) {

    390. 

  390.                   check = 1;

    391. 

  391.                   for(y=0; y<NNZ; y++) {

    392. 

  392.                      if(colSmall[rowStartsSmall[i]+y] == (LtoG[valGtoL][x])) {   // for checking if column overlap occurs or not

    393. 

  393.                         check = 0;

    394. 

  394.                      }

    395. 

  395.                   }

    396. 

  396.                   if (check) {

    397. 

  397.                      colSmall[rowStartsSmall[i]+NNZ] = (LtoG[valGtoL][x]);

    398. 

  398.                      NNZ++;

    399. 

  399.                   }

    400. 

  400.                }

    401. 

  401.             }   

    402. 

  402.          }

    403. 

  403.          

    404. 

  404.          for(k=1; k<NNZ; k++) {           // sorting the column vector values

    405. 

  405.             for(m=1; m<NNZ; m++) {        // for each row from smaller to bigger

    406. 

  406.                if(colSmall[rowStartsSmall[i]+m] < colSmall[rowStartsSmall[i]+m-1]) {

    407. 

  407.                   temp = colSmall[rowStartsSmall[i]+m];

    408. 

  408.                   colSmall[rowStartsSmall[i]+m] = colSmall[rowStartsSmall[i]+m-1];

    409. 

  409.                   colSmall[rowStartsSmall[i]+m-1] = temp;

    410. 

  410.                }

    411. 

  411.             }   

    412. 

  412.          }

    413. 

  413. 

    414. 

  414.       }      

    415. 

  415.    }

    416. 

  416.    

    417. 

  417.    // creation of col from colSmall

    418. 

  418.    m=0;

    419. 

  419.    for (i=0; i<NN; i++) {

    420. 

  420.       for (k=0; k<(4*NN); k=k+NN) {

    421. 

  421.          for (j=rowStartsSmall[i]; j<rowStartsSmall[i+1]; j++) {

    422. 

  422.             col[m]= colSmall[j] + k;

    423. 

  423.             m++;

    424. 

  424.          }

    425. 

  425.       }  

    426. 

  426.    }   

    427. 

  427.    

    428. 

  428.    for (i=0; i<rowStarts[NN]*3; i++) {

    429. 

  429.       col[i+rowStarts[NN]]=col[i];

    430. 

  430.    }

    431. 

  431. 

    432. 

  432. 

    433. 

  433.    ////----------------------CONTROL---------------------------------------

    434. 

  434. 

    435. 

  435.    //cout << endl;

    436. 

  436.    //for (i=0; i<5; i++) { 

    437. 

  437.    //   for (j=rowStartsSmall[i]; j<rowStartsSmall[i+1]; j++) {

    438. 

  438.    //      cout << colSmall[j] << " " ;   

    439. 

  439.    //   }

    440. 

  440.    //   cout << endl;

    441. 

  441.    //}

    442. 

  442.    //cout << endl;

    443. 

  443.    //cout << endl;

    444. 

  444.    //for (i=0; i<5; i++) { 

    445. 

  445.    //   for (j=rowStarts[i]; j<rowStarts[i+1]; j++) {

    446. 

  446.    //      cout << col[j] << " " ;   

    447. 

  447.    //   }

    448. 

  448.    //   cout << endl;

    449. 

  449.    //}

    450. 

  450. 

    451. 

  451.    ////----------------------CONTROL---------------------------------------

    452. 

  452. 

    453. 

  453. 

    454. 

  454.    NNZ = rowStarts[Ndof];

    455. 

  455. 

    456. 

  456.    //initializing val vector

    457. 

  457.    val = new real[rowStarts[Ndof]];

    458. 

  458. 

    459. 

  459.    for(i=0; i<rowStarts[Ndof]; i++) {

    460. 

  460.       val[i] = 0;

    461. 

  461.    }

    462. 

  462. 

    463. 

  463.    for (i = 0; i<NN; i++) {   //deleting the unnecessary arrays for future

    464. 

  464.       delete GtoL[i];

    465. 

  465.    }

    466. 

  466.    delete [] GtoL;

    467. 

  467. 

    468. 

  468. } // End of function compressedSparseRowStorage()

    469. 

  469. 

    470. 

  470. 

    471. 

  471. 

    472. 

  472. 

    473. 

  473. //------------------------------------------------------------------------------

    474. 

  474. void gaussQuad()

    475. 

  475. //------------------------------------------------------------------------------

    476. 

  476. {

    477. 

  477.    // Generates the NGP-point Gauss quadrature points and weights.

    478. 

  478. 

    479. 

  479.    GQpoint = new double*[NGP];

    480. 

  480.    for (int i=0; i<NGP; i++) {

    481. 

  481.       GQpoint[i] = new double[3];     

    482. 

  482.    }

    483. 

  483.    GQweight = new double[NGP];

    484. 

  484. 

    485. 

  485.    if (eType == 3) {    // 3D Hexahedron element

    486. 

  486.    

    487. 

  487.       GQpoint[0][0] = -sqrt(1.0/3);   GQpoint[0][1] = -sqrt(1.0/3);  GQpoint[0][2] = -sqrt(1.0/3); 

    488. 

  488.       GQpoint[1][0] = sqrt(1.0/3);    GQpoint[1][1] = -sqrt(1.0/3);  GQpoint[1][2] = -sqrt(1.0/3);

    489. 

  489.       GQpoint[2][0] = sqrt(1.0/3);    GQpoint[2][1] = sqrt(1.0/3);   GQpoint[2][2] = -sqrt(1.0/3); 

    490. 

  490.       GQpoint[3][0] = -sqrt(1.0/3);   GQpoint[3][1] = sqrt(1.0/3);   GQpoint[3][2] = -sqrt(1.0/3); 

    491. 

  491.       GQpoint[4][0] = -sqrt(1.0/3);   GQpoint[4][1] = -sqrt(1.0/3);  GQpoint[4][2] = sqrt(1.0/3); 

    492. 

  492.       GQpoint[5][0] = sqrt(1.0/3);    GQpoint[5][1] = -sqrt(1.0/3);  GQpoint[5][2] = sqrt(1.0/3);

    493. 

  493.       GQpoint[6][0] = sqrt(1.0/3);    GQpoint[6][1] = sqrt(1.0/3);   GQpoint[6][2] = sqrt(1.0/3); 

    494. 

  494.       GQpoint[7][0] = -sqrt(1.0/3);   GQpoint[7][1] = sqrt(1.0/3);   GQpoint[7][2] = sqrt(1.0/3);      

    495. 

  495.       GQweight[0] = 1.0;

    496. 

  496.       GQweight[1] = 1.0;

    497. 

  497.       GQweight[2] = 1.0;

    498. 

  498.       GQweight[3] = 1.0;

    499. 

  499.       GQweight[4] = 1.0;

    500. 

  500.       GQweight[5] = 1.0;

    501. 

  501.       GQweight[6] = 1.0;

    502. 

  502.       GQweight[7] = 1.0;

    503. 

  503.        

    504. 

  504.    }

    505. 

  505.    else if (eType == 1) {       // Quadrilateral element  3 ile 4 ün yeri deðiþmeyecek mi sor!

    506. 

  506.       if (NGP == 1) {      // One-point quadrature

    507. 

  507.          GQpoint[0][0] = 0.0;            GQpoint[0][1] = 0.0;

    508. 

  508.          GQweight[0] = 4.0;

    509. 

  509.       } else if (NGP == 4) {  // Four-point quadrature

    510. 

  510.          GQpoint[0][0] = -sqrt(1.0/3);   GQpoint[0][1] = -sqrt(1.0/3);

    511. 

  511.          GQpoint[1][0] = sqrt(1.0/3);    GQpoint[1][1] = -sqrt(1.0/3);

    512. 

  512.          GQpoint[2][0] = -sqrt(1.0/3);   GQpoint[2][1] = sqrt(1.0/3);

    513. 

  513.          GQpoint[3][0] = sqrt(1.0/3);    GQpoint[3][1] = sqrt(1.0/3);

    514. 

  514.          GQweight[0] = 1.0;

    515. 

  515.          GQweight[1] = 1.0;

    516. 

  516.          GQweight[2] = 1.0;

    517. 

  517.          GQweight[3] = 1.0;

    518. 

  518.       } else if (NGP == 9) {  // Nine-point quadrature

    519. 

  519.          GQpoint[0][0] = -sqrt(3.0/5);  GQpoint[0][1] = -sqrt(3.0/5);

    520. 

  520.          GQpoint[1][0] = 0.0;           GQpoint[1][1] = -sqrt(3.0/5);

    521. 

  521.          GQpoint[2][0] = sqrt(3.0/5);   GQpoint[2][1] = -sqrt(3.0/5);

    522. 

  522.          GQpoint[3][0] = -sqrt(3.0/5);  GQpoint[3][1] = 0.0;

    523. 

  523.          GQpoint[4][0] = 0.0;           GQpoint[4][1] = 0.0;

    524. 

  524.          GQpoint[5][0] = sqrt(3.0/5);   GQpoint[5][1] = 0.0;

    525. 

  525.          GQpoint[6][0] = -sqrt(3.0/5);  GQpoint[6][1] = sqrt(3.0/5);

    526. 

  526.          GQpoint[7][0] = 0.0;           GQpoint[7][1] = sqrt(3.0/5);

    527. 

  527.          GQpoint[8][0] = sqrt(3.0/5);   GQpoint[8][1] = sqrt(3.0/5);

    528. 

  528. 

    529. 

  529.          GQweight[0] = 5.0/9 * 5.0/9;

    530. 

  530.          GQweight[1] = 8.0/9 * 5.0/9;

    531. 

  531.          GQweight[2] = 5.0/9 * 5.0/9;

    532. 

  532.          GQweight[3] = 5.0/9 * 8.0/9;

    533. 

  533.          GQweight[4] = 8.0/9 * 8.0/9;

    534. 

  534.          GQweight[5] = 5.0/9 * 8.0/9;

    535. 

  535.          GQweight[6] = 5.0/9 * 5.0/9;

    536. 

  536.          GQweight[7] = 8.0/9 * 5.0/9;

    537. 

  537.          GQweight[8] = 5.0/9 * 5.0/9;

    538. 

  538.       } else if (NGP == 16) { // Sixteen-point quadrature

    539. 

  539.          GQpoint[0][0] = -0.8611363116;   GQpoint[0][1] = -0.8611363116;

    540. 

  540.          GQpoint[1][0] = -0.3399810435;   GQpoint[1][1] = -0.8611363116;

    541. 

  541.          GQpoint[2][0] =  0.3399810435;   GQpoint[2][1] = -0.8611363116;

    542. 

  542.          GQpoint[3][0] =  0.8611363116;   GQpoint[3][1] = -0.8611363116;

    543. 

  543.          GQpoint[4][0] = -0.8611363116;   GQpoint[4][1] = -0.3399810435;

    544. 

  544.          GQpoint[5][0] = -0.3399810435;   GQpoint[5][1] = -0.3399810435;

    545. 

  545.          GQpoint[6][0] =  0.3399810435;   GQpoint[6][1] = -0.3399810435;

    546. 

  546.          GQpoint[7][0] =  0.8611363116;   GQpoint[7][1] = -0.3399810435;

    547. 

  547.          GQpoint[8][0] = -0.8611363116;   GQpoint[8][1] =  0.3399810435;

    548. 

  548.          GQpoint[9][0] = -0.3399810435;   GQpoint[9][1] =  0.3399810435;

    549. 

  549.          GQpoint[10][0]=  0.3399810435;   GQpoint[10][1]=  0.3399810435;

    550. 

  550.          GQpoint[11][0]=  0.8611363116;   GQpoint[11][1]=  0.3399810435;

    551. 

  551.          GQpoint[12][0]= -0.8611363116;   GQpoint[12][1]=  0.8611363116;

    552. 

  552.          GQpoint[13][0]= -0.3399810435;   GQpoint[13][1]=  0.8611363116;

    553. 

  553.          GQpoint[14][0]=  0.3399810435;   GQpoint[14][1]=  0.8611363116;

    554. 

  554.          GQpoint[15][0]=  0.8611363116;   GQpoint[15][1]=  0.8611363116;

    555. 

  555. 

    556. 

  556.          GQweight[0] = 0.3478548451 * 0.3478548451;

    557. 

  557.          GQweight[1] = 0.3478548451 * 0.6521451548;

    558. 

  558.          GQweight[2] = 0.3478548451 * 0.6521451548;

    559. 

  559.          GQweight[3] = 0.3478548451 * 0.3478548451;

    560. 

  560.          GQweight[4] = 0.6521451548 * 0.3478548451;

    561. 

  561.          GQweight[5] = 0.6521451548 * 0.6521451548;

    562. 

  562.          GQweight[6] = 0.6521451548 * 0.6521451548;

    563. 

  563.          GQweight[7] = 0.6521451548 * 0.3478548451;

    564. 

  564.          GQweight[8] = 0.6521451548 * 0.3478548451;

    565. 

  565.          GQweight[9] = 0.6521451548 * 0.6521451548;

    566. 

  566.          GQweight[10] = 0.6521451548 * 0.6521451548;

    567. 

  567.          GQweight[11] = 0.6521451548 * 0.3478548451;

    568. 

  568.          GQweight[12] = 0.3478548451 * 0.3478548451;

    569. 

  569.          GQweight[13] = 0.3478548451 * 0.6521451548;

    570. 

  570.          GQweight[14] = 0.3478548451 * 0.6521451548;

    571. 

  571.          GQweight[15] = 0.3478548451 * 0.3478548451;

    572. 

  572.       }

    573. 

  573.    } else if (eType == 2) {  // Triangular element

    574. 

  574.       if (NGP == 1) {          // One-point quadrature

    575. 

  575.          GQpoint[0][0] = 1.0/3;  GQpoint[0][1] = 1.0/3;

    576. 

  576.          GQweight[0] = 0.5;

    577. 

  577.       } else if (NGP == 3) {   // Two-point quadrature

    578. 

  578.          GQpoint[0][0] = 0.5;   GQpoint[0][1] = 0.0;

    579. 

  579.          GQpoint[1][0] = 0.0;   GQpoint[1][1] = 0.5;

    580. 

  580.          GQpoint[2][0] = 0.5;   GQpoint[2][1] = 0.5;

    581. 

  581. 

    582. 

  582.          GQweight[0] = 1.0/6;

    583. 

  583.          GQweight[1] = 1.0/6;

    584. 

  584.          GQweight[2] = 1.0/6;

    585. 

  585.       } else if (NGP == 4) {   // Four-point quadrature

    586. 

  586.          GQpoint[0][0] = 1.0/3;   GQpoint[0][1] = 1.0/3;

    587. 

  587.          GQpoint[1][0] = 0.6;     GQpoint[1][1] = 0.2;

    588. 

  588.          GQpoint[2][0] = 0.2;     GQpoint[2][1] = 0.6;

    589. 

  589.          GQpoint[3][0] = 0.2;     GQpoint[3][1] = 0.2;

    590. 

  590. 

    591. 

  591.          GQweight[0] = -27.0/96;

    592. 

  592.          GQweight[1] = 25.0/96;

    593. 

  593.          GQweight[2] = 25.0/96;

    594. 

  594.          GQweight[3] = 25.0/96;

    595. 

  595.       } else if (NGP == 7) {  // Seven-point quadrature

    596. 

  596.          GQpoint[0][0] = 1.0/3;               GQpoint[0][1] = 1.0/3;

    597. 

  597.          GQpoint[1][0] = 0.059715871789770;   GQpoint[1][1] = 0.470142064105115;

    598. 

  598.          GQpoint[2][0] = 0.470142064105115;   GQpoint[2][1] = 0.059715871789770;

    599. 

  599.          GQpoint[3][0] = 0.470142064105115;   GQpoint[3][1] = 0.470142064105115;

    600. 

  600.          GQpoint[4][0] = 0.101286507323456;   GQpoint[4][1] = 0.797426985353087;

    601. 

  601.          GQpoint[5][0] = 0.101286507323456;   GQpoint[5][1] = 0.101286507323456;

    602. 

  602.          GQpoint[6][0] = 0.797426985353087;   GQpoint[6][1] = 0.101286507323456;

    603. 

  603. 

    604. 

  604.          GQweight[0] = 0.225 / 2;

    605. 

  605.          GQweight[1] = 0.132394152788 / 2;

    606. 

  606.          GQweight[2] = 0.132394152788 / 2;

    607. 

  607.          GQweight[3] = 0.132394152788 / 2;

    608. 

  608.          GQweight[4] = 0.125939180544 / 2;

    609. 

  609.          GQweight[5] = 0.125939180544 / 2;

    610. 

  610.          GQweight[6] = 0.125939180544 / 2;

    611. 

  611.       }

    612. 

  612.    }

    613. 

  613. 

    614. 

  614. } // End of function gaussQuad()

    615. 

  615. 

    616. 

  616. 

    617. 

  617. 

    618. 

  618. 

    619. 

  619. //------------------------------------------------------------------------------

    620. 

  620. void calcShape()

    621. 

  621. //------------------------------------------------------------------------------

    622. 

  622. {

    623. 

  623.    // Calculates the values of the shape functions and their derivatives with

    624. 

  624.    // respect to ksi and eta at GQ points.

    625. 

  625. 

    626. 

  626.    double ksi, eta, zeta;

    627. 

  627.    

    628. 

  628.    if (eType == 3) { // 3D Hexahedron Element

    629. 

  629.       if (NENp == 8) {

    630. 

  630.       

    631. 

  631.          Sp = new double*[8];

    632. 

  632.          for (int i=0; i<NGP; i++) {

    633. 

  633.             Sp[i] = new double[NGP];

    634. 

  634.          }

    635. 

  635.          

    636. 

  636.          DSp = new double **[3];	

    637. 

  637.          for (int i=0; i<3; i++) {

    638. 

  638.             DSp[i] = new double *[8];

    639. 

  639.             for (int j=0; j<8; j++) {

    640. 

  640.                DSp[i][j] = new double[NGP];

    641. 

  641.             }

    642. 

  642.          }     

    643. 

  643.          

    644. 

  644.          for (int k = 0; k<NGP; k++) {

    645. 

  645.             ksi = GQpoint[k][0];

    646. 

  646.             eta = GQpoint[k][1];

    647. 

  647.             zeta = GQpoint[k][2];

    648. 

  648.             

    649. 

  649.             Sp[0][k] = 0.25*(1-ksi)*(1-eta)*(1-zeta);

    650. 

  650.             Sp[1][k] = 0.25*(1+ksi)*(1-eta)*(1-zeta);

    651. 

  651.             Sp[2][k] = 0.25*(1+ksi)*(1+eta)*(1-zeta);

    652. 

  652.             Sp[3][k] = 0.25*(1-ksi)*(1+eta)*(1-zeta);   

    653. 

  653.             Sp[4][k] = 0.25*(1-ksi)*(1-eta)*(1+zeta);

    654. 

  654.             Sp[5][k] = 0.25*(1+ksi)*(1-eta)*(1+zeta);

    655. 

  655.             Sp[6][k] = 0.25*(1+ksi)*(1+eta)*(1+zeta);

    656. 

  656.             Sp[7][k] = 0.25*(1-ksi)*(1+eta)*(1+zeta); 

    657. 

  657. 

    658. 

  658.             DSp[0][0][k] = -0.25*(1-eta)*(1-zeta);  // ksi derivative of the 1st shape funct. at k-th GQ point.

    659. 

  659.             DSp[1][0][k] = -0.25*(1-ksi)*(1-zeta);  // eta derivative of the 1st shape funct. at k-th GQ point.  

    660. 

  660.             DSp[2][0][k] = -0.25*(1-ksi)*(1-eta);   // zeta derivative of the 1st shape funct. at k-th GQ point.

    661. 

  661.             DSp[0][1][k] =  0.25*(1-eta)*(1-zeta);   

    662. 

  662.             DSp[1][1][k] = -0.25*(1+ksi)*(1-zeta);

    663. 

  663.             DSp[2][1][k] = -0.25*(1+ksi)*(1-eta);         

    664. 

  664.             DSp[0][2][k] =  0.25*(1+eta)*(1-zeta);   

    665. 

  665.             DSp[1][2][k] =  0.25*(1+ksi)*(1-zeta);

    666. 

  666.             DSp[2][2][k] = -0.25*(1+ksi)*(1+eta);         

    667. 

  667.             DSp[0][3][k] = -0.25*(1+eta)*(1-zeta);   

    668. 

  668.             DSp[1][3][k] =  0.25*(1-ksi)*(1-zeta);

    669. 

  669.             DSp[2][3][k] = -0.25*(1-ksi)*(1+eta); 

    670. 

  670.             DSp[0][4][k] = -0.25*(1-eta)*(1+zeta);  // ksi derivative of the 1st shape funct. at k-th GQ point.

    671. 

  671.             DSp[1][4][k] = -0.25*(1-ksi)*(1+zeta);  // eta derivative of the 1st shape funct. at k-th GQ point.  

    672. 

  672.             DSp[2][4][k] = 0.25*(1-ksi)*(1-eta);   // zeta derivative of the 1st shape funct. at k-th GQ point.

    673. 

  673.             DSp[0][5][k] =  0.25*(1-eta)*(1+zeta);   

    674. 

  674.             DSp[1][5][k] = -0.25*(1+ksi)*(1+zeta);

    675. 

  675.             DSp[2][5][k] = 0.25*(1+ksi)*(1-eta);         

    676. 

  676.             DSp[0][6][k] =  0.25*(1+eta)*(1+zeta);   

    677. 

  677.             DSp[1][6][k] =  0.25*(1+ksi)*(1+zeta);

    678. 

  678.             DSp[2][6][k] = 0.25*(1+ksi)*(1+eta);         

    679. 

  679.             DSp[0][7][k] = -0.25*(1+eta)*(1+zeta);   

    680. 

  680.             DSp[1][7][k] =  0.25*(1-ksi)*(1+zeta);

    681. 

  681.             DSp[2][7][k] = 0.25*(1-ksi)*(1+eta); 

    682. 

  682.          }

    683. 

  683.       }

    684. 

  684.       

    685. 

  685.       if (NENv == 8) {

    686. 

  686.       

    687. 

  687.          Sv = new double*[8];

    688. 

  688.          for (int i=0; i<NGP; i++) {

    689. 

  689.             Sv[i] = new double[NGP];

    690. 

  690.          }

    691. 

  691.       

    692. 

  692.          DSv = new double **[3];	

    693. 

  693.          for (int i=0; i<3; i++) {

    694. 

  694.             DSv[i] = new double *[8];

    695. 

  695.             for (int j=0; j<8; j++) {

    696. 

  696.                DSv[i][j] = new double[NGP];

    697. 

  697.             }

    698. 

  698.          }

    699. 

  699.          

    700. 

  700.          for (int i=0; i<NGP; i++) {

    701. 

  701.             Sv[i] = Sp[i];

    702. 

  702.          }

    703. 

  703.          

    704. 

  704.          for (int i=0; i<3; i++) {

    705. 

  705.             for (int j=0; j<8; j++) {

    706. 

  706.                DSv[i][j] = DSp[i][j];

    707. 

  707.             }

    708. 

  708.          }         

    709. 

  709.                

    710. 

  710.       }

    711. 

  711.    }   

    712. 

  712.       

    713. 

  713.    for (int i = 0; i<8; i++) {   //deleting the unnecessary arrays for future

    714. 

  714.       delete GQpoint[i];

    715. 

  715.    }

    716. 

  716.    delete [] GQpoint;

    717. 

  717.       

    718. 

  718. } // End of function calcShape()

    719. 

  719. 

    720. 

  720. 

    721. 

  721. 

    722. 

  722. 

    723. 

  723. //------------------------------------------------------------------------------

    724. 

  724. void calcJacobian()

    725. 

  725. //------------------------------------------------------------------------------

    726. 

  726. {

    727. 

  727.    // Calculates the Jacobian matrix, its inverse and determinant for all

    728. 

  728.    // elements at all GQ points. Also evaluates and stores derivatives of shape

    729. 

  729.    // functions wrt global coordinates x and y.

    730. 

  730. 

    731. 

  731.    int e, i, j, k, x, iG; 

    732. 

  732.    double **e_coord;

    733. 

  733.    double **Jacob, **invJacob;

    734. 

  734.    double temp;

    735. 

  735.    

    736. 

  736.    if (eType == 3) {

    737. 

  737.    

    738. 

  738.       e_coord = new double*[NENp];

    739. 

  739.    

    740. 

  740.       for (i=0; i<NENp; i++) {

    741. 

  741.          e_coord[i] = new double[3];

    742. 

  742.       }

    743. 

  743.    

    744. 

  744.       Jacob = new double*[3];

    745. 

  745.       invJacob = new double*[3];

    746. 

  746.       for (i=0; i<3; i++) {

    747. 

  747.          Jacob[i] = new double[3];

    748. 

  748.          invJacob[i] = new double[3];

    749. 

  749.       }

    750. 

  750.       

    751. 

  751.       detJacob = new double*[NE];

    752. 

  752.       for (i=0; i<NE; i++) {

    753. 

  753.          detJacob[i] = new double[NGP];

    754. 

  754.       }

    755. 

  755.       

    756. 

  756.       gDSp = new double***[NE];

    757. 

  757. 

    758. 

  758.       for (i=0; i<NE; i++) {

    759. 

  759.          gDSp[i] = new double**[3];

    760. 

  760.          for(j=0; j<3; j++) {

    761. 

  761.             gDSp[i][j] = new double*[NENv];

    762. 

  762.             for(k=0; k<NENp; k++) {

    763. 

  763.                gDSp[i][j][k] = new double[NGP];

    764. 

  764.             }

    765. 

  765.          }	

    766. 

  766.       }

    767. 

  767.       

    768. 

  768.       gDSv = new double***[NE];

    769. 

  769. 

    770. 

  770.       for (i=0; i<NE; i++) {

    771. 

  771.          gDSv[i] = new double**[3];

    772. 

  772.          for(j=0; j<3; j++) {

    773. 

  773.             gDSv[i][j] = new double*[NENv];

    774. 

  774.             for(k=0; k<NENv; k++) {

    775. 

  775.                gDSv[i][j][k] = new double[NGP];

    776. 

  776.             }

    777. 

  777.          }	

    778. 

  778.       }      

    779. 

  779.    

    780. 

  780.       for (e = 0; e<NE; e++){

    781. 

  781.          // To calculate Jacobian matrix of an element we need e_coord matrix of

    782. 

  782.          // size NEN*2. Each row of it stores x, y & z coordinates of the nodes of

    783. 

  783.          // an element.

    784. 

  784.          for (i = 0; i<NENp; i++){

    785. 

  785.             iG = LtoG[e][i];

    786. 

  786.             e_coord[i][0] = coord[iG][0]; 

    787. 

  787.             e_coord[i][1] = coord[iG][1];

    788. 

  788.             e_coord[i][2] = coord[iG][2];

    789. 

  789.          }

    790. 

  790.       

    791. 

  791.          // For each GQ point calculate 3*3 Jacobian matrix, its inverse and its

    792. 

  792.          // determinant. Also calculate derivatives of shape functions wrt global

    793. 

  793.          // coordinates x and y & z. These are the derivatives that we'll use in

    794. 

  794.          // evaluating K and F integrals. 

    795. 

  795.    	

    796. 

  796.          for (k = 0; k<NGP; k++) {

    797. 

  797.             for (i = 0; i<3; i++) {

    798. 

  798.                for (j = 0; j<3; j++) {

    799. 

  799.                   temp = 0;

    800. 

  800.                   for (x = 0; x<NENp; x++) {

    801. 

  801.                      temp = DSp[i][x][k] * e_coord[x][j]+temp;

    802. 

  802.                   }

    803. 

  803.                   Jacob[i][j] = temp;

    804. 

  804.                }

    805. 

  805.             }

    806. 

  806.             

    807. 

  807.             invJacob[0][0] = Jacob[1][1]*Jacob[2][2]-Jacob[2][1]*Jacob[1][2];

    808. 

  808.             invJacob[0][1] = -(Jacob[0][1]*Jacob[2][2]-Jacob[0][2]*Jacob[2][1]);

    809. 

  809.             invJacob[0][2] = Jacob[0][1]*Jacob[1][2]-Jacob[1][1]*Jacob[0][2];

    810. 

  810.             invJacob[1][0] =  -(Jacob[1][0]*Jacob[2][2]-Jacob[1][2]*Jacob[2][0]);

    811. 

  811.             invJacob[1][1] = Jacob[2][2]*Jacob[0][0]-Jacob[2][0]*Jacob[0][2];

    812. 

  812.             invJacob[1][2] = -(Jacob[1][2]*Jacob[0][0]-Jacob[1][0]*Jacob[0][2]);

    813. 

  813.             invJacob[2][0] = Jacob[1][0]*Jacob[2][1]-Jacob[2][0]*Jacob[1][1];

    814. 

  814.             invJacob[2][1] = -(Jacob[2][1]*Jacob[0][0]-Jacob[2][0]*Jacob[0][1]);

    815. 

  815.             invJacob[2][2] = Jacob[1][1]*Jacob[0][0]-Jacob[1][0]*Jacob[0][1];

    816. 

  816. 

    817. 

  817.             detJacob[e][k] = Jacob[0][0]*(Jacob[1][1]*Jacob[2][2]-Jacob[2][1]*Jacob[1][2]) +

    818. 

  818.                             Jacob[0][1]*(Jacob[1][2]*Jacob[2][0]-Jacob[1][0]*Jacob[2][2]) +

    819. 

  819.                             Jacob[0][2]*(Jacob[1][0]*Jacob[2][1]-Jacob[1][1]*Jacob[2][0]);

    820. 

  820.             

    821. 

  821.             for (i = 0; i<3; i++){

    822. 

  822.                for (j = 0; j<3; j++){

    823. 

  823.                   invJacob[i][j] = invJacob[i][j]/detJacob[e][k];

    824. 

  824.                }    

    825. 

  825.             }

    826. 

  826.          

    827. 

  827.             for (i = 0; i<3; i++){

    828. 

  828.                for (j = 0; j<8; j++){

    829. 

  829.                   temp = 0;

    830. 

  830.                   for (x = 0; x<3; x++){

    831. 

  831.                      temp = invJacob[i][x] * DSp[x][j][k]+temp;

    832. 

  832.                   }

    833. 

  833.                   gDSp[e][i][j][k] = temp;

    834. 

  834.                   gDSv[e][i][j][k] = gDSp[e][i][j][k];

    835. 

  835.                }

    836. 

  836.             }

    837. 

  837.          }

    838. 

  838.       }

    839. 

  839.    } 

    840. 

  840.    

    841. 

  841. }  // End of function calcJacobian()

    842. 

  842. 

    843. 

  843. 

    844. 

  844. 

    845. 

  845. 

    846. 

  846. //------------------------------------------------------------------------------

    847. 

  847. void initGlobalSysVariables()

    848. 

  848. //------------------------------------------------------------------------------

    849. 

  849. {

    850. 

  850.    // Allocating the memory for stiffness matrices and force vectors

    851. 

  851.    

    852. 

  852.    int i;

    853. 

  853.    

    854. 

  854.    F = new real[Ndof];

    855. 

  855.    K = new real*[Ndof];	

    856. 

  856.    for (i=0; i<Ndof; i++) {

    857. 

  857.       K[i] = new real[Ndof];

    858. 

  858.    }

    859. 

  859.    

    860. 

  860.    Fe = new double[NEU];

    861. 

  861.    Ke = new double*[NEU];	

    862. 

  862.    for (i=0; i<NEU; i++) {

    863. 

  863.       Ke[i] = new double[NEU];

    864. 

  864.    }

    865. 

  865.    

    866. 

  866.    Fe_1 = new double[NENv];

    867. 

  867.    Fe_2 = new double[NENv];

    868. 

  868.    Fe_3 = new double[NENv];

    869. 

  869.    Fe_4 = new double[NENp];

    870. 

  870.    

    871. 

  871.    Ke_11 = new double*[NENv];

    872. 

  872.    Ke_12 = new double*[NENv];

    873. 

  873. 	Ke_13 = new double*[NENv];

    874. 

  874.    Ke_14 = new double*[NENv];

    875. 

  875.    // Ke_21 is the transpose of Ke_12

    876. 

  876.    Ke_22 = new double*[NENv];

    877. 

  877.    Ke_23 = new double*[NENv];

    878. 

  878. 	Ke_24 = new double*[NENv]; 

    879. 

  879.    // Ke_31 is the transpose of Ke_13

    880. 

  880.    // Ke_32 is the transpose of Ke_23

    881. 

  881. 	Ke_33 = new double*[NENv];

    882. 

  882. 	Ke_34 = new double*[NENv];

    883. 

  883.    Ke_41 = new double*[NENv];

    884. 

  884.    Ke_42 = new double*[NENv];

    885. 

  885. 	Ke_43 = new double*[NENv];

    886. 

  886. 	Ke_44 = new double*[NENp];

    887. 

  887.    

    888. 

  888.    Ke_11_add = new double*[NENv];

    889. 

  889.    Ke_12_add = new double*[NENv];

    890. 

  890.    Ke_13_add = new double*[NENv];

    891. 

  891.    Ke_14_add = new double*[NENv];

    892. 

  892.    Ke_22_add = new double*[NENv];

    893. 

  893.    Ke_23_add = new double*[NENv];

    894. 

  894.    Ke_24_add = new double*[NENv];

    895. 

  895.    Ke_33_add = new double*[NENv];

    896. 

  896.    Ke_34_add = new double*[NENv];

    897. 

  897.    Ke_41_add = new double*[NENp];

    898. 

  898.    Ke_42_add = new double*[NENp];

    899. 

  899.    Ke_43_add = new double*[NENp];

    900. 

  900.    Ke_44_add = new double*[NENp];

    901. 

  901.    for (i=0; i<NENv; i++) {

    902. 

  902.       Ke_11[i] = new double[NENv];

    903. 

  903.       Ke_12[i] = new double[NENv];

    904. 

  904.       Ke_13[i] = new double[NENv];

    905. 

  905.       Ke_14[i] = new double[NENp];

    906. 

  906.       Ke_22[i] = new double[NENv];

    907. 

  907.       Ke_23[i] = new double[NENv];

    908. 

  908.       Ke_24[i] = new double[NENp];

    909. 

  909.       Ke_33[i] = new double[NENv];

    910. 

  910.       Ke_34[i] = new double[NENp];

    911. 

  911.       

    912. 

  912.       Ke_11_add[i] = new double[NENv];

    913. 

  913.       Ke_12_add[i] = new double[NENv];

    914. 

  914.       Ke_13_add[i] = new double[NENv];

    915. 

  915.       Ke_14_add[i] = new double[NENp];

    916. 

  916.       Ke_22_add[i] = new double[NENv];

    917. 

  917.       Ke_23_add[i] = new double[NENv];

    918. 

  918.       Ke_24_add[i] = new double[NENp];

    919. 

  919.       Ke_33_add[i] = new double[NENv];

    920. 

  920.       Ke_34_add[i] = new double[NENp];

    921. 

  921.    }

    922. 

  922.    for (i=0; i<NENp; i++) {

    923. 

  923.       Ke_41[i] = new double[NENv];

    924. 

  924.       Ke_42[i] = new double[NENv];

    925. 

  925. 	   Ke_43[i] = new double[NENv];

    926. 

  926.       Ke_44[i] = new double[NENp]; 

    927. 

  927.       Ke_41_add[i] = new double[NENp];

    928. 

  928.       Ke_42_add[i] = new double[NENp];

    929. 

  929.       Ke_43_add[i] = new double[NENp];

    930. 

  930.       Ke_44_add[i] = new double[NENp];

    931. 

  931.    }   

    932. 

  932.    

    933. 

  933.    uNodal = new double[NENv];

    934. 

  934.    vNodal = new double[NENv];

    935. 

  935.    wNodal = new double[NENv];

    936. 

  936.    

    937. 

  937. }   

    938. 

  938. 

    939. 

  939. 

    940. 

  940. 

    941. 

  941. 

    942. 

  942. //------------------------------------------------------------------------------

    943. 

  943. void calcGlobalSys()

    944. 

  944. //------------------------------------------------------------------------------

    945. 

  945. {

    946. 

  946.    // Calculates Ke and Fe one by one for each element and assembles them into

    947. 

  947.    // the global K and F.

    948. 

  948. 

    949. 

  949.    int e, i, j, k, m, n, node;

    950. 

  950.    double Tau;

    951. 

  951.    //   double x, y, axy, fxy;

    952. 

  952. 

    953. 

  953.    // Initialize the arrays

    954. 

  954. 

    955. 

  955.    for (i=0; i<Ndof; i++) {

    956. 

  956.       F[i] = 0;

    957. 

  957.       for (j=0; j<Ndof; j++) {

    958. 

  958.          K[i][j] = 0;

    959. 

  959. 	   }

    960. 

  960.    }

    961. 

  961.    

    962. 

  962.    // Calculating the elemental stiffness matrix(Ke) and force vector(Fe)

    963. 

  963. 

    964. 

  964.    for (e = 0; e<NE; e++) {

    965. 

  965.       // Intitialize Ke and Fe to zero.

    966. 

  966.       for (i=0; i<NEU; i++) {

    967. 

  967.          Fe[i] = 0;

    968. 

  968.          for (j=0; j<NEU; j++) {

    969. 

  969.             Ke[i][j] = 0;

    970. 

  970.          }

    971. 

  971.       }

    972. 

  972. 

    973. 

  973.       for (i=0; i<NENv; i++) {

    974. 

  974.          Fe_1[i] = 0;

    975. 

  975.          Fe_2[i] = 0;

    976. 

  976.          Fe_3[i] = 0;

    977. 

  977.          Fe_4[i] = 0;

    978. 

  978.          for (j=0; j<NENv; j++) {

    979. 

  979.             Ke_11[i][j] = 0;

    980. 

  980.             Ke_12[i][j] = 0;

    981. 

  981.             Ke_13[i][j] = 0;

    982. 

  982.             Ke_22[i][j] = 0;

    983. 

  983.             Ke_23[i][j] = 0;

    984. 

  984.             Ke_33[i][j] = 0;

    985. 

  985.          }

    986. 

  986.          for (j=0; j<NENp; j++) {

    987. 

  987.             Ke_14[i][j] = 0;

    988. 

  988.             Ke_24[i][j] = 0;

    989. 

  989.             Ke_34[i][j] = 0;

    990. 

  990.          }         

    991. 

  991.       }   

    992. 

  992.       for (i=0; i<NENp; i++) {

    993. 

  993.          for (j=0; j<NENp; j++) {

    994. 

  994.             Ke_41[i][j] = 0;

    995. 

  995.             Ke_42[i][j] = 0;

    996. 

  996.             Ke_43[i][j] = 0;

    997. 

  997.             Ke_44[i][j] = 0;

    998. 

  998.          }

    999. 

  999.       }

    1000. 

  1000.       

    1001. 

  1001.       for (i=0; i<NENv; i++) {

    1002. 

  1002.          uNodal[i] = u[LtoG[e][i]]; 

    1003. 

  1003.          vNodal[i] = u[LtoG[e][i]+NN];

    1004. 

  1004.          wNodal[i] = u[LtoG[e][i]+2*NN];

    1005. 

  1005.       }

    1006. 

  1006.       

    1007. 

  1007.       for (k = 0; k<NGP; k++) {   // Gauss quadrature loop

    1008. 

  1008.             

    1009. 

  1009.          u0 = 0;

    1010. 

  1010.          v0 = 0;

    1011. 

  1011.          w0 = 0;

    1012. 

  1012.          

    1013. 

  1013.          for (i=0; i<NENv; i++) {

    1014. 

  1014.             u0 = u0 + Sp[i][k] * uNodal[i];

    1015. 

  1015.             v0 = v0 + Sp[i][k] * vNodal[i];

    1016. 

  1016.             w0 = w0 + Sp[i][k] * wNodal[i];            

    1017. 

  1017.          }

    1018. 

  1018.          

    1019. 

  1019.          for (i=0; i<NENv; i++) {

    1020. 

  1020.             for (j=0; j<NENv; j++) {

    1021. 

  1021.                Ke_11_add[i][j] = 0;

    1022. 

  1022.                Ke_12_add[i][j] = 0;

    1023. 

  1023.                Ke_13_add[i][j] = 0;

    1024. 

  1024.                Ke_22_add[i][j] = 0;

    1025. 

  1025.                Ke_23_add[i][j] = 0;

    1026. 

  1026.                Ke_33_add[i][j] = 0;

    1027. 

  1027.             }

    1028. 

  1028.             for (j=0; j<NENp; j++) {

    1029. 

  1029.                Ke_14_add[i][j] = 0;

    1030. 

  1030.                Ke_24_add[i][j] = 0;

    1031. 

  1031.                Ke_34_add[i][j] = 0;

    1032. 

  1032.             }         

    1033. 

  1033.          }

    1034. 

  1034.          for (i=0; i<NENp; i++) {

    1035. 

  1035.             for (j=0; j<NENp; j++) {

    1036. 

  1036.                Ke_41_add[i][j] = 0;

    1037. 

  1037.                Ke_42_add[i][j] = 0;

    1038. 

  1038.                Ke_43_add[i][j] = 0;

    1039. 

  1039.                Ke_44_add[i][j] = 0;

    1040. 

  1040.             }

    1041. 

  1041.          }

    1042. 

  1042. 

    1043. 

  1043.          for (i=0; i<NENv; i++) {

    1044. 

  1044.             for (j=0; j<NENv; j++) {

    1045. 

  1045.                Ke_11_add[i][j] = Ke_11_add[i][j] + viscosity * (

    1046. 

  1046.                      2 * gDSv[e][0][i][k] * gDSv[e][0][j][k] + 

    1047. 

  1047.                      gDSv[e][1][i][k] * gDSv[e][1][j][k] +

    1048. 

  1048.                      gDSv[e][2][i][k] * gDSv[e][2][j][k]) +

    1049. 

  1049.                      density * Sv[i][k] * (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);

    1050. 

  1050.                      //density * Sv[i][k];

    1051. 

  1051.                      

    1052. 

  1052.                Ke_12_add[i][j] = Ke_12_add[i][j] + viscosity * gDSv[e][1][i][k] * gDSv[e][0][j][k];

    1053. 

  1053.                

    1054. 

  1054.                Ke_13_add[i][j] = Ke_13_add[i][j] + viscosity * gDSv[e][2][i][k] * gDSv[e][0][j][k];

    1055. 

  1055.                

    1056. 

  1056.                Ke_22_add[i][j] = Ke_22_add[i][j] + viscosity * (

    1057. 

  1057.                      gDSv[e][0][i][k] * gDSv[e][0][j][k] + 

    1058. 

  1058.                      2 * gDSv[e][1][i][k] * gDSv[e][1][j][k] +

    1059. 

  1059.                      gDSv[e][2][i][k] * gDSv[e][2][j][k]) +

    1060. 

  1060.                      density * Sv[i][k] * (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);

    1061. 

  1061.                      

    1062. 

  1062.                Ke_23_add[i][j] = Ke_23_add[i][j] + viscosity * gDSv[e][2][i][k] * gDSv[e][1][j][k];

    1063. 

  1063.                

    1064. 

  1064.                Ke_33_add[i][j] = Ke_33_add[i][j] + viscosity * (

    1065. 

  1065.                      gDSv[e][0][i][k] * gDSv[e][0][j][k] + 

    1066. 

  1066.                      gDSv[e][1][i][k] * gDSv[e][1][j][k] +

    1067. 

  1067.                      2 * gDSv[e][2][i][k] * gDSv[e][2][j][k]) +

    1068. 

  1068.                      density * Sv[i][k] * (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);

    1069. 

  1069.                      

    1070. 

  1070.             }

    1071. 

  1071.             for (j=0; j<NENp; j++) {

    1072. 

  1072.                Ke_14_add[i][j] = Ke_14_add[i][j] - gDSv[e][0][i][k] * Sp[j][k];

    1073. 

  1073.                Ke_24_add[i][j] = Ke_24_add[i][j] - gDSv[e][1][i][k] * Sp[j][k];

    1074. 

  1074.                Ke_34_add[i][j] = Ke_34_add[i][j] - gDSv[e][2][i][k] * Sp[j][k];

    1075. 

  1075.             }         

    1076. 

  1076.          }

    1077. 

  1077.          

    1078. 

  1078.          for (i=0; i<NENv; i++) {

    1079. 

  1079.             for (j=0; j<NENp; j++) {

    1080. 

  1080.                Ke_41_add[j][i] = Ke_14_add[i][j] ;

    1081. 

  1081.                Ke_42_add[j][i] = Ke_24_add[i][j] ;

    1082. 

  1082.                Ke_43_add[j][i] = Ke_34_add[i][j] ;

    1083. 

  1083.             }

    1084. 

  1084.          }

    1085. 

  1085. 

    1086. 

  1086.   // Apply GLS stabilization for linear elements with NENv = NENp

    1087. 

  1087.     

    1088. 

  1088.          Tau = ((1.0/12.0)*2) / viscosity;  // GLS parameter

    1089. 

  1089. 

    1090. 

  1090.          for (i=0; i<NENv; i++) {

    1091. 

  1091.             for (j=0; j<NENv; j++) {

    1092. 

  1092.               Ke_11_add[i][j] = Ke_11_add[i][j] + Tau * density * density *

    1093. 

  1093.                              (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) *

    1094. 

  1094.                              (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);                             

    1095. 

  1095.              

    1096. 

  1096.               Ke_22_add[i][j] = Ke_22_add[i][j] + Tau * density * density * 

    1097. 

  1097.                              (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) *

    1098. 

  1098.                              (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);   

    1099. 

  1099.                              

    1100. 

  1100.               Ke_33_add[i][j] = Ke_33_add[i][j] + Tau * density * density * 

    1101. 

  1101.                              (u0 * gDSv[e][0][i][k] + v0 * gDSv[e][1][i][k] + w0 * gDSv[e][2][i][k]) *

    1102. 

  1102.                              (u0 * gDSv[e][0][j][k] + v0 * gDSv[e][1][j][k] + w0 * gDSv[e][2][j][k]);                                

    1103. 

  1103.             } 

    1104. 

  1104.          } 

    1105. 

  1105. 

    1106. 

  1106.          for (i=0; i<NENv; i++) {

    1107. 

  1107.             for (j=0; j<NENp; j++) {

    1108. 

  1108.               Ke_14_add[i][j] = Ke_14_add[i][j] + Tau * density * 

    1109. 

  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. 

  1110.              

    1111. 

  1111.               Ke_24_add[i][j] = Ke_24_add[i][j] + Tau * density * 

    1112. 

  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. 

  1113.                              

    1114. 

  1114.               Ke_34_add[i][j] = Ke_34_add[i][j] + Tau * density * 

    1115. 

  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. 

  1116.             } 

    1117. 

  1117.          } 

    1118. 

  1118. 

    1119. 

  1119.          for (i=0; i<NENp; i++) {

    1120. 

  1120.             for (j=0; j<NENv; j++) {

    1121. 

  1121.               Ke_41_add[i][j] = Ke_41_add[i][j] - Tau * density * 

    1122. 

  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. 

  1123.              

    1124. 

  1124.               Ke_42_add[i][j] = Ke_42_add[i][j] - Tau * density * 

    1125. 

  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. 

  1126.                              

    1127. 

  1127.               Ke_43_add[i][j] = Ke_43_add[i][j] - Tau * density * 

    1128. 

  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. 

  1129.             } 

    1130. 

  1130.          } 

    1131. 

  1131.          

    1132. 

  1132.          for (i=0; i<NENp; i++) {

    1133. 

  1133.             for (j=0; j<NENp; j++) {

    1134. 

  1134.                Ke_44_add[i][j] = Ke_44_add[i][j] - Tau *

    1135. 

  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. 

  1136.             }

    1137. 

  1137.          }

    1138. 

  1138. 

    1139. 

  1139.   // Apply GLS stabilization for linear elements with NENv = NENp

    1140. 

  1140.             

    1141. 

  1141.          for (i=0; i<NENv; i++) {

    1142. 

  1142.             for (j=0; j<NENv; j++) {            

    1143. 

  1143.                Ke_11[i][j] += Ke_11_add[i][j] * detJacob[e][k] * GQweight[k];

    1144. 

  1144.                Ke_12[i][j] += Ke_12_add[i][j] * detJacob[e][k] * GQweight[k];

    1145. 

  1145.                Ke_13[i][j] += Ke_13_add[i][j] * detJacob[e][k] * GQweight[k];

    1146. 

  1146.                Ke_22[i][j] += Ke_22_add[i][j] * detJacob[e][k] * GQweight[k];

    1147. 

  1147.                Ke_23[i][j] += Ke_23_add[i][j] * detJacob[e][k] * GQweight[k];

    1148. 

  1148.                Ke_33[i][j] += Ke_33_add[i][j] * detJacob[e][k] * GQweight[k];

    1149. 

  1149.             }

    1150. 

  1150.             for (j=0; j<NENp; j++) {               

    1151. 

  1151.                Ke_14[i][j] += Ke_14_add[i][j] * detJacob[e][k] * GQweight[k];

    1152. 

  1152.                Ke_24[i][j] += Ke_24_add[i][j] * detJacob[e][k] * GQweight[k];

    1153. 

  1153.                Ke_34[i][j] += Ke_34_add[i][j] * detJacob[e][k] * GQweight[k];

    1154. 

  1154.             }

    1155. 

  1155.          }

    1156. 

  1156.          for (i=0; i<NENv; i++) {

    1157. 

  1157.             for (j=0; j<NENv; j++) {

    1158. 

  1158.                Ke_41[i][j] += Ke_41_add[i][j] * detJacob[e][k] * GQweight[k];

    1159. 

  1159.                Ke_42[i][j] += Ke_42_add[i][j] * detJacob[e][k] * GQweight[k];

    1160. 

  1160.                Ke_43[i][j] += Ke_43_add[i][j] * detJacob[e][k] * GQweight[k];

    1161. 

  1161.                Ke_44[i][j] += Ke_44_add[i][j] * detJacob[e][k] * GQweight[k];

    1162. 

  1162.             }

    1163. 

  1163.          }

    1164. 

  1164. 

    1165. 

  1165.       }   // End GQ loop  

    1166. 

  1166. 

    1167. 

  1167.       // Assembly of Fe

    1168. 

  1168. 

    1169. 

  1169.       i=0;

    1170. 

  1170.       for (j=0; j<NENv; j++) {

    1171. 

  1171.          Fe[i]=Fe_1[j];

    1172. 

  1172.          i++;

    1173. 

  1173.       }

    1174. 

  1174.       for (j=0; j<NENv; j++) {

    1175. 

  1175.          Fe[i]=Fe_2[j];

    1176. 

  1176.          i++;

    1177. 

  1177.       }

    1178. 

  1178.       for (j=0; j<NENv; j++) {

    1179. 

  1179.          Fe[i]=Fe_3[j];

    1180. 

  1180.          i++;

    1181. 

  1181.       }

    1182. 

  1182.       for (j=0; j<NENp; j++) {

    1183. 

  1183.          Fe[i]=Fe_4[j];

    1184. 

  1184.          i++;

    1185. 

  1185.       }         

    1186. 

  1186. 

    1187. 

  1187.       // Assembly of Ke

    1188. 

  1188. 

    1189. 

  1189.       i=0;

    1190. 

  1190.       j=0;    

    1191. 

  1191.       for (m=0; m<NENv; m++) {

    1192. 

  1192.          j=0;

    1193. 

  1193.          for (n=0; n<NENv; n++) {       

    1194. 

  1194.             Ke[i][j] =  Ke_11[m][n];

    1195. 

  1195.             j++;

    1196. 

  1196.          }

    1197. 

  1197.          i++;

    1198. 

  1198.       }

    1199. 

  1199.       

    1200. 

  1200.       i=i-NENv;

    1201. 

  1201.       for (m=0; m<NENv; m++) {

    1202. 

  1202.          j=0;

    1203. 

  1203.          for (n=0; n<NENv; n++) {       

    1204. 

  1204.             Ke[i][j+NENv] =  Ke_12[m][n];

    1205. 

  1205.             j++;

    1206. 

  1206.          }

    1207. 

  1207.          i++;

    1208. 

  1208.       }

    1209. 

  1209. 

    1210. 

  1210.       i=i-NENv;

    1211. 

  1211.       for (m=0; m<NENv; m++) {

    1212. 

  1212.          j=0;

    1213. 

  1213.          for (n=0; n<NENv; n++) {       

    1214. 

  1214.             Ke[i][j+2*NENv] =  Ke_13[m][n];

    1215. 

  1215.             j++;

    1216. 

  1216.          }

    1217. 

  1217.          i++;

    1218. 

  1218.       }

    1219. 

  1219.       

    1220. 

  1220.       i=i-NENv;

    1221. 

  1221.       for (m=0; m<NENv; m++) {

    1222. 

  1222.          j=0;

    1223. 

  1223.          for (n=0; n<NENp; n++) {       

    1224. 

  1224.             Ke[i][j+3*NENv] =  Ke_14[m][n];

    1225. 

  1225.             j++;

    1226. 

  1226.          }

    1227. 

  1227.          i++;

    1228. 

  1228.       }

    1229. 

  1229.       

    1230. 

  1230.       for (m=0; m<NENv; m++) {

    1231. 

  1231.          j=0;

    1232. 

  1232.          for (n=0; n<NENv; n++) {  

    1233. 

  1233.             Ke[i][j] =  Ke_12[n][m];

    1234. 

  1234.             j++;

    1235. 

  1235.          }

    1236. 

  1236.          i++;

    1237. 

  1237.       }

    1238. 

  1238. 

    1239. 

  1239.       i=i-NENv;

    1240. 

  1240.       for (m=0; m<NENv; m++) {

    1241. 

  1241.          j=0;

    1242. 

  1242.          for (n=0; n<NENv; n++) {       

    1243. 

  1243.             Ke[i][j+NENv] =  Ke_22[m][n];

    1244. 

  1244.             j++;

    1245. 

  1245.          }

    1246. 

  1246.          i++;

    1247. 

  1247.       }

    1248. 

  1248. 

    1249. 

  1249.       i=i-NENv;

    1250. 

  1250.       for (m=0; m<NENv; m++) {

    1251. 

  1251.          j=0;

    1252. 

  1252.          for (n=0; n<NENv; n++) {       

    1253. 

  1253.             Ke[i][j+2*NENv] =  Ke_23[m][n];

    1254. 

  1254.             j++;

    1255. 

  1255.          }

    1256. 

  1256.          i++;

    1257. 

  1257.       }

    1258. 

  1258. 

    1259. 

  1259.       i=i-NENv;

    1260. 

  1260.       for (m=0; m<NENv; m++) {

    1261. 

  1261.          j=0;

    1262. 

  1262.          for (n=0; n<NENp; n++) {       

    1263. 

  1263.             Ke[i][j+3*NENv] =  Ke_24[m][n];

    1264. 

  1264.             j++;

    1265. 

  1265.          }

    1266. 

  1266.          i++;

    1267. 

  1267.       }

    1268. 

  1268. 

    1269. 

  1269.       for (m=0; m<NENv; m++) {

    1270. 

  1270.          j=0;

    1271. 

  1271.          for (n=0; n<NENv; n++) {       

    1272. 

  1272.             Ke[i][j] =  Ke_13[n][m];

    1273. 

  1273.             j++;

    1274. 

  1274.          }

    1275. 

  1275.          i++;

    1276. 

  1276.       }      

    1277. 

  1277. 

    1278. 

  1278.       i=i-NENv;

    1279. 

  1279.       for (m=0; m<NENv; m++) {

    1280. 

  1280.          j=0;

    1281. 

  1281.          for (n=0; n<NENv; n++) {       

    1282. 

  1282.             Ke[i][j+NENv] =  Ke_23[n][m];

    1283. 

  1283.             j++;

    1284. 

  1284.          }

    1285. 

  1285.          i++;

    1286. 

  1286.       } 

    1287. 

  1287. 

    1288. 

  1288.       i=i-NENv;

    1289. 

  1289.       for (m=0; m<NENv; m++) {

    1290. 

  1290.          j=0;

    1291. 

  1291.          for (n=0; n<NENv; n++) {       

    1292. 

  1292.             Ke[i][j+2*NENv] =  Ke_33[m][n];

    1293. 

  1293.             j++;

    1294. 

  1294.          }

    1295. 

  1295.          i++;

    1296. 

  1296.       }

    1297. 

  1297. 

    1298. 

  1298.       i=i-NENv;

    1299. 

  1299.       for (m=0; m<NENv; m++) {

    1300. 

  1300.          j=0;

    1301. 

  1301.          for (n=0; n<NENp; n++) {       

    1302. 

  1302.             Ke[i][j+3*NENv] =  Ke_34[m][n];

    1303. 

  1303.             j++;

    1304. 

  1304.          }

    1305. 

  1305.          i++;

    1306. 

  1306.       }

    1307. 

  1307. 

    1308. 

  1308.       j=0;

    1309. 

  1309.       for (m=0; m<NENv; m++) {

    1310. 

  1310.          j=0;

    1311. 

  1311.          for (n=0; n<NENp; n++) {       

    1312. 

  1312.             Ke[i][j] =  Ke_41[m][n];

    1313. 

  1313.             j++;

    1314. 

  1314.          }

    1315. 

  1315.          i++;

    1316. 

  1316.       }

    1317. 

  1317.       

    1318. 

  1318.       i=i-NENv;

    1319. 

  1319.       for (m=0; m<NENv; m++) {

    1320. 

  1320.          j=0;

    1321. 

  1321.          for (n=0; n<NENp; n++) {       

    1322. 

  1322.             Ke[i][j+NENv] =  Ke_42[m][n];

    1323. 

  1323.             j++;

    1324. 

  1324.          }

    1325. 

  1325.          i++;

    1326. 

  1326.       } 

    1327. 

  1327. 

    1328. 

  1328.       i=i-NENv;

    1329. 

  1329.       for (m=0; m<NENv; m++) {

    1330. 

  1330.          j=0;

    1331. 

  1331.          for (n=0; n<NENp; n++) {       

    1332. 

  1332.             Ke[i][j+2*NENv] =  Ke_43[m][n];

    1333. 

  1333.             j++;

    1334. 

  1334.          }

    1335. 

  1335.          i++;

    1336. 

  1336.       } 

    1337. 

  1337. 

    1338. 

  1338.       i=i-NENv;

    1339. 

  1339.       for (m=0; m<NENp; m++) {

    1340. 

  1340.          j=0;

    1341. 

  1341.          for (n=0; n<NENp; n++) {       

    1342. 

  1342.             Ke[i][j+3*NENv] = Ke_44[m][n];

    1343. 

  1343.             j++;

    1344. 

  1344.          }

    1345. 

  1345.          i++;

    1346. 

  1346.       }

    1347. 

  1347. 

    1348. 

  1348.       assemble(e, Ke, Fe);  // sending Ke & Fe for assembly

    1349. 

  1349. 

    1350. 

  1350. 

    1351. 

  1351.    }   // End element loop

    1352. 

  1352.  

    1353. 

  1353. } // End of function calcGlobalSys()

    1354. 

  1354. 

    1355. 

  1355. 

    1356. 

  1356. 

    1357. 

  1357. 

    1358. 

  1358. //------------------------------------------------------------------------------

    1359. 

  1359. void assemble(int e, double **Ke, double *Fe)

    1360. 

  1360. //------------------------------------------------------------------------------

    1361. 

  1361. {

    1362. 

  1362.    // Inserts Ke and Fe into proper locations of K and F.

    1363. 

  1363. 

    1364. 

  1364.    int i, j, m, n, iG, jG; 

    1365. 

  1365.    

    1366. 

  1366.    // Inserts [Ke] and {Fe} into proper locations of [K] and {F} by the help of LtoG array. 

    1367. 

  1367.    

    1368. 

  1368.    for (i = 0; i<NEU; i++) {

    1369. 

  1369.          

    1370. 

  1370.       iG = LtoG[e][i%NENv] + NN*(i/NENv);

    1371. 

  1371.       F[iG] = F[iG] + Fe[i];

    1372. 

  1372.       

    1373. 

  1373.       for (j = 0; j<NEU; j++) {   

    1374. 

  1374.          jG = LtoG[e][j%NENv] + NN*(j/NENv);

    1375. 

  1375.          K[iG][jG] = K[iG][jG] + Ke[i][j];

    1376. 

  1376.       }

    1377. 

  1377.    }   

    1378. 

  1378. 

    1379. 

  1379.    // Assembly process for compressed sparse storage 

    1380. 

  1380.    // KeKMapSmall creation

    1381. 

  1381. 

    1382. 

  1382.    int shiftRow, shiftCol;

    1383. 

  1383.    int *nodeData, p, q, k;

    1384. 

  1384.    // int *nodeDataEBC;

    1385. 

  1385.     

    1386. 

  1386.    // nodeDataEBC = new int[NENv];     // Elemental node data(LtoG data) (modified with EBC)    // BC implementationu sonraya býrak ineff ama bir anda zor

    1387. 

  1387.    nodeData = new int[NENv];           // Stores sorted LtoG data

    1388. 

  1388.    for(k=0; k<NENv; k++) {

    1389. 

  1389.       nodeData[k] = (LtoG[e][k]);      // Takes node data from LtoG

    1390. 

  1390.    } 

    1391. 

  1391. 

    1392. 

  1392.    // for(i=0; i<NEN; i++) {

    1393. 

  1393.       // nodeData[i]=nodeDataEBC[i];

    1394. 

  1394.       // if (GtoL[nodeDataEBC[i]][0] == -1) {

    1395. 

  1395.          // val[rowStarts[nodeDataEBC[i]]] = 1*bigNumber;      // Must be fixed, val[x] = 1 repeating for every element that contains the node!

    1396. 

  1396.          // nodeDataEBC[i] = -1;

    1397. 

  1397.       // }

    1398. 

  1398.    // }

    1399. 

  1399. 

    1400. 

  1400.    KeKMapSmall = new int*[NENv];        // NENv X NENv                                          

    1401. 

  1401.    for(j=0; j<NENv; j++) {              // Stores map data between K elemental and value vector

    1402. 

  1402.       KeKMapSmall[j] = new int[NENv];

    1403. 

  1403.    }

    1404. 

  1404. 

    1405. 

  1405.    for(i=0; i<NENv; i++) {

    1406. 

  1406.       for(j=0; j<NENv; j++) {

    1407. 

  1407.          q=0;

    1408. 

  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. 

  1409.             if(colSmall[p] == nodeData[j]) {                                         // Selection process of the KeKMapSmall data from the col vector

    1410. 

  1410.                KeKMapSmall[i][j] = q; 

    1411. 

  1411.                break;

    1412. 

  1412.             }

    1413. 

  1413.             q++;

    1414. 

  1414.          }

    1415. 

  1415.       }

    1416. 

  1416.    }

    1417. 

  1417. 

    1418. 

  1418.    // Creating val vector

    1419. 

  1419.    for(shiftRow=1; shiftRow<5; shiftRow++) {

    1420. 

  1420.       for(shiftCol=1; shiftCol<5; shift…
    1421.

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