80cf2adc871b4a156cee6a41692e482a0894c11c
[2dsolver.git] / src / solver.c
1 /* Aquí entrará todo el tema de la construcción de la matriz de rigidez
2 y demás.
3
4 Para acelerar los cálculos, omitiremos los productos de la matriz de
5 ordenación B en la generación de la matriz de rigidez, colocando las
6 matrices elementales directamente sobre la matriz global.
7
8 */
9
10 #include <2dsolve.h>
11
12 struct polynomial
13 {
14 float a, b, c;
15 };
16
17 /* Muy feo, pero funcional */
18 INLINE void
19 __get_basis (struct polynomial *out, struct point *center, struct point *p2, struct point *p3)
20 {
21 double a_data[] = { center->x, center->y, 1.0,
22 p2->x, p2->y, 1.0,
23 p3->x, p3->y, 1.0 };
24
25 double b_data[] = { 1.0, 0.0, 0.0 };
26 int s;
27
28 gsl_matrix_view m
29 = gsl_matrix_view_array (a_data, 3, 3);
30
31 gsl_vector_view b
32 = gsl_vector_view_array (b_data, 3);
33
34 gsl_vector *x = gsl_vector_alloc (3);
35
36 gsl_permutation * p = gsl_permutation_alloc (3);
37
38 gsl_linalg_LU_decomp (&m.matrix, p, &s);
39
40 gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
41
42 out->a = gsl_vector_get (x, 0);
43 out->b = gsl_vector_get (x, 1);
44 out->c = gsl_vector_get (x, 2);
45
46 gsl_permutation_free (p);
47 gsl_vector_free (x);
48 }
49
50 INLINE void
51 get_basis (struct polynomial *out, struct triangle *triangle, int node)
52 {
53 __get_basis (out, triangle->nodes[node % 3],
54 triangle->nodes[(node + 1) % 3],
55 triangle->nodes[(node + 2) % 3]);
56 }
57
58 INLINE float
59 triangle_evaluate (struct problem2d *problem,
60 struct polynomial *pu, struct polynomial *pv,
61 float x, float y)
62 {
63 float a_eval[2][2];
64 float a0;
65
66 problem->x = x;
67 problem->y = y;
68
69 /* Cálculo de la matriz A */
70 a_eval[0][0] = jit_eval (problem->a[A_11]);
71 a_eval[0][1] = jit_eval (problem->a[A_12]);
72 a_eval[1][0] = jit_eval (problem->a[A_21]);
73 a_eval[1][1] = jit_eval (problem->a[A_22]);
74
75 a0 = jit_eval (problem->a0);
76
77 /* DP{v_node}.x * (DP{u_node}.x * A11 + DP{u_node}.y * A12) +
78 DP{v_node}.y * (DP{u_node}.x * A21 + DP{u_node}.y * A22) +
79 + a0 * P{v_node} * P{u_node}; */
80
81 /* Nabla(P) = (a, b) */
82
83 /* Evaluación */
84 return pv->a * (pu->a * a_eval[0][0] + pu->b * a_eval[0][1]) +
85 pv->b * (pu->a * a_eval[1][0] + pu->b * a_eval[1][1]) +
86 a0 * (pv->a * x + pv->b * y + pv->c) * (pu->a * x + pu->b * y + pu->c);
87
88 }
89
90 INLINE float
91 triangle_integrate (struct problem2d *problem, struct triangle *triangle, struct polynomial *u_p, struct polynomial *v_p)
92 {
93 float b1x, b1y;
94 float b2x, b2y;
95 float b3x, b3y;
96 float I1, I2, I3;
97
98 float area;
99
100
101
102
103 /*
104 P = ax + by + c
105 DP = (a, b) */
106 /* Aquí tenemos que evaluar lo siguiente:
107
108 _
109 |
110 | (DPt x A x DP + a0 * Pt x P) dT
111 _|
112 T
113
114 Donde:
115 DP -> 2 x 3
116 P -> 1 x 3
117
118 Estamos calculando el punto (v_node, u_node) de la matriz elemental. Esto
119 quiere decir, que estamos calculando la integral de:
120
121 f(x, y) = (DP{v_node})t x A x (DP{u_node}) + a0 * P{v_node} * P{u_node}
122
123
124 Con el número entre llaves la columna de la matriz. Entonces:
125
126 f(x, y) = (DP{v_node})t x (DP{u_node}.x * A11 + DP{u_node}.y * A12,
127 DP{u_node}.x * A21 + DP{u_node}.y * A22) +
128 + a0 * P{v_node} * P{u_node} =
129 = DP{v_node}.x * (DP{u_node}.x * A11 + DP{u_node}.y * A12) +
130 DP{v_node}.y * (DP{u_node}.x * A21 + DP{u_node}.y * A22) +
131 + a0 * P{v_node} * P{u_node};
132
133 de esta función se va a hacer mediante la fórmula de cuadratura:
134
135 I ~ |T| / 3 * (f(b12) + f(b23) + f(b31))
136
137 La cual es exacta para polinomios de orden 2 (como es nuestro caso).
138
139 Con |T| el área del triángulo, y bXY el punto medio del lado XY del
140 triángulo. Estos puntos medios se pueden calcular como una simple media,
141 el área del triángulo se puede calcular como la mitad del paralelogramo
142 formado por dos de sus lados.
143
144 Ahora nos queda calcular las expresiones de estos polinomios de base
145 centrados en v_node y u_node, y podremos integrar.
146
147 Veamos como construir P. Si P está centrado en ux, uy:
148
149 P(x, y) = a + bx + cx (elemento P1 triangular)
150 P(ux, uy) = 1
151 P(ux, uy)
152
153 */
154
155 /* Cálculo de los puntos medio */
156 b1x = 0.5 * (triangle->nodes[0]->x + triangle->nodes[1]->x);
157 b1y = 0.5 * (triangle->nodes[0]->y + triangle->nodes[1]->y);
158
159 b2x = 0.5 * (triangle->nodes[1]->x + triangle->nodes[2]->x);
160 b2y = 0.5 * (triangle->nodes[1]->y + triangle->nodes[2]->y);
161
162 b3x = 0.5 * (triangle->nodes[2]->x + triangle->nodes[0]->x);
163 b3y = 0.5 * (triangle->nodes[2]->y + triangle->nodes[0]->y);
164
165 /* Cálculo del área. Esto se puede hacer directamente a través de un determinante */
166 area = 0.5 * fabs ((triangle->nodes[1]->x - triangle->nodes[0]->x) * (triangle->nodes[2]->y - triangle->nodes[0]->y) - (triangle->nodes[1]->y - triangle->nodes[0]->y) * (triangle->nodes[2]->x - triangle->nodes[0]->x));
167
168 /* Cálculo de las integrales parciales */
169 I1 = triangle_evaluate (problem, u_p, v_p, b1x, b1y);
170 I2 = triangle_evaluate (problem, u_p, v_p, b2x, b2y);
171 I3 = triangle_evaluate (problem, u_p, v_p, b3x, b3y);
172
173 /* Cuadratura */
174 return area / 0.3 * (I1 + I2 + I3);
175 }
176
177 INLINE float
178 triangle_integrate_second (struct problem2d *problem, struct triangle *triangle, struct polynomial *v_p)
179 {
180 float bx, by;
181 float area;
182 float integral;
183 float f;
184
185 /* Aquí bastará con un trapecio */
186
187
188 /* La función a integrar es:
189 P * f (+ términos de flujo) */
190
191 /* La cuadratura del trapecio (en su equivalente 2D) es:
192 |T| * P * f (baricentro) */
193
194 /* Calculemos el baricentro */
195 bx = (triangle->nodes[0]->x + triangle->nodes[1]->x + triangle->nodes[2]->x) / 3.0;
196 by = (triangle->nodes[0]->y + triangle->nodes[1]->y + triangle->nodes[2]->y) / 3.0;
197
198 /* Área del triángulo, paralelogramo de los lados entre dos. Esto no es más que un determinante */
199 area = 0.5 * fabs ((triangle->nodes[1]->x - triangle->nodes[0]->x) * (triangle->nodes[2]->y - triangle->nodes[0]->y) - (triangle->nodes[1]->y - triangle->nodes[0]->y) * (triangle->nodes[2]->x - triangle->nodes[0]->x));
200
201 problem->x = bx;
202 problem->y = by;
203 f = jit_eval (problem->f);
204
205 integral = area * (v_p->a * bx + v_p->b * by + v_p->c) * f;
206
207 return integral;
208 }
209
210 int
211 build_matrix_and_vector (struct problem2d *problem, gsl_matrix **A, gsl_vector **b)
212 {
213 int i, j, k;
214 float integral;
215 double prev;
216 struct polynomial u_p, v_p;
217
218 /* Generamos una matriz de rigidez que va a tener NxN, con
219 N = número de nodos. */
220
221 *A = gsl_matrix_calloc (problem->geometry->point_count,
222 problem->geometry->point_count);
223
224 *b = gsl_vector_calloc (problem->geometry->point_count);
225
226 for (i = 0; i < problem->geometry->triangle_count; i++)
227 {
228 /* Esta parte sería la que calcularía la matriz de rigidez elemental
229 para este triángulo. */
230
231 /* Por cada triángulo, vamos a integrar contra cada función básica cen-
232 trada en el nodo triangle_list[i]->nodes[j]. Aquí "j" vendría
233 a indexar la fila de la matriz elemental, es la que indexa los
234 "vh" en la forma bilineal a(uh, vh) */
235
236 for (j = 0; j < TRIANGLE_NODES; j++)
237 {
238 get_basis (&v_p, problem->geometry->triangle_list[i], j);
239
240 /* Podemos ir construyendo el segundo miembro: */
241 prev = gsl_vector_get (*b, problem->geometry->triangle_list[i]->nodes[j]->index);
242
243 integral = triangle_integrate_second (problem, problem->geometry->triangle_list[i], &v_p);
244
245 gsl_vector_set (*b, problem->geometry->triangle_list[i]->nodes[j]->index, prev + integral);
246
247 /* Y por cada una de ellas, integraremos en todo el triángulo cuando
248 uh es una función básica centrada en triangle_list[i]->nodes[k].
249 Esto no hace más que indexar los uh en a(uh, vh) */
250
251 for (k = 0; k < TRIANGLE_NODES; k++)
252 {
253 get_basis (&u_p, problem->geometry->triangle_list[i], k);
254
255 /* Nótese que el lugar correspondiente en la matriz de rigidez
256 para esta integral es (nodes[j]->index, nodes[k]->index) */
257
258 prev = gsl_matrix_get (*A, problem->geometry->triangle_list[i]->nodes[j]->index,
259 problem->geometry->triangle_list[i]->nodes[k]->index);
260
261 integral = triangle_integrate (problem, problem->geometry->triangle_list[i], &u_p, &v_p);
262
263 gsl_matrix_set (*A, problem->geometry->triangle_list[i]->nodes[j]->index,
264 problem->geometry->triangle_list[i]->nodes[k]->index,
265 prev + integral);
266 }
267 }
268 }
269
270 return 0;
271 }
272
273