adds test for yac_fdef_grid_cloud

src/core/basic_grid.h

@@ -136,6 +136,10 @@ void yac_basic_grid_to_file_parallel(
 void yac_basic_grid_compute_cell_areas(
   struct yac_basic_grid * grid, double * cell_areas);
 
+void yac_rotate_coordinates(
+  yac_coordinate_pointer coordinates, size_t num_coordinates,
+  double north_pole[3]);
+
 // YAC PUBLIC HEADER STOP
 
 #endif // BASIC_GRID_H

src/core/grid_reg2d_rot.c

@@ -7,6 +7,43 @@
 #include "geometry.h"
 #include "utils_common.h"
 
+void yac_rotate_coordinates(
+  yac_coordinate_pointer coordinates, size_t num_coordinates,
+  double north_pole[3]) {
+
+  // compute angle between original and new north pol
+  struct sin_cos_angle rot_angle =
+    get_vector_angle_2(north_pole, (double[]){0.0, 0.0, 1.0});
+
+  YAC_ASSERT_F(
+    compare_angles(rot_angle, SIN_COS_TOL) > 0,
+    "ERROR(yac_rotate_coordinates): "
+    "new north pole is the original north pole "
+    "(new north pole (%.3lf; %.3lf; %.3lf); angle (sin: %e, cos: %e)",
+    north_pole[0], north_pole[1], north_pole[2], rot_angle.sin, rot_angle.cos);
+
+  YAC_ASSERT_F(
+    compare_angles(rot_angle, SIN_COS_M_PI_2) <= 0,
+    "ERROR(yac_rotate_coordinates): "
+    "new north pole is on the southern hemisphere "
+    "(new north pole (%.3lf; %.3lf; %.3lf); angle (sin: %e, cos: %e)",
+    north_pole[0], north_pole[1], north_pole[2], rot_angle.sin, rot_angle.cos);
+
+  // compute rotatation axis
+  // (has an angle of 90 degree between original and new north pole)
+  double rot_axis[3];
+  crossproduct_kahan((double[]){0.0, 0.0, 1.0}, north_pole, rot_axis);
+  normalise_vector(rot_axis);
+
+  for (size_t i = 0; i < num_coordinates; ++i) {
+    double temp[3];
+    rotate_vector2(rot_axis, rot_angle, coordinates[i], temp);
+    coordinates[i][0] = temp[0];
+    coordinates[i][1] = temp[1];
+    coordinates[i][2] = temp[2];
+  }
+}
+
 static struct yac_basic_grid_data yac_generate_basic_grid_data_reg_2d_rot_(
   size_t nbr_vertices[2], int cyclic[2],
   double *lon_vertices, double *lat_vertices,
@@ -26,41 +63,16 @@ static struct yac_basic_grid_data yac_generate_basic_grid_data_reg_2d_rot_(
     (num_cells_2d[0] + (cyclic[0]?0:1)) * num_cells_2d[1] +
     num_cells_2d[0] * (num_cells_2d[1] + 1);
 
-  // compute angle between original and new north pol
-  double north_pol[3];
-  LLtoXYZ_ptr(north_pol_lon, north_pol_lat, north_pol);
-  struct sin_cos_angle rot_angle =
-    get_vector_angle_2(north_pol, (double[]){0.0, 0.0, 1.0});
-
-  YAC_ASSERT_F(
-    compare_angles(rot_angle, SIN_COS_TOL) > 0,
-    "ERROR(yac_generate_basic_grid_data_reg_2d_rot): "
-    "new north pole is the original north pole "
-    "(new north pole (%.3lf; %.3lf); angle (sin: %e, cos: %e)",
-    north_pol_lon, north_pol_lat, rot_angle.sin, rot_angle.cos);
-
-  YAC_ASSERT_F(
-    compare_angles(rot_angle, SIN_COS_M_PI_2) <= 0,
-    "ERROR(yac_generate_basic_grid_data_reg_2d_rot): "
-    "new north pole is on the southern hemisphere "
-    "(new north pole (%.3lf; %.3lf); angle (sin: %e, cos: %e)",
-    north_pol_lon, north_pol_lat, rot_angle.sin, rot_angle.cos);
-
-  // compute rotatation axis
-  // (has an angle of 90 degree between original and new north pole)
-  double rot_axis[3];
-  crossproduct_kahan((double[]){0.0, 0.0, 1.0}, north_pol, rot_axis);
-  normalise_vector(rot_axis);
-
   yac_coordinate_pointer vertex_coordinates =
     xmalloc(num_vertices * sizeof(*vertex_coordinates));
-  for (size_t i = 0, k = 0; i < num_vertices_2d[1]; ++i) {
-    for (size_t j = 0; j < num_vertices_2d[0]; ++j, ++k) {
-      double temp[3];
-      LLtoXYZ_ptr(lon_vertices[j], lat_vertices[i], temp);
-      rotate_vector2(rot_axis, rot_angle, temp, &(vertex_coordinates[k][0]));
-    }
-  }
+  for (size_t i = 0, k = 0; i < num_vertices_2d[1]; ++i)
+    for (size_t j = 0; j < num_vertices_2d[0]; ++j, ++k)
+      LLtoXYZ_ptr(lon_vertices[j], lat_vertices[i], vertex_coordinates[k]);
+
+  // rotate all coordinates according to provided north pole
+  double north_pole[3];
+  LLtoXYZ_ptr(north_pol_lon, north_pol_lat, north_pole);
+  yac_rotate_coordinates(vertex_coordinates, num_vertices, north_pole);
 
   enum yac_edge_type * edge_type = xmalloc(num_edges * sizeof(*edge_type));
   for (size_t i = 0; i < num_edges; ++i) edge_type[i] = YAC_GREAT_CIRCLE_EDGE;

tests/mci/test_def_grid.inc

@@ -502,6 +502,34 @@ program test_def_grid
   print *, ' def_grid returned grid_id ', grid_id
   call test ( grid_id /= -99 )
 
+  ! cloud grid
+
+  !  5       6
+  !
+  !
+  !
+  !  3       4
+  !
+  !
+  !
+  !  1       2
+
+  x_vertices(1) = -0.5; y_vertices(1) =  1.0
+  x_vertices(2) =  0.5; y_vertices(2) =  1.0
+  x_vertices(3) = -0.5; y_vertices(3) =  0.0
+  x_vertices(4) =  0.5; y_vertices(4) =  0.0
+  x_vertices(5) = -0.5; y_vertices(5) = -1.0
+  x_vertices(6) =  0.5; y_vertices(6) = -1.0
+
+  call yac_fdef_grid ( 'grid_cloud',          &
+                       nbr_vertices,          &
+                       x_vertices,            &
+                       y_vertices,            &
+                       grid_id)
+
+  print *, ' def_grid returned grid_id ', grid_id
+  call test ( grid_id /= -99 )
+
   call yac_ffinalize
 
   call stop_test