CoCalc -- refine

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_fuber.c
⁵⁶³⁶⁶⁵ views
1
#include "typedef.h"
2
#include "tools.h"
3
#include "matrix.h"
4

5
void vertex_standard(v)
6
vertex_TYP *v;
7
{
8
  int i,j;
9
  int w1;
10
 
11
  if(v->dim>0)
12
  {
13
    w1 = 0;
14
    i=0;
15
    while(i<v->dim && v->v[i] == 0)
16
       i++;
17
    if(i<v->dim)
18
      w1 = v->v[i];
19
    for(j=i+1;j<v->dim;j++)
20
    {
21
     if(v->v[j] != 0)
22
       w1 = GGT(w1, v->v[j]);
23
    }
24
    if(w1 < 0)
25
      w1 = -w1;
26
    if(w1 != 0)
27
    {
28
       for(j=i;j<v->dim;j++)
29
         v->v[j] /=w1;
30
    }
31
  }
32
}
33

34
void umsortieren(F,old_no, count, test, l)
35
fund_domain *F;
36
int old_no, *test, l;
37
int *count;
38
{
39
  int i;
40
  int d=0;
41
  for(i=0;i<old_no;i++)
42
  {
43
    if(count[i] < 0)
44
      free_vertex_fuber(&(F->vert[i]));
45
  }
46
  for(i=0;d<F->vert_no;i++)
47
  {
48
    while(d<F->vert_no && F->vert[d] == NULL)
49
      d++;
50
    if(d != i)
51
    {
52
      F->vert[i] = F->vert[d];
53
      F->vert[d] = NULL;
54
    }
55
    d++;
56
  }
57
  while(F->vert[i-1] == NULL)
58
    i--;
59
  F->vert_no = i;
60

61
  d=0;
62
  for(i=0;i<F->wall_no;i++)
63
  {
64
    if(test[i] == -1)
65
      free_wall_fuber(&(F->wall[i]));
66
  }
67
  for(i=0;d<F->wall_no;i++)
68
  {
69
    while(d<F->wall_no && F->wall[d] == NULL)
70
      d++;
71
    if(d!= i)
72
    {
73
      F->wall[i] = F->wall[d];
74
      F->wall[d] = NULL;
75
    }
76
    d++;
77
  }
78
  F->wall_no = i;
79
}
80

81

82
void renumerate(v,test,l)
83
vertex_TYP *v;
84
int *test,l;
85
{
86
  int i,a;
87
  for(i=0;i<v->wall_no;i++)
88
  {
89
     a = v->wall[i];
90
     v->wall[i] = test[a];
91
  }
92
}
93

94

95
void wallAdd_vertex(i, v)
96
int i;
97
vertex_TYP *v;
98
{
99
  if(v->wall_SIZE == 0)
100
  {
101
    v->wall = (int *)malloc(extsize1 *sizeof(int));
102
    v->wall_SIZE = extsize1;
103
  }
104
  if(v->wall_no == v->wall_SIZE)
105
  {
106
    v->wall_SIZE += extsize1;
107
    v->wall = (int *)realloc(v->wall,v->wall_SIZE *sizeof(int));
108
  }
109
  v->wall[v->wall_no] = i;
110
  v->wall_no++;
111
}
112

113

114
void streichen(v, test)
115
vertex_TYP *v;
116
int *test;
117
{
118
 int i,d;
119

120
 d=0;
121
 for(i=0; d<v->wall_no;i++)
122
 {
123
   while(d< v->wall_no && test[v->wall[d]] == -1)
124
     d++;
125
   v->wall[i] = v->wall[d];
126
   d++;
127
 }
128
 v->wall_no = i;
129
}
130
  
131

132
static int is_element(v, w)
133
vertex_TYP *v;
134
int w;
135
{
136
  int o,u,t;
137
  o=w;
138
  u=0;
139
  if(o>= v->wall_no)
140
    o=v->wall_no-1;
141
  while(o>u)
142
  {
143
    t=(o+u)/2;
144
    if(v->wall[t] < w)
145
      u=t+1;
146
    else
147
      o=t;
148
  } 
149
  if(w == v->wall[u])
150
    return(TRUE);
151
  return(FALSE);
152
}
153

