21 current_vertices(init_vertices), current_vertex_order(init_vertex_order),
22 current_delete_size(init_delete_size), current_delete2_size(init_delete2_size),
23 ed(new int*[current_vertices]), nu(new int[current_vertices]),
24 pts(new double[3*current_vertices]), mem(new int[current_vertex_order]),
25 mec(new int[current_vertex_order]), mep(new int*[current_vertex_order]),
26 ds(new int[current_delete_size]), stacke(ds+current_delete_size),
27 ds2(new int[current_delete2_size]), stacke2(ds2+current_delete_size),
28 current_marginal(init_marginal), marg(new int[current_marginal]) {
31 mem[i]=init_n_vertices;
mec[i]=0;
32 mep[i]=
new int[init_n_vertices*((i<<1)+1)];
34 mem[3]=init_3_vertices;
mec[3]=0;
35 mep[3]=
new int[init_3_vertices*7];
37 mem[i]=init_n_vertices;
mec[i]=0;
38 mep[i]=
new int[init_n_vertices*((i<<1)+1)];
46 delete [] ds2;
delete [] ds;
47 delete []
mep;
delete []
mec;
48 delete []
mem;
delete []
pts;
49 delete []
nu;
delete []
ed;
55 template<
class vc_
class>
59 while(current_vertices<vb->
p) add_memory_vertices(vc);
70 for(j=0;j<
mec[i]*(2*i+1);j++)
mep[i][j]=vb->
mep[i][j];
71 for(j=0;j<
mec[i]*(2*i+1);j+=2*i+1)
ed[
mep[i][j+2*i]]=
mep[i]+j;
73 for(i=0;i<
p;i++)
nu[i]=vb->
nu[i];
74 for(i=0;i<3*p;i++)
pts[i]=vb->
pts[i];
85 for(j=0;j<c.
mec[i]*i;j++)
mne[i][j]=0;
86 for(j=0;j<c.
mec[i];j++)
ne[c.
mep[i][(2*i+1)*j+2*i]]=
mne[i]+(j*i);
98 for(j=0;j<c.
mec[i]*i;j++)
mne[i][j]=c.
mne[i][j];
99 for(j=0;j<c.
mec[i];j++)
ne[c.
mep[i][(2*i+1)*j+2*i]]=
mne[i]+(j*i);
108 while(ptsp<
pts+3*
p) {
109 *(ptsp++)=x;*(ptsp++)=y;*(ptsp++)=z;
124 template<
class vc_
class>
125 void voronoicell_base::add_memory(vc_class &vc,
int i,
int *stackp2) {
128 vc.n_allocate(i,init_n_vertices);
129 mep[i]=
new int[init_n_vertices*s];
130 mem[i]=init_n_vertices;
131 #if VOROPP_VERBOSE >=2
132 fprintf(stderr,
"Order %d vertex memory created\n",i);
137 if(
mem[i]>max_n_vertices) voro_fatal_error(
"Point memory allocation exceeded absolute maximum",
VOROPP_MEMORY_ERROR);
138 #if VOROPP_VERBOSE >=2
139 fprintf(stderr,
"Order %d vertex memory scaled up to %d\n",i,
mem[i]);
143 vc.n_allocate_aux1(i);
148 vc.n_set_to_aux1_offset(k,m);
151 for(dsp=ds2;dsp<stackp2;dsp++) {
152 if(
ed[*dsp]==
mep[i]+j) {
154 vc.n_set_to_aux1_offset(*dsp,m);
159 #if VOROPP_VERBOSE >=3
160 fputs(
"Relocated dangling pointer",stderr);
163 for(k=0;k<s;k++,j++) l[j]=
mep[i][j];
164 for(k=0;k<i;k++,m++) vc.n_copy_to_aux1(i,m);
168 vc.n_switch_to_aux1(i);
177 template<
class vc_
class>
178 void voronoicell_base::add_memory_vertices(vc_class &vc) {
180 if(i>max_vertices) voro_fatal_error(
"Vertex memory allocation exceeded absolute maximum",
VOROPP_MEMORY_ERROR);
181 #if VOROPP_VERBOSE >=2
182 fprintf(stderr,
"Vertex memory scaled up to %d\n",i);
188 vc.n_add_memory_vertices(i);
192 ppts=
new double[3*i];
203 template<
class vc_
class>
204 void voronoicell_base::add_memory_vorder(vc_class &vc) {
206 if(i>max_vertex_order) voro_fatal_error(
"Vertex order memory allocation exceeded absolute maximum",
VOROPP_MEMORY_ERROR);
207 #if VOROPP_VERBOSE >=2
208 fprintf(stderr,
"Vertex order memory scaled up to %d\n",i);
212 delete []
mem;mem=p1;
219 vc.n_add_memory_vorder(i);
220 current_vertex_order=i;
226 void voronoicell_base::add_memory_ds(
int *&stackp) {
229 #if VOROPP_VERBOSE >=2
233 while(dsp<stackp) *(dsnp++)=*(dsp++);
234 delete [] ds;ds=dsn;stackp=dsnp;
241 void voronoicell_base::add_memory_ds2(
int *&stackp2) {
244 #if VOROPP_VERBOSE >=2
248 while(dsp<stackp2) *(dsnp++)=*(dsp++);
249 delete [] ds2;ds2=dsn;stackp2=dsnp;
259 mec[3]=
p=8;xmin*=2;xmax*=2;ymin*=2;ymax*=2;zmin*=2;zmax*=2;
269 *q=1;q[1]=4;q[2]=2;q[3]=2;q[4]=1;q[5]=0;q[6]=0;
270 q[7]=3;q[8]=5;q[9]=0;q[10]=2;q[11]=1;q[12]=0;q[13]=1;
271 q[14]=0;q[15]=6;q[16]=3;q[17]=2;q[18]=1;q[19]=0;q[20]=2;
272 q[21]=2;q[22]=7;q[23]=1;q[24]=2;q[25]=1;q[26]=0;q[27]=3;
273 q[28]=6;q[29]=0;q[30]=5;q[31]=2;q[32]=1;q[33]=0;q[34]=4;
274 q[35]=4;q[36]=1;q[37]=7;q[38]=2;q[39]=1;q[40]=0;q[41]=5;
275 q[42]=7;q[43]=2;q[44]=4;q[45]=2;q[46]=1;q[47]=0;q[48]=6;
276 q[49]=5;q[50]=3;q[51]=6;q[52]=2;q[53]=1;q[54]=0;q[55]=7;
277 *
ed=q;
ed[1]=q+7;
ed[2]=q+14;
ed[3]=q+21;
278 ed[4]=q+28;
ed[5]=q+35;
ed[6]=q+42;
ed[7]=q+49;
296 *q=2;q[1]=5;q[2]=3;q[3]=4;q[4]=0;q[5]=0;q[6]=0;q[7]=0;q[8]=0;
297 q[9]=2;q[10]=4;q[11]=3;q[12]=5;q[13]=2;q[14]=2;q[15]=2;q[16]=2;q[17]=1;
298 q[18]=0;q[19]=4;q[20]=1;q[21]=5;q[22]=0;q[23]=3;q[24]=0;q[25]=1;q[26]=2;
299 q[27]=0;q[28]=5;q[29]=1;q[30]=4;q[31]=2;q[32]=3;q[33]=2;q[34]=1;q[35]=3;
300 q[36]=0;q[37]=3;q[38]=1;q[39]=2;q[40]=3;q[41]=3;q[42]=1;q[43]=1;q[44]=4;
301 q[45]=0;q[46]=2;q[47]=1;q[48]=3;q[49]=1;q[50]=3;q[51]=3;q[52]=1;q[53]=5;
302 *
ed=q;
ed[1]=q+9;
ed[2]=q+18;
ed[3]=q+27;
ed[4]=q+36;
ed[5]=q+45;
312 void voronoicell_base::init_tetrahedron_base(
double x0,
double y0,
double z0,
double x1,
double y1,
double z1,
double x2,
double y2,
double z2,
double x3,
double y3,
double z3) {
320 *q=1;q[1]=3;q[2]=2;q[3]=0;q[4]=0;q[5]=0;q[6]=0;
321 q[7]=0;q[8]=2;q[9]=3;q[10]=0;q[11]=2;q[12]=1;q[13]=1;
322 q[14]=0;q[15]=3;q[16]=1;q[17]=2;q[18]=2;q[19]=1;q[20]=2;
323 q[21]=0;q[22]=1;q[23]=2;q[24]=1;q[25]=2;q[26]=1;q[27]=3;
324 *
ed=q;
ed[1]=q+7;
ed[2]=q+14;
ed[3]=q+21;
333 for(i=0;i<
p;i++)
for(j=0;j<
nu[i];j++)
if(
ed[
ed[i][j]][
ed[i][nu[i]+j]]!=i)
334 printf(
"Relational error at point %d, edge %d.\n",i,j);
344 for(i=0;i<
p;i++)
for(j=1;j<
nu[i];j++)
for(k=0;k<j;k++)
if(
ed[i][j]==
ed[i][k])
345 printf(
"Duplicate edges: (%d,%d) and (%d,%d) [%d]\n",i,j,i,k,
ed[i][j]);
351 for(i=0;i<
p;i++)
for(j=0;j<
nu[i];j++) {
367 template<
class vc_
class>
368 inline bool voronoicell_base::search_for_outside_edge(vc_class &vc,
int &up) {
369 int i,lp,lw,*j(ds2),*stackp2(ds2);
374 for(i=0;i<
nu[
up];i++) {
377 if(lw==-1)
return true;
378 else if(lw==0) add_to_stack(vc,lp,stackp2);
388 template<
class vc_
class>
389 inline void voronoicell_base::add_to_stack(vc_class &vc,
int lp,
int *&stackp2) {
390 for(
int *k(ds2);k<stackp2;k++)
if(*k==lp)
return;
391 if(stackp2==stacke2) add_memory_ds2(stackp2);
403 template<
class vc_
class>
405 int count=0,i,j,k,lp=
up,cp,qp,rp,*stackp(ds),*stackp2(ds2),*dsp;
406 int us=0,ls=0,qs,iqs,cs,uw,qw,lw;
408 double u,l,r,q;
bool complicated_setup=
false,new_double_edge=
false,double_edge=
false;
411 n_marg=0;px=x;py=y;pz=z;prsq=rsq;
438 if(++count>=
p)
throw true;
440 for(us=0;us<ls;us++) {
465 complicated_setup=
true;
466 }
else complicated_setup=
false;
475 if(us==nu[up])
return true;
479 if(++count>=
p)
throw true;
481 for(us=0;us<qs;us++) {
494 if(us==nu[up])
return true;
499 up=qp;us=
ed[lp][nu[lp]+ls];u=q;
500 complicated_setup=
false;
503 complicated_setup=
true;
510 complicated_setup=
true;
518 #if VOROPP_VERBOSE >=1
519 fputs(
"Bailed out of convex calculation\n",stderr);
522 for(qp=0;qp<
p;qp++) {
528 for(us=0;us<nu[qp];us++) {
538 complicated_setup=
true;
540 complicated_setup=
false;
550 for(ls=0;ls<nu[qp];ls++) {
560 complicated_setup=
true;
562 complicated_setup=
false;
564 us=
ed[lp][nu[lp]+ls];
573 complicated_setup=
true;
577 if(qp==p)
return qw==-1?
true:
false;
583 if(
p==current_vertices) add_memory_vertices(vc);
585 if(complicated_setup) {
592 if(!search_for_outside_edge(vc,up))
return false;
620 if(i==nu[up])
return false;
643 if(j==nu[up]&&i==1&&rp==0) {
651 while (nu[
p]>=current_vertex_order) add_memory_vorder(vc);
652 if(
mec[nu[
p]]==mem[nu[
p]]) add_memory(vc,nu[p],stackp2);
653 vc.n_set_pointer(p,nu[p]);
687 if(i==0)
return true;
719 while(nu[
p]>=current_vertex_order) add_memory_vorder(vc);
720 if(
mec[nu[
p]]==mem[nu[
p]]) add_memory(vc,nu[p],stackp2);
725 vc.n_set_pointer(p,nu[p]);
755 vc.n_copy(
p,k,up,qs);
757 }
else vc.n_copy(
p,0,up,qs);
760 if(stackp2==stacke2) add_memory_ds2(stackp2);
772 ed[qp][nu[qp]<<1]=-
p;
780 if(stackp==stacke) add_memory_ds(stackp);
790 if(
mec[3]==mem[3]) add_memory(vc,3,stackp2);
791 vc.n_set_pointer(
p,3);
793 vc.n_copy(
p,1,up,us);
794 vc.n_copy(
p,2,lp,ls);
812 while(qp!=up||qs!=us) {
827 if(stackp==stacke) add_memory_ds(stackp);
837 if(
p==current_vertices) add_memory_vertices(vc);
843 if(
mec[3]==mem[3]) add_memory(vc,3,stackp2);
844 ls=
ed[qp][qs+nu[qp]];
845 vc.n_set_pointer(
p,3);
847 vc.n_copy(
p,1,qp,qs);
848 vc.n_copy(
p,2,lp,ls);
869 if(
p==current_vertices) add_memory_vertices(vc);
874 qs=
ed[qp][nu[qp]+qs];
893 j=-
ed[qp][nu[qp]<<1];
899 new_double_edge=
false;
918 i=-
ed[lp][nu[lp]<<1];
923 if(
ed[i][nu[i]-1]==j) {
924 new_double_edge=
true;
926 }
else new_double_edge=
false;
935 if(j==rp&&lp==up&&
ed[qp][nu[qp]+qs]==us) {
936 new_double_edge=
true;
938 }
else new_double_edge=
false;
940 }
else new_double_edge=
false;
953 i=-
ed[lp][nu[lp]<<1];
955 new_double_edge=
true;
957 }
else new_double_edge=
false;
958 }
else new_double_edge=
false;
965 while(k>=current_vertex_order) add_memory_vorder(vc);
966 if(
mec[k]==mem[k]) add_memory(vc,k,stackp2);
979 edp=
mep[k]+((k<<1)+1)*
mec[k]++;
984 edp[k+i]=
ed[j][nu[j]+i];
992 edd=
mep[nu[j]]+((nu[j]<<1)+1)*--
mec[nu[j]];
994 for(lw=0;lw<=(nu[j]<<1);lw++)
ed[j][lw]=edd[lw];
995 vc.n_set_aux2_copy(j,nu[j]);
996 vc.n_copy_pointer(edd[nu[j]<<1],j);
997 ed[edd[nu[j]<<1]]=
ed[j];
1005 vc.n_set_pointer(
p,k);
1008 if(stackp2==stacke2) add_memory_ds2(stackp2);
1013 ed[qp][nu[qp]<<1]=-
p;
1027 ed[cp][nu[cp]+cs]=i;
1034 while(i<(new_double_edge?k:k-1)) {
1036 lp=
ed[qp][qs];ls=
ed[qp][nu[qp]+qs];
1037 vc.n_copy(j,i,qp,qs);
1041 ed[lp][nu[lp]+ls]=i;
1048 vc.n_copy(j,new_double_edge?0:cs,qp,qs);
1052 double_edge=new_double_edge;
1059 ed[cp][nu[cp]+cs]=0;
1066 if(
ed[j][nu[j]]!=-1) {
1069 }
else *dsp=*(--stackp);
1074 for(dsp=ds2;dsp<stackp2;dsp++) {
1077 if(
ed[j][nu[j]]!=-1) {
1079 if(stackp==stacke) add_memory_ds(stackp);
1085 for(dsp=ds;dsp<stackp;dsp++) {
1087 for(edp=
ed[cp];edp<
ed[cp]+nu[cp];edp++) {
1089 if(qp!=-1&&
ed[qp][nu[qp]]!=-1) {
1090 if(stackp==stacke) {
1092 add_memory_ds(stackp);
1105 while(
ed[
p][nu[
p]]==-1) {
1107 edp=
ed[
p];edd=(
mep[j]+((j<<1)+1)*--
mec[j]);
1108 while(edp<
ed[
p]+(j<<1)+1) *(edp++)=*(edd++);
1109 vc.n_set_aux2_copy(
p,j);
1110 vc.n_copy_pointer(
ed[
p][(j<<1)],
p);
1125 while(edp<
ed[up]+(j<<1)+1) *(edp++)=*(edd++);
1126 vc.n_set_aux2_copy(up,j);
1127 vc.n_copy_pointer(
ed[up][j<<1],up);
1128 vc.n_copy_pointer(up,
p);
1134 for(i=0;i<nu[
up];i++) ed[ed[up][i]][ed[up][nu[up]+i]]=up;
1143 return collapse_order2(vc);
1156 template<
class vc_
class>
1157 inline bool voronoicell_base::collapse_order2(vc_class &vc) {
1158 if(!collapse_order1(vc))
return false;
1164 j=
mep[2][5*i];k=
mep[2][5*i+1];
1166 #if VOROPP_VERBOSE >=1
1167 fputs(
"Order two vertex joins itself",stderr);
1173 for(l=0;l<nu[j];l++) {
1174 if(
ed[j][l]==k)
break;
1180 a=
mep[2][5*i+2];b=
mep[2][5*i+3];i=
mep[2][5*i+4];
1187 if(!delete_connection(vc,j,a,
false))
return false;
1188 if(!delete_connection(vc,k,b,
true))
return false;
1199 for(k=0;k<nu[
p];k++)
ed[
ed[
p][k]][
ed[
p][nu[
p]+k]]=i;
1200 vc.n_copy_pointer(i,
p);
1207 if(!collapse_order1(vc))
return false;
1217 template<
class vc_
class>
1218 inline bool voronoicell_base::collapse_order1(vc_class &vc) {
1222 #if VOROPP_VERBOSE >=1
1223 fputs(
"Order one collapse\n",stderr);
1226 j=
mep[1][3*i];k=
mep[1][3*i+1];
1228 if(!delete_connection(vc,j,k,
false))
return false;
1236 for(k=0;k<nu[
p];k++)
ed[
ed[
p][k]][
ed[
p][nu[
p]+k]]=i;
1237 vc.n_copy_pointer(i,
p);
1252 template<
class vc_
class>
1253 inline bool voronoicell_base::delete_connection(vc_class &vc,
int j,
int k,
bool hand) {
1255 int i=nu[j]-1,l,*edp,*edd,m;
1256 #if VOROPP_VERBOSE >=1
1258 fputs(
"Zero order vertex formed\n",stderr);
1262 if(
mec[i]==mem[i]) add_memory(vc,i,ds2);
1264 for(l=0;l<q;l++) vc.n_copy_aux1(j,l);
1266 vc.n_copy_aux1_shift(j,l);
1269 edp=
mep[i]+((i<<1)+1)*
mec[i]++;
1273 edp[l+i]=
ed[j][l+nu[j]];
1284 edd=
mep[nu[j]]+((nu[j]<<1)+1)*--
mec[nu[j]];
1285 for(l=0;l<=(nu[j]<<1);l++)
ed[j][l]=edd[l];
1286 vc.n_set_aux2_copy(j,nu[j]);
1287 vc.n_set_to_aux2(edd[nu[j]<<1]);
1288 vc.n_set_to_aux1(j);
1289 ed[edd[nu[j]<<1]]=edd;
1300 const double fe=1/48.0;
1303 double ux,uy,uz,vx,vy,vz,wx,wy,wz;
1306 uy=pts[1]-pts[3*i+1];
1307 uz=pts[2]-pts[3*i+2];
1308 for(j=0;j<nu[i];j++) {
1314 vy=pts[3*k+1]-pts[1];
1315 vz=pts[3*k+2]-pts[2];
1316 m=
ed[k][l];
ed[k][l]=-1-m;
1320 wy=pts[3*m+1]-pts[1];
1321 wz=pts[3*m+2]-pts[2];
1322 vol+=ux*vy*wz+uy*vz*wx+uz*vx*wy-uz*vy*wx-uy*vx*wz-ux*vz*wy;
1323 k=m;l=n;vx=wx;vy=wy;vz=wz;
1324 m=
ed[k][l];
ed[k][l]=-1-m;
1340 double ux,uy,uz,vx,vy,vz,wx,wy,wz;
1341 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1347 m=
ed[k][l];
ed[k][l]=-1-m;
1351 uy=pts[3*k+1]-pts[3*i+1];
1352 uz=pts[3*k+2]-pts[3*i+2];
1353 vx=pts[3*m]-pts[3*i];
1354 vy=pts[3*m+1]-pts[3*i+1];
1355 vz=pts[3*m+2]-pts[3*i+2];
1359 area+=sqrt(wx*wx+wy*wy+wz*wz);
1361 m=
ed[k][l];
ed[k][l]=-1-m;
1363 v.push_back(0.125*area);
1375 double ux,uy,uz,vx,vy,vz,wx,wy,wz;
1376 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1381 m=
ed[k][l];
ed[k][l]=-1-m;
1385 uy=pts[3*k+1]-pts[3*i+1];
1386 uz=pts[3*k+2]-pts[3*i+2];
1387 vx=pts[3*m]-pts[3*i];
1388 vy=pts[3*m+1]-pts[3*i+1];
1389 vz=pts[3*m+2]-pts[3*i+2];
1393 area+=sqrt(wx*wx+wy*wy+wz*wz);
1395 m=
ed[k][l];
ed[k][l]=-1-m;
1409 double tvol,vol=0;cx=cy=cz=0;
1411 double ux,uy,uz,vx,vy,vz,wx,wy,wz;
1414 uy=pts[1]-pts[3*i+1];
1415 uz=pts[2]-pts[3*i+2];
1416 for(j=0;j<nu[i];j++) {
1422 vy=pts[3*k+1]-pts[1];
1423 vz=pts[3*k+2]-pts[2];
1424 m=
ed[k][l];
ed[k][l]=-1-m;
1428 wy=pts[3*m+1]-pts[1];
1429 wz=pts[3*m+2]-pts[2];
1430 tvol=ux*vy*wz+uy*vz*wx+uz*vx*wy-uz*vy*wx-uy*vx*wz-ux*vz*wy;
1432 cx+=(wx+vx-ux)*tvol;
1433 cy+=(wy+vy-uy)*tvol;
1434 cz+=(wz+vz-uz)*tvol;
1435 k=m;l=n;vx=wx;vy=wy;vz=wz;
1436 m=
ed[k][l];
ed[k][l]=-1-m;
1442 if(vol>tolerance_sq) {
1444 cx=cx*vol+0.5*(*pts);
1445 cy=cy*vol+0.5*
pts[1];
1446 cz=cz*vol+0.5*
pts[2];
1455 double r,s,*ptsp=
pts+3,*ptse=
pts+3*
p;
1458 s=*ptsp*(*ptsp);ptsp++;
1459 s+=*ptsp*(*ptsp);ptsp++;
1460 s+=*ptsp*(*ptsp);ptsp++;
1470 double dis=0,dx,dy,dz;
1471 for(i=0;i<
p-1;i++)
for(j=0;j<nu[i];j++) {
1475 dy=pts[3*k+1]-pts[3*i+1];
1476 dz=pts[3*k+2]-pts[3*i+2];
1477 dis+=sqrt(dx*dx+dy*dy+dz*dz);
1488 int i,j,k;
double *ptsp=
pts,*pt2;
1489 char posbuf1[128],posbuf2[128];
1490 for(i=0;i<
p;i++,ptsp+=3) {
1491 sprintf(posbuf1,
"%g,%g,%g",x+*ptsp*0.5,y+ptsp[1]*0.5,z+ptsp[2]*0.5);
1492 fprintf(fp,
"sphere{<%s>,r}\n",posbuf1);
1493 for(j=0;j<nu[i];j++) {
1497 sprintf(posbuf2,
"%g,%g,%g",x+*pt2*0.5,y+0.5*pt2[1],z+0.5*pt2[2]);
1498 if(strcmp(posbuf1,posbuf2)!=0) fprintf(fp,
"cylinder{<%s>,<%s>,r}\n",posbuf1,posbuf2);
1509 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1512 fprintf(fp,
"%g %g %g\n",x+0.5*
pts[3*i],y+0.5*
pts[3*i+1],z+0.5*
pts[3*i+2]);
1515 ed[k][
ed[l][nu[l]+m]]=-1-l;
1518 fprintf(fp,
"%g %g %g\n",x+0.5*
pts[3*k],y+0.5*
pts[3*k+1],z+0.5*
pts[3*k+2]);
1519 }
while (search_edge(l,m,k));
1526 inline bool voronoicell_base::search_edge(
int l,
int &m,
int &k) {
1527 for(m=0;m<nu[l];m++) {
1529 if(k>=0)
return true;
1545 fprintf(fp,
"mesh2 {\nvertex_vectors {\n%d\n",
p);
1546 for(i=0;i<
p;i++,ptsp+=3) fprintf(fp,
",<%g,%g,%g>\n",x+*ptsp*0.5,y+ptsp[1]*0.5,z+ptsp[2]*0.5);
1547 fprintf(fp,
"}\nface_indices {\n%d\n",(p-2)<<1);
1548 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1553 m=
ed[k][l];
ed[k][l]=-1-m;
1556 fprintf(fp,
",<%d,%d,%d>\n",i,k,m);
1558 m=
ed[k][l];
ed[k][l]=-1-m;
1562 fputs(
"}\ninside_vector <0,0,1>\n}\n",fp);
1574 for(i=0;i<
p;i++)
for(j=0;j<nu[i];j++) {
1575 if(
ed[i][j]>=0) voro_fatal_error(
"Edge reset routine found a previously untested edge",
VOROPP_INTERNAL_ERROR);
1576 ed[i][j]=-1-
ed[i][j];
1590 inline int voronoicell_base::m_test(
int n,
double &ans) {
1591 double *pp=
pts+n+(n<<1);
1595 if(ans<-tolerance2) {
1597 }
else if(ans>tolerance2) {
1600 return check_marginal(n,ans);
1615 int voronoicell_base::check_marginal(
int n,
double &ans) {
1617 for(i=0;i<n_marg;i+=2)
if(marg[i]==n)
return marg[i+1];
1618 if(n_marg==current_marginal) {
1619 current_marginal<<=1;
1620 if(current_marginal>max_marginal)
1621 voro_fatal_error(
"Marginal case buffer allocation exceeded absolute maximum",
VOROPP_MEMORY_ERROR);
1622 #if VOROPP_VERBOSE >=2
1623 fprintf(stderr,
"Marginal cases buffer scaled up to %d\n",i);
1625 int *pmarg=
new int[current_marginal];
1626 for(
int j=0;j<n_marg;j++) pmarg[j]=marg[j];
1631 marg[n_marg++]=ans>tolerance?1:(ans<-tolerance?-1:0);
1632 return marg[n_marg-1];
1642 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1644 if(k>=0) normals_search(v,i,j,k);
1660 inline void voronoicell_base::normals_search(std::vector<double> &v,
int i,
int j,
int k) {
1663 double ux,uy,uz,vx,vy,vz,wx,wy,wz,wmag;
1665 m=
ed[k][l];
ed[k][l]=-1-m;
1667 uy=pts[3*m+1]-pts[3*k+1];
1668 uz=pts[3*m+2]-pts[3*k+2];
1671 if(ux*ux+uy*uy+uz*uz>tolerance_sq) {
1674 k=m;m=
ed[k][l];
ed[k][l]=-1-m;
1675 vx=pts[3*m]-pts[3*k];
1676 vy=pts[3*m+1]-pts[3*k+1];
1677 vz=pts[3*m+2]-pts[3*k+2];
1684 wmag=wx*wx+wy*wy+wz*wz;
1688 if(wmag>tolerance_sq) {
1692 v.push_back(wx*wmag);
1693 v.push_back(wy*wmag);
1694 v.push_back(wz*wmag);
1700 k=m;m=
ed[k][l];
ed[k][l]=-1-m;
1723 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1746 for(
int i=0;i<
p;i++) v[i]=nu[i];
1753 fprintf(fp,
"%d",*nu);
1754 for(
int *nup=nu+1;nup<nu+
p;nup++) fprintf(fp,
" %d",*nup);
1763 for(
int i=0;i<3*
p;i+=3) {
1765 v[i+1]=*(ptsp++)*0.5;
1766 v[i+2]=*(ptsp++)*0.5;
1774 fprintf(fp,
"(%g,%g,%g)",*pts*0.5,pts[1]*0.5,pts[2]*0.5);
1775 for(
double *ptsp=pts+3;ptsp<pts+3*
p;ptsp+=3) fprintf(fp,
" (%g,%g,%g)",*ptsp*0.5,ptsp[1]*0.5,ptsp[2]*0.5);
1787 for(
int i=0;i<3*
p;i+=3) {
1788 v[i]=x+*(ptsp++)*0.5;
1789 v[i+1]=y+*(ptsp++)*0.5;
1790 v[i+2]=z+*(ptsp++)*0.5;
1800 fprintf(fp,
"(%g,%g,%g)",x+*pts*0.5,y+pts[1]*0.5,z+pts[2]*0.5);
1801 for(
double *ptsp=pts+3;ptsp<pts+3*
p;ptsp+=3) fprintf(fp,
" (%g,%g,%g)",x+*ptsp*0.5,y+ptsp[1]*0.5,z+ptsp[2]*0.5);
1810 double dx,dy,dz,perim;
1811 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1814 dx=pts[3*k]-pts[3*i];
1815 dy=pts[3*k+1]-pts[3*i+1];
1816 dz=pts[3*k+2]-pts[3*i+2];
1817 perim=sqrt(dx*dx+dy*dy+dz*dz);
1822 dx=pts[3*m]-pts[3*k];
1823 dy=pts[3*m+1]-pts[3*k+1];
1824 dz=pts[3*m+2]-pts[3*k+2];
1825 perim+=sqrt(dx*dx+dy*dy+dz*dz);
1830 v.push_back(0.5*perim);
1840 int i,j,k,l,m,vp(0),vn;
1842 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1869 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1894 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
1907 if((
unsigned int) q>=v.size()) v.resize(q+1,0);
1921 double g=x*pts[3*
up]+y*pts[3*up+1]+z*pts[3*up+2];
1922 if(g<rsq)
return plane_intersects_track(x,y,z,rsq,g);
1935 double g=x*pts[3*
up]+y*pts[3*up+1]+z*pts[3*up+2];
1937 int ca=1,cc=
p>>3,mp=1;
1940 m=x*pts[3*mp]+y*pts[3*mp+1]+z*pts[3*mp+2];
1942 if(m>rsq)
return true;
1947 return plane_intersects_track(x,y,z,rsq,g);
1959 inline bool voronoicell_base::plane_intersects_track(
double x,
double y,
double z,
double rsq,
double g) {
1960 int count=0,ls,us,tp;
1964 for(us=0;us<nu[
up];us++) {
1966 t=x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2];
1972 #if VOROPP_VERBOSE >=1
1973 fputs(
"Bailed out of convex calculation",stderr);
1975 for(tp=0;tp<p;tp++) if(x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2]>rsq)
return true;
1982 for(us=0;us<ls;us++) {
1984 g=x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2];
1991 g=x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2];
1995 if(us==nu[up])
return false;
1997 ls=
ed[
up][nu[
up]+us];up=tp;t=g;
2008 int edges=0,*nup=
nu;
2009 while(nup<nu+
p) edges+=*(nup++);
2024 char *fmp=(
const_cast<char*
>(format));
2025 std::vector<int> vi;
2026 std::vector<double> vd;
2033 case 'i': fprintf(fp,
"%d",i);
break;
2034 case 'x': fprintf(fp,
"%g",x);
break;
2035 case 'y': fprintf(fp,
"%g",y);
break;
2036 case 'z': fprintf(fp,
"%g",z);
break;
2037 case 'q': fprintf(fp,
"%g %g %g",x,y,z);
break;
2038 case 'r': fprintf(fp,
"%g",r);
break;
2041 case 'w': fprintf(fp,
"%d",
p);
break;
2057 voro_print_vector(vi,fp);
2059 case 'a':
face_orders(vi);voro_print_vector(vi,fp);
break;
2060 case 'f':
face_areas(vd);voro_print_vector(vd,fp);
break;
2063 voro_print_face_vertices(vi,fp);
2066 voro_print_positions(vd,fp);
2069 voro_print_vector(vi,fp);
2073 case 'v': fprintf(fp,
"%g",
volume());
break;
2077 fprintf(fp,
"%g %g %g",cx,cy,cz);
2082 fprintf(fp,
"%g %g %g",x+cx,y+cy,z+cz);
2086 case 0: fmp--;
break;
2090 default: putc(
'%',fp);putc(*fmp,fp);
2092 }
else putc(*fmp,fp);
2106 init_base(xmin,xmax,ymin,ymax,zmin,zmax);
2108 *q=-5;q[1]=-3;q[2]=-1;
2109 q[3]=-5;q[4]=-2;q[5]=-3;
2110 q[6]=-5;q[7]=-1;q[8]=-4;
2111 q[9]=-5;q[10]=-4;q[11]=-2;
2112 q[12]=-6;q[13]=-1;q[14]=-3;
2113 q[15]=-6;q[16]=-3;q[17]=-2;
2114 q[18]=-6;q[19]=-4;q[20]=-1;
2115 q[21]=-6;q[22]=-2;q[23]=-4;
2116 *
ne=q;
ne[1]=q+3;
ne[2]=q+6;
ne[3]=q+9;
2117 ne[4]=q+12;
ne[5]=q+15;
ne[6]=q+18;
ne[7]=q+21;
2130 *q=-5;q[1]=-6;q[2]=-7;q[3]=-8;
2131 q[4]=-1;q[5]=-2;q[6]=-3;q[7]=-4;
2132 q[8]=-6;q[9]=-5;q[10]=-2;q[11]=-1;
2133 q[12]=-8;q[13]=-7;q[14]=-4;q[15]=-3;
2134 q[16]=-5;q[17]=-8;q[18]=-3;q[19]=-2;
2135 q[20]=-7;q[21]=-6;q[22]=-1;q[23]=-4;
2136 *
ne=q;
ne[1]=q+4;
ne[2]=q+8;
ne[3]=q+12;
ne[4]=q+16;
ne[5]=q+20;
2147 void voronoicell_neighbor::init_tetrahedron(
double x0,
double y0,
double z0,
double x1,
double y1,
double z1,
double x2,
double y2,
double z2,
double x3,
double y3,
double z3) {
2148 init_tetrahedron_base(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3);
2150 *q=-4;q[1]=-3;q[2]=-2;
2151 q[3]=-3;q[4]=-4;q[5]=-1;
2152 q[6]=-4;q[7]=-2;q[8]=-1;
2153 q[9]=-2;q[10]=-3;q[11]=-1;
2154 *
ne=q;
ne[1]=q+3;
ne[2]=q+6;
ne[3]=q+9;
2161 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
2170 if(
ne[k][l]!=q) fprintf(stderr,
"Facet error at (%d,%d)=%d, started from (%d,%d)=%d\n",k,l,
ne[k][l],i,j,q);
2184 for(i=0;i<3;i++)
mne[i]=
new int[init_n_vertices*i];
2185 mne[3]=
new int[init_3_vertices*3];
2192 for(
int i=current_vertex_order-1;i>=0;i--)
if(mem[i]>0)
delete []
mne[i];
2201 for(i=1;i<
p;i++)
for(j=0;j<nu[i];j++) {
2204 v.push_back(
ne[i][j]);
2223 for(
int i=0;i<
p;i++,ptsp+=3) {
2224 printf(
"%d %d ",i,nu[i]);
2225 for(j=0;j<nu[i];j++) printf(
" %d",
ed[i][j]);
2227 while(j<(nu[i]<<1)) printf(
" %d",
ed[i][j]);
2228 printf(
" %d",
ed[i][j]);
2230 printf(
" %g %g %g %p",*ptsp,ptsp[1],ptsp[2],(
void*)
ed[i]);
2231 if(
ed[i]>=
mep[nu[i]]+
mec[nu[i]]*((nu[i]<<1)+1)) puts(
" Memory error");
2241 while(j<nu[i]-1) printf(
"%d,",
ne[i][j++]);
2242 printf(
"%d)",
ne[i][j]);
2243 }
else printf(
" ()");
void init_tetrahedron(double x0, double y0, double z0, double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3)
void output_vertex_orders(FILE *fp=stdout)
Extension of the voronoicell_base class to represent a Voronoi cell with neighbor information...
int cycle_down(int a, int p)
void centroid(double &cx, double &cy, double &cz)
bool plane_intersects(double x, double y, double z, double rsq)
void face_areas(std::vector< double > &v)
#define VOROPP_MEMORY_ERROR
void init_tetrahedron_base(double x0, double y0, double z0, double x1, double y1, double z1, double x2, double y2, double z2, double x3, double y3, double z3)
void draw_gnuplot(double x, double y, double z, FILE *fp=stdout)
virtual void neighbors(std::vector< int > &v)
void construct_relations()
void init(double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
virtual void print_edges_neighbors(int i)
Header file for the voronoicell and related classes.
Master configuration file for setting various compile-time options.
bool nplane(vc_class &vc, double x, double y, double z, double rsq, int p_id)
void init_octahedron(double l)
virtual void neighbors(std::vector< int > &v)
void check_memory_for_copy(vc_class &vc, voronoicell_base *vb)
#define VOROPP_INTERNAL_ERROR
void vertices(std::vector< double > &v)
void face_perimeters(std::vector< double > &v)
virtual void print_edges_neighbors(int i)
void output_vertices(FILE *fp=stdout)
void operator=(voronoicell &c)
void draw_pov(double x, double y, double z, FILE *fp=stdout)
double max_radius_squared()
void face_orders(std::vector< int > &v)
int cycle_up(int a, int p)
void face_vertices(std::vector< int > &v)
void init_octahedron_base(double l)
bool plane_intersects_guess(double x, double y, double z, double rsq)
Header file for the small helper functions.
double total_edge_distance()
void copy(voronoicell_base *vb)
void face_freq_table(std::vector< int > &v)
void output_custom(const char *format, FILE *fp=stdout)
A class representing a single Voronoi cell.
Extension of the voronoicell_base class to represent a Voronoi cell without neighbor information...
void draw_pov_mesh(double x, double y, double z, FILE *fp=stdout)
void init_base(double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
void vertex_orders(std::vector< int > &v)
void translate(double x, double y, double z)
void normals(std::vector< double > &v)