CoCalc -- classifyhpel.hpp

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/ElmerGUI/netgen/libsrc/meshing/classifyhpel.hpp
³²⁰⁶ views
1
HPREF_ELEMENT_TYPE ClassifyTet(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
2
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
3
		 INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint)
4
{
5
  int ep1(0), ep2(0), ep3(0), ep4(0), cp1(0), cp2(0), cp3(0), cp4(0), fp1, fp2, fp3, fp4;
6
  int isedge1(0), isedge2(0), isedge3(0), isedge4(0), isedge5(0), isedge6(0);
7
  int isfedge1, isfedge2, isfedge3, isfedge4, isfedge5, isfedge6;
8
  int isface1(0), isface2(0), isface3(0), isface4(0);
9

10
  HPREF_ELEMENT_TYPE type = HP_NONE; 
11
  
12

13
  int debug = 0;
14
  for (int j = 0;j < 4; j++)
15
    {
16
      if (el.pnums[j] == 444) debug++;
17
      if (el.pnums[j] == 115) debug++;
18
      if (el.pnums[j] == 382) debug++;
19
      if (el.pnums[j] == 281) debug++;
20
    }
21
  if (debug < 4) debug = 0;
22
  
23

24

25
  for (int j = 0; j < 4; j++)
26
    for (int k = 0; k < 4; k++)
27
      {
28
	if (j == k) continue;
29
	if (type) break;
30
	
31
	int pi3 = 0;
32
	while (pi3 == j || pi3 == k) pi3++;
33
	int pi4 = 6 - j - k - pi3;
34
	
35
	// preserve orientation
36
	int sort[4];
37
	sort[0] = j; sort[1] = k; sort[2] = pi3; sort[3] = pi4;
38
	int cnt = 0;
39
	for (int jj = 0; jj < 4; jj++)
40
	  for (int kk = 0; kk < 3; kk++)
41
	    if (sort[kk] > sort[kk+1])
42
	      {
43
		cnt++;
44
		Swap (sort[kk], sort[kk+1]); 
45
	      }
46
	if (cnt % 2 == 1) Swap (pi3, pi4);
47
	
48
	ep1 = edgepoint.Test (el.pnums[j]);
49
	ep2 = edgepoint.Test (el.pnums[k]);
50
	ep3 = edgepoint.Test (el.pnums[pi3]);
51
	ep4 = edgepoint.Test (el.pnums[pi4]);
52
	
53
	cp1 = cornerpoint.Test (el.pnums[j]);
54
	cp2 = cornerpoint.Test (el.pnums[k]);
55
	cp3 = cornerpoint.Test (el.pnums[pi3]);
56
	cp4 = cornerpoint.Test (el.pnums[pi4]);
57
	
58
	isedge1 = edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[k]));
59
	isedge2 = edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi3]));
60
	isedge3 = edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi4]));
61
	isedge4 = edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi3]));
62
	isedge5 = edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi4]));
63
	isedge6 = edges.Used (INDEX_2::Sort (el.pnums[pi3], el.pnums[pi4]));
64
	
65
	if (debug)
66
	  {
67
	    cout << "debug" << endl;
68
	    *testout  << "debug" << endl;
69
	    *testout << "ep = " << ep1 << ep2 << ep3 << ep4 << endl;
70
	    *testout << "cp = " << cp1 << cp2 << cp3 << cp4 << endl;
71
	    *testout << "edge = " << isedge1 << isedge2 << isedge3 << isedge4 << isedge5 << isedge6 << endl;
72
	  }
73

74

75
	isface1 = isface2 = isface3 = isface4 = 0;
76
	for (int l = 0; l < 4; l++)
77
	  {
78
	    INDEX_3 i3(0,0,0);
79
	    switch (l)
80
	      {
81
		case 0: i3.I1() = el.pnums[k]; i3.I1() = el.pnums[pi3]; i3.I1() = el.pnums[pi4]; break;
82
		case 1: i3.I1() = el.pnums[j]; i3.I1() = el.pnums[pi3]; i3.I1() = el.pnums[pi4]; break;
83
		case 2: i3.I1() = el.pnums[j]; i3.I1() = el.pnums[k]; i3.I1() = el.pnums[pi4]; break;
84
		case 3: i3.I1() = el.pnums[j]; i3.I1() = el.pnums[k]; i3.I1() = el.pnums[pi3]; break;
85
	      }
86
	    i3.Sort();
87
	    if (faces.Used (i3))
88
	      {
89
		int domnr = faces.Get(i3);
90
		if (domnr == -1 || domnr == el.GetIndex())
91
		  {
92
		    switch (l)
93
		      {
94
		      case 0: isface1 = 1; break;
95
		      case 1: isface2 = 1; break;
96
		      case 2: isface3 = 1; break;
97
		      case 3: isface4 = 1; break;
98
		      }
99
		  }
100
	      }
101
	  }
102
	/*
103
	  isface1 = faces.Used (INDEX_3::Sort (el.pnums[k], el.pnums[pi3], el.pnums[pi4]));
104
	  isface2 = faces.Used (INDEX_3::Sort (el.pnums[j], el.pnums[pi3], el.pnums[pi4]));
105
	  isface3 = faces.Used (INDEX_3::Sort (el.pnums[j], el.pnums[k], el.pnums[pi4]));
106
	  isface4 = faces.Used (INDEX_3::Sort (el.pnums[j], el.pnums[k], el.pnums[pi3]));
107
	*/
108
	
109
	isfedge1 = isfedge2 = isfedge3 = isfedge4 = isfedge5 = isfedge6 = 0;
110
	for (int l = 0; l < 6; l++)
111
	  {
112
	    INDEX_2 i2(0,0);
113
	    switch (l)
114
	      {
115
		case 0: i2.I1() = el.pnums[j]; i2.I2() = el[k]; break;
116
		case 1: i2.I1() = el.pnums[j]; i2.I2() = el.pnums[pi3]; break;
117
		case 2: i2.I1() = el.pnums[j]; i2.I2() = el.pnums[pi4]; break;
118
		case 3: i2.I1() = el.pnums[k]; i2.I2() = el.pnums[pi3]; break;
119
		case 4: i2.I1() = el.pnums[k]; i2.I2() = el.pnums[pi4]; break;
120
		case 5: i2.I1() = el.pnums[pi3]; i2.I2() = el.pnums[pi4]; break;
121
	      }
122
	    i2.Sort();
123
	    if (face_edges.Used (i2))
124
	      {
125
		int domnr = face_edges.Get(i2);
126
		if (domnr == -1 || domnr == el.GetIndex())
127
		  {
128
		    switch (l)
129
		      {
130
		      case 0: isfedge1 = 1; break;
131
		      case 1: isfedge2 = 1; break;
132
		      case 2: isfedge3 = 1; break;
133
		      case 3: isfedge4 = 1; break;
134
		      case 4: isfedge5 = 1; break;
135
		      case 5: isfedge6 = 1; break;
136
		      }
137
		  }
138
	      }
139
	  }
140
	/*
141
	  isfedge1 = face_edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[k]));
142
	  isfedge2 = face_edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi3]));
143
	  isfedge3 = face_edges.Used (INDEX_2::Sort (el.pnums[j], el.pnums[pi4]));
144
	  isfedge4 = face_edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi3]));
145
	  isfedge5 = face_edges.Used (INDEX_2::Sort (el.pnums[k], el.pnums[pi4]));
146
	  isfedge6 = face_edges.Used (INDEX_2::Sort (el.pnums[pi3], el.pnums[pi4]));
147
	*/
148
	
149
	fp1 = fp2 = fp3 = fp4 = 0;
150
	for (int l = 0; l < 4; l++)
151
	  {
152
	    int pti(0);
153
	    switch (l)
154
	      {
155
	      case 0: pti = el.pnums[j]; break;
156
	      case 1: pti = el.pnums[k]; break;
157
	      case 2: pti = el.pnums[pi3]; break;
158
	      case 3: pti = el.pnums[pi4]; break;
159
	      }
160
	    int domnr = facepoint[pti];
161
	    if (domnr == -1 || domnr == el.GetIndex())
162
	      {
163
		switch (l)
164
		  {
165
		  case 0: fp1 = 1; break;
166
		  case 1: fp2 = 1; break;
167
		  case 2: fp3 = 1; break;
168
		  case 3: fp4 = 1; break;
169
		  }
170
	      }
171
	  }
172
	
173
	/*
174
	  fp1 = facepoint[el.pnums[j]] != 0;
175
	  fp2 = facepoint[el.pnums[k]] != 0;
176
	  fp3 = facepoint[el.pnums[pi3]] != 0;
177
	  fp4 = facepoint[el.pnums[pi4]] != 0;
178
	*/
179
	
180
	
181
	switch (isface1+isface2+isface3+isface4)
