1
/*============================================================================
2
* Structure describing a mesh section's subdivision into simpler elements
4
* This is mostly useful to replace polygons or polyhedra by simpler
5
* elements such as triangles, tetrahedra, and prisms upon data export.
6
*============================================================================*/
9
This file is part of Code_Saturne, a general-purpose CFD tool.
11
Copyright (C) 1998-2011 EDF S.A.
13
This program is free software; you can redistribute it and/or modify it under
14
the terms of the GNU General Public License as published by the Free Software
15
Foundation; either version 2 of the License, or (at your option) any later
18
This program is distributed in the hope that it will be useful, but WITHOUT
19
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
23
You should have received a copy of the GNU General Public License along with
24
this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25
Street, Fifth Floor, Boston, MA 02110-1301, USA.
28
/*----------------------------------------------------------------------------*/
33
* For a polygon with n vertices (and thus n edges), a compact way of
34
* storing a tesselation is described here:
35
* The vertices of each of the polygon's triangles are defined in
36
* terms of the polgon's local vertex numbers. As the local number
37
* of resulting triangles should remain small, these numbers may
38
* be represented by bit fields. for example, with 32 bit integers,
39
* we may use bits 1-10, 11-20, and 21-30 for each of a triangle's
40
* vertices (allowing for polygons with up to 1024 vertices, which
41
* vastly exceeds the number of vertices we expect for polygons
42
* (usually in the range of 5-10).
43
* This requires n-2 values per polygon, while the expanded
44
* connectivity requires (n - 2) * 3 values.
46
* As polyhedra are defined by their polygonal faces, this encoding
49
* For a quadrangle, a better way is to use a single diagonal_flag, equal to
50
* O if the quadrangle is split along its diagonal (1, 3), or
51
* 1 if the quadrangle is split along its diagonal (2, 4).
53
* In the fvm_tesselation_t structure, the encoding[] array corresponds
54
* to the opposite_vertex_num[] array (of the same size as vertex_num[])
55
* for polygons and polyhedra, and to a diagonal_flag[] array of size
56
* n_elements for quadrangles.
58
* If FVM model were ever extended to elements with more than 2 vertices
59
* per edge, we would need to replace n_vertices by n_edges in some of
60
* this file's function's logic, but the biggest difficulty would
61
* be in adding extra vertices along new edges without duplicating them:
62
* this would only require some extra processing for a given section,
63
* (including some communication between processors), but would
64
* involve difficulties mainly when using multiple tesselation structures
65
* for multiple mesh sections, as extra vertices along section boundaries
66
* would need to be shared, involving extra book-keeping (and not sharing
67
* those vertices would make the corresponding elements non-conforming,
68
* so it would not be an option). As in most known data models, general
69
* simple polygons or polyhedra only have two vertices per edge,
70
* restricting tesselation to linear element types would not be a problem
71
* as regards its main intended use, which is to allow replacement of
72
* polygons and/or polyhedra upon export to some formats or for use with
73
* some tools which do not support these types of elements.
76
/*----------------------------------------------------------------------------*/
78
#if defined(HAVE_CONFIG_H)
79
#include "cs_config.h"
82
/*----------------------------------------------------------------------------
83
* Standard C library headers
84
*----------------------------------------------------------------------------*/
92
/*----------------------------------------------------------------------------
94
*----------------------------------------------------------------------------*/
96
#include <bft_error.h>
98
#include <bft_printf.h>
100
/*----------------------------------------------------------------------------
102
*----------------------------------------------------------------------------*/
104
#include "fvm_config_defs.h"
105
#include "fvm_defs.h"
106
#include "fvm_io_num.h"
107
#include "fvm_parall.h"
108
#include "fvm_triangulate.h"
110
/*----------------------------------------------------------------------------
111
* Header for the current file
112
*----------------------------------------------------------------------------*/
114
#include "fvm_tesselation.h"
116
/*----------------------------------------------------------------------------*/
121
} /* Fake brace to force Emacs auto-indentation back to column 0 */
123
#endif /* __cplusplus */
125
/*=============================================================================
127
*============================================================================*/
129
/* Geometric primitives */
131
#undef _CROSS_PRODUCT_3D
132
#undef _DOT_PRODUCT_3D
135
#define _CROSS_PRODUCT_3D(cross_v1_v2, v1, v2) ( \
136
cross_v1_v2[0] = v1[1]*v2[2] - v1[2]*v2[1], \
137
cross_v1_v2[1] = v1[2]*v2[0] - v1[0]*v2[2], \
138
cross_v1_v2[2] = v1[0]*v2[1] - v1[1]*v2[0] )
140
#define _DOT_PRODUCT_3D(v1, v2) ( \
141
v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])
143
#define _MODULE_3D(v) \
144
sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])
146
/* Tesselation encoding */
148
#undef _ENCODING_BITS
149
#undef _DECODE_TRIANGLE_VERTICES
151
#define _ENCODING_BITS (sizeof(fvm_tesselation_encoding_t)*8/3)
152
#define _DECODE_TRIANGLE_VERTICES(vertex_list, encoding, decoding_mask) (\
153
vertex_list[0] = encoding & decoding_mask[0], \
154
vertex_list[1] = (encoding & decoding_mask[1]) >> _ENCODING_BITS, \
155
vertex_list[2] = (encoding & decoding_mask[2]) >> (_ENCODING_BITS*2) )
157
/*============================================================================
159
*============================================================================*/
161
typedef unsigned fvm_tesselation_encoding_t;
163
/*----------------------------------------------------------------------------
164
* Structure defining a tesselation of a mesh section.
165
*----------------------------------------------------------------------------*/
167
struct _fvm_tesselation_t {
169
/* Parent section information */
170
/*----------------------------*/
172
fvm_element_t type; /* Element types */
174
fvm_lnum_t n_elements; /* Number of elements */
176
int dim; /* Spatial dimension */
178
int entity_dim; /* Entity dimension */
180
int stride; /* Element size for regular elements
181
(0 for polygons and polyhedra) */
183
fvm_lnum_t n_faces; /* Number of faces defining polyhedra */
185
/* Pointers to shared vertex coordinates */
187
const fvm_coord_t *vertex_coords; /* pointer to vertex coordinates
189
x1, y1, z1, x2, y2, z2, ...) */
191
const fvm_lnum_t *parent_vertex_num; /* Local numbers (1 to n) of local
192
vertices in the parent mesh.
193
This array is present only when non
194
"trivial" (i.e. not 1, 2, ..., n). */
196
/* Pointers to shared connectivity arrays */
198
const fvm_lnum_t *face_index; /* polyhedron -> faces index (O to n-1);
199
dimension [n_elements + 1] */
200
const fvm_lnum_t *face_num; /* polyhedron -> face numbers (1 to n, signed,
201
> 0 for outwards pointing face normal
202
< 0 for inwards pointing face normal);
203
dimension: [face_index[n_elements]] */
205
const fvm_lnum_t *vertex_index; /* polygon face -> vertices index (O to n-1);
206
dimension: [n_cell_faces + 1] */
208
const fvm_lnum_t *vertex_num; /* vertex numbers (1 to n);
209
dimension: connectivity_size */
211
/* Pointer to shared global element numbers */
213
const fvm_io_num_t *global_element_num;
215
/* Basic information */
216
/*-------------------*/
218
int n_sub_types; /* Number of sub-element types
220
fvm_element_t sub_type[2]; /* Sub-element types */
222
fvm_lnum_t n_sub_max[2]; /* Maximum number of sub-elements of
223
each type per element */
224
fvm_lnum_t n_sub_max_glob[2]; /* Global maximum number of
225
sub-elements of each type
227
fvm_lnum_t n_sub[2]; /* Number of sub-elements of
229
fvm_gnum_t n_sub_glob[2]; /* Global number of sub-elements
232
const fvm_tesselation_encoding_t *encoding; /* Compact tesselation
236
fvm_tesselation_encoding_t *_encoding; /* encoding if owner,
239
const fvm_lnum_t *sub_elt_index[2]; /* index of sub-elements of each
240
given type associated with
241
each element (0 to n-1); */
242
fvm_lnum_t *_sub_elt_index[2]; /* sub_elt_index if owner,
247
/*============================================================================
248
* Static global variables
249
*============================================================================*/
251
/*============================================================================
252
* Private function definitions
253
*============================================================================*/
255
/*----------------------------------------------------------------------------
256
* Compute coordinates of added vertices for a tesselation of a polyhedron.
259
* ts <-- tesselation structure
260
* vertex_coords --> coordinates of added vertices
261
* n_vertices_tot --> total number of vertex references (or NULL)
262
* element_id <-- element index
263
*----------------------------------------------------------------------------*/
266
_added_vertex_coords(const fvm_tesselation_t *ts,
267
fvm_coord_t vertex_coords[3],
269
fvm_lnum_t element_id)
271
fvm_lnum_t n_vertices;
272
fvm_lnum_t i, j, k, face_id, vertex_id;
274
fvm_lnum_t _n_vertices_tot = 0;
275
fvm_coord_t cell_center[3] = {0., 0., 0.};
276
fvm_coord_t cell_center_divisor = 0.;
278
/* Loop on polyhedron's faces */
279
/*----------------------------*/
281
for (i = ts->face_index[element_id];
282
i < ts->face_index[element_id + 1];
285
fvm_coord_t sign, v1[3], v2[3];
287
fvm_coord_t vertices_center[3] = {0., 0., 0.};
288
fvm_coord_t face_center[3] = {0., 0., 0.};
289
fvm_coord_t face_normal[3] = {0., 0., 0.};
290
fvm_coord_t triangle_center[3] = {0., 0., 0.};
291
fvm_coord_t triangle_normal[3] = {0., 0., 0.};
292
fvm_coord_t face_surface = 0.;
293
fvm_coord_t triangle_surface = 0.;
294
const fvm_coord_t *current_coords = NULL;
296
face_id = FVM_ABS(ts->face_num[i]) - 1;
298
n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id];
300
_n_vertices_tot += n_vertices;
302
/* First loop on face vertices to estimate face center location */
304
for (j = 0; j < n_vertices; j++) {
306
vertex_id = ts->vertex_num[ts->vertex_index[face_id] + j] - 1;
308
if (ts->parent_vertex_num != NULL)
310
= ts->vertex_coords + ((ts->parent_vertex_num[vertex_id] - 1) * 3);
312
current_coords = ts->vertex_coords + (vertex_id * 3);
314
for (k = 0; k < 3; k++)
315
vertices_center[k] += current_coords[k];
319
for (k = 0; k < 3; k++)
320
vertices_center[k] /= (fvm_coord_t)n_vertices;
322
/* At this stage, current_coords points
323
to the last face vertice's coords */
325
for (k = 0; k < 3; k++) {
326
v1[k] = current_coords[k] - vertices_center[k];
327
triangle_center[k] = current_coords[k] + vertices_center[k];
330
/* Second loop on face vertices to estimate face normal */
332
for (j = 0; j < n_vertices; j++) {
334
vertex_id = ts->vertex_num[ts->vertex_index[face_id] + j] - 1;
336
if (ts->parent_vertex_num != NULL)
338
= ts->vertex_coords + ((ts->parent_vertex_num[vertex_id] - 1) * 3);
340
current_coords = ts->vertex_coords + (vertex_id * 3);
342
for (k = 0; k < 3; k++) {
343
v2[k] = current_coords[k] - vertices_center[k];
344
triangle_center[k] = (triangle_center[k] + current_coords[k])
348
_CROSS_PRODUCT_3D(triangle_normal, v1, v2);
350
for (k = 0; k < 3; k++)
351
face_normal[k] += triangle_normal[k];
353
triangle_surface = _MODULE_3D(triangle_normal) * 0.5;
355
if (_DOT_PRODUCT_3D(triangle_normal, face_normal) > 0)
360
/* Add triangle contributions to face center and normal */
362
face_surface += triangle_surface * sign;
363
for (k = 0; k < 3; k++)
364
face_center[k] += triangle_center[k] * triangle_surface * sign;
366
/* Prepare for next face vertex */
368
for (k = 0; k < 3; k++) {
370
triangle_center[k] = current_coords[k] + vertices_center[k];
373
} /* End of second loop on face vertices */
375
/* If first normal was not oriented like the face normal, correct sign */
377
if (face_surface < 0) {
378
face_surface = - face_surface;
379
for (k = 0; k < 3; k++)
380
face_center[k] = -face_center[k];
383
/* We should now divide the face_center[] values by face_surface to
384
finish, but the contribution of face center to element center is:
385
face_center[] * face_surface;
386
which the face_center[] array contains, so we use those directly */
388
cell_center_divisor += face_surface;
389
for (k = 0; k < 3; k++)
390
cell_center[k] += face_center[k];
392
} /* End of loop on polyhedron's faces */
394
for (k = 0; k < 3; k++)
395
vertex_coords[k] = cell_center[k] / cell_center_divisor;
397
if (n_vertices_tot != NULL)
398
*n_vertices_tot = _n_vertices_tot;
401
/*----------------------------------------------------------------------------
402
* Solve Ax = B for a 4x4 system.
406
* b <-- right hand side
407
* x --> solution of Ax = b
409
* returns: 0 if solution was found, 1 in case of error (zero pivot)
410
*----------------------------------------------------------------------------*/
413
_solve_ax_b_4(double a[4][4],
417
int i, j, k, k_pivot;
419
double abs_pivot, abs_a_ki, factor;
420
double _a[4][4], _b[4];
421
double _a_swap[4], _b_swap;
423
const double _epsilon = 1.0e-15;
425
/* Copy array and RHS (avoid modifying a and b) */
427
for (i = 0; i < 4; i++) {
428
for (j = 0; j < 4; j++) {
434
/* Forward elimination */
436
for (i = 0; i < 4-1; i++) {
438
/* Seek best pivot */
441
abs_pivot = FVM_ABS(_a[i][i]);
443
for (k = i+1; k < 4; k++) {
444
abs_a_ki = FVM_ABS(_a[k][i]);
445
if (abs_a_ki > abs_pivot) {
446
abs_pivot = abs_a_ki;
451
/* Swap if necessary */
455
for (j = 0; j < 4; j++) {
456
_a_swap[j] = _a[i][j];
457
_a[i][j] = _a[k_pivot][j];
458
_a[k_pivot][j] = _a_swap[j];
463
_b[k_pivot] = _b_swap;
467
if (abs_pivot < _epsilon)
470
/* Eliminate values */
472
for (k = i+1; k < 4; k++) {
474
factor = _a[k][i] / _a[i][i];
477
for (j = i+1; j < 4; j++) {
478
_a[k][j] -= _a[i][j]*factor;
480
_b[k] -= _b[i]*factor;
488
x[3] = _b[3] / _a[3][3];
489
x[2] = (_b[2] - _a[2][3]*x[3]) / _a[2][2];
490
x[1] = (_b[1] - _a[1][3]*x[3] - _a[1][2]*x[2]) / _a[1][1];
491
x[0] = (_b[0] - _a[0][3]*x[3] - _a[0][2]*x[2] - _a[0][1]*x[1]) / _a[0][0];
496
/*----------------------------------------------------------------------------
497
* Compute coordinates of added vertices for a tesselation of a polyhedron.
499
* The heap and vertex_list arrays should be large enough to contain
500
* all vertex references of all the polyhedron's faces.
503
* ts <-- tesselation structure
504
* heap --- working array
505
* vertex_list --> list of polyhedron's vertices
506
* vertex_list_size --> list size
507
* element_id <-- element index
508
*----------------------------------------------------------------------------*/
511
_polyhedron_vertices(const fvm_tesselation_t *ts,
513
fvm_lnum_t vertex_list[],
514
int *vertex_list_size,
515
fvm_lnum_t element_id)
517
fvm_lnum_t n_vertices;
518
fvm_lnum_t i, j, face_id, vertex_id;
522
/* Loop on polyhedron's faces */
523
/*----------------------------*/
525
for (i = ts->face_index[element_id];
526
i < ts->face_index[element_id + 1];
529
/* Loop on face vertices */
530
/*-----------------------*/
532
face_id = FVM_ABS(ts->face_num[i]) - 1;
534
n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id];
536
for (j = 0; j < n_vertices; j++) {
540
int hole = heap_size;
541
int parent = (hole - 1) / 2;
543
vertex_id = ts->vertex_num[ts->vertex_index[face_id] + j] - 1;
545
while (hole > 0 && heap[parent] > vertex_id) {
546
heap[hole] = heap[parent];
548
parent = (parent - 1) / 2;
551
heap[hole] = vertex_id;
558
/* Now remove entries from heap */
563
vertex_list[i++] = heap[0];
565
while (heap_size > 0) {
567
if (heap[0] != vertex_list[i-1])
568
vertex_list[i++] = heap[0];
572
/* Insert the last item in the heap whose root is empty */
576
int child = 1; /* left child of 0 */
578
vertex_id = heap[heap_size];
580
while (child < heap_size) {
582
if (child < (heap_size - 1) && heap[child] > heap[child+1])
583
child++; /* child is the smaller child of start */
585
if (vertex_id < heap[child])
586
break; /* item stays in start */
589
heap[start] = heap[child];
596
heap[start] = vertex_id;
602
*vertex_list_size = i;
605
/*----------------------------------------------------------------------------
606
* Compute field values at added vertices for a tesselation of polyhedra.
608
* One additional vertex is added near the center of each polyhedra.
609
* For element types other than polyhedra, there is no need for added
610
* vertices, so this function returns immediately.
613
* this_tesselation <-- tesselation structure
614
* vertex_coords <-- coordinates of added vertices
615
* src_dim <-- dimension of source data
616
* src_dim_shift <-- source data dimension shift (start index)
617
* dest_dim <-- destination data dimension (1 if non interlaced)
618
* start_id <-- added vertices start index
619
* end_id <-- added vertices past the end index
620
* src_interlace <-- indicates if source data is interlaced
621
* src_datatype <-- source data type (float, double, or int)
622
* dest_datatype <-- destination data type (float, double, or int)
623
* n_parent_lists <-- number of parent lists (if parent_num != NULL)
624
* parent_num_shift <-- parent number to value array index shifts;
625
* size: n_parent_lists
626
* parent_num <-- if n_parent_lists > 0, parent entity numbers
627
* src_data <-- array of source arrays (at least one, with one per
628
* source dimension if non interlaced, times one per
629
* parent list if multiple parent lists, with
630
* x_parent_1, y_parent_1, ..., x_parent_2, ...) order
631
* dest_data --> destination buffer
632
*----------------------------------------------------------------------------*/
635
_vertex_field_of_real_values(const fvm_tesselation_t *this_tesselation,
641
fvm_interlace_t src_interlace,
642
fvm_datatype_t src_datatype,
643
fvm_datatype_t dest_datatype,
645
const fvm_lnum_t parent_num_shift[],
646
const fvm_lnum_t parent_num[],
647
const void *const src_data[],
648
void *const dest_data)
650
int n_vertices_tot, pl, vertex_list_size;
651
fvm_lnum_t i, j, parent_id, src_id, vertex_id;
653
int _src_dim_shift = src_dim_shift;
654
int max_list_size = 128;
656
fvm_lnum_t _heap[128];
657
fvm_lnum_t _vertex_list[128];
658
fvm_lnum_t *heap = _heap;
659
fvm_lnum_t *vertex_list = _vertex_list;
661
const fvm_tesselation_t *ts = this_tesselation;
662
const fvm_coord_t *current_coords = NULL;
664
if (ts->type != FVM_CELL_POLY)
667
/* Main loop on polyhedra */
668
/*------------------------*/
670
for (i = start_id; i < end_id; i++) {
672
for (j = 0; j < dest_dim; j++) { /* Loop on destination dimension */
674
fvm_coord_t vertex_coords[3];
676
double interpolated_value, v_x, v_y, v_z;
679
double a[4][4] = {{0., 0., 0., 0.},
683
double b[4] = {0., 0., 0., 0.};
685
_added_vertex_coords(ts,
690
/* Allocate or resize working arrays if too small */
692
if (n_vertices_tot > max_list_size) {
698
BFT_REALLOC(heap, max_list_size, fvm_lnum_t);
699
BFT_REALLOC(vertex_list, max_list_size, fvm_lnum_t);
702
/* Obtain list of polyhedron's vertices */
704
_polyhedron_vertices(ts,
710
/* Build matrix for least squares */
712
for (j = 0; j < vertex_list_size; j++) {
714
vertex_id = vertex_list[j];
716
if (ts->parent_vertex_num != NULL)
718
= ts->vertex_coords + ((ts->parent_vertex_num[vertex_id] - 1) * 3);
720
current_coords = ts->vertex_coords + (vertex_id * 3);
722
v_x = current_coords[0];
723
v_y = current_coords[1];
724
v_z = current_coords[2];
726
/* Position in source data for current vertex value */
728
if (n_parent_lists == 0) {
730
parent_id = vertex_id;
732
else if (parent_num != NULL) {
733
for (parent_id = parent_num[vertex_id] - 1, pl = n_parent_lists - 1;
734
parent_id < parent_num_shift[pl];
737
parent_id -= parent_num_shift[pl];
740
for (parent_id = vertex_id, pl = n_parent_lists - 1 ;
741
parent_id < parent_num_shift[pl] ;
744
parent_id -= parent_num_shift[pl];
747
if (src_interlace == FVM_INTERLACE) {
748
src_id = parent_id * src_dim + _src_dim_shift;
751
pl = src_dim*pl + _src_dim_shift;
755
/* Now convert based on source datatype */
757
switch(src_datatype) {
759
v_f = ((const float *const *const)src_data)[pl][src_id];
762
v_f = ((const double *const *const)src_data)[pl][src_id];
768
a[0][0] += v_x * v_x;
769
a[0][1] += v_x * v_y;
770
a[0][2] += v_x * v_z;
773
a[1][1] += v_y * v_y;
774
a[1][2] += v_y * v_z;
777
a[2][2] += v_z * v_z;
787
} /* End of loop on element vertices */
789
/* Matrix is symmetric */
800
/* Find hyperplane coefficients and compute added value */
802
if (_solve_ax_b_4(a, b, coeff) == 0) {
803
interpolated_value = ( coeff[0] * vertex_coords[0]
804
+ coeff[1] * vertex_coords[1]
805
+ coeff[2] * vertex_coords[2]
809
interpolated_value = v_f; /* last encountered value */
812
/* Now convert based on destination datatype */
814
switch(dest_datatype) {
816
((float *const)dest_data)[i] = interpolated_value;
819
((double *const)dest_data)[i] = interpolated_value;
825
} /* End of loop on destination dimension */
827
} /* End of main loop on polyhedra */
831
BFT_FREE(vertex_list);
835
/*----------------------------------------------------------------------------
836
* Initialize decoding mask for polygon triangulations.
839
* decoding_mask <-- decoding mask for triangle
840
*----------------------------------------------------------------------------*/
843
_init_decoding_mask(fvm_tesselation_encoding_t decoding_mask[3])
847
for (i = 0; i < _ENCODING_BITS; i++)
848
decoding_mask[0] = (decoding_mask[0] << 1) + 1;
849
decoding_mask[1] = decoding_mask[0] << _ENCODING_BITS;
850
decoding_mask[2] = decoding_mask[1] << _ENCODING_BITS;
853
/*----------------------------------------------------------------------------
854
* Tesselate polygons (2D elements or 3D element faces) of a mesh section
855
* referred to by an fvm_tesselation_t structure. Quadrangles are not
859
* this_tesselation <-> partially initialized tesselation structure
860
* dim <-- spatial dimension
861
* vertex_coords <-- associated vertex coordinates array
862
* parent_vertex_num <-- optional indirection to vertex coordinates
863
* error_count --> number of triangulation errors counter (optional)
864
*----------------------------------------------------------------------------*/
867
_tesselate_polygons(fvm_tesselation_t *this_tesselation,
869
const fvm_coord_t vertex_coords[],
870
const fvm_lnum_t parent_vertex_num[],
871
fvm_lnum_t *error_count)
874
fvm_lnum_t n_vertices, n_triangles, n_quads, n_elements;
875
fvm_lnum_t n_vertices_max, n_triangles_max;
876
fvm_lnum_t i, j, k, vertex_id, encoding_id;
877
fvm_tesselation_encoding_t encoding_sub[3];
879
fvm_gnum_t n_g_elements_tot[2] = {0, 0}; /* Global new elements count */
880
fvm_lnum_t n_elements_tot[2] = {0, 0}; /* New triangles/quadrangles */
881
fvm_lnum_t n_g_elements_max[2] = {0, 0}; /* Global max triangles/quadrangles */
882
fvm_lnum_t n_elements_max[2] = {0, 0}; /* Max triangles/quadrangles */
883
fvm_lnum_t *triangle_vertices = NULL;
884
fvm_triangulate_state_t *state = NULL;
886
fvm_tesselation_t *ts = this_tesselation;
890
if (error_count != NULL)
893
if (ts->type != FVM_CELL_POLY)
894
n_elements = ts->n_elements;
896
n_elements = ts->n_faces;
898
/* Count expected local numbers of triangles and quadrangles */
903
if (ts->vertex_index != NULL) {
904
for (i = 0 ; i < n_elements ; i++) {
905
n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i];
907
n_elements_tot[1] += 1;
909
n_elements_tot[0] += n_vertices - 2;
910
if (n_vertices > n_vertices_max)
911
n_vertices_max = n_vertices;
917
n_triangles_max = n_vertices_max - 2;
919
if (n_vertices_max > 2<<_ENCODING_BITS)
920
bft_error(__FILE__, __LINE__, 0,
921
_("Encoding of tesselation impossible:\n"
922
"maximum number of vertices per polygon: %u\n"
923
"maximum integer encoded on %d/3 = %d bits: %ld\n"),
924
(unsigned)n_triangles_max,
925
sizeof(fvm_tesselation_encoding_t)*8,
926
_ENCODING_BITS, (long)(2<<_ENCODING_BITS));
928
/* Destroy previous tesselation description if present */
931
if (ts->_encoding != NULL)
932
BFT_FREE(ts->_encoding);
934
/* Allocate memory and state variables */
935
/*-------------------------------------*/
937
if (n_vertices_max > 4) {
938
BFT_MALLOC(ts->_encoding,
939
ts->vertex_index[n_elements] - n_elements*2,
940
fvm_tesselation_encoding_t);
941
ts->encoding = ts->_encoding;
942
BFT_MALLOC(triangle_vertices, (n_vertices_max - 2) * 3, fvm_lnum_t);
943
state = fvm_triangulate_state_create(n_vertices_max);
946
n_elements_tot[0] = 0; n_elements_tot[1] = 0; /* reset */
948
/* Main loop on section face elements */
949
/*------------------------------------*/
951
for (i = 0 ; i < n_elements ; i++) {
955
n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i];
956
vertex_id = ts->vertex_index[i];
958
/* We calculate the encoding index base from the polygon's
959
connectivity index base, knowing that for a polygon
960
with n vertices, we have at most n-2 triangles
961
(exactly n-2 when no error occurs) */
963
encoding_id = ts->vertex_index[i] - (i*2);
970
/* If face must be subdivided */
972
if (n_vertices > 4) {
974
n_triangles = fvm_triangulate_polygon(dim,
980
FVM_TRIANGULATE_ELT_DEF,
984
if (n_triangles != (n_vertices - 2) && error_count != NULL)
987
/* Encode local triangle connectivity */
989
for (j = 0; j < n_triangles; j++) {
991
for (k = 0; k < 3; k++)
993
= ( ((fvm_tesselation_encoding_t)(triangle_vertices[j*3 + k] - 1))
994
<< (_ENCODING_BITS * k));
996
ts->_encoding[encoding_id + j]
997
= encoding_sub[0] | encoding_sub[1] | encoding_sub[2];
1001
/* In case of incomplete tesselation due to errors,
1002
blank unused encoding values */
1004
for (j = n_triangles; j < (n_vertices - 2); j++)
1005
ts->_encoding[encoding_id + j] = 0;
1007
n_elements_tot[0] += n_triangles;
1011
/* Otherwise, tesselation trivial or not necessary for this face */
1015
if (ts->_encoding != NULL) {
1016
for (j = 0; j < (n_vertices - 2); j++)
1017
ts->_encoding[encoding_id + j] = 0;
1020
if (n_vertices == 4) {
1021
n_elements_tot[1] += 1;
1025
else if (n_vertices == 3) {
1026
n_elements_tot[0] += 1;
1032
if (n_triangles > n_elements_max[0])
1033
n_elements_max[0] = n_triangles;
1035
if (n_quads > n_elements_max[1])
1036
n_elements_max[1] = n_triangles;
1038
} /* End of loop on elements */
1040
/* Free memory and state variables */
1042
if (n_vertices_max > 4) {
1043
BFT_FREE(triangle_vertices);
1044
state = fvm_triangulate_state_destroy(state);
1047
/* Update tesselation structure info */
1049
for (type_id = 0; type_id < 2; type_id++) {
1050
n_g_elements_tot[type_id] = n_elements_tot[type_id];
1051
n_g_elements_max[type_id] = n_elements_max[type_id];
1054
fvm_parall_counter(n_g_elements_tot, 2);
1055
fvm_parall_counter_max(n_g_elements_max, 2);
1057
for (type_id = 0; type_id < 2; type_id++) {
1058
if (n_g_elements_tot[type_id] > 0) {
1060
ts->sub_type[ts->n_sub_types] = FVM_FACE_TRIA;
1062
ts->sub_type[ts->n_sub_types] = FVM_FACE_QUAD;
1063
ts->n_sub_max[ts->n_sub_types] = n_elements_max[type_id];
1064
ts->n_sub_max_glob[ts->n_sub_types] = n_g_elements_max[type_id];
1065
ts->n_sub[ts->n_sub_types] = n_elements_tot[type_id];
1066
ts->n_sub_glob[ts->n_sub_types] = n_elements_tot[type_id];
1067
ts->n_sub_types += 1;
1073
/*----------------------------------------------------------------------------
1074
* Build indexes and global counts for polygons of a mesh section referred
1075
* to by an fvm_tesselation_t structure.
1078
* this_tesselation <-> partially initialized tesselation structure
1079
* global_count <-- indicates if global counts should be built
1080
*----------------------------------------------------------------------------*/
1083
_count_and_index_sub_polygons(fvm_tesselation_t *this_tesselation,
1086
int sub_type_id, type_id;
1087
fvm_lnum_t n_vertices, n_triangles, n_elements;
1088
fvm_lnum_t i, j, encoding_id;
1090
fvm_lnum_t *n_sub_elements[2] = {NULL, NULL};
1092
fvm_tesselation_t *ts = this_tesselation;
1094
/* Initial checks */
1096
if (ts->vertex_index == NULL)
1099
/* Initialization */
1101
n_elements = ts->n_elements;
1103
/* Count expected local numbers of triangles and quadrangles */
1105
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) {
1108
switch(ts->sub_type[sub_type_id]) {
1119
BFT_MALLOC(ts->_sub_elt_index[sub_type_id], n_elements + 1, fvm_lnum_t);
1121
for (i = 0 ; i < n_elements + 1 ; i++)
1122
ts->_sub_elt_index[sub_type_id][i] = 0;
1124
/* The index array will first be used to hold n_sub_elements[],
1125
of size n_elements. To simplify future conversion to index,
1126
we shift shuch that n_sub_elements[i] corresponds to index[i+1]. */
1128
n_sub_elements[type_id] = ts->_sub_elt_index[sub_type_id] + 1;
1132
/* Loop on elements to populate n_sub_elements[].
1133
Note that each n_sub_elements[] array has been initialized with zeroes,
1134
as it is mapped to a ts->sub_elt_index[] thus initialized. */
1136
for (i = 0 ; i < n_elements ; i++) {
1138
n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i];
1141
if (n_vertices == 3) {
1142
n_sub_elements[0][i] = 1;
1146
else { /* if (n_vertices > 3) */
1148
encoding_id = ts->vertex_index[i] - (i*2);
1150
for (j = 0; j < (n_vertices - 2); j++) {
1151
if (ts->encoding != NULL) {
1152
if (ts->encoding[encoding_id + j] != 0)
1156
if (n_triangles > 0)
1157
n_sub_elements[0][i] = n_triangles;
1158
else if (n_vertices == 4)
1159
n_sub_elements[1][i] = 1;
1160
assert(n_triangles > 0 || n_vertices == 4);
1165
/* Now compute max global number if necessary */
1167
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++)
1168
ts->n_sub_glob[sub_type_id] = ts->n_sub[sub_type_id];
1170
if (global_count == true && ts->global_element_num != NULL) {
1172
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) {
1173
if (ts->_sub_elt_index[sub_type_id] != NULL) {
1174
fvm_lnum_t * _n_sub_elements = ts->_sub_elt_index[sub_type_id] + 1;
1175
if (n_elements == 0) _n_sub_elements = NULL;
1176
ts->n_sub_glob[sub_type_id]
1177
= fvm_io_num_global_sub_size(ts->global_element_num,
1181
ts->n_sub_glob[sub_type_id] = 0;
1186
/* Finally, build index from n_sub_elements_glob */
1188
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) {
1190
fvm_lnum_t *sub_elt_index = ts->_sub_elt_index[sub_type_id];
1191
for (i = 0 ; i < n_elements ; i++)
1192
sub_elt_index[i+1] = sub_elt_index[i] + sub_elt_index[i+1];
1193
ts->sub_elt_index[sub_type_id] = ts->_sub_elt_index[sub_type_id];
1199
/*----------------------------------------------------------------------------
1200
* Count resulting number of tetrahedra and pyramids when tesselating
1201
* polyhedra of a mesh section referred to by an fvm_tesselation_t
1202
* structure. This requires that the mesh section's polygonal faces have
1203
* already been tesselated (i.e. _tesselate_polygons has been called).
1206
* this_tesselation <-> partially initialized tesselation structure
1207
* error_count --> number of tesselation errors counter (optional)
1208
* global_count <-- indicates if global counts should be built
1209
*----------------------------------------------------------------------------*/
1212
_count_and_index_sub_polyhedra(fvm_tesselation_t *this_tesselation,
1213
fvm_lnum_t *error_count,
1216
int sub_type_id, type_id;
1217
fvm_lnum_t n_vertices, n_triangles, n_elements;
1218
fvm_lnum_t n_tetras, n_pyrams;
1219
fvm_lnum_t i, j, k, face_id, encoding_id;
1221
fvm_gnum_t n_g_elements_tot[2] = {0, 0}; /* Global new elements count */
1222
fvm_lnum_t n_elements_tot[2] = {0, 0}; /* New tetrahedra/pyramids */
1223
fvm_lnum_t *n_sub_elements[2] = {NULL, NULL};
1224
fvm_lnum_t n_g_elements_max[2] = {0, 0}; /* Global max tetrahedra/pyramids */
1225
fvm_lnum_t n_elements_max[2] = {0, 0}; /* Max tetrahedra/pyramids */
1227
fvm_tesselation_t *ts = this_tesselation;
1229
/* Initial checks */
1231
assert(ts->vertex_index != NULL);
1233
/* Initialization */
1235
if (error_count != NULL)
1238
n_elements = ts->n_elements;
1240
/* Count expected total and local numbers of tetrahedra and pyramids;
1241
at this stage, the tesselation has been initialized as a
1242
polygon tesselation; adding a "center" vertex, triangle faces
1243
lead to tetrahedra, and quadrangles to pyramids */
1245
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) {
1248
switch(ts->sub_type[sub_type_id]) {
1259
BFT_MALLOC(ts->_sub_elt_index[sub_type_id], n_elements + 1, fvm_lnum_t);
1261
for (i = 0 ; i < n_elements + 1 ; i++)
1262
ts->_sub_elt_index[sub_type_id][i] = 0;
1264
/* The index array will first be used to hold n_sub_elements[],
1265
of size n_elements. To simplify future conversion to index,
1266
we shift shuch that n_sub_elements[i] corresponds to index[i+1]. */
1268
n_sub_elements[type_id] = ts->_sub_elt_index[sub_type_id] + 1;
1272
/* Counting loop on polyhedra elements */
1273
/*-------------------------------------*/
1275
for (i = 0 ; i < n_elements ; i++) {
1280
for (j = ts->face_index[i]; /* Loop on element faces */
1281
j < ts->face_index[i+1];
1284
face_id = FVM_ABS(ts->face_num[j]) - 1;
1286
n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id];
1288
if (n_vertices == 3)
1291
else { /* if (n_vertices > 3) */
1295
encoding_id = ts->vertex_index[face_id] - (face_id*2);
1297
for (k = encoding_id; k < (encoding_id + n_vertices - 2); k++) {
1298
if (ts->encoding != NULL) {
1299
if (ts->encoding[k] != 0)
1304
if (error_count != NULL && n_triangles < n_vertices - 2)
1307
if (n_triangles > 0)
1308
n_tetras += n_triangles;
1309
else if (n_vertices == 4)
1314
} /* End of loop on element faces */
1316
n_elements_tot[0] += n_tetras;
1317
n_elements_tot[1] += n_pyrams;
1319
if (n_tetras > n_elements_max[0])
1320
n_elements_max[0] = n_tetras;
1322
if (n_pyrams > n_elements_max[1])
1323
n_elements_max[1] = n_pyrams;
1325
if (n_sub_elements[0] != NULL)
1326
n_sub_elements[0][i] = n_tetras;
1328
if (n_sub_elements[1] != NULL)
1329
n_sub_elements[1][i] = n_pyrams;
1331
} /* End of loop on elements */
1333
/* Update tesselation structure info */
1335
for (type_id = 0; type_id < 2; type_id++) {
1336
n_g_elements_tot[type_id] = n_elements_tot[type_id];
1337
n_g_elements_max[type_id] = n_elements_max[type_id];
1340
fvm_parall_counter(n_g_elements_tot, 2);
1341
fvm_parall_counter_max(n_g_elements_max, 2);
1343
ts->n_sub_types = 0;
1345
for (type_id = 0; type_id < 2; type_id++) {
1346
if (n_g_elements_tot[type_id] > 0) {
1348
ts->sub_type[ts->n_sub_types] = FVM_CELL_TETRA;
1350
ts->sub_type[ts->n_sub_types] = FVM_CELL_PYRAM;
1351
ts->n_sub_max[ts->n_sub_types] = n_elements_max[type_id];
1352
ts->n_sub_max_glob[ts->n_sub_types] = n_g_elements_max[type_id];
1353
ts->n_sub[ts->n_sub_types] = n_elements_tot[type_id];
1354
ts->n_sub_glob[ts->n_sub_types] = n_elements_tot[type_id];
1355
ts->n_sub_types += 1;
1359
/* Now compute max global number if necessary */
1361
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++)
1362
ts->n_sub_glob[sub_type_id] = ts->n_sub[sub_type_id];
1364
if (global_count == true && ts->global_element_num != NULL) {
1366
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) {
1367
if (ts->_sub_elt_index[sub_type_id] != NULL) {
1368
fvm_lnum_t * _n_sub_elements = ts->_sub_elt_index[sub_type_id] + 1;
1369
if (n_elements == 0) _n_sub_elements = NULL;
1370
ts->n_sub_glob[sub_type_id]
1371
= fvm_io_num_global_sub_size(ts->global_element_num,
1375
ts->n_sub_glob[sub_type_id] = 0;
1380
/* Finally, build index from n_sub_elements_glob */
1382
for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) {
1384
fvm_lnum_t *sub_elt_index = ts->_sub_elt_index[sub_type_id];
1385
for (i = 0 ; i < n_elements ; i++)
1386
sub_elt_index[i+1] = sub_elt_index[i] + sub_elt_index[i+1];
1387
ts->sub_elt_index[sub_type_id] = ts->_sub_elt_index[sub_type_id];
1393
#if defined(HAVE_MPI)
1395
/*----------------------------------------------------------------------------
1396
* Update global_num_end when limiting partial connectivity range end index.
1399
* this_tesselation <-- tesselation structure
1400
* end_id <-- initial end of subset in parent section
1401
* global_num_end <-> past the end (maximum + 1) parent element
1402
* global number (reduced on return if required
1404
* comm <-- associated MPI communicator
1407
* index corresponding to end of decoded range
1408
*----------------------------------------------------------------------------*/
1411
_expand_limit_g(const fvm_tesselation_t *this_tesselation,
1413
fvm_gnum_t *global_num_end,
1416
fvm_gnum_t local_max, global_max;
1418
const fvm_gnum_t *global_element_num
1419
= fvm_io_num_get_global_num(this_tesselation->global_element_num);
1421
/* Check if the maximum id returned on some ranks leads to
1422
a lower global_num_end than initially required
1423
(due to local buffer being full) */
1425
if (end_id < this_tesselation->n_elements) {
1426
if (global_element_num != NULL)
1427
local_max = global_element_num[end_id];
1429
local_max = end_id + 1;
1432
local_max = *global_num_end;
1434
MPI_Allreduce(&local_max, &global_max, 1, FVM_MPI_GNUM, MPI_MIN, comm);
1436
if (global_max < *global_num_end)
1437
*global_num_end = global_max;
1440
#endif /* defined(HAVE_MPI) */
1442
#if defined(HAVE_MPI)
1444
/*----------------------------------------------------------------------------
1445
* Decode polygons tesselation to a connectivity buffer.
1447
* To avoid requiring huge buffers and computing unneeded element
1448
* connectivities when exporting data in slices, this function may decode
1449
* a partial connectivity range, starting at polygon index start_id and ending
1450
* either when the indicated buffer size is attained, or the global element
1451
* number corresponding to a given polygon exceeds a given value.
1452
* It returns the effective polygon index end.
1455
* this_tesselation <-- tesselation structure
1456
* connect_type <-- destination element type
1457
* start_id <-- start index of polygons subset in parent section
1458
* buffer_limit <-- maximum number of sub-elements of destination
1459
* element type allowable for sub_element_idx[]
1460
* and vertex_num[] buffers
1461
* global_num_end <-- past the end (maximum + 1) parent element
1463
* global_vertex_num <-- global vertex numbering
1464
* vertex_num --> sub-element (global) vertex connectivity
1467
* polygon index end corresponding to decoded range
1468
*----------------------------------------------------------------------------*/
1471
_decode_polygons_tesselation_g(const fvm_tesselation_t *this_tesselation,
1472
fvm_element_t connect_type,
1473
fvm_lnum_t start_id,
1474
fvm_lnum_t buffer_limit,
1475
fvm_gnum_t global_num_end,
1476
const fvm_io_num_t *global_vertex_num,
1477
fvm_gnum_t vertex_num[])
1479
fvm_lnum_t n_vertices;
1480
fvm_lnum_t i, j, k, l, vertex_id, encoding_id;
1482
fvm_tesselation_encoding_t decoding_mask[3] = {0, 0, 0};
1484
fvm_lnum_t n_sub_tot = 0;
1485
const fvm_tesselation_t *ts = this_tesselation;
1486
fvm_lnum_t retval = start_id;
1488
const fvm_gnum_t *_global_vertex_num
1489
= fvm_io_num_get_global_num(global_vertex_num);
1490
const fvm_gnum_t *global_element_num
1491
= fvm_io_num_get_global_num(ts->global_element_num);
1493
_init_decoding_mask(decoding_mask);
1495
/* Main loop on polygons */
1496
/*-----------------------*/
1498
for (i = 0, j = start_id ;
1499
j < this_tesselation->n_elements;
1502
if ( global_element_num != NULL
1503
&& global_element_num[j] >= global_num_end)
1506
n_vertices = ts->vertex_index[j+1] - ts->vertex_index[j];
1507
vertex_id = ts->vertex_index[j];
1508
encoding_id = ts->vertex_index[j] - (j*2);
1510
/* Sub-elements (triangles) connectivity */
1511
/*---------------------------------------*/
1513
if ( connect_type == FVM_FACE_TRIA
1514
&& n_vertices > 3 && ts->encoding != NULL) {
1516
/* Loop on triangles */
1518
for (k = 0; k < (n_vertices - 2); k++) {
1520
if (ts->encoding[encoding_id + k] != 0) {
1522
/* Fill connectivity array */
1523
/* Return previous element's end index if buffer size reached */
1525
if (n_sub_tot >= buffer_limit)
1528
/* Fill connectivity array */
1530
_DECODE_TRIANGLE_VERTICES(tv,
1531
ts->encoding[encoding_id + k],
1534
if (_global_vertex_num != NULL) {
1535
for (l = 0; l < 3; l++)
1536
vertex_num[n_sub_tot*3 + l]
1537
= _global_vertex_num[ts->vertex_num[vertex_id + tv[l]] - 1];
1540
for (l = 0; l < 3; l++)
1541
vertex_num[n_sub_tot*3 + l]
1542
= ts->vertex_num[vertex_id + tv[l]];
1549
} /* End of loop on polygon vertices */
1553
/* Non-tesselated elements connectivity */
1554
/*--------------------------------------*/
1556
else if ( (connect_type == FVM_FACE_TRIA && n_vertices == 3)
1557
|| (connect_type == FVM_FACE_QUAD && n_vertices == 4)) {
1559
/* Check that polygon is not subdivided */
1562
if (ts->encoding != NULL) {
1563
for (k = 0; k < (n_vertices - 2); k++)
1564
if (ts->encoding[encoding_id + k] > 0)
1568
if (k == n_vertices - 2 || ts->encoding == NULL) {
1570
/* Return previous element's end index if buffer size reached */
1572
if (n_sub_tot >= buffer_limit)
1575
/* Fill connectivity array */
1577
if (_global_vertex_num != NULL) {
1578
for (k = 0; k < n_vertices; k++)
1579
vertex_num[n_sub_tot*n_vertices +k]
1580
= _global_vertex_num[ts->vertex_num[vertex_id + k] - 1];
1583
for (k = 0; k < n_vertices; k++)
1584
vertex_num[n_sub_tot*n_vertices +k]
1585
= ts->vertex_num[vertex_id + k];
1592
} /* End of case where polygon is not tesselated */
1594
retval = j+1; /* If end of buffer has not already been reached */
1596
} /* End of loop on polygons */
1601
#endif /* defined(HAVE_MPI) */
1603
/*----------------------------------------------------------------------------
1604
* Decode polygons tesselation to a connectivity buffer.
1606
* To avoid requiring huge buffers and computing unneeded element
1607
* connectivities this function may decode a partial connectivity range,
1608
* starting at polygon index start_id and ending either when the indicated
1609
* buffer size or the last polygon is attained.
1610
* It returns the effective polygon index end.
1613
* this_tesselation <-- tesselation structure
1614
* connect_type <-- destination element type
1615
* start_id <-- start index of polygons subset in parent section
1616
* buffer_limit <-- maximum number of sub-elements of destination
1617
* element type allowable for vertex_num[] buffer
1618
* vertex_num --> sub-element (global) vertex connectivity
1621
* polygon index end corresponding to decoded range
1622
*----------------------------------------------------------------------------*/
1625
_decode_polygons_tesselation_l(const fvm_tesselation_t *this_tesselation,
1626
fvm_element_t connect_type,
1627
fvm_lnum_t start_id,
1628
fvm_lnum_t buffer_limit,
1629
fvm_lnum_t vertex_num[])
1631
fvm_lnum_t n_vertices;
1632
fvm_lnum_t i, j, k, vertex_id, encoding_id;
1634
fvm_tesselation_encoding_t decoding_mask[3] = {0, 0, 0};
1636
fvm_lnum_t n_sub_tot = 0;
1637
const fvm_tesselation_t *ts = this_tesselation;
1638
fvm_lnum_t retval = start_id;
1640
_init_decoding_mask(decoding_mask);
1642
/* Main loop on polygons */
1643
/*-----------------------*/
1645
for (i = start_id ; i < this_tesselation->n_elements; i++) {
1647
n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i];
1648
vertex_id = ts->vertex_index[i];
1649
encoding_id = ts->vertex_index[i] - (i*2);
1651
/* Sub-elements (triangles) connectivity */
1652
/*---------------------------------------*/
1654
if ( connect_type == FVM_FACE_TRIA
1655
&& n_vertices > 3 && ts->encoding != NULL) {
1657
/* Loop on polygon vertices */
1659
for (j = 0; j < (n_vertices - 2); j++) {
1661
if (ts->encoding[encoding_id + j] != 0) {
1663
/* Fill connectivity array */
1664
/* Return previous element's end index if buffer size reached */
1666
if (n_sub_tot >= buffer_limit)
1669
/* Fill connectivity array */
1671
_DECODE_TRIANGLE_VERTICES(tv,
1672
ts->encoding[encoding_id + j],
1675
for (k = 0; k < 3; k++)
1676
vertex_num[n_sub_tot*3 + k] = ts->vertex_num[vertex_id + tv[k]];
1682
} /* End of loop on polygon vertices */
1686
/* Non-tesselated elements connectivity */
1687
/*--------------------------------------*/
1689
else if ( (connect_type == FVM_FACE_TRIA && n_vertices == 3)
1690
|| (connect_type == FVM_FACE_QUAD && n_vertices == 4)) {
1692
/* Check that polygon is not subdivided */
1695
if (ts->encoding != NULL) {
1696
for (j = 0; j < (n_vertices - 2); j++)
1697
if (ts->encoding[encoding_id + j] != 0)
1701
if (j == n_vertices - 2 || ts->encoding == NULL) {
1703
/* Return previous element's end index if buffer size reached */
1705
if (n_sub_tot >= buffer_limit)
1708
/* Fill connectivity array */
1710
for (j = 0; j < n_vertices; j++)
1711
vertex_num[n_sub_tot*n_vertices +j] = ts->vertex_num[vertex_id + j];
1717
} /* End of case where polygon is not tesselated */
1719
retval = i+1; /* If end of buffer has not already been reached */
1721
} /* End of loop on polygons */
1726
#if defined(HAVE_MPI)
1728
/*----------------------------------------------------------------------------
1729
* Decode polyhedra tesselation to a connectivity buffer.
1731
* To avoid requiring huge buffers and computing unneeded element
1732
* connectivities when exporting data in slices, this function may decode a
1733
* partial connectivity range, starting at polyhedron index start_id and ending
1734
* either when the indicated buffer size is attained, or the global element
1735
* number corresponding to a given polyhedron exceeds a given value.
1736
* It returns the effective polyhedron index end.
1739
* this_tesselation <-- tesselation structure
1740
* connect_type <-- destination element type
1741
* start_id <-- start index of polyhedra subset in parent section
1742
* buffer_limit <-- maximum number of sub-elements of destination
1743
* element type allowable for vertex_num[] buffer
1744
* global_num_end <-- past the end (maximum + 1) parent element
1746
* extra_vertex_base <-- starting number for added vertices
1747
* global_vertex_num <-- global vertex numbering
1748
* vertex_num --> sub-element (global) vertex connectivity
1751
* polyhedron index end corresponding to decoded range
1752
*----------------------------------------------------------------------------*/
1755
_decode_polyhedra_tesselation_g(const fvm_tesselation_t *this_tesselation,
1756
fvm_element_t connect_type,
1757
fvm_lnum_t start_id,
1758
fvm_lnum_t buffer_limit,
1759
fvm_gnum_t global_num_end,
1760
fvm_gnum_t extra_vertex_base_num,
1761
const fvm_io_num_t *global_vertex_num,
1762
fvm_gnum_t vertex_num[])
1765
fvm_lnum_t n_vertices;
1766
fvm_lnum_t i, j, k, l, m, base_dest_id, face_id, vertex_id, encoding_id;
1768
fvm_tesselation_encoding_t decoding_mask[3] = {0, 0, 0};
1770
fvm_lnum_t n_sub_tot = 0;
1771
const fvm_tesselation_t *ts = this_tesselation;
1772
fvm_lnum_t retval = start_id;
1774
const fvm_gnum_t *_global_vertex_num
1775
= fvm_io_num_get_global_num(global_vertex_num);
1776
const fvm_gnum_t *global_element_num
1777
= fvm_io_num_get_global_num(ts->global_element_num);
1779
_init_decoding_mask(decoding_mask);
1781
/* Main loop on polyhedra */
1782
/*------------------------*/
1784
for (i = 0, j = start_id ;
1785
j < this_tesselation->n_elements;
1788
if ( global_element_num != NULL
1789
&& global_element_num[j] >= global_num_end)
1792
for (k = ts->face_index[j]; /* Loop on element faces */
1793
k < ts->face_index[j+1];
1796
face_id = FVM_ABS(ts->face_num[k]) - 1;
1798
n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id];
1799
vertex_id = ts->vertex_index[face_id];
1800
encoding_id = ts->vertex_index[face_id] - (face_id*2);
1802
/* Invert base orientation as it points outwards in polyhedron
1803
definition, but inwards in tetrahedron and pyramid reference
1806
if (ts->face_num[k] > 0)
1811
/* Sub-elements (tetrahedra) connectivity */
1812
/*----------------------------------------*/
1814
if ( connect_type == FVM_CELL_TETRA
1815
&& n_vertices > 3 && ts->encoding != NULL) {
1817
/* Loop on face vertices */
1819
for (l = 0; l < (n_vertices - 2); l++) {
1821
if (ts->encoding[encoding_id + l] != 0) {
1823
/* Return previous element's end index if buffer size reached */
1825
if (n_sub_tot >= buffer_limit)
1829
base_dest_id = n_sub_tot*4;
1831
base_dest_id = n_sub_tot*4 + 2;
1833
/* Fill connectivity array */
1835
_DECODE_TRIANGLE_VERTICES(tv,
1836
ts->encoding[encoding_id + l],
1839
if (_global_vertex_num != NULL) {
1840
for (m = 0; m < 3; m++)
1841
vertex_num[base_dest_id + orient*m]
1842
= _global_vertex_num[ts->vertex_num[vertex_id + tv[m]] - 1];
1845
for (m = 0; m < 3; m++)
1846
vertex_num[base_dest_id + orient*m]
1847
= ts->vertex_num[vertex_id + tv[m]];
1850
/* Add extra "top" vertex */
1852
if (global_element_num != NULL)
1853
vertex_num[n_sub_tot*4 + 3] = extra_vertex_base_num
1854
+ global_element_num[j] - 1;
1856
vertex_num[n_sub_tot*4 + 3] = extra_vertex_base_num + j;
1866
/* Non-tesselated faces connectivity */
1867
/*-----------------------------------*/
1869
else if ( (connect_type == FVM_CELL_TETRA && n_vertices == 3)
1870
|| (connect_type == FVM_CELL_PYRAM && n_vertices == 4)) {
1872
/* Check that face is not subdivided */
1875
if (ts->encoding != NULL) {
1876
for (l = 0; l < (n_vertices - 2); l++)
1877
if (ts->encoding[encoding_id + l] != 0)
1881
if (l == n_vertices - 2 || ts->encoding == NULL) {
1883
fvm_lnum_t stride = n_vertices + 1;
1885
/* Return previous element's end index if buffer size reached */
1887
if (n_sub_tot >= buffer_limit)
1891
base_dest_id = n_sub_tot*stride;
1893
base_dest_id = n_sub_tot*stride + n_vertices-1;
1895
/* Fill connectivity array */
1897
if (_global_vertex_num != NULL) {
1898
for (l = 0; l < n_vertices; l++)
1899
vertex_num[base_dest_id + l*orient]
1900
= _global_vertex_num[ts->vertex_num[vertex_id + l] - 1];
1903
for (l = 0; l < n_vertices; l++)
1904
vertex_num[base_dest_id + l*orient]
1905
= ts->vertex_num[vertex_id + l];
1908
/* Add extra "top" vertex */
1910
if (global_element_num != NULL)
1911
vertex_num[n_sub_tot*stride + n_vertices]
1912
= extra_vertex_base_num + global_element_num[j] - 1;
1914
vertex_num[n_sub_tot*stride + n_vertices]
1915
= extra_vertex_base_num + j;
1921
} /* End of case where face is not tesselated */
1923
} /* End of loop on element faces */
1925
retval = j+1; /* If end of buffer has not already been reached */
1927
} /* End of main loop on polyhedra */
1932
#endif /* defined(HAVE_MPI) */
1934
/*----------------------------------------------------------------------------
1935
* Decode polyhedra tesselation to a connectivity buffer.
1937
* To avoid requiring huge buffers, this function may decode a partial
1938
* connectivity range, starting at polyhedron index start_id and ending either
1939
* when the indicated buffer size or the last polyhedra is attained.
1940
* It returns the effective polyhedron index end.
1943
* this_tesselation <-- tesselation structure
1944
* connect_type <-- destination element type
1945
* start_id <-- start index of polyhedra subset in parent section
1946
* buffer_limit <-- maximum number of sub-elements of destination
1947
* element type allowable for vertex_num[] buffer
1948
* extra_vertex_base <-- starting number for added vertices
1949
* vertex_num --> sub-element (global) vertex connectivity
1952
* polyhedron index end corresponding to decoded range
1953
*----------------------------------------------------------------------------*/
1956
_decode_polyhedra_tesselation_l(const fvm_tesselation_t *this_tesselation,
1957
fvm_element_t connect_type,
1958
fvm_lnum_t start_id,
1959
fvm_lnum_t buffer_limit,
1960
fvm_lnum_t extra_vertex_base_num,
1961
fvm_lnum_t vertex_num[])
1964
fvm_lnum_t n_vertices;
1965
fvm_lnum_t i, j, k, l, m, base_dest_id, face_id, vertex_id, encoding_id;
1967
fvm_tesselation_encoding_t decoding_mask[3] = {0, 0, 0};
1969
fvm_lnum_t n_sub_tot = 0;
1970
const fvm_tesselation_t *ts = this_tesselation;
1971
fvm_lnum_t retval = start_id;
1973
_init_decoding_mask(decoding_mask);
1975
/* Main loop on polyhedra */
1976
/*------------------------*/
1978
for (i = 0, j = start_id ;
1979
j < this_tesselation->n_elements;
1982
for (k = ts->face_index[j]; /* Loop on element faces */
1983
k < ts->face_index[j+1];
1986
face_id = FVM_ABS(ts->face_num[k]) - 1;
1988
n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id];
1989
vertex_id = ts->vertex_index[face_id];
1990
encoding_id = ts->vertex_index[face_id] - (face_id*2);
1992
/* Invert base orientation as it points outwards in polyhedron
1993
definition, but inwards in tetrahedron and pyramid reference
1996
if (ts->face_num[k] > 0)
2001
/* Sub-elements (tetrahedra) connectivity */
2002
/*----------------------------------------*/
2004
if ( connect_type == FVM_CELL_TETRA
2005
&& n_vertices > 3 && ts->encoding != NULL) {
2007
/* Loop on face vertices */
2009
for (l = 0; l < (n_vertices - 2); l++) {
2011
if (ts->encoding[encoding_id + l] != 0) {
2013
/* Return previous element's end index if buffer size reached */
2015
if (n_sub_tot >= buffer_limit)
2019
base_dest_id = n_sub_tot*4;
2021
base_dest_id = n_sub_tot*4 + 2;
2023
/* Fill connectivity array */
2025
_DECODE_TRIANGLE_VERTICES(tv,
2026
ts->encoding[encoding_id + l],
2029
for (m = 0; m < 3; m++)
2030
vertex_num[base_dest_id + orient*m]
2031
= ts->vertex_num[vertex_id + tv[m]];
2033
/* Add extra "top" vertex */
2035
vertex_num[n_sub_tot*4 + 3] = extra_vertex_base_num + j;
2045
/* Non-tesselated faces connectivity */
2046
/*-----------------------------------*/
2048
else if ( (connect_type == FVM_CELL_TETRA && n_vertices == 3)
2049
|| (connect_type == FVM_CELL_PYRAM && n_vertices == 4)) {
2051
/* Check that face is not subdivided */
2054
if (ts->encoding != NULL) {
2055
for (l = 0; l < (n_vertices - 2); l++)
2056
if (ts->encoding[encoding_id + l] != 0)
2060
if (l == n_vertices - 2 || ts->encoding == NULL) {
2062
fvm_lnum_t stride = n_vertices + 1;
2064
/* Return previous element's end index if buffer size reached */
2066
if (n_sub_tot >= buffer_limit)
2070
base_dest_id = n_sub_tot*stride;
2072
base_dest_id = n_sub_tot*stride + n_vertices-1;
2074
/* Fill connectivity array */
2076
for (l = 0; l < n_vertices; l++)
2077
vertex_num[base_dest_id + l*orient]
2078
= ts->vertex_num[vertex_id + l];
2080
/* Add extra "top" vertex */
2082
vertex_num[n_sub_tot*stride + n_vertices] = extra_vertex_base_num + j;
2088
} /* End of case where face is not tesselated */
2090
} /* End of loop on element faces */
2092
retval = j+1; /* If end of buffer has not already been reached */
2094
} /* End of main loop on polyhedra */
2099
/*============================================================================
2100
* Public function definitions
2101
*============================================================================*/
2103
/*----------------------------------------------------------------------------
2104
* Creation of a mesh section's subdivision into simpler elements.
2106
* The structure contains pointers to the mesh section's connectivity,
2107
* (passed as arguments), which is not copied. This structure should thus
2108
* always be destroyed before the mesh section to which it relates.
2110
* Unused connectivity array arguments should be set to NULL (such as
2111
* face_index[] and face_num[] for 2D or regular (strided) elements,
2112
* and vertex_index[] for strided elements.
2114
* At this stage, the structure does not yet contain tesselation information.
2117
* element_type <-- type of elements considered
2118
* n_elements <-- number of elements
2119
* face_index <-- polyhedron -> faces index (O to n-1)
2120
* dimension [n_elements + 1]
2121
* face_num <-- element -> face numbers (1 to n, signed,
2122
* > 0 for outwards pointing face normal
2123
* < 0 for inwards pointing face normal);
2124
* dimension: [face_index[n_elements]], or NULL
2125
* vertex_index <-- element face -> vertices index (O to n-1);
2126
* dimension: [n_cell_faces + 1], [n_elements + 1],
2127
* or NULL depending on face_index and vertex_index
2128
* vertex_num <-- element -> vertex connectivity (1 to n)
2129
* global_element_num <-- global element numbers (NULL in serial mode)
2132
* pointer to created mesh section tesselation structure
2133
*----------------------------------------------------------------------------*/
2136
fvm_tesselation_create(fvm_element_t element_type,
2137
fvm_lnum_t n_elements,
2138
const fvm_lnum_t face_index[],
2139
const fvm_lnum_t face_num[],
2140
const fvm_lnum_t vertex_index[],
2141
const fvm_lnum_t vertex_num[],
2142
const fvm_io_num_t *global_element_num)
2146
int entity_dim = 0, stride = 0;
2147
fvm_tesselation_t *this_tesselation;
2149
/* First, check if element type is handled and initialize
2152
switch (element_type) {
2169
/* Now, create structure */
2171
BFT_MALLOC(this_tesselation, 1, fvm_tesselation_t);
2173
/* Assign parent mesh section info */
2175
this_tesselation->type = element_type;
2176
this_tesselation->n_elements = n_elements;
2177
this_tesselation->dim = 0;
2178
this_tesselation->entity_dim = entity_dim;
2180
this_tesselation->stride = stride;
2181
this_tesselation->n_faces = 0;
2183
this_tesselation->vertex_coords = NULL;
2184
this_tesselation->parent_vertex_num = NULL;
2186
this_tesselation->face_index = face_index;
2187
this_tesselation->face_num = face_num;
2188
this_tesselation->vertex_index = vertex_index;
2189
this_tesselation->vertex_num = vertex_num;
2191
this_tesselation->global_element_num = global_element_num;
2193
/* Check argument consistency */
2195
if (face_index != NULL || face_num != NULL) {
2196
if (element_type != FVM_CELL_POLY)
2197
bft_error(__FILE__, __LINE__, 0,
2198
_("Incoherent connectivity for tesselation:\n"
2199
"Connectivity face_index or face_num non NULL,\n"
2200
"but element type != FVM_CELL_POLY"));
2203
else if (vertex_index != NULL) {
2204
if (element_type != FVM_FACE_POLY)
2205
bft_error(__FILE__, __LINE__, 0,
2206
_("Incoherent connectivity for tesselation:\n"
2207
"Connectivy vertex_index non NULL,\n"
2208
"but element type != FVM_FACE_POLY"));
2211
/* Compute number of polyhedron faces */
2213
if (n_elements > 0 && face_index != NULL) {
2214
fvm_lnum_t j, k, face_id;
2215
fvm_lnum_t max_face_id = 0;
2216
for (j = 0; j < n_elements; j++) {
2217
for (k = face_index[j]; k < face_index[j+1]; k++) {
2218
face_id = FVM_ABS(face_num[k]) - 1;
2219
if (face_id > max_face_id)
2220
max_face_id = face_id;
2223
this_tesselation->n_faces = max_face_id + 1;
2226
/* Set tesselation structures to zero */
2228
this_tesselation->n_sub_types = 0;
2230
for (i = 0; i < FVM_TESSELATION_N_SUB_TYPES_MAX; i++) {
2231
this_tesselation->n_sub_max[i] = 0;
2232
this_tesselation->n_sub_max_glob[i] = 0;
2233
this_tesselation->n_sub[i] = 0;
2234
this_tesselation->n_sub_glob[i] =0;
2235
this_tesselation->sub_type[i] = FVM_N_ELEMENT_TYPES;
2238
this_tesselation->encoding = NULL;
2239
this_tesselation->_encoding = NULL;
2241
for (i = 0; i < FVM_TESSELATION_N_SUB_TYPES_MAX; i++) {
2242
this_tesselation->sub_elt_index[i] = NULL;
2243
this_tesselation->_sub_elt_index[i] = NULL;
2246
return (this_tesselation);
2249
/*----------------------------------------------------------------------------
2250
* Destruction of a nodal mesh polygon splitting representation structure.
2253
* this_tesselation <-> pointer to structure that should be destroyed
2257
*----------------------------------------------------------------------------*/
2260
fvm_tesselation_destroy(fvm_tesselation_t * this_tesselation)
2264
if (this_tesselation->_encoding != NULL)
2265
BFT_FREE(this_tesselation->_encoding);
2267
for (i = 0; i < this_tesselation->n_sub_types; i++) {
2268
if (this_tesselation->_sub_elt_index[i] != NULL)
2269
BFT_FREE(this_tesselation->_sub_elt_index[i]);
2271
BFT_FREE(this_tesselation);
2276
/*----------------------------------------------------------------------------
2277
* Tesselate a mesh section referred to by an fvm_tesselation_t structure.
2280
* this_tesselation <-> partially initialized tesselation structure
2281
* dim <-- spatial dimension
2282
* vertex_coords <-- associated vertex coordinates array
2283
* parent_vertex_num <-- optional indirection to vertex coordinates
2284
* error_count --> number of elements with a tesselation error
2285
* counter (optional)
2286
*----------------------------------------------------------------------------*/
2289
fvm_tesselation_init(fvm_tesselation_t *this_tesselation,
2291
const fvm_coord_t vertex_coords[],
2292
const fvm_lnum_t parent_vertex_num[],
2293
fvm_lnum_t *error_count)
2295
assert(this_tesselation != NULL);
2297
this_tesselation->dim = dim;
2299
this_tesselation->vertex_coords = vertex_coords;
2300
this_tesselation->parent_vertex_num = parent_vertex_num;
2302
switch(this_tesselation->type) {
2305
_tesselate_polygons(this_tesselation,
2310
_count_and_index_sub_polyhedra(this_tesselation,
2316
_tesselate_polygons(this_tesselation,
2321
_count_and_index_sub_polygons(this_tesselation,
2326
bft_error(__FILE__, __LINE__, 0,
2327
_("Tesselation of element type %s not implemented."),
2328
fvm_elements_type_name[this_tesselation->type]);
2333
/*----------------------------------------------------------------------------
2334
* Reduction of a nodal mesh polygon splitting representation structure;
2335
* only the associations (numberings) necessary to redistribution of fields
2336
* for output are conserved, the full connectivity being no longer useful
2337
* once it has been output.
2340
* this_tesselation <-> pointer to structure that should be reduced
2341
*----------------------------------------------------------------------------*/
2344
fvm_tesselation_reduce(fvm_tesselation_t * this_tesselation)
2346
this_tesselation->stride = 0;
2347
this_tesselation->n_faces = 0;
2349
if (this_tesselation->face_index == NULL) {
2350
this_tesselation->face_num = NULL;
2351
this_tesselation->vertex_index = NULL;
2352
this_tesselation->vertex_num = NULL;
2355
this_tesselation->encoding = NULL;
2356
if (this_tesselation->_encoding != NULL)
2357
BFT_FREE(this_tesselation->_encoding);
2360
/*----------------------------------------------------------------------------
2361
* Return number of parent elements of a tesselation.
2364
* this_tesselation <-- tesselation structure
2367
* number of parent elements
2368
*----------------------------------------------------------------------------*/
2371
fvm_tesselation_n_elements(const fvm_tesselation_t *this_tesselation)
2373
fvm_lnum_t retval = 0;
2375
if (this_tesselation != NULL)
2376
retval = this_tesselation->n_elements;
2381
/*----------------------------------------------------------------------------
2382
* Return global number of added vertices associated with a tesselation.
2385
* this_tesselation <-- tesselation structure
2388
* global number of added vertices associated with the tesselation
2389
*----------------------------------------------------------------------------*/
2392
fvm_tesselation_n_g_vertices_add(const fvm_tesselation_t *this_tesselation)
2394
fvm_gnum_t retval = 0;
2396
assert(this_tesselation != NULL);
2398
if (this_tesselation->type == FVM_CELL_POLY) {
2400
if (this_tesselation->global_element_num != NULL)
2401
retval = fvm_io_num_get_global_count(this_tesselation->global_element_num);
2403
retval = this_tesselation->n_elements;
2410
/*----------------------------------------------------------------------------
2411
* Return (local) number of added vertices associated with a tesselation.
2414
* this_tesselation <-- tesselation structure
2417
* global number of added vertices associated with the tesselation
2418
*----------------------------------------------------------------------------*/
2421
fvm_tesselation_n_vertices_add(const fvm_tesselation_t *this_tesselation)
2423
fvm_gnum_t retval = 0;
2425
assert(this_tesselation != NULL);
2427
if (this_tesselation->type == FVM_CELL_POLY)
2428
retval = this_tesselation->n_elements;
2433
/*----------------------------------------------------------------------------
2434
* Return number of resulting sub-types of a tesselation.
2437
* this_tesselation <-- tesselation structure
2440
* number of resulting sub-types of the tesselation
2441
*----------------------------------------------------------------------------*/
2444
fvm_tesselation_n_sub_types(const fvm_tesselation_t *this_tesselation)
2448
if (this_tesselation != NULL)
2449
retval = this_tesselation->n_sub_types;
2454
/*----------------------------------------------------------------------------
2455
* Return given sub-types of a tesselation.
2458
* this_tesselation <-- tesselation structure
2459
* sub_type_id <-- index of sub-type in tesselation (0 to n-1)
2462
* sub-types of the tesselation with the given index
2463
*----------------------------------------------------------------------------*/
2466
fvm_tesselation_sub_type(const fvm_tesselation_t *this_tesselation,
2469
fvm_element_t retval = FVM_N_ELEMENT_TYPES;
2471
if (this_tesselation == NULL)
2472
retval = FVM_N_ELEMENT_TYPES;
2474
assert(sub_type_id < this_tesselation->n_sub_types);
2475
retval = this_tesselation->sub_type[sub_type_id];
2481
/*----------------------------------------------------------------------------
2482
* Return number of elements of a given sub-type of a tesselation.
2485
* this_tesselation <-- tesselation structure
2486
* sub_type_id <-- index of sub-type in tesselation (0 to n-1)
2489
* sub-types of the tesselation with the given index
2490
*----------------------------------------------------------------------------*/
2493
fvm_tesselation_n_sub_elements(const fvm_tesselation_t *this_tesselation,
2494
fvm_element_t sub_type)
2498
fvm_lnum_t retval = 0;
2500
if (this_tesselation != NULL) {
2501
for (id = 0; id < this_tesselation->n_sub_types; id++) {
2502
if (this_tesselation->sub_type[id] == sub_type) {
2503
retval = this_tesselation->n_sub[id];
2512
/*----------------------------------------------------------------------------
2513
* Obtain the global and maximum number of elements of a given sub-type
2517
* this_tesselation <-- tesselation structure
2518
* sub_type_id <-- index of sub-type in tesselation (0 to n-1)
2519
* n_sub_elements_glob --> global number of sub-elements of the given type
2520
* n_sub_elements_max --> maximum number of sub-elements per element
2521
* of the given type (for all ranks)
2522
*----------------------------------------------------------------------------*/
2525
fvm_tesselation_get_global_size(const fvm_tesselation_t *this_tesselation,
2526
fvm_element_t sub_type,
2527
fvm_gnum_t *n_sub_elements_glob,
2528
fvm_lnum_t *n_sub_elements_max)
2532
if (n_sub_elements_max != NULL)
2533
*n_sub_elements_max = 0;
2535
if (n_sub_elements_glob != NULL)
2536
*n_sub_elements_glob = 0;
2538
if (this_tesselation != NULL) {
2539
for (id = 0; id < this_tesselation->n_sub_types; id++) {
2540
if (this_tesselation->sub_type[id] == sub_type) {
2541
if (n_sub_elements_max != NULL)
2542
*n_sub_elements_max = this_tesselation->n_sub_max_glob[id];
2543
if (n_sub_elements_glob != NULL)
2544
*n_sub_elements_glob = this_tesselation->n_sub_glob[id];
2551
/*----------------------------------------------------------------------------
2552
* Return global numbering of added vertices associated with a tesselation.
2555
* this_tesselation <-- tesselation structure
2558
* pointer to global numbering of added vertices for this tesselation,
2559
* or NULL if no added vertices are present.
2560
*----------------------------------------------------------------------------*/
2562
const fvm_io_num_t *
2563
fvm_tesselation_global_vertex_num(const fvm_tesselation_t *this_tesselation)
2565
const fvm_io_num_t *retval = NULL;
2567
assert(this_tesselation != NULL);
2569
if (this_tesselation->type == FVM_CELL_POLY)
2570
retval = this_tesselation->global_element_num;
2575
/*----------------------------------------------------------------------------
2576
* Compute coordinates of added vertices for a tesselation of polyhedra.
2578
* One additional vertex is added near the center of each polyhedra.
2579
* For element types other than polyhedra, there is no need for added
2580
* vertices, so this function returns immediately.
2583
* this_tesselation <-- tesselation structure
2584
* vertex_coords --> coordinates of added vertices
2585
*----------------------------------------------------------------------------*/
2588
fvm_tesselation_vertex_coords(const fvm_tesselation_t *this_tesselation,
2589
fvm_coord_t vertex_coords[])
2593
if (this_tesselation->type != FVM_CELL_POLY)
2596
for (i = 0; i < this_tesselation->n_elements; i++) {
2598
_added_vertex_coords(this_tesselation, vertex_coords + i*3, NULL, i);
2604
/*----------------------------------------------------------------------------
2605
* Return index of sub-elements associated with each element of a given
2606
* sub-type of a tesselation.
2609
* this_tesselation <-- tesselation structure
2610
* sub_type_id <-- index of sub-type in tesselation (0 to n-1)
2613
* index of sub-elements associated with each element (0 to n-1 numbering)
2614
*----------------------------------------------------------------------------*/
2617
fvm_tesselation_sub_elt_index(const fvm_tesselation_t *this_tesselation,
2618
fvm_element_t sub_type)
2621
const fvm_lnum_t *retval = NULL;
2623
if (this_tesselation != NULL) {
2624
for (id = 0; id < this_tesselation->n_sub_types; id++) {
2625
if (this_tesselation->sub_type[id] == sub_type) {
2626
retval = this_tesselation->sub_elt_index[id];
2635
#if defined(HAVE_MPI)
2637
/*----------------------------------------------------------------------------
2638
* Compute index values corresponding to given range of indices,
2639
* for an element -> sub-element value distribution.
2641
* This index is used mainly to gather a decoded tesselation connectivity
2642
* or element -> sub-element data distribution, expanding the corresponding
2643
* data only on the given range.
2644
* Only the index values in the start_id to end_id range are set by
2645
* this function, starting with index[start_id] = 0.
2648
* this_tesselation <-- tesselation structure
2649
* connect_type <-- destination element type
2650
* stride <-- number of associated values per sub-element
2651
* start_id <-- start index of polyhedra subset in parent section
2652
* buffer_limit <-- maximum number of sub-elements of destination
2653
* element type allowable for vertex_num[] buffer
2654
* global_num_end <-> past the end (maximum + 1) parent element
2655
* global number (reduced on return if required
2656
* by buffer_size limits)
2657
* index --> sub-element index
2658
* comm <-- associated MPI communicator
2661
* polyhedron index end corresponding to decoded range
2662
*----------------------------------------------------------------------------*/
2665
fvm_tesselation_range_index_g(const fvm_tesselation_t *this_tesselation,
2666
fvm_element_t connect_type,
2668
fvm_lnum_t start_id,
2669
fvm_lnum_t buffer_limit,
2670
fvm_gnum_t *global_num_end,
2675
fvm_lnum_t buffer_size = buffer_limit * stride;
2677
const fvm_gnum_t *global_element_num
2678
= fvm_io_num_get_global_num(this_tesselation->global_element_num);
2680
const fvm_lnum_t *sub_element_idx
2681
= fvm_tesselation_sub_elt_index(this_tesselation, connect_type);
2683
/* In parallel mode, global_element_num should exist */
2688
/* Loop on tesselated elements */
2689
/*-----------------------------*/
2691
index[start_id] = 0;
2693
for (i = start_id; i < this_tesselation->n_elements; i++) {
2695
if (global_element_num[i] >= *global_num_end)
2698
index[i+1] = index[i] + ( sub_element_idx[i+1]
2699
- sub_element_idx[i]) * stride;
2701
if (index[i+1] > buffer_size) {
2702
*global_num_end = global_element_num[i];
2708
/* Check if the maximum id returned on some ranks leads to
2709
a lower global_num_end than initially required
2710
(due to local buffer being full) */
2712
_expand_limit_g(this_tesselation,
2720
/*----------------------------------------------------------------------------
2721
* Decode tesselation to a parent element->sub elements index and
2722
* connectivity buffer.
2724
* To avoid requiring huge buffers and computing unneeded element
2725
* connectivities when exporting data in slices, this function may decode
2726
* a partial connectivity range, starting at polygon index start_id and ending
2727
* either when the indicated buffer size is attained, or the global element
2728
* number corresponding to a given polygon exceeds a given value.
2729
* It returns the effective polygon index end.
2732
* this_tesselation <-- tesselation structure
2733
* connect_type <-- destination element type
2734
* start_id <-- start index of polygons subset in parent section
2735
* buffer_limit <-- maximum number of sub-elements of destination
2736
* element type allowable for sub_element_idx[]
2737
* and vertex_num[] buffers
2738
* global_num_end <-> past the end (maximum + 1) parent element
2739
* global number (reduced on return if required
2741
* extra_vertex_base <-- starting number for added vertices
2742
* global_vertex_num <-- global vertex numbering
2743
* vertex_num --> sub-element (global) vertex connectivity
2744
* comm <-- associated MPI communicator
2747
* polygon index corresponding to end of decoded range
2748
*----------------------------------------------------------------------------*/
2751
fvm_tesselation_decode_g(const fvm_tesselation_t *this_tesselation,
2752
fvm_element_t connect_type,
2753
fvm_lnum_t start_id,
2754
fvm_lnum_t buffer_limit,
2755
fvm_gnum_t *global_num_end,
2756
const fvm_io_num_t *global_vertex_num,
2757
fvm_gnum_t extra_vertex_base,
2758
fvm_gnum_t vertex_num[],
2761
fvm_lnum_t retval = 0;
2763
switch(this_tesselation->type) {
2766
retval = _decode_polyhedra_tesselation_g(this_tesselation,
2777
retval = _decode_polygons_tesselation_g(this_tesselation,
2791
/* Check if the maximum id returned on some ranks leads to
2792
a lower global_num_end than initially required
2793
(due to local buffer being full) */
2795
_expand_limit_g(this_tesselation,
2803
#endif /* defined(HAVE_MPI) */
2805
/*----------------------------------------------------------------------------
2806
* Decode tesselation to a connectivity buffer.
2808
* To avoid requiring huge buffers and computing unneeded element
2809
* connectivities, this function may decode a partial connectivity range,
2810
* starting at polygon index start_id and ending either when the indicated
2811
* buffer size or the last polygon is attained.
2812
* It returns the effective polygon index end.
2815
* this_tesselation <-- tesselation structure
2816
* connect_type <-- destination element type
2817
* start_id <-- start index of polygons subset in parent section
2818
* buffer_limit <-- maximum number of sub-elements of destination
2819
* element type allowable for sub_element_idx[]
2820
* and vertex_num[] buffers
2821
* extra_vertex_base <-- starting number for added vertices
2822
* vertex_num --> sub-element (global) vertex connectivity
2825
* polygon index corresponding to end of decoded range
2826
*----------------------------------------------------------------------------*/
2829
fvm_tesselation_decode(const fvm_tesselation_t *this_tesselation,
2830
fvm_element_t connect_type,
2831
fvm_lnum_t start_id,
2832
fvm_lnum_t buffer_limit,
2833
fvm_lnum_t extra_vertex_base,
2834
fvm_lnum_t vertex_num[])
2836
fvm_lnum_t retval = 0;
2838
switch(this_tesselation->type) {
2841
retval = _decode_polyhedra_tesselation_l(this_tesselation,
2850
retval = _decode_polygons_tesselation_l(this_tesselation,
2865
/*----------------------------------------------------------------------------
2866
* Distribute "per element" data from the base elements to their tesselation.
2868
* The same data array is used for input and output, so as to avoid requiring
2869
* excess allocation in typical use cases (extracting data from a parent mesh
2870
* to a buffer and distributing it as per its tesselation).
2871
* The data array should be at least of size:
2872
* [sub_elt_index[end_id] - sub_elt_index[start_id] * size
2875
* this_tesselation <-- tesselation structure
2876
* connect_type <-- destination element type
2877
* start_id <-- start index of elements subset in parent section
2878
* end_id <-- end index of elements subset in parent section
2879
* size <-- data size for each element (sizeof(type)*stride)
2880
* data <-> undistributed data in, distributed data out
2881
*----------------------------------------------------------------------------*/
2884
fvm_tesselation_distribute(const fvm_tesselation_t *this_tesselation,
2885
fvm_element_t connect_type,
2886
fvm_lnum_t start_id,
2892
fvm_lnum_t i, j, k, n_sub;
2896
const fvm_lnum_t *sub_elt_index = NULL;
2898
/* Find index, or return */
2900
if (this_tesselation == NULL)
2903
for (id = 0; id < this_tesselation->n_sub_types; id++) {
2904
if (this_tesselation->sub_type[id] == connect_type) {
2905
sub_elt_index = this_tesselation->sub_elt_index[id];
2909
if (id == this_tesselation->n_sub_types)
2912
/* Distribute data (starting from the end so as not to overwrite
2913
data at the beginning of the array) */
2915
for (i = end_id, j = end_id - start_id - 1; i > start_id; i--, j--) {
2917
src = ((char *)data) + j*size;
2918
dest = ((char *)data) + (sub_elt_index[i-1] - sub_elt_index[start_id])*size;
2919
n_sub = sub_elt_index[i] - sub_elt_index[i-1];
2921
for (k = 0; k < n_sub; k++) {
2922
for (l = 0; l < size; l++)
2923
dest[k*size + l] = src[l];
2929
/*----------------------------------------------------------------------------
2930
* Compute field values at added vertices for a tesselation of polyhedra.
2932
* One additional vertex is added near the center of each polyhedra.
2933
* For element types other than polyhedra, there is no need for added
2934
* vertices, so this function returns immediately.
2937
* this_tesselation <-- tesselation structure
2938
* vertex_coords <-- coordinates of added vertices
2939
* src_dim <-- dimension of source data
2940
* src_dim_shift <-- source data dimension shift (start index)
2941
* dest_dim <-- destination data dimension (1 if non interlaced)
2942
* start_id <-- added vertices start index
2943
* end_id <-- added vertices past the end index
2944
* src_interlace <-- indicates if source data is interlaced
2945
* src_datatype <-- source data type (float, double, or int)
2946
* dest_datatype <-- destination data type (float, double, or int)
2947
* n_parent_lists <-- number of parent lists (if parent_num != NULL)
2948
* parent_num_shift <-- parent number to value array index shifts;
2949
* size: n_parent_lists
2950
* parent_num <-- if n_parent_lists > 0, parent entity numbers
2951
* src_data <-- array of source arrays (at least one, with one per
2952
* source dimension if non interlaced, times one per
2953
* parent list if multiple parent lists, with
2954
* x_parent_1, y_parent_1, ..., x_parent_2, ...) order
2955
* dest_data --> destination buffer
2956
*----------------------------------------------------------------------------*/
2959
fvm_tesselation_vertex_values(const fvm_tesselation_t *this_tesselation,
2963
fvm_lnum_t start_id,
2965
fvm_interlace_t src_interlace,
2966
fvm_datatype_t src_datatype,
2967
fvm_datatype_t dest_datatype,
2969
const fvm_lnum_t parent_num_shift[],
2970
const fvm_lnum_t parent_num[],
2971
const void *const src_data[],
2972
void *const dest_data)
2974
/* If source or destination datatype is not floating-point,
2975
set all return values to zero */
2977
if ( (src_datatype != FVM_DOUBLE && src_datatype != FVM_FLOAT)
2978
|| (dest_datatype != FVM_DOUBLE && dest_datatype != FVM_FLOAT)) {
2980
unsigned char *_dest_data = dest_data;
2982
size_t data_shift = start_id
2983
* (dest_dim * fvm_datatype_size[dest_datatype]);
2984
size_t data_size_c = (end_id - start_id)
2985
* (dest_dim * fvm_datatype_size[dest_datatype]);
2987
memset(_dest_data + data_shift, 0, data_size_c);
2991
/* Otherwise, interpolate values */
2995
_vertex_field_of_real_values(this_tesselation,
3013
/*----------------------------------------------------------------------------
3014
* Dump printout of a mesh section tesselation structure.
3017
* this_tesselation <-- pointer to structure that should be dumped
3018
*----------------------------------------------------------------------------*/
3021
fvm_tesselation_dump(const fvm_tesselation_t *this_tesselation)
3024
fvm_lnum_t n_elements, j, k;
3025
const fvm_lnum_t *idx;
3027
if (this_tesselation == NULL)
3030
/* Global indicators */
3031
/*--------------------*/
3035
"Element type: %s\n"
3036
"Number of elements: %ld\n"
3037
"Spatial dimension: %d\n"
3038
"Entity dimension: %d\n",
3039
fvm_elements_type_name[this_tesselation->type],
3040
(long)this_tesselation->n_elements,
3041
this_tesselation->dim, this_tesselation->entity_dim);
3045
"Number of faces: %d\n",
3046
this_tesselation->stride,
3047
(long)(this_tesselation->n_faces));
3050
"Pointers to shared arrays:\n"
3051
" vertex_coords %p\n"
3052
" parent_vertex_num %p\n"
3054
" vertex_index: %p\n"
3055
" vertex_num: %p\n",
3056
this_tesselation->vertex_coords,
3057
this_tesselation->parent_vertex_num,
3058
this_tesselation->face_index, this_tesselation->face_num,
3059
this_tesselation->vertex_index, this_tesselation->vertex_num);
3062
"Pointers to shared global numbering:\n"
3063
" global_element_num %p\n",
3064
this_tesselation->global_element_num);
3067
/* Basic information */
3068
/*-------------------*/
3071
"Number of sub-entity types: %d\n\n",
3072
this_tesselation->n_sub_types);
3074
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3075
bft_printf("Maximum local number of resulting %s per element: %ld\n",
3076
fvm_elements_type_name[this_tesselation->sub_type[i]],
3077
(long)(this_tesselation->n_sub_max[i]));
3080
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3081
bft_printf("Maximum global number of resulting %s per element: %ld\n",
3082
fvm_elements_type_name[this_tesselation->sub_type[i]],
3083
(long)(this_tesselation->n_sub_max_glob[i]));
3088
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3089
bft_printf("Local number of resulting %s: %ld\n",
3090
fvm_elements_type_name[this_tesselation->sub_type[i]],
3091
(long)(this_tesselation->n_sub[i]));
3094
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3095
bft_printf("Global number of resulting %s: %llu\n",
3096
fvm_elements_type_name[this_tesselation->sub_type[i]],
3097
(unsigned long long)(this_tesselation->n_sub_glob[i]));
3101
"Pointers to shareable arrays:\n"
3103
this_tesselation->encoding);
3105
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3106
if (this_tesselation->sub_elt_index[i] != NULL)
3107
bft_printf(" sub_elt_index[%s]: %p\n",
3108
fvm_elements_type_name[this_tesselation->sub_type[i]],
3109
this_tesselation->sub_elt_index[i]);
3113
"Pointers to local arrays:\n"
3115
this_tesselation->_encoding);
3117
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3118
if (this_tesselation->sub_elt_index[i] != NULL)
3119
bft_printf(" _sub_elt_index[%s]: %p\n",
3120
fvm_elements_type_name[this_tesselation->sub_type[i]],
3121
this_tesselation->_sub_elt_index[i]);
3124
if (this_tesselation->encoding != NULL) {
3126
fvm_tesselation_encoding_t decoding_mask[3] = {0, 0, 0};
3129
_init_decoding_mask(decoding_mask);
3131
if (this_tesselation->type != FVM_FACE_QUAD) {
3132
bft_printf("\nEncoding (local vertex numbers):\n\n");
3133
if (this_tesselation->n_faces > 0)
3134
n_elements = this_tesselation->n_faces;
3136
n_elements = this_tesselation->n_elements;
3137
idx = this_tesselation->vertex_index;
3138
for (j = 0; j < n_elements; j++) {
3139
_DECODE_TRIANGLE_VERTICES(tv,
3140
this_tesselation->encoding[idx[j] - 2*j],
3142
bft_printf("%10d (idx = %10d) %10d %10d %10d\n",
3143
j+1, idx[j], (int)tv[0], (int)tv[1], (int)tv[2]);
3144
for (k = idx[j] -2*j + 1; k < idx[j+1] - 2*j; k++) {
3145
_DECODE_TRIANGLE_VERTICES(tv,
3146
this_tesselation->encoding[k],
3148
bft_printf(" %10d %10d %10d\n",
3149
(int)tv[0], (int)tv[1], (int)tv[2]);
3152
bft_printf(" end (idx = %10d)\n", idx[n_elements]);
3154
else { /* if (this_tesselation->type != FVM_FACE_QUAD) */
3155
bft_printf("\nEncoding (diagonal flag):\n\n");
3156
n_elements = this_tesselation->n_elements;
3157
for (j = 0; j < n_elements; j++)
3158
bft_printf("%10d: %10d\n", j+1, (int)this_tesselation->encoding[j]);
3163
for (i = 0; i < this_tesselation->n_sub_types; i++) {
3164
if (this_tesselation->sub_elt_index[i] != NULL) {
3165
bft_printf("\nSub-element index [%s]:\n\n",
3166
fvm_elements_type_name[this_tesselation->sub_type[i]]);
3167
n_elements = this_tesselation->n_elements;
3168
idx = this_tesselation->sub_elt_index[i];
3169
for (j = 0; j < n_elements; j++)
3170
bft_printf("%10d: idx = %10d\n", j+1, idx[j]);
3171
bft_printf(" end: idx = %10d\n", idx[n_elements]);
3177
/*----------------------------------------------------------------------------*/
3181
#endif /* __cplusplus */
added
removed
src/fvm/fvm_tesselation.c
