bes  Updated for version 3.20.6
ugrid_restrict.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 // NOTE: This file is built only when the gridfields library is linked with
29 // the netcdf_handler (i.e., the handler's build is configured using the
30 // --with-gridfields=... option to the 'configure' script).
31 
32 #include "config.h"
33 
34 #include <limits.h>
35 
36 #include <cstdlib> // used by strtod()
37 #include <cerrno>
38 #include <cmath>
39 #include <iostream>
40 #include <sstream>
41 //#include <cxxabi.h>
42 
43 #include <curl/curl.h>
44 
45 #include <BaseType.h>
46 #include <Int32.h>
47 #include <Str.h>
48 #include <Array.h>
49 #include <Structure.h>
50 #include <Error.h>
51 #include <util.h>
52 #include <escaping.h>
53 
54 #include "BESDebug.h"
55 #include "BESError.h"
56 #include "BESStopWatch.h"
57 #include "util.h"
58 
59 #include "ugrid_utils.h"
60 #include "MeshDataVariable.h"
61 #include "TwoDMeshTopology.h"
62 #include "NDimensionalArray.h"
63 #include <gridfields/GFError.h>
64 
65 #include "ugrid_restrict.h"
66 
67 #ifdef NDEBUG
68 #undef BESDEBUG
69 #define BESDEBUG( x, y )
70 #endif
71 
72 using namespace std;
73 using namespace libdap;
74 
75 namespace ugrid {
76 
84 struct UgridRestrictArgs {
94  locationType dimension;
95 
102  vector<libdap::Array *> rangeVars;
103 
107  string filterExpression;
108 };
109 
115 static void addRangeVar(DDS *dds, libdap::Array *rangeVar, map<string, vector<MeshDataVariable *> *> *rangeVariables)
116 {
117  MeshDataVariable *mdv = new MeshDataVariable();
118  mdv->init(rangeVar);
119  string meshVarName = mdv->getMeshName();
120 
121  BaseType *meshVar = dds->var(meshVarName);
122 
123  if (meshVar == 0) {
124  string msg = "The range variable '" + mdv->getName() + "' references the mesh variable '" + meshVarName
125  + "' which cannot be located in this dataset.";
126  throw Error(malformed_expr, msg);
127  }
128 
129  // Get the rangeVariable vector for this mesh name from the map.
130  vector<MeshDataVariable *> *requestedRangeVarsForMesh;
131  map<string, vector<MeshDataVariable *> *>::iterator mit = rangeVariables->find(meshVarName);
132  if (mit == rangeVariables->end()) {
133  // Not there? Make a new one.
134  BESDEBUG("ugrid",
135  "addRangeVar() - MeshTopology object for '" << meshVarName <<"' does not exist. Getting a new one... " << endl);
136 
137  requestedRangeVarsForMesh = new vector<MeshDataVariable *>();
138  (*rangeVariables)[meshVarName] = requestedRangeVarsForMesh;
139  }
140  else {
141  // Sweet! Found it....
142  BESDEBUG("ugrid",
143  "addRangeVar() - MeshTopology object for '" << meshVarName <<"' exists. Retrieving... " << endl);
144  requestedRangeVarsForMesh = mit->second;
145  }
146 
147  requestedRangeVarsForMesh->push_back(mdv);
148 }
149 
150 
151 string usage(string fnc){
152 
153  string usage = fnc+"(rangeVariable:string, [rangeVariable:string, ... ] condition:string)";
154 
155  return usage;
156 }
157 
161 static UgridRestrictArgs processUgrArgs(string func_name, locationType dimension, int argc, BaseType *argv[])
162 {
163 
164  BESDEBUG("ugrid", "processUgrArgs() - BEGIN" << endl);
165 
166  UgridRestrictArgs args;
167  args.dimension = dimension;
168  BESDEBUG("ugrid", "args.dimension: " << libdap::long_to_string(args.dimension) << endl);
169 
170  args.rangeVars = vector<libdap::Array *>();
171 
172  // Check number of arguments;
173  if (argc < 2)
174  throw Error(malformed_expr,
175  "Wrong number of arguments to ugrid restrict function: " + usage(func_name) + " was passed " + long_to_string(argc)
176  + " argument(s)");
177 
178  BaseType * bt;
179 
180 #if 0
181  // ---------------------------------------------
182  // Process the first arg, which is the rank of the Restriction Clause
183  bt = argv[0];
184  if (bt->type() != dods_int32_c)
185  throw Error(malformed_expr,
186  "Wrong type for first argument, expected DAP Int32. " + ugrSyntax + " was passed a/an " + bt->type_name());
187 
188  args.dimension = (locationType) dynamic_cast<Int32&>(*argv[0]).value();
189  BESDEBUG("ugrid", "args.dimension: " << libdap::long_to_string(args.dimension) << endl);
190 #endif
191 
192  // ---------------------------------------------
193  // Process the last argument, the relational/filter expression used to restrict the ugrid content.
194  bt = argv[argc - 1];
195  if (bt->type() != dods_str_c)
196  throw Error(malformed_expr,
197  "Wrong type for third argument, expected DAP String. " + usage(func_name) + " was passed a/an "
198  + bt->type_name());
199 
200  args.filterExpression = dynamic_cast<Str&>(*bt).value();
201  BESDEBUG("ugrid", "args.filterExpression: '" << args.filterExpression << "' (AS RECEIVED)" << endl);
202 
203  args.filterExpression = www2id(args.filterExpression);
204 
205  BESDEBUG("ugrid", "args.filterExpression: '" << args.filterExpression << "' (URL DECODED)" << endl);
206 
207  // --------------------------------------------------
208  // Process the range variables selected by the user.
209  // We know that argc>=3, because we checked so the
210  // following loop will try to find at least one rangeVar,
211  // and it won't try to process the first or last members
212  // of argv.
213  for (int i = 0; i < (argc - 1); i++) {
214  bt = argv[i];
215  if (bt->type() != dods_array_c)
216  throw Error(malformed_expr,
217  "Wrong type for second argument, expected DAP Array. " + usage(func_name) + " was passed a/an "
218  + bt->type_name());
219 
220  libdap::Array *newRangeVar = dynamic_cast<libdap::Array*>(bt);
221  if (newRangeVar == 0) {
222  throw Error(malformed_expr,
223  "Wrong type for second argument. " + usage(func_name) + " was passed a/an " + bt->type_name());
224  }
225  args.rangeVars.push_back(newRangeVar);
226  }
227 
228  BESDEBUG("ugrid", "processUgrArgs() - END" << endl);
229 
230  return args;
231 }
232 
233 static string arrayState(libdap::Array *dapArray, string indent)
234 {
235 
236  libdap::Array::Dim_iter thisDim;
237  stringstream arrayState;
238  arrayState << endl;
239  arrayState << indent << "dapArray: '" << dapArray->name() << "' ";
240  arrayState << "type: '" << dapArray->type_name() << "' ";
241 
242  arrayState << "read(): " << dapArray->read() << " ";
243  arrayState << "read_p(): " << dapArray->read_p() << " ";
244  arrayState << endl;
245 
246  for (thisDim = dapArray->dim_begin(); thisDim != dapArray->dim_end(); thisDim++) {
247 
248  arrayState << indent << " dim: '" << dapArray->dimension_name(thisDim) << "' ";
249  arrayState << indent << " start: " << dapArray->dimension_start(thisDim) << " ";
250  arrayState << indent << " stride: " << dapArray->dimension_stride(thisDim) << " ";
251  arrayState << indent << " stop: " << dapArray->dimension_stop(thisDim) << " ";
252  arrayState << endl;
253  }
254  return arrayState.str();
255 }
256 
257 static void copyUsingSubsetIndex(libdap::Array *sourceArray, vector<unsigned int> *subsetIndex, void *result)
258 {
259  BESDEBUG("ugrid", "ugrid::copyUsingIndex() - BEGIN" << endl);
260 
261  switch (sourceArray->var()->type()) {
262  case dods_byte_c:
263  sourceArray->value(subsetIndex, (dods_byte *) result);
264  break;
265  case dods_uint16_c:
266  sourceArray->value(subsetIndex, (dods_uint16 *) result);
267  break;
268  case dods_int16_c:
269  sourceArray->value(subsetIndex, (dods_int16 *) result);
270  break;
271  case dods_uint32_c:
272  sourceArray->value(subsetIndex, (dods_uint32 *) result);
273  break;
274  case dods_int32_c:
275  sourceArray->value(subsetIndex, (dods_int32 *) result);
276  break;
277 
278  case dods_float32_c:
279  sourceArray->value(subsetIndex, (dods_float32 *) result);
280  break;
281  case dods_float64_c:
282  sourceArray->value(subsetIndex, (dods_float64 *) result);
283  break;
284 
285  default:
286  throw InternalErr(__FILE__, __LINE__, "ugrid::hgr5::copyUsingSubsetIndex() - Unknown DAP type encountered.");
287  }
288 
289  BESDEBUG("ugrid", "ugrid::copyUsingIndex() - END" << endl);
290 }
291 
292 // This is only used for the ugrid2 BESDEBUG lines.
293 static string vectorToString(vector<unsigned int> *index)
294 {
295  BESDEBUG("ugrid", "indexToString() - BEGIN"<< endl);
296  BESDEBUG("ugrid", "indexToString() - index.size(): " << libdap::long_to_string(index->size()) << endl);
297  stringstream s;
298  s << "[";
299 
300  for (unsigned int i = 0; i < index->size(); ++i) {
301  s << ((i > 0) ? ", " : " ");
302  s << (*index)[i];
303  }
304  s << "]";
305  BESDEBUG("ugrid", "indexToString() - END"<< endl);
306 
307  return s.str();
308 }
309 
313 static void rDAWorker(MeshDataVariable *mdv, libdap::Array::Dim_iter thisDim, vector<unsigned int> *slab_subset_index,
314  NDimensionalArray *results)
315 {
316  libdap::Array *dapArray = mdv->getDapArray();
317 
318  // For real data, this output is huge. jhrg 4/15/15
319  BESDEBUG("ugrid2",
320  "rDAWorker() - slab_subset_index" << vectorToString(slab_subset_index) << " size: " << slab_subset_index->size() << endl);
321 
322  // The locationCoordinateDimension is the dimension of the array that is associated with the ugrid "rank" - e.g. it is the
323  // dimension that ties the variable to the 'nodes' (rank 0) or 'edges' (rank 1) or 'faces' (rank 2) of the ugrid.
324  libdap::Array::Dim_iter locationCoordinateDim = mdv->getLocationCoordinateDimension();
325 
326  BESDEBUG("ugrid", "rdaWorker() - thisDim: '" << dapArray->dimension_name(thisDim) << "'" << endl);
327  BESDEBUG("ugrid",
328  "rdaWorker() - locationCoordinateDim: '" << dapArray->dimension_name(locationCoordinateDim) << "'" << endl);
329  // locationCoordinateDim is, e.g., 'nodes'. jhrg 10/25/13
330  if (thisDim != locationCoordinateDim) {
331  unsigned int start = dapArray->dimension_start(thisDim, true);
332  unsigned int stride = dapArray->dimension_stride(thisDim, true);
333  unsigned int stop = dapArray->dimension_stop(thisDim, true);
334 
335  BESDEBUG("ugrid", "rdaWorker() - array state: " << arrayState(dapArray, " "));
336 
337  for (unsigned int dimIndex = start; dimIndex <= stop; dimIndex += stride) {
338  dapArray->add_constraint(thisDim, dimIndex, 1, dimIndex);
339  rDAWorker(mdv, thisDim + 1, slab_subset_index, results);
340  }
341 
342  // Reset the constraint for this dimension.
343  dapArray->add_constraint(thisDim, start, stride, stop);
344  }
345  else {
346  BESDEBUG("ugrid", "rdaWorker() - Found locationCoordinateDim" << endl);
347 
348  if ((thisDim + 1) != dapArray->dim_end()) {
349  string msg =
350  "rDAWorker() - The location coordinate dimension is not the last dimension in the array. Hyperslab subsetting of this dimension is not supported.";
351  BESDEBUG("ugrid", msg << endl);
352  throw Error(malformed_expr, msg);
353  }
354 
355  BESDEBUG("ugrid", "rdaWorker() - Arrived At Last Dimension" << endl);
356 
357  dapArray->set_read_p(false);
358 
359  BESDEBUG("ugrid", "rdaWorker() - dap array: " << arrayState(dapArray, " "));
360 
361  vector<unsigned int> lastDimHyperSlabLocation;
362  NDimensionalArray::retrieveLastDimHyperSlabLocationFromConstrainedArrray(dapArray, &lastDimHyperSlabLocation);
363 
364  BESDEBUG("ugrid",
365  "rdaWorker() - lastDimHyperSlabLocation: " << NDimensionalArray::vectorToIndices(&lastDimHyperSlabLocation) << endl);
366 
367  // unused. 4/7/14 jhrg. unsigned int elementCount;
368 
369  void *slab;
370  //results->getLastDimensionHyperSlab(&lastDimHyperSlabLocation, &slab, &elementCount);
371  results->getNextLastDimensionHyperSlab(&slab);
372 
373  dapArray->read();
374 
375  copyUsingSubsetIndex(dapArray, slab_subset_index, slab);
376  }
377 }
378 
386 static libdap::Array *restrictRangeVariableByOneDHyperSlab(MeshDataVariable *mdv,
387  vector<unsigned int> *slab_subset_index)
388 {
389 
390  long restrictedSlabSize = slab_subset_index->size();
391 
392  BESDEBUG("ugrid2",
393  "restrictRangeVariableByOneDHyperSlab() - slab_subset_index"<< vectorToString(slab_subset_index) << " size: " << libdap::long_to_string(restrictedSlabSize) << endl);
394 
395  libdap::Array *sourceDapArray = mdv->getDapArray();
396 
397  BESDEBUG("ugrid",
398  "restrictDapArrayByOneDHyperSlab() - locationCoordinateDim: '" << sourceDapArray->dimension_name(mdv->getLocationCoordinateDimension()) << "'" << endl);
399 
400  // We want the manipulate the Array's Dimensions so that only a single dimensioned slab of the location coordinate dimension
401  // is read at a time. We need to cache the original constrained dimensions so that we can build the correct collection of
402  // location coordinate dimension slabs.
403 
404  // Now we need to compute the shape of the final subset result array from the source range variable array and the slab subset.
405  // What's the shape of the source array with any constraint applied?
406  vector<unsigned int> resultArrayShape(sourceDapArray->dimensions(true));
407  NDimensionalArray::computeConstrainedShape(sourceDapArray, &resultArrayShape);
408 
409  stringstream msg;
410  for (unsigned int i = 0; i < resultArrayShape.size(); i++) {
411  msg << "[" << resultArrayShape[i] << "]";
412  }
413  BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - Constrained source array shape" << msg.str() << endl);
414  msg.str("");
415 
416  // Now, we know that the result array has a last dimension slab size determined by the slab_subset_index (whichwas made by the ugrid sub-setting),
417  // so we make the result array shape reflect that
418 
419  resultArrayShape[sourceDapArray->dimensions(true) - 1] = restrictedSlabSize; // jhrg 10/25/13 resultArrayShape.last() = ...
420  libdap::Type dapType = sourceDapArray->var()->type();
421 
422  BESDEBUG("ugrid",
423  "restrictDapArrayByOneDHyperSlab() - UGrid restricted HyperSlab has "<< restrictedSlabSize << " elements." << endl);
424  BESDEBUG("ugrid",
425  "restrictDapArrayByOneDHyperSlab() - Array is of type '"<< libdap::type_name(dapType) << "'" << endl);
426 
427  for (unsigned int i = 0; i < resultArrayShape.size(); i++) {
428  msg << "[" << resultArrayShape[i] << "]";
429  }
430  BESDEBUG("ugrid", "restrictDapArrayByOneDHyperSlab() - resultArrayShape" << msg.str() << endl);
431  msg.str(std::string());
432 
433  // Now we make a new NDimensionalArray instance that we will use to hold the results.
434  NDimensionalArray *result = new NDimensionalArray(&resultArrayShape, dapType);
435 
436  // And we pass that along with other stuff into the recursive rDAWorker that's going to go get all the stuff
437  rDAWorker(mdv, sourceDapArray->dim_begin(), slab_subset_index, result);
438 
439  // And now that the recursion we grab have the NDimensionalArray cough up the rteuslt as a libdap::Array
440  libdap::Array *resultDapArray = result->getArray(sourceDapArray);
441 
442  // Delete the NDimensionalArray
443  delete result;
444 
445  // Return the result as libdap:Array
446  return resultDapArray;
447 }
448 
468 void ugrid_restrict(string func_name, locationType location, int argc, BaseType *argv[], DDS &dds, BaseType **btpp)
469 {
470  try { // This top level try block is used to catch gridfields library errors.
471 
472  BESStopWatch sw;
473  if (BESISDEBUG(TIMING_LOG)) sw.start("ugrid::ugrid_restrict()", "[function_invocation]");
474 
475  BESDEBUG("ugrid", "ugrid_restrict() - BEGIN" << endl);
476 
477  string info = string("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
478  + "<function name=\"ugrid_restrict\" version=\"0.1\">\n" + "Server function for Unstructured grid operations.\n"
479  + "usage: " + usage(func_name) + "\n" + "</function>";
480 
481  if (argc == 0) {
482  Str *response = new Str("info");
483  response->set_value(info);
484  *btpp = response;
485  return;
486  }
487 
488  // Process and QC the arguments
489  UgridRestrictArgs args = processUgrArgs(func_name, location, argc, argv);
490 
491  // Each range variable is associated with a "mesh" i.e. a mesh topology variable. Since there may be more than one mesh in a
492  // dataset, and the user may request more than one range variable for each mesh we need to sift through the list of requested
493  // range variables and organize them by mesh topology variable name.
494  map<string, vector<MeshDataVariable *> *> *meshToRangeVarsMap = new map<string, vector<MeshDataVariable *> *>();
495 
496  // For every Range variable in the arguments list, locate it and ingest it.
497  vector<libdap::Array *>::iterator it;
498  for (it = args.rangeVars.begin(); it != args.rangeVars.end(); ++it) {
499  addRangeVar(&dds, *it, meshToRangeVarsMap);
500  }
501  BESDEBUG("ugrid", "ugrid_restrict() - The user requested "<< args.rangeVars.size() << " range data variables." << endl);
502  BESDEBUG("ugrid",
503  "ugrid_restrict() - The user's request referenced "<< meshToRangeVarsMap->size() << " mesh topology variables." << endl);
504 
505  // ----------------------------------
506  // OK, so up to this point we have not read any data from the data set, but we have QC'd the inputs and verified that
507  // it looks like the request is consistent with the semantics of the dataset.
508  // Now it's time to read some data and pack it into the GridFields library...
509 
510  // TODO This returns a single structure but it would make better sense to the
511  // world if it could return a vector of objects and have them appear at the
512  // top level of the DDS.
513  // FIXME fix the names of the variables in the mesh_topology attributes
514  // If the server side function can be made to return a DDS or a collection of BaseType's then the
515  // names won't change and the original mesh_topology variable and it's metadata will be valid
516  Structure *dapResult = new Structure("ugr_result_unwrap");
517 
518  // Now we need to grab an top level metadata (attriubutes) and copy them into the dapResult Structure
519  // Add any global attributes to the netcdf file
520 #if 1
521  AttrTable &globals = dds.get_attr_table();
522  BESDEBUG("ugrid", "ugrid_restrict() - Copying Global Attributes" << endl << globals << endl);
523 
524  AttrTable *newDatasetAttrTable = new AttrTable(globals);
525  dapResult->set_attr_table(*newDatasetAttrTable);
526  BESDEBUG("ugrid", "ugrid_restrict() - Result Structure attrs: " << endl << dapResult->get_attr_table() << endl);
527 
528 
529 #endif
530 
531  // Since we only want each ugrid structure to appear in the results one time (cause otherwise we might be trying to add
532  // the same variables with the same names to the result multiple times.) we grind on this by iterating over the
533  // names of the mesh topology names.
534  map<string, vector<MeshDataVariable *> *>::iterator mit;
535  for (mit = meshToRangeVarsMap->begin(); mit != meshToRangeVarsMap->end(); ++mit) {
536 
537  string meshVariableName = mit->first;
538  vector<MeshDataVariable *> *requestedRangeVarsForMesh = mit->second;
539 
540  // Building the restricted TwoDMeshTopology without adding any range variables and then converting the result
541  // Grid field to Dap Objects should return all of the Ugrid structural stuff - mesh variable, node coordinate variables,
542  // face and edge coordinate variables if present.
543  BESDEBUG("ugrid",
544  "ugrid_restrict() - Adding restricted mesh_topology structure for mesh '" << meshVariableName << "' to DAP response." << endl);
545 
546  TwoDMeshTopology *tdmt = new TwoDMeshTopology();
547  tdmt->init(meshVariableName, &dds);
548 
549  tdmt->buildBasicGfTopology();
550  tdmt->addIndexVariable(node);
551  tdmt->addIndexVariable(face);
552  tdmt->applyRestrictOperator(args.dimension, args.filterExpression);
553 
554  // 3: because there are nodes (rank = 0), edges (rank = 1), and faces (rank = 2). jhrg 10/25/13
555  vector<vector<unsigned int> *> location_subset_indices(3);
556 
557  long nodeResultSize = tdmt->getResultGridSize(node);
558  BESDEBUG("ugrid", "ugrid_restrict() - there are "<< nodeResultSize << " nodes in the subset." << endl);
559  vector<unsigned int> node_subset_index(nodeResultSize);
560  if (nodeResultSize > 0) {
561  tdmt->getResultIndex(node, &node_subset_index[0]);
562  }
563  location_subset_indices[node] = &node_subset_index;
564 
565  BESDEBUG("ugrid2", "ugrid_restrict() - node_subset_index"<< vectorToString(&node_subset_index) << endl);
566 
567  long faceResultSize = tdmt->getResultGridSize(face);
568  BESDEBUG("ugrid", "ugrid_restrict() - there are "<< faceResultSize << " faces in the subset." << endl);
569  vector<unsigned int> face_subset_index(faceResultSize);
570  if (faceResultSize > 0) {
571  tdmt->getResultIndex(face, &face_subset_index[0]);
572  }
573  location_subset_indices[face] = &face_subset_index;
574  BESDEBUG("ugrid2", "ugrid_restrict() - face_subset_index: "<< vectorToString(&face_subset_index) << endl);
575 
576  // This gets all the stuff that's attached to the grid - which at this point does not include the range variables but does include the
577  // index variable. good enough for now but need to drop the index....
578  vector<BaseType *> dapResults;
579  tdmt->convertResultGridFieldStructureToDapObjects(&dapResults);
580 
581  BESDEBUG("ugrid",
582  "ugrid_restrict() - Restriction of mesh_topology '"<< tdmt->getMeshVariable()->name() << "' structure completed." << endl);
583 
584  // now that we have the mesh topology variable we are going to look at each of the requested
585  // range variables (aka MeshDataVariable instances) and we're going to subset that using the
586  // gridfields library and add its subset version to the results.
587  vector<MeshDataVariable *>::iterator rvit;
588  for (rvit = requestedRangeVarsForMesh->begin(); rvit != requestedRangeVarsForMesh->end(); rvit++) {
589  MeshDataVariable *mdv = *rvit;
590 
591  BESDEBUG("ugrid",
592  "ugrid_restrict() - Processing MeshDataVariable '"<< mdv->getName() << "' associated with rank/location: "<< mdv->getGridLocation() << endl);
593 
594  tdmt->setLocationCoordinateDimension(mdv);
595 
600  libdap::Array *restrictedRangeVarArray = restrictRangeVariableByOneDHyperSlab(mdv,
601  location_subset_indices[mdv->getGridLocation()]);
602 
603  BESDEBUG("ugrid",
604  "ugrid_restrict() - Adding resulting dapArray '"<< restrictedRangeVarArray->name() << "' to dapResults." << endl);
605 
606  dapResults.push_back(restrictedRangeVarArray);
607  }
608 
609  delete tdmt;
610 
611  BESDEBUG("ugrid", "ugrid_restrict() - Adding GF::GridField results to DAP structure " << dapResult->name() << endl);
612 
613  for (vector<BaseType *>::iterator i = dapResults.begin(); i != dapResults.end(); ++i) {
614  BESDEBUG("ugrid",
615  "ugrid_restrict() - Adding variable "<< (*i)->name() << " to DAP structure " << dapResult->name() << endl);
616  dapResult->add_var_nocopy(*i);
617  }
618  }
619 
620  *btpp = dapResult;
621 
622  // TODO replace with auto_ptr. jhrg 4/17/15
623 
624  BESDEBUG("ugrid", "ugrid_restrict() - Releasing maps and vectors..." << endl);
625  for (mit = meshToRangeVarsMap->begin(); mit != meshToRangeVarsMap->end(); ++mit) {
626  vector<MeshDataVariable *> *requestedRangeVarsForMesh = mit->second;
627  vector<MeshDataVariable *>::iterator rvit;
628  for (rvit = requestedRangeVarsForMesh->begin(); rvit != requestedRangeVarsForMesh->end(); rvit++) {
629  MeshDataVariable *mdv = *rvit;
630  delete mdv;
631  }
632  delete requestedRangeVarsForMesh;
633  }
634  delete meshToRangeVarsMap;
635 
636  BESDEBUG("ugrid", "ugrid_restrict() - END" << endl);
637  }
638  catch (GFError &gfe) {
639  throw BESError(gfe.get_message(), gfe.get_error_type(), gfe.get_file(), gfe.get_line());
640  }
641 
642  return;
643 }
644 
645 
646 
647 
648 
667 void ugnr(int argc, BaseType *argv[], DDS &dds, BaseType **btpp) {
668  ugrid_restrict("ugnr",node,argc,argv,dds,btpp);
669 }
670 
671 
672 
691 void uger(int argc, BaseType *argv[], DDS &dds, BaseType **btpp) {
692  ugrid_restrict("uger",edge,argc,argv,dds,btpp);
693 }
694 
695 
714 void ugfr(int argc, BaseType *argv[], DDS &dds, BaseType **btpp) {
715  ugrid_restrict("ugfr",face,argc,argv,dds,btpp);
716 }
717 
718 
719 
720 
721 } // namespace ugrid_restrict
BESStopWatch::start
virtual bool start(std::string name)
Definition: BESStopWatch.cc:58
libdap::NDimensionalArray
Definition: NDimensionalArray.h:39
Type
Type
Type of JSON value.
Definition: cmr_module/rapidjson/rapidjson.h:603
libdap
Definition: BESDapFunctionResponseCache.h:35
BESStopWatch
Definition: BESStopWatch.h:55
Error
BESError
Abstract exception class for the BES with basic string message.
Definition: BESError.h:58