XdmfArray.cpp 63.5 KB
Newer Older
Kenneth Leiter's avatar
Kenneth Leiter committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
/*****************************************************************************/
/*                                    XDMF                                   */
/*                       eXtensible Data Model and Format                    */
/*                                                                           */
/*  Id : XdmfArray.cpp                                                       */
/*                                                                           */
/*  Author:                                                                  */
/*     Kenneth Leiter                                                        */
/*     kenneth.leiter@arl.army.mil                                           */
/*     US Army Research Laboratory                                           */
/*     Aberdeen Proving Ground, MD                                           */
/*                                                                           */
/*     Copyright @ 2011 US Army Research Laboratory                          */
/*     All Rights Reserved                                                   */
/*     See Copyright.txt for details                                         */
/*                                                                           */
/*     This software is distributed WITHOUT ANY WARRANTY; without            */
/*     even the implied warranty of MERCHANTABILITY or FITNESS               */
/*     FOR A PARTICULAR PURPOSE.  See the above copyright notice             */
/*     for more information.                                                 */
/*                                                                           */
/*****************************************************************************/
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
23

24
#include <boost/tokenizer.hpp>
25
#include <boost/assign.hpp>
26
#include <limits>
27
#include <sstream>
28
#include <utility>
29 30
#include <stack>
#include <math.h>
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
31
#include "XdmfArray.hpp"
32
#include "XdmfArrayType.hpp"
33
#include "XdmfHDF5Controller.hpp"
34
#include "XdmfHeavyDataController.hpp"
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
35
#include "XdmfVisitor.hpp"
36
#include "XdmfError.hpp"
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
37

38
std::string XdmfArray::mSupportedOperations = "|#()";
39 40
const std::string XdmfArray::mValidVariableChars =
  "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ1234567890_:.";
41
const std::string XdmfArray::mValidDigitChars = "-1234567890.";
42 43 44
// List the priorities for the operations, based on the order of operations
// The index of the corresponding operation in validOperationChars
// is the same as the index of its priority in this array
45
int XdmfArray::mOperationPriority [8] = {2, 1, 0, 0};
46 47 48 49
// The higher the value, the earlier the operation is
// evaluated in the order of operations
// With the exception of parenthesis which are evaluated
// as soon as the closing parenthesis is found
50

51 52 53 54 55 56 57 58
// Note, it doesn't handle overloaded functions well.
// Will generate errors unless overload methods are typecast.
std::map<std::string, shared_ptr<XdmfArray> (*)(std::vector<shared_ptr<XdmfArray> >)>
  XdmfArray::arrayFunctions = 
    boost::assign::map_list_of ("SUM", (shared_ptr<XdmfArray> (*)(std::vector<shared_ptr<XdmfArray> >))XdmfArray::sum)
      ("AVE", (shared_ptr<XdmfArray> (*)(std::vector<shared_ptr<XdmfArray> >))XdmfArray::ave);
std::map<char, shared_ptr<XdmfArray> (*)(shared_ptr<XdmfArray>, shared_ptr<XdmfArray>)>
  XdmfArray::operations = boost::assign::map_list_of ('|', XdmfArray::chunk) ('#', XdmfArray::interlace);
59 60 61



62 63 64 65 66
XDMF_CHILDREN_IMPLEMENTATION(XdmfArray,
                             XdmfHeavyDataController,
                             HeavyDataController,
                             Name)

67
class XdmfArray::Clear : public boost::static_visitor<void> {
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
68
public:
69

70 71 72 73 74 75
  Clear(XdmfArray * const array) :
    mArray(array)
  {
  }

  void
76
  operator()(const boost::blank & array) const
77
  {
78
    return;
79 80 81 82
  }

  template<typename T>
  void
83
  operator()(const shared_ptr<std::vector<T> > & array) const
84
  {
85
    array->clear();
86
  }
87 88 89 90 91 92 93 94 95 96

  template<typename T>
  void
  operator()(const boost::shared_array<const T> & array) const
  {
    mArray->internalizeArrayPointer();
    boost::apply_visitor(*this,
                         mArray->mArray);
  }

97
private: 
98
  XdmfArray * const mArray;
99 100
};

101
class XdmfArray::Erase : public boost::static_visitor<void> {
102 103
public:

104 105 106
  Erase(XdmfArray * const array,
        const unsigned int index) :
    mArray(array),
107 108 109
    mIndex(index)
  {
  }
110

111
  void
112
  operator()(const boost::blank & array) const
113 114 115 116
  {
    return;
  }

117 118
  template<typename T>
  void
119
  operator()(const shared_ptr<std::vector<T> > & array) const
120 121 122
  {
    array->erase(array->begin() + mIndex);
  }
123

124 125
  template<typename T>
  void
126
  operator()(const boost::shared_array<const T> & array) const
127 128 129 130 131 132
  {
    mArray->internalizeArrayPointer();
    boost::apply_visitor(*this,
                         mArray->mArray);
  }

133 134
private:

135
  XdmfArray * const mArray;
136
  const unsigned int mIndex;
137 138
};

