1 /* pdnmesh - a 2D finite element solver 2 * 3 Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net> 4 This program is free software; you can redistribute it and/or modify 5 it under the terms of the GNU General Public License as published by 6 the Free Software Foundation; either version 2 of the License, or 7 (at your option) any later version. 8 9 This program is distributed in the hope that it will be useful, 10 but WITHOUT ANY WARRANTY; without even the implied warranty of 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 GNU General Public License for more details. 13 14 You should have received a copy of the GNU General Public License 15 along with this program; if not, write to the Free Software 16 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 17 $Id: types.h,v 1.81 2005/04/24 05:32:24 sarod Exp $ 18 */ 19 20 #ifndef TYPES_H 21 #define TYPES_H 1 22 #ifdef __cplusplus 23 extern "C" { 24 #endif 25 #ifndef WIN32 26 #include <gtk/gtk.h> 27 #include <gdk/gdkkeysyms.h> 28 #include <gtk/gtkgl.h> 29 #endif /* WIN32 */ 30 #ifdef WIN32 31 #include <windows.h> 32 #define PACKAGE "pdnmesh" 33 #define VERSION "0.2.1" 34 //#define USE_MKL /* only set when using intel MKL for windows */ 35 #define PACKAGE_BUGREPORT "pdnmesh-bugs" 36 #define YY_NO_UNISTD_H /* do not include unistd.h */ 37 #ifndef isblank 38 #define isblank(x) (((x)<33)||((x)>126)) 39 #endif/*!isblank*/ 40 #endif/*!WIN32*/ 41 #include <GL/gl.h> 42 #include <GL/glu.h> 43 44 /* need this to remove a warning */ 45 extern float sqrtf(float x); 46 47 #define ONE_OVER_TWELVE 0.0833333333333333333333f 48 #define ONE_OVER_SIX 0.166666666666666666667f 49 #define ONE_OVER_THREE 0.33333333333333333333f 50 51 #ifndef M_PI 52 #define M_PI 3.14159265358979323846 /* pi */ 53 #endif /* M_PI */ 54 55 #define TOL 1e-9 56 /* 57 * bitree.c 58 */ 59 #define RESIZE 100 /* number to resize the mesh by realloc() */ 60 61 typedef enum { RED, BLACK } nodecolour; /* red-black tree */ 62 typedef enum { VAR, /* varible point */ 63 FX, /* fixed */ 64 P_UNKNOWN /* unknown - initial imported point */ 65 } pointvalue; /* point type */ 66 typedef enum { NU, /* Neumann edge */ 67 DR, /* Dirichlet boundary */ 68 E_UNKNOWN } etype; 69 /* coordinate precision */ 70 #define MY_DOUBLE double 71 #ifndef MY_DOUBLE 72 #define MY_DOUBLE float 73 #endif 74 /* format for MY_DOUBLE */ 75 #define MDF "%lf" 76 #ifndef MY_INT 77 #define MY_INT int 78 #endif 79 #define MIF "%d" 80 81 #define DEFAULT_TITLE "pdnmesh" 82 #define DEFAULT_WIDTH 400 83 #define DEFAULT_HEIGHT 400 84 85 /* for the point tree */ 86 typedef struct Node_x_ { 87 struct Node_x_ * left; 88 struct Node_x_ * right; 89 struct Node_x_ * parent; 90 struct Node_y_ * Ytree; 91 nodecolour colour; 92 93 /* user data */ 94 MY_DOUBLE x; /* X coordinate */ 95 } node_x; 96 97 typedef struct Node_y_ { 98 struct Node_y_ * left; 99 struct Node_y_ * right; 100 struct Node_y_ * parent; 101 102 unsigned int n; /* point number */ 103 nodecolour colour; 104 105 /* user data */ 106 MY_DOUBLE *xp; /* pointer to X coordinate */ 107 MY_DOUBLE y; /* Y coordinate */ 108 MY_DOUBLE *z; /* potential, etc, etc. array, length=degree_of_freedom*/ 109 pointvalue val; /* VAR or FX */ 110 } node_y; 111 112 typedef struct Mesh_{ 113 node_x * Xtree; 114 node_y ** Narray; /* we keep an array for fast access */ 115 unsigned int count; 116 unsigned int maxpoints; /* current size of Narray */ 117 } Mesh; 118 119 typedef struct point_ { 120 MY_DOUBLE x; 121 MY_DOUBLE y; 122 MY_DOUBLE *z; /* potential, etc. etc. array, length=degree_of_freedom*/ 123 /* any other problem data */ 124 pointvalue val; /* fixed or variable? */ 125 } point; 126 127 128 129 /* 130 * parser.y 131 * 132 */ 133 /* data structures for the parsing part */ 134 /* node types of the syntax tree */ 135 typedef enum { CCONST, VARI, OPR } symtype; 136 /* symbol types */ 137 138 /* constants - double */ 139 typedef struct { 140 symtype type; /* type */ 141 MY_DOUBLE value; /* value */ 142 } node_const; 143 144 /* variables - lower case letters */ 145 typedef struct { 146 symtype type; /* type */ 147 int d; /* identifier (ASCII) of variable to point to array symval[] */ 148 } node_var; 149 150 /* operators */ 151 typedef struct { 152 symtype type; /* type */ 153 int operatr; /* type of operator */ 154 int nops; /* no of operands */ 155 union exp_nodetype_* operand[1]; /* array of operands - similar nodes */ 156 } node_opr; 157 158 typedef union exp_nodetype_ { 159 symtype type; /* either const, var or operand */ 160 node_const constnt; 161 node_var var; 162 node_opr opr; 163 } exp_nodetype; 164 165 166 /* function to evaluate parse treee */ 167 /* input are triangle coordinates and root pf parse tree */ 168 extern MY_DOUBLE 169 ex( exp_nodetype *p, MY_INT p1 ); 170 /* free the parse tree */ 171 extern void 172 free_parse_tree( exp_nodetype *p); 173 /* 174 * * lexer.l 175 * * 176 * */ 177 /* parse equation in string */ 178 extern int 179 equation_parse(char *string, exp_nodetype **bdata); 180 181 182 /* 183 * poly.c 184 */ 185 186 /*******************/ 187 /* edge */ 188 typedef struct edge_ { 189 int p1; /* point 1 */ 190 int p2; /* point 2 */ 191 /* more data */ 192 etype type; /* Dirichlet or Neumann edge */ 193 } poly_edge; 194 195 /* doubly connected edge list for polygon */ 196 typedef struct elist_{ 197 struct elist_* next; 198 struct elist_* prev; 199 poly_edge data; 200 } poly_elist; 201 202 203 /* polygon */ 204 typedef struct poly_{ 205 poly_elist* head; 206 int nedges; /* no of edges */ 207 /* more data */ 208 int hollow; /* 0: hollow boundary- remove triangles 1: keep triangles */ 209 MY_DOUBLE rho; /* current density */ 210 /* we can use an expression for rho as well */ 211 /* however if rhoexp=0 we use the static value */ 212 exp_nodetype *rhoexp; /* expression for rho */ 213 char *rhostr; /* remember the string expression for rho as well */ 214 215 MY_DOUBLE mu; /* permeability */ 216 exp_nodetype *muexp; /* parsed expression */ 217 char *mustr; /* string */ 218 219 MY_DOUBLE epsilon; /* permittivity (we keep this separate from mu) */ 220 exp_nodetype *epsexp; /* parsed expression */ 221 char *epsstr; /* string */ 222 } polygon; 223 224 /* all boundaries */ 225 /* boundaries are closed polygons */ 226 /* we can have open boundaries, but for them, 227 * only the edges exist in edge array */ 228 typedef struct bound_{ 229 polygon *poly_array; /* array of polygons */ 230 int npoly; /* number of polygons */ 231 } boundary; 232 233 /* the following are for dxf.c */ 234 235 typedef struct dxf_int_{ 236 int n; 237 struct dxf_int_ *next; 238 struct dxf_int_ *prev; 239 } dxf_int; 240 241 /* dxf polygon: edges as integers */ 242 typedef struct dxf_poly_{ 243 dxf_int* head; /* doubly linked list,sorted in number order */ 244 int nedges; /* no of edges */ 245 /* more data */ 246 int hollow; /* 0: hollow boundary- remove triangles, 1: keep triangles */ 247 /* strings for mu, epsilon, rho */ 248 char *mu; /* permeability */ 249 char *eps; /* permeability */ 250 char *rho; /* current density */ 251 /* we can use an expression for rho as well */ 252 } dxf_polygon; 253 254 255 /* list of polygons */ 256 typedef struct poly_list_{ 257 dxf_polygon *data; 258 struct poly_list_ *next; 259 struct poly_list_ *prev; 260 int number; /* 0 is the outer boundary */ 261 } poly_list; 262 263 /* boundary as a list of polygons */ 264 typedef struct bound_list_{ 265 poly_list * head; 266 int count; 267 } boundary_list; 268 269 /* 270 * dag.c 271 */ 272 #define DAG_STATUS_OK 1 273 274 typedef struct dag_node_ { 275 struct dag_node_list_* list; 276 void *rec; /* data record */ 277 } DAG_node; 278 279 typedef struct dag_node_list_ { 280 struct dag_node_list_ *next; 281 struct dag_node_ *child; 282 } DAG_node_list; 283 284 typedef struct dag_node_queue_ { 285 struct dag_node_queue_ *next; 286 struct dag_node_queue_ *prev; 287 struct dag_node_ *child; 288 } DAG_node_queue; 289 290 typedef struct dag_tree_ { 291 /* function for testing inclusion */ 292 int (*comp_INC)(const void *rec1,const void *rec2, const void *Data); 293 /* function for printing record data */ 294 void (*print_record)(const void *rec); 295 296 297 /* root node */ 298 DAG_node *root; 299 300 /* the following used for the queue of leaves */ 301 /* head of linked list of leaf nodes */ 302 DAG_node_queue *qhead; 303 /* auxilliary poiter used to traverse the list */ 304 DAG_node_queue *qtail; 305 /* sentinel in the queue */ 306 DAG_node_queue *sentinel; 307 } DAG_tree; 308 309 typedef enum {LIVE, HIDDEN,DEAD} triangle_status; 310 /* LIVE - leaf in dag, active triangle 311 HIDDEN - leaf in dag, but belongs to outer plane, not plotted 312 DEAD - internal node in dag 313 */ 314 315 typedef struct triag_ { 316 MY_INT p0,p1,p2; /* points in ccw order */ 317 MY_INT n; /* unique triangle number */ 318 DAG_node *dag_p; /* pointer to host node in DAG */ 319 /* neighbouring triangles */ 320 int t0,t1,t2; /* t0->p0, t1->p2, t2->p2 */ 321 triangle_status status; /* status flag to indicate triangle still in mesh */ 322 /* boundary */ 323 int boundary; 324 } triangle; 325 326 327 /* 328 * rbt.c 329 */ 330 331 332 /* colors */ 333 typedef enum { rbtBLACK, rbtRED } RBT_node_color; 334 335 /* error codes */ 336 #define RBT_STATUS_OK 0 337 #define RBT_STATUS_KEY_FOUND 1 338 #define RBT_STATUS_KEY_NOT_FOUND 2 339 340 /* generic tree node */ 341 typedef struct node_ { 342 struct node_* left; 343 struct node_* right; 344 struct node_* parent; 345 RBT_node_color color; 346 void *rec; /* record for data, including the rec for order */ 347 } RBT_node_type; 348 349 350 351 /* generic tree */ 352 typedef struct tree_ { 353 /* function for comparison of recs == */ 354 int (*comp_EQ) (const void *rec1, const void *rec2); 355 /* function for comparison of recs < */ 356 int (*comp_LT) (const void *rec1, const void *rec2); 357 /* function to print node data */ 358 void (*print_record)(const void *rec); 359 /* function to free memory used by records */ 360 void (*free_record)(void *rec); 361 362 363 /* root node */ 364 RBT_node_type *root; 365 } RBT_tree; 366 367 368 369 370 /* 371 * bitree.c 372 */ 373 /* (re) initialise the mesh */ 374 extern void 375 BIT_init(Mesh *Ms, int maxpoints); 376 377 /* the following will insert the given point into the mesh if is already is */ 378 /* not there, else it will do nothing but in both cases, the point number will be returned */ 379 extern unsigned int 380 BIT_insert(Mesh *Ms, point p); 381 382 /* free the mesh */ 383 void 384 BIT_free(Mesh *Ms); 385 386 extern Mesh M; 387 extern int maxpoints; 388 /* macros for accessing nodes */ 389 /* first to get x cordinate */ 390 #define Mx(i,M) *(((M).Narray[(i)])->xp) 391 /* to get y coordinate */ 392 #define My(i,M) ((M).Narray[(i)])->y 393 /* to get potential */ 394 #define Mz(i,M) ((M).Narray[(i)])->z[0] 395 /* to get any value of z */ 396 #define Mzz(i,j,M) ((M).Narray[(i)])->z[(j)] 397 398 /* get fixed or variable property */ 399 #define Mval(i,M) ((M).Narray[(i)])->val 400 /* to get position relative to a boundary */ 401 #define Mpos(i,M) ((M).Narray[(i)])->pos 402 403 404 #define DETERM(x1,y1,x2,y2,x3,y3) \ 405 ((x1)*(y2-y3)+(x2)*(y3-y1)+(x3)*(y1-y2)) 406 407 #ifndef ABS 408 #define ABS(x) \ 409 ((x)>=0? (x):-(x)) 410 #endif 411 412 #ifndef LENGTH 413 #define LENGTH(p1,p2,M) \ 414 ((Mx(p1,M)-Mx(p2,M))*(Mx(p1,M)-Mx(p2,M))+(My(p1,M)-My(p2,M))*(My(p1,M)-My(p2,M))) 415 #endif 416 417 #ifndef LENGTH1 418 #define LENGTH1(p,x,y,M) \ 419 ((Mx(p,M)-(x))*(Mx(p,M)-(x))+(My(p,M)-(y))*(My(p,M)-(y))) 420 #endif 421 422 423 #define IS_ZERO(x) \ 424 (ABS(x) <= TOL) 425 426 /* pdnmesh.c */ 427 428 extern int degree_of_freedom; /* degree of freedom for a point */ 429 430 extern int requested_degree_of_freedom; /* user request, for eigenvalue problems */ 431 432 extern int max_iteration_limit; /* limit for iterative solution */ 433 434 /* input cord file name */ 435 extern char * input_cord_filename; 436 437 extern int 438 insert_point_to_mesh(point p); 439 /* generate mesh */ 440 extern void 441 generate_mesh(const char *filename); 442 /* initial buildup of mesh */ 443 extern int 444 solve_problem_from_scratch(const char *filename); 445 446 /* free all allocated variables */ 447 extern void 448 free_everything(void); 449 450 extern void 451 initialize_global_variables(int argc, char *argv[]); 452 453 454 455 /*********************/ 456 /* 457 * input.c 458 */ 459 /* globals */ 460 extern poly_edge *edge_array; 461 extern int nedges; 462 extern boundary bound; /* boundaries */ 463 extern MY_DOUBLE g_xoff,g_yoff; /* x,y offsets */ 464 extern MY_DOUBLE g_xscale,g_yscale; /* x,y scales */ 465 extern MY_DOUBLE g_xrange,g_yrange; /* x,y ranges, used for contour plot */ 466 467 468 /* read input file */ 469 extern int 470 read_input_file(point **input_points, const char *infile); 471 472 473 /* 474 * poly.c 475 */ 476 477 /* poly.c */ 478 /* initialize a polygon */ 479 extern void 480 init_poly(polygon *poly); 481 /* insert an edge to a polygon */ 482 extern int 483 insert_edge_to_poly(polygon *poly,poly_edge *e); 484 /* check point is inside the polygon */ 485 /* returns 1 if true */ 486 /* note the points, and the edges has to be in ccw order */ 487 extern int 488 inside_poly(polygon *poly,int p); 489 /* check if point p is inside polygon poly */ 490 int 491 inside_poly_point(polygon *poly,point p); 492 /* sanity check polygon */ 493 /* arrange edges in correct adjacent order */ 494 extern void 495 arrange_polgon_edges(polygon *poly); 496 497 /* function to remove intersection of the triangles 498 with polygon edges */ 499 extern int 500 remove_polygon_mesh_intersections(polygon *poly); 501 502 /* copies polygon a to b, returns edges copied */ 503 extern int 504 copy_polygon(polygon *a,polygon *b); 505 506 /* destroy a polygon */ 507 void 508 destroy_polygon(polygon *a); 509 510 511 /* routines to handle a list of polygons */ 512 extern void 513 init_boundary_list(boundary_list *b); 514 extern void 515 destroy_boundary_list (boundary_list * b); 516 517 /* insert a polygon */ 518 extern int 519 insert_poly_to_boundary(boundary_list *b, dxf_polygon *p); 520 /* remove a polygon from list */ 521 extern void 522 remove_poly_from_boundary(boundary_list *b, int number); 523 /* get polygon from list, given by number */ 524 /* returns 0 if not found */ 525 extern dxf_polygon * 526 lookup_poly_from_boundary(boundary_list *b, int number); 527 528 /* destroy a polygon */ 529 extern void 530 destroy_dxf_polygon (dxf_polygon * a); 531 532 /* insert edge if it is not included in the polygon 533 * else remove it from the polygon */ 534 extern int 535 insert_remove_from_polygon(dxf_polygon *p,int edge_number); 536 /* sort all edges in polygons of the 537 * boundaries so that they are adjacent */ 538 extern int 539 sort_edges_in_polygons(boundary_list *b); 540 541 /* 542 * dag.c 543 */ 544 545 extern int 546 DAG_init(DAG_tree *t, int (*comp_INC)(const void *rec1,const void *rec2 ,const void* MM/* Mesh */), 547 void (*print_record)(const void *rec)); 548 /* we return the dag node with piggyback data record */ 549 /* because we need it for cross referencing */ 550 extern void * 551 DAG_insert(DAG_tree *t, DAG_node *parent, void *rec); 552 /* returns node pointer where data is found */ 553 /* assertion - nodes at same level cover disjoint regions */ 554 /* nodes at higher level cover same region as nodes at lower level */ 555 extern void * 556 DAG_find(DAG_tree *t, void *datarec, const void *MM/*Mesh*/); 557 /* traverse and print */ 558 extern void 559 DAG_traverse(DAG_tree *t); 560 /* add a cross link between two nodes */ 561 extern void * 562 DAG_link(DAG_tree *t, DAG_node *parent, DAG_node *child); 563 /* traverse the leaf list and prune it */ 564 /* that is, remove nodes that are not leaves */ 565 /* returns the next leaf node data pointer */ 566 extern void * 567 DAG_traverse_prune_list(DAG_tree *t); 568 569 /* reset queue for a new traversal */ 570 extern void 571 DAG_traverse_list_reset(DAG_tree *t); 572 573 /* BIG NOTE: we do not free the record here 574 * * it has to be freed using the RBT */ 575 /* free DAG */ 576 void 577 DAG_free(DAG_tree *t); 578 579 580 /* these are for reverse traversal, from tail to head */ 581 extern void * 582 DAG_reverse_traverse_prune_list(DAG_tree *t); 583 584 extern void 585 DAG_reverse_traverse_list_reset(DAG_tree *t); 586 587 /* 588 * graphics.c 589 */ 590 /* globals needed */ 591 extern double zoom_x1,zoom_x2,zoom_y1,zoom_y2; 592 /* current transformation matrix */ 593 extern GLdouble current_mat[16]; 594 extern int ntriangles,mouse_responce_flag,windowid; 595 extern DAG_tree dt; 596 extern triangle *chosen_triangle; 597 598 #ifndef WIN32 599 extern GtkWidget *global_status_bar,*global_status_label; 600 #endif/* WIN32 */ 601 602 603 /* global to update status bar */ 604 extern void 605 update_status_bar(char *msg); 606 607 /* global rendering routine */ 608 extern void 609 global_render(void); 610 611 /* option flags for plotting */ 612 extern int plot_mesh, plot_cont, plot_fill, plot_grad, plot_legend; 613 /* flag to detect existing mesh: 1 is mesh exists, 0 otherwise */ 614 extern int mesh_already_present; 615 616 /* menu items */ 617 #define MENU_0 0 618 #define MENU_FULL 1 619 #define MENU_ORIG 2 620 #define MENU_EXIT 3 621 #define MENU_INFO 4 622 #define MENU_ZOOM_START 5 623 #define MENU_ZOOM_END 6 624 #define MENU_ZOOM_ALL 7 625 #define MENU_ZOOM_BACK 8 626 #define MENU_SPLIT 9 627 /* the following to modify point/edge properties */ 628 #define MENU_EDIT_POINT 10 629 #define MENU_EDIT_EDGE 11 630 #define MENU_ADD_EDGE 12 631 /* followind to spint triangles */ 632 #define MENU_SPLIT_TRIANGLE 13 633 634 /* flags for plotting */ 635 #define PLOT_MESH 0 636 #define PLOT_CONT 1 637 #define PLOT_GRAD 2 638 #define PLOT_FILL 3 639 640 /* flags for drawing modes */ 641 #define DRAW_OUTPUT 0 642 #define DRAW_INPUT 1 643 #define DRAW_DXF 3 644 645 #ifndef WIN32 646 /* macro to update gui under long computations */ 647 #define UPDATE_GUI\ 648 while(gtk_events_pending())\ 649 gtk_main_iteration(); 650 /* 651 * graphics.c 652 */ 653 /* main window creation */ 654 extern GtkWidget* 655 create_window(GdkGLConfig *glconfig); 656 657 /* openGL framebuffer configuratiion */ 658 extern GdkGLConfig * 659 configure_gl(void); 660 #endif/* WIN32 */ 661 662 /* 663 * rbt.c 664 */ 665 666 extern RBT_tree rt; 667 /* intialize tree */ 668 /* returns RBT_STATUS_OK if no error */ 669 extern int 670 RBT_init(RBT_tree *t, int (*comp_EQ) (const void *rec1, const void *rec2), 671 int (*comp_LT) (const void *rec1, const void *rec2), 672 void (*print_record) (const void *rec), 673 void (*free_record)(void *rec)); 674 675 /* insert */ 676 /* returns RBT_STATUS_OK if inserted, RBT_STATUS_KEY_FOUND if key exists */ 677 extern int 678 RBT_insert(void *rec, RBT_tree *t); 679 /* delete */ 680 /* returns RBT_STATUS_OK if key deleted, RBT_STATUS_KEY_NOT_FOUND if not found */ 681 extern int 682 RBT_delete(void * rec,RBT_tree *t); 683 /* search */ 684 /* returns RBT_STATUS_KEY_FOUND if found, RBT_STATUS_KEY_NOT_FOUND if key not found*/ 685 extern void * 686 RBT_find(void *rec, RBT_tree *t); 687 /* traverse */ 688 extern void 689 RBT_traverse(RBT_tree *t); 690 691 /* free memory used by a tree */ 692 /* makes t=NIL */ 693 void 694 RBT_free(RBT_tree *t); 695 696 697 /* 698 * subdivide.c 699 */ 700 /* subdivision of a triangle */ 701 extern int 702 subdivide_this_triangle(triangle *tg); 703 704 /* subdivides all live triangles */ 705 extern int 706 subdivide_all_triangles(void); 707 708 extern int 709 refine_mesh_with_iterations(MY_DOUBLE priority); 710 extern int 711 refine_this_triangle(triangle *tg, int override, MY_DOUBLE priority); 712 713 /* parameters for modifying the mesh */ 714 extern MY_DOUBLE g_badness_limit; 715 extern MY_DOUBLE g_area_floor; 716 extern MY_DOUBLE g_area_ceil; 717 extern MY_DOUBLE g_gradient_limit; 718 719 /* 720 * solve.c 721 */ 722 typedef enum { POISSON, POISSON_SPARSE, HELMHOLTZ, HELMHOLTZ_INHOMO, HELMHOLTZ_FREQ, HELMHOLTZ_BETA } eq_type; /* type of equation to solve */ 723 /* HELMHOLTZ: homogeneous wave eqn. HELMHOLTZ_INHOMO: inhomogeneous 724 HELMHOLTZ_FREQ: given beta, find cutoff 725 HELMHOLTZ_BETA: given freq. find beta */ 726 extern MY_DOUBLE g_maxpot,g_minpot; 727 728 /* flag to use ARPACK solver */ 729 extern int use_arpack_solver; 730 /*extern MY_INT *renumber,*re_renumber; */ 731 732 extern int 733 re_number_and_solve_helmholtz(void); 734 extern int 735 re_number_and_solve_helmholtz_sparse(Mesh Ms, DAG_tree Dt); 736 737 extern int 738 re_number_and_solve_poisson(void); 739 740 741 /* 742 * contour.c 743 */ 744 /* function to plot the contour */ 745 extern void 746 plot_contour_all(Mesh Ms); 747 extern void 748 plot_contour_all_in_3d(Mesh Ms); 749 /* printing contour plot */ 750 void 751 print_contour_all(const char* filename,int color_flag,Mesh Ms); 752 /* printing mesh as eps*/ 753 extern void 754 print_mesh_eps(const char* filename,Mesh MM); 755 /* printing mesh as an ASCII file */ 756 extern void 757 print_mesh_ascii(const char* filename,Mesh MM); 758 /* printing potentials as ACII file*/ 759 extern void 760 print_potential(const char* filename, Mesh MM); 761 /* plotting gradient */ 762 extern void 763 plot_gradient(Mesh MM); 764 /* plotting legend */ 765 extern void 766 display_legend_window(void); 767 /* a function to print the gradient of potential */ 768 extern void 769 print_gradient(const char *filename,Mesh Ms); 770 771 /* to get a color to plot */ 772 extern int 773 get_colour(float* redp, float* greenp, float *bluep, int level, int max_levels); 774 775 extern int current_plotting_contour; /* number in z array */ 776 extern int contour_levels; /* no of contour levels */ 777 /* 778 * eig.c 779 */ 780 /* inverse iteration for finding eigenvectors 781 * given the eigenvalues */ 782 #define EIGENVECTOR_ITERATION_LIMIT 15 783 /* type of equation POISSON, HELMHOLTZ */ 784 extern eq_type solve_equation; 785 extern MY_DOUBLE *eig_array; /* array to store solved eigenvalues */ 786 /* find the eigenvector(s) given (an) eigenvalue(s) */ 787 extern void 788 find_eigenvectors(MY_DOUBLE **a, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT m, MY_INT n); 789 790 791 /* 792 * dxf.c 793 */ 794 /* read DXF file */ 795 extern int 796 read_dxf_file_and_triangulate(const char *filename); 797 798 #ifndef WIN32 799 /* display edit window */ 800 extern void 801 display_edit_window(const char* title); 802 803 /* generic dialog function */ 804 extern void 805 generic_dialog(const gchar *msg,const gchar *title); 806 #endif/*WIN32*/ 807 /* array of poins in the DXF file */ 808 extern point * dxf_point_array; 809 extern int dxf_points; 810 811 812 /* 813 * output.c 814 */ 815 /* output choice */ 816 typedef enum {OUT_CONT, OUT_CONT_COLOR, OUT_CONT_COLOR_LEGEND, OUT_GRAD, OUT_MESH, OUT_MESH_EPS, OUT_POTENTIAL} outfile_type; 817 extern outfile_type output_file_option; 818 819 extern void 820 generate_output_files(void); 821 822 823 /* 824 * wave.c 825 */ 826 /* cutoff frequency and propagation constant used 827 in waveguide problems */ 828 extern MY_DOUBLE global_wave_k0, global_wave_beta; 829 /* if 1 use full matrix solver, else use half the matrix, for symmetric 830 positive definite eigen equation */ 831 extern int use_full_eigensolver; 832 833 extern int 834 buildup_hash_table_and_equation(Mesh Ms, DAG_tree t); 835 836 /* 837 * integrate.c 838 */ 839 /* the following function was sent by Werner Hoch */ 840 /* calc the integral seperate for each dirichlet edge */ 841 /* input int con - eigenmode of potential */ 842 extern void 843 calc_integral(int con); 844 845 846 /* 847 * glist.c 848 */ 849 850 /* generic list structure */ 851 typedef struct glist_node_ { 852 struct glist_node_ *next; 853 struct glist_node_ *prev; 854 void *data; 855 } glist_node; 856 857 typedef struct glist_ { 858 glist_node *head; 859 glist_node *tail; 860 glist_node *iter; /* used to iterate the list */ 861 /* destructor function for each list data */ 862 void (*list_data_free)(void *); 863 /* function to compare two data records */ 864 int (*list_data_comp)(const void *d1, const void *d2); 865 /* data size */ 866 size_t datasize; 867 } glist; 868 869 /* intialize a generic doubly linked list */ 870 /* datasize: size of list data 871 * list_data_free: function to free list data 872 */ 873 extern void 874 glist_init(glist *L,size_t datasize, void (*list_data_free)(void *)); 875 /* initialize a generic sorted doubly linked list */ 876 /* datasize: size of list data 877 * list_data_free: function to free list data 878 * list_data_comp: compare two data items: return 0 if match, else non zero 879 */ 880 extern void 881 glist_sorted_init(glist *L,size_t datasize, void (*list_data_free)(void *), int (*list_data_comp)(const void *d1, const void *d2)); 882 883 /* insert data to list */ 884 extern void * 885 glist_insert(glist *L, const void *data ); 886 /* insert data to a sorted list */ 887 extern void * 888 glist_sorted_insert(glist *L, const void *data); 889 /* delete a complete list */ 890 extern void 891 glist_delete(glist *L); 892 /* initialize list traversalt from head to tail */ 893 extern void 894 glist_set_iter_forward(glist *L); 895 /* traverse from head to tail */ 896 extern void * 897 glist_iter_forward(glist *L); 898 /* initialize list traversal from tail to head */ 899 extern void 900 glist_set_iter_backward(glist *L); 901 /* traverse list from tail to head */ 902 extern void * 903 glist_iter_backward(glist *L); 904 /* check if the list is empty */ 905 /* return 1 if empty, 0 if not */ 906 extern int 907 glist_empty(glist *L); 908 /* remove first element from the list (at head)*/ 909 /* and return its address */ 910 extern void * 911 glist_remove(glist *L); 912 /* remove element from the list if it exists */ 913 /* returns 0 if not deleted (not exist), else 1*/ 914 extern int 915 glist_sorted_remove(glist *L, const void *data); 916 /* look for element from the list if it exists */ 917 /* returns 0 if not found (not exist), else the address*/ 918 extern void * 919 glist_sorted_lookup(glist *L, const void *data); 920 921 922 /* 923 * spmat.c 924 */ 925 /* sparse matrix routines */ 926 typedef struct SPMat_ { 927 int N; /* size is N by N */ 928 RBT_tree **rt; /* array of Balanced Trees for each row */ 929 glist **list; /* array of linked lists to access each row 930 element serially */ 931 unsigned int T; /* no. of non zero elements */ 932 } SPMat; 933 934 /* intialize a matrix */ 935 extern int 936 SPM_init( SPMat *sM, int N ); /* N - N by N sparse matrix M */ 937 938 /* add z to i,j location of M */ 939 /* initially assumed zero */ 940 /* i,j has to be within [0,N-1] */ 941 extern int 942 SPM_add( SPMat *sM, int i, int j, MY_DOUBLE z); 943 944 /* myltiply i,j location of M by z */ 945 /* initially assumed zero */ 946 /* i,j has to be within [0,N-1] */ 947 extern int 948 SPM_mult( SPMat *sM, int i, int j, MY_DOUBLE z); 949 950 /* returns size of Matrix, or value of N */ 951 extern int 952 SPM_dim(SPMat *sM); 953 /* returns filled size of M, or value of N */ 954 955 extern unsigned long int 956 SPM_size(SPMat *sM); 957 958 /* delete matrix M */ 959 extern int 960 SPM_destroy(SPMat *sM); 961 962 963 /* copy row i of M to vector v */ 964 /* v size N by 1 */ 965 extern int 966 SPM_row(SPMat *sM, MY_INT i, MY_DOUBLE *v); 967 968 /* copy column i of M to vector v */ 969 /* v size N by 1 */ 970 /* note column access is more expensive than row access */ 971 extern int 972 SPM_col(SPMat *sM, MY_INT i, MY_DOUBLE *v); 973 974 /* returns the value at location (i,j) 975 of the matrix or zero if none found */ 976 extern MY_DOUBLE 977 SPM_e(SPMat *sM, MY_INT i, MY_INT j); 978 979 /* replaces the value at location (i,j) 980 of the matrix with z */ 981 extern int 982 SPM_eq(SPMat *sM, MY_INT i, MY_INT j,MY_DOUBLE z); 983 984 /* Matrix vector product A.x*/ 985 /* A-sparse matrix*/ 986 /* x- vector*/ 987 /* v- result = Ax */ 988 /* Note: call SPM_build_list() before this is called */ 989 extern int 990 SPM_vec(SPMat *sM, MY_DOUBLE *x, MY_DOUBLE *v); 991 992 /* copy matrix A to matrix B */ 993 /* B should be properly intialized to same size as A 994 before copying. if no error, 0 returned. else -1 returned */ 995 extern int 996 SPM_copy(SPMat *A, SPMat *B); 997 998 /* add matrix A multiplied by scalar c to matrix B */ 999 /* B = B + A x c */ 1000 extern int 1001 SPM_add_scalar(SPMat *A, SPMat *B, MY_DOUBLE c); 1002 1003 /* build linked lists of each row elements for serial access */ 1004 /* returns 0 if no error */ 1005 extern int 1006 SPM_build_lists(SPMat *sM); 1007 1008 /* build compressed sparse row storage format arrays 1009 from the given sparse matrix. 1010 AA - size Nz x 1 : array of non zero elements 1011 AJ - size Nz x 1: array of column indices 1012 AI - size N x 1: array of indices for starting point of rows 1013 returns the final filled position of array AA 1014 */ 1015 extern int 1016 SPM_build_CSR(SPMat *A, double *AA, int* JA, int*IA); 1017 1018 /* create a transpose of a matrix 1019 by visiting the lists */ 1020 /* B = A^T */ 1021 extern int 1022 SPM_transpose(SPMat *A, SPMat *B); 1023 1024 /* 1025 * sparse_solve.c 1026 */ 1027 1028 /* sparse solver using conjugate gradients for 1029 Poisson equation */ 1030 extern int 1031 re_number_and_solve_poisson_sparse(void); 1032 1033 /* 1034 * eig_lapack.c 1035 */ 1036 1037 /* solves the symmetrix, generalized eigenvalue problem 1038 * (P-lambda. T)x= 0 1039 * P, T = N by N matrices, only lower triangle 1040 * x = eigenvectors to be found size (N x M) 1041 * using LAPACK 1042 */ 1043 #if defined (USE_LAPACK) || defined (USE_MKL) 1044 /* symmetric definite eigensolver */ 1045 extern void 1046 solve_helmholtz_lapack(MY_DOUBLE **P,MY_DOUBLE **T, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT N, MY_INT G); 1047 extern void 1048 solve_helmholtz_lapack_sparse(SPMat *P,SPMat *T, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT N, MY_INT G); 1049 1050 /* symmetric non definite eigensolver */ 1051 extern void 1052 solve_generalized_eig_lapack(SPMat P,SPMat T, MY_DOUBLE **x, MY_DOUBLE *ev, MY_INT N, MY_INT G); 1053 #endif /* !USE_LAPACK || USE_MKL */ 1054 1055 1056 1057 /* 1058 * hash. c 1059 */ 1060 /* data structure for hash table */ 1061 typedef struct htbl_ { 1062 int positions; /* no. of positions in table */ 1063 int size; /* no of items inserted */ 1064 1065 /* hash functions used */ 1066 unsigned int (*h) (const void *key); 1067 /* function for comparison */ 1068 int (*match) (const void *key1, const void *key2); 1069 /* function for deletion */ 1070 void (*destroy) (void *data); 1071 1072 /* array of linked lists */ 1073 glist *table; 1074 } htbl; 1075 1076 1077 /* function interface */ 1078 /* intialize hash table */ 1079 extern void 1080 htbl_init(htbl *tbl, int hashtable_size, size_t datasize, unsigned int (*h)(const void *key), 1081 int (*match)(const void *key1, const void *key2), 1082 void (*destroy)(void *data)); 1083 /* destroy hash table */ 1084 extern void 1085 htbl_destroy(htbl* tbl); 1086 /* insert data */ 1087 /* if data inserted return 0, else return 1 (data already there) */ 1088 extern int 1089 htbl_insert(htbl *tbl, const void *data); 1090 /* remove data */ 1091 /* if data is removed return 0, else return 1 */ 1092 extern int 1093 htbl_remove(htbl *tbl,const void *data); 1094 /* lookup data and return address*/ 1095 /* if data not found, null is returned */ 1096 extern void * 1097 htbl_lookup(const htbl *tbl,const void *data); 1098 1099 /* 1100 * surfmesh.c 1101 */ 1102 1103 #ifndef WIN32 1104 /* display 3d view */ 1105 extern void 1106 display_mesh_grid(GtkWidget *widget, gpointer data); 1107 1108 /* to keep track of surf mesh window from main window */ 1109 extern GtkWidget *global_surf_mesh_window; 1110 #endif /* ! WIN32 */ 1111 1112 /* 1113 * eig_arpack.c 1114 */ 1115 extern int arpack_max_iterations; 1116 extern int arpack_subspace_dimension; 1117 1118 /* find generalized eigenvalues of the pencil (A,B) using ARPACK 1119 store eigenvectors in E - size N x G 1120 and eigenvalues in r - size G x 1 1121 G - the no. of eigenvalues to find 1122 lower_triangular - if 1, matrices A, B are lower triangular 1123 */ 1124 extern int 1125 arpack_find_eigs(SPMat *A, SPMat *B, MY_DOUBLE **v, MY_DOUBLE *r, MY_INT G, int lower_triangular); 1126 /* dummy drivers when not using ARPACK */ 1127 #ifndef USE_ARPACK 1128 extern void dsaupd_(int *ido, char *bmat, int *n, char *which, 1129 int *nev, double *tol, double *resid, int *ncv, 1130 double *v, int *ldv, int *param, int *ipntr, 1131 double *workd, double *workl, int *lworkl, 1132 int *info, unsigned long bmat_len, unsigned long which_len); 1133 extern void dseupd_(int *rvec, char *all, int *select_, double *d, 1134 double *v1, int *ldv1, double *sigma, 1135 char *bmat, int *n, char *which, int *nev, 1136 double *tol, double *resid, int *ncv, double *v, 1137 int *ldv, int *iparam, int *ipntr, double *workd, 1138 double *workl, int *lworkl, int *ierr); 1139 #endif /* USE_ARPACK */ 1140 1141 1142 #ifdef __cplusplus 1143 } /* closing brace for extern "C" */ 1144 #endif 1145 #endif /* ! TYPES_H*/ 1146