182
	  {
183
	  case 0:
184
	    {
185
	      isedge1 |= isfedge1;
186
	      isedge2 |= isfedge2;
187
	      isedge3 |= isfedge3;
188
	      isedge4 |= isfedge4;
189
	      isedge5 |= isfedge5;
190
	      isedge6 |= isfedge6;
191
	      
192
	      ep1 |= fp1;
193
	      ep2 |= fp2;
194
	      ep3 |= fp3;
195
	      ep4 |= fp4;
196
	      
197
	      switch (isedge1+isedge2+isedge3+isedge4+isedge5+isedge6)
198
		{
199
		case 0:
200
		  {		
201
		    if (!ep1 && !ep2 && !ep3 && !ep4)
202
		      type = HP_TET;
203
				
204
		    if (ep1 && !ep2 && !ep3 && !ep4)
205
		      type = HP_TET_0E_1V;
206
		    
207
		    if (ep1 && ep2 && !ep3 && !ep4)
208
		      type = HP_TET_0E_2V;
209
		    
210
		    if (ep1 && ep2 && ep3 && !ep4)
211
		      type = HP_TET_0E_3V;
212
		    
213
		    if (ep1 && ep2 && ep3 && ep4)
214
		      type = HP_TET_0E_4V;
215
		    
216
		    break;
217
		  }
218
		  
219
		case 1:
220
		  {
221
		    if (!isedge1) break;
222
		    
223
		    if (!cp1 && !cp2 && !ep3 && !ep4)
224
		      type = HP_TET_1E_0V;
225
		    
226
		    if (cp1 && !cp2 && !ep3 && !ep4)
227
		      type = HP_TET_1E_1VA;
228
		    
229
		    if (!cp1 && !cp2 && !ep3 && ep4)
230
		      type = HP_TET_1E_1VB;
231
		    
232
		    if (cp1 && cp2 && !ep3 && !ep4)
233
		      type = HP_TET_1E_2VA;
234
		    
235
		    if (cp1 && !cp2 && ep3 && !ep4)
236
		      type = HP_TET_1E_2VB;
237
		    
238
		    if (cp1 && !cp2 && !ep3 && ep4)
239
		      type = HP_TET_1E_2VC;
240
		    
241
		    if (!cp1 && !cp2 && ep3 && ep4)
242
		      type = HP_TET_1E_2VD;
243
		    
244
		    if (cp1 && cp2 && ep3 && !ep4)
245
		      type = HP_TET_1E_3VA;
246
		    
247
		    if (cp1 && !cp2 && ep3 && ep4)
248
		      type = HP_TET_1E_3VB;
249
		    
250
		    if (cp1 && cp2 && ep3 && ep4)
251
		      type = HP_TET_1E_4V;
252
		    
253
		    break;
254
		  }
255
		case 2:
256
		  {
257
		    if (isedge1 && isedge2)
258
		      {
259
			if (!cp2 && !cp3 && !ep4)
260
			  type = HP_TET_2EA_0V;
261
			
262
			if (cp2 && !cp3 && !ep4)
263
			  type = HP_TET_2EA_1VA;
264
			if (!cp2 && cp3 && !ep4)
265
			  type = HP_TET_2EA_1VB;
266
			
267
			if (!cp2 && !cp3 && ep4)
268
			  type = HP_TET_2EA_1VC;
269
			
270
			if (cp2 && cp3 && !ep4)
271
			  type = HP_TET_2EA_2VA;
272
			if (cp2 && !cp3 && ep4)
273
			  type = HP_TET_2EA_2VB;
274
			if (!cp2 && cp3 && ep4)
275
			  type = HP_TET_2EA_2VC;
276
			
277
			if (cp2 && cp3 && ep4)
278
			  type = HP_TET_2EA_3V;
279
		      }
280
		    if (isedge1 && isedge6)
281
		      {
282
			if (!cp1 && !cp2 && !cp3 && !cp4)
283
			  type = HP_TET_2EB_0V;
284
			if (cp1 && !cp2 && !cp3 && !cp4)
285
			  type = HP_TET_2EB_1V;
286
			if (cp1 && cp2 && !cp3 && !cp4)
287
			  type = HP_TET_2EB_2VA;
288
			if (cp1 && !cp2 && cp3 && !cp4)
289
			  type = HP_TET_2EB_2VB;
290
			if (cp1 && !cp2 && !cp3 && cp4)
291
			  type = HP_TET_2EB_2VC;
292
			if (cp1 && cp2 && cp3 && !cp4)
293
			  type = HP_TET_2EB_3V;
294
			if (cp1 && cp2 && cp3 && cp4)
295
			  type = HP_TET_2EB_4V;
296
		      }
297
		    break;
298
		  }
299
		case 3:
300
		  {
301
		    if (isedge1 && isedge2 && isedge3)
302
		      {
303
			if (!cp2 && !cp3 && !cp4)
304
			  type = HP_TET_3EA_0V;
305
			if (cp2 && !cp3 && !cp4)
306
			  type = HP_TET_3EA_1V;
307
			if (cp2 && cp3 && !cp4)
308
			  type = HP_TET_3EA_2V;
309
			if (cp2 && cp3 && cp4)
310
			  type = HP_TET_3EA_3V;
311
		      }
312
		    if (isedge1 && isedge3 && isedge4)
313
		      {
314
			if (!cp3 && !cp4)
315
			  type = HP_TET_3EB_0V;
316
			if (cp3 && !cp4)
317
				      type = HP_TET_3EB_1V;
318
			if (cp3 && cp4)
319
			  type = HP_TET_3EB_2V;
320
		      }
321
		    if (isedge1 && isedge2 && isedge5)
322
		      {
323
			if (!cp3 && !cp4)
324
			  type = HP_TET_3EC_0V;
325
			if (cp3 && !cp4)
326
			  type = HP_TET_3EC_1V;
327
			if (cp3 && cp4)
328
			  type = HP_TET_3EC_2V;
329
		      }
330
		    break;
331
		  }
332
		}
333
	      break;
334
	    }
335
	    
336
	    
337
	    
338
	  case 1:  // one singular face
339
	    {
340
	      if (!isface1) break;
341
	      
342
	      switch (isfedge1+isfedge2+isfedge3+isedge4+isedge5+isedge6)
343
		{
344
		case 0:
345
		  {
346
		    if (!fp1 && !ep2 && !ep3 && !ep4)
347
		      type = HP_TET_1F_0E_0V;
348
		    if (fp1 && !ep2 && !ep3 && !ep4)
349
		      type = HP_TET_1F_0E_1VB;
350
		    if (!fp1 && ep2 && !ep3 & !ep4)
351
		      type = HP_TET_1F_0E_1VA;
352
		    break;
353
		  }
354
		case 1:
355
		  {
356
		    if (isfedge1)
357
		      {
358
			if (!ep1 && !ep3 && !ep4)
359
			  type = HP_TET_1F_1EA_0V;
360
		      }
361
		    if (isedge4) // V1-V3
362
		      {
363
			if (!ep1 && !cp2 && !cp3 && !ep4)
364
			  type = HP_TET_1F_1EB_0V;
365
		      }
366
		    break;
367
		  }
368
		}
369
	      break;
370
	    }
371
	    
372
	    
373
	  case 2:  // two singular faces
374
	    {
375
	      if (!isface1 || !isface2) break;
376
	      
377
	      switch (isfedge1+isedge2+isedge3+isedge4+isedge5)
378
		{
379
		case 0:
380
		  {
381
		    if (!ep1 && !ep2 && !cp3 && !cp4)
382
		      type = HP_TET_2F_0E_0V;
383
		    break;
384
		  }
385
		}
386
	      break;
387
	    }
388
	    
389
	    
390
	  }
391
	
392
	if (type != HP_NONE)
393
	  {
394
	    int pnums[4]; 
395
	    pnums[0] = el.pnums[j];
396
	    pnums[1] = el.pnums[k];
397
	    pnums[2] = el.pnums[pi3];
398
	    pnums[3] = el.pnums[pi4];
399
	    for(k=0;k<4;k++) el.pnums[k] = pnums[k]; 
400
	    break;
401
	  }
402
      }
403
  
404
  
405
  if (debug) cout << "type = " << type << endl;
406

407
  if (type == HP_NONE)
408
    {
409
      //     cnt_undef++;
410
      (*testout) << "undefined element" << endl
411
		 << "cp = " << cp1 << cp2 << cp3 << cp4 << endl
412
		 << "ep = " << ep1 << ep2 << ep3 << ep4 << endl
413
		 << "isedge = " << isedge1 << isedge2 << isedge3 
414
		 << isedge4 << isedge5 << isedge6 << endl
415
		 << "isface = " << isface1 << isface2 << isface3 << isface4 << endl;
416
      cout << "undefined element !!! " << endl;
417

418
      
419
    }
420
  return(type); 
421
}
422

423

424

425
HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
426
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
427
		 INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint)