139
class XdmfArray::GetArrayType :
140
  public boost::static_visitor<shared_ptr<const XdmfArrayType> > {
141 142
public:

143 144
  GetArrayType(const shared_ptr<XdmfHeavyDataController> & heavyDataController) :
    mHeavyDataController(heavyDataController)
145 146 147
  {
  }

148
  shared_ptr<const XdmfArrayType>
149 150 151 152 153
  getArrayType(const char * const) const
  {
    return XdmfArrayType::Int8();
  }

154
  shared_ptr<const XdmfArrayType>
155 156 157 158 159
  getArrayType(const short * const) const
  {
    return XdmfArrayType::Int16();
  }

160
  shared_ptr<const XdmfArrayType>
161 162 163 164 165
  getArrayType(const int * const) const
  {
    return XdmfArrayType::Int32();
  }

166
  shared_ptr<const XdmfArrayType>
167 168 169 170 171
  getArrayType(const long * const) const
  {
    return XdmfArrayType::Int64();
  }

172
  shared_ptr<const XdmfArrayType>
173 174 175 176 177
  getArrayType(const float * const) const
  {
    return XdmfArrayType::Float32();
  }

178
  shared_ptr<const XdmfArrayType>
179 180 181 182 183
  getArrayType(const double * const) const
  {
    return XdmfArrayType::Float64();
  }

184
  shared_ptr<const XdmfArrayType>
185 186 187 188 189
  getArrayType(const unsigned char * const) const
  {
    return XdmfArrayType::UInt8();
  }

190
  shared_ptr<const XdmfArrayType>
191 192 193 194 195
  getArrayType(const unsigned short * const) const
  {
    return XdmfArrayType::UInt16();
  }

196
  shared_ptr<const XdmfArrayType>
197 198 199 200 201
  getArrayType(const unsigned int * const) const
  {
    return XdmfArrayType::UInt32();
  }

Kenneth Leiter's avatar
Kenneth Leiter committed
202 203 204 205 206 207
  shared_ptr<const XdmfArrayType>
  getArrayType(const std::string * const) const
  {
    return XdmfArrayType::String();
  }

208
  shared_ptr<const XdmfArrayType>
209
  operator()(const boost::blank & array) const
210
  {
211 212
    if(mHeavyDataController) {
      return mHeavyDataController->getType();
213 214 215 216
    }
    return XdmfArrayType::Uninitialized();
  }

217
  template<typename T>
218 219
  shared_ptr<const XdmfArrayType>
  operator()(const shared_ptr<std::vector<T> > & array) const
220 221 222 223 224
  {
    return this->getArrayType(&(array.get()->operator[](0)));
  }

  template<typename T>
225
  shared_ptr<const XdmfArrayType>
226 227 228 229
  operator()(const boost::shared_array<const T> & array) const
  {
    return this->getArrayType(array.get());
  }
230 231 232 233

private:

  const shared_ptr<XdmfHeavyDataController> mHeavyDataController;
234 235 236 237 238
};

class XdmfArray::GetCapacity : public boost::static_visitor<unsigned int> {
public:

239 240 241
  GetCapacity()
  {
  }
242

243
  unsigned int
244
  operator()(const boost::blank & array) const
245 246 247 248
  {
    return 0;
  }

249 250
  template<typename T>
  unsigned int
251
  operator()(const shared_ptr<std::vector<T> > & array) const
252 253 254
  {
    return array->capacity();
  }
255 256 257

  template<typename T>
  unsigned int
258
  operator()(const boost::shared_array<const T> & array) const
259 260 261
  {
    return 0;
  }
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
262 263
};

264
class XdmfArray::GetValuesPointer :
265
  public boost::static_visitor<const void *> {
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
266 267
public:

268 269 270 271
  GetValuesPointer()
  {
  }

272
  const void *
273
  operator()(const boost::blank & array) const
274 275 276 277
  {
    return NULL;
  }

278
  template<typename T>
279
  const void *
280
  operator()(const shared_ptr<std::vector<T> > & array) const
281 282 283 284 285
  {
    return &array->operator[](0);
  }

  template<typename T>
286
  const void *
287 288 289 290
  operator()(const boost::shared_array<const T> & array) const
  {
    return array.get();
  }
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
291 292
};

