1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5 //
6 // Contributor(s):
7 // Christophe Trophime
8 //
9
10 #include <math.h>
11 #include "GetDPConfig.h"
12 #include "Pos_Element.h"
13 #include "GeoData.h"
14 #include "Get_Geometry.h"
15 #include "Get_DofOfElement.h"
16 #include "Cal_Value.h"
17 #include "MallocUtils.h"
18 #include "Message.h"
19 #include "ProDefine.h"
20
21 #if defined(HAVE_GMSH)
22 #include <gmsh/BasisFactory.h>
23 #include <gmsh/nodalBasis.h>
24 #endif
25
26 extern struct CurrentData Current ;
27
28 /* ------------------------------------------------------------------------ */
29 /* Create/Destroy/Compare */
30 /* ------------------------------------------------------------------------ */
31
Alloc_PostElement(struct PostElement * PostElement)32 void Alloc_PostElement(struct PostElement * PostElement)
33 {
34 PostElement->NumNodes = (int *) Malloc(PostElement->NbrNodes * sizeof(int)) ;
35 /* allocate as much as possible in one step */
36 PostElement->u = (double *) Malloc(6 * PostElement->NbrNodes * sizeof(double)) ;
37 PostElement->v = &PostElement->u[PostElement->NbrNodes] ;
38 PostElement->w = &PostElement->u[2*PostElement->NbrNodes] ;
39 PostElement->x = &PostElement->u[3*PostElement->NbrNodes] ;
40 PostElement->y = &PostElement->u[4*PostElement->NbrNodes] ;
41 PostElement->z = &PostElement->u[5*PostElement->NbrNodes] ;
42 PostElement->Value = (struct Value *)
43 Malloc(PostElement->NbrNodes * sizeof(struct Value)) ;
44 }
45
Create_PostElement(int Index,int Type,int NbrNodes,int Depth)46 struct PostElement * Create_PostElement(int Index, int Type,
47 int NbrNodes, int Depth)
48 {
49 struct PostElement * PostElement ;
50
51 PostElement = (struct PostElement *) Malloc(sizeof(struct PostElement)) ;
52 PostElement->Index = Index ;
53 PostElement->Type = Type ;
54 PostElement->Depth = Depth ;
55 PostElement->NbrNodes = NbrNodes ;
56 if(NbrNodes > 0) Alloc_PostElement(PostElement);
57
58 return PostElement ;
59 }
60
Destroy_PostElement(struct PostElement * PostElement)61 void Destroy_PostElement(struct PostElement * PostElement)
62 {
63 if(PostElement->NbrNodes > 0){
64 Free(PostElement->NumNodes) ;
65 if(PostElement->u) Free(PostElement->u); /* normal case */
66 else if(PostElement->x) Free(PostElement->x); /* partial copy */
67 Free(PostElement->Value) ;
68 }
69 Free(PostElement) ;
70 }
71
NodeCopy_PostElement(struct PostElement * PostElement)72 struct PostElement * NodeCopy_PostElement(struct PostElement *PostElement)
73 {
74 struct PostElement * Copy ;
75 int i ;
76
77 Copy = (struct PostElement *) Malloc(sizeof(struct PostElement)) ;
78 Copy->Index = PostElement->Index ;
79 Copy->Type = PostElement->Type ;
80 Copy->Depth = PostElement->Depth ;
81 Copy->NbrNodes = PostElement->NbrNodes ;
82 if(Copy->NbrNodes > 0){
83 Alloc_PostElement(Copy);
84 for(i = 0 ; i < Copy->NbrNodes ; i++){
85 Copy->NumNodes[i] = PostElement->NumNodes[i];
86 Copy->u[i] = PostElement->u[i] ;
87 Copy->v[i] = PostElement->v[i] ;
88 Copy->w[i] = PostElement->w[i] ;
89 Copy->x[i] = PostElement->x[i] ;
90 Copy->y[i] = PostElement->y[i] ;
91 Copy->z[i] = PostElement->z[i] ;
92 }
93 }
94
95 return Copy ;
96 }
97
PartialCopy_PostElement(struct PostElement * PostElement)98 struct PostElement * PartialCopy_PostElement(struct PostElement *PostElement)
99 {
100 struct PostElement * Copy ;
101 int i ;
102
103 Copy = (struct PostElement *) Malloc(sizeof(struct PostElement)) ;
104 Copy->Index = PostElement->Index ;
105 Copy->Type = PostElement->Type ;
106 Copy->Depth = PostElement->Depth ;
107 Copy->NbrNodes = PostElement->NbrNodes ;
108 if(Copy->NbrNodes > 0){
109 Copy->NumNodes = NULL ;
110 Copy->u = Copy->v = Copy->w = NULL ;
111 /* allocate as much as possible in one step */
112 Copy->x = (double *) Malloc(3* Copy->NbrNodes * sizeof(double)) ;
113 Copy->y = &Copy->x[Copy->NbrNodes];
114 Copy->z = &Copy->x[2 * Copy->NbrNodes];
115 Copy->Value = (struct Value *) Malloc(Copy->NbrNodes * sizeof(struct Value)) ;
116 for(i = 0 ; i < Copy->NbrNodes ; i++){
117 Copy->x[i] = PostElement->x[i] ;
118 Copy->y[i] = PostElement->y[i] ;
119 Copy->z[i] = PostElement->z[i] ;
120 Cal_CopyValue(&PostElement->Value[i], &Copy->Value[i]);
121 }
122 }
123
124 return Copy ;
125 }
126
127 /* 2 PostElements never have the same barycenter unless they are identical */
128
fcmp_PostElement(const void * a,const void * b)129 int fcmp_PostElement(const void *a, const void *b)
130 {
131 struct PostElement *PE1, *PE2 ;
132 double s1, s2, TOL=Current.GeoData->CharacteristicLength * 1.e-12 ;
133 int i;
134
135 PE1 = *(struct PostElement**)a; PE2 = *(struct PostElement**)b;
136
137 if(PE1->NbrNodes != PE2->NbrNodes)
138 return PE1->NbrNodes - PE2->NbrNodes;
139
140 s1 = s2 = 0.0 ;
141 for(i=0;i<PE1->NbrNodes;i++){ s1 += PE1->x[i]; s2 += PE2->x[i]; }
142 if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1;
143
144 s1 = s2 = 0.0 ;
145 for(i=0;i<PE1->NbrNodes;i++){ s1 += PE1->y[i]; s2 += PE2->y[i]; }
146 if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1;
147
148 s1 = s2 = 0.0 ;
149 for(i=0;i<PE1->NbrNodes;i++){ s1 += PE1->z[i]; s2 += PE2->z[i]; }
150 if(s1-s2 > TOL) return 1; else if(s1-s2 < -TOL) return -1;
151
152 return 0;
153 }
154
fcmp_PostElement_v0(const void * a,const void * b)155 int fcmp_PostElement_v0(const void *a, const void *b)
156 {
157 return (int)( (*(struct PostElement**)a)->v[0] -
158 (*(struct PostElement**)b)->v[0] ) ;
159 }
160
fcmp_PostElement_absu0(const void * a,const void * b)161 int fcmp_PostElement_absu0(const void *a, const void *b)
162 {
163 return (int)( fabs((*(struct PostElement**)b)->u[0]) -
164 fabs((*(struct PostElement**)a)->u[0]) ) ;
165 }
166
167 /* ------------------------------------------------------------------------ */
168 /* C u t _ P o s t E l e m e n t */
169 /* ------------------------------------------------------------------------ */
170
Cut_PostElement(struct PostElement * PE,struct Geo_Element * GE,List_T * PE_L,int Index,int Depth,int Skin,int DecomposeInSimplex)171 void Cut_PostElement(struct PostElement * PE, struct Geo_Element * GE,
172 List_T * PE_L, int Index, int Depth, int Skin,
173 int DecomposeInSimplex)
174 {
175 struct Element E ;
176 struct PostElement * C[8] ;
177
178 double u01, u02, u03, u12, u13, u23 ;
179 double v01, v02, v03, v12, v13, v23 ;
180 double w01, w02, w03, w12, w13, w23 ;
181 int i, j, NbCut = 0 ;
182
183 /* Recursive division */
184
185 if(PE->Depth < Depth){
186
187 switch(PE->Type){
188
189 case POINT_ELEMENT :
190 Message::Error("Impossible to divide a Point recursively");
191 break;
192
193 case LINE :
194 case LINE_2 :
195 case LINE_3 :
196 case LINE_4 :
197 u01 = .5 * (PE->u[0] + PE->u[1]);
198 v01 = .5 * (PE->v[0] + PE->v[1]);
199 w01 = .5 * (PE->w[0] + PE->w[1]);
200
201 C[0] = Create_PostElement(Index, LINE, 2, PE->Depth) ;
202 C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
203 C[0]->u[1] = u01 ; C[0]->v[1] = v01 ; C[0]->w[1] = w01 ;
204
205 C[1] = PE ;
206 C[1]->u[0] = u01 ; C[1]->v[0] = v01 ; C[1]->w[0] = w01 ;
207
208 NbCut = 2 ;
209 break;
210
211 case TRIANGLE :
212 case TRIANGLE_2 :
213 case TRIANGLE_3 :
214 case TRIANGLE_4 :
215 u01 = .5 * (PE->u[0] + PE->u[1]); u02 = .5 * (PE->u[0] + PE->u[2]);
216 v01 = .5 * (PE->v[0] + PE->v[1]); v02 = .5 * (PE->v[0] + PE->v[2]);
217 w01 = .5 * (PE->w[0] + PE->w[1]); w02 = .5 * (PE->w[0] + PE->w[2]);
218
219 u12 = .5 * (PE->u[1] + PE->u[2]);
220 v12 = .5 * (PE->v[1] + PE->v[2]);
221 w12 = .5 * (PE->w[1] + PE->w[2]);
222
223 C[0] = Create_PostElement(Index, TRIANGLE, 3, PE->Depth) ;
224 C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
225 C[0]->u[1] = u01 ; C[0]->v[1] = v01 ; C[0]->w[1] = w01 ;
226 C[0]->u[2] = u02 ; C[0]->v[2] = v02 ; C[0]->w[2] = w02 ;
227
228 C[1] = Create_PostElement(Index, TRIANGLE, 3, PE->Depth) ;
229 C[1]->u[0] = u01 ; C[1]->v[0] = v01 ; C[1]->w[0] = w01 ;
230 C[1]->u[1] = PE->u[1] ; C[1]->v[1] = PE->v[1] ; C[1]->w[1] = PE->w[1] ;
231 C[1]->u[2] = u12 ; C[1]->v[2] = v12 ; C[1]->w[2] = w12 ;
232
233 C[2] = Create_PostElement(Index, TRIANGLE, 3, PE->Depth) ;
234 C[2]->u[0] = u02 ; C[2]->v[0] = v02 ; C[2]->w[0] = w02 ;
235 C[2]->u[1] = u12 ; C[2]->v[1] = v12 ; C[2]->w[1] = w12 ;
236 C[2]->u[2] = PE->u[2] ; C[2]->v[2] = PE->v[2] ; C[2]->w[2] = PE->w[2] ;
237
238 C[3] = PE ;
239 C[3]->u[0] = u01 ; C[3]->v[0] = v01 ; C[3]->w[0] = w01 ;
240 C[3]->u[1] = u12 ; C[3]->v[1] = v12 ; C[3]->w[1] = w12 ;
241 C[3]->u[2] = u02 ; C[3]->v[2] = v02 ; C[3]->w[2] = w02 ;
242
243 NbCut = 4 ;
244 break;
245
246 case QUADRANGLE :
247 case QUADRANGLE_2 :
248 case QUADRANGLE_3 :
249 case QUADRANGLE_4 :
250 u01 = .5 * (PE->u[0] + PE->u[1]); u12 = .5 * (PE->u[1] + PE->u[2]);
251 v01 = .5 * (PE->v[0] + PE->v[1]); v12 = .5 * (PE->v[1] + PE->v[2]);
252 w01 = .5 * (PE->w[0] + PE->w[1]); w12 = .5 * (PE->w[1] + PE->w[2]);
253
254 u23 = .5 * (PE->u[2] + PE->u[3]); u03 = .5 * (PE->u[0] + PE->u[3]);
255 v23 = .5 * (PE->v[2] + PE->v[3]); v03 = .5 * (PE->v[0] + PE->v[3]);
256 w23 = .5 * (PE->w[2] + PE->w[3]); w03 = .5 * (PE->w[0] + PE->w[3]);
257
258 u02 = .5 * (PE->u[0] + PE->u[2]);
259 v02 = .5 * (PE->v[0] + PE->v[2]);
260 w02 = .5 * (PE->w[0] + PE->w[2]);
261
262 C[0] = Create_PostElement(Index, QUADRANGLE, 4, PE->Depth) ;
263 C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
264 C[0]->u[1] = u01 ; C[0]->v[1] = v01 ; C[0]->w[1] = w01 ;
265 C[0]->u[2] = u02 ; C[0]->v[2] = v02 ; C[0]->w[2] = w02 ;
266 C[0]->u[3] = u03 ; C[0]->v[3] = v03 ; C[0]->w[3] = w03 ;
267
268 C[1] = Create_PostElement(Index, QUADRANGLE, 4, PE->Depth) ;
269 C[1]->u[0] = u01 ; C[1]->v[0] = v01 ; C[1]->w[0] = w01 ;
270 C[1]->u[1] = PE->u[1] ; C[1]->v[1] = PE->v[1] ; C[1]->w[1] = PE->w[1] ;
271 C[1]->u[2] = u12 ; C[1]->v[2] = v12 ; C[1]->w[2] = w12 ;
272 C[1]->u[3] = u02 ; C[1]->v[3] = v02 ; C[1]->w[3] = w02 ;
273
274 C[2] = Create_PostElement(Index, QUADRANGLE, 4, PE->Depth) ;
275 C[2]->u[0] = u02 ; C[2]->v[0] = v02 ; C[2]->w[0] = w02 ;
276 C[2]->u[1] = u12 ; C[2]->v[1] = v12 ; C[2]->w[1] = w12 ;
277 C[2]->u[2] = PE->u[2] ; C[2]->v[2] = PE->v[2] ; C[2]->w[2] = PE->w[2] ;
278 C[2]->u[3] = u23 ; C[2]->v[3] = v23 ; C[2]->w[3] = w23 ;
279
280 C[3] = Create_PostElement(Index, QUADRANGLE, 4, PE->Depth) ;
281 C[3]->u[0] = u03 ; C[3]->v[0] = v03 ; C[3]->w[0] = w03 ;
282 C[3]->u[1] = u02 ; C[3]->v[1] = v02 ; C[3]->w[1] = w02 ;
283 C[3]->u[2] = u23 ; C[3]->v[2] = v23 ; C[3]->w[2] = w23 ;
284 C[3]->u[3] = PE->u[3] ; C[3]->v[3] = PE->v[3] ; C[3]->w[3] = PE->w[3] ;
285
286 NbCut = 4 ;
287 break;
288
289
290 case TETRAHEDRON :
291 case TETRAHEDRON_2 :
292 case TETRAHEDRON_3 :
293 case TETRAHEDRON_4 :
294 u01 = .5 * (PE->u[0] + PE->u[1]); u02 = .5 * (PE->u[0] + PE->u[2]);
295 v01 = .5 * (PE->v[0] + PE->v[1]); v02 = .5 * (PE->v[0] + PE->v[2]);
296 w01 = .5 * (PE->w[0] + PE->w[1]); w02 = .5 * (PE->w[0] + PE->w[2]);
297
298 u03 = .5 * (PE->u[0] + PE->u[3]); u12 = .5 * (PE->u[1] + PE->u[2]);
299 v03 = .5 * (PE->v[0] + PE->v[3]); v12 = .5 * (PE->v[1] + PE->v[2]);
300 w03 = .5 * (PE->w[0] + PE->w[3]); w12 = .5 * (PE->w[1] + PE->w[2]);
301
302 u13 = .5 * (PE->u[1] + PE->u[3]); u23 = .5 * (PE->u[2] + PE->u[3]);
303 v13 = .5 * (PE->v[1] + PE->v[3]); v23 = .5 * (PE->v[2] + PE->v[3]);
304 w13 = .5 * (PE->w[1] + PE->w[3]); w23 = .5 * (PE->w[2] + PE->w[3]);
305
306 C[0] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
307 C[0]->u[0] = PE->u[0] ; C[0]->v[0] = PE->v[0] ; C[0]->w[0] = PE->w[0] ;
308 C[0]->u[1] = u01 ; C[0]->v[1] = v01 ; C[0]->w[1] = w01 ;
309 C[0]->u[2] = u02 ; C[0]->v[2] = v02 ; C[0]->w[2] = w02 ;
310 C[0]->u[3] = u03 ; C[0]->v[3] = v03 ; C[0]->w[3] = w03 ;
311
312 C[1] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
313 C[1]->u[0] = PE->u[1] ; C[1]->v[0] = PE->v[1] ; C[1]->w[0] = PE->w[1] ;
314 C[1]->u[1] = u01 ; C[1]->v[1] = v01 ; C[1]->w[1] = w01 ;
315 C[1]->u[2] = u12 ; C[1]->v[2] = v12 ; C[1]->w[2] = w12 ;
316 C[1]->u[3] = u13 ; C[1]->v[3] = v13 ; C[1]->w[3] = w13 ;
317
318 C[2] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
319 C[2]->u[0] = PE->u[2] ; C[2]->v[0] = PE->v[2] ; C[2]->w[0] = PE->w[2] ;
320 C[2]->u[1] = u02 ; C[2]->v[1] = v02 ; C[2]->w[1] = w02 ;
321 C[2]->u[2] = u12 ; C[2]->v[2] = v12 ; C[2]->w[2] = w12 ;
322 C[2]->u[3] = u23 ; C[2]->v[3] = v23 ; C[2]->w[3] = w23 ;
323
324 C[3] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
325 C[3]->u[0] = PE->u[3] ; C[3]->v[0] = PE->v[3] ; C[3]->w[0] = PE->w[3] ;
326 C[3]->u[1] = u03 ; C[3]->v[1] = v03 ; C[3]->w[1] = w03 ;
327 C[3]->u[2] = u13 ; C[3]->v[2] = v13 ; C[3]->w[2] = w13 ;
328 C[3]->u[3] = u23 ; C[3]->v[3] = v23 ; C[3]->w[3] = w23 ;
329
330 C[4] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
331 C[4]->u[0] = u01 ; C[4]->v[0] = v01 ; C[4]->w[0] = w01 ;
332 C[4]->u[1] = u02 ; C[4]->v[1] = v02 ; C[4]->w[1] = w02 ;
333 C[4]->u[2] = u03 ; C[4]->v[2] = v03 ; C[4]->w[2] = w03 ;
334 C[4]->u[3] = u23 ; C[4]->v[3] = v23 ; C[4]->w[3] = w23 ;
335
336 C[5] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
337 C[5]->u[0] = u01 ; C[5]->v[0] = v01 ; C[5]->w[0] = w01 ;
338 C[5]->u[1] = u02 ; C[5]->v[1] = v02 ; C[5]->w[1] = w02 ;
339 C[5]->u[2] = u12 ; C[5]->v[2] = v12 ; C[5]->w[2] = w12 ;
340 C[5]->u[3] = u23 ; C[5]->v[3] = v23 ; C[5]->w[3] = w23 ;
341
342 C[6] = Create_PostElement(Index, TETRAHEDRON, 4, PE->Depth) ;
343 C[6]->u[0] = u01 ; C[6]->v[0] = v01 ; C[6]->w[0] = w01 ;
344 C[6]->u[1] = u12 ; C[6]->v[1] = v12 ; C[6]->w[1] = w12 ;
345 C[6]->u[2] = u13 ; C[6]->v[2] = v13 ; C[6]->w[2] = w13 ;
346 C[6]->u[3] = u23 ; C[6]->v[3] = v23 ; C[6]->w[3] = w23 ;
347
348 C[7] = PE ;
349 C[7]->u[0] = u01 ; C[7]->v[0] = v01 ; C[7]->w[0] = w01 ;
350 C[7]->u[1] = u03 ; C[7]->v[1] = v03 ; C[7]->w[1] = w03 ;
351 C[7]->u[2] = u13 ; C[7]->v[2] = v13 ; C[7]->w[2] = w13 ;
352 C[7]->u[3] = u23 ; C[7]->v[3] = v23 ; C[7]->w[3] = w23 ;
353
354 NbCut = 8 ;
355 break ;
356
357 default :
358 Message::Error("Recursive division not implemented for Quadrangles, Hexahedra, "
359 "Prisms and Pyramids") ;
360 break;
361 }
362
363 for(i = 0 ; i < NbCut ; i++){
364 C[i]->Depth ++ ;
365 for(j = 0 ; j < C[i]->NbrNodes ; j++) C[i]->NumNodes[j] = -1 ;
366 Cut_PostElement(C[i], GE, PE_L, Index, Depth, Skin, DecomposeInSimplex);
367 }
368
369 }
370 else{
371
372 Get_InitDofOfElement(&E) ;
373 E.GeoElement = GE ;
374 E.Num = E.GeoElement->Num ;
375 E.Type = E.GeoElement->Type ;
376 E.Region = E.GeoElement->Region ;
377 Get_NodesCoordinatesOfElement(&E) ;
378
379 for(i = 0 ; i < PE->NbrNodes ; i++){
380 if( Skin == 0 && PE->Depth == 1 &&
381 ( DecomposeInSimplex == 0 || E.GeoElement->Type == LINE ||
382 E.GeoElement->Type == TRIANGLE || E.GeoElement->Type == TETRAHEDRON ) ){
383 PE->x[i] = E.x[i] ;
384 PE->y[i] = E.y[i] ;
385 PE->z[i] = E.z[i] ;
386 }
387 else{
388 Get_BFGeoElement(&E, PE->u[i], PE->v[i], PE->w[i]) ;
389
390 PE->x[i] = PE->y[i] = PE->z[i] = 0. ;
391 for (j = 0 ; j < E.GeoElement->NbrNodes ; j++) {
392 PE->x[i] += E.x[j] * E.n[j] ;
393 PE->y[i] += E.y[j] * E.n[j] ;
394 PE->z[i] += E.z[j] * E.n[j] ;
395 }
396 }
397
398 }
399
400 List_Add(PE_L, &PE);
401 }
402 }
403
404 /* ------------------------------------------------------------------------ */
405 /* F i l l _ P o s t E l e m e n t */
406 /* ------------------------------------------------------------------------ */
407
Create_HighOrderPostElement(struct Geo_Element * GE,int Index)408 static struct PostElement *Create_HighOrderPostElement(struct Geo_Element * GE,
409 int Index)
410 {
411 struct PostElement * PE = NULL;
412 #if defined(HAVE_GMSH)
413 const nodalBasis *basis = BasisFactory::getNodalBasis(GetDP2Gmsh(GE->Type));
414 if(basis){
415 int NbrNodes = basis->points.size1();
416 PE = Create_PostElement(Index, GE->Type, NbrNodes, 1) ;
417 for(int i = 0; i < NbrNodes; i++){
418 PE->NumNodes[i] = GE->NumNodes[i] ;
419 int dim = basis->points.size2();
420 PE->u[i] = (dim >= 1) ? basis->points(i, 0) : 0;
421 PE->v[i] = (dim >= 2) ? basis->points(i, 1) : 0;
422 PE->w[i] = (dim >= 3) ? basis->points(i, 2) : 0;
423 }
424 }
425 #endif
426 return PE;
427 }
428
429 #define POS_CUT_FILL Cut_PostElement(PE, GE, PE_L, Index, Depth, 0, DecomposeInSimplex)
430 #define POS_CUT_SKIN Cut_PostElement(PE, GE, PE_L, Index, Depth, 1, DecomposeInSimplex)
431
Fill_PostElement(struct Geo_Element * GE,List_T * PE_L,int Index,int Depth,int Skin,int DecomposeInSimplex,int HighOrder,int Gauss)432 void Fill_PostElement(struct Geo_Element * GE, List_T * PE_L,
433 int Index, int Depth, int Skin,
434 int DecomposeInSimplex, int HighOrder, int Gauss)
435 {
436 struct PostElement * PE ;
437
438 if(Gauss > 0){
439
440 Depth = 0;
441 int error;
442 void (*f)(int, int, double *, double *, double *, double *);
443 Get_FunctionForDefine(FunctionForGauss, GE->Type, &error, (void (**)())&f);
444 if(!error){
445 double dummy;
446 for(int i = 0; i < Gauss; i++){
447 PE = Create_PostElement(Index, POINT_ELEMENT, 1, 0) ;
448 f(Gauss, i, &PE->u[0], &PE->v[0], &PE->w[0], &dummy);
449 POS_CUT_FILL ;
450 }
451 }
452
453 }
454 else if(!Depth){
455
456 PE = Create_PostElement(Index, POINT_ELEMENT, 1, 0) ;
457 switch(GE->Type){
458 case POINT_ELEMENT :
459 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ; break ;
460 case LINE :
461 case LINE_2 :
462 case LINE_3 :
463 case LINE_4 :
464 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ; break ;
465 case TRIANGLE :
466 case TRIANGLE_2 :
467 case TRIANGLE_3 :
468 case TRIANGLE_4 :
469 PE->u[0] = 1./3.; PE->v[0] = 1./3.; PE->w[0] = 0. ; break ;
470 case QUADRANGLE :
471 case QUADRANGLE_2 :
472 case QUADRANGLE_2_8N :
473 case QUADRANGLE_3 :
474 case QUADRANGLE_4 :
475 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ; break ;
476 case TETRAHEDRON :
477 case TETRAHEDRON_2 :
478 case TETRAHEDRON_3 :
479 case TETRAHEDRON_4 :
480 PE->u[0] = 0.25 ; PE->v[0] = 0.25 ; PE->w[0] = 0.25 ; break ;
481 case HEXAHEDRON :
482 case HEXAHEDRON_2 :
483 case HEXAHEDRON_2_20N :
484 case HEXAHEDRON_3 :
485 case HEXAHEDRON_4 :
486 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ; break ;
487 case PRISM :
488 case PRISM_2 :
489 case PRISM_2_15N :
490 case PRISM_3 :
491 case PRISM_4 :
492 PE->u[0] = 1./3.; PE->v[0] = 1./3.; PE->w[0] = 0. ; break ;
493 case PYRAMID :
494 case PYRAMID_2 :
495 case PYRAMID_2_13N :
496 case PYRAMID_3 :
497 //case PYRAMID_4 :
498 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 1./3.; break ;
499 }
500 POS_CUT_FILL ;
501
502 }
503 else{
504
505 if(!Skin){
506
507 PE = NULL;
508
509 switch(GE->Type){
510
511 case POINT_ELEMENT :
512 PE = Create_PostElement(Index, POINT_ELEMENT, 1, 1) ; /* node 1 */
513 PE->NumNodes[0] = GE->NumNodes[0] ;
514 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
515 POS_CUT_FILL ;
516 break ;
517
518 case LINE : case LINE_2 : case LINE_3 : case LINE_4 :
519 if(HighOrder && GE->Type != LINE)
520 PE = Create_HighOrderPostElement(GE, Index);
521 if(!PE){
522 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 1 2 */
523 PE->NumNodes[0] = GE->NumNodes[0] ;
524 PE->NumNodes[1] = GE->NumNodes[1] ;
525 PE->u[0] =-1. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
526 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
527 }
528 POS_CUT_FILL ;
529 break ;
530
531 case TRIANGLE : case TRIANGLE_2 : case TRIANGLE_3 : case TRIANGLE_4 :
532 if(HighOrder && GE->Type != TRIANGLE)
533 PE = Create_HighOrderPostElement(GE, Index);
534 if(!PE){
535 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 2 3 */
536 PE->NumNodes[0] = GE->NumNodes[0] ;
537 PE->NumNodes[1] = GE->NumNodes[1] ;
538 PE->NumNodes[2] = GE->NumNodes[2] ;
539 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
540 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
541 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
542 }
543 POS_CUT_FILL ;
544 break ;
545
546 case QUADRANGLE : case QUADRANGLE_2 : case QUADRANGLE_2_8N :
547 case QUADRANGLE_3 : case QUADRANGLE_4 :
548 if(HighOrder && GE->Type != QUADRANGLE)
549 PE = Create_HighOrderPostElement(GE, Index);
550 if(!PE){
551 if(DecomposeInSimplex){
552 PE = Create_PostElement(Index, TRIANGLE, 3, 1); /* nodes 1 2 4 */
553 PE->NumNodes[0] = GE->NumNodes[0] ;
554 PE->NumNodes[1] = GE->NumNodes[1] ;
555 PE->NumNodes[2] = GE->NumNodes[3] ;
556 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
557 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
558 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
559 POS_CUT_FILL;
560 PE = Create_PostElement(Index, TRIANGLE, 3, 1); /* nodes 2 3 4 */
561 PE->NumNodes[0] = GE->NumNodes[1] ;
562 PE->NumNodes[1] = GE->NumNodes[2] ;
563 PE->NumNodes[2] = GE->NumNodes[3] ;
564 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
565 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
566 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
567 }
568 else{
569 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 2 3 4 */
570 PE->NumNodes[0] = GE->NumNodes[0] ;
571 PE->NumNodes[1] = GE->NumNodes[1] ;
572 PE->NumNodes[2] = GE->NumNodes[2] ;
573 PE->NumNodes[3] = GE->NumNodes[3] ;
574 PE->u[0] = -1. ; PE->v[0] = -1. ; PE->w[0] = 0. ;
575 PE->u[1] = 1. ; PE->v[1] = -1. ; PE->w[1] = 0. ;
576 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
577 PE->u[3] = -1. ; PE->v[3] = 1. ; PE->w[3] = 0. ;
578 }
579 }
580 POS_CUT_FILL ;
581 break ;
582
583 case TETRAHEDRON : case TETRAHEDRON_2 : case TETRAHEDRON_3 : case TETRAHEDRON_4 :
584 if(HighOrder && GE->Type != TETRAHEDRON)
585 PE = Create_HighOrderPostElement(GE, Index);
586 if(!PE){
587 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1) ; /* nodes 1 2 3 4 */
588 PE->NumNodes[0] = GE->NumNodes[0] ;
589 PE->NumNodes[1] = GE->NumNodes[1] ;
590 PE->NumNodes[2] = GE->NumNodes[2] ;
591 PE->NumNodes[3] = GE->NumNodes[3] ;
592 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
593 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
594 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
595 PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
596 }
597 POS_CUT_FILL;
598 break ;
599
600 case HEXAHEDRON : case HEXAHEDRON_2 : case HEXAHEDRON_2_20N :
601 case HEXAHEDRON_3 :case HEXAHEDRON_4 :
602 if(HighOrder && GE->Type != HEXAHEDRON)
603 PE = Create_HighOrderPostElement(GE, Index);
604 if(!PE){
605 if(DecomposeInSimplex){
606 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 2 3 6 */
607 PE->NumNodes[0] = GE->NumNodes[0] ;
608 PE->NumNodes[1] = GE->NumNodes[1] ;
609 PE->NumNodes[2] = GE->NumNodes[2] ;
610 PE->NumNodes[3] = GE->NumNodes[5] ;
611 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
612 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] =-1. ;
613 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
614 PE->u[3] = 1. ; PE->v[3] =-1. ; PE->w[3] = 1. ;
615 POS_CUT_FILL;
616 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 3 6 7 */
617 PE->NumNodes[0] = GE->NumNodes[0] ;
618 PE->NumNodes[1] = GE->NumNodes[2] ;
619 PE->NumNodes[2] = GE->NumNodes[5] ;
620 PE->NumNodes[3] = GE->NumNodes[6] ;
621 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
622 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
623 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ;
624 PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
625 POS_CUT_FILL;
626 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 5 6 7 */
627 PE->NumNodes[0] = GE->NumNodes[0] ;
628 PE->NumNodes[1] = GE->NumNodes[4] ;
629 PE->NumNodes[2] = GE->NumNodes[5] ;
630 PE->NumNodes[3] = GE->NumNodes[6] ;
631 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
632 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 1. ;
633 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ;
634 PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
635 POS_CUT_FILL;
636 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 3 4 7 */
637 PE->NumNodes[0] = GE->NumNodes[0] ;
638 PE->NumNodes[1] = GE->NumNodes[2] ;
639 PE->NumNodes[2] = GE->NumNodes[3] ;
640 PE->NumNodes[3] = GE->NumNodes[6] ;
641 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
642 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
643 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
644 PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
645 POS_CUT_FILL;
646 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 5 7 8 */
647 PE->NumNodes[0] = GE->NumNodes[0] ;
648 PE->NumNodes[1] = GE->NumNodes[4] ;
649 PE->NumNodes[2] = GE->NumNodes[6] ;
650 PE->NumNodes[3] = GE->NumNodes[7] ;
651 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
652 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 1. ;
653 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
654 PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
655 POS_CUT_FILL;
656 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 4 7 8 */
657 PE->NumNodes[0] = GE->NumNodes[0] ;
658 PE->NumNodes[1] = GE->NumNodes[3] ;
659 PE->NumNodes[2] = GE->NumNodes[6] ;
660 PE->NumNodes[3] = GE->NumNodes[7] ;
661 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
662 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
663 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
664 PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
665 }
666 else{
667 PE = Create_PostElement(Index, HEXAHEDRON, 8, 1) ; /* nodes 1 2 3 4 5 6 7 8 */
668 PE->NumNodes[0] = GE->NumNodes[0] ;
669 PE->NumNodes[1] = GE->NumNodes[1] ;
670 PE->NumNodes[2] = GE->NumNodes[2] ;
671 PE->NumNodes[3] = GE->NumNodes[3] ;
672 PE->NumNodes[4] = GE->NumNodes[4] ;
673 PE->NumNodes[5] = GE->NumNodes[5] ;
674 PE->NumNodes[6] = GE->NumNodes[6] ;
675 PE->NumNodes[7] = GE->NumNodes[7] ;
676 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
677 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] =-1. ;
678 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
679 PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] =-1. ;
680 PE->u[4] =-1. ; PE->v[4] =-1. ; PE->w[4] = 1. ;
681 PE->u[5] = 1. ; PE->v[5] =-1. ; PE->w[5] = 1. ;
682 PE->u[6] = 1. ; PE->v[6] = 1. ; PE->w[6] = 1. ;
683 PE->u[7] =-1. ; PE->v[7] = 1. ; PE->w[7] = 1. ;
684 }
685 }
686 POS_CUT_FILL;
687 break ;
688
689 case PRISM : case PRISM_2 : case PRISM_2_15N : case PRISM_3 : case PRISM_4 :
690 if(HighOrder && GE->Type != PRISM)
691 PE = Create_HighOrderPostElement(GE, Index);
692 if(!PE){
693 if(DecomposeInSimplex){
694 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 2 3 5 */
695 PE->NumNodes[0] = GE->NumNodes[0] ;
696 PE->NumNodes[1] = GE->NumNodes[1] ;
697 PE->NumNodes[2] = GE->NumNodes[2] ;
698 PE->NumNodes[3] = GE->NumNodes[4] ;
699 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
700 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] =-1. ;
701 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
702 PE->u[3] = 1. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
703 POS_CUT_FILL;
704 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 3 5 6 */
705 PE->NumNodes[0] = GE->NumNodes[0] ;
706 PE->NumNodes[1] = GE->NumNodes[2] ;
707 PE->NumNodes[2] = GE->NumNodes[4] ;
708 PE->NumNodes[3] = GE->NumNodes[5] ;
709 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
710 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
711 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
712 PE->u[3] = 0. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
713 POS_CUT_FILL;
714 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 4 5 6 */
715 PE->NumNodes[0] = GE->NumNodes[0] ;
716 PE->NumNodes[1] = GE->NumNodes[3] ;
717 PE->NumNodes[2] = GE->NumNodes[4] ;
718 PE->NumNodes[3] = GE->NumNodes[5] ;
719 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
720 PE->u[1] = 0. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
721 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
722 PE->u[3] = 0. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
723 }
724 else{
725 PE = Create_PostElement(Index, PRISM, 6, 1) ; /* nodes 1 2 3 4 5 6 */
726 PE->NumNodes[0] = GE->NumNodes[0] ;
727 PE->NumNodes[1] = GE->NumNodes[1] ;
728 PE->NumNodes[2] = GE->NumNodes[2] ;
729 PE->NumNodes[3] = GE->NumNodes[3] ;
730 PE->NumNodes[4] = GE->NumNodes[4] ;
731 PE->NumNodes[5] = GE->NumNodes[5] ;
732 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
733 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] =-1. ;
734 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
735 PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
736 PE->u[4] = 1. ; PE->v[4] = 0. ; PE->w[4] = 1. ;
737 PE->u[5] = 0. ; PE->v[5] = 1. ; PE->w[5] = 1. ;
738 }
739 }
740 POS_CUT_FILL;
741 break ;
742
743 case PYRAMID : case PYRAMID_2 : case PYRAMID_2_13N :
744 case PYRAMID_3 : // case PYRAMID_4 :
745 if(HighOrder && GE->Type != PYRAMID)
746 PE = Create_HighOrderPostElement(GE, Index);
747 if(!PE){
748 if(DecomposeInSimplex){
749 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 1 2 4 5 */
750 PE->NumNodes[0] = GE->NumNodes[0] ;
751 PE->NumNodes[1] = GE->NumNodes[1] ;
752 PE->NumNodes[2] = GE->NumNodes[3] ;
753 PE->NumNodes[3] = GE->NumNodes[4] ;
754 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
755 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
756 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
757 PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
758 POS_CUT_FILL;
759 PE = Create_PostElement(Index, TETRAHEDRON, 4, 1); /* nodes 2 3 4 5 */
760 PE->NumNodes[0] = GE->NumNodes[1] ;
761 PE->NumNodes[1] = GE->NumNodes[2] ;
762 PE->NumNodes[2] = GE->NumNodes[3] ;
763 PE->NumNodes[3] = GE->NumNodes[4] ;
764 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
765 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
766 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
767 PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
768 }
769 else{
770 PE = Create_PostElement(Index, PYRAMID, 5, 1) ; /* nodes 1 2 3 4 5 */
771 PE->NumNodes[0] = GE->NumNodes[0] ;
772 PE->NumNodes[1] = GE->NumNodes[1] ;
773 PE->NumNodes[2] = GE->NumNodes[2] ;
774 PE->NumNodes[3] = GE->NumNodes[3] ;
775 PE->NumNodes[4] = GE->NumNodes[4] ;
776 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
777 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
778 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
779 PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] = 0. ;
780 PE->u[4] = 0. ; PE->v[4] = 0. ; PE->w[4] = 1. ;
781 }
782 }
783 POS_CUT_FILL;
784 break ;
785
786 }
787 }
788 else { /* Skin: facets oriented with normals pointing outwards */
789
790 switch(GE->Type){
791
792 case TRIANGLE :
793 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 1 2 */
794 PE->NumNodes[0] = GE->NumNodes[0] ;
795 PE->NumNodes[1] = GE->NumNodes[1] ;
796 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
797 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
798 POS_CUT_SKIN ;
799
800 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 2 3 */
801 PE->NumNodes[0] = GE->NumNodes[1] ;
802 PE->NumNodes[1] = GE->NumNodes[2] ;
803 PE->u[0] = 1. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
804 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
805 POS_CUT_SKIN ;
806
807 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 3 1 */
808 PE->NumNodes[0] = GE->NumNodes[2] ;
809 PE->NumNodes[1] = GE->NumNodes[0] ;
810 PE->u[0] = 0. ; PE->v[0] = 1. ; PE->w[0] = 0. ;
811 PE->u[1] = 0. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
812 POS_CUT_SKIN ;
813 break ;
814
815 case QUADRANGLE :
816 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 1 2 */
817 PE->NumNodes[0] = GE->NumNodes[0] ;
818 PE->NumNodes[1] = GE->NumNodes[1] ;
819 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
820 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
821 POS_CUT_SKIN ;
822
823 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 2 3 */
824 PE->NumNodes[0] = GE->NumNodes[1] ;
825 PE->NumNodes[1] = GE->NumNodes[2] ;
826 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
827 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
828 POS_CUT_SKIN ;
829
830 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 3 4 */
831 PE->NumNodes[0] = GE->NumNodes[2] ;
832 PE->NumNodes[1] = GE->NumNodes[3] ;
833 PE->u[0] = 1. ; PE->v[0] = 1. ; PE->w[0] = 0. ;
834 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
835 POS_CUT_SKIN ;
836
837 PE = Create_PostElement(Index, LINE, 2, 1) ; /* nodes 4 1 */
838 PE->NumNodes[0] = GE->NumNodes[3] ;
839 PE->NumNodes[1] = GE->NumNodes[0] ;
840 PE->u[0] =-1. ; PE->v[0] = 1. ; PE->w[0] = 0. ;
841 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
842 POS_CUT_SKIN ;
843 break ;
844
845 case TETRAHEDRON :
846 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 2 4 */
847 PE->NumNodes[0] = GE->NumNodes[0] ;
848 PE->NumNodes[1] = GE->NumNodes[1] ;
849 PE->NumNodes[2] = GE->NumNodes[3] ;
850 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
851 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 0. ;
852 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
853 POS_CUT_SKIN;
854
855 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 3 2 */
856 PE->NumNodes[0] = GE->NumNodes[0] ;
857 PE->NumNodes[1] = GE->NumNodes[2] ;
858 PE->NumNodes[2] = GE->NumNodes[1] ;
859 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
860 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
861 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 0. ;
862 POS_CUT_SKIN;
863
864 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 4 3 */
865 PE->NumNodes[0] = GE->NumNodes[0] ;
866 PE->NumNodes[1] = GE->NumNodes[3] ;
867 PE->NumNodes[2] = GE->NumNodes[2] ;
868 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
869 PE->u[1] = 0. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
870 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
871 POS_CUT_SKIN;
872
873 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 2 3 4 */
874 PE->NumNodes[0] = GE->NumNodes[1] ;
875 PE->NumNodes[1] = GE->NumNodes[2] ;
876 PE->NumNodes[2] = GE->NumNodes[3] ;
877 PE->u[0] = 1. ; PE->v[0] = 0. ; PE->w[0] = 0. ;
878 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
879 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
880 POS_CUT_SKIN;
881 break ;
882
883 case HEXAHEDRON :
884 if(DecomposeInSimplex){
885 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 4 2 */
886 PE->NumNodes[0] = GE->NumNodes[0] ;
887 PE->NumNodes[1] = GE->NumNodes[3] ;
888 PE->NumNodes[2] = GE->NumNodes[1] ;
889 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
890 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
891 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] =-1. ;
892 POS_CUT_SKIN;
893
894 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 2 4 3 */
895 PE->NumNodes[0] = GE->NumNodes[1] ;
896 PE->NumNodes[1] = GE->NumNodes[3] ;
897 PE->NumNodes[2] = GE->NumNodes[2] ;
898 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
899 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
900 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
901 POS_CUT_SKIN;
902
903 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 5 6 8 */
904 PE->NumNodes[0] = GE->NumNodes[4] ;
905 PE->NumNodes[1] = GE->NumNodes[5] ;
906 PE->NumNodes[2] = GE->NumNodes[7] ;
907 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
908 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 1. ;
909 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
910 POS_CUT_SKIN;
911
912 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 6 7 8 */
913 PE->NumNodes[0] = GE->NumNodes[5] ;
914 PE->NumNodes[1] = GE->NumNodes[6] ;
915 PE->NumNodes[2] = GE->NumNodes[7] ;
916 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
917 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 1. ;
918 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
919 POS_CUT_SKIN;
920
921 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 5 4 */
922 PE->NumNodes[0] = GE->NumNodes[0] ;
923 PE->NumNodes[1] = GE->NumNodes[4] ;
924 PE->NumNodes[2] = GE->NumNodes[3] ;
925 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
926 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 1. ;
927 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
928 POS_CUT_SKIN;
929
930 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 5 8 4 */
931 PE->NumNodes[0] = GE->NumNodes[4] ;
932 PE->NumNodes[1] = GE->NumNodes[7] ;
933 PE->NumNodes[2] = GE->NumNodes[3] ;
934 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
935 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] = 1. ;
936 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
937 POS_CUT_SKIN;
938
939 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 2 3 6 */
940 PE->NumNodes[0] = GE->NumNodes[1] ;
941 PE->NumNodes[1] = GE->NumNodes[2] ;
942 PE->NumNodes[2] = GE->NumNodes[5] ;
943 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
944 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
945 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ;
946 POS_CUT_SKIN;
947
948 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 6 3 7 */
949 PE->NumNodes[0] = GE->NumNodes[5] ;
950 PE->NumNodes[1] = GE->NumNodes[2] ;
951 PE->NumNodes[2] = GE->NumNodes[6] ;
952 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
953 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
954 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
955 POS_CUT_SKIN;
956
957 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 5 1 6 */
958 PE->NumNodes[0] = GE->NumNodes[4] ;
959 PE->NumNodes[1] = GE->NumNodes[0] ;
960 PE->NumNodes[2] = GE->NumNodes[5] ;
961 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
962 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] =-1. ;
963 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ;
964 POS_CUT_SKIN;
965
966 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 6 1 2 */
967 PE->NumNodes[0] = GE->NumNodes[5] ;
968 PE->NumNodes[1] = GE->NumNodes[0] ;
969 PE->NumNodes[2] = GE->NumNodes[1] ;
970 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
971 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] =-1. ;
972 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] =-1. ;
973 POS_CUT_SKIN;
974
975 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 8 7 4 */
976 PE->NumNodes[0] = GE->NumNodes[7] ;
977 PE->NumNodes[1] = GE->NumNodes[6] ;
978 PE->NumNodes[2] = GE->NumNodes[3] ;
979 PE->u[0] =-1. ; PE->v[0] = 1. ; PE->w[0] = 1. ;
980 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 1. ;
981 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
982 POS_CUT_SKIN;
983
984 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 7 3 4 */
985 PE->NumNodes[0] = GE->NumNodes[6] ;
986 PE->NumNodes[1] = GE->NumNodes[2] ;
987 PE->NumNodes[2] = GE->NumNodes[3] ;
988 PE->u[0] = 1. ; PE->v[0] = 1. ; PE->w[0] = 1. ;
989 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
990 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
991 POS_CUT_SKIN;
992 }
993 else{
994 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 2 6 5 */
995 PE->NumNodes[0] = GE->NumNodes[0] ;
996 PE->NumNodes[1] = GE->NumNodes[1] ;
997 PE->NumNodes[2] = GE->NumNodes[5] ;
998 PE->NumNodes[3] = GE->NumNodes[4] ;
999 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
1000 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] =-1. ;
1001 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 1. ;
1002 PE->u[3] =-1. ; PE->v[3] =-1. ; PE->w[3] = 1. ;
1003 POS_CUT_SKIN;
1004
1005 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 4 3 2 */
1006 PE->NumNodes[0] = GE->NumNodes[0] ;
1007 PE->NumNodes[1] = GE->NumNodes[3] ;
1008 PE->NumNodes[2] = GE->NumNodes[2] ;
1009 PE->NumNodes[3] = GE->NumNodes[1] ;
1010 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
1011 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
1012 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
1013 PE->u[3] = 1. ; PE->v[3] =-1. ; PE->w[3] =-1. ;
1014 POS_CUT_SKIN;
1015
1016 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 5 8 4 */
1017 PE->NumNodes[0] = GE->NumNodes[0] ;
1018 PE->NumNodes[1] = GE->NumNodes[4] ;
1019 PE->NumNodes[2] = GE->NumNodes[7] ;
1020 PE->NumNodes[3] = GE->NumNodes[3] ;
1021 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
1022 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 1. ;
1023 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1024 PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] =-1. ;
1025 POS_CUT_SKIN;
1026
1027 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 2 3 7 6 */
1028 PE->NumNodes[0] = GE->NumNodes[1] ;
1029 PE->NumNodes[1] = GE->NumNodes[2] ;
1030 PE->NumNodes[2] = GE->NumNodes[6] ;
1031 PE->NumNodes[3] = GE->NumNodes[5] ;
1032 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] =-1. ;
1033 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
1034 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1035 PE->u[3] = 1. ; PE->v[3] =-1. ; PE->w[3] = 1. ;
1036 POS_CUT_SKIN;
1037
1038 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 3 4 8 7 */
1039 PE->NumNodes[0] = GE->NumNodes[2] ;
1040 PE->NumNodes[1] = GE->NumNodes[3] ;
1041 PE->NumNodes[2] = GE->NumNodes[7] ;
1042 PE->NumNodes[3] = GE->NumNodes[6] ;
1043 PE->u[0] = 1. ; PE->v[0] = 1. ; PE->w[0] =-1. ;
1044 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
1045 PE->u[2] =-1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1046 PE->u[3] = 1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
1047 POS_CUT_SKIN;
1048
1049 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 5 6 7 8 */
1050 PE->NumNodes[0] = GE->NumNodes[4] ;
1051 PE->NumNodes[1] = GE->NumNodes[5] ;
1052 PE->NumNodes[2] = GE->NumNodes[6] ;
1053 PE->NumNodes[3] = GE->NumNodes[7] ;
1054 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 1. ;
1055 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 1. ;
1056 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1057 PE->u[3] =-1. ; PE->v[3] = 1. ; PE->w[3] = 1. ;
1058 POS_CUT_SKIN;
1059 }
1060 break ;
1061
1062 case PRISM :
1063 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 3 2 */
1064 PE->NumNodes[0] = GE->NumNodes[0] ;
1065 PE->NumNodes[1] = GE->NumNodes[2] ;
1066 PE->NumNodes[2] = GE->NumNodes[1] ;
1067 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1068 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
1069 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] =-1. ;
1070 POS_CUT_SKIN;
1071
1072 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 4 5 6 */
1073 PE->NumNodes[0] = GE->NumNodes[3] ;
1074 PE->NumNodes[1] = GE->NumNodes[4] ;
1075 PE->NumNodes[2] = GE->NumNodes[5] ;
1076 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] = 1. ;
1077 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
1078 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1079 POS_CUT_SKIN;
1080
1081 if(DecomposeInSimplex){
1082 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 2 5 */
1083 PE->NumNodes[0] = GE->NumNodes[0] ;
1084 PE->NumNodes[1] = GE->NumNodes[1] ;
1085 PE->NumNodes[2] = GE->NumNodes[4] ;
1086 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1087 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] =-1. ;
1088 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1089 POS_CUT_SKIN;
1090
1091 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 5 4 */
1092 PE->NumNodes[0] = GE->NumNodes[0] ;
1093 PE->NumNodes[1] = GE->NumNodes[4] ;
1094 PE->NumNodes[2] = GE->NumNodes[3] ;
1095 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1096 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
1097 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1098 POS_CUT_SKIN;
1099
1100 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 6 3 */
1101 PE->NumNodes[0] = GE->NumNodes[0] ;
1102 PE->NumNodes[1] = GE->NumNodes[5] ;
1103 PE->NumNodes[2] = GE->NumNodes[2] ;
1104 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1105 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] = 1. ;
1106 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] =-1. ;
1107 POS_CUT_SKIN;
1108
1109 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 4 6 */
1110 PE->NumNodes[0] = GE->NumNodes[0] ;
1111 PE->NumNodes[1] = GE->NumNodes[3] ;
1112 PE->NumNodes[2] = GE->NumNodes[5] ;
1113 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1114 PE->u[1] = 0. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
1115 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1116 POS_CUT_SKIN;
1117
1118 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 2 3 5 */
1119 PE->NumNodes[0] = GE->NumNodes[1] ;
1120 PE->NumNodes[1] = GE->NumNodes[2] ;
1121 PE->NumNodes[2] = GE->NumNodes[4] ;
1122 PE->u[0] = 1. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1123 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
1124 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1125 POS_CUT_SKIN;
1126
1127 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 3 6 5 */
1128 PE->NumNodes[0] = GE->NumNodes[2] ;
1129 PE->NumNodes[1] = GE->NumNodes[5] ;
1130 PE->NumNodes[2] = GE->NumNodes[4] ;
1131 PE->u[0] = 0. ; PE->v[0] = 1. ; PE->w[0] =-1. ;
1132 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] = 1. ;
1133 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1134 POS_CUT_SKIN;
1135 }
1136 else{
1137 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 2 5 4 */
1138 PE->NumNodes[0] = GE->NumNodes[0] ;
1139 PE->NumNodes[1] = GE->NumNodes[1] ;
1140 PE->NumNodes[2] = GE->NumNodes[4] ;
1141 PE->NumNodes[3] = GE->NumNodes[3] ;
1142 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1143 PE->u[1] = 1. ; PE->v[1] = 0. ; PE->w[1] =-1. ;
1144 PE->u[2] = 1. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1145 PE->u[3] = 0. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
1146 POS_CUT_SKIN;
1147
1148 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 4 6 3 */
1149 PE->NumNodes[0] = GE->NumNodes[0] ;
1150 PE->NumNodes[1] = GE->NumNodes[3] ;
1151 PE->NumNodes[2] = GE->NumNodes[5] ;
1152 PE->NumNodes[3] = GE->NumNodes[2] ;
1153 PE->u[0] = 0. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1154 PE->u[1] = 0. ; PE->v[1] = 0. ; PE->w[1] = 1. ;
1155 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1156 PE->u[3] = 0. ; PE->v[3] = 1. ; PE->w[3] =-1. ;
1157 POS_CUT_SKIN;
1158
1159 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 2 3 6 5 */
1160 PE->NumNodes[0] = GE->NumNodes[1] ;
1161 PE->NumNodes[1] = GE->NumNodes[2] ;
1162 PE->NumNodes[2] = GE->NumNodes[5] ;
1163 PE->NumNodes[3] = GE->NumNodes[4] ;
1164 PE->u[0] = 1. ; PE->v[0] = 0. ; PE->w[0] =-1. ;
1165 PE->u[1] = 0. ; PE->v[1] = 1. ; PE->w[1] =-1. ;
1166 PE->u[2] = 0. ; PE->v[2] = 1. ; PE->w[2] = 1. ;
1167 PE->u[3] = 1. ; PE->v[3] = 0. ; PE->w[3] = 1. ;
1168 POS_CUT_SKIN;
1169 }
1170 break ;
1171
1172 case PYRAMID :
1173 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 2 5 */
1174 PE->NumNodes[0] = GE->NumNodes[0] ;
1175 PE->NumNodes[1] = GE->NumNodes[1] ;
1176 PE->NumNodes[2] = GE->NumNodes[4] ;
1177
1178 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
1179 PE->u[1] = 1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
1180 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1181
1182 POS_CUT_SKIN;
1183
1184 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 2 3 5 */
1185 PE->NumNodes[0] = GE->NumNodes[1] ;
1186 PE->NumNodes[1] = GE->NumNodes[2] ;
1187 PE->NumNodes[2] = GE->NumNodes[4] ;
1188
1189 PE->u[0] = 1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
1190 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
1191 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1192
1193 POS_CUT_SKIN;
1194
1195 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 3 4 5 */
1196 PE->NumNodes[0] = GE->NumNodes[2] ;
1197 PE->NumNodes[1] = GE->NumNodes[3] ;
1198 PE->NumNodes[2] = GE->NumNodes[4] ;
1199
1200 PE->u[0] = 1. ; PE->v[0] = 1. ; PE->w[0] = 0. ;
1201 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
1202 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1203
1204 POS_CUT_SKIN;
1205
1206 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 4 1 5 */
1207 PE->NumNodes[0] = GE->NumNodes[3] ;
1208 PE->NumNodes[1] = GE->NumNodes[0] ;
1209 PE->NumNodes[2] = GE->NumNodes[4] ;
1210
1211 PE->u[0] =-1. ; PE->v[0] = 1. ; PE->w[0] = 0. ;
1212 PE->u[1] =-1. ; PE->v[1] =-1. ; PE->w[1] = 0. ;
1213 PE->u[2] = 0. ; PE->v[2] = 0. ; PE->w[2] = 1. ;
1214
1215 POS_CUT_SKIN;
1216
1217 if(DecomposeInSimplex){
1218 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 3 2 */
1219 PE->NumNodes[0] = GE->NumNodes[0] ;
1220 PE->NumNodes[1] = GE->NumNodes[2] ;
1221 PE->NumNodes[2] = GE->NumNodes[1] ;
1222
1223 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
1224 PE->u[1] = 1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
1225 PE->u[2] = 1. ; PE->v[2] =-1. ; PE->w[2] = 0. ;
1226
1227 POS_CUT_SKIN;
1228
1229 PE = Create_PostElement(Index, TRIANGLE, 3, 1) ; /* nodes 1 4 3 */
1230 PE->NumNodes[0] = GE->NumNodes[0] ;
1231 PE->NumNodes[1] = GE->NumNodes[3] ;
1232 PE->NumNodes[2] = GE->NumNodes[2] ;
1233
1234 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
1235 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
1236 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
1237
1238 POS_CUT_SKIN;
1239 }
1240 else{
1241 PE = Create_PostElement(Index, QUADRANGLE, 4, 1) ; /* nodes 1 4 3 2 */
1242 PE->NumNodes[0] = GE->NumNodes[0] ;
1243 PE->NumNodes[1] = GE->NumNodes[3] ;
1244 PE->NumNodes[2] = GE->NumNodes[2] ;
1245 PE->NumNodes[3] = GE->NumNodes[1] ;
1246
1247 PE->u[0] =-1. ; PE->v[0] =-1. ; PE->w[0] = 0. ;
1248 PE->u[1] =-1. ; PE->v[1] = 1. ; PE->w[1] = 0. ;
1249 PE->u[2] = 1. ; PE->v[2] = 1. ; PE->w[2] = 0. ;
1250 PE->u[3] = 1. ; PE->v[3] =-1. ; PE->w[3] = 0. ;
1251
1252 POS_CUT_SKIN;
1253 }
1254 break ;
1255
1256 }
1257
1258 }
1259
1260 }
1261 }
1262
1263 #undef POS_CUT_FILL
1264 #undef POS_CUT_SKIN
1265
1266 /* ------------------------------------------------------------------------ */
1267 /* S o r t B y C o n n e c t i v i t y */
1268 /* ------------------------------------------------------------------------ */
1269
Compare_PostElement_Node(struct PostElement * PE1,int n1,struct PostElement * PE2,int n2)1270 int Compare_PostElement_Node(struct PostElement * PE1, int n1,
1271 struct PostElement * PE2, int n2)
1272 {
1273 double TOL=Current.GeoData->CharacteristicLength * 1.e-8;
1274 if ( (fabs(PE1->x[n1] - PE2->x[n2]) < TOL) &&
1275 (fabs(PE1->y[n1] - PE2->y[n2]) < TOL) &&
1276 (fabs(PE1->z[n1] - PE2->z[n2]) < TOL) )
1277 return 1 ;
1278 return 0 ;
1279 }
1280
Sort_PostElement_Connectivity(List_T * PostElement_L)1281 void Sort_PostElement_Connectivity(List_T *PostElement_L)
1282 {
1283 int ii, jj, start, end, iPost, NbrPost ;
1284 struct PostElement *PE, *PE2 ;
1285
1286 NbrPost = List_Nbr(PostElement_L) ;
1287
1288 /*
1289 u[0] = 1 if the element is in the ordered list, with natural orientation
1290 -1 if the element is in the ordered list, but with opposite orientation
1291 0 if the element is not in the list
1292 v[0] = relative index (to the first element) in the ordered list
1293 */
1294
1295 for(ii = 0 ; ii < NbrPost ; ii++){
1296 PE = *(struct PostElement**)List_Pointer(PostElement_L, ii);
1297 if(PE->NbrNodes != 2){
1298 Message::Error("Connectivity sorting impossible for %d-noded elements",
1299 PE->NbrNodes) ;
1300 return;
1301 }
1302 PE->u[0] = 0. ;
1303 }
1304
1305 PE = *(struct PostElement**)List_Pointer(PostElement_L, 0);
1306 PE->u[0] = 1. ; PE->v[0] = 0. ;
1307
1308 iPost = 1 ;
1309 while(iPost < NbrPost){
1310 for(ii = 0 ; ii < NbrPost ; ii++){
1311 PE = *(struct PostElement**)List_Pointer(PostElement_L, ii);
1312 if(PE->u[0]){
1313 if(PE->u[0] > 0){
1314 start = 0 ; end = 1 ;
1315 }
1316 else{
1317 start = 1 ; end = 0 ;
1318 }
1319 for(jj = 0 ; jj < NbrPost ; jj++){
1320 if(jj != ii){
1321 PE2 = *(struct PostElement**)List_Pointer(PostElement_L, jj);
1322 if(!PE2->u[0]){
1323 if(Compare_PostElement_Node(PE, end, PE2, 0)){
1324 PE2->u[0] = 1. ; PE2->v[0] = PE->v[0] + 1. ; iPost++ ;
1325 }
1326 else if (Compare_PostElement_Node(PE, start, PE2, 0)){
1327 PE2->u[0] = -1. ; PE2->v[0] = PE->v[0] - 1. ; iPost++ ;
1328 }
1329 else if (Compare_PostElement_Node(PE, start, PE2, 1)){
1330 PE2->u[0] = 1. ; PE2->v[0] = PE->v[0] - 1. ; iPost++ ;
1331 }
1332 else if (Compare_PostElement_Node(PE, end, PE2, 1)){
1333 PE2->u[0] = -1. ; PE2->v[0] = PE->v[0] + 1. ; iPost++ ;
1334 }
1335 }
1336 }
1337 }
1338 }
1339 }
1340 List_Sort(PostElement_L, fcmp_PostElement_absu0) ;
1341 }
1342 List_Sort(PostElement_L, fcmp_PostElement_v0) ;
1343 }
1344