/* pdnmesh - a 2D finite element solver Copyright (C) 2001-2005 Sarod Yatawatta This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA $Id: integrate.c,v 1.7 2005/03/12 09:21:05 sarod Exp $ */ #ifdef HAVE_CONFIG_H #include #endif #include #include "types.h" static int on_line(poly_edge e, MY_INT p); /* the following function was sent by Werner Hoch */ /* calc the integral separate for each dirichlet edge */ /* input int con - eigenmode of potential */ #ifdef WIN32 void calc_integral(int con) { MY_INT n, al_p1, al_p2, al_p3; MY_DOUBLE yy23,yy31,yy12,xx23,xx31,xx12,del,a,b,ez; MY_DOUBLE dx,dy; MY_DOUBLE q=0.0; /* charge */ triangle *t; poly_edge e; /* examine each edge, take the integral in the distance of a half trianlge */ for ( n=0; n < nedges; n++) { #ifdef DEBUG printf("calc_integral: edge "MIF"\n",n); #endif /* only Dirichlet is for interest */ if (edge_array[n].type==DR) { e.p1 = edge_array[n].p1; e.p2 = edge_array[n].p2; #ifdef DEBUG printf("calc_integral: edge is Dirichlet\n"); #endif q = 0.0; /* examine each triangle */ DAG_traverse_list_reset(&dt); t=DAG_traverse_prune_list(&dt); while(t) { if ( t && t->status==LIVE ) { al_p1= on_line( e, t->p0); al_p2= on_line( e, t->p1); al_p3= on_line( e, t->p2); if (( al_p1 || al_p2 || al_p3)) { /* triangle has to touch the edge */ /* calculate gradient */ /* phi =a*x+b*y+c for each triangle */ /* need to find a and b only */ /* using cramer's rule for this */ yy23 = My(t->p1,M)-My(t->p2,M); yy31 = My(t->p2,M)-My(t->p0,M); yy12 = My(t->p0,M)-My(t->p1,M); xx23 = Mx(t->p1,M)-Mx(t->p2,M); xx31 = Mx(t->p2,M)-Mx(t->p0,M); xx12 = Mx(t->p0,M)-Mx(t->p1,M); del = Mx(t->p0,M)*yy23 + Mx(t->p1,M)*yy31 + Mx(t->p2,M)*yy12; if (!(IS_ZERO(del))) { del = 1.0/del; a = del*(Mzz(t->p0,con,M)*yy23 + Mzz(t->p1,con,M)*yy31 + Mzz(t->p2,con,M)*yy12); /* b = del*(Mx(t->p0)*(Mzz(t->p1,con)-Mzz(t->p2,con)) +Mx(t->p1)*(Mzz(t->p2,con)-Mzz(t->p0,con)) +Mx(t->p2)*(Mzz(t->p0,con)-Mzz(t->p1,con))); */ b = -del*(Mzz(t->p0,con,M)*yy23 + Mzz(t->p1,con,M)*yy31 + Mzz(t->p2,con,M)*yy12); if ( al_p1 && al_p2 && !al_p3 ) { dx = xx12*0.5; dy = yy12*0.5; } else if ( al_p1 && !al_p2 && al_p3 ) { dx = xx31*0.5; dy = yy31*0.5; } else if ( !al_p1 && al_p2 && al_p3 ) { dx = xx23*0.5; dy = yy23*0.5; } else if ( al_p1 && !al_p2 && !al_p3 ) { dx = - xx23*0.5; dy = - yy23*0.5; } else if ( !al_p1 && al_p2 && !al_p3 ) { dx = - xx31*0.5; dy = - yy31*0.5; } else if ( !al_p1 && !al_p2 && al_p3 ) { dx = - xx12*0.5; dy = - yy12*0.5; } else { /* this should never happen */ fprintf(stderr,"%s: %d: three colinear points?",__FILE__,__LINE__); dx=dy=0.0; } /* we calc the cross product (dx, dy, 0) X (a, b, 0) = (0, 0, ez) */ ez = (a*dy -b*dx)*bound.poly_array[t->boundary].mu; q += ez; #ifdef DEBUG printf("triangle "MIF", dx: "MDF", dy: "MDF", ez: "MDF"\n", t->n, dx, dy, ez); printf(" points on line: p1: %d, p2: %d, p3: %d\n",al_p1,al_p2,al_p3); printf(" p1.x: "MDF" p1.y: "MDF"\n",Mx(t->p0,M), My(t->p0,M)); printf(" p2.x: "MDF" p2.y: "MDF"\n",Mx(t->p1,M), My(t->p1,M)); printf(" p3.x: "MDF" p3.y: "MDF"\n",Mx(t->p2,M), My(t->p2,M)); #endif } } } t=DAG_traverse_prune_list(&dt); } printf("(edge "MIF" charge="MDF")\n",n,q); } } printf("\n"); /*return q; */ } #endif /* ! WIN32 */ #ifndef WIN32 /* calc the integral separate for each dirichlet edge */ /* input int con - eigenmode of potential */ /* display output on a gtk window */ void calc_integral( int con ) { GtkWidget *window, *box, *hbox, *llabel,*rlabel,*button; GtkWidget *scrolled_window,*outer_box; gchar buf[256]; /* variables for calculation of integral */ MY_INT n, al_p1, al_p2, al_p3; MY_DOUBLE yy23,yy31,yy12,xx23,xx31,xx12,del,a,b,ez; MY_DOUBLE dx,dy; MY_DOUBLE q=0.0; /* charge */ triangle *t; poly_edge e; #ifdef DEBUG printf ("calc_integral:\n"); #endif window=gtk_window_new(GTK_WINDOW_TOPLEVEL); g_signal_connect_swapped(GTK_OBJECT(window),"delete_event", G_CALLBACK(gtk_widget_destroy),(gpointer)window); gtk_window_set_title(GTK_WINDOW(window),"Integration"); gtk_container_set_border_width(GTK_CONTAINER(window),0); gtk_widget_set_size_request(window,300,100); outer_box=gtk_vbox_new(FALSE,0); gtk_container_add(GTK_CONTAINER(window),outer_box); gtk_container_set_border_width(GTK_CONTAINER(outer_box),10); gtk_widget_show(outer_box); /* scrolled window */ scrolled_window=gtk_scrolled_window_new(NULL,NULL); gtk_scrolled_window_set_policy(GTK_SCROLLED_WINDOW(scrolled_window), GTK_POLICY_AUTOMATIC,GTK_POLICY_ALWAYS); gtk_box_pack_start(GTK_BOX(outer_box),scrolled_window,TRUE,TRUE,0); gtk_widget_show(scrolled_window); /* inner box within the scrolled window */ box=gtk_vbox_new(FALSE,0); gtk_scrolled_window_add_with_viewport( GTK_SCROLLED_WINDOW(scrolled_window),box); gtk_container_set_border_width(GTK_CONTAINER(box),10); gtk_widget_show(box); /********** actual calculation of integral here */ /* examine each edge, take the integral in the distance of a half trianlge */ for ( n=0; n < nedges; n++) { #ifdef DEBUG printf("calc_integral: edge "MIF"\n",n); #endif /* only Dirichlet is for interest */ if (edge_array[n].type==DR) { e.p1 = edge_array[n].p1; e.p2 = edge_array[n].p2; #ifdef DEBUG printf("calc_integral: edge is Dirichlet\n"); #endif q = 0.0; /* examine each triangle */ DAG_traverse_list_reset(&dt); t=DAG_traverse_prune_list(&dt); while(t) { if ( t && t->status==LIVE ) { al_p1= on_line( e, t->p0); al_p2= on_line( e, t->p1); al_p3= on_line( e, t->p2); if (( al_p1 || al_p2 || al_p3)) { /* triangle has to touch the edge */ /* calculate gradient */ /* phi =a*x+b*y+c for each triangle */ /* need to find a and b only */ /* using cramer's rule for this */ yy23 = My(t->p1,M)-My(t->p2,M); yy31 = My(t->p2,M)-My(t->p0,M); yy12 = My(t->p0,M)-My(t->p1,M); xx23 = Mx(t->p1,M)-Mx(t->p2,M); xx31 = Mx(t->p2,M)-Mx(t->p0,M); xx12 = Mx(t->p0,M)-Mx(t->p1,M); del = Mx(t->p0,M)*yy23 + Mx(t->p1,M)*yy31 + Mx(t->p2,M)*yy12; if (!(IS_ZERO(del))) { del = 1.0/del; a = del*(Mzz(t->p0,con,M)*yy23 + Mzz(t->p1,con,M)*yy31 + Mzz(t->p2,con,M)*yy12); /* b = del*(Mx(t->p0)*(Mzz(t->p1,con)-Mzz(t->p2,con)) +Mx(t->p1)*(Mzz(t->p2,con)-Mzz(t->p0,con)) +Mx(t->p2)*(Mzz(t->p0,con)-Mzz(t->p1,con))); */ b = -del*(Mzz(t->p0,con,M)*yy23 + Mzz(t->p1,con,M)*yy31 + Mzz(t->p2,con,M)*yy12); if ( al_p1 && al_p2 && !al_p3 ) { dx = xx12*0.5; dy = yy12*0.5; } else if ( al_p1 && !al_p2 && al_p3 ) { dx = xx31*0.5; dy = yy31*0.5; } else if ( !al_p1 && al_p2 && al_p3 ) { dx = xx23*0.5; dy = yy23*0.5; } else if ( al_p1 && !al_p2 && !al_p3 ) { dx = - xx23*0.5; dy = - yy23*0.5; } else if ( !al_p1 && al_p2 && !al_p3 ) { dx = - xx31*0.5; dy = - yy31*0.5; } else if ( !al_p1 && !al_p2 && al_p3 ) { dx = - xx12*0.5; dy = - yy12*0.5; } else { /* this should never happen */ fprintf(stderr,"%s: %d: three colinear points?",__FILE__,__LINE__); dx=dy=0.0; } /* we calc the cross product (dx, dy, 0) X (a, b, 0) = (0, 0, ez) */ ez = (a*dy -b*dx)*bound.poly_array[t->boundary].mu; q += ez; #ifdef DEBUG printf("triangle "MIF", dx: "MDF", dy: "MDF", ez: "MDF"\n", t->n, dx, dy, ez); printf(" points on line: p1: %d, p2: %d, p3: %d\n",al_p1,al_p2,al_p3); printf(" p1.x: "MDF" p1.y: "MDF"\n",Mx(t->p0,M), My(t->p0,M)); printf(" p2.x: "MDF" p2.y: "MDF"\n",Mx(t->p1,M), My(t->p1,M)); printf(" p3.x: "MDF" p3.y: "MDF"\n",Mx(t->p2,M), My(t->p2,M)); #endif } } } t=DAG_traverse_prune_list(&dt); } hbox=gtk_hbox_new(FALSE,0); gtk_container_set_border_width(GTK_CONTAINER(hbox),0); gtk_box_pack_start(GTK_BOX(box),hbox,TRUE,FALSE,0); gtk_widget_show(hbox); /* printf("(edge "MIF" charge="MDF")\n",n,q); */ llabel=gtk_label_new(""); gtk_box_pack_start(GTK_BOX(hbox),llabel,FALSE,FALSE,0); gtk_widget_show(llabel); gtk_label_set_justify(GTK_LABEL(llabel),GTK_JUSTIFY_LEFT); rlabel=gtk_label_new(""); gtk_box_pack_start(GTK_BOX(hbox),rlabel,FALSE,FALSE,0); gtk_widget_show(rlabel); gtk_label_set_justify(GTK_LABEL(rlabel),GTK_JUSTIFY_LEFT); g_snprintf(buf, 256, "edge\t "MIF"\t ",n); gtk_label_set_markup(GTK_LABEL(llabel),buf); g_snprintf(buf, 256, "charge= "MDF"",q); gtk_label_set_markup(GTK_LABEL(rlabel),buf); } } /********** end actual calculation of integral */ button=gtk_button_new_with_label("OK"); g_signal_connect_swapped(G_OBJECT(button),"clicked",G_CALLBACK(gtk_widget_destroy), G_OBJECT(window)); gtk_box_pack_start(GTK_BOX(outer_box),button,FALSE,FALSE,0); gtk_widget_show(button); gtk_widget_show(window); } #endif /* WIN32 */ /* a function to check whether a point p is on-line an edge e */ /* source taken from subdivide.c is_on_edge */ static int on_line(poly_edge e, MY_INT p) { MY_DOUBLE dx1,dx2,dy1,dy2; dx1 = Mx(e.p1,M)-Mx(p,M); #ifdef DEBUG printf("dx1 "MDF" ",dx1); #endif dy1 = My(e.p1,M)-My(p,M); #ifdef DEBUG printf("dy1 "MDF" ",dy1); #endif dx2 = Mx(p,M)-Mx(e.p2,M); #ifdef DEBUG printf("dx2 "MDF" ",dx2); #endif dy2 = My(p,M)-My(e.p2,M); #ifdef DEBUG printf("dy2 "MDF"\n",dy2); #endif /* first the coincident case */ if ( ( IS_ZERO(dx1) && IS_ZERO(dy1) ) || ( IS_ZERO(dx2)&& IS_ZERO(dy2) )) { #ifdef DEBUG printf("******************case1\n"); #endif return(1); } /* now to check parrallbility */ /* we allow a toelarance, otherwise unpredictable */ if ( IS_ZERO(dx1*dy2 - dx2*dy1)) { #ifdef DEBUG printf("__________________case2\n"); #endif /* now check whether point is inbetween */ /* this again gave problems, but when I changed */ /* tha and to an or, and > to a < it worked */ if ( (dx1*dx2 >= TOL && dy1*dy2 >= TOL ) || ( (ABS(dx1) < TOL) && (ABS(dx2) < TOL ) && dy1*dy2 >= TOL ) || ( (ABS(dy1) < TOL ) && (ABS(dy2) < TOL ) && dx1*dx2 >= TOL ) ) return(1); return(0); } return(0); }