293
class XdmfArray::GetValuesString : public boost::static_visitor<std::string> {
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
294 295
public:

296 297 298 299 300
  GetValuesString(const int arrayPointerNumValues) :
    mArrayPointerNumValues(arrayPointerNumValues)
  {
  }

301
  template<typename T, typename U>
302
  std::string
303
  getValuesString(const T * const array,
304 305
                  const int numValues) const
  {
306
    const int lastIndex = numValues - 1;
307 308 309 310 311

    if(lastIndex < 0) {
      return "";
    }

312
    std::stringstream toReturn;
313
    toReturn.precision(std::numeric_limits<U>::digits10 + 2);
314
    for(int i=0; i<lastIndex; ++i) {
315
      toReturn << (U)array[i] << " ";
316
    }
317
    toReturn << (U)array[lastIndex];
318 319 320
    return toReturn.str();
  }

321 322 323 324 325 326 327
  std::string
  getValuesString(const char * const array,
                  const int numValues) const
  {
    return getValuesString<char, int>(array, numValues);
  }

328 329 330 331
  std::string
  getValuesString(const unsigned char * const array,
                  const int numValues) const
  {
332
    return getValuesString<unsigned char, int>(array, numValues);
333 334 335 336 337 338 339
  }

  template<typename T>
  std::string
  getValuesString(const T * const array,
                  const int numValues) const
  {
340
    return getValuesString<T, T>(array, numValues);
341 342
  }

343
  std::string
344
  operator()(const boost::blank & array) const
345 346 347 348
  {
    return "";
  }

349 350
  template<typename T>
  std::string
351
  operator()(const shared_ptr<std::vector<T> > & array) const
352 353 354 355 356 357 358 359 360 361
  {
    return getValuesString(&(array->operator[](0)), array->size());
  }

  template<typename T>
  std::string
  operator()(const boost::shared_array<const T> & array) const
  {
    return getValuesString(array.get(), mArrayPointerNumValues);
  }
362 363 364

private:

365
  const unsigned int mArrayPointerNumValues;
366 367
};

368 369 370
class XdmfArray::InsertArray : public boost::static_visitor<void> {
public:

371 372
  InsertArray(XdmfArray * const array,
              const unsigned int startIndex,
373 374 375
              const unsigned int valuesStartIndex,
              const unsigned int numValues,
              const unsigned int arrayStride,
376
              const unsigned int valuesStride,
377 378 379
              std::vector<unsigned int> & dimensions,
              const shared_ptr<const XdmfArray> & arrayToCopy) :
    mArray(array),
380 381 382 383
    mStartIndex(startIndex),
    mValuesStartIndex(valuesStartIndex),
    mNumValues(numValues),
    mArrayStride(arrayStride),
384
    mValuesStride(valuesStride),
385 386
    mDimensions(dimensions),
    mArrayToCopy(arrayToCopy)
387 388 389 390
  {
  }

  void
391
  operator()(const boost::blank & array) const
392 393 394 395 396 397 398 399 400
  {
    mArray->initialize(mArrayToCopy->getArrayType());
    boost::apply_visitor(*this,
                         mArray->mArray);
  }

  template<typename T>
  void
  operator()(const shared_ptr<std::vector<T> > & array) const
401
  {
402
    unsigned int size = mStartIndex + (mNumValues - 1) * mArrayStride + 1;
403 404
    if(array->size() < size) {
      array->resize(size);
405
      mDimensions.clear();
406
    }
407 408 409 410 411 412 413 414 415 416 417 418 419 420
    mArrayToCopy->getValues(mValuesStartIndex,
                            &(array->operator[](mStartIndex)),
                            mNumValues,
                            mValuesStride,
                            mArrayStride);
  }

  template<typename T>
  void
  operator()(const boost::shared_array<const T> & array) const
  {
    mArray->internalizeArrayPointer();
    boost::apply_visitor(*this,
                         mArray->mArray);
421
  }
422 423 424

private:

425
  XdmfArray * const mArray;
426 427 428 429 430
  const unsigned int mStartIndex;
  const unsigned int mValuesStartIndex;
  const unsigned int mNumValues;
  const unsigned int mArrayStride;
  const unsigned int mValuesStride;
431
  std::vector<unsigned int> & mDimensions;
432
  const shared_ptr<const XdmfArray> mArrayToCopy;
433 434
};

435
class XdmfArray::InternalizeArrayPointer : public boost::static_visitor<void> {
436 437
public:

438 439 440 441 442
  InternalizeArrayPointer(XdmfArray * const array) :
    mArray(array)
  {
  }

443
  void
444
  operator()(const boost::blank & array) const
445 446 447 448 449 450
  {
    return;
  }

