Added remap_method_conserv.cc

src/Makefile.am

@@ -228,6 +228,8 @@ libcdo_la_SOURCES =		after_dvtrans.cc          \
 				remap_distwgt.cc          \
 				remap_grid_cell_search.cc \
 				remap_grid_cell_search.h  \
+				remap_method_conserv.cc   \
+				remap_method_conserv.h    \
 				remap_point_search.cc     \
 				remap_scrip_io.cc         \
 				remap_search_latbins.cc   \

src/lib/yac/test/remap_method_conserv.cc

+#include "remap_method_conserv.h"
+
+void
+cdo_compute_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell)
+{
+  auto &overlapCells = search.overlapCells;
+
+  // Do the clipping and get the cell for the overlapping area
+  yac_cell_clipping(numSearchCells, search.gridCells.data(), gridCell.yacGridCell, overlapCells.data());
+
+  // Get the partial areas for the overlapping regions
+  for (size_t i = 0; i < numSearchCells; ++i) search.partialAreas[i] = gridcell_area(overlapCells[i]);
+
+#ifdef VERBOSE
+  for (size_t i = 0; i < numSearchCells; ++i) cdo_print("overlap area : %lf", search.partialAreas[i]);
+#endif
+}
+
+static double
+get_edge_direction(const double *ref_corner, double *corner_a, double *corner_b)
+{
+  double edge_norm[3];
+  crossproduct_kahan(corner_a, corner_b, edge_norm);
+  normalise_vector(edge_norm);
+
+  // sine of the angle between the edge and the reference corner
+  double angle = edge_norm[0] * ref_corner[0] + edge_norm[1] * ref_corner[1] + edge_norm[2] * ref_corner[2];
+
+  // if the reference corner is directly on the edge
+  // (for small angles sin(x)==x)
+  if (fabs(angle) < yac_angle_tol) return 0.0;
+
+  return copysign(1.0, angle);
+}
+
+void
+cdo_compute_concave_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell)
+{
+  auto targetCell = gridCell.yacGridCell;
+  auto &overlapAreas = search.partialAreas;
+  auto &overlapCells = search.overlapCells;
+  auto &sourceCells = search.gridCells;
+
+  // common node point to all partial target cells
+  double *baseCorner = targetCell.coordinates_xyz[0];
+
+  // the triangulation algorithm only works for cells that only have great circle edges
+  enum yac_edge_type edgeTypes[3] = { YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE };
+  double coordinates_xyz[3][3] = { { -1.0, -1.0, -1.0 }, { -1.0, -1.0, -1.0 }, { -1.0, -1.0, -1.0 } };
+  for (int k = 0; k < 3; ++k) coordinates_xyz[0][k] = baseCorner[k];
+
+  // data structure to hold the triangles of the target cell
+  yac_grid_cell partialCell;
+  partialCell.array_size = 3;
+  partialCell.num_corners = 3;
+  partialCell.coordinates_xyz = coordinates_xyz;
+  partialCell.edge_type = edgeTypes;
+
+  // Do the clipping and get the cell for the overlapping area
+
+  for (size_t i = 0; i < numSearchCells; ++i) overlapAreas[i] = 0.0;
+
+  // for all triangles of the target cell
+  // (triangles a formed by first corner of the target cells and each edge of
+  //  the cell; the first and last edge of the cell already, contain the
+  //  first corner, therefore we can skip them)
+  for (size_t cornerIdx = 1; cornerIdx < targetCell.num_corners - 1; ++cornerIdx)
+    {
+      auto cornerA = targetCell.coordinates_xyz[cornerIdx];
+      auto cornerB = targetCell.coordinates_xyz[(cornerIdx + 1)];
+
+      // if the current edge has a length of zero
+      if (points_are_identically(cornerA, cornerB)) continue;
+
+      auto edgeDirection = get_edge_direction(baseCorner, cornerA, cornerB);
+
+      for (int k = 0; k < 3; ++k) partialCell.coordinates_xyz[1][k] = cornerA[k];
+      for (int k = 0; k < 3; ++k) partialCell.coordinates_xyz[2][k] = cornerB[k];
+
+      // clip the current target cell triangle with all source cells
+      yac_cell_clipping(numSearchCells, sourceCells.data(), partialCell, overlapCells.data());
+
+      // Get the partial areas for the overlapping regions as sum over the partial target cells.
+      for (size_t i = 0; i < numSearchCells; ++i)
+        {
+          if (overlapCells[i].num_corners == 0) continue;
+
+          overlapAreas[i] += gridcell_area(overlapCells[i]) * edgeDirection;
+        }
+    }
+
+  for (size_t i = 0; i < numSearchCells; ++i)
+    {
+      if (overlapAreas[i] < 0.0) overlapAreas[i] = -overlapAreas[i];
+    }
+
+#ifdef VERBOSE
+  for (size_t i = 0; i < numSearchCells; ++i) cdo_print("overlap area %zu: %lf", i, partialAreas[i]);
+#endif
+}