428
{
429

430
  HPREF_ELEMENT_TYPE type = HP_NONE;
431
  
432
  int p[6];
433
  for(int m=1;m<=6;m++)
434
    {
435
      int point_sing[6]={0,0,0,0,0,0}; 
436
      int face_sing[5]={0,0,0,0,0};
437
      int edge_sing[9]={0,0,0,0,0,0,0,0,0}; 
438
      
439
      if(m<4)
440
	{ 
441
	  p[0]= m; p[1]=m%3+1; p[2]=(m%3+1)%3+1;
442
	  for(int l=3;l<6;l++) p[l]=p[l-3]+3;  
443
	}
444
      else
445
	{
446
	  p[0] = m; p[1]=(m%3+1)%3+4; p[2]=m%3+4;
447
	  for(int l=3;l<6;l++) p[l]=p[l-3]-3; 
448
	}
449
      
450
      for(int j=0;j<6;j++) 
451
	{ 
452
	  if(cornerpoint.Test(el.PNum(p[j])))  { point_sing[p[j]-1]=3;}
453
	  else if(edgepoint.Test(el.PNum(p[j]))) point_sing[p[j]-1]=2;
454
	  else if (facepoint[el.PNum(p[j])] == -1 || facepoint[el.PNum(p[j])] == el.GetIndex())
455
	    point_sing[p[j]-1] = 1;  
456
	}
457
      
458
      const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (PRISM);
459
      for(int k=0;k<9;k++)
460
	{
461
	  INDEX_2 i2 = INDEX_2 :: Sort(el.PNum(p[eledges[k][0]-1]),el.PNum(p[eledges[k][1]-1])); 
462
	  if (edges.Used(i2)) edge_sing[k] = 2;
463
	  else edge_sing[k] = face_edges.Used(i2);
464
	}
465
      
466
      const ELEMENT_FACE * elfaces  = MeshTopology::GetFaces (PRISM);
467
      for (int k=0;k<5;k++)
468
	{
469
	  INDEX_3 i3; 
470
	  
471
	  if(k<2) 
472
	    i3 = INDEX_3::Sort(el.pnums[p[elfaces[k][0]-1]-1], el.pnums[p[elfaces[k][1]-1]-1], 
473
			       el.pnums[p[elfaces[k][2]-1]-1]); 
474
	  else 
475
	    { 
476
	      INDEX_4  i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]-1], el.pnums[p[elfaces[k][1]-1]-1], el.pnums[p[elfaces[k][2]-1]-1],el.pnums[p[elfaces[k][3]-1]-1]); 
477
	      i4.Sort();
478
	      i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3()); 
479
	    }
480
	  
481
	  if (faces.Used (i3))
482
	    {
483
	      int domnr = faces.Get(i3); 
484
	      if (domnr == -1 || domnr == el.GetIndex())
485
		face_sing[k] = 1; 
486
	      
487
	    } 
488
	} 
489
      if (face_sing[1] > face_sing[0]) {m=m+2; continue;}  
490
      
491
      
492
      //int cp = 0;  
493
      
494
      int qfsing = face_sing[2] + face_sing[3] + face_sing[4];
495
      int tfsing = face_sing[0] + face_sing[1]; 
496
      int evsing = edge_sing[6] + edge_sing[7] + edge_sing[8];
497
      int ehsing = edge_sing[0] + edge_sing[1] + edge_sing[2] + edge_sing[3] + edge_sing[4] + edge_sing[5];
498
      
499
      if (qfsing + tfsing + evsing + ehsing == 0)  
500
	{ type = HP_PRISM;  break;}
501
      
502
      HPREF_ELEMENT_TYPE types[] = {HP_NONE,HP_NONE,HP_NONE};   
503
      
504
      int fb = (1-face_sing[4])* face_sing[3] * (face_sing[2] + face_sing[3]) + 3*face_sing[4]*face_sing[3]*face_sing[2];  
505
      int sve[3] = {edge_sing[7] , edge_sing[8], edge_sing[6]}; 
506
      
507
            
508
      if(fb!=qfsing) continue; 
509
      
510
      
511
      switch(fb)
512
	{ 
513
	case 0: 
514
	  if (evsing == 0 && ehsing==3*tfsing) 
515
	    {
516
	      types[0] = HP_PRISM; 
517
	      types[1] = HP_PRISM_1FA_0E_0V;   
518
	      types[2] = HP_PRISM_2FA_0E_0V; 
519
	    } 
520
	  if(evsing > 0 &&  sve[0] == evsing) // 1 vertical edge 1-4 
521
	    { 
522
	      types[0] = HP_PRISM_SINGEDGE;
523
	      types[1] = HP_PRISM_1FA_1E_0V;
524
	      types[2] = HP_PRISM_2FA_1E_0V;   
525
	    }
526
	  
527
	  if(sve[0] > 0 && sve[1] > 0 && sve[2] == 0)
528
	    {
529
	      types[0] = HP_PRISM_SINGEDGE_V12;
530
	      types[1] = HP_PRISM_1FA_2E_0V; 
531
	      types[2] = HP_PRISM_2FA_2E_0V; 
532
	    }
533
	  if(sve[0] > 0 && sve[1] > 0 && sve[2] > 0) 
534
	    {
535
	      types[0] = HP_PRISM_3E_0V;
536
	      types[1] = HP_PRISM_1FA_3E_0V;
537
	      types[2] = HP_PRISM_2FA_3E_0V;
538
	      
539
	      if ( edge_sing[0] > 1 && edge_sing[2] > 1 &&  
540
		   edge_sing[4] > 1 && edge_sing[5] > 1 && tfsing==0)
541
		types[0] = HP_PRISM_3E_4EH; 
542
	    }
543
	  
544
	  break;
545
	case 1:
546
	  if(sve[0] <= 1 && sve[1] <= 1)  
547
	    if(sve[2]==0)
548
	      { 
549
		types[0] = HP_PRISM_1FB_0E_0V;
550
		types[1] = HP_PRISM_1FA_1FB_0E_0V;
551
		types[2] = HP_PRISM_2FA_1FB_0E_0V;  
552
	      }
553
	    else
554
	      { 
555
		types[0] = HP_PRISM_1FB_1EC_0V;
556
		types[1] = HP_PRISM_1FA_1FB_1EC_0V;
557
		types[2] = HP_PRISM_2FA_1FB_1EC_0V; 
558
	      }
559
	  
560
	  if(sve[0] > 1 && sve[2] >= 1 && sve[1] <= 1)
561
	    { 
562
	      types[0] = HP_PRISM_1FB_2EB_0V;  
563
	      types[1] = HP_PRISM_1FA_1FB_2EB_0V;
564
	      types[2] = HP_PRISM_2FA_1FB_2EB_0V; 
565
	    }
566
	  
567
	  if(sve[0] > 1 && sve[1] <= 1 && sve[2] == 0) // ea && !eb  
568
	    {
569
	      types[0] = HP_PRISM_1FB_1EA_0V;
570
	      types[1] = HP_PRISM_1FA_1FB_1EA_0V;
571
	      types[2] = HP_PRISM_2FA_1FB_1EA_0V; 
572
	    } 
573
	  
574
	  if(sve[0] <= 1 && sve[1] > 1 && sve[2] == 0)
575
	    types[1] = HP_PRISM_1FA_1FB_1EB_0V; 
576
	  
577
	  if(sve[0] > 1 && sve[1]>1) 
578
	    if(sve[2] == 0)  // ea && eb 
579
	      {
580
		types[0] = HP_PRISM_1FB_2EA_0V;
581
		types[1] = HP_PRISM_1FA_1FB_2EA_0V;
582
		types[2] = HP_PRISM_2FA_1FB_2EA_0V; 
583
	      }
584
	  if(sve[0] <= 1 && sve[1] > 1 && sve[2] >0)
585
	    types[1] = HP_PRISM_1FA_1FB_2EC_0V; 
586
	  
587
	  if(sve[0] > 1 && sve[1] > 1 && sve[2] >= 1) //sve[2] can also be a face-edge  
588
	    {
589
	      types[0] = HP_PRISM_1FB_3E_0V;  
590
	      types[1] = HP_PRISM_1FA_1FB_3E_0V; 
591
	      types[2] = HP_PRISM_2FA_1FB_3E_0V; 
592
	    } 
593
	  
594
	  break;  
595
	  
596
	case 2:
597
	  if(sve[0] <= 1) 
598
	    cout << " **** WARNING: Edge between to different singular faces should be marked singular " << endl; 
599
		      
600
	  if(sve[1] <= 1)   
601
	    if(sve[2] <=1) 
602
	      { 
603
		types[0] = HP_PRISM_2FB_0E_0V; 
604
		types[1] = HP_PRISM_1FA_2FB_0E_0V;
605
		types[2] = HP_PRISM_2FA_2FB_0E_0V;
606
	      }
607
	    else
608
	      { 
609
		types[0] = HP_PRISM_2FB_1EC_0V; 
610
		types[1] = HP_PRISM_1FA_2FB_1EC_0V; 
611
		types[2] = HP_PRISM_2FA_2FB_1EC_0V;   
612
	      }
613
	  else
614
	    if(sve[2] <= 1) 
615
	      types[1] = HP_PRISM_1FA_2FB_1EB_0V; 
616
	    else
617
	      { 
618
		types[0] = HP_PRISM_2FB_3E_0V; 
619
		types[1] = HP_PRISM_1FA_2FB_3E_0V; 
620
		types[2] = HP_PRISM_2FA_2FB_3E_0V; 
621
	      }
622
	  
623
	  break;
624
	  
625
	case 3: 
626
	  types[0] = HP_PRISM_3FB_0V; 
627
	  types[1] = HP_PRISM_1FA_3FB_0V; 
628
	  types[2] = HP_PRISM_2FA_3FB_0V; 
629
	  break;
630
	}
631
      type = types[tfsing];
632
      
633
         
634
      if(type != HP_NONE)  
635
	break;
636
    }
637
	 
638
  /*
639
   *testout << " Prism with pnums " << endl; 
640
   for(int j=0;j<6;j++) *testout << el.pnums[j] << "\t"; 
641
   *testout << endl; 
642
   */