  template<typename T>
  void
451
  operator()(const shared_ptr<std::vector<T> > & array) const
452 453 454 455
  {
    return;
  }

456 457 458 459
  template<typename T>
  void
  operator()(const boost::shared_array<const T> & array) const
  {
460 461 462 463
    const T * const pointer = array.get();
    shared_ptr<std::vector<T> > newArray(new std::vector<T>(pointer,
                                                            pointer + mArray->mArrayPointerNumValues));
    mArray->mArray = newArray;
464 465
    mArray->mArrayPointerNumValues = 0;
  }
466 467 468

private:

469
  XdmfArray * const mArray;
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
470 471
};

472 473 474 475 476 477 478 479
class XdmfArray::IsInitialized : public boost::static_visitor<bool> {
public:

  IsInitialized()
  {
  }

  bool
480
  operator()(const boost::blank &) const
481 482 483 484 485 486
  {
    return false;
  }

  template<typename T>
  bool
487
  operator()(const shared_ptr<std::vector<T> > &) const
488 489 490 491 492 493
  {
    return true;
  }

  template<typename T>
  bool
494
  operator()(const T &) const
495 496 497 498 499
  {
    return true;
  }
};

500
class XdmfArray::Reserve : public boost::static_visitor<void> {
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
501 502
public:

503 504 505 506
  Reserve(XdmfArray * const array,
          const unsigned int size):
    mArray(array),
    mSize(size)
507 508 509 510
  {
  }

  void
511
  operator()(const boost::blank & array) const
512
  {
513
    mArray->mTmpReserveSize = mSize;
514
  }
515

516 517 518
  template<typename T>
  void
  operator()(shared_ptr<std::vector<T> > & array) const
519
  {
520
    array->reserve(mSize);
521
  }
522

523 524
  template<typename T>
  void
525
  operator()(const boost::shared_array<const T> & array) const
526
  {
527 528 529
    mArray->internalizeArrayPointer();
    boost::apply_visitor(*this,
                         mArray->mArray);
530
  }
531 532 533

private:

534
  XdmfArray * const mArray;
535
  const unsigned int mSize;
536 537
};

538 539 540
class XdmfArray::Size : public boost::static_visitor<unsigned int> {
public:

541 542 543 544 545 546
  Size(const XdmfArray * const array) :
    mArray(array)
  {
  }

  unsigned int
547
  operator()(const boost::blank & array) const
548
  {
549
    if(mArray->mHeavyDataControllers.size()>0) {
550
      int total = 0;
551
      for (unsigned int i = 0; i < mArray->mHeavyDataControllers.size(); ++i) {
552 553
        total += mArray->mHeavyDataControllers[i]->getSize();
      }
554
      return total;
555 556
    }
    return 0;
557
  }
558

559 560
  template<typename T>
  unsigned int
561
  operator()(const shared_ptr<std::vector<T> > & array) const
562 563 564
  {
    return array->size();
  }
565 566 567

  template<typename T>
  unsigned int
568
  operator()(const boost::shared_array<const T> & array) const
569 570 571 572 573 574
  {
    return mArray->mArrayPointerNumValues;
  }

private:

575
  const XdmfArray * const mArray; 
576 577
};

578
shared_ptr<XdmfArray>
579
XdmfArray::New()
580
{
581
  shared_ptr<XdmfArray> p(new XdmfArray());
582
  return p;
583 584
}

Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
585
XdmfArray::XdmfArray() :
586 587 588
  mArrayPointerNumValues(0),
  mName(""),
  mTmpReserveSize(0)
Ken Leiter (Civ ARL/CISD) kleiter's avatar
ENH:  
Ken Leiter (Civ ARL/CISD) kleiter committed
589 590 591 592 593 594 595
{
}

XdmfArray::~XdmfArray()
{
}

596
const std::string XdmfArray::ItemTag = "DataItem";
597

598
int
599 600
XdmfArray::addFunction(std::string name,
                       shared_ptr<XdmfArray>(*functionref)(std::vector<shared_ptr<XdmfArray> >))
601
{
602
  // Check to ensure that the name has valid characters
603
  for (unsigned int i = 0; i < name.size(); ++i) {
604
    // If the character is not found in the list of valid characters
605
    if (mValidVariableChars.find(name[i]) == std::string::npos) {
606
      // Then throw an error
607 608 609 610 611 612 613 614 615 616 617
      try {
        XdmfError::message(XdmfError::FATAL,
                           "Error: Function Name Contains Invalid Character(s)");
      }
      catch (XdmfError e) {
        throw e;
      }
    }
  }
  size_t origsize = arrayFunctions.size();
  arrayFunctions[name] = functionref;
618
  // If no new functions were added
619
  if (origsize == arrayFunctions.size()) {
620
    // Toss a warning, it's nice to let people know that they're doing this
621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645
    try {
      XdmfError::message(XdmfError::WARNING,
                         "Warning: Function Overwritten");
    }
    catch (XdmfError e) {
      throw e;
    }
  }
  return arrayFunctions.size();
}

