1 #include "Global.h"
2 #include "Nef.h"
3
4 #ifndef CEQ_Nmax
5 #define CEQ_Nmax EQUA_Nmax
6 #endif
7
8 /* =============== Typedefs and Headers =================== */
9
10 int GLZ_Make_Trian_NF(Long X[][VERT_Nmax], int *n, int *nv,
11 GL_Long G[][POLY_Dmax]); /* current=best */
12
13 void Poly_Sym(PolyPointList *_P, VertexNumList *_V, EqList *_F, int *sym_num,
14 int V_perm[][VERT_Nmax]);
15
16 /* =============== local Typedefs and Headers =================== */
17
18 typedef struct {
19 int ne;
20 Equation e[CEQ_Nmax];
21 } CEqList;
22
23 typedef struct {
24 int nv; /* #vertices of face */
25 int v[VERT_Nmax]; /* vertices of face */
26 } VList;
27
28 typedef struct {
29 int Nv; /* #vertices */
30 int nf; /* #facets */
31 VList *vl; /* vertices of facets */
32 } FVList;
33
34 typedef struct {
35 int f;
36 int v;
37 } Step;
38
39 typedef struct {
40 int s[VERT_Nmax];
41 } V_Flag;
42
43 typedef struct {
44 int m[FACE_Nmax];
45 } M_Rank;
46
47 typedef struct {
48 int M[POLY_Dmax][POLY_Dmax];
49 int d;
50 int codim;
51 } MMatrix;
52
53 typedef struct {
54 int d;
55 int nv;
56 Long X[POLY_Dmax][VERT_Nmax];
57 } XMatrix;
58
59 typedef struct {
60 int d;
61 GL_Long G[POLY_Dmax][POLY_Dmax];
62 } GMatrix;
63
64 typedef struct {
65 int Vp[SYM_Nmax][VERT_Nmax];
66 int ns;
67 } SYM;
68
69 typedef struct {
70 int A[POLY_Dmax];
71 int m;
72 int M;
73 } Subset;
74
75 int GLZ_Start_Simplex(PolyPointList *_P, VertexNumList *_V, CEqList *_C);
76
77
78 /* =============== End of Typedefs and Headers =================== */
79
Print_FVl(FVList * _FVl,const char * comment)80 void Print_FVl(FVList *_FVl, const char *comment){
81 int i,j;
82
83 fprintf(outFILE,"%s\n",comment);
84 for (i = 0; i < _FVl->nf; i++){
85 for (j = 0; j < _FVl->vl[i].nv; j++)
86 fprintf(outFILE, " %d ", _FVl->vl[i].v[j]);
87 fprintf(outFILE, "nv: %d\n",_FVl->vl[i].nv);fflush(0);
88 }
89 }
90
Print_M(MMatrix * _M,int * _nf,const char * comment)91 void Print_M(MMatrix *_M, int *_nf, const char *comment){
92 int i, j, k;
93
94 fprintf(outFILE,"%s\n",comment);
95 for (i = 0; i < *_nf; i++){
96 fprintf(outFILE, "\n\n facet %d:\n", i);
97 for(j = 0; j < _M[i].codim; j++){
98 fprintf(outFILE, "m[%d]: ",j);
99 for(k = 0; k < _M[i].d; k++)
100 fprintf(outFILE, "%d ",_M[i].M[k][j]);
101 fprintf(outFILE,"\n");
102 }
103 }
104 }
105
Dir_Product(PartList * _PTL,VertexNumList * _V,PolyPointList * _P)106 void Dir_Product(PartList *_PTL, VertexNumList *_V, PolyPointList *_P){
107
108 int i, j, k, d, rank;
109
110 CEqList CEtemp;
111 VertexNumList Vtemp;
112
113 PolyPointList *_PV = (PolyPointList *) malloc(sizeof(PolyPointList));
114 if(_PV == NULL) Die("Unable to alloc space for PolyPointList _PV");
115
116 _PV->n = _P->n;
117 for(i = 0; i < _PTL->n; i++){
118 rank = 0;
119 for(j = 0; j < _PTL->codim; j++){
120 _PV->np = 0;
121 for(k = 0; k < _PTL->nv; k++)
122 if(_PTL->S[i][k] == j){
123 for(d = 0; d < _P->n; d++)
124 _PV->x[_PV->np][d] = _P->x[_V->v[k]][d];
125 _PV->np += 1;
126 }
127 for(d = 0; d < _P->n; d++)
128 _PV->x[_PV->np][d] = 0;
129 _PV->np += 1;
130 rank += GLZ_Start_Simplex(_PV, &Vtemp, &CEtemp);
131 }
132 if(rank == _P->n)
133 _PTL->DirProduct[i] = 1;
134 else
135 _PTL->DirProduct[i] = 0;
136 }
137 free(_PV);
138 }
139
Set_To_Vlist(int * _S,VertexNumList * _V,PolyPointList * _P,PolyPointList * _PV,Subset * _Pset)140 void Set_To_Vlist(int *_S, VertexNumList *_V, PolyPointList *_P, PolyPointList
141 *_PV, Subset *_Pset){
142 int i, d;
143
144 _PV->n = _P->n; _PV->np = 0;
145 for(i = 0; i < _V->nv; i++)
146 if( _Pset->A[_S[i]] == 1){
147 for(d = 0; d < _P->n; d++)
148 _PV->x[_PV->np][d] = _P->x[_V->v[i]][d];
149 _PV->np ++;
150 }
151 for(d = 0; d < _P->n; d++)
152 _PV->x[_PV->np][d] = 0;
153 _PV->np += 1;
154 }
155
New_Set(int i,Subset * _Pset,Subset * _CPset)156 int New_Set(int i, Subset *_Pset, Subset *_CPset){
157
158 if(_Pset->A[i] == 0){
159 _Pset->A[i] = 1; _Pset->m++;
160 _CPset->A[i] = 0; _CPset->m--;
161 return 1;
162 }
163 else
164 return 0;
165 }
166
Old_Set(int i,Subset * _Pset,Subset * _CPset)167 void Old_Set(int i, Subset *_Pset, Subset *_CPset){
168 _Pset->A[i] = 0; _Pset->m--;
169 _CPset->A[i] = 1; _CPset->m++;
170 }
171
Select_Set(int * _S,VertexNumList * _V,PolyPointList * _P,PolyPointList * _PV,CEqList * _CEtemp,VertexNumList * _Vtemp,int * _DirFlag,Subset * _Pset,Subset * _CPset,int nmax_old)172 void Select_Set(int *_S, VertexNumList *_V, PolyPointList *_P, PolyPointList *_PV,
173 CEqList *_CEtemp, VertexNumList *_Vtemp, int * _DirFlag, Subset
174 *_Pset, Subset *_CPset, int nmax_old){
175 int rank, nmax;
176
177 Set_To_Vlist(_S, _V, _P, _PV, _Pset);
178 rank = GLZ_Start_Simplex(_PV, _Vtemp, _CEtemp);
179 Set_To_Vlist(_S, _V, _P, _PV, _CPset);
180 rank += GLZ_Start_Simplex(_PV, _Vtemp, _CEtemp);
181
182 /*{int i; for(i=0;i<_Pset->M;i++) printf(" %d",_Pset->A[i]); printf("\n");}*/
183 if(rank == _P->n)
184 *_DirFlag = 1;
185 else{
186 if(_Pset->m < (_Pset->M/2 + (_Pset->M % 2))){
187 nmax = (nmax_old + 1);
188 while(!*_DirFlag && (nmax < _Pset->M)){
189 /*printf("aaa");*/
190 if (New_Set(nmax, _Pset, _CPset)){
191 Select_Set(_S, _V, _P, _PV, _CEtemp, _Vtemp, _DirFlag, _Pset, _CPset,
192 nmax);
193 Old_Set(nmax, _Pset, _CPset);
194 }
195 nmax++;
196 }
197 }
198 }
199 }
200
REC_Dir_Product(PartList * _PTL,VertexNumList * _V,PolyPointList * _P)201 void REC_Dir_Product(PartList *_PTL, VertexNumList *_V, PolyPointList *_P){
202
203 int i, j, k;
204 CEqList *_CEtemp;
205 VertexNumList *_Vtemp;
206 PolyPointList *_PV;
207 Subset Pset, CPset;
208
209 _PV = (PolyPointList *) malloc(sizeof(PolyPointList));
210 if(_PV == NULL) Die("Unable to alloc space for PolyPointList _PV");
211 _CEtemp = (CEqList *) malloc(sizeof(CEqList));
212 if(_CEtemp == NULL) Die("Unable to alloc space for CEqList _CEtemp");
213 _Vtemp = (VertexNumList *) malloc(sizeof(VertexNumList));
214 if(_Vtemp == NULL) Die("Unable to alloc space for VertexNumList _Vtemp");
215
216 for(i = 0; i < _PTL->n; i++){
217 _PTL->DirProduct[i] = 0; j = 0;
218 /*printf("\n*********************************\n");*/
219 while((j < _PTL->codim) && !_PTL->DirProduct[i]){
220 for(k = 0; k < _PTL->codim; k++){
221 if(k == j){
222 Pset.A[k] = 1; CPset.A[k] = 0;
223 }
224 else{
225 Pset.A[k] = 0; CPset.A[k] = 1;
226 }
227 }
228 Pset.m = 1; Pset.M = _PTL->codim;
229 CPset.m = (_PTL->codim - 1); CPset.M = _PTL->codim;
230 Select_Set(_PTL->S[i], _V, _P, _PV, _CEtemp, _Vtemp, &_PTL->DirProduct[i],
231 &Pset, &CPset, 0);
232 j++;
233 }
234 }
235 free(_Vtemp); free(_CEtemp); free(_PV);
236 }
237
COMP_S(int SA[],int SB[],int * _nv)238 int COMP_S(int SA[], int SB[], int *_nv){
239
240 int eq_flag = 1, i = 0;
241
242 while((i < *_nv) && eq_flag){
243 if((SA[i] > SB[i]) || (SA[i] < SB[i]))
244 eq_flag = 0;
245 i++;
246 }
247 if(eq_flag)
248 return 0;
249 else{
250 if(SA[i-1] > SB[i-1])
251 return 1;
252 else
253 return -1;
254 }
255 }
256
NForm_S(int S[],int * _nv)257 void NForm_S(int S[], int *_nv){
258
259 int s[POLY_Dmax], d=0, n_flag, i, j;
260
261 for(i = 0; i < *_nv; i++){
262 n_flag = 1; j = 0;
263 while((j < i) && n_flag){
264 if(S[j] == S[i])
265 n_flag = 0;
266 j++;
267 }
268 if(n_flag){
269 s[S[i]] = d;
270 d++;
271 }
272 }
273 for(i = 0; i < *_nv; i++)
274 S[i] = s[S[i]];
275 }
276
Bisection_PTL(PartList * _PTL,int s[],int S[])277 int Bisection_PTL(PartList *_PTL, int s[], int S[]){
278
279 int min_pos = -1, max_pos = _PTL->n, pos = 1, c = 1;
280
281 while((max_pos - min_pos ) > 1){
282 pos = (max_pos + min_pos) / 2;
283 c = COMP_S(S, _PTL->S[s[pos]], &_PTL->nv);
284 if (c == 1)
285 min_pos = pos;
286 else if (c == -1)
287 max_pos = pos;
288 else
289 min_pos = max_pos;
290 }
291 assert (c == 0);
292 return s[pos];
293 }
294
Bubble_PTL(PartList * _PTL,int s[])295 void Bubble_PTL(PartList *_PTL, int s[]){
296
297 int i, j, diff=0;
298
299 for(i = 0; i < _PTL->n; i++)
300 s[i] = i;
301 for(i = 0; i < _PTL->n - 1; ++i)
302 for(j = _PTL->n - 1; j > i; --j)
303 if(COMP_S(_PTL->S[s[j-1]], _PTL->S[s[j]], &_PTL->nv) == 1)
304 swap(&s[j-1], &s[j]);
305 for(i = 1; i < _PTL->n; i++){
306 if(COMP_S(_PTL->S[s[i-1]], _PTL->S[s[i]], &_PTL->nv) == 0)
307 diff += 1;
308 else
309 s[i - diff] = s[i];
310 }
311 _PTL->n = _PTL->n - diff;
312 }
313
Remove_Sym(SYM * _VP,PartList * _PTL,PartList * _S_PTL)314 void Remove_Sym(SYM *_VP, PartList *_PTL, PartList *_S_PTL){
315
316 int *_s, *_s_sym, *_p;
317 int n=1, i, j, k;
318 PartList *_SYM_PTL;
319
320 _s = (int*) calloc(_PTL->n, sizeof(int));
321 assert(_s != NULL);
322 _s_sym = (int*) calloc(_VP->ns, sizeof(int));
323 assert(_s_sym != NULL);
324 _p = (int*) calloc(_PTL->n, sizeof(int));
325 assert(_p != NULL);
326 _SYM_PTL = (PartList*) malloc(sizeof(PartList));
327 assert(_SYM_PTL != NULL);
328
329 if(Nef_Max < SYM_Nmax) Die("\nNeed Nef_Max >= SYM_Nmax!!!\n");
330 for(i = 0; i < _PTL->n; i++){
331 _p[i] = 0; _s[i] = 0;
332 }
333 for(i = 0; i < _PTL->n; i++)
334 NForm_S(_PTL->S[i], &_PTL->nv);
335 Bubble_PTL(_PTL, _s);
336 for(i = 0; i < _PTL->n; i++){
337 if(_p[_s[i]] == 0){
338 _p[_s[i]] = n;
339 for(j = 0; j < _VP->ns; j++){
340 for(k = 0; k < _PTL->nv; k++)
341 _SYM_PTL->S[j][k] = _PTL->S[_s[i]][_VP->Vp[j][k]];
342 }
343 _SYM_PTL->n = _VP->ns; _SYM_PTL->nv = _PTL->nv;
344 for(j = 0; j < _VP->ns; j++)
345 NForm_S(_SYM_PTL->S[j], &_SYM_PTL->nv);
346 Bubble_PTL(_SYM_PTL, _s_sym);
347 for(j = 1; j < _VP->ns; j++)
348 if(COMP_S(_SYM_PTL->S[_s_sym[j]], _SYM_PTL->S[_s_sym[j-1]], &_PTL->nv) != 0)
349 _p[Bisection_PTL(_PTL, _s, _SYM_PTL->S[_s_sym[j]])] = n;
350 n++;
351 }
352 }
353 n--; _S_PTL->nv = _PTL->nv; _S_PTL->n = 0; _S_PTL->codim = _PTL->codim;
354 while(n > 0){
355 i = 0;
356 while(_p[_s[i]] != n)
357 i++;
358 for(j = 0; j < _PTL->nv; j ++)
359 _S_PTL->S[_S_PTL->n][j] = _PTL->S[_s[i]][j];
360 _S_PTL->n += 1;
361 n--;
362 }
363 free(_s);free(_s_sym);free(_p);free(_SYM_PTL);
364 }
365
M_TO_MM(MMatrix * _M,MMatrix * _MM,GMatrix * _G,int * _nf)366 void M_TO_MM(MMatrix *_M, MMatrix *_MM, GMatrix *_G, int *_nf){
367
368 int i, l, c, j;
369
370 for(i = 0; i < *_nf; i++){
371 for(l = 0; l < _M[i].d; l++)
372 for(c = 0; c < _M[i].codim; c++){
373 _MM[i].M[l][c] = 0;
374 for(j = 0; j < _M[i].d; j++)
375 _MM[i].M[l][c] += _G[i].G[j][l] * _M[i].M[j][c];
376 }
377 _MM[i].d = _M[i].d; _MM[i].codim =_M[i].codim;
378 }
379 }
380
Convex_Check(MMatrix * _M,GMatrix * _G,XMatrix * _X,int S[],FVList * _FVl,NEF_Flags * _F)381 int Convex_Check(MMatrix *_M, GMatrix *_G, XMatrix *_X, int S[], FVList *_FVl,
382 NEF_Flags *_F){
383 MMatrix *_MM;
384 int c_flag = 1, i, j, k, l, d, IP;
385
386 if (_F->noconvex) return 1;
387
388 _MM = (MMatrix *) calloc(_FVl->nf, sizeof(MMatrix));
389 assert(_MM != NULL);
390
391 M_TO_MM(_M, _MM, _G, &_FVl->nf);
392
393 i = 0;
394 while((i < _FVl->nf) && c_flag){
395 j = 0;
396 while((j < _FVl->nf) && c_flag){
397 if(i != j){
398 k = 0;
399 while((k < _MM[i].codim) && c_flag){
400 l = 0;
401 while((l < _X[i].nv) && c_flag){
402 if(S[_FVl->vl[i].v[l]] == k)
403 IP = 1;
404 else
405 IP = 0;
406 for(d = 0; d < _MM[i].d; d++)
407 IP += - _X[i].X[d][l] * _MM[j].M[d][k];
408 if (IP < 0)
409 c_flag = 0;
410 l++;
411 }
412 k++;
413 }
414 }
415 j++;
416 }
417 i++;
418 }
419 if (c_flag && _F->Test)
420 Print_M(_MM, &_FVl->nf, "M-Matrix:");
421 free(_MM);
422 return c_flag;
423 }
424
Codim_Check(int S[],int * _codim,int * _Nv)425 int Codim_Check(int S[], int *_codim, int *_Nv){
426
427 int gp[POLY_Dmax] = {0}, gp_flag = 1, i, j=0;
428
429 while( (j < *_codim) && gp_flag){
430 i = 0;
431 while( (i < *_Nv) && (gp[j] == 0)){
432 if (S[i] == j)
433 gp[j] = 1;
434 i++;
435 }
436 gp_flag = gp[j];
437 j++;
438 }
439 return gp_flag;
440 }
441
Fix_M(Step * _step,XMatrix * _Y,MMatrix * _M,int S[],M_Rank * _MR,FVList * _FVl)442 int Fix_M(Step *_step, XMatrix *_Y, MMatrix *_M, int S[], M_Rank *_MR,
443 FVList *_FVl){
444
445 int f_flag=1, i=0, j, m, IP;
446
447 while(f_flag && (i!=_M[_step->f].codim)){
448 if(i == S[_FVl->vl[_step->f].v[_step->v]])
449 IP = 1;
450 else
451 IP = 0;
452 for(j = 0; j < _MR->m[_step->f]; j++)
453 IP += -_Y[_step->f].X[j][_step->v] * _M[_step->f].M[j][i];
454 m = (IP / _Y[_step->f].X[_MR->m[_step->f]][_step->v]);
455 if((IP - m * _Y[_step->f].X[_MR->m[_step->f]][_step->v]) != 0)
456 f_flag=0;
457 else
458 _M[_step->f].M[_MR->m[_step->f]][i] = m;
459 i++;
460 }
461 return f_flag;
462 }
463
Check_Consistence(Step * _step,XMatrix * _Y,MMatrix * _M,int S[],M_Rank * _MR,FVList * _FVl)464 int Check_Consistence(Step *_step, XMatrix *_Y, MMatrix *_M, int S[],
465 M_Rank *_MR, FVList *_FVl){
466
467 int c_flag=1, i=0, j, IP;
468
469 while(c_flag && (i!=_M[_step->f].codim)){
470 IP = 0;
471 for(j = 0; j < _MR->m[_step->f]; j++)
472 IP += _Y[_step->f].X[j][_step->v] * _M[_step->f].M[j][i];
473 if(i == S[_FVl->vl[_step->f].v[_step->v]]){
474 if(IP != 1)
475 c_flag = 0;
476 }
477 else{
478 if(IP != 0)
479 c_flag = 0;
480 }
481 i++;
482 }
483 return c_flag;
484 }
485
New_V(V_Flag * _VF,int * _i)486 int New_V(V_Flag *_VF, int *_i){
487
488 int new_flag = 0;
489
490 if (!_VF->s[*_i]){
491 _VF->s[*_i] = 1;
492 new_flag = 1;
493 }
494 return new_flag;
495 }
496
Next_Step(FVList * _FVl,Step * _step)497 int Next_Step(FVList *_FVl, Step *_step){
498
499 int step_flag = 1;
500
501 if(_step->v == _FVl->vl[_step->f].nv-1){
502 _step->f += 1;
503 _step->v = 0;
504 }
505 else
506 _step->v += 1;
507 if(_step->f == _FVl->nf)
508 step_flag = 0;
509 return step_flag;
510 }
511
New_VFlag(V_Flag * _VF,int * _n)512 void New_VFlag(V_Flag *_VF, /* int *_Nv, */ int *_n){
513
514 _VF->s[*_n] = 1;
515 }
516
Old_VFlag(V_Flag * _VF,int * _n)517 void Old_VFlag(V_Flag *_VF, /* int *_Nv, */ int *_n){
518
519 _VF->s[*_n] = 0;
520 }
521
Raise_M_Rank(M_Rank * _MR,int * _facet)522 void Raise_M_Rank(M_Rank *_MR, int *_facet){
523
524 _MR->m[*_facet] += 1;
525 }
526
Lower_M_Rank(M_Rank * _MR,int * _facet)527 void Lower_M_Rank(M_Rank *_MR, int *_facet){
528
529 _MR->m[*_facet] -= 1;
530 }
531
Zero_M_Rank(M_Rank * _MR,int * _facet)532 void Zero_M_Rank(M_Rank *_MR, int *_facet){
533 _MR->m[*_facet] = 0;
534 }
535
Initial_Conditions(MMatrix * _M,XMatrix * _Y,M_Rank * _MR,Step * _step,FVList * _FVl,V_Flag * _VF,int S[],int * _codim,int * _dim,PartList * _PTL)536 void Initial_Conditions(MMatrix *_M, XMatrix *_Y, M_Rank *_MR, Step *_step,
537 FVList *_FVl, V_Flag *_VF, int S[], int *_codim,
538 int *_dim, PartList *_PTL){
539
540 int i;
541 for(i = 0; i < _FVl->nf; i++){
542 _M[i].codim = *_codim;
543 _M[i].d = *_dim;
544 _MR->m[i] = 0;
545 }
546 _step->f = 0; _step->v = 0;
547 for(i = 0; i < _FVl->Nv; i++)
548 _VF->s[i] = 0;
549 _VF->s[_FVl->vl[0].v[0]] = 1;
550 _MR->m[0] = 1;
551 assert((_Y[0].X[0][0] == 1) || (_Y[0].X[0][0] == -1));
552 for(i = 0; i < *_codim; i++)
553 _M[0].M[0][i] = 0;
554 _M[0].M[0][0] = _Y[0].X[0][0];
555 S[_FVl->vl[0].v[0]] = 0;
556 _PTL->n = 0;
557 _PTL->nv = _FVl->Nv;
558 _PTL->codim = *_codim;
559 }
560
INCI_To_VList(INCI * _X,VList * _Vl,int * _N)561 void INCI_To_VList( INCI *_X, VList *_Vl, int *_N){
562 /* INCI X -> {3, 5, 7, ...} ... List of Vertices in
563 VertexNumList of corresponding face */
564
565 int i, n=0;
566 INCI Y = *_X;
567
568 for (i = 0; i < *_N; i++) {
569 if (INCI_M2(Y)){
570 _Vl->v[n] = *_N-i-1;
571 n++;
572 }
573 Y = INCI_D2(Y);
574 }
575 _Vl->nv = n;
576 }
577
INCI_To_FVList(FaceInfo * _I,PolyPointList * _P,FVList * _FVl)578 void INCI_To_FVList(FaceInfo *_I, PolyPointList *_P, FVList *_FVl){
579 /* FaceInfo I -> {3, 5, 7, ...}, {1, 4, 8, ...}, ...
580 Lists of Vertices in VertexNumList of all facets */
581
582 int i;
583
584 for (i = 0; i < _I->nf[_P->n - 1]; i++)
585 INCI_To_VList( &_I->v[_P->n - 1][i], &_FVl->vl[i], &_I->nf[0]);
586 _FVl->nf = _I->nf[_P->n - 1];
587 _FVl->Nv = _I->nf[0];
588 }
589
Sort_FVList(FVList * _FVl,FVList * _FVl_new,int f[])590 void Sort_FVList(FVList *_FVl, FVList *_FVl_new, int f[]){
591 /* gives f[] = LIST OF SORTED FACETS IN _FVl:
592 f[0]: facet with max # of Vertices
593 f[k != 0]: facet with max # of Vertices in f[0],...,f[k-1]
594 _FVl->vl[k] -> _FVl->vl[f[k]] =
595 {v in f[0],...,f[k-1]}U{v not in f[0],...,f[k-1]} =: FVl_new->vl[k] */
596
597 int i, j, k, l, n, Vequal, vequal, newfacet, newvertex,
598 v[VERT_Nmax], w[VERT_Nmax], nv, nw;
599
600 f[0] = 0;
601 for (i = 1; i < _FVl->nf; i++) /* first facet */
602 if (_FVl->vl[i].nv > _FVl->vl[f[0]].nv)
603 f[0] = i;
604 for(i = 0; i < _FVl->vl[f[0]].nv; i++)
605 v[i] = _FVl->vl[f[0]].v[i];
606 nv = _FVl->vl[f[0]].nv;
607 for (n = 1; n < _FVl->nf; n++){ /* next facet */
608 Vequal = 0;
609 for (i = 0; i < _FVl->nf; i++){
610 newfacet = 1;
611 k = 0;
612 while((k < n) && newfacet){
613 if (f[k] == i)
614 newfacet = 0;
615 k++;
616 }
617 if(newfacet){
618 vequal = 0;
619 for(j = 0; j < _FVl->vl[i].nv; j++){
620 newvertex = 1;
621 k = 0;
622 while((k < nv) && newvertex){
623 if (v[k] == _FVl->vl[i].v[j])
624 newvertex = 0;
625 k++;
626 }
627 if(!newvertex)
628 vequal++;
629 }
630 if(vequal >= Vequal){
631 Vequal = vequal;
632 f[n] = i;
633 }
634 }
635 }
636 nw = 0;
637 for(j = 0; j < _FVl->vl[f[n]].nv; j++){
638 newvertex = 1;
639 k = 0;
640 while((k < nv) && newvertex){
641 if (v[k] == _FVl->vl[f[n]].v[j])
642 newvertex = 0;
643 k++;
644 }
645 if(newvertex){
646 v[nv] = _FVl->vl[f[n]].v[j];
647 nv++;
648 }
649 else{
650 w[nw] = _FVl->vl[f[n]].v[j];
651 nw++;
652 }
653 }
654 for(l = 0; l < nw; l++)
655 _FVl->vl[f[n]].v[l] = w[l];
656 for(l = 0; l < _FVl->vl[f[n]].nv - nw; l++)
657 _FVl->vl[f[n]].v[nw+l] = v[nv-l-1];
658 }
659 for(i = 0; i < _FVl->nf; i++){
660 for(j = 0; j <_FVl->vl[f[i]].nv; j++)
661 _FVl_new->vl[i].v[j] = _FVl->vl[f[i]].v[j];
662 _FVl_new->vl[i].nv = _FVl->vl[f[i]].nv;
663 }
664 _FVl_new->nf = _FVl->nf;
665 _FVl_new->Nv = _FVl->Nv;
666 }
667
Make_Matrix(XMatrix * _X,XMatrix * _Y,VList * _Vl,PolyPointList * _P,VertexNumList * _V)668 void Make_Matrix(XMatrix * _X, XMatrix * _Y, VList *_Vl,
669 PolyPointList *_P, VertexNumList *_V){
670 /* Vl -> Matrix X */
671
672 int d, j;
673
674 for (d=0; d < _P->n; d++)
675 for (j=0; j < _Vl->nv; j++){
676 _X->X[d][j] = _P->x[_V->v[_Vl->v[j]]][d];
677 _Y->X[d][j] = _P->x[_V->v[_Vl->v[j]]][d];
678 }
679 _X->d = _P->n; _X->nv = _Vl->nv;
680 _Y->d = _P->n; _Y->nv = _Vl->nv;
681 }
682
Copy_PTL(PartList * _IN_PTL,PartList * _OUT_PTL)683 void Copy_PTL(PartList *_IN_PTL, PartList *_OUT_PTL){
684
685 int i, j;
686 for(i = 0; i < _IN_PTL->n; i++)
687 for(j = 0; j < _IN_PTL->nv; j++)
688 _OUT_PTL->S[i][j] = _IN_PTL->S[i][j];
689 _OUT_PTL->n = _IN_PTL->n; _OUT_PTL->nv = _IN_PTL->nv;
690 _OUT_PTL->codim = _IN_PTL->codim;
691 }
692
Select_Sv(int S[],V_Flag * _VF,MMatrix * _M,GMatrix * _G,XMatrix * _X,XMatrix * _Y,M_Rank * _MR,FVList * _FVl,Step step,PartList * _PTL,NEF_Flags * _F)693 void Select_Sv(int S[], V_Flag *_VF, MMatrix *_M, GMatrix *_G, XMatrix *_X,
694 XMatrix *_Y, M_Rank * _MR, FVList *_FVl, Step step,
695 PartList *_PTL, NEF_Flags* _F){
696
697 int i;
698 if(Next_Step(_FVl, &step)){
699 if(step.v == 0)
700 Zero_M_Rank(_MR, &step.f);
701 if(!New_V(_VF, &_FVl->vl[step.f].v[step.v])){
702 if (_F->Test){
703 char c;
704 fprintf(outFILE, "\nold vertex f:%d v:%d\n",step.f,step.v);
705 scanf("%c",&c);
706 }
707 if((_MR->m[step.f] == _M[step.f].d) ||
708 (_Y[step.f].X[_MR->m[step.f]][step.v] == 0)){
709 if(Check_Consistence(&step, _Y, _M, S, _MR, _FVl)){
710 if (_F->Test)
711 fprintf(outFILE, "f:%d v:%d consistent\n",step.f,step.v);
712 Select_Sv(S, _VF, _M, _G, _X, _Y, _MR, _FVl, step, _PTL, _F);
713 }
714 }
715 else{
716 if(Fix_M(&step, _Y, _M, S, _MR, _FVl)){
717 if (_F->Test)
718 fprintf(outFILE, "f:%d v:%d fixed\n",step.f,step.v);
719 Raise_M_Rank(_MR, &step.f);
720 Select_Sv(S, _VF, _M, _G, _X, _Y, _MR, _FVl, step, _PTL, _F);
721 }
722 }
723 }
724 else{
725 if (_F->Test){
726 fprintf(outFILE, "\nnew vertex f:%d v:%d",step.f,step.v);
727 }
728 New_VFlag(_VF, /* &_FVl->Nv,*/ &_FVl->vl[step.f].v[step.v]);
729 for(i = 0; i < _M[step.f].codim; i++){
730 S[_FVl->vl[step.f].v[step.v]] = i;
731 if (_F->Test)
732 fprintf(outFILE, "to partition: %d\n",i);
733 if((_MR->m[step.f] == _M[step.f].d) ||
734 (_Y[step.f].X[_MR->m[step.f]][step.v] == 0)){
735 if(Check_Consistence(&step, _Y, _M, S, _MR, _FVl)){
736 if (_F->Test)
737 fprintf(outFILE, "f:%d v:%d consistent\n",step.f,step.v);
738 Select_Sv(S, _VF, _M, _G, _X, _Y, _MR, _FVl, step, _PTL, _F);
739 }
740 }
741 else{
742 if(Fix_M(&step, _Y, _M, S, _MR, _FVl)){
743 if (_F->Test)
744 fprintf(outFILE, "f:%d v:%d fixed\n",step.f,step.v);
745 Raise_M_Rank(_MR, &step.f);
746 Select_Sv(S, _VF, _M, _G, _X, _Y, _MR, _FVl, step, _PTL, _F);
747 Lower_M_Rank(_MR, &step.f);
748 }
749 }
750 }
751 Old_VFlag(_VF, /* &_FVl->Nv, */ &_FVl->vl[step.f].v[step.v]);
752 }
753 }
754 else{
755 if(_F->Test){
756 fprintf(outFILE, "\n********************************************\n");
757 for(i = 0; i < _FVl->Nv; i++)
758 fprintf(outFILE, " %d ", S[i]);
759 fprintf(outFILE, "\n");
760 fprintf(outFILE, "**********************************************\n");
761 fflush(0);
762 }
763 if (Codim_Check(S, &_M[step.f-1].codim, &_FVl->Nv))
764 if (Convex_Check(_M, _G, _X, S, _FVl, _F)){
765 assert(_PTL->n < Nef_Max);
766 for(i = 0; i < _FVl->Nv; i++)
767 _PTL->S[_PTL->n][i] = S[i];
768 _PTL->n += 1;
769 }
770 }
771 }
772
part_nef(PolyPointList * _P,VertexNumList * _V,EqList * _E,PartList * _OUT_PTL,int * _codim,NEF_Flags * _F)773 void part_nef(PolyPointList *_P, VertexNumList *_V, EqList *_E,
774 PartList *_OUT_PTL, int *_codim, NEF_Flags* _F){
775
776 FaceInfo I;
777 FVList FVl_temp, FVl;
778 XMatrix *_X, *_Y;
779 MMatrix *_M;
780 GMatrix *_G;
781 Step step;
782 V_Flag VF;
783 M_Rank MR;
784 int i, f[FACE_Nmax], S[VERT_Nmax] = {0};
785
786 PartList *_PTL = (PartList *) malloc(sizeof(PartList));
787 assert(_PTL != NULL);
788
789 Make_Incidence(_P, _V, _E, &I);
790
791 FVl.vl = (VList *) calloc(I.nf[_P->n - 1], sizeof(VList));
792 assert(FVl.vl != NULL);
793
794 if (_F->Sort){
795 FVl_temp.vl = (VList *) calloc(I.nf[_P->n - 1], sizeof(VList));
796 assert(FVl_temp.vl != NULL);
797 INCI_To_FVList(&I, _P, &FVl_temp);
798 Sort_FVList(&FVl_temp, &FVl, f);
799 free(FVl_temp.vl);
800 }
801 else
802 INCI_To_FVList(&I, _P, &FVl);
803 if (_F->Test){
804 Print_VL(_P, _V, "Vertices of P:");
805 Print_FVl(&FVl, "Facets/Vertices:");
806 }
807 _X = (XMatrix *) calloc(FVl.nf, sizeof(XMatrix));
808 assert(_X != NULL);
809 _Y = (XMatrix *) calloc(FVl.nf, sizeof(XMatrix));
810 assert(_Y != NULL);
811 _M = (MMatrix *) calloc(FVl.nf, sizeof(MMatrix));
812 assert(_M != NULL);
813 _G = (GMatrix *) calloc(FVl.nf, sizeof(GMatrix));
814 assert(_G != NULL);
815
816 for(i = 0; i < FVl.nf; i++){
817 Make_Matrix(&_X[i], &_Y[i], &FVl.vl[i], _P, _V);
818 GLZ_Make_Trian_NF(_Y[i].X, &_P->n, &FVl.vl[i].nv, _G[i].G);
819 }
820 Initial_Conditions(_M, _Y, &MR, &step, &FVl, &VF, S, _codim, &_P->n, _PTL);
821 Select_Sv(S, &VF, _M, _G, _X, _Y, &MR, &FVl, step, _PTL, _F);
822 free(_X); free(_Y); free(_M); free(_G); free(FVl.vl);
823 if(_F->Sym){
824 SYM *_VP = (SYM *) malloc(sizeof(SYM));
825 assert(_VP != NULL);
826
827 Poly_Sym(_P, _V, _E, &_VP->ns, _VP->Vp);
828 Remove_Sym(_VP, _PTL, _OUT_PTL);
829 free(_VP);
830 }
831 else
832 Copy_PTL(_PTL, _OUT_PTL);
833 /*Dir_Product(_OUT_PTL, _V, _P);*/
834 if(*_codim > 1)
835 REC_Dir_Product(_OUT_PTL, _V, _P);
836 free(_PTL);
837 }
838