diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc index b87c17987..d1ea5addf 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc @@ -8,8 +8,10 @@ * nor does it submit to any jurisdiction. */ +#include #include #include +#include // for mkdir #include "ConservativeSphericalPolygonInterpolation.h" @@ -22,6 +24,7 @@ #include "atlas/mesh/actions/BuildNode2CellConnectivity.h" #include "atlas/meshgenerator.h" #include "atlas/parallel/mpi/mpi.h" +#include "atlas/parallel/omp/omp.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" #include "atlas/runtime/Trace.h" @@ -31,6 +34,8 @@ #include "eckit/log/Bytes.h" +#define PRINT_BAD_POLYGONS 0 + namespace atlas { namespace interpolation { namespace method { @@ -40,8 +45,99 @@ using util::ConvexSphericalPolygon; namespace { +template +std::string to_json(const It& begin, const It& end, int precision = 0) { + std::stringstream ss; + ss << "[\n"; + size_t size = std::distance(begin,end); + size_t c=0; + for( auto it = begin; it != end; ++it, ++c ) { + ss << " " << it->json(precision); + if( c < size-1 ) { + ss << ",\n"; + } + } + ss << "\n]"; + return ss.str(); +} + +template +std::string to_json(const ConvexSphericalPolygonContainer& polygons, int precision = 0) { + return to_json(polygons.begin(),polygons.end(),precision); +} + + MethodBuilder __builder("conservative-spherical-polygon"); +template +std::string polygons_to_json(const It& begin, const It& end, int precision = 0) { + std::stringstream ss; + ss << "[\n"; + size_t size = std::distance(begin,end); + size_t c=0; + for( auto it = begin; it != end; ++it, ++c ) { + ss << " " << it->json(precision); + if( c < size-1 ) { + ss << ",\n"; + } + } + ss << "\n]"; + return ss.str(); +} + +template +std::string polygons_to_json(const ConvexSphericalPolygonContainer& polygons, int precision = 0) { + return polygons_to_json(polygons.begin(),polygons.end(),precision); +} + +void dump_polygons_to_json( const ConvexSphericalPolygon& t_csp, + double pointsSameEPS, + const std::vector>& source_polygons, + const std::vector& source_polygons_considered_indices, + const std::string folder, + const std::string name) { + std::vector csp_arr{ t_csp }; + std::vector csp_arr_intersecting {t_csp}; + std::vector intersections; + int count = 1; + for( auto& s_idx : source_polygons_considered_indices ) { + auto s_csp = std::get<0>(source_polygons[s_idx]); + csp_arr.emplace_back( s_csp ); + std::fstream file_plg_debug(folder + name + "_" + std::to_string(count++) + ".debug", std::ios::out); + ConvexSphericalPolygon iplg = t_csp.intersect(s_csp, &file_plg_debug, pointsSameEPS); + file_plg_debug.close(); + if (iplg) { + if( iplg.area() > 0. ) { + csp_arr_intersecting.emplace_back( s_csp ); + intersections.emplace_back( iplg ); + } + } + } + double tot = 0.; + Log::info().indent(); + std::fstream file_info(folder + name + ".info", std::ios::out); + file_info << "List of intersection weights: " << std::endl; + count = 1; + for( auto& iplg : intersections ) { + csp_arr.emplace_back( iplg ); + csp_arr_intersecting.emplace_back( iplg ); + tot += iplg.area() / t_csp.area(); + file_info << "\t" << count++ << ": " << iplg.area() / t_csp.area() << std::endl; + } + file_info << std::endl << name + ": " << 100. * tot << " % covered."<< std::endl << std::endl; + file_info << "Target polygon + candidate source polygons + " << intersections.size() << " intersections in file:" << std::endl; + file_info << "\t" << folder + name + ".candidates\n" << std::endl; + std::fstream file_plg(folder + name + ".candidates", std::ios::out); + file_plg << polygons_to_json(csp_arr, 16); + file_plg.close(); + file_info << "Target polygon + intersecting source polygon + " << intersections.size() << " intersections in file:" << std::endl; + file_info << "\t" << folder + name + ".intersections\n" << std::endl; + file_plg.open(folder + name + ".intersections", std::ios::out); + file_plg << polygons_to_json(csp_arr_intersecting, 16); + file_plg.close(); + Log::info().unindent(); +} + constexpr double unit_sphere_area() { // 4*pi*r^2 with r=1 return 4. * M_PI; @@ -49,7 +145,7 @@ constexpr double unit_sphere_area() { template size_t memory_of(const std::vector& vector) { - return sizeof(T) * vector.size(); + return sizeof(T) * vector.capacity(); } template size_t memory_of(const std::vector>& vector_of_vector) { @@ -66,7 +162,7 @@ size_t memory_of( for (const auto& params : vector_of_params) { mem += memory_of(params.cell_idx); mem += memory_of(params.centroids); - mem += memory_of(params.src_weights); + mem += memory_of(params.weights); mem += memory_of(params.tgt_weights); } return mem; @@ -108,9 +204,6 @@ void sort_and_accumulate_triplets(std::vector& triplets) } } - -} // namespace - int inside_vertices(const ConvexSphericalPolygon& plg1, const ConvexSphericalPolygon& plg2, int& pout) { int points_in = 0; pout = 0; @@ -132,6 +225,8 @@ int inside_vertices(const ConvexSphericalPolygon& plg1, const ConvexSphericalPol return points_in; } +} // namespace + ConservativeSphericalPolygonInterpolation::ConservativeSphericalPolygonInterpolation(const Config& config): Method(config) { config.get("validate", validate_ = false); @@ -321,8 +416,9 @@ std::vector ConservativeSphericalPolygonInterpolation::get_node_neighbour // Create polygons for cell-centred data. Here, the polygons are mesh cells ConservativeSphericalPolygonInterpolation::CSPolygonArray -ConservativeSphericalPolygonInterpolation::get_polygons_celldata(Mesh& mesh) const { +ConservativeSphericalPolygonInterpolation::get_polygons_celldata(FunctionSpace fs) const { CSPolygonArray cspolygons; + auto mesh = extract_mesh(fs); const idx_t n_cells = mesh.cells().size(); cspolygons.resize(n_cells); const auto& cell2node = mesh.cells().node_connectivity(); @@ -331,7 +427,11 @@ ConservativeSphericalPolygonInterpolation::get_polygons_celldata(Mesh& mesh) con const auto& cell_flags = array::make_view(mesh.cells().flags()); const auto& cell_part = array::make_view(mesh.cells().partition()); std::vector pts_ll; + const int fs_halo = functionspace::CellColumns(fs).halo().size(); for (idx_t cell = 0; cell < n_cells; ++cell) { + if( cell_halo(cell) > fs_halo ) { + continue; + } int halo_type = cell_halo(cell); const idx_t n_nodes = cell2node.cols(cell); pts_ll.clear(); @@ -354,12 +454,13 @@ ConservativeSphericalPolygonInterpolation::get_polygons_celldata(Mesh& mesh) con // (cell_centre, edge_centre, cell_vertex, edge_centre) // additionally, subcell-to-node and node-to-subcells mapping are computed ConservativeSphericalPolygonInterpolation::CSPolygonArray -ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std::vector& csp2node, +ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(FunctionSpace fs, std::vector& csp2node, std::vector>& node2csp, std::array& errors) const { CSPolygonArray cspolygons; csp2node.clear(); node2csp.clear(); + auto mesh = extract_mesh(fs); node2csp.resize(mesh.nodes().size()); const auto nodes_ll = array::make_view(mesh.nodes().lonlat()); const auto& cell2node = mesh.cells().node_connectivity(); @@ -376,9 +477,19 @@ ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std eckit::geometry::Sphere::convertSphericalToCartesian(1., p_ll, p_xyz); return p_xyz; }; + const auto node_halo = array::make_view(mesh.nodes().halo()); + const auto node_ghost = array::make_view(mesh.nodes().ghost()); + const auto node_flags = array::make_view(mesh.nodes().flags()); + const auto node_part = array::make_view(mesh.nodes().partition()); + idx_t cspol_id = 0; // subpolygon enumeration errors = {0., 0.}; // over/undershoots in creation of subpolygons + const int fs_halo = functionspace::NodeColumns(fs).halo().size(); + for (idx_t cell = 0; cell < mesh.cells().size(); ++cell) { + if( cell_halo(cell) > fs_halo ) { + continue; + } ATLAS_ASSERT(cell < cell2node.rows()); const idx_t n_nodes = cell2node.cols(cell); ATLAS_ASSERT(n_nodes > 2); @@ -409,6 +520,7 @@ ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std PointLonLat cell_ll = xyz2ll(cell_mid); double loc_csp_area_shoot = ConvexSphericalPolygon(pts_ll).area(); // get ConvexSphericalPolygon for each valid edge + int halo_type; for (int inode = 0; inode < pts_idx.size(); inode++) { int inode_n = next_index(inode, pts_idx.size()); idx_t node = cell2node(cell, inode); @@ -429,18 +541,22 @@ ConservativeSphericalPolygonInterpolation::get_polygons_nodedata(Mesh& mesh, std subpol_pts_ll[1] = xyz2ll(iedge_mid); subpol_pts_ll[2] = pts_ll[inode_n]; subpol_pts_ll[3] = xyz2ll(jedge_mid); - int halo_type = cell_halo(cell); - if (util::Bitflags::view(cell_flags(cell)).check(util::Topology::PERIODIC) and - cell_part(cell) == mpi::rank()) { + halo_type = node_halo(node_n); + + if (util::Bitflags::view(node_flags(node_n)).check(util::Topology::PERIODIC) and + node_part(node_n) == mpi::rank()) { halo_type = -1; } + ConvexSphericalPolygon cspi(subpol_pts_ll); loc_csp_area_shoot -= cspi.area(); cspolygons.emplace_back(cspi, halo_type); cspol_id++; } - errors[0] += std::abs(loc_csp_area_shoot); - errors[1] = std::max(std::abs(loc_csp_area_shoot), errors[1]); + if (halo_type == 0) { + errors[0] += std::abs(loc_csp_area_shoot); + errors[1] = std::max(std::abs(loc_csp_area_shoot), errors[1]); + } } ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&errors[0], 1, eckit::mpi::sum()); @@ -465,7 +581,7 @@ void ConservativeSphericalPolygonInterpolation::do_setup_impl(const Grid& src_gr tgt_fs_ = functionspace::CellColumns(tgt_mesh_, option::halo(0)); } else { - tgt_fs_ = functionspace::NodeColumns(tgt_mesh_, option::halo(0)); + tgt_fs_ = functionspace::NodeColumns(tgt_mesh_, option::halo(1)); } } sharable_data_->tgt_fs_ = tgt_fs_; @@ -561,48 +677,68 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr tgt_mesh_ = extract_mesh(tgt_fs_); { - // we need src_halo_size >= 2, whereas tgt_halo_size >= 0 is enough - int src_halo_size = 0; - src_mesh_.metadata().get("halo", src_halo_size); - ATLAS_ASSERT(src_halo_size > 1); + // we need src_halo_size >= 2, tgt_halo_size >= 0 for CellColumns + // if target is NodeColumns, we need: + // tgt_halo_size >= 1 and + // src_halo_size large enough to cover the the target halo cells in the first row + int halo_size = 0; + src_mesh_.metadata().get("halo", halo_size); + if (halo_size < 2) { + Log::info() << "WARNING The halo size on source mesh should be at least 2.\n"; + } + if (not tgt_cell_data_) { + Log::info() << "WARNING The source cells should cover the first row of the target halos.\n"; + } + halo_size = 0; + tgt_mesh_.metadata().get("halo", halo_size); + if (not tgt_cell_data_ and halo_size == 0) { + Log::info() << "WARNING The halo size on target mesh should be at least 1 for the target NodeColumns.\n"; + } } CSPolygonArray src_csp; CSPolygonArray tgt_csp; - std::array errors = {0., 0.}; if (compute_cache) { + std::array errors = {0., 0.}; ATLAS_TRACE("Get source polygons"); StopWatch stopwatch; stopwatch.start(); if (src_cell_data_) { - src_csp = get_polygons_celldata(src_mesh_); + src_csp = get_polygons_celldata(src_fs_); } else { src_csp = - get_polygons_nodedata(src_mesh_, sharable_data_->src_csp2node_, sharable_data_->src_node2csp_, errors); + get_polygons_nodedata(src_fs_, sharable_data_->src_csp2node_, sharable_data_->src_node2csp_, errors); } stopwatch.stop(); sharable_data_->timings.source_polygons_assembly = stopwatch.elapsed(); - } - remap_stat_.errors[Statistics::Errors::SRC_PLG_L1] = errors[0]; - remap_stat_.errors[Statistics::Errors::SRC_PLG_LINF] = errors[1]; - if (compute_cache) { + remap_stat_.counts[Statistics::Counts::SRC_PLG] = src_csp.size(); + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_L1] = errors[0]; + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_LINF] = errors[1]; + + errors = {0., 0.}; ATLAS_TRACE("Get target polygons"); - StopWatch stopwatch; stopwatch.start(); if (tgt_cell_data_) { - tgt_csp = get_polygons_celldata(tgt_mesh_); + tgt_csp = get_polygons_celldata(tgt_fs_); } else { tgt_csp = - get_polygons_nodedata(tgt_mesh_, sharable_data_->tgt_csp2node_, sharable_data_->tgt_node2csp_, errors); + get_polygons_nodedata(tgt_fs_, sharable_data_->tgt_csp2node_, sharable_data_->tgt_node2csp_, errors); } stopwatch.stop(); sharable_data_->timings.target_polygons_assembly = stopwatch.elapsed(); + remap_stat_.counts[Statistics::Counts::TGT_PLG] = tgt_csp.size(); + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_L1] = errors[0]; + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_LINF] = errors[1]; + } + else { + remap_stat_.counts[Statistics::Counts::SRC_PLG] = -1111; + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::SRC_SUBPLG_LINF] = -1111.; + remap_stat_.counts[Statistics::Counts::TGT_PLG] = -1111; + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::TGT_SUBPLG_LINF] = -1111.; } - remap_stat_.counts[Statistics::Counts::SRC_PLG] = src_csp.size(); - remap_stat_.counts[Statistics::Counts::TGT_PLG] = tgt_csp.size(); - remap_stat_.errors[Statistics::Errors::TGT_PLG_L1] = errors[0]; - remap_stat_.errors[Statistics::Errors::TGT_PLG_LINF] = errors[1]; n_spoints_ = src_fs_.size(); n_tpoints_ = tgt_fs_.size(); @@ -611,82 +747,124 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr intersect_polygons(src_csp, tgt_csp); auto& src_points_ = sharable_data_->src_points_; - auto& tgt_points_ = sharable_data_->tgt_points_; src_points_.resize(n_spoints_); - tgt_points_.resize(n_tpoints_); sharable_data_->src_areas_.resize(n_spoints_); auto& src_areas_v = sharable_data_->src_areas_; - if (src_cell_data_) { - for (idx_t spt = 0; spt < n_spoints_; ++spt) { - const auto& s_csp = std::get<0>(src_csp[spt]); - src_points_[spt] = s_csp.centroid(); - src_areas_v[spt] = s_csp.area(); - } - } - else { - auto& src_node2csp_ = sharable_data_->src_node2csp_; - const auto lonlat = array::make_view(src_mesh_.nodes().lonlat()); - for (idx_t spt = 0; spt < n_spoints_; ++spt) { - if (src_node2csp_[spt].size() == 0) { - // this is a node to which no subpolygon is associated - // maximal twice per mesh we end here, and that is only when mesh has nodes on poles - auto p = PointLonLat{lonlat(spt, 0), lonlat(spt, 1)}; - eckit::geometry::Sphere::convertSphericalToCartesian(1., p, src_points_[spt]); - } - else { - // .. in the other case, start computing the barycentre - src_points_[spt] = PointXYZ{0., 0., 0.}; + { + ATLAS_TRACE("Store src_areas and src_point"); + if (src_cell_data_) { + for (idx_t spt = 0; spt < n_spoints_; ++spt) { + const auto& s_csp = std::get<0>(src_csp[spt]); + src_points_[spt] = s_csp.centroid(); + src_areas_v[spt] = s_csp.area(); } - src_areas_v[spt] = 0.; - for (idx_t isubcell = 0; isubcell < src_node2csp_[spt].size(); ++isubcell) { - idx_t subcell = src_node2csp_[spt][isubcell]; - const auto& s_csp = std::get<0>(src_csp[subcell]); - src_areas_v[spt] += s_csp.area(); - src_points_[spt] = src_points_[spt] + PointXYZ::mul(s_csp.centroid(), s_csp.area()); + } + else { + auto& src_node2csp_ = sharable_data_->src_node2csp_; + const auto lonlat = array::make_view(src_mesh_.nodes().lonlat()); + for (idx_t spt = 0; spt < n_spoints_; ++spt) { + if (src_node2csp_[spt].size() == 0) { + // this is a node to which no subpolygon is associated + // maximal twice per mesh we end here, and that is only when mesh has nodes on poles + auto p = PointLonLat{lonlat(spt, 0), lonlat(spt, 1)}; + eckit::geometry::Sphere::convertSphericalToCartesian(1., p, src_points_[spt]); + } + else { + // .. in the other case, start computing the barycentre + src_points_[spt] = PointXYZ{0., 0., 0.}; + } + src_areas_v[spt] = 0.; + for (idx_t isubcell = 0; isubcell < src_node2csp_[spt].size(); ++isubcell) { + idx_t subcell = src_node2csp_[spt][isubcell]; + const auto& s_csp = std::get<0>(src_csp[subcell]); + src_areas_v[spt] += s_csp.area(); + src_points_[spt] = src_points_[spt] + PointXYZ::mul(s_csp.centroid(), s_csp.area()); + } + double src_point_norm = PointXYZ::norm(src_points_[spt]); + if (src_point_norm == 0.) { + ATLAS_DEBUG_VAR(src_point_norm); + ATLAS_DEBUG_VAR(src_points_[spt]); + ATLAS_DEBUG_VAR(src_node2csp_[spt].size()); + for (idx_t isubcell = 0; isubcell < src_node2csp_[spt].size(); ++isubcell) { + idx_t subcell = src_node2csp_[spt][isubcell]; + ATLAS_DEBUG_VAR(subcell); + const auto& s_csp = std::get<0>(src_csp[subcell]); + s_csp.print(Log::info()); + Log::info() << std::endl; + src_areas_v[spt] += s_csp.area(); + ATLAS_DEBUG_VAR(s_csp.area()); + ATLAS_DEBUG_VAR(s_csp.centroid()); + src_points_[spt] = src_points_[spt] + PointXYZ::mul(s_csp.centroid(), s_csp.area()); + } + Log::info().flush(); + // something went wrong, improvise + src_point_norm = 1.; + } + src_points_[spt] = PointXYZ::div(src_points_[spt], src_point_norm); } - double src_point_norm = PointXYZ::norm(src_points_[spt]); - ATLAS_ASSERT(src_point_norm > 0.); - src_points_[spt] = PointXYZ::div(src_points_[spt], src_point_norm); } } + auto& tgt_points_ = sharable_data_->tgt_points_; + tgt_points_.resize(n_tpoints_); sharable_data_->tgt_areas_.resize(n_tpoints_); auto& tgt_areas_v = sharable_data_->tgt_areas_; - if (tgt_cell_data_) { - for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { - const auto& t_csp = std::get<0>(tgt_csp[tpt]); - tgt_points_[tpt] = t_csp.centroid(); - tgt_areas_v[tpt] = t_csp.area(); - } - } - else { - auto& tgt_node2csp_ = sharable_data_->tgt_node2csp_; - const auto lonlat = array::make_view(tgt_mesh_.nodes().lonlat()); - for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { - if (tgt_node2csp_[tpt].size() == 0) { - // this is a node to which no subpolygon is associated - // maximal twice per mesh we end here, and that is only when mesh has nodes on poles - auto p = PointLonLat{lonlat(tpt, 0), lonlat(tpt, 1)}; - eckit::geometry::Sphere::convertSphericalToCartesian(1., p, tgt_points_[tpt]); - } - else { - // .. in the other case, start computing the barycentre - tgt_points_[tpt] = PointXYZ{0., 0., 0.}; + { + ATLAS_TRACE("Store src_areas and src_point"); + if (tgt_cell_data_) { + for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { + const auto& t_csp = std::get<0>(tgt_csp[tpt]); + tgt_points_[tpt] = t_csp.centroid(); + tgt_areas_v[tpt] = t_csp.area(); } - tgt_areas_v[tpt] = 0.; - for (idx_t isubcell = 0; isubcell < tgt_node2csp_[tpt].size(); ++isubcell) { - idx_t subcell = tgt_node2csp_[tpt][isubcell]; - const auto& t_csp = std::get<0>(tgt_csp[subcell]); - tgt_areas_v[tpt] += t_csp.area(); - tgt_points_[tpt] = tgt_points_[tpt] + PointXYZ::mul(t_csp.centroid(), t_csp.area()); + } + else { + auto& tgt_node2csp_ = sharable_data_->tgt_node2csp_; + const auto lonlat = array::make_view(tgt_mesh_.nodes().lonlat()); + for (idx_t tpt = 0; tpt < n_tpoints_; ++tpt) { + if (tgt_node2csp_[tpt].size() == 0) { + // this is a node to which no subpolygon is associated + // maximal twice per mesh we end here, and that is only when mesh has nodes on poles + auto p = PointLonLat{lonlat(tpt, 0), lonlat(tpt, 1)}; + eckit::geometry::Sphere::convertSphericalToCartesian(1., p, tgt_points_[tpt]); + } + else { + // .. in the other case, start computing the barycentre + tgt_points_[tpt] = PointXYZ{0., 0., 0.}; + } + tgt_areas_v[tpt] = 0.; + for (idx_t isubcell = 0; isubcell < tgt_node2csp_[tpt].size(); ++isubcell) { + idx_t subcell = tgt_node2csp_[tpt][isubcell]; + const auto& t_csp = std::get<0>(tgt_csp[subcell]); + tgt_areas_v[tpt] += t_csp.area(); + tgt_points_[tpt] = tgt_points_[tpt] + PointXYZ::mul(t_csp.centroid(), t_csp.area()); + } + double tgt_point_norm = PointXYZ::norm(tgt_points_[tpt]); + ATLAS_ASSERT(tgt_point_norm > 0.); + if (tgt_point_norm == 0.) { + ATLAS_DEBUG_VAR(tgt_point_norm); + ATLAS_DEBUG_VAR(tgt_points_[tpt]); + ATLAS_DEBUG_VAR(tgt_node2csp_[tpt].size()); + for (idx_t isubcell = 0; isubcell < tgt_node2csp_[tpt].size(); ++isubcell) { + idx_t subcell = tgt_node2csp_[tpt][isubcell]; + ATLAS_DEBUG_VAR(subcell); + const auto& t_csp = std::get<0>(tgt_csp[subcell]); + t_csp.print(Log::info()); + Log::info() << std::endl; + tgt_areas_v[tpt] += t_csp.area(); + ATLAS_DEBUG_VAR(t_csp.area()); + ATLAS_DEBUG_VAR(t_csp.centroid()); + tgt_points_[tpt] = tgt_points_[tpt] + PointXYZ::mul(t_csp.centroid(), t_csp.area()); + } + Log::info().flush(); + // something went wrong, improvise + tgt_point_norm = 1.; + } + tgt_points_[tpt] = PointXYZ::div(tgt_points_[tpt], tgt_point_norm); } - double tgt_point_norm = PointXYZ::norm(tgt_points_[tpt]); - ATLAS_ASSERT(tgt_point_norm > 0.); - tgt_points_[tpt] = PointXYZ::div(tgt_points_[tpt], tgt_point_norm); } } } - if (not matrix_free_) { StopWatch stopwatch; stopwatch.start(); @@ -718,30 +896,6 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr } } -namespace { -// needed for intersect_polygons only, merely for detecting duplicate points -struct ComparePointXYZ { - bool operator()(const PointXYZ& f, const PointXYZ& s) const { - // eps = ConvexSphericalPolygon::EPS which is the threshold when two points are "same" - double eps = 1e4 * std::numeric_limits::epsilon(); - if (f[0] < s[0] - eps) { - return true; - } - else if (std::abs(f[0] - s[0]) < eps) { - if (f[1] < s[1] - eps) { - return true; - } - else if (std::abs(f[1] - s[1]) < eps) { - if (f[2] < s[2] - eps) { - return true; - } - } - } - return false; - } -}; -} // namespace - void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolygonArray& src_csp, const CSPolygonArray& tgt_csp) { ATLAS_TRACE(); @@ -751,8 +905,10 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg util::KDTree kdt_search; kdt_search.reserve(tgt_csp.size()); double max_tgtcell_rad = 0.; + const int tgt_halo_intersection_depth = (tgt_cell_data_ ? 0 : 1); // if target NodeColumns, one target halo required for subcells around target nodes for (idx_t jcell = 0; jcell < tgt_csp.size(); ++jcell) { - if (std::get<1>(tgt_csp[jcell]) == 0) { + auto tgt_halo_type = std::get<1>(tgt_csp[jcell]); + if (tgt_halo_type <= tgt_halo_intersection_depth) { // and tgt_halo_type != -1) { const auto& t_csp = std::get<0>(tgt_csp[jcell]); kdt_search.insert(t_csp.centroid(), jcell); max_tgtcell_rad = std::max(max_tgtcell_rad, t_csp.radius()); @@ -767,18 +923,38 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg StopWatch stopwatch_polygon_intersections; stopwatch_src_already_in.start(); - std::set src_cent; - auto polygon_point = [](const ConvexSphericalPolygon& pol) { - PointXYZ p{0., 0., 0.}; - for (int i = 0; i < pol.size(); i++) { - p = p + pol[i]; - } - p /= pol.size(); - return p; + + + // needed for intersect_polygons only, merely for detecting duplicate points + // Treshold at which points are considered same + double compare_pointxyz_eps = 1.e8 * std::numeric_limits::epsilon(); + if (::getenv("ATLAS_COMPAREPOINTXYZ_EPS_FACTOR")) { + compare_pointxyz_eps = std::atof(::getenv("ATLAS_COMPAREPOINTXYZ_EPS_FACTOR")) * std::numeric_limits::epsilon(); + } + + auto compare_pointxyz = [eps=compare_pointxyz_eps] (const PointXYZ& f, const PointXYZ& s) -> bool { + if (f[0] < s[0] - eps) { + return true; + } + else if (std::abs(f[0] - s[0]) < eps) { + if (f[1] < s[1] - eps) { + return true; + } + else if (std::abs(f[1] - s[1]) < eps) { + if (f[2] < s[2] - eps) { + return true; + } + } + } + return false; }; + + std::set src_cent(compare_pointxyz); auto src_already_in = [&](const PointXYZ& halo_cent) { if (src_cent.find(halo_cent) == src_cent.end()) { - src_cent.insert(halo_cent); + atlas_omp_critical{ + src_cent.insert(halo_cent); + } return false; } return true; @@ -807,79 +983,114 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg tgt_iparam.resize(tgt_csp.size()); } + // the worst target polygon coverage for analysis of intersection + std::pair worst_tgt_overcover; + std::pair worst_tgt_undercover; + worst_tgt_overcover.first = -1; + worst_tgt_overcover.second = -1.; + worst_tgt_undercover.first = -1; + worst_tgt_undercover.second = M_PI; + + // NOTE: polygon vertex points at distance < pointsSameEPS will be replaced with one point + constexpr double pointsSameEPS = 5.e6 * std::numeric_limits::epsilon(); + eckit::Channel blackhole; - eckit::ProgressTimer progress("Intersecting polygons ", src_csp.size(), " cell", double(10), - src_csp.size() > 50 ? Log::info() : blackhole); - for (idx_t scell = 0; scell < src_csp.size(); ++scell, ++progress) { - stopwatch_src_already_in.start(); - if (src_already_in(polygon_point(std::get<0>(src_csp[scell])))) { + eckit::ProgressTimer progress("Intersecting polygons ", src_csp.size() / atlas_omp_get_max_threads(), " (cell/thread)", double(10), + src_csp.size() / atlas_omp_get_max_threads() > 50 ? Log::info() : blackhole); + atlas_omp_parallel_for (idx_t scell = 0; scell < src_csp.size(); ++scell) { + if ( atlas_omp_get_thread_num() == 0 ) { + ++progress; + stopwatch_src_already_in.start(); + } + bool already_in = src_already_in((std::get<0>(src_csp[scell]).centroid())); + if ( atlas_omp_get_thread_num() == 0 ) { stopwatch_src_already_in.stop(); - continue; } - stopwatch_src_already_in.stop(); - - const auto& s_csp = std::get<0>(src_csp[scell]); - const double s_csp_area = s_csp.area(); - double src_cover_area = 0.; - - stopwatch_kdtree_search.start(); - auto tgt_cells = kdt_search.closestPointsWithinRadius(s_csp.centroid(), s_csp.radius() + max_tgtcell_rad); - stopwatch_kdtree_search.stop(); - for (idx_t ttcell = 0; ttcell < tgt_cells.size(); ++ttcell) { - auto tcell = tgt_cells[ttcell].payload(); - const auto& t_csp = std::get<0>(tgt_csp[tcell]); - stopwatch_polygon_intersections.start(); - ConvexSphericalPolygon csp_i = s_csp.intersect(t_csp); - double csp_i_area = csp_i.area(); - stopwatch_polygon_intersections.stop(); - if (validate_) { - // check zero area intersections with inside_vertices - int pout; - if (inside_vertices(s_csp, t_csp, pout) > 2 && csp_i.area() < 3e-16) { - dump_intersection(s_csp, tgt_csp, tgt_cells); + if (not already_in) { + const auto& s_csp = std::get<0>(src_csp[scell]); + const double s_csp_area = s_csp.area(); + double src_cover_area = 0.; + + stopwatch_kdtree_search.start(); + auto tgt_cells = kdt_search.closestPointsWithinRadius(s_csp.centroid(), s_csp.radius() + max_tgtcell_rad); + stopwatch_kdtree_search.stop(); + for (idx_t ttcell = 0; ttcell < tgt_cells.size(); ++ttcell) { + auto tcell = tgt_cells[ttcell].payload(); + const auto& t_csp = std::get<0>(tgt_csp[tcell]); + if( atlas_omp_get_thread_num() == 0 ) { + stopwatch_polygon_intersections.start(); + } + ConvexSphericalPolygon csp_i = s_csp.intersect(t_csp, nullptr, pointsSameEPS); + double csp_i_area = csp_i.area(); + if( atlas_omp_get_thread_num() == 0 ) { + stopwatch_polygon_intersections.stop(); } - } - if (csp_i_area > 0.) { if (validate_) { - tgt_iparam[tcell].cell_idx.emplace_back(scell); - tgt_iparam[tcell].tgt_weights.emplace_back(csp_i_area); + int pout; + // TODO: this can be removed soon + if (inside_vertices(s_csp, t_csp, pout) > 2 && csp_i.area() < 3e-16) { + dump_intersection("Zero area intersections with inside_vertices", s_csp, tgt_csp, tgt_cells); + } + // TODO: assuming intersector search works fine, this should be move under "if (csp_i_area > 0)" + atlas_omp_critical { + tgt_iparam[tcell].cell_idx.emplace_back(scell); + tgt_iparam[tcell].tgt_weights.emplace_back(csp_i_area); + } } - src_iparam_[scell].cell_idx.emplace_back(tcell); - src_iparam_[scell].src_weights.emplace_back(csp_i_area); - double target_weight = csp_i_area / t_csp.area(); - src_iparam_[scell].tgt_weights.emplace_back(target_weight); - src_iparam_[scell].centroids.emplace_back(csp_i.centroid()); - src_cover_area += csp_i_area; - ATLAS_ASSERT(target_weight < 1.1); - ATLAS_ASSERT(csp_i_area / s_csp_area < 1.1); - } - } - const double src_cover_err = std::abs(s_csp_area - src_cover_area); - const double src_cover_err_percent = 100. * src_cover_err / s_csp_area; - if (src_cover_err_percent > 0.1 and std::get<1>(src_csp[scell]) == 0) { - // HACK: source cell at process boundary will not be covered by target cells, skip them - // TODO: mark these source cells beforehand and compute error in them among the processes - - if (validate_) { - if (mpi::size() == 1) { - Log::info() << "WARNING src cell covering error : " << src_cover_err_percent << "%\n"; - dump_intersection(s_csp, tgt_csp, tgt_cells); + if (csp_i_area > 0) { + src_iparam_[scell].cell_idx.emplace_back(tcell); + src_iparam_[scell].weights.emplace_back(csp_i_area); + double target_weight = csp_i_area / t_csp.area(); + src_iparam_[scell].tgt_weights.emplace_back(target_weight); // TODO: tgt_weights vector should be removed for the sake of highres + if (order_ == 2 or not matrix_free_ or not matrixAllocated()) { + src_iparam_[scell].centroids.emplace_back(csp_i.centroid()); + } + src_cover_area += csp_i_area; + + if (validate_) { + // this check is concerned with accuracy of polygon intersections + if (target_weight > 1.1) { + dump_intersection("Intersection larger than target", s_csp, tgt_csp, tgt_cells); + } + if (csp_i_area / s_csp_area > 1.1) { + dump_intersection("Intersection larger than source", s_csp, tgt_csp, tgt_cells); + } + } } } - area_coverage[TOTAL_SRC] += src_cover_err; - area_coverage[MAX_SRC] = std::max(area_coverage[MAX_SRC], src_cover_err); - } - if (src_iparam_[scell].cell_idx.size() == 0) { - num_pol[SRC_NONINTERSECT]++; - } - if (normalise_intersections_ && src_cover_err_percent < 1.) { - double wfactor = s_csp.area() / (src_cover_area > 0. ? src_cover_area : 1.); - for (idx_t i = 0; i < src_iparam_[scell].src_weights.size(); i++) { - src_iparam_[scell].src_weights[i] *= wfactor; - src_iparam_[scell].tgt_weights[i] *= wfactor; + const double src_cover_err = std::abs(s_csp_area - src_cover_area); + const double src_cover_err_percent = 100. * src_cover_err / s_csp_area; + if (src_cover_err_percent > 0.1 and std::get<1>(src_csp[scell]) == 0) { + // NOTE: source cell at process boundary will not be covered by target cells, skip them + // TODO: mark these source cells beforehand and compute error in them among the processes + if (validate_ and mpi::size() == 1) { + dump_intersection("Source cell not exactly covered", s_csp, tgt_csp, tgt_cells); + if (statistics_intersection_) { + atlas_omp_critical{ + area_coverage[TOTAL_SRC] += src_cover_err; + area_coverage[MAX_SRC] = std::max(area_coverage[MAX_SRC], src_cover_err); + } + } + } } - } - num_pol[SRC_TGT_INTERSECT] += src_iparam_[scell].src_weights.size(); + if (src_iparam_[scell].cell_idx.size() == 0 and statistics_intersection_) { + atlas_omp_critical{ + num_pol[SRC_NONINTERSECT]++; + } + } + if (normalise_intersections_ && src_cover_err_percent < 1.) { + double wfactor = s_csp.area() / (src_cover_area > 0. ? src_cover_area : 1.); + for (idx_t i = 0; i < src_iparam_[scell].weights.size(); i++) { + src_iparam_[scell].weights[i] *= wfactor; + src_iparam_[scell].tgt_weights[i] *= wfactor; + } + } + if (statistics_intersection_) { + atlas_omp_critical{ + num_pol[SRC_TGT_INTERSECT] += src_iparam_[scell].weights.size(); + } + } + } // already in } timings.polygon_intersections = stopwatch_polygon_intersections.elapsed(); timings.target_kdtree_search = stopwatch_kdtree_search.elapsed(); @@ -888,48 +1099,151 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg num_pol[TGT] = tgt_csp.size(); ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(num_pol.data(), num_pol.size(), eckit::mpi::sum()); - mpi::comm().allReduceInPlace(area_coverage.data(), area_coverage.size(), eckit::mpi::max()); } remap_stat_.counts[Statistics::Counts::INT_PLG] = num_pol[SRC_TGT_INTERSECT]; remap_stat_.counts[Statistics::Counts::UNCVR_SRC] = num_pol[SRC_NONINTERSECT]; - remap_stat_.errors[Statistics::Errors::GEO_L1] = area_coverage[TOTAL_SRC]; - remap_stat_.errors[Statistics::Errors::GEO_LINF] = area_coverage[MAX_SRC]; - - double geo_err_l1 = 0.; - double geo_err_linf = 0.; - for (idx_t scell = 0; scell < src_csp.size(); ++scell) { - const int cell_flag = std::get<1>(src_csp[scell]); - if (cell_flag == -1 or cell_flag > 0) { - // skip periodic & halo cells - continue; + + const std::string polygon_intersection_folder = "polygon_intersection/"; + if (validate_ && mpi::rank() == 0) { + if (mkdir(polygon_intersection_folder.c_str(), 0777) != 0) { + Log::info() << "WARNING Polygon intersection relevant information in is the folder \e[1mpolygon_intersection\e[0m." << std::endl; } - double diff_cell = std::get<0>(src_csp[scell]).area(); - for (idx_t icell = 0; icell < src_iparam_[scell].src_weights.size(); ++icell) { - diff_cell -= src_iparam_[scell].src_weights[icell]; + else { + Log::info() << "WARNING Could not create the folder \e[1mpolygon_intersection\e[0m." << std::endl; } - geo_err_l1 += std::abs(diff_cell); - geo_err_linf = std::max(geo_err_linf, std::abs(diff_cell)); } - ATLAS_TRACE_MPI(ALLREDUCE) { - mpi::comm().allReduceInPlace(geo_err_l1, eckit::mpi::sum()); - mpi::comm().allReduceInPlace(geo_err_linf, eckit::mpi::max()); - } - remap_stat_.errors[Statistics::Errors::GEO_L1] = geo_err_l1 / unit_sphere_area(); - remap_stat_.errors[Statistics::Errors::GEO_LINF] = geo_err_linf; if (validate_) { + double geo_err_l1 = 0.; + double geo_err_linf = -1.; + for (idx_t scell = 0; scell < src_csp.size(); ++scell) { + if (std::get<1>(src_csp[scell]) != 0 ) { + // skip periodic & halo cells + continue; + } + double diff_cell = std::get<0>(src_csp[scell]).area(); + for (idx_t icell = 0; icell < src_iparam_[scell].weights.size(); ++icell) { + diff_cell -= src_iparam_[scell].weights[icell]; + } + geo_err_l1 += std::abs(diff_cell); + geo_err_linf = std::max(geo_err_linf, std::abs(diff_cell)); + } + ATLAS_TRACE_MPI(ALLREDUCE) { + mpi::comm().allReduceInPlace(geo_err_l1, eckit::mpi::sum()); + mpi::comm().allReduceInPlace(geo_err_linf, eckit::mpi::max()); + } + if (mpi::size() == 1) { + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_L1] = geo_err_l1 / unit_sphere_area(); + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF] = geo_err_linf; + } + else { + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF] = -1111.; + } + + std::fstream polygon_intersection_info("polygon_intersection", std::ios::out); + + geo_err_l1 = 0.; + geo_err_linf = -1.; for (idx_t tcell = 0; tcell < tgt_csp.size(); ++tcell) { + if (std::get<1>(tgt_csp[tcell]) != 0) { + // skip periodic & halo cells + continue; + } const auto& t_csp = std::get<0>(tgt_csp[tcell]); + + double tgt_cover_area = 0.; const auto& tiparam = tgt_iparam[tcell]; +#if 0 + // dump polygons in json format + if( tcell == 27609 ) { + dump_polygons_to_json(t_csp, src_csp, tiparam.cell_idx, 1.e-16); + } +#endif for (idx_t icell = 0; icell < tiparam.cell_idx.size(); ++icell) { tgt_cover_area += tiparam.tgt_weights[icell]; } - const double tgt_cover_err_percent = 100. * std::abs(t_csp.area() - tgt_cover_area) / t_csp.area(); - if (tgt_cover_err_percent > 0.1 and std::get<1>(tgt_csp[tcell]) == 0) { - Log::info() << "WARNING tgt cell covering error : " << tgt_cover_err_percent << " %\n"; - dump_intersection(t_csp, src_csp, tiparam.cell_idx); + /* + // TODO: normalise to target cell + double normm = t_csp.area() / (tgt_cover_area > 0. ? tgt_cover_area : t_csp.area()); + for (idx_t icell = 0; icell < tiparam.cell_idx.size(); ++icell) { + idx_t scell = tiparam.cell_idx[icell]; + auto siparam = src_iparam_[scell]; + size_t tgt_intersectors = siparam.cell_idx.size(); + for (idx_t sicell = 0; sicell < tgt_intersectors; sicell++ ) { + if (siparam.cell_idx[icell] == tcell) {; + siparam.weights[icell] *= normm; + siparam.tgt_weights[icell] *= normm; + } + } } + */ + double diff_cell = tgt_cover_area - t_csp.area(); + geo_err_l1 += std::abs(diff_cell); + geo_err_linf = std::max(geo_err_linf, std::abs(diff_cell)); + const double tgt_cover_err = 100. * diff_cell / t_csp.area(); + if (worst_tgt_overcover.second < tgt_cover_err) { + worst_tgt_overcover.second = tgt_cover_err;; + worst_tgt_overcover.first = tcell; + } + if (worst_tgt_undercover.second > tgt_cover_err) { + worst_tgt_undercover.second = tgt_cover_err;; + worst_tgt_undercover.first = tcell; + } + if (std::abs(tgt_cover_err) > 0.5) { + PointLonLat centre_ll; + eckit::geometry::Sphere::convertCartesianToSpherical(1., t_csp.centroid(), centre_ll); + polygon_intersection_info << "WARNING tgt cell " << tcell << " over-covering: \e[1m" << tgt_cover_err << "\e[0m %, cell-centre: " + << centre_ll <<"\n"; + polygon_intersection_info << "source indices: " << tiparam.cell_idx << std::endl; + dump_intersection("Target cell not exaclty covered", t_csp, src_csp, tiparam.cell_idx); + //dump_polygons_to_json(t_csp, src_csp, tiparam.cell_idx, "bad_polygon", 1.e-16); + } + } + ATLAS_TRACE_MPI(ALLREDUCE) { + mpi::comm().allReduceInPlace(geo_err_l1, eckit::mpi::sum()); + mpi::comm().allReduceInPlace(geo_err_linf, eckit::mpi::max()); + } + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_L1] = geo_err_l1 / unit_sphere_area(); + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_LINF] = geo_err_linf; + } + else { + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF] = -1111.; + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_L1] = -1111.; + remap_stat_.errors[Statistics::Errors::TGT_INTERSECTPLG_LINF] = -1111.; + } + + if (validate_) { + std::vector first(mpi::comm().size()); + std::vector second(mpi::comm().size()); + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGather(worst_tgt_overcover.first, first.begin(), first.end()); + mpi::comm().allGather(worst_tgt_overcover.second, second.begin(), second.end()); + } + auto max_over = std::max_element(second.begin(), second.end()); + auto rank_over = std::distance(second.begin(), max_over); + Log::info() << "WARNING The worst target polygon over-coveraging: \e[1m" + << *max_over + << "\e[0m %. For details, check the file: worst_target_cell_overcover.info " << std::endl; + if (rank_over == mpi::rank()) { + auto tcell = worst_tgt_overcover.first; + dump_polygons_to_json(std::get<0>(tgt_csp[tcell]), pointsSameEPS, src_csp, tgt_iparam[tcell].cell_idx, polygon_intersection_folder, "worst_target_cell_overcover"); + } + ATLAS_TRACE_MPI(ALLGATHER) { + mpi::comm().allGather(worst_tgt_undercover.first, first.begin(), first.end()); + mpi::comm().allGather(worst_tgt_undercover.second, second.begin(), second.end()); + } + auto min_under = std::min_element(second.begin(), second.end()); + auto rank_under = std::distance(second.begin(), min_under); + Log::info() << "WARNING The worst target polygon under-coveraging: \e[1m" + << *min_under + << "\e[0m %. For details, check the file: worst_target_cell_undercover.info " << std::endl; + + if (rank_under == mpi::rank()) { + auto tcell = worst_tgt_undercover.first; + dump_polygons_to_json(std::get<0>(tgt_csp[tcell]), pointsSameEPS, src_csp, tgt_iparam[tcell].cell_idx, polygon_intersection_folder, "worst_target_cell_undercover"); } } } @@ -943,7 +1257,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 // determine the size of array of triplets used to define the sparse matrix if (src_cell_data_) { for (idx_t scell = 0; scell < n_spoints_; ++scell) { - triplets_size += src_iparam_[scell].centroids.size(); + triplets_size += src_iparam_[scell].cell_idx.size(); } } else { @@ -962,7 +1276,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 if (src_cell_data_ && tgt_cell_data_) { for (idx_t scell = 0; scell < n_spoints_; ++scell) { const auto& iparam = src_iparam_[scell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; triplets.emplace_back(tcell, scell, iparam.tgt_weights[icell]); } @@ -974,8 +1288,9 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { const idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; + ATLAS_ASSERT(tcell < n_tpoints_); triplets.emplace_back(tcell, snode, iparam.tgt_weights[icell]); } } @@ -985,11 +1300,15 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 auto& tgt_csp2node_ = data_->tgt_csp2node_; for (idx_t scell = 0; scell < n_spoints_; ++scell) { const auto& iparam = src_iparam_[scell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; + if (tnode >= n_tpoints_) { + Log::info() << "tnode, n_tpoints = " << tnode << ", " << n_tpoints_ << std::endl; + ATLAS_ASSERT(false); + } double inv_node_weight = (tgt_areas_v[tnode] > 0. ? 1. / tgt_areas_v[tnode] : 0.); - triplets.emplace_back(tnode, scell, iparam.src_weights[icell] * inv_node_weight); + triplets.emplace_back(tnode, scell, iparam.weights[icell] * inv_node_weight); } } } @@ -1000,16 +1319,42 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { const idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; + ATLAS_ASSERT(tnode < n_tpoints_); double inv_node_weight = (tgt_areas_v[tnode] > 0. ? 1. / tgt_areas_v[tnode] : 0.); - triplets.emplace_back(tnode, snode, iparam.src_weights[icell] * inv_node_weight); + triplets.emplace_back(tnode, snode, iparam.weights[icell] * inv_node_weight); + if( tnode >= n_tpoints_) { + Log::info() << tnode << " = tnode, " << n_tpoints_ << " = n_tpoints\n";; + ATLAS_ASSERT(false); + } } } } } sort_and_accumulate_triplets(triplets); // Very expensive!!! (90% of this routine). We need to avoid it + +// if (validate_) { +// std::vector weight_sum(n_tpoints_); +// for( auto& triplet : triplets ) { +// weight_sum[triplet.row()] += triplet.value(); +// } +// if (order_ == 1) { +// // first order should not give overshoots +// double eps = 1e4 * std::numeric_limits::epsilon(); +// for( auto& triplet : triplets ) { +// if (triplet.value() > 1. + eps or triplet.value() < -eps) { +// Log::info() << "target point " << triplet.row() << " weight: " << triplet.value() << std::endl; +// } +// } +// } +// for( size_t row=0; row < n_tpoints_; ++row ) { +// if (std::abs(weight_sum[row] - 1.) > 1e-11) { +// Log::info() << "target weight in row " << row << " differs from 1 by " << std::abs(weight_sum[row] - 1.) << std::endl; +// } +// } +// } return Matrix(n_tpoints_, n_spoints_, triplets); } @@ -1026,28 +1371,30 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 const auto src_halo = array::make_view(src_mesh_.cells().halo()); for (idx_t scell = 0; scell < n_spoints_; ++scell) { const auto nb_cells = get_cell_neighbours(src_mesh_, scell); - triplets_size += (2 * nb_cells.size() + 1) * src_iparam_[scell].centroids.size(); + triplets_size += (2 * nb_cells.size() + 1) * src_iparam_[scell].cell_idx.size(); } triplets.reserve(triplets_size); for (idx_t scell = 0; scell < n_spoints_; ++scell) { - const auto nb_cells = get_cell_neighbours(src_mesh_, scell); const auto& iparam = src_iparam_[scell]; - if (iparam.centroids.size() == 0 && not src_halo(scell)) { + if (iparam.cell_idx.size() == 0 && not src_halo(scell)) { continue; } /* // better conservation after Kritsikis et al. (2017) + // NOTE: ommited here at cost of conservation due to more involved implementation in parallel + // TEMP: use barycentre computed based on the source polygon's vertices, rather then + // the numerical barycentre based on barycentres of the intersections with target cells. + // TODO: for a given source cell, collect the centroids of all its intersections with target cells + // to compute the numerical barycentre of the cell bases on intersection. PointXYZ Cs = {0., 0., 0.}; - for ( idx_t icell = 0; icell < iparam.centroids.size(); ++icell ) { - Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.src_weights[icell] ); + for ( idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell ) { + Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.weights[icell] ); } - const double Cs_norm = PointXYZ::norm( Cs ); - ATLAS_ASSERT( Cs_norm > 0. ); - Cs = PointXYZ::div( Cs, Cs_norm ); */ const PointXYZ& Cs = src_points_[scell]; // compute gradient from cells double dual_area_inv = 0.; std::vector Rsj; + const auto nb_cells = get_cell_neighbours(src_mesh_, scell); Rsj.resize(nb_cells.size()); for (idx_t j = 0; j < nb_cells.size(); ++j) { idx_t nj = next_index(j, nb_cells.size()); @@ -1071,15 +1418,15 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } // assemble the matrix std::vector Aik; - Aik.resize(iparam.centroids.size()); - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + Aik.resize(iparam.cell_idx.size()); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const PointXYZ& Csk = iparam.centroids[icell]; const PointXYZ Csk_Cs = Csk - Cs; Aik[icell] = Csk_Cs - PointXYZ::mul(Cs, PointXYZ::dot(Cs, Csk_Cs)); Aik[icell] = PointXYZ::mul(Aik[icell], iparam.tgt_weights[icell] * dual_area_inv); } if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const idx_t tcell = iparam.cell_idx[icell]; for (idx_t j = 0; j < nb_cells.size(); ++j) { idx_t nj = next_index(j, nb_cells.size()); @@ -1093,11 +1440,11 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } else { auto& tgt_csp2node_ = data_->tgt_csp2node_; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; double inv_node_weight = (tgt_areas_v[tnode] > 0.) ? 1. / tgt_areas_v[tnode] : 0.; - double csp2node_coef = iparam.src_weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; + double csp2node_coef = iparam.weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; for (idx_t j = 0; j < nb_cells.size(); ++j) { idx_t nj = next_index(j, nb_cells.size()); idx_t sj = nb_cells[j]; @@ -1118,7 +1465,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 const auto nb_nodes = get_node_neighbours(src_mesh_, snode); for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { idx_t subcell = src_node2csp_[snode][isubcell]; - triplets_size += (2 * nb_nodes.size() + 1) * src_iparam_[subcell].centroids.size(); + triplets_size += (2 * nb_nodes.size() + 1) * src_iparam_[subcell].cell_idx.size(); } } triplets.reserve(triplets_size); @@ -1130,14 +1477,14 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 for ( idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell ) { idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - for ( idx_t icell = 0; icell < iparam.centroids.size(); ++icell ) { - Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.src_weights[icell] ); + for ( idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell ) { + Cs = Cs + PointXYZ::mul( iparam.centroids[icell], iparam.weights[icell] ); } + const double Cs_norm = PointXYZ::norm( Cs ); + ATLAS_ASSERT( Cs_norm > 0. ); + Cs = PointXYZ::div( Cs, Cs_norm ); } - const double Cs_norm = PointXYZ::norm( Cs ); - ATLAS_ASSERT( Cs_norm > 0. ); - Cs = PointXYZ::div( Cs, Cs_norm ); -*/ + */ const PointXYZ& Cs = src_points_[snode]; // compute gradient from nodes double dual_area_inv = 0.; @@ -1168,19 +1515,19 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 for (idx_t isubcell = 0; isubcell < src_node2csp_[snode].size(); ++isubcell) { idx_t subcell = src_node2csp_[snode][isubcell]; const auto& iparam = src_iparam_[subcell]; - if (iparam.centroids.size() == 0) { + if (iparam.cell_idx.size() == 0) { continue; } std::vector Aik; - Aik.resize(iparam.centroids.size()); - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + Aik.resize(iparam.cell_idx.size()); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const PointXYZ& Csk = iparam.centroids[icell]; const PointXYZ Csk_Cs = Csk - Cs; Aik[icell] = Csk_Cs - PointXYZ::mul(Cs, PointXYZ::dot(Cs, Csk_Cs)); Aik[icell] = PointXYZ::mul(Aik[icell], iparam.tgt_weights[icell] * dual_area_inv); } if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { const idx_t tcell = iparam.cell_idx[icell]; for (idx_t j = 0; j < nb_nodes.size(); ++j) { idx_t nj = next_index(j, nb_nodes.size()); @@ -1194,11 +1541,11 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } else { auto& tgt_csp2node_ = data_->tgt_csp2node_; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; double inv_node_weight = (tgt_areas_v[tnode] > 1e-15) ? 1. / tgt_areas_v[tnode] : 0.; - double csp2node_coef = iparam.src_weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; + double csp2node_coef = iparam.weights[icell] / iparam.tgt_weights[icell] * inv_node_weight; for (idx_t j = 0; j < nb_nodes.size(); ++j) { idx_t nj = next_index(j, nb_nodes.size()); idx_t sj = nb_nodes[j]; @@ -1218,6 +1565,16 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 return Matrix(n_tpoints_, n_spoints_, triplets); } +void ConservativeSphericalPolygonInterpolation::do_execute(const FieldSet& src_fields, FieldSet& tgt_fields, + Metadata& metadata) const { + std::vector md_array; + md_array.resize( src_fields.size() ); + for (int i = 0; i < src_fields.size(); i++) { // TODO: memory-wise should we openmp this? + do_execute(src_fields[i], tgt_fields[i], md_array[i]); + } + metadata = md_array[0]; // TODO: reduce metadata of a set of variables to a single metadata of a variable? +} + void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_field, Field& tgt_field, Metadata& metadata) const { ATLAS_TRACE("ConservativeMethod::do_execute()"); @@ -1245,8 +1602,8 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } for (idx_t scell = 0; scell < src_vals.size(); ++scell) { const auto& iparam = src_iparam_[scell]; - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { - tgt_vals(iparam.cell_idx[icell]) += iparam.src_weights[icell] * src_vals(scell); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { + tgt_vals(iparam.cell_idx[icell]) += iparam.weights[icell] * src_vals(scell); } } for (idx_t tcell = 0; tcell < tgt_vals.size(); ++tcell) { @@ -1260,16 +1617,13 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } else if (order_ == 2) { if (matrix_free_) { - ATLAS_TRACE("matrix_free_order_2"); - const auto& src_iparam_ = data_->src_iparam_; - const auto& tgt_areas_v = data_->tgt_areas_; - if (not src_cell_data_ or not tgt_cell_data_) { ATLAS_NOTIMPLEMENTED; } - + ATLAS_TRACE("matrix_free_order_2"); + const auto& src_iparam_ = data_->src_iparam_; + const auto& tgt_areas_v = data_->tgt_areas_; auto& src_points_ = data_->src_points_; - const auto src_vals = array::make_view(src_field); auto tgt_vals = array::make_view(tgt_field); const auto halo = array::make_view(src_mesh_.cells().halo()); @@ -1277,55 +1631,42 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel tgt_vals(tcell) = 0.; } for (idx_t scell = 0; scell < src_vals.size(); ++scell) { - if (halo(scell)) { + const auto& iparam = src_iparam_[scell]; + if (iparam.cell_idx.size() == 0 or halo(scell)) { continue; } - const auto& iparam = src_iparam_[scell]; - const PointXYZ& P = src_points_[scell]; + const PointXYZ& Cs = src_points_[scell]; PointXYZ grad = {0., 0., 0.}; PointXYZ src_barycenter = {0., 0., 0.}; - auto src_neighbour_cells = get_cell_neighbours(src_mesh_, scell); - double dual_area = 0.; - for (idx_t nb_id = 0; nb_id < src_neighbour_cells.size(); ++nb_id) { - idx_t nnb_id = next_index(nb_id, src_neighbour_cells.size()); - idx_t ncell = src_neighbour_cells[nb_id]; - idx_t nncell = src_neighbour_cells[nnb_id]; - const auto& Pn = src_points_[ncell]; - const auto& Pnn = src_points_[nncell]; - if (ncell != scell && nncell != scell) { - double val = 0.5 * (src_vals(ncell) + src_vals(nncell)) - src_vals(scell); - auto csp = ConvexSphericalPolygon({Pn, Pnn, P}); - if (csp.area() < std::numeric_limits::epsilon()) { - csp = ConvexSphericalPolygon({Pn, P, Pnn}); - } - auto NsNsj = ConvexSphericalPolygon::GreatCircleSegment(P, Pn); - val *= (NsNsj.inLeftHemisphere(Pnn, -1e-16) ? -1 : 1); - dual_area += std::abs(csp.area()); - grad = grad + PointXYZ::mul(PointXYZ::cross(Pn, Pnn), val); - } - else if (ncell != scell) { - ATLAS_NOTIMPLEMENTED; - //double val = 0.5 * ( src_vals( ncell ) - src_vals( scell ) ); - //grad = grad + PointXYZ::mul( PointXYZ::cross( Pn, P ), val ); + auto nb_cells = get_cell_neighbours(src_mesh_, scell); + double dual_area_inv = 0.; + for (idx_t j = 0; j < nb_cells.size(); ++j) { + idx_t nj = next_index(j, nb_cells.size()); + idx_t sj = nb_cells[j]; + idx_t nsj = nb_cells[nj]; + const auto& Csj = src_points_[sj]; + const auto& Cnsj = src_points_[nsj]; + double val = 0.5 * (src_vals(nj) + src_vals(nsj)) - src_vals(j); + if (ConvexSphericalPolygon::GreatCircleSegment(Cs, Csj).inLeftHemisphere(Cnsj, -1e-16)) { + dual_area_inv += ConvexSphericalPolygon({Cs, Csj, Cnsj}).area(); } - else if (nncell != scell) { - ATLAS_NOTIMPLEMENTED; - //double val = 0.5 * ( src_vals( nncell ) - src_vals( scell ) ); - //grad = grad + PointXYZ::mul( PointXYZ::cross( P, Pnn ), val ); + else { + //val *= -1; + dual_area_inv += ConvexSphericalPolygon({Cs, Cnsj, Csj}).area(); } + grad = grad + PointXYZ::mul(PointXYZ::cross(Csj, Cnsj), val); } - if (dual_area > std::numeric_limits::epsilon()) { - grad = PointXYZ::div(grad, dual_area); - } - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { - src_barycenter = src_barycenter + PointXYZ::mul(iparam.centroids[icell], iparam.src_weights[icell]); + dual_area_inv = ((dual_area_inv > 0.) ? 1. / dual_area_inv : 0.); + grad = PointXYZ::mul(grad, dual_area_inv); + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { + src_barycenter = src_barycenter + PointXYZ::mul(iparam.centroids[icell], iparam.weights[icell]); } src_barycenter = PointXYZ::div(src_barycenter, PointXYZ::norm(src_barycenter)); grad = grad - PointXYZ::mul(src_barycenter, PointXYZ::dot(grad, src_barycenter)); ATLAS_ASSERT(std::abs(PointXYZ::dot(grad, src_barycenter)) < 1e-14); - for (idx_t icell = 0; icell < iparam.centroids.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.cell_idx.size(); ++icell) { tgt_vals(iparam.cell_idx[icell]) += - iparam.src_weights[icell] * + iparam.weights[icell] * (src_vals(scell) + PointXYZ::dot(grad, iparam.centroids[icell] - src_barycenter)); } } @@ -1340,6 +1681,10 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } stopwatch.stop(); + { + ATLAS_TRACE("halo exchange target"); + tgt_field.haloExchange(); + } auto remap_stat = remap_stat_; if (statistics_conservation_) { @@ -1369,7 +1714,7 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel else { auto& src_node2csp_ = data_->src_node2csp_; for (idx_t spt = 0; spt < src_vals.size(); ++spt) { - if (src_node_ghost(spt) or src_areas_v[spt] < 1e-14) { + if (src_node_halo(spt) or src_node_ghost(spt)) { continue; } err_remap_cons += src_vals(spt) * src_areas_v[spt]; @@ -1386,24 +1731,28 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel } else { for (idx_t tpt = 0; tpt < tgt_vals.size(); ++tpt) { - if (tgt_node_ghost(tpt)) { + if (tgt_node_halo(tpt) or tgt_node_ghost(tpt)) { continue; } err_remap_cons -= tgt_vals(tpt) * tgt_areas_v[tpt]; } } ATLAS_TRACE_MPI(ALLREDUCE) { mpi::comm().allReduceInPlace(&err_remap_cons, 1, eckit::mpi::sum()); } - remap_stat.errors[Statistics::Errors::REMAP_CONS] = std::sqrt(std::abs(err_remap_cons) / unit_sphere_area()); + remap_stat.errors[Statistics::Errors::REMAP_CONS] = err_remap_cons / unit_sphere_area(); metadata.set("conservation_error", remap_stat.errors[Statistics::Errors::REMAP_CONS]); } if (statistics_intersection_) { - metadata.set("polygons.source", remap_stat.counts[Statistics::SRC_PLG]); - metadata.set("polygons.target", remap_stat.counts[Statistics::TGT_PLG]); - metadata.set("polygons.intersections", remap_stat.counts[Statistics::INT_PLG]); + metadata.set("polygons.source", remap_stat.counts[Statistics::Counts::SRC_PLG]); + metadata.set("polygons.target", remap_stat.counts[Statistics::Counts::TGT_PLG]); + metadata.set("polygons.intersections", remap_stat.counts[Statistics::Counts::INT_PLG]); metadata.set("polygons.uncovered_source", remap_stat.counts[Statistics::UNCVR_SRC]); - metadata.set("source_area_error.L1", remap_stat.errors[Statistics::Errors::GEO_L1]); - metadata.set("source_area_error.Linf", remap_stat.errors[Statistics::Errors::GEO_LINF]); + if (validate_) { + metadata.set("source_area_error.L1", remap_stat.errors[Statistics::Errors::SRC_INTERSECTPLG_L1]); + metadata.set("source_area_error.Linf", remap_stat.errors[Statistics::Errors::SRC_INTERSECTPLG_LINF]); + metadata.set("target_area_error.L1", remap_stat.errors[Statistics::Errors::TGT_INTERSECTPLG_L1]); + metadata.set("target_area_error.Linf", remap_stat.errors[Statistics::Errors::TGT_INTERSECTPLG_LINF]); + } } if (statistics_intersection_ || statistics_conservation_) { @@ -1430,15 +1779,16 @@ void ConservativeSphericalPolygonInterpolation::do_execute(const Field& src_fiel metadata.set("memory.src_node2csp", memory_of(data_->src_node2csp_)); metadata.set("memory.tgt_node2csp", memory_of(data_->tgt_node2csp_)); metadata.set("memory.src_iparam", memory_of(data_->src_iparam_)); - - tgt_field.set_dirty(); } void ConservativeSphericalPolygonInterpolation::print(std::ostream& out) const { out << "ConservativeMethod{"; out << "order:" << order_; - out << ", source:" << (src_cell_data_ ? "cells" : "nodes"); - out << ", target:" << (tgt_cell_data_ ? "cells" : "nodes"); + int halo = 0; + src_mesh_.metadata().get("halo", halo); + out << ", source:" << (src_cell_data_ ? "cells(" : "nodes(") << src_mesh_.grid().name() << ",halo=" << halo << ")"; + tgt_mesh_.metadata().get("halo", halo); + out << ", target:" << (tgt_cell_data_ ? "cells(" : "nodes(") << tgt_mesh_.grid().name() << ",halo=" << halo << ")"; out << ", normalise_intersections:" << normalise_intersections_; out << ", matrix_free:" << matrix_free_; out << ", statistics.intersection:" << statistics_intersection_; @@ -1506,7 +1856,7 @@ void ConservativeSphericalPolygonInterpolation::setup_stat() const { remap_stat_.tgt_area_sum = src_tgt_sums[1]; geo_create_err = std::abs(src_tgt_sums[0] - src_tgt_sums[1]) / unit_sphere_area(); - remap_stat_.errors[Statistics::Errors::GEO_DIFF] = geo_create_err; + remap_stat_.errors[Statistics::Errors::SRCTGT_INTERSECTPLG_DIFF] = geo_create_err; } Field ConservativeSphericalPolygonInterpolation::Statistics::diff(const Interpolation& interpolation, @@ -1543,29 +1893,34 @@ Field ConservativeSphericalPolygonInterpolation::Statistics::diff(const Interpol double diff = src_vals(spt) * src_areas_v[spt]; const auto& iparam = src_iparam_[spt]; if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; if (tgt_cell_halo(tcell) < 1) { - diff -= tgt_vals(iparam.cell_idx[icell]) * iparam.src_weights[icell]; + diff -= tgt_vals(tcell) * iparam.weights[icell]; } } } else { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; - if (tgt_node_halo(tnode) < 1) { - diff -= tgt_vals(tnode) * iparam.src_weights[icell]; + if (not tgt_node_ghost(tnode)) { + diff -= tgt_vals(tnode) * iparam.weights[icell]; } } } - diff_vals(spt) = std::abs(diff) / src_areas_v[spt]; + if (src_areas_v[spt] > 0.) { + diff_vals(spt) = diff / src_areas_v[spt]; + } + else { + Log::info() << " at cell " << spt << " cell-area: " << src_areas_v[spt] << std::endl; + diff_vals(spt) = std::numeric_limits::max(); + } } } else { for (idx_t spt = 0; spt < src_vals.size(); ++spt) { - if (src_node_ghost(spt) or src_areas_v[spt] < 1e-14) { - diff_vals(spt) = 0.; + if (src_node_ghost(spt) or src_node_halo(spt) != 0) { continue; } double diff = src_vals(spt) * src_areas_v[spt]; @@ -1573,38 +1928,50 @@ Field ConservativeSphericalPolygonInterpolation::Statistics::diff(const Interpol for (idx_t subcell = 0; subcell < node2csp.size(); ++subcell) { const auto& iparam = src_iparam_[node2csp[subcell]]; if (tgt_cell_data_) { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { - diff -= tgt_vals(iparam.cell_idx[icell]) * iparam.src_weights[icell]; + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { + idx_t tcell = iparam.cell_idx[icell]; + if (tgt_cell_halo(tcell) < 1) { + diff -= tgt_vals(tcell) * iparam.weights[icell]; + } } } else { - for (idx_t icell = 0; icell < iparam.src_weights.size(); ++icell) { + for (idx_t icell = 0; icell < iparam.weights.size(); ++icell) { idx_t tcell = iparam.cell_idx[icell]; idx_t tnode = tgt_csp2node_[tcell]; - diff -= tgt_vals(tnode) * iparam.src_weights[icell]; + if (not tgt_node_ghost(tnode)) { + diff -= tgt_vals(tnode) * iparam.weights[icell]; + } } } } - diff_vals(spt) = std::abs(diff) / src_areas_v[spt]; + if (src_areas_v[spt] > 0.) { + diff_vals(spt) = diff / src_areas_v[spt]; + } + else { + diff_vals(spt) = std::numeric_limits::max(); + Log::info() << " at cell " << spt << " cell-area: " << src_areas_v[spt] << std::endl; + } } } return diff; } -void ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Interpolation& interpolation, - const Field target, - std::function func) { +ConservativeSphericalPolygonInterpolation::Metadata ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Interpolation& interpolation, + const Field target, + std::function func) { auto tgt_vals = array::make_view(target); auto cachable_data_ = ConservativeSphericalPolygonInterpolation::Cache(interpolation).get(); auto tgt_mesh_ = extract_mesh(cachable_data_->src_fs_); auto tgt_cell_data_ = extract_mesh(cachable_data_->tgt_fs_); const auto tgt_cell_halo = array::make_view(tgt_mesh_.cells().halo()); const auto tgt_node_ghost = array::make_view(tgt_mesh_.nodes().ghost()); + const auto tgt_node_halo = array::make_view(tgt_mesh_.nodes().halo()); const auto& tgt_areas_v = cachable_data_->tgt_areas_; + auto& tgt_points_ = cachable_data_->tgt_points_; double err_remap_l2 = 0.; double err_remap_linf = 0.; - auto& tgt_points_ = cachable_data_->tgt_points_; if (tgt_cell_data_) { size_t ncells = std::min(tgt_vals.size(), tgt_mesh_.cells().size()); for (idx_t tpt = 0; tpt < ncells; ++tpt) { @@ -1623,7 +1990,7 @@ void ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Inter else { size_t nnodes = std::min(tgt_vals.size(), tgt_mesh_.nodes().size()); for (idx_t tpt = 0; tpt < nnodes; ++tpt) { - if (tgt_node_ghost(tpt)) { + if (tgt_node_ghost(tpt) or tgt_node_halo(tpt)) { continue; } auto p = tgt_points_[tpt]; @@ -1638,77 +2005,71 @@ void ConservativeSphericalPolygonInterpolation::Statistics::accuracy(const Inter mpi::comm().allReduceInPlace(&err_remap_l2, 1, eckit::mpi::sum()); mpi::comm().allReduceInPlace(&err_remap_linf, 1, eckit::mpi::max()); } - this->errors[Statistics::Errors::REMAP_L2] = std::sqrt(err_remap_l2 / unit_sphere_area()); - this->errors[Statistics::Errors::REMAP_LINF] = err_remap_linf; + errors[Statistics::Errors::REMAP_L2] = std::sqrt(err_remap_l2 / unit_sphere_area()); + errors[Statistics::Errors::REMAP_LINF] = err_remap_linf; + Metadata metadata; + metadata.set("errors.REMAP_L2", errors[Statistics::Errors::REMAP_L2]); + metadata.set("errors.REMAP_LINF", errors[Statistics::Errors::REMAP_LINF]); + return metadata; } -auto debug_intersection = [](const ConvexSphericalPolygon& plg_1, const ConvexSphericalPolygon& plg_2, - const ConvexSphericalPolygon& iplg, const ConvexSphericalPolygon& jplg) { - const double intersection_comm_err = std::abs(iplg.area() - jplg.area()) / (plg_1.area() > 0 ? plg_1.area() : 1.); - Log::info().indent(); - if (intersection_comm_err > 1e-6) { - Log::info() << "PLG_1 : " << std::setprecision(10) << plg_1 << "\n"; - Log::info() << "area(PLG_1) : " << plg_1.area() << "\n"; - Log::info() << "PLG_2 :" << plg_2 << "\n"; - Log::info() << "PLG_12 : " << iplg << "\n"; - Log::info() << "PLG_21 : " << jplg << "\n"; - Log::info() << "area(PLG_12 - PLG_21) : " << intersection_comm_err << "\n"; - Log::info() << "area(PLG_21) : " << jplg.area() << "\n"; - //ATLAS_ASSERT( false, "SRC.intersect.TGT =/= TGT.intersect.SRC."); - } - int pout; - int pin = inside_vertices(plg_1, plg_2, pout); - if (pin > 2 && iplg.area() < 3e-16) { - Log::info() << " pin : " << pin << ", pout :" << pout << ", total vertices : " << plg_2.size() << "\n"; - Log::info() << "PLG_2 :" << plg_2 << "\n"; - Log::info() << "PLG_12 : " << iplg << "\n"; - Log::info() << "area(PLG_12) : " << iplg.area() << "\n\n"; - //ATLAS_ASSERT( false, "SRC must intersect TGT." ); - } - Log::info().unindent(); -}; - -void ConservativeSphericalPolygonInterpolation::dump_intersection(const ConvexSphericalPolygon& plg_1, +void ConservativeSphericalPolygonInterpolation::dump_intersection(const std::string msg, + const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const std::vector& plg_2_idx_array) const { +#if PRINT_BAD_POLYGONS double plg_1_coverage = 0.; + std::vector int_areas; + int_areas.resize(plg_2_idx_array.size()); for (int i = 0; i < plg_2_idx_array.size(); ++i) { const auto plg_2_idx = plg_2_idx_array[i]; const auto& plg_2 = std::get<0>(plg_2_array[plg_2_idx]); - auto iplg = plg_1.intersect(plg_2); - auto jplg = plg_2.intersect(plg_1); - debug_intersection(plg_1, plg_2, iplg, jplg); + auto iplg = plg_1.intersect(plg_2, false, 1.e5 * std::numeric_limits::epsilon()); plg_1_coverage += iplg.area(); + int_areas[i] = iplg.area(); } - Log::info().indent(); if (std::abs(plg_1.area() - plg_1_coverage) > 0.01 * plg_1.area()) { - Log::info() << "Polygon coverage incomplete. Printing polygons." << std::endl; - Log::info() << "Polygon 1 : "; + Log::info().indent(); + Log::info() << msg << ", Polygon_1.area: " << plg_1.area() << ", covered: " << plg_1_coverage << std::endl; + Log::info() << "Polygon_1 : "; + Log::info().precision(18); plg_1.print(Log::info()); Log::info() << std::endl << "Printing " << plg_2_idx_array.size() << " covering polygons -->" << std::endl; Log::info().indent(); + plg_1_coverage = plg_1.area(); for (int i = 0; i < plg_2_idx_array.size(); ++i) { const auto plg_2_idx = plg_2_idx_array[i]; const auto& plg_2 = std::get<0>(plg_2_array[plg_2_idx]); - Log::info() << "Polygon " << i + 1 << " : "; plg_2.print(Log::info()); + Log::info() << "\ninters:"; + plg_1_coverage -= int_areas[i]; + auto iplg = plg_1.intersect(plg_2, false, 1.e5 * std::numeric_limits::epsilon()); + Log::info().indent(); + iplg.print(Log::info()); + Log::info() << ", inters.-area: " << iplg.area() << ", plg1-area left: " << plg_1_coverage << std::endl; + Log::info().unindent(); Log::info() << std::endl; } + Log::info() << std::endl; + Log::info().unindent(); Log::info().unindent(); } - Log::info().unindent(); +#endif } template -void ConservativeSphericalPolygonInterpolation::dump_intersection(const ConvexSphericalPolygon& plg_1, +void ConservativeSphericalPolygonInterpolation::dump_intersection(const std::string msg, + const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const TargetCellsIDs& plg_2_idx_array) const { +#if PRINT_BAD_POLYGONS std::vector idx_array; idx_array.resize(plg_2_idx_array.size()); for (int i = 0; i < plg_2_idx_array.size(); ++i) { idx_array[i] = plg_2_idx_array[i].payload(); } - dump_intersection(plg_1, plg_2_array, idx_array); + dump_intersection(msg, plg_1, plg_2_array, idx_array); +#endif } ConservativeSphericalPolygonInterpolation::Cache::Cache(std::shared_ptr entry): @@ -1749,20 +2110,16 @@ void ConservativeSphericalPolygonInterpolation::Data::print(std::ostream& out) c } void ConservativeSphericalPolygonInterpolation::Statistics::fillMetadata(Metadata& metadata) { - // counts - metadata.set("counts.SRC_PLG", counts[SRC_PLG]); - metadata.set("counts.TGT_PLG", counts[TGT_PLG]); - metadata.set("counts.INT_PLG", counts[INT_PLG]); - metadata.set("counts.UNCVR_SRC", counts[UNCVR_SRC]); - // errors - metadata.set("errors.SRC_PLG_L1", errors[SRC_PLG_L1]); - metadata.set("errors.SRC_PLG_LINF", errors[SRC_PLG_LINF]); - metadata.set("errors.TGT_PLG_L1", errors[TGT_PLG_L1]); - metadata.set("errors.TGT_PLG_LINF", errors[TGT_PLG_LINF]); - metadata.set("errors.GEO_L1", errors[GEO_L1]); - metadata.set("errors.GEO_LINF", errors[GEO_LINF]); - metadata.set("errors.GEO_DIFF", errors[GEO_DIFF]); + metadata.set("errors.SRC_SUBPLG_L1", errors[SRC_SUBPLG_L1]); + metadata.set("errors.SRC_SUBPLG_LINF", errors[SRC_SUBPLG_LINF]); + metadata.set("errors.TGT_SUBPLG_L1", errors[TGT_SUBPLG_L1]); + metadata.set("errors.TGT_SUBPLG_LINF", errors[TGT_SUBPLG_LINF]); + metadata.set("errors.SRC_INTERSECTPLG_L1", errors[SRC_INTERSECTPLG_L1]); + metadata.set("errors.SRC_INTERSECTPLG_LINF", errors[SRC_INTERSECTPLG_LINF]); + metadata.set("errors.TGT_INTERSECTPLG_L1", errors[TGT_INTERSECTPLG_L1]); + metadata.set("errors.TGT_INTERSECTPLG_LINF", errors[TGT_INTERSECTPLG_LINF]); + metadata.set("errors.SRCTGT_INTERSECTPLG_DIFF", errors[SRCTGT_INTERSECTPLG_DIFF]); metadata.set("errors.REMAP_CONS", errors[REMAP_CONS]); metadata.set("errors.REMAP_L2", errors[REMAP_L2]); metadata.set("errors.REMAP_LINF", errors[REMAP_LINF]); @@ -1774,20 +2131,16 @@ ConservativeSphericalPolygonInterpolation::Statistics::Statistics() { } ConservativeSphericalPolygonInterpolation::Statistics::Statistics(const Metadata& metadata): Statistics() { - // counts - metadata.get("counts.SRC_PLG", counts[SRC_PLG]); - metadata.get("counts.TGT_PLG", counts[TGT_PLG]); - metadata.get("counts.INT_PLG", counts[INT_PLG]); - metadata.get("counts.UNCVR_SRC", counts[UNCVR_SRC]); - // errors - metadata.get("errors.SRC_PLG_L1", errors[SRC_PLG_L1]); - metadata.get("errors.SRC_PLG_LINF", errors[SRC_PLG_LINF]); - metadata.get("errors.TGT_PLG_L1", errors[TGT_PLG_L1]); - metadata.get("errors.TGT_PLG_LINF", errors[TGT_PLG_LINF]); - metadata.get("errors.GEO_L1", errors[GEO_L1]); - metadata.get("errors.GEO_LINF", errors[GEO_LINF]); - metadata.get("errors.GEO_DIFF", errors[GEO_DIFF]); + metadata.get("errors.SRC_SUBPLG_L1", errors[SRC_SUBPLG_L1]); + metadata.get("errors.SRC_SUBPLG_LINF", errors[SRC_SUBPLG_LINF]); + metadata.get("errors.TGT_SUBPLG_L1", errors[TGT_SUBPLG_L1]); + metadata.get("errors.TGT_SUBPLG_LINF", errors[TGT_SUBPLG_LINF]); + metadata.get("errors.SRC_INTERSECTPLG_L1", errors[SRC_INTERSECTPLG_L1]); + metadata.get("errors.SRC_INTERSECTPLG_LINF", errors[SRC_INTERSECTPLG_LINF]); + metadata.get("errors.TGT_INTERSECTPLG_L1", errors[TGT_INTERSECTPLG_L1]); + metadata.get("errors.TGT_INTERSECTPLG_LINF", errors[TGT_INTERSECTPLG_LINF]); + metadata.get("errors.SRCTGT_INTERSECTPLG_DIFF", errors[SRCTGT_INTERSECTPLG_DIFF]); metadata.get("errors.REMAP_CONS", errors[REMAP_CONS]); metadata.get("errors.REMAP_L2", errors[REMAP_L2]); metadata.get("errors.REMAP_LINF", errors[REMAP_LINF]); diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h index bcab9fdba..3dd79099b 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h @@ -25,10 +25,8 @@ class ConservativeSphericalPolygonInterpolation : public Method { struct InterpolationParameters { // one polygon intersection std::vector cell_idx; // target cells used for intersection std::vector centroids; // intersection cell centroids - std::vector src_weights; // intersection cell areas - - // TODO: tgt_weights can be computed on the fly - std::vector tgt_weights; // (intersection cell areas) / (target cell area) + std::vector weights; // intersection cell areas + std::vector tgt_weights; // (intersection cell areas) / (target cell area) // TODO: tgt_weights vector should be removed for the sake of highres }; private: @@ -55,7 +53,6 @@ class ConservativeSphericalPolygonInterpolation : public Method { std::vector> src_node2csp_; std::vector> tgt_node2csp_; - // Timings struct Timings { double source_polygons_assembly{0}; @@ -102,18 +99,20 @@ class ConservativeSphericalPolygonInterpolation : public Method { std::array counts; enum Errors { - SRC_PLG_L1 = 0, // index, over/undershoot in source subpolygon creation - SRC_PLG_LINF, - TGT_PLG_L1, // index, over/untershoot in target subpolygon creation - TGT_PLG_LINF, - GEO_L1, // index, cumulative area mismatch in polygon intersections - GEO_LINF, // index, like GEO_L1 but in L_infinity norm - GEO_DIFF, // index, difference in earth area coverages + SRC_SUBPLG_L1 = 0, // index, \sum_{cell of mesh} {cell.area - \sum_{subpol of cell} subpol.area} + SRC_SUBPLG_LINF, // index, \max_-||- + TGT_SUBPLG_L1, // see above + TGT_SUBPLG_LINF, // see above + SRC_INTERSECTPLG_L1, // index, \sum_{cell of src-mesh} {cell.area - \sum_{intersect of cell} intersect.area} + SRC_INTERSECTPLG_LINF, // index, \max_-||- + TGT_INTERSECTPLG_L1, // see above + TGT_INTERSECTPLG_LINF, // see above + SRCTGT_INTERSECTPLG_DIFF, // index, 1/(unit_sphere.area) ( \sum_{scell} scell.area - \sum{tcell} tcell.area ) REMAP_CONS, // index, error in mass conservation REMAP_L2, // index, error accuracy for given analytical function REMAP_LINF // index, like REMAP_L2 but in L_infinity norm }; - std::array errors; + std::array errors; double tgt_area_sum; double src_area_sum; @@ -123,8 +122,8 @@ class ConservativeSphericalPolygonInterpolation : public Method { Statistics(); Statistics(const Metadata&); - void accuracy(const Interpolation& interpolation, const Field target, - std::function func); + Metadata accuracy(const Interpolation& interpolation, const Field target, + std::function func); // compute difference field of source and target mass @@ -139,6 +138,7 @@ class ConservativeSphericalPolygonInterpolation : public Method { void do_setup(const FunctionSpace& src_fs, const FunctionSpace& tgt_fs) override; void do_setup(const Grid& src_grid, const Grid& tgt_grid, const interpolation::Cache&) override; void do_execute(const Field& src_field, Field& tgt_field, Metadata&) const override; + void do_execute(const FieldSet& src_fields, FieldSet& tgt_fields, Metadata&) const override; void print(std::ostream& out) const override; @@ -160,17 +160,17 @@ class ConservativeSphericalPolygonInterpolation : public Method { void intersect_polygons(const CSPolygonArray& src_csp, const CSPolygonArray& tgt_scp); Matrix compute_1st_order_matrix(); Matrix compute_2nd_order_matrix(); - void dump_intersection(const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, + void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const std::vector& plg_2_idx_array) const; template - void dump_intersection(const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, + void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const TargetCellsIDs& plg_2_idx_array) const; std::vector sort_cell_edges(Mesh& mesh, idx_t cell_id) const; std::vector sort_node_edges(Mesh& mesh, idx_t cell_id) const; - std::vector get_cell_neighbours(Mesh& mesh, idx_t jcell) const; - std::vector get_node_neighbours(Mesh& mesh, idx_t jcell) const; - CSPolygonArray get_polygons_celldata(Mesh& mesh) const; - CSPolygonArray get_polygons_nodedata(Mesh& mesh, std::vector& csp2node, + std::vector get_cell_neighbours(Mesh&, idx_t jcell) const; + std::vector get_node_neighbours(Mesh&, idx_t jcell) const; + CSPolygonArray get_polygons_celldata(FunctionSpace) const; + CSPolygonArray get_polygons_nodedata(FunctionSpace, std::vector& csp2node, std::vector>& node2csp, std::array& errors) const; diff --git a/src/atlas/util/ConvexSphericalPolygon.cc b/src/atlas/util/ConvexSphericalPolygon.cc index fe266f87f..40b3584c7 100644 --- a/src/atlas/util/ConvexSphericalPolygon.cc +++ b/src/atlas/util/ConvexSphericalPolygon.cc @@ -19,9 +19,9 @@ #include "atlas/util/ConvexSphericalPolygon.h" #include "atlas/util/CoordinateEnums.h" #include "atlas/util/NormaliseLongitude.h" +#include "atlas/library/FloatingPointExceptions.h" -#define DEBUG_OUTPUT 0 -#define DEBUG_OUTPUT_DETAIL 0 +// #define DEBUG_OUTPUT 1 namespace atlas { namespace util { @@ -32,14 +32,6 @@ namespace { constexpr double EPS = std::numeric_limits::epsilon(); constexpr double EPS2 = EPS * EPS; -constexpr double TOL = 1.e4 * EPS; // two points considered "same" -constexpr double TOL2 = TOL * TOL; - -enum IntersectionType -{ - NO_INTERSECT = -100, - OVERLAP -}; double distance2(const PointXYZ& p1, const PointXYZ& p2) { double dx = p2[0] - p1[0]; @@ -52,38 +44,35 @@ double norm2(const PointXYZ& p) { return p[0] * p[0] + p[1] * p[1] + p[2] * p[2]; } -bool approx_eq(const double& v1, const double& v2, const double& tol) { - return std::abs(v1 - v2) < tol; +bool approx_eq(const double& v1, const double& v2) { + return std::abs(v1 - v2) <= EPS; } -bool approx_eq(const PointXYZ& v1, const PointXYZ& v2, const double& tol) { +bool approx_eq(const PointXYZ& v1, const PointXYZ& v2) { //return approx_eq( v1[0], v2[0], t ) && approx_eq( v1[1], v2[1], t ) && approx_eq( v1[2], v2[2], t ); - return distance2(v1, v2) < tol * tol; -} - -bool approx_eq_null(const PointXYZ& v1, const double& tol) { - //return approx_eq( v1[0], 0., t ) && approx_eq( v1[1], 0., t ) && approx_eq( v1[2], 0., t ); - double n = v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]; - return n < tol * tol; + return distance2(v1, v2) <= EPS2; } void lonlat2xyz(const PointLonLat& lonlat, PointXYZ& xyz) { eckit::geometry::Sphere::convertSphericalToCartesian(1., lonlat, xyz); } -void xyz2lonlat(const PointXYZ xyz, PointLonLat& lonlat) { +void xyz2lonlat(const PointXYZ& xyz, PointLonLat& lonlat) { eckit::geometry::Sphere::convertCartesianToSpherical(1., xyz, lonlat); } -double norm_max(const PointXYZ& p, const PointXYZ& q) { - double n01 = std::max(std::abs(p[0] - q[0]), std::abs(p[1] - q[1])); - return std::max(n01, std::abs(p[2] - q[2])); +PointLonLat xyz2lonlat(const PointXYZ& xyz) { + PointLonLat lonlat; + eckit::geometry::Sphere::convertCartesianToSpherical(1., xyz, lonlat); + return lonlat; } -template +#if 0 +// NOTE: StackVector is not used +template struct StackVector { private: - using Wrapped = std::array; + using Wrapped = std::array; public: using reference = typename Wrapped::reference; @@ -101,7 +90,7 @@ struct StackVector { #endif void push_back(const T& value) { wrapped_[size_++] = value; - ATLAS_ASSERT(size_ < ConvexSphericalPolygon::MAX_SIZE); + ATLAS_ASSERT(size_ < MAX_SIZE); } size_t size() const { return size_; } @@ -109,38 +98,7 @@ struct StackVector { size_t size_{0}; Wrapped wrapped_; }; - -struct PolygonEdgeIntersection { - static constexpr int BEGIN = 1; - static constexpr int END = 2; - static constexpr int INSIDE = 3; - - PolygonEdgeIntersection(const ConvexSphericalPolygon& polygon, int edge_index, const PointXYZ& point) { - auto matches = [](const PointXYZ& p1, const PointXYZ& p2) { - return (distance2(p1, p2) < 1e-16); - // We would like this to be TOL2 instead, but that gives bad results - }; - - ATLAS_ASSERT(edge_index >= 0); - ATLAS_ASSERT(edge_index < polygon.size()); - - const int node_index = edge_index; - - if (matches(point, polygon[node_index])) { - location = BEGIN; - } - else if (matches(point, polygon[polygon.next(node_index)])) { - location = END; - } - else { - location = INSIDE; - } - } - bool isPointAtBegin() const { return location == BEGIN; } - bool isPointAtEnd() const { return location == END; } - bool isPointInside() const { return location == INSIDE; } - int location; -}; +#endif } // namespace @@ -153,13 +111,14 @@ ConvexSphericalPolygon::ConvexSphericalPolygon(const PointLonLat points[], size_ size_t isp = 1; for (size_t i = 1; i < size_ - 1; ++i) { lonlat2xyz(points[i], sph_coords_[isp]); - if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { +// Log::info() << " d : " << PointXYZ::distance(sph_coords_[isp], sph_coords_[isp - 1]) << std::endl; + if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { continue; } ++isp; } lonlat2xyz(points[size_ - 1], sph_coords_[isp]); - if (approx_eq(sph_coords_[isp], sph_coords_[0], TOL) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { + if (approx_eq(sph_coords_[isp], sph_coords_[0]) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { } else { ++isp; @@ -178,13 +137,13 @@ ConvexSphericalPolygon::ConvexSphericalPolygon(const PointXYZ points[], size_t s size_t isp = 1; for (size_t i = 1; i < size_ - 1; ++i) { sph_coords_[isp] = points[i]; - if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { + if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { continue; } ++isp; } sph_coords_[isp] = points[size_ - 1]; - if (approx_eq(sph_coords_[isp], sph_coords_[0], TOL) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { + if (approx_eq(sph_coords_[isp], sph_coords_[0]) or approx_eq(sph_coords_[isp], sph_coords_[isp - 1])) { } else { ++isp; @@ -196,74 +155,8 @@ ConvexSphericalPolygon::ConvexSphericalPolygon(const PointXYZ points[], size_t s } } - -void ConvexSphericalPolygon::simplify() { - ATLAS_ASSERT(size_ < MAX_SIZE); - if (size_ < 3) { - size_ = 0; - valid_ = false; - return; - } - idx_t isp = 0; - idx_t i = 0; - idx_t j; - idx_t k; - bool search_3pts = true; - auto& points = sph_coords_; - for (; i < size_ && search_3pts; ++i) { - const PointXYZ& P0 = points[i]; - for (j = i + 1; j < size_ && search_3pts; ++j) { - const PointXYZ& P1 = points[j]; - if (approx_eq(P0, P1, 1.e-10)) { - continue; - } - for (k = j + 1; k < size_ && search_3pts; ++k) { - const PointXYZ& P2 = points[k]; - if (approx_eq(P1, P2, TOL) or approx_eq(P0, P2, TOL)) { - continue; - } - if (GreatCircleSegment{P0, P1}.inLeftHemisphere(P2, -EPS)) { - sph_coords_[isp++] = P0; - sph_coords_[isp++] = P1; - sph_coords_[isp++] = P2; - search_3pts = false; - } - } - } - } - if (search_3pts) { - valid_ = false; - size_ = 0; - return; - } - for (; k < size_ - 1; ++k) { - if (approx_eq(points[k], sph_coords_[isp - 1], TOL) or - (not GreatCircleSegment{sph_coords_[isp - 2], sph_coords_[isp - 1]}.inLeftHemisphere(points[k], -EPS))) { - continue; - } - sph_coords_[isp] = points[k]; - isp++; - } - const PointXYZ& Pl2 = sph_coords_[isp - 2]; - const PointXYZ& Pl1 = sph_coords_[isp - 1]; - const PointXYZ& P0 = sph_coords_[0]; - const PointXYZ& P = points[size_ - 1]; - if ((not approx_eq(P, P0, EPS)) and (not approx_eq(P, Pl1, EPS)) and - GreatCircleSegment{Pl2, Pl1}.inLeftHemisphere(P, -EPS) and - GreatCircleSegment{Pl1, P}.inLeftHemisphere(P0, -EPS)) { - sph_coords_[isp] = P; - ++isp; - } - size_ = isp; - valid_ = size_ > 2; - - computed_area_ = false; - computed_radius_ = false; - computed_centroid_ = false; -} - void ConvexSphericalPolygon::compute_centroid() const { - const auto triangles = triangulate(radius()); + const auto triangles = triangulate(); area_ = triangles.area(); computed_area_ = true; @@ -288,9 +181,9 @@ bool ConvexSphericalPolygon::validate() { int nni = next(ni); const PointXYZ& P = sph_coords_[i]; const PointXYZ& nextP = sph_coords_[ni]; - ATLAS_ASSERT(std::abs(PointXYZ::dot(P, P) - 1.) < 10. * EPS); - ATLAS_ASSERT(not approx_eq(P, PointXYZ::mul(nextP, -1.), TOL)); - valid_ = valid_ && GreatCircleSegment{P, nextP}.inLeftHemisphere(sph_coords_[nni], -EPS); + ATLAS_ASSERT(std::abs(PointXYZ::dot(P, P) - 1.) < 10 * EPS); + ATLAS_ASSERT(not approx_eq(P, PointXYZ::mul(nextP, -1.))); + valid_ = valid_ && GreatCircleSegment{P, nextP}.inLeftHemisphere(sph_coords_[nni], -0.5*EPS); } } return valid_; @@ -300,8 +193,16 @@ bool ConvexSphericalPolygon::equals(const ConvexSphericalPolygon& plg, const dou if (size_ == 0 and plg.size_ == 0) { return true; } - if ((not plg.valid_) || (not valid_) || size_ != plg.size()) { - Log::info() << " ConvexSphericalPolygon::equals == not compatible\n"; + if (not valid_) { + Log::info() << " ConvexSphericalPolygon::equals : this polygon is not valid\n"; + return false; + } + if (not plg.valid_) { + Log::info() << " ConvexSphericalPolygon::equals : other polygon passed as argument is not valid\n"; + return false; + } + if (size_ != plg.size()) { + Log::info() << " ConvexSphericalPolygon::equals : incompatible sizes: " << size_ << " != " << plg.size() << "\n"; return false; } int offset = 0; @@ -313,7 +214,7 @@ bool ConvexSphericalPolygon::equals(const ConvexSphericalPolygon& plg, const dou } } if (offset == size_) { - Log::info() << "ConvexSphericalPolygon::equals == no point equal\n"; + Log::info() << "ConvexSphericalPolygon::equals : no point equal\n"; return false; } @@ -321,82 +222,38 @@ bool ConvexSphericalPolygon::equals(const ConvexSphericalPolygon& plg, const dou int idx = (offset + j) % size_; auto dist2 = distance2(plg.sph_coords_[j], sph_coords_[idx]); if (dist2 > le2) { - Log::info() << " ConvexSphericalPolygon::equals == point distance " << std::sqrt(dist2) << "\n"; + Log::info() << " ConvexSphericalPolygon::equals : point distance " << std::sqrt(dist2) << " < " << le2 << "(= "<< deg_prec << " deg)" << "\n"; return false; } } return true; } -// note: unit sphere! -// I. Todhunter (1886), Paragr. 99 -ConvexSphericalPolygon::SubTriangles ConvexSphericalPolygon::triangulate(const double cell_radius) const { +// cf. Folke Eriksson, "On the Measure of Solid Angles", Mathematics Magazine, Vol. 63, No. 3, pp. 184-187 (1990) +ConvexSphericalPolygon::SubTriangles ConvexSphericalPolygon::triangulate() const { SubTriangles triangles; if (size_ < 3) { return triangles; } - size_t itri{0}; - if (cell_radius < 1.e-6) { // plane area - for (int i = 1; i < size_ - 1; i++) { - const PointXYZ pl = sph_coords_[i] - sph_coords_[0]; - const PointXYZ pr = sph_coords_[i + 1] - sph_coords_[0]; - triangles[itri].centroid = PointXYZ::normalize(sph_coords_[0] + sph_coords_[i] + sph_coords_[i + 1]); - triangles[itri].area = 0.5 * PointXYZ::norm(PointXYZ::cross(pl, pr)); - ++itri; - } - } - else { // spherical area - const PointXYZ& a = sph_coords_[0]; - for (size_t i = 1; i < size_ - 1; i++) { - const PointXYZ& b = sph_coords_[i]; - const PointXYZ& c = sph_coords_[i + 1]; - auto ab = PointXYZ::cross(a, b); - auto bc = PointXYZ::cross(b, c); - auto ca = PointXYZ::cross(c, a); - const double ab_norm = PointXYZ::norm(ab); - const double bc_norm = PointXYZ::norm(bc); - const double ca_norm = PointXYZ::norm(ca); - if (ab_norm < EPS or bc_norm < EPS or ca_norm < EPS) { - continue; - } - double abc = -PointXYZ::dot(ab, bc) / (ab_norm * bc_norm); - double bca = -PointXYZ::dot(bc, ca) / (bc_norm * ca_norm); - double cab = -PointXYZ::dot(ca, ab) / (ca_norm * ab_norm); - if (abc <= -1.) { - abc = M_PI; - } - else if (abc < 1.) { - abc = std::acos(abc); - } - else { - abc = 0.; - } - if (bca <= -1.) { - bca = M_PI; - } - else if (bca < 1.) { - bca = std::acos(bca); - } - else { - bca = 0.; - } - if (cab <= -1.) { - cab = M_PI; - } - else if (cab < 1.) { - cab = std::acos(cab); - } - else { - cab = 0.; - } - triangles[itri].centroid = PointXYZ::normalize(a + b + c); - triangles[itri].area = abc + bca + cab - M_PI; - ++itri; - } + const PointXYZ& a = sph_coords_[0]; + for (size_t i = 1; i < size_ - 1; i++) { + const PointXYZ& b = sph_coords_[i]; + const PointXYZ& c = sph_coords_[i + 1]; + triangles[itri].centroid = PointXYZ::normalize(a + b + c); +// if (PointXYZ::distance(a, b) + PointXYZ::distance(b, c) < 1e-10) { +// triangles[itri].area = 0.5 * PointXYZ::norm(PointXYZ::cross(b - a, c - b)); +// } +// else { + auto abc = PointXYZ::dot(a, b) + PointXYZ::dot(b, c) + PointXYZ::dot(c, a); + auto a_bc = PointXYZ::dot(a, PointXYZ::cross(b, c)); + triangles[itri].area = 2. * std::atan(std::abs(a_bc) / (1. + abc)); +// } + ++itri; } triangles.size() = itri; return triangles; + } @@ -408,95 +265,93 @@ double ConvexSphericalPolygon::SubTriangles::area() const { return area; } -// @return lowest point id of this polygon's segment intersecting [s1,s2)) -int ConvexSphericalPolygon::intersect(const int start, const GreatCircleSegment& s, PointXYZ& I) const { - for (int i = start; i < size_; i++) { - const int id0 = i; - const int id1 = next(i); - const GreatCircleSegment p(sph_coords_[id0], sph_coords_[id1]); - I = s.intersect(p); - if (I[0] == 0 && I[1] == 0 && I[2] == 0) { - // intersection not on [p1,p2) - continue; - } - if (I[0] == 1 && I[1] == 1) { - return OVERLAP; - } - return id0; - } - return NO_INTERSECT; -} - - -void ConvexSphericalPolygon::clip(const GreatCircleSegment& great_circle) { +void ConvexSphericalPolygon::clip(const GreatCircleSegment& great_circle, std::ostream* out, double pointsSameEPS) { ATLAS_ASSERT(valid_); - ATLAS_ASSERT(distance2(great_circle.first(), great_circle.second()) > TOL2); - + ATLAS_ASSERT(not approx_eq(great_circle.first(), great_circle.second())); + std::vector clipped_sph_coords; + clipped_sph_coords.reserve(ConvexSphericalPolygon::MAX_SIZE); auto invalidate_this_polygon = [&]() { size_ = 0; valid_ = false; area_ = 0.; }; - - // Count and mark all vertices to be possibly considered in clipped polygon - StackVector vertex_in(size_); - int num_vertices_in = 0; +#if DEBUG_OUTPUT + int add_point_num = 0; +#endif + bool first_in = great_circle.inLeftHemisphere(sph_coords_[0], -1.5 * EPS, out); for (int i = 0; i < size_; i++) { - vertex_in[i] = great_circle.inLeftHemisphere(sph_coords_[i], -10. * TOL); - num_vertices_in += vertex_in[i] ? 1 : 0; - } - - PointXYZ i1; - const int f1 = intersect(0, great_circle, i1); - const bool segment_only_touches_last_point = (f1 == size_ - 1); - if (f1 == OVERLAP || f1 == NO_INTERSECT || segment_only_touches_last_point) { - if (num_vertices_in < 3) { - invalidate_this_polygon(); - } - return; - } - PolygonEdgeIntersection intersection_1(*this, f1, i1); - - PointXYZ i2; // second intersection point - auto start2 = [&]() { return f1 + 1 + (intersection_1.isPointAtEnd() ? 1 : 0); }; - const int f2 = intersect(start2(), great_circle, i2); - if (f2 == OVERLAP || f2 == NO_INTERSECT) { - if (num_vertices_in < 3) { - invalidate_this_polygon(); + int in = (i+1) % size_; + bool second_in = great_circle.inLeftHemisphere(sph_coords_[in], -1.5 * EPS, out); +#if DEBUG_OUTPUT + if (out) { + out->precision(18); + *out << " ** first: " << xyz2lonlat(sph_coords_[i]) << ", in ? " << first_in << std::endl; + *out << " second: " << xyz2lonlat(sph_coords_[in]) << ", in ? " << second_in << std::endl; } - return; - } - PolygonEdgeIntersection intersection_2(*this, f2, i2); - - // Create new vector of clipped coordinates - StackVector clipped_sph_coords; - { - auto keep_vertex = [&](int index) { clipped_sph_coords.push_back(sph_coords_[index]); }; - auto insert_point = [&](const PointXYZ& p) { clipped_sph_coords.push_back(p); }; - - for (int i = 0; i <= f1; i++) { - if (vertex_in[i]) { - keep_vertex(i); +#endif + if (first_in and second_in) { + clipped_sph_coords.emplace_back(sph_coords_[in]); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) both_in add second: " << xyz2lonlat(sph_coords_[in]) << std::endl; } +#endif } - if ((not vertex_in[f1] and intersection_1.isPointAtBegin()) or - (not vertex_in[next(f1)] and intersection_1.isPointAtEnd()) or intersection_1.isPointInside()) { - insert_point(i1); + else if (not first_in and not second_in) { + // continue to update first_in } - for (int i = f1 + 1; i <= f2; i++) { - if (vertex_in[i]) { - keep_vertex(i); + else { + const GreatCircleSegment segment(sph_coords_[i], sph_coords_[in]); + PointXYZ ip = great_circle.intersect(segment, out, pointsSameEPS); +#if DEBUG_OUTPUT + if (out) { + *out << " ip : " << xyz2lonlat(ip) << std::endl; } - } - if ((not vertex_in[f2] and intersection_2.isPointAtBegin()) or - (not vertex_in[next(f2)] and intersection_2.isPointAtEnd()) or intersection_2.isPointInside()) { - insert_point(i2); - } - for (int i = f2 + 1; i < size_; i++) { - if (vertex_in[i]) { - keep_vertex(i); +#endif + if (ip[0] == 1 and ip[1] == 1 and ip[2] == 1) { + // consider the segments parallel +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) ip=(1,1,1) add second: " + << xyz2lonlat(sph_coords_[in]) << std::endl; + } +#endif + clipped_sph_coords.emplace_back(sph_coords_[in]); + first_in = second_in; + continue; + } + if (second_in) { + int inn = (in+1) % size_; + const GreatCircleSegment segment_n(sph_coords_[in], sph_coords_[inn]); + if (segment.inLeftHemisphere(ip, -1.5 * EPS, out) and + segment_n.inLeftHemisphere(ip, -1.5 * EPS, out) and + (PointXYZ::distance(ip, sph_coords_[in]) > pointsSameEPS)) { + clipped_sph_coords.emplace_back(ip); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) second_in add ip: " << xyz2lonlat(ip) << std::endl; + } +#endif + } + clipped_sph_coords.emplace_back(sph_coords_[in]); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) second_in add second: " << xyz2lonlat(sph_coords_[in]) << std::endl; + } +#endif + } + else { + if (PointXYZ::distance(ip, sph_coords_[i]) > pointsSameEPS) { + clipped_sph_coords.emplace_back(ip); +#if DEBUG_OUTPUT + if (out) { + *out << " ((" << ++add_point_num << ")) first_in add ip: " << xyz2lonlat(ip) << std::endl; + } +#endif + } } } + first_in = second_in; } // Update polygon @@ -516,22 +371,68 @@ void ConvexSphericalPolygon::clip(const GreatCircleSegment& great_circle) { // intersect a polygon with this polygon // @param[in] pol clipping polygon // @param[out] intersecting polygon -ConvexSphericalPolygon ConvexSphericalPolygon::intersect(const ConvexSphericalPolygon& plg) const { - ConvexSphericalPolygon intersection = *this; +ConvexSphericalPolygon ConvexSphericalPolygon::intersect(const ConvexSphericalPolygon& plg, std::ostream* out, double pointsSameEPS) const { + + bool fpe_disabled = atlas::library::disable_floating_point_exception(FE_INVALID); + auto restore_fpe = [fpe_disabled] { + if (fpe_disabled) { + atlas::library::enable_floating_point_exception(FE_INVALID); + } + }; + + // the larger area polygon is the intersector + ConvexSphericalPolygon intersection; + ConvexSphericalPolygon intersector; + std::string intor_id = "P1"; + std::string inted_id = "P2"; + + bool this_intersector = (area() > plg.area()); + if (approx_eq(area(), plg.area())) { + PointXYZ dc = centroid() - plg.centroid(); + if (dc[0] > 0. or (dc[0] == 0. and (dc[1] > 0. or (dc[1] == 0. and dc[2] > 0.)))) { + this_intersector = true; + } + } + if (this_intersector) { + intersector = *this; + intersection = plg; + } + else { + intor_id = "P2"; + inted_id = "P1"; + intersector = plg; + intersection = *this; + } +#if DEBUG_OUTPUT + if (out) { + *out << inted_id << " : "; + print(*out); + *out << std::endl << intor_id << " : "; + plg.print(*out); + *out << std::endl; + } +#endif if (intersection.valid_) { - for (size_t i = 0; i < plg.size_; i++) { - const PointXYZ& s1 = plg.sph_coords_[i]; - const PointXYZ& s2 = plg.sph_coords_[(i != plg.size_ - 1) ? i + 1 : 0]; - intersection.clip(GreatCircleSegment(s1, s2)); + for (size_t i = 0; i < intersector.size_; i++) { + const PointXYZ& s1 = intersector.sph_coords_[i]; + const PointXYZ& s2 = intersector.sph_coords_[(i != intersector.size_ - 1) ? i + 1 : 0]; +#if DEBUG_OUTPUT + if (out) { + *out << std::endl << "Clip with [" << intor_id << "_" << i << ", " << intor_id << "_" + << (i+1) % intersector.size_ << "]" << std::endl; + } +#endif + intersection.clip(GreatCircleSegment(s1, s2), out, pointsSameEPS); if (not intersection.valid_) { + restore_fpe(); return intersection; } } } - intersection.simplify(); intersection.computed_area_ = false; intersection.computed_radius_ = false; intersection.computed_centroid_ = false; + restore_fpe(); return intersection; } @@ -548,72 +449,64 @@ void ConvexSphericalPolygon::print(std::ostream& out) const { out << "}"; } +std::string ConvexSphericalPolygon::json(int precision) const { + std::stringstream ss; + if( precision ) { + ss.precision(16); + } + ss << "["; + for (size_t i = 0; i < size(); ++i) { + if (i > 0) { + ss << ","; + } + PointLonLat ip_ll; + xyz2lonlat(sph_coords_[i], ip_ll); + ss << "[" << ip_ll.lon() << "," << ip_ll.lat() << "]"; + } + ss << "]"; + return ss.str(); +} + + double ConvexSphericalPolygon::compute_radius() const { double radius{0.}; if (valid_) { - PointXYZ centroid; - centroid = sph_coords_[0]; - size_t isp = 1; - for (size_t i = 1; i < size_; ++i) { - if (approx_eq(sph_coords_[isp], sph_coords_[isp - 1], TOL)) { - continue; - } - centroid = centroid + sph_coords_[isp]; - ++isp; + if (not computed_centroid_) { + compute_centroid(); } - centroid = PointXYZ::div(centroid, PointXYZ::norm(centroid)); - for (size_t i = 0; i < size_; ++i) { - radius = std::max(radius, PointXYZ::distance(sph_coords_[i], centroid)); + radius = std::max(radius, PointXYZ::distance(sph_coords_[i], centroid_)); } } return radius; } -bool ConvexSphericalPolygon::GreatCircleSegment::contains(const PointXYZ& p) const { - /* - * @brief Point-on-segment test on great circle segments - * @param[in] P given point in (x,y,z) coordinates - * @return - */ - constexpr double eps_large = 1.e3 * EPS; - - // Case where p is one of the endpoints - double pp1n2 = distance2(p, p1_); - double pp2n2 = distance2(p, p2_); - if (pp1n2 < EPS2 or pp2n2 < EPS2) { - return true; - } - - PointXYZ p12 = cross(); - double p12n2 = norm2(p12); - double p12n = std::sqrt(p12n2); - p12 /= p12n; - if (std::abs(PointXYZ::dot(p, p12)) > eps_large) { - return false; - } - double pp = PointXYZ::distance(p1_, p2_); - double pp1 = PointXYZ::distance(p, p1_); - double pp2 = PointXYZ::distance(p, p2_); - return (std::min(pp - pp1, pp - pp2) > -eps_large); -} - -PointXYZ ConvexSphericalPolygon::GreatCircleSegment::intersect(const GreatCircleSegment& p) const { +// 'this' great circle's intersection with the segment 'p': [p.first(), p.second()) +PointXYZ ConvexSphericalPolygon::GreatCircleSegment::intersect(const GreatCircleSegment& p, std::ostream* out, double pointsSameEPS) const { const auto& s = *this; PointXYZ sp = PointXYZ::cross(s.cross(), p.cross()); double sp_norm = PointXYZ::norm(sp); - if (sp_norm > EPS) { + bool gcircles_distinct = (sp_norm > EPS); +#if DEBUG_OUTPUT + if (out) { + *out << " Great circles distinct ? " << sp_norm << " > " << EPS << " ? " + << gcircles_distinct << std::endl; + } +#endif + if (gcircles_distinct) { sp /= sp_norm; - if (p.contains(sp)) { + auto sp_2 = sp * -1.; + double d = distance2(p.first(), p.second()); + double d1 = std::max(distance2(sp, p.first()), distance2(sp, p.second())); + double d2 = std::max(distance2(sp_2, p.first()), distance2(sp_2, p.second())); + if (d1 < d2) { return sp; } - sp *= -1.; - if (p.contains(sp)) { - return sp; + else { + return sp_2; } - return PointXYZ(0, 0, 0); } else { return PointXYZ(1, 1, 1); diff --git a/src/atlas/util/ConvexSphericalPolygon.h b/src/atlas/util/ConvexSphericalPolygon.h index 9a5783fa0..7b3338fd5 100644 --- a/src/atlas/util/ConvexSphericalPolygon.h +++ b/src/atlas/util/ConvexSphericalPolygon.h @@ -16,6 +16,9 @@ #include "atlas/util/Point.h" #include "atlas/util/Polygon.h" +#include "atlas/runtime/Log.h" + +#define DEBUG_OUTPUT 1 namespace atlas { namespace util { @@ -31,17 +34,20 @@ class ConvexSphericalPolygon { class GreatCircleSegment { public: GreatCircleSegment(const PointXYZ& p1, const PointXYZ& p2): p1_(p1), p2_(p2), cross_(PointXYZ::cross(p1, p2)) {} - bool contains(const PointXYZ&) const; - - PointXYZ intersect(const GreatCircleSegment&) const; - // Hemisphere is defined by segment when walking from first() to second() - // Positive offset: distance into left hemisphere, e.g. to exclude segment itself with tolerance - // Negative offset: distance into right hemisphere, e.g. to include segment with tolerance - bool inLeftHemisphere(const PointXYZ& P, const double offset = 0.) const { - return (PointXYZ::dot(cross(), P) > offset); + // For a given segment, the "left" hemisphere is defined on the left of the segment when walking from first() to second() + inline bool inLeftHemisphere(const PointXYZ& P, const double offset = 0., std::ostream* out = NULL) const { +#if DEBUG_OUTPUT + if (out) { + *out << " inLeftHemi: " <= " << offset << " ? " + << (PointXYZ::dot(cross(), P) >= offset) << std::endl; + } +#endif + return (PointXYZ::dot(cross(), P) >= offset); // has to have = included } + PointXYZ intersect(const GreatCircleSegment&, std::ostream* f = NULL, double pointsSameEPS = std::numeric_limits::epsilon()) const; + const PointXYZ& first() const { return p1_; } const PointXYZ& second() const { return p2_; } @@ -96,7 +102,7 @@ class ConvexSphericalPolygon { return radius_; } - ConvexSphericalPolygon intersect(const ConvexSphericalPolygon& pol) const; + ConvexSphericalPolygon intersect(const ConvexSphericalPolygon& pol, std::ostream* f = nullptr, double pointsEqualEPS = std::numeric_limits::epsilon()) const; /* * @brief check if two spherical polygons area equal @@ -107,6 +113,8 @@ class ConvexSphericalPolygon { void print(std::ostream&) const; + std::string json(int precision = 0) const; + friend std::ostream& operator<<(std::ostream& out, const ConvexSphericalPolygon& p) { p.print(out); return out; @@ -136,29 +144,15 @@ class ConvexSphericalPolygon { double compute_radius() const; - SubTriangles triangulate(const double radius) const; + SubTriangles triangulate() const; - void clip(const GreatCircleSegment&); + void clip(const GreatCircleSegment&, std::ostream* f = nullptr, double pointsSameEPS = std::numeric_limits::epsilon()); /* * @return true:polygon is convex */ bool validate(); - /* - * @brief Segment-sph_polygon intersection - * @param[in] s1, s2 segment endpoints in (x,y,z) coordinates - * @param[in] start start with polygon segments [pol[start],pol[start+1]],... - * @param[out] ip intersection point or nullptr - * @return -1: overlap with one of polygon edges, - * 0: no_intersect, - * 1 + (id of this polygon's segment intersecting [s1,s2)): otherwise - */ - int intersect(const int start, const GreatCircleSegment&, PointXYZ& ip) const; - - /// Makes the polygon convex and skips consequtive node that is too close to previous - void simplify(); - private: std::array sph_coords_; mutable PointXYZ centroid_; diff --git a/src/sandbox/interpolation/atlas-conservative-interpolation.cc b/src/sandbox/interpolation/atlas-conservative-interpolation.cc index 48b6e71eb..6e5199674 100644 --- a/src/sandbox/interpolation/atlas-conservative-interpolation.cc +++ b/src/sandbox/interpolation/atlas-conservative-interpolation.cc @@ -8,19 +8,6 @@ * nor does it submit to any jurisdiction. */ - -// TODO: -// ----- -// * Fix abort encountered with -// mpirun -np 4 atlas-conservative-interpolation --source.grid=O20 --target.grid=H8 --order=2 -// -// QUESTIONS: -// ---------- -// * Why sqrt in ConservativeSphericalPolygon in line -// remap_stat.errors[Statistics::Errors::REMAP_CONS] = std::sqrt(std::abs(err_remap_cons) / unit_sphere_area()); -// used to compute conservation_error - - #include #include #include @@ -35,6 +22,8 @@ #include "atlas/array/MakeView.h" #include "atlas/field.h" #include "atlas/grid.h" +#include "atlas/grid/Spacing.h" +#include "atlas/grid/detail/spacing/CustomSpacing.h" #include "atlas/interpolation/Interpolation.h" #include "atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h" #include "atlas/mesh.h" @@ -73,13 +62,15 @@ class AtlasParallelInterpolation : public AtlasTool { add_option(new eckit::option::Separator("Source/Target options")); add_option(new SimpleOption("source.grid", "source gridname")); + add_option(new SimpleOption("source.partitioner", "source partitioner name (spherical-polygon, lonlat-polygon, brute-force)")); add_option(new SimpleOption("target.grid", "target gridname")); + add_option(new SimpleOption("target.partitioner", "target partitioner name (equal_regions, regular_bands, equal_bands)")); add_option(new SimpleOption("source.functionspace", "source functionspace, to override source grid default")); add_option(new SimpleOption("target.functionspace", "target functionspace, to override target grid default")); add_option(new SimpleOption("source.halo", "default=2")); - add_option(new SimpleOption("target.halo", "default=0")); + add_option(new SimpleOption("target.halo", "default=0 for CellColumns and 1 for NodeColumns")); add_option(new eckit::option::Separator("Interpolation options")); add_option(new SimpleOption("order", "Interpolation order. Supported: 1, 2 (default=1)")); @@ -182,8 +173,22 @@ std::function get_init(const AtlasTool::Args& args) } int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { - auto src_grid = Grid{args.getString("source.grid", "H16")}; - auto tgt_grid = Grid{args.getString("target.grid", "H32")}; + auto get_grid = [](std::string grid_name) { + int grid_number = std::atoi( grid_name.substr(1, grid_name.size()).c_str() ); + if (grid_name.at(0) == 'P') { + Log::info() << "P-grid number: " << grid_number << std::endl; + ATLAS_ASSERT(grid_number > 3); + std::vector y = {90, 89.9999, 0, -90}; + auto xspace = StructuredGrid::XSpace( grid::LinearSpacing(0, 360, grid_number, false) ); + auto yspace = StructuredGrid::YSpace( new grid::spacing::CustomSpacing( y.size(), y.data() ) ); + return StructuredGrid(xspace, yspace); + } + else { + return StructuredGrid{grid_name}; + } + }; + auto src_grid = get_grid(args.getString("source.grid", "H32")); + auto tgt_grid = get_grid(args.getString("target.grid", "H32")); auto create_functionspace = [&](Mesh& mesh, int halo, std::string type) -> FunctionSpace { if (type.empty()) { @@ -199,13 +204,13 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { return functionspace::CellColumns(mesh, option::halo(halo)); } else if (type == "NodeColumns") { - return functionspace::NodeColumns(mesh, option::halo(halo)); + return functionspace::NodeColumns(mesh, option::halo(std::max(1,halo))); } ATLAS_THROW_EXCEPTION("FunctionSpace " << type << " is not recognized."); }; timers.target_setup.start(); - auto tgt_mesh = Mesh{tgt_grid}; + auto tgt_mesh = Mesh{tgt_grid, grid::Partitioner(args.getString("target.partitioner", "regular_bands"))}; auto tgt_functionspace = create_functionspace(tgt_mesh, args.getLong("target.halo", 0), args.getString("target.functionspace", "")); auto tgt_field = tgt_functionspace.createField(); @@ -213,8 +218,8 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { timers.source_setup.start(); auto src_meshgenerator = - MeshGenerator{src_grid.meshgenerator() | option::halo(2) | util::Config("pole_elements", "pentagons")}; - auto src_partitioner = grid::MatchingPartitioner{tgt_mesh}; + MeshGenerator{src_grid.meshgenerator() | option::halo(2) | util::Config("pole_elements", "")}; + auto src_partitioner = grid::MatchingPartitioner{tgt_mesh, util::Config("partitioner",args.getString("source.partitioner", "spherical-polygon"))}; auto src_mesh = src_meshgenerator.generate(src_grid, src_partitioner); auto src_functionspace = create_functionspace(src_mesh, args.getLong("source.halo", 2), args.getString("source.functionspace", "")); @@ -230,14 +235,14 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { for (idx_t n = 0; n < lonlat.shape(0); ++n) { src_view(n) = f(PointLonLat{lonlat(n, LON), lonlat(n, LAT)}); } - src_field.set_dirty(false); + src_field.set_dirty(true); timers.initial_condition.start(); } - timers.interpolation_setup.start(); auto interpolation = Interpolation(option::type("conservative-spherical-polygon") | args, src_functionspace, tgt_functionspace); + Log::info() << interpolation << std::endl; timers.interpolation_setup.stop(); @@ -251,11 +256,13 @@ int AtlasParallelInterpolation::execute(const AtlasTool::Args& args) { using Statistics = interpolation::method::ConservativeSphericalPolygonInterpolation::Statistics; Statistics stats(metadata); if (args.getBool("statistics.accuracy", false)) { - stats.accuracy(interpolation, tgt_field, get_init(args)); + metadata.set( stats.accuracy(interpolation, tgt_field, get_init(args) ) ); } if (args.getBool("statistics.conservation", false)) { // compute difference field src_conservation_field = stats.diff(interpolation, src_field, tgt_field); + src_conservation_field.set_dirty(true); + src_conservation_field.haloExchange(); } } diff --git a/src/tests/interpolation/test_interpolation_conservative.cc b/src/tests/interpolation/test_interpolation_conservative.cc index 92b426ff0..5cfde0448 100644 --- a/src/tests/interpolation/test_interpolation_conservative.cc +++ b/src/tests/interpolation/test_interpolation_conservative.cc @@ -25,6 +25,7 @@ #include "atlas/meshgenerator.h" #include "atlas/option.h" #include "atlas/util/Config.h" +#include "atlas/util/function/VortexRollup.h" #include "tests/AtlasTestEnvironment.h" @@ -36,17 +37,26 @@ using ConservativeMethod = interpolation::method::ConservativeSphericalPolygonIn using Statistics = ConservativeMethod::Statistics; void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function func, - Statistics& remap_stat_1, Statistics& remap_stat_2) { + Statistics& remap_stat_1, Statistics& remap_stat_2, bool src_cell_data, bool tgt_cell_data) { + std::string src_data_type = (src_cell_data ? "CellColumns(" : "NodeColumns("); + std::string tgt_data_type = (tgt_cell_data ? "CellColumns(" : "NodeColumns("); + Log::info() << "+-----------------------\n"; + Log::info() << src_data_type << src_grid.name() << ") --> " << tgt_data_type << tgt_grid.name() <<")\n"; + Log::info() << "+-----------------------\n"; Log::info().indent(); + // setup conservative remap: compute weights, polygon intersection, etc util::Config config("type", "conservative-spherical-polygon"); config.set("order", 1); config.set("validate", true); config.set("statistics.intersection", true); config.set("statistics.conservation", true); + config.set("src_cell_data", src_cell_data); + config.set("tgt_cell_data", tgt_cell_data); auto conservative_interpolation = Interpolation(config, src_grid, tgt_grid); Log::info() << conservative_interpolation << std::endl; + Log::info() << std::endl; // create source field from analytic function "func" const auto& src_fs = conservative_interpolation.source(); @@ -85,6 +95,7 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function 1st order constructing new matrix"); @@ -95,14 +106,16 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function 1st order matrix-free"); cfg.set("matrix_free", true); cfg.set("order", 1); auto interpolation = Interpolation(cfg, src_grid, tgt_grid, cache); Log::info() << interpolation << std::endl; interpolation.execute(src_field, tgt_field); + Log::info() << std::endl; } auto cache_2 = interpolation::Cache{}; { @@ -113,14 +126,16 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function 2nd order matrix-free"); cfg.set("matrix_free", true); cfg.set("order", 2); auto interpolation = Interpolation(cfg, src_grid, tgt_grid, cache); Log::info() << interpolation << std::endl; interpolation.execute(src_field, tgt_field); + Log::info() << std::endl; } { ATLAS_TRACE("cached -> 2nd order using cached matrix"); @@ -129,10 +144,10 @@ void do_remapping_test(Grid src_grid, Grid tgt_grid, std::function tol) { - auto improvement = [](double& e, double& r) { return 100. * (r - e) / r; }; + //auto improvement = [](double& e, double& r) { return 100. * (r - e) / r; }; + auto improvement = [](double& e, double& r) { return r - e; }; double err; // check polygon intersections - err = remap_stat_1.errors[Statistics::Errors::GEO_DIFF]; - Log::info() << "Polygon area computation improvement: " << improvement(err, tol[0]) << " %" << std::endl; + err = remap_stat_1.errors[Statistics::Errors::SRCTGT_INTERSECTPLG_DIFF]; + Log::info() << "Polygon area computation (new < ref) = (" << err << " < " << tol[0] << ")" << std::endl; EXPECT(err < tol[0]); - err = remap_stat_1.errors[Statistics::Errors::GEO_L1]; - Log::info() << "Polygon intersection improvement : " << improvement(err, tol[1]) << " %" << std::endl; + err = remap_stat_1.errors[Statistics::Errors::TGT_INTERSECTPLG_L1]; + Log::info() << "Polygon intersection (new < ref) = (" << err << " < " << tol[1] << ")" << std::endl; EXPECT(err < tol[1]); // check remap accuracy - err = remap_stat_1.errors[Statistics::Errors::REMAP_L2]; - Log::info() << "1st order accuracy improvement : " << improvement(err, tol[2]) << " %" << std::endl; + err = std::abs(remap_stat_1.errors[Statistics::Errors::REMAP_L2]); + Log::info() << "1st order accuracy (new < ref) = (" << err << " < " << tol[2] << ")" << std::endl; EXPECT(err < tol[2]); - err = remap_stat_2.errors[Statistics::Errors::REMAP_L2]; - Log::info() << "2nd order accuracy improvement : " << improvement(err, tol[3]) << " %" << std::endl; + err = std::abs(remap_stat_2.errors[Statistics::Errors::REMAP_L2]); + Log::info() << "2nd order accuracy (new < ref) = (" << err << " < " << tol[3] << ")" << std::endl; EXPECT(err < tol[3]); // check mass conservation - err = remap_stat_1.errors[Statistics::Errors::REMAP_CONS]; - Log::info() << "1st order conservation improvement : " << improvement(err, tol[4]) << " %" << std::endl; + err = std::abs(remap_stat_1.errors[Statistics::Errors::REMAP_CONS]); + Log::info() << "1st order conservation (new < ref) = (" << err << " < " << tol[4] << ")" << std::endl; EXPECT(err < tol[4]); - err = remap_stat_2.errors[Statistics::Errors::REMAP_CONS]; - Log::info() << "2nd order conservation improvement : " << improvement(err, tol[5]) << " %" << std::endl - << std::endl; + err = std::abs(remap_stat_2.errors[Statistics::Errors::REMAP_CONS]); + Log::info() << "2nd order conservation (new < ref) = (" << err << " < " << tol[5] << ")" << std::endl; EXPECT(err < tol[5]); Log::info().unindent(); } CASE("test_interpolation_conservative") { -#if 1 SECTION("analytic constfunc") { auto func = [](const PointLonLat& p) { return 1.; }; Statistics remap_stat_1; Statistics remap_stat_2; - do_remapping_test(Grid("H47"), Grid("H48"), func, remap_stat_1, remap_stat_2); - check(remap_stat_1, remap_stat_2, {1.e-13, 5.e-8, 2.9e-6, 2.9e-6, 5.5e-5, 5.5e-5}); + bool src_cell_data = true; + bool tgt_cell_data = true; + do_remapping_test(Grid("O128"), Grid("H48"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-13, 1.0e-13, 1.0e-13, 1.0e-13, 1.0e-13}); } - SECTION("analytic Y_2^2 as in Jones(1998)") { + SECTION("vortex_rollup") { auto func = [](const PointLonLat& p) { - double cos = std::cos(0.025 * p[0]); - return 2. + cos * cos * std::cos(2 * 0.025 * p[1]); + return util::function::vortex_rollup(p[0], p[1], 0.5); }; Statistics remap_stat_1; Statistics remap_stat_2; - do_remapping_test(Grid("H47"), Grid("H48"), func, remap_stat_1, remap_stat_2); - check(remap_stat_1, remap_stat_2, {1.e-13, 5.e-8, 4.8e-4, 1.1e-4, 8.9e-5, 1.1e-4}); + + bool src_cell_data = true; + bool tgt_cell_data = true; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); + + src_cell_data = true; + tgt_cell_data = false; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); + + src_cell_data = false; + tgt_cell_data = true; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-13, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); + + src_cell_data = false; + tgt_cell_data = false; + do_remapping_test(Grid("O32"), Grid("H24"), func, remap_stat_1, remap_stat_2, src_cell_data, tgt_cell_data); + check(remap_stat_1, remap_stat_2, {1.0e-12, 1.0e-12, 3.0e-3, 6.0e-4, 1.0e-15, 1.0e-8}); } -#endif } } // namespace test diff --git a/src/tests/util/test_convexsphericalpolygon.cc b/src/tests/util/test_convexsphericalpolygon.cc index f4e158855..f7897a29e 100644 --- a/src/tests/util/test_convexsphericalpolygon.cc +++ b/src/tests/util/test_convexsphericalpolygon.cc @@ -2,6 +2,8 @@ #include #include #include +#include +#include #include "atlas/util/ConvexSphericalPolygon.h" #include "atlas/util/Geometry.h" @@ -13,6 +15,7 @@ namespace atlas { namespace test { using ConvexSphericalPolygon = util::ConvexSphericalPolygon; +const double EPS = std::numeric_limits::epsilon(); util::ConvexSphericalPolygon make_polygon(const std::initializer_list& list) { return util::ConvexSphericalPolygon{std::vector(list)}; @@ -29,17 +32,112 @@ util::ConvexSphericalPolygon make_polygon() { return util::ConvexSphericalPolygon{}; } + +template +std::string to_json(const It& begin, const It& end, int precision = 0) { + std::stringstream ss; + ss << "[\n"; + size_t size = std::distance(begin,end); + size_t c=0; + for( auto it = begin; it != end; ++it, ++c ) { + ss << " " << it->json(precision); + if( c < size-1 ) { + ss << ",\n"; + } + } + ss << "\n]"; + return ss.str(); +} + +template +std::string to_json(const ConvexSphericalPolygonContainer& polygons, int precision = 0) { + return to_json(polygons.begin(),polygons.end(),precision); +} +std::string to_json(std::initializer_list&& polygons, int precision = 0) { + return to_json(polygons.begin(),polygons.end(),precision); +} + +void check_intersection(const ConvexSphericalPolygon& plg1, const ConvexSphericalPolygon& plg2, const ConvexSphericalPolygon& iplg_sol, double pointsSameEPS = 5.e6 * EPS, std::ostream* out = nullptr) { + auto iplg = plg1.intersect(plg2, out, pointsSameEPS); + Log::info().indent(); + Log::info() << "plg1 area : " << plg1.area() << "\n"; + Log::info() << "plg2 area : " << plg2.area() << "\n"; + Log::info() << "iplg area : " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg1,plg2,iplg},20) << "\n"; + EXPECT(std::min(plg1.area(), plg2.area()) >= iplg.area()); + EXPECT(iplg.equals(iplg_sol, 1.e-8)); + Log::info().unindent(); +} + CASE("test default constructor") { ConvexSphericalPolygon p; EXPECT(bool(p) == false); } +CASE("Size of ConvexSphericalPolygon") { + // This test illustrates that ConvexSphericalPolygon is allocated on the stack completely, + // as sizeof(ConvexSphericalPolygon) includes space for MAX_SIZE coordinates of type PointXYZ + EXPECT(sizeof(PointXYZ) == sizeof(double) * 3); + size_t expected_size = 0; + expected_size += (1 + ConvexSphericalPolygon::MAX_SIZE) * sizeof(PointXYZ); + expected_size += sizeof(size_t); + expected_size += sizeof(bool); + expected_size += 2 * sizeof(double); + EXPECT(sizeof(ConvexSphericalPolygon) >= expected_size); // greater because compiler may add some padding +} + +CASE("analyse intersect") { + const double du = 0.5; + const double dv = 1.1 * EPS; + const double duc = 0.5 * du; + const double sduc = std::sqrt(1. - 0.25 * du * du); + const double dvc = 1. - 0.5 * dv * dv; + const double sdvc = dv * std::sqrt(1. - 0.25 * dv * dv); + PointXYZ s0p0{sduc, -duc, 0.}; + PointXYZ s0p1{sduc, duc, 0.}; + PointXYZ s1p0{dvc * sduc, -dvc * duc, -sdvc}; + PointXYZ s1p1{dvc * sduc, dvc * duc, sdvc}; + + EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p0 - s1p0), EPS); + EXPECT_APPROX_EQ(du, PointXYZ::norm(s0p0 - s0p1), EPS); + EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p1 - s1p1), EPS); + + ConvexSphericalPolygon::GreatCircleSegment s1(s0p0, s0p1); + ConvexSphericalPolygon::GreatCircleSegment s2(s1p0, s1p1); + + // analytical solution + PointXYZ Isol{1., 0., 0.}; + + // test "intersection" + PointXYZ I12 = s1.intersect(s2); + PointXYZ I21 = s2.intersect(s1); + EXPECT_APPROX_EQ(std::abs(PointXYZ::norm(I12) - 1.), 0., EPS); + EXPECT_APPROX_EQ(PointXYZ::norm(I12 - Isol), 0., EPS); + EXPECT_APPROX_EQ(std::abs(PointXYZ::norm(I21) - 1.), 0., EPS); + EXPECT_APPROX_EQ(PointXYZ::norm(I21 - Isol), 0., EPS); +} + +CASE("test_json_format") { + auto plg = make_polygon({{0., 45.}, {0., 0.}, {90., 0.}, {90., 45.}}); + EXPECT_EQ(plg.json(), "[[0,45],[0,0],[90,0],[90,45]]"); + + std::vector plg_g = { + make_polygon({{0, 60}, {0, 50}, {40, 60}}), //0 + make_polygon({{0, 60}, {0, 50}, {20, 60}}), + make_polygon({{10, 60}, {10, 50}, {30, 60}}), //3 + }; + Log::info() << to_json(plg_g) << std::endl; + +} + CASE("test_spherical_polygon_area") { auto plg1 = make_polygon({{0., 90.}, {0., 0.}, {90., 0.}}); EXPECT_APPROX_EQ(plg1.area(), M_PI_2); auto plg2 = make_polygon({{0., 45.}, {0., 0.}, {90., 0.}, {90., 45.}}); auto plg3 = make_polygon({{0., 90.}, {0., 45.}, {90., 45.}}); - Log::info() << "area diff: " << plg1.area() - plg2.area() - plg3.area() << std::endl; + Log::info().indent(); + EXPECT(plg1.area() - plg2.area() - plg3.area() <= EPS); + Log::info().unindent(); EXPECT_APPROX_EQ(std::abs(plg1.area() - plg2.area() - plg3.area()), 0, 1e-15); } @@ -50,23 +148,23 @@ CASE("test_spherical_polygon_intersection") { std::array plg_f = {make_polygon({{0, 70}, {0, 60}, {40, 60}, {40, 70}}), make_polygon({{0, 90}, {0, 0}, {40, 0}})}; std::array plg_g = { - make_polygon({{0, 60}, {0, 50}, {40, 60}}), //0 - make_polygon({{0, 60}, {0, 50}, {20, 60}}), - make_polygon({{10, 60}, {10, 50}, {30, 60}}), //3 - make_polygon({{40, 80}, {0, 60}, {40, 60}}), - make_polygon({{0, 80}, {0, 60}, {40, 60}}), //5 - make_polygon({{20, 80}, {0, 60}, {40, 60}}), - make_polygon({{20, 70}, {0, 50}, {40, 50}}), //7 - make_polygon({{0, 90}, {0, 60}, {40, 60}}), - make_polygon({{-10, 80}, {-10, 50}, {50, 80}}), //9 - make_polygon({{0, 80}, {0, 50}, {40, 50}, {40, 80}}), - make_polygon({{0, 65}, {20, 55}, {40, 60}, {20, 65}}), //11 - make_polygon({{20, 65}, {0, 60}, {20, 55}, {40, 60}}), - make_polygon({{10, 63}, {20, 55}, {30, 63}, {20, 65}}), //13 - make_polygon({{20, 75}, {0, 70}, {5, 5}, {10, 0}, {20, 0}, {40, 70}}), - make_polygon({{0, 50}, {0, 40}, {5, 45}}), //15 - make_polygon({{0, 90}, {0, 80}, {20, 0}, {40, 80}}), - make_polygon({{0, 65}, {0, 55}, {40, 65}, {40, 75}}), //17 + make_polygon({{0, 60}, {0, 50}, {40, 60}}), // 0 + make_polygon({{0, 60}, {0, 50}, {20, 60}}), // 1 + make_polygon({{10, 60}, {10, 50}, {30, 60}}), // 2 + make_polygon({{40, 80}, {0, 60}, {40, 60}}), // 3 + make_polygon({{0, 80}, {0, 60}, {40, 60}}), // 4 + make_polygon({{20, 80}, {0, 60}, {40, 60}}), // 5 + make_polygon({{20, 70}, {0, 50}, {40, 50}}), // 6 + make_polygon({{0, 90}, {0, 60}, {40, 60}}), // 7 + make_polygon({{-10, 80}, {-10, 50}, {50, 80}}), // 8 + make_polygon({{0, 80}, {0, 50}, {40, 50}, {40, 80}}), // 9 + make_polygon({{0, 65}, {20, 55}, {40, 60}, {20, 65}}), // 10 + make_polygon({{20, 65}, {0, 60}, {20, 55}, {40, 60}}), // 11 + make_polygon({{10, 63}, {20, 55}, {30, 63}, {20, 65}}), // 12 + make_polygon({{20, 75}, {0, 70}, {5, 5}, {10, 0}, {20, 0}, {40, 70}}), // 13 + make_polygon({{0, 50}, {0, 40}, {5, 45}}), // 14 + make_polygon({{0, 90}, {0, 80}, {20, 0}, {40, 80}}), // 15 + make_polygon({{0, 65}, {0, 55}, {40, 65}, {40, 75}}), // 16 }; std::array plg_i = { make_polygon(), //0 @@ -103,154 +201,208 @@ CASE("test_spherical_polygon_intersection") { make_polygon({{0, 50}, {0, 40}, {5, 45}}), make_polygon({{0, 90}, {0, 80}, {20, 0}, {40, 80}}), //32 make_polygon({{0, 65}, {0, 55}, {40, 65}, {40, 75}})}; + Log::info().indent(); for (int i = 0; i < nplg_f; i++) { for (int j = 0; j < nplg_g; j++) { - Log::info() << "\n(" << i * nplg_g + j << ") Intersecting polygon\n " << plg_f[i] << std::endl; - Log::info() << "with polygon\n " << plg_g[j] << std::endl; - auto plg_fg = plg_f[i].intersect(plg_g[j]); - auto plg_gf = plg_g[j].intersect(plg_f[i]); - Log::info() << "got polygon\n "; - if (plg_fg) { - Log::info() << plg_fg << std::endl; - Log::info() << " " << plg_gf << std::endl; - EXPECT(plg_fg.equals(plg_gf)); - EXPECT(plg_fg.equals(plg_i[i * nplg_g + j], 0.1)); - Log::info() << "instead of polygon\n "; - Log::info() << plg_i[i * nplg_g + j] << std::endl; - } - else { - Log::info() << " empty" << std::endl; - Log::info() << "instead of polygon\n "; - Log::info() << plg_i[i * nplg_g + j] << std::endl; + auto plg_fg = plg_f[i].intersect(plg_g[j], nullptr, 5.e6 * std::numeric_limits::epsilon()); + bool polygons_equal = plg_i[i * nplg_g + j].equals(plg_fg, 0.1); + EXPECT(polygons_equal); + if( not polygons_equal or (i==0 && j==10)) { + Log::info() << "\nIntersected the polygon plg_f["<= expected_size); // greater because compiler may add some padding -} +CASE("test_spherical_polygon_intersection_stretched") { + // plg1 is an octant + const auto plg1 = make_polygon({{0, 90}, {0, 0}, {90, 0}}); + Log::info().precision(20); -CASE("analyse intersect") { - const double EPS = std::numeric_limits::epsilon(); - const double du = 0.5; - const double dv = 1.1 * EPS; - const double duc = 0.5 * du; - const double sduc = std::sqrt(1. - 0.25 * du * du); - const double dvc = 1. - 0.5 * dv * dv; - const double sdvc = dv * std::sqrt(1. - 0.25 * dv * dv); - PointXYZ s0p0{sduc, -duc, 0.}; - PointXYZ s0p1{sduc, duc, 0.}; - PointXYZ s1p0{dvc * sduc, -dvc * duc, -sdvc}; - PointXYZ s1p1{dvc * sduc, dvc * duc, sdvc}; + const double delta = 1.1 * EPS; + const double u = 1. - delta; + const double v = std::sqrt(1 - u * u); + const double w = delta; + const double a = sqrt(1 - w * w); - EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p0 - s1p0), EPS); - EXPECT_APPROX_EQ(du, PointXYZ::norm(s0p0 - s0p1), EPS); - EXPECT_APPROX_EQ(dv, PointXYZ::norm(s0p1 - s1p1), EPS); + // plg2 is a stretched polygon + std::vector plg2_points = { + PointXYZ{u, v, 0.}, + PointXYZ{a, 0., w }, + PointXYZ{u, -v, 0.}, + PointXYZ{0., 0., -1.}}; + auto plg2 = util::ConvexSphericalPolygon(plg2_points.data(), plg2_points.size()); + EXPECT(plg2.size() == 4); - ConvexSphericalPolygon::GreatCircleSegment s1(s0p0, s0p1); - ConvexSphericalPolygon::GreatCircleSegment s2(s1p0, s1p1); + auto plg11 = plg1.intersect(plg1); + auto plg12 = plg1.intersect(plg2); + auto plg22 = plg2.intersect(plg2); - // analytical solution - PointXYZ Isol{1., 0., 0.}; + EXPECT(plg11.size() == 3); + EXPECT(plg12.size() == 3); + EXPECT(plg22.size() == 4); - // test "intersection" - PointXYZ I = s1.intersect(s2); - EXPECT_APPROX_EQ(std::abs(PointXYZ::norm(I) - 1.), 0., EPS); - EXPECT_APPROX_EQ(PointXYZ::norm(I - Isol), 0., EPS); + // plg3 is the analytical solution + std::vector plg3_points = { + PointXYZ{u, v, 0.}, + PointXYZ{a, 0, w}, + PointXYZ{1, 0, 0.}}; + auto plg3 = util::ConvexSphericalPolygon(plg3_points.data(), plg3_points.size()); + EXPECT(plg3.size() == 3); + EXPECT(plg12.size() == 3); + auto plg_sol = make_polygon({{1.2074e-06, 0.}, {0., 1.3994e-14}, {0., 0.}}); + EXPECT(plg12.equals(plg_sol, 1e-10)); - // test "contains" - EXPECT(s1.contains(Isol) && s2.contains(Isol)); - EXPECT(s1.contains(I) && s2.contains(I)); + Log::info().indent(); + Log::info() << "delta : " << delta << std::endl; + Log::info() << "plg12.area : " << plg12.area() << std::endl; + Log::info() << "exact intersection area : " << plg3.area() << std::endl; + double error_area = plg12.area() - plg3.area(); + EXPECT(error_area < EPS and error_area >= 0.); + Log::info().unindent(); + Log::info().precision(-1); } -CASE("source_covered") { - const double dd = 0.; - const auto csp0 = make_polygon({{0, 90}, {0, 0}, {90, 0}}); - double dcov1 = csp0.area(); // optimal coverage - double dcov2 = dcov1; // intersection-based coverage - double dcov3 = dcov1; // normalized intersection-based coverage - double darea = 0; // commutative area error in intersection: |area(A^B)-area(B^A)| - - double max_tarea = 0.; - double min_tarea = 0.; - double accumulated_tarea = 0.; - - const int n = 900; - const int m = 900; - const double dlat = 90. / n; +CASE("test_lonlat_pole_problem") { + const auto north_octant = make_polygon({{0, 90}, {0, 0}, {90, 0}}); + const double first_lat = 90. - 1.e+12 * EPS; + const int m = 10000; const double dlon = 90. / m; + std::vector csp(m); + Log::info().indent(); + + ATLAS_TRACE_SCOPE("create polygons") + for (int j = 0; j < m; j++) { + csp[j] = + make_polygon(PointLonLat{0, 90}, PointLonLat{dlon * j, first_lat}, PointLonLat{dlon * (j + 1), first_lat}); + } + + double max_area_overshoot = 0.; + int false_zero = 0; + for (int j = 0; j < m; j++) { + auto csp_i = csp[j].intersect(north_octant); + // intersection area should not be larger than its father polygon's + max_area_overshoot = std::max(max_area_overshoot, csp_i.area() - csp[j].area()); + if (csp_i.area() < EPS) { + false_zero++; + } + } + Log::info() << "False zero area intersection: " << false_zero << std::endl; + Log::info() << "Max area overshoot: " << max_area_overshoot << std::endl; + EXPECT(max_area_overshoot <= m * EPS); + EXPECT(false_zero == 0); + Log::info().unindent(); +} + +CASE("test_thin_elements_area") { + const auto north_octant = make_polygon({{0, 90}, {0, 0}, {90, 0}}); + const auto south_octant = make_polygon({{0,0}, {0, -90},{90, 0}}); + const int n = 3; + const int m = 2500; + const double dlat = 180. / n; + const double dlon = 90. / m; + + ATLAS_ASSERT(n > 1); std::vector csp(n * m); - std::vector tgt_area(n * m); - ATLAS_TRACE_SCOPE("1") + + ATLAS_TRACE_SCOPE("create polygons") for (int j = 0; j < m; j++) { csp[j] = make_polygon(PointLonLat{0, 90}, PointLonLat{dlon * j, 90 - dlat}, PointLonLat{dlon * (j + 1), 90 - dlat}); } - ATLAS_TRACE_SCOPE("2") - for (int i = 1; i < n; i++) { + for (int i = 1; i < n - 1; i++) { for (int j = 0; j < m; j++) { csp[i * m + j] = make_polygon( PointLonLat{dlon * j, 90 - dlat * i}, PointLonLat{dlon * j, 90 - dlat * (i + 1)}, PointLonLat{dlon * (j + 1), 90 - dlat * (i + 1)}, PointLonLat{dlon * (j + 1), 90 - dlat * i}); } } - ConvexSphericalPolygon cspi0; - ConvexSphericalPolygon csp0i; + for (int j = 0; j < m; j++) { + csp[(n - 1) * m + j] = + make_polygon(PointLonLat{dlon * j, 90 - dlat * (n-1)}, PointLonLat{dlon * j, -90.}, + PointLonLat{dlon * (j + 1), 90 - dlat * (n-1)}); + } + + double coverage_north = 0.; // north octant coverage by intersections with "csp"s + double coverage_south = 0.; // south octant coverage by intersections with "csp"s + double coverage_norm = 0.; // area sum of intersections in the north octant normalised to sum up to the area of the north octant + double coverage_csp = 0.; // area sum of all "csp"s + double accumulated_tarea = 0.; + std::vector i_north_area(n * m); - ATLAS_TRACE_SCOPE("3") + ATLAS_TRACE_SCOPE("intersect polygons") for (int i = 0; i < n; i++) { for (int j = 0; j < m; j++) { auto ipoly = i * m + j; - cspi0 = csp[ipoly].intersect(csp0); // intersect small csp[...] with large polygon csp0 - // should be approx csp[...] as well - csp0i = csp0.intersect(csp[ipoly]); // opposite: intersect csp0 with small csp[...] - // should be approx csp[...] as well - double a_csp = csp[ipoly].area(); - double a_cspi0 = cspi0.area(); // should be approx csp[...].area() - double a_csp0i = csp0i.area(); // should approx match a_cspi0 - EXPECT_APPROX_EQ(a_cspi0, a_csp, 1.e-10); - EXPECT_APPROX_EQ(a_csp0i, a_csp, 1.e-10); - darea = std::max(darea, a_cspi0 - a_csp0i); // should remain approx zero - dcov1 -= csp[ipoly].area(); - dcov2 -= a_cspi0; - tgt_area[ipoly] = a_cspi0; - accumulated_tarea += tgt_area[ipoly]; + double a_north_csp_i = csp[ipoly].intersect(north_octant).area(); // intersect narrow csp with the north octant + double a_south_csp_i = csp[ipoly].intersect(south_octant).area(); // intersect narrow csp with the south octant + if (i == 0) { + if (n == 2) { + EXPECT_APPROX_EQ(a_north_csp_i, csp[ipoly].area(), EPS); + } + else { + if (a_north_csp_i > csp[ipoly].area()) { + Log::info() << " error: " << a_north_csp_i - csp[ipoly].area() << std::endl; + } + } + } + coverage_north += a_north_csp_i; + coverage_csp += csp[ipoly].area(); + coverage_south += a_south_csp_i; + i_north_area[ipoly] = a_north_csp_i; } } - // normalize weights - double norm_fac = csp0.area() / accumulated_tarea; - ATLAS_TRACE_SCOPE("4") + // normalise weights of the intersection polygons to sum up to the area of the north octant + double norm_fac = north_octant.area() / coverage_north; for (int i = 0; i < n; i++) { for (int j = 0; j < m; j++) { - tgt_area[i * m + j] *= norm_fac; - dcov3 -= tgt_area[i * m + j]; + coverage_norm += i_north_area[i * m + j] * norm_fac; } } + EXPECT(north_octant.area() - M_PI_2 < EPS); // spherical area error for the north octant + EXPECT(south_octant.area() - M_PI_2 < EPS); // spherical area error for the south octant + EXPECT(north_octant.area() - coverage_north < m * EPS); // error polygon intersec. north + EXPECT(south_octant.area() - coverage_south < m * EPS); // error polygon intersec. north + EXPECT(north_octant.area() - coverage_norm < m * EPS); // error polygon intersec. north + auto err = std::max(north_octant.area() - coverage_north, south_octant.area() - coverage_south); + auto err_normed = north_octant.area() - coverage_norm; + Log::info() << "Octant coverting error : " << err << std::endl; + Log::info() << "Octant coverting error normed : " << err_normed << std::endl; +} - Log::info() << " dlat, dlon : " << dlat << ", " << dlon << "\n"; - Log::info() << " max (commutative_area) : " << darea << "\n"; - Log::info() << " dcov1 : " << dcov1 << "\n"; - Log::info() << " dcov2 : " << dcov2 << "\n"; - Log::info() << " dcov3 : " << dcov3 << "\n"; - Log::info() << " accumulated small polygon area : " << accumulated_tarea << "\n"; - Log::info() << " large polygon area : " << csp0.area() << "\n"; - EXPECT_APPROX_EQ(darea, 0., 1.e-10); - EXPECT_APPROX_EQ(accumulated_tarea, csp0.area(), 1.e-8); +CASE("intesection of epsilon-distorted polygons") { + const double eps = 1e-4; // degrees + const auto plg0 = make_polygon({{42.45283019, 50.48004426}, + {43.33333333, 49.70239033}, + {44.1509434, 50.48004426}, + {43.26923077, 51.2558069}}); + const auto plg1 = make_polygon({{42.45283019, 50.48004426 - eps}, + {43.33333333 + eps, 49.70239033}, + {44.1509434, 50.48004426 + eps}, + {43.26923077 - eps, 51.2558069}}); + const auto iplg_sol = make_polygon({{44.15088878324276, 50.48009332686897}, + {43.68455392823953, 50.89443301919586}, + {43.26918271448949, 51.25576215711414}, + {42.86876285000331, 50.87921047438197}, + {42.45288468219661, 50.47999711267543}, + {42.92307395320301, 50.06869923211562}, + {43.33338148668725, 49.70243705225555}, + {43.72937844034824, 50.08295071539503}}); + check_intersection(plg0, plg1, iplg_sol); } -CASE("edge cases") { +CASE("Edge cases in polygon intersections") { Log::info().precision(20); - SECTION("CS-LFR-256 -> H1280 problem polygon intersection") { + + SECTION("CS-LFR-256 -> H1280") { const auto plg0 = make_polygon({{-23.55468749999994, -41.11286269132660}, {-23.20312500000000, -41.18816845938357}, {-23.20312500000000, -40.83947225425061}, @@ -260,14 +412,15 @@ CASE("edge cases") { {-23.23828125000000, -40.81704944888558}, {-23.27343750000000, -40.77762996221442}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 2e-12); // can not take 1e-15 - EXPECT_EQ(iplg.size(), 3); - EXPECT_EQ(jplg.size(), 3); - EXPECT(iplg.equals(jplg, 5.e-7)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},16) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } - SECTION("CS-LFR-16 -> O32 problem polygon intersection") { + SECTION("CS-LFR-16 -> O32") { const auto plg0 = make_polygon({{174.3750000000001, -16.79832945594544}, {174.3750000000001, -11.19720014633353}, {168.7500000000000, -11.03919441545243}, @@ -276,25 +429,29 @@ CASE("edge cases") { {174.1935483870968, -15.34836475949100}, {177.0967741935484, -15.34836475949100}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-13); // can not take 1e-15 - EXPECT_EQ(iplg.size(), 3); - EXPECT_EQ(jplg.size(), 3); - EXPECT(iplg.equals(jplg, 1.e-11)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},16) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } - SECTION("CS-LFR-2 -> O48 problem polygon intersection") { + SECTION("CS-LFR-2 -> O48") { const auto plg0 = make_polygon({{180, -45}, {180, 0}, {135, 0}, {135, -35.26438968}}); const auto plg1 = make_polygon({{180, -23.31573073}, {180, -25.18098558}, {-177.75, -23.31573073}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-15); EXPECT_EQ(iplg.size(), 0); - EXPECT_EQ(jplg.size(), 0); - EXPECT(iplg.equals(jplg, 1.e-12)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},10) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + EXPECT_APPROX_EQ(iplg.area(), 0.); + Log::info().unindent(); } - SECTION("H128->H256 problem polygon intersection") { + SECTION("H128->H256") { const auto plg0 = make_polygon({{42.45283019, 50.48004426}, {43.33333333, 49.70239033}, {44.15094340, 50.48004426}, @@ -304,12 +461,13 @@ CASE("edge cases") { {44.15094340, 50.48004426}, {43.71428571, 50.86815893}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-15); - EXPECT_EQ(iplg.size(), 4); - EXPECT_EQ(jplg.size(), 4); - EXPECT(iplg.equals(jplg, 1.e-12)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},10) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } SECTION("intesection of epsilon-distorted polygons") { @@ -323,13 +481,104 @@ CASE("edge cases") { {44.1509434, 50.48004426 + eps}, {43.26923077 - eps, 51.2558069}}); auto iplg = plg0.intersect(plg1); - auto jplg = plg1.intersect(plg0); - EXPECT_APPROX_EQ(iplg.area(), jplg.area(), 1e-10); EXPECT_EQ(iplg.size(), 8); - EXPECT_EQ(jplg.size(), 8); - EXPECT(iplg.equals(jplg, 1.e-8)); - Log::info() << "Intersection area difference: " << std::abs(iplg.area() - jplg.area()) << "\n"; + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},10) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); + } + SECTION("one") { // TODO: remove, do not know what example + auto plg0 = make_polygon({{0,90},{67.463999999999998636,89.999899999085130275},{67.5,89.999899999085130275}}); + auto plg1 = make_polygon({{72,85.760587120443801723},{90,85.760587120443801723},{-90,85.760587120443801723}}); + auto iplg = plg0.intersect(plg1); + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg},20) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); } + + SECTION("debug") { // TODO: remove, do not know what example + auto plg0 = make_polygon({{168.599999999999994,-6.88524379355500038},{168.561872909699019,-7.16627415158899961},{168.862876254180634,-7.16627415158899961}}); + auto plg1 = make_polygon({{168.84,-6.84},{168.84,-7.2},{169.2,-7.2},{169.2,-6.84}}); + auto plg2 = make_polygon({{168.48,-6.84},{168.48,-7.2},{168.84,-7.2},{168.84,-6.84}}); + auto iplg = plg0.intersect(plg1); + auto iplg2 = plg0.intersect(plg2); + Log::info().indent(); + Log::info() << "source area: " << plg0.area() << "\n"; + Log::info() << "target area: " << plg1.area() << "\n"; + Log::info() << "inters area: " << iplg.area() << "\n"; + Log::info() << "json: \n" << to_json({plg0,plg1,iplg2,iplg,iplg2},20) << "\n"; + EXPECT(std::min(plg0.area(), plg1.area()) >= iplg.area()); + Log::info().unindent(); + } + + + SECTION("debug2") { // TODO: remove, do not know what example + auto plg0 = make_polygon({{168.599999999999994, -6.88524379355500038}, + {168.561872909699019, -7.16627415158899961}, + {168.862876254180634, -7.16627415158899961}}); + auto plg1 = make_polygon({{168.48, -6.84}, + {168.48, -7.20}, + {168.84, -7.20}, + {168.84, -6.84}}); + const auto iplg_sol = make_polygon({{168.600000000000, -6.885243793555000}, + {168.561872909699, -7.166274151589000}, + {168.840000000000, -7.166281023991750}, + {168.840000000000, -7.141837487873564}}); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O320c_O320n_tcell-2524029") { + auto plg0 = make_polygon({{ 16.497973615369710, 89.85157892074884}, + { 0.000000000000000, 89.87355342974176}, + {-54.000000000000010, 89.78487690721863}, + { -9.000000000000002, 89.84788464490568}}); + auto plg1 = make_polygon({{ 36.000000000000000, 89.78487690721863}, + {-54.000000000000010, 89.78487690721863}, + {-36.000000000000000, 89.78487690721863}}); + const auto iplg_sol = make_polygon(); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O128c_F128c_tcell-77536") { + auto plg0 = make_polygon({{157.5,-16.49119584364},{157.5,-17.192948774143},{158.203125,-17.192948774143},{158.203125,-16.49119584364}}); + auto plg1 = make_polygon({{157.7064220183486,-16.49119584364},{157.5,-17.192948774143},{158.3333333333333,-17.192948774143}}); + const auto iplg_sol = make_polygon({{157.7066386314724, -16.49143953797217}, + {157.7063506671565, -16.49143933951490}, + {157.5000000000000, -17.19294877414300}, + {158.2031250000000, -17.19294877414292}, + {158.2031250000000, -17.04778221576889}}); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O128c_F128c_tcell") { + auto plg0 = make_polygon({{135,-10.505756145244},{135,-11.906523334954},{136.40625,-11.906523334954},{136.40625,-10.505756145244}}); + auto plg1 = make_polygon({{135.7377049180328,-10.505756145244},{135,-11.906523334954},{136.5,-11.906523334954}}); + const auto iplg_sol = make_polygon({{135.7381225488238, -10.50652773728132}, + {135.7373006953013, -10.50652782622945}, + {135.0000000000000, -11.90652333495400}, + {136.4062500000000, -11.90652333495396}, + {136.4062500000000, -11.73509894855489}}); + check_intersection(plg0, plg1, iplg_sol); + } + + SECTION("O128c_F128c_tcell-2") { + auto plg0 = make_polygon({{134.296875,-10.877172064989},{134.296875,-11.578925065131},{135,-11.578925065131},{135,-10.877172064989}}); + auto plg1 = make_polygon({{134.6153846153846,-10.877172064989},{134.2241379310345,-11.578925065131},{135,-11.578925065131}}); + const auto iplg_sol = make_polygon({{134.6154929040914, -10.87737018871372}, + {134.6152744739456, -10.87737016536156}, + {134.2968750000000, -11.44875997313061}, + {134.2968749999999, -11.57892506513095}, + {135.0000000000000, -11.57892506513100}}); + check_intersection(plg0, plg1, iplg_sol); + } + } //-----------------------------------------------------------------------------