643
  
644
  if(type != HP_NONE) 
645
    {
646
      int pnums[6]; 
647
      for(int j=0;j<6;j++) pnums[j] = el.PNum (p[j]);
648
      for(int k=0;k<6;k++) el.pnums[k] = pnums[k]; 
649
    }
650

651
  /* *testout << " Classified Prism with pnums " << endl; 
652
  for(int j=0;j<6;j++) *testout << el.pnums[j] << "\t"; 
653
  *testout << endl; 
654
  */ 
655
  return(type); 
656
}
657

658

659
// #ifdef SABINE 
660
HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
661
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int>  & face_edges, 
662
				INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
663

664
{
665
  HPREF_ELEMENT_TYPE type = HP_NONE; 
666
  
667
  int pnums[3]; 
668
  int p[3];   
669
  
670
  INDEX_3 i3 (el.pnums[0], el.pnums[1], el.pnums[2]);
671
  i3.Sort();
672
  bool sing_face = faces.Used (i3);
673
  
674
  // *testout << " facepoint " << facepoint << endl;  
675
      
676

677
  // Try all rotations of the trig 
678
  for (int j=0;j<3;j++) 
679
    {
680
      int point_sing[3] = {0,0,0}; 
681
      int edge_sing[3] = {0,0,0}; 
682
      // *testout << " actual rotation of trig points " ;  
683
      for(int m=0;m<3;m++) 
684
	{ 
685
	  p[m] = (j+m)%3 +1; // local vertex number
686
	  pnums[m] = el.PNum(p[m]); // global vertex number 
687
	  // *testout << pnums[m] << " \t "; 
688
	}
689
      // *testout << endl ; 
690
      
691
      if(dim == 3) 
692
	{
693
	  // face point 
694
	  for(int k=0;k<3;k++)
695
	    if(!sing_face)
696
	      { 
697
		//	*testout << " fp [" << k << "] = " << facepoint[pnums[k]] << endl;   
698
		//	*testout << " fd.DomainIn()" <<  fd.DomainIn() << endl; 
699
		//	*testout  << " fd.DomainOut()" <<  fd.DomainOut() << endl; 
700
		if( facepoint[pnums[k]]  && (facepoint[pnums[k]] ==-1 || 
701
					     facepoint[pnums[k]] == fd.DomainIn() ||   facepoint[pnums[k]] == fd.DomainOut()))
702
		  point_sing[p[k]-1] = 1; 
703
	      } 
704
	  // if point is on face_edge in next step sing = 2 
705

706
	  /*	  *testout << " pointsing NACH FACEPOints ... FALLS EDGEPOINT UMSETZEN" ; 
707
	  for (int k=0;k<3;k++) *testout << "\t" << point_sing[p[k]-1] ;
708
	  *testout << endl; */
709
	}
710
      
711
      const ELEMENT_EDGE * eledges = MeshTopology::GetEdges(TRIG); 
712
      
713
      if(dim==3)
714
	{
715
	  for(int k=0;k<3;k++) 
716
	    { 
717
	      int ep1=p[eledges[k][0]-1];  
718
	      int ep2=p[eledges[k][1]-1];  
719
	      INDEX_2 i2(el.PNum(ep1),el.PNum(ep2)); 
720
	      
721
	      if(edges.Used(i2)) 
722
		{
723
		  
724
		  edge_sing[k]=2;
725
		  point_sing[ep1-1] = 2; 
726
		  point_sing[ep2-1] = 2; 
727
		}
728
	      else // face_edge? 
729
		{	  
730
		  i2.Sort();  
731
		  if(surf_edges.Used(i2) && surf_edges.Get(i2) != fd.SurfNr()+1)  // edge not face_edge acc. to surface in which trig lies
732
		    if(face_edges.Get(i2)==-1 ||face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut() )
733
		      { 
734
			edge_sing[k]=1;
735
		      } 
736
		    else
737
		      { 
738
			point_sing[ep1-1] = 0; // set to edge_point 
739
			point_sing[ep2-1] = 0; // set to edge_point
740
		      } 
741
		}
742
	      
743
	      /*  *testout << " pointsing NACH edges UND FACEEDGES UMSETZEN ... " ; 
744
	      for (int k=0;k<3;k++) *testout << "\t" << point_sing[p[k]-1] ; 
745
	      *testout << endl;          
746
	      */
747
	    }
748
	}
749
      /*
750
      *testout << " dim " << dim << endl; 
751
      *testout << " edgepoint_dom " << edgepoint_dom << endl; 
752
      */
753
      if(dim==2)
754
	{
755
	  for(int k=0;k<3;k++) 
756
	    { 
757
	      int ep1=p[eledges[k][0]-1];  
758
	      int ep2=p[eledges[k][1]-1];  
759
	     
760
	      INDEX_2 i2(el.PNum(ep1),el.PNum(ep2));  
761
	     
762
	      if(edges.Used(i2)) 
763
		{
764
		  
765
		  if(edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep1-1])) || 
766
		     edgepoint_dom.Used(INDEX_2(-1,pnums[ep1-1])) || 
767
		     edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep2-1])) || 
768
		     edgepoint_dom.Used(INDEX_2(-1,pnums[ep2-1]))) 
769
		    {
770
		      edge_sing[k]=2;
771
		      point_sing[ep1-1] = 2;
772
		      point_sing[ep2-1] = 2; 
773
		    }
774
		}
775
	     
776
	    }
777
	}
778

779
     
780
	 
781
      for(int k=0;k<3;k++) 
782
	if(edgepoint.Test(pnums[k])) //edgepoint, but not member of sing_edge on trig -> cp 
783
	  {
784
	    INDEX_2 i2a=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+1)%3])); 
785
	    INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3])); 
786
	    
787
	    if(!edges.Used(i2a) && !edges.Used(i2b)) 
788
	      point_sing[p[k]-1] = 3; 	
789
	  } 
790
      
791
      for(int k=0;k<3;k++) 
792
	if(cornerpoint.Test(el.PNum(p[k]))) 
793
	  point_sing[p[k]-1] = 3; 
794
      
795

796
     if(edge_sing[0] + edge_sing[1] + edge_sing[2] == 0) 
797
       { 
798
	 int ps = point_sing[0] + point_sing[1] + point_sing[2]; 
799
	 
800
	 if(ps==0) 
801
	   type = HP_TRIG; 
802
	 else if(point_sing[p[0]-1]  && !point_sing[p[1]-1] && !point_sing[p[2]-1])
803
	   type = HP_TRIG_SINGCORNER;
804
	 else if(point_sing[p[0]-1] && point_sing[p[1]-1] && !point_sing[p[2]-1]) 
805
	   type = HP_TRIG_SINGCORNER12; 
806
	 else if(point_sing[p[0]-1] && point_sing[p[1]-1] && point_sing[p[2]-1]) 
807
	   { 
808
	     if(dim==2) type = HP_TRIG_SINGCORNER123_2D; 
809
	     else type = HP_TRIG_SINGCORNER123; 
810
	   } 
811
       } 
812
     else
813
       if (edge_sing[2] && !edge_sing[0] && !edge_sing[1]) //E[2]=(1,2) 
814
      { 
815
	int code = 0; 
816
	if(point_sing[p[0]-1] > edge_sing[2]) code+=1; 
817
	if(point_sing[p[1]-1] > edge_sing[2]) code+=2; 
818
	if(point_sing[p[2]-1]) code+=4; 
819
	
820
	HPREF_ELEMENT_TYPE types[] =
821
	  {
822
	    HP_TRIG_SINGEDGE, 
823
	    HP_TRIG_SINGEDGECORNER1, 
824
	    HP_TRIG_SINGEDGECORNER2,
825
	    HP_TRIG_SINGEDGECORNER12, 
826
	    HP_TRIG_SINGEDGECORNER3, 
827
	    HP_TRIG_SINGEDGECORNER13, 
828
	    HP_TRIG_SINGEDGECORNER23, 
829
	    HP_TRIG_SINGEDGECORNER123, 
830
	  };
831
	type = types[code]; 
832
	
833
      }  // E[0] = [0,2], E[1] =[1,2], E[2] = [0,1]
834
    else 
835
      if(edge_sing[2] && !edge_sing[1] && edge_sing[0])
836
	{
837
	  if(point_sing[p[2]-1] <= edge_sing[0] ) 
838
	    { 
839
	      if(point_sing[p[1]-1]<= edge_sing[2]) type = HP_TRIG_SINGEDGES; 
840
	      else type = HP_TRIG_SINGEDGES2; 
841
	    } 
842
	  else 
843
	    {
844
	      if(point_sing[p[1]-1]<= edge_sing[2]) 
845
		type = HP_TRIG_SINGEDGES3; 
846
	      else type = HP_TRIG_SINGEDGES23; 
847
	    }
848
	}
849
      else if (edge_sing[2] && edge_sing[1] && edge_sing[0])
850
	type = HP_TRIG_3SINGEDGES; 
851
     
852
     //  cout << " run for " <<  j << " gives type " << type << endl; 
853
     //*testout << " run for " <<  j << " gives type " << type << endl; 
854
     if(type!=HP_NONE) break;
855
    }
856
    
857
  for(int k=0;k<3;k++) el[k] = pnums[k]; 
