/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1796 by jfenwick, Wed Sep 17 01:45:46 2008 UTC revision 1859 by gross, Wed Oct 8 03:03:37 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  #include "Data.h"  #include "Data.h"
16    
# Line 109  Data::Data(const Data& inData, Line 108  Data::Data(const Data& inData,
108  Data::Data(const Data& inData,  Data::Data(const Data& inData,
109             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
110  {  {
111      if (inData.isEmpty())
112      {
113        throw DataException("Error - will not interpolate for instances of DataEmpty.");
114      }
115    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
116      m_data=inData.m_data;      m_data=inData.m_data;
117    } else if (inData.isConstant()) { // for a constant function, we just need to use the new function space    } else if (inData.isConstant()) { // for a constant function, we just need to use the new function space
# Line 413  Data::getShapeTuple() const Line 416  Data::getShapeTuple() const
416          throw DataException("Error - illegal Data rank.");          throw DataException("Error - illegal Data rank.");
417    }    }
418  }  }
419    
420    
421    // The different name is needed because boost has trouble with overloaded functions.
422    // It can't work out what type the function is based soley on its name.
423    // There are ways to fix this involving creating function pointer variables for each form
424    // but there doesn't seem to be a need given that the methods have the same name from the python point of view
425    Data*
426    Data::copySelf()
427    {
428       DataAbstract* temp=m_data->deepCopy();
429       return new Data(temp);
430    }
431    
432  void  void
433  Data::copy(const Data& other)  Data::copy(const Data& other)
434  {  {
435    //    DataAbstract* temp=other.m_data->deepCopy();
436    // Perform a deep copy    shared_ptr<DataAbstract> temp_data(temp);
437    {    m_data=temp_data;
     DataExpanded* temp=dynamic_cast<DataExpanded*>(other.m_data.get());  
     if (temp!=0) {  
       //  
       // Construct a DataExpanded copy  
       DataAbstract* newData=new DataExpanded(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
   }  
   {  
     DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());  
     if (temp!=0) {  
       //  
       // Construct a DataTagged copy  
       DataAbstract* newData=new DataTagged(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
   }  
   {  
     DataConstant* temp=dynamic_cast<DataConstant*>(other.m_data.get());  
     if (temp!=0) {  
       //  
       // Construct a DataConstant copy  
       DataAbstract* newData=new DataConstant(*temp);  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
   }  
   {  
     DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());  
     if (temp!=0) {  
       //  
       // Construct a DataEmpty copy  
       DataAbstract* newData=new DataEmpty();  
       shared_ptr<DataAbstract> temp_data(newData);  
       m_data=temp_data;  
       return;  
     }  
   }  
   throw DataException("Error - Copy not implemented for this Data type.");  
438  }  }
439    
440    
441  void  void
442  Data::setToZero()  Data::setToZero()
443  {  {
444      if (isEmpty())
445      {
446         throw DataException("Error - Operations not permitted on instances of DataEmpty.");
447      }
448    {    {
449      DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());      DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
450      if (temp!=0) {      if (temp!=0) {
# Line 493  Data::setToZero() Line 469  Data::setToZero()
469    throw DataException("Error - Data can not be set to zero.");    throw DataException("Error - Data can not be set to zero.");
470  }  }
471    
472    // void
473    // Data::copyWithMask(const Data& other,
474    //                    const Data& mask)
475    // {
476    //   if (other.isEmpty() || mask.isEmpty())
477    //   {
478    //  throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
479    //   }
480    //   Data mask1;
481    //   Data mask2;
482    //   mask1 = mask.wherePositive();
483    //
484    //   mask2.copy(mask1);
485    //   mask1 *= other;
486    //
487    //   mask2 *= *this;
488    //   mask2 = *this - mask2;
489    //   *this = mask1 + mask2;
490    // }
491    
492  void  void
493  Data::copyWithMask(const Data& other,  Data::copyWithMask(const Data& other,
494                     const Data& mask)                     const Data& mask)
495  {  {
496    Data mask1;    // 1. Interpolate if required so all Datas use the same FS as this
497    Data mask2;    // 2. Tag or Expand so that all Data's are the same type
498      // 3. Iterate over the data vectors copying values where mask is >0
499    mask1 = mask.wherePositive();    if (other.isEmpty() || mask.isEmpty())
500    mask2.copy(mask1);    {
501        throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
502      }
503      Data other2(other);
504      Data mask2(mask);
505      FunctionSpace myFS=getFunctionSpace();
506      FunctionSpace oFS=other2.getFunctionSpace();
507      FunctionSpace mFS=mask2.getFunctionSpace();
508      if (oFS!=myFS)
509      {
510         if (other2.probeInterpolation(myFS))
511         {
512        other2=other2.interpolate(myFS);
513         }
514         else
515         {
516        throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
517         }
518      }
519      if (mFS!=myFS)
520      {
521         if (mask2.probeInterpolation(myFS))
522         {
523        mask2=mask2.interpolate(myFS);
524         }
525         else
526         {
527        throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
528         }
529      }
530                // Ensure that all args have the same type
531      if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
532      {
533        this->expand();
534        other2.expand();
535        mask2.expand();
536      }
537      else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
538      {
539        this->tag();
540        other2.tag();
541        mask2.tag();
542      }
543      else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
544      {
545      }
546      else
547      {
548        throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
549      }
550      // Now we iterate over the elements
551      DataVector& self=m_data->getVector();
552      const DataVector& ovec=other2.m_data->getVector();
553      const DataVector& mvec=mask2.m_data->getVector();
554      if ((self.size()!=ovec.size()) || (self.size()!=mvec.size()))
555      {
556        throw DataException("Error - size mismatch in arguments to copyWithMask.");
557      }
558      size_t num_points=self.size();
559      long i;
560      #pragma omp parallel for private(i) schedule(static)
561      for (i=0;i<num_points;++i)
562      {
563        if (mvec[i]>0)
564        {
565           self[i]=ovec[i];
566        }
567      }
568    }
569    
   mask1 *= other;  
   mask2 *= *this;  
   mask2 = *this - mask2;  
570    
   *this = mask1 + mask2;  
 }  
571    
572  bool  bool
573  Data::isExpanded() const  Data::isExpanded() const
# Line 662  Data::probeInterpolation(const FunctionS Line 721  Data::probeInterpolation(const FunctionS
721  Data  Data
722  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
723  {  {
724      if (isEmpty())
725      {
726        throw DataException("Error - operation not permitted on instances of DataEmpty.");
727      }
728    double blocktimer_start = blocktimer_time();    double blocktimer_start = blocktimer_time();
729    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
730      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
# Line 676  Data::gradOn(const FunctionSpace& functi Line 739  Data::gradOn(const FunctionSpace& functi
739  Data  Data
740  Data::grad() const  Data::grad() const
741  {  {
742      if (isEmpty())
743      {
744        throw DataException("Error - operation not permitted on instances of DataEmpty.");
745      }
746    return gradOn(escript::function(getDomain()));    return gradOn(escript::function(getDomain()));
747  }  }
748    
# Line 691  Data::getLength() const Line 758  Data::getLength() const
758    return m_data->getLength();    return m_data->getLength();
759  }  }
760    
 // const DataTypes::ShapeType&  
 // Data::getDataPointShape() const  
 // {  
 //   return getPointDataView().getShape();  
 // }  
   
   
   
   
761  const  const
762  boost::python::numeric::array  boost::python::numeric::array
763  Data:: getValueOfDataPoint(int dataPointNo)  Data:: getValueOfDataPoint(int dataPointNo)
# Line 799  Data::setValueOfDataPointToPyObject(int Line 857  Data::setValueOfDataPointToPyObject(int
857      // this will throw if the value cannot be represented      // this will throw if the value cannot be represented
858      boost::python::numeric::array num_array(py_object);      boost::python::numeric::array num_array(py_object);
859      setValueOfDataPointToArray(dataPointNo,num_array);      setValueOfDataPointToArray(dataPointNo,num_array);
   
   
860  }  }
861    
862  void  void
# Line 822  Data::setValueOfDataPointToArray(int dat Line 878  Data::setValueOfDataPointToArray(int dat
878    }    }
879    //    //
880    // make sure data is expanded:    // make sure data is expanded:
881      //
882    if (!isExpanded()) {    if (!isExpanded()) {
883      expand();      expand();
884    }    }
# Line 1582  Data::calc_minGlobalDataPoint(int& ProcN Line 1639  Data::calc_minGlobalDataPoint(int& ProcN
1639  void  void
1640  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1641  {  {
1642      if (isEmpty())
1643      {
1644        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1645      }
1646    boost::python::dict args;    boost::python::dict args;
1647    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
1648    getDomain().saveDX(fileName,args);    getDomain().saveDX(fileName,args);
# Line 1591  Data::saveDX(std::string fileName) const Line 1652  Data::saveDX(std::string fileName) const
1652  void  void
1653  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1654  {  {
1655      if (isEmpty())
1656      {
1657        throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1658      }
1659    boost::python::dict args;    boost::python::dict args;
1660    args["data"]=boost::python::object(this);    args["data"]=boost::python::object(this);
1661    getDomain().saveVTK(fileName,args);    getDomain().saveVTK(fileName,args);
# Line 2010  Data::setTaggedValueFromCPP(int tagKey, Line 2075  Data::setTaggedValueFromCPP(int tagKey,
2075  int  int
2076  Data::getTagNumber(int dpno)  Data::getTagNumber(int dpno)
2077  {  {
2078      if (isEmpty())
2079      {
2080        throw DataException("Error - operation not permitted on instances of DataEmpty.");
2081      }
2082    return getFunctionSpace().getTagFromDataPointNo(dpno);    return getFunctionSpace().getTagFromDataPointNo(dpno);
2083  }  }
2084    

Legend:
Removed from v.1796  
changed lines
  Added in v.1859

  ViewVC Help
Powered by ViewVC 1.1.26