API tools faq
paste
Advertisement
amber91

Untitled

amber91
Apr 9th, 2017
149
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 9.69 KB | None | 0 0
raw download clone embed print report
  1. dane do txt
  2. 12  100 1200    25  300 100
  3. 0   5   6   1   0 0 1 0 1 1 0 1 0 0
  4. 1   6   7   2   0 1 1 1 1 2 0 2 0 0
  5. 2   7   8   3   0 2 1 2 1 3 0 3 0 0
  6. 3   8   9   4   0 3 1 3 1 4 0 4 1 2
  7. 5   10  11  6   1 0 2 0 2 1 1 1 0 0
  8. 6   11  12  7   1 1 2 1 2 2 1 2 0 0
  9. 7   12  13  8   1 2 2 2 2 3 1 3 0 0
  10. 8   13  14  9   1 3 2 3 2 4 1 4 1 2
  11. 10  15  16  11  2 0 3 0 3 1 2 1 0 0
  12. 11  16  17  12  2 1 3 1 3 2 2 2 0 0
  13. 15  18  19  16  3 0 4 0 4 1 3 1 2 1
  14. 16  19  20  17  3 1 4 1 4 2 3 2 2 1
  15.  
  16.  
  17.  
  18.  
  19. #include <iostream>
  20. #include <iomanip>
  21. #include <math.h>
  22. #include <fstream>
  23.  
  24. using namespace std;
  25.  
  26. const double eps = 1e-12;   // stala przyblizenia zera potrzebna do obliczenia ukladu rownan metoda Gaussa
  27. int licznik = 0;
  28. float t0, tinf, k, alfa, q;
  29. float ksi[4][2] = {
  30.     { -0.5773502692, -0.5773502692 },
  31.     { 0.5773502692, -0.5773502692 },
  32.     { 0.5773502692, 0.5773502692 },
  33.     { -0.5773502692, 0.5773502692 }
  34. };
  35.  
  36.  
  37. float N(int no, float ksi, float ni){
  38.     switch (no)
  39.     {
  40.     case 1:
  41.         return (1 - ksi)*(1 - ni) / 4;
  42.     case 2:
  43.         return (1 + ksi)*(1 - ni) / 4;
  44.     case 3:
  45.         return (1 + ksi)*(1 + ni) / 4;
  46.     case 4:
  47.         return (1 - ksi)*(1 + ni) / 4;
  48.     default:
  49.         break;
  50.     }
  51.  
  52. }
  53.  
  54. float NKsi(int no, float ni){
  55.     switch (no)
  56.     {
  57.     case 1:
  58.         return (-1)*(1 - ni) / 4;
  59.     case 2:
  60.         return (1 - ni) / 4;
  61.     case 3:
  62.         return (1 + ni) / 4;
  63.     case 4:
  64.         return (-1)*(1 + ni) / 4;
  65.     default:
  66.         break;
  67.     }
  68.  
  69. }
  70. float Nni(int no, float ksi){
  71.     switch (no)
  72.     {
  73.     case 1:
  74.         return (1 - ksi)*(-1) / 4;
  75.     case 2:
  76.         return (1 + ksi)*(-1) / 4;
  77.     case 3:
  78.         return (1 + ksi) / 4;
  79.     case 4:
  80.         return (1 - ksi) / 4;
  81.     default:
  82.         break;
  83.     }
  84.  
  85. }
  86. float Rp(float ksi, float ni, float x1, float x2, float x3, float x4){  //Funkcja przejscia elementu lokalnego do globalnego
  87.     return N(1, ksi, ni)*x1 + N(2, ksi, ni)*x2 + N(3, ksi, ni)*x3 + N(4, ksi, ni)*x4;
  88.  
  89. }
  90.  
  91. float J00(float x1, float x2, float x3, float x4, float ni){
  92.     return (NKsi(1, ni)*x1 + NKsi(2, ni)*x2 + NKsi(3, ni)*x3 + NKsi(4, ni)*x4)*1.0f;
  93. }
  94. float J01(float y1, float y2, float y3, float y4, float ni){
  95.     return (NKsi(1, ni)*y1 + NKsi(2, ni)*y2 + NKsi(3, ni)*y3 + NKsi(4, ni)*y4)*1.0f;
  96. }
  97. float J10(int x1, int x2, int x3, int x4, float ksi){
  98.  
  99.     return (Nni(1, ksi)*x1 + Nni(2, ksi)*x2 + Nni(3, ksi)*x3 + Nni(4, ksi)*x4)*1.0f;
  100. }
  101. float J11(float y1, float y2, float y3, float y4, float ksi){
  102.     return (Nni(1, ksi)*y1 + Nni(2, ksi)*y2 + Nni(3, ksi)*y3 + Nni(4, ksi)*y4)*1.0f;
  103. }
  104. float detJfun(float x1, float x2, float x3, float x4){
  105.     return (x1*x4) - (x2*x3);
  106. }
  107. float detJminusFun(float x1, float x2, float x3, float x4, float det){
  108.     return (x1*x4) - (x2*x3);
  109. }
  110. class ElementSkonczony
  111. {
  112. public:
  113.     float x1, x2, x3, x4;
  114.     float y1, y2, y3, y4;       //położenie węzłów elementu skonczonego w ukladzie globalnym
  115.     float t01, t02, t03, t04;       //temperatura poczatkowa w poszczególnych węzłach C
  116.     float K;            //Wspolczynnik przewodzenia ciapla W/m*C
  117.     int NOP[4];         //NOP1, NOP2, NOP3, NOP4- wspolczynniki w odniesieniu do macierzy globalnej
  118.     float K_l[4][4];    //K_l- lokalna macierz wspolczynnikow ukladu rownan
  119.     float F_l[4];       //F_l- lokalny wektor prawej czesci ukladu
  120.     float J[4][4];
  121.     float detJ[4];
  122.     float Jminus[4][4];
  123.     float detJminus[4];
  124.     float pochodnaLokalna[4][2];
  125.     int war = 0;
  126.     int warPow = 0;
  127.     //float JElem[2][2];
  128.     ElementSkonczony(int NOP1n, int NOP2n, int NOP3n, int NOP4n, float x1n, float y1n, float x2n, float y2n, float x3n, float y3n, float x4n, float y4n, float t0, int warn, int warPown){ //Konstruktor elementu skonczonego
  129.         x1 = x1n;
  130.         x2 = x2n;
  131.         x3 = x3n;
  132.         x4 = x4n;
  133.         y1 = y1n;
  134.         y2 = y2n;
  135.         y3 = y3n;
  136.         y4 = y4n;
  137.         NOP[0] = NOP1n;
  138.         NOP[1] = NOP2n;
  139.         NOP[2] = NOP3n;
  140.         NOP[3] = NOP4n;
  141.         this->war = warn;
  142.         this->warPow = warPown;
  143.         t01 = t02 = t03 = t04 = t0;
  144.         for (int i = 0; i<4; ++i)
  145.             for (int j = 0; j<4; ++j){
  146.                 K_l[i][j] = 0;
  147.             }
  148.     }
  149.     void jakobiMatrix(){
  150.         J[0][0] = J00(x1, x2, x3, x4, -0.577);
  151.         J[0][1] = J01(y1, y2, y3, y4, -0.577);
  152.         J[0][2] = J10(x1, x2, x3, x4, -0.577);
  153.         J[0][3] = J11(y1, y2, y3, y4, -0.577);
  154.         this->detJ[0] = detJfun(J[0][0], J[0][1], J[0][2], J[0][3]);
  155.         Jminus[0][0] = J[0][3] / detJ[0];
  156.         Jminus[0][1] = -J[0][1] / detJ[0];
  157.         Jminus[0][2] = -J[0][2] / detJ[0];
  158.         Jminus[0][3] = J[0][0] / detJ[0];
  159.  
  160.         J[1][0] = J00(x1, x2, x3, x4, -0.577);
  161.         J[1][1] = J01(y1, y2, y3, y4, 0.577);
  162.         J[1][2] = J10(x1, x2, x3, x4, -0.577);
  163.         J[1][3] = J11(y1, y2, y3, y4, 0.577);
  164.         this->detJ[1] = detJfun(J[1][0], J[1][1], J[1][2], J[1][3]);
  165.         Jminus[1][0] = J[1][3] / detJ[1];
  166.         Jminus[1][1] = -J[1][1] / detJ[1];
  167.         Jminus[1][2] = -J[1][2] / detJ[1];
  168.         Jminus[1][3] = J[1][0] / detJ[1];
  169.  
  170.         J[2][0] = J00(x1, x2, x3, x4, 0.577);
  171.         J[2][1] = J01(y1, y2, y3, y4, 0.577);
  172.         J[2][2] = J10(x1, x2, x3, x4, 0.577);
  173.         J[2][3] = J11(y1, y2, y3, y4, 0.577);
  174.         this->detJ[2] = detJfun(J[2][0], J[2][1], J[2][2], J[2][3]);
  175.         Jminus[2][0] = J[2][3] / detJ[2];
  176.         Jminus[2][1] = -J[2][1] / detJ[2];
  177.         Jminus[2][2] = -J[2][2] / detJ[2];
  178.         Jminus[2][3] = J[2][0] / detJ[2];
  179.  
  180.         J[3][0] = J00(x1, x2, x3, x4, 0.577);
  181.         J[3][1] = J01(y1, y2, y3, y4, -0.577);
  182.         J[3][2] = J10(x1, x2, x3, x4, 0.577);
  183.         J[3][3] = J11(y1, y2, y3, y4, -0.577);
  184.         this->detJ[3] = detJfun(J[3][0], J[3][1], J[3][2], J[3][3]);
  185.         Jminus[3][0] = J[3][3] / detJ[3];
  186.         Jminus[3][1] = -J[3][1] / detJ[3];
  187.         Jminus[3][2] = -J[3][2] / detJ[3];
  188.         Jminus[3][3] = J[3][0] / detJ[3];
  189.     }
  190.     void liczPochodne(){
  191.         for (int j = 0; j<4; j++){
  192.             for (int i = 0; i<4; i++){//uzupelnianie macierzy lokalnej
  193.                 for (int l = 0; l<4; l++){
  194.                     float dndxl = NKsi(i + 1, ksi[j][1]) * Jminus[0][j] + Nni(i + 1, ksi[j][0]) * Jminus[j][1];
  195.                     float dndyl = NKsi(i + 1, ksi[j][1]) * Jminus[2][j] + Nni(i + 1, ksi[j][0]) * Jminus[j][3];
  196.                     float dndxk = NKsi(i + 1, ksi[j][1]) * Jminus[0][j] + Nni(i + 1, ksi[j][0]) * Jminus[j][1];
  197.                     float dndyk = NKsi(i + 1, ksi[j][1]) * Jminus[2][j] + Nni(i + 1, ksi[j][0]) * Jminus[j][3];
  198.                     K_l[i][l] += k*detJ[i] * (dndxl*dndxk + dndyl*dndyk);
  199.  
  200.                 }
  201.             }
  202.         }
  203.         if (war == 1){ // dodajemy strumien ciepla q
  204.  
  205.             K_l[warPow][warPow] += q / 2;
  206.             K_l[warPow + 1][warPow + 1] += q / 2;
  207.             cout << q << "   ";
  208.         }
  209.         else if (war == 2){ //dodajemy wymniane ciepla alfa
  210.             K_l[warPow][warPow] += 2.0 / 6.0*alfa;
  211.             K_l[warPow + 1][warPow + 1] += 2.0 / 6.0*alfa;
  212.             K_l[warPow + 1][warPow] += 1.0 / 6.0*alfa;
  213.             K_l[warPow][warPow + 1] += 1.0 / 6.0*alfa;
  214.             F_l[warPow] = 1 / 2 * alfa*tinf;
  215.             F_l[warPow + 1] = 1 / 2 * alfa*tinf;
  216.             cout << "   " << alfa;
  217.         }
  218.     }
  219.    
  220.     void macierzLokalnaWyswietl(){
  221.         cout << endl;
  222.  
  223.         for (int i = 0; i<4; ++i){
  224.             for (int j = 0; j<4; ++j){
  225.                 cout << J[i][j] << "  ";
  226.                 if (j % 2 == 1) cout << endl;
  227.             }
  228.             cout << endl;
  229.         }
  230.  
  231.         cout << x1 << "  " << x2 << "  " << x3 << "  " << x4 << endl;
  232.         cout << y1 << "  " << y2 << "  " << y3 << "  " << y4 << endl;
  233.  
  234.         for (int i = 0; i<4; ++i){
  235.             for (int j = 0; j<4; ++j){
  236.                 cout << K_l[i][j] << " ";
  237.             }
  238.             cout << endl;
  239.         }
  240.     }
  241. };
  242.  
  243.  
  244. bool gauss(int n, float ** AB, float * T)
  245. {
  246.     int i, j, k;
  247.     float m, s;
  248.  
  249.     //Eliminacja wspolczynnikow
  250.     for (i = 0; i < n - 1; i++){
  251.         for (j = i + 1; j < n; j++){
  252.             if (fabs(AB[i][i]) < eps) return false;
  253.             m = -AB[j][i] / AB[i][i]; // obliczamy mnoznik
  254.             for (k = i + 1; k <= n; k++){
  255.                 AB[j][k] += m * AB[i][j];
  256.             }
  257.         }
  258.     }
  259.  
  260.     //Wyliczanie niewiadomych
  261.     for (i = n - 1; i >= 0; i--){
  262.         s = AB[i][n];
  263.         for (j = n - 1; j >= i + 1; j--){
  264.             s -= AB[i][j] * T[j];
  265.         }
  266.         if (fabs(AB[i][i]) < eps) return false;
  267.         T[i] = s / AB[i][i];
  268.     }
  269.     return true;
  270. }
  271.  
  272. int main(){
  273.     int NOP1n, NOP2n, NOP3n, NOP4n;
  274.     int ne, nh;
  275.     float x1n, x2n, x3n, x4n, y1n, y2n, y3n, y4n;
  276.     int war, warPow;
  277.     fstream plik;
  278.     plik.open("dane.txt", ios::in); //Otwieranie pliku z danymi
  279.     if (!plik.good())
  280.         return false;
  281.  
  282.     plik >> ne >> t0 >> tinf >> k >> alfa >> q,
  283.         nh = ne + 1;
  284.  
  285.     ElementSkonczony** tabElementowSkonczonych = new ElementSkonczony*[ne]; //Tworzymy tablice elementow skonczonych
  286.  
  287.     for (int i = 0; i<ne; i++){
  288.         plik >> NOP1n >> NOP2n >> NOP3n >> NOP4n >> x1n >> y1n >> x2n >> y2n >> x3n >> y3n >> x4n >> y4n >> war >> warPow;
  289.         tabElementowSkonczonych[i] = new ElementSkonczony(NOP1n, NOP2n, NOP3n, NOP4n, x1n, y1n, x2n, y2n, x3n, y3n, x4n, y4n, t0, war, warPow);
  290.     }
  291.     for (int i = 0; i<ne; i++){
  292.         tabElementowSkonczonych[i]->jakobiMatrix();
  293.     }
  294.  
  295.  
  296.     for (int i = 0; i<ne; i++){
  297.         tabElementowSkonczonych[i]->liczPochodne();
  298.     }
  299.  
  300.     nh = 21;
  301.     float **AB = new float*[nh];        //Macierz AB- pomocnicza macierz do obliczenia ukladu rownan metoda Gaussa
  302.     float** KG = new float*[nh];        //Macierz KG- globalna macierz wspolczynnikow uklady rownan
  303.     float* FG = new float[nh];          //F- wektor prawej czesci ukladu rownan
  304.     float* T = new float[nh];           //T- wektor niewiadomych temperatur
  305.  
  306.     for (int i = 0; i<nh; i++)
  307.         AB[i] = new float[nh + 1];
  308.     for (int i = 0; i<nh; i++){
  309.         T[i] = t0;
  310.         KG[i] = new float[nh];
  311.     }
  312.  
  313.     for (int i = 0; i<nh; i++)
  314.     {
  315.         for (int j = 0; j<nh; j++)
  316.             KG[i][j] = 0;
  317.         FG[i] = 0;
  318.         T[i] = 0;
  319.     }
  320.     for (int i = 0; i<12; ++i){
  321.         for (int k = 0; k<4; k++){ // uzupelnianie macierzy globalnej
  322.             for (int l = 0; l<4; l++){
  323.                 KG[tabElementowSkonczonych[i]->NOP[k]][tabElementowSkonczonych[i]->NOP[l]] += tabElementowSkonczonych[i]->K_l[k][l];
  324.  
  325.             }
  326.             FG[tabElementowSkonczonych[i]->NOP[k]] += tabElementowSkonczonych[i]->F_l[k];
  327.         }
  328.  
  329.     }
  330.  
  331.  
  332.     for (int i = 0; i < nh; i++){
  333.         for (int j = 0; j < nh; j++)
  334.             AB[i][j] = KG[i][j];
  335.     }
  336.     for (int i = 0; i < nh; i++){
  337.         AB[i][nh] = FG[i];
  338.     }
  339.  
  340.     if (gauss(nh, AB, T)){  //Wywolanie funkcji obliczajacej uklad rownan metoda Gaussa
  341.         cout << endl << endl;
  342.         cout << "Niewiadome temperatury w wezlach wynosza : " << endl << endl;
  343.         for (int i = 0; i < nh; i++){
  344.             cout << "t" << i + 1 << " = " << setw(9) << T[i] << endl;
  345.         }
  346.     }
  347.     else{
  348.         cout << "DZIELNIK ZERO\n";
  349.         cout << endl;
  350.        
  351.    
  352.     }
  353.     system("PAUSE");
  354.     return 0;
  355. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!