src/lib/yac/test/remap_method_conserv.h

+#ifndef REMAP_METHOD_CONSERV_H
+#define REMAP_METHOD_CONSERV_H
+
+// #include "remap_grid_cell_search.h"
+#include "cdo_cell_search.h"
+
+extern "C"
+{
+#include "../../lib/yac/clipping.h"
+#include "../../lib/yac/area.h"
+#include "../../lib/yac/geometry.h"
+}
+
+struct CellSearch
+{
+  enum yac_edge_type *edgeType = nullptr;
+  size_t numCellCorners = 0;
+  size_t maxCells = 0;
+  Varray<double> partialAreas;
+  Varray<yac_grid_cell> gridCells;
+  Varray<yac_grid_cell> overlapCells;
+
+  void
+  realloc(size_t numCells)
+  {
+    if (numCells > maxCells)
+      {
+        partialAreas.resize(numCells);
+        overlapCells.resize(numCells);
+        gridCells.resize(numCells);
+
+        for (size_t i = maxCells; i < numCells; ++i)
+          {
+            overlapCells[i].array_size = 0;
+            overlapCells[i].num_corners = 0;
+            overlapCells[i].edge_type = nullptr;
+            overlapCells[i].coordinates_xyz = nullptr;
+
+            gridCells[i].array_size = numCellCorners;
+            gridCells[i].num_corners = numCellCorners;
+            gridCells[i].edge_type = edgeType;
+            gridCells[i].coordinates_xyz = new double[numCellCorners][3];
+          }
+
+        maxCells = numCells;
+      }
+  }
+
+  void
+  free()
+  {
+    for (size_t i = 0; i < maxCells; ++i)
+      {
+        if (overlapCells[i].array_size > 0)
+          {
+            if (overlapCells[i].coordinates_xyz) std::free(overlapCells[i].coordinates_xyz);
+            if (overlapCells[i].edge_type) std::free(overlapCells[i].edge_type);
+          }
+
+        delete[] gridCells[i].coordinates_xyz;
+      }
+
+    varray_free(partialAreas);
+    varray_free(overlapCells);
+    varray_free(gridCells);
+  }
+};
+
+inline double
+gridcell_area(const yac_grid_cell &cell)
+{
+  return (cell.num_corners > 1) ? yac_huiliers_area(cell) : 0.0;
+}
+
+void cdo_compute_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell);
+void cdo_compute_concave_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell);
+
+#endif

src/remap_conserv.cc

@@ -11,74 +11,13 @@
 #include "process_int.h"
 #include "cdo_timer.h"
 #include "remap.h"
+#include "remap_method_conserv.h"
 #include "remap_store_link.h"
 #include "cdo_options.h"
 #include "progress.h"
 #include "cdo_omp.h"
 #include <mpim_grid.h>
 
