43 #include <curl/curl.h>
49 #include <Structure.h>
56 #include "BESStopWatch.h"
59 #include "ugrid_utils.h"
60 #include "MeshDataVariable.h"
61 #include "TwoDMeshTopology.h"
62 #include "NDimensionalArray.h"
63 #include <gridfields/GFError.h>
65 #include "ugrid_restrict.h"
69 #define BESDEBUG( x, y )
84 struct UgridRestrictArgs {
94 locationType dimension;
102 vector<libdap::Array *> rangeVars;
107 string filterExpression;
115 static void addRangeVar(DDS *dds, libdap::Array *rangeVar, map<
string, vector<MeshDataVariable *> *> *rangeVariables)
117 MeshDataVariable *mdv =
new MeshDataVariable();
119 string meshVarName = mdv->getMeshName();
121 BaseType *meshVar = dds->var(meshVarName);
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);
130 vector<MeshDataVariable *> *requestedRangeVarsForMesh;
131 map<string, vector<MeshDataVariable *> *>::iterator mit = rangeVariables->find(meshVarName);
132 if (mit == rangeVariables->end()) {
135 "addRangeVar() - MeshTopology object for '" << meshVarName <<
"' does not exist. Getting a new one... " << endl);
137 requestedRangeVarsForMesh =
new vector<MeshDataVariable *>();
138 (*rangeVariables)[meshVarName] = requestedRangeVarsForMesh;
143 "addRangeVar() - MeshTopology object for '" << meshVarName <<
"' exists. Retrieving... " << endl);
144 requestedRangeVarsForMesh = mit->second;
147 requestedRangeVarsForMesh->push_back(mdv);
151 string usage(
string fnc){
153 string usage = fnc+
"(rangeVariable:string, [rangeVariable:string, ... ] condition:string)";
161 static UgridRestrictArgs processUgrArgs(
string func_name, locationType dimension,
int argc, BaseType *argv[])
164 BESDEBUG(
"ugrid",
"processUgrArgs() - BEGIN" << endl);
166 UgridRestrictArgs args;
167 args.dimension = dimension;
168 BESDEBUG(
"ugrid",
"args.dimension: " << libdap::long_to_string(args.dimension) << endl);
170 args.rangeVars = vector<libdap::Array *>();
174 throw Error(malformed_expr,
175 "Wrong number of arguments to ugrid restrict function: " + usage(func_name) +
" was passed " + long_to_string(argc)
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());
188 args.dimension = (locationType)
dynamic_cast<Int32&
>(*argv[0]).value();
189 BESDEBUG(
"ugrid",
"args.dimension: " << libdap::long_to_string(args.dimension) << endl);
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 "
200 args.filterExpression =
dynamic_cast<Str&
>(*bt).value();
201 BESDEBUG(
"ugrid",
"args.filterExpression: '" << args.filterExpression <<
"' (AS RECEIVED)" << endl);
203 args.filterExpression = www2id(args.filterExpression);
205 BESDEBUG(
"ugrid",
"args.filterExpression: '" << args.filterExpression <<
"' (URL DECODED)" << endl);
213 for (
int i = 0; i < (argc - 1); 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 "
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());
225 args.rangeVars.push_back(newRangeVar);
228 BESDEBUG(
"ugrid",
"processUgrArgs() - END" << endl);
233 static string arrayState(libdap::Array *dapArray,
string indent)
236 libdap::Array::Dim_iter thisDim;
237 stringstream arrayState;
239 arrayState << indent <<
"dapArray: '" << dapArray->name() <<
"' ";
240 arrayState <<
"type: '" << dapArray->type_name() <<
"' ";
242 arrayState <<
"read(): " << dapArray->read() <<
" ";
243 arrayState <<
"read_p(): " << dapArray->read_p() <<
" ";
246 for (thisDim = dapArray->dim_begin(); thisDim != dapArray->dim_end(); thisDim++) {
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) <<
" ";
254 return arrayState.str();
257 static void copyUsingSubsetIndex(libdap::Array *sourceArray, vector<unsigned int> *subsetIndex,
void *result)
259 BESDEBUG(
"ugrid",
"ugrid::copyUsingIndex() - BEGIN" << endl);
261 switch (sourceArray->var()->type()) {
263 sourceArray->value(subsetIndex, (dods_byte *) result);
266 sourceArray->value(subsetIndex, (dods_uint16 *) result);
269 sourceArray->value(subsetIndex, (dods_int16 *) result);
272 sourceArray->value(subsetIndex, (dods_uint32 *) result);
275 sourceArray->value(subsetIndex, (dods_int32 *) result);
279 sourceArray->value(subsetIndex, (dods_float32 *) result);
282 sourceArray->value(subsetIndex, (dods_float64 *) result);
286 throw InternalErr(__FILE__, __LINE__,
"ugrid::hgr5::copyUsingSubsetIndex() - Unknown DAP type encountered.");
289 BESDEBUG(
"ugrid",
"ugrid::copyUsingIndex() - END" << endl);
293 static string vectorToString(vector<unsigned int> *index)
295 BESDEBUG(
"ugrid",
"indexToString() - BEGIN"<< endl);
296 BESDEBUG(
"ugrid",
"indexToString() - index.size(): " << libdap::long_to_string(index->size()) << endl);
300 for (
unsigned int i = 0; i < index->size(); ++i) {
301 s << ((i > 0) ?
", " :
" ");
305 BESDEBUG(
"ugrid",
"indexToString() - END"<< endl);
313 static void rDAWorker(MeshDataVariable *mdv, libdap::Array::Dim_iter thisDim, vector<unsigned int> *slab_subset_index,
316 libdap::Array *dapArray = mdv->getDapArray();
320 "rDAWorker() - slab_subset_index" << vectorToString(slab_subset_index) <<
" size: " << slab_subset_index->size() << endl);
324 libdap::Array::Dim_iter locationCoordinateDim = mdv->getLocationCoordinateDimension();
326 BESDEBUG(
"ugrid",
"rdaWorker() - thisDim: '" << dapArray->dimension_name(thisDim) <<
"'" << endl);
328 "rdaWorker() - locationCoordinateDim: '" << dapArray->dimension_name(locationCoordinateDim) <<
"'" << endl);
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);
335 BESDEBUG(
"ugrid",
"rdaWorker() - array state: " << arrayState(dapArray,
" "));
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);
343 dapArray->add_constraint(thisDim, start, stride, stop);
346 BESDEBUG(
"ugrid",
"rdaWorker() - Found locationCoordinateDim" << endl);
348 if ((thisDim + 1) != dapArray->dim_end()) {
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);
355 BESDEBUG(
"ugrid",
"rdaWorker() - Arrived At Last Dimension" << endl);
357 dapArray->set_read_p(
false);
359 BESDEBUG(
"ugrid",
"rdaWorker() - dap array: " << arrayState(dapArray,
" "));
361 vector<unsigned int> lastDimHyperSlabLocation;
362 NDimensionalArray::retrieveLastDimHyperSlabLocationFromConstrainedArrray(dapArray, &lastDimHyperSlabLocation);
365 "rdaWorker() - lastDimHyperSlabLocation: " << NDimensionalArray::vectorToIndices(&lastDimHyperSlabLocation) << endl);
371 results->getNextLastDimensionHyperSlab(&slab);
375 copyUsingSubsetIndex(dapArray, slab_subset_index, slab);
386 static libdap::Array *restrictRangeVariableByOneDHyperSlab(MeshDataVariable *mdv,
387 vector<unsigned int> *slab_subset_index)
390 long restrictedSlabSize = slab_subset_index->size();
393 "restrictRangeVariableByOneDHyperSlab() - slab_subset_index"<< vectorToString(slab_subset_index) <<
" size: " << libdap::long_to_string(restrictedSlabSize) << endl);
395 libdap::Array *sourceDapArray = mdv->getDapArray();
398 "restrictDapArrayByOneDHyperSlab() - locationCoordinateDim: '" << sourceDapArray->dimension_name(mdv->getLocationCoordinateDimension()) <<
"'" << endl);
406 vector<unsigned int> resultArrayShape(sourceDapArray->dimensions(
true));
407 NDimensionalArray::computeConstrainedShape(sourceDapArray, &resultArrayShape);
410 for (
unsigned int i = 0; i < resultArrayShape.size(); i++) {
411 msg <<
"[" << resultArrayShape[i] <<
"]";
413 BESDEBUG(
"ugrid",
"restrictDapArrayByOneDHyperSlab() - Constrained source array shape" << msg.str() << endl);
419 resultArrayShape[sourceDapArray->dimensions(
true) - 1] = restrictedSlabSize;
423 "restrictDapArrayByOneDHyperSlab() - UGrid restricted HyperSlab has "<< restrictedSlabSize <<
" elements." << endl);
425 "restrictDapArrayByOneDHyperSlab() - Array is of type '"<< libdap::type_name(dapType) <<
"'" << endl);
427 for (
unsigned int i = 0; i < resultArrayShape.size(); i++) {
428 msg <<
"[" << resultArrayShape[i] <<
"]";
430 BESDEBUG(
"ugrid",
"restrictDapArrayByOneDHyperSlab() - resultArrayShape" << msg.str() << endl);
431 msg.str(std::string());
437 rDAWorker(mdv, sourceDapArray->dim_begin(), slab_subset_index, result);
440 libdap::Array *resultDapArray = result->getArray(sourceDapArray);
446 return resultDapArray;
468 void ugrid_restrict(
string func_name, locationType location,
int argc, BaseType *argv[], DDS &dds, BaseType **btpp)
473 if (BESISDEBUG(TIMING_LOG)) sw.
start(
"ugrid::ugrid_restrict()",
"[function_invocation]");
475 BESDEBUG(
"ugrid",
"ugrid_restrict() - BEGIN" << endl);
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>";
482 Str *response =
new Str(
"info");
483 response->set_value(info);
489 UgridRestrictArgs args = processUgrArgs(func_name, location, argc, argv);
494 map<string, vector<MeshDataVariable *> *> *meshToRangeVarsMap =
new map<string, vector<MeshDataVariable *> *>();
497 vector<libdap::Array *>::iterator it;
498 for (it = args.rangeVars.begin(); it != args.rangeVars.end(); ++it) {
499 addRangeVar(&dds, *it, meshToRangeVarsMap);
501 BESDEBUG(
"ugrid",
"ugrid_restrict() - The user requested "<< args.rangeVars.size() <<
" range data variables." << endl);
503 "ugrid_restrict() - The user's request referenced "<< meshToRangeVarsMap->size() <<
" mesh topology variables." << endl);
516 Structure *dapResult =
new Structure(
"ugr_result_unwrap");
521 AttrTable &globals = dds.get_attr_table();
522 BESDEBUG(
"ugrid",
"ugrid_restrict() - Copying Global Attributes" << endl << globals << endl);
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);
534 map<string, vector<MeshDataVariable *> *>::iterator mit;
535 for (mit = meshToRangeVarsMap->begin(); mit != meshToRangeVarsMap->end(); ++mit) {
537 string meshVariableName = mit->first;
538 vector<MeshDataVariable *> *requestedRangeVarsForMesh = mit->second;
544 "ugrid_restrict() - Adding restricted mesh_topology structure for mesh '" << meshVariableName <<
"' to DAP response." << endl);
546 TwoDMeshTopology *tdmt =
new TwoDMeshTopology();
547 tdmt->init(meshVariableName, &dds);
549 tdmt->buildBasicGfTopology();
550 tdmt->addIndexVariable(node);
551 tdmt->addIndexVariable(face);
552 tdmt->applyRestrictOperator(args.dimension, args.filterExpression);
555 vector<vector<unsigned int> *> location_subset_indices(3);
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]);
563 location_subset_indices[node] = &node_subset_index;
565 BESDEBUG(
"ugrid2",
"ugrid_restrict() - node_subset_index"<< vectorToString(&node_subset_index) << endl);
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]);
573 location_subset_indices[face] = &face_subset_index;
574 BESDEBUG(
"ugrid2",
"ugrid_restrict() - face_subset_index: "<< vectorToString(&face_subset_index) << endl);
578 vector<BaseType *> dapResults;
579 tdmt->convertResultGridFieldStructureToDapObjects(&dapResults);
582 "ugrid_restrict() - Restriction of mesh_topology '"<< tdmt->getMeshVariable()->name() <<
"' structure completed." << endl);
587 vector<MeshDataVariable *>::iterator rvit;
588 for (rvit = requestedRangeVarsForMesh->begin(); rvit != requestedRangeVarsForMesh->end(); rvit++) {
589 MeshDataVariable *mdv = *rvit;
592 "ugrid_restrict() - Processing MeshDataVariable '"<< mdv->getName() <<
"' associated with rank/location: "<< mdv->getGridLocation() << endl);
594 tdmt->setLocationCoordinateDimension(mdv);
600 libdap::Array *restrictedRangeVarArray = restrictRangeVariableByOneDHyperSlab(mdv,
601 location_subset_indices[mdv->getGridLocation()]);
604 "ugrid_restrict() - Adding resulting dapArray '"<< restrictedRangeVarArray->name() <<
"' to dapResults." << endl);
606 dapResults.push_back(restrictedRangeVarArray);
611 BESDEBUG(
"ugrid",
"ugrid_restrict() - Adding GF::GridField results to DAP structure " << dapResult->name() << endl);
613 for (vector<BaseType *>::iterator i = dapResults.begin(); i != dapResults.end(); ++i) {
615 "ugrid_restrict() - Adding variable "<< (*i)->name() <<
" to DAP structure " << dapResult->name() << endl);
616 dapResult->add_var_nocopy(*i);
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;
632 delete requestedRangeVarsForMesh;
634 delete meshToRangeVarsMap;
636 BESDEBUG(
"ugrid",
"ugrid_restrict() - END" << endl);
638 catch (GFError &gfe) {
639 throw BESError(gfe.get_message(), gfe.get_error_type(), gfe.get_file(), gfe.get_line());
667 void ugnr(
int argc, BaseType *argv[], DDS &dds, BaseType **btpp) {
668 ugrid_restrict(
"ugnr",node,argc,argv,dds,btpp);
691 void uger(
int argc, BaseType *argv[], DDS &dds, BaseType **btpp) {
692 ugrid_restrict(
"uger",edge,argc,argv,dds,btpp);
714 void ugfr(
int argc, BaseType *argv[], DDS &dds, BaseType **btpp) {
715 ugrid_restrict(
"ugfr",face,argc,argv,dds,btpp);