858
  /*if(type != HP_NONE) 
859
   {
860
     
861
      cout << " TRIG with pnums " << pnums[0] << "\t"  << 
862
	pnums[1] << "\t"  << pnums[2] << endl; 
863
      cout << " type "  << type << endl; 
864
      }
865
  */
866
  return(type);
867
}
868
#ifdef HPREF_OLD 
869
HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
870
				BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
871
				INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
872
{
873
  HPREF_ELEMENT_TYPE type = HP_NONE; 
874
  
875
  int pnums[3]; 
876
	      
877
  INDEX_3 i3 (el.pnums[0], el.pnums[1], el.pnums[2]);
878
  i3.Sort();
879
  bool sing_face = faces.Used (i3);
880
   
881
  
882
  for (int j = 1; j <= 3; j++)
883
    {
884
      int ep1 = edgepoint.Test (el.PNumMod (j));
885
      int ep2 = edgepoint.Test (el.PNumMod (j+1));
886
      int ep3 = edgepoint.Test (el.PNumMod (j+2));
887
      
888
      if (dim == 2)
889
	{
890
	  // JS, Dec 11
891
	  ep1 = edgepoint_dom.Used (INDEX_2 (fd.SurfNr(), el.PNumMod(j))) ||
892
	    edgepoint_dom.Used (INDEX_2 (-1, el.PNumMod(j)));
893
	  ep2 = edgepoint_dom.Used (INDEX_2 (fd.SurfNr(), el.PNumMod(j+1))) ||
894
	    edgepoint_dom.Used (INDEX_2 (-1, el.PNumMod(j+1)));
895
	  ep3 = edgepoint_dom.Used (INDEX_2 (fd.SurfNr(), el.PNumMod(j+2))) ||
896
	    edgepoint_dom.Used (INDEX_2 (-1, el.PNumMod(j+2)));
897
	  /*
898
	  ep1 = edgepoint_dom.Used (INDEX_2 (el.index, el.PNumMod(j)));
899
	  ep2 = edgepoint_dom.Used (INDEX_2 (el.index, el.PNumMod(j+1)));
900
	  ep3 = edgepoint_dom.Used (INDEX_2 (el.index, el.PNumMod(j+2)));
901
	  */
902
	  // ep3 = edgepoint_dom.Used (INDEX_2 (mesh.SurfaceElement(i).GetIndex(), el.PNumMod(j+2)));
903
	}
904
      
905
      
906
      
907
      int cp1 = cornerpoint.Test (el.PNumMod (j));
908
      int cp2 = cornerpoint.Test (el.PNumMod (j+1));
909
      int cp3 = cornerpoint.Test (el.PNumMod (j+2));
910
      
911
      ep1 |= cp1;
912
      ep2 |= cp2;
913
      ep3 |= cp3;
914
      
915
      
916
      // (*testout) << "cp = " << cp1 << cp2 << cp3 << ", ep = " << ep1 << ep2 << ep3 << endl;
917

918
      int p[3] = { el.PNumMod (j), el.PNumMod (j+1), el.PNumMod (j+2)};
919
      if(ep1)
920
	{
921
	  INDEX_2 i2a=INDEX_2::Sort(p[0], p[1]); 
922
	  INDEX_2 i2b=INDEX_2::Sort(p[0], p[2]); 
923
	  if(!edges.Used(i2a) && !edges.Used(i2b)) 
924
	    cp1 = 1; 
925
	}
926
      if(ep2)
927
	{
928
	  INDEX_2 i2a=INDEX_2::Sort(p[1], p[0]); 
929
	  INDEX_2 i2b=INDEX_2::Sort(p[1], p[2]); 
930
	  if(!edges.Used(i2a) && !edges.Used(i2b)) 
931
	    cp2 = 1; 
932
	}
933
      if(ep3)
934
	{
935
	  INDEX_2 i2a=INDEX_2::Sort(p[2], p[0]); 
936
	  INDEX_2 i2b=INDEX_2::Sort(p[2], p[1]); 
937
	  if(!edges.Used(i2a) && !edges.Used(i2b)) 
938
	    cp3= 1; 
939
	}		      
940
      
941
      
942
      int isedge1=0, isedge2=0, isedge3=0; 
943
      if(dim == 3 )
944
	{
945
	  INDEX_2 i2;
946
	  i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
947
	  isedge1 = edges.Used (i2);
948
	  i2.Sort();
949
	  if(surf_edges.Used(i2) &&  surf_edges.Get(i2)   != fd.SurfNr()+1 && 
950
	     (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
951
	    {
952
	      isedge1=1;
953
	      ep1 = 1; ep2=1;
954
	    }
955
	  
956
	  i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
957
	  isedge2 = edges.Used (i2);
958
	  i2.Sort();
959
	  if(surf_edges.Used(i2) &&  surf_edges.Get(i2)   != fd.SurfNr()+1 &&
960
	     (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
961
	    {
962
	      isedge2=1;
963
	      ep2 = 1; ep3=1;
964
	    }
965
	  i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
966
	  isedge3 = edges.Used (i2);
967
	  i2.Sort();
968
	  if(surf_edges.Used(i2) &&  surf_edges.Get(i2)   != fd.SurfNr()+1 && 
969
	     (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
970
	    {
971
	      isedge3=1;
972
	      ep1 = 1; ep3=1;
973
	    }
974
	  
975
	  // cout << " isedge " << isedge1 << " \t " << isedge2 << " \t " << isedge3 << endl;  
976
	
977
	  if (!sing_face)
978
	{
979
	  /*
980
	    if (!isedge1)  { cp1 |= ep1; cp2 |= ep2; }
981
	    if (!isedge2)  { cp2 |= ep2; cp3 |= ep3; }
982
	    if (!isedge3)  { cp3 |= ep3; cp1 |= ep1; }
983
	  */
984
	  ep1 |= facepoint [el.PNumMod(j)] != 0;
985
	  ep2 |= facepoint [el.PNumMod(j+1)] != 0;
986
	  ep3 |= facepoint [el.PNumMod(j+2)] != 0;
987
	  
988
	  
989
	  isedge1 |= face_edges.Used (INDEX_2::Sort (el.PNumMod(j), el.PNumMod(j+1)));
990
	  isedge2 |= face_edges.Used (INDEX_2::Sort (el.PNumMod(j+1), el.PNumMod(j+2)));
991
	  isedge3 |= face_edges.Used (INDEX_2::Sort (el.PNumMod(j+2), el.PNumMod(j+3)));
992
	}
993
	}
994
      
995
      if(dim ==2) 
996
	{ 
997
	  INDEX_2 i2;
998
	  i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
999
	  i2.Sort();
1000
	  isedge1 = edges.Used (i2);
1001
	  if(isedge1)
1002
	    {
1003
	      ep1 = 1; ep2=1;
1004
	    }
1005
	  
1006
	  i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
1007
	  i2.Sort();
1008
	  isedge2 = edges.Used (i2);
1009
	  if(isedge2)
1010
	    {
1011
	      ep2 = 1; ep3=1;
1012
	    }
1013
	  i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
1014
	  i2.Sort();
1015
	  isedge3 = edges.Used (i2);
1016
	  if(isedge3)
1017
	    {
1018
	      ep1 = 1; ep3=1;
1019
	    }
1020
	  
1021
	  
1022
	}
1023
      
1024
		  
1025
		  /*
1026
		  cout << " used " << face_edges.Used (INDEX_2::Sort (el.PNumMod(j), el.PNumMod(j+1))) << endl; 
1027

1028
		  cout << " isedge " << isedge1 << " \t " << isedge2 << " \t " << isedge3 << endl; 
1029
		  cout << " ep " << ep1 << "\t" << ep2 << " \t " << ep3 << endl; 
1030
		  cout << " cp " << cp1 << "\t" << cp2 << " \t " << cp3 << endl; 
1031
		  */
1032
		  
1033

1034
      
1035
      if (isedge1 + isedge2 + isedge3 == 0)
1036
	{
1037
	  if (!ep1 && !ep2 && !ep3)
1038
	    type = HP_TRIG;
1039
	  
1040
	  if (ep1 && !ep2 && !ep3)
1041
	    type = HP_TRIG_SINGCORNER;
1042
	  
1043
	  if (ep1 && ep2 && !ep3)
1044
	    type = HP_TRIG_SINGCORNER12;
1045
	  
1046
	  if (ep1 && ep2 && ep3)
1047
	    {
1048
	      if (dim == 2)
1049
			    type = HP_TRIG_SINGCORNER123_2D;
1050
	      else
1051
		type = HP_TRIG_SINGCORNER123;
1052
	    }
1053
	  
1054
	  if (type != HP_NONE)
1055
	    {
1056
	      pnums[0] = el.PNumMod (j);
1057
	      pnums[1] = el.PNumMod (j+1);
1058
	      pnums[2] = el.PNumMod (j+2);
1059
	      break;
1060
	    }
1061
	}
1062
      
1063
      if (isedge1 && !isedge2 && !isedge3)
1064
	{
1065
	  int code = 0;
1066
	  if (cp1) code += 1;
1067
	  if (cp2) code += 2;
1068
	  if (ep3) code += 4;
1069
	  
1070
	  HPREF_ELEMENT_TYPE types[] =
1071
	    {
1072
	      HP_TRIG_SINGEDGE, 
1073
	      HP_TRIG_SINGEDGECORNER1, 
1074
	      HP_TRIG_SINGEDGECORNER2,
1075
	      HP_TRIG_SINGEDGECORNER12, 
1076
	      HP_TRIG_SINGEDGECORNER3, 
1077
	      HP_TRIG_SINGEDGECORNER13, 
1078
	      HP_TRIG_SINGEDGECORNER23, 
1079
	      HP_TRIG_SINGEDGECORNER123, 
1080
	    };
1081
	  type = types[code];
1082
	  pnums[0] = el.PNumMod (j);
1083
	  pnums[1] = el.PNumMod (j+1);
1084
	  pnums[2] = el.PNumMod (j+2);
1085
	  break;
1086
	}
1087
      
1088
      
1089
      if (isedge1 && !isedge2 && isedge3)
1090
	{
1091
	  if (!cp3)
1092
	    {
1093
	      if (!cp2) type = HP_TRIG_SINGEDGES;
1094
	      else      type = HP_TRIG_SINGEDGES2;
1095
	    }
1096
	  else
1097
	    { 
1098
	      if (!cp2) type = HP_TRIG_SINGEDGES3;
1099
	      else      type = HP_TRIG_SINGEDGES23;
1100
	    }
1101
	  
1102
	  pnums[0] = el.PNumMod (j);
1103
	  pnums[1] = el.PNumMod (j+1);
1104
	  pnums[2] = el.PNumMod (j+2);
1105
	  break;
1106
	}
1107
       
1108
      if (isedge1 && isedge2 && isedge3)
1109
	{
1110
	  type = HP_TRIG_3SINGEDGES;
1111
	  pnums[0] = el.PNumMod (j);
1112
	  pnums[1] = el.PNumMod (j+1);
1113
	  pnums[2] = el.PNumMod (j+2);
1114
	  break;
1115
	}
1116
    }
1117
  
1118
  for(int k=0;k<3;k++) el[k] = pnums[k]; 
1119
  /*if(type != HP_NONE) 
1120
   {
1121
     
1122
      cout << " TRIG with pnums " << pnums[0] << "\t"  << 
1123
	pnums[1] << "\t"  << pnums[2] << endl; 
1124
      cout << " type "  << type << endl; 
1125
      }
1126
  */
1127
  return(type);
1128
}
1129
#endif
1130
HPREF_ELEMENT_TYPE ClassifyQuad(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
1131
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int>  & face_edges, 
1132
				INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint, int dim, const FaceDescriptor & fd)
1133
{
1134
  HPREF_ELEMENT_TYPE type = HP_NONE; 
1135
  
1136
  int ep1(-1), ep2(-1), ep3(-1), ep4(-1), cp1(-1), cp2(-1), cp3(-1), cp4(-1);
1137
	      int isedge1, isedge2, isedge3, isedge4;
1138

1139
	      for (int j = 1; j <= 4; j++)
1140
		{
1141
		  ep1 = edgepoint.Test (el.PNumMod (j));
1142
		  ep2 = edgepoint.Test (el.PNumMod (j+1));
1143
		  ep3 = edgepoint.Test (el.PNumMod (j+2));
1144
		  ep4 = edgepoint.Test (el.PNumMod (j+3));
1145

1146
		  if (dim == 2)
1147
		    {
1148
		      ep1 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j)));
1149
		      ep2 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j+1)));
1150
		      ep3 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j+2)));
