CoCalc -- refine_polyeder.c

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
gap4r8 / pkg / carat / carat-2.1b1 / functions / Polyeder / refine_polyeder.c
⁵⁶³⁶⁰³ views
1
#include"typedef.h"
2
#include"matrix.h"
3
#include"longtools.h"
4
#include"polyeder.h"
5
/**************************************************************************\
6
@---------------------------------------------------------------------------
7
@---------------------------------------------------------------------------
8
@ FILE: refine_polyeder.c
9
@---------------------------------------------------------------------------
10
@---------------------------------------------------------------------------
11
@
12
\**************************************************************************/
13

14

15
static int g_option;
16

17
static void re_sort(F,old_no, count, test, l)
18
polyeder_TYP *F;
19
int old_no, *test, l;
20
int *count;
21
{
22
  int i;
23
  int d=0;
24
  for(i=0;i<old_no;i++)
25
  {
26
    if(count[i]< 0)
27
      free_vertex(&(F->vert[i]));
28
  }
29
  for(i=0;d<F->vert_no;i++)
30
  {
31
    while(d<F->vert_no && F->vert[d] == NULL)
32
      d++;
33
    if(d != i)
34
    {
35
      F->vert[i] = F->vert[d];
36
      F->vert[d] = NULL;
37
    }
38
    d++;
39
  }
40
  while(F->vert[i-1] == NULL)
41
    i--;
42
  F->vert_no = i;
43

44
  d=0;
45
  for(i=0;i<F->wall_no;i++)
46
  {
47
    if(test[i] == -1)
48
      free_wall(&(F->wall[i]));
49
  }
50
  for(i=0;d<F->wall_no;i++)
51
  {
52
    while(d<F->wall_no && F->wall[d] == NULL)
53
      d++;
54
    if(d!= i)
55
    {
56
      F->wall[i] = F->wall[d];
57
      F->wall[d] = NULL;
58
    }
59
    d++;
60
  }
61
  F->wall_no = i;
62
}
63

64

65
static void renumerate(v,test)
66
vertex_TYP *v;
67
int *test;
68
{
69
  int i,a;
70
  for(i=0;i<v->wall_no;i++)
71
  {
72
     a = v->wall[i];
73
     v->wall[i] = test[a];
74
  }
75
}
76

77

78
static void add_wall_to_vertex(i, v)
79
int i;
80
vertex_TYP *v;
81
{
82
  if(v->wall_SIZE == 0)
83
  {
84
    v->wall = (int *)malloc(extsize1 *sizeof(int));
85
    v->wall_SIZE = extsize1;
86
  }
87
  if(v->wall_no == v->wall_SIZE)
88
  {
89
    v->wall_SIZE += extsize1;
90
    v->wall = (int *)realloc(v->wall,v->wall_SIZE *sizeof(int));
91
  }
92
  v->wall[v->wall_no] = i;
93
  v->wall_no++;
94
}
95

96

97
static void delete_walls_from_vertex(v, test)
98
vertex_TYP *v;
99
int *test;
100
{
101
 int i,d;
102

103
 d=0;
104
 for(i=0; d<v->wall_no;i++)
105
 {
106
   while(d< v->wall_no && test[v->wall[d]] == -1)
107
     d++;
108
   v->wall[i] = v->wall[d];
109
   d++;
110
 }
111
 v->wall_no = i;
112
}
113
  
114

115

116
static void add_wall_to_polyeder(F, w)
117
polyeder_TYP *F;
118
wall_TYP *w;
119
{
120
  if(F->wall_SIZE == 0)
121
  {
122
    F->wall = (wall_TYP **)malloc(extsize1 *sizeof(wall_TYP *));
123
    F->wall_SIZE = extsize1;
124
  }
125
  if(F->wall_no == F->wall_SIZE)
126
  {
127
    F->wall_SIZE += extsize1;
128
    F->wall = (wall_TYP **)realloc(F->wall,F->wall_SIZE *sizeof(wall_TYP *));
129
  }
130
  F->wall[F->wall_no] = copy_wall(w);
131
  F->wall_no++;
132
}
133

134

135

136
static void add_vertex_to_polyeder(F, w)
137
polyeder_TYP *F;
138
vertex_TYP *w;
139
{
140
  if(F->vert_SIZE == 0)
141
  {
142
    F->vert = (vertex_TYP **)malloc(extsize1 *sizeof(vertex_TYP *));
143
    F->vert_SIZE = extsize1;
144
  }
145
  if(F->vert_no == F->vert_SIZE)
146
  {
147
    F->vert_SIZE += extsize1;
148
    F->vert = (vertex_TYP **)realloc(F->vert,F->vert_SIZE *sizeof(vertex_TYP *));
149
  }
150
  F->vert[F->vert_no] = w;
151
  F->vert_no++;
152
}
153

154

