Skip to content

Commit

Permalink
Merge branch 'release/0.32.1'
Browse files Browse the repository at this point in the history 
* release/0.32.1:
  Version 0.32.1
  Remove deprecated use of sprintf
  Added (lon, lat) to (alpha, beta) transforms to cubed sphere projection (#116)
  • Loading branch information
@wdeconinck
wdeconinck committed Feb 9, 2023
2 parents 60b00c0 + 531b2cf commit f367722
Show file tree
Hide file tree
Showing 15 changed files with 235 additions and 81 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,10 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html

## [Unreleased]

## [0.32.1] - 2023-02-09
### Added
- Added (lon, lat) to (alpha, beta) transforms to cubed sphere projection

## [0.32.0] - 2023-01-23
### Added
- Added BlockStructuredColumns FunctionSpace
Expand Down Expand Up @@ -428,6 +432,7 @@ This project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html
## 0.13.0 - 2018-02-16

[Unreleased]: https://github.com/ecmwf/atlas/compare/master...develop
[0.32.1]: https://github.com/ecmwf/atlas/compare/0.32.0...0.32.1
[0.32.0]: https://github.com/ecmwf/atlas/compare/0.31.1...0.32.0
[0.31.1]: https://github.com/ecmwf/atlas/compare/0.31.0...0.31.1
[0.31.0]: https://github.com/ecmwf/atlas/compare/0.30.0...0.31.0
Expand Down
2 changes: 1 addition & 1 deletion VERSION
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
0.32.0
0.32.1
103 changes: 49 additions & 54 deletions src/atlas/interpolation/method/cubedsphere/CellFinder.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,105 +26,100 @@ CellFinder::CellFinder(const Mesh& mesh, const util::Config& config): mesh_{mesh
projection_ = &(csGrid.cubedSphereProjection());

// Get views to cell data.
const auto xyView = array::make_view<double, 2>(mesh_.cells().field("xy"));
const auto tijView = array::make_view<idx_t, 2>(mesh_.cells().field("tij"));
const auto lonlatView = array::make_view<double, 2>(mesh_.cells().field("lonlat"));
const auto haloView = array::make_view<int, 1>(mesh_.cells().halo());

// make points and payloads vectors.
auto points = std::array<std::vector<Point2>, 6>{};
auto payloads = std::array<std::vector<idx_t>, 6>{};
auto points = std::vector<PointLonLat>{};
auto payloads = std::vector<idx_t>{};

// Iterate over cells.
auto halo = config.getInt("halo", 0);
for (idx_t i = 0; i < mesh_.cells().size(); ++i) {
if (haloView(i) <= halo) {
const auto t = tijView(i, 0);
const auto xy = PointXY(xyView(i, XX), xyView(i, YY));

points[size_t(t)].push_back(projection_->xy2alphabeta(xy, t));
payloads[size_t(t)].push_back(i);
points.push_back(PointLonLat(lonlatView(i, LON), lonlatView(i, LAT)));
payloads.push_back(i);
}
}

// build cell kd-trees.
for (size_t t = 0; t < 6; ++t) {
trees_[t].build(points[t], payloads[t]);
}
// build cell kd-tree.
tree_.build(points, payloads);
}

CellFinder::Cell CellFinder::getCell(const PointXY& xy, size_t listSize, double edgeEpsilon, double epsilon) const {
CellFinder::Cell CellFinder::getCell(const PointLonLat& lonlat, size_t listSize, double edgeEpsilon, double epsilon) const {
// Convert xy to alphabeta and t;
const auto& tiles = projection_->getCubedSphereTiles();
const auto t = tiles.indexFromXY(xy.data());
const auto alphabeta = projection_->xy2alphabeta(xy, t);


// Get mesh nodes and connectivity table.
const auto nodeXyView = array::make_view<double, 2>(mesh_.nodes().xy());
const auto& nodeConnectivity = mesh_.cells().node_connectivity();
const auto nodeXyView = array::make_view<double, 2>(mesh_.nodes().xy());
const auto& nodeConnectivity = mesh_.cells().node_connectivity();

// Get four nearest cell-centres to xy.
const auto& tree = trees_[size_t(t)];
if (tree.size() == 0) {
return Cell{{}, Intersect{}.fail()};
if (tree_.size() == 0) {
return Cell{-1, {}, Intersect{}.fail()};
}

const auto valueList = tree.closestPoints(alphabeta, std::min(listSize, tree.size()));
const auto valueList = tree_.closestPoints(lonlat, std::min(listSize, tree_.size()));

// Get view to cell flags.
const auto cellFlagsView = array::make_view<int, 1>(mesh_.cells().flags());
const auto cellTijView = array::make_view<idx_t, 2>(mesh_.cells().field("tij"));

for (const auto& value : valueList) {
const auto i = value.payload();

// Ignore invalid cells.
if (!(cellFlagsView(i) & Topology::INVALID)) {
const auto row = nodeConnectivity.row(i);
if ((cellFlagsView(i) & Topology::INVALID)) {
break;
}

if (row.size() == 4) {
auto quadAlphabeta = std::array<Vector2D, 4>{};
for (size_t k = 0; k < 4; ++k) {
const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
quadAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
}
const auto t = cellTijView(i);
const auto row = nodeConnectivity.row(i);

const auto quad =
element::Quad2D(quadAlphabeta[0], quadAlphabeta[1], quadAlphabeta[2], quadAlphabeta[3]);
if (row.size() == 4) {
auto quadAlphabeta = std::array<Vector2D, 4>{};
for (size_t k = 0; k < 4; ++k) {
const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
quadAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
}

const auto quad =
element::Quad2D(quadAlphabeta[0], quadAlphabeta[1], quadAlphabeta[2], quadAlphabeta[3]);

const auto isect = quad.localRemap(alphabeta, edgeEpsilon, epsilon);
const auto alphabeta = projection_->lonlat2alphabeta(lonlat, t);
const auto isect = quad.localRemap(alphabeta, edgeEpsilon, epsilon);

if (isect) {
return Cell{{row(0), row(1), row(2), row(3)}, isect};
}
if (isect) {
return Cell{i, {row(0), row(1), row(2), row(3)}, isect};
}
}
else {
// Cell is triangle.
auto triagAlphabeta = std::array<Vector2D, 3>{};
for (size_t k = 0; k < 3; ++k) {
const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
triagAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
}
else {
// Cell is triangle.
auto triagAlphabeta = std::array<Vector2D, 3>{};
for (size_t k = 0; k < 3; ++k) {
const auto xyNode = PointXY(nodeXyView(row(k), XX), nodeXyView(row(k), YY));
triagAlphabeta[k] = Vector2D(projection_->xy2alphabeta(xyNode, t).data());
}

const auto triag = element::Triag2D(triagAlphabeta[0], triagAlphabeta[1], triagAlphabeta[2]);
const auto triag = element::Triag2D(triagAlphabeta[0], triagAlphabeta[1], triagAlphabeta[2]);

const auto isect = triag.intersects(alphabeta, edgeEpsilon, epsilon);
const auto alphabeta = projection_->lonlat2alphabeta(lonlat, t);
const auto isect = triag.intersects(alphabeta, edgeEpsilon, epsilon);

if (isect) {
return Cell{{row(0), row(1), row(2)}, isect};
}
if (isect) {
return Cell{i, {row(0), row(1), row(2)}, isect};
}
}
}

// Couldn't find a cell.
return Cell{{}, Intersect{}.fail()};
return Cell{-1, {}, Intersect{}.fail()};
}

CellFinder::Cell CellFinder::getCell(const PointLonLat& lonlat, size_t listSize, double edgeEpsilon,
CellFinder::Cell CellFinder::getCell(const PointXY& xy, size_t listSize, double edgeEpsilon,
double epsilon) const {
auto xy = PointXY(lonlat.data());
projection_->lonlat2xy(xy.data());
return getCell(xy, listSize, edgeEpsilon, epsilon);
const auto lonlat = projection_->lonlat(xy);
return getCell(lonlat, listSize, edgeEpsilon, epsilon);
}


Expand Down
4 changes: 3 additions & 1 deletion src/atlas/interpolation/method/cubedsphere/CellFinder.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "atlas/interpolation/method/Intersect.h"
#include "atlas/mesh/Mesh.h"
#include "atlas/util/Config.h"
#include "atlas/util/Geometry.h"
#include "atlas/util/KDTree.h"
#include "atlas/util/Point.h"

Expand All @@ -36,6 +37,7 @@ using namespace util;
class CellFinder {
public:
struct Cell {
idx_t idx;
std::vector<idx_t> nodes;
Intersect isect;
};
Expand All @@ -57,7 +59,7 @@ class CellFinder {
private:
Mesh mesh_{};
const projection::detail::CubedSphereProjectionBase* projection_{};
std::array<util::IndexKDTree2D, 6> trees_{};
util::IndexKDTree tree_{Geometry{}};
};

} // namespace cubedsphere
Expand Down
8 changes: 8 additions & 0 deletions src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,10 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
auto weights = std::vector<Triplet>{};
const auto ghostView = array::make_view<int, 1>(target_.ghost());
const auto lonlatView = array::make_view<double, 2>(target_.lonlat());
const auto tijView = array::make_view<idx_t, 2>(ncSource.mesh().cells().field("tij"));

// Make vector of tile indices for each target point (needed for vector field interpolation).
std::vector<idx_t> tileIndex{};

for (idx_t i = 0; i < target_.size(); ++i) {
if (!ghostView(i)) {
Expand All @@ -66,6 +70,7 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
std::to_string(i) + ".");
}

tileIndex.push_back(tijView(cell.idx, 0));
const auto& isect = cell.isect;
const auto& j = cell.nodes;

Expand Down Expand Up @@ -95,6 +100,9 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp
// fill sparse matrix and return.
Matrix A(target_.size(), source_.size(), weights);
setMatrix(A);

// Add tile index metadata to target.
target_->metadata().set("tile index", tileIndex);
}

void CubedSphereBilinear::print(std::ostream&) const {
Expand Down
4 changes: 2 additions & 2 deletions src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ namespace method {
/// accross cubed sphere tiles in (alpha, beta) coordinates.
/// Adding int "halo" (default 0) to the config controls the amount of
/// halo to consider when seraching for interpolation polygons. Adding
/// int "list_size" (defualt 4) will change the number of cells
/// int "list_size" (defualt 8) will change the number of cells
/// returned by the internal kd-tree search. Increasing both values
/// may resolve errors if setup method fails to find cells.
/// The automatic halo exchange in the execute method can be disabled
Expand All @@ -46,7 +46,7 @@ class CubedSphereBilinear : public Method {
FunctionSpace target_;

int halo_{0};
int listSize_{4};
int listSize_{8};
bool halo_exchange_{true};
};

Expand Down
10 changes: 6 additions & 4 deletions src/atlas/output/detail/GmshIO.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -221,7 +221,8 @@ std::vector<int> get_levels(int nlev, const Metadata& gmsh_options) {
std::string field_lev(const Field& field, int jlev) {
if (field.levels()) {
char str[6] = {0, 0, 0, 0, 0, 0};
std::sprintf(str, "[%03d]", jlev);
auto str_len = std::snprintf(str, sizeof(str), "[%03d]", jlev);
ATLAS_ASSERT(str_len == 5);
return std::string(str);
}
else {
Expand Down Expand Up @@ -388,8 +389,9 @@ void write_field_nodes(const Metadata& gmsh_options, const functionspace::NoFunc

// ----------------------------------------------------------------------------

void print_field_lev(char field_lev[], int jlev) {
std::sprintf(field_lev, "[%03d]", jlev);
void print_field_lev(char field_lev[], size_t size, int jlev) {
ATLAS_ASSERT(size > 5);
auto str_len = std::snprintf(field_lev, size, "[%03d]", jlev);
}

/* unused
Expand Down Expand Up @@ -433,7 +435,7 @@ void write_field_nodes(const Metadata& gmsh_options, const functionspace::Struct
char field_lev[6] = {0, 0, 0, 0, 0, 0};

if (field.levels()) {
print_field_lev(field_lev, jlev);
print_field_lev(field_lev, sizeof(field_lev), jlev);
}

out << "$NodeData\n";
Expand Down
Loading

0 comments on commit f367722

Please sign in to comment.