1151
		      ep4 = edgepoint_dom.Used (INDEX_2 (el.GetIndex(), el.PNumMod(j+3)));
1152
		    }
1153

1154
		  cp1 = cornerpoint.Test (el.PNumMod (j));
1155
		  cp2 = cornerpoint.Test (el.PNumMod (j+1));
1156
		  cp3 = cornerpoint.Test (el.PNumMod (j+2));
1157
		  cp4 = cornerpoint.Test (el.PNumMod (j+3));
1158

1159
		  ep1 |= cp1;
1160
		  ep2 |= cp2;
1161
		  ep3 |= cp3;
1162
		  ep4 |= cp4;
1163
		
1164
		  int p[4] = { el.PNumMod (j), el.PNumMod (j+1), el.PNumMod (j+2), el.PNumMod(j+4)};
1165
		  //int epp[4] = { ep1, ep2, ep3, ep4}; 
1166
		  int cpp[4] = { cp1, cp2, cp3, cp4};
1167
		  for(int k=0;k<0;k++)
1168
		    {
1169
		      INDEX_2 i2a=INDEX_2::Sort(p[k], p[(k+1)%4]); 
1170
		      INDEX_2 i2b=INDEX_2::Sort(p[k], p[(k-1)%4]); 
1171
		      if(!edges.Used(i2a) && !edges.Used(i2b)) 
1172
			cpp[k] = 1; 
1173
		    }
1174
		  cp1= cpp[0]; cp2=cpp[1]; cp3=cpp[2]; cp4=cpp[3];
1175
		  
1176

1177
		  if(dim ==3) 
