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