Trisurf Monte Carlo simulator
Samo Penic
2014-06-15 7d62fc765adc11bbcb7db0716ac5dc19145dca80
commit | author | age
0af0cf 1 #include<general.h>
SP 2 #include<cross-section.h>
3 #include<coord.h>
a61c00 4 #include<cairo/cairo.h>
0af0cf 5 /** @brief Calculates cross-section of vesicle with plane.
SP 6  *
7  * Function returns points of cross-section of vesicle with plane. Plane is described with equation $ax+by+cz+d=0$. Algorithm extracts coordinates of each vertex of a vesicle and then:
8  *
9  * if a distance of point to plane (given by equation $D=\frac{ax_0+by_0+cz_0+d}{\sqrt{a^2+b^2+c^2}}$, where $x_0$, $y_0$ and $z_0$ are coordinates of a given vertex) is less than maximal allowed distance between vertices {\tt sqrt(vesicle->dmax)} than vertex is a candidate for crossection calculation.
10  *
11  */
d27c07 12 ts_coord_list *get_crossection_with_plane(ts_vesicle *vesicle,ts_double a,ts_double b,ts_double c, ts_double d){
0af0cf 13
SP 14
ad4e70 15     ts_uint i, j,k,l;
0af0cf 16     ts_double pp,Dsq; // distance from the plane squared
ad4e70 17     ts_double ppn1; // distance from the plane squared of a neighbor
0af0cf 18     ts_vertex *vtx;
a61c00 19     ts_uint ntria=0; // number triangles
SP 20     ts_triangle *tria[2]; // list of triangles
0af0cf 21     ts_coord_list *pts=init_coord_list();
d27c07 22     for(i=0;i<vesicle->vlist->n;i++){
0af0cf 23         vtx=vesicle->vlist->vtx[i];
SP 24
25         pp=vtx->x*a+vtx->y*b+vtx->z*c+d;
26         Dsq=pp*pp/(a*a+b*b+c*c);
27         if(Dsq<vesicle->dmax){
28             for(j=0;j<vtx->neigh_no;j++){
a61c00 29                 ppn1=vtx->neigh[j]->x*a+vtx->neigh[j]->y*b+vtx->neigh[j]->z*c+d;
7d62fc 30                 if(pp*ppn1<=0.0){ //the combination of vertices are good candidates for a crossection
a61c00 31                     //find triangle that belongs to the two vertices
SP 32                     ntria=0;
33                     for(k=0;k<vtx->tristar_no;k++){
ad4e70 34                         if(vtx->tristar[k]->vertex[0]==vtx && ( vtx->tristar[k]->vertex[1]==vtx->neigh[j] || vtx->tristar[k]->vertex[2]==vtx->neigh[j]) ){ 
a61c00 35                             //triangle found.
SP 36                             tria[ntria]=vtx->tristar[k];
37                             ntria++;
38                         }
39                     }
7d62fc 40                     if(ntria==0) continue; // no need to continue
a61c00 41                     //find the two intersections (in general) to form a intersection line
ad4e70 42                     for(l=0;l<ntria;l++){
SP 43                     //we add intersection line between two points for each of the triangles found above.
44                         add_crosssection_point(pts,a,b,c,d,tria[l]->vertex[0], tria[l]->vertex[1]);
45                         add_crosssection_point(pts,a,b,c,d,tria[l]->vertex[0], tria[l]->vertex[2]);
46                         add_crosssection_point(pts,a,b,c,d,tria[l]->vertex[1], tria[l]->vertex[2]);
a61c00 47
SP 48                     }
0af0cf 49                 }
SP 50             }
51         }
52     }
53     return pts;
54 }
55
56
ad4e70 57
SP 58 ts_bool add_crosssection_point(ts_coord_list *pts, ts_double a, ts_double b, ts_double c, ts_double d, ts_vertex *vtx1, ts_vertex *vtx2){
59     ts_double pp=vtx1->x*a+vtx1->y*b+vtx1->z*c+d;
60     ts_double pp2=vtx2->x*a+vtx1->y*b+vtx2->z*c+d;
7d62fc 61     if(pp*pp2<=0.0){
ad4e70 62     ts_double u=pp/(a*(vtx1->x-vtx2->x)+b*(vtx1->y-vtx2->y)+c*(vtx1->z-vtx2->z));
SP 63     add_coord(pts, vtx1->x+u*(vtx2->x - vtx1->x),
64         vtx1->y+u*(vtx2->y - vtx1->y),        
65         vtx1->z+u*(vtx2->z - vtx1->z),
66         TS_COORD_CARTESIAN);
67     
68     return TS_SUCCESS;
69     } else {
70     return TS_FAIL;
71     }
72 }
73
a61c00 74 /** Saves calculated crossection as a png image */
SP 75 ts_bool crossection_to_png(ts_coord_list *pts, char *filename){
76
77     cairo_surface_t *surface;
78     cairo_t *cr;
79     ts_uint i;
ad4e70 80     surface = cairo_image_surface_create (CAIRO_FORMAT_RGB24, 1800, 1800);
a61c00 81     cr = cairo_create (surface);
SP 82     cairo_rectangle(cr, 0.0, 0.0, 1800,1800);
83     cairo_set_source_rgb(cr, 0.3, 0.3, 0.3);
84     cairo_fill(cr);
85     cairo_set_line_width (cr, 5.0/30.0);
324ad7 86     cairo_set_line_cap(cr, CAIRO_LINE_CAP_ROUND);
SP 87     cairo_set_line_join(cr, CAIRO_LINE_JOIN_ROUND); 
a61c00 88     cairo_translate(cr, 900,900);
SP 89     cairo_scale (cr, 30, 30);
90     cairo_set_source_rgb (cr, 1.0, 1.0, 1.0);
91
92     for(i=0;i<pts->n;i+=2){
93         cairo_move_to(cr, pts->coord[i]->e1, pts->coord[i]->e2);
94         cairo_line_to(cr, pts->coord[i+1]->e1, pts->coord[i+1]->e2);
95     }
96     cairo_stroke(cr);
97     cairo_surface_write_to_png (surface,filename);
98     cairo_surface_finish (surface);    
99     return TS_SUCCESS;
100 }
ad4e70 101
SP 102
103 ts_bool save_crossection_snapshot(ts_coord_list *pts, ts_uint timestepno){
104     char filename[255];
105         sprintf(filename,"timestep_%.6u.png",timestepno);
106     crossection_to_png(pts,filename);
107     return TS_SUCCESS;
108 }