/*============================================================================
* 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 <fvm_config.h>
/*----------------------------------------------------------------------------
* Standard C library headers
*----------------------------------------------------------------------------*/
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#if defined(FVM_HAVE_MPI)
#include <mpi.h>
#endif
/*----------------------------------------------------------------------------
* BFT library headers
*----------------------------------------------------------------------------*/
#include <bft_error.h>
#include <bft_mem.h>
#include <bft_printf.h>
/*----------------------------------------------------------------------------
* Local headers
*----------------------------------------------------------------------------*/
#include <fvm_config_defs.h>
#include <fvm_defs.h>
#include <fvm_nodal.h>
#include <fvm_nodal_priv.h>
#include <fvm_parall.h>
/*----------------------------------------------------------------------------
* Header for the current file
*----------------------------------------------------------------------------*/
#include <fvm_locator.h>
/*----------------------------------------------------------------------------*/
#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 <math.h> 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 */
syntax highlighted by Code2HTML, v. 0.9.1