1178
		    { 
1179
		      INDEX_2 i2;
1180
		      i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
1181
		      // i2.Sort();
1182
		      isedge1 = edges.Used (i2);
1183
		      i2.Sort();
1184
		      if(surf_edges.Used(i2)  &&  surf_edges.Get(i2)   != fd.SurfNr()+1 && 
1185
			 (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
1186
			{
1187
			  isedge1=1;
1188
			  ep1 = 1; ep2=1;
1189
			}
1190
		      i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
1191
		      // i2.Sort();
1192
		      isedge2 = edges.Used (i2);
1193
		      i2.Sort();
1194
		      if(surf_edges.Used(i2)  &&  surf_edges.Get(i2)   != fd.SurfNr()+1 && 
1195
			 (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
1196
			{  
1197
			  isedge2=1;
1198
			  ep2=1; ep3=1;
1199
			}
1200
		      i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
1201
		      // i2.Sort();
1202
		      isedge3 = edges.Used (i2); 
1203
		      i2.Sort();
1204
		      if(surf_edges.Used(i2)   &&  surf_edges.Get(i2)   != fd.SurfNr()+1 && (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
1205
			{
1206
			  isedge3=1;
1207
			  ep3=1; ep4=1;
1208
			}
1209
		      i2 = INDEX_2(el.PNumMod (j+3), el.PNumMod (j+4));
1210
		      // i2.Sort();
1211
		      isedge4 = edges.Used (i2);
1212
		      i2.Sort();
1213
		      if(surf_edges.Used(i2)  &&  surf_edges.Get(i2)   != fd.SurfNr()+1 && 
1214
			 (face_edges.Get(i2) == -1 || face_edges.Get(i2) == fd.DomainIn() || face_edges.Get(i2) == fd.DomainOut()) ) 
1215
			{ 
1216
			  isedge4=1;
1217
			  ep4=1; ep1=1;
1218
			} 
1219
		    
1220

1221
//MH***********************************************************************************************************
1222
		      if(ep1)
1223
			if(edgepoint.Test(p[0]))
1224
			{
1225
			  INDEX_2 i2a=INDEX_2::Sort(p[0], p[1]); 
1226
	  		  INDEX_2 i2b=INDEX_2::Sort(p[0], p[3]); 
1227
	  		  if(!edges.Used(i2a) && !edges.Used(i2b)) 
1228
	    		    cp1 = 1; 
1229
			}
1230
		      if(ep2)
1231
			if(edgepoint.Test(p[1]))
1232
			{
1233
			  INDEX_2 i2a=INDEX_2::Sort(p[0], p[1]); 
1234
	  		  INDEX_2 i2b=INDEX_2::Sort(p[1], p[2]); 
1235
	  		  if(!edges.Used(i2a) && !edges.Used(i2b)) 
1236
	    		    cp2 = 1; 
1237
			}
1238
		      if(ep3)
1239
			if(edgepoint.Test(p[2]))
1240
			{
1241
			  INDEX_2 i2a=INDEX_2::Sort(p[2], p[1]); 
1242
	  		  INDEX_2 i2b=INDEX_2::Sort(p[3], p[2]); 
1243
	  		  if(!edges.Used(i2a) && !edges.Used(i2b)) 
1244
	    		    cp3 = 1; 
1245
			}
1246
		      if(ep4)
1247
			if(edgepoint.Test(p[3]))
1248
			{
1249
			  INDEX_2 i2a=INDEX_2::Sort(p[0], p[3]); 
1250
	  		  INDEX_2 i2b=INDEX_2::Sort(p[3], p[2]); 
1251
	  		  if(!edges.Used(i2a) && !edges.Used(i2b)) 
1252
	    		    cp4 = 1; 
1253
			}
1254
//MH*****************************************************************************************************************************
1255
		}
1256
		  else
1257
		    { 
1258
		      INDEX_2 i2;
1259
		      i2 = INDEX_2(el.PNumMod (j), el.PNumMod (j+1));
1260
		      i2.Sort();
1261
		      isedge1 = edges.Used (i2);
1262
		      if(isedge1)
1263
			{
1264
			  ep1 = 1; ep2=1;
1265
			}
1266
		      i2 = INDEX_2(el.PNumMod (j+1), el.PNumMod (j+2));
1267
		      i2.Sort();
1268
		      isedge2 = edges.Used (i2);
1269
		      if(isedge2)
1270
			{  
1271
			  ep2=1; ep3=1;
1272
			}
1273
		      i2 = INDEX_2(el.PNumMod (j+2), el.PNumMod (j+3));
1274
		      i2.Sort();
1275
		      isedge3 = edges.Used (i2); 
1276
		      
1277
		      if(isedge3)
1278
			{
1279
			  ep3=1; ep4=1;
1280
			}
1281
		      i2 = INDEX_2(el.PNumMod (j+3), el.PNumMod (j+4));
1282
		      i2.Sort();
1283
		      isedge4 = edges.Used (i2);
1284
		      if(isedge4)
1285
			{ 
1286
			  ep4=1; ep1=1;
1287
			} 
1288
		    }
1289

1290
		  int sumcp = cp1 + cp2 + cp3 + cp4;
1291
		  int sumep = ep1 + ep2 + ep3 + ep4;
1292
		  int sumedge = isedge1 + isedge2 + isedge3 + isedge4;
1293

1294
		  switch (sumedge)
1295
		    {
1296
		    case 0:
1297
		      {
1298
			switch (sumep)
1299
			  {
1300
			  case 0: 
1301
			    type = HP_QUAD; 
1302
			    break;
1303
			  case 1: 
1304
			    if (ep1) type = HP_QUAD_SINGCORNER;
1305
			    break; 
1306
			  case 2:
1307
			    {
1308
			      if (ep1 && ep2) type = HP_QUAD_0E_2VA;
1309
			      if (ep1 && ep3) type = HP_QUAD_0E_2VB;
1310
			      break;
1311
			    }
1312
			  case 3: 
1313
			    if (!ep4) type = HP_QUAD_0E_3V; 
1314
			    break; 
1315
			  case 4: 
1316
			    type = HP_QUAD_0E_4V; 
1317
			    break; 
1318
			  }
1319
			break;
1320
		      }
1321
		    case 1:
1322
		      {
1323
			if (isedge1)
1324
			  {
1325
			    switch (cp1+cp2+ep3+ep4)
1326
			      {
1327
			      case 0: 
1328
				type = HP_QUAD_SINGEDGE; 
1329
				break;
1330
			      case 1:
1331
				{
1332
				  if (cp1) type = HP_QUAD_1E_1VA;
1333
				  if (cp2) type = HP_QUAD_1E_1VB;
1334
				  if (ep3) type = HP_QUAD_1E_1VC;
1335
				  if (ep4) type = HP_QUAD_1E_1VD; 
1336
				  break; 
1337
				}
1338
			      case 2:
1339
				{
1340
				  if (cp1 && cp2) type = HP_QUAD_1E_2VA; 
1341
				  if (cp1 && ep3) type = HP_QUAD_1E_2VB; 
1342
				  if (cp1 && ep4) type = HP_QUAD_1E_2VC; 
1343
				  if (cp2 && ep3) type = HP_QUAD_1E_2VD; 
1344
				  if (cp2 && ep4) type = HP_QUAD_1E_2VE; 
1345
				  if (ep3 && ep4) type = HP_QUAD_1E_2VF; 
1346
				  break; 
1347
				}
1348
			      case 3:
1349
				{
1350
				  if (cp1 && cp2 && ep3) type = HP_QUAD_1E_3VA;
1351
				  if (cp1 && cp2 && ep4) type = HP_QUAD_1E_3VB;
1352
				  if (cp1 && ep3 && ep4) type = HP_QUAD_1E_3VC;
1353
				  if (cp2 && ep3 && ep4) type = HP_QUAD_1E_3VD;
1354
				  break;
1355
				}
1356
			      case 4:
1357
				{
1358
				  type = HP_QUAD_1E_4V; 
1359
				  break;
1360
				}
1361
			      }
1362
			  }
1363
			break;
1364
		      }
1365
		    case 2:
1366
		      {
1367
			if (isedge1 && isedge4)
1368
			  {
1369
			    if (!cp2 && !ep3 && !cp4)
1370
			      type = HP_QUAD_2E;
1371
			  
1372
			    if (cp2 && !ep3 && !cp4)
1373
			      type = HP_QUAD_2E_1VA;
1374
			    if (!cp2 && ep3 && !cp4)
1375
			      type = HP_QUAD_2E_1VB;
1376
			    if (!cp2 && !ep3 && cp4)
1377
			      type = HP_QUAD_2E_1VC;
1378

1379
			    if (cp2 && ep3 && !cp4)
1380
			      type = HP_QUAD_2E_2VA;
1381
			    if (cp2 && !ep3 && cp4)
1382
			      type = HP_QUAD_2E_2VB;
1383
			    if (!cp2 && ep3 && cp4)
1384
			      type = HP_QUAD_2E_2VC;
1385

1386
			    if (cp2 && ep3 && cp4)
1387
			      type = HP_QUAD_2E_3V;
1388
			  }
1389
			if (isedge1 && isedge3)
1390
			  {
1391
			    switch (sumcp)
1392
			      {
1393
			      case 0: 
1394
				type = HP_QUAD_2EB_0V; break;
1395
			      case 1:
1396
				{
1397
				  if (cp1) type = HP_QUAD_2EB_1VA; 
1398
				  if (cp2) type = HP_QUAD_2EB_1VB; 
1399
				  break;
1400
				}
1401
			      case 2:
1402
				{
1403
				  if (cp1 && cp2) { type = HP_QUAD_2EB_2VA; }
1404
				  if (cp1 && cp3) { type = HP_QUAD_2EB_2VB; }
1405
				  if (cp1 && cp4) { type = HP_QUAD_2EB_2VC; }
1406
				  if (cp2 && cp4) { type = HP_QUAD_2EB_2VD; }
1407
				  break;
1408
				}
1409
			      case 3:
1410
				{
1411
				  if (cp1 && cp2 && cp3) { type = HP_QUAD_2EB_3VA; }
1412
				  if (cp1 && cp2 && cp4) { type = HP_QUAD_2EB_3VB; }
1413
				  break;
1414
				}
1415
			      case 4:
1416
				{
1417
				  type = HP_QUAD_2EB_4V; break;
1418
				}
1419
			      }
1420
			  }
1421
			break;
1422
		      }
1423

1424
		    case 3:
1425
		      {
1426
			if (isedge1 && isedge2 && isedge4)
1427
			  {
1428
			    if (!cp3 && !cp4) type = HP_QUAD_3E;
1429
			    if (cp3 && !cp4) type = HP_QUAD_3E_3VA;
1430
			    if (!cp3 && cp4) type = HP_QUAD_3E_3VB;
1431
			    if (cp3 && cp4) type = HP_QUAD_3E_4V;
1432
			  }
1433
			break;
1434
		      }
1435

1436
		    case 4:
1437
		      {
1438
			type = HP_QUAD_4E;
1439
			break;
1440
		      }
1441
		    }
1442

1443
		  if (type != HP_NONE)
1444
		    {
1445
		      int pnums[4]; 
1446
		      pnums[0] = el.PNumMod (j); 
1447
		      pnums[1] = el.PNumMod (j+1);
1448
		      pnums[2] = el.PNumMod (j+2); 
1449
		      pnums[3] = el.PNumMod (j+3);
1450
		      for (int k=0;k<4;k++) el[k] = pnums[k]; 
1451
	
1452
		      /*  cout << " QUAD with pnums " << pnums[0] << "\t"  << 
1453
			pnums[1] << "\t"  << pnums[2] << "\t"  << pnums[3] 
1454
			<< endl << " of type " << type << endl; */
1455
		   		      
1456
		      break;
1457
		    }
1458
		}
1459
	      if (type == HP_NONE)
1460
		{
1461
		  (*testout) << "undefined element" << endl
1462
			     << "cp = " << cp1 << cp2 << cp3 << cp4 << endl
1463
			     << "ep = " << ep1 << ep2 << ep3 << ep4 << endl
1464
			     << "isedge = " << isedge1 << isedge2 << isedge3 
1465
			     << isedge4 << endl;
1466
		}
1467
	    
1468

1469
	      return(type);  
1470
}	    
1471

1472

1473
HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
1474
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
1475
		 INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint)
1476
{
1477
  HPREF_ELEMENT_TYPE type = HP_NONE;
1478
  
1479
  // implementation only for HP_HEX_1F_0E_0V
1480
  //                         HP_HEX_1FA_1FB_0E_0V
1481
  //                         HP_HEX 
1482
  // up to now other cases are refine dummies 
1483
  
1484
  // indices of bot,top-faces combinations
1485
  int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; 
1486
  int p[8]; 
1487
  const ELEMENT_FACE * elfaces  = MeshTopology::GetFaces (HEX);
1488
  const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (HEX);
1489
  
1490
  for(int m=0;m<6 && type == HP_NONE;m++) 
1491
    for(int j=0;j<4 && type == HP_NONE;j++) 
1492
      { 
1493
	int point_sing[8]={0,0,0,0,0,0,0,0}; 
1494
	int face_sing[6] = {0,0,0,0,0,0};
1495
	int edge_sing[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
1496
	int spoint=0, sface=0, sedge=0; 
1497
	for(int l=0;l<4;l++) 
1498
	  {
1499
	    p[l] = elfaces[index[m][0]][(4-j-l)%4]; 
1500
	    p[l+4] = elfaces[index[m][1]][(j+l)%4];
1501
	  }
1502
	
1503
	for(int l=0;l<8;l++) 
1504
	  if(cornerpoint.Test(el.PNum(p[l])))  
1505
	    { 
1506
	      point_sing[p[l]-1]=3;
1507
	      spoint++; 
1508
	    }
1509
	  else if(edgepoint.Test(el.PNum(p[l]))) point_sing[p[l]-1]=2;
1510
	  else if (facepoint[el.PNum(p[l])] == -1 || facepoint[el.PNum(p[l])] == el.GetIndex())
1511
	    point_sing[p[l]-1] = 1;   
1512
	
1513
	for(int k=0;k<12;k++)
1514
	{
1515
	  INDEX_2 i2 = INDEX_2 :: Sort(el.PNum(p[eledges[k][0]-1]),el.PNum(p[eledges[k][1]-1])); 
1516
	  if (edges.Used(i2)) 
1517
	    { 
1518
	      edge_sing[k] = 2;
1519
	      sedge++; 
1520
	    }
1521
	  else edge_sing[k] = face_edges.Used(i2);
1522
	}
1523
	
1524
	for (int k=0;k<6;k++)
1525
	{
1526
	  INDEX_3 i3; 
1527
	  
1528
	
1529
	  INDEX_4  i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]-1], el.pnums[p[elfaces[k][1]-1]-1], el.pnums[p[elfaces[k][2]-1]-1],el.pnums[p[elfaces[k][3]-1]-1]); 
1530
	  i4.Sort();
1531
	  i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3()); 
1532
	  
1533
	  if (faces.Used (i3))
1534
	    {
1535
	      
1536
	      int domnr = faces.Get(i3); 
1537
	      if (domnr == -1 || domnr == el.GetIndex())
1538
		{
1539
		  face_sing[k] = 1;
1540
		  sface++;
1541
		}
1542
	      
1543
	    } 
1544
	} 
1545
	
1546
	if(!sface && !sedge && !spoint) type = HP_HEX; 
1547
	if(!sedge && !spoint) 
1548
	  {
1549
	    if(face_sing[0] && face_sing[2] && sface==2) 
1550
	      type = HP_HEX_1FA_1FB_0E_0V; 
1551
	    if (face_sing[0] && sface==1)  
1552
	      type = HP_HEX_1F_0E_0V; 
1553
	  }
1554
	
1555
	el.type=type; 
1556

1557
	if(type != HP_NONE) 
1558
	  {
1559
	    int pnums[8]; 
1560
	    for(int l=0;l<8;l++) pnums[l] = el[p[l]-1];
1561
	    for(int l=0;l<8;l++) el[l] = pnums[l];
1562
	    /* cout << " HEX with pnums " << pnums[0] << "\t"  << 
1563
	      pnums[1] << "\t"  << pnums[2] << "\t"  << pnums[3] << "\t"  << 
1564
	      pnums[4] << "\t"  <<  pnums[5] << endl << " of type " << type << endl; */
1565
	    break; 
1566
	  }
1567
      }
1568
  
1569
  return (type); 
1570

1571
}
1572

