/*============================================================================ * Locate points in a nodal representation associated with a mesh *============================================================================*/ /* 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) 2005-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 */ /*----------------------------------------------------------------------------*/ #include /*---------------------------------------------------------------------------- * Standard C library headers *----------------------------------------------------------------------------*/ #include #include #include #include #include #if defined(FVM_HAVE_MPI) #include #endif /*---------------------------------------------------------------------------- * 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 back Emacs auto-indentation back to column 0 */ #endif #endif /* __cplusplus */ /*============================================================================ * Local macro definitions *============================================================================*/ /* In case does not define HUGE_VAL, use a "safe" value */ #if !defined(HUGE_VAL) #define HUGE_VAL 1.0e+17 #endif /*============================================================================ * Type definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Structure defining a locator *----------------------------------------------------------------------------*/ struct _fvm_locator_t { /* Basic information */ /*-------------------*/ double tolerance; /* Associated tolerance */ int dim; /* Spatial dimension */ #if defined(FVM_HAVE_MPI) MPI_Comm comm; /* Associated MPI communicator */ #endif int n_ranks; /* Number of MPI ranks of distant location */ int start_rank; /* First MPI rank of distant location */ int n_intersects; /* Number of intersecting distant ranks */ int *intersect_rank; /* List of intersecting distant ranks */ double *intersect_extents; /* List of intersecting distant extents */ fvm_lnum_t *local_points_idx; /* Start index of local points per rank (size: n_intersects + 1)*/ fvm_lnum_t *distant_points_idx; /* Start index of distant points per rank (size: n_intersects + 1)*/ fvm_lnum_t *local_point_index; /* Local point index for data received (with blocs starting at local_points_idx[] indexes, 0 to n-1 numbering) */ fvm_lnum_t *distant_point_location; /* Location of distant points by parent element number (with blocs starting at distant_points_idx[] indexes) */ fvm_coord_t *distant_point_coords; /* Coordinates of distant points (with blocs starting at distant_points_idx[]*dim indexes) */ fvm_lnum_t n_interior; /* Number of local points located */ fvm_lnum_t *interior_list; /* List (1 to n numbering) of points located */ fvm_lnum_t n_exterior; /* Number of local points not located */ fvm_lnum_t *exterior_list; /* List (1 to n numbering) of points not located */ }; /*============================================================================ * Static global variables *============================================================================*/ /*============================================================================ * Private function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Test if two extents intersect * * parameters: * dim <-- spatial (coordinates) dimension * extents_1 <-- extents: x_min, y_min, ..., x_max, y_max, ... * size: dim*2 * extents_2 <-- extents: x_min, y_min, ..., x_max, y_max, ... * size: dim*2 * * returns: * true if extents intersect, false otherwise *----------------------------------------------------------------------------*/ inline static _Bool _intersect_extents(const int dim, const double extents_1[], const double extents_2[]) { int i; _Bool retval = true; for (i = 0 ; i < dim ; i++) { if ( (extents_1[i] > extents_2[i + dim]) || (extents_2[i] > extents_1[i + dim])) { retval = false; break; } } return retval; } /*---------------------------------------------------------------------------- * Test if a point is within given extents * * parameters: * dim <-- spatial (coordinates) dimension * coords <-- coordinates: x, y, ... * size: dim * extents <-- extents: x_min, y_min, ..., x_max, y_max, ... * size: dim*2 * * returns: * true if point lies within extents, false otherwise *----------------------------------------------------------------------------*/ inline static _Bool _within_extents(const int dim, const double coords[], const double extents[]) { int i; _Bool retval = true; for (i = 0 ; i < dim ; i++) { if ( (coords[i] < extents[i]) || (coords[i] > extents[i + dim])) { retval = false; break; } } return retval; } /*---------------------------------------------------------------------------- * Update element extents with a given vertex * * parameters: * dim <-- spatial (coordinates) dimension * vertex_id <-- vertex index (0 to n-1) * parent_vertex_num <-- pointer to parent vertex numbers (or NULL) * vertex_coords <-- pointer to vertex coordinates * elt_extents <-> extents associated with element: * x_min, y_min, ..., x_max, y_max, ... (size: 2*dim) * elt_initialized <-> are extents already initialized for this vertex * (for all element vertices except the first) ? *----------------------------------------------------------------------------*/ inline static void _update_elt_extents(const int dim, const fvm_lnum_t vertex_id, const fvm_lnum_t *const parent_vertex_num, const fvm_coord_t vertex_coords[], double elt_extents[], _Bool *const elt_initialized) { fvm_lnum_t i, coord_idx; if (parent_vertex_num == NULL) coord_idx = vertex_id; else coord_idx = parent_vertex_num[vertex_id] - 1; if (*elt_initialized == false) { for (i = 0 ; i < dim ; i++) { elt_extents[i] = vertex_coords[(coord_idx * dim) + i]; elt_extents[i + dim] = vertex_coords[(coord_idx * dim) + i]; } *elt_initialized = true; } else { for (i = 0 ; i < dim ; i++) { if (elt_extents[i] > vertex_coords[(coord_idx * dim) + i]) elt_extents[i] = vertex_coords[(coord_idx * dim) + i]; if (elt_extents[i + dim] < vertex_coords[(coord_idx * dim) + i]) elt_extents[i + dim] = vertex_coords[(coord_idx * dim) + i]; } } } /*---------------------------------------------------------------------------- * Adjust element extents with search tolerance and update global extents * * parameters: * dim <-- spatial (coordinates) dimension * elt_dim <-- element dimension * tolerance <-- addition to local extents of each element: * extent = base_extent * (1 + tolerance) * elt_extents <-> extents associated with element: * x_min, y_min, ..., x_max, y_max, ... (size: 2*dim) * extents <-> mesh or section extents, to be updated: * x_min, y_min, z_min, x_max, y_max, z_max *----------------------------------------------------------------------------*/ inline static void _elt_extents_finalize(const int dim, const int elt_dim, const double tolerance, double *const restrict elt_extents, double *const restrict extents) { int i; double delta[3]; for (i = 0 ; i < dim ; i++) delta[i] = (elt_extents[i+dim] - elt_extents[i]) * tolerance; if (elt_dim < 3) { double delta_max = delta[0]; /* for 1d or 2d elements, ensure */ for (i = 0 ; i < dim ; i++) { /* search extent "thickness" */ if (delta[i] > delta_max) delta_max = delta[i]; } for (i = 0 ; i < dim ; i++) delta[i] = delta_max; } for (i = 0 ; i < dim ; i++) { elt_extents[i] = elt_extents[i] - delta[i]; elt_extents[i+dim] = elt_extents[i+dim] + delta[i]; } for (i = 0 ; i < dim ; i++) { if (elt_extents[i] < extents[i]) extents[i] = elt_extents[i]; if (elt_extents[i+dim] > extents[i+dim]) extents[i+dim] = elt_extents[i+dim]; } } /*---------------------------------------------------------------------------- * Compute extents of a nodal mesh representation section * * parameters: * this_section <-- pointer to section structure * dim <-- spatial (coordinates) dimension * parent_vertex_num <-- pointer to parent vertex numbers (or NULL) * vertex_coords <-- pointer to vertex coordinates * tolerance <-- addition to local extents of each element: * extent = base_extent * (1 + tolerance) * extents <-> extents associated with section: * x_min, y_min, ..., x_max, y_max, ... (size: 2*dim) *----------------------------------------------------------------------------*/ static void _nodal_section_extents(const fvm_nodal_section_t *const this_section, const int dim, const fvm_lnum_t *const parent_vertex_num, const fvm_coord_t vertex_coords[], const double tolerance, double extents[]) { fvm_lnum_t i, j, k, face_id, vertex_id; double elt_extents[6]; /* initialize extents in case section is empty */ for (i = 0 ; i < dim ; i++) { extents[i] = HUGE_VAL; extents[i + dim] = -HUGE_VAL; } /* Extents for polyhedra */ if (this_section->face_index != NULL) { for (i = 0 ; i < this_section->n_elements ; i++) { _Bool elt_initialized = false; for (j = this_section->face_index[i] ; j < this_section->face_index[i + 1] ; j++) { face_id = FVM_ABS(this_section->face_num[j]) - 1; for (k = this_section->vertex_index[face_id] ; k < this_section->vertex_index[face_id + 1] ; k++) { vertex_id = this_section->vertex_num[k] - 1; _update_elt_extents(dim, vertex_id, parent_vertex_num, vertex_coords, elt_extents, &elt_initialized); } } _elt_extents_finalize(dim, 3, tolerance, elt_extents, extents); } } /* Extents for polygons */ else if (this_section->vertex_index != NULL) { fvm_lnum_t n_faces = (this_section->n_faces > 0) ? this_section->n_faces : this_section->n_elements; for (i = 0 ; i < n_faces ; i++) { _Bool elt_initialized = false; for (j = this_section->vertex_index[i] ; j < this_section->vertex_index[i + 1] ; j++) { vertex_id = this_section->vertex_num[j] - 1; _update_elt_extents(dim, vertex_id, parent_vertex_num, vertex_coords, elt_extents, &elt_initialized); } _elt_extents_finalize(dim, 2, tolerance, elt_extents, extents); } } /* Extents for regular elements */ else { for (i = 0 ; i < this_section->n_elements ; i++) { _Bool elt_initialized = false; for (j = 0 ; j < this_section->stride ; j++) { vertex_id = this_section->vertex_num[i*this_section->stride + j] - 1; _update_elt_extents(dim, vertex_id, parent_vertex_num, vertex_coords, elt_extents, &elt_initialized); } _elt_extents_finalize(dim, this_section->entity_dim, tolerance, elt_extents, extents); } } } /*---------------------------------------------------------------------------- * Compute extents of a nodal mesh representation * * parameters: * this_nodal <-- pointer to mesh representation structure * tolerance <-- addition to local extents of each element: * extent = base_extent * (1 + tolerance) * extents <-> extents associated with mesh: * x_min, y_min, ..., x_max, y_max, ... (size: 2*dim) *----------------------------------------------------------------------------*/ static void _nodal_extents(const fvm_nodal_t *const this_nodal, const double tolerance, double extents[]) { int i, j; int dim; double section_extents[6]; if (this_nodal == NULL) return; dim = this_nodal->dim; /* initialize extents in case mesh is empty or dim < 3 */ for (i = 0 ; i < dim ; i++) { extents[i] = HUGE_VAL; extents[i + dim] = -HUGE_VAL; } /* Compute extents */ for (i = 0 ; i < this_nodal->n_sections ; i++) { _nodal_section_extents(this_nodal->sections[i], this_nodal->dim, this_nodal->parent_vertex_num, this_nodal->vertex_coords, tolerance, section_extents); for (j = 0 ; j < this_nodal->dim ; j++) { if (section_extents[j] < extents[j]) extents[j] = section_extents[j]; if (section_extents[j+dim] > extents[j+dim]) extents[j+dim] = section_extents[j+dim]; } } } /*---------------------------------------------------------------------------- * Compute extents of a point set * * parameters: * dim <-- space dimension of points to locate * n_points <-- number of points to locate * point_index <-- optional indirection array to point_coords * (1 to n_points numbering) * point_coords <-- coordinates of points to locate * (dimension: dim * n_points) * extents --> extents associated with mesh: * x_min, y_min, ..., x_max, y_max, ... (size: 2*dim) *----------------------------------------------------------------------------*/ static void _point_extents(const int dim, const fvm_lnum_t n_points, const fvm_lnum_t point_index[], const fvm_coord_t point_coords[], double extents[]) { int i; fvm_lnum_t j, coord_idx; /* initialize extents in case mesh is empty or dim < 3 */ for (i = 0 ; i < dim ; i++) { extents[i] = HUGE_VAL; extents[i + dim] = -HUGE_VAL; } /* Compute extents */ if (point_index != NULL) { for (j = 0 ; j < n_points ; j++) { coord_idx = point_index[j] - 1; for (i = 0 ; i < dim ; i++) { if (extents[i] > point_coords[(coord_idx * dim) + i]) extents[i] = point_coords[(coord_idx * dim) + i]; if (extents[i + dim] < point_coords[(coord_idx * dim) + i]) extents[i + dim] = point_coords[(coord_idx * dim) + i]; } } } else { for (coord_idx = 0 ; coord_idx < n_points ; coord_idx++) { for (i = 0 ; i < dim ; i++) { if (extents[i] > point_coords[(coord_idx * dim) + i]) extents[i] = point_coords[(coord_idx * dim) + i]; if (extents[i + dim] < point_coords[(coord_idx * dim) + i]) extents[i + dim] = point_coords[(coord_idx * dim) + i]; } } } } /*---------------------------------------------------------------------------- * Updates the location[] and distance[] arrays associated with a set * of points for points that are in a given element extent, based only * on this extent only (temporary, unoptimzed location). * * parameters: * elt_num <-- number of element corresponding to extents * extents <-> extents associated with element: * x_min, y_min, ..., x_max, y_max, ... (size: 2*dim) * dim <-- spatial dimension * n_points <-- number of points to locate * point_coords <-- point coordinates * location <-> number of element containing or closest to each * point (size: n_points) * distance <-> distance from point to element indicated by * location[]: < 0 if unlocated, 0 - 1 if inside, * > 1 if outside (size: n_points) *----------------------------------------------------------------------------*/ static void _locate_by_extents(const fvm_lnum_t elt_num, const int dim, const double extents[], const fvm_lnum_t n_points, const fvm_coord_t point_coords[], fvm_lnum_t location[], float distance[]) { fvm_lnum_t i, j; /* For now, we base a minimal location test on the element extents */ /* The behavior is quadradic, nothing is optimized yet */ for (i = 0 ; i < n_points ; i++) { double elt_coord_max = -1; double elt_coord = -1; for (j = 0 ; j < dim ; j++) { double cur_coord = point_coords[i*dim + j]; elt_coord = (cur_coord - 0.5*(extents[j+dim] + extents[j])) / ( 0.5*(extents[j+dim] - extents[j])); elt_coord = FVM_ABS(elt_coord); if (elt_coord > elt_coord_max) elt_coord_max = elt_coord; } if ( (distance[i] < 0 && elt_coord_max < 1) || elt_coord_max < distance[i]) { location[i] = elt_num; distance[i] = elt_coord_max; } } } /*---------------------------------------------------------------------------- * Find elements in a given section containing points: updates the * location[] and distance[] arrays associated with a set of points * for points that are in an element of this section, or closer to one * than to previously encountered elements. * * parameters: * this_section <-- pointer to mesh section representation structure * dim <-- spatial dimension * parent_vertex_num <-- pointer to parent vertex numbers (or NULL) * vertex_coords <-- pointer to vertex coordinates * tolerance <-- associated tolerance * n_points <-- number of points to locate * point_coords <-- point coordinates * location <-> number of element containing or closest to each * point (size: n_points) * distance <-> distance from point to element indicated by * location[]: < 0 if unlocated, 0 - 1 if inside, * > 1 if outside (size: n_points) *----------------------------------------------------------------------------*/ static void _nodal_section_locate(const fvm_nodal_section_t *const this_section, const int dim, const fvm_lnum_t *const parent_vertex_num, const fvm_coord_t vertex_coords[], const double tolerance, const fvm_lnum_t n_points, const fvm_coord_t point_coords[], fvm_lnum_t location[], float distance[]) { fvm_lnum_t i, j, k, face_id, vertex_id, elt_num; double elt_extents[6], section_extents[6]; /* Initialize section extents (argument to _update_elt_extents(), not used here) */ for (i = 0 ; i < dim ; i++) { section_extents[i] = HUGE_VAL; section_extents[i + dim] = -HUGE_VAL; } /* If section contains polyhedra */ if (this_section->face_index != NULL) { for (i = 0 ; i < this_section->n_elements ; i++) { _Bool elt_initialized = false; /* Compute extents */ for (j = this_section->face_index[i] ; j < this_section->face_index[i + 1] ; j++) { face_id = FVM_ABS(this_section->face_num[j]) - 1; for (k = this_section->vertex_index[face_id] ; k < this_section->vertex_index[face_id + 1] ; k++) { vertex_id = this_section->vertex_num[k] - 1; _update_elt_extents(dim, vertex_id, parent_vertex_num, vertex_coords, elt_extents, &elt_initialized); } } _elt_extents_finalize(dim, 3, tolerance, elt_extents, section_extents); if (this_section->parent_element_num != NULL) elt_num = this_section->parent_element_num[i]; else elt_num = i + 1; _locate_by_extents(elt_num, dim, elt_extents, n_points, point_coords, location, distance); } } /* Extents for polygons */ else if (this_section->vertex_index != NULL) { fvm_lnum_t n_faces = (this_section->n_faces > 0) ? this_section->n_faces : this_section->n_elements; for (i = 0 ; i < n_faces ; i++) { _Bool elt_initialized = false; for (j = this_section->vertex_index[i] ; j < this_section->vertex_index[i + 1] ; j++) { vertex_id = this_section->vertex_num[j] - 1; _update_elt_extents(dim, vertex_id, parent_vertex_num, vertex_coords, elt_extents, &elt_initialized); } _elt_extents_finalize(dim, 2, tolerance, elt_extents, section_extents); if (this_section->parent_element_num != NULL) elt_num = this_section->parent_element_num[i]; else elt_num = i + 1; _locate_by_extents(elt_num, dim, elt_extents, n_points, point_coords, location, distance); } } /* Extents for regular elements */ else { for (i = 0 ; i < this_section->n_elements ; i++) { _Bool elt_initialized = false; for (j = 0 ; j < this_section->stride ; j++) { vertex_id = this_section->vertex_num[i*this_section->stride + j] - 1; _update_elt_extents(dim, vertex_id, parent_vertex_num, vertex_coords, elt_extents, &elt_initialized); } _elt_extents_finalize(dim, this_section->entity_dim, tolerance, elt_extents, section_extents); if (this_section->parent_element_num != NULL) elt_num = this_section->parent_element_num[i]; else elt_num = i + 1; _locate_by_extents(elt_num, dim, elt_extents, n_points, point_coords, location, distance); } } } /*---------------------------------------------------------------------------- * Find elements in a given nodal mesh containing points: updates the * location[] and distance[] arrays associated with a set of points * for points that are in an element of this mesh, or closer to one * than to previously encountered elements. * * parameters: * this_nodal <-- pointer to nodal mesh representation structure * tolerance <-- associated tolerance * n_points <-- number of points to locate * point_coords <-- point coordinates * location <-> number of element containing or closest to each * point (size: n_points) * distance <-> distance from point to element indicated by * location[]: < 0 if unlocated, 0 - 1 if inside, * > 1 if outside (size: n_points) *----------------------------------------------------------------------------*/ static void _nodal_locate(const fvm_nodal_t *const this_nodal, const double tolerance, const fvm_lnum_t n_points, const fvm_coord_t point_coords[], fvm_lnum_t location[], float distance[]) { int i; if (this_nodal == NULL) return; /* Compute extents */ for (i = 0 ; i < this_nodal->n_sections ; i++) { const fvm_nodal_section_t *this_section = this_nodal->sections[i]; if (this_section->entity_dim == this_nodal->dim) _nodal_section_locate(this_section, this_nodal->dim, this_nodal->parent_vertex_num, this_nodal->vertex_coords, tolerance, n_points, point_coords, location, distance); } } #if defined(FVM_HAVE_MPI) /*---------------------------------------------------------------------------- * Prepare locator for use with a given nodal mesh representation and * point set. * * parameters: * this_locator <-> pointer to locator structure * this_nodal <-- pointer to mesh representation structure * dim <-- spatial dimension * n_points <-- number of points to locate * point_index <-- optional indirection array to point_coords * (1 to n_points numbering) * point_coords <-- coordinates of points to locate * (dimension: dim * n_points) *----------------------------------------------------------------------------*/ static void _locate_all_distant(fvm_locator_t *const this_locator, const fvm_nodal_t *const this_nodal, const int dim, const fvm_lnum_t n_points, const fvm_lnum_t point_index[], const fvm_coord_t point_coords[]) { int i, k, stride; int dist_rank, dist_index; fvm_lnum_t j; fvm_lnum_t n_coords_loc, n_coords_dist, n_interior, n_exterior; fvm_lnum_t coord_idx, start_idx; fvm_lnum_t *location_loc, *location_dist; fvm_lnum_t *location, *location_rank_id; fvm_lnum_t *send_index, *send_location; fvm_lnum_t *location_count, *location_shift; fvm_coord_t *coords_dist, *send_coords; float *distance; float *distance_dist, *distance_loc; double extents[6]; MPI_Status status; /* Initialization */ stride = dim * 2; BFT_MALLOC(send_coords, n_points * dim, fvm_coord_t); BFT_MALLOC(send_index, n_points, fvm_lnum_t); BFT_MALLOC(location, n_points, fvm_lnum_t); BFT_MALLOC(location_rank_id, n_points, fvm_lnum_t); BFT_MALLOC(distance, n_points, float); for (j = 0 ; j < n_points ; j++) { location[j] = -1; location_rank_id[j] = -1; distance[j] = -1.0; } /* First loop on possibly intersecting distant ranks */ /*---------------------------------------------------*/ for (i = 0 ; i < this_locator->n_intersects ; i++) { dist_index = i; /* Ordering (communication schema) not yet optimized */ dist_rank = this_locator->intersect_rank[dist_index]; /* Prepare and send coords that should fit in each send buffer */ /* Reset buffers for current intersect rank */ n_coords_loc = 0; for (k = 0 ; k < stride ; k++) extents[k] = this_locator->intersect_extents[dist_index*stride + k]; /* Build partial buffer */ for (j = 0 ; j < n_points ; j++) { if (point_index != NULL) coord_idx = point_index[j] - 1; else coord_idx = j; if (_within_extents(dim, &(point_coords[dim*coord_idx]), extents) == true) { send_index[n_coords_loc] = j; for (k = 0 ; k < dim ; k++) send_coords[n_coords_loc*dim + k] = point_coords[dim*coord_idx + k]; n_coords_loc += 1; } } /* Send then receive partial buffer */ dist_rank = this_locator->intersect_rank[dist_index]; MPI_Sendrecv(&n_coords_loc, 1, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, &n_coords_dist, 1, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); BFT_MALLOC(coords_dist, n_coords_dist*dim, fvm_coord_t); MPI_Sendrecv(send_coords, (int)(n_coords_loc*dim), FVM_MPI_COORD, dist_rank, FVM_MPI_TAG, coords_dist, (int)(n_coords_dist*dim), FVM_MPI_COORD, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); /* Now locate received coords on local rank */ BFT_MALLOC(location_dist, n_coords_dist, fvm_lnum_t); BFT_MALLOC(distance_dist, n_coords_dist, float); for (j = 0 ; j < n_coords_dist ; j++) { location_dist[j] = -1; distance_dist[j] = -1.0; } _nodal_locate(this_nodal, this_locator->tolerance, n_coords_dist, coords_dist, location_dist, distance_dist); BFT_FREE(coords_dist); /* Exchange location return information with distant rank */ BFT_MALLOC(location_loc, n_coords_loc, fvm_lnum_t); BFT_MALLOC(distance_loc, n_coords_loc, float); MPI_Sendrecv(location_dist, (int)n_coords_dist, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, location_loc, (int)n_coords_loc, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); MPI_Sendrecv(distance_dist, (int)n_coords_dist, MPI_FLOAT, dist_rank, FVM_MPI_TAG, distance_loc, (int)n_coords_loc, MPI_FLOAT, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); BFT_FREE(location_dist); BFT_FREE(distance_dist); /* Now update location information */ for (j = 0 ; j < n_coords_loc ; j++) { fvm_lnum_t l = send_index[j]; if ( (distance_loc[j] > -0.1) && (distance_loc[j] < distance[l] || distance[l] < -0.1)) { location_rank_id[l] = i; location[l] = location_loc[j]; distance[l] = distance_loc[j]; } } BFT_FREE(location_loc); BFT_FREE(distance_loc); } /* Reorganize localization information */ /*-------------------------------------*/ /* Now that localization is done, the location[] array contains either -1 if a point was not localized, or a local index (associated with the corresponding rank); the distance[] array is not needed anymore now that all comparisons have been done */ BFT_FREE(distance); BFT_MALLOC(location_shift, this_locator->n_intersects, fvm_lnum_t); BFT_MALLOC(location_count, this_locator->n_intersects, fvm_lnum_t); for (i = 0 ; i < this_locator->n_intersects ; i++) location_count[i] = 0; n_exterior = 0; for (j = 0 ; j < n_points ; j++) { if (location_rank_id[j] > -1) location_count[location_rank_id[j]] += 1; else n_exterior += 1; } this_locator->n_interior = n_points - n_exterior; BFT_MALLOC(this_locator->interior_list, this_locator->n_interior, fvm_lnum_t); this_locator->n_exterior = n_exterior; BFT_MALLOC(this_locator->exterior_list, this_locator->n_exterior, fvm_lnum_t); if (this_locator->n_intersects > 0) location_shift[0] = 0; for (i = 1 ; i < this_locator->n_intersects ; i++) location_shift[i] = location_shift[i-1] + location_count[i-1]; for (i = 0 ; i < this_locator->n_intersects ; i++) location_count[i] = 0; /* send_index[] will now contain information for all blocks */ for (j = 0 ; j < n_points ; j++) send_index[j] = -1; BFT_MALLOC(send_location, n_points, fvm_lnum_t); n_interior = 0; n_exterior = 0; for (j = 0 ; j < n_points ; j++) { const int l_rank = location_rank_id[j]; if (l_rank > -1) { send_index[location_shift[l_rank] + location_count[l_rank]] = j; location_count[l_rank] += 1; this_locator->interior_list[n_interior] = j + 1; n_interior += 1; } else { this_locator->exterior_list[n_exterior] = j + 1; n_exterior += 1; } } /* Second loop on possibly intersecting distant ranks */ /*----------------------------------------------------*/ /* Count and organize total number of local and distant points */ BFT_MALLOC(this_locator->local_points_idx, this_locator->n_intersects + 1, fvm_lnum_t); BFT_MALLOC(this_locator->distant_points_idx, this_locator->n_intersects + 1, fvm_lnum_t); this_locator->local_points_idx[0] = 0; this_locator->distant_points_idx[0] = 0; for (i = 0 ; i < this_locator->n_intersects ; i++) { dist_index = i; /* Ordering (communication schema) not yet optimized */ dist_rank = this_locator->intersect_rank[dist_index]; n_coords_loc = location_count[i]; this_locator->local_points_idx[i+1] = this_locator->local_points_idx[i] + n_coords_loc; MPI_Sendrecv(&n_coords_loc, 1, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, &n_coords_dist, 1, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); this_locator->distant_points_idx[i+1] = this_locator->distant_points_idx[i] + n_coords_dist; } /* Third loop on possibly intersecting distant ranks */ /*----------------------------------------------------*/ BFT_MALLOC(this_locator->local_point_index, this_locator->local_points_idx[this_locator->n_intersects], fvm_lnum_t); BFT_MALLOC(this_locator->distant_point_location, this_locator->distant_points_idx[this_locator->n_intersects], fvm_lnum_t); BFT_MALLOC(this_locator->distant_point_coords, this_locator->distant_points_idx[this_locator->n_intersects] * dim, fvm_coord_t); for (i = 0 ; i < this_locator->n_intersects ; i++) { dist_index = i; /* Ordering (communication schema) not yet optimized */ dist_rank = this_locator->intersect_rank[dist_index]; n_coords_loc = this_locator->local_points_idx[i+1] - this_locator->local_points_idx[i]; n_coords_dist = this_locator->distant_points_idx[i+1] - this_locator->distant_points_idx[i]; start_idx = this_locator->local_points_idx[i]; for (j = 0 ; j < n_coords_loc ; j++) { coord_idx = send_index[location_shift[i] + j]; this_locator->local_point_index[start_idx + j] = coord_idx; send_location[j] = location[coord_idx]; if (point_index != NULL) { for (k = 0 ; k < dim ; k++) send_coords[j*dim + k] = point_coords[dim*(point_index[coord_idx] - 1) + k]; } else { for (k = 0 ; k < dim ; k++) send_coords[j*dim + k] = point_coords[dim*coord_idx + k]; } } MPI_Sendrecv(send_location, (int)n_coords_loc, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, (this_locator->distant_point_location + this_locator->distant_points_idx[i]), (int)n_coords_dist, FVM_MPI_LNUM, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); MPI_Sendrecv(send_coords, (int)(n_coords_loc*dim), FVM_MPI_COORD, dist_rank, FVM_MPI_TAG, (this_locator->distant_point_coords + (this_locator->distant_points_idx[i]*dim)), (int)(n_coords_dist*dim), FVM_MPI_COORD, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); } BFT_FREE(location_count); BFT_FREE(location_shift); BFT_FREE(send_index); BFT_FREE(send_location); BFT_FREE(send_coords); BFT_FREE(location_rank_id); BFT_FREE(location); } #endif /* defined(FVM_HAVE_MPI) */ /*---------------------------------------------------------------------------- * Prepare locator for use with a given nodal mesh representation and * point set. * * parameters: * this_locator <-> pointer to locator structure * this_nodal <-- pointer to mesh representation structure * dim <-- spatial dimension * n_points <-- number of points to locate * point_index <-- optional indirection array to point_coords * (1 to n_points numbering) * point_coords <-- coordinates of points to locate * (dimension: dim * n_points) *----------------------------------------------------------------------------*/ static void _locate_all_local(fvm_locator_t *const this_locator, const fvm_nodal_t *const this_nodal, const int dim, const fvm_lnum_t n_points, const fvm_lnum_t point_index[], const fvm_coord_t point_coords[]) { int k, stride; fvm_lnum_t j; fvm_lnum_t n_coords, n_interior, n_exterior, coord_idx; fvm_lnum_t *location; fvm_lnum_t location_count; fvm_coord_t *coords; float *distance; double extents[6]; /* Initialization */ stride = this_nodal->dim * 2; BFT_MALLOC(coords, n_points * this_nodal->dim, fvm_coord_t); BFT_MALLOC(location, n_points, fvm_lnum_t); BFT_MALLOC(distance, n_points, float); for (j = 0 ; j < n_points ; j++) { location[j] = -1; distance[j] = -1.0; } /* Initialize localization information */ /*-------------------------------------*/ n_coords = 0; for (k = 0 ; k < stride ; k++) extents[k] = this_locator->intersect_extents[k]; /* Build partial buffer */ for (j = 0 ; j < n_points ; j++) { if (point_index != NULL) coord_idx = point_index[j] - 1; else coord_idx = j; if (_within_extents(dim, &(point_coords[dim*coord_idx]), extents) == true) { for (k = 0 ; k < dim ; k++) coords[n_coords*dim + k] = point_coords[dim*coord_idx + k]; n_coords += 1; } } /* Now locate coords */ BFT_MALLOC(location, n_coords, fvm_lnum_t); BFT_MALLOC(distance, n_coords, float); for (j = 0 ; j < n_coords ; j++) { location[j] = -1; distance[j] = -1.0; } _nodal_locate(this_nodal, this_locator->tolerance, n_coords, coords, location, distance); /* Reorganize localization information */ /*-------------------------------------*/ /* Now that localization is done, the location[] array contains either -1 if a point was not localized, or a local index; the distance[] array is not needed anymore now that all comparisons have been done */ BFT_FREE(distance); location_count = 0; n_exterior = 0; for (j = 0 ; j < n_points ; j++) { if (location[j] > -1) location_count += 1; else n_exterior += 1; } this_locator->n_interior = n_points - n_exterior; BFT_MALLOC(this_locator->interior_list, this_locator->n_interior, fvm_lnum_t); this_locator->n_exterior = n_exterior; BFT_MALLOC(this_locator->exterior_list, this_locator->n_exterior, fvm_lnum_t); /* Organize total number of "local" and "distant" points */ BFT_MALLOC(this_locator->local_points_idx, 2, fvm_lnum_t); BFT_MALLOC(this_locator->distant_points_idx, 2, fvm_lnum_t); this_locator->local_points_idx[0] = 0; this_locator->local_points_idx[1] = location_count; this_locator->distant_points_idx[0] = 0; this_locator->distant_points_idx[1] = n_coords; assert(location_count == n_coords); /* If this is false, review _exchange_point_var_local() */ BFT_MALLOC(this_locator->local_point_index, location_count, fvm_lnum_t); BFT_MALLOC(this_locator->distant_point_location, location_count, fvm_lnum_t); BFT_MALLOC(this_locator->distant_point_coords, n_coords * dim, fvm_coord_t); location_count = 0; n_interior = 0; n_exterior = 0; for (j = 0 ; j < n_points ; j++) { if (location[j] > -1) { this_locator->local_point_index[location_count] = j; this_locator->distant_point_location[location_count] = location[j]; for (k = 0 ; k < dim ; k++) { this_locator->distant_point_coords[location_count*dim + k] = point_coords[j*dim + k]; } location_count += 1; this_locator->interior_list[n_interior] = j + 1; n_interior += 1; } else this_locator->exterior_list[n_exterior] = j + 1; n_exterior += 1; } BFT_FREE(location); BFT_FREE(coords); } #if defined(FVM_HAVE_MPI) /*---------------------------------------------------------------------------- * Distribute variable defined on distant points to processes owning * the original points (i.e. distant processes). * * The exchange is symmetric if both variables are defined, receive * only if distant_var is NULL, or send only if local_var is NULL. * * parameters: * this_locator <-- pointer to locator structure * distant_var <-- variable defined on distant points (ready to send) * local_var --> variable defined on local points (received) * datatype <-- variable type * stride <-- dimension (1 for scalar, 3 for interlaced vector) *----------------------------------------------------------------------------*/ static void _exchange_point_var_distant(fvm_locator_t *const this_locator, void *const distant_var, void *const local_var, const MPI_Datatype datatype, const size_t stride) { int i, send_count, recv_count, size; int dist_rank, dist_index; int send_flag, recv_flag; fvm_lnum_t n_points_loc, n_points_loc_max, n_points_dist; size_t send_idx; void *send_ptr; void *recvbuf; MPI_Aint extent; MPI_Status status; /* Check extent of datatype */ MPI_Type_extent(datatype, &extent); MPI_Type_size(datatype, &size); if (extent != size) bft_error(__FILE__, __LINE__, 0, _("_exchange_point_var() is not implemented for use with\n" "MPI datatypes associated with structures using padding\n" "(for which size != extent).")); /* Initialization */ n_points_loc_max = 0; for (i = 0 ; i < this_locator->n_intersects ; i++) { n_points_loc = this_locator->local_points_idx[i+1] - this_locator->local_points_idx[i]; if (n_points_loc > n_points_loc_max) n_points_loc_max = n_points_loc; } BFT_MALLOC(recvbuf, n_points_loc_max*size*stride, char); /* Loop on possibly intersecting distant ranks */ /*---------------------------------------------*/ for (i = 0 ; i < this_locator->n_intersects ; i++) { dist_index = i; /* Ordering (communication schema) not yet optimized */ dist_rank = this_locator->intersect_rank[dist_index]; n_points_loc = this_locator->local_points_idx[i+1] - this_locator->local_points_idx[i]; n_points_dist = this_locator->distant_points_idx[i+1] - this_locator->distant_points_idx[i]; if (distant_var != NULL && n_points_dist > 0) send_flag = 1; else send_flag = 0; MPI_Sendrecv(&send_flag, 1, MPI_INT, dist_rank, FVM_MPI_TAG, &recv_flag, 1, MPI_INT, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); if ( (recv_flag == 1 && (local_var == NULL || n_points_loc == 0)) || (recv_flag == 0 && n_points_loc > 0)) bft_error(__FILE__, __LINE__, 0, _("Incoherent arguments to different instances in " "_exchange_point_var().\n" "Send and receive operations do not match " "(dist_rank = %d\n)\n"), dist_rank); send_idx = this_locator->distant_points_idx[i] * stride*size; send_count = n_points_dist * stride * send_flag; if (recv_flag > 0) recv_count = n_points_loc*stride; else recv_count = 0; /* Exchange information */ if (distant_var != NULL) send_ptr = (void *)(((char *)distant_var) + send_idx); else send_ptr = NULL; MPI_Sendrecv(send_ptr, send_count, datatype, dist_rank, FVM_MPI_TAG, recvbuf, recv_count, datatype, dist_rank, FVM_MPI_TAG, this_locator->comm, &status); { int k; size_t l; const size_t nbytes = stride*size; for (k = 0; k < n_points_loc; k++) { char *local_var_p = (char *)local_var + ( this_locator->local_point_index + this_locator->local_points_idx[i])[k]*nbytes; const char *recvbuf_p = (const char *)recvbuf + k*nbytes; for (l = 0; l < nbytes; l++) local_var_p[l] = recvbuf_p[l]; } } } /* End of loop on possibly intersecting ranks */ BFT_FREE(recvbuf); } #endif /* defined(FVM_HAVE_MPI) */ /*---------------------------------------------------------------------------- * Distribute variable defined on "distant points" to the original ("local") * points. * * The exchange is symmetric if both variables are defined, receive * only if distant_var is NULL, or send only if local_var is NULL. * * parameters: * this_locator <-- pointer to locator structure * distant_var <-- variable defined on distant points (ready to send) * local_var --> variable defined on local points (received) * type_size <-- sizeof (float or double) variable type * stride <-- dimension (1 for scalar, 3 for interlaced vector) *----------------------------------------------------------------------------*/ static void _exchange_point_var_local(fvm_locator_t *const this_locator, void *const distant_var, void *const local_var, const size_t type_size, const size_t stride) { fvm_lnum_t i; size_t j; fvm_lnum_t n_points_loc, n_points_dist; const size_t nbytes = stride*type_size; /* Initialization */ n_points_loc = this_locator->local_points_idx[1] - this_locator->local_points_idx[0]; n_points_dist = this_locator->distant_points_idx[1] - this_locator->distant_points_idx[0]; assert(n_points_loc == n_points_dist); /* Exchange information */ for (i = 0; i < n_points_loc; i++) { char *local_var_p = (char *)local_var + (this_locator->local_point_index[i])*nbytes; const char *distant_var_p = (const char *)distant_var + i*nbytes; for (j = 0; j < nbytes; j++) local_var_p[j] = distant_var_p[j]; } } /*============================================================================ * Public function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Creation of a locator structure. * * Note that depending on the choice of ranks of the associated communicator, * distant ranks may in fact be truly distant or not. If n_ranks = 1 and * start_rank is equal to the current rank in the communicator, the locator * will work only locally. * * parameters: * tolerance <-- addition to local extents of each element: * extent = base_extent * (1 + tolerance) * comm <-- associated MPI communicator * n_ranks <-- number of MPI ranks associated with distant location * start_rank <-- first MPI rank associated with distant location * * returns: * pointer to locator *----------------------------------------------------------------------------*/ #if defined(FVM_HAVE_MPI) fvm_locator_t * fvm_locator_create(const double tolerance, const MPI_Comm comm, const int n_ranks, const int start_rank) #else fvm_locator_t * fvm_locator_create(const double tolerance) #endif { fvm_locator_t *this_locator; BFT_MALLOC(this_locator, 1, fvm_locator_t); this_locator->tolerance = tolerance; this_locator->dim = 0; #if defined(FVM_HAVE_MPI) this_locator->comm = comm; this_locator->n_ranks = n_ranks; this_locator->start_rank = start_rank; #else this_locator->n_ranks = 1; this_locator->start_rank = 0; #endif this_locator->n_intersects = 0 ; this_locator->intersect_rank = NULL; this_locator->intersect_extents = NULL; this_locator->local_points_idx = NULL; this_locator->distant_points_idx = NULL; this_locator->local_point_index = NULL; this_locator->distant_point_location = NULL; this_locator->distant_point_coords = NULL; this_locator->n_interior = 0; this_locator->interior_list = NULL; this_locator->n_exterior = 0; this_locator->exterior_list = NULL; return this_locator; } /*---------------------------------------------------------------------------- * Destruction of a locator structure. * * parameters: * this_locator --> locator to destroy * * returns: * NULL pointer *----------------------------------------------------------------------------*/ fvm_locator_t * fvm_locator_destroy(fvm_locator_t * this_locator) { if (this_locator != NULL) { BFT_FREE(this_locator->local_points_idx); BFT_FREE(this_locator->distant_points_idx); BFT_FREE(this_locator->local_point_index); BFT_FREE(this_locator->distant_point_location); BFT_FREE(this_locator->distant_point_coords); BFT_FREE(this_locator->intersect_rank); BFT_FREE(this_locator->intersect_extents); BFT_FREE(this_locator->interior_list); BFT_FREE(this_locator->exterior_list); BFT_FREE(this_locator); } return NULL; } /*---------------------------------------------------------------------------- * Prepare locator for use with a given nodal mesh representation. * * parameters: * this_locator <-> pointer to locator structure * this_nodal <-- pointer to (containing) mesh representation structure * dim <-- space dimension of points to locate * n_points <-- number of points to locate * point_index <-- optional indirection array to point_coords * (1 to n_points numbering) * point_coords <-- coordinates of points to locate * (dimension: dim * n_points) *----------------------------------------------------------------------------*/ void fvm_locator_set_nodal(fvm_locator_t *const this_locator, const fvm_nodal_t *const this_nodal, const int dim, const fvm_lnum_t n_points, const fvm_lnum_t point_index[], const fvm_coord_t point_coords[]) { int i; int stride2; double tolerance; double extents[12]; #if defined(FVM_HAVE_MPI) int j; int stride4; int comm_rank, comm_size; int locdim[2], globdim[2]; int n_intersects; int *intersect_rank; double *recvbuf; #endif int mpi_flag = 0; /* Compute extents */ this_locator->dim = dim; tolerance = FVM_MAX(this_locator->tolerance, 0.1); _nodal_extents(this_nodal, tolerance, extents); _point_extents(dim, n_points, point_index, point_coords, extents + 2*dim); /* Release information if previously present */ if (this_locator->intersect_rank != NULL) BFT_FREE(this_locator->intersect_rank); if (this_locator->intersect_extents != NULL) BFT_FREE(this_locator->intersect_extents); if (this_locator->local_points_idx != NULL) BFT_FREE(this_locator->local_points_idx); if (this_locator->distant_points_idx != NULL) BFT_FREE(this_locator->distant_points_idx); if (this_locator->local_point_index != NULL) BFT_FREE(this_locator->local_point_index); if (this_locator->distant_point_location != NULL) BFT_FREE(this_locator->distant_point_location); if (this_locator->distant_point_coords != NULL) BFT_FREE(this_locator->distant_point_coords); if (this_locator->interior_list != NULL) BFT_FREE(this_locator->interior_list); if (this_locator->exterior_list != NULL) BFT_FREE(this_locator->exterior_list); /* Prepare locator (MPI version) */ /*-------------------------------*/ #if defined(FVM_HAVE_MPI) MPI_Initialized(&mpi_flag); if (mpi_flag && this_locator->comm == MPI_COMM_NULL) mpi_flag = 0; if (mpi_flag) { /* Check that at least one of the local or distant nodal meshes is non-NULL, and at least one of the local or distant point sets is non null */ if (this_nodal != NULL) locdim[0] = this_nodal->dim; else locdim[0] = -1; if (n_points > 0) locdim[1] = dim; else locdim[1] = -1; MPI_Allreduce(locdim, globdim, 2, MPI_INT, MPI_MAX, this_locator->comm); if (globdim[0] < 0 || globdim[1] < 0) return; else if (this_nodal != NULL && globdim[1] != this_nodal->dim) bft_error(__FILE__, __LINE__, 0, _("Locator trying to use distant space dimension %d\n" "with local mesh dimension %d\n"), globdim[1], this_nodal->dim); else if (this_nodal == NULL && globdim[0] != dim) bft_error(__FILE__, __LINE__, 0, _("Locator trying to use local space dimension %d\n" "with distant mesh dimension %d\n"), dim, globdim[0]); /* Exchange extent information */ MPI_Comm_rank(this_locator->comm, &comm_rank); MPI_Comm_size(this_locator->comm, &comm_size); stride2 = dim * 2; /* Stride for one type of extent */ stride4 = dim * 4; /* Stride for element and vertex extents, end-to-end */ BFT_MALLOC(recvbuf, stride4*comm_size, double); MPI_Allgather(extents, stride4, MPI_DOUBLE, recvbuf, stride4, MPI_DOUBLE, this_locator->comm); /* Count and mark possible overlaps */ n_intersects = 0; BFT_MALLOC(intersect_rank, this_locator->n_ranks, int); for (i = 0 ; i < this_locator->n_ranks ; i++) { j = this_locator->start_rank + i; if ( (_intersect_extents(dim, extents + (2*dim), recvbuf + (j*stride4)) == true) || (_intersect_extents(dim, extents, recvbuf + (j*stride4) + (2*dim)) == true)) { intersect_rank[n_intersects] = j; n_intersects += 1; } } this_locator->n_intersects = n_intersects; BFT_MALLOC(this_locator->intersect_rank, this_locator->n_intersects, int); BFT_MALLOC(this_locator->intersect_extents, this_locator->n_intersects * stride2, double); for (i = 0 ; i < this_locator->n_intersects ; i++) { this_locator->intersect_rank[i] = intersect_rank[i]; /* Copy only distant element (and not point) extents */ for (j = 0 ; j < stride2 ; j++) this_locator->intersect_extents[i*stride2 + j] = recvbuf[intersect_rank[i]*stride4 + j]; } /* Free temporary memory */ BFT_FREE(intersect_rank); BFT_FREE(recvbuf); _locate_all_distant(this_locator, this_nodal, dim, n_points, point_index, point_coords); } #endif /* Prepare locator (local version) */ /*---------------------------------*/ if (!mpi_flag) { if (this_nodal == NULL || n_points == 0) return; stride2 = this_nodal->dim * 2; /* Count and mark possible overlaps */ if (_intersect_extents(this_nodal->dim, extents, extents + (2*dim)) == true) { this_locator->n_intersects = 1; BFT_MALLOC(this_locator->intersect_rank, 1, int); BFT_MALLOC(this_locator->intersect_extents, stride2, double); this_locator->intersect_rank[0] = 0; for (i = 0 ; i < stride2 ; i++) this_locator->intersect_extents[i] = extents[i]; _locate_all_local(this_locator, this_nodal, dim, n_points, point_index, point_coords); } } /* Update local_point_index values */ /*---------------------------------*/ if (this_locator->n_interior > 0) { fvm_lnum_t *reduced_index; BFT_MALLOC(reduced_index, n_points, fvm_lnum_t); for (i = 0; i < n_points; i++) reduced_index[i] = -1; assert( this_locator->local_points_idx[this_locator->n_intersects] == this_locator->n_interior); for (i = 0; i < this_locator->n_interior; i++) reduced_index[this_locator->interior_list[i] - 1] = i; /* Update this_locator->local_point_index[] so that it refers to an index in a dense [0, this_locator->n_interior] subset of the local points */ for (i = 0; i < this_locator->n_interior; i++) this_locator->local_point_index[i] = reduced_index[this_locator->local_point_index[i]]; BFT_FREE(reduced_index); } /* If an initial point index was given, update this_locator->interior_list and this_locator->exterior_list so that they refer to the same point set as that initial index (and not to an index within the selected point set) */ if (point_index != NULL) { for (i = 0; i < this_locator->n_interior; i++) this_locator->interior_list[i] = point_index[this_locator->interior_list[i] - 1]; for (i = 0; i < this_locator->n_exterior; i++) this_locator->exterior_list[i] = point_index[this_locator->exterior_list[i] - 1]; } } /*---------------------------------------------------------------------------- * Return number of distant points after locator initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * number of distant points. *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_locator_get_n_dist_points(const fvm_locator_t *const this_locator) { fvm_lnum_t retval = 0; if (this_locator != NULL) { if (this_locator->n_intersects != 0) retval = this_locator->distant_points_idx[this_locator->n_intersects]; } return retval; } /*---------------------------------------------------------------------------- * Return an array of local element numbers containing (or nearest to) * each distant point after locator initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * local element numbers associated with distant points (1 to n numbering). *----------------------------------------------------------------------------*/ const fvm_lnum_t * fvm_locator_get_dist_locations(const fvm_locator_t *const this_locator) { const fvm_lnum_t * retval = NULL; if (this_locator != NULL) { if (this_locator->n_ranks != 0) retval = this_locator->distant_point_location; } return retval; } /*---------------------------------------------------------------------------- * Return an array of coordinates of each distant point after * locator initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * coordinate array associated with distant points (interlaced). *----------------------------------------------------------------------------*/ const fvm_coord_t * fvm_locator_get_dist_coords(const fvm_locator_t *const this_locator) { const fvm_coord_t * retval = NULL; if (this_locator != NULL) { if (this_locator->n_intersects != 0) retval = this_locator->distant_point_coords; } return retval; } /*---------------------------------------------------------------------------- * Return number of points located after locator initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * number of points located. *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_locator_get_n_interior(const fvm_locator_t *const this_locator) { fvm_lnum_t retval = 0; if (this_locator != NULL) retval = this_locator->n_interior; return retval; } /*---------------------------------------------------------------------------- * Return list of points not located after locator initialization. * This list defines a subset of the point set used at initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * list of points not located (1 to n numbering). *----------------------------------------------------------------------------*/ const fvm_lnum_t * fvm_locator_get_interior_list(const fvm_locator_t *const this_locator) { return this_locator->interior_list; } /*---------------------------------------------------------------------------- * Return number of points not located after locator initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * number of points not located. *----------------------------------------------------------------------------*/ fvm_lnum_t fvm_locator_get_n_exterior(const fvm_locator_t *const this_locator) { return this_locator->n_exterior; } /*---------------------------------------------------------------------------- * Return list of points not located after locator initialization. * This list defines a subset of the point set used at initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * list of points not located (1 to n numbering). *----------------------------------------------------------------------------*/ const fvm_lnum_t * fvm_locator_get_exterior_list(const fvm_locator_t *const this_locator) { return this_locator->exterior_list; } /*---------------------------------------------------------------------------- * Discard list of points not located after locator initialization. * This list defines a subset of the point set used at initialization. * * parameters: * this_locator <-- pointer to locator structure * * returns: * list of points not located (1 to n numbering). *----------------------------------------------------------------------------*/ void fvm_locator_discard_exterior(fvm_locator_t *const this_locator) { this_locator->n_exterior = 0; BFT_FREE(this_locator->exterior_list); } /*---------------------------------------------------------------------------- * Distribute variable defined on distant points to processes owning * the original points (i.e. distant processes). * * The exchange is symmetric if both variables are defined, receive * only if distant_var is NULL, or send only if local_var is NULL. * * parameters: * this_locator <-- pointer to locator structure * distant_var <-- variable defined on distant points (ready to send) * local_var --> variable defined on local points (received) * type_size <-- sizeof (float or double) variable type * stride <-- dimension (1 for scalar, 3 for interlaced vector) *----------------------------------------------------------------------------*/ void fvm_locator_exchange_point_var(fvm_locator_t *const this_locator, void *const distant_var, void *const local_var, const size_t type_size, const size_t stride) { int mpi_flag = 0; #if defined(FVM_HAVE_MPI) MPI_Initialized(&mpi_flag); if (mpi_flag && this_locator->comm == MPI_COMM_NULL) mpi_flag = 0; if (mpi_flag) { MPI_Datatype datatype = MPI_DATATYPE_NULL; if (type_size == sizeof(double)) datatype = MPI_DOUBLE; else if (type_size == sizeof(float)) datatype = MPI_FLOAT; else bft_error(__FILE__, __LINE__, 0, _("type_size passed to fvm_locator_exchange_point_var() does\n" "not correspond to double or float.")); assert (datatype != MPI_DATATYPE_NULL); _exchange_point_var_distant(this_locator, distant_var, local_var, datatype, stride); } #endif /* defined(FVM_HAVE_MPI) */ if (!mpi_flag) _exchange_point_var_local(this_locator, distant_var, local_var, type_size, stride); } /*---------------------------------------------------------------------------- * Dump printout of a locator structure. * * parameters: * this_locator <-- pointer to structure that should be dumped *----------------------------------------------------------------------------*/ void fvm_locator_dump(const fvm_locator_t *this_locator) { int i; fvm_lnum_t j, k; const fvm_lnum_t *idx, *index, *loc; const fvm_coord_t *coords; if (this_locator == NULL) return; /* Basic information */ /*-------------------*/ bft_printf(_("\n" "Locator:\n\n" "Spatial dimension: %d\n" "Tolerance: %f\n" "Number of ranks of distant location: %d\n" "First rank of distant location: %d\n" "Number of intersecting distant ranks: %d\n"), this_locator->dim, this_locator->tolerance, this_locator->n_ranks, this_locator->start_rank, this_locator->n_intersects); #if defined(FVM_HAVE_MPI) if (this_locator->comm != MPI_COMM_NULL) bft_printf(_("\n" "Associated MPI communicator: %ld\n"), (long)(this_locator->comm)); #endif /* Arrays indexed by rank */ /*------------------------*/ for (i = 0; i < this_locator->n_intersects; i++) { bft_printf(_("\n" " Intersection %d with distant rank %d\n\n"), i+1, this_locator->intersect_rank[i]); bft_printf(_(" Distant rank extents:\n")); k = i * (this_locator->dim) * 2; for (j = 0; j < this_locator->dim; j++) bft_printf(_(" [%12.5e, %12.5e]\n"), this_locator->intersect_extents[k + j], this_locator->intersect_extents[k + this_locator->dim + j]); } if (this_locator->n_interior > 0) { bft_printf(_("\n Local point index (for receiving):\n\n")); idx = this_locator->local_points_idx; index = this_locator->local_point_index; for (i = 0; i < this_locator->n_intersects; i++) { if (idx[i+1] > idx[i]) { bft_printf("%6d (idx = %10d) %10d\n", i + 1, idx[i], index[idx[i]]); for (j = idx[i] + 1; j < idx[i + 1]; j++) bft_printf(" %10d\n", index[j]); } else { bft_printf("%6d (idx = %10d)\n", i + 1, idx[i]); } bft_printf(_(" end (idx = %10d)\n"), idx[this_locator->n_intersects]); } } if (this_locator->distant_points_idx != NULL) { idx = this_locator->distant_points_idx; loc = this_locator->distant_point_location; coords = this_locator->distant_point_coords; if (idx[this_locator->n_intersects] > 0) bft_printf(_("\n Distant point location:\n\n")); for (i = 0; i < this_locator->n_intersects; i++) { if (idx[i+1] > idx[i]) { if (this_locator->dim == 1) { bft_printf("%6d (idx = %10d) %10d [%12.5e]\n", i + 1, this_locator->intersect_rank[i], idx[i], loc[idx[i]], coords[idx[i]]); for (j = idx[i] + 1; j < idx[i + 1]; j++) bft_printf(" %10d [%12.5e]\n", loc[j], coords[j]); } else if (this_locator->dim == 2) { bft_printf("%6d (idx = %10d) %10d [%12.5e, %12.5e]\n", i + 1, idx[i], loc[idx[i]], coords[2*idx[i]], coords[2*idx[i]+1]); for (j = idx[i] + 1; j < idx[i + 1]; j++) bft_printf(" %10d [%12.5e, %12.5e]\n", loc[j], coords[2*j], coords[2*j+1]); } else if (this_locator->dim == 3) { bft_printf("%6d (idx = %10d) %10d [%12.5e, %12.5e, %12.5e]\n", i + 1, idx[i], loc[idx[i]], coords[3*idx[i]], coords[3*idx[i]+1], coords[3*idx[i]]+2); for (j = idx[i] + 1; j < idx[i + 1]; j++) bft_printf(" " "%10d [%12.5e, %12.5e, %12.5e]\n", loc[j], coords[3*j], coords[3*j+1], coords[3*j+2]); } } /* if (idx[i+1] > idx[i]) */ } if (idx[this_locator->n_intersects] > 0) bft_printf(_(" end (idx = %10d)\n"), idx[this_locator->n_intersects]); } /* Local arrays */ /*--------------*/ bft_printf(_("\n" " Number of local points successfully located: %d\n\n"), this_locator->n_interior); for (j = 0; j < this_locator->n_interior; j++) bft_printf(" %10d\n", this_locator->interior_list[j]); if (this_locator->n_interior > 0) bft_printf("\n"); bft_printf(_(" Number of local points not located: %d\n"), this_locator->n_exterior); for (j = 0; j < this_locator->n_exterior; j++) bft_printf(" %10d\n", this_locator->exterior_list[j]); if (this_locator->n_exterior > 0) bft_printf("\n"); } /*----------------------------------------------------------------------------*/ #ifdef __cplusplus } #endif /* __cplusplus */