First partially-implemented functional release
authorGonzalo J. Carracedo <BatchDrake@gmail.com>
Mon, 23 Apr 2012 16:06:37 +0000 (18:06 +0200)
committerGonzalo J. Carracedo <BatchDrake@gmail.com>
Mon, 23 Apr 2012 16:06:37 +0000 (18:06 +0200)
configure.in
src/2dsolve.h
src/Makefile.am
src/config.h [new file with mode: 0644]
src/display.c [new file with mode: 0644]
src/geometry.c
src/main.c
src/problem.c [new file with mode: 0644]
src/solver.c

index 5a88cb1..8bbb6e3 100644 (file)
@@ -85,12 +85,23 @@ GLOBAL_LDFLAGS="$GLOBAL_LDFLAGS -lz"
 
 
 dnl Macro snippets imported from dependency `gsl'
-AC_CHECK_LIB([m],[cos])
-AC_CHECK_LIB([gslcblas],[cblas_dgemm])
-AC_CHECK_LIB([gsl],[gsl_blas_dgemm])
+AC_CHECK_LIB([m],[cos],, [ AC_MSG_ERROR(error: libm not found)])
+AC_CHECK_LIB([gslcblas],[cblas_dgemm],, [ AC_MSG_ERROR(error: libgslcblas not found (you may need to install GSL devel libraries))])
+AC_CHECK_LIB([gsl],[gsl_blas_dgemm],, [ AC_MSG_ERROR(error: libgsl not found (you may need to install GSL devel libraries))])
+
+AC_CHECK_HEADER(gsl/gsl_vector.h,
+                ,
+                [ AC_MSG_ERROR(error: GSL include files not found (you may need to install GSL devel libraries))])
+AC_CHECK_HEADER(gsl/gsl_matrix.h,
+                ,
+                [ AC_MSG_ERROR(error: GSL include files not found (you may need to install GSL devel libraries))])
+AC_CHECK_HEADER(gsl/gsl_linalg.h,
+                ,
+                [ AC_MSG_ERROR(error: GSL include files not found (you may need to install GSL devel libraries))])
+
+GLOBAL_LDFLAGS="$GLOBAL_LDFLAGS -lgsl -lgslcblas -lm"
+
 
-GLOBAL_CFLAGS="$GLOBAL_CFLAGS $GSL_CFLAGS"
-GLOBAL_LDFLAGS="$GLOBAL_LDFLAGS $GSL_LIBS"
 
 AC_SUBST(GSL_CFLAGS)
 AC_SUBST(GSL_LIBS)
index b031c23..b4572f6 100644 (file)
@@ -23,6 +23,9 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_linalg.h>
 
+#define PROBLEM_INCOMPLETE_NO_MATRIX 1
+#define PROBLEM_INCOMPLETE_NO_VECTOR 2
+#define PROBLEM_INCOMPLETE_NO_BOUNDS 3
 
 #define SCREEN_WIDTH  640
 #define SCREEN_HEIGHT 480
@@ -54,6 +57,15 @@ struct geometry_info
   PTR_LIST (struct triangle, triangle);
 };
 
+#define BOUNDARY_CONDITION_DIRICHLET 0
+#define BOUNDARY_CONDITION_NEUMANN   1
+
+struct boundary_condition
+{
+  int type;
+  float value;
+};
+
 #define A_11 0
 #define A_12 1
 #define A_21 2
@@ -65,6 +77,12 @@ struct problem2d
   struct jit_expr *a[4];
   struct jit_expr *a0, *f, *g;
   
+  gsl_matrix *A;
+  gsl_vector *b;
+  gsl_vector *solution;
+  
+  PTR_LIST (struct boundary_condition, boundary);
+  
   float x, y;
 };
 
@@ -73,8 +91,14 @@ struct problem2d
 struct edge_desc *edge_desc_new (struct point *, struct point *);
 struct point *point_new (float, float);
 struct triangle *triangle_new (struct point *, struct point *, struct point *);
-int build_matrix_and_vector (struct problem2d *, gsl_matrix **, gsl_vector **);
+int build_matrix_and_vector (struct problem2d *);
 struct geometry_info *geometry_load_from_matlab (const char *, const char *, const char *, const char *);
 void geometry_destroy (struct geometry_info *);
-
+int export_matrix_and_vector (struct problem2d *, FILE *);
+void display_geometry (display_t *, struct geometry_info *, float);
+struct problem2d *problem_get_sample_from_matlab (const char *);
+int problem_lock (struct problem2d *);
+int problem_is_incomplete (struct problem2d *);
+int problem_solve (struct problem2d *);
+void display_solution (display_t *, struct problem2d *, float);
 #endif /* _MAIN_APP_INCLUDE_H */
index 30a70d4..43098d5 100644 (file)
@@ -7,4 +7,4 @@ bin_PROGRAMS = 2dsolve
 
 2dsolve_LDADD = ../gsl/libgslbinding.a ../matlab/libmatlab.a ../jit/libjit.a ../sim-static/libsim.a ../util/libutil.a
 
-2dsolve_SOURCES = geometry.c main.c 2dsolve.h solver.c triangle.c triangle.h
+2dsolve_SOURCES = display.c geometry.c main.c 2dsolve.h problem.c solver.c triangle.c triangle.h
diff --git a/src/config.h b/src/config.h
new file mode 100644 (file)
index 0000000..6f54bf8
--- /dev/null
@@ -0,0 +1,122 @@
+/* src/config.h.  Generated from config.h.in by configure.  */
+/* src/config.h.in.  Generated from configure.in by autoheader.  */
+
+/* Define to 1 if you have the <dlfcn.h> header file. */
+#define HAVE_DLFCN_H 1
+
+/* Define to 1 if you have the `fork' function. */
+#define HAVE_FORK 1
+
+/* This compiler supports inline functions */
+#define HAVE_INLINE 1
+
+/* Define to 1 if you have the <inttypes.h> header file. */
+#define HAVE_INTTYPES_H 1
+
+/* Define to 1 if you have the `gsl' library (-lgsl). */
+#define HAVE_LIBGSL 1
+
+/* Define to 1 if you have the `gslcblas' library (-lgslcblas). */
+#define HAVE_LIBGSLCBLAS 1
+
+/* Define to 1 if you have the `m' library (-lm). */
+#define HAVE_LIBM 1
+
+/* Define to 1 if you have the `z' library (-lz). */
+#define HAVE_LIBZ 1
+
+/* Define to 1 if your system has a GNU libc compatible `malloc' function, and
+   to 0 otherwise. */
+#define HAVE_MALLOC 1
+
+/* Define to 1 if you have the <memory.h> header file. */
+#define HAVE_MEMORY_H 1
+
+/* Define to 1 if you have the <stdint.h> header file. */
+#define HAVE_STDINT_H 1
+
+/* Define to 1 if you have the <stdlib.h> header file. */
+#define HAVE_STDLIB_H 1
+
+/* Define to 1 if you have the <strings.h> header file. */
+#define HAVE_STRINGS_H 1
+
+/* Define to 1 if you have the <string.h> header file. */
+#define HAVE_STRING_H 1
+
+/* Define to 1 if you have the <sys/stat.h> header file. */
+#define HAVE_SYS_STAT_H 1
+
+/* Define to 1 if you have the <sys/types.h> header file. */
+#define HAVE_SYS_TYPES_H 1
+
+/* Define to 1 if you have the <unistd.h> header file. */
+#define HAVE_UNISTD_H 1
+
+/* Define to 1 if you have the `vfork' function. */
+#define HAVE_VFORK 1
+
+/* Define to 1 if you have the <vfork.h> header file. */
+/* #undef HAVE_VFORK_H */
+
+/* Define to 1 if `fork' works. */
+#define HAVE_WORKING_FORK 1
+
+/* Define to 1 if `vfork' works. */
+#define HAVE_WORKING_VFORK 1
+
+/* Define to the sub-directory in which libtool stores uninstalled libraries.
+   */
+#define LT_OBJDIR ".libs/"
+
+/* Define to 1 if your C compiler doesn't accept -c and -o together. */
+/* #undef NO_MINUS_C_MINUS_O */
+
+/* Name of package */
+#define PACKAGE "2dsolve"
+
+/* Define to the address where bug reports for this package should be sent. */
+#define PACKAGE_BUGREPORT ""
+
+/* Define to the full name of this package. */
+#define PACKAGE_NAME "2dsolve"
+
+/* Define to the full name and version of this package. */
+#define PACKAGE_STRING "2dsolve 0.1"
+
+/* Define to the one symbol short name of this package. */
+#define PACKAGE_TARNAME "2dsolve"
+
+/* Define to the home page for this package. */
+#define PACKAGE_URL ""
+
+/* Define to the version of this package. */
+#define PACKAGE_VERSION "0.1"
+
+/* Define to 1 if you have the ANSI C header files. */
+#define STDC_HEADERS 1
+
+/* Define to 1 if you can safely include both <sys/time.h> and <time.h>. */
+#define TIME_WITH_SYS_TIME 1
+
+/* Version number of package */
+#define VERSION "0.1"
+
+/* Define to 1 if `lex' declares `yytext' as a `char *' by default, not a
+   `char[]'. */
+/* #undef YYTEXT_POINTER */
+
+/* Define to `__inline__' or `__inline' if that's what the C compiler
+   calls it, or to nothing if 'inline' is not supported under any name.  */
+#ifndef __cplusplus
+/* #undef inline */
+#endif
+
+/* Define to rpl_malloc if the replacement function should be used. */
+/* #undef malloc */
+
+/* Define to `int' if <sys/types.h> does not define. */
+/* #undef pid_t */
+
+/* Define as `fork' if `vfork' does not work. */
+/* #undef vfork */
diff --git a/src/display.c b/src/display.c
new file mode 100644 (file)
index 0000000..e4e3f27
--- /dev/null
@@ -0,0 +1,238 @@
+#include <2dsolve.h>
+
+#define SX(x, zoom) ((int) (((x) * (zoom)) + SCREEN_WIDTH / 2))
+#define SY(y, zoom) ((int) (-((y) * (zoom)) + SCREEN_HEIGHT / 2))
+
+int triangle_select;
+
+/*
+int
+khandler (int key, struct display_info *disp, struct event_info *event)
+{
+  if (event->state)
+    return;
+    
+  if (event->code == SDLK_n)
+    triangle_select++;
+  else if (event->code == SDLK_p && triangle_select)
+    triangle_select--;
+  else
+    return 0;
+    
+  display_geometry (disp, curr_info, curr_zoom);
+  display_refresh (disp);
+}
+*/
+
+/* Vamos a recurrir a isométricas:
+
+                 Z
+                 ^
+                 |
+                 |
+                 |
+                 |
+                 |
+                 |
+               _-O-_
+            _-       -_
+         _-             -_
+      _-                   -_
+  X <                         > Y
+  
+  La idea será calcular las proyecciones de los tres ejes
+  sobre la pantalla, y multiplicar cada punto 3D por
+  estos ejes proyectados.
+  
+  Las coordenadas 2D de estos ejes son:
+  
+  X = (-cos (30), -sen (30)) = (-sqrt(3) / 2, -1/2)
+  Y = (cos (30), -sen (30))  = (sqrt(3) / 2, -1/2)
+  Z = (0, 1) 
+*/
+
+  
+#define SQRT_3_2  0.86603
+         
+static void
+line3d (display_t *disp, float x1, float y1, float z1,
+                         float x2, float y2, float z2,
+                         Uint32 color, float zoom)
+{
+  float p1_x;
+  float p1_y;
+  
+  float p2_x;
+  float p2_y;
+  
+  p1_x = SQRT_3_2 * (y1 - x1);
+  p1_y = - 0.5 * (x1 + y1) + z1;
+  
+  p2_x = SQRT_3_2 * (y2 - x2);
+  p2_y = - 0.5 * (x2 + y2) + z2;
+  
+  line (disp, 
+        SX (p1_x, zoom), SY (p1_y, zoom),
+        SX (p2_x, zoom), SY (p2_y, zoom),
+        color);             
+}
+
+static void
+mark3d (display_t *disp, float x1, float y1, float z1,
+                         Uint32 color, float zoom)
+{
+  float p1_x;
+  float p1_y;
+  
+  p1_x = SQRT_3_2 * (y1 - x1);
+  p1_y = - 0.5 * (x1 + y1) + z1;
+  
+  mark (disp, 
+        SX (p1_x, zoom), SY (p1_y, zoom),
+        2, color, color);             
+}
+
+void
+display_solution (display_t *disp, struct problem2d *problem, float zoom)
+{
+  int i, j;
+  
+  float x, y, z;
+  float nx, ny, nz;
+  
+  if (problem->solution == NULL)
+  {
+    disputs (disp, "Cannot paint solution: system not solved\n");
+    return;
+  }
+  
+  line3d (disp, 0, 0, 0, 1, 0, 0, OPAQUE (0xffffff), zoom);
+  line3d (disp, 0, 0, 0, 0, 1, 0, OPAQUE (0xffffff), zoom);
+  line3d (disp, 0, 0, 0, 0, 0, 1, OPAQUE (0xffffff), zoom);
+  
+  for (i = 0; i < problem->geometry->triangle_count; i++)
+  {
+    for (j = 0; j < TRIANGLE_NODES; j++)
+    {
+      x = problem->geometry->triangle_list[i]->nodes[j]->x;
+      y = problem->geometry->triangle_list[i]->nodes[j]->y;
+      z = 10.0 * gsl_vector_get (problem->solution, problem->geometry->triangle_list[i]->nodes[j]->index);
+      
+      nx = problem->geometry->triangle_list[i]->nodes[(j + 1) % TRIANGLE_NODES]->x;
+      ny = problem->geometry->triangle_list[i]->nodes[(j + 1) % TRIANGLE_NODES]->y;
+      nz = 10.0 * gsl_vector_get (problem->solution, problem->geometry->triangle_list[i]->nodes[(j + 1) % TRIANGLE_NODES]->index);
+      
+      line3d (disp, x, y, z, nx, ny, nz, OPAQUE (calc_color (400 * (z + nz))), zoom);
+      mark3d (disp, x, y, z, OPAQUE (calc_color (800 * z)), zoom);
+      // mark3d (disp, nx, ny, nz, OPAQUE (calc_color (8000 * nz)), zoom);
+    }
+  }
+}
+
+/*
+void
+display_geometry (display_t *disp, struct geometry_info *info, float zoom)
+{
+  int i;
+  int select;
+  
+  select = triangle_select % info->triangle_count;
+  
+  clear (disp, OPAQUE (0));
+  
+  display_printf (disp, 0, 0, OPAQUE (0xffa500), OPAQUE (0), "Triangle select: %5d/%d",
+    triangle_select, info->triangle_count);
+    
+  for (i = 0; i < info->triangle_count; i++)
+  {
+    line (disp, SX (info->triangle_list[i]->nodes[0]->x, zoom), 
+                SY (info->triangle_list[i]->nodes[0]->y, zoom),
+                SX (info->triangle_list[i]->nodes[1]->x, zoom), 
+                SY (info->triangle_list[i]->nodes[1]->y, zoom),
+                OPAQUE (0xff0000));
+                
+    line (disp, SX (info->triangle_list[i]->nodes[1]->x, zoom), 
+                SY (info->triangle_list[i]->nodes[1]->y, zoom),
+                SX (info->triangle_list[i]->nodes[2]->x, zoom), 
+                SY (info->triangle_list[i]->nodes[2]->y, zoom),
+                OPAQUE (0xff0000));
+                
+    line (disp, SX (info->triangle_list[i]->nodes[2]->x, zoom), 
+                SY (info->triangle_list[i]->nodes[2]->y, zoom),
+                SX (info->triangle_list[i]->nodes[0]->x, zoom), 
+                SY (info->triangle_list[i]->nodes[0]->y, zoom),
+                OPAQUE (0xff0000));
+  }
+  
+  line (disp, SX (FT_x (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
+              SY (FT_y (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
+              SX (FT_x (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
+              SY (FT_y (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
+              OPAQUE (0x00ffff));
+  
+  line (disp, SX (FT_x (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
+              SY (FT_y (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
+              SX (FT_x (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
+              SY (FT_y (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
+              OPAQUE (0x00ffff));
+  
+  line (disp, SX (FT_x (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
+              SY (FT_y (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
+              SX (FT_x (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
+              SY (FT_y (info->triangle_list[select]->nodes[0], 
+                        info->triangle_list[select]->nodes[1], 
+                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
+              OPAQUE (0x00ffff));
+  
+  for (i = 0; i < info->edge_count; i++)
+    line (disp, SX (info->edge_list[i]->p1->x, zoom), 
+                SY (info->edge_list[i]->p1->y, zoom),
+                SX (info->edge_list[i]->p2->x, zoom), 
+                SY (info->edge_list[i]->p2->y, zoom),
+                OPAQUE (0xffa500));
+  
+  for (i = 0; i < info->point_count; i++)
+    mark (disp, SX (info->point_list[i]->x, zoom),
+                SY (info->point_list[i]->y, zoom),
+                2, OPAQUE (0x00ff00), 0);
+}
+*/
+
+void
+prepare_display (display_t *disp)
+{
+  /* display_register_key_handler (disp, SDLK_n, khandler);
+  display_register_key_handler (disp, SDLK_p, khandler);*/
+  
+    /*
+  int i, j;
+  
+  for (i = 0; i < A->size1; i++)
+    for (j = 0; j < A->size2; j++)
+      pset_abs (disp, i, j, OPAQUE (calc_color (15 * gsl_matrix_get (A, i, j))));
+  */
+}
+
+
index 633640c..574409c 100644 (file)
@@ -37,8 +37,17 @@ edge_desc_new (struct point *p1, struct point *p2)
   
   edge_desc = xmalloc (sizeof (struct edge_desc));
   
