First partially-implemented functional release
[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 int
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 if (gsl_linalg_LU_decomp (&m.matrix, p, &s))
39 {
40 fprintf (stderr, "__get_basis: cannot build basic polynomial for triangle [%d], %d, %d: LU decomposition failed\n",
41 center->index, p2->index, p3->index);
42
43 gsl_permutation_free (p);
44 gsl_vector_free (x);
45
46 return -1;
47 }
48
49
50 if (gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x))
51 {
52 fprintf (stderr, "__get_basis: cannot build basic polynomial for triangle [%d], %d, %d: LU solver failed\n",
53 center->index, p2->index, p3->index);
54 gsl_permutation_free (p);
55 gsl_vector_free (x);
56
57 return -1;
58 }
59
60 out->a = gsl_vector_get (x, 0);
61 out->b = gsl_vector_get (x, 1);
62 out->c = gsl_vector_get (x, 2);
63
64 gsl_permutation_free (p);
65 gsl_vector_free (x);
66
67 return 0;
68 }
69
70 INLINE int
71 get_basis (struct polynomial *out, struct triangle *triangle, int node)
72 {
73 return __get_basis (out, triangle->nodes[node % 3],
74 triangle->nodes[(node + 1) % 3],
75 triangle->nodes[(node + 2) % 3]);
76 }
77
78 INLINE float
79 triangle_evaluate (struct problem2d *problem,
80 struct polynomial *pu, struct polynomial *pv,
81 float x, float y)
82 {
83 float a_eval[2][2];
84 float a0;
85
86 problem->x = x;
87 problem->y = y;
88
89 /* Cálculo de la matriz A */
90 a_eval[0][0] = jit_eval (problem->a[A_11]);
91 a_eval[0][1] = jit_eval (problem->a[A_12]);
92 a_eval[1][0] = jit_eval (problem->a[A_21]);
93 a_eval[1][1] = jit_eval (problem->a[A_22]);
94
95 a0 = jit_eval (problem->a0);
96
97 /* DP{v_node}.x * (DP{u_node}.x * A11 + DP{u_node}.y * A12) +
98 DP{v_node}.y * (DP{u_node}.x * A21 + DP{u_node}.y * A22) +
99 + a0 * P{v_node} * P{u_node}; */
100
101 /* Nabla(P) = (a, b) */
102
103 /* Evaluación */
104 return pv->a * (pu->a * a_eval[0][0] + pu->b * a_eval[0][1]) +
105 pv->b * (pu->a * a_eval[1][0] + pu->b * a_eval[1][1]) +
106 a0 * (pv->a * x + pv->b * y + pv->c) * (pu->a * x + pu->b * y + pu->c);
107
108 }
109
110 INLINE float
111 triangle_integrate (struct problem2d *problem, struct triangle *triangle, struct polynomial *u_p, struct polynomial *v_p)
112 {
113 float b1x, b1y;
114 float b2x, b2y;
115 float b3x, b3y;
116 float I1, I2, I3;
117
118 float area;
119
120
121
122
123 /*
124 P = ax + by + c
125 DP = (a, b) */
126 /* Aquí tenemos que evaluar lo siguiente:
127
128 _
129 |
130 | (DPt x A x DP + a0 * Pt x P) dT
131 _|
132 T
133
134 Donde:
135 DP -> 2 x 3
136 P -> 1 x 3
137
138 Estamos calculando el punto (v_node, u_node) de la matriz elemental. Esto
139 quiere decir, que estamos calculando la integral de:
140
141 f(x, y) = (DP{v_node})t x A x (DP{u_node}) + a0 * P{v_node} * P{u_node}
142
143
144 Con el número entre llaves la columna de la matriz. Entonces:
145
146 f(x, y) = (DP{v_node})t x (DP{u_node}.x * A11 + DP{u_node}.y * A12,
147 DP{u_node}.x * A21 + DP{u_node}.y * A22) +
148 + a0 * P{v_node} * P{u_node} =
149 = DP{v_node}.x * (DP{u_node}.x * A11 + DP{u_node}.y * A12) +
150 DP{v_node}.y * (DP{u_node}.x * A21 + DP{u_node}.y * A22) +
151 + a0 * P{v_node} * P{u_node};
152
153 de esta función se va a hacer mediante la fórmula de cuadratura:
154
155 I ~ |T| / 3 * (f(b12) + f(b23) + f(b31))
156
157 La cual es exacta para polinomios de orden 2 (como es nuestro caso).
158
159 Con |T| el área del triángulo, y bXY el punto medio del lado XY del
160 triángulo. Estos puntos medios se pueden calcular como una simple media,
161 el área del triángulo se puede calcular como la mitad del paralelogramo
162 formado por dos de sus lados.
163
164 Ahora nos queda calcular las expresiones de estos polinomios de base
165 centrados en v_node y u_node, y podremos integrar.
166
167 Veamos como construir P. Si P está centrado en ux, uy:
168
169 P(x, y) = a + bx + cx (elemento P1 triangular)
170 P(ux, uy) = 1
171 P(ux, uy)
172
173 */
174
175 /* Cálculo de los puntos medio */
176 b1x = 0.5 * (triangle->nodes[0]->x + triangle->nodes[1]->x);
177 b1y = 0.5 * (triangle->nodes[0]->y + triangle->nodes[1]->y);
178
179 b2x = 0.5 * (triangle->nodes[1]->x + triangle->nodes[2]->x);
180 b2y = 0.5 * (triangle->nodes[1]->y + triangle->nodes[2]->y);
181
182 b3x = 0.5 * (triangle->nodes[2]->x + triangle->nodes[0]->x);
183 b3y = 0.5 * (triangle->nodes[2]->y + triangle->nodes[0]->y);
184
185 /* Cálculo del área. Esto se puede hacer directamente a través de un determinante */
186 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));
187
188 /* Cálculo de las integrales parciales */
189 I1 = triangle_evaluate (problem, u_p, v_p, b1x, b1y);
190 I2 = triangle_evaluate (problem, u_p, v_p, b2x, b2y);
191 I3 = triangle_evaluate (problem, u_p, v_p, b3x, b3y);
192
193 /* Cuadratura */
194 return area / 0.3 * (I1 + I2 + I3);
195 }
196
197 INLINE float
198 triangle_integrate_second (struct problem2d *problem, struct triangle *triangle, struct polynomial *v_p)
199 {
200 float bx, by;
201 float area;
202 float integral;
203 float f;
204
205 /* Aquí bastará con un trapecio */
206
207
208 /* La función a integrar es:
209 P * f (+ términos de flujo) */
210
211 /* La cuadratura del trapecio (en su equivalente 2D) es:
212 |T| * P * f (baricentro) */
213
214 /* Calculemos el baricentro */
215 bx = (triangle->nodes[0]->x + triangle->nodes[1]->x + triangle->nodes[2]->x) / 3.0;
216 by = (triangle->nodes[0]->y + triangle->nodes[1]->y + triangle->nodes[2]->y) / 3.0;
217
218 /* Área del triángulo, paralelogramo de los lados entre dos. Esto no es más que un determinante */
219 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));
220
221 problem->x = bx;
222 problem->y = by;
223 f = jit_eval (problem->f);
224
225 integral = area * (v_p->a * bx + v_p->b * by + v_p->c) * f;
226
227 return integral;
228 }
229
230 #define xfprintf(fp, fmt, arg...) \
231 do \
232 { \
233 if (fprintf (fp, fmt, ##arg) < 0) \
234 return -1; \
235 } \
236 while (0)
237
238 int
239 export_matrix_and_vector (struct problem2d *problem, FILE *fp)
240 {
241 int i, j;
242
243 xfprintf (fp, "A = [");
244
245 for (i = 0; i < problem->A->size1; i++)
246 {
247 if (i)
248 xfprintf (fp, ";\n");
249
250 for (j = 0; j < problem->A->size2; j++)
251 xfprintf (fp, "%s%lf", j ? ", " : "", gsl_matrix_get (problem->A, i, j));
252 }
253
254 xfprintf (fp, "];\n");
255
256 xfprintf (fp, "b = [");
257
258 for (i = 0; i < problem->b->size; i++)
259 xfprintf (fp, "%s%lf", i ? "; " : "", gsl_vector_get (problem->b, i));
260
261 xfprintf (fp, "];\n");
262 }
263
264 int
265 build_matrix_and_vector (struct problem2d *problem)
266 {
267 int i, j, k;
268 float integral;
269 double prev;
270 struct polynomial u_p, v_p;
271
272 /* Generamos una matriz de rigidez que va a tener NxN, con
273 N = número de nodos. */
274
275 problem->A = gsl_matrix_calloc (problem->geometry->point_count,
276 problem->geometry->point_count);
277
278 problem->b = gsl_vector_calloc (problem->geometry->point_count);
279
280 for (i = 0; i < problem->geometry->triangle_count; i++)
281 {
282 /* Esta parte sería la que calcularía la matriz de rigidez elemental
283 para este triángulo. */
284
285 /* Por cada triángulo, vamos a integrar contra cada función básica cen-
286 trada en el nodo triangle_list[i]->nodes[j]. Aquí "j" vendría
287 a indexar la fila de la matriz elemental, es la que indexa los
288 "vh" en la forma bilineal a(uh, vh) */
289
290 for (j = 0; j < TRIANGLE_NODES; j++)
291 {
292 if (get_basis (&v_p, problem->geometry->triangle_list[i], j) == -1)
293 return -1;
294
295 /* Podemos ir construyendo el segundo miembro: */
296 prev = gsl_vector_get (problem->b, problem->geometry->triangle_list[i]->nodes[j]->index);
297
298 integral = triangle_integrate_second (problem, problem->geometry->triangle_list[i], &v_p);
299
300 gsl_vector_set (problem->b, problem->geometry->triangle_list[i]->nodes[j]->index, prev + integral);
301
302 /* Y por cada una de ellas, integraremos en todo el triángulo cuando
303 uh es una función básica centrada en triangle_list[i]->nodes[k].
304 Esto no hace más que indexar los uh en a(uh, vh) */
305
306 for (k = 0; k < TRIANGLE_NODES; k++)
307 {
308 if (get_basis (&u_p, problem->geometry->triangle_list[i], k) == -1)
309 return -1;
310
311 /* Nótese que el lugar correspondiente en la matriz de rigidez
312 para esta integral es (nodes[j]->index, nodes[k]->index) */
313
314 prev = gsl_matrix_get (problem->A, problem->geometry->triangle_list[i]->nodes[j]->index,
315 problem->geometry->triangle_list[i]->nodes[k]->index);
316
317 integral = triangle_integrate (problem, problem->geometry->triangle_list[i], &u_p, &v_p);
318
319 gsl_matrix_set (problem->A, problem->geometry->triangle_list[i]->nodes[j]->index,
320 problem->geometry->triangle_list[i]->nodes[k]->index,
321 prev + integral);
322 }
323 }
324 }
325
326 return 0;
327 }
328
329 /* TODO: añadir estados */
330 int
331 problem_solve (struct problem2d *problem)
332 {
333 gsl_permutation * p;
334 problem->solution = gsl_vector_alloc (problem->b->size);
335 int s;
336
337 p = gsl_permutation_alloc (problem->A->size1);
338
339 if (gsl_linalg_LU_decomp (problem->A, p, &s))
340 {
341 fprintf (stderr, "problem_solve: cannot perform LU decomposition!\n");
342
343 gsl_permutation_free (p);
344 return -1;
345 }
346
347 if (gsl_linalg_LU_solve (problem->A, p, problem->b, problem->solution))
348 {
349 fprintf (stderr, "problem_solve: cannot solve via LU!\n");
350 gsl_permutation_free (p);
351
352 return -1;
353 }
354
355 gsl_permutation_free (p);
356
357
358 return 0;
359 }
360
361