OOFEM 3.0
Loading...
Searching...
No Matches
refinedmesh.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "refinedmesh.h"
36#include "domain.h"
37#include "element.h"
38#include "node.h"
39#include "refinedelement.h"
40
41namespace oofem {
42#define EDGE_ELEM 1
43#define FACE_ELEM 2
44#define QUAD_ELEM 3
45#define TETRA_ELEM 4
46#define HEXA_ELEM 5
47
48
49#define CUMUL_EDGES ( fe_edges )
50#define CUMUL_FACES ( fe_edges + fe_faces )
51#define CUMUL_QUADS ( fe_edges + fe_faces + fe_quads )
52#define CUMUL_TETRAS ( fe_edges + fe_faces + fe_quads + fe_tetras )
53#define CUMUL_HEXAS ( fe_edges + fe_faces + fe_quads + fe_tetras + fe_hexas )
54
55
56/* < : C style */
57#define elem_type(ID) ( ( ( ID ) <= CUMUL_EDGES ) ? EDGE_ELEM : \
58 ( ( ID ) <= CUMUL_FACES ) ? FACE_ELEM : \
59 ( ( ID ) <= CUMUL_QUADS ) ? QUAD_ELEM : \
60 ( ( ID ) <= CUMUL_TETRAS ) ? TETRA_ELEM : \
61 ( ( ID ) <= CUMUL_HEXAS ) ? HEXA_ELEM : 0 )
62
63#define is_edge(ID) ( ( elem_type(ID) == EDGE_ELEM ) ? 1 : 0 )
64#define is_face(ID) ( ( elem_type(ID) == FACE_ELEM ) ? 1 : 0 )
65#define is_quad(ID) ( ( elem_type(ID) == QUAD_ELEM ) ? 1 : 0 )
66#define is_tetra(ID) ( ( elem_type(ID) == TETRA_ELEM ) ? 1 : 0 )
67#define is_hexa(ID) ( ( elem_type(ID) == HEXA_ELEM ) ? 1 : 0 )
68
69/* <, >= : C style */
70/*
71 * #define is_edge(ID) (((ID) <= CUMUL_EDGES && (ID) > 0) ? 1 : 0)
72 * #define is_face(ID) (((ID) <= CUMUL_FACES && (ID) > CUMUL_EDGES) ? 1 : 0)
73 * #define is_quad(ID) (((ID) <= CUMUL_QUADS && (ID) > CUMUL_FACES) ? 1 : 0)
74 * #define is_tetra(ID) (((ID) <= CUMUL_TETRAS && (ID) > CUMUL_QUADS) ? 1 : 0)
75 * #define is_hexa(ID) (((ID) <= CUMUL_HEXAS && (ID) > CUMUL_TETRAS) ? 1 : 0)
76 */
77
78#define global_edge_id(ID) ( ( ID ) )
79#define global_face_id(ID) ( ( ID ) + CUMUL_EDGES )
80#define global_quad_id(ID) ( ( ID ) + CUMUL_FACES )
81#define global_tetra_id(ID) ( ( ID ) + CUMUL_QUADS )
82#define global_hexa_id(ID) ( ( ID ) + CUMUL_TETRAS )
83
84
85#define local_edge_id(ID) ( ( ID ) )
86#define local_face_id(ID) ( ( ID ) -CUMUL_EDGES )
87#define local_quad_id(ID) ( ( ID ) -CUMUL_FACES )
88#define local_tetra_id(ID) ( ( ID ) -CUMUL_QUADS )
89#define local_hexa_id(ID) ( ( ID ) -CUMUL_TETRAS )
90
91
92#define matrix_2d(ARRAY, U, V) ( ARRAY ) [ ( level + 2 ) * ( V ) + ( U ) ]
93#define matrix_3d(ARRAY, U, V, W) ( ARRAY ) [ ( level + 2 ) * ( level + 2 ) * ( W ) + ( level + 2 ) * ( V ) + ( U ) ]
94
95
96#define error_message(MSG) { OOFEM_LOG_RELEVANT( "refineMeshGlobally: %s\n", ( MSG ) ); return ( -1 ); }
97
98
99/* I do my own local and global numbering;
100 * the reason is to mimic t3d output file numbering;
101 *
102 * local numbering: each type of elem (edge, face, quad, tetra, hexa) is numbered separately
103 * but respecting the relative position in global element list
104 * global numbering: local numbering with element types ordered as edge face quad tetra hexa
105 *
106 * example:
107 * no: 1 2 3 4 5 6 7 8 9
108 * tp: t q h f t f h q f
109 *
110 * local numbering: faces 1(4) 2(6) 3(9)
111 * quads 1(2) 2(8)
112 * tetras 1(1) 2(5)
113 * hexas 1(3) 2(7)
114 * global numbering: 1(4) 2(6) 3(9) 4(2) 5(8) 6(1) 7(5) 8(3) 9(7) */
115
116
117
118int
119RefinedMesh :: refineMeshGlobally(Domain *d, int level, std :: vector< RefinedElement > &refinedElementList)
120{
121 int jj, kk = 0, j1, j2, j3, j4, type;
122 long node, elem, nd, pos = 0, p, number, size;
123 long i, j, k, m, n, node1, node2, node3, node4, elem1, elem2;
124 long refine_nodes, fine_node_id, refine_node_id;
125 long mesh_edge_id, fine_edge_id, fine_quad_id, fine_hexa_id;
126 long *tmp_array_start [ 4 ], *tmp_array_end [ 4 ];
127 long *tmp_array_cen [ 6 ], *tmp_array [ 12 ];
128 int swap [ 6 ], flag;
129
130 Element *element = NULL;
131 IntArray *connectivity = NULL, *boundary = NULL;
132
133 long fe_nodes = 0, fe_elems = 0;
134 long fe_edges = 0, fe_faces = 0, fe_quads = 0, fe_tetras = 0, fe_hexas = 0, fe_pyrams = 0, fe_wedges = 0;
135 long edge_id = 0, face_id = 0, quad_id = 0, tetra_id = 0, hexa_id = 0; // pyram_id = 0, wedge_id = 0;
136 long mesh_edges = 0, mesh_faces = 0, mesh_quads = 0;
137 long fine_nodes = 0, fine_edges = 0, fine_quads = 0, fine_hexas = 0;
138
139 fe_node_rec *fe_node = NULL, *fe_node_array = NULL;
140 fe_edge_rec *fe_edge = NULL, *fe_edge_array = NULL;
141 fe_face_rec *fe_face = NULL, *fe_face_array = NULL;
142 fe_quad_rec *fe_quad = NULL, *fe_quad_array = NULL;
143 fe_tetra_rec *fe_tetra = NULL, *fe_tetra_array = NULL;
144 fe_hexa_rec *fe_hexa = NULL, *fe_hexa_array = NULL;
145
146 fe_node_rec *fine_node = NULL, *fine_node_array = NULL;
147 /*
148 * fe_node_rec *refine_node = NULL, *refine_node_array = NULL;
149 */
150 fe_node_rec *nd1 = NULL, *nd2 = NULL, *nd3 = NULL, *nd4 = NULL;
151
152 mesh_edge_rec *mesh_edge = NULL, *mesh_edge_array = NULL;
153 mesh_face_rec *mesh_face = NULL, *mesh_face_array = NULL, *mesh_fc [ 4 ];
154 mesh_quad_rec *mesh_quad = NULL, *mesh_quad_array = NULL, *mesh_qd [ 6 ];
155
156#ifdef COMPLETE_DATA_STRUCTURE
157 tmp_face_rec *tmp_face = NULL, *tmp_face_array = NULL, *tmp_fc = NULL;
158 tmp_quad_rec *tmp_quad = NULL, *tmp_quad_array = NULL, *tmp_qd = NULL;
159#endif
160
161 tmp_tetra_rec *tmp_tetra = NULL, *tmp_tetra_array = NULL, *tmp_tet = NULL;
162 tmp_hexa_rec *tmp_hexa = NULL, *tmp_hexa_array = NULL, *tmp_hx = NULL;
163
164 fine_edge_rec *fine_edge = NULL, *fine_edge_array = NULL;
165 fine_quad_rec *fine_quad = NULL, *fine_quad_array = NULL;
166 fine_hexa_rec *fine_hexa = NULL, *fine_hexa_array = NULL;
167
168 long *node_num_elems = NULL, *node_con_elems = NULL;
169 long *node_num_nodes = NULL, *node_con_nodes = NULL;
170
171 /* caution: side and face ordering on face and tetra is not consistent with T3D */
172
173 short face_ed_nd [ 3 ] [ 2 ] = { { 0, 1 }, { 1, 2 }, { 2, 0 } };
174 short quad_ed_nd [ 4 ] [ 2 ] = { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 } };
175
176 /* note: ordering of nodes on tetra faces is important (normal: inner, outer, outer, outer) */
177
178 short tetra_fc_nd [ 4 ] [ 3 ] = { { 0, 1, 2 }, { 0, 1, 3 }, { 1, 2, 3 }, { 2, 0, 3 } };
179
180 /* note: ordering of nodes on hexa faces is important (normal: inner, outer, outer, outer, outer, inner) */
181
182 short hexa_fc_nd [ 6 ] [ 4 ] = { { 0, 1, 2, 3 }, { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 3, 7, 6 }, { 3, 0, 4, 7 }, { 7, 6, 5, 4 } };
183
184 short quad_con_nd [ 4 ] [ 2 ] = { { 1, 3 }, { 2, 0 }, { 3, 1 }, { 0, 2 } };
185 short hexa_con_nd [ 8 ] [ 3 ] = { { 1, 3, 4 }, { 2, 0, 5 }, { 3, 1, 6 }, { 0, 2, 7 }, { 7, 5, 0 }, { 4, 6, 1 }, { 5, 7, 2 }, { 6, 4, 3 } };
186
187 fe_nodes = d->giveNumberOfDofManagers();
188 fe_elems = d->giveNumberOfElements();
189
190 for ( i = 0; i < fe_elems; i++ ) {
191 element = d->giveElement(i + 1);
192 switch ( element->giveGeometryType() ) {
193 case EGT_line_1:
194 fe_edges++;
195 break;
196 case EGT_triangle_1:
197 fe_faces++;
198 break;
199 case EGT_quad_1:
200 fe_quads++;
201 break;
202 case EGT_tetra_1:
203 fe_tetras++;
204 break;
205 case EGT_hexa_1:
206 fe_hexas++;
207 break;
208 default:
209 error_message("Not supported element type");
210 break;
211 }
212 }
213
214 if ( fe_nodes == 0 ) {
215 error_message("Not enough nodes found");
216 }
217
218 if ( fe_nodes != 0 ) {
219 if ( ( fe_node_array = ( fe_node_rec * ) calloc( fe_nodes, sizeof( fe_node_rec ) ) ) == NULL ) {
220 error_message("Memory allocation error");
221 }
222 }
223
224 if ( fe_edges != 0 ) {
225 if ( ( fe_edge_array = ( fe_edge_rec * ) calloc( fe_edges, sizeof( fe_edge_rec ) ) ) == NULL ) {
226 error_message("Memory allocation error");
227 }
228 }
229
230 if ( fe_faces != 0 ) {
231 if ( ( fe_face_array = ( fe_face_rec * ) calloc( fe_faces, sizeof( fe_face_rec ) ) ) == NULL ) {
232 error_message("Memory allocation error");
233 }
234 }
235
236 if ( fe_quads != 0 ) {
237 if ( ( fe_quad_array = ( fe_quad_rec * ) calloc( fe_quads, sizeof( fe_quad_rec ) ) ) == NULL ) {
238 error_message("Memory allocation error");
239 }
240 }
241
242 if ( fe_tetras != 0 ) {
243 if ( ( fe_tetra_array = ( fe_tetra_rec * ) calloc( fe_tetras, sizeof( fe_tetra_rec ) ) ) == NULL ) {
244 error_message("Memory allocation error");
245 }
246 }
247
248 if ( fe_hexas != 0 ) {
249 if ( ( fe_hexa_array = ( fe_hexa_rec * ) calloc( fe_hexas, sizeof( fe_hexa_rec ) ) ) == NULL ) {
250 error_message("Memory allocation error");
251 }
252 }
253
254 fe_node = fe_node_array;
255 for ( i = 0; i < fe_nodes; i++, fe_node++ ) {
256 fe_node->id = i + 1;
257 }
258
259 edge_id = face_id = quad_id = tetra_id = hexa_id = 0;
260 for ( i = 0; i < fe_elems; i++ ) {
261 element = d->giveElement(i + 1);
262 switch ( element->giveGeometryType() ) {
263 case EGT_line_1:
264 fe_edge = & ( fe_edge_array [ edge_id++ ] );
265 for ( j = 0; j < 2; j++ ) {
266 nd = element->giveNode(j + 1)->giveNumber();
267 if ( nd <= 0 || nd > fe_nodes ) {
268 error_message("Invalid edge nodes");
269 }
270
271 fe_edge->node [ j ] = & ( fe_node_array [ nd - 1 ] );
272 }
273
274 break;
275 case EGT_triangle_1:
276 fe_face = & ( fe_face_array [ face_id++ ] );
277 for ( j = 0; j < 3; j++ ) {
278 nd = element->giveNode(j + 1)->giveNumber();
279 if ( nd <= 0 || nd > fe_nodes ) {
280 error_message("Invalid face nodes");
281 }
282
283 fe_face->node [ j ] = & ( fe_node_array [ nd - 1 ] );
284 }
285
286 break;
287 case EGT_quad_1:
288 fe_quad = & ( fe_quad_array [ quad_id++ ] );
289 for ( j = 0; j < 4; j++ ) {
290 nd = element->giveNode(j + 1)->giveNumber();
291 if ( nd <= 0 || nd > fe_nodes ) {
292 error_message("Invalid quad nodes");
293 }
294
295 fe_quad->node [ j ] = & ( fe_node_array [ nd - 1 ] );
296 }
297
298 break;
299 case EGT_tetra_1:
300 fe_tetra = & ( fe_tetra_array [ tetra_id++ ] );
301 for ( j = 0; j < 4; j++ ) {
302 nd = element->giveNode(j + 1)->giveNumber();
303 if ( nd <= 0 || nd > fe_nodes ) {
304 error_message("Invalid tetra nodes");
305 }
306
307 fe_tetra->node [ j ] = & ( fe_node_array [ nd - 1 ] );
308 }
309
310 break;
311 case EGT_hexa_1:
312 fe_hexa = & ( fe_hexa_array [ hexa_id++ ] );
313 for ( j = 0; j < 8; j++ ) {
314 nd = element->giveNode(j + 1)->giveNumber();
315 if ( nd <= 0 || nd > fe_nodes ) {
316 error_message("Invalid hexa nodes");
317 }
318
319 fe_hexa->node [ j ] = & ( fe_node_array [ nd - 1 ] );
320 }
321
322 break;
323 default:
324 error_message("Not supported element type");
325 break;
326 }
327 }
328
329 /* count number of elements incident to nodes */
330
331 if ( ( node_num_elems = ( long * ) calloc( fe_nodes + 1, sizeof( long ) ) ) == NULL ) {
332 error_message("Memory allocation error");
333 }
334
335 fe_edge = fe_edge_array;
336 for ( i = 0; i < fe_edges; i++, fe_edge++ ) {
337 for ( j = 0; j < 2; j++ ) {
338 nd = fe_edge->node [ j ]->id;
339 node_num_elems [ nd - 1 ]++;
340 }
341 }
342
343 fe_face = fe_face_array;
344 for ( i = 0; i < fe_faces; i++, fe_face++ ) {
345 for ( j = 0; j < 3; j++ ) {
346 nd = fe_face->node [ j ]->id;
347 node_num_elems [ nd - 1 ]++;
348 }
349 }
350
351 fe_quad = fe_quad_array;
352 for ( i = 0; i < fe_quads; i++, fe_quad++ ) {
353 for ( j = 0; j < 4; j++ ) {
354 nd = fe_quad->node [ j ]->id;
355 node_num_elems [ nd - 1 ]++;
356 }
357 }
358
359 fe_tetra = fe_tetra_array;
360 for ( i = 0; i < fe_tetras; i++, fe_tetra++ ) {
361 for ( j = 0; j < 4; j++ ) {
362 nd = fe_tetra->node [ j ]->id;
363 node_num_elems [ nd - 1 ]++;
364 }
365 }
366
367 fe_hexa = fe_hexa_array;
368 for ( i = 0; i < fe_hexas; i++, fe_hexa++ ) {
369 for ( j = 0; j < 8; j++ ) {
370 nd = fe_hexa->node [ j ]->id;
371 node_num_elems [ nd - 1 ]++;
372 }
373 }
374
375 /* recalculate the number of elements to current addresses */
376
377 pos = 0;
378 for ( i = 0; i < fe_nodes; i++ ) {
379 number = node_num_elems [ i ];
380 node_num_elems [ i ] = pos;
381 pos += number;
382 }
383
384 node_num_elems [ fe_nodes ] = size = pos;
385 if ( ( node_con_elems = ( long * ) calloc( size, sizeof( long ) ) ) == NULL ) {
386 error_message("Memory allocation error");
387 }
388
389 /* store element numbers incident to nodes */
390
391 fe_edge = fe_edge_array;
392 for ( i = 0; i < fe_edges; i++, fe_edge++ ) {
393 for ( j = 0; j < 2; j++ ) {
394 nd = fe_edge->node [ j ]->id;
395 node_con_elems [ node_num_elems [ nd - 1 ]++ ] = i + 1;
396 }
397 }
398
399 fe_face = fe_face_array;
400 for ( i = 0; i < fe_faces; i++, fe_face++ ) {
401 for ( j = 0; j < 3; j++ ) {
402 nd = fe_face->node [ j ]->id;
403 node_con_elems [ node_num_elems [ nd - 1 ]++ ] = i + 1 + fe_edges;
404 }
405 }
406
407 fe_quad = fe_quad_array;
408 for ( i = 0; i < fe_quads; i++, fe_quad++ ) {
409 for ( j = 0; j < 4; j++ ) {
410 nd = fe_quad->node [ j ]->id;
411 node_con_elems [ node_num_elems [ nd - 1 ]++ ] = i + 1 + fe_edges + fe_faces;
412 }
413 }
414
415 fe_tetra = fe_tetra_array;
416 for ( i = 0; i < fe_tetras; i++, fe_tetra++ ) {
417 for ( j = 0; j < 4; j++ ) {
418 nd = fe_tetra->node [ j ]->id;
419 node_con_elems [ node_num_elems [ nd - 1 ]++ ] = i + 1 + fe_edges + fe_faces + fe_quads;
420 }
421 }
422
423 fe_hexa = fe_hexa_array;
424 for ( i = 0; i < fe_hexas; i++, fe_hexa++ ) {
425 for ( j = 0; j < 8; j++ ) {
426 nd = fe_hexa->node [ j ]->id;
427 node_con_elems [ node_num_elems [ nd - 1 ]++ ] = i + 1 + fe_edges + fe_faces + fe_quads + fe_tetras;
428 }
429 }
430
431 /* recalculate the addresses to address of the first element */
432
433 pos = 0;
434 for ( i = 0; i < fe_nodes; i++ ) {
435 number = node_num_elems [ i ] - pos;
436 node_num_elems [ i ] = pos;
437 pos += number;
438 }
439
440#ifdef DEBUG
441 for ( i = 0; i < fe_nodes; i++ ) {
442 OOFEM_LOG_DEBUG("node %ld: %ld:", i + 1, node_num_elems [ i + 1 ] - node_num_elems [ i ]);
443 for ( j = node_num_elems [ i ]; j < node_num_elems [ i + 1 ]; j++ ) {
444 switch ( elem_type(node_con_elems [ j ]) ) {
445 case EDGE_ELEM:
446 OOFEM_LOG_DEBUG( " %ld (e)", local_edge_id(node_con_elems [ j ]) );
447 break;
448 case FACE_ELEM:
449 OOFEM_LOG_DEBUG( " %ld (f)", local_face_id(node_con_elems [ j ]) );
450 break;
451 case QUAD_ELEM:
452 OOFEM_LOG_DEBUG( " %ld (q)", local_quad_id(node_con_elems [ j ]) );
453 break;
454 case TETRA_ELEM:
455 OOFEM_LOG_DEBUG( " %ld (t)", local_tetra_id(node_con_elems [ j ]) );
456 break;
457 case HEXA_ELEM:
458 OOFEM_LOG_DEBUG( " %ld (h)", local_hexa_id(node_con_elems [ j ]) );
459 break;
460 default:
461 OOFEM_LOG_DEBUG(" %ld (?)", node_con_elems [ j ]);
462 }
463 }
464
465 OOFEM_LOG_DEBUG("\n");
466 }
467
468 OOFEM_LOG_DEBUG("\n");
469#endif
470
471 /* count number of nodes incident to nodes */
472
473 if ( ( node_num_nodes = ( long * ) calloc( fe_nodes + 1, sizeof( long ) ) ) == NULL ) {
474 error_message("Memory allocation error");
475 }
476
477 fe_node = fe_node_array;
478 for ( i = 0; i < fe_nodes; i++, fe_node++ ) {
479 node = i + 1;
480 for ( j = node_num_elems [ i ]; j < node_num_elems [ node ]; j++ ) {
481 elem = node_con_elems [ j ];
482 type = elem_type(elem); /* macro */
483 switch ( type ) {
484 case EDGE_ELEM:
485 elem = local_edge_id(elem); /* macro */
486 fe_edge = & ( fe_edge_array [ elem - 1 ] );
487 for ( k = 0; k < 2; k++ ) {
488 nd = fe_edge->node [ k ]->id;
489 if ( nd != node && nd > 0 ) {
490 fe_edge->node [ k ]->id = -nd;
491 node_num_nodes [ i ]++;
492 }
493 }
494
495 break;
496 case FACE_ELEM:
497 elem = local_face_id(elem); /* macro */
498 fe_face = & ( fe_face_array [ elem - 1 ] );
499 for ( k = 0; k < 3; k++ ) {
500 nd = fe_face->node [ k ]->id;
501 if ( nd != node && nd > 0 ) {
502 fe_face->node [ k ]->id = -nd;
503 node_num_nodes [ i ]++;
504 }
505 }
506
507 break;
508 case QUAD_ELEM:
509 elem = local_quad_id(elem); /* macro */
510 fe_quad = & ( fe_quad_array [ elem - 1 ] );
511 for ( k = 0; k < 4; k++ ) {
512 if ( fe_node == fe_quad->node [ k ] ) {
513 for ( m = 0; m < 2; m++ ) {
514 nd = fe_quad->node [ quad_con_nd [ k ] [ m ] ]->id;
515 if ( nd > 0 ) {
516 fe_quad->node [ quad_con_nd [ k ] [ m ] ]->id = -nd;
517 node_num_nodes [ i ]++;
518 }
519 }
520
521 break;
522 }
523 }
524
525 break;
526 case TETRA_ELEM:
527 elem = local_tetra_id(elem); /* macro */
528 fe_tetra = & ( fe_tetra_array [ elem - 1 ] );
529 for ( k = 0; k < 4; k++ ) {
530 nd = fe_tetra->node [ k ]->id;
531 if ( nd != node && nd > 0 ) {
532 fe_tetra->node [ k ]->id = -nd;
533 node_num_nodes [ i ]++;
534 }
535 }
536
537 break;
538 case HEXA_ELEM:
539 elem = local_hexa_id(elem); /* macro */
540 fe_hexa = & ( fe_hexa_array [ elem - 1 ] );
541 for ( k = 0; k < 8; k++ ) {
542 if ( fe_node == fe_hexa->node [ k ] ) {
543 for ( m = 0; m < 3; m++ ) {
544 nd = fe_hexa->node [ hexa_con_nd [ k ] [ m ] ]->id;
545 if ( nd > 0 ) {
546 fe_hexa->node [ hexa_con_nd [ k ] [ m ] ]->id = -nd;
547 node_num_nodes [ i ]++;
548 }
549 }
550
551 break;
552 }
553 }
554
555 break;
556 default:
557 error_message("Invalid element number");
558 break;
559 }
560 }
561
562 /* restore connected node id to positive value */
563
564 for ( j = node_num_elems [ i ]; j < node_num_elems [ node ]; j++ ) {
565 elem = node_con_elems [ j ];
566 type = elem_type(elem); /* macro */
567 switch ( type ) {
568 case EDGE_ELEM:
569 elem = local_edge_id(elem); /* macro */
570 fe_edge = & ( fe_edge_array [ elem - 1 ] );
571 for ( k = 0; k < 2; k++ ) {
572 fe_edge->node [ k ]->id = abs(fe_edge->node [ k ]->id);
573 }
574
575 break;
576 case FACE_ELEM:
577 elem = local_face_id(elem); /* macro */
578 fe_face = & ( fe_face_array [ elem - 1 ] );
579 for ( k = 0; k < 3; k++ ) {
580 fe_face->node [ k ]->id = abs(fe_face->node [ k ]->id);
581 }
582
583 break;
584 case QUAD_ELEM:
585 elem = local_quad_id(elem); /* macro */
586 fe_quad = & ( fe_quad_array [ elem - 1 ] );
587 for ( k = 0; k < 4; k++ ) {
588 fe_quad->node [ k ]->id = abs(fe_quad->node [ k ]->id);
589 }
590
591 break;
592 case TETRA_ELEM:
593 elem = local_tetra_id(elem); /* macro */
594 fe_tetra = & ( fe_tetra_array [ elem - 1 ] );
595 for ( k = 0; k < 4; k++ ) {
596 fe_tetra->node [ k ]->id = abs(fe_tetra->node [ k ]->id);
597 }
598
599 break;
600 case HEXA_ELEM:
601 elem = local_hexa_id(elem); /* macro */
602 fe_hexa = & ( fe_hexa_array [ elem - 1 ] );
603 for ( k = 0; k < 8; k++ ) {
604 fe_hexa->node [ k ]->id = abs(fe_hexa->node [ k ]->id);
605 }
606
607 break;
608 default:
609 error_message("Invalid element number");
610 break;
611 }
612 }
613 }
614
615 /* recalculate the number of nodes to current addresses */
616
617 pos = 0;
618 for ( i = 0; i < fe_nodes; i++ ) {
619 number = node_num_nodes [ i ];
620 node_num_nodes [ i ] = pos;
621 pos += number;
622 }
623
624 node_num_nodes [ fe_nodes ] = size = pos;
625 if ( ( node_con_nodes = ( long * ) calloc( size, sizeof( long ) ) ) == NULL ) {
626 error_message("Memory allocation error");
627 }
628
629 mesh_edges = size / 2;
630
631 /* store nodes incident to nodes */
632
633 fe_node = fe_node_array;
634 for ( i = 0; i < fe_nodes; i++, fe_node++ ) {
635 node = i + 1;
636 for ( j = node_num_elems [ i ]; j < node_num_elems [ node ]; j++ ) {
637 elem = node_con_elems [ j ];
638 type = elem_type(elem); /* macro */
639 switch ( type ) {
640 case EDGE_ELEM:
641 elem = local_edge_id(elem); /* macro */
642 fe_edge = & ( fe_edge_array [ elem - 1 ] );
643 for ( k = 0; k < 2; k++ ) {
644 nd = fe_edge->node [ k ]->id;
645 if ( nd != node && nd > 0 ) {
646 fe_edge->node [ k ]->id = -nd;
647 node_con_nodes [ node_num_nodes [ i ]++ ] = nd;
648 }
649 }
650
651 break;
652 case FACE_ELEM:
653 elem = local_face_id(elem); /* macro */
654 fe_face = & ( fe_face_array [ elem - 1 ] );
655 for ( k = 0; k < 3; k++ ) {
656 nd = fe_face->node [ k ]->id;
657 if ( nd != node && nd > 0 ) {
658 fe_face->node [ k ]->id = -nd;
659 node_con_nodes [ node_num_nodes [ i ]++ ] = nd;
660 }
661 }
662
663 break;
664 case QUAD_ELEM:
665 elem = local_quad_id(elem); /* macro */
666 fe_quad = & ( fe_quad_array [ elem - 1 ] );
667 for ( k = 0; k < 4; k++ ) {
668 if ( fe_node == fe_quad->node [ k ] ) {
669 for ( m = 0; m < 2; m++ ) {
670 nd = fe_quad->node [ quad_con_nd [ k ] [ m ] ]->id;
671 if ( nd > 0 ) {
672 fe_quad->node [ quad_con_nd [ k ] [ m ] ]->id = -nd;
673 node_con_nodes [ node_num_nodes [ i ]++ ] = nd;
674 }
675 }
676
677 break;
678 }
679 }
680
681 break;
682 case TETRA_ELEM:
683 elem = local_tetra_id(elem); /* macro */
684 fe_tetra = & ( fe_tetra_array [ elem - 1 ] );
685 for ( k = 0; k < 4; k++ ) {
686 nd = fe_tetra->node [ k ]->id;
687 if ( nd != node && nd > 0 ) {
688 fe_tetra->node [ k ]->id = -nd;
689 node_con_nodes [ node_num_nodes [ i ]++ ] = nd;
690 }
691 }
692
693 break;
694 case HEXA_ELEM:
695 elem = local_hexa_id(elem); /* macro */
696 fe_hexa = & ( fe_hexa_array [ elem - 1 ] );
697 for ( k = 0; k < 8; k++ ) {
698 if ( fe_node == fe_hexa->node [ k ] ) {
699 for ( m = 0; m < 3; m++ ) {
700 nd = fe_hexa->node [ hexa_con_nd [ k ] [ m ] ]->id;
701 if ( nd > 0 ) {
702 fe_hexa->node [ hexa_con_nd [ k ] [ m ] ]->id = -nd;
703 node_con_nodes [ node_num_nodes [ i ]++ ] = nd;
704 }
705 }
706
707 break;
708 }
709 }
710
711 break;
712 default:
713 error_message("Invalid element number");
714 break;
715 }
716 }
717
718 /* restore connected node id to positive value */
719
720 for ( j = node_num_elems [ i ]; j < node_num_elems [ node ]; j++ ) {
721 elem = node_con_elems [ j ];
722 type = elem_type(elem); /* macro */
723 switch ( type ) {
724 case EDGE_ELEM:
725 elem = local_edge_id(elem); /* macro */
726 fe_edge = & ( fe_edge_array [ elem - 1 ] );
727 for ( k = 0; k < 2; k++ ) {
728 fe_edge->node [ k ]->id = abs(fe_edge->node [ k ]->id);
729 }
730
731 break;
732 case FACE_ELEM:
733 elem = local_face_id(elem); /* macro */
734 fe_face = & ( fe_face_array [ elem - 1 ] );
735 for ( k = 0; k < 3; k++ ) {
736 fe_face->node [ k ]->id = abs(fe_face->node [ k ]->id);
737 }
738
739 break;
740 case QUAD_ELEM:
741 elem = local_quad_id(elem); /* macro */
742 fe_quad = & ( fe_quad_array [ elem - 1 ] );
743 for ( k = 0; k < 4; k++ ) {
744 fe_quad->node [ k ]->id = abs(fe_quad->node [ k ]->id);
745 }
746
747 break;
748 case TETRA_ELEM:
749 elem = local_tetra_id(elem); /* macro */
750 fe_tetra = & ( fe_tetra_array [ elem - 1 ] );
751 for ( k = 0; k < 4; k++ ) {
752 fe_tetra->node [ k ]->id = abs(fe_tetra->node [ k ]->id);
753 }
754
755 break;
756 case HEXA_ELEM:
757 elem = local_hexa_id(elem); /* macro */
758 fe_hexa = & ( fe_hexa_array [ elem - 1 ] );
759 for ( k = 0; k < 8; k++ ) {
760 fe_hexa->node [ k ]->id = abs(fe_hexa->node [ k ]->id);
761 }
762
763 break;
764 default:
765 error_message("Invalid element number");
766 break;
767 }
768 }
769 }
770
771 /* recalculate the addresses to address of the first node */
772
773 pos = 0;
774 for ( i = 0; i < fe_nodes; i++ ) {
775 number = node_num_nodes [ i ] - pos;
776 node_num_nodes [ i ] = pos;
777 pos += number;
778 }
779
780#ifdef DEBUG
781 for ( i = 0; i < fe_nodes; i++ ) {
782 OOFEM_LOG_DEBUG("node %ld: %ld:", i + 1, node_num_nodes [ i + 1 ] - node_num_nodes [ i ]);
783 for ( j = node_num_nodes [ i ]; j < node_num_nodes [ i + 1 ]; j++ ) {
784 OOFEM_LOG_DEBUG(" %ld", node_con_nodes [ j ]);
785 }
786
787 OOFEM_LOG_DEBUG("\n");
788 }
789
790 OOFEM_LOG_DEBUG("\n");
791#endif
792
793#ifdef COMPLETE_DATA_STRUCTURE
794 if ( fe_faces != 0 ) {
795 if ( ( tmp_face_array = ( tmp_face_rec * ) calloc( fe_faces, sizeof( tmp_face_rec ) ) ) == NULL ) {
796 error_message("Memory allocation error");
797 }
798 }
799
800 if ( fe_quads != 0 ) {
801 if ( ( tmp_quad_array = ( tmp_quad_rec * ) calloc( fe_quads, sizeof( tmp_quad_rec ) ) ) == NULL ) {
802 error_message("Memory allocation error");
803 }
804 }
805
806#endif
807
808 if ( fe_tetras != 0 ) {
809 if ( ( tmp_tetra_array = ( tmp_tetra_rec * ) calloc( fe_tetras, sizeof( tmp_tetra_rec ) ) ) == NULL ) {
810 error_message("Memory allocation error");
811 }
812 }
813
814 if ( fe_hexas != 0 ) {
815 if ( ( tmp_hexa_array = ( tmp_hexa_rec * ) calloc( fe_hexas, sizeof( tmp_hexa_rec ) ) ) == NULL ) {
816 error_message("Memory allocation error");
817 }
818 }
819
820#ifdef COMPLETE_DATA_STRUCTURE
821
822 /* setup 2D ngb elems */
823
824 fe_face = fe_face_array;
825 for ( i = 0; i < fe_faces; i++, fe_face++ ) {
826 elem1 = global_face_id(i) + 1; /* macro */
827 for ( j = 0; j < 3; j++ ) {
828 j1 = face_ed_nd [ j ] [ 0 ];
829 j2 = face_ed_nd [ j ] [ 1 ];
830 node1 = fe_face->node [ j1 ]->id;
831 node2 = fe_face->node [ j2 ]->id;
832 for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
833 elem2 = node_con_elems [ k ];
834 if ( elem2 == elem1 ) {
835 continue;
836 }
837
838 if ( is_face(elem2) == 1 || is_quad(elem2) == 1 ) { /* macro */
839 for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
840 if ( elem2 == node_con_elems [ m ] ) {
841 /* this error message is generated to prevent multiple creation of the non-manifold edge;
842 * therefore do not break the outer for cycle to detect the non-manifold edge */
843
844 if ( ( & ( tmp_face_array [ i ] ) )->ngb_elem_id [ j ] != 0 ) {
845 error_message("Domain is not 2-manifold");
846 }
847
848 ( & ( tmp_face_array [ i ] ) )->ngb_elem_id [ j ] = elem2;
849 break;
850 }
851 }
852 }
853 }
854 }
855 }
856
857 fe_quad = fe_quad_array;
858 for ( i = 0; i < fe_quads; i++, fe_quad++ ) {
859 elem1 = global_quad_id(i) + 1; /* macro */
860 for ( j = 0; j < 4; j++ ) {
861 j1 = quad_ed_nd [ j ] [ 0 ];
862 j2 = quad_ed_nd [ j ] [ 1 ];
863 node1 = fe_quad->node [ j1 ]->id;
864 node2 = fe_quad->node [ j2 ]->id;
865 for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
866 elem2 = node_con_elems [ k ];
867 if ( elem2 == elem1 ) {
868 continue;
869 }
870
871 if ( is_face(elem2) == 1 || is_quad(elem2) == 1 ) { /* macro */
872 for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
873 if ( elem2 == node_con_elems [ m ] ) {
874 /* this error message is generated to prevent multiple creation of the non-manifold edge;
875 * therefore do not break the outer for cycle to detect the non-manifold edge */
876
877 if ( ( & ( tmp_quad_array [ i ] ) )->ngb_elem_id [ j ] != 0 ) {
878 error_message("Domain is not 2-manifold");
879 }
880
881 ( & ( tmp_quad_array [ i ] ) )->ngb_elem_id [ j ] = elem2;
882 break;
883 }
884 }
885 }
886 }
887 }
888 }
889
890 #ifdef DEBUG
891 tmp_face = tmp_face_array;
892 for ( i = 0; i < fe_faces; i++, tmp_face++ ) {
893 OOFEM_LOG_DEBUG("face %ld: ", i + 1);
894 for ( j = 0; j < 3; j++ ) {
895 switch ( elem_type(tmp_face->ngb_elem_id [ j ]) ) {
896 case FACE_ELEM:
897 OOFEM_LOG_DEBUG( " %ld (f)", local_face_id(tmp_face->ngb_elem_id [ j ]) );
898 break;
899 case QUAD_ELEM:
900 OOFEM_LOG_DEBUG( " %ld (q)", local_quad_id(tmp_face->ngb_elem_id [ j ]) );
901 break;
902 }
903 }
904
905 OOFEM_LOG_DEBUG("\n");
906 }
907
908 tmp_quad = tmp_quad_array;
909 for ( i = 0; i < fe_quads; i++, tmp_quad++ ) {
910 OOFEM_LOG_DEBUG("quad %ld: ", i + 1);
911 for ( j = 0; j < 4; j++ ) {
912 switch ( elem_type(tmp_quad->ngb_elem_id [ j ]) ) {
913 case FACE_ELEM:
914 OOFEM_LOG_DEBUG( " %ld (f)", local_face_id(tmp_quad->ngb_elem_id [ j ]) );
915 break;
916 case QUAD_ELEM:
917 OOFEM_LOG_DEBUG( " %ld (q)", local_quad_id(tmp_quad->ngb_elem_id [ j ]) );
918 break;
919 }
920 }
921
922 OOFEM_LOG_DEBUG("\n");
923 }
924
925 OOFEM_LOG_DEBUG("\n");
926 #endif
927
928#endif
929
930 /* setup 3D ngb elems */
931
932 fe_tetra = fe_tetra_array;
933 for ( i = 0; i < fe_tetras; i++, fe_tetra++ ) {
934 elem1 = global_tetra_id(i) + 1; /* macro */
935 mesh_faces += 4;
936 for ( j = 0; j < 4; j++ ) {
937 j1 = tetra_fc_nd [ j ] [ 0 ];
938 j2 = tetra_fc_nd [ j ] [ 1 ];
939 j3 = tetra_fc_nd [ j ] [ 2 ];
940 node1 = fe_tetra->node [ j1 ]->id;
941 node2 = fe_tetra->node [ j2 ]->id;
942 node3 = fe_tetra->node [ j3 ]->id;
943 for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
944 elem2 = node_con_elems [ k ];
945 if ( elem2 == elem1 ) {
946 continue;
947 }
948
949 if ( is_tetra(elem2) == 0 ) {
950 continue; /* macro */
951 }
952
953 for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
954 if ( elem2 != node_con_elems [ m ] ) {
955 continue;
956 }
957
958 for ( n = node_num_elems [ node3 - 1 ]; n < node_num_elems [ node3 ]; n++ ) {
959 if ( elem2 == node_con_elems [ n ] ) {
960 ( & ( tmp_tetra_array [ i ] ) )->ngb_elem_id [ j ] = elem2;
961 if ( elem1 > elem2 ) {
962 mesh_faces--;
963 }
964
965 k = node_num_elems [ node1 ];
966 m = node_num_elems [ node2 ];
967 break;
968 }
969 }
970 }
971 }
972 }
973 }
974
975 fe_hexa = fe_hexa_array;
976 for ( i = 0; i < fe_hexas; i++, fe_hexa++ ) {
977 elem1 = global_hexa_id(i) + 1; /* macro */
978 mesh_quads += 6;
979 for ( j = 0; j < 6; j++ ) {
980 j1 = hexa_fc_nd [ j ] [ 0 ];
981 j2 = hexa_fc_nd [ j ] [ 1 ];
982 j3 = hexa_fc_nd [ j ] [ 2 ];
983 j4 = hexa_fc_nd [ j ] [ 3 ];
984 node1 = fe_hexa->node [ j1 ]->id;
985 node2 = fe_hexa->node [ j2 ]->id;
986 node3 = fe_hexa->node [ j3 ]->id;
987 node4 = fe_hexa->node [ j4 ]->id;
988 for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
989 elem2 = node_con_elems [ k ];
990 if ( elem2 == elem1 ) {
991 continue;
992 }
993
994 if ( is_hexa(elem2) == 0 ) {
995 continue; /* macro */
996 }
997
998 for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
999 if ( elem2 != node_con_elems [ m ] ) {
1000 continue;
1001 }
1002
1003 for ( n = node_num_elems [ node3 - 1 ]; n < node_num_elems [ node3 ]; n++ ) {
1004 if ( elem2 != node_con_elems [ n ] ) {
1005 continue;
1006 }
1007
1008 for ( p = node_num_elems [ node4 - 1 ]; p < node_num_elems [ node4 ]; p++ ) {
1009 if ( elem2 == node_con_elems [ p ] ) {
1010 ( & ( tmp_hexa_array [ i ] ) )->ngb_elem_id [ j ] = elem2;
1011 if ( elem1 > elem2 ) {
1012 mesh_quads--;
1013 }
1014
1015 k = node_num_elems [ node1 ];
1016 m = node_num_elems [ node2 ];
1017 n = node_num_elems [ node3 ];
1018 break;
1019 }
1020 }
1021 }
1022 }
1023 }
1024 }
1025 }
1026
1027#ifdef DEBUG
1028 tmp_tetra = tmp_tetra_array;
1029 for ( i = 0; i < fe_tetras; i++, tmp_tetra++ ) {
1030 OOFEM_LOG_DEBUG("tetra %ld: ", i + 1);
1031 for ( j = 0; j < 4; j++ ) {
1032 if ( tmp_tetra->ngb_elem_id [ j ] == 0 ) {
1033 continue;
1034 }
1035
1036 OOFEM_LOG_DEBUG( " %ld (t)", local_tetra_id(tmp_tetra->ngb_elem_id [ j ]) );
1037 }
1038
1039 OOFEM_LOG_DEBUG("\n");
1040 }
1041
1042 tmp_hexa = tmp_hexa_array;
1043 for ( i = 0; i < fe_hexas; i++, tmp_hexa++ ) {
1044 OOFEM_LOG_DEBUG("hexa %ld: ", i + 1);
1045 for ( j = 0; j < 6; j++ ) {
1046 if ( tmp_hexa->ngb_elem_id [ j ] == 0 ) {
1047 continue;
1048 }
1049
1050 OOFEM_LOG_DEBUG( " %ld (h)", local_hexa_id(tmp_hexa->ngb_elem_id [ j ]) );
1051 }
1052
1053 OOFEM_LOG_DEBUG("\n");
1054 }
1055
1056 OOFEM_LOG_DEBUG("\n");
1057#endif
1058
1059 if ( mesh_edges != 0 ) {
1060 if ( ( mesh_edge_array = ( mesh_edge_rec * ) calloc( mesh_edges, sizeof( mesh_edge_rec ) ) ) == NULL ) {
1061 error_message("Memory allocation error");
1062 }
1063 }
1064
1065 if ( mesh_faces != 0 ) {
1066 if ( ( mesh_face_array = ( mesh_face_rec * ) calloc( mesh_faces, sizeof( mesh_face_rec ) ) ) == NULL ) {
1067 error_message("Memory allocation error");
1068 }
1069 }
1070
1071 if ( mesh_quads != 0 ) {
1072 if ( ( mesh_quad_array = ( mesh_quad_rec * ) calloc( mesh_quads, sizeof( mesh_quad_rec ) ) ) == NULL ) {
1073 error_message("Memory allocation error");
1074 }
1075 }
1076
1077 fine_nodes = mesh_edges + mesh_faces + mesh_quads + fe_faces + fe_quads + fe_tetras + fe_pyrams + fe_wedges + fe_hexas;
1078
1079 fine_edges = fe_edges * 2;
1080 fine_quads = fe_faces * 3 + fe_quads * 4;
1081 fine_hexas = fe_tetras * 4 + fe_pyrams * 8 + fe_wedges * 8 + fe_hexas * 8;
1082
1083 if ( fine_nodes != 0 ) {
1084 if ( ( fine_node_array = ( fe_node_rec * ) calloc( fine_nodes, sizeof( fe_node_rec ) ) ) == NULL ) {
1085 error_message("Memory allocation error");
1086 }
1087 }
1088
1089 if ( fine_edges != 0 ) {
1090 if ( ( fine_edge_array = ( fine_edge_rec * ) calloc( fine_edges, sizeof( fine_edge_rec ) ) ) == NULL ) {
1091 error_message("Memory allocation error");
1092 }
1093 }
1094
1095 if ( fine_quads != 0 ) {
1096 if ( ( fine_quad_array = ( fine_quad_rec * ) calloc( fine_quads, sizeof( fine_quad_rec ) ) ) == NULL ) {
1097 error_message("Memory allocation error");
1098 }
1099 }
1100
1101 if ( fine_hexas != 0 ) {
1102 if ( ( fine_hexa_array = ( fine_hexa_rec * ) calloc( fine_hexas, sizeof( fine_hexa_rec ) ) ) == NULL ) {
1103 error_message("Memory allocation error");
1104 }
1105 }
1106
1107 fine_node_id = fe_nodes;
1108 fine_node = fine_node_array;
1109
1110 mesh_edge_id = 0;
1111 mesh_edge = mesh_edge_array;
1112
1113 if ( fe_edges != 0 ) {
1114 /* ensure that mesh_edges for fe_edges are created first */
1115
1116 fe_edge = fe_edge_array;
1117 for ( i = 0; i < fe_edges; i++, fe_edge++ ) {
1118 node1 = ( nd1 = fe_edge->node [ 0 ] )->id;
1119 node2 = ( nd2 = fe_edge->node [ 1 ] )->id;
1120
1121 fine_node->id = ++fine_node_id;
1122
1123 /* note: do not swap nodes, keep tham in the same order as on fe_edge !!! */
1124
1125 mesh_edge->node [ 0 ] = nd1;
1126 mesh_edge->node [ 1 ] = nd2;
1127
1128 mesh_edge->mid_node = fine_node++;
1129 mesh_edge_id++;
1130
1131 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1132 if ( node_con_nodes [ k ] == node2 ) {
1133 node_con_nodes [ k ] = -mesh_edge_id;
1134 break;
1135 }
1136 }
1137
1138 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1139 if ( node_con_nodes [ k ] == node1 ) {
1140 node_con_nodes [ k ] = -mesh_edge_id;
1141 break;
1142 }
1143 }
1144
1145 mesh_edge++;
1146 }
1147 }
1148
1149#ifdef COMPLETE_DATA_STRUCTURE
1150
1151 /* create mesh edges */
1152
1153 fe_face = fe_face_array;
1154 tmp_face = tmp_face_array;
1155 for ( i = 0; i < fe_faces; i++, fe_face++, tmp_face++ ) {
1156 elem1 = global_face_id(i) + 1; /* macro */
1157 for ( j = 0; j < 3; j++ ) {
1158 elem2 = tmp_face->ngb_elem_id [ j ];
1159
1160 if ( elem2 == 0 || elem1 < elem2 ) {
1161 j1 = face_ed_nd [ j ] [ 0 ];
1162 j2 = face_ed_nd [ j ] [ 1 ];
1163 node1 = ( nd1 = fe_face->node [ j1 ] )->id;
1164 node2 = ( nd2 = fe_face->node [ j2 ] )->id;
1165
1166 fine_node->id = ++fine_node_id;
1167
1168 /* note: ensure that starting node has smaller id */
1169
1170 if ( node1 < node2 ) {
1171 mesh_edge->node [ 0 ] = nd1;
1172 mesh_edge->node [ 1 ] = nd2;
1173 } else {
1174 mesh_edge->node [ 0 ] = nd2;
1175 mesh_edge->node [ 1 ] = nd1;
1176 }
1177
1178 mesh_edge->mid_node = fine_node++;
1179 mesh_edge_id++;
1180
1181 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1182 node = node_con_nodes [ k ];
1183 if ( node == node2 ) {
1184 node_con_nodes [ k ] = -mesh_edge_id;
1185 break;
1186 }
1187 }
1188
1189 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1190 node = node_con_nodes [ k ];
1191 if ( node == node1 ) {
1192 node_con_nodes [ k ] = -mesh_edge_id;
1193 break;
1194 }
1195 }
1196
1197 tmp_face->edge [ j ] = mesh_edge;
1198
1199 if ( elem2 != 0 ) {
1200 if ( is_face(elem2) == 1 ) {
1201 elem = local_face_id(elem2); /* macro */
1202 tmp_fc = & ( tmp_face_array [ elem - 1 ] );
1203 for ( k = 0; k < 3; k++ ) {
1204 if ( tmp_fc->ngb_elem_id [ k ] == elem1 ) {
1205 tmp_fc->edge [ k ] = mesh_edge;
1206 break;
1207 }
1208 }
1209 }
1210
1211 if ( is_quad(elem2) == 1 ) {
1212 elem = local_quad_id(elem2); /* macro */
1213 tmp_qd = & ( tmp_quad_array [ elem - 1 ] );
1214 for ( k = 0; k < 4; k++ ) {
1215 if ( tmp_qd->ngb_elem_id [ k ] == elem1 ) {
1216 tmp_qd->edge [ k ] = mesh_edge;
1217 break;
1218 }
1219 }
1220 }
1221 }
1222
1223 mesh_edge++;
1224 }
1225 }
1226 }
1227
1228 fe_quad = fe_quad_array;
1229 tmp_quad = tmp_quad_array;
1230 for ( i = 0; i < fe_quads; i++, fe_quad++, tmp_quad++ ) {
1231 elem1 = global_quad_id(i) + 1; /* macro */
1232 for ( j = 0; j < 4; j++ ) {
1233 elem2 = tmp_quad->ngb_elem_id [ j ];
1234
1235 if ( elem2 == 0 || elem1 < elem2 ) {
1236 j1 = quad_ed_nd [ j ] [ 0 ];
1237 j2 = quad_ed_nd [ j ] [ 1 ];
1238 node1 = ( nd1 = fe_quad->node [ j1 ] )->id;
1239 node2 = ( nd2 = fe_quad->node [ j2 ] )->id;
1240
1241 fine_node->id = ++fine_node_id;
1242
1243 /* note: ensure that starting node has smaller id */
1244
1245 if ( node1 < node2 ) {
1246 mesh_edge->node [ 0 ] = nd1;
1247 mesh_edge->node [ 1 ] = nd2;
1248 } else {
1249 mesh_edge->node [ 0 ] = nd2;
1250 mesh_edge->node [ 1 ] = nd1;
1251 }
1252
1253 mesh_edge->mid_node = fine_node++;
1254 mesh_edge_id++;
1255
1256 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1257 node = node_con_nodes [ k ];
1258 if ( node == node2 ) {
1259 node_con_nodes [ k ] = -mesh_edge_id;
1260 break;
1261 }
1262 }
1263
1264 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1265 node = node_con_nodes [ k ];
1266 if ( node == node1 ) {
1267 node_con_nodes [ k ] = -mesh_edge_id;
1268 break;
1269 }
1270 }
1271
1272 tmp_quad->edge [ j ] = mesh_edge;
1273
1274 if ( elem2 != 0 ) {
1275 if ( is_face(elem2) == 1 ) {
1276 elem = local_face_id(elem2); /* macro */
1277 tmp_fc = & ( tmp_face_array [ elem - 1 ] );
1278 for ( k = 0; k < 3; k++ ) {
1279 if ( tmp_fc->ngb_elem_id [ k ] == elem1 ) {
1280 tmp_fc->edge [ k ] = mesh_edge;
1281 break;
1282 }
1283 }
1284 }
1285
1286 if ( is_quad(elem2) == 1 ) {
1287 elem = local_quad_id(elem2); /* macro */
1288 tmp_qd = & ( tmp_quad_array [ elem - 1 ] );
1289 for ( k = 0; k < 4; k++ ) {
1290 if ( tmp_qd->ngb_elem_id [ k ] == elem1 ) {
1291 tmp_qd->edge [ k ] = mesh_edge;
1292 break;
1293 }
1294 }
1295 }
1296 }
1297
1298 mesh_edge++;
1299 }
1300 }
1301 }
1302
1303#endif
1304
1305 /* create mesh edges (only on remaining 3D mesh if COMPLETE_DATA_STRUCTURE is defined) */
1306
1307 for ( i = 0; i < fe_nodes; i++ ) {
1308 node = i + 1;
1309 for ( j = node_num_nodes [ i ]; j < node_num_nodes [ node ]; j++ ) {
1310 /* note: negative nd means that there is already mesh_edge number */
1311
1312 nd = node_con_nodes [ j ];
1313 if ( nd < 0 ) {
1314 node_con_nodes [ j ] = -nd;
1315 continue;
1316 }
1317
1318#ifdef DEBUG
1319 if ( node > nd ) {
1320 error_message("Unexpected situation");
1321 }
1322
1323#endif
1324
1325 nd1 = & ( fe_node_array [ node - 1 ] );
1326 nd2 = & ( fe_node_array [ nd - 1 ] );
1327
1328 fine_node->id = ++fine_node_id;
1329
1330 /* note: nd1 should be smaller than nd2 */
1331
1332 mesh_edge->node [ 0 ] = nd1;
1333 mesh_edge->node [ 1 ] = nd2;
1334 mesh_edge->mid_node = fine_node++;
1335 mesh_edge++;
1336
1337 node_con_nodes [ j ] = ++mesh_edge_id;
1338
1339 for ( k = node_num_nodes [ nd - 1 ]; k < node_num_nodes [ nd ]; k++ ) {
1340 if ( node_con_nodes [ k ] == node ) {
1341 node_con_nodes [ k ] = -mesh_edge_id;
1342 break;
1343 }
1344 }
1345 }
1346 }
1347
1348 /* create mesh faces */
1349
1350 mesh_face = mesh_face_array;
1351
1352 fe_tetra = fe_tetra_array;
1353 tmp_tetra = tmp_tetra_array;
1354 for ( i = 0; i < fe_tetras; i++, fe_tetra++, tmp_tetra++ ) {
1355 elem1 = global_tetra_id(i) + 1; /* macro */
1356 for ( j = 0; j < 4; j++ ) {
1357 elem2 = tmp_tetra->ngb_elem_id [ j ];
1358
1359 if ( elem2 == 0 || elem1 < elem2 ) {
1360 j1 = tetra_fc_nd [ j ] [ 0 ];
1361 j2 = tetra_fc_nd [ j ] [ 1 ];
1362 j3 = tetra_fc_nd [ j ] [ 2 ];
1363 nd1 = fe_tetra->node [ j1 ];
1364 nd2 = fe_tetra->node [ j2 ];
1365 nd3 = fe_tetra->node [ j3 ];
1366
1367 fine_node->id = ++fine_node_id;
1368
1369 mesh_face->node [ 0 ] = nd1;
1370 mesh_face->node [ 1 ] = nd2;
1371 mesh_face->node [ 2 ] = nd3;
1372 mesh_face->mid_node = fine_node++;
1373
1374 tmp_tetra->face [ j ] = mesh_face;
1375
1376 if ( elem2 != 0 ) {
1377 elem = local_tetra_id(elem2); /* macro */
1378 tmp_tet = & ( tmp_tetra_array [ elem - 1 ] );
1379 for ( k = 0; k < 4; k++ ) {
1380 if ( tmp_tet->ngb_elem_id [ k ] == elem1 ) {
1381 tmp_tet->face [ k ] = mesh_face;
1382 break;
1383 }
1384 }
1385 }
1386
1387 mesh_face++;
1388 }
1389 }
1390 }
1391
1392 /* create mesh quads */
1393
1394 mesh_quad = mesh_quad_array;
1395
1396 fe_hexa = fe_hexa_array;
1397 tmp_hexa = tmp_hexa_array;
1398 for ( i = 0; i < fe_hexas; i++, fe_hexa++, tmp_hexa++ ) {
1399 elem1 = global_hexa_id(i) + 1; /* macro */
1400 for ( j = 0; j < 6; j++ ) {
1401 elem2 = tmp_hexa->ngb_elem_id [ j ];
1402
1403 if ( elem2 == 0 || elem1 < elem2 ) {
1404 j1 = hexa_fc_nd [ j ] [ 0 ];
1405 j2 = hexa_fc_nd [ j ] [ 1 ];
1406 j3 = hexa_fc_nd [ j ] [ 2 ];
1407 j4 = hexa_fc_nd [ j ] [ 3 ];
1408 nd1 = fe_hexa->node [ j1 ];
1409 nd2 = fe_hexa->node [ j2 ];
1410 nd3 = fe_hexa->node [ j3 ];
1411 nd4 = fe_hexa->node [ j4 ];
1412
1413 fine_node->id = ++fine_node_id;
1414
1415 mesh_quad->node [ 0 ] = nd1;
1416 mesh_quad->node [ 1 ] = nd2;
1417 mesh_quad->node [ 2 ] = nd3;
1418 mesh_quad->node [ 3 ] = nd4;
1419 mesh_quad->mid_node = fine_node++;
1420
1421 tmp_hexa->quad [ j ] = mesh_quad;
1422
1423 if ( elem2 != 0 ) {
1424 elem = local_hexa_id(elem2); /* macro */
1425 tmp_hx = & ( tmp_hexa_array [ elem - 1 ] );
1426 for ( k = 0; k < 6; k++ ) {
1427 if ( tmp_hx->ngb_elem_id [ k ] == elem1 ) {
1428 tmp_hx->quad [ k ] = mesh_quad;
1429 break;
1430 }
1431 }
1432 }
1433
1434 mesh_quad++;
1435 }
1436 }
1437 }
1438
1439 refine_nodes = 0;
1440 refine_nodes += mesh_edges * level * 2;
1441 refine_nodes += mesh_faces * ( level * 3 + level * level * 3 );
1442 refine_nodes += mesh_quads * ( level * 4 + level * level * 4 );
1443 refine_nodes += fe_faces * ( level * 3 + level * level * 3 );
1444 refine_nodes += fe_quads * ( level * 4 + level * level * 4 );
1445 refine_nodes += fe_tetras * ( level * 4 + level * level * 6 + level * level * level * 4 );
1446 refine_nodes += fe_hexas * ( level * 6 + level * level * 12 + level * level * level * 8 );
1447
1448 refine_node_id = fe_nodes + fine_nodes;
1449
1450 /* refine mesh edges */
1451
1452 mesh_edge = mesh_edge_array;
1453 for ( i = 0; i < mesh_edges; i++, mesh_edge++ ) {
1454 for ( j = 0; j < 2; j++ ) {
1455 if ( ( mesh_edge->fine_id [ j ] = ( int * ) calloc( level + 2, sizeof( int ) ) ) == NULL ) {
1456 error_message("Memory allocation error");
1457 }
1458
1459 mesh_edge->fine_id [ j ] [ 0 ] = mesh_edge->node [ j ]->id;
1460 for ( k = 1; k < level + 1; k++ ) {
1461 mesh_edge->fine_id [ j ] [ k ] = ++refine_node_id;
1462 }
1463
1464 mesh_edge->fine_id [ j ] [ level + 1 ] = mesh_edge->mid_node->id;
1465 }
1466 }
1467
1468 /*
1469 * o-------o-------o o
1470 | | | / \
1471 | | | / \
1472 | | | / \
1473 * o-------o-------o o o
1474 | | | / \ / \
1475 | cen | / o \
1476 | | | / | \
1477 * o-start-o--end--o o-------o-------o
1478 */
1479
1480 for ( j = 0; j < 4; j++ ) {
1481 if ( ( tmp_array_start [ j ] = ( long * ) calloc( level + 2, sizeof( long ) ) ) == NULL ) {
1482 error_message("Memory allocation error");
1483 }
1484
1485 if ( ( tmp_array_end [ j ] = ( long * ) calloc( level + 2, sizeof( long ) ) ) == NULL ) {
1486 error_message("Memory allocation error");
1487 }
1488
1489 if ( ( tmp_array_cen [ j ] = ( long * ) calloc( level + 2, sizeof( long ) ) ) == NULL ) {
1490 error_message("Memory allocation error");
1491 }
1492 }
1493
1494 /* refine mesh faces */
1495
1496 mesh_face = mesh_face_array;
1497 for ( i = 0; i < mesh_faces; i++, mesh_face++ ) {
1498 for ( j = 0; j < 3; j++ ) {
1499 j1 = face_ed_nd [ j ] [ 0 ];
1500 j2 = face_ed_nd [ j ] [ 1 ];
1501
1502 node1 = mesh_face->node [ j1 ]->id;
1503 node2 = mesh_face->node [ j2 ]->id;
1504
1505 /* I do rely on the fact that mesh_edge starts with node with smaller id */
1506
1507 if ( node1 < node2 ) {
1508 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1509 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1510 if ( mesh_edge->node [ 1 ]->id == node2 ) {
1511 for ( m = 0; m < level + 2; m++ ) {
1512 tmp_array_start [ j ] [ m ] = mesh_edge->fine_id [ 0 ] [ m ];
1513 tmp_array_end [ j ] [ m ] = mesh_edge->fine_id [ 1 ] [ m ];
1514 }
1515
1516 break;
1517 }
1518 }
1519 } else {
1520 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1521 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1522 if ( mesh_edge->node [ 1 ]->id == node1 ) {
1523 for ( m = 0; m < level + 2; m++ ) {
1524 tmp_array_start [ j ] [ m ] = mesh_edge->fine_id [ 1 ] [ m ];
1525 tmp_array_end [ j ] [ m ] = mesh_edge->fine_id [ 0 ] [ m ];
1526 }
1527
1528 break;
1529 }
1530 }
1531 }
1532
1533 tmp_array_cen [ j ] [ 0 ] = mesh_edge->mid_node->id;
1534 for ( m = 1; m < level + 1; m++ ) {
1535 tmp_array_cen [ j ] [ m ] = ++refine_node_id;
1536 }
1537
1538 tmp_array_cen [ j ] [ level + 1 ] = mesh_face->mid_node->id;
1539 }
1540
1541 for ( j = 0; j < 3; j++ ) {
1542 if ( ( mesh_face->fine_id [ j ] = ( int * ) calloc( ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
1543 error_message("Memory allocation error");
1544 }
1545
1546 if ( ( jj = j - 1 ) < 0 ) {
1547 jj = 2;
1548 }
1549
1550 j1 = face_ed_nd [ jj ] [ 0 ];
1551 j2 = face_ed_nd [ jj ] [ 1 ];
1552
1553 pos = 0;
1554 for ( k = 0; k < level + 2; k++ ) {
1555 mesh_face->fine_id [ j ] [ pos++ ] = tmp_array_start [ j2 ] [ k ];
1556 }
1557
1558 for ( k = 1; k < level + 1; k++ ) {
1559 mesh_face->fine_id [ j ] [ pos++ ] = tmp_array_end [ j1 ] [ k ];
1560 for ( m = 1; m < level + 1; m++ ) {
1561 mesh_face->fine_id [ j ] [ pos++ ] = ++refine_node_id;
1562 }
1563
1564 mesh_face->fine_id [ j ] [ pos++ ] = tmp_array_cen [ j2 ] [ k ];
1565 }
1566
1567 for ( k = 0; k < level + 2; k++ ) {
1568 mesh_face->fine_id [ j ] [ pos++ ] = tmp_array_cen [ j1 ] [ k ];
1569 }
1570 }
1571 }
1572
1573 /* refine mesh quads */
1574
1575 mesh_quad = mesh_quad_array;
1576 for ( i = 0; i < mesh_quads; i++, mesh_quad++ ) {
1577 for ( j = 0; j < 4; j++ ) {
1578 j1 = quad_ed_nd [ j ] [ 0 ];
1579 j2 = quad_ed_nd [ j ] [ 1 ];
1580
1581 node1 = mesh_quad->node [ j1 ]->id;
1582 node2 = mesh_quad->node [ j2 ]->id;
1583
1584 /* I do rely on the fact that mesh_edge starts with node with smaller id */
1585
1586 if ( node1 < node2 ) {
1587 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1588 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1589 if ( mesh_edge->node [ 1 ]->id == node2 ) {
1590 for ( m = 0; m < level + 2; m++ ) {
1591 tmp_array_start [ j ] [ m ] = mesh_edge->fine_id [ 0 ] [ m ];
1592 tmp_array_end [ j ] [ m ] = mesh_edge->fine_id [ 1 ] [ m ];
1593 }
1594
1595 break;
1596 }
1597 }
1598 } else {
1599 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1600 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1601 if ( mesh_edge->node [ 1 ]->id == node1 ) {
1602 for ( m = 0; m < level + 2; m++ ) {
1603 tmp_array_start [ j ] [ m ] = mesh_edge->fine_id [ 1 ] [ m ];
1604 tmp_array_end [ j ] [ m ] = mesh_edge->fine_id [ 0 ] [ m ];
1605 }
1606
1607 break;
1608 }
1609 }
1610 }
1611
1612 tmp_array_cen [ j ] [ 0 ] = mesh_edge->mid_node->id;
1613 for ( m = 1; m < level + 1; m++ ) {
1614 tmp_array_cen [ j ] [ m ] = ++refine_node_id;
1615 }
1616
1617 tmp_array_cen [ j ] [ level + 1 ] = mesh_quad->mid_node->id;
1618 }
1619
1620 for ( j = 0; j < 4; j++ ) {
1621 if ( ( mesh_quad->fine_id [ j ] = ( int * ) calloc( ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
1622 error_message("Memory allocation error");
1623 }
1624
1625 if ( ( jj = j - 1 ) < 0 ) {
1626 jj = 3;
1627 }
1628
1629 j1 = quad_ed_nd [ jj ] [ 0 ];
1630 j2 = quad_ed_nd [ jj ] [ 1 ];
1631
1632 pos = 0;
1633 for ( k = 0; k < level + 2; k++ ) {
1634 mesh_quad->fine_id [ j ] [ pos++ ] = tmp_array_start [ j2 ] [ k ];
1635 }
1636
1637 for ( k = 1; k < level + 1; k++ ) {
1638 mesh_quad->fine_id [ j ] [ pos++ ] = tmp_array_end [ j1 ] [ k ];
1639 for ( m = 1; m < level + 1; m++ ) {
1640 mesh_quad->fine_id [ j ] [ pos++ ] = ++refine_node_id;
1641 }
1642
1643 mesh_quad->fine_id [ j ] [ pos++ ] = tmp_array_cen [ j2 ] [ k ];
1644 }
1645
1646 for ( k = 0; k < level + 2; k++ ) {
1647 mesh_quad->fine_id [ j ] [ pos++ ] = tmp_array_cen [ j1 ] [ k ];
1648 }
1649 }
1650 }
1651
1652#ifdef COMPLETE_DATA_STRUCTURE
1653 if ( node_num_nodes != NULL ) {
1654 free(node_num_nodes);
1655 }
1656
1657 if ( node_con_nodes != NULL ) {
1658 free(node_con_nodes);
1659 }
1660
1661#endif
1662
1663 /* create fine edges */
1664
1665 /* this is really not necessary because everything is stored in mesh_edge;
1666 * however, mesh_edge_array is going to be freed */
1667
1668 fine_edge_id = 0;
1669
1670 mesh_edge = mesh_edge_array;
1671 for ( i = 0; i < fe_edges; i++, mesh_edge++ ) {
1672 for ( j = 0; j < 2; j++ ) {
1673 fine_edge = & ( fine_edge_array [ fine_edge_id++ ] );
1674 if ( ( fine_edge->fine_id = ( int * ) calloc( ( level + 2 ), sizeof( int ) ) ) == NULL ) {
1675 error_message("Memory allocation error");
1676 }
1677
1678 for ( k = 0; k < level + 2; k++ ) {
1679 fine_edge->fine_id [ k ] = mesh_edge->fine_id [ j ] [ k ];
1680 }
1681 }
1682 }
1683
1684 /* create fine quads */
1685
1686 fine_quad_id = 0;
1687
1688 fe_face = fe_face_array;
1689 for ( i = 0; i < fe_faces; i++, fe_face++, fine_node++ ) {
1690 fine_node->id = ++fine_node_id;
1691
1692#ifdef COMPLETE_DATA_STRUCTURE
1693 tmp_face = & tmp_face_array [ i ];
1694#endif
1695
1696 for ( j = 0; j < 3; j++ ) {
1697 j1 = face_ed_nd [ j ] [ 0 ];
1698 j2 = face_ed_nd [ j ] [ 1 ];
1699
1700 node1 = fe_face->node [ j1 ]->id;
1701 node2 = fe_face->node [ j2 ]->id;
1702
1703#ifdef COMPLETE_DATA_STRUCTURE
1704 mesh_edge = tmp_face->edge [ j ];
1705#endif
1706
1707 /* I do rely on the fact that mesh_edge starts with node with smaller id */
1708
1709 if ( node1 < node2 ) {
1710#ifndef COMPLETE_DATA_STRUCTURE
1711 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1712 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1713 if ( mesh_edge->node [ 1 ]->id == node2 ) {
1714 break;
1715 }
1716 }
1717
1718#endif
1719
1720 for ( k = 0; k < level + 2; k++ ) {
1721 tmp_array_start [ j ] [ k ] = mesh_edge->fine_id [ 0 ] [ k ];
1722 tmp_array_end [ j ] [ k ] = mesh_edge->fine_id [ 1 ] [ k ];
1723 }
1724 } else {
1725#ifndef COMPLETE_DATA_STRUCTURE
1726 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1727 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1728 if ( mesh_edge->node [ 1 ]->id == node1 ) {
1729 break;
1730 }
1731 }
1732
1733#endif
1734
1735 for ( k = 0; k < level + 2; k++ ) {
1736 tmp_array_start [ j ] [ k ] = mesh_edge->fine_id [ 1 ] [ k ];
1737 tmp_array_end [ j ] [ k ] = mesh_edge->fine_id [ 0 ] [ k ];
1738 }
1739 }
1740
1741 tmp_array_cen [ j ] [ 0 ] = mesh_edge->mid_node->id;
1742 for ( k = 1; k < level + 1; k++ ) {
1743 tmp_array_cen [ j ] [ k ] = ++refine_node_id;
1744 }
1745
1746 tmp_array_cen [ j ] [ level + 1 ] = fine_node->id;
1747 }
1748
1749 for ( j = 0; j < 3; j++ ) {
1750 fine_quad = & ( fine_quad_array [ fine_quad_id++ ] );
1751 if ( ( fine_quad->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
1752 error_message("Memory allocation error");
1753 }
1754
1755 if ( ( jj = j - 1 ) < 0 ) {
1756 jj = 2;
1757 }
1758
1759 j1 = face_ed_nd [ jj ] [ 0 ];
1760 j2 = face_ed_nd [ jj ] [ 1 ];
1761
1762 pos = 0;
1763 for ( k = 0; k < level + 2; k++ ) {
1764 fine_quad->fine_id [ pos++ ] = tmp_array_start [ j2 ] [ k ];
1765 }
1766
1767 for ( k = 1; k < level + 1; k++ ) {
1768 fine_quad->fine_id [ pos++ ] = tmp_array_end [ j1 ] [ k ];
1769 for ( m = 1; m < level + 1; m++ ) {
1770 fine_quad->fine_id [ pos++ ] = ++refine_node_id;
1771 }
1772
1773 fine_quad->fine_id [ pos++ ] = tmp_array_cen [ j2 ] [ k ];
1774 }
1775
1776 for ( k = 0; k < level + 2; k++ ) {
1777 fine_quad->fine_id [ pos++ ] = tmp_array_cen [ j1 ] [ k ];
1778 }
1779 }
1780 }
1781
1782 fe_quad = fe_quad_array;
1783 for ( i = 0; i < fe_quads; i++, fe_quad++, fine_node++ ) {
1784 fine_node->id = ++fine_node_id;
1785
1786#ifdef COMPLETE_DATA_STRUCTURE
1787 tmp_quad = & ( tmp_quad_array [ i ] );
1788#endif
1789
1790 for ( j = 0; j < 4; j++ ) {
1791 j1 = quad_ed_nd [ j ] [ 0 ];
1792 j2 = quad_ed_nd [ j ] [ 1 ];
1793
1794 node1 = fe_quad->node [ j1 ]->id;
1795 node2 = fe_quad->node [ j2 ]->id;
1796
1797#ifdef COMPLETE_DATA_STRUCTURE
1798 mesh_edge = tmp_quad->edge [ j ];
1799#endif
1800
1801 /* I do rely on the fact that mesh_edge starts with node with smaller id */
1802
1803 if ( node1 < node2 ) {
1804#ifndef COMPLETE_DATA_STRUCTURE
1805 for ( k = node_num_nodes [ node1 - 1 ]; k < node_num_nodes [ node1 ]; k++ ) {
1806 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1807 if ( mesh_edge->node [ 1 ]->id == node2 ) {
1808 break;
1809 }
1810 }
1811
1812#endif
1813
1814 for ( k = 0; k < level + 2; k++ ) {
1815 tmp_array_start [ j ] [ k ] = mesh_edge->fine_id [ 0 ] [ k ];
1816 tmp_array_end [ j ] [ k ] = mesh_edge->fine_id [ 1 ] [ k ];
1817 }
1818 } else {
1819#ifndef COMPLETE_DATA_STRUCTURE
1820 for ( k = node_num_nodes [ node2 - 1 ]; k < node_num_nodes [ node2 ]; k++ ) {
1821 mesh_edge = & ( mesh_edge_array [ node_con_nodes [ k ] - 1 ] );
1822 if ( mesh_edge->node [ 1 ]->id == node1 ) {
1823 break;
1824 }
1825 }
1826
1827#endif
1828
1829 for ( k = 0; k < level + 2; k++ ) {
1830 tmp_array_start [ j ] [ k ] = mesh_edge->fine_id [ 1 ] [ k ];
1831 tmp_array_end [ j ] [ k ] = mesh_edge->fine_id [ 0 ] [ k ];
1832 }
1833 }
1834
1835 tmp_array_cen [ j ] [ 0 ] = mesh_edge->mid_node->id;
1836 for ( k = 1; k < level + 1; k++ ) {
1837 tmp_array_cen [ j ] [ k ] = ++refine_node_id;
1838 }
1839
1840 tmp_array_cen [ j ] [ level + 1 ] = fine_node->id;
1841 }
1842
1843 for ( j = 0; j < 4; j++ ) {
1844 fine_quad = & ( fine_quad_array [ fine_quad_id++ ] );
1845 if ( ( fine_quad->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
1846 error_message("Memory allocation error");
1847 }
1848
1849 if ( ( jj = j - 1 ) < 0 ) {
1850 jj = 3;
1851 }
1852
1853 j1 = quad_ed_nd [ jj ] [ 0 ];
1854 j2 = quad_ed_nd [ jj ] [ 1 ];
1855
1856 pos = 0;
1857 for ( k = 0; k < level + 2; k++ ) {
1858 fine_quad->fine_id [ pos++ ] = tmp_array_start [ j2 ] [ k ];
1859 }
1860
1861 for ( k = 1; k < level + 1; k++ ) {
1862 fine_quad->fine_id [ pos++ ] = tmp_array_end [ j1 ] [ k ];
1863 for ( m = 1; m < level + 1; m++ ) {
1864 fine_quad->fine_id [ pos++ ] = ++refine_node_id;
1865 }
1866
1867 fine_quad->fine_id [ pos++ ] = tmp_array_cen [ j2 ] [ k ];
1868 }
1869
1870 for ( k = 0; k < level + 2; k++ ) {
1871 fine_quad->fine_id [ pos++ ] = tmp_array_cen [ j1 ] [ k ];
1872 }
1873 }
1874 }
1875
1876 for ( j = 0; j < 4; j++ ) {
1877 free(tmp_array_start [ j ]);
1878 free(tmp_array_end [ j ]);
1879 free(tmp_array_cen [ j ]);
1880 }
1881
1882#ifdef COMPLETE_DATA_STRUCTURE
1883 if ( tmp_face_array != NULL ) {
1884 free(tmp_face_array);
1885 }
1886
1887 if ( tmp_quad_array != NULL ) {
1888 free(tmp_quad_array);
1889 }
1890
1891#else
1892 if ( node_num_nodes != NULL ) {
1893 free(node_num_nodes);
1894 }
1895
1896 if ( node_con_nodes != NULL ) {
1897 free(node_con_nodes);
1898 }
1899
1900#endif
1901
1902 mesh_edge = mesh_edge_array;
1903 for ( i = 0; i < mesh_edges; i++, mesh_edge++ ) {
1904 for ( j = 0; j < 2; j++ ) {
1905 if ( mesh_edge->fine_id [ j ] != NULL ) {
1906 free(mesh_edge->fine_id [ j ]);
1907 }
1908 }
1909 }
1910
1911 if ( mesh_edge_array != NULL ) {
1912 free(mesh_edge_array);
1913 }
1914
1915 /* create fine hexas */
1916
1917 for ( j = 0; j < 6; j++ ) {
1918 if ( ( tmp_array_cen [ j ] = ( long * ) calloc( level + 2, sizeof( long ) ) ) == NULL ) {
1919 error_message("Memory allocation error");
1920 }
1921 }
1922
1923 for ( j = 0; j < 12; j++ ) {
1924 if ( ( tmp_array [ j ] = ( long * ) calloc( ( level + 2 ) * ( level + 2 ), sizeof( long ) ) ) == NULL ) {
1925 error_message("Memory allocation error");
1926 }
1927 }
1928
1929 fine_hexa_id = 0;
1930
1931 fe_tetra = fe_tetra_array;
1932 tmp_tetra = tmp_tetra_array;
1933 for ( i = 0; i < fe_tetras; i++, fe_tetra++, tmp_tetra++, fine_node++ ) {
1934 fine_node->id = ++fine_node_id;
1935
1936 for ( j = 0; j < 4; j++ ) {
1937 tmp_array_cen [ j ] [ 0 ] = tmp_tetra->face [ j ]->mid_node->id;
1938 for ( k = 1; k < level + 1; k++ ) {
1939 tmp_array_cen [ j ] [ k ] = ++refine_node_id;
1940 }
1941
1942 tmp_array_cen [ j ] [ level + 1 ] = fine_node->id;
1943 }
1944
1945 /* mark swapped faces (see ordering of nodes on tetra faces above);
1946 * bottom face is swapped if normal is outer;
1947 * side faces are swapped if normal is inner */
1948
1949 for ( j = 0; j < 4; j++ ) {
1950 mesh_face = mesh_fc [ j ] = tmp_tetra->face [ j ];
1951 for ( k = 0; k < 3; k++ ) {
1952 if ( mesh_face->node [ k ] == fe_tetra->node [ tetra_fc_nd [ j ] [ 0 ] ] ) {
1953 if ( ( kk = k + 1 ) == 3 ) {
1954 kk = 0;
1955 }
1956
1957 if ( mesh_face->node [ kk ] == fe_tetra->node [ tetra_fc_nd [ j ] [ 1 ] ] ) {
1958 swap [ j ] = 0;
1959 } else {
1960 swap [ j ] = 1;
1961 }
1962
1963 break;
1964 }
1965 }
1966 }
1967
1968 /* form inner tmp fine quads "perpendicular" to bottom face;
1969 * each one is oriented u: from center of bottom edge to center of bottom face
1970 * v: from center of bottom edge to center of side face */
1971
1972 for ( j = 0; j < 3; j++ ) {
1973 pos = 0;
1974
1975 p = ( level + 2 ) * ( level + 1 );
1976 mesh_face = mesh_fc [ 0 ];
1977 #pragma GCC diagnostic push
1978 // swap[0] possibly uninitialized
1979 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
1980 for ( k = 0; k < 3; k++ ) {
1981 if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
1982 if ( swap [ 0 ] == 1 ) {
1983 kk = k;
1984 } else {
1985 if ( ( kk = k + 1 ) == 3 ) {
1986 kk = 0;
1987 }
1988 }
1989
1990 break;
1991 }
1992 }
1993 #pragma GCC diagnostic pop
1994
1995 for ( m = 0; m < level + 2; m++ ) {
1996 tmp_array [ j ] [ pos++ ] = mesh_face->fine_id [ kk ] [ p++ ];
1997 }
1998
1999 if ( level != 0 ) {
2000 p = ( level + 2 ) * ( level + 1 );
2001 mesh_face = mesh_fc [ j + 1 ];
2002 for ( k = 0; k < 3; k++ ) {
2003 if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2004 if ( swap [ j + 1 ] == 1 ) {
2005 kk = k;
2006 } else {
2007 if ( ( kk = k + 1 ) == 3 ) {
2008 kk = 0;
2009 }
2010 }
2011
2012 break;
2013 }
2014 }
2015
2016 for ( n = 1; n < level + 1; n++ ) {
2017 tmp_array [ j ] [ pos++ ] = mesh_face->fine_id [ kk ] [ ++p ];
2018 for ( m = 0; m < level; m++ ) {
2019 tmp_array [ j ] [ pos++ ] = ++refine_node_id;
2020 }
2021
2022 tmp_array [ j ] [ pos++ ] = tmp_array_cen [ 0 ] [ n ];
2023 }
2024 }
2025
2026 for ( m = 0; m < level + 2; m++ ) {
2027 tmp_array [ j ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ m ];
2028 }
2029 }
2030
2031 /* form inner tmp fine quads "parallel" to bottom face;
2032 * each one is oriented as fine quads on bottom face if not swapped (this is 0-1-2) */
2033
2034 for ( j = 0; j < 3; j++ ) {
2035 j1 = face_ed_nd [ j ] [ 0 ];
2036 j2 = face_ed_nd [ j ] [ 1 ];
2037
2038 pos = 0;
2039
2040 p = ( level + 2 ) * ( level + 1 );
2041 mesh_face = mesh_fc [ j + 1 ];
2042 for ( k = 0; k < 3; k++ ) {
2043 if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2044 if ( swap [ j + 1 ] == 1 ) {
2045 kk = k;
2046 } else {
2047 if ( ( kk = k + 1 ) == 3 ) {
2048 kk = 0;
2049 }
2050 }
2051
2052 break;
2053 }
2054 }
2055
2056 for ( m = 0; m < level + 2; m++ ) {
2057 tmp_array [ j + 3 ] [ pos++ ] = mesh_face->fine_id [ kk ] [ p++ ];
2058 }
2059
2060 if ( ( jj = j ) == 0 ) {
2061 jj = 3;
2062 }
2063
2064 if ( level != 0 ) {
2065 p = ( level + 2 ) * ( level + 1 );
2066 mesh_face = mesh_fc [ jj ];
2067 for ( k = 0; k < 3; k++ ) {
2068 if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2069 if ( swap [ jj ] == 1 ) {
2070 kk = k;
2071 } else {
2072 if ( ( kk = k + 1 ) == 3 ) {
2073 kk = 0;
2074 }
2075 }
2076
2077 break;
2078 }
2079 }
2080
2081 for ( n = 1; n < level + 1; n++ ) {
2082 tmp_array [ j + 3 ] [ pos++ ] = mesh_face->fine_id [ kk ] [ ++p ];
2083 for ( m = 0; m < level; m++ ) {
2084 tmp_array [ j + 3 ] [ pos++ ] = ++refine_node_id;
2085 }
2086
2087 tmp_array [ j + 3 ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ n ];
2088 }
2089 }
2090
2091 for ( m = 0; m < level + 2; m++ ) {
2092 tmp_array [ j + 3 ] [ pos++ ] = tmp_array_cen [ jj ] [ m ];
2093 }
2094 }
2095
2096 /* form bottom fine hexas */
2097
2098 /* filling 3d matrix in a sequential way would be too difficult
2099 * with respect to switching among boundary and inner sources;
2100 * therefore it is much easier to traverse boundary sources
2101 * and localize them into 3d matrix;
2102 * this is done using macros matrix_3d and possible matrix_2d;
2103 * when traversing whole boudary source it is better to use direct
2104 * access instead of macro matrix_2d */
2105
2106 /* note: I ignore that some entries of 3d matrix are rewritten several times;
2107 * in order to enable direct access */
2108
2109 for ( j = 0; j < 3; j++ ) {
2110 fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2111 if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2112 error_message("Memory allocation error");
2113 }
2114
2115 mesh_face = mesh_fc [ 0 ];
2116 for ( k = 0; k < 3; k++ ) {
2117 if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2118 break;
2119 }
2120 }
2121
2122 if ( swap [ 0 ] == 1 ) {
2123 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2124 for ( m = 0; m < level + 2; m++, pos++ ) {
2125 matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2126 }
2127 }
2128 } else {
2129 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2130 for ( m = 0; m < level + 2; m++, pos++ ) {
2131 matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2132 }
2133 }
2134 }
2135
2136 mesh_face = mesh_fc [ j + 1 ];
2137 for ( k = 0; k < 3; k++ ) {
2138 if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2139 break;
2140 }
2141 }
2142
2143 if ( swap [ j + 1 ] == 1 ) {
2144 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2145 for ( m = 0; m < level + 2; m++, pos++ ) {
2146 matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2147 }
2148 }
2149 } else {
2150 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2151 for ( m = 0; m < level + 2; m++, pos++ ) {
2152 matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2153 }
2154 }
2155 }
2156
2157 if ( ( jj = j ) == 0 ) {
2158 jj = 3;
2159 }
2160
2161 mesh_face = mesh_fc [ jj ];
2162 for ( k = 0; k < 3; k++ ) {
2163 if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2164 break;
2165 }
2166 }
2167
2168 if ( swap [ jj ] == 1 ) {
2169 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2170 for ( m = 0; m < level + 2; m++, pos++ ) {
2171 matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2172 }
2173 }
2174 } else {
2175 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2176 for ( m = 0; m < level + 2; m++, pos++ ) {
2177 matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2178 }
2179 }
2180 }
2181
2182 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2183 for ( m = 0; m < level + 2; m++, pos++ ) {
2184 matrix_3d(fine_hexa->fine_id, level + 1, m, n) = tmp_array [ j ] [ pos ]; /* macro */
2185 matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ jj - 1 ] [ pos ]; /* macro */
2186 matrix_3d(fine_hexa->fine_id, m, n, level + 1) = tmp_array [ j + 3 ] [ pos ]; /* macro */
2187 }
2188 }
2189
2190 if ( level != 0 ) {
2191 for ( k = 1; k < level + 1; k++ ) {
2192 for ( n = 1; n < level + 1; n++ ) {
2193 for ( m = 1; m < level + 1; m++ ) {
2194 matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2195 }
2196 }
2197 }
2198 }
2199 }
2200
2201 /* form top hexa */
2202
2203 fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2204 if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2205 error_message("Memory allocation error");
2206 }
2207
2208 mesh_face = mesh_fc [ 1 ];
2209 for ( k = 0; k < 3; k++ ) {
2210 if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2211 break;
2212 }
2213 }
2214
2215 if ( swap [ 1 ] == 1 ) {
2216 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2217 for ( m = 0; m < level + 2; m++, pos++ ) {
2218 matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2219 }
2220 }
2221 } else {
2222 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2223 for ( m = 0; m < level + 2; m++, pos++ ) {
2224 matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2225 }
2226 }
2227 }
2228
2229 mesh_face = mesh_fc [ 2 ];
2230 for ( k = 0; k < 3; k++ ) {
2231 if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2232 break;
2233 }
2234 }
2235
2236 if ( swap [ 2 ] == 1 ) {
2237 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2238 for ( m = 0; m < level + 2; m++, pos++ ) {
2239 matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2240 }
2241 }
2242 } else {
2243 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2244 for ( m = 0; m < level + 2; m++, pos++ ) {
2245 matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2246 }
2247 }
2248 }
2249
2250 mesh_face = mesh_fc [ 3 ];
2251 for ( k = 0; k < 3; k++ ) {
2252 if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2253 break;
2254 }
2255 }
2256
2257 if ( swap [ 3 ] == 1 ) {
2258 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2259 for ( m = 0; m < level + 2; m++, pos++ ) {
2260 matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2261 }
2262 }
2263 } else {
2264 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2265 for ( m = 0; m < level + 2; m++, pos++ ) {
2266 matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2267 }
2268 }
2269 }
2270
2271 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2272 for ( m = 0; m < level + 2; m++, pos++ ) {
2273 matrix_3d(fine_hexa->fine_id, level + 1, n, m) = tmp_array [ 0 + 3 ] [ pos ]; /* macro */
2274 matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ 2 + 3 ] [ pos ]; /* macro */
2275 matrix_3d(fine_hexa->fine_id, n, m, level + 1) = tmp_array [ 1 + 3 ] [ pos ]; /* macro */
2276 }
2277 }
2278
2279 if ( level != 0 ) {
2280 for ( k = 1; k < level + 1; k++ ) {
2281 for ( n = 1; n < level + 1; n++ ) {
2282 for ( m = 1; m < level + 1; m++ ) {
2283 matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2284 }
2285 }
2286 }
2287 }
2288 }
2289
2290 fe_hexa = fe_hexa_array;
2291 tmp_hexa = tmp_hexa_array;
2292 for ( i = 0; i < fe_hexas; i++, fe_hexa++, tmp_hexa++, fine_node++ ) {
2293 fine_node->id = ++fine_node_id;
2294
2295 for ( j = 0; j < 6; j++ ) {
2296 tmp_array_cen [ j ] [ 0 ] = tmp_hexa->quad [ j ]->mid_node->id;
2297 for ( k = 1; k < level + 1; k++ ) {
2298 tmp_array_cen [ j ] [ k ] = ++refine_node_id;
2299 }
2300
2301 tmp_array_cen [ j ] [ level + 1 ] = fine_node->id;
2302 }
2303
2304 /* mark swapped faces (see ordering of nodes on hexa faces above);
2305 * bottom and top faces are swapped if normal is outer;
2306 * side faces are swapped if normal is inner */
2307
2308 for ( j = 0; j < 6; j++ ) {
2309 mesh_quad = mesh_qd [ j ] = tmp_hexa->quad [ j ];
2310 for ( k = 0; k < 4; k++ ) {
2311 if ( mesh_quad->node [ k ] == fe_hexa->node [ hexa_fc_nd [ j ] [ 0 ] ] ) {
2312 if ( ( kk = k + 1 ) == 4 ) {
2313 kk = 0;
2314 }
2315
2316 if ( mesh_quad->node [ kk ] == fe_hexa->node [ hexa_fc_nd [ j ] [ 1 ] ] ) {
2317 swap [ j ] = 0;
2318 } else {
2319 swap [ j ] = 1;
2320 }
2321
2322 break;
2323 }
2324 }
2325 }
2326
2327 /* form inner bottom tmp fine quads "perpendicular" to bottom face;
2328 * each one is oriented u: from center of bottom edge to center of bottom face
2329 * v: from center of bottom edge to center of side face */
2330
2331 for ( j = 0; j < 4; j++ ) {
2332 pos = 0;
2333
2334 p = ( level + 2 ) * ( level + 1 );
2335 mesh_quad = mesh_qd [ 0 ];
2336 #pragma GCC diagnostic push
2337 // swap[0] possibly uninitialized
2338 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
2339 for ( k = 0; k < 4; k++ ) {
2340 if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2341 if ( swap [ 0 ] == 1 ) {
2342 kk = k;
2343 } else {
2344 if ( ( kk = k + 1 ) == 4 ) {
2345 kk = 0;
2346 }
2347 }
2348
2349 break;
2350 }
2351 }
2352 #pragma GCC diagnostic pop
2353
2354 for ( m = 0; m < level + 2; m++ ) {
2355 tmp_array [ j ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ p++ ];
2356 }
2357
2358 if ( level != 0 ) {
2359 p = ( level + 2 ) * ( level + 1 );
2360 mesh_quad = mesh_qd [ j + 1 ];
2361 for ( k = 0; k < 4; k++ ) {
2362 if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2363 if ( swap [ j + 1 ] == 1 ) {
2364 kk = k;
2365 } else {
2366 if ( ( kk = k + 1 ) == 4 ) {
2367 kk = 0;
2368 }
2369 }
2370
2371 break;
2372 }
2373 }
2374
2375 for ( n = 1; n < level + 1; n++ ) {
2376 tmp_array [ j ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ ++p ];
2377 for ( m = 0; m < level; m++ ) {
2378 tmp_array [ j ] [ pos++ ] = ++refine_node_id;
2379 }
2380
2381 tmp_array [ j ] [ pos++ ] = tmp_array_cen [ 0 ] [ n ];
2382 }
2383 }
2384
2385 for ( m = 0; m < level + 2; m++ ) {
2386 tmp_array [ j ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ m ];
2387 }
2388 }
2389
2390 /* form top tmp fine quads "perpendicular" to top face;
2391 * each one is oriented u: from center of top edge to center of top face
2392 * v: from center of top edge to center of side face */
2393
2394 for ( j = 0; j < 4; j++ ) {
2395 pos = 0;
2396
2397 p = ( level + 2 ) * ( level + 1 );
2398 mesh_quad = mesh_qd [ 5 ];
2399 #pragma GCC diagnostic push
2400 // swap[5] possibly uninitialized
2401 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
2402 for ( k = 0; k < 4; k++ ) {
2403 if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2404 if ( swap [ 5 ] == 0 ) {
2405 kk = k;
2406 } else {
2407 if ( ( kk = k + 1 ) == 4 ) {
2408 kk = 0;
2409 }
2410 }
2411
2412 break;
2413 }
2414 }
2415 #pragma GCC diagnostic pop
2416
2417 for ( m = 0; m < level + 2; m++ ) {
2418 tmp_array [ j + 4 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ p++ ];
2419 }
2420
2421 if ( level != 0 ) {
2422 p = ( level + 2 ) * ( level + 1 );
2423 mesh_quad = mesh_qd [ j + 1 ];
2424 for ( k = 0; k < 4; k++ ) {
2425 if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2426 if ( swap [ j + 1 ] == 0 ) {
2427 kk = k;
2428 } else {
2429 if ( ( kk = k + 1 ) == 4 ) {
2430 kk = 0;
2431 }
2432 }
2433
2434 break;
2435 }
2436 }
2437
2438 for ( n = 1; n < level + 1; n++ ) {
2439 tmp_array [ j + 4 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ ++p ];
2440 for ( m = 0; m < level; m++ ) {
2441 tmp_array [ j + 4 ] [ pos++ ] = ++refine_node_id;
2442 }
2443
2444 tmp_array [ j + 4 ] [ pos++ ] = tmp_array_cen [ 5 ] [ n ];
2445 }
2446 }
2447
2448 for ( m = 0; m < level + 2; m++ ) {
2449 tmp_array [ j + 4 ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ m ];
2450 }
2451 }
2452
2453 /* form inner tmp fine quads "parallel" to bottom and top face;
2454 * each one is oriented as fine quads on bottom face if not swapped (this is 0-1-2-3) */
2455
2456 for ( j = 0; j < 4; j++ ) {
2457 pos = 0;
2458
2459 p = ( level + 2 ) * ( level + 1 );
2460 mesh_quad = mesh_qd [ j + 1 ];
2461 for ( k = 0; k < 4; k++ ) {
2462 if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2463 if ( swap [ j + 1 ] == 1 ) {
2464 kk = k;
2465 } else {
2466 if ( ( kk = k + 1 ) == 4 ) {
2467 kk = 0;
2468 }
2469 }
2470
2471 break;
2472 }
2473 }
2474
2475 for ( m = 0; m < level + 2; m++ ) {
2476 tmp_array [ j + 8 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ p++ ];
2477 }
2478
2479 if ( ( jj = j ) == 0 ) {
2480 jj = 4;
2481 }
2482
2483 if ( level != 0 ) {
2484 p = ( level + 2 ) * ( level + 1 );
2485 mesh_quad = mesh_qd [ jj ];
2486 for ( k = 0; k < 4; k++ ) {
2487 if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2488 if ( swap [ jj ] == 1 ) {
2489 kk = k;
2490 } else {
2491 if ( ( kk = k + 1 ) == 4 ) {
2492 kk = 0;
2493 }
2494 }
2495
2496 break;
2497 }
2498 }
2499
2500 for ( n = 1; n < level + 1; n++ ) {
2501 tmp_array [ j + 8 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ ++p ];
2502 for ( m = 0; m < level; m++ ) {
2503 tmp_array [ j + 8 ] [ pos++ ] = ++refine_node_id;
2504 }
2505
2506 tmp_array [ j + 8 ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ n ];
2507 }
2508 }
2509
2510 for ( m = 0; m < level + 2; m++ ) {
2511 tmp_array [ j + 8 ] [ pos++ ] = tmp_array_cen [ jj ] [ m ];
2512 }
2513 }
2514
2515 /* form bottom fine hexas */
2516
2517 /* filling 3d matrix in a sequential way would be too difficult
2518 * with respect to switching among boundary and inner sources;
2519 * therefore it is much easier to traverse boundary sources
2520 * and localize them into 3d matrix;
2521 * this is done using macros matrix_3d and possible matrix_2d;
2522 * when traversing whole boudary source it is better to use direct
2523 * access instead of macro matrix_2d */
2524
2525 /* note: I ignore that some entries of 3d matrix are rewritten several times;
2526 * in order to enable direct access */
2527
2528 for ( j = 0; j < 4; j++ ) {
2529 fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2530 if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2531 error_message("Memory allocation error");
2532 }
2533
2534 mesh_quad = mesh_qd [ 0 ];
2535 for ( k = 0; k < 4; k++ ) {
2536 if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2537 break;
2538 }
2539 }
2540
2541 if ( swap [ 0 ] == 1 ) {
2542 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2543 for ( m = 0; m < level + 2; m++, pos++ ) {
2544 matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2545 }
2546 }
2547 } else {
2548 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2549 for ( m = 0; m < level + 2; m++, pos++ ) {
2550 matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2551 }
2552 }
2553 }
2554
2555 mesh_quad = mesh_qd [ j + 1 ];
2556 for ( k = 0; k < 4; k++ ) {
2557 if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2558 break;
2559 }
2560 }
2561
2562 if ( swap [ j + 1 ] == 1 ) {
2563 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2564 for ( m = 0; m < level + 2; m++, pos++ ) {
2565 matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2566 }
2567 }
2568 } else {
2569 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2570 for ( m = 0; m < level + 2; m++, pos++ ) {
2571 matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2572 }
2573 }
2574 }
2575
2576 if ( ( jj = j ) == 0 ) {
2577 jj = 4;
2578 }
2579
2580 mesh_quad = mesh_qd [ jj ];
2581 for ( k = 0; k < 4; k++ ) {
2582 if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2583 break;
2584 }
2585 }
2586
2587 if ( swap [ jj ] == 1 ) {
2588 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2589 for ( m = 0; m < level + 2; m++, pos++ ) {
2590 matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2591 }
2592 }
2593 } else {
2594 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2595 for ( m = 0; m < level + 2; m++, pos++ ) {
2596 matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2597 }
2598 }
2599 }
2600
2601 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2602 for ( m = 0; m < level + 2; m++, pos++ ) {
2603 matrix_3d(fine_hexa->fine_id, level + 1, m, n) = tmp_array [ j ] [ pos ]; /* macro */
2604 matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ jj - 1 ] [ pos ]; /* macro */
2605 matrix_3d(fine_hexa->fine_id, m, n, level + 1) = tmp_array [ j + 8 ] [ pos ]; /* macro */
2606 }
2607 }
2608
2609 if ( level != 0 ) {
2610 for ( k = 1; k < level + 1; k++ ) {
2611 for ( n = 1; n < level + 1; n++ ) {
2612 for ( m = 1; m < level + 1; m++ ) {
2613 matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2614 }
2615 }
2616 }
2617 }
2618 }
2619
2620 /* form top fine hexas */
2621
2622 for ( j = 0; j < 4; j++ ) {
2623 fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2624 if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2625 error_message("Memory allocation error");
2626 }
2627
2628 mesh_quad = mesh_qd [ 5 ];
2629 for ( k = 0; k < 4; k++ ) {
2630 if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2631 break;
2632 }
2633 }
2634
2635 if ( swap [ 5 ] == 1 ) {
2636 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2637 for ( m = 0; m < level + 2; m++, pos++ ) {
2638 matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2639 }
2640 }
2641 } else {
2642 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2643 for ( m = 0; m < level + 2; m++, pos++ ) {
2644 matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2645 }
2646 }
2647 }
2648
2649 mesh_quad = mesh_qd [ j + 1 ];
2650 for ( k = 0; k < 4; k++ ) {
2651 if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2652 break;
2653 }
2654 }
2655
2656 if ( swap [ j + 1 ] == 1 ) {
2657 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2658 for ( m = 0; m < level + 2; m++, pos++ ) {
2659 matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2660 }
2661 }
2662 } else {
2663 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2664 for ( m = 0; m < level + 2; m++, pos++ ) {
2665 matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2666 }
2667 }
2668 }
2669
2670 if ( ( jj = j ) == 0 ) {
2671 jj = 4;
2672 }
2673
2674 mesh_quad = mesh_qd [ jj ];
2675 for ( k = 0; k < 4; k++ ) {
2676 if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2677 break;
2678 }
2679 }
2680
2681 if ( swap [ jj ] == 1 ) {
2682 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2683 for ( m = 0; m < level + 2; m++, pos++ ) {
2684 matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2685 }
2686 }
2687 } else {
2688 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2689 for ( m = 0; m < level + 2; m++, pos++ ) {
2690 matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2691 }
2692 }
2693 }
2694
2695 for ( pos = 0, n = 0; n < level + 2; n++ ) {
2696 for ( m = 0; m < level + 2; m++, pos++ ) {
2697 matrix_3d(fine_hexa->fine_id, level + 1, m, n) = tmp_array [ jj - 1 + 4 ] [ pos ]; /* macro */
2698 matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ j + 4 ] [ pos ]; /* macro */
2699 matrix_3d(fine_hexa->fine_id, n, m, level + 1) = tmp_array [ j + 8 ] [ pos ]; /* macro */
2700 }
2701 }
2702
2703 if ( level != 0 ) {
2704 for ( k = 1; k < level + 1; k++ ) {
2705 for ( n = 1; n < level + 1; n++ ) {
2706 for ( m = 1; m < level + 1; m++ ) {
2707 matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2708 }
2709 }
2710 }
2711 }
2712 }
2713 }
2714
2715 for ( j = 0; j < 6; j++ ) {
2716 free(tmp_array_cen [ j ]);
2717 free(tmp_array [ j ]);
2718 free(tmp_array [ j + 6 ]);
2719 }
2720
2721 if ( fe_node_array != NULL ) {
2722 free(fe_node_array);
2723 }
2724
2725 if ( fe_edge_array != NULL ) {
2726 free(fe_edge_array);
2727 }
2728
2729 if ( fe_face_array != NULL ) {
2730 free(fe_face_array);
2731 }
2732
2733 if ( fe_quad_array != NULL ) {
2734 free(fe_quad_array);
2735 }
2736
2737 if ( fe_tetra_array != NULL ) {
2738 free(fe_tetra_array);
2739 }
2740
2741 if ( fe_hexa_array != NULL ) {
2742 free(fe_hexa_array);
2743 }
2744
2745 mesh_face = mesh_face_array;
2746 for ( i = 0; i < mesh_faces; i++, mesh_face++ ) {
2747 for ( j = 0; j < 3; j++ ) {
2748 if ( mesh_face->fine_id [ j ] != NULL ) {
2749 free(mesh_face->fine_id [ j ]);
2750 }
2751 }
2752 }
2753
2754 mesh_quad = mesh_quad_array;
2755 for ( i = 0; i < mesh_quads; i++, mesh_quad++ ) {
2756 for ( j = 0; j < 4; j++ ) {
2757 if ( mesh_quad->fine_id [ j ] != NULL ) {
2758 free(mesh_quad->fine_id [ j ]);
2759 }
2760 }
2761 }
2762
2763 if ( mesh_face_array != NULL ) {
2764 free(mesh_face_array);
2765 }
2766
2767 if ( mesh_quad_array != NULL ) {
2768 free(mesh_quad_array);
2769 }
2770
2771 /*
2772 * if(refine_nodes != 0){
2773 * if((refine_node_array = (fe_node_rec *)calloc(refine_nodes, sizeof(fe_node_rec))) == NULL)
2774 * error_message("Memory allocation error");
2775 *
2776 * fine_edge = fine_edge_array;
2777 * for(i = 0; i < fine_edges; i++, fine_edge++){
2778 * pos = 0;
2779 * for(m = 0; m < level + 2; m++, pos++){
2780 * node_id = fine_edge -> fine_id[pos];
2781 * if(node_id <= fe_nodes + fine_nodes)continue;
2782 *
2783 * refine_node = &(refine_node_array[node_id - fe_nodes - fine_nodes - 1]);
2784 * refine_node -> id = node_id;
2785 * }
2786 * }
2787 *
2788 * fine_quad = fine_quad_array;
2789 * for(i = 0; i < fine_quads; i++, fine_quad++){
2790 * pos = 0;
2791 * for(n = 0; n < level + 2; n++){
2792 * for(m = 0; m < level + 2; m++, pos++){
2793 * node_id = fine_quad -> fine_id[pos];
2794 * if(node_id <= fe_nodes + fine_nodes)continue;
2795 *
2796 * refine_node = &(refine_node_array[node_id - fe_nodes - fine_nodes - 1]);
2797 * refine_node -> id = node_id;
2798 * }
2799 * }
2800 * }
2801 *
2802 * fine_hexa = fine_hexa_array;
2803 * for(i = 0; i < fine_hexas; i++, fine_hexa++){
2804 * pos = 0;
2805 * for(k = 0; k < level + 2; k++){
2806 * for(n = 0; n < level + 2; n++){
2807 * for(m = 0; m < level + 2; m++, pos++){
2808 * node_id = fine_hexa -> fine_id[pos];
2809 * if(node_id <= fe_nodes + fine_nodes)continue;
2810 *
2811 * refine_node = &(refine_node_array[node_id - fe_nodes - fine_nodes - 1]);
2812 * refine_node -> id = node_id;
2813 * }
2814 * }
2815 * }
2816 * }
2817 * }
2818 */
2819
2820 if ( fine_node_array != NULL ) {
2821 free(fine_node_array);
2822 }
2823
2824 /*
2825 * if(refine_node_array != NULL)free(refine_node_array);
2826 */
2827
2828 edge_id = face_id = quad_id = tetra_id = hexa_id = 0;
2829 for ( i = 0; i < fe_elems; i++ ) {
2830 element = d->giveElement(i + 1);
2831 RefinedElement &refinedElement = refinedElementList.at(i + 1);
2832 boundary = refinedElement.giveBoundaryFlagArray();
2833
2834 switch ( element->giveGeometryType() ) {
2835 case EGT_line_1:
2836 pos = edge_id * 2;
2837 fine_edge = & ( fine_edge_array [ pos ] );
2838 for ( j = 0; j < 2; j++, fine_edge++ ) {
2839 connectivity = refinedElement.giveFineNodeArray(j + 1);
2840 connectivity->resize(level + 2);
2841 for ( k = 0; k < level + 2; k++ ) {
2842 connectivity->at(k + 1) = fine_edge->fine_id [ k ];
2843 }
2844 }
2845
2846 for ( j = 0; j < 2; j++ ) {
2847 /* note: I cannot use fe_edge (it is already released)
2848 *
2849 * node = fe_edge -> node[j] -> id;
2850 */
2851
2852 node = ( & ( fine_edge_array [ pos + j ] ) )->fine_id [ 0 ];
2853 if ( node_num_elems [ node ] - node_num_elems [ node - 1 ] == 1 ) {
2854 boundary->at(j + 1) = 1;
2855 } else {
2856 boundary->at(j + 1) = 0;
2857 }
2858 }
2859
2860 edge_id++;
2861 break;
2862
2863 case EGT_triangle_1:
2864 pos = face_id * 3;
2865 fine_quad = & ( fine_quad_array [ pos ] );
2866 for ( j = 0; j < 3; j++, fine_quad++ ) {
2867 connectivity = refinedElement.giveFineNodeArray(j + 1);
2868 connectivity->resize( ( level + 2 ) * ( level + 2 ) );
2869 for ( k = 0; k < ( level + 2 ) * ( level + 2 ); k++ ) {
2870 connectivity->at(k + 1) = fine_quad->fine_id [ k ];
2871 }
2872 }
2873
2874 /* note: I am not using tmp_face (and its ngb_elem_id) because tmp_face exists
2875 * only if COMPLETE_DATA_STRUCTURE is required but that is supported for 2-manifolds only */
2876
2877 elem1 = global_face_id(face_id) + 1;
2878 for ( j = 0; j < 3; j++ ) {
2879 boundary->at(j + 1) = 1;
2880
2881 j1 = face_ed_nd [ j ] [ 0 ];
2882 j2 = face_ed_nd [ j ] [ 1 ];
2883
2884 /* note: I cannot use fe_face (it is already released)
2885 *
2886 * node1 = fe_face -> node[j1] -> id;
2887 * node2 = fe_face -> node[j2] -> id;
2888 */
2889
2890 node1 = ( & ( fine_quad_array [ pos + j1 ] ) )->fine_id [ 0 ];
2891 node2 = ( & ( fine_quad_array [ pos + j2 ] ) )->fine_id [ 0 ];
2892
2893 flag = 1;
2894 for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
2895 elem2 = node_con_elems [ k ];
2896 if ( elem2 == elem1 ) {
2897 continue;
2898 }
2899
2900 if ( is_edge(elem2) == 0 ) { /* macro */
2901 for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
2902 if ( elem2 == node_con_elems [ m ] ) {
2903 boundary->at(j + 1) = flag = 0;
2904 break;
2905 }
2906 }
2907
2908 if ( flag == 0 ) {
2909 break;
2910 }
2911 }
2912 }
2913 }
2914
2915 face_id++;
2916 break;
2917
2918 case EGT_quad_1:
2919 pos = fe_faces * 3 + quad_id * 4;
2920 fine_quad = & ( fine_quad_array [ pos ] );
2921 for ( j = 0; j < 4; j++, fine_quad++ ) {
2922 connectivity = refinedElement.giveFineNodeArray(j + 1);
2923 connectivity->resize( ( level + 2 ) * ( level + 2 ) );
2924 for ( k = 0; k < ( level + 2 ) * ( level + 2 ); k++ ) {
2925 connectivity->at(k + 1) = fine_quad->fine_id [ k ];
2926 }
2927 }
2928
2929 /* note: I am not using tmp_quad (and its ngb_elem_id) because tmp_quad exists
2930 * only if COMPLETE_DATA_STRUCTURE is required but that is supported for 2-manifolds only */
2931
2932 elem1 = global_quad_id(quad_id) + 1;
2933 for ( j = 0; j < 4; j++ ) {
2934 boundary->at(j + 1) = 1;
2935
2936 j1 = quad_ed_nd [ j ] [ 0 ];
2937 j2 = quad_ed_nd [ j ] [ 1 ];
2938
2939 /* note: I cannot use fe_quad (it is already released)
2940 *
2941 * node1 = fe_quad -> node[j1] -> id;
2942 * node2 = fe_quad -> node[j2] -> id;
2943 */
2944
2945 node1 = ( & ( fine_quad_array [ pos + j1 ] ) )->fine_id [ 0 ];
2946 node2 = ( & ( fine_quad_array [ pos + j2 ] ) )->fine_id [ 0 ];
2947
2948 flag = 1;
2949 for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
2950 elem2 = node_con_elems [ k ];
2951 if ( elem2 == elem1 ) {
2952 continue;
2953 }
2954
2955 if ( is_edge(elem2) == 0 ) { /* macro */
2956 for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
2957 if ( elem2 == node_con_elems [ m ] ) {
2958 boundary->at(j + 1) = flag = 0;
2959 break;
2960 }
2961 }
2962
2963 if ( flag == 0 ) {
2964 break;
2965 }
2966 }
2967 }
2968 }
2969
2970 quad_id++;
2971 break;
2972
2973 case EGT_tetra_1:
2974 pos = tetra_id * 4;
2975 fine_hexa = & ( fine_hexa_array [ pos ] );
2976 for ( j = 0; j < 4; j++, fine_hexa++ ) {
2977 connectivity = refinedElement.giveFineNodeArray(j + 1);
2978 connectivity->resize( ( level + 2 ) * ( level + 2 ) * ( level + 2 ) );
2979 for ( k = 0; k < ( level + 2 ) * ( level + 2 ) * ( level + 2 ); k++ ) {
2980 connectivity->at(k + 1) = fine_hexa->fine_id [ k ];
2981 }
2982 }
2983
2984 tmp_tetra = & ( tmp_tetra_array [ tetra_id ] );
2985 for ( j = 0; j < 4; j++ ) {
2986 if ( tmp_tetra->ngb_elem_id [ j ] != 0 ) {
2987 boundary->at(j + 1) = 0;
2988 } else {
2989 boundary->at(j + 1) = 1;
2990 }
2991 }
2992
2993 tetra_id++;
2994 break;
2995
2996 case EGT_hexa_1:
2997 pos = fe_tetras * 4 + hexa_id * 8;
2998 fine_hexa = & ( fine_hexa_array [ pos ] );
2999 for ( j = 0; j < 8; j++, fine_hexa++ ) {
3000 connectivity = refinedElement.giveFineNodeArray(j + 1);
3001 connectivity->resize( ( level + 2 ) * ( level + 2 ) * ( level + 2 ) );
3002 for ( k = 0; k < ( level + 2 ) * ( level + 2 ) * ( level + 2 ); k++ ) {
3003 connectivity->at(k + 1) = fine_hexa->fine_id [ k ];
3004 }
3005 }
3006
3007 tmp_hexa = & ( tmp_hexa_array [ hexa_id ] );
3008 for ( j = 0; j < 6; j++ ) {
3009 if ( tmp_hexa->ngb_elem_id [ j ] != 0 ) {
3010 boundary->at(j + 1) = 0;
3011 } else {
3012 boundary->at(j + 1) = 1;
3013 }
3014 }
3015
3016 hexa_id++;
3017 break;
3018
3019 default:
3020 error_message("Not supported element type");
3021 break;
3022 }
3023 }
3024
3025 /* the following 4 arrays cannot be released earlier because they are used to recover boundary flag */
3026
3027 if ( node_num_elems != NULL ) {
3028 free(node_num_elems);
3029 }
3030
3031 if ( node_con_elems != NULL ) {
3032 free(node_con_elems);
3033 }
3034
3035 if ( tmp_tetra_array != NULL ) {
3036 free(tmp_tetra_array);
3037 }
3038
3039 if ( tmp_hexa_array != NULL ) {
3040 free(tmp_hexa_array);
3041 }
3042
3043 fine_edge = fine_edge_array;
3044 for ( i = 0; i < fine_edges; i++, fine_edge++ ) {
3045 if ( fine_edge->fine_id != NULL ) {
3046 free(fine_edge->fine_id);
3047 }
3048 }
3049
3050 fine_quad = fine_quad_array;
3051 for ( i = 0; i < fine_quads; i++, fine_quad++ ) {
3052 if ( fine_quad->fine_id != NULL ) {
3053 free(fine_quad->fine_id);
3054 }
3055 }
3056
3057 fine_hexa = fine_hexa_array;
3058 for ( i = 0; i < fine_hexas; i++, fine_hexa++ ) {
3059 if ( fine_hexa->fine_id != NULL ) {
3060 free(fine_hexa->fine_id);
3061 }
3062 }
3063
3064 if ( fine_edge_array != NULL ) {
3065 free(fine_edge_array);
3066 }
3067
3068 if ( fine_quad_array != NULL ) {
3069 free(fine_quad_array);
3070 }
3071
3072 if ( fine_hexa_array != NULL ) {
3073 free(fine_hexa_array);
3074 }
3075
3076 this->nodes = fe_nodes + fine_nodes + refine_nodes;
3077 this->edges = fine_edges * ( level + 1 );
3078 this->quads = fine_quads * ( level + 1 ) * ( level + 1 );
3079 this->hexas = fine_hexas * ( level + 1 ) * ( level + 1 ) * ( level + 1 );
3080 this->elems = this->edges + this->quads + this->hexas;
3081
3082 return ( 0 );
3083}
3084} // end namespace oofem
oofem::Domain::giveNumberOfElements
int giveNumberOfElements() const
Returns number of elements in domain.
Definition domain.h:463
oofem::Domain::giveNumberOfDofManagers
int giveNumberOfDofManagers() const
Returns number of dof managers in domain.
Definition domain.h:461
oofem::Domain::giveElement
Element * giveElement(int n)
Definition domain.C:165
oofem::Element::giveNode
Node * giveNode(int i) const
Definition element.h:629
oofem::Element::giveGeometryType
virtual Element_Geometry_Type giveGeometryType() const =0
oofem::FEMComponent::giveNumber
int giveNumber() const
Definition femcmpnn.h:104
oofem::IntArray::resize
void resize(int n)
Definition intarray.C:73
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
oofem::RefinedElement::giveBoundaryFlagArray
IntArray * giveBoundaryFlagArray(void)
Definition refinedelement.h:61
oofem::RefinedElement::giveFineNodeArray
IntArray * giveFineNodeArray(int node)
Definition refinedelement.C:81
OOFEM_LOG_DEBUG
#define OOFEM_LOG_DEBUG(...)
Definition logger.h:144
TETRA_ELEM
#define TETRA_ELEM
Definition refinedmesh.C:45
QUAD_ELEM
#define QUAD_ELEM
Definition refinedmesh.C:44
matrix_3d
#define matrix_3d(ARRAY, U, V, W)
Definition refinedmesh.C:93
is_hexa
#define is_hexa(ID)
Definition refinedmesh.C:67
global_hexa_id
#define global_hexa_id(ID)
Definition refinedmesh.C:82
global_face_id
#define global_face_id(ID)
Definition refinedmesh.C:79
error_message
#define error_message(MSG)
Definition refinedmesh.C:96
is_tetra
#define is_tetra(ID)
Definition refinedmesh.C:66
is_face
#define is_face(ID)
Definition refinedmesh.C:64
is_quad
#define is_quad(ID)
Definition refinedmesh.C:65
global_tetra_id
#define global_tetra_id(ID)
Definition refinedmesh.C:81
FACE_ELEM
#define FACE_ELEM
Definition refinedmesh.C:43
local_quad_id
#define local_quad_id(ID)
Definition refinedmesh.C:87
elem_type
#define elem_type(ID)
Definition refinedmesh.C:57
local_tetra_id
#define local_tetra_id(ID)
Definition refinedmesh.C:88
HEXA_ELEM
#define HEXA_ELEM
Definition refinedmesh.C:46
local_hexa_id
#define local_hexa_id(ID)
Definition refinedmesh.C:89
local_edge_id
#define local_edge_id(ID)
Definition refinedmesh.C:85
EDGE_ELEM
#define EDGE_ELEM
Definition refinedmesh.C:42
local_face_id
#define local_face_id(ID)
Definition refinedmesh.C:86
global_quad_id
#define global_quad_id(ID)
Definition refinedmesh.C:80
is_edge
#define is_edge(ID)
Definition refinedmesh.C:63
oofem::RefinedMesh::fe_edge_rec
Definition refinedmesh.h:77
oofem::RefinedMesh::fe_edge_rec::node
fe_node_rec * node[2]
Definition refinedmesh.h:78
oofem::RefinedMesh::fe_face_rec
Definition refinedmesh.h:81
oofem::RefinedMesh::fe_face_rec::node
fe_node_rec * node[3]
Definition refinedmesh.h:82
oofem::RefinedMesh::fe_hexa_rec
Definition refinedmesh.h:93
oofem::RefinedMesh::fe_hexa_rec::node
fe_node_rec * node[8]
Definition refinedmesh.h:94
oofem::RefinedMesh::fe_node_rec
Definition refinedmesh.h:73
oofem::RefinedMesh::fe_node_rec::id
int id
Definition refinedmesh.h:74
oofem::RefinedMesh::fe_quad_rec
Definition refinedmesh.h:85
oofem::RefinedMesh::fe_quad_rec::node
fe_node_rec * node[4]
Definition refinedmesh.h:86
oofem::RefinedMesh::fe_tetra_rec
Definition refinedmesh.h:89
oofem::RefinedMesh::fe_tetra_rec::node
fe_node_rec * node[4]
Definition refinedmesh.h:90
oofem::RefinedMesh::fine_edge_rec
Definition refinedmesh.h:60
oofem::RefinedMesh::fine_edge_rec::fine_id
int * fine_id
Definition refinedmesh.h:61
oofem::RefinedMesh::fine_hexa_rec
Definition refinedmesh.h:68
oofem::RefinedMesh::fine_hexa_rec::fine_id
int * fine_id
Definition refinedmesh.h:69
oofem::RefinedMesh::fine_quad_rec
Definition refinedmesh.h:64
oofem::RefinedMesh::fine_quad_rec::fine_id
int * fine_id
Definition refinedmesh.h:65
oofem::RefinedMesh::mesh_edge_rec
Definition refinedmesh.h:99
oofem::RefinedMesh::mesh_edge_rec::mid_node
fe_node_rec * mid_node
Definition refinedmesh.h:101
oofem::RefinedMesh::mesh_edge_rec::node
fe_node_rec * node[2]
Definition refinedmesh.h:100
oofem::RefinedMesh::mesh_edge_rec::fine_id
int * fine_id[2]
Definition refinedmesh.h:102
oofem::RefinedMesh::mesh_face_rec
Definition refinedmesh.h:105
oofem::RefinedMesh::mesh_face_rec::mid_node
fe_node_rec * mid_node
Definition refinedmesh.h:107
oofem::RefinedMesh::mesh_face_rec::node
fe_node_rec * node[3]
Definition refinedmesh.h:106
oofem::RefinedMesh::mesh_face_rec::fine_id
int * fine_id[3]
Definition refinedmesh.h:108
oofem::RefinedMesh::mesh_quad_rec
Definition refinedmesh.h:111
oofem::RefinedMesh::mesh_quad_rec::fine_id
int * fine_id[4]
Definition refinedmesh.h:114
oofem::RefinedMesh::mesh_quad_rec::node
fe_node_rec * node[4]
Definition refinedmesh.h:112
oofem::RefinedMesh::mesh_quad_rec::mid_node
fe_node_rec * mid_node
Definition refinedmesh.h:113
oofem::RefinedMesh::tmp_hexa_rec
Definition refinedmesh.h:137
oofem::RefinedMesh::tmp_hexa_rec::quad
mesh_quad_rec * quad[6]
Definition refinedmesh.h:139
oofem::RefinedMesh::tmp_hexa_rec::ngb_elem_id
int ngb_elem_id[6]
Definition refinedmesh.h:138
oofem::RefinedMesh::tmp_tetra_rec
Definition refinedmesh.h:132
oofem::RefinedMesh::tmp_tetra_rec::ngb_elem_id
int ngb_elem_id[4]
Definition refinedmesh.h:133
oofem::RefinedMesh::tmp_tetra_rec::face
mesh_face_rec * face[4]
Definition refinedmesh.h:134
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011