-  edge_desc->p1 = p1;
-  edge_desc->p2 = p2;
+  /* Este truco nos permite asociar ejes a nodos. */
+  if (p1->index < p2->index)
+  {
+    edge_desc->p1 = p1;
+    edge_desc->p2 = p2;
+  }
+  else
+  {
+    edge_desc->p1 = p2;
+    edge_desc->p2 = p1;
+  }
   
   return edge_desc;
 }
index 202a18e..ed5ec83 100644 (file)
 
 #include <2dsolve.h>
 
-#define SX(x, zoom) ((int) (((x - 0.05) * (zoom)) + SCREEN_WIDTH / 2))
-#define SY(y, zoom) ((int) (-((y + 0.45) * (zoom)) + SCREEN_HEIGHT / 2))
-
-int triangle_select;
-struct geometry_info *curr_info;
-float curr_zoom;
-
-void display_geometry (display_t *, struct geometry_info *, float);
-
-int
-khandler (int key, struct display_info *disp, struct event_info *event)
-{
-  if (event->state)
-    return;
-    
-  if (event->code == SDLK_n)
-    triangle_select++;
-  else if (event->code == SDLK_p && triangle_select)
-    triangle_select--;
-  else
-    return 0;
-    
-  display_geometry (disp, curr_info, curr_zoom);
-  display_refresh (disp);
-}
-
-void
-display_geometry (display_t *disp, struct geometry_info *info, float zoom)
-{
-  int i;
-  int select;
-  
-  select = triangle_select % info->triangle_count;
-  
-  clear (disp, OPAQUE (0));
-  
-  display_printf (disp, 0, 0, OPAQUE (0xffa500), OPAQUE (0), "Triangle select: %5d/%d",
-    triangle_select, info->triangle_count);
-    
-  for (i = 0; i < info->triangle_count; i++)
-  {
-    line (disp, SX (info->triangle_list[i]->nodes[0]->x, zoom), 
-                SY (info->triangle_list[i]->nodes[0]->y, zoom),
-                SX (info->triangle_list[i]->nodes[1]->x, zoom), 
-                SY (info->triangle_list[i]->nodes[1]->y, zoom),
-                OPAQUE (0xff0000));
-                
-    line (disp, SX (info->triangle_list[i]->nodes[1]->x, zoom), 
-                SY (info->triangle_list[i]->nodes[1]->y, zoom),
-                SX (info->triangle_list[i]->nodes[2]->x, zoom), 
-                SY (info->triangle_list[i]->nodes[2]->y, zoom),
-                OPAQUE (0xff0000));
-                
-    line (disp, SX (info->triangle_list[i]->nodes[2]->x, zoom), 
-                SY (info->triangle_list[i]->nodes[2]->y, zoom),
-                SX (info->triangle_list[i]->nodes[0]->x, zoom), 
-                SY (info->triangle_list[i]->nodes[0]->y, zoom),
-                OPAQUE (0xff0000));
-  }
-  
-  line (disp, SX (FT_x (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
-              SY (FT_y (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
-              SX (FT_x (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
-              SY (FT_y (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
-              OPAQUE (0x00ffff));
-  
-  line (disp, SX (FT_x (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
-              SY (FT_y (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 1.0, 0.0), zoom),
-              SX (FT_x (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
-              SY (FT_y (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
-              OPAQUE (0x00ffff));
-  
-  line (disp, SX (FT_x (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
-              SY (FT_y (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 1.0), zoom),
-              SX (FT_x (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
-              SY (FT_y (info->triangle_list[select]->nodes[0], 
-                        info->triangle_list[select]->nodes[1], 
-                        info->triangle_list[select]->nodes[2], 0.0, 0.0), zoom),
-              OPAQUE (0x00ffff));
-  
-  for (i = 0; i < info->edge_count; i++)
-    line (disp, SX (info->edge_list[i]->p1->x, zoom), 
-                SY (info->edge_list[i]->p1->y, zoom),
-                SX (info->edge_list[i]->p2->x, zoom), 
-                SY (info->edge_list[i]->p2->y, zoom),
-                OPAQUE (0xffa500));
-  
-  for (i = 0; i < info->point_count; i++)
-    mark (disp, SX (info->point_list[i]->x, zoom),
-                SY (info->point_list[i]->y, zoom),
-                2, OPAQUE (0x00ff00), 0);
-}
-
-struct problem2d *
-get_default_problem (struct geometry_info *geom)
-{
-  struct problem2d *new;
-  
-  new = xmalloc (sizeof (struct problem2d));
-  
-  memset (new, 0, sizeof (struct problem2d));
-  
-  new->geometry = geom;
-  jit_compile (new->a[A_11] = jit_expr_new ("1"));
-  jit_compile (new->a[A_22] = jit_expr_new ("1"));
-  jit_compile (new->a[A_21] = jit_expr_new ("0"));
-  jit_compile (new->a[A_12] = jit_expr_new ("0"));
-  jit_compile (new->a0      = jit_expr_new ("1"));
-  jit_compile (new->f       = jit_expr_new ("9.8"));
-  
-  return new;
-}
-
 int
 main (int argc, char *argv[], char *envp[])
 {
-  int i, j;
-  
   display_t *disp;
   struct problem2d *problem;
-  gsl_matrix *A;
-  gsl_vector *b;
-  
-  curr_zoom = 200.0;
   
   if (argc != 2)
   {
@@ -164,32 +23,39 @@ main (int argc, char *argv[], char *envp[])
     exit (EXIT_FAILURE);
   }
   
-  if ((curr_info = geometry_load_from_matlab (argv[1], "p", "e", "t")) == NULL)
+  if ((problem = problem_get_sample_from_matlab (argv[1])) == NULL)
   {
-    fprintf (stderr, "%s: couldn't load geometry from `%s'\n", argv[0], argv[1]);
+    fprintf (stderr, "Cannot load problem from MATLAB\n");
     exit (EXIT_FAILURE);
-  } 
+  }
   
-  problem = get_default_problem (curr_info);
+  if (build_matrix_and_vector (problem) == -1)
+  {
+    fprintf (stderr, "Cannot load stiffness matrix\n");
+    exit (EXIT_FAILURE);
+  }
   
-  build_matrix_and_vector (problem, &A, &b);
+  problem_force_dirichlet (problem, 0.0);
   
-  if ((disp = display_new (SCREEN_WIDTH, SCREEN_HEIGHT)) == NULL)
+  problem_lock (problem);
+  
+  if (problem_solve (problem) == -1)
+  {
+    fprintf (stderr, "Solver failed\n");
     exit (EXIT_FAILURE);
+  }
   
-  display_register_key_handler (disp, SDLK_n, khandler);
-  display_register_key_handler (disp, SDLK_p, khandler);
-    
-  display_geometry (disp, curr_info, curr_zoom);
+  // export_matrix_and_vector (problem, stdout);
   
-  for (i = 0; i < A->size1; i++)
-    for (j = 0; j < A->size2; j++)
-      pset_abs (disp, i, j, OPAQUE (calc_color (15 * gsl_matrix_get (A, i, j))));
+  if ((disp = display_new (SCREEN_WIDTH, SCREEN_HEIGHT)) == NULL)
+    exit (EXIT_FAILURE);
       
+  display_solution (disp, problem, 200.0);
+  
   display_refresh (disp);
   
   display_end (disp);
-  geometry_destroy (curr_info);
+  
   return 0;
 }
 
diff --git a/src/problem.c b/src/problem.c
new file mode 100644 (file)
index 0000000..f74db00
--- /dev/null
@@ -0,0 +1,140 @@
+#include <2dsolve.h>
+
+struct boundary_condition *
+boundary_condition_new (int type, float value)
+{
+  struct boundary_condition *new;
+  
+  new = xmalloc (sizeof (struct boundary_condition));
+  
+  new->type = type;
+  new->value = value;
+  
+  return new;
+}
+
+struct problem2d *
+problem_get_default (struct geometry_info *geom)
+{
+  struct problem2d *new;
+  
+  new = xmalloc (sizeof (struct problem2d));
+  
+  memset (new, 0, sizeof (struct problem2d));
+  
+  new->geometry = geom;
+  
+  jit_compile (new->a[A_11] = jit_expr_new ("1"));
+  jit_compile (new->a[A_22] = jit_expr_new ("1"));
+  jit_compile (new->a[A_21] = jit_expr_new ("0"));
+  jit_compile (new->a[A_12] = jit_expr_new ("0"));
+  jit_compile (new->a0      = jit_expr_new ("1"));
+  jit_compile (new->f       = jit_expr_new ("-4.9"));
+  
+  new->boundary_count = geom->edge_count;
+  new->boundary_list = xmalloc (geom->edge_count * sizeof (struct boundary_condition *));
+  
+  memset (new->boundary_list, 0, geom->edge_count * sizeof (struct boundary_condition *));
+  
+  return new;
+}
+
+void
+problem_force_dirichlet (struct problem2d *problem, float value)
+{
+  int i;
+  
+  for (i = 0; i < problem->boundary_count; i++)
+    if (problem->boundary_list[i] == NULL)
+      problem->boundary_list[i] = boundary_condition_new (BOUNDARY_CONDITION_DIRICHLET, value);
+}
+
+struct problem2d *
+problem_get_sample_from_matlab (const char *path)
+{
+  struct problem2d *problem;
+  struct geometry_info *curr_info;
+  
+  if ((curr_info = geometry_load_from_matlab (path, "p", "e", "t")) == NULL)
+  {
+    fprintf (stderr, "cannot open %s: %s\n", path, strerror (errno));
+    return NULL;
+  }
+  
+  return problem_get_default (curr_info);
+}
+
+#define REALLY_HUGE 1e30
+
+int
+problem_is_incomplete (struct problem2d *problem)
+{
+  int i;
+  
+  if (problem->A == NULL)
+    return PROBLEM_INCOMPLETE_NO_MATRIX;
+    
+  if (problem->b == NULL)
+    return PROBLEM_INCOMPLETE_NO_VECTOR;
+    
+  for (i = 0; i < problem->boundary_count; i++)
+    if (problem->boundary_list[i] == NULL)
+      return PROBLEM_INCOMPLETE_NO_BOUNDS;
+      
+  return 0;
+}
+
+static void
+lock_the_right_way (struct problem2d *problem, int var, float value)
+{
+  int i;
+  
+  for (i = 0; i < problem->A->size2; i++)
+    gsl_matrix_set (problem->A, var, i, i == var);
+  
+  gsl_vector_set (problem->b, var, value);
+}
+
+int
+problem_lock (struct problem2d *problem)
+{
+  int code;
+  int i;
+  
+  code = problem_is_incomplete (problem);
+  
+  if (code == PROBLEM_INCOMPLETE_NO_MATRIX ||
+      code == PROBLEM_INCOMPLETE_NO_VECTOR)
+  {
+    fprintf (stderr, "Cannot lock boundary conditions: stiffness matrix/second member uninitialized\n");
+    return -1;
+  }
+  else if (code == PROBLEM_INCOMPLETE_NO_BOUNDS)
+  {
+    fprintf (stderr, "Cannot lock boundary conditions: some boundaries are undefined\n");
+    return -1;
+  }
+  
+  for (i = 0; i < problem->boundary_count; i++)
+  {
+    if (problem->boundary_list[i]->type != BOUNDARY_CONDITION_DIRICHLET)
+    {
+      fprintf (stderr, "Cannot understand Neumann conditions (yet)\n");
+      return -1;
+    }
+    
+    
+    lock_the_right_way (problem, problem->geometry->edge_list[i]->p1->index, problem->boundary_list[i]->value);
+    lock_the_right_way (problem, problem->geometry->edge_list[i]->p2->index, problem->boundary_list[i]->value);
+    
+    /*
+    gsl_matrix_set (problem->A, problem->geometry->edge_list[i]->p1->index, problem->geometry->edge_list[i]->p1->index, REALLY_HUGE);
+    gsl_vector_set (problem->b, problem->geometry->edge_list[i]->p1->index, REALLY_HUGE * problem->boundary_list[i]->value);
+    
+    gsl_matrix_set (problem->A, problem->geometry->edge_list[i]->p2->index, problem->geometry->edge_list[i]->p2->index, REALLY_HUGE);
+    gsl_vector_set (problem->b, problem->geometry->edge_list[i]->p2->index, REALLY_HUGE * problem->boundary_list[i]->value);
+    */
+    
+  }
+}
+
index 80cf2ad..c317f72 100644 (file)
@@ -15,7 +15,7 @@ struct polynomial
 };
 
 /* Muy feo, pero funcional */
-INLINE void
+INLINE int
 __get_basis (struct polynomial *out, struct point *center, struct point *p2, struct point *p3)
 {
   double a_data[] = { center->x, center->y, 1.0,
@@ -35,9 +35,27 @@ __get_basis (struct polynomial *out, struct point *center, struct point *p2, str
 
   gsl_permutation * p = gsl_permutation_alloc (3);
 
-  gsl_linalg_LU_decomp (&m.matrix, p, &s);
+  if (gsl_linalg_LU_decomp (&m.matrix, p, &s))
+  {
+    fprintf (stderr, "__get_basis: cannot build basic polynomial for triangle [%d], %d, %d: LU decomposition failed\n", 
+      center->index, p2->index, p3->index);
+    
+    gsl_permutation_free (p);
+    gsl_vector_free (x);
+    
+    return -1;
+  }
+  
 
-  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
+  if (gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x))
+  {
+    fprintf (stderr, "__get_basis: cannot build basic polynomial for triangle [%d], %d, %d: LU solver failed\n", 
+      center->index, p2->index, p3->index);
+    gsl_permutation_free (p);
+    gsl_vector_free (x);
+    
+    return -1;
+  }
 
   out->a = gsl_vector_get (x, 0);
   out->b = gsl_vector_get (x, 1);
@@ -45,14 +63,16 @@ __get_basis (struct polynomial *out, struct point *center, struct point *p2, str
   
   gsl_permutation_free (p);
   gsl_vector_free (x);
+  
+  return 0;
 }
 
-INLINE void
+INLINE int
 get_basis (struct polynomial *out, struct triangle *triangle, int node)
 {
-  __get_basis (out, triangle->nodes[node % 3], 
-                    triangle->nodes[(node + 1) % 3], 
-                    triangle->nodes[(node + 2) % 3]);
+  return __get_basis (out, triangle->nodes[node % 3], 
+                       triangle->nodes[(node + 1) % 3], 
+                       triangle->nodes[(node + 2) % 3]);
 }
 
 INLINE float
@@ -207,8 +227,42 @@ triangle_integrate_second (struct problem2d *problem, struct triangle *triangle,
   return integral;
 }
 
+#define xfprintf(fp, fmt, arg...)     \
+  do                                  \
+  {                                   \
+    if (fprintf (fp, fmt, ##arg) < 0) \
+      return -1;                      \
+  }                                   \
+  while (0)
+  
+int
+export_matrix_and_vector (struct problem2d *problem, FILE *fp)
+{
+  int i, j;
+  
+  xfprintf (fp, "A = [");
+  
+  for (i = 0; i < problem->A->size1; i++)
+  {
+    if (i)
+      xfprintf (fp, ";\n");
+      
+    for (j = 0; j < problem->A->size2; j++)
+      xfprintf (fp, "%s%lf", j ? ", " : "", gsl_matrix_get (problem->A, i, j));
+  }
+  
+  xfprintf (fp, "];\n");
+  
+  xfprintf (fp, "b = [");
+  
+  for (i = 0; i < problem->b->size; i++)
+    xfprintf (fp, "%s%lf", i ? "; " : "", gsl_vector_get (problem->b, i));
+  
+  xfprintf (fp, "];\n");
+}
+
 int
-build_matrix_and_vector (struct problem2d *problem, gsl_matrix **A, gsl_vector **b)
+build_matrix_and_vector (struct problem2d *problem)
 {
   int i, j, k;
   float integral;
@@ -218,10 +272,10 @@ build_matrix_and_vector (struct problem2d *problem, gsl_matrix **A, gsl_vector *
   /* Generamos una matriz de rigidez que va a tener NxN, con
      N = número de nodos. */
      
-  *A = gsl_matrix_calloc (problem->geometry->point_count,
-                          problem->geometry->point_count);
+  problem->A = gsl_matrix_calloc (problem->geometry->point_count,
+                                  problem->geometry->point_count);
                           
-  *b = gsl_vector_calloc (problem->geometry->point_count);
+  problem->b = gsl_vector_calloc (problem->geometry->point_count);
   
   for (i = 0; i < problem->geometry->triangle_count; i++)
   {
@@ -235,14 +289,15 @@ build_matrix_and_vector (struct problem2d *problem, gsl_matrix **A, gsl_vector *
        
     for (j = 0; j < TRIANGLE_NODES; j++)
     {
-      get_basis (&v_p,  problem->geometry->triangle_list[i], j);
+      if (get_basis (&v_p,  problem->geometry->triangle_list[i], j) == -1)
+        return -1;
       
       /* Podemos ir construyendo el segundo miembro: */
-      prev = gsl_vector_get (*b, problem->geometry->triangle_list[i]->nodes[j]->index);
+      prev = gsl_vector_get (problem->b, problem->geometry->triangle_list[i]->nodes[j]->index);
       
       integral = triangle_integrate_second (problem, problem->geometry->triangle_list[i], &v_p);
       
-      gsl_vector_set (*b, problem->geometry->triangle_list[i]->nodes[j]->index, prev + integral);
+      gsl_vector_set (problem->b, problem->geometry->triangle_list[i]->nodes[j]->index, prev + integral);
       
       /* Y por cada una de ellas, integraremos en todo el triángulo cuando
          uh es una función básica centrada en triangle_list[i]->nodes[k].
@@ -250,19 +305,20 @@ build_matrix_and_vector (struct problem2d *problem, gsl_matrix **A, gsl_vector *
       
       for (k = 0; k < TRIANGLE_NODES; k++)
       {
-        get_basis (&u_p,  problem->geometry->triangle_list[i], k);
+        if (get_basis (&u_p,  problem->geometry->triangle_list[i], k) == -1)
+          return -1;
         
         /* Nótese que el lugar correspondiente en la matriz de rigidez
            para esta integral es (nodes[j]->index, nodes[k]->index) */
         
-        prev = gsl_matrix_get (*A, problem->geometry->triangle_list[i]->nodes[j]->index,
+        prev = gsl_matrix_get (problem->A, problem->geometry->triangle_list[i]->nodes[j]->index,
                                problem->geometry->triangle_list[i]->nodes[k]->index);
                                
         integral = triangle_integrate (problem, problem->geometry->triangle_list[i], &u_p, &v_p);
 
-        gsl_matrix_set (*A, problem->geometry->triangle_list[i]->nodes[j]->index,
-                               problem->geometry->triangle_list[i]->nodes[k]->index,
-                               prev + integral);
+        gsl_matrix_set (problem->A, problem->geometry->triangle_list[i]->nodes[j]->index,
+                        problem->geometry->triangle_list[i]->nodes[k]->index,
+                        prev + integral);
       }
     }
   }
@@ -270,4 +326,36 @@ build_matrix_and_vector (struct problem2d *problem, gsl_matrix **A, gsl_vector *
   return 0;
 }
 
+/* TODO: añadir estados */
+int
+problem_solve (struct problem2d *problem)
+{ 
+  gsl_permutation * p;
+  problem->solution = gsl_vector_alloc (problem->b->size);
+  int s;
+  
+  p = gsl_permutation_alloc (problem->A->size1);
+
+  if (gsl_linalg_LU_decomp (problem->A, p, &s))
+  {
+    fprintf (stderr, "problem_solve: cannot perform LU decomposition!\n");
+
+    gsl_permutation_free (p);    
+    return -1;
+  }
+  
+  if (gsl_linalg_LU_solve (problem->A, p, problem->b, problem->solution))
+  {
+    fprintf (stderr, "problem_solve: cannot solve via LU!\n");
+    gsl_permutation_free (p);
+    
+    return -1;
+  }
+
+  gsl_permutation_free (p);
+
+  
+  return 0;
+}
+