diff --git a/src/Infrastructure/Mesh/include/ESMCI_Mesh_Vector_Regrid.h b/src/Infrastructure/Mesh/include/ESMCI_Mesh_Vector_Regrid.h new file mode 100644 index 0000000000..923202a090 --- /dev/null +++ b/src/Infrastructure/Mesh/include/ESMCI_Mesh_Vector_Regrid.h @@ -0,0 +1,66 @@ +// $Id$ +// +// Earth System Modeling Framework +// Copyright 2002-2022, University Corporation for Atmospheric Research, +// Massachusetts Institute of Technology, Geophysical Fluid Dynamics +// Laboratory, University of Michigan, National Centers for Environmental +// Prediction, Los Alamos National Laboratory, Argonne National Laboratory, +// NASA Goddard Space Flight Center. +// Licensed under the University of Illinois-NCSA License. +// +//============================================================================== + +//------------------------------------------------------------------------------ +// INCLUDES +//------------------------------------------------------------------------------ + +#ifndef ESMCI_Mesh_Vector_Regrid_H +#define ESMCI_Mesh_Vector_Regrid_H + +#include + +#include "ESMCI_Macros.h" +#include "ESMCI_F90Interface.h" +#include "ESMCI_LogErr.h" +#include "ESMCI_VM.h" +#include "VM/include/ESMC_VM.h" +#include "ESMCI_Array.h" + +#include +#include +#include +#include +#include +#include + +#include "Mesh/include/ESMCI_Mesh.h" +#include "ESMCI_CoordSys.h" +#include +#include "Mesh/include/Legacy/ESMCI_ParEnv.h" +#include + +//----------------------------------------------------------------------------- + // leave the following line as-is; it will insert the cvs ident string + // into the object file for tracking purposes. +// static const char *const version = "$Id$"; +//----------------------------------------------------------------------------- + +using namespace ESMCI; + +// TODO: move below +void get_vec_dims_for_vectorRegrid(ESMCI::Array &array, int &num_vec_dims, int *vec_dims_undist_seqind); + + + + +void create_vector_sparse_mat_from_reg_sparse_mat(int num_entries, int *iientries, double *factors, + int num_vec_dims, int *src_vec_dims_undist_seqind, int *dst_vec_dims_undist_seqind, + Mesh *src_mesh, PointList *src_pl, + Mesh *dst_mesh, PointList *dst_pl, + int &num_entries_vec, int *&iientries_vec, double *&factors_vec); + + + + + +#endif // ESMCI_Mesh_Vector_Regrid_H diff --git a/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C b/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C index 4114238f36..e36d07c0cd 100644 --- a/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C +++ b/src/Infrastructure/Mesh/src/ESMCI_Mesh_Regrid_Glue.C @@ -38,9 +38,9 @@ #include "Mesh/include/ESMCI_MathUtil.h" #include "Mesh/include/Legacy/ESMCI_Phedra.h" #include "Mesh/include/ESMCI_Mesh_Regrid_Glue.h" +#include "Mesh/include/ESMCI_Mesh_Vector_Regrid.h" #include "Mesh/include/Legacy/ESMCI_MeshMerge.h" #include "Mesh/include/ESMCI_Mesh_GToM_Glue.h" -#include "Mesh/include/ESMCI_DInfo.h" #include @@ -59,8 +59,6 @@ using namespace ESMCI; - - // prototypes from below static bool all_mesh_node_ids_in_wmat(PointList *pointlist, WMat &wts, int *missing_id); static bool all_mesh_elem_ids_in_wmat(Mesh *mesh, WMat &wts, int *missing_id); @@ -84,39 +82,6 @@ void CpMeshDataToArray(Grid &grid, int staggerLoc, ESMCI::Mesh &mesh, ESMCI::Arr void CpMeshElemDataToArray(Grid &grid, int staggerloc, ESMCI::Mesh &mesh, ESMCI::Array &array, MEField<> *dataToArray); void PutElemAreaIntoArray(Grid &grid, int staggerLoc, ESMCI::Mesh &mesh, ESMCI::Array &array); -// Get information about vector dimensions for vectorRegrid capability -// TODO: move below -static void _get_vec_dims_for_vectorRegrid(ESMCI::Array &array, int &num_vec_dims, int *vec_dims_undist_seqind) { - - // Get undisributed dimension sizes - const int *undistLBound=array.getUndistLBound(); - const int *undistUBound=array.getUndistUBound(); - - // We error checked earlier that there is only one undist bound, but double check here that these aren't NULL - ThrowRequire(undistLBound != NULL); - ThrowRequire(undistUBound != NULL); - - // Calculate size of vector dims - num_vec_dims=undistUBound[0]-undistLBound[0]+1; - - // Fill undistibuted seqinds for vector dims - // (This is easy now, but will be harder as the number of undist. dims gets larger than 1.) - for (int i=0; i { - public: - CoordFromIdEntry_just_id_less() {} - bool operator()(const CoordFromIdEntry &l, const CoordFromIdEntry &r) { - return l.id < r.id; - } - }; - - - - // Data - bool is_searchable; - std::vector searchable; // After committing this will contain a sorted list for searching - - - - -public: - - // Create empty - CoordFromId(): is_searchable(false) { - - } - - - // Add points to query structure from Mesh - void add(Mesh *mesh, MeshObj::id_type obj_type); - - // Add points to query structure from PointList - void add(PointList *pl); - - // Make searchable - void make_searchable(); - - // This method makes the structure searchable for a particular set - // of possibly non-local ids. In this case, for convenience, the ids - // are passed in in the form of the weight matrix index list. sord - // give whether the ids from the source (sord=0) or destination (sord=1) should - // be used. - void make_searchable(int num_entries, int *iientries, int sord); - - // Search - // Returns true if id has been found. In that case fills coords_out with coords of - // point for that id, coords_out must be of size >=3. - bool search(int search_id, double *coords_out) { - - // Find id in searchable list - std::vector::iterator ei; // break into two lines, so easier to read - ei = std::lower_bound(searchable.begin(), - searchable.end(), - CoordFromIdEntry(search_id, 0.0, 0.0, 0.0), - CoordFromIdEntry_just_id_less()); - - - // If within list - if (ei != searchable.end()) { - CoordFromIdEntry &lb_cfie = *ei; - - // If the ids match, then we've found an answer - if (lb_cfie.id == search_id) { - - // Get coords from entry - lb_cfie.coord.get_coords(coords_out); - - // Report success - return true; - } - } - - // Didn't find id, so report that - return false; - } - -}; - - - - -void CoordFromId::add(Mesh *mesh, MeshObj::id_type obj_type) { - - // If already searchable, don't allow more to be added - if (is_searchable) { - Throw() << "CoordFromId object has already been made searchable, more points can't be added."; - } - - // Get spatial dimension - int sdim=mesh->spatial_dim(); - - // Add based on obj_type - if (obj_type == MeshObj::NODE) { - - // Reserve to the correct size - searchable.reserve(mesh->num_nodes()); - - // Get coordinate data - MEField<> *node_coords=mesh->GetField("coordinates"); - ThrowRequire(node_coords != NULL); - - // Add - if (sdim == 2) { - Mesh::iterator ni = mesh->node_begin(), ne = mesh->node_end(); - for (; ni != ne; ++ni) { - MeshObj &node = *ni; - - // Skip if not local - if (!GetAttr(node).is_locally_owned()) continue; - - // Get id - int id=node.get_id(); - - // Get pointer to coords - double *c = node_coords->data(node); - - // Add to list - searchable.push_back(CoordFromIdEntry(id,c[0],c[1],0.0)); - } - } else if (sdim == 3) { - Mesh::iterator ni = mesh->node_begin(), ne = mesh->node_end(); - for (; ni != ne; ++ni) { - MeshObj &node = *ni; - - // Skip if not local - if (!GetAttr(node).is_locally_owned()) continue; - - // Get id - int id=node.get_id(); - - // Get pointer to coords - double *c = node_coords->data(node); - - // Add to list - searchable.push_back(CoordFromIdEntry(id,c[0],c[1],c[2])); - } - } else { - Throw() << "Meshes of spatial dim= "<num_elems()); - - // Get element coordinate data - MEField<> *elem_coords=mesh->GetField("elem_coordinates"); - if (elem_coords == NULL) Throw() << "Vector regridding not supported on Mesh elements when the elements don't have center coordinates."; - - // Add based on spatial dim - if (sdim == 2) { - Mesh::iterator ei = mesh->elem_begin(), ee = mesh->elem_end(); - for (; ei != ee; ++ei) { - MeshObj &elem = *ei; - - // Skip if not local - if (!GetAttr(elem).is_locally_owned()) continue; - - // Get id - int id=elem.get_id(); - - // Get pointer to coords - double *c = elem_coords->data(elem); - - // Add to list - searchable.push_back(CoordFromIdEntry(id,c[0],c[1],0.0)); - } - } else if (sdim == 3) { - Mesh::iterator ei = mesh->elem_begin(), ee = mesh->elem_end(); - for (; ei != ee; ++ei) { - MeshObj &elem = *ei; - - // Skip if not local - if (!GetAttr(elem).is_locally_owned()) continue; - - // Get id - int id=elem.get_id(); - - // Get pointer to coords - double *c = elem_coords->data(elem); - - // Add to list - searchable.push_back(CoordFromIdEntry(id,c[0],c[1],c[2])); - } - } else { - Throw() << "Meshes of spatial dim= "<get_coord_dim(); - - // Reserve to the correct size - searchable.reserve(pl->get_curr_num_pts()); - - // Add points based on spatial dim - if (sdim == 2) { - // Loop adding ids and points - int num_pts = pl->get_curr_num_pts(); - for (int i=0; iget_id(i); - - // Get point coords - const double *c = pl->get_coord_ptr(i); - - // Add to list - searchable.push_back(CoordFromIdEntry(id,c[0],c[1],0.0)); - } - } else if (sdim == 3) { - // Loop adding ids and points - int num_pts = pl->get_curr_num_pts(); - for (int i=0; iget_id(i); - - // Get point coords - const double *c = pl->get_coord_ptr(i); - - // Add to list - searchable.push_back(CoordFromIdEntry(id,c[0],c[1],c[2])); - } - } else { - Throw() << "Geometries of spatial dim= "< 1)) { - Throw() << "Value= "< search_ids; - - // Reserve a reasonable size - search_ids.reserve(num_entries); - - // Get list of ids on this PET - int pos=sord; - for (int i=0; i< num_entries; i++) { - search_ids.push_back(iientries[pos]); - pos += 2; // Advance to next pair of entries - } - - // Sort to allow unique to work - std::sort(search_ids.begin(),search_ids.end()); - - // Unique ids to make search and comm more efficient - search_ids.erase(std::unique(search_ids.begin(), search_ids.end()), search_ids.end()); - - - //// Set up class for distributed search - - // Allocate array where the results will go - Coord *results = new Coord[search_ids.size()]; - - // Declare distributed search structure - { - DInfo di_search; - - // Reserve space - di_search.reserve(searchable.size()); - - // Add search items - for (CoordFromIdEntry cfie : searchable) { - di_search.add(cfie.id,cfie.coord); - } - - // Commit to make searchable - di_search.commit(); - - - //// Do distributed search on local set of ids - - // Declare not found info value - Coord bad_value(0.0,0.0,0.0); - - // Do search - // (Because of true will return error if id not found) - if (search_ids.empty()) { - di_search.search(0, NULL, true, bad_value, results); - } else { - di_search.search(search_ids.size(), search_ids.data(), true, bad_value, results); - } - - } // End of block to get rid of search struct - - - //// Use distributed results to reconstruct search vector - - // Empty current search vector - searchable.clear(); - - // Add new distributed objects to search vector - for (auto i=0; i +#include +#include + + //------------------------------------------------------------------------------ +//BOP +// !DESCRIPTION: +// +// + //EOP +//------------------------------------------------------------------------- + + +using namespace ESMCI; + + + +// Get information about vector dimensions for vectorRegrid capability +void get_vec_dims_for_vectorRegrid(ESMCI::Array &array, int &num_vec_dims, int *vec_dims_undist_seqind) { + + // Get undisributed dimension sizes + const int *undistLBound=array.getUndistLBound(); + const int *undistUBound=array.getUndistUBound(); + + // We error checked earlier that there is only one undist bound, but double check here that these aren't NULL + ThrowRequire(undistLBound != NULL); + ThrowRequire(undistUBound != NULL); + + // Calculate size of vector dims + num_vec_dims=undistUBound[0]-undistLBound[0]+1; + + // Fill undistibuted seqinds for vector dims + // (This is easy now, but will be harder as the number of undist. dims gets larger than 1.) + for (int i=0; i { + public: + CoordFromIdEntry_just_id_less() {} + bool operator()(const CoordFromIdEntry &l, const CoordFromIdEntry &r) { + return l.id < r.id; + } + }; + + + + // Data + bool is_searchable; + std::vector searchable; // After committing this will contain a sorted list for searching + + + + +public: + + // Create empty + CoordFromId(): is_searchable(false) { + + } + + + // Add points to query structure from Mesh + void add(Mesh *mesh, MeshObj::id_type obj_type); + + // Add points to query structure from PointList + void add(PointList *pl); + + // Make searchable + void make_searchable(); + + // This method makes the structure searchable for a particular set + // of possibly non-local ids. In this case, for convenience, the ids + // are passed in in the form of the weight matrix index list. sord + // give whether the ids from the source (sord=0) or destination (sord=1) should + // be used. + void make_searchable(int num_entries, int *iientries, int sord); + + // Search + // Returns true if id has been found. In that case fills coords_out with coords of + // point for that id, coords_out must be of size >=3. + bool search(int search_id, double *coords_out) { + + // Find id in searchable list + std::vector::iterator ei; // break into two lines, so easier to read + ei = std::lower_bound(searchable.begin(), + searchable.end(), + CoordFromIdEntry(search_id, 0.0, 0.0, 0.0), + CoordFromIdEntry_just_id_less()); + + + // If within list + if (ei != searchable.end()) { + CoordFromIdEntry &lb_cfie = *ei; + + // If the ids match, then we've found an answer + if (lb_cfie.id == search_id) { + + // Get coords from entry + lb_cfie.coord.get_coords(coords_out); + + // Report success + return true; + } + } + + // Didn't find id, so report that + return false; + } + +}; + + + + +void CoordFromId::add(Mesh *mesh, MeshObj::id_type obj_type) { + + // If already searchable, don't allow more to be added + if (is_searchable) { + Throw() << "CoordFromId object has already been made searchable, more points can't be added."; + } + + // Get spatial dimension + int sdim=mesh->spatial_dim(); + + // Add based on obj_type + if (obj_type == MeshObj::NODE) { + + // Reserve to the correct size + searchable.reserve(mesh->num_nodes()); + + // Get coordinate data + MEField<> *node_coords=mesh->GetField("coordinates"); + ThrowRequire(node_coords != NULL); + + // Add + if (sdim == 2) { + Mesh::iterator ni = mesh->node_begin(), ne = mesh->node_end(); + for (; ni != ne; ++ni) { + MeshObj &node = *ni; + + // Skip if not local + if (!GetAttr(node).is_locally_owned()) continue; + + // Get id + int id=node.get_id(); + + // Get pointer to coords + double *c = node_coords->data(node); + + // Add to list + searchable.push_back(CoordFromIdEntry(id,c[0],c[1],0.0)); + } + } else if (sdim == 3) { + Mesh::iterator ni = mesh->node_begin(), ne = mesh->node_end(); + for (; ni != ne; ++ni) { + MeshObj &node = *ni; + + // Skip if not local + if (!GetAttr(node).is_locally_owned()) continue; + + // Get id + int id=node.get_id(); + + // Get pointer to coords + double *c = node_coords->data(node); + + // Add to list + searchable.push_back(CoordFromIdEntry(id,c[0],c[1],c[2])); + } + } else { + Throw() << "Meshes of spatial dim= "<num_elems()); + + // Get element coordinate data + MEField<> *elem_coords=mesh->GetField("elem_coordinates"); + if (elem_coords == NULL) Throw() << "Vector regridding not supported on Mesh elements when the elements don't have center coordinates."; + + // Add based on spatial dim + if (sdim == 2) { + Mesh::iterator ei = mesh->elem_begin(), ee = mesh->elem_end(); + for (; ei != ee; ++ei) { + MeshObj &elem = *ei; + + // Skip if not local + if (!GetAttr(elem).is_locally_owned()) continue; + + // Get id + int id=elem.get_id(); + + // Get pointer to coords + double *c = elem_coords->data(elem); + + // Add to list + searchable.push_back(CoordFromIdEntry(id,c[0],c[1],0.0)); + } + } else if (sdim == 3) { + Mesh::iterator ei = mesh->elem_begin(), ee = mesh->elem_end(); + for (; ei != ee; ++ei) { + MeshObj &elem = *ei; + + // Skip if not local + if (!GetAttr(elem).is_locally_owned()) continue; + + // Get id + int id=elem.get_id(); + + // Get pointer to coords + double *c = elem_coords->data(elem); + + // Add to list + searchable.push_back(CoordFromIdEntry(id,c[0],c[1],c[2])); + } + } else { + Throw() << "Meshes of spatial dim= "<get_coord_dim(); + + // Reserve to the correct size + searchable.reserve(pl->get_curr_num_pts()); + + // Add points based on spatial dim + if (sdim == 2) { + // Loop adding ids and points + int num_pts = pl->get_curr_num_pts(); + for (int i=0; iget_id(i); + + // Get point coords + const double *c = pl->get_coord_ptr(i); + + // Add to list + searchable.push_back(CoordFromIdEntry(id,c[0],c[1],0.0)); + } + } else if (sdim == 3) { + // Loop adding ids and points + int num_pts = pl->get_curr_num_pts(); + for (int i=0; iget_id(i); + + // Get point coords + const double *c = pl->get_coord_ptr(i); + + // Add to list + searchable.push_back(CoordFromIdEntry(id,c[0],c[1],c[2])); + } + } else { + Throw() << "Geometries of spatial dim= "< 1)) { + Throw() << "Value= "< search_ids; + + // Reserve a reasonable size + search_ids.reserve(num_entries); + + // Get list of ids on this PET + int pos=sord; + for (int i=0; i< num_entries; i++) { + search_ids.push_back(iientries[pos]); + pos += 2; // Advance to next pair of entries + } + + // Sort to allow unique to work + std::sort(search_ids.begin(),search_ids.end()); + + // Unique ids to make search and comm more efficient + search_ids.erase(std::unique(search_ids.begin(), search_ids.end()), search_ids.end()); + + + //// Set up class for distributed search + + // Allocate array where the results will go + Coord *results = new Coord[search_ids.size()]; + + // Declare distributed search structure + { + DInfo di_search; + + // Reserve space + di_search.reserve(searchable.size()); + + // Add search items + for (CoordFromIdEntry cfie : searchable) { + di_search.add(cfie.id,cfie.coord); + } + + // Commit to make searchable + di_search.commit(); + + + //// Do distributed search on local set of ids + + // Declare not found info value + Coord bad_value(0.0,0.0,0.0); + + // Do search + // (Because of true will return error if id not found) + if (search_ids.empty()) { + di_search.search(0, NULL, true, bad_value, results); + } else { + di_search.search(search_ids.size(), search_ids.data(), true, bad_value, results); + } + + } // End of block to get rid of search struct + + + //// Use distributed results to reconstruct search vector + + // Empty current search vector + searchable.clear(); + + // Add new distributed objects to search vector + for (auto i=0; i