155
static vertex_TYP *gis_neighbour(i,j,F)
156
int i,j;
157
polyeder_TYP *F;
158
{
159
  int k,l,m;
160
  int a, u, tester1;
161
  int anz;
162
  vertex_TYP *erg;
163
  matrix_TYP *A;
164

165

166
  k=F->vert[j]->wall_no;
167
  if(k>F->vert[i]->wall_no)
168
    k=F->vert[i]->wall_no;
169
  erg = init_vertex(F->vert[0]->dim, k+1);
170
  anz=0;
171
  u=0;
172
  for(k=0;k<F->vert[i]->wall_no;k++)
173
  {
174
    a = F->vert[i]->wall[k];
175
    for(l=u; l<F->vert[j]->wall_no && F->vert[j]->wall[l]< a; l++);
176
    u=l;
177
    if(u<F->vert[j]->wall_no && F->vert[j]->wall[u] == a)
178
       { erg->wall[anz] = a; anz++;}
179
  }
180
  erg->wall_no = anz;
181
  if(anz < F->vert[0]->dim-2)
182
   { free_vertex(&erg); erg=NULL; return(NULL);}
183

184
  A = init_mat(anz,F->vert[0]->dim, "l");
185
  for(k=0;k<anz;k++)
186
    for(l=0; l<F->vert[0]->dim;l++)
187
      A->array.SZ[k][l] = F->wall[erg->wall[k]]->gl[l];
188
  a = long_row_gauss(A);
189
  free_mat(A);
190
  if(a < F->vert[0]->dim-2)
191
   { free_vertex(&erg); erg=NULL; return(NULL);}
192
  return(erg);
193
}
194

195

196

197
static vertex_TYP *is_neighbour(i,j,F, vertex_no)
198
int i,j;
199
polyeder_TYP *F;
200
int vertex_no;
201
{
202
  int k,l,m;
203
  int a, u, tester1;
204
  int anz;
205
  vertex_TYP *erg;
206

207

208
  k=F->vert[j]->wall_no;
209
  if(k>F->vert[i]->wall_no)
210
    k=F->vert[i]->wall_no;
211
  erg = init_vertex(F->vert[0]->dim, k+1);
212
  anz=0;
213
  u=0;
214
  for(k=0;k<F->vert[i]->wall_no;k++)
215
  {
216
    a = F->vert[i]->wall[k];
217
    for(l=u; l<F->vert[j]->wall_no && F->vert[j]->wall[l]< a; l++);
218
    u=l;
219
    if(u<F->vert[j]->wall_no && F->vert[j]->wall[u] == a)
220
       { erg->wall[anz] = a; anz++;}
221
  }
222
  erg->wall_no = anz;
223
  if(anz < F->vert[0]->dim-2)
224
   { free_vertex(&erg); erg=NULL; return(NULL);}
225

226
  for(k=0; k<vertex_no ;k++)
227
  {
228
   if(k != i && k != j)
229
   {
230
     tester1 = TRUE;
231
        u=0;
232
     for(l=0; l<anz && tester1 == TRUE; l++)
233
     {
234
        a = erg->wall[l];
235
        m=u;
236
        while(m<F->vert[k]->wall_no && F->vert[k]->wall[m] <a)
237
           m++;
238
        if(m == F->vert[k]->wall_no || F->vert[k]->wall[m] != a)
239
           tester1 = FALSE;
240
        u=m;
241
     }
242
     if(tester1 == TRUE)
243
       { free_vertex(&erg); erg=NULL; return(NULL);}
244
   }
245
  }
246
  return(erg);
247
}
248

249

250

251

