bes  Updated for version 3.20.6
TwoDMeshTopology.cc
1 // -*- mode: c++; c-basic-offset:4 -*-
2 
3 // This file is part of libdap, A C++ implementation of the OPeNDAP Data
4 // Access Protocol.
5 
6 // Copyright (c) 2002,2003,2011,2012 OPeNDAP, Inc.
7 // Authors: Nathan Potter <ndp@opendap.org>
8 // James Gallagher <jgallagher@opendap.org>
9 // Scott Moe <smeest1@gmail.com>
10 // Bill Howe <billhowe@cs.washington.edu>
11 //
12 // This library is free software; you can redistribute it and/or
13 // modify it under the terms of the GNU Lesser General Public
14 // License as published by the Free Software Foundation; either
15 // version 2.1 of the License, or (at your option) any later version.
16 //
17 // This library is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 // Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License along with this library; if not, write to the Free Software
24 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 //
26 // You can contact OPeNDAP, Inc. at PO Box 112, Saunderstown, RI. 02874-0112.
27 
28 #include "config.h"
29 
30 #include <sstream>
31 #include <vector>
32 #include <algorithm>
33 #include <functional>
34 
35 #include <gridfields/type.h>
36 #include <gridfields/gridfield.h>
37 #include <gridfields/grid.h>
38 #include <gridfields/implicit0cells.h>
39 #include <gridfields/gridfieldoperator.h>
40 #include <gridfields/restrict.h>
41 #include <gridfields/refrestrict.h>
42 #include <gridfields/array.h>
43 #include <gridfields/GFError.h>
44 
45 #include "BaseType.h"
46 #include "Int32.h"
47 #include "Float64.h"
48 #include "Array.h"
49 #include "util.h"
50 
51 #include "ugrid_utils.h"
52 //#include "NDimensionalArray.h"
53 #include "MeshDataVariable.h"
54 #include "TwoDMeshTopology.h"
55 
56 #include "BESDebug.h"
57 #include "BESError.h"
58 
59 #ifdef NDEBUG
60 #undef BESDEBUG
61 #define BESDEBUG( x, y )
62 #endif
63 
64 using namespace std;
65 using namespace libdap;
66 using namespace ugrid;
67 
68 namespace ugrid {
69 
70 /* not used. faceCoordinateNames(0), */
71 TwoDMeshTopology::TwoDMeshTopology() :
72  d_meshVar(0), nodeCoordinateArrays(0), nodeCount(0), faceNodeConnectivityArray(0), faceCount(0), faceCoordinateArrays(
73  0), gridTopology(0), d_inputGridField(0), resultGridField(0), fncCellArray(0), _initialized(false)
74 {
75  rangeDataArrays = new vector<MeshDataVariable *>();
76  sharedIntArrays = new vector<int *>();
77  sharedFloatArrays = new vector<float *>();
78 }
79 
80 TwoDMeshTopology::~TwoDMeshTopology()
81 {
82  BESDEBUG("ugrid", "~TwoDMeshTopology() - BEGIN" << endl);
83  BESDEBUG("ugrid", "~TwoDMeshTopology() - (" << this << ")" << endl);
84 
85  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::GridField 'resultGridField'." << endl);
86  delete resultGridField;
87 
88  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::GridField 'inputGridField'." << endl);
89  delete d_inputGridField;
90 
91  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::Grid 'gridTopology'." << endl);
92  delete gridTopology;
93 
94  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::Arrays..." << endl);
95  for (vector<GF::Array *>::iterator it = gfArrays.begin(); it != gfArrays.end(); ++it) {
96  GF::Array *gfa = *it;
97  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF:Array '" << gfa->getName() << "'" << endl);
98  delete gfa;
99  }
100  BESDEBUG("ugrid", "~TwoDMeshTopology() - GF::Arrays deleted." << endl);
101 
102  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting sharedIntArrays..." << endl);
103  for (vector<int *>::iterator it = sharedIntArrays->begin(); it != sharedIntArrays->end(); ++it) {
104  int *i = *it;
105  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting integer array '" << i << "'" << endl);
106  delete[] i;
107  }
108  delete sharedIntArrays;
109 
110  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting sharedFloatArrays..." << endl);
111  for (vector<float *>::iterator it = sharedFloatArrays->begin(); it != sharedFloatArrays->end(); ++it) {
112  float *f = *it;
113  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting float array '" << f << "'" << endl);
114  delete[] f;
115  }
116  delete sharedFloatArrays;
117 
118  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting range data vector" << endl);
119  delete rangeDataArrays;
120  BESDEBUG("ugrid", "~TwoDMeshTopology() - Range data vector deleted." << endl);
121 
122  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting vector of node coordinate arrays." << endl);
123  delete nodeCoordinateArrays;
124 
125  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting vector of face coordinate arrays." << endl);
126  delete faceCoordinateArrays;
127 
128  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting face node connectivity cell array (GF::Node's)." << endl);
129  delete[] fncCellArray;
130 
131  BESDEBUG("ugrid", "~TwoDMeshTopology() - END" << endl);
132 }
133 
137 void TwoDMeshTopology::init(string meshVarName, DDS *dds)
138 {
139  if (_initialized) return;
140 
141  d_meshVar = dds->var(meshVarName);
142 
143  if (!d_meshVar) throw Error("Unable to locate variable: " + meshVarName);
144 
145  dimension = getAttributeValue(d_meshVar, UGRID_TOPOLOGY_DIMENSION);
146  if (dimension.empty()) dimension = getAttributeValue(d_meshVar, UGRID_DIMENSION);
147 
148  if (dimension.empty()) {
149  string msg = "ugr5(): The mesh topology variable '" + d_meshVar->name()
150  + "' is missing the required attribute named '" + UGRID_TOPOLOGY_DIMENSION + "'";
151  BESDEBUG("ugrid", "TwoDMeshTopology::init() - " << msg);
152  throw Error(msg);
153  }
154 
155  // Retrieve the node coordinate arrays for the mesh
156  ingestNodeCoordinateArrays(d_meshVar, dds);
157 
158  // Why would we retrieve this? I think Bill said that this needs to be recomputed after a restrict operation.
159  // @TODO Verify that Bill actually said that this needs to be recomputed.
160  // Retrieve the face coordinate arrays (if any) for the mesh
161  ingestFaceCoordinateArrays(d_meshVar, dds);
162 
163  // Inspect and QC the face node connectivity array for the mesh
164  ingestFaceNodeConnectivityArray(d_meshVar, dds);
165 
166  // Load d_meshVar with some data value - needed because this variable
167  // will be returned as part of the result and DAP does not allow
168  // empty variables. Since this code is designed to work with the UGrid
169  // specification being developed as an extension to (or in conjunction
170  // with) CF, I am assuming the UGrids will always be netCDF files and
171  // that the dummy mesh variable will always get a value when read()
172  // is called because netCDF guarantees that reading missing values
173  // returns the 'fill value'.
174  // See https://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c/Fill-Values.html
175  // jhrg 4/15/15
176  try {
177  d_meshVar->read(); // read() sets read_p to true
178  }
179  catch (Error &e) {
180  throw Error(malformed_expr,
181  "ugr5(): While trying to read the UGrid mesh variable, an error occurred: " + e.get_error_message());
182  }
183  catch (std::exception &e) {
184  throw Error(malformed_expr,
185  string("ugr5(): While trying to read the UGrid mesh variable, an error occurred: ") + e.what());
186  }
187 
188  _initialized = true;
189 }
190 
191 void TwoDMeshTopology::setNodeCoordinateDimension(MeshDataVariable *mdv)
192 {
193  BESDEBUG("ugrid", "TwoDMeshTopology::setNodeCoordinateDimension() - BEGIN" << endl);
194  libdap::Array *dapArray = mdv->getDapArray();
195  libdap::Array::Dim_iter ait1;
196 
197  for (ait1 = dapArray->dim_begin(); ait1 != dapArray->dim_end(); ++ait1) {
198 
199  string dimName = dapArray->dimension_name(ait1);
200 
201  if (dimName.compare(nodeDimensionName) == 0) { // are the names the same?
202  BESDEBUG("ugrid",
203  "TwoDMeshTopology::setNodeCoordinateDimension() - Found dimension name matching nodeDimensionName '"<< nodeDimensionName << "'" << endl);
204  int size = dapArray->dimension_size(ait1, true);
205  if (size == nodeCount) { // are they the same size?
206  BESDEBUG("ugrid",
207  "TwoDMeshTopology::setNodeCoordinateDimension() - Dimension sizes match (" << libdap::long_to_string(nodeCount) << ") - DONE" << endl);
208  // Yay! We found the node coordinate dimension
209  mdv->setLocationCoordinateDimension(ait1);
210  return;
211  }
212  }
213  }
214 
215  throw Error(
216  "Unable to determine the node coordinate dimension of the range variable '" + mdv->getName()
217  + "' The node dimension is named '" + nodeDimensionName + "' with size "
218  + libdap::long_to_string(nodeCount));
219 
220 }
221 
222 void TwoDMeshTopology::setFaceCoordinateDimension(MeshDataVariable *mdv)
223 {
224 
225  libdap::Array *dapArray = mdv->getDapArray();
226  libdap::Array::Dim_iter ait1;
227 
228  for (ait1 = dapArray->dim_begin(); ait1 != dapArray->dim_end(); ++ait1) {
229  string dimName = dapArray->dimension_name(ait1);
230 
231  if (dimName.compare(faceDimensionName) == 0) { // are the names the same?
232  int size = dapArray->dimension_size(ait1, true);
233  if (size == faceCount) { // are they the same size?
234  // Yay! We found the coordinate dimension
235  mdv->setLocationCoordinateDimension(ait1);
236  return;
237  }
238  }
239 
240  }
241  throw Error(
242  "Unable to determine the face coordinate dimension of the range variable '" + mdv->getName()
243  + "' The face coordinate dimension is named '" + faceDimensionName + "' with size "
244  + libdap::long_to_string(faceCount));
245 
246 }
247 
248 void TwoDMeshTopology::setLocationCoordinateDimension(MeshDataVariable *mdv)
249 {
250 
251  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - BEGIN" << endl);
252 
253  // libdap::Array *dapArray = mdv->getDapArray();
254 
255  string locstr;
256  switch (mdv->getGridLocation()) {
257 
258  case node: {
259  BESDEBUG("ugrid",
260  "TwoDMeshTopology::setLocationCoordinateDimension() - Checking Node variable '"<< mdv->getName() << "'" << endl);
261  // Locate and set the MDV's node coordinate dimension.
262  setNodeCoordinateDimension(mdv);
263  locstr = "node";
264  }
265  break;
266 
267  case face: {
268  BESDEBUG("ugrid",
269  "TwoDMeshTopology::setLocationCoordinateDimension() - Checking Face variable '"<< mdv->getName() << "'" << endl);
270  // Locate and set the MDV's face coordinate dimension.
271  setFaceCoordinateDimension(mdv);
272  locstr = "face";
273  }
274  break;
275 
276  case edge:
277  default: {
278  string msg = "TwoDMeshTopology::setLocationCoordinateDimension() - Unknown/Unsupported location value '"
279  + libdap::long_to_string(mdv->getGridLocation()) + "'";
280  BESDEBUG("ugrid", msg << endl);
281  throw Error(msg);
282  }
283  break;
284  }
285  BESDEBUG("ugrid",
286  "TwoDMeshTopology::setLocationCoordinateDimension() - MeshDataVariable '"<< mdv->getName() << "' is a "<< locstr <<" variable," << endl);
287  BESDEBUG("ugrid",
288  "TwoDMeshTopology::setLocationCoordinateDimension() - MeshDataVariable '"<< mdv->getName() << "' location coordinate dimension is '" << /*dapArray*/mdv->getDapArray()->dimension_name(mdv->getLocationCoordinateDimension()) << "'" << endl);
289 
290  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - DONE" << endl);
291 
292 }
293 
301 void TwoDMeshTopology::ingestFaceNodeConnectivityArray(libdap::BaseType *meshTopology, libdap::DDS *dds)
302 {
303 
304  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Locating FNCA" << endl);
305 
306  string face_node_connectivity_var_name;
307  AttrTable at = meshTopology->get_attr_table();
308 
309  AttrTable::Attr_iter iter_fnc = at.simple_find(UGRID_FACE_NODE_CONNECTIVITY);
310  if (iter_fnc != at.attr_end()) {
311  face_node_connectivity_var_name = at.get_attr(iter_fnc, 0);
312  }
313  else {
314  throw Error(
315  "Could not locate the " UGRID_FACE_NODE_CONNECTIVITY " attribute in the " UGRID_MESH_TOPOLOGY " variable! "
316  "The mesh_topology variable is named " + meshTopology->name());
317  }
318  BESDEBUG("ugrid",
319  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Located the '" << UGRID_FACE_NODE_CONNECTIVITY << "' attribute." << endl);
320 
321  // Find the variable using the name
322 
323  BaseType *btp = dds->var(face_node_connectivity_var_name);
324 
325  if (btp == 0)
326  throw Error(
327  "Could not locate the " UGRID_FACE_NODE_CONNECTIVITY " variable named '" + face_node_connectivity_var_name
328  + "'! " + "The mesh_topology variable is named " + meshTopology->name());
329 
330  // Is it an array?
331  libdap::Array *fncArray = dynamic_cast<libdap::Array*>(btp);
332  if (fncArray == 0) {
333  throw Error(malformed_expr,
334  "Face Node Connectivity variable '" + face_node_connectivity_var_name
335  + "' is not an Array type. It's an instance of " + btp->type_name());
336  }
337 
338  BESDEBUG("ugrid",
339  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Located FNC Array '" << fncArray->name() << "'." << endl);
340 
341  // It's got to have exactly 2 dimensions - [max#nodes_per_face][#faces]
342  int numDims = fncArray->dimensions(true);
343  if (numDims != 2) {
344  throw Error(malformed_expr,
345  "Face Node Connectivity variable '" + face_node_connectivity_var_name
346  + "' Must have two (2) dimensions. It has " + libdap::long_to_string(numDims));
347  }
348 
349  BESDEBUG("ugrid",
350  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - FNC Array '" << fncArray->name() << "' has two (2) dimensions." << endl);
351 
352  // We just need to have both dimensions handy, so get the first and the second
353  libdap::Array::Dim_iter firstDim = fncArray->dim_begin();
354  libdap::Array::Dim_iter secondDim = fncArray->dim_begin(); // same as the first for a moment...
355  secondDim++; // now it's second!
356 
357  if (faceDimensionName.empty()) {
358  // By now we know it only has two dimensions, but since there is no promise that they'll be in a particular order
359  // we punt: We'll assume that smallest of the two is in fact the nodes per face and the larger the face index dimensions.
360  int sizeFirst = fncArray->dimension_size(firstDim, true);
361  int sizeSecond = fncArray->dimension_size(secondDim, true);
362  BESDEBUG("ugrid",
363  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - sizeFirst: "<< libdap::long_to_string(sizeFirst) << endl);
364  BESDEBUG("ugrid",
365  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - sizeSecond: "<< libdap::long_to_string(sizeSecond) << endl);
366 
367  if (sizeFirst < sizeSecond) {
368  BESDEBUG("ugrid",
369  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - FNC Array first dimension is smaller than second." << endl);
370  fncNodesDim = firstDim;
371  fncFacesDim = secondDim;
372  }
373  else {
374  BESDEBUG("ugrid",
375  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - FNC Array second dimension is smaller than first." << endl);
376  fncNodesDim = secondDim;
377  fncFacesDim = firstDim;
378  }
379 
380  BESDEBUG("ugrid",
381  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - fncFacesDim meshVarName: '" << fncArray->dimension_name(fncFacesDim) << "'" << endl);
382  BESDEBUG("ugrid",
383  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - fncFacesDim size: '" << fncArray->dimension_size(fncFacesDim,true) << "'" << endl);
384 
385  faceDimensionName = fncArray->dimension_name(fncFacesDim);
386  }
387  else {
388  // There is already a faceDimensionName defined - possibly from loading face coordinate variables.
389 
390  // Does it match the name of the first or second dimensions of the fncArray? It better!
391  if (faceDimensionName.compare(fncArray->dimension_name(firstDim)) == 0) {
392  fncNodesDim = secondDim;
393  fncFacesDim = firstDim;
394  }
395  else if (faceDimensionName.compare(fncArray->dimension_name(secondDim)) == 0) {
396  fncNodesDim = firstDim;
397  fncFacesDim = secondDim;
398  }
399  else {
400  string msg = "The face coordinate dimension of the Face Node Connectivity variable '"
401  + face_node_connectivity_var_name + "' Has dimension name.'" + fncArray->dimension_name(fncFacesDim)
402  + "' which does not match the existing face coordinate dimension meshVarName '" + faceDimensionName
403  + "'";
404  BESDEBUG("ugrid", msg << endl);
405  throw Error(msg);
406  }
407  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face dimension names match." << endl);
408 
409  }
410 
411  BESDEBUG("ugrid",
412  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face dimension name: '" << faceDimensionName << "'" << endl);
413 
414  // Check to see if faceCount is initialized and do so if needed
415  if (faceCount == 0) {
416  faceCount = fncArray->dimension_size(fncFacesDim, true);
417  BESDEBUG("ugrid",
418  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face count: "<< libdap::long_to_string(faceCount) << endl);
419  }
420  else {
421  // Make sure the face counts match.
422 
423  BESDEBUG("ugrid",
424  "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face count: "<< libdap::long_to_string(faceCount) << endl);
425  if (faceCount != fncArray->dimension_size(fncFacesDim, true)) {
426  string msg = "The faces dimension of the Face Node Connectivity variable '"
427  + face_node_connectivity_var_name + "' Has size "
428  + libdap::long_to_string(fncArray->dimension_size(fncFacesDim, true))
429  + " which does not match the existing face count of " + libdap::long_to_string(faceCount);
430  BESDEBUG("ugrid", msg << endl);
431  throw Error(msg);
432  }
433  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face counts match!" << endl);
434  }
435 
436  faceNodeConnectivityArray = fncArray;
437 
438  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityArray() - Got FCNA '"+fncArray->name()+"'" << endl);
439 
440 }
446 void TwoDMeshTopology::ingestFaceCoordinateArrays(libdap::BaseType *meshTopology, libdap::DDS *dds)
447 {
448  BESDEBUG("ugrid",
449  "TwoDMeshTopology::ingestFaceCoordinateArrays() - BEGIN Gathering face coordinate arrays..." << endl);
450 
451  if (faceCoordinateArrays == 0) faceCoordinateArrays = new vector<libdap::Array *>();
452 
453  faceCoordinateArrays->clear();
454 
455  string face_coordinates;
456  AttrTable at = meshTopology->get_attr_table();
457 
458  AttrTable::Attr_iter iter_nodeCoors = at.simple_find(UGRID_FACE_COORDINATES);
459  if (iter_nodeCoors != at.attr_end()) {
460  face_coordinates = at.get_attr(iter_nodeCoors, 0);
461 
462  BESDEBUG("ugrid",
463  "TwoDMeshTopology::ingestFaceCoordinateArrays() - Located '"<< UGRID_FACE_COORDINATES << "' attribute." << endl);
464 
465  // Split the node_coordinates string up on spaces
466  vector<string> faceCoordinateNames = split(face_coordinates, ' ');
467 
468  // Find each variable in the resulting list
469  vector<string>::iterator coorName_it;
470  for (coorName_it = faceCoordinateNames.begin(); coorName_it != faceCoordinateNames.end(); ++coorName_it) {
471  string faceCoordinateName = *coorName_it;
472 
473  BESDEBUG("ugrid",
474  "TwoDMeshTopology::ingestFaceCoordinateArrays() - Processing face coordinate '"<< faceCoordinateName << "'." << endl);
475 
476  //Now that we have the name of the coordinate variable get it from the DDS!!
477  BaseType *btp = dds->var(faceCoordinateName);
478  if (btp == 0)
479  throw Error(
480  "Could not locate the " UGRID_FACE_COORDINATES " variable named '" + faceCoordinateName + "'! "
481  + "The mesh_topology variable is named " + meshTopology->name());
482 
483  libdap::Array *newFaceCoordArray = dynamic_cast<libdap::Array*>(btp);
484  if (newFaceCoordArray == 0) {
485  throw Error(malformed_expr,
486  "Face coordinate variable '" + faceCoordinateName + "' is not an Array type. It's an instance of "
487  + btp->type_name());
488  }
489 
490  BESDEBUG("ugrid",
491  "TwoDMeshTopology::ingestFaceCoordinateArrays() - Found face coordinate array '"<< faceCoordinateName << "'." << endl);
492 
493  // Coordinate arrays MUST be single dimensioned.
494  if (newFaceCoordArray->dimensions(true) != 1) {
495  throw Error(malformed_expr,
496  "Face coordinate variable '" + faceCoordinateName
497  + "' has more than one dimension. That's just not allowed. It has "
498  + long_to_string(newFaceCoordArray->dimensions(true)) + " dimensions.");
499  }
500  BESDEBUG("ugrid",
501  "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face coordinate array '"<< faceCoordinateName << "' has a single dimension." << endl);
502 
503  // Make sure this node coordinate variable has the same size and meshVarName as all the others on the list - error if not true.
504  string dimName = newFaceCoordArray->dimension_name(newFaceCoordArray->dim_begin());
505  int dimSize = newFaceCoordArray->dimension_size(newFaceCoordArray->dim_begin(), true);
506 
507  BESDEBUG("ugrid",
508  "TwoDMeshTopology::ingestFaceCoordinateArrays() - dimName: '"<< dimName << "' dimSize: " << libdap::long_to_string(dimSize) << endl);
509 
510  if (faceDimensionName.empty()) {
511  faceDimensionName = dimName;
512  }
513  BESDEBUG("ugrid",
514  "TwoDMeshTopology::ingestFaceCoordinateArrays() - faceDimensionName: '"<< faceDimensionName << "' " << endl);
515 
516  if (faceDimensionName.compare(dimName) != 0) {
517  throw Error(
518  "The face coordinate array '" + faceCoordinateName + "' has the named dimension '" + dimName
519  + "' which differs from the expected dimension meshVarName '" + faceDimensionName
520  + "'. The mesh_topology variable is named " + meshTopology->name());
521  }
522  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face dimension names match." << endl);
523 
524  if (faceCount == 0) {
525  faceCount = dimSize;
526  }
527  BESDEBUG("ugrid",
528  "TwoDMeshTopology::ingestFaceCoordinateArrays() - faceCount: "<< libdap::long_to_string(faceCount) << endl);
529 
530  if (faceCount != dimSize) {
531  throw Error(
532  "The face coordinate array '" + faceCoordinateName + "' has a dimension size of "
533  + libdap::long_to_string(dimSize) + " which differs from the the expected size of "
534  + libdap::long_to_string(faceCount) + " The mesh_topology variable is named "
535  + meshTopology->name());
536  }
537  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face counts match." << endl);
538 
539  // Add variable to faceCoordinateArrays.
540  faceCoordinateArrays->push_back(newFaceCoordArray);
541  BESDEBUG("ugrid",
542  "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face coordinate array '"<< faceCoordinateName << "' ingested." << endl);
543  }
544  BESDEBUG("ugrid",
545  "TwoDMeshTopology::ingestFaceCoordinateArrays() - Located "<< libdap::long_to_string(faceCoordinateArrays->size()) << " face coordinate arrays." << endl);
546 
547  }
548  else {
549  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - No Face Coordinates Found." << endl);
550  }
551 
552  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - DONE" << endl);
553 
554 }
555 
561 void TwoDMeshTopology::ingestNodeCoordinateArrays(libdap::BaseType *meshTopology, libdap::DDS *dds)
562 {
563  BESDEBUG("ugrid",
564  "TwoDMeshTopology::ingestNodeCoordinateArrays() - BEGIN Gathering node coordinate arrays..." << endl);
565 
566  string node_coordinates;
567  AttrTable at = meshTopology->get_attr_table();
568 
569  AttrTable::Attr_iter iter_nodeCoors = at.simple_find(UGRID_NODE_COORDINATES);
570  if (iter_nodeCoors != at.attr_end()) {
571  node_coordinates = at.get_attr(iter_nodeCoors, 0);
572  }
573  else {
574  throw Error("Could not locate the " UGRID_NODE_COORDINATES " attribute in the " UGRID_MESH_TOPOLOGY
575  " variable! The mesh_topology variable is named " + meshTopology->name());
576  }
577  BESDEBUG("ugrid",
578  "TwoDMeshTopology::ingestNodeCoordinateArrays() - Located '"<< UGRID_NODE_COORDINATES << "' attribute." << endl);
579 
580  if (nodeCoordinateArrays == 0) nodeCoordinateArrays = new vector<libdap::Array *>();
581 
582  nodeCoordinateArrays->clear();
583 
584  // Split the node_coordinates string up on spaces
585  // TODO make this work on situations where multiple spaces in the node_coorindates string doesn't hose the split()
586  vector<string> nodeCoordinateNames = split(node_coordinates, ' ');
587 
588  // Find each variable in the resulting list
589  vector<string>::iterator coorName_it;
590  for (coorName_it = nodeCoordinateNames.begin(); coorName_it != nodeCoordinateNames.end(); ++coorName_it) {
591  string nodeCoordinateName = *coorName_it;
592 
593  BESDEBUG("ugrid",
594  "TwoDMeshTopology::ingestNodeCoordinateArrays() - Processing node coordinate '"<< nodeCoordinateName << "'." << endl);
595 
596  //Now that we have the name of the coordinate variable get it from the DDS!!
597  BaseType *btp = dds->var(nodeCoordinateName);
598  if (btp == 0)
599  throw Error(
600  "Could not locate the " UGRID_NODE_COORDINATES " variable named '" + nodeCoordinateName + "'! "
601  + "The mesh_topology variable is named " + meshTopology->name());
602 
603  libdap::Array *newNodeCoordArray = dynamic_cast<libdap::Array*>(btp);
604  if (newNodeCoordArray == 0) {
605  throw Error(malformed_expr,
606  "Node coordinate variable '" + nodeCoordinateName + "' is not an Array type. It's an instance of "
607  + btp->type_name());
608  }
609  BESDEBUG("ugrid",
610  "TwoDMeshTopology::ingestNodeCoordinateArrays() - Found node coordinate array '"<< nodeCoordinateName << "'." << endl);
611 
612  // Coordinate arrays MUST be single dimensioned.
613  if (newNodeCoordArray->dimensions(true) != 1) {
614  throw Error(malformed_expr,
615  "Node coordinate variable '" + nodeCoordinateName
616  + "' has more than one dimension. That's just not allowed. It has "
617  + long_to_string(newNodeCoordArray->dimensions(true)) + " dimensions.");
618  }
619  BESDEBUG("ugrid",
620  "TwoDMeshTopology::ingestNodeCoordinateArrays() - Node coordinate array '"<< nodeCoordinateName << "' has a single dimension." << endl);
621 
622  // Make sure this node coordinate variable has the same size and meshVarName as all the others on the list - error if not true.
623  string dimName = newNodeCoordArray->dimension_name(newNodeCoordArray->dim_begin());
624  int dimSize = newNodeCoordArray->dimension_size(newNodeCoordArray->dim_begin(), true);
625 
626  BESDEBUG("ugrid",
627  "TwoDMeshTopology::ingestNodeCoordinateArrays() - dimName: '"<< dimName << "' dimSize: " << libdap::long_to_string(dimSize) << endl);
628 
629  if (nodeDimensionName.empty()) {
630  nodeDimensionName = dimName;
631  }
632  BESDEBUG("ugrid",
633  "TwoDMeshTopology::ingestNodeCoordinateArrays() - nodeDimensionName: '"<< nodeDimensionName << "' " << endl);
634  if (nodeDimensionName.compare(dimName) != 0) {
635  throw Error(
636  "The node coordinate array '" + nodeCoordinateName + "' has the named dimension '" + dimName
637  + "' which differs from the expected dimension meshVarName '" + nodeDimensionName
638  + "'. The mesh_topology variable is named " + meshTopology->name());
639  }
640  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Node dimension names match." << endl);
641 
642  if (nodeCount == 0) {
643  nodeCount = dimSize;
644  }
645  BESDEBUG("ugrid",
646  "TwoDMeshTopology::ingestNodeCoordinateArrays() - nodeCount: "<< libdap::long_to_string(nodeCount) << endl);
647 
648  if (nodeCount != dimSize) {
649  throw Error(
650  "The node coordinate array '" + nodeCoordinateName + "' has a dimension size of "
651  + libdap::long_to_string(dimSize) + " which differs from the the expected size of "
652  + libdap::long_to_string(nodeCount) + " The mesh_topology variable is named "
653  + meshTopology->name());
654  }
655  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Node counts match." << endl);
656 
657  // Add variable to nodeCoordinateArrays vector.
658  nodeCoordinateArrays->push_back(newNodeCoordArray);
659  BESDEBUG("ugrid",
660  "TwoDMeshTopology::ingestNodeCoordinateArrays() - Coordinate array '"<< nodeCoordinateName << "' ingested." << endl);
661 
662  }
663 
664  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - DONE" << endl);
665 
666 }
667 
668 void TwoDMeshTopology::buildBasicGfTopology()
669 {
670 
671  BESDEBUG("ugrid",
672  "TwoDMeshTopology::buildBasicGfTopology() - Building GridFields objects for mesh_topology variable "<< getMeshVariable()->name() << endl);
673  // Start building the Grid for the GridField operation.
674  BESDEBUG("ugrid",
675  "TwoDMeshTopology::buildGridFieldsTopology() - Constructing new GF::Grid for "<< meshVarName() << endl);
676  gridTopology = new GF::Grid(meshVarName());
677 
678  // 1) Make the implicit nodes - same size as the node coordinate arrays
679  BESDEBUG("ugrid",
680  "TwoDMeshTopology::buildGridFieldsTopology() - Building and adding implicit range Nodes to the GF::Grid" << endl);
681  GF::AbstractCellArray *nodes = new GF::Implicit0Cells(nodeCount);
682  // Attach the implicit nodes to the grid at rank 0
683  gridTopology->setKCells(nodes, node);
684 
685  // @TODO Do I need to add implicit k-cells for faces (rank 2) if I plan to add range data on faces later?
686  // Apparently not...
687 
688  // Attach the Mesh to the grid.
689  // Get the face node connectivity cells (i think these correspond to the GridFields K cells of Rank 2)
690  // FIXME Read this array once! It is read again below..
691  BESDEBUG("ugrid",
692  "TwoDMeshTopology::buildGridFieldsTopology() - Building face node connectivity Cell array from the DAP version" << endl);
693  GF::CellArray *faceNodeConnectivityCells = getFaceNodeConnectivityCells();
694 
695  // Attach the Mesh to the grid at rank 2
696  // This 2 stands for rank 2, or faces.
697  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Attaching Cell array to GF::Grid" << endl);
698  gridTopology->setKCells(faceNodeConnectivityCells, face);
699 
700  // The Grid is complete. Now we make a GridField from the Grid
701  BESDEBUG("ugrid",
702  "TwoDMeshTopology::buildGridFieldsTopology() - Construct new GF::GridField from GF::Grid" << endl);
703  d_inputGridField = new GF::GridField(gridTopology);
704  // TODO Question for Bill: Can we delete the GF::Grid (tdmt->gridTopology) here?
705 
706  // We read and add the coordinate data (using GridField->addAttribute()) to the GridField at
707  // grid dimension/rank/dimension 0 (a.k.a. node)
708  vector<libdap::Array *>::iterator ncit;
709  for (ncit = nodeCoordinateArrays->begin(); ncit != nodeCoordinateArrays->end(); ++ncit) {
710  libdap::Array *nca = *ncit;
711  BESDEBUG("ugrid",
712  "TwoDMeshTopology::buildGridFieldsTopology() - Adding node coordinate "<< nca->name() << " to GF::GridField at rank 0" << endl);
713  GF::Array *gfa = extractGridFieldArray(nca, sharedIntArrays, sharedFloatArrays);
714  gfArrays.push_back(gfa);
715  d_inputGridField->AddAttribute(node, gfa);
716  }
717 
718  // We read and add the coordinate data (using GridField->addAttribute() to the GridField at
719  // grid dimension/rank/dimension 0 (a.k.a. node)
720  for (ncit = faceCoordinateArrays->begin(); ncit != faceCoordinateArrays->end(); ++ncit) {
721  libdap::Array *fca = *ncit;
722  BESDEBUG("ugrid",
723  "TwoDMeshTopology::buildGridFieldsTopology() - Adding face coordinate "<< fca->name() << " to GF::GridField at rank " << face << endl);
724  GF::Array *gfa = extractGridFieldArray(fca, sharedIntArrays, sharedFloatArrays);
725  gfArrays.push_back(gfa);
726  d_inputGridField->AddAttribute(face, gfa);
727  }
728 }
729 
730 int TwoDMeshTopology::getResultGridSize(locationType dim)
731 {
732  return resultGridField->Size(dim);
733 }
734 
745 GF::Node *TwoDMeshTopology::getFncArrayAsGFCells(libdap::Array *fncVar)
746 {
747  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFCells() - BEGIN" << endl);
748 
749  int nodesPerFace = fncVar->dimension_size(fncNodesDim, true);
750  int faceCount = fncVar->dimension_size(fncFacesDim, true);
751 
752  GF::Node *cells = 0; // return the result in this var
753 
754  if (fncVar->dim_begin() == fncNodesDim) {
755  // This dataset/file stores the face-node connectivity array as a
756  // 3xN, but gridfields needs that information in an Nx3; twiddle
757  BESDEBUG("ugrid", "Reorganizing the data from teh DAP FNC Array for GF." << endl);
758 
759  cells = new GF::Node[faceCount * nodesPerFace];
760  GF::Node *temp_nodes = 0;
761  // optimize; extractArray uses an extra copy in this case
762  if (fncVar->type() == dods_int32_c || fncVar->type() == dods_uint32_c) {
763  temp_nodes = new GF::Node[faceCount * nodesPerFace];
764  fncVar->value(temp_nodes);
765  }
766  else {
767  // cover the odd case when the FNC Array is not a int/uint.
768  temp_nodes = ugrid::extractArray<GF::Node>(fncVar);
769  }
770 
771  for (int fIndex = 0; fIndex < faceCount; fIndex++) {
772  for (int nIndex = 0; nIndex < nodesPerFace; nIndex++) {
773  cells[nodesPerFace * fIndex + nIndex] = *(temp_nodes + (fIndex + (faceCount * nIndex)));
774  }
775  }
776 
777  delete[] temp_nodes;
778  }
779  else {
780  // This dataset/file stores the face-node connectivity array as an Nx3 which is what
781  // gridfields expects. We can use the libdap::BaseType::value() method to slurp
782  // up the stuff.
783 
784  if (fncVar->type() == dods_int32_c || fncVar->type() == dods_uint32_c) {
785  cells = new GF::Node[faceCount * nodesPerFace];
786  fncVar->value(cells);
787  }
788  else {
789  // cover the odd case when the FNC Array is not a int/uint.
790  cells = ugrid::extractArray<GF::Node>(fncVar);
791  }
792  }
793 
794  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFCells() - DONE" << endl);
795  return cells;
796 }
797 
802 int TwoDMeshTopology::getStartIndex(libdap::Array *array)
803 {
804  AttrTable &at = array->get_attr_table();
805  AttrTable::Attr_iter start_index_iter = at.simple_find(UGRID_START_INDEX);
806  if (start_index_iter != at.attr_end()) {
807  BESDEBUG("ugrid", "TwoDMeshTopology::getStartIndex() - Found the "<< UGRID_START_INDEX <<" attribute." << endl);
808  AttrTable::entry *start_index_entry = *start_index_iter;
809  if (start_index_entry->attr->size() == 1) {
810  string val = (*start_index_entry->attr)[0];
811  BESDEBUG("TwoDMeshTopology::getStartIndex", "value: " << val << endl);
812  stringstream buffer(val);
813  // what happens if string cannot be converted to an integer?
814  int start_index;
815  buffer >> start_index;
816  return start_index;
817  }
818  else {
819  throw Error(malformed_expr,
820  "Index origin attribute exists, but either no value supplied, or more than one value supplied.");
821  }
822  }
823 
824  return 0;
825 }
826 
831 GF::CellArray *TwoDMeshTopology::getFaceNodeConnectivityCells()
832 {
833  BESDEBUG("ugrid",
834  "TwoDMeshTopology::getFaceNodeConnectivityCells() - Building face node connectivity Cell array from the Array '" << faceNodeConnectivityArray->name() << "'" << endl);
835 
836  int nodesPerFace = faceNodeConnectivityArray->dimension_size(fncNodesDim);
837  int total_size = nodesPerFace * faceCount;
838 
839  BESDEBUG("ugrid",
840  "TwoDMeshTopology::getFaceNodeConnectivityCells() - Converting FNCArray to GF::Node array." << endl);
841  fncCellArray = getFncArrayAsGFCells(faceNodeConnectivityArray);
842 
843  // adjust for the start_index (cardinal or ordinal array access)
844  int startIndex = getStartIndex(faceNodeConnectivityArray);
845  if (startIndex != 0) {
846  BESDEBUG("ugrid",
847  "TwoDMeshTopology::getFaceNodeConnectivityCells() - Applying startIndex to GF::Node array." << endl);
848  for (int j = 0; j < total_size; j++) {
849  fncCellArray[j] -= startIndex;
850  }
851  }
852  // Create the cell array
853  GF::CellArray *rankTwoCells = new GF::CellArray(fncCellArray, faceCount, nodesPerFace);
854 
855  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityCells() - DONE" << endl);
856  return rankTwoCells;
857 
858 }
859 
860 void TwoDMeshTopology::applyRestrictOperator(locationType loc, string filterExpression)
861 {
862 
863  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - BEGIN" << endl);
864 
865  // I think this function could be done with just the following single line:
866  // resultGridField = GF::RefRestrictOp::Restrict(filterExpression,loc,d_inputGridField);
867 
868  // Build the restriction operator
869  BESDEBUG("ugrid",
870  "TwoDMeshTopology::applyRestrictOperator() - Constructing new GF::RestrictOp using user "<< "supplied 'dimension' value and filter expression combined with the GF:GridField " << endl);
871  GF::RestrictOp op = GF::RestrictOp(filterExpression, loc, d_inputGridField);
872 
873  // Apply the operator and get the result;
874  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - Applying GridField operator." << endl);
875  GF::GridField *resultGF = op.getResult();
876 
877  resultGridField = resultGF;
878  BESDEBUG("ugrid",
879  "TwoDMeshTopology::applyRestrictOperator() - GridField operator applied and result obtained." << endl);
880 
881  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - END" << endl);
882 }
883 
884 void TwoDMeshTopology::convertResultGridFieldStructureToDapObjects(vector<BaseType *> *results)
885 {
886  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - BEGIN" << endl);
887 
888  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Normalizing Grid." << endl);
889  resultGridField->GetGrid()->normalize();
890 
891  BESDEBUG("ugrid",
892  "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - resultGridField->MaxRank(): "<< resultGridField->MaxRank() << endl);
893 
894  if (resultGridField->MaxRank() < 0) {
895  throw BESError("Oops! The ugrid constraint expression resulted in an empty response.", BES_SYNTAX_USER_ERROR,
896  __FILE__, __LINE__);
897  }
898 
899  // Add the node coordinate arrays to the results.
900  BESDEBUG("ugrid",
901  "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Converting the node coordinate arrays to DAP arrays." << endl);
902  vector<libdap::Array *>::iterator it;
903  for (it = nodeCoordinateArrays->begin(); it != nodeCoordinateArrays->end(); ++it) {
904  libdap::Array *sourceCoordinateArray = *it;
905  libdap::Array *resultCoordinateArray = getGFAttributeAsDapArray(sourceCoordinateArray, node, resultGridField);
906  results->push_back(resultCoordinateArray);
907  }
908 
909 #if 1
910  // Add the face coordinate arrays to the results.
911  BESDEBUG("ugrid",
912  "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Converting the face coordinate arrays to DAP arrays." << endl);
913  for (it = faceCoordinateArrays->begin(); it != faceCoordinateArrays->end(); ++it) {
914  libdap::Array *sourceCoordinateArray = *it;
915  libdap::Array *resultCoordinateArray = getGFAttributeAsDapArray(sourceCoordinateArray, face, resultGridField);
916  results->push_back(resultCoordinateArray);
917  }
918 #endif
919 
920  // Add the new face node connectivity array - make sure it has the same attributes as the original.
921  BESDEBUG("ugrid",
922  "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Adding the new face node connectivity array to the response." << endl);
923  libdap::Array *resultFaceNodeConnectivityDapArray = getGridFieldCellArrayAsDapArray(resultGridField,
924  faceNodeConnectivityArray);
925  results->push_back(resultFaceNodeConnectivityDapArray);
926 
927  results->push_back(getMeshVariable()->ptr_duplicate());
928 
929  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - END" << endl);
930 }
931 
937 libdap::Array *TwoDMeshTopology::getGridFieldCellArrayAsDapArray(GF::GridField *resultGridField,
938  libdap::Array *sourceFcnArray)
939 {
940  BESDEBUG("ugrid", "TwoDMeshTopology::getGridFieldCellArrayAsDapArray() - BEGIN" << endl);
941 
942  // Get the rank 2 k-cells from the GridField object.
943  GF::CellArray* gfCellArray = (GF::CellArray*) (resultGridField->GetGrid()->getKCells(2));
944 
945  // This is a vector of size N holding vectors of size 3
946  vector<vector<int> > nodes2 = gfCellArray->makeArrayInts();
947 
948  libdap::Array *resultFncDapArray = new libdap::Array(sourceFcnArray->name(), new Int32(sourceFcnArray->name()));
949 
950  // Is the sourceFcnArray a Nx3 (follows the ugrid 0.9 spec) or 3xN - both
951  // commonly appear. Make the resultFncDapArray match the source's organization
952  // modulo that 'N' is a different value now given that the ugrid has been
953  // subset. jhrg 4/17/15
954  bool follows_ugrid_09_spec = true;
955  libdap::Array::Dim_iter di = sourceFcnArray->dim_begin();
956  if (di->size == 3) {
957  follows_ugrid_09_spec = false;
958 
959  resultFncDapArray->append_dim(3, di->name);
960  ++di;
961  resultFncDapArray->append_dim(nodes2.size(), di->name);
962  }
963  else {
964  resultFncDapArray->append_dim(nodes2.size(), di->name);
965  ++di;
966  resultFncDapArray->append_dim(3, di->name);
967  }
968 
969  // Copy the attributes of the template array to our new array.
970  resultFncDapArray->set_attr_table(sourceFcnArray->get_attr_table());
971  resultFncDapArray->reserve_value_capacity(3 * nodes2.size());
972 
973  // adjust for the start_index (cardinal or ordinal array access)
974  int startIndex = getStartIndex(sourceFcnArray);
975 
976  // Now transfer the index information returned by Gridfields to the newly
977  // made DAP array. Since GF returns those indexes as an Nx3 using zero-based
978  // indexing, we may have to transform the values.
979  if (!follows_ugrid_09_spec) {
980  // build a new set of values and set them as the value of the result fnc array.
981  vector<dods_int32> node_data(3 * nodes2.size());
982  vector<dods_int32>::iterator ndi = node_data.begin();
983 
984  if (startIndex == 0) {
985  for (unsigned int i = 0; i < 3; ++i) {
986  for (unsigned int n = 0; n < nodes2.size(); ++n) {
987  *ndi++ = nodes2[n][i];
988  }
989  }
990  }
991  else {
992  for (unsigned int i = 0; i < 3; ++i) {
993  for (unsigned int n = 0; n < nodes2.size(); ++n) {
994  *ndi++ = nodes2[n][i] + startIndex;
995  }
996  }
997  }
998 
999  resultFncDapArray->set_value(node_data, node_data.size());
1000  }
1001  else {
1002  // build a new set of values and set them as the value of the result fnc array.
1003  vector<dods_int32> node_data(nodes2.size() * 3);
1004  vector<dods_int32>::iterator ndi = node_data.begin();
1005 
1006  if (startIndex == 0) {
1007  for (unsigned int n = 0; n < nodes2.size(); ++n) {
1008  for (unsigned int i = 0; i < 3; ++i) {
1009  *ndi++ = nodes2[n][i];
1010  }
1011  }
1012  }
1013  else {
1014  for (unsigned int n = 0; n < nodes2.size(); ++n) {
1015  for (unsigned int i = 0; i < 3; ++i) {
1016  *ndi++ = nodes2[n][i] + startIndex;
1017  }
1018  }
1019  }
1020 
1021  resultFncDapArray->set_value(node_data, node_data.size());
1022  }
1023 
1024  BESDEBUG("ugrid", "TwoDMeshTopology::getGridFieldCellArrayAsDapArray() - DONE" << endl);
1025 
1026  return resultFncDapArray;
1027 }
1028 
1039 static string copySizeOneDimensions(libdap::Array *sourceArray, libdap::Array *dapArray)
1040 {
1041  for (libdap::Array::Dim_iter srcArrIt = sourceArray->dim_begin(); srcArrIt != sourceArray->dim_end(); ++srcArrIt) {
1042  // Get the original dimension size
1043  int dimSize = sourceArray->dimension_size(srcArrIt, true);
1044  string dimName = sourceArray->dimension_name(srcArrIt);
1045 
1046  // Preserve single dimensions
1047  if (dimSize == 1)
1048  dapArray->append_dim(dimSize, dimName);
1049  else
1050  return dimName;
1051  }
1052 
1053  return "";
1054 }
1055 
1056 #if 0
1057 // the original version jhrg 4/15/15
1058 
1059 static libdap::Array::Dim_iter copySizeOneDimensions(libdap::Array *sourceArray, libdap::Array *dapArray)
1060 {
1061  libdap::Array::Dim_iter dataDimension;
1062  for (libdap::Array::Dim_iter srcArrIt = sourceArray->dim_begin(); srcArrIt != sourceArray->dim_end(); ++srcArrIt) {
1063 
1064  // Get the original dimension size
1065  int dimSize = sourceArray->dimension_size(srcArrIt, true);
1066  string dimName = sourceArray->dimension_name(srcArrIt);
1067 
1068  // Preserve single dimensions
1069  if (dimSize == 1) {
1070  BESDEBUG("ugrid", "TwoDMeshTopology::copySizeOneDimensions() - Adding size one dimension '"<< dimName << "' from source array into result." << endl);
1071  dapArray->append_dim(dimSize, dimName);
1072  }
1073  else {
1074  BESDEBUG("ugrid", "TwoDMeshTopology::copySizeOneDimensions() - Located data dimension '"<< dimName << "' in source array." << endl);
1075  dataDimension = srcArrIt;
1076  }
1077  }
1078 
1079  BESDEBUG("ugrid", "TwoDMeshTopology::copySizeOneDimensions() - Returning dimension iterator pointing to '"<< sourceArray->dimension_name(dataDimension) << "'." << endl);
1080  return dataDimension;
1081 }
1082 #endif
1083 
1088 libdap::Array *TwoDMeshTopology::getGFAttributeAsDapArray(libdap::Array *templateArray, locationType rank,
1089  GF::GridField *resultGridField)
1090 {
1091  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - BEGIN" << endl);
1092 
1093  // The result variable is assumed to be bound to the GridField at 'rank'
1094  // Try to get the Attribute from 'rank' with the same name as the source array
1095  BESDEBUG("ugrid",
1096  "TwoDMeshTopology::getGFAttributeAsDapArray() - Retrieving rank "<< rank << " GF::GridField Attribute '" << templateArray->name() << "'" << endl);
1097  GF::Array* gfa = 0;
1098  gfa = resultGridField->GetAttribute(rank, templateArray->name());
1099 
1100  libdap::Array *dapArray;
1101  BaseType *templateVar = templateArray->var();
1102  string dimName;
1103 
1104  switch (templateVar->type()) {
1105  case dods_byte_c:
1106  case dods_uint16_c:
1107  case dods_int16_c:
1108  case dods_uint32_c:
1109  case dods_int32_c: {
1110  // Get the data
1111  vector<dods_int32> GF_ints = gfa->makeArray();
1112 
1113  // Make a DAP array to put the data into.
1114  dapArray = new libdap::Array(templateArray->name(), new libdap::Int32(templateVar->name()) /*&tmpltVar*/);
1115 
1116  // copy the dimensions whose size is "1" from the source array to the result.
1117  dimName = copySizeOneDimensions(templateArray, dapArray);
1118 
1119  // Add the result dimension
1120  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Adding dimension " << dimName << endl);
1121 
1122  dapArray->append_dim(GF_ints.size(), dimName);
1123 
1124  // Add the data
1125  dapArray->set_value(GF_ints, GF_ints.size());
1126  break;
1127  }
1128  case dods_float32_c:
1129  case dods_float64_c: {
1130  // Get the data
1131  vector<dods_float64> GF_floats = gfa->makeArrayf();
1132 
1133  // Make a DAP array to put the data into.
1134  dapArray = new libdap::Array(templateArray->name(), new libdap::Float64(templateVar->name()) /*&tmpltVar*/);
1135 
1136  // copy the dimensions whose size is "1" from the source array to the result.
1137  dimName = copySizeOneDimensions(templateArray, dapArray);
1138 
1139  // Add the result dimension
1140  dapArray->append_dim(GF_floats.size(), dimName);
1141 
1142  // Add the data
1143  dapArray->set_value(GF_floats, GF_floats.size());
1144  break;
1145  }
1146  default:
1147  throw InternalErr(__FILE__, __LINE__, "Unknown DAP type encountered when converting to gridfields array");
1148  }
1149 
1150  // Copy the source objects attributes.
1151  BESDEBUG("ugrid",
1152  "TwoDMeshTopology::getGFAttributeAsDapArray() - Copying libdap::Attribute's from template array " << templateArray->name() << endl);
1153  dapArray->set_attr_table(templateArray->get_attr_table());
1154 
1155  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - DONE" << endl);
1156 
1157  return dapArray;
1158 }
1159 
1164 void TwoDMeshTopology::getResultGFAttributeValues(string attrName, libdap::Type dapType, locationType rank,
1165  void *target)
1166 {
1167  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - BEGIN" << endl);
1168 
1169  // The result variable is assumed to be bound to the GridField at 'rank'
1170  // Try to get the Attribute from 'rank' with the same name as the source array
1171  BESDEBUG("ugrid",
1172  "TwoDMeshTopology::getResultGFAttributeValues() - Retrieving GF::GridField Attribute '" << attrName << "'" << endl);
1173 
1174  GF::Array *gfa = 0;
1175  if (resultGridField->IsAttribute(rank, attrName)) {
1176  gfa = resultGridField->GetAttribute(rank, attrName);
1177  }
1178  else {
1179  string msg = "Oddly, the requested attribute '" + attrName + "' associated with rank "
1180  + libdap::long_to_string(rank) + " does not appear in the resultGridField object! \n"
1181  + "resultGridField->MaxRank(): " + libdap::long_to_string(resultGridField->MaxRank());
1182 
1183  throw InternalErr(__FILE__, __LINE__, "ERROR - Unable to locate requested GridField attribute. " + msg);
1184  }
1185 
1186  switch (dapType) {
1187  case dods_byte_c:
1188  case dods_uint16_c:
1189  case dods_int16_c:
1190  case dods_uint32_c:
1191  case dods_int32_c: {
1192  // Get the data
1193  BESDEBUG("ugrid",
1194  "TwoDMeshTopology::getResultGFAttributeValues() - Retrieving GF::Array as libdap::Array of libdap::Int32." << endl);
1195  vector<dods_int32> GF_ints = gfa->makeArray();
1196 #if 0
1197  // Removing this reduced the runtime by ~1s for one test case
1198  // (from real 0m6.285s; user 0m5.816s to real 0m5.136s; user 0m4.808s)
1199  // jhrg 4/15/15
1200  stringstream s;
1201  for(unsigned int i=0; i<GF_ints.size(); i++) {
1202  s << "GF_ints[" << i << "]: " << GF_ints[i] << endl;
1203  }
1204  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieved GF_ints: "<< endl << s.str());
1205 #endif
1206 
1207  BESDEBUG("ugrid",
1208  "TwoDMeshTopology::getResultGFAttributeValues() - Copying GF result to target memory" << endl);
1209  memcpy(target, &GF_ints[0], GF_ints.size() * sizeof(dods_int32));
1210  break;
1211  }
1212  case dods_float32_c:
1213  case dods_float64_c: {
1214  // Get the data
1215  BESDEBUG("ugrid",
1216  "TwoDMeshTopology::getResultGFAttributeValues() - Retrieving GF::Array as libdap::Array of libdap::Float64." << endl);
1217  vector<dods_float64> GF_floats = gfa->makeArrayf();
1218 #if 0
1219  stringstream s;
1220  for(unsigned int i=0; i<GF_floats.size(); i++) {
1221  s << "GF_ints[" << i << "]: " << GF_floats[i] << endl;
1222  }
1223  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieved GF_floats: "<< endl << s.str());
1224 #endif
1225  BESDEBUG("ugrid",
1226  "TwoDMeshTopology::getResultGFAttributeValues() - Copying GF result to target memory" << endl);
1227  memcpy(target, &GF_floats[0], GF_floats.size() * sizeof(dods_float64));
1228  break;
1229 
1230  }
1231  default:
1232  throw InternalErr(__FILE__, __LINE__,
1233  "Unknown DAP type encountered when converting to gridfields result values");
1234  }
1235 
1236  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - END" << endl);
1237 }
1238 
1239 static string getIndexVariableName(locationType location)
1240 {
1241  switch (location) {
1242 
1243  case node:
1244  return "node_index";
1245 
1246  case face:
1247  return "face_index";
1248 
1249  case edge:
1250  default:
1251  break;
1252  }
1253 
1254  string msg = "ugr5(): Unknown/Unsupported location value '" + libdap::long_to_string(location) + "'";
1255  BESDEBUG("ugrid", "TwoDMeshTopology::getIndexVariableName() - " << msg << endl);
1256  throw Error(malformed_expr, msg);
1257 }
1258 
1259 int TwoDMeshTopology::getInputGridSize(locationType location)
1260 {
1261  switch (location) {
1262  case node:
1263  return nodeCount;
1264 
1265  case face:
1266  return faceCount;
1267 
1268  case edge:
1269  default:
1270  break;
1271  }
1272 
1273  string msg = "ugr5(): Unknown/Unsupported location value '" + libdap::long_to_string(location) + "'";
1274  BESDEBUG("ugrid", "TwoDMeshTopology::getInputGridSize() - " << msg << endl);
1275  throw Error(malformed_expr, msg);
1276 }
1277 
1281 void TwoDMeshTopology::addIndexVariable(locationType location)
1282 {
1283  int size = getInputGridSize(location);
1284  string name = getIndexVariableName(location);
1285 
1286  BESDEBUG("ugrid",
1287  "TwoDMeshTopology::addIndexVariable() - Adding index variable '" << name << "' size: " << libdap::long_to_string(size) << " at rank " << libdap::long_to_string(location) << endl);
1288 
1289  GF::Array *indexArray = newGFIndexArray(name, size, sharedIntArrays);
1290  d_inputGridField->AddAttribute(location, indexArray);
1291  gfArrays.push_back(indexArray);
1292 }
1293 
1297 void TwoDMeshTopology::getResultIndex(locationType location, void *target)
1298 {
1299  string name = getIndexVariableName(location);
1300  getResultGFAttributeValues(name, dods_int32_c, location, target);
1301 }
1302 
1303 } // namespace ugrid
ugrid::MeshDataVariable
Definition: MeshDataVariable.h:39
ugrid::TwoDMeshTopology::addIndexVariable
void addIndexVariable(locationType location)
Adds an index variable at the gridfields rank as indicated by the passed locationType.
Definition: TwoDMeshTopology.cc:1281
Type
Type
Type of JSON value.
Definition: cmr_module/rapidjson/rapidjson.h:603
libdap
Definition: BESDapFunctionResponseCache.h:35
ugrid::TwoDMeshTopology::getResultGFAttributeValues
void getResultGFAttributeValues(string attrName, libdap::Type type, locationType rank, void *target)
Definition: TwoDMeshTopology.cc:1164
ugrid::TwoDMeshTopology::init
void init(string meshVarName, libdap::DDS *dds)
Definition: TwoDMeshTopology.cc:137
Error
BESError
Abstract exception class for the BES with basic string message.
Definition: BESError.h:58