shared_ptr<XdmfArray>
XdmfArray::ave(std::vector<shared_ptr<XdmfArray> > values)
{
        double total = sum(values)->getValue<double>(0);;
        int totalSize = 0;
        for (unsigned int i = 0; i < values.size(); ++i)
        {
                totalSize += values[i]->getSize();
        }
        shared_ptr<XdmfArray> returnArray = XdmfArray::New();
        returnArray->insert(0, total/totalSize);
        return returnArray;
}

646 647
void
XdmfArray::clear()
648
{
649
  boost::apply_visitor(Clear(this), 
650 651
                       mArray);
  mDimensions.clear();
652 653
}

654 655 656
shared_ptr<XdmfArray>
XdmfArray::chunk(shared_ptr<XdmfArray> val1, shared_ptr<XdmfArray> val2)
{
657 658
  // Join chunk (add the new array to the end of the first one)
  // Joins into new array and returns it
659 660 661 662 663 664
  shared_ptr<XdmfArray> returnArray = XdmfArray::New();
  returnArray->insert(0, val1, 0, val1->getSize(),  1, 1);
  returnArray->insert(val1->getSize(), val2, 0, val2->getSize(), 1, 1);
  return returnArray;
}

665 666
void
XdmfArray::erase(const unsigned int index)
667
{
668 669 670 671
  boost::apply_visitor(Erase(this,
                             index),
                       mArray);
  mDimensions.clear();
672 673
}

674 675 676 677
// This function and the functions it depends on are reimplemented
// in XdmfCore.i in order to properly interface with python.
// There are differences between the versions,
// but the overall algorithm remains mostly the same.
678
shared_ptr<XdmfArray>
679 680 681
XdmfArray::evaluateExpression(std::string expression,
                              std::map<std::string,
                              shared_ptr<XdmfArray> > variables)