1573
HPREF_ELEMENT_TYPE ClassifySegm(HPRefElement & hpel, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
1574
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
1575
		 INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint)
1576
{
1577
  
1578
  int cp1 = cornerpoint.Test (hpel[0]);
1579
  int cp2 = cornerpoint.Test (hpel[1]);
1580
  
1581
  INDEX_2 i2;
1582
  i2 = INDEX_2(hpel[0], hpel[1]);
1583
  i2.Sort();
1584
  if (!edges.Used (i2))
1585
    {
1586
      cp1 = edgepoint.Test (hpel[0]);
1587
      cp2 = edgepoint.Test (hpel[1]);
1588
    }
1589
  
1590
  if(!edges.Used(i2) && !face_edges.Used(i2))
1591
    {
1592
      if(facepoint[hpel[0]]!=0) cp1=1; 
1593
      if(facepoint[hpel[1]]!=0) cp2=1; 
1594
    }
1595
  
1596
  if(edges.Used(i2) && !face_edges.Used(i2))
1597
    {
1598
      if(facepoint[hpel[0]]) cp1 = 1; 
1599
	    if(facepoint[hpel[1]]) cp2 = 1; 
1600
    }
1601
  
1602
  if (!cp1 && !cp2)
1603
    hpel.type = HP_SEGM;
1604
  else if (cp1 && !cp2)
1605
    hpel.type = HP_SEGM_SINGCORNERL;
1606
  else if (!cp1 && cp2)
1607
    hpel.type = HP_SEGM_SINGCORNERR;
1608
  else
1609
    hpel.type = HP_SEGM_SINGCORNERS;
1610
  
1611
  // cout << " SEGM found with " << hpel[0] << " \t" << hpel[1] << endl << " of type " << hpel.type << endl; 
1612
  return(hpel.type) ; 
1613
} 
1614

1615

1616
HPREF_ELEMENT_TYPE ClassifyPyramid(HPRefElement & el, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom, 
1617
		 BitArray & cornerpoint, BitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges, 
1618
		 INDEX_2_HASHTABLE<int> & surf_edges, ARRAY<int, PointIndex::BASE> & facepoint)
1619
{
1620
  HPREF_ELEMENT_TYPE type = HP_NONE;
1621
  
1622
  // implementation only for HP_PYRAMID
1623
  //                         HP_PYRAMID_0E_1V
1624
  //                         HP_PYRAMID_EDGES
1625
  //                         HP_PYRAMID_1FB_0E_1VA
1626
  // up to now other cases are refine dummies 
1627
  
1628
  // indices of bot,top-faces combinations
1629
  // int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; 
1630

1631
  const ELEMENT_FACE * elfaces  = MeshTopology::GetFaces (PYRAMID);
1632
  const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (PYRAMID);
1633
  
1634
  int point_sing[5]={0,0,0,0,0}; 
1635
  int face_sing[5] = {0,0,0,0,0};
1636
  int edge_sing[8] = {0,0,0,0,0,0,0,0};
1637
  
1638
  int spoint=0, sedge=0, sface=0; 
1639
   
1640
  for(int m=0;m<4 && type == HP_NONE;m++) 
1641
    {
1642
      int p[5] = {m%4, m%4+1, m%4+2, m%4+3, 4}; 
1643

1644
      for(int l=0;l<5;l++) 
1645
	{
1646
	  if(cornerpoint.Test(el.pnums[p[l]]))  
1647
	    point_sing[l]=3;
1648
	  
1649
	     else if(edgepoint.Test(el.pnums[p[l]]))
1650
	    point_sing[l]=2;
1651
	  
1652
	  else if (facepoint[el.pnums[p[l]]] == -1 || facepoint[el.pnums[p[l]]] == el.GetIndex())
1653
	    point_sing[l] = 1;   
1654
	  
1655
	  spoint += point_sing[l]; 
1656
	}
1657
      
1658
      for(int k=0;k<8;k++)
1659
	{
1660
	  INDEX_2 i2 = INDEX_2 :: Sort(el.pnums[p[eledges[k][0]-1]],
1661
				       el.pnums[p[eledges[k][1]-1]]); 
1662
	  if (edges.Used(i2)) 
1663
	    edge_sing[k] = 2;
1664
	  else 
1665
	    edge_sing[k] = face_edges.Used(i2);
1666
	  
1667
	  sedge += edge_sing[k]; 
1668
	}
1669
  
1670
      for (int k=0;k<5;k++)
1671
	{
1672
	  INDEX_3 i3; 
1673
	  INDEX_4  i4 = INDEX_4(el.pnums[p[elfaces[k][0]-1]], el.pnums[p[elfaces[k][1]-1]], el.pnums[p[elfaces[k][2]-1]],
1674
				el.pnums[p[elfaces[k][3]-1]]); 
1675
	  i4.Sort();
1676
	  i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3()); 
1677
	  
1678
	  if (faces.Used (i3))
1679
	    {
1680
	      
1681
	      int domnr = faces.Get(i3); 
1682
	      if (domnr == -1 || domnr == el.GetIndex())
1683
		face_sing[k] = 1;
1684
	    } 
1685
	  sface +=face_sing[k]; 
1686
	} 
1687
  
1688
      if(!sface && !spoint && !sedge) return(HP_PYRAMID); 
1689
      
1690
      if(!sface && !sedge && point_sing[p[0]] == spoint) 
1691
	type = HP_PYRAMID_0E_1V; 
1692
      
1693
      if(!sface && edge_sing[0] + edge_sing[2] == sedge && 
1694
	 spoint == point_sing[0] + point_sing[1] + point_sing[3]) 
1695
	type = HP_PYRAMID_EDGES; 
1696
      
1697
      if(sface && sface == face_sing[0] && spoint == point_sing[4] + 2)
1698
	type = HP_PYRAMID_1FB_0E_1VA; 
1699
      
1700
      
1701
      if(type != HP_NONE) 
1702
	{ 
1703
	  int pnums[8]; 
1704
	  for(int l=0;l<5;l++) pnums[l] = el[p[l]];
1705
	  for(int l=0;l<5;l++) el[l] = pnums[l];
1706
	  el.type=type; 
1707
	  break; 
1708
	} 
1709
    }
1710
  
1711
  return (type); 
1712
  
1713
}
1714

1715
Product

Resources

Company
