diff --git a/src/syrthes-kernel/include/syr_cfd_point_location.h b/src/syrthes-kernel/include/syr_cfd_point_location.h index 676649e..6b5dd35 100755 --- a/src/syrthes-kernel/include/syr_cfd_point_location.h +++ b/src/syrthes-kernel/include/syr_cfd_point_location.h @@ -88,50 +88,29 @@ syr_cfd_point_location_extents(const void *mesh, * than to previously encountered elements. * * parameters: - * mesh <-- pointer to 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, - * and > 1 if outside a volume element, or absolute - * distance to a surface element (size: n_points) + * mesh <-- pointer to mesh representation structure + * tolerance_base <-- associated base tolerance (used for bounding + * box check only, not for location test) + * tolerance_fraction <-- associated fraction of element bounding boxes + * added to tolerance + * n_points <-- number of points to locate + * point_coords <-- point coordinates + * point_tag <-- optional point tag (unused here) + * 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, + * and > 1 if outside a volume element, or absolute + * distance to a surface element (size: n_points) *----------------------------------------------------------------------------*/ void syr_cfd_point_location_contain(const void *mesh, - double tolerance, - ple_lnum_t n_points, - const ple_coord_t point_coords[], - ple_lnum_t location[], - float distance[]); - -/*---------------------------------------------------------------------------- - * Find elements in a given mesh closest to points: updates the - * location[] and distance[] arrays associated with a set of points - * for points that are closer to an element of this mesh than to previously - * encountered elements. - * - * This function currently only handles elements of lower dimension than - * the spatial dimension. - * - * parameters: - * mesh <-- pointer to mesh representation structure - * 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, or absolute - * distance to a surface element (size: n_points) - *----------------------------------------------------------------------------*/ - -void -syr_cfd_point_location_closest(const void *mesh, + float tolerance_base, + float tolerance_fraction, ple_lnum_t n_points, const ple_coord_t point_coords[], + const ple_lnum_t point_tag[], ple_lnum_t location[], float distance[]); diff --git a/src/syrthes-kernel/src/syr_cfd_coupling.c b/src/syrthes-kernel/src/syr_cfd_coupling.c index 4396102..e7ac550 100644 --- a/src/syrthes-kernel/src/syr_cfd_coupling.c +++ b/src/syrthes-kernel/src/syr_cfd_coupling.c @@ -32,9 +32,9 @@ /*=========================================================================== * Main API functions for coupling between Syrthes and Code_Saturne - * AUTHORS : J. Bonelle, Y Fournier + * AUTHORS : J. Bonelle, Y Fournier * - * Library: Syrthes Copyright EDF 2008-2010 + * Library: Syrthes Copyright EDF 2008-2014 *===========================================================================*/ #if defined(_SYRTHES_CFD_) @@ -150,7 +150,8 @@ struct _syr_cfd_coupling_t { int conservativity; /* Conservativity forcing flag */ int allow_nearest; /* Allow nearest flag */ - double tolerance; /* Tolerance for location */ + double tolerance[2]; /* Tolerance for location + 0: absolute, 1: relative*/ }; @@ -378,8 +379,8 @@ _init_mpi_comm(syr_cfd_coupling_t *cfd_coupling, { int mpi_flag = 0; - int scheme_min = 1; - int scheme_max = 1; + int scheme_min = 2; + int scheme_max = 2; char volume_flag = ' '; char boundary_flag = ' '; @@ -408,7 +409,7 @@ _init_mpi_comm(syr_cfd_coupling_t *cfd_coupling, int distant_range[2] = {-1, -1}; ple_coupling_mpi_intracomm_create(MPI_COMM_WORLD, - _syr_coupling_loc_comm, + _syr_coupling_loc_comm, cfd_coupling->root_dist_rank, &(cfd_coupling->comm), local_range, @@ -453,53 +454,55 @@ _init_mpi_comm(syr_cfd_coupling_t *cfd_coupling, /* Leave id 17 open for volume conservativity flag */ - /* For newer versions of Code_Saturne (3.2 +), exchange scheme version information - is transferred, with a min an max scheme number so as to allow forward/backward - compatibility across some scheme changes: both codes should follow the maximum - common handled scheme. Scheme numbers are only to be increased when additional + /* Exchange scheme version information, with a min and max scheme number so as to + allow forward/backward compatibility across some scheme changes: + both codes should follow the maximum common handled scheme. + Scheme numbers are only to be increased when additional or different data are passed. Changes not breaking a given scheme (such as better error handling by passing additional messages instead of aborting immediately do not require increasing its number. As the first versioned scheme is the current version, a blank in the scheme version character places is assumed to be scheme 1. */ - if (op_name_recv[18] == ' ' || op_name_recv[18] == '\1') - scheme_min = 1; + if (op_name_recv[18] == ' ' || op_name_recv[18] == '\2') + scheme_min = 2; else scheme_min = (int)(op_name_recv[18]); - if (op_name_recv[19] == ' ' || op_name_recv[19] == '\1') - scheme_max = 1; + if (op_name_recv[19] == ' ' || op_name_recv[19] == '\2') + scheme_max = 2; else scheme_max = (int)(op_name_recv[19]); - /* For Code_Saturne 3.2 and later, tolerance information is exchanged also */ + /* Tolerance information is exchanged also */ if (op_name_recv[20] == '0') cfd_coupling->allow_nearest = 0; else cfd_coupling->allow_nearest = 1; - cfd_coupling->tolerance = 0.1; /* default */ + cfd_coupling->tolerance[0] = 0.; /* default */ + cfd_coupling->tolerance[1] = 0.1; /* default */ if (op_name_recv[21] == '(') { int i = 0, j = 22; char tolerance_str[32] = ""; while (j < 33) { - if (op_name_recv[j] == ')') + if (op_name_recv[j] == ')') break; - else + else tolerance_str[i++] = op_name_recv[j++]; - } + } tolerance_str[i++] = '\0'; - cfd_coupling->tolerance = atof(tolerance_str); - if (cfd_coupling->tolerance <= 0) - cfd_coupling->tolerance = 0.1; + cfd_coupling->tolerance[1] = atof(tolerance_str); + if (cfd_coupling->tolerance[1] <= 0) + cfd_coupling->tolerance[1] = 0.1; } ple_printf(_(" communication scheme version: %d\n" " conservativity flag: %d\n" " allow nonmatching: %d\n" - " tolerance: %-6.2g\n"), + " tolerance: %-6.2g %-6.2g\n"), scheme_max, cfd_coupling->conservativity, - cfd_coupling->allow_nearest, cfd_coupling->tolerance); + cfd_coupling->allow_nearest, + cfd_coupling->tolerance[0], cfd_coupling->tolerance[1]); fflush(stdout); /* Check compatibility */ @@ -1101,6 +1104,56 @@ _interpolation_destroy(syr_cfd_interpolation_t **ip) } /*---------------------------------------------------------------------------- + * Check if coupling location is complete + * + * arguments: + * coupl <-> coupling helper structure + * ip <-> coupling interpolation helper structure + + * returns: + * 1 if location is incomplete, 0 if it is complete + *----------------------------------------------------------------------------*/ + +static int +_is_location_incomplete(syr_cfd_coupling_t *coupl, + syr_cfd_interpolation_t *ip) +{ + /* Local variables */ + + int location_incomplete = 0; + int flag_glob = 0; + + char op_name_send[32 + 1]; + char op_name_recv[32 + 1]; + + if (ple_locator_get_n_exterior(ip->locator) > 0) + location_incomplete = 1; + +#if defined(_SYRTHES_MPI_) + { + int n_ranks = 1; + MPI_Comm_size(_syr_coupling_loc_comm, &n_ranks); + if (n_ranks > 1) { + int flag_loc = location_incomplete; + MPI_Allreduce(&flag_loc, &location_incomplete, 1, MPI_INT, MPI_MAX, + _syr_coupling_loc_comm); + } + } +#endif + + if (location_incomplete > 0) + strcpy(op_name_send, "coupling:location:incomplete"); + else + strcpy(op_name_send, "coupling:location:ok"); + + _exchange_sync(coupl, op_name_send, op_name_recv); + if (!strcmp(op_name_recv, "coupling:location:incomplete")) + location_incomplete = 1; + + return location_incomplete; +} + +/*---------------------------------------------------------------------------- * Initialize interpolation structure. * * Local elements are used as a location basis for values of distant @@ -1142,8 +1195,6 @@ _interpolation_init(syr_cfd_coupling_t *coupl, ple_coord_t *elt_centers = NULL; float *dist_to_cfd = NULL; - ple_mesh_elements_closest_t *locate_on_closest = NULL; - /* Check if cells are coupled, and return if this is not the case */ if (n_coupled_elts > 0) { @@ -1166,16 +1217,12 @@ _interpolation_init(syr_cfd_coupling_t *coupl, /* Initialization */ - if (elt_dim == dim - 1 && coupl->allow_nearest == 1) - locate_on_closest = syr_cfd_point_location_closest; - #if defined(_SYRTHES_MPI_) - ip->locator = ple_locator_create(coupl->tolerance, - coupl->comm, + ip->locator = ple_locator_create(coupl->comm, coupl->n_dist_ranks, coupl->root_dist_rank); #else - ip->locator = ple_locator_create(tolerance); + ip->locator = ple_locator_create(); #endif if (dim == 3) { @@ -1223,14 +1270,46 @@ _interpolation_init(syr_cfd_coupling_t *coupl, ple_locator_set_mesh(ip->locator, elts, + NULL, + coupl->tolerance[0], + coupl->tolerance[1], dim, n_coupled_elts, NULL, + NULL, elt_centers, dist_to_cfd, syr_cfd_point_location_extents, - syr_cfd_point_location_contain, - locate_on_closest); + syr_cfd_point_location_contain); + + flag_loc = _is_location_incomplete(coupl, ip); + + if (coupl->allow_nearest) { + + float tolerance = coupl->tolerance[1]; + + while (flag_loc > 0) { + + tolerance *= 4; + + ple_locator_extend_search(ip->locator, + elts, + NULL, + 0, + tolerance, + n_coupled_elts, + NULL, + NULL, + elt_centers, + dist_to_cfd, + syr_cfd_point_location_extents, + syr_cfd_point_location_contain); + + } + + flag_loc = _is_location_incomplete(coupl, ip); + + } if (elt_dim == dim - 1) { ple_locator_exchange_point_var(ip->locator, @@ -1547,7 +1626,7 @@ syr_cfd_coupling_finalize(void) syr_cfd_coupling_t * syr_cfd_coupling_create(int dim, const char *app_name, - syr_cfd_coupling_type_t type) + syr_cfd_coupling_type_t type) { /* local variables */ int len; @@ -1585,7 +1664,7 @@ syr_cfd_coupling_create(int dim, coupling->cell_ip = _interpolation_create(); } else if (type == SYR_NO_COUPLING) - ple_printf(_("\n Warning: No coupling type defined\n")); + ple_printf(_("\n Warning: No coupling type defined\n")); else ple_error(__FILE__, __LINE__, 0, _("\n Unknown coupling type.\n --> abort.")); @@ -1594,7 +1673,8 @@ syr_cfd_coupling_create(int dim, coupling->conservativity = 0; coupling->allow_nearest = 1; - coupling->tolerance = 0.1; + coupling->tolerance[0] = 0.; + coupling->tolerance[1] = 0.1; return coupling; } diff --git a/src/syrthes-kernel/src/syr_cfd_point_location.c b/src/syrthes-kernel/src/syr_cfd_point_location.c index e2e803c..c7f3db7 100644 --- a/src/syrthes-kernel/src/syr_cfd_point_location.c +++ b/src/syrthes-kernel/src/syr_cfd_point_location.c @@ -1704,59 +1704,6 @@ _mesh_locate_3d(const syr_cfd_mesh_t *mesh, } /*---------------------------------------------------------------------------- - * Find elements in a given section closest to 3d 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: - * mesh <-- pointer to mesh representation structure - * point_coords <-- point coordinates - * n_point_ids <-- number of points to locate - * point_id <-- ids of points to locate - * 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, or absolute - * distance to a surface element (size: n_points) - *----------------------------------------------------------------------------*/ - -static void -_mesh_closest_3d(const syr_cfd_mesh_t *mesh, - const ple_coord_t point_coords[], - int n_point_ids, - const int point_ids[], - ple_lnum_t location[], - float distance[]) -{ - int i; - - const int stride = syr_cfd_mesh_n_vertices_element[mesh->element_type]; - - if (mesh->n_elements == 0) - return; - - if (mesh->element_type == SYR_CFD_TETRA) - ple_error(__FILE__, __LINE__, 0, - _("Locating volume elements closest to points not handled yet")); - - else if (mesh->element_type == SYR_CFD_TRIA) { - - for (i = 0; i < mesh->n_elements; i++) - _locate_on_triangle_3d(i + 1, - mesh->vertex_num + i*stride, - mesh->vertex_coords, - point_coords, - n_point_ids, - point_ids, - (HUGE_VAL / 4.), - location, - distance); - - } -} - -/*---------------------------------------------------------------------------- * Find elements in a given section containing 2d 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 @@ -1846,60 +1793,6 @@ _mesh_locate_2d(const syr_cfd_mesh_t *mesh, } /* End of loop on elements */ } -/*---------------------------------------------------------------------------- - * Find elements in a given section closest to 3d 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: - * mesh <-- pointer to mesh representation structure - * point_coords <-- point coordinates - * n_point_ids <-- number of points to locate - * point_id <-- ids of points to locate - * 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, or absolute - * distance to a line element (size: n_points) - *----------------------------------------------------------------------------*/ - -static void -_mesh_closest_2d(const syr_cfd_mesh_t *mesh, - const ple_coord_t point_coords[], - int n_point_ids, - const int point_id[], - ple_lnum_t location[], - float distance[]) -{ - int i; - - const int stride = syr_cfd_mesh_n_vertices_element[mesh->element_type]; - - /* Return immediately if nothing to do for this rank */ - - if (mesh->n_elements == 0 || mesh->element_type != SYR_CFD_EDGE) - return; - - /* Main loop on elements */ - - for (i = 0; i < mesh->n_elements; i++) { - - /* Locate on edge */ - - _locate_on_edge_2d(i + 1, - mesh->vertex_num + i*stride, - mesh->vertex_coords, - point_coords, - n_point_ids, - point_id, - -1.0, - location, - distance); - - } /* End of loop on elements */ -} - /*============================================================================ * Public function definitions *============================================================================*/ @@ -2031,23 +1924,29 @@ syr_cfd_point_location_extents(const void *mesh, * than to previously encountered elements. * * parameters: - * mesh <-- pointer to 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, - * and > 1 if outside a volume element, or absolute - * distance to a surface element (size: n_points) + * mesh <-- pointer to mesh representation structure + * tolerance_base <-- associated base tolerance (used for bounding + * box check only, not for location test) + * tolerance_fraction <-- associated fraction of element bounding boxes + * added to tolerance + * n_points <-- number of points to locate + * point_coords <-- point coordinates + * point_tag <-- optional point tag (unused here) + * 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, + * and > 1 if outside a volume element, or absolute + * distance to a surface element (size: n_points) *----------------------------------------------------------------------------*/ void syr_cfd_point_location_contain(const void *mesh, - double tolerance, + float tolerance_base, + float tolerance_fraction, ple_lnum_t n_points, const ple_coord_t point_coords[], + const ple_lnum_t point_tag[], ple_lnum_t location[], float distance[]) { @@ -2074,7 +1973,7 @@ syr_cfd_point_location_contain(const void *mesh, /* Locate for all sections */ _mesh_locate_3d(m, - tolerance, + tolerance_fraction, point_coords, &octree, points_in_extents, @@ -2091,7 +1990,7 @@ syr_cfd_point_location_contain(const void *mesh, _quadtree_t quadtree = _build_quadtree(n_points, point_coords); _mesh_locate_2d(m, - tolerance, + tolerance_fraction, point_coords, &quadtree, points_in_extents, @@ -2105,72 +2004,6 @@ syr_cfd_point_location_contain(const void *mesh, } } -/*---------------------------------------------------------------------------- - * Find elements in a given mesh closest to points: updates the - * location[] and distance[] arrays associated with a set of points - * for points that are closer to an element of this mesh than to previously - * encountered elements. - * - * This function currently only handles elements of lower dimension than - * the spatial dimension. - * - * parameters: - * mesh <-- pointer to mesh representation structure - * 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, or absolute - * distance to a surface element (size: n_points) - *----------------------------------------------------------------------------*/ - -void -syr_cfd_point_location_closest(const void *mesh, - ple_lnum_t n_points, - const ple_coord_t point_coords[], - ple_lnum_t location[], - float distance[]) -{ - if (mesh == NULL) - return; - - else { - - int j; - int *point_ids = NULL; - - const syr_cfd_mesh_t *m = mesh; - - PLE_MALLOC(point_ids, n_points, int); - for (j = 0; j < n_points; j++) - point_ids[j] = j; - - /* Use brute force for closest 3d point location */ - - if (m->dim == 3) - _mesh_closest_3d(m, - point_coords, - n_points, - point_ids, - location, - distance); - - /* Use brute force for closest 2d point location */ - - else if (m->dim == 2) - _mesh_closest_2d(m, - point_coords, - n_points, - point_ids, - location, - distance); - - if (point_ids != NULL) - PLE_FREE(point_ids); - } -} - /*----------------------------------------------------------------------------*/ #undef _DOT_PRODUCT