682
{
683 684 685
  std::stack<shared_ptr<XdmfArray> > valueStack;
  std::stack<char> operationStack;

686 687
  // String is parsed left to right
  // Elements of the same priority are evaluated right to left
688
  for (unsigned int i = 0; i < expression.size(); ++i) {
689 690 691
    // Found to be a digit
    if (mValidDigitChars.find(expression[i]) != std::string::npos) {
      // Progress until a non-digit is found
692
      int valueStart = i;
693 694 695 696
      if (i + 1 < expression.size()) {
        while (mValidDigitChars.find(expression[i+1]) != std::string::npos) {
          i++;
        }
697
      }
698
      // Push back to the value stack
699
      shared_ptr<XdmfArray> valueArray = XdmfArray::New();
700
      // Use this to convert to double
701
      valueArray->insert(0, atof(expression.substr(valueStart, i + 1 - valueStart).c_str()));
702 703
      valueStack.push(valueArray);
    }
704
    else if (mValidVariableChars.find(expression[i]) != std::string::npos) {
705
      // Found to be a variable
706
      int valueStart = i;
707
      // Progress until a nonvariable value is found
708 709 710 711
      if (i+1 < expression.size()){ 
        while (mValidVariableChars.find(expression[i+1]) != std::string::npos) {
          i++;
        }
712
      }
713 714 715 716 717
      // Convert to equivalent
      if (variables.find(expression.substr(valueStart, i + 1 - valueStart))
          == variables.end()) {
        if (arrayFunctions.find(expression.substr(valueStart, i + 1 - valueStart))
            == arrayFunctions.end()) {
718 719
          try {
            XdmfError::message(XdmfError::FATAL,
720 721
                               "Error: Invalid Variable in evaluateExpression "
                               + expression.substr(valueStart, i + 1 - valueStart));
722 723 724 725
          }
          catch (XdmfError e) {
            throw e;
          }
726 727 728
        }
        else {
          std::string currentFunction = expression.substr(valueStart, i + 1 - valueStart);
729
          // Check if next character is an open parenthesis
730 731 732 733
          if (i+1 >= expression.size()) {
            if (expression[i+1] != '(') {
              try {
                XdmfError::message(XdmfError::FATAL,
734 735
                                   "Error: No values supplied to function "
                                   + expression.substr(valueStart, i + 1 - valueStart));
736 737 738 739
              }
              catch (XdmfError e) {
                throw e;
              }
740
            }
741
          }
742
          // If it is grab the string between paranthesis
743 744 745

          if (i + 2 >= expression.size()) {
            XdmfError::message(XdmfError::FATAL,
746 747
                               "Error: Missing closing parethesis to function "
                               + expression.substr(valueStart, i + 1 - valueStart));
748
          }
749 750 751 752 753 754 755 756 757 758 759 760 761 762
          i = i + 2;
          valueStart = i;
          int numOpenParenthesis = 0;
          while ((expression[i] != ')' || numOpenParenthesis) && i < expression.size()) {
            if (expression[i] == '(') {
              numOpenParenthesis++;
            }
            else if (expression[i] == ')') {
              numOpenParenthesis--;
            }
            i++;
          }
          std::string functionParameters = expression.substr(valueStart, i - valueStart);
          std::vector<shared_ptr<XdmfArray> > parameterVector;
763
          // Split that string at commas
764 765 766 767
          size_t parameterSplit = 0;
          while (parameterSplit != std::string::npos) {
            parameterSplit = 0;
            parameterSplit = functionParameters.find_first_of(",", parameterSplit);
768
            // Feed the substrings to the parse function
769 770 771 772 773 774 775 776 777 778 779 780
            if (parameterSplit == std::string::npos) {
              parameterVector.push_back(evaluateExpression(functionParameters, variables));
            }
            else {
              parameterVector.push_back(evaluateExpression(functionParameters.substr(0, parameterSplit), variables));
              functionParameters = functionParameters.substr(parameterSplit+1);
            }
          }
          valueStack.push(evaluateFunction(parameterVector, currentFunction));
        }
      }
      else {
781
        // Push equivalent to value stack
782 783 784
        valueStack.push(variables.find(expression.substr(valueStart, i + 1 - valueStart))->second);
      }
    }
785
    else if (mSupportedOperations.find(expression[i]) != std::string::npos) {
786 787
      // Found to be an operation
      // Ppop operations off the stack until one of a lower or equal importance is found
788 789
      if (operationStack.size() > 0) {
        if (expression[i] == ')') {
790
          // To close a parenthesis pop off all operations until another parentheis is found
791
          while (operationStack.size() > 0 && operationStack.top() != '(') {
792
            // Must be at least two values for this loop to work properly
793
            if (valueStack.size() < 2) {
794 795 796 797 798 799 800
              try {
                XdmfError::message(XdmfError::FATAL,
                                   "Error: Not Enough Values in evaluateExpression");
              }
              catch (XdmfError e) {
                throw e;
              }
801 802 803 804 805 806 807 808 809 810 811 812 813
            }
            else {
              shared_ptr<XdmfArray> val2 = valueStack.top();
              valueStack.pop();
              shared_ptr<XdmfArray> val1 = valueStack.top();
              valueStack.pop();
              valueStack.push(evaluateOperation(val1, val2, operationStack.top()));
              operationStack.pop();
            }
          }
          operationStack.pop();
        }
        else if (expression[i] == '(') {
814 815 816
          // Just add it if it's a start parenthesis
          // Nothing happens here in that case
          // Addition happens after the if statement
817 818 819 820
        }
        else {
          int operationLocation = getOperationPriority(expression[i]);
          int topOperationLocation = getOperationPriority(operationStack.top());
821
          // See order of operations to determine importance
822
          while (operationStack.size() > 0 && operationLocation < topOperationLocation) {
823
            // Must be at least two values for this loop to work properly
824
            if (valueStack.size() < 2) {
825 826 827 828 829 830 831
              try {
                XdmfError::message(XdmfError::FATAL,
                                   "Error: Not Enough Values in evaluateExpression");
              }
              catch (XdmfError e) {
                throw e;
              }
832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848
            }
            else {
              shared_ptr<XdmfArray> val2 = valueStack.top();
              valueStack.pop();
              shared_ptr<XdmfArray> val1 = valueStack.top();
              valueStack.pop();
              valueStack.push(evaluateOperation(val1, val2, operationStack.top()));
              operationStack.pop();
              if (operationStack.size() == 0) {
                break;
              }
              topOperationLocation = getOperationPriority(operationStack.top());
            }
          }
        }
      }
      if (expression[i] != ')') {
849
        // Add the operation to the operation stack
850 851 852
        operationStack.push(expression[i]);
      }
    }
853
    // If not a value or operation the character is ignored
854
  }
855

856
  // Empty what's left in the stacks before finishing
857
  while (valueStack.size() > 1 && operationStack.size() > 0) {
858 859
    if (valueStack.size() < 2) {
      // Must be at least two values for this loop to work properly
860 861 862 863 864 865 866
      try {
        XdmfError::message(XdmfError::FATAL,
                           "Error: Not Enough Values in evaluateExpression");
      }
      catch (XdmfError e) {
        throw e;
      }
867 868 869
    }
    else {
      if(operationStack.top() == '(') {
870 871 872 873 874 875 876
        try {
          XdmfError::message(XdmfError::WARNING,
                             "Warning: Unpaired Parenthesis");
        }
        catch (XdmfError e) {
          throw e;
        }
877 878 879 880 881 882 883
      }
      else {
        shared_ptr<XdmfArray> val2 = valueStack.top();
        valueStack.pop();
        shared_ptr<XdmfArray> val1 = valueStack.top();
        valueStack.pop();
        if (operationStack.size() == 0) {
884 885 886 887 888 889 890
          try {
            XdmfError::message(XdmfError::FATAL,
                               "Error: Not Enough Operators in evaluateExpression");
          }
          catch (XdmfError e) {
            throw e;
          }
891 892 893 894 895 896 897 898
        }
        else {
          valueStack.push(evaluateOperation(val1, val2, operationStack.top()));
          operationStack.pop();
        }
      }
    }
  }
899

900
  // Throw error if there's extra operations
901
  if (operationStack.size() > 0) {
902 903 904 905 906 907 908
    try {
      XdmfError::message(XdmfError::WARNING,
                         "Warning: Left Over Operators in evaluateExpression");
    }
    catch (XdmfError e) {
      throw e;
    }
909
  }
910

911
  if (valueStack.size() > 1) {
912 913 914 915 916 917 918
    try {
      XdmfError::message(XdmfError::WARNING,
                         "Warning: Left Over Values in evaluateExpression");
    }
    catch (XdmfError e) {
      throw e;
    }
919
  }
920

921
  return valueStack.top();
922 923 924
}

