?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Strip extra spaces from code.
[voro++.git] / branches / 2d / src / v_connect.cc
blobbb98fee8101fb8db67482da023fb890dec2b22f5
1
2
3
4 #include "v_connect.hh"
5 #include <math.h>
6 #include <stdio.h>
7
8 namespace voro{
9
10 void v_connect::import(FILE *fp){
11 bool neg_label=false,boundary_track=false,start=false;
12 char *buf(new char[512]);
13 int i=0,id;
14 double x, y,pad=.05;
15
16 while(fgets(buf,512,fp)!=NULL) {
17
18 if(strcmp(buf,"#Start\n")==0||strcmp(buf,"# Start\n")==0) {
19
20 // Check that two consecutive start tokens haven't been
21 // encountered
22 if(boundary_track) voro_fatal_error("File import error - two consecutive start tokens found",VOROPP_FILE_ERROR);
23 start=true;boundary_track=true;
24
25 } else if(strcmp(buf,"#End\n")==0||strcmp(buf,"# End\n")==0||
26 strcmp(buf,"#End")==0||strcmp(buf,"# End")==0) {
27
28 // Check that two consecutive end tokens haven't been
29 // encountered
30 if(start) voro_fatal_error("File import error - end token immediately after start token",VOROPP_FILE_ERROR);
31 if(!boundary_track) voro_fatal_error("File import error - found end token without start token",VOROPP_FILE_ERROR);
32 vbd[i-1]|=2;boundary_track=false;
33 } else {
34 if(!boundary_track && bd==-1) bd=i;
35 // Try and read three entries from the line
36 if(sscanf(buf,"%d %lg %lg",&id,&x,&y)!=3) voro_fatal_error("File import error #1",VOROPP_FILE_ERROR);
37 vid.push_back(id);
38 vpos.push_back(x);
39 vpos.push_back(y);
40 vbd.push_back(start?1:0);
41 i++;
42
43 // Determine bounds
44 if(id<0) neg_label=true;
45 if(id>mid) mid=id;
46 if(x<minx) minx=x;
47 if(x>maxx) maxx=x;
48 if(y<miny) miny=y;
49 if(y>maxy) maxy=y;
50
51 start=false;
52 }
53 }
54
55 if(boundary_track) voro_fatal_error("File import error - boundary not finished",VOROPP_FILE_ERROR);
56 if(!feof(fp)) voro_fatal_error("File import error #2",VOROPP_FILE_ERROR);
57 delete [] buf;
58
59 // Add small amount of padding to container bounds
60 double dx=maxx-minx,dy=maxy-miny;
61 minx-=pad*dx;maxx+=pad*dx;dx+=2*pad*dx;
62 miny-=pad*dy;maxy+=pad*dy;dy+=2*pad*dy;
63
64 // Guess the optimal computationl grid, aiming at eight particles per
65 // grid square
66 double lscale=sqrt(8.0*dx*dy/i);
67 nx=(int)(dx/lscale)+1,ny=(int)(dy/lscale)+1;
68 ng=i;
69 gen_to_vert= new vector<int>[i];
70 gen_to_ed= new vector<int>[i];
71 gen_to_gen_e= new vector<int>[i];
72 gen_to_gen_v=new vector<int>[i];
73 generator_is_vertex=new int[i];
74 for(int j=0;j<ng;j++) generator_is_vertex[j]=-1;
75 mp=new int[mid+1];
76 for(int j=0;j<ng;j++) mp[vid[j]]=j;
77 }
78 // Assemble vert_to_gen,gen_to_vert,vert_to_ed,ed_to_vert
79 void v_connect::assemble_vertex(){
80 bool arrange=false;
81 int cv,lv,cvi,lvi,fvi,j,id,pcurrent_vertices=init_vertices,vert_size=0,g1,g2,g3,gl1,gl2,gl3,i,pne=0;
82 int *pmap;
83 int *ped_to_vert=new int[2*current_edges];
84 bool seencv=false,seenlv=false;
85 vector<int> gens;
86 vector<int> *pvert_to_gen= new vector<int>[pcurrent_vertices];
87 vector<int> *pgen_to_vert=new vector<int>[ng];
88 vector<int> problem_verts;
89 vector<int> problem_verts21;
90 vector<int> problem_verts32;
91 vector<int> problem_gen_to_vert;
92 double gx1,gy1,gx2,gy2;
93 double *pvertl=new double[2*pcurrent_vertices];
94 unsigned int *globvertc=new unsigned int[2*(ng+1)*ng*ng];
95
96 for(int i=0;i<2*((ng+1)*ng*ng);i++){
97 globvertc[i]=0;
98 }
99
100 cout << "2.1" << endl;
101
102 double x,y,vx,vy;
103 // Create container
104 container_boundary_2d con(minx,maxx,miny,maxy,nx,ny,false,false,16);
105
106 // Import data
107 for(j=0;j<vid.size();j++) {
108 if(vbd[j]&1) con.start_boundary();
109 con.put(vid[j],vpos[2*j],vpos[2*j+1]);
110 if(vbd[j]&2) con.end_boundary();
111 }
112
113 // Carry out all of the setup prior to computing any Voronoi cells
114 con.setup();
115 con.full_connect_on();
116 voronoicell_nonconvex_neighbor_2d c;
117 c_loop_all_2d cl(con);
118
119 //Compute Voronoi Cells, adding vertices to potential data structures (p*)
120 if(cl.start()) do if(con.compute_cell(c,cl)){
121 cv=0;
122 lv=c.ed[2*cv+1];
123 id=cl.pid();
124 x=vpos[2*mp[id]];
125 y=vpos[2*mp[id]+1];
126 groom_vertexg(c);
127 do{
128 gens=c.vertexg[cv];
129 vx=.5*c.pts[2*cv]+x;
130 vy=.5*c.pts[2*cv+1]+y;
131 if(gens.size()==1){
132 seencv=false;
133 if(pcurrent_vertices==vert_size){
134 pcurrent_vertices<<=1;
135 add_memory_vector(pvert_to_gen,pcurrent_vertices);
136 add_memory_array(pvertl,pcurrent_vertices);
137
138 }
139 pgen_to_vert[gens[0]].push_back(vert_size);
140 pvertl[2*vert_size]=vx;
141 pvertl[2*vert_size+1]=vy;
142 pvert_to_gen[vert_size]=gens;
143 cvi=vert_size;
144 if(cv==0) fvi=cvi;
145 vert_size++;
146 }else if(gens.size()==2){
147 gx1=vpos[2*mp[gens[0]]]-(c.pts[2*cv]*.5+x);gy1=vpos[2*mp[gens[0]]+1]-(c.pts[2*cv+1]*.5+y);
148 gx2=vpos[2*mp[gens[1]]]-(c.pts[2*cv]*.5+x);gy2=vpos[2*mp[gens[1]]+1]-(c.pts[2*cv+1]*.5+y);
149 if((((gx1*gy2)-(gy1*gx2))>-tolerance) && ((gx1*gy2)-(gy1*gx2))<tolerance){
150 if((vbd[mp[gens[0]]]|2)!=0 && vbd[mp[gens[1]]]==1){
151
152 g1=gens[1]; g2=gens[0];
153 }else if((vbd[mp[gens[1]]]|2)!=0 && vbd[mp[gens[0]]]==1){
154
155 g1=gens[0]; g2=gens[1];
156 }else{
157 if(mp[gens[0]]<mp[gens[1]]){
158 g1=gens[1];g2=gens[0];
159 }else{
160 g1=gens[0];g2=gens[1];
161 }
162 }
163 }else if(((gx1*gy2)-(gy1*gx2))>0){
164 g1=gens[0]; g2=gens[1];
165 }else{
166 g1=gens[1]; g2=gens[0];
167 }
168 if(globvertc[2*((ng*ng*g1)+(ng*g2))]!=0){
169 seencv=true;
170 globvertc[2*((ng*ng*g1)+(ng*g2))]+=1;
171 cvi=globvertc[2*((ng*ng*g1)+(ng*g2))+1];
172 pgen_to_vert[id].push_back(cvi);
173 if(cv==0) fvi=cvi;
174 }else{
175 seencv=false;
176 if(pcurrent_vertices==vert_size){
177 pcurrent_vertices<<=1;
178 add_memory_vector(pvert_to_gen,pcurrent_vertices);
179 add_memory_array(pvertl,pcurrent_vertices);
180
181 }
182 pgen_to_vert[id].push_back(vert_size);
183 pvertl[2*vert_size]=vx;
184 pvertl[2*vert_size+1]=vy;
185 pvert_to_gen[vert_size]=gens;
186 globvertc[2*((ng*ng*g1)+(ng*g2))]=(unsigned int)1;
187 globvertc[2*((ng*ng*g1)+(ng*g2))+1]=vert_size;
188 cvi=vert_size;
189 if(cv==0) fvi=cvi;
190 vert_size++;
191 }
192 }else{
193 arrange_cc_x_to_gen(gens,vx,vy);
194 g1=gens[0];g2=gens[1];g3=gens[2];
195 if(g1<g2 && g1<g3){
196 gl1=g1;gl2=g2;gl3=g3;
197 }else if(g2<g1 && g2<g3){
198 gl1=g2;gl2=g3;gl3=g1;
199 }else{
200 gl1=g3;gl2=g1;gl3=g2;
201 }
202 if(globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)]!=0){
203 seencv=true;
204 globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)]+=1;
205 cvi=globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)+1];
206 pgen_to_vert[id].push_back(cvi);
207 if(cv==0)fvi=cvi;
208 }else{
209 seencv=false;
210 if(pcurrent_vertices==vert_size){
211 pcurrent_vertices<<=1;
212 add_memory_vector(pvert_to_gen,pcurrent_vertices);
213 add_memory_array(pvertl,pcurrent_vertices);
214
215 }
216 pgen_to_vert[id].push_back(vert_size);
217 pvertl[2*vert_size]=vx;
218 pvertl[2*vert_size+1]=vy;
219 pvert_to_gen[vert_size]=gens;
220 globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)]=(unsigned int)1;
221 globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)+1]=vert_size;
222 cvi=vert_size;
223 if(cv==0)fvi=cvi;
224 vert_size++;
225 }
226 }
227
228 //add edges to potential edge structures
229 if(cv!=0){
230 if(pne==current_edges){
231 current_edges<<=1;
232 add_memory_array(ped_to_vert,current_edges);
233 }
234 if(cvi<lvi){
235 ped_to_vert[2*pne]=cvi;
236 ped_to_vert[2*pne+1]=lvi;
237 }else{
238 ped_to_vert[2*pne]=lvi;
239 ped_to_vert[2*pne+1]=cvi;
240 }
241 pne++;
242 }
243
244 lv=cv;
245 lvi=cvi;
246 seenlv=seencv;
247 cv=c.ed[2*lv];
248
249
250 }while(cv!=0);
251 //deal with final edge (last vertex to first vertex)
252 if(pne==current_edges){
253 current_edges<<=1;
254 add_memory_array(ped_to_vert,current_edges);
255 }
256 if(fvi<lvi){
257 ped_to_vert[2*pne]=fvi;
258 ped_to_vert[2*pne+1]=lvi;
259 }else{
260 ped_to_vert[2*pne]=lvi;
261 ped_to_vert[2*pne+1]=fvi;
262 }
263 pne++;
264
265 }while (cl.inc());
266 cout << "2.2" << endl;
267 //Add non-problem vertices(connectivity<=3) to class variables, add problem vertices to problem vertice data structures
268 pmap=new int[pcurrent_vertices];
269 j=0;
270 int v1=0;
271 for(int i=0;i<vert_size;i++){
272 gens=pvert_to_gen[i];
273 if(gens.size()==1){
274 if(j==current_vertices){
275 add_memory_vertices();
276 }
277 generator_is_vertex[gens[0]]=j;
278 vertex_is_generator[j]=gens[0];
279 vertl[2*j]=pvertl[2*i];
280 vertl[2*j+1]=pvertl[2*i+1];
281 vert_to_gen[j]=pvert_to_gen[i];
282 pmap[i]=j;
283 j++;
284 v1++;
285 }
286 else if(gens.size()==2){
287 gx1=vpos[2*mp[gens[0]]]-(pvertl[2*i]);gy1=vpos[2*mp[gens[0]]+1]-(pvertl[2*i+1]);
288 gx2=vpos[2*mp[gens[1]]]-(pvertl[2*i]);gy2=vpos[2*mp[gens[1]]+1]-(pvertl[2*i+1]);
289 if((((gx1*gy2)-(gy1*gx2))>-tolerance) && ((gx1*gy2)-(gy1*gx2))<tolerance){
290 if((vbd[mp[gens[0]]]|2)!=0 && vbd[mp[gens[1]]]==1){
291 g1=gens[1]; g2=gens[0];
292 }else if((vbd[mp[gens[1]]]|2)!=0 && vbd[mp[gens[0]]]==1){
293 g1=gens[0]; g2=gens[1];
294 }else{
295 if(mp[gens[0]]<mp[gens[1]]){
296 g1=gens[1];g2=gens[0];
297 }else{
298 g1=gens[0];g2=gens[1];
299 }
300 }
301 }else if(((gx1*gy2)-(gy1*gx2))>0){
302 g1=gens[0]; g2=gens[1];
303 }else{
304 g1=gens[1]; g2=gens[0];
305 }
306 if(globvertc[2*((ng*ng*g1)+(ng*g2))]!=2){
307 problem_verts21.push_back(i);
308 }
309 else{
310 if(j==current_vertices){
311 add_memory_vertices();
312 }
313 vertl[2*j]=pvertl[2*i];
314 vertl[2*j+1]=pvertl[2*i+1];
315 vert_to_gen[j]=pvert_to_gen[i];
316 pmap[i]=j;
317 j++;
318 }
319 }else{
320 g1=gens[0];g2=gens[1];g3=gens[2];
321 if(g1<g2 && g1<g3){
322 gl1=g1;gl2=g2;gl3=g3;
323 }else if(g2<g1 && g2<g3){
324 gl1=g2;gl2=g3;gl3=g1;
325 }else{
326 gl1=g3;gl2=g1;gl3=g2;
327 }
328 if(globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)]!=3){
329 if(globvertc[2*((ng*ng*gl1)+(ng*gl2)+gl3)]==1) problem_verts.push_back(i);
330 else problem_verts32.push_back(i);
331 }else{
332 if(j==current_vertices){
333 add_memory_vertices();
334 }
335 vertl[2*j]=pvertl[2*i];
336 vertl[2*j+1]=pvertl[2*i+1];
337 vert_to_gen[j]=pvert_to_gen[i];
338 pmap[i]=j;
339 j++;
340 }
341 }
342
343
344 }
345 cout << "2.3\n problemverts21.size()=" << problem_verts21.size() << "\nproblem_verts32.size()=" << problem_verts32.size() << "\nproblem_verts.size()=" << problem_verts.size() << endl;
346 nv=j;
347 degenerate_vertices=j;
348 //deal with problem verts
349 while(problem_verts21.size()>problem_verts32.size()){
350 for(int i=0;i<problem_verts21.size();i++){
351 for(int j=i+1;j<problem_verts21.size();j++){
352 one_in_common(pvert_to_gen[problem_verts21[i]],pvert_to_gen[problem_verts21[j]],g1);
353 if(g1==-1) continue;
354 gens=pvert_to_gen[problem_verts21[i]];
355 gens.push_back(not_this_one(pvert_to_gen[problem_verts21[j]],g1));
356 for(int k=0;k<problem_verts.size();k++){
357 if(contain_same_elements(gens,pvert_to_gen[problem_verts[k]])){
358 if(nv==current_vertices) add_memory_vertices();
359 vertl[2*nv]=pvertl[2*problem_verts[k]];
360 vertl[2*nv+1]=pvertl[2*problem_verts[k]+1];
361 vert_to_gen[nv]=pvert_to_gen[problem_verts[k]];
362 pmap[problem_verts21[i]]=nv;
363 pmap[problem_verts21[j]]=nv;
364 pmap[problem_verts[k]]=nv;
365 nv++;
366 problem_gen_to_vert.push_back(problem_verts21[i]);
367 problem_gen_to_vert.push_back(problem_verts21[j]);
368 problem_verts21.erase(problem_verts21.begin()+i);
369 problem_verts21.erase(problem_verts21.begin()+(j-1));
370 problem_verts.erase(problem_verts.begin()+k);
371 break;
372 }
373 }
374 }
375 }
376 }
377 cout << "part 2" << endl;
378 while(problem_verts32.size()>0){
379 for(int i=0;i<problem_verts32.size();i++){
380 gens=pvert_to_gen[problem_verts32[i]];
381 for(int j=0;j<problem_verts21.size();j++){
382 if(subset(pvert_to_gen[problem_verts21[j]],gens)){
383 if(nv==current_vertices) add_memory_vertices();
384 vertl[2*nv]=pvertl[2*problem_verts32[i]];
385 vertl[2*nv+1]=pvertl[2*problem_verts32[i]+1];
386 vert_to_gen[nv]=gens;
387 pmap[problem_verts32[i]]=nv;
388 pmap[problem_verts21[j]]=nv;
389 nv++;
390 problem_gen_to_vert.push_back(problem_verts21[j]);
391 problem_verts21.erase(problem_verts21.begin()+j);
392 problem_verts32.erase(problem_verts32.begin()+i);
393 break;
394 }
395 }
396 }
397 }
398
399
400 double standard,distance;
401 cout << "part3" << endl;
402 while(problem_verts.size()>0){
403
404 if(nv==current_vertices) add_memory_vertices();
405 gens=pvert_to_gen[problem_verts[0]];
406 vx=pvertl[2*problem_verts[0]];vy=pvertl[2*problem_verts[0]+1];
407 standard=pow(vx-vpos[2*mp[gens[0]]],2)+pow(vy-vpos[2*mp[gens[0]]+1],2);
408 g1=gens[0];g2=gens[1];
409 pmap[problem_verts[0]]=nv;
410 problem_verts.erase(problem_verts.begin());
411 i=0;
412 while(true){
413 g3=not_these_two(pvert_to_gen[problem_verts[i]],g1,g2);
414 if(g3!=-1){
415 distance=pow(vx-vpos[2*mp[g3]],2)+pow(vy-vpos[2*mp[g3]+1],2);
416 }
417 if(contains_two(pvert_to_gen[problem_verts[i]],g1,g2) &&
418 (distance<(standard+tolerance) && (distance>(standard-tolerance)))){
419 if(g3==gens[0]){
420 pmap[problem_verts[i]]=nv;
421 problem_verts.erase(problem_verts.begin()+i);
422 break;
423 }
424 if(g3!=gens[2]) gens.push_back(g3);
425 pmap[problem_verts[i]]=nv;
426 g2=g3;
427 g1=pvert_to_gen[problem_verts[i]][0];
428 if(problem_verts.size()>1){
429 problem_verts.erase(problem_verts.begin()+i);
430 }else{
431 break;
432 }
433 i=0;
434 }else{
435 i++;
436 }
437 }
438
439 vertl[2*nv]=vx;
440 vertl[2*nv+1]=vy;
441 arrange_cc_x_to_gen(gens,vx,vy);
442 vert_to_gen[nv]=gens;
443 nv++;
444
445 }
446 delete [] pvert_to_gen;
447 delete [] pvertl;
448 delete [] globvertc;
449 cout << "2.4" << endl;
450 //assemble edge data structures
451 ed_to_vert=new int[2*pne];
452 vert_to_ed=new vector<int>[nv];
453 unsigned int* globedgec=new unsigned int[nv*nv];
454 for(int i=0;i<(nv*nv);i++){
455 globedgec[i]=0;
456 }
457 for(int i=0;i<pne;i++){
458 g1=pmap[ped_to_vert[2*i]];g2=pmap[ped_to_vert[2*i+1]];
459 if(g2<g1){
460 g2^=g1;
461 g1^=g2;
462 g2^=g1;
463 }
464 if(globedgec[(nv*g1+g2)]!=0){
465 continue;
466 }else{
467 globedgec[(nv*g1+g2)]=1;
468 ed_to_vert[2*ne]=g1;
469 ed_to_vert[2*ne+1]=g2;
470 vert_to_ed[g1].push_back(ne);
471 vert_to_ed[g2].push_back(ne);
472 ne++;
473 }
474 }
475 for(int i=0;i<ng;i++){
476 for(int k=0;k<pgen_to_vert[i].size();k++){
477 gen_to_vert[i].push_back(pmap[pgen_to_vert[i][k]]);
478 if(contains(problem_gen_to_vert,pgen_to_vert[i][k])) arrange=true;
479 }
480 if(arrange) arrange_cc_gen_to_vert(gen_to_vert[i],vpos[2*mp[i]],vpos[2*mp[i]+1]);
481 arrange=false;
482 }
483
484 delete [] pgen_to_vert;
485 delete [] pmap;
486 delete [] ped_to_vert;
487 delete [] globedgec;
488 cout << "out" << endl;
489 }
490 //assemble gen_to_gen , gen_to_ed , ed_to_gen
491 void v_connect::assemble_gen_ed(){
492 cout << "gen_ed 1" << endl;
493 ed_to_gen=new vector<int>[ne];
494 //though neither ed_on_bd or vert_on_bd are modified during this method, they are initialized here in case the user
495 //does not require boundary information and will not be calling assemble_boundary
496 ed_on_bd=new vector<int>[ne];
497 vert_on_bd=new vector<int>[nv];
498 vector<int> gens;
499 int g1,g2,v1,v2,j,vi;
500 double gx1,gy1,gx2,gy2,vx1,vy1,vx2,vy2;
501 cout << "gen_ed 2" << endl;
502 for(int i=0;i<ne;i++){
503
504 v1=ed_to_vert[2*i];v2=ed_to_vert[2*i+1];
505 two_in_common(vert_to_gen[v1],vert_to_gen[v2],g1,g2);
506 if(g1!=-1 && g2!=-1){
507 gen_to_gen_e[g1].push_back(g2);
508 gen_to_gen_e[g2].push_back(g1);
509 for(int k=0;k<gen_to_vert[g1].size();k++){
510 if(gen_to_vert[g1][k]==v1){
511 vi=k;
512 break;
513 }
514 }
515 if((vi!=(gen_to_vert[g1].size()-1) && gen_to_vert[g1][vi+1]==v2) ||
516 (vi==(gen_to_vert[g1].size()-1) && gen_to_vert[g1][0]==v2)){
517 ed_to_gen[i].push_back(g1);
518 ed_to_gen[i].push_back(g2);
519 gen_to_ed[g1].push_back(i);
520 gen_to_ed[g2].push_back(~i);
521 }else{
522 ed_to_gen[i].push_back(g2);
523 ed_to_gen[i].push_back(g1);
524 gen_to_ed[g2].push_back(i);
525 gen_to_ed[g1].push_back(~i);
526 }
527 }else{
528 if(g1==-1) continue;
529 ed_to_gen[i].push_back(g1);
530 vx1=vertl[2*v1];vy1=vertl[2*v1+1];
531 vx2=vertl[2*v2];vy2=vertl[2*v2+1];
532 gx1=vpos[2*mp[g1]];gy1=vpos[2*mp[g1]+1];
533 for(int k=0;k<gen_to_vert[g1].size();k++){
534 if(gen_to_vert[g1][k]==v1){
535 vi=k;
536 break;
537 }
538 }
539 if((vi!=(gen_to_vert[g1].size()-1) && gen_to_vert[g1][vi+1]==v2) ||
540 (vi==(gen_to_vert[g1].size()-1) && gen_to_vert[g1][0]==v2)) gen_to_ed[g1].push_back(i);
541 else gen_to_ed[g1].push_back(~i);
542 }
543 }
544 cout << "gen_ed 3" << endl;
545 for(int i=0;i<nv;i++){
546 if(vert_to_ed[i].size()<=3) continue;
547 else{
548 gens=vert_to_gen[i];
549 for(int j=0;j<gens.size();j++){
550 for(int k=0;k<gens.size();k++){
551 if(!contains(gen_to_gen_e[gens[j]],gens[k])) gen_to_gen_v[gens[j]].push_back(gens[k]);
552
553 }
554 }
555 }
556 }
557 cout << "gen_ed 4" << endl;
558 //arrange gen_to_ed gen_to_gen_e gen_to_gen_v counterclockwise
559 for(int i=0;i<ng;i++){
560 gx1=vpos[2*mp[i]];gy1=vpos[2*mp[i]+1];
561 arrange_cc_gen_to_ed(gen_to_ed[i]);
562 arrange_cc_x_to_gen(gen_to_gen_e[i],gx1,gy1);
563 arrange_cc_x_to_gen(gen_to_gen_v[i],gx1,gy1);
564 }
565 cout << "gen_ed 5" << endl;
566 //arrange vert_to_ed cc
567 for(int i=0;i<nv;i++){
568 gx1=vertl[2*i];gy1=vertl[2*i+1];
569 arrange_cc_vert_to_ed(vert_to_ed[i],gx1,gy1,i);
570 }
571
572 }
573
574 //assemble vert_on_bd and ed_on_bd as well as side edge information if neccessary.
575 void v_connect::assemble_boundary(){
576 bool begun=false;
577 int i=0,cg,ng,fg,lv,cv,nv,ev,ei,j;
578 while(true){
579 if(vbd[i]==1){
580 begun=true;
581 fg=mp[i];
582 cg=fg;
583 ng=mp[i+1];
584 ev=gen_to_vert[ng][0];
585 }else if(vbd[i]==2){
586 begun=false;
587 cg=mp[i];
588 ng=fg;
589 ev=gen_to_vert[ng][0];
590 }else{
591 if(!begun) break;
592 cg=mp[i];
593 ng=mp[i+1];
594 ev=gen_to_vert[ng][0];
595 }
596 cv=gen_to_vert[cg][0];
597 nv=gen_to_vert[cg][1];
598 one_in_common(vert_to_ed[cv],vert_to_ed[nv],ei);
599 vert_on_bd[cv].push_back(cg);
600 vert_on_bd[cv].push_back(ng);
601 vert_on_bd[nv].push_back(cg);
602 vert_on_bd[nv].push_back(ng);
603 ed_on_bd[ei].push_back(cg);
604 ed_on_bd[ei].push_back(ng);
605 lv=cv;
606 cv=nv;
607 while(nv!=ev){
608 j=0;
609 while(true){
610 if(ed_to_vert[2*vert_to_ed[cv][j]]==lv) break;
611 else if(ed_to_vert[2*vert_to_ed[cv][j]+1]==lv) break;
612 j++;
613 }
614 if(j==(vert_to_ed[cv].size()-1)) ei=vert_to_ed[cv][0];
615 else ei=vert_to_ed[cv][j+1];
616 if(ed_to_vert[2*ei]==cv) nv=ed_to_vert[2*ei+1];
617 else nv=ed_to_vert[2*ei];
618 vert_on_bd[nv].push_back(cg);
619 vert_on_bd[nv].push_back(ng);
620 ed_on_bd[ei].push_back(cg);
621 ed_on_bd[ei].push_back(ng);
622 lv=cv;
623 cv=nv;
624 }
625 i++;
626 }
627
628 }
629
630
631
632 //(x,y)=my coordinates
633 //(vx,vy)=vertex coordinates
634 vector<int> v_connect::groom_vertexg_help(double x,double y,double vx, double vy,vector<int> &g){
635 if(g.size()<2) return g;
636 bool rightside=false;
637 int g0=g[0],g1,bestg,besti;
638 double d1;
639 double standard=pow((vpos[2*mp[g0]]-vx),2)+pow((vpos[2*mp[g0]+1]-vy),2);
640 double gx0,gy0,gx1,gy1,best,current;
641 vector<int> newg,temp;
642 temp.push_back(g0);
643 for(int i=1;i<g.size();i++){
644 g1=g[i];
645 if(contains(temp,g1)) continue;
646 d1=pow((vpos[2*mp[g1]]-vx),2)+pow((vpos[2*mp[g1]+1]-vy),2);
647 if(d1<(standard+tolerance)){
648 temp.push_back(g1);
649 }
650 }
651 if(temp.size()<3) return temp;
652 gx0=vpos[2*mp[g0]]; gy0=vpos[2*mp[g0]+1];
653 newg.push_back(g0);
654 //find right hand
655 g1=temp[1];
656 gx1=vpos[2*mp[g1]]; gy1=vpos[2*mp[g1]+1];
657 if(cross_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy)>0){
658 rightside=true;
659 }
660 best=dot_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy);
661 bestg=g1;
662 besti=1;
663 for(int i=2;i<temp.size();i++){
664 g1=temp[i];
665 if(contains(newg,g1)) continue;
666 gx1=vpos[2*mp[g1]]; gy1=vpos[2*mp[g1]+1];
667 if(cross_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy)>0){
668 if(!rightside){
669 rightside=true;
670 best=dot_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy);
671 bestg=g1;
672 besti=i;
673 }else{
674 current=dot_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy);
675 if(current>best){
676 best=current;
677 bestg=g1;
678 besti=i;
679 }
680 }
681 }else{
682 if(rightside) continue;
683 else{
684 current=dot_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy);
685 if(current<best){
686 best=current;
687 bestg=g1;
688 besti=i;
689 }
690 }
691 }
692 }
693 if(!contains(newg,bestg)) newg.push_back(bestg);
694
695 //find left hand
696 rightside=false;
697 g1=temp[1];
698 gx1=vpos[2*mp[g1]]; gy1=vpos[2*mp[g1]+1];
699 if(cross_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy)<0){
700 rightside=true;
701 }
702 best=dot_product(gx0-vx,gy0-vy,gx1-vy,gy1-vy);
703 bestg=g1;
704 for(int i=2;i<temp.size();i++){
705 g1=temp[i];
706 if(contains(newg,g1)) continue;
707 gx1=vpos[2*mp[g1]]; gy1=vpos[2*mp[g1]+1];
708 if(cross_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy)<0){
709 if(!rightside){
710 rightside=true;
711 best=dot_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy);
712 bestg=g1;
713 }else{
714 current=dot_product(gx0-vx,gy0-vy,gx1-vx,gy1-vy);
715 if(current>best){
716 best=current;
717 bestg=g1;
718 }
719 }
720 }else{
721 if(rightside) continue;
722 else{
723 current=dot_product(gx0-vx,gy0-vx,gx1-vx,gy1-vy);
724 if(current<best){
725 best=current;
726 bestg=g1;
727 }
728 }
729 }
730 }
731 if(!contains(newg,bestg)) newg.push_back(bestg);
732 if(newg.size()<3) return groom_vertexg_help2(x,y,vx,vy,g);
733 return newg;
734
735 }
736 vector<int> v_connect::groom_vertexg_help2(double x,double y,double vx,double vy,vector<int> &g){
737 if(g.size()<2) return g;
738 int m0=g[0],m1=g[1],m2,p,i=1;
739 double d0=pow((vpos[2*mp[m0]]-x),2)+pow((vpos[2*mp[m0]+1]-y),2);
740 double d1=pow((vpos[2*mp[m1]]-vx),2)+pow((vpos[2*mp[m1]+1]-vy),2);
741 double d2;
742 double standard=pow((vpos[2*mp[m0]]-vx),2)+pow((vpos[2*mp[m0]+1]-vy),2);
743 double dp,dcompare, temp;
744 vector<int> newg;
745 while(d1>=standard+tolerance){
746 if(i==g.size()-1){
747 newg.push_back(m0);
748 return newg;
749 }
750 i++;
751 m1=g[i];
752 d1=pow((vpos[2*mp[m1]]-vx),2)+pow((vpos[2*mp[m1]+1]-vy),2);
753 }
754 if(i==g.size()-1){
755 newg.push_back(m0);
756 newg.push_back(m1);
757 return newg;
758 }
759 i++;
760 m2=g[i];
761 d2=pow((vpos[2*mp[m2]]-vx),2)+pow((vpos[2*mp[m2]+1]-vy),2);
762 while(d2>=standard+tolerance){
763 if(i==g.size()-1){
764 newg.push_back(m0);
765 newg.push_back(m1);
766 return newg;
767 }
768 i++;
769 m2=g[i];
770 d2=pow((vpos[2*mp[m2]]-vx),2)+pow((vpos[2*mp[m2]+1]-vy),2);
771 }
772 if(i==g.size()-1){
773 newg.push_back(m0);
774 newg.push_back(m1);
775 newg.push_back(m2);
776 return newg;
777 }
778 i++;
779 d1=pow((vpos[2*mp[m1]]-x),2)+pow((vpos[2*mp[m1]+1]-y),2);
780 d2=pow((vpos[2*mp[m2]]-x),2)+pow((vpos[2*mp[m2]+1]-y),2);
781 if(d0<d2 && d2<d1){
782 temp=d1;
783 d1=d2;
784 d2=temp;
785 m2^=m1;
786 m1^=m2;
787 m2^=m1;
788 }else if(d1<d0 && d0<d2){
789 temp=d1;
790 d1=d0;
791 d0=temp;
792 m0^=m1;
793 m1^=m0;
794 m0^=m1;
795 }else if(d1<d2 && d2<d0){
796 temp=d1;
797 d1=d0;
798 d0=temp;
799 m0^=m1;
800 m1^=m0;
801 m0^=m1;
802 temp=d1;
803 d1=d2;
804 d2=temp;
805 m2^=m1;
806 m1^=m2;
807 m2^=m1;
808 }else if(d2<d1 && d1<d0){
809 temp=d2;
810 d2=d0;
811 d0=temp;
812 m0^=m2;
813 m2^=m0;
814 m0^=m2;
815 }else if(d2<d0 && d0<d1){
816 temp=d2;
817 d2=d0;
818 d0=temp;
819 m0^=m2;
820 m2^=m0;
821 m0^=m2;
822 temp=d1;
823 d1=d2;
824 d2=temp;
825 m2^=m1;
826 m1^=m2;
827 m2^=m1;
828 }
829 for(int j=i;j<g.size();j++){
830 p=g[j];
831 dcompare=pow((vpos[2*mp[p]]-vx),2)+pow((vpos[2*mp[p]+1]-vy),2);
832 if(dcompare<=(standard+tolerance)){
833 dp=pow((vpos[2*mp[p]]-x),2)+pow((vpos[2*mp[p]+1]-y),2);
834 if(dp<d2){
835 temp=d2;
836 d2=dp;
837 dp=temp;
838 p^=m2;
839 m2^=p;
840 p^=m2;
841 if(d2<d1){
842 temp=d1;
843 d1=d2;
844 d2=temp;
845 m2^=m1;
846 m1^=m2;
847 m2^=m1;
848 if(d1<d0){
849 temp=d1;
850 d1=d0;
851 d0=temp;
852 m0^=m1;
853 m1^=m0;
854 m0^=m1;
855 }
856 }
857 }
858 }
859 }
860 newg.push_back(m0);
861 newg.push_back(m1);
862 newg.push_back(m2);
863 return newg;
864 }
865
866 void v_connect::arrange_cc_x_to_gen(vector<int> &list,double cx,double cy){
867 if(list.size()==0) return;
868 bool wrongside;
869 int g1,ng,ni;
870 double x1,y1,x2,y2,best,current;
871 vector<int> newlist;
872 vector<int> potential;
873 newlist.push_back(list[0]);
874 x1=vpos[2*mp[list[0]]];y1=vpos[2*mp[list[0]]+1];
875 list.erase(list.begin());
876 while(list.size()>0){
877 wrongside=true;
878 for(int i=0;i<list.size();i++){
879 g1=list[i];
880 x2=vpos[2*mp[g1]];y2=vpos[2*mp[g1]+1];
881 if(cross_product(cx-x1,cy-y1,cx-x2,cy-y2)>=0){
882 current=dot_product(cx-x1,cy-y1,cx-x2,cy-y2);
883 if(wrongside){
884 ng=g1;
885 ni=i;
886 best=current;
887 wrongside=false;
888 }else if(current>best){
889 best=current;
890 ng=g1;
891 ni=i;
892 }
893 }else{
894 if(!wrongside) continue;
895 current=dot_product(cx-x1,cy-y1,cx-x2,cy-y2);
896 if(i==0){
897 best=current;
898 ng=g1;
899 ni=i;
900 }else if(current<best){
901 best=current;
902 ng=g1;
903 ni=i;
904 }
905 }
906 }
907 newlist.push_back(ng);
908 list.erase(list.begin()+ni);
909 }
910 list=newlist;
911 }
912
913 void v_connect::arrange_cc_gen_to_vert(vector<int> &list,double cx,double cy){
914 if(list.size()==0) return;
915 bool wrongside;
916 int g1,ng,ni;
917 double x1,y1,x2,y2,best,current;
918 vector<int> newlist;
919 vector<int> potential;
920 newlist.push_back(list[0]);
921 x1=vertl[2*list[0]];y1=vertl[2*list[0]+1];
922 list.erase(list.begin());
923 while(list.size()>0){
924 wrongside=true;
925 for(int i=0;i<list.size();i++){
926 g1=list[i];
927 x2=vertl[2*g1];y2=vertl[2*g1+1];
928 if(cross_product(cx-x1,cy-y1,cx-x2,cy-y2)>=0){
929 current=dot_product(cx-x1,cy-y1,cx-x2,cy-y2);
930 if(wrongside){
931 ng=g1;
932 ni=i;
933 best=current;
934 wrongside=false;
935 }else if(current>best){
936 best=current;
937 ng=g1;
938 ni=i;
939 }
940 }else{
941 if(!wrongside) continue;
942 current=dot_product(cx-x1,cy-y1,cx-x2,cy-y2);
943 if(i==0){
944 best=current;
945 ng=g1;
946 ni=i;
947 }else if(current<best){
948 best=current;
949 ng=g1;
950 ni=i;
951 }
952 }
953 }
954 newlist.push_back(ng);
955 list.erase(list.begin()+ni);
956 }
957 list=newlist;
958 }
959
960 void v_connect::arrange_cc_gen_to_ed(vector<int> &list){
961 vector<int> newlist;
962 int v1,v2,ed,i=0;
963 newlist.push_back(list[0]);
964 list.erase(list.begin());
965 if(newlist[0]<0){
966 ed=~newlist[0];
967 v1=ed_to_vert[2*ed];
968 }else{
969 ed=newlist[0];
970 v1=ed_to_vert[2*ed+1];
971 }
972 while(list.size()>0){
973 if(list[i]>=0){
974 v2=ed_to_vert[2*list[i]];
975 if(v2==v1){
976 v1=ed_to_vert[2*list[i]+1];
977 newlist.push_back(list[i]);
978 list.erase(list.begin()+i);
979 i=0;
980 }else i++;
981 }else{
982 v2=ed_to_vert[2*(~list[i])+1];
983
984 if(v2==v1){
985 v1=ed_to_vert[2*(~list[i])];
986 newlist.push_back(list[i]);
987 list.erase(list.begin()+i);
988 i=0;
989 }else i++;
990 }
991 }
992 list=newlist;
993 }
994
995 void v_connect::arrange_cc_vert_to_ed(vector<int> &list,double cx, double cy,int id){
996
997 if(list.size()==0) return;
998 bool wrongside;
999 int g1,ng,ni,index;
1000 double x1,y1,x2,y2,best,current;
1001 vector<int> newlist;
1002 vector<int> potential;
1003 newlist.push_back(list[0]);
1004 if(ed_to_vert[2*list[0]]==id) index=ed_to_vert[2*list[0]+1];
1005 else index=ed_to_vert[2*list[0]];
1006 x1=vertl[2*index];y1=vertl[2*index+1];
1007 list.erase(list.begin());
1008 while(list.size()>0){
1009 wrongside=true;
1010 for(int i=0;i<list.size();i++){
1011 g1=list[i];
1012 if(ed_to_vert[2*g1]==id) index=ed_to_vert[2*g1+1];
1013 else index=ed_to_vert[2*g1];
1014 x2=vertl[2*index];y2=vertl[2*index+1];
1015 if(cross_product(cx-x1,cy-y1,cx-x2,cy-y2)>=0){
1016 current=dot_product(cx-x1,cy-y1,cx-x2,cy-y2);
1017 if(wrongside){
1018 ng=g1;
1019 ni=i;
1020 best=current;
1021 wrongside=false;
1022 }else if(current>best){
1023 best=current;
1024 ng=g1;
1025 ni=i;
1026 }
1027 }else{
1028 if(!wrongside) continue;
1029 current=dot_product(cx-x1,cy-y1,cx-x2,cy-y2);
1030 if(i==0){
1031 best=current;
1032 ng=g1;
1033 ni=i;
1034 }else if(current<best){
1035 best=current;
1036 ng=g1;
1037 ni=i;
1038 }
1039 }
1040 }
1041 newlist.push_back(ng);
1042 list.erase(list.begin()+ni);
1043 }
1044 list=newlist;
1045 }
1046
1047 void v_connect::draw_gnu(FILE *fp){
1048 int vert;
1049 for(int i=0;i<ng;i++){
1050 fprintf(fp,"# cell number %i\n",i);
1051 for(int j=0;j<gen_to_vert[i].size();j++){
1052 vert=gen_to_vert[i][j];
1053 fprintf(fp,"%g %g\n",vertl[2*vert],vertl[2*vert+1]);
1054 }
1055 fprintf(fp,"%g %g\n",vertl[2*gen_to_vert[i][0]],vertl[2*gen_to_vert[i][0]+1]);
1056 fprintf(fp,"\n");
1057 }
1058 }
1059
1060
1061
1062
1063 void v_connect::draw_vtg_gnu(FILE *fp){
1064 double vx, vy, gx, gy;
1065 int l2;
1066 for(int i=0;i<ng;i++){
1067 gx=vpos[2*mp[i]];gy=vpos[2*mp[i]+1];
1068 for(int k=0;k<gen_to_vert[i].size();k++){
1069 vx=vertl[2*gen_to_vert[i][k]];vy=vertl[2*gen_to_vert[i][k]+1];
1070 fprintf(fp, "%g %g\n %g %g\n\n\n",vx,vy,gx,gy);
1071 }
1072 }
1073 }
1074
1075 void v_connect::draw_gen_gen(FILE *fp){
1076 int g2;
1077 double gx1,gy1,gx2,gy2;
1078 for(int i=0;i<ng;i++){
1079 gx1=vpos[2*mp[i]];gy1=vpos[2*mp[i]+1];
1080 for(int k=0;k<gen_to_gen_e[i].size();k++){
1081 g2=gen_to_gen_e[i][k];
1082 gx2=vpos[2*mp[g2]];gy2=vpos[2*mp[g2]+1];
1083 fprintf(fp, "%g %g\n %g %g\n\n\n",gx1,gy1,gx2,gy2);
1084 }
1085 for(int k=0;k<gen_to_gen_v[i].size();k++){
1086 g2=gen_to_gen_v[i][k];
1087 gx2=vpos[2*mp[g2]];gy2=vpos[2*mp[g2]+1];
1088 fprintf(fp, "%g %g\n %g %g \n\n\n",gx1,gy1,gx2,gy2);
1089 }
1090 }
1091 }
1092
1093 void v_connect::label_vertices(FILE *fp){
1094 double vx,vy;
1095 for(int i=0;i<nv;i++){
1096 vx=vertl[2*i];vy=vertl[2*i+1];
1097 fprintf(fp,"set label '%i' at %g,%g point lt 2 pt 3 ps 2 offset -3,3\n",i,vx,vy);
1098 }
1099 }
1100
1101 void v_connect::label_generators(FILE *fp){
1102 double gx,gy;
1103 for(int i=0;i<vid.size();i++){
1104 gx=vpos[2*i];gy=vpos[2*i+1];
1105 fprintf(fp,"set label '%i' at %g,%g point lt 4 pt 4 ps 2 offset 3,-3\n",vid[i],gx,gy);
1106 }
1107 }
1108
1109 void v_connect::label_edges(FILE *fp){
1110 double ex,ey,vx1,vy1,vx2,vy2;
1111 for(int i=0;i<ne;i++){
1112 vx1=vertl[2*ed_to_vert[2*i]];vy1=vertl[2*ed_to_vert[2*i]+1];
1113 vx2=vertl[2*ed_to_vert[2*i+1]];vy2=vertl[2*ed_to_vert[2*i+1]+1];
1114 ex=(vx1+vx2)/2; ey=(vy1+vy2)/2;
1115 fprintf(fp,"set label '%i' at %g,%g point lt 3 pt 1 ps 2 offset 0,-3\n",i,ex,ey);
1116 }
1117 }
1118
1119 void v_connect::label_centroids(FILE *fp){
1120 double x,y;
1121 for(int i=0;i<ng;i++){
1122 centroid(i,x,y);
1123 fprintf(fp,"set label '%i' at %g,%g point lt 5 pt 5 ps 2 offset 0,3\n",i,x,y);
1124 }
1125 }
1126
1127 void v_connect::print_gen_to_ed_table(FILE *fp){
1128 fprintf(fp,"generator to edge connectivity, arranged counterclockwise. Negative edge number means edge with reverse orientation\n\n");
1129 for(int i=0;i<ng;i++){
1130 fprintf(fp,"generator %i\n", i);
1131 for(int k=0;k<gen_to_ed[i].size();k++){
1132 if(gen_to_ed[i][k]>=0) fprintf(fp,"\t %i\n",gen_to_ed[i][k]);
1133 else fprintf(fp,"\t -%i\n",~gen_to_ed[i][k]);
1134 }
1135 fprintf(fp,"\n\n");
1136 }
1137 }
1138
1139 void v_connect::print_gen_to_vert_table(FILE *fp){
1140 fprintf(fp,"generator to vertex connectivity, arranged conterclockwise\n\n");
1141 for(int i=0;i<ng;i++){
1142 fprintf(fp,"generator %i\n",i);
1143 for(int k=0;k<gen_to_vert[i].size();k++){
1144 fprintf(fp,"\t %i\n",gen_to_vert[i][k]);
1145 }
1146 }
1147 }
1148
1149 void v_connect::print_vert_to_gen_table(FILE *fp){
1150 fprintf(fp,"vertex to generator connectivity, arranged counterclockwise\n\n");
1151 for(int i=0;i<nv;i++){
1152 fprintf(fp,"vertex %i\n",i);
1153 for(int k=0;k<vert_to_gen[i].size();k++){
1154 fprintf(fp,"\t %i\n",vert_to_gen[i][k]);
1155 }
1156 }
1157 }
1158
1159 void v_connect::print_ed_to_gen_table(FILE *fp){
1160 fprintf(fp,"edge to generator connectivity, arranged left-side, right-side\n\n");
1161 for(int i=0;i<ne;i++){
1162 fprintf(fp,"ed %i\n",i);
1163 for(int k=0;k<ed_to_gen[i].size();k++){
1164 fprintf(fp,"\t %i\n",ed_to_gen[i][k]);
1165 }
1166 }
1167 }
1168
1169 void v_connect::print_vert_to_ed_table(FILE *fp){
1170 fprintf(fp,"vert to edge connectivity, arranged cc\n\n");
1171 for(int i=0;i<nv;i++){
1172 fprintf(fp,"vert %i\n",i);
1173 for(int k=0;k<vert_to_ed[i].size();k++){
1174 fprintf(fp,"\t %i\n",vert_to_ed[i][k]);
1175 }
1176 }
1177 }
1178
1179 void v_connect::print_vert_boundary(FILE *fp){
1180 fprintf(fp,"vertex on boundary\n\n");
1181 for(int i=0;i<nv;i++){
1182 fprintf(fp,"\nvert %i\n",i);
1183 if(vert_on_bd[i].size()==0) fprintf(fp,"\t vertex not on bound\n");
1184 else{
1185 fprintf(fp,"\tvertex on bound");
1186 for(int k=0;k<vert_on_bd[i].size();k++) fprintf(fp,"\t %i",vert_on_bd[i][k]);
1187 }
1188 }
1189 }
1190
1191 void v_connect::print_ed_boundary(FILE *fp){
1192 fprintf(fp,"edge on boundar \n\n");
1193 for(int i=0;i<ne;i++){
1194 fprintf(fp,"\nedge %i\n",i);
1195 if(ed_on_bd[i].size()==0) fprintf(fp,"\t edge not on bound\n");
1196 else{
1197 fprintf(fp,"\t edge on bound");
1198 for(int k=0;k<ed_on_bd[i].size();k++) fprintf(fp,"\t %i",ed_on_bd[i][k]);
1199 }
1200 }
1201 }
1202
1203 void v_connect::ascii_output(FILE *fp){
1204 fprintf(fp,"\n# General Mesh Data\nNnp\t%i\nNel\t%i\nNel_tri3\t0\nNel_poly2d\t%i\nNel_quad4\t0\nNel_hexh8\t0\nNel_poly3d\t0\nNel_pyr5\t0\nNel_tet4\t0\nNel_wedget\t0\nNdim\t2\nNnd_sets\t0\nNsd_sets\t0\nendi\n",nv,ng,ng);
1205 fprintf(fp,"# element data: global id, block id, number of nodes, nodes\n");
1206 for(int i=0;i<ng;i++){
1207 fprintf(fp,"%i\t1\t%i",i,gen_to_vert[i].size());
1208 for(int j=0;j<gen_to_vert[i].size();j++){
1209 fprintf(fp,"\t%i",gen_to_vert[i][j]);
1210 }
1211 fprintf(fp,"\n");
1212 }
1213 fprintf(fp,"#\n#nodal data:global id, xcoord, ycoord\n#\n");
1214 for(int i=0;i<nv;i++){
1215 fprintf(fp,"%i\t%g\t%g\n",i,vertl[2*i],vertl[2*i+1]);
1216 }
1217 }
1218
1219
1220
1221 void v_connect::add_memory_vertices(){
1222 double *vertle(vertl+2*current_vertices);
1223 int *vertex_is_generatore(vertex_is_generator+current_vertices);
1224 current_vertices<<=1;
1225 cout << "2.2.1" << endl;
1226 //copy vertl
1227 double *nvertl(new double[2*current_vertices]),*nvp(nvertl),*vp(vertl);
1228 while(vp<vertle) *(nvp++)=*(vp++);
1229 delete [] vertl;vertl=nvertl;
1230 cout << "2.2.2" << endl;
1231 //copy vert_to_gen, vert_to_ed
1232 vector<int> *nvert_to_gen(new vector<int>[current_vertices]);
1233 for(int i=0;i<(current_vertices>>1);i++){
1234 nvert_to_gen[i]=vert_to_gen[i];
1235 }
1236 delete [] vert_to_gen;vert_to_gen=nvert_to_gen;
1237 cout << "2.2.3" << endl;
1238 //copy vertex_is_generator
1239 int *nvertex_is_generator(new int[current_vertices]),*nvig(nvertex_is_generator),*vig(vertex_is_generator),*nvige(nvertex_is_generator+current_vertices);
1240 while(vig<vertex_is_generatore) *(nvig++)=*(vig++);
1241 cout << "inbetween" << endl;
1242 while(nvig<nvige) *(nvig++)=-1;
1243 cout << "2.2.4" << endl;
1244 delete [] vertex_is_generator; vertex_is_generator=nvertex_is_generator;
1245 }
1246
1247 //returns the signed area of the cell corresponding to generator g
1248 double v_connect::signed_area(int g){
1249 double area=0,vx1,vy1,vx2,vy2;
1250 vx1=vertl[2*gen_to_vert[g][0]];
1251 vy1=vertl[2*gen_to_vert[g][0]+1];
1252 for(int i=0;i<gen_to_vert[g].size();i++){
1253 if(i==gen_to_vert[g].size()-1){
1254 vx2=vertl[2*gen_to_vert[g][0]];
1255 vy2=vertl[2*gen_to_vert[g][0]+1];
1256 }else{
1257 vx2=vertl[2*gen_to_vert[g][i+1]];
1258 vy2=vertl[2*gen_to_vert[g][i+1]+1];
1259 }
1260 area+=((vx1*vy2)-(vx2*vy1));
1261 vx1=vx2;
1262 vy1=vy2;
1263 }
1264 area*=.5;
1265 return area;
1266 }
1267
1268 //returns the centroid of the cell corresponding to generator g in x,y
1269 void v_connect::centroid(int g,double &x,double &y){
1270 double area,vx1,vy1,vx2,vy2,cp;
1271 x=0; y=0; area=signed_area(g);
1272 vx1=vertl[2*gen_to_vert[g][0]];
1273 vy1=vertl[2*gen_to_vert[g][0]+1];
1274 for(int i=0;i<gen_to_vert[g].size();i++){
1275 if(i==gen_to_vert[g].size()-1){
1276 vx2=vertl[2*gen_to_vert[g][0]];
1277 vy2=vertl[2*gen_to_vert[g][0]+1];
1278 }else{
1279 vx2=vertl[2*gen_to_vert[g][i+1]];
1280 vy2=vertl[2*gen_to_vert[g][i+1]+1];
1281 }
1282 cp=cross_product(vx1,vy1,vx2,vy2);
1283 x+=((vx1+vx2)*cp);
1284 y+=((vy1+vy2)*cp);
1285 vx1=vx2;
1286 vy1=vy2;
1287 }
1288 x=((1/(6*area))*x);
1289 y=((1/(6*area))*y);
1290 }
1291 //outputs lloyds allgorithms in files lloyds# until max distance from generator to centroid < epsilon
1292 void v_connect::lloyds(double epsilon){
1293 double gx,gy,cx,cy,max,current;
1294 int j=0;
1295 char *outfn1(new char[1000]);
1296 sprintf(outfn1,"lloyds");
1297 FILE *fp=safe_fopen(outfn1,"w");
1298 do{
1299 cout << "iteration" << j << endl;
1300 max=0;
1301 char *outfn2(new char[1000]);
1302 sprintf(outfn2,"lloyds%i",j);
1303 cout << "blah" << endl;
1304 draw_gnu(outfn2);
1305 cout << "yeah" << endl;
1306 delete[] outfn2;
1307 if(bd!=-1){
1308 for(int i=bd;i<ng;i++){
1309 cout << "i=" << i << endl;
1310 gx=vpos[2*mp[i]]; gy=vpos[2*mp[i]+1];
1311 centroid(i,cx,cy);
1312 current=pow(cx-gx,2)+pow(cy-gy,2);
1313 if(current>max) max=current;
1314 vpos[2*mp[i]]=cx;
1315 vpos[2*mp[i]+1]=cy;
1316 }
1317 }
1318 fprintf(fp,"iteration %i, max=%f, epsilon=%f",j,max,epsilon);
1319 nv=0;
1320 ne=0;
1321 cout << "1" << endl;
1322 current_vertices=init_vertices;
1323 current_edges=init_vertices;
1324 delete [] gen_to_gen_e;
1325 gen_to_gen_e=new vector<int>[ng];
1326 cout << "2" << endl;
1327 delete [] gen_to_gen_v;
1328 gen_to_gen_v=new vector<int>[ng];
1329 cout << "3" << endl;
1330 delete [] gen_to_ed;
1331 gen_to_ed=new vector<int>[ng];
1332 cout << "4" << endl;
1333 delete [] gen_to_vert;
1334 gen_to_vert=new vector<int>[ng];
1335 cout << "5" << endl;
1336 delete [] vertl;
1337 vertl=new double[2*init_vertices];
1338 cout << "6" << endl;
1339 delete [] vert_to_gen;
1340 vert_to_gen=new vector<int>[init_vertices];
1341 cout << "7" << endl;
1342 delete [] vertex_is_generator;
1343 cout << "7.1" << endl;
1344 vertex_is_generator=new int[init_vertices];
1345 cout << "7.2" << endl;
1346 for(int k=0;k<init_vertices;k++) vertex_is_generator[k]=-1;
1347 cout << "8" << endl;
1348 delete [] generator_is_vertex;
1349 generator_is_vertex=new int[ng];
1350 for(int k=0;k<ng;k++) generator_is_vertex[k]=-1;
1351 cout << "9" << endl;
1352 delete [] ed_to_vert;
1353 delete [] vert_to_ed;
1354 delete [] ed_to_gen;
1355 delete [] ed_on_bd;
1356 delete [] vert_on_bd;
1357 cout << "1..." << endl;
1358 assemble_vertex();
1359 cout << "2..." << endl;
1360 assemble_gen_ed();
1361 cout << "3..." << endl;
1362 assemble_boundary();
1363 cout << "4..." << endl;
1364 j++;
1365 }while(max>epsilon);
1366 fclose(fp);
1367 delete [] outfn1;
1368 }
1369
1370 void v_connect::draw_median_mesh(FILE *fp){
1371 double ex,ey,cx,cy;
1372 int e;
1373 for(int i=0;i<ng;i++){
1374 if(generator_is_vertex[i]!=-1) continue;
1375 centroid(i,cx,cy);
1376 for(int j=0;j<gen_to_ed[i].size();j++){
1377 e=gen_to_ed[i][j];
1378 if(e<0) e=~e;
1379 ex=(vertl[2*ed_to_vert[2*e]] + vertl[2*ed_to_vert[2*e+1]])/2;
1380 ey=(vertl[2*ed_to_vert[2*e]+1] + vertl[2*ed_to_vert[2*e+1]+1])/2;
1381 fprintf(fp,"\n %g %g\n %g %g \n\n\n",ex,ey,cx,cy);
1382 }
1383 }
1384
1385 }
1386
1387 void v_connect::draw_closest_generator(FILE *fp,double x,double y){
1388 int cg,pg,ng=0;
1389 double gx,gy,best,bestl,current;
1390 gx=vpos[2*mp[0]]; gy=vpos[2*mp[0]+1];
1391 best = pow(gx-x,2)+pow(gy-y,2);
1392 cout << "best=" << best << endl;
1393 fprintf(fp,"%g %g\n",gx,gy);
1394 do{
1395 cg=ng;
1396 for(int i=0;i<gen_to_gen_e[cg].size();i++){
1397 pg=gen_to_gen_e[cg][i];
1398 gx=vpos[2*mp[pg]]; gy=vpos[2*mp[pg]+1];
1399 current=pow(gx-x,2)+pow(gy-y,2);
1400 cout << "current=" << current << endl;
1401 if(current<best){
1402 cout << "changed" << endl;
1403 best=current;
1404 ng=pg;
1405 }
1406 }
1407 fprintf(fp,"%g %g\n",vpos[2*mp[ng]],vpos[2*mp[ng]+1]);
1408 }while(cg!=ng);
1409 }
1410
1411
1412
1413 }
Voro++ is a open source software library for the computation of the Voronoi diagram