OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
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 - 2013 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 
41 namespace 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 
118 int
119 RefinedMesh :: 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  for ( k = 0; k < 3; k++ ) {
1978  if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
1979  if ( swap [ 0 ] == 1 ) {
1980  kk = k;
1981  } else {
1982  if ( ( kk = k + 1 ) == 3 ) {
1983  kk = 0;
1984  }
1985  }
1986 
1987  break;
1988  }
1989  }
1990 
1991  for ( m = 0; m < level + 2; m++ ) {
1992  tmp_array [ j ] [ pos++ ] = mesh_face->fine_id [ kk ] [ p++ ];
1993  }
1994 
1995  if ( level != 0 ) {
1996  p = ( level + 2 ) * ( level + 1 );
1997  mesh_face = mesh_fc [ j + 1 ];
1998  for ( k = 0; k < 3; k++ ) {
1999  if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2000  if ( swap [ j + 1 ] == 1 ) {
2001  kk = k;
2002  } else {
2003  if ( ( kk = k + 1 ) == 3 ) {
2004  kk = 0;
2005  }
2006  }
2007 
2008  break;
2009  }
2010  }
2011 
2012  for ( n = 1; n < level + 1; n++ ) {
2013  tmp_array [ j ] [ pos++ ] = mesh_face->fine_id [ kk ] [ ++p ];
2014  for ( m = 0; m < level; m++ ) {
2015  tmp_array [ j ] [ pos++ ] = ++refine_node_id;
2016  }
2017 
2018  tmp_array [ j ] [ pos++ ] = tmp_array_cen [ 0 ] [ n ];
2019  }
2020  }
2021 
2022  for ( m = 0; m < level + 2; m++ ) {
2023  tmp_array [ j ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ m ];
2024  }
2025  }
2026 
2027  /* form inner tmp fine quads "parallel" to bottom face;
2028  * each one is oriented as fine quads on bottom face if not swapped (this is 0-1-2) */
2029 
2030  for ( j = 0; j < 3; j++ ) {
2031  j1 = face_ed_nd [ j ] [ 0 ];
2032  j2 = face_ed_nd [ j ] [ 1 ];
2033 
2034  pos = 0;
2035 
2036  p = ( level + 2 ) * ( level + 1 );
2037  mesh_face = mesh_fc [ j + 1 ];
2038  for ( k = 0; k < 3; k++ ) {
2039  if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2040  if ( swap [ j + 1 ] == 1 ) {
2041  kk = k;
2042  } else {
2043  if ( ( kk = k + 1 ) == 3 ) {
2044  kk = 0;
2045  }
2046  }
2047 
2048  break;
2049  }
2050  }
2051 
2052  for ( m = 0; m < level + 2; m++ ) {
2053  tmp_array [ j + 3 ] [ pos++ ] = mesh_face->fine_id [ kk ] [ p++ ];
2054  }
2055 
2056  if ( ( jj = j ) == 0 ) {
2057  jj = 3;
2058  }
2059 
2060  if ( level != 0 ) {
2061  p = ( level + 2 ) * ( level + 1 );
2062  mesh_face = mesh_fc [ jj ];
2063  for ( k = 0; k < 3; k++ ) {
2064  if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2065  if ( swap [ jj ] == 1 ) {
2066  kk = k;
2067  } else {
2068  if ( ( kk = k + 1 ) == 3 ) {
2069  kk = 0;
2070  }
2071  }
2072 
2073  break;
2074  }
2075  }
2076 
2077  for ( n = 1; n < level + 1; n++ ) {
2078  tmp_array [ j + 3 ] [ pos++ ] = mesh_face->fine_id [ kk ] [ ++p ];
2079  for ( m = 0; m < level; m++ ) {
2080  tmp_array [ j + 3 ] [ pos++ ] = ++refine_node_id;
2081  }
2082 
2083  tmp_array [ j + 3 ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ n ];
2084  }
2085  }
2086 
2087  for ( m = 0; m < level + 2; m++ ) {
2088  tmp_array [ j + 3 ] [ pos++ ] = tmp_array_cen [ jj ] [ m ];
2089  }
2090  }
2091 
2092  /* form bottom fine hexas */
2093 
2094  /* filling 3d matrix in a sequential way would be too difficult
2095  * with respect to switching among boundary and inner sources;
2096  * therefore it is much easier to traverse boundary sources
2097  * and localize them into 3d matrix;
2098  * this is done using macros matrix_3d and possible matrix_2d;
2099  * when traversing whole boudary source it is better to use direct
2100  * access instead of macro matrix_2d */
2101 
2102  /* note: I ignore that some entries of 3d matrix are rewritten several times;
2103  * in order to enable direct access */
2104 
2105  for ( j = 0; j < 3; j++ ) {
2106  fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2107  if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2108  error_message("Memory allocation error");
2109  }
2110 
2111  mesh_face = mesh_fc [ 0 ];
2112  for ( k = 0; k < 3; k++ ) {
2113  if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2114  break;
2115  }
2116  }
2117 
2118  if ( swap [ 0 ] == 1 ) {
2119  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2120  for ( m = 0; m < level + 2; m++, pos++ ) {
2121  matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2122  }
2123  }
2124  } else {
2125  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2126  for ( m = 0; m < level + 2; m++, pos++ ) {
2127  matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2128  }
2129  }
2130  }
2131 
2132  mesh_face = mesh_fc [ j + 1 ];
2133  for ( k = 0; k < 3; k++ ) {
2134  if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2135  break;
2136  }
2137  }
2138 
2139  if ( swap [ j + 1 ] == 1 ) {
2140  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2141  for ( m = 0; m < level + 2; m++, pos++ ) {
2142  matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2143  }
2144  }
2145  } else {
2146  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2147  for ( m = 0; m < level + 2; m++, pos++ ) {
2148  matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2149  }
2150  }
2151  }
2152 
2153  if ( ( jj = j ) == 0 ) {
2154  jj = 3;
2155  }
2156 
2157  mesh_face = mesh_fc [ jj ];
2158  for ( k = 0; k < 3; k++ ) {
2159  if ( mesh_face->node [ k ] == fe_tetra->node [ j ] ) {
2160  break;
2161  }
2162  }
2163 
2164  if ( swap [ jj ] == 1 ) {
2165  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2166  for ( m = 0; m < level + 2; m++, pos++ ) {
2167  matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2168  }
2169  }
2170  } else {
2171  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2172  for ( m = 0; m < level + 2; m++, pos++ ) {
2173  matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2174  }
2175  }
2176  }
2177 
2178  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2179  for ( m = 0; m < level + 2; m++, pos++ ) {
2180  matrix_3d(fine_hexa->fine_id, level + 1, m, n) = tmp_array [ j ] [ pos ]; /* macro */
2181  matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ jj - 1 ] [ pos ]; /* macro */
2182  matrix_3d(fine_hexa->fine_id, m, n, level + 1) = tmp_array [ j + 3 ] [ pos ]; /* macro */
2183  }
2184  }
2185 
2186  if ( level != 0 ) {
2187  for ( k = 1; k < level + 1; k++ ) {
2188  for ( n = 1; n < level + 1; n++ ) {
2189  for ( m = 1; m < level + 1; m++ ) {
2190  matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2191  }
2192  }
2193  }
2194  }
2195  }
2196 
2197  /* form top hexa */
2198 
2199  fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2200  if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2201  error_message("Memory allocation error");
2202  }
2203 
2204  mesh_face = mesh_fc [ 1 ];
2205  for ( k = 0; k < 3; k++ ) {
2206  if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2207  break;
2208  }
2209  }
2210 
2211  if ( swap [ 1 ] == 1 ) {
2212  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2213  for ( m = 0; m < level + 2; m++, pos++ ) {
2214  matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2215  }
2216  }
2217  } else {
2218  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2219  for ( m = 0; m < level + 2; m++, pos++ ) {
2220  matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2221  }
2222  }
2223  }
2224 
2225  mesh_face = mesh_fc [ 2 ];
2226  for ( k = 0; k < 3; k++ ) {
2227  if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2228  break;
2229  }
2230  }
2231 
2232  if ( swap [ 2 ] == 1 ) {
2233  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2234  for ( m = 0; m < level + 2; m++, pos++ ) {
2235  matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2236  }
2237  }
2238  } else {
2239  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2240  for ( m = 0; m < level + 2; m++, pos++ ) {
2241  matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2242  }
2243  }
2244  }
2245 
2246  mesh_face = mesh_fc [ 3 ];
2247  for ( k = 0; k < 3; k++ ) {
2248  if ( mesh_face->node [ k ] == fe_tetra->node [ 3 ] ) {
2249  break;
2250  }
2251  }
2252 
2253  if ( swap [ 3 ] == 1 ) {
2254  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2255  for ( m = 0; m < level + 2; m++, pos++ ) {
2256  matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2257  }
2258  }
2259  } else {
2260  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2261  for ( m = 0; m < level + 2; m++, pos++ ) {
2262  matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_face->fine_id [ k ] [ pos ]; /* macro */
2263  }
2264  }
2265  }
2266 
2267  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2268  for ( m = 0; m < level + 2; m++, pos++ ) {
2269  matrix_3d(fine_hexa->fine_id, level + 1, n, m) = tmp_array [ 0 + 3 ] [ pos ]; /* macro */
2270  matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ 2 + 3 ] [ pos ]; /* macro */
2271  matrix_3d(fine_hexa->fine_id, n, m, level + 1) = tmp_array [ 1 + 3 ] [ pos ]; /* macro */
2272  }
2273  }
2274 
2275  if ( level != 0 ) {
2276  for ( k = 1; k < level + 1; k++ ) {
2277  for ( n = 1; n < level + 1; n++ ) {
2278  for ( m = 1; m < level + 1; m++ ) {
2279  matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2280  }
2281  }
2282  }
2283  }
2284  }
2285 
2286  fe_hexa = fe_hexa_array;
2287  tmp_hexa = tmp_hexa_array;
2288  for ( i = 0; i < fe_hexas; i++, fe_hexa++, tmp_hexa++, fine_node++ ) {
2289  fine_node->id = ++fine_node_id;
2290 
2291  for ( j = 0; j < 6; j++ ) {
2292  tmp_array_cen [ j ] [ 0 ] = tmp_hexa->quad [ j ]->mid_node->id;
2293  for ( k = 1; k < level + 1; k++ ) {
2294  tmp_array_cen [ j ] [ k ] = ++refine_node_id;
2295  }
2296 
2297  tmp_array_cen [ j ] [ level + 1 ] = fine_node->id;
2298  }
2299 
2300  /* mark swapped faces (see ordering of nodes on hexa faces above);
2301  * bottom and top faces are swapped if normal is outer;
2302  * side faces are swapped if normal is inner */
2303 
2304  for ( j = 0; j < 6; j++ ) {
2305  mesh_quad = mesh_qd [ j ] = tmp_hexa->quad [ j ];
2306  for ( k = 0; k < 4; k++ ) {
2307  if ( mesh_quad->node [ k ] == fe_hexa->node [ hexa_fc_nd [ j ] [ 0 ] ] ) {
2308  if ( ( kk = k + 1 ) == 4 ) {
2309  kk = 0;
2310  }
2311 
2312  if ( mesh_quad->node [ kk ] == fe_hexa->node [ hexa_fc_nd [ j ] [ 1 ] ] ) {
2313  swap [ j ] = 0;
2314  } else {
2315  swap [ j ] = 1;
2316  }
2317 
2318  break;
2319  }
2320  }
2321  }
2322 
2323  /* form inner bottom tmp fine quads "perpendicular" to bottom face;
2324  * each one is oriented u: from center of bottom edge to center of bottom face
2325  * v: from center of bottom edge to center of side face */
2326 
2327  for ( j = 0; j < 4; j++ ) {
2328  pos = 0;
2329 
2330  p = ( level + 2 ) * ( level + 1 );
2331  mesh_quad = mesh_qd [ 0 ];
2332  for ( k = 0; k < 4; k++ ) {
2333  if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2334  if ( swap [ 0 ] == 1 ) {
2335  kk = k;
2336  } else {
2337  if ( ( kk = k + 1 ) == 4 ) {
2338  kk = 0;
2339  }
2340  }
2341 
2342  break;
2343  }
2344  }
2345 
2346  for ( m = 0; m < level + 2; m++ ) {
2347  tmp_array [ j ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ p++ ];
2348  }
2349 
2350  if ( level != 0 ) {
2351  p = ( level + 2 ) * ( level + 1 );
2352  mesh_quad = mesh_qd [ j + 1 ];
2353  for ( k = 0; k < 4; k++ ) {
2354  if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2355  if ( swap [ j + 1 ] == 1 ) {
2356  kk = k;
2357  } else {
2358  if ( ( kk = k + 1 ) == 4 ) {
2359  kk = 0;
2360  }
2361  }
2362 
2363  break;
2364  }
2365  }
2366 
2367  for ( n = 1; n < level + 1; n++ ) {
2368  tmp_array [ j ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ ++p ];
2369  for ( m = 0; m < level; m++ ) {
2370  tmp_array [ j ] [ pos++ ] = ++refine_node_id;
2371  }
2372 
2373  tmp_array [ j ] [ pos++ ] = tmp_array_cen [ 0 ] [ n ];
2374  }
2375  }
2376 
2377  for ( m = 0; m < level + 2; m++ ) {
2378  tmp_array [ j ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ m ];
2379  }
2380  }
2381 
2382  /* form top tmp fine quads "perpendicular" to top face;
2383  * each one is oriented u: from center of top edge to center of top face
2384  * v: from center of top edge to center of side face */
2385 
2386  for ( j = 0; j < 4; j++ ) {
2387  pos = 0;
2388 
2389  p = ( level + 2 ) * ( level + 1 );
2390  mesh_quad = mesh_qd [ 5 ];
2391  for ( k = 0; k < 4; k++ ) {
2392  if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2393  if ( swap [ 5 ] == 0 ) {
2394  kk = k;
2395  } else {
2396  if ( ( kk = k + 1 ) == 4 ) {
2397  kk = 0;
2398  }
2399  }
2400 
2401  break;
2402  }
2403  }
2404 
2405  for ( m = 0; m < level + 2; m++ ) {
2406  tmp_array [ j + 4 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ p++ ];
2407  }
2408 
2409  if ( level != 0 ) {
2410  p = ( level + 2 ) * ( level + 1 );
2411  mesh_quad = mesh_qd [ j + 1 ];
2412  for ( k = 0; k < 4; k++ ) {
2413  if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2414  if ( swap [ j + 1 ] == 0 ) {
2415  kk = k;
2416  } else {
2417  if ( ( kk = k + 1 ) == 4 ) {
2418  kk = 0;
2419  }
2420  }
2421 
2422  break;
2423  }
2424  }
2425 
2426  for ( n = 1; n < level + 1; n++ ) {
2427  tmp_array [ j + 4 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ ++p ];
2428  for ( m = 0; m < level; m++ ) {
2429  tmp_array [ j + 4 ] [ pos++ ] = ++refine_node_id;
2430  }
2431 
2432  tmp_array [ j + 4 ] [ pos++ ] = tmp_array_cen [ 5 ] [ n ];
2433  }
2434  }
2435 
2436  for ( m = 0; m < level + 2; m++ ) {
2437  tmp_array [ j + 4 ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ m ];
2438  }
2439  }
2440 
2441  /* form inner tmp fine quads "parallel" to bottom and top face;
2442  * each one is oriented as fine quads on bottom face if not swapped (this is 0-1-2-3) */
2443 
2444  for ( j = 0; j < 4; j++ ) {
2445  pos = 0;
2446 
2447  p = ( level + 2 ) * ( level + 1 );
2448  mesh_quad = mesh_qd [ j + 1 ];
2449  for ( k = 0; k < 4; k++ ) {
2450  if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2451  if ( swap [ j + 1 ] == 1 ) {
2452  kk = k;
2453  } else {
2454  if ( ( kk = k + 1 ) == 4 ) {
2455  kk = 0;
2456  }
2457  }
2458 
2459  break;
2460  }
2461  }
2462 
2463  for ( m = 0; m < level + 2; m++ ) {
2464  tmp_array [ j + 8 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ p++ ];
2465  }
2466 
2467  if ( ( jj = j ) == 0 ) {
2468  jj = 4;
2469  }
2470 
2471  if ( level != 0 ) {
2472  p = ( level + 2 ) * ( level + 1 );
2473  mesh_quad = mesh_qd [ jj ];
2474  for ( k = 0; k < 4; k++ ) {
2475  if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2476  if ( swap [ jj ] == 1 ) {
2477  kk = k;
2478  } else {
2479  if ( ( kk = k + 1 ) == 4 ) {
2480  kk = 0;
2481  }
2482  }
2483 
2484  break;
2485  }
2486  }
2487 
2488  for ( n = 1; n < level + 1; n++ ) {
2489  tmp_array [ j + 8 ] [ pos++ ] = mesh_quad->fine_id [ kk ] [ ++p ];
2490  for ( m = 0; m < level; m++ ) {
2491  tmp_array [ j + 8 ] [ pos++ ] = ++refine_node_id;
2492  }
2493 
2494  tmp_array [ j + 8 ] [ pos++ ] = tmp_array_cen [ j + 1 ] [ n ];
2495  }
2496  }
2497 
2498  for ( m = 0; m < level + 2; m++ ) {
2499  tmp_array [ j + 8 ] [ pos++ ] = tmp_array_cen [ jj ] [ m ];
2500  }
2501  }
2502 
2503  /* form bottom fine hexas */
2504 
2505  /* filling 3d matrix in a sequential way would be too difficult
2506  * with respect to switching among boundary and inner sources;
2507  * therefore it is much easier to traverse boundary sources
2508  * and localize them into 3d matrix;
2509  * this is done using macros matrix_3d and possible matrix_2d;
2510  * when traversing whole boudary source it is better to use direct
2511  * access instead of macro matrix_2d */
2512 
2513  /* note: I ignore that some entries of 3d matrix are rewritten several times;
2514  * in order to enable direct access */
2515 
2516  for ( j = 0; j < 4; j++ ) {
2517  fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2518  if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2519  error_message("Memory allocation error");
2520  }
2521 
2522  mesh_quad = mesh_qd [ 0 ];
2523  for ( k = 0; k < 4; k++ ) {
2524  if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2525  break;
2526  }
2527  }
2528 
2529  if ( swap [ 0 ] == 1 ) {
2530  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2531  for ( m = 0; m < level + 2; m++, pos++ ) {
2532  matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2533  }
2534  }
2535  } else {
2536  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2537  for ( m = 0; m < level + 2; m++, pos++ ) {
2538  matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2539  }
2540  }
2541  }
2542 
2543  mesh_quad = mesh_qd [ j + 1 ];
2544  for ( k = 0; k < 4; k++ ) {
2545  if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2546  break;
2547  }
2548  }
2549 
2550  if ( swap [ j + 1 ] == 1 ) {
2551  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2552  for ( m = 0; m < level + 2; m++, pos++ ) {
2553  matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2554  }
2555  }
2556  } else {
2557  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2558  for ( m = 0; m < level + 2; m++, pos++ ) {
2559  matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2560  }
2561  }
2562  }
2563 
2564  if ( ( jj = j ) == 0 ) {
2565  jj = 4;
2566  }
2567 
2568  mesh_quad = mesh_qd [ jj ];
2569  for ( k = 0; k < 4; k++ ) {
2570  if ( mesh_quad->node [ k ] == fe_hexa->node [ j ] ) {
2571  break;
2572  }
2573  }
2574 
2575  if ( swap [ jj ] == 1 ) {
2576  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2577  for ( m = 0; m < level + 2; m++, pos++ ) {
2578  matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2579  }
2580  }
2581  } else {
2582  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2583  for ( m = 0; m < level + 2; m++, pos++ ) {
2584  matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2585  }
2586  }
2587  }
2588 
2589  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2590  for ( m = 0; m < level + 2; m++, pos++ ) {
2591  matrix_3d(fine_hexa->fine_id, level + 1, m, n) = tmp_array [ j ] [ pos ]; /* macro */
2592  matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ jj - 1 ] [ pos ]; /* macro */
2593  matrix_3d(fine_hexa->fine_id, m, n, level + 1) = tmp_array [ j + 8 ] [ pos ]; /* macro */
2594  }
2595  }
2596 
2597  if ( level != 0 ) {
2598  for ( k = 1; k < level + 1; k++ ) {
2599  for ( n = 1; n < level + 1; n++ ) {
2600  for ( m = 1; m < level + 1; m++ ) {
2601  matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2602  }
2603  }
2604  }
2605  }
2606  }
2607 
2608  /* form top fine hexas */
2609 
2610  for ( j = 0; j < 4; j++ ) {
2611  fine_hexa = & ( fine_hexa_array [ fine_hexa_id++ ] );
2612  if ( ( fine_hexa->fine_id = ( int * ) calloc( ( level + 2 ) * ( level + 2 ) * ( level + 2 ), sizeof( int ) ) ) == NULL ) {
2613  error_message("Memory allocation error");
2614  }
2615 
2616  mesh_quad = mesh_qd [ 5 ];
2617  for ( k = 0; k < 4; k++ ) {
2618  if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2619  break;
2620  }
2621  }
2622 
2623  if ( swap [ 5 ] == 1 ) {
2624  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2625  for ( m = 0; m < level + 2; m++, pos++ ) {
2626  matrix_3d(fine_hexa->fine_id, n, m, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2627  }
2628  }
2629  } else {
2630  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2631  for ( m = 0; m < level + 2; m++, pos++ ) {
2632  matrix_3d(fine_hexa->fine_id, m, n, 0) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2633  }
2634  }
2635  }
2636 
2637  mesh_quad = mesh_qd [ j + 1 ];
2638  for ( k = 0; k < 4; k++ ) {
2639  if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2640  break;
2641  }
2642  }
2643 
2644  if ( swap [ j + 1 ] == 1 ) {
2645  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2646  for ( m = 0; m < level + 2; m++, pos++ ) {
2647  matrix_3d(fine_hexa->fine_id, 0, m, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2648  }
2649  }
2650  } else {
2651  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2652  for ( m = 0; m < level + 2; m++, pos++ ) {
2653  matrix_3d(fine_hexa->fine_id, 0, n, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2654  }
2655  }
2656  }
2657 
2658  if ( ( jj = j ) == 0 ) {
2659  jj = 4;
2660  }
2661 
2662  mesh_quad = mesh_qd [ jj ];
2663  for ( k = 0; k < 4; k++ ) {
2664  if ( mesh_quad->node [ k ] == fe_hexa->node [ j + 4 ] ) {
2665  break;
2666  }
2667  }
2668 
2669  if ( swap [ jj ] == 1 ) {
2670  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2671  for ( m = 0; m < level + 2; m++, pos++ ) {
2672  matrix_3d(fine_hexa->fine_id, n, 0, m) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2673  }
2674  }
2675  } else {
2676  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2677  for ( m = 0; m < level + 2; m++, pos++ ) {
2678  matrix_3d(fine_hexa->fine_id, m, 0, n) = mesh_quad->fine_id [ k ] [ pos ]; /* macro */
2679  }
2680  }
2681  }
2682 
2683  for ( pos = 0, n = 0; n < level + 2; n++ ) {
2684  for ( m = 0; m < level + 2; m++, pos++ ) {
2685  matrix_3d(fine_hexa->fine_id, level + 1, m, n) = tmp_array [ jj - 1 + 4 ] [ pos ]; /* macro */
2686  matrix_3d(fine_hexa->fine_id, m, level + 1, n) = tmp_array [ j + 4 ] [ pos ]; /* macro */
2687  matrix_3d(fine_hexa->fine_id, n, m, level + 1) = tmp_array [ j + 8 ] [ pos ]; /* macro */
2688  }
2689  }
2690 
2691  if ( level != 0 ) {
2692  for ( k = 1; k < level + 1; k++ ) {
2693  for ( n = 1; n < level + 1; n++ ) {
2694  for ( m = 1; m < level + 1; m++ ) {
2695  matrix_3d(fine_hexa->fine_id, m, n, k) = ++refine_node_id; /* macro */
2696  }
2697  }
2698  }
2699  }
2700  }
2701  }
2702 
2703  for ( j = 0; j < 6; j++ ) {
2704  free(tmp_array_cen [ j ]);
2705  free(tmp_array [ j ]);
2706  free(tmp_array [ j + 6 ]);
2707  }
2708 
2709  if ( fe_node_array != NULL ) {
2710  free(fe_node_array);
2711  }
2712 
2713  if ( fe_edge_array != NULL ) {
2714  free(fe_edge_array);
2715  }
2716 
2717  if ( fe_face_array != NULL ) {
2718  free(fe_face_array);
2719  }
2720 
2721  if ( fe_quad_array != NULL ) {
2722  free(fe_quad_array);
2723  }
2724 
2725  if ( fe_tetra_array != NULL ) {
2726  free(fe_tetra_array);
2727  }
2728 
2729  if ( fe_hexa_array != NULL ) {
2730  free(fe_hexa_array);
2731  }
2732 
2733  mesh_face = mesh_face_array;
2734  for ( i = 0; i < mesh_faces; i++, mesh_face++ ) {
2735  for ( j = 0; j < 3; j++ ) {
2736  if ( mesh_face->fine_id [ j ] != NULL ) {
2737  free(mesh_face->fine_id [ j ]);
2738  }
2739  }
2740  }
2741 
2742  mesh_quad = mesh_quad_array;
2743  for ( i = 0; i < mesh_quads; i++, mesh_quad++ ) {
2744  for ( j = 0; j < 4; j++ ) {
2745  if ( mesh_quad->fine_id [ j ] != NULL ) {
2746  free(mesh_quad->fine_id [ j ]);
2747  }
2748  }
2749  }
2750 
2751  if ( mesh_face_array != NULL ) {
2752  free(mesh_face_array);
2753  }
2754 
2755  if ( mesh_quad_array != NULL ) {
2756  free(mesh_quad_array);
2757  }
2758 
2759  /*
2760  * if(refine_nodes != 0){
2761  * if((refine_node_array = (fe_node_rec *)calloc(refine_nodes, sizeof(fe_node_rec))) == NULL)
2762  * error_message("Memory allocation error");
2763  *
2764  * fine_edge = fine_edge_array;
2765  * for(i = 0; i < fine_edges; i++, fine_edge++){
2766  * pos = 0;
2767  * for(m = 0; m < level + 2; m++, pos++){
2768  * node_id = fine_edge -> fine_id[pos];
2769  * if(node_id <= fe_nodes + fine_nodes)continue;
2770  *
2771  * refine_node = &(refine_node_array[node_id - fe_nodes - fine_nodes - 1]);
2772  * refine_node -> id = node_id;
2773  * }
2774  * }
2775  *
2776  * fine_quad = fine_quad_array;
2777  * for(i = 0; i < fine_quads; i++, fine_quad++){
2778  * pos = 0;
2779  * for(n = 0; n < level + 2; n++){
2780  * for(m = 0; m < level + 2; m++, pos++){
2781  * node_id = fine_quad -> fine_id[pos];
2782  * if(node_id <= fe_nodes + fine_nodes)continue;
2783  *
2784  * refine_node = &(refine_node_array[node_id - fe_nodes - fine_nodes - 1]);
2785  * refine_node -> id = node_id;
2786  * }
2787  * }
2788  * }
2789  *
2790  * fine_hexa = fine_hexa_array;
2791  * for(i = 0; i < fine_hexas; i++, fine_hexa++){
2792  * pos = 0;
2793  * for(k = 0; k < level + 2; k++){
2794  * for(n = 0; n < level + 2; n++){
2795  * for(m = 0; m < level + 2; m++, pos++){
2796  * node_id = fine_hexa -> fine_id[pos];
2797  * if(node_id <= fe_nodes + fine_nodes)continue;
2798  *
2799  * refine_node = &(refine_node_array[node_id - fe_nodes - fine_nodes - 1]);
2800  * refine_node -> id = node_id;
2801  * }
2802  * }
2803  * }
2804  * }
2805  * }
2806  */
2807 
2808  if ( fine_node_array != NULL ) {
2809  free(fine_node_array);
2810  }
2811 
2812  /*
2813  * if(refine_node_array != NULL)free(refine_node_array);
2814  */
2815 
2816  edge_id = face_id = quad_id = tetra_id = hexa_id = 0;
2817  for ( i = 0; i < fe_elems; i++ ) {
2818  element = d->giveElement(i + 1);
2819  RefinedElement &refinedElement = refinedElementList.at(i + 1);
2820  boundary = refinedElement.giveBoundaryFlagArray();
2821 
2822  switch ( element->giveGeometryType() ) {
2823  case EGT_line_1:
2824  pos = edge_id * 2;
2825  fine_edge = & ( fine_edge_array [ pos ] );
2826  for ( j = 0; j < 2; j++, fine_edge++ ) {
2827  connectivity = refinedElement.giveFineNodeArray(j + 1);
2828  connectivity->resize(level + 2);
2829  for ( k = 0; k < level + 2; k++ ) {
2830  connectivity->at(k + 1) = fine_edge->fine_id [ k ];
2831  }
2832  }
2833 
2834  for ( j = 0; j < 2; j++ ) {
2835  /* note: I cannot use fe_edge (it is already released)
2836  *
2837  * node = fe_edge -> node[j] -> id;
2838  */
2839 
2840  node = ( & ( fine_edge_array [ pos + j ] ) )->fine_id [ 0 ];
2841  if ( node_num_elems [ node ] - node_num_elems [ node - 1 ] == 1 ) {
2842  boundary->at(j + 1) = 1;
2843  } else {
2844  boundary->at(j + 1) = 0;
2845  }
2846  }
2847 
2848  edge_id++;
2849  break;
2850 
2851  case EGT_triangle_1:
2852  pos = face_id * 3;
2853  fine_quad = & ( fine_quad_array [ pos ] );
2854  for ( j = 0; j < 3; j++, fine_quad++ ) {
2855  connectivity = refinedElement.giveFineNodeArray(j + 1);
2856  connectivity->resize( ( level + 2 ) * ( level + 2 ) );
2857  for ( k = 0; k < ( level + 2 ) * ( level + 2 ); k++ ) {
2858  connectivity->at(k + 1) = fine_quad->fine_id [ k ];
2859  }
2860  }
2861 
2862  /* note: I am not using tmp_face (and its ngb_elem_id) because tmp_face exists
2863  * only if COMPLETE_DATA_STRUCTURE is required but that is supported for 2-manifolds only */
2864 
2865  elem1 = global_face_id(face_id) + 1;
2866  for ( j = 0; j < 3; j++ ) {
2867  boundary->at(j + 1) = 1;
2868 
2869  j1 = face_ed_nd [ j ] [ 0 ];
2870  j2 = face_ed_nd [ j ] [ 1 ];
2871 
2872  /* note: I cannot use fe_face (it is already released)
2873  *
2874  * node1 = fe_face -> node[j1] -> id;
2875  * node2 = fe_face -> node[j2] -> id;
2876  */
2877 
2878  node1 = ( & ( fine_quad_array [ pos + j1 ] ) )->fine_id [ 0 ];
2879  node2 = ( & ( fine_quad_array [ pos + j2 ] ) )->fine_id [ 0 ];
2880 
2881  flag = 1;
2882  for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
2883  elem2 = node_con_elems [ k ];
2884  if ( elem2 == elem1 ) {
2885  continue;
2886  }
2887 
2888  if ( is_edge(elem2) == 0 ) { /* macro */
2889  for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
2890  if ( elem2 == node_con_elems [ m ] ) {
2891  boundary->at(j + 1) = flag = 0;
2892  break;
2893  }
2894  }
2895 
2896  if ( flag == 0 ) {
2897  break;
2898  }
2899  }
2900  }
2901  }
2902 
2903  face_id++;
2904  break;
2905 
2906  case EGT_quad_1:
2907  pos = fe_faces * 3 + quad_id * 4;
2908  fine_quad = & ( fine_quad_array [ pos ] );
2909  for ( j = 0; j < 4; j++, fine_quad++ ) {
2910  connectivity = refinedElement.giveFineNodeArray(j + 1);
2911  connectivity->resize( ( level + 2 ) * ( level + 2 ) );
2912  for ( k = 0; k < ( level + 2 ) * ( level + 2 ); k++ ) {
2913  connectivity->at(k + 1) = fine_quad->fine_id [ k ];
2914  }
2915  }
2916 
2917  /* note: I am not using tmp_quad (and its ngb_elem_id) because tmp_quad exists
2918  * only if COMPLETE_DATA_STRUCTURE is required but that is supported for 2-manifolds only */
2919 
2920  elem1 = global_quad_id(quad_id) + 1;
2921  for ( j = 0; j < 4; j++ ) {
2922  boundary->at(j + 1) = 1;
2923 
2924  j1 = quad_ed_nd [ j ] [ 0 ];
2925  j2 = quad_ed_nd [ j ] [ 1 ];
2926 
2927  /* note: I cannot use fe_quad (it is already released)
2928  *
2929  * node1 = fe_quad -> node[j1] -> id;
2930  * node2 = fe_quad -> node[j2] -> id;
2931  */
2932 
2933  node1 = ( & ( fine_quad_array [ pos + j1 ] ) )->fine_id [ 0 ];
2934  node2 = ( & ( fine_quad_array [ pos + j2 ] ) )->fine_id [ 0 ];
2935 
2936  flag = 1;
2937  for ( k = node_num_elems [ node1 - 1 ]; k < node_num_elems [ node1 ]; k++ ) {
2938  elem2 = node_con_elems [ k ];
2939  if ( elem2 == elem1 ) {
2940  continue;
2941  }
2942 
2943  if ( is_edge(elem2) == 0 ) { /* macro */
2944  for ( m = node_num_elems [ node2 - 1 ]; m < node_num_elems [ node2 ]; m++ ) {
2945  if ( elem2 == node_con_elems [ m ] ) {
2946  boundary->at(j + 1) = flag = 0;
2947  break;
2948  }
2949  }
2950 
2951  if ( flag == 0 ) {
2952  break;
2953  }
2954  }
2955  }
2956  }
2957 
2958  quad_id++;
2959  break;
2960 
2961  case EGT_tetra_1:
2962  pos = tetra_id * 4;
2963  fine_hexa = & ( fine_hexa_array [ pos ] );
2964  for ( j = 0; j < 4; j++, fine_hexa++ ) {
2965  connectivity = refinedElement.giveFineNodeArray(j + 1);
2966  connectivity->resize( ( level + 2 ) * ( level + 2 ) * ( level + 2 ) );
2967  for ( k = 0; k < ( level + 2 ) * ( level + 2 ) * ( level + 2 ); k++ ) {
2968  connectivity->at(k + 1) = fine_hexa->fine_id [ k ];
2969  }
2970  }
2971 
2972  tmp_tetra = & ( tmp_tetra_array [ tetra_id ] );
2973  for ( j = 0; j < 4; j++ ) {
2974  if ( tmp_tetra->ngb_elem_id [ j ] != 0 ) {
2975  boundary->at(j + 1) = 0;
2976  } else {
2977  boundary->at(j + 1) = 1;
2978  }
2979  }
2980 
2981  tetra_id++;
2982  break;
2983 
2984  case EGT_hexa_1:
2985  pos = fe_tetras * 4 + hexa_id * 8;
2986  fine_hexa = & ( fine_hexa_array [ pos ] );
2987  for ( j = 0; j < 8; j++, fine_hexa++ ) {
2988  connectivity = refinedElement.giveFineNodeArray(j + 1);
2989  connectivity->resize( ( level + 2 ) * ( level + 2 ) * ( level + 2 ) );
2990  for ( k = 0; k < ( level + 2 ) * ( level + 2 ) * ( level + 2 ); k++ ) {
2991  connectivity->at(k + 1) = fine_hexa->fine_id [ k ];
2992  }
2993  }
2994 
2995  tmp_hexa = & ( tmp_hexa_array [ hexa_id ] );
2996  for ( j = 0; j < 6; j++ ) {
2997  if ( tmp_hexa->ngb_elem_id [ j ] != 0 ) {
2998  boundary->at(j + 1) = 0;
2999  } else {
3000  boundary->at(j + 1) = 1;
3001  }
3002  }
3003 
3004  hexa_id++;
3005  break;
3006 
3007  default:
3008  error_message("Not supported element type");
3009  break;
3010  }
3011  }
3012 
3013  /* the following 4 arrays cannot be released earlier because they are used to recover boundary flag */
3014 
3015  if ( node_num_elems != NULL ) {
3016  free(node_num_elems);
3017  }
3018 
3019  if ( node_con_elems != NULL ) {
3020  free(node_con_elems);
3021  }
3022 
3023  if ( tmp_tetra_array != NULL ) {
3024  free(tmp_tetra_array);
3025  }
3026 
3027  if ( tmp_hexa_array != NULL ) {
3028  free(tmp_hexa_array);
3029  }
3030 
3031  fine_edge = fine_edge_array;
3032  for ( i = 0; i < fine_edges; i++, fine_edge++ ) {
3033  if ( fine_edge->fine_id != NULL ) {
3034  free(fine_edge->fine_id);
3035  }
3036  }
3037 
3038  fine_quad = fine_quad_array;
3039  for ( i = 0; i < fine_quads; i++, fine_quad++ ) {
3040  if ( fine_quad->fine_id != NULL ) {
3041  free(fine_quad->fine_id);
3042  }
3043  }
3044 
3045  fine_hexa = fine_hexa_array;
3046  for ( i = 0; i < fine_hexas; i++, fine_hexa++ ) {
3047  if ( fine_hexa->fine_id != NULL ) {
3048  free(fine_hexa->fine_id);
3049  }
3050  }
3051 
3052  if ( fine_edge_array != NULL ) {
3053  free(fine_edge_array);
3054  }
3055 
3056  if ( fine_quad_array != NULL ) {
3057  free(fine_quad_array);
3058  }
3059 
3060  if ( fine_hexa_array != NULL ) {
3061  free(fine_hexa_array);
3062  }
3063 
3064  this->nodes = fe_nodes + fine_nodes + refine_nodes;
3065  this->edges = fine_edges * ( level + 1 );
3066  this->quads = fine_quads * ( level + 1 ) * ( level + 1 );
3067  this->hexas = fine_hexas * ( level + 1 ) * ( level + 1 ) * ( level + 1 );
3068  this->elems = this->edges + this->quads + this->hexas;
3069 
3070  return ( 0 );
3071 }
3072 } // end namespace oofem
oofem::RefinedMesh::fine_hexa_rec
Definition: refinedmesh.h:68
is_edge
#define is_edge(ID)
Definition: refinedmesh.C:63
local_quad_id
#define local_quad_id(ID)
Definition: refinedmesh.C:87
oofem::RefinedMesh::mesh_edge_rec::mid_node
fe_node_rec * mid_node
Definition: refinedmesh.h:101
oofem::Domain
Class and object Domain.
Definition: domain.h:115
oofem::RefinedMesh::fe_tetra_rec::node
fe_node_rec * node[4]
Definition: refinedmesh.h:90
HEXA_ELEM
#define HEXA_ELEM
Definition: refinedmesh.C:46
oofem::RefinedMesh::mesh_face_rec
Definition: refinedmesh.h:105
oofem::RefinedMesh::mesh_quad_rec::node
fe_node_rec * node[4]
Definition: refinedmesh.h:112
oofem::Domain::giveNumberOfDofManagers
int giveNumberOfDofManagers() const
Returns number of dof managers in domain.
Definition: domain.h:432
matrix_3d
#define matrix_3d(ARRAY, U, V, W)
Definition: refinedmesh.C:93
global_quad_id
#define global_quad_id(ID)
Definition: refinedmesh.C:80
oofem::RefinedMesh::fine_quad_rec::fine_id
int * fine_id
Definition: refinedmesh.h:65
global_face_id
#define global_face_id(ID)
Definition: refinedmesh.C:79
oofem::RefinedMesh::fe_edge_rec::node
fe_node_rec * node[2]
Definition: refinedmesh.h:78
TETRA_ELEM
#define TETRA_ELEM
Definition: refinedmesh.C:45
oofem::RefinedElement::giveFineNodeArray
IntArray * giveFineNodeArray(int node)
Definition: refinedelement.C:81
oofem::RefinedMesh::mesh_quad_rec::mid_node
fe_node_rec * mid_node
Definition: refinedmesh.h:113
oofem::Element
Abstract base class for all finite elements.
Definition: element.h:145
oofem::RefinedMesh::fe_quad_rec
Definition: refinedmesh.h:85
oofem::RefinedMesh::fe_hexa_rec
Definition: refinedmesh.h:93
oofem::Domain::giveNumberOfElements
int giveNumberOfElements() const
Returns number of elements in domain.
Definition: domain.h:434
oofem::RefinedMesh::fine_edge_rec
Definition: refinedmesh.h:60
OOFEM_LOG_DEBUG
#define OOFEM_LOG_DEBUG(...)
Definition: logger.h:128
oofem::IntArray
Class implementing an array of integers.
Definition: intarray.h:61
oofem::IntArray::at
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
is_tetra
#define is_tetra(ID)
Definition: refinedmesh.C:66
local_edge_id
#define local_edge_id(ID)
Definition: refinedmesh.C:85
oofem::RefinedMesh::mesh_quad_rec::fine_id
int * fine_id[4]
Definition: refinedmesh.h:114
oofem::RefinedMesh::fine_quad_rec
Definition: refinedmesh.h:64
QUAD_ELEM
#define QUAD_ELEM
Definition: refinedmesh.C:44
is_quad
#define is_quad(ID)
Definition: refinedmesh.C:65
local_face_id
#define local_face_id(ID)
Definition: refinedmesh.C:86
oofem::RefinedMesh::mesh_edge_rec::fine_id
int * fine_id[2]
Definition: refinedmesh.h:102
oofem::RefinedMesh::tmp_tetra_rec::face
mesh_face_rec * face[4]
Definition: refinedmesh.h:134
oofem::RefinedMesh::tmp_hexa_rec::quad
mesh_quad_rec * quad[6]
Definition: refinedmesh.h:139
oofem::Domain::giveElement
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
oofem::RefinedMesh::fe_quad_rec::node
fe_node_rec * node[4]
Definition: refinedmesh.h:86
oofem::RefinedMesh::mesh_edge_rec::node
fe_node_rec * node[2]
Definition: refinedmesh.h:100
oofem::RefinedMesh::mesh_face_rec::mid_node
fe_node_rec * mid_node
Definition: refinedmesh.h:107
oofem::RefinedMesh::refineMeshGlobally
int refineMeshGlobally(Domain *d, int level, std::vector< RefinedElement > &refinedElementList)
Definition: refinedmesh.C:119
oofem::RefinedMesh::fe_edge_rec
Definition: refinedmesh.h:77
local_tetra_id
#define local_tetra_id(ID)
Definition: refinedmesh.C:88
oofem::IntArray::resize
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
global_hexa_id
#define global_hexa_id(ID)
Definition: refinedmesh.C:82
is_face
#define is_face(ID)
Definition: refinedmesh.C:64
elem_type
#define elem_type(ID)
Definition: refinedmesh.C:57
oofem::RefinedMesh::fe_face_rec::node
fe_node_rec * node[3]
Definition: refinedmesh.h:82
oofem::RefinedMesh::fine_edge_rec::fine_id
int * fine_id
Definition: refinedmesh.h:61
local_hexa_id
#define local_hexa_id(ID)
Definition: refinedmesh.C:89
oofem::RefinedMesh::fe_node_rec
Definition: refinedmesh.h:73
FACE_ELEM
#define FACE_ELEM
Definition: refinedmesh.C:43
oofem::RefinedMesh::mesh_quad_rec
Definition: refinedmesh.h:111
oofem::RefinedMesh::tmp_tetra_rec::ngb_elem_id
int ngb_elem_id[4]
Definition: refinedmesh.h:133
is_hexa
#define is_hexa(ID)
Definition: refinedmesh.C:67
oofem::RefinedMesh::fine_hexa_rec::fine_id
int * fine_id
Definition: refinedmesh.h:69
oofem::RefinedMesh::fe_tetra_rec
Definition: refinedmesh.h:89
oofem::RefinedMesh::mesh_face_rec::fine_id
int * fine_id[3]
Definition: refinedmesh.h:108
oofem::RefinedMesh::tmp_tetra_rec
Definition: refinedmesh.h:132
oofem::Element::giveGeometryType
virtual Element_Geometry_Type giveGeometryType() const
Returns the element geometry type.
Definition: element.C:1529
oofem::RefinedMesh::tmp_hexa_rec::ngb_elem_id
int ngb_elem_id[6]
Definition: refinedmesh.h:138
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::FEMComponent::giveNumber
int giveNumber() const
Definition: femcmpnn.h:107
oofem::RefinedMesh::fe_node_rec::id
int id
Definition: refinedmesh.h:74
oofem::Element::giveNode
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
oofem::RefinedMesh::mesh_face_rec::node
fe_node_rec * node[3]
Definition: refinedmesh.h:106
oofem::RefinedMesh::tmp_hexa_rec
Definition: refinedmesh.h:137
error_message
#define error_message(MSG)
Definition: refinedmesh.C:96
oofem::RefinedMesh::mesh_edge_rec
Definition: refinedmesh.h:99
oofem::RefinedMesh::fe_face_rec
Definition: refinedmesh.h:81
oofem::RefinedElement::giveBoundaryFlagArray
IntArray * giveBoundaryFlagArray(void)
Definition: refinedelement.h:61
oofem::RefinedMesh::fe_hexa_rec::node
fe_node_rec * node[8]
Definition: refinedmesh.h:94
EDGE_ELEM
#define EDGE_ELEM
Definition: refinedmesh.C:42
global_tetra_id
#define global_tetra_id(ID)
Definition: refinedmesh.C:81
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011