shared_ptr<XdmfArray>
925 926 927
XdmfArray::evaluateOperation(shared_ptr<XdmfArray> val1,
                             shared_ptr<XdmfArray> val2,
                             char operation)
928
{
929
  if (operations.find(operation) != operations.end()) {
930 931 932
    return (*(shared_ptr<XdmfArray>(*)(shared_ptr<XdmfArray>,
                                       shared_ptr<XdmfArray>))operations[operation])(val1,
                                                                                     val2);
933 934 935 936
  }
  else {
    return shared_ptr<XdmfArray>();
  }
937 938 939
}

int
940 941 942 943
XdmfArray::addOperation(char newoperator,
                        shared_ptr<XdmfArray>(*operationref)(shared_ptr<XdmfArray>,
                                                             shared_ptr<XdmfArray>),
                        int priority)
944
{
945
  if (newoperator == '(' || newoperator == ')') {
946 947 948 949 950 951 952
    try {
      XdmfError::message(XdmfError::FATAL,
                         "Error: Parenthesis can not be redefined");
    }
    catch (XdmfError e) {
      throw e;
    }
953
  }
954
  // Give warning if the operation already exists
955
  size_t origsize = operations.size();
956
  // Place reference in the associated location
957
  operations[newoperator] = operationref;
958 959
  // It's nice to let people know they're doing this
  // So they don't get surprised about changes in behavior
960
  if (origsize == operations.size()) {
961 962 963 964 965 966 967
    try {
      XdmfError::message(XdmfError::WARNING,
                         "Warning: Function Overwritten");
    }
    catch (XdmfError e) {
      throw e;
    }
968 969
    // Overwrite the existing info for that operation
    // Add the priority to the specified location in the priority array
970 971 972 973
    size_t priorityLocation = mSupportedOperations.find(newoperator);
    mOperationPriority[priorityLocation] = priority;
  }
  else {
974 975 976 977
    // Create new operation
    // As long as the operation isn't a valid function character
    if (mValidVariableChars.find(newoperator) != std::string::npos
        || mValidDigitChars.find(newoperator) != std::string::npos) {
978 979 980 981 982 983 984
      try {
        XdmfError::message(XdmfError::FATAL,
                           "Error: Operation Overlaps with Variables");
      }
      catch (XdmfError e) {
        throw e;
      }
985
    }
986
    else {
987 988
      // Build the operation
      // Add operation to the supported character string
989
      mSupportedOperations.push_back(newoperator);
990
      unsigned int priorityArraySize = sizeof(mOperationPriority)/sizeof(int);
991
      // First check to see if the priority array is large enough
992
      if (mSupportedOperations.size()-1 > priorityArraySize) {
993
        // If it isn't make it bigger, double size should be fine
994 995
        int newArray [priorityArraySize*2];
        std::copy(mOperationPriority, mOperationPriority+(priorityArraySize-1), newArray);
996
	*mOperationPriority = *newArray;
997 998 999 1000 1001 1002
      }
      size_t priorityLocation = mSupportedOperations.find(newoperator);
      mOperationPriority[priorityLocation] = priority;
    }
  }
  return operations.size();
1003 1004
}