252
/**************************************************************************\
253
@---------------------------------------------------------------------------
254
@ int refine_polyeder(F, h)
255
@ polyeder_TYP *F;
256
@ wall_TYP *h;
257
@ 
258
@  calculates the intersection of the linear Polyeder 'F' with the halfspace
259
@  defined by the inequality
260
@              H :=  h->gl * X^{tr} >= 0
261
@  and stores it in F.
262
@  If the intersection of F and H is already F, 'refine_polyeder returns 0,
263
@  otherwise 1.
264
@  All vertices of the intersecting polyeder are calculated.
265
@  If the intersection has not full dimension, i.e the new Polyeder is
266
@  non-degenerate, F->is_denerate is set to 0.
267
@
268
@  CAUTION: The result is correct only if the new Polyeder is non-degenerate.
269
@
270
@---------------------------------------------------------------------------
271
@
272
\**************************************************************************/
273
int refine_polyeder(F, h)
274
polyeder_TYP *F;
275
wall_TYP *h;
276
{
277
 int i,j,k;
278
 int l;
279
 int old_no;
280
 int *count;
281
 int waste;
282
 int tester, *test;
283
 int anz;
284
 int p,n, z;
285
 vertex_TYP *v;
286
 int *wall_test, *necessary;
287

288
 
289

290
 p=0;n=0;
291
 old_no = F->vert_no;
292
 count = (int *)malloc(F->vert_no *sizeof(int *));
293
 for(i=0;i<F->vert_no;i++)
294
  count[i] = 0;
295
 for(i=0;i<old_no;i++)
296
 {
297
     count[i] = wall_times_vertex(h,F->vert[i]);
298
    if(count[i] > 0)
299
      p++;
300
    if(count[i] < 0)
301
      n++;
302
    if(count[i] == 0)
303
      z++;
304
 }
305
 if(n == 0)
306
 {
307
   free(count);
308
   return(0);
309
 }
310
 if(p == 0)
311
 {
312
   F->is_degenerate = TRUE;
313
 }
314

315

316
 g_option = FALSE;
317
 if(is_option('g') == TRUE)
318
    g_option = TRUE;
319
 wall_test = (int *)malloc(F->wall_no *sizeof(int));
320
 for(i=0;i<F->wall_no;i++)
321
     wall_test[i] = 2;
322
 for(i=0;i<old_no;i++)
323
 {
324
   if(count[i]>0)
325
   {
326
      for(j=0;j<F->vert[i]->wall_no;j++)
327
      {
328
        if(wall_test[F->vert[i]->wall[j]] == 2)
329
          wall_test[F->vert[i]->wall[j]] = 1;
330
        if(wall_test[F->vert[i]->wall[j]] == -1)
331
          wall_test[F->vert[i]->wall[j]] = 0;
332
      }
333
   }
334
   if(count[i] < 0)
335
   {
336
      for(j=0;j<F->vert[i]->wall_no;j++)
337
      {
338
        if(wall_test[F->vert[i]->wall[j]] == 2)
339
          wall_test[F->vert[i]->wall[j]] = -1;
340
        if(wall_test[F->vert[i]->wall[j]] == 1)
341
          wall_test[F->vert[i]->wall[j]] = 0;
342
      }
343
   }
344
 }
345
 necessary = (int *)malloc(F->vert_no *sizeof(int));
346
 for(i=0;i<old_no;i++)
347
   necessary[i] = FALSE;
348
 for(i=0;i<old_no;i++)
349
 {
350
  if(count[i] != 0)
351
  {
352
   k=0;
353
   for(j=0;j<F->vert[i]->wall_no && necessary[i] == FALSE;j++)
354
   {
355
     if( wall_test[F->vert[i]->wall[j]] == 0)
356
       k++;
357
     if(k>= F->vert[i]->dim-2)
358
     {
359
       if(count[i] > 0)
360
         necessary[i] = 1;
361
       if(count[i] < 0)
362
         necessary[i] = -1;
363
       if(count[i] == 0)
364
         necessary[i] = 0;
365
     }
366
   }
367
  }
368
 }
369

370
 if(F->vert[0]->dim == 2)
371
 {
372
   for(i=0;i<old_no;i++)
373
   {
374
     if(count[i] > 0)
375
       necessary[i] = 1;
376
     if(count[i] < 0)
377
       necessary[i] = -1;
378
     if(count[i] == 0)
379
       necessary[i] = 0;
380
   }
381
 }
382

383

384
 for(i=0;i<old_no;i++)
385
 {
386
  if(necessary[i] == -1)
387
  {  
388
    for(j=0;j<old_no;j++)
389
    {
390
     if(necessary[j] == 1)
391
     {
392
        if(g_option == TRUE)
393
          v = gis_neighbour(i, j, F);
394
        else
395
          v = is_neighbour(i, j, F, old_no);
396
        if(v != NULL)
397
        {
398
          for(k=0; k<v->dim; k++)
399
            v->v[k] = count[j] * F->vert[i]->v[k] - count[i] *F->vert[j]->v[k];
400
          normal_vertex(v);
401
          add_vertex_to_polyeder(F, v);
402
        }
403
     }
404
    }
405
  }
406
 }
407
  add_wall_to_polyeder(F, h);
408
  test = (int *)malloc((F->wall_no) *sizeof(int));
409
  l=0;
410
  for(i=0; i<F->wall_no-1;i++)
411
  {
412
    tester = 0;
413
    for(j=0;j<old_no && tester == 0; j++)
414
    {
415
     if(count[j] > 0)
416
       if(is_vertex_of_wallno(F->vert[j], i) == TRUE)
417
          tester = 1;
418
    }
419
    if(tester == 0)
420
     test[i] = -1;
421
    else
422
     {test[i] = l; l++;}
423
  }
424
  test[F->wall_no-1] = l;
425
  l++;
426
  if(test[0] == -1)
427
    F->is_closed = TRUE;
428

429
  for(i=0;i<old_no; i++)
430
  {
431
    if(count[i] == 0)
432
      add_wall_to_vertex(F->wall_no-1, F->vert[i]);
433
  }
434

435
  for(i=0;i<old_no;i++)
436
  {
437
     if(count[i]==0)
438
       delete_walls_from_vertex(F->vert[i],test);
439
     if(count[i]>=0)
440
       renumerate(F->vert[i],test);
441
  }
442
  for(i=old_no; i<F->vert_no;i++)
443
  {
444
     F->vert[i]->wall[F->vert[i]->wall_no] = F->wall_no-1;
445
     F->vert[i]->wall_no++;
446
     renumerate(F->vert[i],test);
447
  }
448

449
  re_sort(F, old_no, count, test,l);
450
  
451
   free(count);
452
   free(test);
453
   free(necessary);
454
   free(wall_test);
455
  return(1);
456
}
457

458
Product

Resources

Company
