diff --git a/src/IO/H5/ExtendConnectivityHelpers.cpp b/src/IO/H5/ExtendConnectivityHelpers.cpp index a15501f4d28d..a1d8e1ade461 100644 --- a/src/IO/H5/ExtendConnectivityHelpers.cpp +++ b/src/IO/H5/ExtendConnectivityHelpers.cpp @@ -13,6 +13,7 @@ #include #include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Structure/SegmentId.hpp" #include "IO/H5/VolumeData.hpp" #include "NumericalAlgorithms/Spectral/Basis.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" @@ -358,6 +359,482 @@ generate_new_connectivity( return build_connectivity_by_hexahedron(ordered_x, ordered_y, ordered_z, block_number); } + +// NEW FUNCTIONS FOR NEW UPCOMING VERSION OF EXTEND CONNECTIVITY + +template +std::pair, std::array> +compute_index_and_refinement_for_element(const std::string& element_grid_name) { + // Computes the refinements and indieces for an element + // when given the grid names for that particular element. This function CANNOT + // be given all the gridnames across the whole domain, only for the one + // element. Returns a std::pair where the first entry contains the all the + // indices for the element and the second entry contains the refinement in + // every dimension of the element. + + std::array indices_of_element = {}; + std::array h_ref_of_element = {}; + + // some string indexing gymnastics + size_t grid_points_previous_start_position = 0; + size_t grid_points_start_position = 0; + size_t grid_points_end_position = 0; + + size_t h_ref_previous_start_position = 0; + size_t h_ref_start_position = 0; + size_t h_ref_end_position = 0; + + for (size_t j = 0; j < SpatialDim; ++j) { + grid_points_start_position = + element_grid_name.find('I', grid_points_previous_start_position + 1); + if (j == SpatialDim - 1) { + grid_points_end_position = + element_grid_name.find(')', grid_points_start_position); + } else { + grid_points_end_position = + element_grid_name.find(',', grid_points_start_position); + } + + std::stringstream element_index_substring(element_grid_name.substr( + grid_points_start_position + 1, + grid_points_end_position - grid_points_start_position - 1)); + size_t current_element_index = 0; + element_index_substring >> current_element_index; + gsl::at(indices_of_element, j) = current_element_index; + grid_points_previous_start_position = grid_points_start_position; + + h_ref_start_position = + element_grid_name.find('L', h_ref_previous_start_position + 1); + h_ref_end_position = element_grid_name.find('I', h_ref_start_position); + std::stringstream element_grid_name_substring(element_grid_name.substr( + h_ref_start_position + 1, + h_ref_end_position - h_ref_start_position - 1)); + size_t current_element_h_ref = 0; + element_grid_name_substring >> current_element_h_ref; + gsl::at(h_ref_of_element, j) = current_element_h_ref; + h_ref_previous_start_position = h_ref_start_position; + } + + // pushes back the all indices for an element to the "indices" vector. + return std::pair{indices_of_element, h_ref_of_element}; +} + +// std::pair>, +// std::vector>> +template +std::vector< + std::pair, std::array>> +compute_indices_and_refinements_for_elements( + const std::vector& block_grid_names) { + // Computes the refinements and indieces for all the elements in a given block + // when given the grid names for that particular block. This function CANNOT + // be given all the gridnames across the whole domain, only for the one block. + // Returns a std::pair where the first entry contains the all the indices for + // all elements in a block and the second entry contains all the refinements + // in every dimension of every element within the block. + + std::vector< + std::pair, std::array>> + indices_and_h_refs = {}; + + // std::vector> indices = {}; + // std::vector> h_ref = {}; + + for (const auto& element_grid_name : block_grid_names) { + // pushes back the all indices for an element to the "indices" vector. + indices_and_h_refs.push_back( + compute_index_and_refinement_for_element( + element_grid_name)); + } + + return indices_and_h_refs; +} + +template +std::vector> compute_block_segment_ids( + const std::vector, + std::array>>& + indices_and_refinements_for_elements) { + // Creates a std::vector of the segmentIds of each element in the block that + // is the same block as the block for which the refinement and indices are + // calculated. + std::vector> block_segment_ids = {}; + std::array segment_ids_of_current_element{}; + for (size_t i = 0; i < indices_and_refinements_for_elements.size(); ++i) { + for (size_t j = 0; j < SpatialDim; ++j) { + SegmentId current_segment_id( + indices_and_refinements_for_elements[i].second[j], + indices_and_refinements_for_elements[i].first[j]); + gsl::at(segment_ids_of_current_element, j) = current_segment_id; + } + block_segment_ids.push_back(segment_ids_of_current_element); + } + + return block_segment_ids; +} + +template +std::vector> element_logical_coordinates( + const Mesh& element_mesh) { + // Computes the ELCs of one element. Only needs to take in the element mesh of + // that one element + + std::array grid_point_coordinates{}; + std::vector> element_logical_coordinates = {}; + const auto element_logical_coordinates_tensor = + logical_coordinates(element_mesh); + for (size_t i = 0; i < element_logical_coordinates_tensor.get(0).size(); + ++i) { + for (size_t j = 0; j < SpatialDim; ++j) { + gsl::at(grid_point_coordinates, j) = + element_logical_coordinates_tensor.get(j)[i]; + } + + element_logical_coordinates.push_back(grid_point_coordinates); + } + + return element_logical_coordinates; +} + +bool share_endpoints(const SegmentId& segment_id_1, + const SegmentId& segment_id_2) { + // Returns true if segment_id_1 and segment_id_2 touch in any way (including + // overlapping). Otherwise returns false + const double upper_1 = segment_id_1.endpoint(Side::Upper); + const double lower_1 = segment_id_1.endpoint(Side::Lower); + const double upper_2 = segment_id_2.endpoint(Side::Upper); + const double lower_2 = segment_id_2.endpoint(Side::Lower); + + if (upper_1 == upper_2 || upper_1 == lower_2 || lower_1 == upper_2 || + lower_1 == lower_2) { + return true; + } + + return false; +} + +template +std::vector> find_neighbors( + const std::array& element_of_interest, + std::vector>& all_elements) { + // Identifies all neighbors(face, edge, and corner) in one big vector. + // This does not differentiate between types of neighbors, just that they are + // neighbors. Takes in the element of interest and the rest of the elements in + // the block in the form of SegmentIds. This function also identifies the + // element of interest as it's own neighbor. This is removed in + // "neighbor_direction". + + for (size_t i = 0; i < SpatialDim; ++i) { + const SegmentId current_segment_id = gsl::at(element_of_interest, i); + // std::cout << "size: " << all_elements.size() << '\n'; + // lambda identifies the indicies of the non neighbors and stores it in + // not_neighbors + const auto not_neighbors = std::remove_if( + all_elements.begin(), all_elements.end(), + [&i, ¤t_segment_id]( + std::array element_to_compare) { + // only need touches since if they overlap, they also touch + bool touches = share_endpoints(current_segment_id, + gsl::at(element_to_compare, i)); + // Filters out neighbors with higher refinement + if (current_segment_id.refinement_level() > + gsl::at(element_to_compare, i).refinement_level()) { + touches = false; + } + return !touches; + }); + // removes all std::array that aren't neighbors + all_elements.erase(not_neighbors, all_elements.end()); + } + + // std::cout << "size: " << all_elements.size() << '\n'; + return all_elements; +} + +template +std::array neighbor_direction( + const std::array& element_of_interest, + const std::array& element_to_compare) { + // identifies the direction vector of the element to compare relative to the + // element of interest. It is known that element_to_compare is a neighbor + + std::array direction{}; + for (size_t i = 0; i < SpatialDim; ++i) { + if (gsl::at(element_of_interest, i).endpoint(Side::Upper) == + gsl::at(element_to_compare, i).endpoint(Side::Lower)) { + gsl::at(direction, i) = 1; + } + if (gsl::at(element_of_interest, i).endpoint(Side::Lower) == + gsl::at(element_to_compare, i).endpoint(Side::Upper)) { + gsl::at(direction, i) = -1; + } + if (gsl::at(element_of_interest, i).endpoint(Side::Upper) == + gsl::at(element_to_compare, i).endpoint(Side::Upper) || + gsl::at(element_of_interest, i).endpoint(Side::Lower) == + gsl::at(element_to_compare, i).endpoint(Side::Lower)) { + gsl::at(direction, i) = 0; + } + } + + return direction; +} + +template +std::vector< + std::pair, std::array>> +compute_neighbors_with_direction( + const std::array& element_of_interest, + const std::vector>& all_neighbors) { + // outputs a std::array of length SpatialDim of element segmentIds. First are + // face neighbors, then edge neighbors, then corner neighbors. Then it also + // returns a std::array of the direction vector of each neighbor as well in a + // std::pair. Takes in the element of interest and the list of all neighbors. + std::vector< + std::pair, std::array>> + neighbors_with_direction = {}; + + for (const auto& neighbor : all_neighbors) { + std::array direction = + neighbor_direction(element_of_interest, neighbor); + std::pair, std::array> + neighbor_with_direction{neighbor, direction}; + int direction_magnitude = 0; + + // Why are we doing this loop for every neighbour just to be able to + // eliminate itself + for (size_t j = 0; j < SpatialDim; ++j) { + direction_magnitude += abs(gsl::at(direction, j)); + } + // Filter out itself fromt the neighbors + if (direction_magnitude != 0) { + neighbors_with_direction.push_back(neighbor_with_direction); + } + } + + return neighbors_with_direction; +} + +template +std::vector> +block_logical_coordinates_for_element( + const std::vector>& + element_logical_coordinates, + const std::pair, + std::array>& + element_indices_and_refinements, + const bool& sort_flag) { + // Generates the BLC of a particular element given that elements ELC and it's + // indices and refinements as defined from its grid name. + + std::vector> BLCs_for_element = {}; + for (const auto& grid_point : element_logical_coordinates) { + std::array BLCs_of_gridpoint{}; + for (size_t j = 0; j < SpatialDim; ++j) { + int number_of_elements = + two_to_the(gsl::at(element_indices_and_refinements.second, j)); + double shift = -1 + (2 * static_cast(gsl::at( + element_indices_and_refinements.first, j)) + + 1) / + static_cast(number_of_elements); + gsl::at(BLCs_of_gridpoint, j) = + gsl::at(grid_point, j) / number_of_elements + shift; + } + BLCs_for_element.push_back(BLCs_of_gridpoint); + } + + // sort the BLC by increasing z, then y, then x + if (sort_flag) { + std::sort(BLCs_for_element.begin(), BLCs_for_element.end(), + [](const auto& lhs, const auto& rhs) { + return std::lexicographical_compare(lhs.begin(), lhs.end(), + rhs.begin(), rhs.end()); + }); + } + + return BLCs_for_element; +} + +template +std::string grid_name_reconstruction( + const std::array& element, + const std::vector& block_grid_names) { + // Gets the block number. This can be CHANGED to be passed in as an argument + // later. + std::string element_grid_name = + block_grid_names[0].substr(0, block_grid_names[0].find('(', 0) + 1); + for (size_t i = 0; i < SpatialDim; ++i) { + element_grid_name += + "L" + std::to_string(gsl::at(element, i).refinement_level()) + "I" + + std::to_string(gsl::at(element, i).index()); + if (i < SpatialDim - 1) { + element_grid_name += ","; + } else if (i == SpatialDim - 1) { + element_grid_name += ")]"; + } + } + return element_grid_name; +} + +template +std::vector> compute_element_BLCs( + const std::array& element, + const std::vector& block_grid_names, + const std::vector>& block_extents, + const std::vector>& block_bases, + const std::vector>& block_quadratures, + const std::vector, + std::array>>& + indices_and_refinements_for_elements) { + // grid_name_reconstruction requires block_grid_names for the block number. + // This should be CHANGED!!! later when we loop over the blocks to take in + // the index of the block. + + const std::string element_grid_name = + grid_name_reconstruction(element, block_grid_names); + // Find the index of the element of interest within grid_names so we can + // find it later in extents, bases, and quadratures + const auto found_element_grid_name = + alg::find(block_grid_names, element_grid_name); + const auto element_index = static_cast( + std::distance(block_grid_names.begin(), found_element_grid_name)); + + // Construct the mesh for the element of interest. mesh_for_grid finds the + // index internally for us. + const Mesh element_mesh = h5::mesh_for_grid( + element_grid_name, block_grid_names, block_extents, block_bases, + block_quadratures); + + // Compute the element logical coordinates for the element of interest + std::vector> element_ELCs = + element_logical_coordinates(element_mesh); + + // Access the indices and refinements for the element of interest and change + // container type. Was orginally a std::pair, + // std::vector<...>>. We index into the vectors and reconstruct the pair. + const std::pair, + std::array>& + element_indices_and_refinements{ + indices_and_refinements_for_elements[element_index].first, + indices_and_refinements_for_elements[element_index].second}; + + // Compute BLC for the element of interest + const std::vector>& element_BLCs = + block_logical_coordinates_for_element( + element_ELCs, element_indices_and_refinements, true); + + return element_BLCs; +} + +// Need offsets to determine how to index element for gridpoints +template +std::pair gridpoints_BLCs_dim_offsets( + const std::vector>& + element_gridpoints_BLCs) { + // Comparison function to check if BLC at index dimension is equal + auto is_equal = [element_gridpoints_BLCs]( + std::array gridpoint_BLCs, + size_t index) { + return gridpoint_BLCs[index] == element_gridpoints_BLCs[0][index]; + }; + + // Find first value at which is_equal is not true in both dimensions + // Add check that iterator did not return last, i.e. couldn't determine offset + size_t y_offset_index = std::distance( + element_gridpoints_BLCs.begin(), + std::find_if_not(element_gridpoints_BLCs.begin(), + element_gridpoints_BLCs.end(), + std::bind(is_equal, std::placeholders::_1, 1))); + size_t x_offset_index = std::distance( + element_gridpoints_BLCs.begin(), + std::find_if_not(element_gridpoints_BLCs.begin(), + element_gridpoints_BLCs.end(), + std::bind(is_equal, std::placeholders::_1, 0))); + + return std::make_pair(x_offset_index, y_offset_index); +} + +// Returns the BLCs and directions of all the neighbors in +// neighbors_with_direction +template +std::vector>, + std::array>> +compute_neighbor_BLCs_and_directions( + const std::vector& block_grid_names, + const std::vector>& block_extents, + const std::vector>& block_bases, + const std::vector>& block_quadratures, + const std::vector, + std::array>>& + neighbors_with_direction, + const std::vector, + std::array>>& + indices_and_refinements_for_elements) { + // Will store all neighbor BLC's and their direction vectors + std::vector>, + std::array>> + neighbor_BLCs_and_directions = {}; + // Need to loop over all the neighbors. First by neighbor type then the + // number of neighbor in each type. + for (size_t k = 0; k < neighbors_with_direction.size(); ++k) { + // Reconstruct the grid name for the neighboring element + // std::cout << "Neighbor grid name: " << '\n'; + // std::cout << "Neighbor BLCs: " << '\n'; + const std::vector> neighbor_BLCs = + compute_element_BLCs(neighbors_with_direction[k].first, + block_grid_names, block_extents, block_bases, + block_quadratures, + indices_and_refinements_for_elements); + + // Access the direction vector. Lives inside neighbors_with_direction + // datastructure + const std::array neighbor_direction = + gsl::at(neighbors_with_direction, k).second; + // std::cout << "direction vector: " << + // neighbors_with_direction[k].second[0] + // << ", " << neighbors_with_direction[k].second[1] << ", " + // << neighbors_with_direction[k].second[2] << '\n'; + + std::pair>, + std::array> + neighbor_BLCs_and_direction{neighbor_BLCs, neighbor_direction}; + neighbor_BLCs_and_directions.push_back(neighbor_BLCs_and_direction); + } + + return neighbor_BLCs_and_directions; +} + +template +std::vector>, + std::array>> +compute_neighbor_info( + const std::array& element_of_interest, + std::vector>& neighbor_segment_ids, + const std::vector& block_grid_names, + const std::vector>& block_extents, + const std::vector>& block_bases, + const std::vector>& block_quadratures, + const std::vector, + std::array>>& + indices_and_refinements_for_elements) { + const std::vector> all_neighbors = + find_neighbors(element_of_interest, neighbor_segment_ids); + + // Gives all neighbors sorted by type. First is face, then edge, then + // corner. Also gives the direction vector in the pair. + const std::vector< + std::pair, std::array>> + neighbors_with_direction = + compute_neighbors_with_direction(element_of_interest, all_neighbors); + + // Stores all neighbor BLCs and directions in neighbor_BLCs_and_directions + const std::vector>, + std::array>> + neighbor_BLCs_and_directions = compute_neighbor_BLCs_and_directions( + block_grid_names, block_extents, block_bases, block_quadratures, + neighbors_with_direction, indices_and_refinements_for_elements); + + return neighbor_BLCs_and_directions; +} } // namespace namespace h5::detail { @@ -482,13 +959,182 @@ std::vector extend_connectivity( return new_connectivity; } +// NEW MAIN EXTEND CONNECTIVITY FUNCTIONS + +// Write new connectivity connections given a std::vector of observation ids +template +std::vector> extend_connectivity_by_block( + const std::vector& block_grid_names, + const std::vector>& block_extents, + const std::vector>& block_bases, + const std::vector>& block_quadratures) { + std::vector> block_connectivity; + + // std::pair of the indices(first entry) and refinements (second entry) for + // the entire block. + // const std::pair>, + // std::vector>>& + const std::vector, + std::array>>& + indices_and_refinements_for_elements = + compute_indices_and_refinements_for_elements( + block_grid_names); + + // Segment Ids for every element in the block. Each element has SpatialDim + // SegmentIds, one for each dimension. + const std::vector>& block_segment_ids = + compute_block_segment_ids( + indices_and_refinements_for_elements); + + for (size_t i = 0; i < block_segment_ids.size(); ++i) { + // Need to copy block_segment_ids to a new container since it will be + // altered. + std::vector> neighbor_segment_ids = + block_segment_ids; + // Identify the element I want to find the neighbors of. + const std::array& element_of_interest = + gsl::at(block_segment_ids, i); + + // Need the grid name of the element of interest to reverse search it in the + // vector of all grid names to get its position in that vector + // std::cout << "Element of interest grid name: " << '\n'; + // std::cout << "Element of interest BLCs: " << '\n'; + std::vector> element_of_interest_BLCs = + compute_element_BLCs(element_of_interest, block_grid_names, + block_extents, block_bases, block_quadratures, + indices_and_refinements_for_elements); + + std::pair elem_of_int_dim_offsets = + gridpoints_BLCs_dim_offsets(element_of_interest_BLCs); + + // Stores all neighbor BLCs and directions in neighbor_info + std::vector>, + std::array>> + neighbour_info = compute_neighbor_info( + element_of_interest, neighbor_segment_ids, block_grid_names, + block_extents, block_bases, block_quadratures, + indices_and_refinements_for_elements); + + /* To be added: Determining the connectivity between the element of + interest and the neighbour + */ + } + + return block_connectivity; +} + +/* Write new connectivity connections given the grid_names, bases, quadratures, + and extents +*/ +template +std::vector new_extend_connectivity( + std::vector& grid_names, + std::vector>& bases, + std::vector>& quadratures, + std::vector>& extents) { + auto [number_of_blocks, block_number_for_each_element, + sorted_element_indices] = compute_and_organize_block_info(grid_names); + + const auto sorted_grid_names = + sort_by_block(sorted_element_indices, grid_names); + const auto sorted_extents = sort_by_block(sorted_element_indices, extents); + const auto sorted_bases = sort_by_block(sorted_element_indices, bases); + const auto sorted_quadratures = + sort_by_block(sorted_element_indices, quadratures); + + // Create an unordered_map to be used to associate a grid point's block + // number and coordinates as an array to the its label + // (B#, grid_point_coord_array) -> grid_point_number + std::unordered_map< + std::pair>, size_t, + boost::hash>>> + block_and_grid_point_map; + + // Counter for the grid points when filling the unordered_map. Grid points + // are labelled by positive integers, so we are numbering them with this + // counter as we associate them to the key (B#, grid_point_coord_array) in + // the unordered map. + size_t grid_point_number = 0; + + for (size_t element_index = 0; + element_index < block_number_for_each_element.size(); ++element_index) { + auto element_mesh = mesh_for_grid( + grid_names[element_index], grid_names, extents, bases, quadratures); + auto element_logical_coordinates_tensor = logical_coordinates(element_mesh); + + std::vector> element_logical_coordinates; + element_logical_coordinates.reserve( + element_logical_coordinates_tensor.get(0).size()); + + for (size_t k = 0; k < element_logical_coordinates_tensor.get(0).size(); + ++k) { + std::array logical_coords_element_increment = {}; + for (size_t l = 0; l < SpatialDim; l++) { + gsl::at(logical_coords_element_increment, l) = + element_logical_coordinates_tensor.get(l)[k]; + } + element_logical_coordinates.push_back(logical_coords_element_increment); + } + + const std::pair, + std::array>& element_inds_and_refs = + compute_index_and_refinement_for_element( + grid_names[element_index]); + + std::vector> block_logical_coordinates = + block_logical_coordinates_for_element( + element_logical_coordinates, element_inds_and_refs, false); + + // Stores (B#, grid_point_coord_array) -> grid_point_number in an + // unordered_map and grid_point_coord_array by block + for (size_t k = 0; k < block_logical_coordinates.size(); ++k) { + std::pair> block_and_grid_point( + block_number_for_each_element[element_index], + block_logical_coordinates[k]); + block_and_grid_point_map.insert( + std::pair>, size_t>( + block_and_grid_point, grid_point_number)); + grid_point_number += 1; + } + } + + std::vector new_connectivity; + + /* Formula for extended connectivity size assumes face, edge, and corner + neighbours are being connected. Since only face neighbours is complete as + of now, reserving memory for the connectivity vector is skipped for now. + */ + + // size_t total_expected_connectivity = 0; + // new_connectivity.reserve(total_expected_connectivity); + + for (size_t j = 0; j < sorted_element_indices.size(); ++j) { + std::vector> block_connectivity = + extend_connectivity_by_block( + sorted_grid_names[j], sorted_extents[j], sorted_bases[j], + sorted_quadratures[j]); + + for (const std::array& it : block_connectivity) { + new_connectivity.push_back( + block_and_grid_point_map[std::make_pair(j, it)]); + } + } + + return new_connectivity; +} + #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template std::vector h5::detail::extend_connectivity( \ - std::vector & grid_names, \ - std::vector> & bases, \ - std::vector> & quadratures, \ +#define INSTANTIATE(_, data) \ + template std::vector h5::detail::extend_connectivity( \ + std::vector & grid_names, \ + std::vector> & bases, \ + std::vector> & quadratures, \ + std::vector> & extents); \ + template std::vector h5::detail::new_extend_connectivity( \ + std::vector & grid_names, \ + std::vector> & bases, \ + std::vector> & quadratures, \ std::vector> & extents); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) diff --git a/src/IO/H5/ExtendConnectivityHelpers.hpp b/src/IO/H5/ExtendConnectivityHelpers.hpp index dc8b62dc8c7d..7cd6170dbf23 100644 --- a/src/IO/H5/ExtendConnectivityHelpers.hpp +++ b/src/IO/H5/ExtendConnectivityHelpers.hpp @@ -21,6 +21,19 @@ std::vector extend_connectivity( std::vector>& bases, std::vector>& quadratures, std::vector>& extents); -} +template +std::vector> extend_connectivity_by_block( + const std::vector& block_grid_names, + const std::vector>& block_extents, + const std::vector>& block_bases, + const std::vector>& block_quadratures); + +template +std::vector new_extend_connectivity( + std::vector& grid_names, + std::vector>& bases, + std::vector>& quadratures, + std::vector>& extents); +} /// \endcond