Commit fedb2260 authored by Uwe Schulzweida's avatar Uwe Schulzweida

Isosurface: renamed lev1 to levels.

parent e7d3a6bb
Pipeline #5590 passed with stages
in 14 minutes and 43 seconds
...@@ -24,37 +24,37 @@ ...@@ -24,37 +24,37 @@
#include "interpol.h" #include "interpol.h"
static void static void
isosurface(double isoval, long nlev1, const Varray<double> &lev1, const FieldVector &field3D, Field &field2D) isosurface(double isoval, int nlevels, const Varray<double> &levels, const FieldVector &field3D, Field &field2D)
{ {
const auto gridsize = gridInqSize(field3D[0].grid); const auto gridsize = gridInqSize(field3D[0].grid);
const auto missval = field3D[0].missval; const auto missval = field3D[0].missval;
auto &data2D = field2D.vec_d; auto &data2D = field2D.vec_d;
auto nmiss = field3D[0].nmiss; auto nmiss = field3D[0].nmiss;
for (long k = 1; k < nlev1; ++k) nmiss += field3D[k].nmiss; for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
data2D[i] = missval; data2D[i] = missval;
for (long k = 0; k < (nlev1 - 1); ++k) for (int k = 0; k < (nlevels - 1); ++k)
{ {
const double val1 = field3D[k].vec_d[i]; const auto val1 = field3D[k].vec_d[i];
const double val2 = field3D[k + 1].vec_d[i]; const auto val2 = field3D[k + 1].vec_d[i];
if (nmiss) if (nmiss)
{ {
const bool lmiss1 = DBL_IS_EQUAL(val1, missval); const bool lmiss1 = DBL_IS_EQUAL(val1, missval);
const bool lmiss2 = DBL_IS_EQUAL(val2, missval); const bool lmiss2 = DBL_IS_EQUAL(val2, missval);
if (lmiss1 && lmiss2) continue; if (lmiss1 && lmiss2) continue;
if (lmiss1 && IS_EQUAL(isoval, val2)) data2D[i] = lev1[k + 1]; if (lmiss1 && IS_EQUAL(isoval, val2)) data2D[i] = levels[k + 1];
if (lmiss2 && IS_EQUAL(isoval, val1)) data2D[i] = lev1[k]; if (lmiss2 && IS_EQUAL(isoval, val1)) data2D[i] = levels[k];
if (lmiss1 || lmiss2) continue; if (lmiss1 || lmiss2) continue;
} }
if ((isoval >= val1 && isoval <= val2) || (isoval >= val2 && isoval <= val1)) if ((isoval >= val1 && isoval <= val2) || (isoval >= val2 && isoval <= val1))
{ {
data2D[i] = IS_EQUAL(val1, val2) ? lev1[k] : intlin(isoval, lev1[k], val1, lev1[k + 1], val2); data2D[i] = IS_EQUAL(val1, val2) ? levels[k] : intlin(isoval, levels[k], val1, levels[k + 1], val2);
break; break;
} }
} }
...@@ -64,20 +64,20 @@ isosurface(double isoval, long nlev1, const Varray<double> &lev1, const FieldVec ...@@ -64,20 +64,20 @@ isosurface(double isoval, long nlev1, const Varray<double> &lev1, const FieldVec
} }
static void static void
layerValueMin(long nlev1, const Varray<double> &lev1, const FieldVector &field3D, Field &field2D) layerValueMin(int nlevels, const FieldVector &field3D, Field &field2D)
{ {
const auto gridsize = gridInqSize(field3D[0].grid); const auto gridsize = gridInqSize(field3D[0].grid);
const auto missval = field3D[0].missval; const auto missval = field3D[0].missval;
auto &data2D = field2D.vec_d; auto &data2D = field2D.vec_d;
auto nmiss = field3D[0].nmiss; auto nmiss = field3D[0].nmiss;
for (long k = 1; k < nlev1; ++k) nmiss += field3D[k].nmiss; for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
data2D[i] = missval; data2D[i] = missval;
for (long k = 0; k < nlev1; ++k) for (int k = 0; k < nlevels; ++k)
{ {
const auto val = field3D[k].vec_d[i]; const auto val = field3D[k].vec_d[i];
if (!DBL_IS_EQUAL(val, missval)) if (!DBL_IS_EQUAL(val, missval))
...@@ -92,20 +92,20 @@ layerValueMin(long nlev1, const Varray<double> &lev1, const FieldVector &field3D ...@@ -92,20 +92,20 @@ layerValueMin(long nlev1, const Varray<double> &lev1, const FieldVector &field3D
} }
static void static void
layerValueMax(long nlev1, const Varray<double> &lev1, const FieldVector &field3D, Field &field2D) layerValueMax(int nlevels, const FieldVector &field3D, Field &field2D)
{ {
const auto gridsize = gridInqSize(field3D[0].grid); const auto gridsize = gridInqSize(field3D[0].grid);
const auto missval = field3D[0].missval; const auto missval = field3D[0].missval;
auto &data2D = field2D.vec_d; auto &data2D = field2D.vec_d;
auto nmiss = field3D[0].nmiss; auto nmiss = field3D[0].nmiss;
for (long k = 1; k < nlev1; ++k) nmiss += field3D[k].nmiss; for (int k = 1; k < nlevels; ++k) nmiss += field3D[k].nmiss;
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
data2D[i] = missval; data2D[i] = missval;
for (long k = nlev1 - 1; k >= 0; --k) for (int k = nlevels - 1; k >= 0; --k)
{ {
const auto val = field3D[k].vec_d[i]; const auto val = field3D[k].vec_d[i];
if (!DBL_IS_EQUAL(val, missval)) if (!DBL_IS_EQUAL(val, missval))
...@@ -128,9 +128,6 @@ positiveIsDown(int zaxisID) ...@@ -128,9 +128,6 @@ positiveIsDown(int zaxisID)
void * void *
Isosurface(void *process) Isosurface(void *process)
{ {
int nrecs;
int zaxisID1 = -1;
cdoInitialize(process); cdoInitialize(process);
// clang-format off // clang-format off
...@@ -160,13 +157,12 @@ Isosurface(void *process) ...@@ -160,13 +157,12 @@ Isosurface(void *process)
const auto taxisID2 = taxisDuplicate(taxisID1); const auto taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2); vlistDefTaxis(vlistID2, taxisID2);
int zaxisID1 = -1;
const auto nzaxis = vlistNzaxis(vlistID1); const auto nzaxis = vlistNzaxis(vlistID1);
int nlevels = 0; for (int i = 0; i < nzaxis; i++)
int i;
for (i = 0; i < nzaxis; i++)
{ {
const auto zaxisID = vlistZaxis(vlistID1, i); const auto zaxisID = vlistZaxis(vlistID1, i);
nlevels = zaxisInqSize(zaxisID); const auto nlevels = zaxisInqSize(zaxisID);
if (zaxisInqType(zaxisID) != ZAXIS_HYBRID && zaxisInqType(zaxisID) != ZAXIS_HYBRID_HALF) if (zaxisInqType(zaxisID) != ZAXIS_HYBRID && zaxisInqType(zaxisID) != ZAXIS_HYBRID_HALF)
if (nlevels > 1) if (nlevels > 1)
{ {
...@@ -174,20 +170,20 @@ Isosurface(void *process) ...@@ -174,20 +170,20 @@ Isosurface(void *process)
break; break;
} }
} }
if (i == nzaxis) cdoAbort("No processable variable found!"); if (zaxisID1 == -1) cdoAbort("No processable variable found!");
const auto nlev1 = nlevels; const auto nlevels = zaxisInqSize(zaxisID1);
Varray<double> lev1(nlev1); Varray<double> levels(nlevels);
cdoZaxisInqLevels(zaxisID1, lev1.data()); cdoZaxisInqLevels(zaxisID1, levels.data());
const auto lpositive = !(lev1[0] < 0.0 && lev1[nlev1 - 1] < 0.0); const auto lpositive = !(levels[0] < 0.0 && levels[nlevels - 1] < 0.0);
const auto lreverse = (lev1[0] > lev1[nlev1 - 1]); const auto lreverse = (levels[0] > levels[nlevels - 1]);
auto bottomValueFunc = lreverse ? layerValueMax : layerValueMin; auto bottomValueFunc = lreverse ? layerValueMax : layerValueMin;
auto topValueFunc = lreverse ? layerValueMin : layerValueMax; auto topValueFunc = lreverse ? layerValueMin : layerValueMax;
if (lpositive && positiveIsDown(zaxisID1)) std::swap(bottomValueFunc, topValueFunc); if (lpositive && positiveIsDown(zaxisID1)) std::swap(bottomValueFunc, topValueFunc);
const auto zaxisIDsfc = zaxisFromName("surface"); const auto zaxisIDsfc = zaxisFromName("surface");
for (i = 0; i < nzaxis; i++) for (int i = 0; i < nzaxis; i++)
if (zaxisID1 == vlistZaxis(vlistID1, i)) vlistChangeZaxisIndex(vlistID2, i, zaxisIDsfc); if (zaxisID1 == vlistZaxis(vlistID1, i)) vlistChangeZaxisIndex(vlistID2, i, zaxisIDsfc);
const auto streamID2 = cdoOpenWrite(1); const auto streamID2 = cdoOpenWrite(1);
...@@ -207,8 +203,11 @@ Isosurface(void *process) ...@@ -207,8 +203,11 @@ Isosurface(void *process)
for (int varID = 0; varID < nvars; varID++) lvar3D[varID] = (varList1[varID].zaxisID == zaxisID1); for (int varID = 0; varID < nvars; varID++) lvar3D[varID] = (varList1[varID].zaxisID == zaxisID1);
int tsID = 0; int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID))) while (true)
{ {
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
taxisCopyTimestep(taxisID2, taxisID1); taxisCopyTimestep(taxisID2, taxisID1);
cdoDefTimestep(streamID2, tsID); cdoDefTimestep(streamID2, tsID);
...@@ -232,9 +231,9 @@ Isosurface(void *process) ...@@ -232,9 +231,9 @@ Isosurface(void *process)
{ {
field2.init(varList1[varID]); field2.init(varList1[varID]);
// clang-format off // clang-format off
if (operatorID == ISOSURFACE) isosurface(isoval, nlev1, lev1, vars1[varID], field2); if (operatorID == ISOSURFACE) isosurface(isoval, nlevels, levels, vars1[varID], field2);
else if (operatorID == BOTTOMVALUE) bottomValueFunc(nlev1, lev1, vars1[varID], field2); else if (operatorID == BOTTOMVALUE) bottomValueFunc(nlevels, vars1[varID], field2);
else if (operatorID == TOPVALUE) topValueFunc(nlev1, lev1, vars1[varID], field2); else if (operatorID == TOPVALUE) topValueFunc(nlevels, vars1[varID], field2);
// clang-format on // clang-format on
cdoDefRecord(streamID2, varID, 0); cdoDefRecord(streamID2, varID, 0);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment