OPeNDAP Hyrax Back End Server (BES)  Updated for version 3.8.3
TwoDMeshTopology.cc
Go to the documentation of this file.
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 
33 #include <gridfields/type.h>
34 #include <gridfields/gridfield.h>
35 #include <gridfields/grid.h>
36 #include <gridfields/implicit0cells.h>
37 #include <gridfields/gridfieldoperator.h>
38 #include <gridfields/restrict.h>
39 #include <gridfields/refrestrict.h>
40 #include <gridfields/array.h>
41 
42 
43 #include "BaseType.h"
44 #include "Int32.h"
45 #include "Float64.h"
46 #include "Array.h"
47 #include "util.h"
48 
49 
50 #include "ugrid_utils.h"
51 #include "NDimensionalArray.h"
52 #include "MeshDataVariable.h"
53 #include "TwoDMeshTopology.h"
54 
55 #include "BESDebug.h"
56 
57 using namespace std;
58 using std::endl;
59 using namespace libdap;
60 using namespace ugrid;
61 
62 namespace ugrid {
63 
64 
65 TwoDMeshTopology::TwoDMeshTopology():
66  d_myMeshVar(0),
67  nodeCoordinateArrays(0),
68  nodeCount(0),
69  faceNodeConnectivityArray(0),
70  faceCount(0),
71  faceCoordinateNames(0),
72  faceCoordinateArrays(0),
73  gridTopology(0),
74  d_inputGridField(0),
75  resultGridField(0),
76  fncCellArray(0),
77  _initialized(false)
78 
79 {
80  rangeDataArrays = new vector<MeshDataVariable *>();
81  sharedIntArrays = new vector<int *>();
82  sharedFloatArrays = new vector<float *>();
83 
84 }
85 
86 
88 {
89  BESDEBUG("ugrid", "~TwoDMeshTopology() - BEGIN" << endl);
90  BESDEBUG("ugrid", "~TwoDMeshTopology() - (" << this << ")" << endl);
91 
92  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::GridField 'resultGridField'." << endl);
93  delete resultGridField;
94 
95 
96 
97  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::GridField 'inputGridField'." << endl);
98  delete d_inputGridField;
99 
100  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::Grid 'gridTopology'." << endl);
101  delete gridTopology;
102 
103 
104  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF::Arrays..." << endl);
105  for (vector<GF::Array *>::iterator it = gfArrays.begin(); it != gfArrays.end(); ++it) {
106  GF::Array *gfa = *it;
107  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting GF:Array '" << gfa->getName() << "'" << endl);
108  delete gfa;
109  }
110  BESDEBUG("ugrid", "~TwoDMeshTopology() - GF::Arrays deleted." << endl);
111 
112 
113 
114  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting sharedIntArrays..." << endl);
115  for (vector<int *>::iterator it = sharedIntArrays->begin(); it != sharedIntArrays->end(); ++it) {
116  int *i = *it;
117  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting integer array '" << i << "'" << endl);
118  delete [] i;
119  }
120  delete sharedIntArrays;
121 
122  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting sharedFloatArrays..." << endl);
123  for (vector<float *>::iterator it = sharedFloatArrays->begin(); it != sharedFloatArrays->end(); ++it) {
124  float *f = *it;
125  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting float array '" << f << "'" << endl);
126  delete [] f;
127  }
128  delete sharedFloatArrays;
129 
130  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting range data vector" << endl);
131  delete rangeDataArrays;
132  BESDEBUG("ugrid", "~TwoDMeshTopology() - Range data vector deleted." << endl);
133 
134 
135  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting vector of node coordinate arrays." << endl);
136  delete nodeCoordinateArrays;
137 
138  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting vector of face coordinate arrays." << endl);
139  delete faceCoordinateArrays;
140 
141  BESDEBUG("ugrid", "~TwoDMeshTopology() - Deleting face node connectivity cell array (GF::Node's)." << endl);
142  delete[] fncCellArray;
143 
144 
145  BESDEBUG("ugrid", "~TwoDMeshTopology() - END" << endl);
146 }
147 
148 
152 void TwoDMeshTopology::init(string meshVarName, DDS *dds)
153 {
154 
155  if(_initialized)
156  return;
157 
158  BaseType *bt = dds->var(meshVarName);
159 
160  if(bt)
161  d_myMeshVar = bt;
162  else
163  throw Error("Unable to locate variable: "+meshVarName);
164 
165 
166  BaseType *meshVar = getMeshVariable();
167 
168  dimension = getAttributeValue(meshVar,UGRID_TOPOLOGY_DIMENSION);
169 
170  if (dimension.empty()) {
171  string msg = "TwoDMeshTopology::init() - The mesh topology variable '" + meshVar->name() + "' is missing the required attribute named '" + UGRID_TOPOLOGY_DIMENSION + "'";
172  BESDEBUG("ugrid",msg );
173  throw Error( msg);
174  }
175 
176  // Retrieve the node coordinate arrays for the mesh
177  ingestNodeCoordinateArrays(meshVar, dds);
178 
179 
180  // Why would we retrieve this? I think Bill said that this needs to be recomputed after a restrict operation.
181  // @TODO Verify that Bill actually said that this needs to be recomputed.
182  // Retrieve the face coordinate arrays (if any) for the mesh
183  ingestFaceCoordinateArrays(meshVar,dds);
184 
185 
186  // Inspect and QC the face node connectivity array for the mesh
187  ingestFaceNodeConnectivityArray(meshVar,dds);
188 
189 
190 
191  _initialized = true;
192 
193 }
194 
195 
196 
197 
198 void TwoDMeshTopology::setNodeCoordinateDimension(MeshDataVariable *mdv)
199 {
200  BESDEBUG("ugrid", "TwoDMeshTopology::setNodeCoordinateDimension() - BEGIN" << endl);
201  libdap::Array *dapArray = mdv->getDapArray();
202  libdap::Array::Dim_iter ait1;
203 
204 
205  for (ait1 = dapArray->dim_begin(); ait1 != dapArray->dim_end(); ++ait1) {
206 
207  string dimName = dapArray->dimension_name(ait1);
208 
209  if(dimName.compare(nodeDimensionName) == 0){ // are the names the same?
210  BESDEBUG("ugrid", "TwoDMeshTopology::setNodeCoordinateDimension() - Found dimension name matching nodeDimensionName '"<< nodeDimensionName << "'" << endl);
211  int size = dapArray->dimension_size(ait1,true);
212  if(size == nodeCount){ // are they the same size?
213  BESDEBUG("ugrid", "TwoDMeshTopology::setNodeCoordinateDimension() - Dimension sizes match (" << libdap::long_to_string(nodeCount) << ") - DONE" << endl);
214  // Yay! We found the node coordinate dimension
216  return;
217  }
218  }
219 
220  }
221 
222  throw Error(
223  "Unable to determine the node coordinate dimension of the range variable '"
224  + mdv->getName() + "' The node dimension is named '"+ nodeDimensionName
225  + "' with size " + libdap::long_to_string(nodeCount)
226  );
227 
228 
229 }
230 
231 
232 void TwoDMeshTopology::setFaceCoordinateDimension(MeshDataVariable *mdv)
233 {
234 
235 
236  libdap::Array *dapArray = mdv->getDapArray();
237  libdap::Array::Dim_iter ait1;
238 
239  for (ait1 = dapArray->dim_begin(); ait1 != dapArray->dim_end(); ++ait1) {
240  string dimName = dapArray->dimension_name(ait1);
241 
242  if(dimName.compare(faceDimensionName) == 0){ // are the names the same?
243  int size = dapArray->dimension_size(ait1,true);
244  if(size == faceCount){ // are they the same size?
245  // Yay! We found the coordinate dimension
247  return;
248  }
249  }
250 
251  }
252  throw Error(
253  "Unable to determine the face coordinate dimension of the range variable '"
254  + mdv->getName() + "' The face coordinate dimension is named '"+ faceDimensionName
255  + "' with size " + libdap::long_to_string(faceCount)
256  );
257 
258 
259 
260 }
261 
262 
264 
265  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - BEGIN" << endl);
266 
267  libdap::Array *dapArray = mdv->getDapArray();
268 
269  string locstr;
270  switch(mdv->getGridLocation()){
271 
272  case node:
273  {
274  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - Checking Node variable '"<< mdv->getName() << "'" << endl);
275  // Locate and set the MDV's node coordinate dimension.
276  setNodeCoordinateDimension(mdv);
277  locstr = "node";
278  }
279  break;
280 
281  case face:
282  {
283  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - Checking Face variable '"<< mdv->getName() << "'" << endl);
284  // Locate and set the MDV's face coordinate dimension.
285  setFaceCoordinateDimension(mdv);
286  locstr = "face";
287  }
288  break;
289 
290  case edge:
291  default:
292  {
293  string msg = "TwoDMeshTopology::setLocationCoordinateDimension() - Unknown/Unsupported location value '" +
294  libdap::long_to_string(mdv->getGridLocation()) + "'";
295  BESDEBUG("ugrid", msg << endl);
296  throw Error( msg );
297  }
298  break;
299  }
300  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - MeshDataVariable '"<< mdv->getName() <<
301  "' is a "<< locstr <<" variable," << endl);
302  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - MeshDataVariable '"<< mdv->getName() <<
303  "' location coordinate dimension is '" << dapArray->dimension_name(mdv->getLocationCoordinateDimension()) << "'" << endl);
304 
305  BESDEBUG("ugrid", "TwoDMeshTopology::setLocationCoordinateDimension() - DONE" << endl);
306 
307 }
308 
309 
310 #if 0 // Unused crufty code from previous efforts
311 void TwoDMeshTopology::addDataVariable(MeshDataVariable *mdv)
312 {
313 
314  BESDEBUG("ugrid", "TwoDMeshTopology::addDataVariable() - BEGIN" << endl);
316  BESDEBUG("ugrid", "TwoDMeshTopology::addDataVariable() - Adding MeshDataVariable '"<< mdv->getName() << "' " << endl);
317  rangeDataArrays->push_back(mdv);
318  BESDEBUG("ugrid", "TwoDMeshTopology::addDataVariable() - DONE" << endl);
319 
320 }
321 #endif
322 
323 
324 
329 void TwoDMeshTopology::ingestFaceNodeConnectivityArray(libdap::BaseType *meshTopology, libdap::DDS *dds)
330 {
331 
332  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Locating FNCA" << endl);
333 
334  string face_node_connectivity_var_name;
335  AttrTable at = meshTopology->get_attr_table();
336 
337  AttrTable::Attr_iter iter_fnc = at.simple_find(UGRID_FACE_NODE_CONNECTIVITY);
338  if (iter_fnc != at.attr_end()) {
339  face_node_connectivity_var_name = at.get_attr(iter_fnc, 0);
340  }
341  else {
342  throw Error("Could not locate the " UGRID_FACE_NODE_CONNECTIVITY " attribute in the " UGRID_MESH_TOPOLOGY " variable! "
343  "The mesh_topology variable is named " + meshTopology->name());
344  }
345  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Located the '" << UGRID_FACE_NODE_CONNECTIVITY << "' attribute." << endl);
346 
347  // Find the variable using the name
348 
349  BaseType *btp = dds->var(face_node_connectivity_var_name);
350 
351  if(btp==0)
352  throw Error("Could not locate the " UGRID_FACE_NODE_CONNECTIVITY " variable named '" + face_node_connectivity_var_name + "'! "+
353  "The mesh_topology variable is named "+meshTopology->name());
354 
355 
356  // Is it an array?
357  libdap::Array *fncArray = dynamic_cast<libdap::Array*>(btp);
358  if(fncArray == 0) {
359  throw Error(malformed_expr,"Face Node Connectivity variable '"+face_node_connectivity_var_name+"' is not an Array type. It's an instance of " + btp->type_name());
360  }
361 
362  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Located FNC Array '" << fncArray->name() << "'." << endl);
363 
364  // It's got to have exactly 2 dimensions - [max#nodes_per_face][#faces]
365  int numDims = fncArray->dimensions(true) ;
366  if( numDims != 2){
367  throw Error(malformed_expr,"Face Node Connectivity variable '"+face_node_connectivity_var_name
368  +"' Must have two (2) dimensions. It has " + libdap::long_to_string(numDims));
369  }
370 
371  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - FNC Array '" << fncArray->name() << "' has two (2) dimensions." << endl);
372 
373 
374  // We just need to have both dimensions handy, so get the first and the second
375  libdap::Array::Dim_iter firstDim = fncArray->dim_begin();
376  libdap::Array::Dim_iter secondDim = fncArray->dim_begin(); // same as the first for a moment...
377  secondDim++; // now it's second!
378 
379 
380  if(faceDimensionName.empty()){
381  // By now we know it only has two dimensions, but since there is no promise that they'll be in a particular order
382  // we punt: We'll assume that smallest of the two is in fact the nodes per face and the larger the face index dimensions.
383  int sizeFirst = fncArray->dimension_size(firstDim,true);
384  int sizeSecond = fncArray->dimension_size(secondDim,true);
385  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - sizeFirst: "<< libdap::long_to_string(sizeFirst) << endl);
386  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - sizeSecond: "<< libdap::long_to_string(sizeSecond) << endl);
387 
388  if(sizeFirst < sizeSecond){
389  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - FNC Array first dimension is smaller than second." << endl);
390  fncNodesDim = firstDim;
391  fncFacesDim = secondDim;
392  }
393  else {
394  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - FNC Array second dimension is smaller than first." << endl);
395  fncNodesDim = secondDim;
396  fncFacesDim = firstDim;
397  }
398 
399  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - fncFacesDim name: '" << fncArray->dimension_name(fncFacesDim) << "'" << endl);
400  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - fncFacesDim size: '" << fncArray->dimension_size(fncFacesDim,true) << "'" << endl);
401 
402  faceDimensionName = fncArray->dimension_name(fncFacesDim);
403  }
404  else {
405  // There is already a faceDimensionName defined - possibly from loading face coordinate variables.
406 
407  // Does it match the name of the first or second dimensions of the fncArray? It better!
408  if(faceDimensionName.compare(fncArray->dimension_name(firstDim)) == 0){
409  fncNodesDim = secondDim;
410  fncFacesDim = firstDim;
411  }
412  else if(faceDimensionName.compare(fncArray->dimension_name(secondDim)) == 0){
413  fncNodesDim = firstDim;
414  fncFacesDim = secondDim;
415  }
416  else {
417  string msg = "The face coordinate dimension of the Face Node Connectivity variable '"+face_node_connectivity_var_name
418  +"' Has dimension name.'"+ fncArray->dimension_name(fncFacesDim) + "' which does not match the existing face coordinate dimension name '"
419  + faceDimensionName + "'";
420  BESDEBUG("ugrid", msg << endl);
421  throw Error(msg);
422  }
423  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face dimension names match." << endl);
424 
425  }
426 
427  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face dimension name: '" << faceDimensionName << "'" << endl);
428 
429  // Check to see if faceCount is initialized and do so if needed
430  if(faceCount==0) {
431  faceCount = fncArray->dimension_size(fncFacesDim,true);
432  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face count: "<< libdap::long_to_string(faceCount) << endl);
433  }
434  else {
435  // Make sure the face counts match.
436 
437  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face count: "<< libdap::long_to_string(faceCount) << endl);
438  if(faceCount!=fncArray->dimension_size(fncFacesDim,true) ){
439  string msg = "The faces dimension of the Face Node Connectivity variable '"+face_node_connectivity_var_name
440  +"' Has size "+ libdap::long_to_string(fncArray->dimension_size(fncFacesDim,true))+ " which does not match the existing face count of "
441  + libdap::long_to_string(faceCount);
442  BESDEBUG("ugrid", msg << endl);
443  throw Error(msg);
444  }
445  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceNodeConnectivityArray() - Face counts match!" << endl);
446  }
447 
448  faceNodeConnectivityArray = fncArray;
449 
450  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityArray() - Got FCNA '"+fncArray->name()+"'" << endl);
451 
452 
453 }
459 void TwoDMeshTopology::ingestFaceCoordinateArrays(libdap::BaseType *meshTopology, libdap::DDS *dds)
460 {
461  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - BEGIN Gathering face coordinate arrays..." << endl);
462 
463  if(faceCoordinateArrays==0)
464  faceCoordinateArrays = new vector<libdap::Array *>();
465 
466  faceCoordinateArrays->clear();
467 
468  string face_coordinates;
469  AttrTable at = meshTopology->get_attr_table();
470 
471  AttrTable::Attr_iter iter_nodeCoors = at.simple_find(UGRID_FACE_COORDINATES);
472  if (iter_nodeCoors != at.attr_end()) {
473  face_coordinates = at.get_attr(iter_nodeCoors, 0);
474 
475  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Located '"<< UGRID_FACE_COORDINATES << "' attribute." << endl);
476 
477 
478 
479  // Split the node_coordinates string up on spaces
480  vector<string> faceCoordinateNames = split(face_coordinates, ' ');
481 
482  // Find each variable in the resulting list
483  vector<string>::iterator coorName_it;
484  for (coorName_it = faceCoordinateNames.begin();
485  coorName_it != faceCoordinateNames.end(); ++coorName_it) {
486  string faceCoordinateName = *coorName_it;
487 
488  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Processing face coordinate '"<< faceCoordinateName << "'." << endl);
489 
490  //Now that we have the name of the coordinate variable get it from the DDS!!
491  BaseType *btp = dds->var(faceCoordinateName);
492  if (btp == 0)
493  throw Error(
494  "Could not locate the " UGRID_FACE_COORDINATES " variable named '"
495  + faceCoordinateName + "'! "
496  + "The mesh_topology variable is named "
497  + meshTopology->name());
498 
499  libdap::Array *newFaceCoordArray = dynamic_cast<libdap::Array*>(btp);
500  if (newFaceCoordArray == 0) {
501  throw Error(malformed_expr,
502  "Face coordinate variable '" + faceCoordinateName
503  + "' is not an Array type. It's an instance of "
504  + btp->type_name());
505  }
506 
507 
508  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Found face coordinate array '"<< faceCoordinateName << "'." << endl);
509 
510  // Coordinate arrays MUST be single dimensioned.
511  if(newFaceCoordArray->dimensions(true) != 1){
512  throw Error(malformed_expr,
513  "Face coordinate variable '" + faceCoordinateName
514  + "' has more than one dimension. That's just not allowed. It has "
515  + long_to_string(newFaceCoordArray->dimensions(true)) + " dimensions."
516  );
517  }
518  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face coordinate array '"<< faceCoordinateName << "' has a single dimension." << endl);
519 
520 
521  // Make sure this node coordinate variable has the same size and name as all the others on the list - error if not true.
522  string dimName = newFaceCoordArray->dimension_name(newFaceCoordArray->dim_begin());
523  int dimSize = newFaceCoordArray->dimension_size(newFaceCoordArray->dim_begin(), true);
524 
525  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - dimName: '"<< dimName << "' dimSize: "
526  << libdap::long_to_string(dimSize) << endl);
527 
528  if(faceDimensionName.empty()){
529  faceDimensionName = dimName;
530  }
531  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - faceDimensionName: '"<< faceDimensionName << "' " << endl);
532 
533  if(faceDimensionName.compare(dimName)!=0 ){
534  throw Error(
535  "The face coordinate array '" + faceCoordinateName
536  + "' has the named dimension '"+ dimName + "' which differs from the expected dimension name '"+faceDimensionName
537  +"'. The mesh_topology variable is named " + meshTopology->name());
538  }
539  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face dimension names match." << endl);
540 
541 
542  if(faceCount == 0){
543  faceCount = dimSize;
544  }
545  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - faceCount: "<< libdap::long_to_string(faceCount) << endl);
546 
547  if(faceCount!=dimSize){
548  throw Error(
549  "The face coordinate array '" + faceCoordinateName
550  + "' has a dimension size of " + libdap::long_to_string(dimSize) + " which differs from the the expected size of "
551  + libdap::long_to_string(faceCount) + " The mesh_topology variable is named "
552  + meshTopology->name());
553  }
554  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face counts match." << endl);
555 
556  // Add variable to faceCoordinateArrays.
557  faceCoordinateArrays->push_back(newFaceCoordArray);
558  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Face coordinate array '"<< faceCoordinateName << "' ingested." << endl);
559  }
560  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - Located "<< libdap::long_to_string(faceCoordinateArrays->size()) << " face coordinate arrays." << endl);
561 
562  } else {
563  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - No Face Coordinates Found." << endl);
564  }
565 
566 
567  BESDEBUG("ugrid", "TwoDMeshTopology::ingestFaceCoordinateArrays() - DONE" << endl);
568 
569 }
570 
571 
577 void TwoDMeshTopology::ingestNodeCoordinateArrays(libdap::BaseType *meshTopology, libdap::DDS *dds)
578 {
579  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - BEGIN Gathering node coordinate arrays..." << endl);
580 
581  string node_coordinates;
582  AttrTable at = meshTopology->get_attr_table();
583 
584  AttrTable::Attr_iter iter_nodeCoors = at.simple_find(UGRID_NODE_COORDINATES);
585  if (iter_nodeCoors != at.attr_end()) {
586  node_coordinates = at.get_attr(iter_nodeCoors, 0);
587  } else {
588  throw Error(
589  "Could not locate the " UGRID_NODE_COORDINATES " attribute in the " UGRID_MESH_TOPOLOGY
590  " variable! The mesh_topology variable is named " + meshTopology->name());
591  }
592  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Located '"<< UGRID_NODE_COORDINATES << "' attribute." << endl);
593 
594  if(nodeCoordinateArrays==0)
595  nodeCoordinateArrays = new vector<libdap::Array *>();
596 
597  nodeCoordinateArrays->clear();
598 
599  // Split the node_coordinates string up on spaces
600  // TODO make this work on situations where multiple spaces in the node_coorindates string doesn't hose the split()
601  vector<string> nodeCoordinateNames = split(node_coordinates, ' ');
602 
603  // Find each variable in the resulting list
604  vector<string>::iterator coorName_it;
605  for (coorName_it = nodeCoordinateNames.begin();
606  coorName_it != nodeCoordinateNames.end(); ++coorName_it) {
607  string nodeCoordinateName = *coorName_it;
608 
609  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Processing node coordinate '"<< nodeCoordinateName << "'." << endl);
610 
611  //Now that we have the name of the coordinate variable get it from the DDS!!
612  BaseType *btp = dds->var(nodeCoordinateName);
613  if (btp == 0)
614  throw Error(
615  "Could not locate the " UGRID_NODE_COORDINATES " variable named '"
616  + nodeCoordinateName + "'! "
617  + "The mesh_topology variable is named "
618  + meshTopology->name());
619 
620  libdap::Array *newNodeCoordArray = dynamic_cast<libdap::Array*>(btp);
621  if (newNodeCoordArray == 0) {
622  throw Error(malformed_expr,
623  "Node coordinate variable '" + nodeCoordinateName
624  + "' is not an Array type. It's an instance of "
625  + btp->type_name());
626  }
627  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Found node coordinate array '"<< nodeCoordinateName << "'." << endl);
628 
629  // Coordinate arrays MUST be single dimensioned.
630  if(newNodeCoordArray->dimensions(true) != 1){
631  throw Error(malformed_expr,
632  "Node coordinate variable '" + nodeCoordinateName
633  + "' has more than one dimension. That's just not allowed. It has "
634  + long_to_string(newNodeCoordArray->dimensions(true)) + " dimensions."
635  );
636  }
637  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Node coordinate array '"<< nodeCoordinateName << "' has a single dimension." << endl);
638 
639 
640  // Make sure this node coordinate variable has the same size and name as all the others on the list - error if not true.
641  string dimName = newNodeCoordArray->dimension_name(newNodeCoordArray->dim_begin());
642  int dimSize = newNodeCoordArray->dimension_size(newNodeCoordArray->dim_begin(), true);
643 
644  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - dimName: '"<< dimName << "' dimSize: "
645  << libdap::long_to_string(dimSize) << endl);
646 
647  if(nodeDimensionName.empty()){
648  nodeDimensionName = dimName;
649  }
650  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - nodeDimensionName: '"<< nodeDimensionName << "' " << endl);
651  if(nodeDimensionName.compare(dimName)!=0 ){
652  throw Error(
653  "The node coordinate array '" + nodeCoordinateName
654  + "' has the named dimension '"+ dimName + "' which differs from the expected dimension name '"+nodeDimensionName
655  +"'. The mesh_topology variable is named " + meshTopology->name());
656  }
657  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Node dimension names match." << endl);
658 
659 
660  if(nodeCount == 0){
661  nodeCount = dimSize;
662  }
663  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - nodeCount: "<< libdap::long_to_string(nodeCount) << endl);
664 
665  if(nodeCount!=dimSize){
666  throw Error(
667  "The node coordinate array '" + nodeCoordinateName
668  + "' has a dimension size of " + libdap::long_to_string(dimSize) + " which differs from the the expected size of "
669  + libdap::long_to_string(nodeCount) + " The mesh_topology variable is named "
670  + meshTopology->name());
671  }
672  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Node counts match." << endl);
673 
674 
675  // Add variable to nodeCoordinateArrays vector.
676  nodeCoordinateArrays->push_back(newNodeCoordArray);
677  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - Coordinate array '"<< nodeCoordinateName << "' ingested." << endl);
678 
679 
680  }
681 
682  BESDEBUG("ugrid", "TwoDMeshTopology::ingestNodeCoordinateArrays() - DONE" << endl);
683 
684 }
685 
686 
687 
689 
690  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Building GridFields objects for mesh_topology variable "<< getMeshVariable()->name() << endl);
691  // Start building the Grid for the GridField operation.
692  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Constructing new GF::Grid for "<< name() << endl);
693  gridTopology = new GF::Grid(name());
694 
695  // 1) Make the implicit nodes - same size as the node coordinate arrays
696  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Building and adding implicit range Nodes to the GF::Grid" << endl);
697  GF::AbstractCellArray *nodes = new GF::Implicit0Cells(nodeCount);
698  // Attach the implicit nodes to the grid at rank 0
699  gridTopology->setKCells(nodes, node);
700 
701  // @TODO Do I need to add implicit k-cells for faces (rank 2) if I plan to add range data on faces later?
702  // Apparently not...
703 
704 
705  // Attach the Mesh to the grid.
706  // Get the face node connectivity cells (i think these correspond to the GridFields K cells of Rank 2)
707  // FIXME Read this array once! It is read again below..
708  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Building face node connectivity Cell array from the DAP version" << endl);
709  GF::CellArray *faceNodeConnectivityCells = getFaceNodeConnectivityCells();
710 
711 
712  // Attach the Mesh to the grid at rank 2
713  // This 2 stands for rank 2, or faces.
714  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Attaching Cell array to GF::Grid" << endl);
715  gridTopology->setKCells(faceNodeConnectivityCells, face);
716 
717  // The Grid is complete. Now we make a GridField from the Grid
718  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Construct new GF::GridField from GF::Grid" << endl);
719  d_inputGridField = new GF::GridField(gridTopology);
720  // TODO Question for Bill: Can we delete the GF::Grid (tdmt->gridTopology) here?
721 
722 
723  // We read and add the coordinate data (using GridField->addAttribute() to the GridField at
724  // grid dimension/rank/dimension 0 (a.k.a. node)
725  vector<libdap::Array *>::iterator ncit;
726  for (ncit = nodeCoordinateArrays->begin(); ncit != nodeCoordinateArrays->end(); ++ncit) {
727  libdap::Array *nca = *ncit;
728  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Adding node coordinate "<< nca->name() << " to GF::GridField at rank 0" << endl);
729  GF::Array *gfa = extractGridFieldArray(nca,sharedIntArrays,sharedFloatArrays);
730  gfArrays.push_back(gfa);
731  d_inputGridField->AddAttribute(node, gfa);
732  }
733 
734  // We read and add the coordinate data (using GridField->addAttribute() to the GridField at
735  // grid dimension/rank/dimension 0 (a.k.a. node)
736  for (ncit = faceCoordinateArrays->begin(); ncit != faceCoordinateArrays->end(); ++ncit) {
737  libdap::Array *fca = *ncit;
738  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Adding face coordinate "<< fca->name() << " to GF::GridField at rank " << face << endl);
739  GF::Array *gfa = extractGridFieldArray(fca,sharedIntArrays,sharedFloatArrays);
740  gfArrays.push_back(gfa);
741  d_inputGridField->AddAttribute(face, gfa);
742  }
743 }
744 
745 
746 
747 #if 0 // Unused crufty code from previous efforts
748 
751 void TwoDMeshTopology::buildRestrictedGfTopology(locationType loc, string filterExpression){
753  applyRestrictOperator(loc, filterExpression);
754 }
755 #endif
756 
758  return resultGridField->Size(dim);
759 }
760 
761 
762 
763 #if 0 // Unused crufty code from previous efforts
764 
767 void TwoDMeshTopology::buildGridFieldsTopology()
768 {
769 
771 
772  // For each range data variable associated with this MeshTopology read and add the data to the GridField
773  // At the appropriate rank.
774  for (vector<MeshDataVariable *>::iterator mdv_it = rangeDataArrays->begin(); mdv_it != rangeDataArrays->end(); ++mdv_it) {
775  MeshDataVariable *mdVar = *mdv_it;
776  GF::Array *gfa = extractGridFieldArray(mdVar->getDapArray(),sharedIntArrays,sharedFloatArrays);
777  BESDEBUG("ugrid", "TwoDMeshTopology::buildGridFieldsTopology() - Adding mesh data variable '"<< mdVar->getName() <<"' to GF::GridField at rank "<< mdVar->getGridLocation() << endl);
778  d_inputGridField->AddAttribute(mdVar->getGridLocation(), gfa);
779  gfArrays.push_back(gfa);
780  }
781 
782 }
783 #endif
784 
785 
786 
797 GF::Node *TwoDMeshTopology::getFncArrayAsGFCells(libdap::Array *fncVar)
798 {
799 
800  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFNodes() - BEGIN" << endl);
801 
802 
803  int nodesPerFace = fncVar->dimension_size(fncNodesDim,true);
804  int faceCount = fncVar->dimension_size(fncFacesDim,true);
805 
806 
807  // interpret the array data as triangles
808  GF::Node *cells = new GF::Node[ faceCount * nodesPerFace ];
809 
810  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFNodes() - Reading DAP data into GF::Node array." << endl);
811  GF::Node *temp_nodes = ugrid::extractArray<GF::Node>(fncVar);
812 
813  // Reorganize the cell ids so that cellids contains
814  // the cells in three consecutive values (0,1,2; 3,4,5; ...).
815  // The the values from fncVar now in cellids2 and ar organized
816  // as 0,N,2N; 1,1+N,1+2N; ...
817 
818  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFNodes() - Re-packing and copying GF::Node array to result." << endl);
819  for (int fIndex = 0; fIndex < faceCount; fIndex++) {
820  for(int nIndex=0; nIndex<nodesPerFace ; nIndex++){
821  cells[nodesPerFace * fIndex + nIndex] = temp_nodes[fIndex + (faceCount * nIndex)];
822 
823  }
824  }
825 
826  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFNodes() - Deleting intermediate GF::Node array." << endl);
827  delete [] temp_nodes;
828 
829  BESDEBUG("ugrid", "TwoDMeshTopology::getFncArrayAsGFNodes() - DONE" << endl);
830 
831  return cells;
832 }
833 
838 int TwoDMeshTopology::getStartIndex(libdap::Array *array)
839 {
840  AttrTable &at = array->get_attr_table();
841  AttrTable::Attr_iter start_index_iter = at.simple_find(UGRID_START_INDEX);
842  if (start_index_iter != at.attr_end()) {
843  BESDEBUG("ugrid", "TwoDMeshTopology::getStartIndex() - Found the "<< UGRID_START_INDEX <<" attribute." << endl);
844  AttrTable::entry *start_index_entry = *start_index_iter;
845  if (start_index_entry->attr->size() == 1) {
846  string val = (*start_index_entry->attr)[0];
847  BESDEBUG("TwoDMeshTopology::getStartIndex", "value: " << val << endl);
848  stringstream buffer(val);
849  // what happens if string cannot be converted to an integer?
850  int start_index;
851  ;
852 
853  buffer >> start_index;
854  return start_index;
855  } else {
856  throw Error(malformed_expr,
857  "Index origin attribute exists, but either no value supplied, or more than one value supplied.");
858  }
859  }
860  return 0;
861 }
862 
866 GF::CellArray *TwoDMeshTopology::getFaceNodeConnectivityCells()
867 {
868 
869  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityCells() - Building face node connectivity Cell " <<
870  "array from the Array "<< faceNodeConnectivityArray->name() << endl);
871 
872  int nodesPerFace = faceNodeConnectivityArray->dimension_size(fncNodesDim);
873  int total_size = nodesPerFace * faceCount;
874 
875  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityCells() - Converting FNCArray to GF::Node array." << endl);
876  fncCellArray = getFncArrayAsGFCells(faceNodeConnectivityArray);
877 
878 
879  // adjust for the start_index (cardinal or ordinal array access)
880  int startIndex = getStartIndex(faceNodeConnectivityArray);
881  if (startIndex != 0) {
882  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityCells() - Applying startIndex to GF::Node array." << endl);
883  for (int j = 0; j < total_size; j++) {
884  fncCellArray[j] -= startIndex;
885  }
886  }
887  // Create the cell array
888  GF::CellArray *rankTwoCells = new GF::CellArray(fncCellArray, faceCount, nodesPerFace);
889 
890 
891  BESDEBUG("ugrid", "TwoDMeshTopology::getFaceNodeConnectivityCells() - DONE" << endl);
892  return rankTwoCells;
893 
894 }
895 
896 
897 
898 void TwoDMeshTopology::applyRestrictOperator(locationType loc, string filterExpression)
899 {
900 
901  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - BEGIN" << endl);
902 
903  // I think this function could be done with just the following single line:
904  // resultGridField = GF::RefRestrictOp::Restrict(filterExpression,loc,d_inputGridField);
905 
906  // Build the restriction operator;
907  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - Constructing new GF::RestrictOp using user "<<
908  "supplied 'dimension' value and filter expression combined with the GF:GridField " << endl);
909  GF::RestrictOp op = GF::RestrictOp(filterExpression, loc, d_inputGridField);
910 
911 
912  // Apply the operator and get the result;
913  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - Applying GridField operator." << endl);
914  GF::GridField *resultGF = op.getResult();
915  resultGridField = resultGF;
916  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - GridField operator applied and result obtained." << endl);
917  BESDEBUG("ugrid", "TwoDMeshTopology::applyRestrictOperator() - END" << endl);
918 }
919 
920 
921 #if 0 // Unused crufty code from previous efforts
922 
925 static void rDAWorker(
926  libdap::Array *dapArray,
927  libdap::Array::Dim_iter thisDim,
928  libdap::Array::Dim_iter locationCoordinateDim,
929  locationType gridLocation,
930  NDimensionalArray *results) {
931 
932  if(thisDim==dapArray->dim_end()){
933 
934  // GF::Array *gfa = extractGridFieldArray(dapArray,sharedIntArrays,sharedFloatArrays);
935 
936  //resultGridField->Bind(gridLocation,gfa);
937  //resultGridField->GetGrid()->normalize();
938 
939  //libdap::Array *resultRangeVar = getGFAttributeAsDapArray( dapArray, gridLocation, resultGridField);
940 
941  //resultGridField->unBind(gridLocation,gfa);
942 
943 
944  }
945  else {
946  libdap::Array::Dim_iter nextDim = thisDim;
947  nextDim++; // now it's the next dimension.
948 
949  if(thisDim==locationCoordinateDim){
950 
951  if(nextDim!=dapArray->dim_end()){
952  string msg = "rDAWorker() - The location coordinate dimension is not the last dimension in the array. Hyperslab subsetting of this dimension is not supported.";
953  BESDEBUG("ugrid", msg << endl);
954  throw Error(malformed_expr, msg);
955 
956  }
957 
958 
959  rDAWorker(dapArray, nextDim, locationCoordinateDim, gridLocation, results);
960  }
961  else {
962 
963 
964  unsigned int start = dapArray->dimension_start(thisDim, true);
965  unsigned int stride = dapArray->dimension_stride(thisDim, true);
966  unsigned int stop = dapArray->dimension_stop(thisDim, true);
967 
968  for(unsigned int dimIndex=start; dimIndex<=stop ; dimIndex+=stride){
969  dapArray->add_constraint(thisDim,dimIndex,1,dimIndex);
970  rDAWorker(dapArray, nextDim, locationCoordinateDim, gridLocation, results);
971  }
972 
973  // Reset the constraint for this dimension.
974  dapArray->add_constraint(thisDim,start,stride,stop);
975  }
976 
977  }
978 
979 
980 }
981 
989 libdap::Array *TwoDMeshTopology::restrictDapArray(libdap::Array *dapArray, libdap::Array::Dim_iter locationCoordinateDim, locationType gridLocation){
990 
991 
992  // We want the manipulate the Array's Dimensions so that only a single dimensioned slab of the location coordinate dimension
993  // is read at a time. We need to cache the original constrained dimensions so that we can build the correct collection of
994  // location coordinate dimension slabs.
995 
996  // What's the shape of the requested range variable array?
997  vector<unsigned int> arrayShape(dapArray->dimensions(true));
998  NDimensionalArray::computeConstrainedShape(dapArray, &arrayShape );
999 
1000  // Now,
1001  arrayShape[dapArray->dimensions(true) - 1] = resultGridField->GetGrid()->countKCells(node);
1002 
1003  NDimensionalArray *result = new NDimensionalArray(&arrayShape, dapArray->var()->type());
1004 
1005  rDAWorker(dapArray, dapArray->dim_begin(), locationCoordinateDim, gridLocation, result);
1006 
1007  libdap::Array *myResult = 0;
1008 
1009  return myResult;
1010 
1011 }
1012 #endif
1013 
1014 #if 0 // Unused crufty code from previous efforts
1015 void TwoDMeshTopology::restrictRange( vector<BaseType *> *results){
1016 
1017 
1018  BESDEBUG("ugrid", "TwoDMeshTopology::restrictRange() - BEGIN" << endl);
1019 
1020 
1021  BESDEBUG("ugrid", "TwoDMeshTopology::restrictRange() - Normalizing Grid." << endl);
1022  resultGridField->GetGrid()->normalize();
1023 
1024  // Add the coordinate node arrays to the response.
1025  BESDEBUG("ugrid", "TwoDMeshTopology::restrictRange() - Adding the coordinate node arrays to the response." << endl);
1026  vector<libdap::Array *>::iterator it;
1027  for (it = nodeCoordinateArrays->begin(); it != nodeCoordinateArrays->end(); ++it) {
1028  libdap::Array *sourceCoordinateArray = *it;
1029  libdap::Array *resultCoordinateArray = getGFAttributeAsDapArray(sourceCoordinateArray, node, resultGridField );
1030  results->push_back(resultCoordinateArray);
1031  }
1032  // Add the new face node connectivity array - make sure it has the same attributes as the original.
1033  BESDEBUG("ugrid", "TwoDMeshTopology::restrictRange() - Adding the new face node connectivity array to the response." << endl);
1034  libdap::Array *resultFaceNodeConnectivityDapArray = getGridFieldCellArrayAsDapArray(resultGridField, faceNodeConnectivityArray);
1035  results->push_back(resultFaceNodeConnectivityDapArray);
1036 
1037 
1038  // Add the range variable data arrays to the response.
1039  BESDEBUG("ugrid", "TwoDMeshTopology::restrictRange() - Processing range variable data arrays." << endl);
1040  vector<MeshDataVariable *>::iterator mdvIt;
1041  for (mdvIt = rangeDataArrays->begin(); mdvIt != rangeDataArrays->end(); ++mdvIt) {
1042  MeshDataVariable *mdv = *mdvIt;
1043 
1044 
1045  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldToDapObjects() - Processing MeshDataVariable '"<< mdv->getName() << "' array type is "<< mdv->getDapArray()->type_name() << endl);
1046 
1047  libdap::Array *resultRangeVar = restrictDapArray(mdv->getDapArray(), mdv->getLocationCoordinateDimension(), mdv->getGridLocation());
1048 
1049  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldToDapObjects() - Retrieved restricted DAP Array for MeshDataVariable '"<< mdv->getName() << "'" << endl);
1050  results->push_back(resultRangeVar);
1051  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldToDapObjects() - DAP Array added to result" << endl);
1052 
1053  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldToDapObjects() - DAP Array added to result" << endl);
1054  }
1055 
1056 
1057 
1058  BESDEBUG("ugrid", "TwoDMeshTopology::restrictRange() - END" << endl);
1059 
1060 }
1061 #endif
1062 
1063 
1064 
1066 {
1067  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - BEGIN" << endl);
1068 
1069  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Normalizing Grid." << endl);
1070  resultGridField->GetGrid()->normalize();
1071 
1072  // Add the node coordinate arrays to the results.
1073  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Converting the node coordinate arrays to DAP arrays." << endl);
1074  vector<libdap::Array *>::iterator it;
1075  for (it = nodeCoordinateArrays->begin(); it != nodeCoordinateArrays->end(); ++it) {
1076  libdap::Array *sourceCoordinateArray = *it;
1077  libdap::Array *resultCoordinateArray = getGFAttributeAsDapArray(sourceCoordinateArray, node, resultGridField );
1078  results->push_back(resultCoordinateArray);
1079  }
1080 
1081 #if 1
1082  // Add the face coordinate arrays to the results.
1083  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Converting the face coordinate arrays to DAP arrays." << endl);
1084  for (it = faceCoordinateArrays->begin(); it != faceCoordinateArrays->end(); ++it) {
1085  libdap::Array *sourceCoordinateArray = *it;
1086  libdap::Array *resultCoordinateArray = getGFAttributeAsDapArray(sourceCoordinateArray, face, resultGridField );
1087  results->push_back(resultCoordinateArray);
1088  }
1089 #endif
1090 
1091 
1092  // Add the new face node connectivity array - make sure it has the same attributes as the original.
1093  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - Adding the new face node connectivity array to the response." << endl);
1094  libdap::Array *resultFaceNodeConnectivityDapArray = getGridFieldCellArrayAsDapArray(resultGridField, faceNodeConnectivityArray);
1095  results->push_back(resultFaceNodeConnectivityDapArray);
1096 
1097  results->push_back(getMeshVariable()->ptr_duplicate());
1098 
1099  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldStructureToDapObjects() - END" << endl);
1100 }
1101 
1102 
1103 
1104 #if 0 // Unused crufty code from previous efforts
1105 
1108 void TwoDMeshTopology::convertResultRangeVarsToDapObjects(vector<BaseType *> *results)
1109 {
1110 
1111  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultRangeVarsToDapObjects() - BEGIN" << endl);
1112 
1113 
1114  // Add the range variable data arrays to the response.
1115  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultRangeVarsToDapObjects() - Adding the range variable data arrays to the response." << endl);
1116  vector<MeshDataVariable *>::iterator mdvIt;
1117  for (mdvIt = rangeDataArrays->begin(); mdvIt != rangeDataArrays->end(); ++mdvIt) {
1118  MeshDataVariable *mdv = *mdvIt;
1119  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultRangeVarsToDapObjects() - Processing MeshDataVariable '"<< mdv->getName() << "' array type is "<< mdv->getDapArray()->type_name() << endl);
1120  libdap::Array *resultRangeVar = getGFAttributeAsDapArray( mdv->getDapArray(), mdv->getGridLocation(), resultGridField);
1121  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultRangeVarsToDapObjects() - Retrieved restricted DAP Array for MeshDataVariable '"<< mdv->getName() << "'" << endl);
1122  results->push_back(resultRangeVar);
1123  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultRangeVarsToDapObjects() - DAP Array added to result" << endl);
1124  }
1125 
1126  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultRangeVarsToDapObjects() - END" << endl);
1127 
1128 }
1129 #endif
1130 
1131 #if 0 // Unused crufty code from previous efforts
1132 
1135 void TwoDMeshTopology::convertResultGridFieldToDapObjects(vector<BaseType *> *results)
1136 {
1137 
1138  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldToDapObjects() - BEGIN" << endl);
1139 
1141  convertResultRangeVarsToDapObjects(results);
1142 
1143  BESDEBUG("ugrid", "TwoDMeshTopology::convertResultGridFieldToDapObjects() - END" << endl);
1144 
1145 }
1146 #endif
1147 
1148 
1149 
1155 libdap::Array *TwoDMeshTopology::getNewFncDapArray(libdap::Array *templateArray, int N)
1156 {
1157 
1158  // Is the template array a 2D array?
1159  int dimCount = templateArray->dimensions(true);
1160  if (dimCount != 2)
1161  throw Error(
1162  "Expected a 2 dimensional array. The array '"
1163  + templateArray->name() + "' has "
1164  + libdap::long_to_string(dimCount) + " dimensions.");
1165 
1166  // Is the template array really 3xN?
1167  libdap::Array::Dim_iter di = templateArray->dim_begin();
1168  if (di->c_size != 3) {
1169  string msg =
1170  "Expected a 2 dimensional array with shape of 3xN! The array '"
1171  + templateArray->name() + "' has a first "
1172  + "dimension of size " + libdap::long_to_string(di->c_size);
1173  BESDEBUG("ugrid", "TwoDMeshTopology::getNewFcnDapArray() - "<< msg << endl);
1174  templateArray->print_decl(cerr, " ", true,false, true);
1175 
1176  throw Error(malformed_expr, msg);
1177  }
1178 
1179  // Get a new template variable for our new array (should be just like the template for the source array)
1180  //BaseType *arrayTemplate = getDapVariableInstance(templateArray->var(0)->name(),templateArray->var(0)->type());
1181  Int32 tmplt(templateArray->name());
1182  libdap::Array *newArray = new libdap::Array(templateArray->name(),&tmplt);
1183 
1184  //Add the first dimension (size 3 same same as template array's first dimension)
1185  newArray->append_dim(3, di->name);
1186 
1187  // Add the second dimension to the result array, but use only the name from the template array's
1188  // second dimension. The size will be from the passed parameter N
1189  di++;
1190  newArray->append_dim(N, di->name);
1191 
1192  // Copy the attributes of the template array to our new array.
1193  newArray->set_attr_table(templateArray->get_attr_table());
1194 
1195  // make the new array big enough to hold all the values.
1196  newArray->reserve_value_capacity(3 * N);
1197 
1198 #if 0 // This is debugging - keep it around we might need it...
1199  cerr<<"getNewFcnDapArray() -"<<endl<<endl;
1200  cerr << "Newly minted Array: "<< endl;
1201  newArray->print_val(cerr);
1202  cerr<<endl<<endl;
1203 
1204 #endif
1205 
1206  return newArray;
1207 
1208 }
1209 
1215 libdap::Array *TwoDMeshTopology::getGridFieldCellArrayAsDapArray(GF::GridField *resultGridField, libdap::Array *sourceFcnArray)
1216 {
1217 
1218  BESDEBUG("ugrid", "TwoDMeshTopology::getGridFieldCellArrayAsDapArray() - BEGIN" << endl);
1219 
1220  // Get the rank 2 k-cells from the GridField object.
1221  GF::CellArray* gfCellArray = (GF::CellArray*) (resultGridField->GetGrid()->getKCells(2));
1222 
1223  // This is a vector of size N holding vectors of size 3
1224  vector<vector<int> > nodes2 = gfCellArray->makeArrayInts();
1225 
1226  libdap::Array *resultFncDapArray = getNewFncDapArray(sourceFcnArray, nodes2.size());
1227 
1228  // Make a vector to hold the re-packed cell nodes.
1229  vector<dods_int32> rowMajorNodes;
1230 
1231 
1232  // adjust for the start_index (cardinal or ordinal array access)
1233  int startIndex = getStartIndex(sourceFcnArray);
1234  if (startIndex != 0) {
1235  BESDEBUG("ugrid", "TwoDMeshTopology::getGridFieldCellArrayAsDapArray() - Updating locations to match source FaceNodeConnectivity array." << endl);
1236  }
1237 
1238 
1239  // Re-pack the mesh nodes.
1240  for (unsigned int firstDim = 0; firstDim < 3; firstDim++) {
1241  for (unsigned int secondDim = 0; secondDim < nodes2.size();
1242  secondDim++) {
1243  dods_int32 val = nodes2.at(secondDim).at(firstDim) + startIndex;
1244  rowMajorNodes.push_back(val);
1245  }
1246  }
1247 
1248 #if 0 // This is debugging - keep it around we might need it.
1249  cerr << "getGridFieldCellArrayAsDapArray() - rowMajorNodes: " << endl << "{";
1250  for (unsigned int j=0; j < rowMajorNodes.size(); j++) {
1251  dods_int32 val = rowMajorNodes.at(j);
1252  cerr << val << ", ";
1253  }
1254  cerr << "}" << endl;
1255 #endif
1256 
1257  // Add them to the DAP array.
1258  resultFncDapArray->set_value(rowMajorNodes, rowMajorNodes.size());
1259 
1260 #if 0 // This is debugging - keep it around we might need it.
1261  cerr << "getGridFieldCellArrayAsDapArray() - DAP Array: "<< endl;
1262  resultFcnDapArray->print_val(cerr);
1263 #endif
1264  BESDEBUG("ugrid", "TwoDMeshTopology::getGridFieldCellArrayAsDapArray() - DONE" << endl);
1265 
1266  return resultFncDapArray;
1267 
1268 }
1269 
1270 
1271 
1272 static libdap::Array::Dim_iter copySizeOneDimensions(libdap::Array *sourceArray, libdap::Array *dapArray){
1273 
1274  libdap::Array::Dim_iter dataDimension;
1275  for (libdap::Array::Dim_iter srcArrIt = sourceArray->dim_begin(); srcArrIt != sourceArray->dim_end() ;++srcArrIt) {
1276 
1277  // Get the original dimension size
1278  int dimSize = sourceArray->dimension_size(srcArrIt,true);
1279  string dimName = sourceArray->dimension_name(srcArrIt);
1280 
1281  // Preserve single dimensions
1282  if(dimSize==1){
1283  BESDEBUG("ugrid", "TwoDMeshTopology::copySizeOneDimensions() - Adding size one dimension '"<< dimName << "' from source array into result." << endl);
1284  dapArray->append_dim(dimSize, dimName);
1285  }
1286  else {
1287  BESDEBUG("ugrid", "TwoDMeshTopology::copySizeOneDimensions() - Located data dimension '"<< dimName << "' in source array." << endl);
1288 
1289  dataDimension = srcArrIt;
1290  }
1291  }
1292  BESDEBUG("ugrid", "TwoDMeshTopology::copySizeOneDimensions() - Returning dimension iterator pointing to '"<< sourceArray->dimension_name(dataDimension) << "'." << endl);
1293  return dataDimension;
1294 }
1295 
1296 
1297 
1302 libdap::Array *TwoDMeshTopology::getGFAttributeAsDapArray(libdap::Array *templateArray, locationType rank, GF::GridField *resultGridField)
1303 {
1304 
1305  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - BEGIN" << endl);
1306 
1307  // The result variable is assumed to be bound to the GridField at 'rank'
1308  // Try to get the Attribute from 'rank' with the same name as the source array
1309  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Retrieving GF::GridField Attribute '" <<
1310  templateArray->name() << "'" << endl);
1311  GF::Array* gfa = resultGridField->GetAttribute(rank, templateArray->name());
1312 
1313  libdap::Array *dapArray;
1314  BaseType *templateVar = templateArray->var();
1315  string dimName;
1316 
1317  switch (templateVar->type()) {
1318  case dods_byte_c:
1319  case dods_uint16_c:
1320  case dods_int16_c:
1321  case dods_uint32_c:
1322  case dods_int32_c: {
1323  // Get the data
1324  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Retrieving GF::Array as libdap::Array of libdap::Int32." << endl);
1325  vector<dods_int32> GF_ints = gfa->makeArray();
1326  // Make a DAP array to put the data into.
1327  libdap::Int32 tmpltVar(templateVar->name());
1328  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - libdap::Int32 template variable name: " << tmpltVar.name() << endl);
1329  dapArray = new libdap::Array(templateArray->name(), &tmpltVar);
1330 
1331  // copy the dimensions whose size is "1" from the source array to the result.
1332  libdap::Array::Dim_iter dataDimension = copySizeOneDimensions(templateArray, dapArray);
1333 
1334  // Add the result dimension
1335  dimName = templateArray->dimension_name(dataDimension);
1336  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Adding dimension " << dimName << endl);
1337 
1338  dapArray->append_dim(GF_ints.size(), dimName);
1339 
1340  // Add the data
1341  dapArray->set_value(GF_ints, GF_ints.size());
1342  break;
1343  }
1344  case dods_float32_c:
1345  case dods_float64_c: {
1346  // Get the data
1347  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Retrieving GF::Array as libdap::Array of libdap::Float64." << endl);
1348  vector<dods_float64> GF_floats = gfa->makeArrayf();
1349  // Make a DAP array to put the data into.
1350  libdap::Float64 tmpltVar(templateVar->name());
1351  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - libdap::Float64 template variable name: " << tmpltVar.name() << endl);
1352  dapArray = new libdap::Array(templateArray->name(), &tmpltVar);
1353 
1354  // copy the dimensions whose size is "1" from the source array to the result.
1355  libdap::Array::Dim_iter dataDimension = copySizeOneDimensions(templateArray, dapArray);
1356 
1357  // Add the result dimension
1358  dimName = templateArray->dimension_name(dataDimension);
1359  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Adding dimension " << dimName << endl);
1360  dapArray->append_dim(GF_floats.size(), dimName);
1361 
1362  // Add the data
1363  dapArray->set_value(GF_floats, GF_floats.size());
1364  break;
1365  }
1366  default:
1367  throw InternalErr(__FILE__, __LINE__,
1368  "Unknown DAP type encountered when converting to gridfields array");
1369  }
1370 
1371  // Copy the source objects attributes.
1372  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - Copying libdap::Attribute's from template array " << templateArray->name() << endl);
1373  dapArray->set_attr_table(templateArray->get_attr_table());
1374 
1375  BESDEBUG("ugrid", "TwoDMeshTopology::getGFAttributeAsDapArray() - DONE" << endl);
1376 
1377  return dapArray;
1378 }
1379 
1380 
1381 
1382 
1383 #if 0
1384 
1388 void TwoDMeshTopology::getResultGFAttributeValues(libdap::Array *templateArray, locationType rank, void *target){
1389 
1390  getResultGFAttributeValues(templateArray->name(),templateArray->var()->type(),rank,target);
1391 }
1392 #endif
1393 
1394 
1399 void TwoDMeshTopology::getResultGFAttributeValues(string attrName, libdap::Type dapType, locationType rank, void *target)
1400 {
1401  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - BEGIN" << endl);
1402 
1403  // The result variable is assumed to be bound to the GridField at 'rank'
1404  // Try to get the Attribute from 'rank' with the same name as the source array
1405  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieving GF::GridField Attribute '" <<
1406  attrName << "'" << endl);
1407 
1408  GF::Array *gfa = 0;
1409  if(resultGridField->IsAttribute(rank,attrName)){
1410  gfa = resultGridField->GetAttribute(rank, attrName);
1411  }
1412  else {
1413  string msg = "Oddly, the requested attribute '"+ attrName +"' associated with rank "+ libdap::long_to_string(rank) +
1414  " does not appear in the resultGridField object! \n" +
1415  "resultGridField->MaxRank(): " + libdap::long_to_string(resultGridField->MaxRank());
1416 
1417  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - ERROR! " << msg << endl);
1418  throw InternalErr(__FILE__, __LINE__,
1419  "ERROR - Unable to locate requested GridField attribute. "+ msg);
1420  }
1421 
1422  switch (dapType) {
1423  case dods_byte_c:
1424  case dods_uint16_c:
1425  case dods_int16_c:
1426  case dods_uint32_c:
1427  case dods_int32_c: {
1428  // Get the data
1429  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieving GF::Array as libdap::Array of libdap::Int32." << endl);
1430  vector<dods_int32> GF_ints = gfa->makeArray();
1431  stringstream s;
1432  for(unsigned int i=0; i<GF_ints.size() ; i++){
1433  s << "GF_ints[" << i << "]: " << GF_ints[i] << endl;
1434  }
1435  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieved GF_ints: "<< endl << s.str());
1436 
1437 
1438 
1439 
1440  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Copying GF result to target memory" << endl);
1441  memcpy(target, &GF_ints[0], GF_ints.size()*sizeof(dods_int32));
1442  break;
1443  }
1444  case dods_float32_c:
1445  case dods_float64_c: {
1446  // Get the data
1447  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieving GF::Array as libdap::Array of libdap::Float64." << endl);
1448  vector<dods_float64> GF_floats = gfa->makeArrayf();
1449 
1450  stringstream s;
1451  for(unsigned int i=0; i<GF_floats.size() ; i++){
1452  s << "GF_ints[" << i << "]: " << GF_floats[i] << endl;
1453  }
1454  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Retrieved GF_floats: "<< endl << s.str());
1455  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - Copying GF result to target memory" << endl);
1456  memcpy(target, &GF_floats[0], GF_floats.size()*sizeof(dods_float64));
1457  break;
1458 
1459  }
1460  default:
1461  throw InternalErr(__FILE__, __LINE__,
1462  "Unknown DAP type encountered when converting to gridfields result values");
1463  }
1464  BESDEBUG("ugrid", "TwoDMeshTopology::getResultGFAttributeValues() - END" << endl);
1465 
1466 }
1467 
1468 
1469 static string getIndexVariableName(locationType location){
1470  switch(location){
1471 
1472  case node:
1473  {
1474  return "node_index";
1475  }
1476  break;
1477 
1478  case face:
1479  {
1480  return "face_index";
1481  }
1482  break;
1483 
1484  case edge:
1485  default:
1486  break;
1487  }
1488  string msg = "TwoDMeshTopology::setLocationCoordinateDimension() - Unknown/Unsupported location value '" +
1489  libdap::long_to_string(location) + "'";
1490  BESDEBUG("ugrid", msg << endl);
1491  throw Error( msg );
1492 
1493 }
1494 
1495 
1497 
1498  switch(location){
1499  case node:
1500  {
1501  return nodeCount;
1502  }
1503  break;
1504 
1505  case face:
1506  {
1507  return faceCount;
1508  }
1509  break;
1510 
1511  case edge:
1512  default:
1513  break;
1514  }
1515 
1516  string msg = "TwoDMeshTopology::setLocationCoordinateDimension() - Unknown/Unsupported location value '" +
1517  libdap::long_to_string(location) + "'";
1518  BESDEBUG("ugrid", msg << endl);
1519  throw Error( msg );
1520 
1521 }
1522 
1523 
1528 {
1529 
1530  int size = getInputGridSize(location);
1531  string name = getIndexVariableName(location);
1532 
1533  BESDEBUG("ugrid", "TwoDMeshTopology::addIndexVariable() - Adding index variable '" << name <<
1534  "' size: " << libdap::long_to_string(size) << " at rank " << libdap::long_to_string(location) << endl);
1535 
1536  GF::Array *indexArray = newGFIndexArray(name,size,sharedIntArrays);
1537  d_inputGridField->AddAttribute(location, indexArray);
1538  gfArrays.push_back(indexArray);
1539 }
1540 
1541 
1542 
1543 
1548 {
1549  string name = getIndexVariableName(location);
1550  getResultGFAttributeValues(name,dods_int32_c,location,target);
1551 
1552 }
1553 
1554 
1555 
1556 
1557 
1558 } // namespace ugrid
void applyRestrictOperator(locationType loc, string filterExpression)
libdap::BaseType * getMeshVariable()
void getResultIndex(locationType location, void *target)
libdap::Array::Dim_iter getLocationCoordinateDimension()
#define UGRID_FACE_NODE_CONNECTIVITY
Definition: ugrid_utils.h:49
void init(string meshVarName, libdap::DDS *dds)
only call this from the constructor?? Seems like the thing to do, but then the constructor will be po...
void setLocationCoordinateDimension(MeshDataVariable *mdv)
STL namespace.
GF::Array * extractGridFieldArray(libdap::Array *a, vector< int * > *sharedIntArrays, vector< float * > *sharedFloatArrays)
Extract data from a DAP array and return those values in a gridfields array.
Definition: ugrid_utils.cc:287
locationType getGridLocation()
void getResultGFAttributeValues(string attrName, libdap::Type type, locationType rank, void *target)
Retrieves a single dimensional GF attribute array from a GF::GridField and places the data into DAP a...
GF::Array * newGFIndexArray(string name, long size, vector< int * > *sharedIntArrays)
Definition: ugrid_utils.cc:265
static class NCMLUtil overview
libdap::Array * getDapArray()
#define UGRID_MESH_TOPOLOGY
Definition: ugrid_utils.h:47
void addIndexVariable(locationType location)
Adds an index variable at the gridfields rank as indicated by the passed locationType.
int getInputGridSize(locationType location)
long computeConstrainedShape(libdap::Array *a, std::vector< unsigned int > *shape)
Compute the constrained shape of the Array and return it in a vector.
Definition: fojson_utils.cc:61
locationType
Definition: LocationType.h:31
#define UGRID_FACE_COORDINATES
Definition: ugrid_utils.h:66
vector< string > & split(const string &s, char delim, vector< string > &elems)
Splits the string on the passed char.
Definition: ugrid_utils.cc:436
void convertResultGridFieldStructureToDapObjects(vector< libdap::BaseType * > *results)
void setLocationCoordinateDimension(libdap::Array::Dim_iter cdim)
int getResultGridSize(locationType location)
#define UGRID_START_INDEX
Definition: ugrid_utils.h:58
#define BESDEBUG(x, y)
macro used to send debug information to the debug stream
Definition: BESDebug.h:64
#define UGRID_TOPOLOGY_DIMENSION
Definition: ugrid_utils.h:51
string getAttributeValue(BaseType *bt, string aName)
Definition: ugrid_utils.cc:456
#define UGRID_NODE_COORDINATES
Definition: ugrid_utils.h:48