/*============================================================================ * Structure describing a mesh section's subdivision into simpler elements * * This is mostly useful to replace polygons or polyhedra by simpler * elements such as triangles, tetrahedra, and prisms upon data export. *============================================================================*/ /* This file is part of the "Finite Volume Mesh" library, intended to provide finite volume mesh and associated fields I/O and manipulation services. Copyright (C) 2006 EDF This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */ /*----------------------------------------------------------------------------*/ /* * Notes: * * For a polygon with n vertices (and thus n edges), a compact way of * storing a tesselation is described here: * For a given vertex index i, a triangle is defined by * {vertex_num[i], vertex_num[i+1], opposite_vertex_num[i]} * when all values are > 0. When opposite_vertex_num[i] = 0, * no triangle is defined (indicating that the corresponding triangle * has already been defined, as we have (n - 2) triangles per polygon * (or less in case of tesselation errors), or that the polygon is * not tesselated (i.e. a quadrangle, if allowed). * This requires n values per polygon, while the full triangle * connectivity requires (n - 2) * 3 values. * * As polyhedra are defined by their polygonal faces, this encoding * is also used. * * For a quadrangle, a better way is to use a single diagonal_flag, equal to * O if the quadrangle is split along its diagonal (1, 3), or * 1 if the quadrangle is split along its diagonal (2, 4). * * In the fvm_tesselation_t structure, the encoding[] array corresponds * to the opposite_vertex_num[] array (of the same size as vertex_num[]) * for polygons and polyhedra, and to a diagonal_flag[] array of size * n_elements for quadrangles. * * If FVM model were ever extended to elements with more than 2 vertices * per edge, we would need to replace n_vertices by n_edges in some of * this file's function's logic, but the biggest difficulty would * be in adding extra vertices along new edges without duplicating them: * this would only require some extra processing for a given section, * (including some communication between processors), but would * involve difficulties mainly when using multiple tesselation structures * for multiple mesh sections, as extra vertices along section boundaries * would need to be shared, involving extra book-keeping (and not sharing * those vertices would make the corresponding elements non-conforming, * so it would not be an option). As in most known data models, general * simple polygons or polyhedra only have two vertices per edge, * restricting tesselation to linear element types would not be a problem * as regards its main intended use, which is to allow replacement of * polygons and/or polyhedra upon export to some formats or for use with * some tools which do not support these types of elements. */ /*---------------------------------------------------------------------------- * Standard C library headers *----------------------------------------------------------------------------*/ #include #include #include #include #include /*---------------------------------------------------------------------------- * BFT library headers *----------------------------------------------------------------------------*/ #include #include #include /*---------------------------------------------------------------------------- * Local headers *----------------------------------------------------------------------------*/ #include #include #include #include #include /*---------------------------------------------------------------------------- * Header for the current file *----------------------------------------------------------------------------*/ #include /*----------------------------------------------------------------------------*/ #ifdef __cplusplus extern "C" { #if 0 } /* Fake brace to force Emacs auto-indentation back to column 0 */ #endif #endif /* __cplusplus */ /*============================================================================= * Macro definitions *============================================================================*/ #undef _CROSS_PRODUCT_3D #undef _DOT_PRODUCT_3D #undef _MODULE_3D #define _CROSS_PRODUCT_3D(cross_v1_v2, v1, v2) ( \ cross_v1_v2[0] = v1[1]*v2[2] - v1[2]*v2[1], \ cross_v1_v2[1] = v1[2]*v2[0] - v1[0]*v2[2], \ cross_v1_v2[2] = v1[0]*v2[1] - v1[1]*v2[0] ) #define _DOT_PRODUCT_3D(v1, v2) ( \ v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]) #define _MODULE_3D(v) \ sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]) /*============================================================================ * Type definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Structure defining a tesselation of a mesh section. *----------------------------------------------------------------------------*/ struct _fvm_tesselation_t { /* Parent section information */ /*----------------------------*/ fvm_element_t type; /* Element types */ fvm_lnum_t n_elements; /* Number of elements */ int dim; /* Spatial dimension */ int entity_dim; /* Entity dimension */ int stride; /* Element size for regular elements (0 for polygons and polyhedra) */ fvm_lnum_t n_faces; /* Number of faces defining polyhedra */ /* Pointers to shared vertex coordinates */ const fvm_coord_t *vertex_coords; /* pointer to vertex coordinates (always interlaced: x1, y1, z1, x2, y2, z2, ...) */ const fvm_lnum_t *parent_vertex_num; /* Local numbers (1 to n) of local vertices in the parent mesh. This array is present only when non "trivial" (i.e. not 1, 2, ..., n). */ /* Pointers to shared connectivity arrays */ const fvm_lnum_t *face_index; /* polyhedron -> faces index (O to n-1); dimension [n_elements + 1] */ const fvm_lnum_t *face_num; /* polyhedron -> face numbers (1 to n, signed, > 0 for outwards pointing face normal < 0 for inwards pointing face normal); dimension: [face_index[n_elements]] */ const fvm_lnum_t *vertex_index; /* polygon face -> vertices index (O to n-1); dimension: [n_cell_faces + 1] */ const fvm_lnum_t *vertex_num; /* vertex numbers (1 to n); dimension: connectivity_size */ /* Pointer to shared global element numbers */ const fvm_io_num_t *global_element_num; /* Basic information */ /*-------------------*/ int n_sub_types; /* Number of sub-element types (currently max 2) */ fvm_element_t sub_type[2]; /* Sub-element types */ fvm_lnum_t n_sub_max[2]; /* Maximum number of sub-elements of each type per element */ fvm_lnum_t n_sub_max_glob[2]; /* Global maximum number of sub-elements of each type per element */ fvm_lnum_t n_sub[2]; /* Number of sub-elements of each type */ fvm_gnum_t n_sub_glob[2]; /* Global number of sub-elements of each type */ const fvm_lnum_t *encoding; /* Compact tesselation encoding array (see notes above) */ fvm_lnum_t *_encoding; /* encoding if owner, NULL if shared */ const fvm_lnum_t *sub_elt_index[2]; /* index of sub-elements of each given type associated with each element (0 to n-1); */ fvm_lnum_t *_sub_elt_index[2]; /* sub_elt_index if owner, NULL if shared */ }; /*============================================================================ * Static global variables *============================================================================*/ /*============================================================================ * Private function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Compute coordinates of added vertices for a tesselation of a polyhedron. * * parameters: * ts <-- tesselation structure * vertex_coords --> coordinates of added vertices * n_vertices_tot --> total number of vertex references (or NULL) * element_id <-- element index *----------------------------------------------------------------------------*/ static inline void _added_vertex_coords(const fvm_tesselation_t *ts, fvm_coord_t vertex_coords[3], int *n_vertices_tot, fvm_lnum_t element_id) { fvm_lnum_t n_vertices; fvm_lnum_t i, j, k, face_id, vertex_id; fvm_lnum_t _n_vertices_tot = 0; fvm_coord_t cell_center[3] = {0., 0., 0.}; fvm_coord_t cell_center_divisor = 0.; /* Loop on polyhedron's faces */ /*----------------------------*/ for (i = ts->face_index[element_id]; i < ts->face_index[element_id + 1]; i++) { fvm_coord_t sign, v1[3], v2[3]; fvm_coord_t vertices_center[3] = {0., 0., 0.}; fvm_coord_t face_center[3] = {0., 0., 0.}; fvm_coord_t face_normal[3] = {0., 0., 0.}; fvm_coord_t triangle_center[3] = {0., 0., 0.}; fvm_coord_t triangle_normal[3] = {0., 0., 0.}; fvm_coord_t face_surface = 0.; fvm_coord_t triangle_surface = 0.; const fvm_coord_t *current_coords = NULL; face_id = FVM_ABS(ts->face_num[i]) - 1; n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id]; _n_vertices_tot += n_vertices; /* First loop on face vertices to estimate face center location */ for (j = 0; j < n_vertices; j++) { vertex_id = ts->vertex_num[ts->vertex_index[face_id] + j] - 1; if (ts->parent_vertex_num != NULL) current_coords = ts->vertex_coords + ((ts->parent_vertex_num[vertex_id] - 1) * 3); else current_coords = ts->vertex_coords + (vertex_id * 3); for (k = 0; k < 3; k++) vertices_center[k] += current_coords[k]; } for (k = 0; k < 3; k++) vertices_center[k] /= (fvm_coord_t)n_vertices; /* At this stage, current_coords points to the last face vertice's coords */ for (k = 0; k < 3; k++) { v1[k] = current_coords[k] - vertices_center[k]; triangle_center[k] = current_coords[k] + vertices_center[k]; } /* Second loop on face vertices to estimate face normal */ for (j = 0; j < n_vertices; j++) { vertex_id = ts->vertex_num[ts->vertex_index[face_id] + j] - 1; if (ts->parent_vertex_num != NULL) current_coords = ts->vertex_coords + ((ts->parent_vertex_num[vertex_id] - 1) * 3); else current_coords = ts->vertex_coords + (vertex_id * 3); for (k = 0; k < 3; k++) { v2[k] = current_coords[k] - vertices_center[k]; triangle_center[k] = (triangle_center[k] + current_coords[k]) * (1./3.); } _CROSS_PRODUCT_3D(triangle_normal, v1, v2); for (k = 0; k < 3; k++) face_normal[k] += triangle_normal[k]; triangle_surface = _MODULE_3D(triangle_normal) * 0.5; if (_DOT_PRODUCT_3D(triangle_normal, face_normal) > 0) sign = 1.; else sign = -1.; /* Add triangle contributions to face center and normal */ face_surface += triangle_surface * sign; for (k = 0; k < 3; k++) face_center[k] += triangle_center[k] * triangle_surface * sign; /* Prepare for next face vertex */ for (k = 0; k < 3; k++) { v1[k] = v2[k]; triangle_center[k] = current_coords[k] + vertices_center[k]; } } /* End of second loop on face vertices */ /* If first normal was not oriented like the face normal, correct sign */ if (face_surface < 0) { face_surface = - face_surface; for (k = 0; k < 3; k++) face_center[k] = -face_center[k]; } /* We should now divide the face_center[] values by face_surface to finish, but the contribution of face center to element center is: face_center[] * face_surface; which the face_center[] array contains, so we use those directly */ cell_center_divisor += face_surface; for (k = 0; k < 3; k++) cell_center[k] += face_center[k]; } /* End of loop on polyhedron's faces */ for (k = 0; k < 3; k++) vertex_coords[k] = cell_center[k] / cell_center_divisor; if (n_vertices_tot != NULL) *n_vertices_tot = _n_vertices_tot; } /*---------------------------------------------------------------------------- * Solve Ax = B for a 4x4 system. * * parameters: * a <-- matrix * b <-- right hand side * x --> solution of Ax = b * * returns: 0 if solution was found, 1 in case of error (zero pivot) *----------------------------------------------------------------------------*/ static int _solve_ax_b_4(double a[4][4], double *restrict b, double *restrict x) { int i, j, k, k_pivot; double abs_pivot, abs_a_ki, factor; double _a[4][4], _b[4]; double _a_swap[4], _b_swap; const double _epsilon = 1.0e-15; /* Copy array and RHS (avoid modifying a and b) */ for (i = 0; i < 4; i++) { for (j = 0; j < 4; j++) { _a[i][j] = a[i][j]; } _b[i] = b[i]; } /* Forward elimination */ for (i = 0; i < 4-1; i++) { /* Seek best pivot */ k_pivot = i; abs_pivot = FVM_ABS(_a[i][i]); for (k = i+1; k < 4; k++) { abs_a_ki = FVM_ABS(_a[k][i]); if (abs_a_ki > abs_pivot) { abs_pivot = abs_a_ki; k_pivot = k; } } /* Swap if necessary */ if (k_pivot != i) { for (j = 0; j < 4; j++) { _a_swap[j] = _a[i][j]; _a[i][j] = _a[k_pivot][j]; _a[k_pivot][j] = _a_swap[j]; } _b_swap = _b[i]; _b[i] = _b[k_pivot]; _b[k_pivot] = _b_swap; } if (abs_pivot < _epsilon) return 1; /* Eliminate values */ for (k = i+1; k < 4; k++) { factor = _a[k][i] / _a[i][i]; _a[k][i] = 0.0; for (j = i+1; j < 4; j++) { _a[k][j] -= _a[i][j]*factor; } _b[k] -= _b[i]*factor; } } /* Solve system */ x[3] = _b[3] / _a[3][3]; x[2] = (_b[2] - _a[2][3]*x[3]) / _a[2][2]; x[1] = (_b[1] - _a[1][3]*x[3] - _a[1][2]*x[2]) / _a[1][1]; x[0] = (_b[0] - _a[0][3]*x[3] - _a[0][2]*x[2] - _a[0][1]*x[1]) / _a[0][0]; return 0; } /*---------------------------------------------------------------------------- * Compute coordinates of added vertices for a tesselation of a polyhedron. * * The heap and vertex_list arrays should be large enough to contain * all vertex references of all the polyhedron's faces. * * parameters: * ts <-- tesselation structure * heap --- working array * vertex_list --> list of polyhedron's vertices * vertex_list_size --> list size * element_id <-- element index *----------------------------------------------------------------------------*/ static inline void _polyhedron_vertices(const fvm_tesselation_t *ts, fvm_lnum_t heap[], fvm_lnum_t vertex_list[], int *vertex_list_size, fvm_lnum_t element_id) { fvm_lnum_t n_vertices; fvm_lnum_t i, j, face_id, vertex_id; int heap_size = 0; /* Loop on polyhedron's faces */ /*----------------------------*/ for (i = ts->face_index[element_id]; i < ts->face_index[element_id + 1]; i++) { /* Loop on face vertices */ /*-----------------------*/ face_id = FVM_ABS(ts->face_num[i]) - 1; n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id]; for (j = 0; j < n_vertices; j++) { /* Add to heap */ int hole = heap_size; int parent = (hole - 1) / 2; vertex_id = ts->vertex_num[ts->vertex_index[face_id] + j] - 1; while (hole > 0 && heap[parent] > vertex_id) { heap[hole] = heap[parent]; hole = parent; parent = (parent - 1) / 2; } heap[hole] = vertex_id; heap_size++; } } /* Now remove entries from heap */ i = 0; if (heap_size > 0) vertex_list[i++] = heap[0]; while (heap_size > 0) { if (heap[0] != vertex_list[i-1]) vertex_list[i++] = heap[0]; heap_size--; /* Insert the last item in the heap whose root is empty */ { int start = 0; int child = 1; /* left child of 0 */ vertex_id = heap[heap_size]; while (child < heap_size) { if (child < (heap_size - 1) && heap[child] > heap[child+1]) child++; /* child is the smaller child of start */ if (vertex_id < heap[child]) break; /* item stays in start */ else { heap[start] = heap[child]; start = child; child = 2*start + 1; } } heap[start] = vertex_id; } } *vertex_list_size = i; } /*---------------------------------------------------------------------------- * Compute field values at added vertices for a tesselation of polyhedra. * * One additional vertex is added near the center of each polyhedra. * For element types other than polyhedra, there is no need for added * vertices, so this function returns immediately. * * parameters: * this_tesselation <-- tesselation structure * vertex_coords <-- coordinates of added vertices * src_dim <-- dimension of source data * src_dim_shift <-- source data dimension shift (start index) * dest_dim <-- destination data dimension (1 if non interlaced) * start_id <-- added vertices start index * end_id <-- added vertices past the end index * src_interlace <-- indicates if source data is interlaced * src_datatype <-- source data type (float, double, or int) * dest_datatype <-- destination data type (float, double, or int) * n_parent_lists <-- number of parent lists (if parent_num != NULL) * parent_num_shift <-- parent number to value array index shifts; * size: n_parent_lists * parent_num <-- if n_parent_lists > 0, parent entity numbers * src_data <-- array of source arrays (at least one, with one per * source dimension if non interlaced, times one per * parent list if multiple parent lists, with * x_parent_1, y_parent_1, ..., x_parent_2, ...) order * dest_data --> destination buffer *----------------------------------------------------------------------------*/ static void _vertex_field_of_real_values(const fvm_tesselation_t *this_tesselation, int src_dim, int src_dim_shift, int dest_dim, fvm_lnum_t start_id, fvm_lnum_t end_id, fvm_interlace_t src_interlace, fvm_datatype_t src_datatype, fvm_datatype_t dest_datatype, int n_parent_lists, const fvm_lnum_t parent_num_shift[], const fvm_lnum_t parent_num[], const void *const src_data[], void *const dest_data) { int n_vertices_tot, pl, vertex_list_size; fvm_lnum_t i, j, parent_id, src_id, vertex_id; int _src_dim_shift = src_dim_shift; int max_list_size = 128; fvm_lnum_t _heap[128]; fvm_lnum_t _vertex_list[128]; fvm_lnum_t *heap = _heap; fvm_lnum_t *vertex_list = _vertex_list; const fvm_tesselation_t *ts = this_tesselation; const fvm_coord_t *current_coords = NULL; if (ts->type != FVM_CELL_POLY) return; /* Main loop on polyhedra */ /*------------------------*/ for (i = start_id; i < end_id; i++) { for (j = 0; j < dest_dim; j++) { /* Loop on destination dimension */ fvm_coord_t vertex_coords[3]; double coeff[4]; double interpolated_value, v_x, v_y, v_z; double v_f = 0.; double a[4][4] = {{0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}, {0., 0., 0., 0.}}; double b[4] = {0., 0., 0., 0.}; _added_vertex_coords(ts, vertex_coords, &n_vertices_tot, i); /* Allocate or resize working arrays if too small */ if (n_vertices_tot > max_list_size) { max_list_size *= 2; if (heap == _heap) { heap = NULL; vertex_list = NULL; } BFT_REALLOC(heap, max_list_size, fvm_lnum_t); BFT_REALLOC(vertex_list, max_list_size, fvm_lnum_t); } /* Obtain list of polyhedron's vertices */ _polyhedron_vertices(ts, heap, vertex_list, &vertex_list_size, i); /* Build matrix for least squares */ for (j = 0; j < vertex_list_size; j++) { vertex_id = vertex_list[j]; if (ts->parent_vertex_num != NULL) current_coords = ts->vertex_coords + ((ts->parent_vertex_num[vertex_id] - 1) * 3); else current_coords = ts->vertex_coords + (vertex_id * 3); v_x = current_coords[0]; v_y = current_coords[1]; v_z = current_coords[2]; /* Position in source data for current vertex value */ if (n_parent_lists == 0) { pl = 0; parent_id = vertex_id; } else if (parent_num != NULL) { for (parent_id = parent_num[vertex_id] - 1, pl = n_parent_lists - 1; parent_id < parent_num_shift[pl]; pl--); assert(pl > -1); parent_id -= parent_num_shift[pl]; } else { for (parent_id = vertex_id, pl = n_parent_lists - 1 ; parent_id < parent_num_shift[pl] ; pl--); assert(pl > -1); parent_id -= parent_num_shift[pl]; } if (src_interlace == FVM_INTERLACE) { src_id = parent_id * src_dim + _src_dim_shift; } else { pl = src_dim*pl + _src_dim_shift; src_id = parent_id; } /* Now convert based on source datatype */ switch(src_datatype) { case FVM_FLOAT: v_f = ((const float *const *const)src_data)[pl][src_id]; break; case FVM_DOUBLE: v_f = ((const double *const *const)src_data)[pl][src_id]; break; default: assert(0); } a[0][0] += v_x * v_x; a[0][1] += v_x * v_y; a[0][2] += v_x * v_z; a[0][3] += v_x; a[1][1] += v_y * v_y; a[1][2] += v_y * v_z; a[1][3] += v_y; a[2][2] += v_z * v_z; a[2][3] += v_z; a[3][3] += 1.; b[0] += v_x * v_f; b[1] += v_y * v_f; b[2] += v_z * v_f; b[3] += v_f; } /* End of loop on element vertices */ /* Matrix is symmetric */ a[1][0] = a[0][1]; a[2][0] = a[0][2]; a[3][0] = a[0][3]; a[2][1] = a[1][2]; a[3][1] = a[1][3]; a[3][2] = a[2][3]; /* Find hyperplane coefficients and compute added value */ if (_solve_ax_b_4(a, b, coeff) == 0) { interpolated_value = ( coeff[0] * vertex_coords[0] + coeff[1] * vertex_coords[1] + coeff[2] * vertex_coords[2] + coeff[3]); } else { interpolated_value = v_f; /* last encountered value */ } /* Now convert based on destination datatype */ switch(dest_datatype) { case FVM_FLOAT: ((float *const)dest_data)[i] = interpolated_value; break; case FVM_DOUBLE: ((double *const)dest_data)[i] = interpolated_value; break; default: assert(0); } } /* End of loop on destination dimension */ } /* End of main loop on polyhedra */ if (heap != _heap) { BFT_FREE(heap); BFT_FREE(vertex_list); } } /*---------------------------------------------------------------------------- * Tesselate polygons (2D elements or 3D element faces) of a mesh section * referred to by an fvm_tesselation_t structure. Quadrangles are not * tesselated. * * parameters: * this_tesselation <-> partially initialized tesselation structure * dim <-- spatial dimension * vertex_coords <-- associated vertex coordinates array * parent_vertex_num <-- optional indirection to vertex coordinates * error_count --> number of triangulation errors counter (optional) *----------------------------------------------------------------------------*/ static void _tesselate_polygons(fvm_tesselation_t *this_tesselation, int dim, const fvm_coord_t vertex_coords[], const fvm_lnum_t parent_vertex_num[], fvm_lnum_t *error_count) { int type_id; fvm_lnum_t n_vertices, n_triangles, n_quads, n_elements; fvm_lnum_t n_vertices_max, n_triangles_max, n_unassigned_triangles; fvm_lnum_t i, j, k, l, v1, v2, vertex_id; fvm_gnum_t n_g_elements_tot[2] = {0, 0}; /* Global new elements count */ fvm_lnum_t n_elements_tot[2] = {0, 0}; /* New triangles/quadrangles */ fvm_lnum_t n_g_elements_max[2] = {0, 0}; /* Global max triangles/quadrangles */ fvm_lnum_t n_elements_max[2] = {0, 0}; /* Max triangles/quadrangles */ fvm_lnum_t *triangle_vertices = NULL; fvm_triangulate_state_t *state = NULL; fvm_tesselation_t *ts = this_tesselation; /* Initialization */ if (error_count != NULL) *error_count = 0; if (ts->type != FVM_CELL_POLY) n_elements = ts->n_elements; else n_elements = ts->n_faces; /* Count expected local numbers of triangles and quadrangles */ n_vertices_max = 0; n_triangles_max = 0; if (ts->vertex_index != NULL) { for (i = 0 ; i < n_elements ; i++) { n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i]; if (n_vertices == 4) n_elements_tot[1] += 1; else n_elements_tot[0] += n_vertices - 2; if (n_vertices > n_vertices_max) n_vertices_max = n_vertices; } } else return; n_triangles_max = n_vertices_max - 2; /* Destroy previous tesselation description if present */ ts->encoding = NULL; if (ts->_encoding != NULL) BFT_FREE(ts->_encoding); /* Allocate memory and state variables */ /*-------------------------------------*/ if (n_vertices_max > 4) { BFT_MALLOC(ts->_encoding, ts->vertex_index[n_elements], fvm_lnum_t); ts->encoding = ts->_encoding; BFT_MALLOC(triangle_vertices, (n_vertices_max - 2) * 3, fvm_lnum_t); state = fvm_triangulate_state_create(n_vertices_max); } n_elements_tot[0] = 0; n_elements_tot[1] = 0; /* reset */ /* Main loop on section face elements */ /*------------------------------------*/ for (i = 0 ; i < n_elements ; i++) { n_triangles = 0; n_quads = 0; n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i]; vertex_id = ts->vertex_index[i]; if (n_vertices == 4) type_id = 1; else type_id = 0; /* If face must be subdivided */ if (n_vertices > 4) { n_triangles = fvm_triangulate_polygon(dim, n_vertices, vertex_coords, parent_vertex_num, ( ts->vertex_num + vertex_id), triangle_vertices, state); if (n_triangles != (n_vertices - 2) && error_count != NULL) *error_count += 1; n_unassigned_triangles = n_triangles; for (j = 0; j < n_vertices; j++) { /* Vertices of corresponding edge */ v1 = ts->vertex_num[vertex_id + j]; v2 = ts->vertex_num[vertex_id + ((j+1) % n_vertices)]; /* Base initialization, in case no remaining triangle can be associated with this edge */ ts->_encoding[vertex_id + j] = -1; k = 0; while (k < n_unassigned_triangles) { for (l = 0; l < 3; l++) { if ( (triangle_vertices[k*3 + (l % 3)] == v1) && (triangle_vertices[k*3 + ((l+1) % 3)] == v2)) { ts->_encoding[vertex_id + j] = triangle_vertices[k*3 + ((l+2) % 3)]; break; } } /* Once found, remove triangle definition to avoid using it twice, (we should only have n-2 triangles for n vertices), and break from loop on triangles */ if (l < 3) { for (l = 0; l < 3; l++) triangle_vertices[k*3 + l] = triangle_vertices[(n_unassigned_triangles - 1)*3 + l]; n_unassigned_triangles--; break; } k++; } /* End of loop on unassigned triangles */ } /* End of loop on polygon vertices */ /* This situation may occur when a triangle shares no edge with the polygon boundary (in which case a hole is left). This will be corrected through a new encoding scheme. */ /* assert(n_unassigned_triangles == 0); */ n_elements_tot[0] += n_triangles; } /* Otherwise, tesselation trivial or not necessary for this face */ else { if (ts->_encoding != NULL) { for (j = ts->vertex_index[i]; j < ts->vertex_index[i+1]; j++) ts->_encoding[j] = -1; } if (n_vertices == 4) { n_elements_tot[1] += 1; n_quads = 1; } else if (n_vertices == 3) { n_elements_tot[0] += 1; n_triangles = 1; } } if (n_triangles > n_elements_max[0]) n_elements_max[0] = n_triangles; if (n_quads > n_elements_max[1]) n_elements_max[1] = n_triangles; } /* End of loop on elements */ /* Free memory and state variables */ if (n_vertices_max > 4) { BFT_FREE(triangle_vertices); state = fvm_triangulate_state_destroy(state); } /* Update tesselation structure info */ for (type_id = 0; type_id < 2; type_id++) { n_g_elements_tot[type_id] = n_elements_tot[type_id]; n_g_elements_max[type_id] = n_elements_max[type_id]; } fvm_parall_counter(n_g_elements_tot, 2); fvm_parall_counter_max(n_g_elements_max, 2); for (type_id = 0; type_id < 2; type_id++) { if (n_g_elements_tot[type_id] > 0) { if (type_id == 0) ts->sub_type[ts->n_sub_types] = FVM_FACE_TRIA; else ts->sub_type[ts->n_sub_types] = FVM_FACE_QUAD; ts->n_sub_max[ts->n_sub_types] = n_elements_max[type_id]; ts->n_sub_max_glob[ts->n_sub_types] = n_g_elements_max[type_id]; ts->n_sub[ts->n_sub_types] = n_elements_tot[type_id]; ts->n_sub_glob[ts->n_sub_types] = n_elements_tot[type_id]; ts->n_sub_types += 1; } } } /*---------------------------------------------------------------------------- * Build indexes and global counts for polygons of a mesh section referred * to by an fvm_tesselation_t structure. * * parameters: * this_tesselation <-> partially initialized tesselation structure * global_count <-- indicates if global counts should be built *----------------------------------------------------------------------------*/ static void _count_and_index_sub_polygons(fvm_tesselation_t *this_tesselation, _Bool global_count) { int sub_type_id, type_id; fvm_lnum_t n_vertices, n_triangles, n_elements; fvm_lnum_t i, j; fvm_lnum_t n_elements_tot[2] = {0, 0}; /* New triangles/quadrangles */ fvm_lnum_t *n_sub_elements[2] = {NULL, NULL}; fvm_tesselation_t *ts = this_tesselation; /* Initial checks */ if (ts->vertex_index == NULL) return; /* Initialization */ n_elements = ts->n_elements; /* Count expected local numbers of triangles and quadrangles */ for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) { type_id = -1; switch(ts->sub_type[sub_type_id]) { case FVM_FACE_TRIA: type_id = 0; break; case FVM_FACE_QUAD: type_id = 1; break; default: assert(0); } n_elements_tot[type_id] = ts->n_sub[sub_type_id]; BFT_MALLOC(ts->_sub_elt_index[sub_type_id], n_elements + 1, fvm_lnum_t); for (i = 0 ; i < n_elements + 1 ; i++) ts->_sub_elt_index[sub_type_id][i] = 0; /* The index array will first be used to hold n_sub_elements[], of size n_elements. To simplify future conversion to index, we shift shuch that n_sub_elements[i] corresponds to index[i+1]. */ n_sub_elements[type_id] = ts->_sub_elt_index[sub_type_id] + 1; } /* Loop on elements to populate n_sub_elements[]. Note that each n_sub_elements[] array has been initialized with zeroes, as it is mapped to a ts->sub_elt_index[] thus initialized. */ for (i = 0 ; i < n_elements ; i++) { n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i]; n_triangles = 0; if (n_vertices == 3) { n_sub_elements[0][i] = 1; n_triangles += 1; } else { /* if (n_vertices > 3) */ for (j = ts->vertex_index[i]; j < ts->vertex_index[i+1]; j++) { if (ts->encoding != NULL) { if (ts->encoding[j] > 0) n_triangles += 1; } } if (n_triangles > 0) n_sub_elements[0][i] = n_triangles; else if (n_vertices == 4) n_sub_elements[1][i] = 1; assert(n_triangles > 0 || n_vertices == 4); } } /* Now compute max global number if necessary */ for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) ts->n_sub_glob[sub_type_id] = ts->n_sub[sub_type_id]; if (global_count == true && ts->global_element_num != NULL) { for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) { if (ts->_sub_elt_index[sub_type_id] != NULL) { fvm_lnum_t * _n_sub_elements = ts->_sub_elt_index[sub_type_id] + 1; if (n_elements == 0) _n_sub_elements = NULL; ts->n_sub_glob[sub_type_id] = fvm_io_num_global_sub_size(ts->global_element_num, _n_sub_elements); } else ts->n_sub_glob[sub_type_id] = 0; } } /* Finally, build index from n_sub_elements_glob */ for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) { fvm_lnum_t *sub_elt_index = ts->_sub_elt_index[sub_type_id]; for (i = 0 ; i < n_elements ; i++) sub_elt_index[i+1] = sub_elt_index[i] + sub_elt_index[i+1]; ts->sub_elt_index[sub_type_id] = ts->_sub_elt_index[sub_type_id]; } } /*---------------------------------------------------------------------------- * Count resulting number of tetrahedra and pyramids when tesselating * polyhedra of a mesh section referred to by an fvm_tesselation_t * structure. This requires that the mesh section's polygonal faces have * already been tesselated (i.e. _tesselate_polygons has been called). * * parameters: * this_tesselation <-> partially initialized tesselation structure * error_count --> number of tesselation errors counter (optional) * global_count <-- indicates if global counts should be built *----------------------------------------------------------------------------*/ static void _count_and_index_sub_polyhedra(fvm_tesselation_t *this_tesselation, fvm_lnum_t *error_count, _Bool global_count) { int sub_type_id, type_id; fvm_lnum_t n_vertices, n_triangles, n_elements; fvm_lnum_t n_tetras, n_pyrams; fvm_lnum_t i, j, k, face_id; fvm_gnum_t n_g_elements_tot[2] = {0, 0}; /* Global new elements count */ fvm_lnum_t n_elements_tot[2] = {0, 0}; /* New tetrahedra/pyramids */ fvm_lnum_t *n_sub_elements[2] = {NULL, NULL}; fvm_lnum_t n_g_elements_max[2] = {0, 0}; /* Global max tetrahedra/pyramids */ fvm_lnum_t n_elements_max[2] = {0, 0}; /* Max tetrahedra/pyramids */ fvm_tesselation_t *ts = this_tesselation; /* Initial checks */ assert(ts->vertex_index != NULL); /* Initialization */ if (error_count != NULL) *error_count = 0; n_elements = ts->n_elements; /* Count expected total and local numbers of tetrahedra and pyramids; at this stage, the tesselation has been initialized as a polygon tesselation; adding a "center" vertex, triangle faces lead to tetrahedra, and quadrangles to pyramids */ for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) { type_id = -1; switch(ts->sub_type[sub_type_id]) { case FVM_FACE_TRIA: type_id = 0; break; case FVM_FACE_QUAD: type_id = 1; break; default: assert(0); } BFT_MALLOC(ts->_sub_elt_index[sub_type_id], n_elements + 1, fvm_lnum_t); for (i = 0 ; i < n_elements + 1 ; i++) ts->_sub_elt_index[sub_type_id][i] = 0; /* The index array will first be used to hold n_sub_elements[], of size n_elements. To simplify future conversion to index, we shift shuch that n_sub_elements[i] corresponds to index[i+1]. */ n_sub_elements[type_id] = ts->_sub_elt_index[sub_type_id] + 1; } /* Counting loop on polyhedra elements */ /*-------------------------------------*/ for (i = 0 ; i < n_elements ; i++) { n_tetras = 0; n_pyrams = 0; for (j = ts->face_index[i]; /* Loop on element faces */ j < ts->face_index[i+1]; j++) { face_id = FVM_ABS(ts->face_num[j]) - 1; n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id]; if (n_vertices == 3) n_tetras += 1; else { /* if (n_vertices > 3) */ n_triangles = 0; for (k = ts->vertex_index[face_id]; k < ts->vertex_index[face_id+1]; k++) { if (ts->encoding != NULL) { if (ts->encoding[k] > 0) n_triangles += 1; } } if (error_count != NULL && n_triangles < n_vertices - 2) *error_count += 1; if (n_triangles > 0) n_tetras += n_triangles; else if (n_vertices == 4) n_pyrams += 1; } } /* End of loop on element faces */ n_elements_tot[0] += n_tetras; n_elements_tot[1] += n_pyrams; if (n_tetras > n_elements_max[0]) n_elements_max[0] = n_tetras; if (n_pyrams > n_elements_max[1]) n_elements_max[1] = n_pyrams; if (n_sub_elements[0] != NULL) n_sub_elements[0][i] = n_tetras; if (n_sub_elements[1] != NULL) n_sub_elements[1][i] = n_pyrams; } /* End of loop on elements */ /* Update tesselation structure info */ for (type_id = 0; type_id < 2; type_id++) { n_g_elements_tot[type_id] = n_elements_tot[type_id]; n_g_elements_max[type_id] = n_elements_max[type_id]; } fvm_parall_counter(n_g_elements_tot, 2); fvm_parall_counter_max(n_g_elements_max, 2); ts->n_sub_types = 0; for (type_id = 0; type_id < 2; type_id++) { if (n_g_elements_tot[type_id] > 0) { if (type_id == 0) ts->sub_type[ts->n_sub_types] = FVM_CELL_TETRA; else ts->sub_type[ts->n_sub_types] = FVM_CELL_PYRAM; ts->n_sub_max[ts->n_sub_types] = n_elements_max[type_id]; ts->n_sub_max_glob[ts->n_sub_types] = n_g_elements_max[type_id]; ts->n_sub[ts->n_sub_types] = n_elements_tot[type_id]; ts->n_sub_glob[ts->n_sub_types] = n_elements_tot[type_id]; ts->n_sub_types += 1; } } /* Now compute max global number if necessary */ for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) ts->n_sub_glob[sub_type_id] = ts->n_sub[sub_type_id]; if (global_count == true && ts->global_element_num != NULL) { for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) { if (ts->_sub_elt_index[sub_type_id] != NULL) { fvm_lnum_t * _n_sub_elements = ts->_sub_elt_index[sub_type_id] + 1; if (n_elements == 0) _n_sub_elements = NULL; ts->n_sub_glob[sub_type_id] = fvm_io_num_global_sub_size(ts->global_element_num, _n_sub_elements); } else ts->n_sub_glob[sub_type_id] = 0; } } /* Finally, build index from n_sub_elements_glob */ for (sub_type_id = 0; sub_type_id < ts->n_sub_types; sub_type_id++) { fvm_lnum_t *sub_elt_index = ts->_sub_elt_index[sub_type_id]; for (i = 0 ; i < n_elements ; i++) sub_elt_index[i+1] = sub_elt_index[i] + sub_elt_index[i+1]; ts->sub_elt_index[sub_type_id] = ts->_sub_elt_index[sub_type_id]; } } #if defined(FVM_HAVE_MPI) /*---------------------------------------------------------------------------- * Update global_num_end when limiting partial connectivity range end index. * * parameters: * this_tesselation <-- tesselation structure * end_id <-- initial end of subset in parent section * global_num_end <-> past the end (maximum + 1) parent element * global number (reduced on return if required * by buffer_size) * comm <-- associated MPI communicator * * returns: * index corresponding to end of decoded range *----------------------------------------------------------------------------*/ static void _expand_limit_g(const fvm_tesselation_t *this_tesselation, fvm_lnum_t end_id, fvm_gnum_t *global_num_end, MPI_Comm comm) { fvm_gnum_t local_max, global_max; const fvm_gnum_t *global_element_num = fvm_io_num_get_global_num(this_tesselation->global_element_num); /* Check if the maximum id returned on some ranks leads to a lower global_num_end than initially required (due to local buffer being full) */ if (end_id < this_tesselation->n_elements) { if (global_element_num != NULL) local_max = global_element_num[end_id]; else local_max = end_id + 1; } else local_max = *global_num_end; MPI_Allreduce(&local_max, &global_max, 1, FVM_MPI_GNUM, MPI_MIN, comm); if (global_max < *global_num_end) *global_num_end = global_max; } #endif /* defined(FVM_HAVE_MPI) */ #if defined(FVM_HAVE_MPI) /*---------------------------------------------------------------------------- * Decode polygons tesselation to a connectivity buffer. * * To avoid requiring huge buffers and computing unneeded element * connectivities when exporting data in slices, this function may decode * a partial connectivity range, starting at polygon index start_id and ending * either when the indicated buffer size is attained, or the global element * number corresponding to a given polygon exceeds a given value. * It returns the effective polygon index end. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of polygons subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for sub_element_idx[] * and vertex_num[] buffers * global_num_end <-- past the end (maximum + 1) parent element * global number * global_vertex_num <-- global vertex numbering * vertex_num --> sub-element (global) vertex connectivity * * returns: * polygon index end corresponding to decoded range *----------------------------------------------------------------------------*/ static fvm_lnum_t _decode_polygons_tesselation_g(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_gnum_t global_num_end, const fvm_io_num_t *global_vertex_num, fvm_gnum_t vertex_num[]) { fvm_lnum_t n_vertices; fvm_lnum_t i, j, k, vertex_id; fvm_lnum_t n_sub_tot = 0; const fvm_tesselation_t *ts = this_tesselation; fvm_lnum_t retval = start_id; const fvm_gnum_t *_global_vertex_num = fvm_io_num_get_global_num(global_vertex_num); const fvm_gnum_t *global_element_num = fvm_io_num_get_global_num(ts->global_element_num); /* Main loop on polygons */ /*-----------------------*/ for (i = 0, j = start_id ; j < this_tesselation->n_elements; i++, j++) { if ( global_element_num != NULL && global_element_num[j] >= global_num_end) break; n_vertices = ts->vertex_index[j+1] - ts->vertex_index[j]; vertex_id = ts->vertex_index[j]; /* Sub-elements (triangles) connectivity */ /*---------------------------------------*/ if ( connect_type == FVM_FACE_TRIA && n_vertices > 3 && ts->encoding != NULL) { /* Loop on polygon vertices */ for (k = 0; k < n_vertices; k++) { if (ts->encoding[vertex_id + k] > 0) { /* Fill connectivity array */ /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return j; /* Fill connectivity array */ if (_global_vertex_num != NULL) { vertex_num[n_sub_tot*3 ] = _global_vertex_num[ts->vertex_num[vertex_id + k] - 1]; vertex_num[n_sub_tot*3 + 1] = _global_vertex_num[ts->vertex_num[ vertex_id + ((k+1) % n_vertices)] - 1]; vertex_num[n_sub_tot*3 + 2] = _global_vertex_num[ts->encoding[vertex_id + k] - 1]; } else { vertex_num[n_sub_tot*3 ] = ts->vertex_num[vertex_id + k]; vertex_num[n_sub_tot*3 + 1] = ts->vertex_num[vertex_id + ((k+1) % n_vertices)]; vertex_num[n_sub_tot*3 + 2] = ts->encoding[vertex_id + k]; } n_sub_tot += 1; } } /* End of loop on polygon vertices */ } /* Non-tesselated elements connectivity */ /*--------------------------------------*/ else if ( (connect_type == FVM_FACE_TRIA && n_vertices == 3) || (connect_type == FVM_FACE_QUAD && n_vertices == 4)) { /* Check that polygon is not subdivided */ k = n_vertices; if (ts->encoding != NULL) { for (k = 0; k < n_vertices; k++) if (ts->encoding[vertex_id + k] > 0) break; } if (k == n_vertices) { /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return j; /* Fill connectivity array */ if (_global_vertex_num != NULL) { for (k = 0; k < n_vertices; k++) vertex_num[n_sub_tot*n_vertices +k] = _global_vertex_num[ts->vertex_num[vertex_id + k] - 1]; } else { for (k = 0; k < n_vertices; k++) vertex_num[n_sub_tot*n_vertices +k] = ts->vertex_num[vertex_id + k]; } n_sub_tot += 1; } } /* End of case where polygon is not tesselated */ retval = j+1; /* If end of buffer has not already been reached */ } /* End of loop on polygons */ return retval; } #endif /* defined(FVM_HAVE_MPI) */ /*---------------------------------------------------------------------------- * Decode polygons tesselation to a connectivity buffer. * * To avoid requiring huge buffers and computing unneeded element * connectivities this function may decode a partial connectivity range, * starting at polygon index start_id and ending either when the indicated * buffer size or the last polygon is attained. * It returns the effective polygon index end. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of polygons subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for vertex_num[] buffer * vertex_num --> sub-element (global) vertex connectivity * * returns: * polygon index end corresponding to decoded range *----------------------------------------------------------------------------*/ static fvm_lnum_t _decode_polygons_tesselation_l(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_lnum_t vertex_num[]) { fvm_lnum_t n_vertices; fvm_lnum_t i, j, vertex_id; fvm_lnum_t n_sub_tot = 0; const fvm_tesselation_t *ts = this_tesselation; fvm_lnum_t retval = start_id; /* Main loop on polygons */ /*-----------------------*/ for (i = start_id ; i < this_tesselation->n_elements; i++) { n_vertices = ts->vertex_index[i+1] - ts->vertex_index[i]; vertex_id = ts->vertex_index[i]; /* Sub-elements (triangles) connectivity */ /*---------------------------------------*/ if ( connect_type == FVM_FACE_TRIA && n_vertices > 3 && ts->encoding != NULL) { /* Loop on polygon vertices */ for (j = 0; j < n_vertices; j++) { if (ts->encoding[vertex_id + j] > 0) { /* Fill connectivity array */ /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return i; /* Fill connectivity array */ vertex_num[n_sub_tot*3 ] = ts->vertex_num[vertex_id + j]; vertex_num[n_sub_tot*3 + 1] = ts->vertex_num[vertex_id + ((j+1) % n_vertices)]; vertex_num[n_sub_tot*3 + 2] = ts->encoding[vertex_id + j]; n_sub_tot += 1; } } /* End of loop on polygon vertices */ } /* Non-tesselated elements connectivity */ /*--------------------------------------*/ else if ( (connect_type == FVM_FACE_TRIA && n_vertices == 3) || (connect_type == FVM_FACE_QUAD && n_vertices == 4)) { /* Check that polygon is not subdivided */ j = n_vertices; if (ts->encoding != NULL) { for (j = 0; j < n_vertices; j++) if (ts->encoding[vertex_id + j] > 0) break; } if (j == n_vertices) { /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return i; /* Fill connectivity array */ for (j = 0; j < n_vertices; j++) vertex_num[n_sub_tot*n_vertices +j] = ts->vertex_num[vertex_id + j]; n_sub_tot += 1; } } /* End of case where polygon is not tesselated */ retval = i+1; /* If end of buffer has not already been reached */ } /* End of loop on polygons */ return retval; } #if defined(FVM_HAVE_MPI) /*---------------------------------------------------------------------------- * Decode polyhedra tesselation to a connectivity buffer. * * To avoid requiring huge buffers and computing unneeded element * connectivities when exporting data in slices, this function may decode a * partial connectivity range, starting at polyhedron index start_id and ending * either when the indicated buffer size is attained, or the global element * number corresponding to a given polyhedron exceeds a given value. * It returns the effective polyhedron index end. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of polyhedra subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for vertex_num[] buffer * global_num_end <-- past the end (maximum + 1) parent element * global number * extra_vertex_base <-- starting number for added vertices * global_vertex_num <-- global vertex numbering * vertex_num --> sub-element (global) vertex connectivity * * returns: * polyhedron index end corresponding to decoded range *----------------------------------------------------------------------------*/ static fvm_lnum_t _decode_polyhedra_tesselation_g(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_gnum_t global_num_end, fvm_gnum_t extra_vertex_base_num, const fvm_io_num_t *global_vertex_num, fvm_gnum_t vertex_num[]) { int orient; fvm_lnum_t n_vertices; fvm_lnum_t i, j, k, l, base_dest_id, face_id, vertex_id; fvm_lnum_t n_sub_tot = 0; const fvm_tesselation_t *ts = this_tesselation; fvm_lnum_t retval = start_id; const fvm_gnum_t *_global_vertex_num = fvm_io_num_get_global_num(global_vertex_num); const fvm_gnum_t *global_element_num = fvm_io_num_get_global_num(ts->global_element_num); /* Main loop on polyhedra */ /*------------------------*/ for (i = 0, j = start_id ; j < this_tesselation->n_elements; i++, j++) { if ( global_element_num != NULL && global_element_num[j] >= global_num_end) break; for (k = ts->face_index[j]; /* Loop on element faces */ k < ts->face_index[j+1]; k++) { face_id = FVM_ABS(ts->face_num[k]) - 1; n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id]; vertex_id = ts->vertex_index[face_id]; /* Invert base orientation as it points outwards in polyhedron definition, but inwards in tetrahedron and pyramid reference connectivity */ if (ts->face_num[k] > 0) orient = -1; else orient = 1; /* Sub-elements (tetrahedra) connectivity */ /*----------------------------------------*/ if ( connect_type == FVM_CELL_TETRA && n_vertices > 3 && ts->encoding != NULL) { /* Loop on face vertices */ for (l = 0; l < n_vertices; l++) { if (ts->encoding[vertex_id + l] > 0) { /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return j; if (orient == 1) base_dest_id = n_sub_tot*4; else base_dest_id = n_sub_tot*4 + 2; /* Fill connectivity array */ if (_global_vertex_num != NULL) { vertex_num[base_dest_id] = _global_vertex_num[ts->vertex_num[vertex_id + l] - 1]; vertex_num[base_dest_id + orient] = _global_vertex_num[ts->vertex_num[ vertex_id + ((l+1) % n_vertices)] - 1]; vertex_num[base_dest_id + 2*orient] = _global_vertex_num[ts->encoding[vertex_id + l] - 1]; } else { vertex_num[base_dest_id] = ts->vertex_num[vertex_id + l]; vertex_num[base_dest_id + orient] = ts->vertex_num[vertex_id + ((l+1) % n_vertices)]; vertex_num[base_dest_id + 2*orient] = ts->encoding[vertex_id + l]; } /* Add extra "top" vertex */ if (global_element_num != NULL) vertex_num[n_sub_tot*4 + 3] = extra_vertex_base_num + global_element_num[j] - 1; else vertex_num[n_sub_tot*4 + 3] = extra_vertex_base_num + j; n_sub_tot += 1; } } } /* Non-tesselated faces connectivity */ /*-----------------------------------*/ else if ( (connect_type == FVM_CELL_TETRA && n_vertices == 3) || (connect_type == FVM_CELL_PYRAM && n_vertices == 4)) { /* Check that face is not subdivided */ l = n_vertices; if (ts->encoding != NULL) { for (l = 0; l < n_vertices; l++) if (ts->encoding[vertex_id + l] > 0) break; } if (l == n_vertices) { fvm_lnum_t stride = n_vertices + 1; /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return j; if (orient == 1) base_dest_id = n_sub_tot*stride; else base_dest_id = n_sub_tot*stride + n_vertices-1; /* Fill connectivity array */ if (_global_vertex_num != NULL) { for (l = 0; l < n_vertices; l++) vertex_num[base_dest_id + l*orient] = _global_vertex_num[ts->vertex_num[vertex_id + l] - 1]; } else { for (l = 0; l < n_vertices; l++) vertex_num[base_dest_id + l*orient] = ts->vertex_num[vertex_id + l]; } /* Add extra "top" vertex */ if (global_element_num != NULL) vertex_num[n_sub_tot*stride + n_vertices] = extra_vertex_base_num + global_element_num[j] - 1; else vertex_num[n_sub_tot*stride + n_vertices] = extra_vertex_base_num + j; n_sub_tot += 1; } } /* End of case were face is not tesselated */ } /* End of loop on element faces */ retval = j+1; /* If end of buffer has not already been reached */ } /* End of main loop on polyhedra */ return retval; } #endif /* defined(FVM_HAVE_MPI) */ /*---------------------------------------------------------------------------- * Decode polyhedra tesselation to a connectivity buffer. * * To avoid requiring huge buffers, this function may decode a partial * connectivity range, starting at polyhedron index start_id and ending either * when the indicated buffer size or the last polyhedra is attained. * It returns the effective polyhedron index end. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of polyhedra subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for vertex_num[] buffer * extra_vertex_base <-- starting number for added vertices * vertex_num --> sub-element (global) vertex connectivity * * returns: * polyhedron index end corresponding to decoded range *----------------------------------------------------------------------------*/ static fvm_lnum_t _decode_polyhedra_tesselation_l(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_lnum_t extra_vertex_base_num, fvm_lnum_t vertex_num[]) { int orient; fvm_lnum_t n_vertices; fvm_lnum_t i, j, k, l, base_dest_id, face_id, vertex_id; fvm_lnum_t n_sub_tot = 0; const fvm_tesselation_t *ts = this_tesselation; fvm_lnum_t retval = start_id; /* Main loop on polyhedra */ /*------------------------*/ for (i = 0, j = start_id ; j < this_tesselation->n_elements; i++, j++) { for (k = ts->face_index[j]; /* Loop on element faces */ k < ts->face_index[j+1]; k++) { face_id = FVM_ABS(ts->face_num[k]) - 1; n_vertices = ts->vertex_index[face_id+1] - ts->vertex_index[face_id]; vertex_id = ts->vertex_index[face_id]; /* Invert base orientation as it points outwards in polyhedron definition, but inwards in tetrahedron and pyramid reference connectivity */ if (ts->face_num[k] > 0) orient = -1; else orient = 1; /* Sub-elements (tetrahedra) connectivity */ /*----------------------------------------*/ if ( connect_type == FVM_CELL_TETRA && n_vertices > 3 && ts->encoding != NULL) { /* Loop on face vertices */ for (l = 0; l < n_vertices; l++) { if (ts->encoding[vertex_id + l] > 0) { /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return j; if (orient == 1) base_dest_id = n_sub_tot*4; else base_dest_id = n_sub_tot*4 + 2; /* Fill connectivity array */ vertex_num[base_dest_id] = ts->vertex_num[vertex_id + l]; vertex_num[base_dest_id + orient] = ts->vertex_num[vertex_id + ((l+1) % n_vertices)]; vertex_num[base_dest_id + 2*orient] = ts->encoding[vertex_id + l]; /* Add extra "top" vertex */ vertex_num[n_sub_tot*4 + 3] = extra_vertex_base_num + j; n_sub_tot += 1; } } } /* Non-tesselated faces connectivity */ /*-----------------------------------*/ else if ( (connect_type == FVM_CELL_TETRA && n_vertices == 3) || (connect_type == FVM_CELL_PYRAM && n_vertices == 4)) { /* Check that face is not subdivided */ l = n_vertices; if (ts->encoding != NULL) { for (l = 0; l < n_vertices; l++) if (ts->encoding[vertex_id + l] > 0) break; } if (l == n_vertices) { fvm_lnum_t stride = n_vertices + 1; /* Return previous element's end index if buffer size reached */ if (n_sub_tot >= buffer_limit) return j; if (orient == 1) base_dest_id = n_sub_tot*stride; else base_dest_id = n_sub_tot*stride + n_vertices-1; /* Fill connectivity array */ for (l = 0; l < n_vertices; l++) vertex_num[base_dest_id + l*orient] = ts->vertex_num[vertex_id + l]; /* Add extra "top" vertex */ vertex_num[n_sub_tot*stride + n_vertices] = extra_vertex_base_num + j; n_sub_tot += 1; } } /* End of case were face is not tesselated */ } /* End of loop on element faces */ retval = j+1; /* If end of buffer has not already been reached */ } /* End of main loop on polyhedra */ return retval; } /*============================================================================ * Public function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Creation of a mesh section's subdivision into simpler elements. * * The structure contains pointers to the mesh section's connectivity, * (passed as arguments), which is not copied. This structure should thus * always be destroyed before the mesh section to which it relates. * * Unused connectivity array arguments should be set to NULL (such as * face_index[] and face_num[] for 2D or regular (strided) elements, * and vertex_index[] for strided elements. * * At this stage, the structure does not yet contain tesselation information. * * parameters: * element_type <-- type of elements considered * n_elements <-- number of elements * face_index <-- polyhedron -> faces index (O to n-1) * dimension [n_elements + 1] * face_num <-- element -> face numbers (1 to n, signed, * > 0 for outwards pointing face normal * < 0 for inwards pointing face normal); * dimension: [face_index[n_elements]], or NULL * vertex_index <-- element face -> vertices index (O to n-1); * dimension: [n_cell_faces + 1], [n_elements + 1], * or NULL depending on face_index and vertex_index * vertex_num <-- element -> vertex connectivity (1 to n) * global_element_num <-- global element numbers (NULL in serial mode) * * returns: * pointer to created mesh section tesselation structure *----------------------------------------------------------------------------*/ fvm_tesselation_t * fvm_tesselation_create(fvm_element_t element_type, fvm_lnum_t n_elements, const fvm_lnum_t face_index[], const fvm_lnum_t face_num[], const fvm_lnum_t vertex_index[], const fvm_lnum_t vertex_num[], const fvm_io_num_t *global_element_num) { int i; int entity_dim = 0, stride = 0; fvm_tesselation_t *this_tesselation; /* First, check if element type is handled and initialize additionnal info */ switch (element_type) { case FVM_FACE_QUAD: entity_dim = 2; stride = 4; break; case FVM_FACE_POLY: entity_dim = 2; stride = 0; break; case FVM_CELL_POLY: entity_dim = 3; stride = 0; break; default: return NULL; } /* Now, create structure */ BFT_MALLOC(this_tesselation, 1, fvm_tesselation_t); /* Assign parent mesh section info */ this_tesselation->type = element_type; this_tesselation->n_elements = n_elements; this_tesselation->dim = 0; this_tesselation->entity_dim = entity_dim; this_tesselation->stride = stride; this_tesselation->n_faces = 0; this_tesselation->vertex_coords = NULL; this_tesselation->parent_vertex_num = NULL; this_tesselation->face_index = face_index; this_tesselation->face_num = face_num; this_tesselation->vertex_index = vertex_index; this_tesselation->vertex_num = vertex_num; this_tesselation->global_element_num = global_element_num; /* Check argument consistency */ if (face_index != NULL || face_num != NULL) { if (element_type != FVM_CELL_POLY) bft_error(__FILE__, __LINE__, 0, _("Incoherent connectivity for tesselation:\n" "Connectivity face_index or face_num non NULL,\n" "but element type != FVM_CELL_POLY")); } else if (vertex_index != NULL) { if (element_type != FVM_FACE_POLY) bft_error(__FILE__, __LINE__, 0, _("Incoherent connectivity for tesselation:\n" "Connectivy vertex_index non NULL,\n" "but element type != FVM_FACE_POLY")); } /* Compute number of polyhedron faces */ if (n_elements > 0 && face_index != NULL) { fvm_lnum_t j, k, face_id; fvm_lnum_t max_face_id = 0; for (j = 0; j < n_elements; j++) { for (k = face_index[j]; k < face_index[j+1]; k++) { face_id = FVM_ABS(face_num[k]) - 1; if (face_id > max_face_id) max_face_id = face_id; } } this_tesselation->n_faces = max_face_id + 1; } /* Set tesselation structures to zero */ this_tesselation->n_sub_types = 0; for (i = 0; i < FVM_TESSELATION_N_SUB_TYPES_MAX; i++) { this_tesselation->n_sub_max[i] = 0; this_tesselation->n_sub_max_glob[i] = 0; this_tesselation->n_sub[i] = 0; this_tesselation->n_sub_glob[i] =0; this_tesselation->sub_type[i] = FVM_N_ELEMENT_TYPES; } this_tesselation->encoding = NULL; this_tesselation->_encoding = NULL; for (i = 0; i < FVM_TESSELATION_N_SUB_TYPES_MAX; i++) { this_tesselation->sub_elt_index[i] = NULL; this_tesselation->_sub_elt_index[i] = NULL; } return (this_tesselation); } /*---------------------------------------------------------------------------- * Destruction of a nodal mesh polygon splitting representation structure. * * parameters: * this_tesselation <-> pointer to structure that should be destroyed * * returns: * NULL pointer *----------------------------------------------------------------------------*/ fvm_tesselation_t * fvm_tesselation_destroy(fvm_tesselation_t * this_tesselation) { int i; if (this_tesselation->_encoding != NULL) BFT_FREE(this_tesselation->_encoding); for (i = 0; i < this_tesselation->n_sub_types; i++) { if (this_tesselation->_sub_elt_index[i] != NULL) BFT_FREE(this_tesselation->_sub_elt_index[i]); } BFT_FREE(this_tesselation); return NULL; } /*---------------------------------------------------------------------------- * Tesselate a mesh section referred to by an fvm_tesselation_t structure. * * parameters: * this_tesselation <-> partially initialized tesselation structure * dim <-- spatial dimension * vertex_coords <-- associated vertex coordinates array * parent_vertex_num <-- optional indirection to vertex coordinates * error_count --> number of elements with a tesselation error * counter (optional) *----------------------------------------------------------------------------*/ void fvm_tesselation_init(fvm_tesselation_t *this_tesselation, int dim, const fvm_coord_t vertex_coords[], const fvm_lnum_t parent_vertex_num[], fvm_lnum_t *error_count) { assert(this_tesselation != NULL); this_tesselation->dim = dim; this_tesselation->vertex_coords = vertex_coords; this_tesselation->parent_vertex_num = parent_vertex_num; switch(this_tesselation->type) { case FVM_CELL_POLY: _tesselate_polygons(this_tesselation, dim, vertex_coords, parent_vertex_num, error_count); _count_and_index_sub_polyhedra(this_tesselation, error_count, true); break; case FVM_FACE_POLY: _tesselate_polygons(this_tesselation, dim, vertex_coords, parent_vertex_num, error_count); _count_and_index_sub_polygons(this_tesselation, true); break; default: bft_error(__FILE__, __LINE__, 0, _("Tesselation of element type %s not implemented."), fvm_elements_type_name[this_tesselation->type]); } } /*---------------------------------------------------------------------------- * Reduction of a nodal mesh polygon splitting representation structure; * only the associations (numberings) necessary to redistribution of fields * for output are conserved, the full connectivity being no longer useful * once it has been output. * * parameters: * this_tesselation <-> pointer to structure that should be reduced *----------------------------------------------------------------------------*/ void fvm_tesselation_reduce(fvm_tesselation_t * this_tesselation) { this_tesselation->stride = 0; this_tesselation->n_faces = 0; if (this_tesselation->face_index == NULL) { this_tesselation->face_num = NULL; this_tesselation->vertex_index = NULL; this_tesselation->vertex_num = NULL; } this_tesselation->encoding = NULL; if (this_tesselation->_encoding != NULL) BFT_FREE(this_tesselation->_encoding); } /*---------------------------------------------------------------------------- * Return number of parent elements of a tesselation. * * parameters: * this_tesselation <-- tesselation structure * * returns: * number of parent elements *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_tesselation_n_elements(const fvm_tesselation_t *this_tesselation) { fvm_lnum_t retval = 0; if (this_tesselation != NULL) retval = this_tesselation->n_elements; return retval; } /*---------------------------------------------------------------------------- * Return global number of added vertices associated with a tesselation. * * parameters: * this_tesselation <-- tesselation structure * * returns: * global number of added vertices associated with the tesselation *----------------------------------------------------------------------------*/ fvm_gnum_t fvm_tesselation_n_g_vertices_add(const fvm_tesselation_t *this_tesselation) { fvm_gnum_t retval = 0; assert(this_tesselation != NULL); if (this_tesselation->type == FVM_CELL_POLY) { if (this_tesselation->global_element_num != NULL) retval = fvm_io_num_get_global_count(this_tesselation->global_element_num); else retval = this_tesselation->n_elements; } return retval; } /*---------------------------------------------------------------------------- * Return (local) number of added vertices associated with a tesselation. * * parameters: * this_tesselation <-- tesselation structure * * returns: * global number of added vertices associated with the tesselation *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_tesselation_n_vertices_add(const fvm_tesselation_t *this_tesselation) { fvm_gnum_t retval = 0; assert(this_tesselation != NULL); if (this_tesselation->type == FVM_CELL_POLY) retval = this_tesselation->n_elements; return retval; } /*---------------------------------------------------------------------------- * Return number of resulting sub-types of a tesselation. * * parameters: * this_tesselation <-- tesselation structure * * returns: * number of resulting sub-types of the tesselation *----------------------------------------------------------------------------*/ int fvm_tesselation_n_sub_types(const fvm_tesselation_t *this_tesselation) { int retval = 0; if (this_tesselation != NULL) retval = this_tesselation->n_sub_types; return retval; } /*---------------------------------------------------------------------------- * Return given sub-types of a tesselation. * * parameters: * this_tesselation <-- tesselation structure * sub_type_id <-- index of sub-type in tesselation (0 to n-1) * * returns: * sub-types of the tesselation with the given index *----------------------------------------------------------------------------*/ fvm_element_t fvm_tesselation_sub_type(const fvm_tesselation_t *this_tesselation, int sub_type_id) { fvm_element_t retval = FVM_N_ELEMENT_TYPES; if (this_tesselation == NULL) retval = 0; else { assert(sub_type_id < this_tesselation->n_sub_types); retval = this_tesselation->sub_type[sub_type_id]; } return retval; } /*---------------------------------------------------------------------------- * Return number of elements of a given sub-type of a tesselation. * * parameters: * this_tesselation <-- tesselation structure * sub_type_id <-- index of sub-type in tesselation (0 to n-1) * * returns: * sub-types of the tesselation with the given index *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_tesselation_n_sub_elements(const fvm_tesselation_t *this_tesselation, fvm_element_t sub_type) { int id; fvm_lnum_t retval = 0; if (this_tesselation != NULL) { for (id = 0; id < this_tesselation->n_sub_types; id++) { if (this_tesselation->sub_type[id] == sub_type) { retval = this_tesselation->n_sub[id]; break; } } } return retval; } /*---------------------------------------------------------------------------- * Obtain the global and maximum number of elements of a given sub-type * of a tesselation. * * parameters: * this_tesselation <-- tesselation structure * sub_type_id <-- index of sub-type in tesselation (0 to n-1) * n_sub_elements_glob --> global number of sub-elements of the given type * n_sub_elements_max --> maximum number of sub-elements per element * of the given type (for all ranks) *----------------------------------------------------------------------------*/ void fvm_tesselation_get_global_size(const fvm_tesselation_t *this_tesselation, fvm_element_t sub_type, fvm_gnum_t *n_sub_elements_glob, fvm_lnum_t *n_sub_elements_max) { int id; if (n_sub_elements_max != NULL) *n_sub_elements_max = 0; if (n_sub_elements_glob != NULL) *n_sub_elements_glob = 0; if (this_tesselation != NULL) { for (id = 0; id < this_tesselation->n_sub_types; id++) { if (this_tesselation->sub_type[id] == sub_type) { if (n_sub_elements_max != NULL) *n_sub_elements_max = this_tesselation->n_sub_max_glob[id]; if (n_sub_elements_glob != NULL) *n_sub_elements_glob = this_tesselation->n_sub_glob[id]; break; } } } } /*---------------------------------------------------------------------------- * Return global numbering of added vertices associated with a tesselation. * * parameters: * this_tesselation <-- tesselation structure * * returns: * pointer to global numbering of added vertices for this tesselation, * or NULL if no added vertices are present. *----------------------------------------------------------------------------*/ const fvm_io_num_t * fvm_tesselation_global_vertex_num(const fvm_tesselation_t *this_tesselation) { const fvm_io_num_t *retval = NULL; assert(this_tesselation != NULL); if (this_tesselation->type == FVM_CELL_POLY) retval = this_tesselation->global_element_num; return retval; } /*---------------------------------------------------------------------------- * Compute coordinates of added vertices for a tesselation of polyhedra. * * One additional vertex is added near the center of each polyhedra. * For element types other than polyhedra, there is no need for added * vertices, so this function returns immediately. * * parameters: * this_tesselation <-- tesselation structure * vertex_coords --> coordinates of added vertices *----------------------------------------------------------------------------*/ void fvm_tesselation_vertex_coords(const fvm_tesselation_t *this_tesselation, fvm_coord_t vertex_coords[]) { fvm_lnum_t i; if (this_tesselation->type != FVM_CELL_POLY) return; for (i = 0; i < this_tesselation->n_elements; i++) { _added_vertex_coords(this_tesselation, vertex_coords + i*3, NULL, i); } } /*---------------------------------------------------------------------------- * Return index of sub-elements associated with each element of a given * sub-type of a tesselation. * * parameters: * this_tesselation <-- tesselation structure * sub_type_id <-- index of sub-type in tesselation (0 to n-1) * * returns: * index of sub-elements associated with each element (0 to n-1 numbering) *----------------------------------------------------------------------------*/ const fvm_lnum_t * fvm_tesselation_sub_elt_index(const fvm_tesselation_t *this_tesselation, fvm_element_t sub_type) { int id; const fvm_lnum_t *retval = NULL; if (this_tesselation != NULL) { for (id = 0; id < this_tesselation->n_sub_types; id++) { if (this_tesselation->sub_type[id] == sub_type) { retval = this_tesselation->sub_elt_index[id]; break; } } } return retval; } #if defined(FVM_HAVE_MPI) /*---------------------------------------------------------------------------- * Compute index values corresponding to given range of indices, * for an element -> sub-element value distribution. * * This index is used mainly to gather a decoded tesselation connectivity * or element -> sub-element data distribution, expanding the corresponding * data only on the given range. * Only the index values in the start_id to end_id range are set by * this function, starting with index[start_id] = 0. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * stride <-- number of associated values per sub-element * start_id <-- start index of polyhedra subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for vertex_num[] buffer * global_num_end <-> past the end (maximum + 1) parent element * global number (reduced on return if required * by buffer_size limits) * index --> sub-element index * comm <-- associated MPI communicator * * returns: * polyhedron index end corresponding to decoded range *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_tesselation_range_index_g(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, int stride, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_gnum_t *global_num_end, fvm_lnum_t index[], MPI_Comm comm) { fvm_lnum_t i; fvm_lnum_t buffer_size = buffer_limit * stride; const fvm_gnum_t *global_element_num = fvm_io_num_get_global_num(this_tesselation->global_element_num); const fvm_lnum_t *sub_element_idx = fvm_tesselation_sub_elt_index(this_tesselation, connect_type); /* In parallel mode, global_element_num should exist */ if (index == NULL) return start_id; /* Loop on tesselated elements */ /*-----------------------------*/ index[start_id] = 0; for (i = start_id; i < this_tesselation->n_elements; i++) { if (global_element_num[i] >= *global_num_end) break; index[i+1] = index[i] + ( sub_element_idx[i+1] - sub_element_idx[i]) * stride; if (index[i+1] > buffer_size) { *global_num_end = global_element_num[i]; break; } } /* Check if the maximum id returned on some ranks leads to a lower global_num_end than initially required (due to local buffer being full) */ _expand_limit_g(this_tesselation, i, global_num_end, comm); return i; } /*---------------------------------------------------------------------------- * Decode tesselation to a parent element->sub elements index and * connectivity buffer. * * To avoid requiring huge buffers and computing unneeded element * connectivities when exporting data in slices, this function may decode * a partial connectivity range, starting at polygon index start_id and ending * either when the indicated buffer size is attained, or the global element * number corresponding to a given polygon exceeds a given value. * It returns the effective polygon index end. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of polygons subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for sub_element_idx[] * and vertex_num[] buffers * global_num_end <-> past the end (maximum + 1) parent element * global number (reduced on return if required * by buffer_limit) * extra_vertex_base <-- starting number for added vertices * global_vertex_num <-- global vertex numbering * vertex_num --> sub-element (global) vertex connectivity * comm <-- associated MPI communicator * * returns: * polygon index corresponding to end of decoded range *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_tesselation_decode_g(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_gnum_t *global_num_end, const fvm_io_num_t *global_vertex_num, fvm_gnum_t extra_vertex_base, fvm_gnum_t vertex_num[], MPI_Comm comm) { fvm_lnum_t retval = 0; switch(this_tesselation->type) { case FVM_CELL_POLY: retval = _decode_polyhedra_tesselation_g(this_tesselation, connect_type, start_id, buffer_limit, *global_num_end, extra_vertex_base, global_vertex_num, vertex_num); break; case FVM_FACE_POLY: retval = _decode_polygons_tesselation_g(this_tesselation, connect_type, start_id, buffer_limit, *global_num_end, global_vertex_num, vertex_num); break; default: assert(0); } /* Check if the maximum id returned on some ranks leads to a lower global_num_end than initially required (due to local buffer being full) */ _expand_limit_g(this_tesselation, retval, global_num_end, comm); return retval; } #endif /* defined(FVM_HAVE_MPI) */ /*---------------------------------------------------------------------------- * Decode tesselation to a connectivity buffer. * * To avoid requiring huge buffers and computing unneeded element * connectivities, this function may decode a partial connectivity range, * starting at polygon index start_id and ending either when the indicated * buffer size or the last polygon is attained. * It returns the effective polygon index end. * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of polygons subset in parent section * buffer_limit <-- maximum number of sub-elements of destination * element type allowable for sub_element_idx[] * and vertex_num[] buffers * extra_vertex_base <-- starting number for added vertices * vertex_num --> sub-element (global) vertex connectivity * * returns: * polygon index corresponding to end of decoded range *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_tesselation_decode(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t buffer_limit, fvm_lnum_t extra_vertex_base, fvm_lnum_t vertex_num[]) { fvm_lnum_t retval = 0; switch(this_tesselation->type) { case FVM_CELL_POLY: retval = _decode_polyhedra_tesselation_l(this_tesselation, connect_type, start_id, buffer_limit, extra_vertex_base, vertex_num); break; case FVM_FACE_POLY: retval = _decode_polygons_tesselation_l(this_tesselation, connect_type, start_id, buffer_limit, vertex_num); break; default: assert(0); } return retval; } /*---------------------------------------------------------------------------- * Distribute "per element" data from the base elements to their tesselation. * * The same data array is used for input and output, so as to avoid requiring * excess allocation in typical use cases (extracting data from a parent mesh * to a buffer and distributing it as per its tesselation). * The data array should be at least of size: * [sub_elt_index[end_id] - sub_elt_index[start_id] * size * * parameters: * this_tesselation <-- tesselation structure * connect_type <-- destination element type * start_id <-- start index of elements subset in parent section * end_id <-- end index of elements subset in parent section * size <-- data size for each element (sizeof(type)*stride) * data <-> undistributed data in, distributed data out *----------------------------------------------------------------------------*/ void fvm_tesselation_distribute(const fvm_tesselation_t *this_tesselation, fvm_element_t connect_type, fvm_lnum_t start_id, fvm_lnum_t end_id, size_t size, void *data) { int id; fvm_lnum_t i, j, k, n_sub; size_t l; char *src, *dest; const fvm_lnum_t *sub_elt_index = NULL; /* Find index, or return */ if (this_tesselation == NULL) return; for (id = 0; id < this_tesselation->n_sub_types; id++) { if (this_tesselation->sub_type[id] == connect_type) { sub_elt_index = this_tesselation->sub_elt_index[id]; break; } } if (id == this_tesselation->n_sub_types) return; /* Distribute data (starting from the end so as not to overwrite data at the beginning of the array) */ for (i = end_id, j = end_id - start_id - 1; i > start_id; i--, j--) { src = ((char *)data) + j*size; dest = ((char *)data) + (sub_elt_index[i-1] - sub_elt_index[start_id])*size; n_sub = sub_elt_index[i] - sub_elt_index[i-1]; for (k = 0; k < n_sub; k++) { for (l = 0; l < size; l++) dest[k*size + l] = src[l]; } } } /*---------------------------------------------------------------------------- * Compute field values at added vertices for a tesselation of polyhedra. * * One additional vertex is added near the center of each polyhedra. * For element types other than polyhedra, there is no need for added * vertices, so this function returns immediately. * * parameters: * this_tesselation <-- tesselation structure * vertex_coords <-- coordinates of added vertices * src_dim <-- dimension of source data * src_dim_shift <-- source data dimension shift (start index) * dest_dim <-- destination data dimension (1 if non interlaced) * start_id <-- added vertices start index * end_id <-- added vertices past the end index * src_interlace <-- indicates if source data is interlaced * src_datatype <-- source data type (float, double, or int) * dest_datatype <-- destination data type (float, double, or int) * n_parent_lists <-- number of parent lists (if parent_num != NULL) * parent_num_shift <-- parent number to value array index shifts; * size: n_parent_lists * parent_num <-- if n_parent_lists > 0, parent entity numbers * src_data <-- array of source arrays (at least one, with one per * source dimension if non interlaced, times one per * parent list if multiple parent lists, with * x_parent_1, y_parent_1, ..., x_parent_2, ...) order * dest_data --> destination buffer *----------------------------------------------------------------------------*/ void fvm_tesselation_vertex_values(const fvm_tesselation_t *this_tesselation, int src_dim, int src_dim_shift, int dest_dim, fvm_lnum_t start_id, fvm_lnum_t end_id, fvm_interlace_t src_interlace, fvm_datatype_t src_datatype, fvm_datatype_t dest_datatype, int n_parent_lists, const fvm_lnum_t parent_num_shift[], const fvm_lnum_t parent_num[], const void *const src_data[], void *const dest_data) { /* If source or destination datatype is not floating-point, set all return values to zero */ if ( (src_datatype != FVM_DOUBLE && src_datatype != FVM_FLOAT) || (dest_datatype != FVM_DOUBLE && dest_datatype != FVM_FLOAT)) { size_t data_size_c = (end_id - start_id) * (dest_dim * fvm_datatype_size[dest_datatype]); memset(dest_data, 0, data_size_c); } /* Otherwise, interpolate values */ else { _vertex_field_of_real_values(this_tesselation, src_dim, src_dim_shift, dest_dim, start_id, end_id, src_interlace, src_datatype, dest_datatype, n_parent_lists, parent_num_shift, parent_num, src_data, dest_data); } } /*---------------------------------------------------------------------------- * Dump printout of a mesh section tesselation structure. * * parameters: * this_tesselation <-- pointer to structure that should be dumped *----------------------------------------------------------------------------*/ void fvm_tesselation_dump(const fvm_tesselation_t *this_tesselation) { int i; fvm_lnum_t n_elements, j, k; const fvm_lnum_t *idx, *num; if (this_tesselation == NULL) return; /* Global indicators */ /*--------------------*/ bft_printf(_("\n" "Tesselation:\n\n" "Element type: %s\n" "Number of elements: %ld\n" "Spatial dimension: %d\n" "Entity dimension: %d\n"), fvm_elements_type_name[this_tesselation->type], (long)this_tesselation->n_elements, this_tesselation->dim, this_tesselation->entity_dim); bft_printf(_("\n" "Stride: %d\n" "Number of faces: %d\n"), this_tesselation->stride, (long)(this_tesselation->n_faces)); bft_printf(_("\n" "Pointers to shared arrays:\n" " vertex_coords %p\n" " parent_vertex_num %p\n" " face_num: %p\n" " vertex_index: %p\n" " vertex_num: %p\n"), this_tesselation->vertex_coords, this_tesselation->parent_vertex_num, this_tesselation->face_index, this_tesselation->face_num, this_tesselation->vertex_index, this_tesselation->vertex_num); bft_printf(_("\n" "Pointers to shared global numbering:\n" " global_element_num %p\n"), this_tesselation->global_element_num); /* Basic information */ /*-------------------*/ bft_printf(_("\n" "Number of sub-entity types: %d\n\n"), this_tesselation->n_sub_types); for (i = 0; i < this_tesselation->n_sub_types; i++) { bft_printf(_("Maximum local number of resulting %s per element: %ld\n"), fvm_elements_type_name[this_tesselation->sub_type[i]], (long)(this_tesselation->n_sub_max[i])); } for (i = 0; i < this_tesselation->n_sub_types; i++) { bft_printf(_("Maximum global number of resulting %s per element: %ld\n"), fvm_elements_type_name[this_tesselation->sub_type[i]], (long)(this_tesselation->n_sub_max_glob[i])); } bft_printf("\n"); for (i = 0; i < this_tesselation->n_sub_types; i++) { bft_printf(_("Local number of resulting %s: %ld\n"), fvm_elements_type_name[this_tesselation->sub_type[i]], (long)(this_tesselation->n_sub[i])); } for (i = 0; i < this_tesselation->n_sub_types; i++) { bft_printf(_("Global number of resulting %s: %lu\n"), fvm_elements_type_name[this_tesselation->sub_type[i]], (unsigned long)(this_tesselation->n_sub_glob[i])); } bft_printf(_("\n" "Pointers to shareable arrays:\n" " encoding: %p\n"), this_tesselation->encoding); for (i = 0; i < this_tesselation->n_sub_types; i++) { if (this_tesselation->sub_elt_index[i] != NULL) bft_printf(_(" sub_elt_index[%s]: %p\n"), fvm_elements_type_name[this_tesselation->sub_type[i]], this_tesselation->sub_elt_index[i]); } bft_printf(_("\n" "Pointers to local arrays:\n" " _encoding: %p\n"), this_tesselation->_encoding); for (i = 0; i < this_tesselation->n_sub_types; i++) { if (this_tesselation->sub_elt_index[i] != NULL) bft_printf(_(" _sub_elt_index[%s]: %p\n"), fvm_elements_type_name[this_tesselation->sub_type[i]], this_tesselation->_sub_elt_index[i]); } if (this_tesselation->encoding != NULL) { if (this_tesselation->type != FVM_FACE_QUAD) { bft_printf(_("\nEncoding (opposite vertex numbers):\n\n")); if (this_tesselation->n_faces > 0) n_elements = this_tesselation->n_faces; else n_elements = this_tesselation->n_elements; idx = this_tesselation->vertex_index; num = this_tesselation->encoding; for (j = 0; j < n_elements; j++) { bft_printf("%10d (idx = %10d) %10d\n", j+1, idx[j], num[idx[j]]); for (k = idx[j] + 1; k < idx[j + 1]; k++) bft_printf(" %10d\n", num[k]); } bft_printf(_(" end (idx = %10d)\n"), idx[n_elements]); } else { /* if (this_tesselation->type != FVM_FACE_QUAD) */ bft_printf(_("\nEncoding (diagonal flag):\n\n")); n_elements = this_tesselation->n_elements; num = this_tesselation->encoding; for (j = 0; j < n_elements; j++) bft_printf("%10d: %10d\n", j+1, num[j]); } } for (i = 0; i < this_tesselation->n_sub_types; i++) { if (this_tesselation->sub_elt_index[i] != NULL) { bft_printf(_("\nSub-element index [%s]:\n\n"), fvm_elements_type_name[this_tesselation->sub_type[i]]); n_elements = this_tesselation->n_elements; idx = this_tesselation->sub_elt_index[i]; for (j = 0; j < n_elements; j++) bft_printf("%10d: idx = %10d\n", j+1, idx[j]); bft_printf(_(" end: idx = %10d\n"), idx[n_elements]); } } } /*----------------------------------------------------------------------------*/ #ifdef __cplusplus } #endif /* __cplusplus */