1005
// This is how you use references to functions
1006
shared_ptr<XdmfArray>
1007 1008
XdmfArray::evaluateFunction(std::vector<shared_ptr<XdmfArray> > valueVector,
                            std::string functionName)
1009
{
1010 1011 1012 1013 1014 1015
  if (arrayFunctions.find(functionName) != arrayFunctions.end()) {
    return (*(shared_ptr<XdmfArray>(*)(std::vector<shared_ptr<XdmfArray> >))arrayFunctions[functionName])(valueVector);
  }
  else {
    return shared_ptr<XdmfArray>();
  }
1016 1017
}

1018 1019
shared_ptr<const XdmfArrayType>
XdmfArray::getArrayType() const
1020
{
1021 1022 1023
  if (mHeavyDataControllers.size()>0) {
    return boost::apply_visitor(GetArrayType(mHeavyDataControllers[0]), 
                                mArray);
1024
  }
1025 1026 1027
  else {
    return boost::apply_visitor(GetArrayType(shared_ptr<XdmfHDF5Controller>()),
                                mArray);
1028
  }
1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044
}

unsigned int
XdmfArray::getCapacity() const
{
  return boost::apply_visitor(GetCapacity(), 
                              mArray);
}

std::vector<unsigned int>
XdmfArray::getDimensions() const
{
  if(mDimensions.size() == 0) {
    if(!this->isInitialized() && mHeavyDataControllers.size() > 0) {
      std::vector<unsigned int> returnDimensions;
      std::vector<unsigned int> tempDimensions;
1045
      // Find the controller with the most dimensions
1046 1047 1048 1049 1050 1051 1052 1053 1054 1055
      int dimControllerIndex = 0;
      unsigned int dimSizeMax = 0;
      unsigned int dimTotal = 0;
      for (unsigned int i = 0; i < mHeavyDataControllers.size(); ++i) {
        dimTotal += mHeavyDataControllers[i]->getSize();
        if (mHeavyDataControllers[i]->getSize() > dimSizeMax) {
          dimSizeMax = mHeavyDataControllers[i]->getSize();
          dimControllerIndex = i;
        }
      }
1056
      // Total up the size of the lower dimensions
1057
      int controllerDimensionSubtotal = 1;
1058 1059 1060
      for (unsigned int i = 0;
           i < mHeavyDataControllers[dimControllerIndex]->getDimensions().size() - 1;
           ++i) {
1061 1062 1063
        returnDimensions.push_back(mHeavyDataControllers[dimControllerIndex]->getDimensions()[i]);
        controllerDimensionSubtotal *= mHeavyDataControllers[dimControllerIndex]->getDimensions()[i];
      }
1064
      // Divide the total contained by the dimensions by the size of the lower dimensions
1065 1066
      returnDimensions.push_back(dimTotal/controllerDimensionSubtotal);
      return returnDimensions;
1067
    }
1068 1069
    const unsigned int size = this->getSize();
    return std::vector<unsigned int>(1, size);
1070
  }
1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088
  return mDimensions;
}

std::string
XdmfArray::getDimensionsString() const
{
  const std::vector<unsigned int> dimensions = this->getDimensions();
  return GetValuesString(dimensions.size()).getValuesString(&dimensions[0],
                                                            dimensions.size());
}

std::map<std::string, std::string>
XdmfArray::getItemProperties() const
{
  std::map<std::string, std::string> arrayProperties;
  if(mHeavyDataControllers.size()>0) {
    arrayProperties.insert(std::make_pair("Format",
                                          mHeavyDataControllers[0]->getName()));
1089
  }
1090 1091
  else {
    arrayProperties.insert(std::make_pair("Format", "XML"));
1092
  }
1093 1094 1095 1096
  arrayProperties.insert(std::make_pair("Dimensions", 
                                        this->getDimensionsString()));
  if(mName.compare("") != 0) {
    arrayProperties.insert(std::make_pair("Name", mName));
1097
  }
1098
  shared_ptr<const XdmfArrayType> type = this->getArrayType();
1099 1100
  type->getProperties(arrayProperties);
  return arrayProperties;
1101 1102
}

1103 1104
std::string
XdmfArray::getItemTag() const
1105
{
1106
  return ItemTag;
1107 1108
}

1109 1110
std::string
XdmfArray::getName() const