-extern "C"
-{
-#include "lib/yac/clipping.h"
-#include "lib/yac/area.h"
-#include "lib/yac/geometry.h"
-}
-
-struct CellSearch
-{
-  enum yac_edge_type *edgeType = nullptr;
-  size_t numCellCorners = 0;
-  size_t maxCells = 0;
-  Varray<double> partialAreas;
-  Varray<struct yac_grid_cell> gridCells;
-  Varray<struct yac_grid_cell> overlapCells;
-};
-
-static void
-cellsearch_realloc(size_t numCells, CellSearch &search)
-{
-  if (numCells > search.maxCells)
-    {
-      search.partialAreas.resize(numCells);
-      search.overlapCells.resize(numCells);
-      search.gridCells.resize(numCells);
-
-      for (size_t i = search.maxCells; i < numCells; ++i)
-        {
-          search.overlapCells[i].array_size = 0;
-          search.overlapCells[i].num_corners = 0;
-          search.overlapCells[i].edge_type = nullptr;
-          search.overlapCells[i].coordinates_xyz = nullptr;
-
-          search.gridCells[i].array_size = search.numCellCorners;
-          search.gridCells[i].num_corners = search.numCellCorners;
-          search.gridCells[i].edge_type = search.edgeType;
-          search.gridCells[i].coordinates_xyz = new double[search.numCellCorners][3];
-        }
-
-      search.maxCells = numCells;
-    }
-}
-
-static void
-cellsearch_free(CellSearch &search)
-{
-  for (size_t i = 0; i < search.maxCells; ++i)
-    {
-      if (search.overlapCells[i].array_size > 0)
-        {
-          if (search.overlapCells[i].coordinates_xyz) free(search.overlapCells[i].coordinates_xyz);
-          if (search.overlapCells[i].edge_type) free(search.overlapCells[i].edge_type);
-        }
-
-      delete[] search.gridCells[i].coordinates_xyz;
-    }
-
-  varray_free(search.partialAreas);
-  varray_free(search.overlapCells);
-  varray_free(search.gridCells);
-}
-
 static void
 boundbox_from_corners1r(size_t ic, size_t nc, const Varray<double> &cornerLon, const Varray<double> &cornerLat, float *boundBox)
 {
@@ -121,114 +60,9 @@ boundbox_from_corners1r(size_t ic, size_t nc, const Varray<double> &cornerLon, c
   */
 }
 
-static inline double
-gridcell_area(const struct yac_grid_cell &cell)
-{
-  return yac_huiliers_area(cell);
-}
-
-static void
-cdo_compute_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell)
-{
-  auto &overlapCells = search.overlapCells;
-
-  // Do the clipping and get the cell for the overlapping area
-  yac_cell_clipping(numSearchCells, search.gridCells.data(), gridCell.yacGridCell, overlapCells.data());
-
-  // Get the partial areas for the overlapping regions
-  for (size_t i = 0; i < numSearchCells; ++i) search.partialAreas[i] = gridcell_area(overlapCells[i]);
-
-#ifdef VERBOSE
-  for (size_t i = 0; i < numSearchCells; ++i) cdo_print("overlap area : %lf", search.partialAreas[i]);
-#endif
-}
-
-static double
-get_edge_direction(const double *ref_corner, double *corner_a, double *corner_b)
-{
-  double edge_norm[3];
-  crossproduct_kahan(corner_a, corner_b, edge_norm);
-  normalise_vector(edge_norm);
-
-  // sine of the angle between the edge and the reference corner
-  double angle = edge_norm[0] * ref_corner[0] + edge_norm[1] * ref_corner[1] + edge_norm[2] * ref_corner[2];
-
-  // if the reference corner is directly on the edge
-  // (for small angles sin(x)==x)
-  if (fabs(angle) < yac_angle_tol) return 0.0;
-
-  return copysign(1.0, angle);
-}
-
-static void
-cdo_compute_concave_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell)
-{
-  auto targetCell = gridCell.yacGridCell;
-  auto &overlapAreas = search.partialAreas;
-  auto &overlapCells = search.overlapCells;
-  auto &sourceCells = search.gridCells;
-
-  // common node point to all partial target cells
-  double *baseCorner = targetCell.coordinates_xyz[0];
-
-  // the triangulation algorithm only works for cells that only have great circle edges
-  enum yac_edge_type edgeTypes[3] = { YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE };
-  double coordinates_xyz[3][3] = { { -1.0, -1.0, -1.0 }, { -1.0, -1.0, -1.0 }, { -1.0, -1.0, -1.0 } };
-  for (int k = 0; k < 3; ++k) coordinates_xyz[0][k] = baseCorner[k];
-
-  // data structure to hold the triangles of the target cell
-  struct yac_grid_cell partialCell;
-  partialCell.array_size = 3;
-  partialCell.num_corners = 3;
-  partialCell.coordinates_xyz = coordinates_xyz;
-  partialCell.edge_type = edgeTypes;
-
-  // Do the clipping and get the cell for the overlapping area
-
-  for (size_t i = 0; i < numSearchCells; ++i) overlapAreas[i] = 0.0;
-
-  // for all triangles of the target cell
-  // (triangles a formed by first corner of the target cells and each edge of
-  //  the cell; the first and last edge of the cell already, contain the
-  //  first corner, therefore we can skip them)
-  for (size_t cornerIdx = 1; cornerIdx < targetCell.num_corners - 1; ++cornerIdx)
-    {
-      auto cornerA = targetCell.coordinates_xyz[cornerIdx];
-      auto cornerB = targetCell.coordinates_xyz[(cornerIdx + 1)];
-
-      // if the current edge has a length of zero
-      if (points_are_identically(cornerA, cornerB)) continue;
-
-      auto edgeDirection = get_edge_direction(baseCorner, cornerA, cornerB);
-
-      for (int k = 0; k < 3; ++k) partialCell.coordinates_xyz[1][k] = cornerA[k];
-      for (int k = 0; k < 3; ++k) partialCell.coordinates_xyz[2][k] = cornerB[k];
-
-      // clip the current target cell triangle with all source cells
-      yac_cell_clipping(numSearchCells, sourceCells.data(), partialCell, overlapCells.data());
-
-      // Get the partial areas for the overlapping regions as sum over the partial target cells.
-      for (size_t i = 0; i < numSearchCells; ++i)
-        {
-          if (overlapCells[i].num_corners == 0) continue;
-
-          overlapAreas[i] += gridcell_area(overlapCells[i]) * edgeDirection;
-        }
-    }
-
-  for (size_t i = 0; i < numSearchCells; ++i)
-    {
-      if (overlapAreas[i] < 0.0) overlapAreas[i] = -overlapAreas[i];
-    }
-
-#ifdef VERBOSE
-  for (size_t i = 0; i < numSearchCells; ++i) cdo_print("overlap area %zu: %lf", i, partialAreas[i]);
-#endif
-}
-
 static void
 set_cell_coordinates_yac(RemapGridType remapGridType, size_t cellIndex, size_t numCorners, const RemapGrid *remapGrid,
-                         const struct yac_grid_cell &yacGridCell, double *x, double *y)
+                         const yac_grid_cell &yacGridCell, double *x, double *y)
 {
   auto storeXY = (x && y);
   auto xyz = yacGridCell.coordinates_xyz;
@@ -256,8 +90,8 @@ set_cell_coordinates_yac(RemapGridType remapGridType, size_t cellIndex, size_t n
     }
   else
     {
-      auto cell_corner_lon = &remapGrid->cell_corner_lon[cellIndex * numCorners];
-      auto cell_corner_lat = &remapGrid->cell_corner_lat[cellIndex * numCorners];
+      const auto cell_corner_lon = &remapGrid->cell_corner_lon[cellIndex * numCorners];
+      const auto cell_corner_lat = &remapGrid->cell_corner_lat[cellIndex * numCorners];
       for (size_t i = 0; i < numCorners; ++i) gcLLtoXYZ(cell_corner_lon[i], cell_corner_lat[i], xyz[i]);
 
       if (storeXY)
@@ -720,7 +554,7 @@ remap_conserv_weights(RemapSearch &remapSearch, RemapVars &rv)
 
       // Create search arrays
 
-      cellsearch_realloc(numSearchCells, cellSearch);
+      cellSearch.realloc(numSearchCells);
 
       if (remapSearch.gcs.xyzCoords)
         set_coordinates_yac(remapSearch.gcs.xyzCoords, numSearchCells, indices, numCellCorners, cellSearch.gridCells);
@@ -773,7 +607,7 @@ remap_conserv_weights(RemapSearch &remapSearch, RemapVars &rv)
   // Finished with all cells: deallocate search arrays
   for (auto ompthID = 0; ompthID < Threading::ompNumThreads; ++ompthID)
     {
-      cellsearch_free(cellSearch2[ompthID]);
+      cellSearch2[ompthID].free();
       gridcell_free_yac(tgtGridCell2[ompthID]);
     }
 
@@ -939,7 +773,7 @@ remap_conserv(const Varray<T1> &srcArray, Varray<T2> &tgtArray, double srcMissva
 
       // Create search arrays
 
-      cellsearch_realloc(numSearchCells, cellSearch);
+      cellSearch.realloc(numSearchCells);
 
       if (remapSearch.gcs.xyzCoords)
         set_coordinates_yac(remapSearch.gcs.xyzCoords, numSearchCells, indices, numCellCorners, cellSearch.gridCells);
@@ -984,7 +818,7 @@ remap_conserv(const Varray<T1> &srcArray, Varray<T2> &tgtArray, double srcMissva
 
   for (auto ompthID = 0; ompthID < Threading::ompNumThreads; ++ompthID)
     {
-      cellsearch_free(cellSearch2[ompthID]);
+      cellSearch2[ompthID].free();
       gridcell_free_yac(tgtGridCell2[ompthID]);
     }
 
@@ -1134,7 +968,7 @@ remap_weights_zonal_mean(int gridID1, int gridID2, Varray2D<size_t> &remapIndice
       search.numCellCorners = nv1;
       search.edgeType = srcEdgeType;
 
-      cellsearch_realloc(numSearchCells, search);
+      search.realloc(numSearchCells);
 
       for (size_t j = 0; j < numSearchCells; ++j)
         {
@@ -1164,7 +998,7 @@ remap_weights_zonal_mean(int gridID1, int gridID2, Varray2D<size_t> &remapIndice
       remapWeights[i2].resize(numWeights);
       for (size_t i = 0; i < numWeights; ++i) remapWeights[i2][i] = partialWeights[i];
 
-      cellsearch_free(search);
+      search.free();
     }
 }
 

src/remap_method_conserv.cc

+#include "remap_method_conserv.h"
+
+void
+cdo_compute_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell)
+{
+  auto &overlapCells = search.overlapCells;
+
+  // Do the clipping and get the cell for the overlapping area
+  yac_cell_clipping(numSearchCells, search.gridCells.data(), gridCell.yacGridCell, overlapCells.data());
+
+  // Get the partial areas for the overlapping regions
+  for (size_t i = 0; i < numSearchCells; ++i) search.partialAreas[i] = gridcell_area(overlapCells[i]);
+
+#ifdef VERBOSE
+  for (size_t i = 0; i < numSearchCells; ++i) cdo_print("overlap area : %lf", search.partialAreas[i]);
+#endif
+}
+
+static double
+get_edge_direction(const double *ref_corner, double *corner_a, double *corner_b)
+{
+  double edge_norm[3];
+  crossproduct_kahan(corner_a, corner_b, edge_norm);
+  normalise_vector(edge_norm);
+
+  // sine of the angle between the edge and the reference corner
+  double angle = edge_norm[0] * ref_corner[0] + edge_norm[1] * ref_corner[1] + edge_norm[2] * ref_corner[2];
+
+  // if the reference corner is directly on the edge
+  // (for small angles sin(x)==x)
+  if (fabs(angle) < yac_angle_tol) return 0.0;
+
+  return copysign(1.0, angle);
+}
+
+void
+cdo_compute_concave_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell)
+{
+  auto targetCell = gridCell.yacGridCell;
+  auto &overlapAreas = search.partialAreas;
+  auto &overlapCells = search.overlapCells;
+  auto &sourceCells = search.gridCells;
+
+  // common node point to all partial target cells
+  double *baseCorner = targetCell.coordinates_xyz[0];
+
+  // the triangulation algorithm only works for cells that only have great circle edges
+  enum yac_edge_type edgeTypes[3] = { YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE, YAC_GREAT_CIRCLE_EDGE };
+  double coordinates_xyz[3][3] = { { -1.0, -1.0, -1.0 }, { -1.0, -1.0, -1.0 }, { -1.0, -1.0, -1.0 } };
+  for (int k = 0; k < 3; ++k) coordinates_xyz[0][k] = baseCorner[k];
+
+  // data structure to hold the triangles of the target cell
+  yac_grid_cell partialCell;
+  partialCell.array_size = 3;
+  partialCell.num_corners = 3;
+  partialCell.coordinates_xyz = coordinates_xyz;
+  partialCell.edge_type = edgeTypes;
+
+  // Do the clipping and get the cell for the overlapping area
+
+  for (size_t i = 0; i < numSearchCells; ++i) overlapAreas[i] = 0.0;
+
+  // for all triangles of the target cell
+  // (triangles a formed by first corner of the target cells and each edge of
+  //  the cell; the first and last edge of the cell already, contain the
+  //  first corner, therefore we can skip them)
+  for (size_t cornerIdx = 1; cornerIdx < targetCell.num_corners - 1; ++cornerIdx)
+    {
+      auto cornerA = targetCell.coordinates_xyz[cornerIdx];
+      auto cornerB = targetCell.coordinates_xyz[(cornerIdx + 1)];
+
+      // if the current edge has a length of zero
+      if (points_are_identically(cornerA, cornerB)) continue;
+
+      auto edgeDirection = get_edge_direction(baseCorner, cornerA, cornerB);
+
+      for (int k = 0; k < 3; ++k) partialCell.coordinates_xyz[1][k] = cornerA[k];
+      for (int k = 0; k < 3; ++k) partialCell.coordinates_xyz[2][k] = cornerB[k];
+
+      // clip the current target cell triangle with all source cells
+      yac_cell_clipping(numSearchCells, sourceCells.data(), partialCell, overlapCells.data());
+
+      // Get the partial areas for the overlapping regions as sum over the partial target cells.
+      for (size_t i = 0; i < numSearchCells; ++i)
+        {
+          if (overlapCells[i].num_corners == 0) continue;
+
+          overlapAreas[i] += gridcell_area(overlapCells[i]) * edgeDirection;
+        }
+    }
+
+  for (size_t i = 0; i < numSearchCells; ++i)
+    {
+      if (overlapAreas[i] < 0.0) overlapAreas[i] = -overlapAreas[i];
+    }
+
+#ifdef VERBOSE
+  for (size_t i = 0; i < numSearchCells; ++i) cdo_print("overlap area %zu: %lf", i, partialAreas[i]);
+#endif
+}

src/remap_method_conserv.h

+#ifndef REMAP_METHOD_CONSERV_H
+#define REMAP_METHOD_CONSERV_H
+
+#include "remap_grid_cell_search.h"
+
+extern "C"
+{
+#include "lib/yac/clipping.h"
+#include "lib/yac/area.h"
+#include "lib/yac/geometry.h"
+}
+
+struct CellSearch
+{
+  enum yac_edge_type *edgeType = nullptr;
+  size_t numCellCorners = 0;
+  size_t maxCells = 0;
+  Varray<double> partialAreas;
+  Varray<yac_grid_cell> gridCells;
+  Varray<yac_grid_cell> overlapCells;
+
+  void
+  realloc(size_t numCells)
+  {
+    if (numCells > maxCells)
+      {
+        partialAreas.resize(numCells);
+        overlapCells.resize(numCells);
+        gridCells.resize(numCells);
+
+        for (size_t i = maxCells; i < numCells; ++i)
+          {
+            overlapCells[i].array_size = 0;
+            overlapCells[i].num_corners = 0;
+            overlapCells[i].edge_type = nullptr;
+            overlapCells[i].coordinates_xyz = nullptr;
+
+            gridCells[i].array_size = numCellCorners;
+            gridCells[i].num_corners = numCellCorners;
+            gridCells[i].edge_type = edgeType;
+            gridCells[i].coordinates_xyz = new double[numCellCorners][3];
+          }
+
+        maxCells = numCells;
+      }
+  }
+
+  void
+  free()
+  {
+    for (size_t i = 0; i < maxCells; ++i)
+      {
+        if (overlapCells[i].array_size > 0)
+          {
+            if (overlapCells[i].coordinates_xyz) std::free(overlapCells[i].coordinates_xyz);
+            if (overlapCells[i].edge_type) std::free(overlapCells[i].edge_type);
+          }
+
+        delete[] gridCells[i].coordinates_xyz;
+      }
+
+    varray_free(partialAreas);
+    varray_free(overlapCells);
+    varray_free(gridCells);
+  }
+};
+
+inline double
+gridcell_area(const yac_grid_cell &cell)
+{
+  return (cell.num_corners > 1) ? yac_huiliers_area(cell) : 0.0;
+}
+
+void cdo_compute_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell);
+void cdo_compute_concave_overlap_areas(size_t numSearchCells, CellSearch &search, const GridCell &gridCell);
+
+#endif