154

155

156
void wallAdd_fund_domain(F, w)
157
fund_domain *F;
158
wall_TYP *w;
159
{
160
  if(F->wall_SIZE == 0)
161
  {
162
    F->wall = (wall_TYP **)malloc(extsize1 *sizeof(wall_TYP *));
163
    F->wall_SIZE = extsize1;
164
  }
165
  if(F->wall_no == F->wall_SIZE)
166
  {
167
    F->wall_SIZE += extsize1;
168
    F->wall = (wall_TYP **)realloc( (int *)F->wall,F->wall_SIZE *sizeof(wall_TYP *));
169
  }
170
  F->wall[F->wall_no] = w;
171
  F->wall_no++;
172
}
173

174

175

176
void vertexAdd_fund_domain(F, w)
177
fund_domain *F;
178
vertex_TYP *w;
179
{
180
  if(F->vert_SIZE == 0)
181
  {
182
    F->vert = (vertex_TYP **)malloc(extsize1 *sizeof(vertex_TYP *));
183
    F->vert_SIZE = extsize1;
184
  }
185
  if(F->vert_no == F->vert_SIZE)
186
  {
187
    F->vert_SIZE += extsize1;
188
    F->vert = (vertex_TYP **)realloc( (int *)F->vert,F->vert_SIZE *sizeof(vertex_TYP *));
189
  }
190
  F->vert[F->vert_no] = w;
191
  F->vert_no++;
192
}
193

194

195
vertex_TYP *gis_neighbour(i,j,F,old_no)
196
int i,j;
197
fund_domain *F;
198
int old_no;
199
{
200
  int k,l,m;
201
  int a, u, tester1;
202
  int anz;
203
  vertex_TYP *erg;
204
  matrix_TYP *A;
205

206
  extern vertex_TYP *init_vertex_fuber();
207
  extern matrix_TYP *init_mat();
208
  extern int row_gauss();
209

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

228
  A = init_mat(anz,F->vert[0]->dim, "");
229
  for(k=0;k<anz;k++)
230
    for(l=0; l<F->vert[0]->dim;l++)
231
      A->array.SZ[k][l] = F->wall[erg->wall[k]]->gl[l];
232
  a = row_gauss(A);
233
  free_mat(A);
234
  if(a < F->vert[0]->dim-2)
235
   { free_vertex_fuber(&erg); erg=NULL; return(NULL);}
236
  return(erg);
237
}
238

239

240

241
vertex_TYP *is_neighbour(i,j,F,old_no)
242
int i,j;
243
fund_domain *F;
244
int old_no;
245
{
246
  int k,l,m;
247
  int a, u, tester1;
248
  int anz;
249
  vertex_TYP *erg;
250

251
  extern vertex_TYP *init_vertex_fuber();
252

253
  k=F->vert[j]->wall_no;
254
  if(k>F->vert[i]->wall_no)
255
    k=F->vert[i]->wall_no;
256
  erg = init_vertex_fuber(F->vert[0]->dim, k+1);
257
  anz=0;
258
  u=0;
259
  for(k=0;k<F->vert[i]->wall_no;k++)
260
  {
261
    a = F->vert[i]->wall[k];
262
    for(l=u; l<F->vert[j]->wall_no && F->vert[j]->wall[l]< a; l++);
263
    u=l;
264
    if(u<F->vert[j]->wall_no && F->vert[j]->wall[u] == a)
265
       { erg->wall[anz] = a; anz++;}
266
  }
267
  erg->wall_no = anz;
268
  if(anz < F->vert[0]->dim-2)
269
   { free_vertex_fuber(&erg); erg=NULL; return(NULL);}
270

271
  for(k=0; k<old_no ;k++)
272
  {
273
   if(k != i && k != j)
274
   {
275
     tester1 = TRUE;
276
        u=0;
277
     for(l=0; l<anz && tester1 == TRUE; l++)
278
     {
279
        a = erg->wall[l];
280
        m=u;
281
        while(m<F->vert[k]->wall_no && F->vert[k]->wall[m] <a)
282
           m++;
283
        if(m == F->vert[k]->wall_no || F->vert[k]->wall[m] != a)
284
           tester1 = FALSE;
285
        u=m;
286
     }
287
     if(tester1 == TRUE)
288
       { free_vertex_fuber(&erg); erg=NULL; return(NULL);}
289
   }
290
  }
291
  return(erg);
292
}
293

