diff --git a/examples/cpp/ex_aniso_lines.cpp b/examples/cpp/ex_aniso_lines.cpp index 2608ead..c71cfc5 100644 --- a/examples/cpp/ex_aniso_lines.cpp +++ b/examples/cpp/ex_aniso_lines.cpp @@ -94,6 +94,9 @@ constexpr comma::CoMMASeedsPoolT seed_order = comma::CoMMASeedsPoolT::BOUNDARY_PRIORITY; constexpr IntT sing_card = 1; constexpr std::optional max_cells_in_line = std::nullopt; +constexpr comma::CoMMACellCouplingT aniso_cell_coupling = + comma::CoMMACellCouplingT::MAX_WEIGHT; +constexpr bool force_line_direction = true; constexpr IntT fc_iter = 1; constexpr comma::CoMMANeighbourhoodT neigh = comma::CoMMANeighbourhoodT::EXTENDED; @@ -144,6 +147,8 @@ void ex_agglomerate_one_level() { AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh ); @@ -188,7 +193,9 @@ void ex_agglomerate_one_level_args() { build_lines, odd_line_length, aniso_thresh, - max_cells_in_line + max_cells_in_line, + aniso_cell_coupling, + force_line_direction ); // Prepare output containers diff --git a/examples/scripts/comma_tools.py b/examples/scripts/comma_tools.py index 62436ab..ac36448 100644 --- a/examples/scripts/comma_tools.py +++ b/examples/scripts/comma_tools.py @@ -13,7 +13,7 @@ import numpy as np import numpy.typing as npt -from CoMMA import AR, Neighbourhood, SeedsPool +from CoMMA import AR, CellCoupling, Neighbourhood, SeedsPool AR_DESCRIPTIONS = { AR.DIAMETER_OVER_RADIUS: "Diameter over radius", @@ -37,6 +37,10 @@ SeedsPool.BOUNDARY_POINT_INIT: "Boundary priority with point initialization", # noqa: E501 SeedsPool.NEIGHBOURHOOD_POINT_INIT: "Neighbourhood priority with point initialization", # noqa: E501 } +CELL_COUPLING_DESCRIPTIONS = { + CellCoupling.MAX_WEIGHT: "Maximum edge weight", + CellCoupling.MIN_DISTANCE: "Minimum distance between cell centers", +} def compute_neighbourhood_wall_distance( diff --git a/examples/scripts/ex_aniso_lines.py b/examples/scripts/ex_aniso_lines.py index 4f1e2a7..d8ab1e7 100644 --- a/examples/scripts/ex_aniso_lines.py +++ b/examples/scripts/ex_aniso_lines.py @@ -23,6 +23,7 @@ from comma_tools import ( AR_DESCRIPTIONS, + CELL_COUPLING_DESCRIPTIONS, NEIGHBOURHOOD_DESCRIPTIONS, SEED_ORDERING_DESCRIPTIONS, assign_anisotropic_line_data_to_cells, @@ -71,6 +72,12 @@ sing_card = 1 # Max cells in an anisotropic line max_cells_in_line = None # Or positive number +# Anisotropic cells coupling. Choices: +# - CellCoupling.MAX_WEIGHT = max edge weight (i.e., area) +# - CellCoupling.MIN_DISTANCE = minimum distance between cell centers +aniso_cell_coupling = CoMMA.CellCoupling.MAX_WEIGHT +# Whether to force aniso lines to have a straight direction +force_line_direction = True # Number of iterations for iterative fine-cell research algorithm fc_iter = 1 @@ -103,6 +110,10 @@ print(f" * seed_ordering={SEED_ORDERING_DESCRIPTIONS[seed_order]}") print(f" * Threshold cardinality for singular cells={sing_card}") print(f" * Max cells in anisotropic line={max_cells_in_line}") +print( + f" * Cell coupling for aniso cells={CELL_COUPLING_DESCRIPTIONS[aniso_cell_coupling]}" +) +print(f" * Force aniso lines direction={force_line_direction}") print(f" * Fine-cell research iterations={fc_iter}") print(" [Output]") renum = renumber_coarse > 1 @@ -189,6 +200,8 @@ AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh_type, ) diff --git a/examples/scripts/ex_default_dualGPy.py b/examples/scripts/ex_default_dualGPy.py index 93b7802..6b7a58d 100644 --- a/examples/scripts/ex_default_dualGPy.py +++ b/examples/scripts/ex_default_dualGPy.py @@ -23,6 +23,7 @@ from comma_tools import ( AR_DESCRIPTIONS, + CELL_COUPLING_DESCRIPTIONS, NEIGHBOURHOOD_DESCRIPTIONS, SEED_ORDERING_DESCRIPTIONS, prepare_meshio_agglomeration_data, @@ -75,6 +76,12 @@ sing_card = 1 # Max cells in an anisotropic line max_cells_in_line = None # Or positive number +# Anisotropic cells coupling. Choices: +# - CellCoupling.MAX_WEIGHT = max edge weight (i.e., area) +# - CellCoupling.MIN_DISTANCE = minimum distance between cell centers +aniso_cell_coupling = CoMMA.CellCoupling.MAX_WEIGHT +# Whether to force aniso lines to have a straight direction +force_line_direction = True # Number of iterations for iterative fine-cell research algorithm fc_iter = 1 @@ -107,6 +114,10 @@ print(f" * seed_ordering={SEED_ORDERING_DESCRIPTIONS[seed_order]}") print(f" * Threshold cardinality for singular cells={sing_card}") print(f" * Max cells in anisotropic line={max_cells_in_line}") +print( + f" * Cell coupling for aniso cells={CELL_COUPLING_DESCRIPTIONS[aniso_cell_coupling]}" +) +print(f" * Force aniso lines direction={force_line_direction}") print(f" * Fine-cell research iterations={fc_iter}") print(" [Output]") renum = renumber_coarse > 1 @@ -181,6 +192,8 @@ AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh_type, ) diff --git a/examples/scripts/ex_rae.py b/examples/scripts/ex_rae.py index b0302c1..66d6214 100644 --- a/examples/scripts/ex_rae.py +++ b/examples/scripts/ex_rae.py @@ -24,6 +24,7 @@ from comma_tools import ( AR_DESCRIPTIONS, + CELL_COUPLING_DESCRIPTIONS, NEIGHBOURHOOD_DESCRIPTIONS, SEED_ORDERING_DESCRIPTIONS, assign_anisotropic_line_data_to_cells, @@ -91,6 +92,12 @@ def limit_line_length(idxs, cells, max_cells_in_line): sing_card = 1 # Max cells in an anisotropic line max_cells_in_line = None # Or positive number +# Anisotropic cells coupling. Choices: +# - CellCoupling.MAX_WEIGHT = max edge weight (i.e., area) +# - CellCoupling.MIN_DISTANCE = minimum distance between cell centers +aniso_cell_coupling = CoMMA.CellCoupling.MAX_WEIGHT +# Whether to force aniso lines to have a straight direction +force_line_direction = True # Number of iterations for iterative fine-cell research algorithm fc_iter = 1 # Number of iteration steps to perform @@ -129,6 +136,10 @@ def limit_line_length(idxs, cells, max_cells_in_line): print(f" * seed_ordering={SEED_ORDERING_DESCRIPTIONS[seed_order]}") print(f" * Threshold cardinality for singular cells={sing_card}") print(f" * Max cells in anisotropic line={max_cells_in_line}") +print( + f" * Cell coupling for aniso cells={CELL_COUPLING_DESCRIPTIONS[aniso_cell_coupling]}" +) +print(f" * Force aniso lines direction={force_line_direction}") print(f" * Fine-cell research iterations={fc_iter}") print(f" * {agglomeration_levels=}") print(" [Output]") @@ -280,6 +291,8 @@ def limit_line_length(idxs, cells, max_cells_in_line): AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh_type, ) diff --git a/examples/scripts/ex_read_mesh.py b/examples/scripts/ex_read_mesh.py index bd5dc7b..894e880 100644 --- a/examples/scripts/ex_read_mesh.py +++ b/examples/scripts/ex_read_mesh.py @@ -26,6 +26,7 @@ from comma_tools import ( AR_DESCRIPTIONS, + CELL_COUPLING_DESCRIPTIONS, NEIGHBOURHOOD_DESCRIPTIONS, SEED_ORDERING_DESCRIPTIONS, prepare_meshio_agglomeration_data, @@ -80,6 +81,12 @@ sing_card = 1 # Max cells in an anisotropic line max_cells_in_line = None # Or positive number +# Anisotropic cells coupling. Choices: +# - CellCoupling.MAX_WEIGHT = max edge weight (i.e., area) +# - CellCoupling.MIN_DISTANCE = minimum distance between cell centers +aniso_cell_coupling = CoMMA.CellCoupling.MAX_WEIGHT +# Whether to force aniso lines to have a straight direction +force_line_direction = True # Number of iterations for iterative fine-cell research algorithm fc_iter = 1 @@ -116,6 +123,10 @@ print(f" * seed_ordering={SEED_ORDERING_DESCRIPTIONS[seed_order]}") print(f" * Threshold cardinality for singular cells={sing_card}") print(f" * Max cells in anisotropic line={max_cells_in_line}") +print( + f" * Cell coupling for aniso cells={CELL_COUPLING_DESCRIPTIONS[aniso_cell_coupling]}" +) +print(f" * Force aniso lines direction={force_line_direction}") print(f" * Fine-cell research iterations={fc_iter}") print(" [Output]") renum = renumber_coarse > 1 @@ -190,6 +201,8 @@ AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh_type, ) diff --git a/examples/scripts/ex_several_agglo_levels.py b/examples/scripts/ex_several_agglo_levels.py index 1f5681a..cd43fc9 100644 --- a/examples/scripts/ex_several_agglo_levels.py +++ b/examples/scripts/ex_several_agglo_levels.py @@ -22,6 +22,7 @@ from comma_tools import ( AR_DESCRIPTIONS, + CELL_COUPLING_DESCRIPTIONS, NEIGHBOURHOOD_DESCRIPTIONS, SEED_ORDERING_DESCRIPTIONS, build_coarse_graph, @@ -75,6 +76,12 @@ sing_card = 1 # Max cells in an anisotropic line max_cells_in_line = None # Or positive number +# Anisotropic cells coupling. Choices: +# - CellCoupling.MAX_WEIGHT = max edge weight (i.e., area) +# - CellCoupling.MIN_DISTANCE = minimum distance between cell centers +aniso_cell_coupling = CoMMA.CellCoupling.MAX_WEIGHT +# Whether to force aniso lines to have a straight direction +force_line_direction = True # Number of iterations for iterative fine-cell research algorithm fc_iter = 1 agglomeration_levels = 3 @@ -109,6 +116,10 @@ print(f" * seed_ordering={SEED_ORDERING_DESCRIPTIONS[seed_order]}") print(f" * Threshold cardinality for singular cells={sing_card}") print(f" * Max cells in anisotropic line={max_cells_in_line}") +print( + f" * Cell coupling for aniso cells={CELL_COUPLING_DESCRIPTIONS[aniso_cell_coupling]}" +) +print(f" * Force aniso lines direction={force_line_direction}") print(f" * Fine-cell research iterations={fc_iter}") print(f" * {agglomeration_levels=}") print(" [Output]") @@ -216,6 +227,8 @@ AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh_type, ) diff --git a/examples/scripts/pyproject.toml b/examples/scripts/pyproject.toml index 9724636..35ac377 100644 --- a/examples/scripts/pyproject.toml +++ b/examples/scripts/pyproject.toml @@ -3,6 +3,8 @@ profile = "black" [tool.ruff] # Same as black line-length = 88 +# Disable line too long since Black is already on +extend-ignore = ["E501"] [tool.pytest.ini_options] addopts = [ "--ignore=deprecated", diff --git a/examples/scripts/test_CoMMA.py b/examples/scripts/test_CoMMA.py index 9f586b3..f9a19fc 100644 --- a/examples/scripts/test_CoMMA.py +++ b/examples/scripts/test_CoMMA.py @@ -27,6 +27,8 @@ def agglomerate_aniso_lines(): neigh_type = CoMMA.Neighbourhood.EXTENDED sing_card = 1 max_cells_in_line = None # Or positive number + aniso_cell_coupling = CoMMA.CellCoupling.MAX_WEIGHT + force_line_direction = True fc_iter = 1 # fmt: off adjMatrix_row_ptr = [ @@ -87,6 +89,8 @@ def agglomerate_aniso_lines(): AR, sing_card, max_cells_in_line, + aniso_cell_coupling, + force_line_direction, fc_iter, neigh_type, ) diff --git a/include/CoMMA/Agglomerator.h b/include/CoMMA/Agglomerator.h index d61c2d1..3e5bf64 100644 --- a/include/CoMMA/Agglomerator.h +++ b/include/CoMMA/Agglomerator.h @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -27,6 +28,7 @@ #include #include "CoMMA/ARComputer.h" +#include "CoMMA/CoMMADefs.h" #include "CoMMA/Coarse_Cell_Container.h" #include "CoMMA/Dual_Graph.h" #include "CoMMA/Neighbourhood.h" @@ -34,45 +36,6 @@ namespace comma { -/** @brief Given the features of a cell, compute its aspect-ratio. This uses - * the diameter and the measure.\n In 3D: - * \f$ AR = \frac{diam_{CC}}{\sqrt{vol_{CC}}} \f$ \n In 2D: - * \f$ AR = \frac{diam_{CC}}{\sqrt[3]{vol_{CC}}} \f$ (Recall that in 2D the - * volume is the surface) \n Generally, for dimension \f$ d >=2 \f$: - * \f$ AR = \frac{diam_{CC}}{\sqrt[d]{vol_{CC}}} \f$, otherwise no root: - * \f$ AR = \frac{diam_{CC}}{vol_{CC}} \f$ - * @tparam T type used for the features - * @tparam dim Dimension (e.g., 2D) - * @param[in] feat Cell features - * @return the aspect-ratio - */ -template -inline RealT compute_AR_max_ov_radius( - const CellFeatures &feat -) { - if constexpr (dim < 2) { - return sqrt(feat._sq_diam) / feat._measure; - } else if constexpr (dim == 2) { - return sqrt(feat._sq_diam) / sqrt(feat._measure); - } else if constexpr (dim == 3) { - return sqrt(feat._sq_diam) / cbrt(feat._measure); - } else { - return sqrt(feat._sq_diam) / pow(feat._measure, 1.0 / dim); - } -} - -/** @brief Given the features of a cell, compute its aspect-ratio. This returns - * the ratio between the diameter and the minimum edge. - * @tparam T type used for the features - * @param[in] feat Cell features - * @return the aspect-ratio - */ -template -inline RealT compute_AR_max_ov_min(const CellFeatures &feat -) { - return feat._diam / feat._min_edge; -} - // How to pass parameters from base class // https://stackoverflow.com/questions/9692675/c-constructor-where-parameters-are-used-by-base-class-constructor // https://stackoverflow.com/questions/120876/what-are-the-rules-for-calling-the-base-class-constructor @@ -245,6 +208,10 @@ class Agglomerator_Anisotropic : * @param[in] max_cells_in_line Maximum number of cells in an anisotropic * line; when this value is reached, all reaming cells are discarded, hence * considered as isotropic + * @param[in] cell_coupling CoMMACellCouplingT indicating the type of + * coupling to consider when building anisotropic lines + * @param[in] force_line_direction Whether a continuous direction of the + * line should be enforced when building it */ Agglomerator_Anisotropic( std::shared_ptr> @@ -261,7 +228,9 @@ class Agglomerator_Anisotropic : const std::vector &priority_weights, const bool build_lines, const bool odd_line_length, - const std::optional max_cells_in_line + const std::optional max_cells_in_line, + CoMMACellCouplingT cell_coupling = CoMMACellCouplingT::MAX_WEIGHT, + const bool force_line_direction = true ) : Agglomerator( graph, cc_graph, seeds_pool, dimension @@ -269,7 +238,9 @@ class Agglomerator_Anisotropic : _should_agglomerate(true), _aniso_neighbours(), _odd_line_length(odd_line_length), - _max_cells_in_line(max_cells_in_line) { + _max_cells_in_line(max_cells_in_line), + _cell_coupling(cell_coupling), + _force_line_direction(force_line_direction) { // for every defined level (1 by default), contains the number of cells // e.g. _l_nb_of_cells[0]= number of cells on finest level // _l_nb_of_cells[1]= number of cells on the first coarse level @@ -368,7 +339,9 @@ class Agglomerator_Anisotropic : // Here we have to consider a special case when we have an odd number // of cells: THIS IS FUNDAMENTAL FOR THE CONVERGENCE OF THE MULTIGRID // ALGORITHM - std::unordered_set s_fc = {*line_it, *(line_it + 1)}; + std::unordered_set s_fc = { + *line_it, *std::next(line_it) + }; n_cells += 2; // We update the neighbours. At this stage, we do not check if it is // or will be agglomerated since there will be a cleaning step after @@ -380,8 +353,8 @@ class Agglomerator_Anisotropic : ); this->_aniso_neighbours.insert( this->_aniso_neighbours.end(), - this->_fc_graph->neighbours_cbegin(*(line_it + 1)), - this->_fc_graph->neighbours_cend(*(line_it + 1)) + this->_fc_graph->neighbours_cbegin(*std::next(line_it)), + this->_fc_graph->neighbours_cend(*std::next(line_it)) ); if ( std::distance(line_it, end) == 3 @@ -390,12 +363,12 @@ class Agglomerator_Anisotropic : && n_cells + 1 == this->_max_cells_in_line.value())) { if (this->_odd_line_length) { // If only three cells left, agglomerate them - s_fc.insert(*(line_it + 2)); + s_fc.insert(*std::next(line_it, 2)); n_cells++; this->_aniso_neighbours.insert( this->_aniso_neighbours.end(), - this->_fc_graph->neighbours_cbegin(*(line_it + 2)), - this->_fc_graph->neighbours_cend(*(line_it + 2)) + this->_fc_graph->neighbours_cbegin(*std::next(line_it, 2)), + this->_fc_graph->neighbours_cend(*std::next(line_it, 2)) ); } line_it++; // Ensure to break the loop after current iteration @@ -571,18 +544,25 @@ class Agglomerator_Anisotropic : std::vector max_weights( this->_fc_graph->_number_of_cells, 0.0 ); + std::vector alt_weights{}; std::vector to_treat(this->_fc_graph->_number_of_cells, false); // Computation of the anisotropic cell, alias of the cells for which the // ratio between the face with maximum area and the face with minimum area // is more than a given threshold. this->_fc_graph->tag_anisotropic_cells( + alt_weights, max_weights, to_treat, aniso_seeds_pool, threshold_anisotropy, priority_weights, + this->_cell_coupling, 0 ); + std::vector &ref_weights = + this->_cell_coupling == CoMMACellCouplingT::MAX_WEIGHT + ? this->_fc_graph->_m_CRS_Values + : alt_weights; if (aniso_seeds_pool.empty()) { std::cout << "CoMMA - No anisotropic cell found. Skipping anisotropic " "agglomeration" @@ -624,40 +604,21 @@ class Agglomerator_Anisotropic : while (!end) { // for the seed (that is updated each time end!= true) we fill the // neighbours and the weights - auto n_it = this->_fc_graph->neighbours_cbegin(seed); - auto w_it = this->_fc_graph->weights_cbegin(seed); - // std::vector of the candidates to continue the line + // Candidates to continue the line CoMMASetOfPairType candidates; // If the line is long enough, we use the direction. Otherwise, we use // the weight. Putting a high-level if to reduce the branching inside // the loop over the neighbours. - if (empty_line) { - for (; n_it != this->_fc_graph->neighbours_cend(seed) - && w_it != this->_fc_graph->weights_cend(seed); - ++n_it, ++w_it) { - if (to_treat[*n_it] && *w_it > 0.90 * max_weights[seed]) { - candidates.emplace(*n_it, *w_it); - } - } // end for loop + if (!this->_force_line_direction || empty_line) { + this->find_cell_candidates_for_line_max_weight( + seed, to_treat, ref_weights, max_weights, candidates + ); } else { // If not an empty line, we check the direction, see // !dot_deviate below - for (; n_it != this->_fc_graph->neighbours_cend(seed) - && w_it != this->_fc_graph->weights_cend(seed); - ++n_it, ++w_it) { - if (to_treat[*n_it] && *w_it > 0.90 * max_weights[seed]) { - std::vector cur_dir(pts_dim); - get_direction( - this->_fc_graph->_centers[seed], - this->_fc_graph->_centers[*n_it], - cur_dir - ); - const auto dot = dot_product(prev_dir, cur_dir); - if (!dot_deviate(dot)) { - candidates.emplace(*n_it, fabs(dot)); - } - } - } // end for loop + this->find_cell_candidates_for_line_direction( + seed, to_treat, ref_weights, max_weights, prev_dir, true, candidates + ); } if (!candidates.empty()) { // Even if we have more than one candidate, we choose just one @@ -677,11 +638,8 @@ class Agglomerator_Anisotropic : } to_treat[seed] = false; empty_line = false; - get_direction( - this->_fc_graph->_centers[old_seed], - this->_fc_graph->_centers[seed], - prev_dir - ); + // Always new_seed - old_seed to be sure to keep the same direction + this->_fc_graph->get_center_direction(seed, old_seed, prev_dir); if (!primal_dir.has_value()) primal_dir = prev_dir; } // 0 candidate, we are at the end of the line or at the end of one @@ -690,26 +648,18 @@ class Agglomerator_Anisotropic : // Before giving up, let's try another thing: Doing the same things as // above with the only difference that we allow to move through any // faces and not only the maximum one but still checking for direction - if (!empty_line) { + if (this->_force_line_direction && !empty_line) { // If not an empty line, we check the direction, see is_parallel // below - for (auto it = this->_fc_graph->neighbours_cbegin(seed); - it != this->_fc_graph->neighbours_cend(seed); - ++it) { - if (to_treat[*it]) { - std::vector cur_dir(pts_dim); - get_direction( - this->_fc_graph->_centers[seed], - this->_fc_graph->_centers[*it], - cur_dir - ); - const auto dot = - dot_product(prev_dir, cur_dir); - if (!dot_deviate(dot)) { - candidates.emplace(*it, fabs(dot)); - } - } - } // end for loop + this->find_cell_candidates_for_line_direction( + seed, + to_treat, + ref_weights, + max_weights, + prev_dir, + false, + candidates + ); if (!candidates.empty()) { // We found one! Keep going! const auto old_seed = seed; @@ -721,11 +671,9 @@ class Agglomerator_Anisotropic : } to_treat[seed] = false; empty_line = false; - get_direction( - this->_fc_graph->_centers[old_seed], - this->_fc_graph->_centers[seed], - prev_dir - ); + // Always new_seed - old_seed to be sure to keep the same + // direction + this->_fc_graph->get_center_direction(seed, old_seed, prev_dir); if (!primal_dir.has_value()) primal_dir = prev_dir; } else { // If we have already looked at the other side of the line or if @@ -734,8 +682,13 @@ class Agglomerator_Anisotropic : if (!opposite_direction_check && seed != primal_seed) { seed = primal_seed; // The check seed != primal_seed should ensure that - // primal_dir != null - prev_dir = primal_dir.value(); + // primal_dir != null. + // Changing direction with - + for (auto xyz = decltype(prev_dir.size()){0}; + xyz < prev_dir.size(); + xyz++) { + prev_dir[xyz] = -primal_dir.value()[xyz]; + } opposite_direction_check = true; } else { end = true; @@ -750,7 +703,12 @@ class Agglomerator_Anisotropic : seed = primal_seed; // The check seed != primal_seed should ensure that // primal_dir != null - prev_dir = primal_dir.value(); + // Changing direction with - + for (auto xyz = decltype(prev_dir.size()){0}; + xyz < prev_dir.size(); + xyz++) { + prev_dir[xyz] = -primal_dir.value()[xyz]; + } opposite_direction_check = true; } else { end = true; @@ -774,6 +732,88 @@ class Agglomerator_Anisotropic : return true; } + /** @brief Tell whether a weight should be considered as high coupling between + * two cells. + * @param[in] weight Value to be examined + * @param[in] ref_weight Reference against which the value will be compared + * @return whether the given value is a mark of high-coupling + */ + inline bool is_high_coupling( + const CoMMAWeightType weight, const CoMMAWeightType ref_weight + ) { + return weight > 0.90 * ref_weight; + } + + /** @brief Find cells which are good candidates to be added to the anisotropic + * line. In order to identify the candidate, the weights are considered. + * @param[in] seed Last cell to be added to the line. + * @param[in] to_treat Vector of booleans telling whether a cell should be + * considered. + * @param[in] weights Weights of the graph. + * @param[in] max_weights Vector holding for each cell the maximum weight + * among its neighbours. + * @param[out] candidates Set of the candidates, which are pairs of cell ID + * and related priority. + */ + inline void find_cell_candidates_for_line_max_weight( + const CoMMAIndexType seed, + const std::vector &to_treat, + const std::vector &weights, + const std::vector &max_weights, + CoMMASetOfPairType &candidates + ) { + const auto offset = this->_fc_graph->_m_CRS_Row_Ptr[seed]; + auto n_it = std::next(this->_fc_graph->_m_CRS_Col_Ind.cbegin(), offset); + auto w_it = std::next(weights.cbegin(), offset); + for (; n_it != this->_fc_graph->neighbours_cend(seed); ++n_it, ++w_it) { + if (to_treat[*n_it] && this->is_high_coupling(*w_it, max_weights[seed])) { + candidates.emplace(*n_it, *w_it); + } + } // end for loop + } + + /** @brief Find cells which are good candidates to be added to the anisotropic + * line. In order to identify the candidate, the direction is compared to the + * last reported one; additionally, if requested, the weights can also be + * considered + * @param[in] seed Last cell to be added to the line. + * @param[in] to_treat Vector of booleans telling whether a cell should be + * considered. + * @param[in] weights Weights of the graph. + * @param[in] max_weights Vector holding for each cell the maximum weight + * among its neighbours. + * @param[in] prev_dir Previous direction. + * @param[in] check_weight Whether to check weight of the potential candidate + * @param[out] candidates Set of the candidates, which are pairs of cell ID + * and related priority. + */ + inline void find_cell_candidates_for_line_direction( + const CoMMAIndexType seed, + const std::vector &to_treat, + const std::vector &weights, + const std::vector &max_weights, + const std::vector &prev_dir, + const bool check_weight, + CoMMASetOfPairType &candidates + ) { + const auto offset = this->_fc_graph->_m_CRS_Row_Ptr[seed]; + auto n_it = std::next(this->_fc_graph->_m_CRS_Col_Ind.cbegin(), offset); + auto w_it = std::next(weights.cbegin(), offset); + const auto pts_dim = this->_fc_graph->_centers[0].size(); + for (; n_it != this->_fc_graph->neighbours_cend(seed); ++n_it, ++w_it) { + if (to_treat[*n_it] + && (!check_weight || this->is_high_coupling(*w_it, max_weights[seed]))) { + std::vector cur_dir(pts_dim); + // Always new_seed - old_seed to be sure to keep the same direction + this->_fc_graph->get_center_direction(*n_it, seed, cur_dir); + const auto dot = dot_product(prev_dir, cur_dir); + // If not > 0 we might change direction 180 degrees + if (dot > 0 && !dot_deviate(dot)) + candidates.emplace(*n_it, dot); + } + } // end for loop + } + /** @brief Neighbours of the anisotropic cells agglomerated. They are used to * update the seeds pool */ @@ -786,6 +826,14 @@ class Agglomerator_Anisotropic : * reached, all reaming cells are discarded, hence considered as isotropic */ std::optional _max_cells_in_line; + + /** Type of cell coupling to consider when building aniso lines */ + CoMMACellCouplingT _cell_coupling; + + /** Whether the force the direction of the anisotropic lines to remain + * straight + */ + bool _force_line_direction; }; /** @brief Agglomerator_Isotropic class is a child class of the Agglomerator diff --git a/include/CoMMA/Args.h b/include/CoMMA/Args.h index c65e093..8c247dc 100644 --- a/include/CoMMA/Args.h +++ b/include/CoMMA/Args.h @@ -197,6 +197,12 @@ class AnisotropicArgs { CoMMAWeightType threshold_anisotropy; /** @brief Maximum number of cells in an anisotropic line */ std::optional max_cells_in_line = std::nullopt; + /** @brief Type of coupling to consider when building lines */ + CoMMACellCouplingT cell_coupling; + /** @brief Whether to force the direction of the anisotropic lines to remain + * straight + */ + bool line_direction; /** @brief Constructor * @param[in] is_anisotropic Whether to consider an anisotropic agglomeration. @@ -210,6 +216,10 @@ class AnisotropicArgs { * cell is considered as anisotropic. Default: 4.0 * @param[in] max_cells_in_line Maximum number of cells in an anisotropic * line. Default: no limit. + * @param[in] cell_coupling Type of coupling to consider when building lines. + * Default: CoMMACellCouplingT::MAX_WEIGHT. + * @param[in] line_direction Whether to force the direction of the + * anisotropic lines to remain straight. Default: true. */ AnisotropicArgs( bool is_anisotropic, @@ -217,14 +227,18 @@ class AnisotropicArgs { bool build_lines = true, bool odd_line_length = true, CoMMAWeightType threshold_anisotropy = 4., - std::optional max_cells_in_line = std::nullopt + std::optional max_cells_in_line = std::nullopt, + CoMMACellCouplingT cell_coupling = CoMMACellCouplingT::MAX_WEIGHT, + bool line_direction = true ) : is_anisotropic(is_anisotropic), anisotropicCompliantCells(anisotropicCompliantCells), build_lines(build_lines), odd_line_length(odd_line_length), threshold_anisotropy(threshold_anisotropy), - max_cells_in_line(max_cells_in_line) {} + max_cells_in_line(max_cells_in_line), + cell_coupling(cell_coupling), + line_direction(line_direction) {} }; } // end namespace comma diff --git a/include/CoMMA/CoMMA.h b/include/CoMMA/CoMMA.h index 0a9e52f..0dec5e0 100644 --- a/include/CoMMA/CoMMA.h +++ b/include/CoMMA/CoMMA.h @@ -105,6 +105,10 @@ constexpr CoMMAIntT iter_agglo_max_iter = 4; * @param[in] max_cells_in_line [Optional] Maximum number of cells in an * anisotropic line; when this value is reached, all reaming cells are * discarded, hence considered as isotropic + * @param[in] aniso_cell_coupling CoMMACellCouplingT indicating the type of + * coupling to consider when building anisotropic lines + * @param[in] force_line_direction Whether to force the direction of the + * anisotropic lines to remain straight * @param[in] fc_choice_iter (optional, default=1) Number of iterations allowed * for the algorithm choosing which fine cell to add next. The cost grows * exponentially, hence use small values. @@ -166,6 +170,8 @@ void agglomerate_one_level( CoMMAAspectRatioT aspect_ratio = CoMMAAspectRatioT::DIAMETER_OVER_RADIUS, CoMMAIntType singular_card_thresh = 1, std::optional max_cells_in_line = std::nullopt, + CoMMACellCouplingT aniso_cell_coupling = CoMMACellCouplingT::MAX_WEIGHT, + bool force_line_direction = true, CoMMAIntType fc_choice_iter = 1, CoMMANeighbourhoodT neighbourhood_type = CoMMANeighbourhoodT::EXTENDED ) { @@ -231,6 +237,17 @@ void agglomerate_one_level( << std::endl; max_cells_in_line = std::nullopt; } + switch (aniso_cell_coupling) { + case CoMMACellCouplingT::MAX_WEIGHT: + case CoMMACellCouplingT::MIN_DISTANCE: + // OK + break; + default: + std::cout << "CoMMA - Warning: Unknown anisotropic cell coupling." + " Switching it to: max weight." + << std::endl; + aniso_cell_coupling = CoMMACellCouplingT::MAX_WEIGHT; + } } else { // Anisotropic lines are provided if (agglomerationLines_Idx.size() < 2 || agglomerationLines.empty()) { @@ -362,7 +379,9 @@ void agglomerate_one_level( priority_weights, build_anisotropic_lines, odd_line_length, - max_cells_in_line + max_cells_in_line, + aniso_cell_coupling, + force_line_direction ); // Agglomerate anisotropic cells only @@ -486,6 +505,8 @@ inline void agglomerate_one_level( agglo.aspect_ratio, agglo.singular_card_thresh, aniso.max_cells_in_line, + aniso.cell_coupling, + aniso.line_direction, agglo.fc_choice_iter, agglo.neighbourhood_type ); diff --git a/include/CoMMA/CoMMADefs.h b/include/CoMMA/CoMMADefs.h index 9254b6f..62a1495 100644 --- a/include/CoMMA/CoMMADefs.h +++ b/include/CoMMA/CoMMADefs.h @@ -41,17 +41,17 @@ enum CoMMANeighbourhoodT : CoMMAIntT { /** @brief Type of seeds pool ordering */ enum CoMMASeedsPoolT : CoMMAIntT { - /** The number of boundary faces has higher priority */ + /** The number of boundary faces has highest priority */ BOUNDARY_PRIORITY = 0, /** The neighbourhood has highest priority (neighbours of coarse cells have - * higher priority) + * priority) */ NEIGHBOURHOOD_PRIORITY = 1, - /** The number of boundary faces has higher priority, and initialize with one + /** The number of boundary faces has highest priority, and initialize with one * point only then let evolve */ BOUNDARY_PRIORITY_ONE_POINT_INIT = 10, - /** The neighbourhood has higher priority, and initialize with one point only + /** The neighbourhood has highest priority, and initialize with one point only * then let evolve */ NEIGHBOURHOOD_PRIORITY_ONE_POINT_INIT = 11 @@ -99,6 +99,14 @@ enum CoMMAAspectRatioT : CoMMAIntT { ALGEBRAIC_PERIMETER_OVER_MEASURE, }; +/** @brief Type of coupling between cells in an anisotropic line */ +enum CoMMACellCouplingT : CoMMAIntT { + /** Maximum edge-weight (i.e., max area) */ + MAX_WEIGHT = 0, + /** Minimum centers distance */ + MIN_DISTANCE = 1 +}; + } // end namespace comma #endif diff --git a/include/CoMMA/Dual_Graph.h b/include/CoMMA/Dual_Graph.h index addfab4..3aa4d1e 100644 --- a/include/CoMMA/Dual_Graph.h +++ b/include/CoMMA/Dual_Graph.h @@ -25,7 +25,9 @@ #include #include +#include "CoMMA/CoMMADefs.h" #include "CoMMA/Seeds_Pool.h" +#include "CoMMA/Util.h" namespace comma { @@ -180,7 +182,7 @@ class Graph { */ inline ContainerIndexConstIt neighbours_cbegin(const CoMMAIndexType &i_c ) const { - return _m_CRS_Col_Ind.cbegin() + _m_CRS_Row_Ptr[i_c]; + return std::next(_m_CRS_Col_Ind.cbegin(), _m_CRS_Row_Ptr[i_c]); } /** @brief Get constant pointer to the element following the last neighbour of @@ -190,7 +192,7 @@ class Graph { */ inline ContainerIndexConstIt neighbours_cend(const CoMMAIndexType &i_c ) const { - return _m_CRS_Col_Ind.cbegin() + _m_CRS_Row_Ptr[i_c + 1]; + return std::next(_m_CRS_Col_Ind.cbegin(), _m_CRS_Row_Ptr[i_c + 1]); } /** @brief Based on the area of the faces composing the cell given as an @@ -218,7 +220,7 @@ class Graph { */ inline ContainerWeightConstIt weights_cbegin(const CoMMAIndexType &i_c ) const { - return _m_CRS_Values.cbegin() + _m_CRS_Row_Ptr[i_c]; + return std::next(_m_CRS_Values.cbegin(), _m_CRS_Row_Ptr[i_c]); } /** @brief Get constant pointer to the element following the last neighbour of @@ -227,7 +229,7 @@ class Graph { * @return The pointer */ inline ContainerWeightConstIt weights_cend(const CoMMAIndexType &i_c) const { - return _m_CRS_Values.cbegin() + _m_CRS_Row_Ptr[i_c + 1]; + return std::next(_m_CRS_Values.cbegin(), _m_CRS_Row_Ptr[i_c + 1]); } /** @brief Check the connectivity of the graph. @@ -694,43 +696,101 @@ class Dual_Graph : public Graph { + this->estimated_boundary_weight(idx_c); } + /** @brief Compute the direction from vertex \p a to vertex \p b and store it + * as unit vector in \p dir. + * @param[in] a Starting point + * @param[in] b End point + * @param[out] dir Unit vector of the direction + * @return the distance from the two points (the norm used for the + * normalization) + * @note This is just a wrapper around get_direction. + */ + inline CoMMAWeightType get_center_direction( + const CoMMAIndexType &a, + const CoMMAIndexType &b, + std::vector &dir + ) { + return get_direction( + this->_centers[a], this->_centers[b], dir + ); + } + /** @brief Tag cells as anisotropic if their aspect-ratio is over a given * threshold and order them according to given priority - * @param[out] max_weights Array of the maximum weight: the biggest area of + * @param[out] alt_weights Alternative version of graph edge weights computed + * if necessary (according to \p cell_coupling) + * @param[out] max_weights Array of the maximum weight: the biggest area of * the faces composing the given fine cell - * @param[out] is_anisotropic Vector of length equal to the total number of + * @param[out] is_anisotropic Vector of length equal to the total number of * cell telling whether a cell is anisotropic - * @param[out] aniso_seeds_pool Container containing the anisotropic cells in + * @param[out] aniso_seeds_pool Container containing the anisotropic cells in * the order they should be considered when computing the lines - * @param[in] threshold_anisotropy Value of the aspect ratio above which a - * cell is considered anisotropic. If negative, all compliant cells are - * considered as anisotropic + * @param[in] threshold_anisotropy Value of the aspect ratio above which a + * cell is considered anisotropic. If negative, all compliant cells are + * considered as anisotropic * @param[in] priority_weights Weights used to set the order telling where to * start agglomerating. The higher the weight, the higher the priority - * @param[in] preserving if 0 does not hit only the BL prism to preserve the + * @param[in] cell_coupling CoMMACellCouplingT indicating the type of + * coupling to consider when building anisotropic lines + * @param[in] preserving if 0 does not hit only the BL prism to preserve the * boundary layer otherwise 2 for 2D or 3 for the 3D to preserve the BL only * in the anisotropic agglomeration */ void tag_anisotropic_cells( + std::vector &alt_weights, ContainerWeightType &max_weights, std::vector &is_anisotropic, std::deque &aniso_seeds_pool, const CoMMAWeightType threshold_anisotropy, const ContainerWeightType &priority_weights, + const CoMMACellCouplingT cell_coupling, const CoMMAIndexType preserving ) { + if (cell_coupling == CoMMACellCouplingT::MIN_DISTANCE) { + // We have to build the new weights. Preparing the storage. + alt_weights.resize(this->_m_CRS_Values.size()); + std::fill(alt_weights.begin(), alt_weights.end(), 0.0); + } CoMMASetOfPairType aniso_w_weights{}; if (threshold_anisotropy < 0) { - for (const CoMMAIndexType i_fc : _s_anisotropic_compliant_cells) { - aniso_w_weights.emplace(i_fc, priority_weights[i_fc]); - max_weights[i_fc] = - *(max_element(this->weights_cbegin(i_fc), this->weights_cend(i_fc))); - } + switch (cell_coupling) { + case CoMMACellCouplingT::MAX_WEIGHT: { + for (const CoMMAIndexType i_fc : _s_anisotropic_compliant_cells) { + aniso_w_weights.emplace(i_fc, priority_weights[i_fc]); + max_weights[i_fc] = *( + max_element(this->weights_cbegin(i_fc), this->weights_cend(i_fc)) + ); + } + break; + } + case CoMMACellCouplingT::MIN_DISTANCE: { + for (const CoMMAIndexType i_fc : _s_anisotropic_compliant_cells) { + aniso_w_weights.emplace(i_fc, priority_weights[i_fc]); + CoMMAWeightType max_ov_dist = 0.0; + const auto offset = this->_m_CRS_Row_Ptr[i_fc]; + auto n_it = std::next(this->_m_CRS_Col_Ind.cbegin(), offset); + auto w_it = std::next(alt_weights.begin(), offset); + for (; n_it != this->neighbours_cend(i_fc); ++n_it, ++w_it) { + const auto dist = 1. + / squared_euclidean_distance( + this->_centers[i_fc], this->_centers[*n_it] + ); + *w_it = dist; + if (dist > max_ov_dist) max_ov_dist = dist; + } // for neigh + max_weights[i_fc] = max_ov_dist; + } + break; + } + default: + break; + } // end switch } else { for (const CoMMAIndexType i_fc : _s_anisotropic_compliant_cells) { + // The max weight is always computed since it's needed for the AR + CoMMAWeightType max_weight = 0.0; CoMMAWeightType min_weight = std::numeric_limits::max(); - CoMMAWeightType max_weight = 0.0; // computation of min_weight, max_weight for the current cell // Process of every faces/Neighbours and compute for the current cell @@ -754,8 +814,6 @@ class Dual_Graph : public Graph { nb_neighbours--; } } - - max_weights[i_fc] = max_weight; // Compute the aspect-ratio and add cell to list if necessary const auto ar = _compute_AR(min_weight, max_weight); // Anisotropy criteria for the line Admissibility @@ -778,6 +836,33 @@ class Dual_Graph : public Graph { break; } // End switch } // End if ar + + // Updating weights according to coupling + switch (cell_coupling) { + case CoMMACellCouplingT::MAX_WEIGHT: { + max_weights[i_fc] = max_weight; + break; + } + case CoMMACellCouplingT::MIN_DISTANCE: { + const auto offset = this->_m_CRS_Row_Ptr[i_fc]; + auto n_it = std::next(this->_m_CRS_Col_Ind.cbegin(), offset); + auto w_it = std::next(alt_weights.begin(), offset); + CoMMAWeightType max_ov_dist = 0.0; + for (; n_it != this->neighbours_cend(i_fc); ++n_it, ++w_it) { + const auto dist = 1. + / squared_euclidean_distance( + this->_centers[i_fc], this->_centers[*n_it] + ); + *w_it = dist; + if (dist > max_ov_dist) max_ov_dist = dist; + } // for neigh + max_weights[i_fc] = max_ov_dist; + break; + } + default: + break; + } // end switch + } // End for compliant cells } // End if threshold diff --git a/src/CoMMA.cpp b/src/CoMMA.cpp index ee833dd..804e454 100644 --- a/src/CoMMA.cpp +++ b/src/CoMMA.cpp @@ -153,6 +153,22 @@ PYBIND11_MODULE(CoMMA, module_handle) { "Algebraic-like perimeter over measure, that is, external weights over cell weight" ) .export_values(); + py::enum_( + module_handle, + "CellCoupling", + "Type of coupling between cells considered when building anisotropic lines" + ) + .value( + "MAX_WEIGHT", + CoMMACellCouplingT::MAX_WEIGHT, + "Maximum edge weight (i.e., max area)" + ) + .value( + "MIN_DISTANCE", + CoMMACellCouplingT::MIN_DISTANCE, + "Minimum centers distance" + ) + .export_values(); // Main function module_handle.def( @@ -191,6 +207,8 @@ PYBIND11_MODULE(CoMMA, module_handle) { CoMMAIntT aspect_ratio, CoMMAIntT singular_card_thresh, optional max_cells_in_line, + CoMMAIntT aniso_cell_coupling, + bool force_line_direction, CoMMAIntT fc_choice_iter, CoMMAIntT type_of_isotropic_agglomeration ) { @@ -220,6 +238,8 @@ PYBIND11_MODULE(CoMMA, module_handle) { static_cast(aspect_ratio), singular_card_thresh, max_cells_in_line, + static_cast(aniso_cell_coupling), + force_line_direction, fc_choice_iter, static_cast(type_of_isotropic_agglomeration) ); @@ -250,6 +270,8 @@ PYBIND11_MODULE(CoMMA, module_handle) { "aspect_ratio"_a = CoMMAAspectRatioT::DIAMETER_OVER_RADIUS, "singular_card_thresh"_a = 1, "max_cells_in_line"_a = std::nullopt, + "aniso_cell_coupling"_a = CoMMACellCouplingT::MAX_WEIGHT, + "force_line_direction"_a = true, "fc_choice_iter"_a = 1, "type_of_isotropic_agglomeration"_a = CoMMANeighbourhoodT::EXTENDED ); diff --git a/tests/DualGraphExamples.h b/tests/DualGraphExamples.h index a463c0b..0f8f4f8 100644 --- a/tests/DualGraphExamples.h +++ b/tests/DualGraphExamples.h @@ -12,6 +12,7 @@ * Any copyright is dedicated to the Public Domain. * https://creativecommons.org/publicdomain/zero/1.0/ */ +#include #include #include "CoMMA/CoMMADefs.h" @@ -20,6 +21,11 @@ using DGEIndexT = comma::CoMMAIndexT; using DGEIntT = comma::CoMMAIntT; using DGEWeightT = comma::CoMMAWeightT; +constexpr DGEWeightT _1ov3 = 1. / 3.; +constexpr DGEWeightT _3_1ov3 = 10. / 3.; +const DGEWeightT sqrt2 = std::sqrt(2.); +const DGEWeightT sqrt101 = std::sqrt(101.); + class BaseDualGEx { public: /** @brief Dimension */ @@ -228,7 +234,7 @@ class DualGEx_hole_w_corners: public BaseDualGEx { // clang-format off /* row_ptr */ {0, 3, 6, 9, 12, 16, 20, 24, 28, 31, 34, 37, 40, 42, 44}, /* col_ind */ {4, 12, 13, 2, 5, 13, 1, 3, 6, 2, 7, 12, 0, 7, 5, 8, 1, 4, 6, 9, 2, 5, 7, 10, 3, 4, 6, 11, 4, 11, 9, 5, 8, 10, 6, 9, 11, 7, 8, 10, 0, 3, 0, 1}, - /* areas */ {6.0, 1.0, 1.0, 1.4142135623730951, 6.0, 1.0, 1.4142135623730951, 1.4142135623730951, 6.0, 1.4142135623730951, 6.0, 1.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 4.0, 1.4142135623730951, 1.4142135623730951, 4.0, 1.4142135623730951, 1.4142135623730951, 4.0, 1.4142135623730951, 1.4142135623730951, 4.0, 1.4142135623730951, 1.4142135623730951, 1.0, 1.0, 1.0, 1.0}, + /* areas */ {6.0, 1.0, 1.0, sqrt2, 6.0, 1.0, sqrt2, sqrt2, 6.0, sqrt2, 6.0, 1.0, 6.0, sqrt2, sqrt2, 4.0, 6.0, sqrt2, sqrt2, 4.0, 6.0, sqrt2, sqrt2, 4.0, 6.0, sqrt2, sqrt2, 4.0, 4.0, sqrt2, sqrt2, 4.0, sqrt2, sqrt2, 4.0, sqrt2, sqrt2, 4.0, sqrt2, sqrt2, 1.0, 1.0, 1.0, 1.0}, /* n_bnd */ {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2}, /* volumes */ {6.00, 6.50, 7.00, 6.50, 5.00, 5.00, 5.00, 5.00, 3.00, 3.00, 3.00, 3.00, 1.00, 1.00}, /* centers */ {{3.5, 0.0}, {-0.25, 3.5}, {-3.5, 0.0}, {-0.25, -3.5}, {2.5, 0.0}, {0.0, 2.5}, {-2.5, 0.0}, @@ -244,7 +250,7 @@ class DualGEx_hole_no_corners : public BaseDualGEx { // clang-format off /* row_ptr */ {0, 3, 6, 9, 12, 16, 20, 24, 28, 31, 34, 37, 40}, /* col_ind */ {3, 1, 4, 0, 2, 5, 1, 3, 6, 0, 2, 7, 0, 7, 5, 8, 1, 4, 6, 9, 2, 5, 7, 10, 3, 4, 6, 11, 4, 11, 9, 5, 8, 10, 6, 9, 11, 7, 8, 10}, - /* areas */ {1.4142135623730951, 1.4142135623730951, 6.0, 1.4142135623730951, 1.4142135623730951, 6.0, 1.4142135623730951, 1.4142135623730951, 6.0, 1.4142135623730951, 1.4142135623730951, 6.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 6.0, 1.4142135623730951, 1.4142135623730951, 4.0, 4.0, 1.4142135623730951, 1.4142135623730951, 4.0, 1.4142135623730951, 1.4142135623730951, 4.0, 1.4142135623730951, 1.4142135623730951, 4.0, 1.4142135623730951, 1.4142135623730951}, + /* areas */ {sqrt2, sqrt2, 6.0, sqrt2, sqrt2, 6.0, sqrt2, sqrt2, 6.0, sqrt2, sqrt2, 6.0, 6.0, sqrt2, sqrt2, 4.0, 6.0, sqrt2, sqrt2, 4.0, 6.0, sqrt2, sqrt2, 4.0, 6.0, sqrt2, sqrt2, 4.0, 4.0, sqrt2, sqrt2, 4.0, sqrt2, sqrt2, 4.0, sqrt2, sqrt2, 4.0, sqrt2, sqrt2}, /* n_bnd */ {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}, /* volumes */ {7.00, 7.00, 7.00, 7.00, 5.00, 5.00, 5.00, 5.00, 3.00, 3.00, 3.00, 3.00}, /* centers */ {{3.5, 0.0}, {0.0, 3.5}, {-3.5, 0.0}, {0.0, -3.5}, {2.5, 0.0}, {0.0, 2.5}, {-2.5, 0.0}, @@ -268,4 +274,51 @@ class DualGEx_T_shaped : public BaseDualGEx { ) {} }; +/** @brief 8 triangles obtained by splitting each rectangle of a 2x2 grid + * stretched (10x) in the x-direction. It aims to simulate a triangular + * boundary layer. + */ +#if 0 +Structure of the mesh, although in reality they are way more stretched in the +x-direction + _________ + | 5/| 7/| + | / | / | + |/_4|/_6| + | 1/| 3/| + | / | / | + |/_0|/_2| +#endif +class DualGEx_Tria_BL : public BaseDualGEx { +public: + DualGEx_Tria_BL() : BaseDualGEx( + // clang-format off + /* row_ptr */ {0, 2, 4, 5, 8, 11, 12, 14, 16}, + /* col_ind */ {1, 3, 0, 4, 3, 0, 2, 6, 1, 5, 7, 4, 3, 7, 4, 6}, + /* areas */ { + /* 0 */ sqrt101, 1.0, + /* 1 */ sqrt101, 10.0, + /* 2 */ sqrt101, + /* 3 */ 1.0, sqrt101, 10.0, + /* 4 */ 10.0, sqrt101, 1.0, + /* 5 */ sqrt101, + /* 6 */ 10.0, sqrt101, + /* 7 */ 1.0, sqrt101 + }, + /* n_bnd */ {1, 1, 2, 0, 0, 2, 1, 1}, + /* volumes */ {5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00, 5.00}, + /* centers */ { + // Bottom left, 0 and 1 + {-_3_1ov3, -2.*_1ov3}, {-2.*_3_1ov3, -_1ov3}, + // Bottom right, 2 and 3 + {+2.*_3_1ov3, -2*_1ov3}, {+_3_1ov3, -_1ov3}, + // Top left, 4 and 5 + {-_3_1ov3, +_1ov3}, {-2.*_3_1ov3, 2.*_1ov3}, + // Top right, 6 and 7 + {2.*_3_1ov3, _1ov3}, {_3_1ov3, 2.*_1ov3} + } + // clang-on + ) {} +}; + #endif diff --git a/tests/test_anisoagglo.cpp b/tests/test_anisoagglo.cpp index aebc46a..0e671ed 100644 --- a/tests/test_anisoagglo.cpp +++ b/tests/test_anisoagglo.cpp @@ -12,23 +12,27 @@ #include #include +#include +#include #include #include #include "CoMMA/Agglomerator.h" #include "CoMMA/CoMMA.h" +#include "CoMMA/CoMMADefs.h" +#include "CoMMA/Util.h" #include "DualGraphExamples.h" #include "test_defs.h" using namespace comma; // NOLINT using namespace std; // NOLINT +using Catch::Matchers::Approx; using Catch::Matchers::Contains; +using Catch::Matchers::Equals; +using Catch::Matchers::SizeIs; const auto EMPTY = Catch::Matchers::IsEmpty(); -#define SING_CARD_THRESH 1 -#define MAX_CELLS_IN_LINE std::nullopt - // For X-Macro where here X = CHECK2CELLS #define CHECK3CELLS(f2c, a, b, c) \ (CHECK2CELLS(f2c, a, b) && CHECK2CELLS(f2c, b, c)) @@ -71,7 +75,9 @@ SCENARIO( Data.weights, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); WHEN( @@ -121,7 +127,9 @@ SCENARIO( Data.weights, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); THEN( @@ -198,6 +206,7 @@ SCENARIO( ); shared_ptr const cc_graph = make_shared(fc_graph, SING_CARD_THRESH); + // Value does not matter since we already provide the lines const CoMMAWeightT aniso_thresh{10000.}; const bool build_lines = false; vector agglomerationLines_Idx{0, 1, 2}; @@ -273,8 +282,11 @@ SCENARIO( Data.weights, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); + THEN( "There is no need to agglomerate anisotropically since no line can be built" ) { @@ -326,7 +338,9 @@ SCENARIO( Data.weights, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); THEN("Only 3 lines are built") { REQUIRE(aniso_agg._nb_lines[0] == 3); @@ -371,7 +385,9 @@ SCENARIO( Data.weights, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); Agglomerator_Biconnected iso_agg( fc_graph, @@ -708,7 +724,9 @@ the line grows vertically wei, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); Agglomerator_Biconnected iso_agg( fc_graph, @@ -814,7 +832,9 @@ the line grows vertically Data.weights, build_lines, ODD_LINE_LENGTH, - MAX_CELLS_IN_LINE + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + FORCE_DIRECTION ); Agglomerator_Biconnected iso_agg( fc_graph, @@ -835,7 +855,7 @@ the line grows vertically ++ref_line, ++read_line) { // I am not sure how to treat this... // NOLINTNEXTLINE - for (auto ref_fc = agglomerationLines.cbegin() + (*(ref_line - 1)); + for (auto ref_fc = agglomerationLines.cbegin() + *std::prev(ref_line); ref_fc != agglomerationLines.cbegin() + (*ref_line); ++ref_fc) { REQUIRE_THAT(**read_line, Contains(*ref_fc)); @@ -1009,8 +1029,265 @@ SCENARIO("Test the anisotropic line computations", "[Anisotropic lines]") { } } } + +#undef CHECK2CELLS +#define CHECK2CELLS(f2c, a, b) ((f2c)[(a)].value() == (f2c)[(b)].value()) + + GIVEN("Triangular mesh stretched in the x-direction representing a BL") { + const DualGEx_Tria_BL Data = DualGEx_Tria_BL(); + shared_ptr const seeds_pool = + make_shared(Data.n_bnd_faces, Data.weights, false); + shared_ptr fc_graph = make_shared( + Data.nb_fc, + Data.adjMatrix_row_ptr, + Data.adjMatrix_col_ind, + Data.adjMatrix_areaValues, + Data.volumes, + Data.centers, + Data.n_bnd_faces, + Data.dim, + Data.anisoCompliantCells + ); + shared_ptr const cc_graph = + make_shared(fc_graph, SING_CARD_THRESH); + const CoMMAWeightT aniso_thresh{-2.}; + vector agglomerationLines_Idx{}; + vector agglomerationLines{}; + const bool build_lines = true; + WHEN("We agglomerate forcing the direction") { + constexpr bool force_direction = true; + Agglomerator_Anisotropic aniso_agg( + fc_graph, + cc_graph, + seeds_pool, + Data.dim, + aniso_thresh, + agglomerationLines_Idx, + agglomerationLines, + Data.weights, + build_lines, + ODD_LINE_LENGTH, + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + force_direction + ); + aniso_agg.agglomerate_one_level(2, 2, 2, Data.weights, false); + aniso_agg.export_anisotropic_lines( + 1, agglomerationLines_Idx, agglomerationLines + ); + THEN("Lines grow diagonally and the bottom right cell is singular") { + const auto &f2c = cc_graph->_fc_2_cc; + const vector ref_lines_idx = {0, 2, 3}; + const vector ref_lines = {0, 1, 2}; + REQUIRE_THAT(agglomerationLines_Idx, Equals(ref_lines_idx)); + REQUIRE_THAT(agglomerationLines, Equals(ref_lines)); + // Line 0 + REQUIRE(CHECK2CELLS(f2c, 0, 1)); + REQUIRE(CHECK2CELLS(f2c, 3, 6)); + // Line 1 + REQUIRE(CHECK3CELLS(f2c, 5, 4, 7)); + // Singular + REQUIRE_FALSE(f2c[2].has_value()); + } + } + WHEN("We agglomerate without forcing the direction") { + constexpr bool force_direction = false; + Agglomerator_Anisotropic aniso_agg( + fc_graph, + cc_graph, + seeds_pool, + Data.dim, + aniso_thresh, + agglomerationLines_Idx, + agglomerationLines, + Data.weights, + build_lines, + ODD_LINE_LENGTH, + MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, + force_direction + ); + aniso_agg.agglomerate_one_level(2, 2, 2, Data.weights, false); + aniso_agg.export_anisotropic_lines( + 1, agglomerationLines_Idx, agglomerationLines + ); + THEN("Lines grow vertically") { + const auto &f2c = cc_graph->_fc_2_cc; + const vector ref_lines_idx = {0, 2, 4}; + const vector ref_lines = {0, 1, 2, 3}; + REQUIRE_THAT(agglomerationLines_Idx, Equals(ref_lines_idx)); + REQUIRE_THAT(agglomerationLines, Equals(ref_lines)); + // Line 0 + REQUIRE(CHECK2CELLS(f2c, 0, 1)); + REQUIRE(CHECK2CELLS(f2c, 4, 5)); + // Line 1 + REQUIRE(CHECK2CELLS(f2c, 2, 3)); + REQUIRE(CHECK2CELLS(f2c, 6, 7)); + } + } + } } #undef CHECK2CELLS + +SCENARIO("Testing cell coupling", "[Anisotropic.CellCoupling]") { + GIVEN("Triangular mesh stretched in the x-direction representing a BL") { + const DualGEx_Tria_BL Data = DualGEx_Tria_BL(); + shared_ptr const seeds_pool = + make_shared(Data.n_bnd_faces, Data.weights, false); + shared_ptr fc_graph = make_shared( + Data.nb_fc, + Data.adjMatrix_row_ptr, + Data.adjMatrix_col_ind, + Data.adjMatrix_areaValues, + Data.volumes, + Data.centers, + Data.n_bnd_faces, + Data.dim, + Data.anisoCompliantCells + ); + // We recreate what's done inside build_anisotropic_lines + std::deque aniso_seeds_pool; + std::vector max_weights(fc_graph->_number_of_cells, 0.0); + std::vector to_treat(fc_graph->_number_of_cells, false); + std::vector alt_weights{}; + WHEN("Asking for maximum weight as coupling and all cells are anisotropic" + ) { + constexpr CoMMAWeightT aniso_thresh = -1.; + fc_graph->tag_anisotropic_cells( + alt_weights, + max_weights, + to_treat, + aniso_seeds_pool, + aniso_thresh, + Data.weights, + CoMMACellCouplingT::MAX_WEIGHT, + 0 + ); + THEN("Alternative weights are not computed") { + REQUIRE_THAT(alt_weights, EMPTY); + } + THEN("Max is computed correctly") { + // Largest side of the triangle + const std::vector ref_max( + fc_graph->_number_of_cells, sqrt101 + ); + REQUIRE_THAT(max_weights, Approx(ref_max).margin(eps)); + } + } + WHEN( + "Asking for maximum weight as coupling and checking anisotropy threshold" + ) { + // The result is same as in the previous one because the threshold is low + // enough to select all cells + constexpr CoMMAWeightT aniso_thresh = 2.; + fc_graph->tag_anisotropic_cells( + alt_weights, + max_weights, + to_treat, + aniso_seeds_pool, + aniso_thresh, + Data.weights, + CoMMACellCouplingT::MAX_WEIGHT, + 0 + ); + THEN("Alternative weights are not computed") { + REQUIRE_THAT(alt_weights, EMPTY); + } + THEN("Max is computed correctly") { + // Largest side of the triangle + const std::vector ref_max( + fc_graph->_number_of_cells, sqrt101 + ); + REQUIRE_THAT(max_weights, Approx(ref_max).margin(eps)); + } + } + WHEN("Asking for minimum distance and all cells are anisotropic") { + constexpr CoMMAWeightT aniso_thresh = -1.; + fc_graph->tag_anisotropic_cells( + alt_weights, + max_weights, + to_treat, + aniso_seeds_pool, + aniso_thresh, + Data.weights, + CoMMACellCouplingT::MIN_DISTANCE, + 0 + ); + const CoMMAWeightT x_dist = _3_1ov3; + const CoMMAWeightT y_dist = _1ov3; + const CoMMAWeightT same_block = 1. / (_sq(x_dist) + _sq(y_dist)); + const CoMMAWeightT side_block = 1. / (_sq(2. * x_dist) + _sq(y_dist)); + const CoMMAWeightT vert_block = 1. / (_sq(x_dist) + _sq(2. * y_dist)); + THEN("Alternative weights are computed") { + REQUIRE_THAT(alt_weights, SizeIs(fc_graph->_m_CRS_Values.size())); + // clang-format off + const std::vector ref_weights = { + /* 0 */ same_block, side_block, + /* 1 */ same_block, vert_block, + /* 2 */ same_block, + /* 3 */ side_block, same_block, vert_block, + /* 4 */ vert_block, same_block, side_block, + /* 5 */ same_block, + /* 6 */ vert_block, same_block, + /* 7 */ side_block, same_block + }; + // clang-format on + REQUIRE_THAT(alt_weights, Approx(ref_weights).margin(eps)); + } + THEN("Max is computed correctly") { + // Largest side of the triangle + const std::vector ref_max( + fc_graph->_number_of_cells, same_block + ); + REQUIRE_THAT(max_weights, Approx(ref_max).margin(eps)); + } + } + WHEN("Asking for minimum distance and checking anisotropy threshold") { + // The result is same as in the previous one because the threshold is low + // enough to select all cells + constexpr CoMMAWeightT aniso_thresh = 2.; + fc_graph->tag_anisotropic_cells( + alt_weights, + max_weights, + to_treat, + aniso_seeds_pool, + aniso_thresh, + Data.weights, + CoMMACellCouplingT::MIN_DISTANCE, + 0 + ); + const CoMMAWeightT x_dist = _3_1ov3; + const CoMMAWeightT y_dist = _1ov3; + const CoMMAWeightT same_block = 1. / (_sq(x_dist) + _sq(y_dist)); + const CoMMAWeightT side_block = 1. / (_sq(2. * x_dist) + _sq(y_dist)); + const CoMMAWeightT vert_block = 1. / (_sq(x_dist) + _sq(2. * y_dist)); + THEN("Alternative weights are computed") { + REQUIRE_THAT(alt_weights, SizeIs(fc_graph->_m_CRS_Values.size())); + // clang-format off + const std::vector ref_weights = { + /* 0 */ same_block, side_block, + /* 1 */ same_block, vert_block, + /* 2 */ same_block, + /* 3 */ side_block, same_block, vert_block, + /* 4 */ vert_block, same_block, side_block, + /* 5 */ same_block, + /* 6 */ vert_block, same_block, + /* 7 */ side_block, same_block + }; + // clang-format on + REQUIRE_THAT(alt_weights, Approx(ref_weights).margin(eps)); + } + THEN("Max is computed correctly") { + // Largest side of the triangle + const std::vector ref_max( + fc_graph->_number_of_cells, same_block + ); + REQUIRE_THAT(max_weights, Approx(ref_max).margin(eps)); + } + } + } +} + #undef CHECK3CELLS #undef CHECK4CELLS diff --git a/tests/test_defs.h b/tests/test_defs.h index 0716c6d..cb76b23 100644 --- a/tests/test_defs.h +++ b/tests/test_defs.h @@ -34,6 +34,13 @@ constexpr bool ODD_LINE_LENGTH = true; constexpr comma::CoMMAAspectRatioT DEFAULT_AR = comma::CoMMAAspectRatioT::DIAMETER_OVER_RADIUS; +constexpr std::optional MAX_CELLS_IN_LINE = std::nullopt; + +constexpr comma::CoMMACellCouplingT CELL_COUPLING_MAX = + comma::CoMMACellCouplingT::MAX_WEIGHT; + +constexpr bool FORCE_DIRECTION = true; + constexpr comma::CoMMAWeightT eps = 1e-10; using SeedsPoolT = comma::Seeds_Pool_Boundary_Priority< diff --git a/tests/test_structure.cpp b/tests/test_structure.cpp index eb8f1d3..7cf18a0 100644 --- a/tests/test_structure.cpp +++ b/tests/test_structure.cpp @@ -18,14 +18,13 @@ #include "CoMMA/Args.h" #include "CoMMA/CoMMA.h" +#include "CoMMA/CoMMADefs.h" #include "DualGraphExamples.h" #include "test_defs.h" using namespace comma; // NOLINT using namespace std; // NOLINT -#define MAX_CELLS_IN_LINE std::nullopt - SCENARIO("Test of a structure", "[structure]") { GIVEN("A simple graph, and we build the Dual Graph") { const DualGEx Data = DualGEx(); @@ -483,6 +482,7 @@ SCENARIO("Test of main function", "[structure]") { Data.centers, Data.weights, Data.anisoCompliantCells, Data.n_bnd_faces, build_lines, aniso, odd_length, aniso_thr, seed, fc2cc, alines_idx, alines, correction, Data.dim, goal_card, min_card, max_card, DEFAULT_AR, SING_CARD_THRESH, MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, FORCE_DIRECTION, 0) ); // Bad iteration number: greater than threshold @@ -492,6 +492,7 @@ SCENARIO("Test of main function", "[structure]") { Data.centers, Data.weights, Data.anisoCompliantCells, Data.n_bnd_faces, build_lines, aniso, odd_length, aniso_thr, seed, fc2cc, alines_idx, alines, correction, Data.dim, goal_card, min_card, max_card, DEFAULT_AR, SING_CARD_THRESH, MAX_CELLS_IN_LINE, + CELL_COUPLING_MAX, FORCE_DIRECTION, comma::iter_agglo_max_iter + 1) ); }