fork download
copy 
  1. #include <iostream>
  2. #include <vector>
  3. #include <cmath>
  4. #include <cassert>
  5.  
  6. using namespace std;
  7.  
  8. // Estructura para representar un nodo
  9. struct Nodo {
  10.     double x, y;  // Coordenadas del nodo
  11.     double fx, fy;  // Fuerzas aplicadas en el nodo (si las hay)
  12.     bool fijo_x, fijo_y; // Condiciones de frontera (si el nodo es fijo en X o Y)
  13. };
  14.  
  15. // Estructura para representar una barra entre dos nodos
  16. struct Barra {
  17.     int nodo1, nodo2;  // Índices de los nodos conectados
  18.     double longitud;   // Longitud de la barra
  19.     double area;       // Área de la sección transversal
  20.     double E;          // Módulo de elasticidad
  21. };
  22.  
  23. // Función para calcular la rigidez de una barra (en 2D)
  24. vector<vector<double>> rigidezBarra(const Barra& barra, const Nodo& n1, const Nodo& n2) {
  25.     double L = barra.longitud;
  26.     double A = barra.area;
  27.     double E = barra.E;
  28.  
  29.     // Cálculo de la matriz de rigidez local de la barra
  30.     double c = (n2.x - n1.x) / L;  // Coseno de la dirección de la barra
  31.     double s = (n2.y - n1.y) / L;  // Seno de la dirección de la barra
  32.  
  33.     vector<vector<double>> K(4, vector<double>(4, 0));  // Matriz de rigidez local (4x4)
  34.  
  35.     // Matriz de rigidez local (considerando 2D)
  36.     double k_local = (A * E) / L;
  37.     K[0][0] = k_local * c * c;
  38.     K[0][1] = k_local * c * s;
  39.     K[1][0] = k_local * c * s;
  40.     K[1][1] = k_local * s * s;
  41.     K[2][2] = k_local * c * c;
  42.     K[2][3] = k_local * c * s;
  43.     K[3][2] = k_local * c * s;
  44.     K[3][3] = k_local * s * s;
  45.  
  46.     return K;
  47. }
  48.  
  49. // Función para ensamblar la matriz de rigidez global
  50. vector<vector<double>> ensamblarRigidezGlobal(const vector<Barra>& barras, const vector<Nodo>& nodos) {
  51.     int n = nodos.size();
  52.     vector<vector<double>> K_global(2 * n, vector<double>(2 * n, 0));  // Matriz de rigidez global (2n x 2n)
  53.  
  54.     // Ensamblaje de la matriz de rigidez global
  55.     for (const Barra& barra : barras) {
  56.         Nodo n1 = nodos[barra.nodo1];
  57.         Nodo n2 = nodos[barra.nodo2];
  58.  
  59.         // Obtener la matriz de rigidez local
  60.         vector<vector<double>> K_local = rigidezBarra(barra, n1, n2);
  61.  
  62.         // Incluir K_local en la matriz de rigidez global
  63.         K_global[2 * barra.nodo1][2 * barra.nodo1] += K_local[0][0];
  64.         K_global[2 * barra.nodo1][2 * barra.nodo1 + 1] += K_local[0][1];
  65.         K_global[2 * barra.nodo1 + 1][2 * barra.nodo1] += K_local[1][0];
  66.         K_global[2 * barra.nodo1 + 1][2 * barra.nodo1 + 1] += K_local[1][1];
  67.  
  68.         K_global[2 * barra.nodo2][2 * barra.nodo2] += K_local[2][2];
  69.         K_global[2 * barra.nodo2][2 * barra.nodo2 + 1] += K_local[2][3];
  70.         K_global[2 * barra.nodo2 + 1][2 * barra.nodo2] += K_local[3][2];
  71.         K_global[2 * barra.nodo2 + 1][2 * barra.nodo2 + 1] += K_local[3][3];
  72.  
  73.         // Ensamblaje de los términos de interacción entre los nodos
  74.         K_global[2 * barra.nodo1][2 * barra.nodo2] += K_local[0][2];
  75.         K_global[2 * barra.nodo1][2 * barra.nodo2 + 1] += K_local[0][3];
  76.         K_global[2 * barra.nodo1 + 1][2 * barra.nodo2] += K_local[1][2];
  77.         K_global[2 * barra.nodo1 + 1][2 * barra.nodo2 + 1] += K_local[1][3];
  78.  
  79.         K_global[2 * barra.nodo2][2 * barra.nodo1] += K_local[2][0];
  80.         K_global[2 * barra.nodo2][2 * barra.nodo1 + 1] += K_local[2][1];
  81.         K_global[2 * barra.nodo2 + 1][2 * barra.nodo1] += K_local[3][0];
  82.         K_global[2 * barra.nodo2 + 1][2 * barra.nodo1 + 1] += K_local[3][1];
  83.     }
  84.  
  85.     return K_global;
  86. }
  87.  
  88. // Función para aplicar condiciones de frontera (restringir desplazamientos)
  89. void aplicarCondicionesFrontera(vector<vector<double>>& K, vector<double>& F, const vector<Nodo>& nodos) {
  90.     int n = nodos.size();
  91.  
  92.     for (int i = 0; i < n; i++) {
  93.         if (nodos[i].fijo_x) {
  94.             // Fijar el desplazamiento en X (columna correspondiente)
  95.             for (int j = 0; j < 2 * n; j++) {
  96.                 K[2 * i][j] = 0;
  97.                 K[j][2 * i] = 0;
  98.             }
  99.             K[2 * i][2 * i] = 1;
  100.             F[2 * i] = 0;  // Fuerza en X es 0
  101.         }
  102.  
  103.         if (nodos[i].fijo_y) {
  104.             // Fijar el desplazamiento en Y (columna correspondiente)
  105.             for (int j = 0; j < 2 * n; j++) {
  106.                 K[2 * i + 1][j] = 0;
  107.                 K[j][2 * i + 1] = 0;
  108.             }
  109.             K[2 * i + 1][2 * i + 1] = 1;
  110.             F[2 * i + 1] = 0;  // Fuerza en Y es 0
  111.         }
  112.     }
  113. }
  114.  
  115. // Función para resolver el sistema de ecuaciones (desplazamientos)
  116. vector<double> resolverDesplazamientos(vector<vector<double>>& K, vector<double>& F) {
  117.     int n = K.size();
  118.  
  119.     // Resolver el sistema Kx = F usando eliminación de Gauss
  120.     vector<double> X(n, 0);  // Desplazamientos inicializados a 0
  121.     vector<vector<double>> A = K;
  122.  
  123.     // Eliminación de Gauss
  124.     for (int i = 0; i < n; i++) {
  125.         // Buscar el máximo en la columna para estabilidad
  126.         int max_row = i;
  127.         for (int k = i + 1; k < n; k++) {
  128.             if (abs(A[k][i]) > abs(A[max_row][i])) {
  129.                 max_row = k;
  130.             }
  131.         }
  132.  
  133.         // Intercambiar filas
  134.         swap(A[i], A[max_row]);
  135.         swap(F[i], F[max_row]);
  136.  
  137.         // Hacer cero las entradas debajo de la diagonal
  138.         for (int j = i + 1; j < n; j++) {
  139.             double factor = A[j][i] / A[i][i];
  140.             for (int k = i; k < n; k++) {
  141.                 A[j][k] -= factor * A[i][k];
  142.             }
  143.             F[j] -= factor * F[i];
  144.         }
  145.     }
  146.  
  147.     // Resolución hacia atrás
  148.     for (int i = n - 1; i >= 0; i--) {
  149.         X[i] = F[i] / A[i][i];
  150.         for (int j = i - 1; j >= 0; j--) {
  151.             F[j] -= A[j][i] * X[i];
  152.         }
  153.     }
  154.  
  155.     return X;
  156. }
  157.  
  158. int main() {
  159.     int num_nodos = 3;
  160.     int num_barras = 3;
  161.  
  162.     // Definir nodos y barras
  163.     vector<Nodo> nodos = {{0, 0, 0, 0, true, true}, {0, 5, 0, 0, false, true}, {2.5, 2.5, 0, -10000, false, false}};
  164.     vector<Barra> barras = {{0, 1, 5.0, 0.01, 210e9}, {1, 2, sqrt(2.5*2.5 + 2.5*2.5), 0.01, 210e9}, {0, 2, sqrt(2.5*2.5 + 2.5*2.5), 0.01, 210e9}};
  165.  
  166.     // Ensamblar la matriz de rigidez global
  167.     vector<vector<double>> K_global = ensamblarRigidezGlobal(barras, nodos);
  168.  
  169.     // Definir el vector de fuerzas (aplicadas en los nodos)
  170.     vector<double> F(2 * nodos.size(), 0);
  171.     F[5] = nodos[2].fy;  // Fuerza de 10 kN en el nodo 3 (Y)
  172.  
  173.     // Aplicar condiciones de frontera
  174.     aplicarCondicionesFrontera(K_global, F, nodos);
  175.  
  176.     // Resolver el sistema de ecuaciones para obtener los desplazamientos
  177.     vector<double> desplazamientos = resolverDesplazamientos(K_global, F);
  178.  
  179.     // Mostrar los resultados de los desplazamientos
  180.     for (int i = 0; i < desplazamientos.size(); ++i) {
  181.         cout << "Desplazamiento en nodo " << i / 2 << " (componente " << i % 2 << "): " 
  182.              << desplazamientos[i] << " m" << endl;
  183.     }
  184.  
  185.     return 0;
  186. }
  187.  
Success #stdin #stdout 0.01s 5292KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Desplazamiento en nodo 0 (componente 0): 0 m
Desplazamiento en nodo 0 (componente 1): 0 m
Desplazamiento en nodo 1 (componente 0): 0 m
Desplazamiento en nodo 1 (componente 1): 0 m
Desplazamiento en nodo 2 (componente 0): 0 m
Desplazamiento en nodo 2 (componente 1): -1.68359e-05 m
https://ideone.com/4DN6rP
language:
C++ (gcc 8.3)
created:
3 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!