294

295

296
int refine_fund_domain(F, h)
297
fund_domain *F;
298
wall_TYP *h;
299
{
300
 int i,j,k;
301
 int l;
302
 int old_no;
303
 int *count;
304
 int waste;
305
 int tester, *test;
306
 int anz;
307
 int p,n;
308
 vertex_TYP *v;
309

310
 extern void umsortieren();
311
 extern void renumerate();
312
 extern void streichen();
313
 extern int wall_times_vertex_fuber();
314
 extern int containes();
315
 extern void  wallAdd_fund_domain();
316
 extern void  vertexAdd_fund_domain();
317
 extern vertex_TYP *is_neighbour();
318
 
319

320
 p=0;n=0;
321
 waste = 0;
322
 old_no = F->vert_no;
323
 count = (int *)malloc(F->vert_no *sizeof(int *));
324
 for(i=0;i<F->vert_no;i++)
325
  count[i] = 0;
326
 for(i=0;i<old_no;i++)
327
 {
328
     count[i] = wall_times_vertex_fuber(h,F->vert[i]);
329
    if(count[i] > 0)
330
      p++;
331
    if(count[i] < 0)
332
      n++;
333
 }
334
 if(n == 0)
335
 {
336
   free(count);
337
   return(0);
338
 }
339
 for(i=0;i<old_no;i++)
340
 {
341
  if(count[i] < 0)
342
  {  
343
    for(j=0;j<old_no;j++)
344
    {
345
     if(count[j] > 0)
346
     {
347
        v = is_neighbour(i, j, F, old_no);
348
        if(v != NULL)
349
        {
350
           for(k=0; k<v->dim; k++)
351
           {
352
             v->v[k] = count[j] * F->vert[i]->v[k];
353
             v->v[k] -= count[i] * F->vert[j]->v[k];
354
           }
355
           vertex_standard(v);
356
           vertexAdd_fund_domain(F, v);
357
        }
358
     }
359
    }
360
  }
361
 }
362
  wallAdd_fund_domain(F, h);
363
  test = (int *)malloc((F->wall_no) *sizeof(int));
364
  l=0;
365
  for(i=0; i<F->wall_no-1;i++)
366
  {
367
    tester = 0;
368
    for(j=0;j<old_no && tester == 0; j++)
369
    {
370
     if(count[j] > 0)
371
       if(is_element(F->vert[j], i) == TRUE)
372
          tester = 1;
373
    }
374
    if(tester == 0)
375
     test[i] = -1;
376
    else
377
     {test[i] = l; l++;}
378
  }
379
  test[F->wall_no-1] = l;
380
  l++;
381

382
  for(i=0;i<old_no; i++)
383
  {
384
    if(count[i] == 0)
385
      wallAdd_vertex(F->wall_no-1, F->vert[i]);
386
  }
387

388
  for(i=0;i<old_no;i++)
389
  {
390
     if(count[i] ==0)
391
       streichen(F->vert[i],test );
392
     if(count[i] >=0)
393
       renumerate(F->vert[i],test,l);
394
  }
395
  for(i=old_no; i<F->vert_no;i++)
396
  {
397
     F->vert[i]->wall[F->vert[i]->wall_no] = F->wall_no-1;
398
     F->vert[i]->wall_no++;
399
     renumerate(F->vert[i],test,l);
400
  }
401

402
  umsortieren(F, old_no, count, test,l);
403
  
404
   free(count);
405
   free(test);
406
  return(1);
407
}
408
/*{{{}}}*/
409

410
Product

Resources

Company
