commit | author | age
|
7958e9
|
1 |
#include<stdlib.h> |
SP |
2 |
#include<math.h> |
|
3 |
#include<stdio.h> |
|
4 |
#include "general.h" |
|
5 |
#include "vertex.h" |
|
6 |
#include "bond.h" |
|
7 |
#include "vesicle.h" |
|
8 |
#include "vertex.h" |
|
9 |
#include "triangle.h" |
|
10 |
#include "initial_distribution.h" |
f74313
|
11 |
#include "energy.h" |
1ab449
|
12 |
#include "poly.h" |
8a6614
|
13 |
#include "io.h" |
dc77e8
|
14 |
#include "sh.h" |
459ff9
|
15 |
#include "shcomplex.h" |
7958e9
|
16 |
|
SP |
17 |
ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){ |
1ab449
|
18 |
ts_fprintf(stdout,"Starting initial_distribution on vesicle with %u shells!...\n",nshell); |
7958e9
|
19 |
ts_bool retval; |
1ab449
|
20 |
ts_uint no_vertices=5*nshell*nshell+2; |
SP |
21 |
ts_vesicle *vesicle=init_vesicle(no_vertices,ncmax1,ncmax2,ncmax3,stepsize); |
|
22 |
vesicle->nshell=nshell; |
|
23 |
//retval = vtx_set_global_values(vesicle); |
|
24 |
retval = pentagonal_dipyramid_vertex_distribution(vesicle->vlist); |
|
25 |
retval = init_vertex_neighbours(vesicle->vlist); |
|
26 |
vesicle->vlist = init_sort_neighbours(vesicle->blist,vesicle->vlist); |
b01cc1
|
27 |
// retval = init_vesicle_bonds(vesicle); // bonds are created in sort_neigh |
1ab449
|
28 |
retval = init_triangles(vesicle); |
SP |
29 |
retval = init_triangle_neighbours(vesicle); |
|
30 |
retval = init_common_vertex_triangle_neighbours(vesicle); |
|
31 |
retval = init_normal_vectors(vesicle->tlist); |
|
32 |
retval = mean_curvature_and_energy(vesicle); |
|
33 |
ts_fprintf(stdout,"initial_distribution finished!\n"); |
41a035
|
34 |
if(retval); |
7958e9
|
35 |
return vesicle; |
SP |
36 |
} |
|
37 |
|
|
38 |
|
1ab449
|
39 |
|
SP |
40 |
ts_vesicle *create_vesicle_from_tape(ts_tape *tape){ |
|
41 |
ts_vesicle *vesicle; |
bcf455
|
42 |
ts_vertex *vtx; |
M |
43 |
|
1ab449
|
44 |
vesicle=initial_distribution_dipyramid(tape->nshell,tape->ncxmax,tape->ncymax,tape->nczmax,tape->stepsize); |
9166cb
|
45 |
vesicle->tape=tape; |
58230a
|
46 |
// Nucleus: |
fe24d2
|
47 |
vesicle->R_nucleus=tape->R_nucleus*tape->R_nucleus; |
M |
48 |
|
|
49 |
vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies; |
bcf455
|
50 |
|
58230a
|
51 |
//Initialize grafted polymers (brush): |
624f81
|
52 |
vesicle->poly_list=init_poly_list(tape->npoly,tape->nmono, vesicle->vlist, vesicle); |
1ab449
|
53 |
vesicle->spring_constant=tape->kspring; |
SP |
54 |
poly_assign_spring_const(vesicle); |
bcf455
|
55 |
|
58230a
|
56 |
//Initialize filaments (polymers inside the vesicle): |
M |
57 |
vesicle->filament_list=init_poly_list(tape->nfil,tape->nfono, NULL, vesicle); |
bcf455
|
58 |
poly_assign_filament_xi(vesicle,tape); |
58230a
|
59 |
|
bcf455
|
60 |
ts_uint i,j; |
M |
61 |
for(i=0;i<vesicle->filament_list->n;i++){ |
|
62 |
for(j=0;j<vesicle->filament_list->poly[i]->blist->n;j++){ |
|
63 |
bond_vector(vesicle->filament_list->poly[i]->blist->bond[j]); |
|
64 |
vesicle->filament_list->poly[i]->blist->bond[j]->bond_length = sqrt(vtx_distance_sq(vesicle->filament_list->poly[i]->blist->bond[j]->vtx1,vesicle->filament_list->poly[i]->blist->bond[j]->vtx2)); |
|
65 |
} |
58230a
|
66 |
} |
bcf455
|
67 |
|
M |
68 |
for(i=0;i<vesicle->filament_list->n;i++){ |
|
69 |
for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){ |
|
70 |
vtx = vesicle->filament_list->poly[i]->vlist->vtx[j]; |
|
71 |
if(vtx->bond_no == 2){ |
|
72 |
vtx->energy = -(vtx->bond[0]->x*vtx->bond[1]->x + vtx->bond[0]->y*vtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length; |
|
73 |
} |
|
74 |
} |
58230a
|
75 |
} |
bcf455
|
76 |
|
ea1cce
|
77 |
for(i=0;i<vesicle->filament_list->n;i++){ |
M |
78 |
vertex_list_assign_id(vesicle->filament_list->poly[i]->vlist,TS_ID_FILAMENT); |
|
79 |
} |
bcf455
|
80 |
|
58230a
|
81 |
// vesicle->spring_constant=tape->kspring; |
M |
82 |
// poly_assign_spring_const(vesicle); |
|
83 |
|
1ab449
|
84 |
|
SP |
85 |
vesicle->nshell=tape->nshell; |
|
86 |
vesicle->dmax=tape->dmax*tape->dmax; /* dmax^2 in the vesicle dmax variable */ |
|
87 |
vesicle->bending_rigidity=tape->xk0; |
|
88 |
vtx_set_global_values(vesicle); /* make xk0 default value for every vertex */ |
|
89 |
ts_fprintf(stdout, "Tape setting: xk0=%e\n",tape->xk0); |
|
90 |
vesicle->stepsize=tape->stepsize; |
|
91 |
vesicle->clist->ncmax[0]=tape->ncxmax; |
|
92 |
vesicle->clist->ncmax[1]=tape->ncymax; |
|
93 |
vesicle->clist->ncmax[2]=tape->nczmax; |
|
94 |
vesicle->clist->max_occupancy=8; /* hard coded max occupancy? */ |
|
95 |
|
|
96 |
vesicle->pressure= tape->pressure; |
|
97 |
vesicle->pswitch=tape->pswitch; |
632960
|
98 |
if(tape->shc>0){ |
459ff9
|
99 |
vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc); |
632960
|
100 |
} |
SP |
101 |
else { |
|
102 |
vesicle->sphHarmonics=NULL; |
|
103 |
} |
1ab449
|
104 |
return vesicle; |
SP |
105 |
|
|
106 |
} |
|
107 |
|
|
108 |
|
|
109 |
|
|
110 |
|
|
111 |
|
7958e9
|
112 |
ts_bool pentagonal_dipyramid_vertex_distribution(ts_vertex_list *vlist){ |
SP |
113 |
/* Some often used relations */ |
|
114 |
const ts_double s1= sin(2.0*M_PI/5.0); |
|
115 |
const ts_double s2= sin(4.0*M_PI/5.0); |
|
116 |
const ts_double c1= cos(2.0*M_PI/5.0); |
|
117 |
const ts_double c2= cos(4.0*M_PI/5.0); |
|
118 |
|
|
119 |
/* Calculates projection lenght of an edge bond to pentagram plane */ |
|
120 |
const ts_double xl0=A0/(2.0*sin(M_PI/5.0)); |
|
121 |
#ifdef TS_DOUBLE_DOUBLE |
|
122 |
const ts_double z0=sqrt(pow(A0,2)-pow(xl0,2)); |
|
123 |
#endif |
|
124 |
#ifdef TS_DOUBLE_FLOAT |
|
125 |
const ts_double z0=sqrtf(powf(A0,2)-powf(xl0,2)); |
|
126 |
#endif |
|
127 |
#ifdef TS_DOUBLE_LONGDOUBLE |
|
128 |
const ts_double z0=sqrtl(powl(A0,2)-powl(xl0,2)); |
|
129 |
#endif |
|
130 |
// const z0=sqrt(A0*A0 -xl0*xl0); /* I could use pow function but if pow is used make a check on the float type. If float then powf, if long double use powl */ |
|
131 |
|
|
132 |
/*placeholder for the pointer to vertex datastructure list... DIRTY: actual pointer points towards invalid address, one position before actual beginning of the list... This is to solve the difference between 1 based indexing in original program in fortran and 0 based indexing in C. All algorithms remain unchanged because of this!*/ |
|
133 |
ts_vertex **vtx=vlist->vtx -1 ; |
|
134 |
|
|
135 |
|
|
136 |
ts_uint nshell=(ts_uint)( sqrt((ts_double)(vlist->n-2)/5)); |
|
137 |
// printf("nshell=%u\n",nshell); |
|
138 |
ts_uint i,n0; // some for loop prereq |
|
139 |
ts_int j,k; |
|
140 |
ts_double dx,dy; // end loop prereq |
|
141 |
|
|
142 |
/* topmost vertex */ |
8f6a69
|
143 |
vtx[1]->x=0.0; |
SP |
144 |
vtx[1]->y=0.0; |
|
145 |
vtx[1]->z=z0*(ts_double)nshell; |
7958e9
|
146 |
|
SP |
147 |
/* starting from to in circular order on pentagrams */ |
|
148 |
for(i=1;i<=nshell;i++){ |
|
149 |
n0=2+5*i*(i-1)/2; //-1 would be for the reason that C index starts from 0 |
8f6a69
|
150 |
vtx[n0]->x=0.0; |
SP |
151 |
vtx[n0]->y=(ts_double)i*xl0; |
|
152 |
vtx[n0+i]->x=vtx[n0]->y*s1; |
|
153 |
vtx[n0+i]->y=vtx[n0]->y*c1; |
|
154 |
vtx[n0+2*i]->x=vtx[n0]->y*s2; |
|
155 |
vtx[n0+2*i]->y=vtx[n0]->y*c2; |
|
156 |
vtx[n0+3*i]->x=-vtx[n0+2*i]->x; |
|
157 |
vtx[n0+3*i]->y=vtx[n0+2*i]->y; |
|
158 |
vtx[n0+4*i]->x=-vtx[n0+i]->x; |
|
159 |
vtx[n0+4*i]->y=vtx[n0+i]->y; |
7958e9
|
160 |
} |
SP |
161 |
|
|
162 |
/* vertexes on the faces of the dipyramid */ |
|
163 |
for(i=1;i<=nshell;i++){ |
|
164 |
n0=2+5*i*(i-1)/2; // -1 would be because of C! |
|
165 |
for(j=1;j<=i-1;j++){ |
8f6a69
|
166 |
dx=(vtx[n0]->x-vtx[n0+4*i]->x)/(ts_double)i; |
SP |
167 |
dy=(vtx[n0]->y-vtx[n0+4*i]->y)/(ts_double)i; |
|
168 |
vtx[n0+4*i+j]->x=(ts_double)j*dx+vtx[n0+4*i]->x; |
|
169 |
vtx[n0+4*i+j]->y=(ts_double)j*dy+vtx[n0+4*i]->y; |
7958e9
|
170 |
} |
SP |
171 |
for(k=0;k<=3;k++){ // I would be worried about zero starting of for |
8f6a69
|
172 |
dx=(vtx[n0+(k+1)*i]->x - vtx[n0+k*i]->x)/(ts_double) i; |
SP |
173 |
dy=(vtx[n0+(k+1)*i]->y - vtx[n0+k*i]->y)/(ts_double) i; |
7958e9
|
174 |
for(j=1; j<=i-1;j++){ |
8f6a69
|
175 |
vtx[n0+k*i+j]->x= (ts_double)j*dx+vtx[n0+k*i]->x; |
SP |
176 |
vtx[n0+k*i+j]->y= (ts_double)j*dy+vtx[n0+k*i]->y; |
7958e9
|
177 |
} |
SP |
178 |
} |
|
179 |
} |
|
180 |
|
|
181 |
for(i=1;i<=nshell;i++){ |
|
182 |
n0= 2+ 5*i*(i-1)/2; |
|
183 |
for(j=0;j<=5*i-1;j++){ |
8f6a69
|
184 |
vtx[n0+j]->z= z0*(ts_double)(nshell-i); // I would be worried about zero starting of for |
7958e9
|
185 |
} |
SP |
186 |
} |
|
187 |
|
|
188 |
/* for botom part of dipyramide we calculate the positions of vertices */ |
|
189 |
for(i=2+5*nshell*(nshell+1)/2;i<=vlist->n;i++){ |
8f6a69
|
190 |
vtx[i]->x=vtx[vlist->n - i +1]->x; |
SP |
191 |
vtx[i]->y=vtx[vlist->n - i +1]->y; |
|
192 |
vtx[i]->z=-vtx[vlist->n - i +1]->z; |
7958e9
|
193 |
} |
SP |
194 |
|
|
195 |
for(i=1;i<=vlist->n;i++){ |
|
196 |
for(j=1;j<=vlist->n;j++){ |
|
197 |
if(i!=j && vtx_distance_sq(vtx[i],vtx[j])<0.001){ |
|
198 |
printf("Vertices %u and %u are the same!\n",i,j); |
|
199 |
} |
|
200 |
} |
|
201 |
} |
|
202 |
return TS_SUCCESS; |
|
203 |
} |
|
204 |
|
|
205 |
|
|
206 |
|
|
207 |
ts_bool init_vertex_neighbours(ts_vertex_list *vlist){ |
|
208 |
ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment. |
|
209 |
const ts_double eps=0.001; //TODO: find out if you can use EPS from math.h |
|
210 |
ts_uint i,j; |
|
211 |
ts_double dist2; // Square of distance of neighbours |
|
212 |
/*this is not required if we zero all data in vertex structure at initialization */ |
|
213 |
/*if we force zeroing at initialization this for loop can safely be deleted */ |
|
214 |
//for(i=1;i<=vlist->n;i++){ |
|
215 |
// vtx[i].neigh_no=0; |
|
216 |
//} |
|
217 |
for(i=1;i<=vlist->n;i++){ |
|
218 |
for(j=1;j<=vlist->n;j++){ |
|
219 |
dist2=vtx_distance_sq(vtx[i],vtx[j]); |
|
220 |
if( (dist2>eps) && (dist2<(A0*A0+eps))){ |
|
221 |
//if it is close enough, but not too much close (solves problem of comparing when i==j) |
|
222 |
vtx_add_neighbour(vtx[i],vtx[j]); |
|
223 |
} |
|
224 |
} |
|
225 |
// printf ("vertex %u ima %u sosedov!\n",i,vtx[i]->data->neigh_no); |
|
226 |
} |
|
227 |
|
|
228 |
return TS_SUCCESS; |
|
229 |
} |
|
230 |
|
b01cc1
|
231 |
// TODO: with new datastructure can be rewritten. Partially it is done, but it is complicated. |
SP |
232 |
ts_vertex_list *init_sort_neighbours(ts_bond_list *blist,ts_vertex_list *vlist){ |
7958e9
|
233 |
ts_vertex **vtx=vlist->vtx -1; // take a look at dipyramid function for comment. |
SP |
234 |
ts_uint i,l,j,jj,jjj,k=0; |
|
235 |
ts_double eps=0.001; // Take a look if EPS from math.h can be used |
|
236 |
|
|
237 |
/*lets initialize memory for temporary vertex_list. Should we write a function instead */ |
b01cc1
|
238 |
ts_vertex_list *tvlist=vertex_list_copy(vlist); |
7958e9
|
239 |
ts_vertex **tvtx=tvlist->vtx -1; /* again to compensate for 0-indexing */ |
SP |
240 |
|
|
241 |
ts_double dist2; // Square of distance of neighbours |
|
242 |
ts_double direct; // Something, dont know what, but could be normal of some kind |
|
243 |
for(i=1;i<=vlist->n;i++){ |
|
244 |
k++; // WHY i IS NOT GOOD?? |
8f6a69
|
245 |
vtx_add_cneighbour(blist,tvtx[k], tvtx[vtx[i]->neigh[0]->idx+1]); //always add 1st |
7958e9
|
246 |
jjj=1; |
SP |
247 |
jj=1; |
8f6a69
|
248 |
for(l=2;l<=vtx[i]->neigh_no;l++){ |
SP |
249 |
for(j=2;j<=vtx[i]->neigh_no;j++){ |
|
250 |
dist2=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]); |
|
251 |
direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]); |
|
252 |
// TODO: check if fabs can be used with all floating point types!! |
7958e9
|
253 |
if( (fabs(dist2-A0*A0)<=eps) && (direct>0.0) && (j!=jjj) ){ |
8f6a69
|
254 |
vtx_add_cneighbour(blist,tvtx[k],tvtx[vtx[i]->neigh[j-1]->idx+1]); |
7958e9
|
255 |
jjj=jj; |
SP |
256 |
jj=j; |
|
257 |
break; |
|
258 |
} |
|
259 |
} |
|
260 |
} |
|
261 |
} |
b01cc1
|
262 |
/* We use the temporary vertex for our main vertices and we abandon main |
SP |
263 |
* vertices, because their neighbours are not correctly ordered */ |
|
264 |
// tvtx=vlist->vtx; |
|
265 |
// vlist->vtx=tvtx; |
|
266 |
// tvlist->vtx=vtx; |
|
267 |
vtx_list_free(vlist); |
|
268 |
/* Let's make a check if the number of bonds is correct */ |
|
269 |
if((blist->n)!=3*(tvlist->n-2)){ |
|
270 |
ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(tvlist->n-2)); |
|
271 |
fatal("Number of bonds is not 3*(no_vertex-2).",4); |
7958e9
|
272 |
} |
SP |
273 |
|
b01cc1
|
274 |
return tvlist; |
7958e9
|
275 |
} |
SP |
276 |
|
|
277 |
|
|
278 |
ts_bool init_vesicle_bonds(ts_vesicle *vesicle){ |
|
279 |
ts_vertex_list *vlist=vesicle->vlist; |
|
280 |
ts_bond_list *blist=vesicle->blist; |
|
281 |
ts_vertex **vtx=vesicle->vlist->vtx - 1; // Because of 0 indexing |
|
282 |
/* lets make correct clockwise ordering of in nearest neighbour list */ |
|
283 |
ts_uint i,j,k; |
|
284 |
for(i=1;i<=vlist->n;i++){ |
|
285 |
for(j=i+1;j<=vlist->n;j++){ |
8f6a69
|
286 |
for(k=0;k<vtx[i]->neigh_no;k++){ // has changed 0 to < instead of 1 and <= |
SP |
287 |
if(vtx[i]->neigh[k]==vtx[j]){ //if addresses matches it is the same |
7958e9
|
288 |
bond_add(blist,vtx[i],vtx[j]); |
SP |
289 |
break; |
|
290 |
} |
|
291 |
} |
|
292 |
} |
|
293 |
} |
|
294 |
/* Let's make a check if the number of bonds is correct */ |
|
295 |
if((blist->n)!=3*(vlist->n-2)){ |
|
296 |
ts_fprintf(stderr,"Number of bonds is %u should be %u!\n", blist->n, 3*(vlist->n-2)); |
|
297 |
fatal("Number of bonds is not 3*(no_vertex-2).",4); |
|
298 |
} |
|
299 |
return TS_SUCCESS; |
|
300 |
} |
|
301 |
|
|
302 |
|
|
303 |
|
|
304 |
ts_bool init_triangles(ts_vesicle *vesicle){ |
|
305 |
ts_uint i,j,jj,k; |
|
306 |
ts_vertex **vtx=vesicle->vlist->vtx -1; // difference between 0 indexing and 1 indexing |
|
307 |
ts_triangle_list *tlist=vesicle->tlist; |
|
308 |
ts_double dist, direct; |
|
309 |
ts_double eps=0.001; // can we use EPS from math.h? |
|
310 |
k=0; |
|
311 |
for(i=1;i<=vesicle->vlist->n;i++){ |
8f6a69
|
312 |
for(j=1;j<=vtx[i]->neigh_no;j++){ |
SP |
313 |
for(jj=1;jj<=vtx[i]->neigh_no;jj++){ |
7958e9
|
314 |
// ts_fprintf(stderr,"%u: (%u,%u) neigh_no=%u ",i,j,jj,vtx[i].neigh_no); |
SP |
315 |
// ts_fprintf(stderr,"%e, %e",vtx[i].neigh[j-1]->x,vtx[i].neigh[jj-1]->x); |
8f6a69
|
316 |
dist=vtx_distance_sq(vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]); |
SP |
317 |
direct=vtx_direct(vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]); |
|
318 |
// TODO: same as above |
|
319 |
if(fabs(dist-A0*A0)<=eps && direct < 0.0 && vtx[i]->neigh[j-1]->idx+1 > i && vtx[i]->neigh[jj-1]->idx+1 >i){ |
|
320 |
triangle_add(tlist,vtx[i],vtx[i]->neigh[j-1],vtx[i]->neigh[jj-1]); |
7958e9
|
321 |
} |
SP |
322 |
} |
|
323 |
} |
|
324 |
} |
|
325 |
/* We check if all triangles have 3 vertices and if the number of triangles |
|
326 |
* matches the theoretical value. |
|
327 |
*/ |
|
328 |
for(i=0;i<tlist->n;i++){ |
|
329 |
k=0; |
|
330 |
for(j=0;j<3;j++){ |
41a035
|
331 |
if(tlist->tria[i]->vertex[j]!=NULL) |
7958e9
|
332 |
k++; |
SP |
333 |
} |
|
334 |
if(k!=3){ |
8f6a69
|
335 |
fatal("Some triangles have less than 3 vertices..",4); |
7958e9
|
336 |
} |
SP |
337 |
} |
|
338 |
if(tlist->n!=2*(vesicle->vlist->n -2)){ |
|
339 |
ts_fprintf(stderr,"The number of triangles is %u but should be %u!\n",tlist->n,2*(vesicle->vlist->n -2)); |
|
340 |
fatal("The number of triangles doesn't match 2*(no_vertex -2).",4); |
|
341 |
} |
|
342 |
return TS_SUCCESS; |
|
343 |
} |
|
344 |
|
|
345 |
|
|
346 |
|
|
347 |
ts_bool init_triangle_neighbours(ts_vesicle *vesicle){ |
|
348 |
ts_uint i,j,nobo; |
|
349 |
ts_vertex *i1,*i2,*i3,*j1,*j2,*j3; |
|
350 |
// ts_vertex **vtx=vesicle->vlist->vtx -1; // difference between 0 indexing and 1 indexing |
|
351 |
ts_triangle_list *tlist=vesicle->tlist; |
|
352 |
ts_triangle **tria=tlist->tria -1; |
|
353 |
nobo=0; |
|
354 |
for(i=1;i<=tlist->n;i++){ |
41a035
|
355 |
i1=tria[i]->vertex[0]; |
SP |
356 |
i2=tria[i]->vertex[1]; |
|
357 |
i3=tria[i]->vertex[2]; |
7958e9
|
358 |
for(j=1;j<=tlist->n;j++){ |
SP |
359 |
if(j==i) continue; |
41a035
|
360 |
j1=tria[j]->vertex[0]; |
SP |
361 |
j2=tria[j]->vertex[1]; |
|
362 |
j3=tria[j]->vertex[2]; |
7958e9
|
363 |
if((i1==j1 && i3==j2) || (i1==j2 && i3==j3) || (i1==j3 && i3==j1)){ |
SP |
364 |
triangle_add_neighbour(tria[i],tria[j]); |
|
365 |
nobo++; |
|
366 |
} |
|
367 |
} |
|
368 |
} |
|
369 |
for(i=1;i<=tlist->n;i++){ |
41a035
|
370 |
i1=tria[i]->vertex[0]; |
SP |
371 |
i2=tria[i]->vertex[1]; |
|
372 |
i3=tria[i]->vertex[2]; |
7958e9
|
373 |
for(j=1;j<=tlist->n;j++){ |
SP |
374 |
if(j==i) continue; |
41a035
|
375 |
j1=tria[j]->vertex[0]; |
SP |
376 |
j2=tria[j]->vertex[1]; |
|
377 |
j3=tria[j]->vertex[2]; |
7958e9
|
378 |
if((i1==j1 && i2==j3) || (i1==j3 && i2==j2) || (i1==j2 && i2==j1)){ |
SP |
379 |
triangle_add_neighbour(tria[i],tria[j]); |
|
380 |
nobo++; |
|
381 |
} |
|
382 |
} |
|
383 |
} |
|
384 |
for(i=1;i<=tlist->n;i++){ |
41a035
|
385 |
i1=tria[i]->vertex[0]; |
SP |
386 |
i2=tria[i]->vertex[1]; |
|
387 |
i3=tria[i]->vertex[2]; |
7958e9
|
388 |
for(j=1;j<=tlist->n;j++){ |
SP |
389 |
if(j==i) continue; |
41a035
|
390 |
j1=tria[j]->vertex[0]; |
SP |
391 |
j2=tria[j]->vertex[1]; |
|
392 |
j3=tria[j]->vertex[2]; |
7958e9
|
393 |
if((i2==j1 && i3==j3) || (i2==j3 && i3==j2) || (i2==j2 && i3==j1)){ |
SP |
394 |
triangle_add_neighbour(tria[i],tria[j]); |
|
395 |
nobo++; |
|
396 |
} |
|
397 |
} |
|
398 |
} |
|
399 |
if(nobo != vesicle->blist->n*2) { |
|
400 |
ts_fprintf(stderr,"Number of triangles= %u, number of bonds= %u\n",nobo/2, vesicle->blist->n); |
|
401 |
fatal("Number of triangle neighbour pairs differs from double the number of bonds!",4); |
|
402 |
} |
|
403 |
return TS_SUCCESS; |
|
404 |
} |
|
405 |
|
|
406 |
|
|
407 |
ts_bool init_common_vertex_triangle_neighbours(ts_vesicle *vesicle){ |
|
408 |
ts_uint i,j,jp,k; |
|
409 |
ts_vertex *k1,*k2,*k3,*k4,*k5; |
|
410 |
ts_vertex **vtx=vesicle->vlist->vtx -1; // difference between 0 indexing and 1 indexing |
|
411 |
ts_triangle_list *tlist=vesicle->tlist; |
|
412 |
ts_triangle **tria=tlist->tria -1; |
|
413 |
|
|
414 |
for(i=1;i<=vesicle->vlist->n;i++){ |
8f6a69
|
415 |
for(j=1;j<=vtx[i]->neigh_no;j++){ |
SP |
416 |
k1=vtx[i]->neigh[j-1]; |
7958e9
|
417 |
jp=j+1; |
8f6a69
|
418 |
if(j == vtx[i]->neigh_no) jp=1; |
SP |
419 |
k2=vtx[i]->neigh[jp-1]; |
7958e9
|
420 |
for(k=1;k<=tlist->n;k++){ // VERY NON-OPTIMAL!!! too many loops (vlist.n * vtx.neigh * tlist.n )! |
41a035
|
421 |
k3=tria[k]->vertex[0]; |
SP |
422 |
k4=tria[k]->vertex[1]; |
|
423 |
k5=tria[k]->vertex[2]; |
7958e9
|
424 |
// ts_fprintf(stderr,"%u %u: k=(%u %u %u)\n",k1,k2,k3,k4,k5); |
SP |
425 |
if((vtx[i]==k3 && k1==k4 && k2==k5) || |
|
426 |
(vtx[i]==k4 && k1==k5 && k2==k3) || |
|
427 |
(vtx[i]==k5 && k1==k3 && k2==k4)){ |
b01cc1
|
428 |
|
SP |
429 |
//TODO: probably something wrong with neighbour distribution. |
|
430 |
// if(vtx[i]==k3 || vtx[i]==k4 || vtx[i]==k5){ |
dac2e5
|
431 |
// if(i==6) ts_fprintf(stdout, "Vtx[%u] > Added to tristar!\n",i); |
7958e9
|
432 |
vertex_add_tristar(vtx[i],tria[k]); |
SP |
433 |
} |
|
434 |
} |
|
435 |
} |
|
436 |
/* ts_fprintf(stderr,"TRISTAR for %u (%u):",i-1,vtx[i].tristar_no); |
|
437 |
for(j=0;j<vtx[i].tristar_no;j++){ |
|
438 |
ts_fprintf(stderr," %u,",vtx[i].tristar[j]->idx); |
|
439 |
} |
|
440 |
ts_fprintf(stderr,"\n"); */ |
|
441 |
} |
|
442 |
return TS_SUCCESS; |
|
443 |
} |
|
444 |
|
|
445 |
|
|
446 |
ts_bool init_normal_vectors(ts_triangle_list *tlist){ |
|
447 |
/* Normals point INSIDE vesicle */ |
|
448 |
ts_uint k; |
|
449 |
ts_triangle **tria=tlist->tria -1; //for 0 indexing |
|
450 |
for(k=1;k<=tlist->n;k++){ |
|
451 |
triangle_normal_vector(tria[k]); |
|
452 |
} |
|
453 |
return TS_SUCCESS; |
|
454 |
} |