/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Diff of /trunk/ripley/src/RipleyDomain.cpp

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

revision 3670 by caltinay, Wed Nov 16 02:21:18 2011 UTC revision 3785 by caltinay, Wed Jan 25 04:00:33 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 12  Line 12 
12  *******************************************************/  *******************************************************/
13    
14  #include <ripley/RipleyDomain.h>  #include <ripley/RipleyDomain.h>
 #include <ripley/IndexList.h>  
 #include <ripley/Util.h>  
 extern "C" {  
 #include "esysUtils/blocktimer.h"  
 }  
 #include <escript/Data.h>  
15  #include <escript/DataFactory.h>  #include <escript/DataFactory.h>
16    #include <escript/FunctionSpaceFactory.h>
17    #include <pasowrap/SystemMatrixAdapter.h>
18    #include <pasowrap/TransportProblemAdapter.h>
19    
 #ifdef USE_NETCDF  
 #include <netcdfcpp.h>  
 #endif  
   
 #ifdef USE_PARMETIS  
 #include <parmetis.h>  
 #endif  
   
 #include <boost/python/import.hpp>  
 #include <boost/python/tuple.hpp>  
20  #include <iomanip>  #include <iomanip>
21    
22  using namespace std;  using namespace std;
23    using paso::SystemMatrixAdapter;
24    using paso::TransportProblemAdapter;
25    
26  namespace ripley {  namespace ripley {
27    
28  //  escript::Domain_ptr RipleyDomain::loadMesh(const string& filename)
 // define the static constants  
 RipleyDomain::FunctionSpaceNamesMapType RipleyDomain::m_functionSpaceTypeNames;  
   
 RipleyDomain::RipleyDomain(const string& name, dim_t numDim, Esys_MPIInfo* mpiInfo) :  
     m_name(name)  
29  {  {
30      setFunctionSpaceTypeNames();      throw RipleyException("loadMesh() not implemented");
     m_mpiInfo = Esys_MPIInfo_getReference(mpiInfo);  
     m_nodes.reset(new NodeFile(numDim, m_mpiInfo));  
     if (numDim==3) {  
         m_elements.reset(new ElementFile(Hex8, m_mpiInfo));  
         m_faceElements.reset(new ElementFile(Rec4, m_mpiInfo));  
     } else if (numDim==2) {  
         m_elements.reset(new ElementFile(Rec4, m_mpiInfo));  
         m_faceElements.reset(new ElementFile(Line2, m_mpiInfo));  
     } else  
         throw RipleyException("RipleyDomain: Invalid number of dimensions");  
     m_points.reset(new ElementFile(Point1, m_mpiInfo));  
     m_fullFullPattern = NULL;  
     m_fullReducedPattern = NULL;  
     m_reducedFullPattern = NULL;  
     m_reducedReducedPattern = NULL;  
31  }  }
32    
33  RipleyDomain::~RipleyDomain()  escript::Domain_ptr RipleyDomain::readMesh(const string& filename)
34  {  {
35      Esys_MPIInfo_free(m_mpiInfo);      throw RipleyException("readMesh() not implemented");
36  }  }
37    
38  int RipleyDomain::getMPISize() const  RipleyDomain::RipleyDomain(dim_t dim) :
39        m_numDim(dim),
40        m_status(0)
41  {  {
42      return m_mpiInfo->size;      m_mpiInfo = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
43  }  }
44    
45  int RipleyDomain::getMPIRank() const  RipleyDomain::~RipleyDomain()
46  {  {
47      return m_mpiInfo->rank;      Esys_MPIInfo_free(m_mpiInfo);
48  }  }
49    
50  void RipleyDomain::MPIBarrier() const  bool RipleyDomain::operator==(const AbstractDomain& other) const
51  {  {
52  #ifdef ESYS_MPI      const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53      MPI_Barrier(m_mpiInfo->comm);      if (o) {
54  #endif          return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55                    && m_elementTags==o->m_elementTags
56                    && m_faceTags==o->m_faceTags);
57        }
58        return false;
59  }  }
60    
61  bool RipleyDomain::onMasterProcessor() const  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62  {  {
63      return m_mpiInfo->rank == 0;      switch (fsType) {
64            case DegreesOfFreedom:
65            case ReducedDegreesOfFreedom:
66            case Nodes:
67            case ReducedNodes:
68            case Elements:
69            case ReducedElements:
70            case FaceElements:
71            case ReducedFaceElements:
72            case Points:
73                return true;
74            default:
75                break;
76        }
77        return false;
78  }  }
79    
80  #ifdef ESYS_MPI  string RipleyDomain::functionSpaceTypeAsString(int fsType) const
 MPI_Comm  
 #else  
 unsigned int  
 #endif  
 RipleyDomain::getMPIComm() const  
81  {  {
82  #ifdef ESYS_MPI      switch (fsType) {
83      return m_mpiInfo->comm;          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84  #else          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85      return 0;          case Nodes: return "Ripley_Nodes";
86  #endif          case ReducedNodes: return "Ripley_ReducedNodes";
87            case Elements: return "Ripley_Elements";
88            case ReducedElements: return "Ripley_ReducedElements";
89            case FaceElements: return "Ripley_FaceElements";
90            case ReducedFaceElements: return "Ripley_ReducedFaceElements";
91            case Points: return "Ripley_Points";
92            default:
93                break;
94        }
95        return "Invalid function space type code";
96  }  }
97    
98  void RipleyDomain::write(const string& filename) const  pair<int,int> RipleyDomain::getDataShape(int fsType) const
99  {  {
100      if (m_mpiInfo->size > 1)      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101          throw RipleyException("write: only single processor runs are supported.");      switch (fsType) {
102            case Nodes:
103      // open file          case ReducedNodes: //FIXME: reduced
104      ofstream f(filename.c_str());              return pair<int,int>(1, getNumNodes());
105      if (!f.is_open()) {          case DegreesOfFreedom:
106          stringstream msg;          case ReducedDegreesOfFreedom: //FIXME: reduced
107          msg << "write: Opening file " << filename << " for writing failed.";              return pair<int,int>(1, getNumDOF());
108          throw RipleyException(msg.str());          case Elements:
109                return pair<int,int>(ptsPerSample, getNumElements());
110            case FaceElements:
111                return pair<int,int>(ptsPerSample/2, getNumFaceElements());
112            case ReducedElements:
113                return pair<int,int>(1, getNumElements());
114            case ReducedFaceElements:
115                return pair<int,int>(1, getNumFaceElements());
116            case Points:
117                return pair<int,int>(1, 0); //FIXME: dirac
118            default:
119                break;
120      }      }
121    
122      // write header      stringstream msg;
123      f << m_name << endl;      msg << "getDataShape(): Unsupported function space type " << fsType
124            << " for " << getDescription();
125      f.precision(15);      throw RipleyException(msg.str());
126      f.setf(ios::scientific, ios::floatfield);  }
   
     // write nodes  
     dim_t numDim = getDim();  
     f << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;  
     for (dim_t i = 0; i < m_nodes->getNumNodes(); i++) {  
         f << m_nodes->getIdVector()[i] << " "  
             << m_nodes->getGlobalDegreesOfFreedom()[i] << " "  
             << m_nodes->getTagVector()[i];  
         for (dim_t j = 0; j < numDim; j++)  
             f << " " << setw(20) << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];  
         f << endl;  
     }  
   
     // write all element types  
     const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };  
   
     for (size_t k=0; k<3; k++) {  
         f << eFiles[k]->getName() << " " << eFiles[k]->getNumElements() << endl;  
         dim_t NN = eFiles[k]->getNumNodes();  
         for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {  
             f << eFiles[k]->getIdVector()[i] << " "  
                 << eFiles[k]->getTagVector()[i];  
             for (dim_t j = 0; j < NN; j++)  
                 f << " " << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];  
             f << endl;  
         }  
     }  
127    
128      // write tags  string RipleyDomain::showTagNames() const
129      if (m_tagMap.size() > 0) {  {
130          f << "Tags" << endl;      stringstream ret;
131          TagMap::const_iterator it;      TagMap::const_iterator it;
132          for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {      for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
133              f << it->first << " " << it->second << endl;          if (it!=m_tagMap.begin()) ret << ", ";
134          }          ret << it->first;
135      }      }
136      f.close();      return ret.str();
 #ifdef Ripley_TRACE  
     cout << "mesh " << m_name << " has been written to file " << filename << endl;  
 #endif  
137  }  }
138    
139  void RipleyDomain::Print_Mesh_Info(const bool full) const  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
140  {  {
141      cout << "Ripley_PrintMesh_Info running on CPU " <<     /*
142          m_mpiInfo->rank << " of " << m_mpiInfo->size << endl;      The idea is to use equivalence classes (i.e. types which can be
143      cout << "\tMesh name '" << m_name << "'" << endl;      interpolated back and forth):
144        class 0: DOF <-> Nodes
145      // write nodes      class 1: ReducedDOF <-> ReducedNodes
146      int numDim = getDim();      class 2: Points
147      cout << "\tNodes: " << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;      class 3: Elements
148      if (full) {      class 4: ReducedElements
149          cout << "\t     Id   Tag  gDOF   gNI grDfI  grNI:  Coordinates" << endl;      class 5: FaceElements
150          cout.precision(15);      class 6: ReducedFaceElements
151          cout.setf(ios::scientific, ios::floatfield);  
152          for (int i = 0; i < m_nodes->getNumNodes(); i++) {      There is also a set of lines. Interpolation is possible down a line but not
153              cout << "\t  " << setw(5) << m_nodes->getIdVector()[i] << " "      between lines.
154                 << setw(5) << m_nodes->getTagVector()[i] << " "      class 0 and 1 belong to all lines so aren't considered.
155                 << setw(5) << m_nodes->getGlobalDegreesOfFreedom()[i] << " "      line 0: class 2
156                 << setw(5) << m_nodes->getGlobalNodesIndex()[i] << " "      line 1: class 3,4
157                 << setw(5) << m_nodes->getGlobalReducedDOFIndex()[i] << " "      line 2: class 5,6
                << setw(5) << m_nodes->getGlobalReducedNodesIndex()[i] << ": ";  
             for (int j = 0; j < numDim; j++)  
                 cout << " " << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];  
             cout << endl;  
         }  
     }  
158    
159      // write all element types      For classes with multiple members (eg class 1) we have vars to record if
160      const char *eNames[3] = { "Elements", "Face elements", "Points" };      there is at least one instance. e.g. hasnodes is true if we have at least
161      const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };      one instance of Nodes.
162        */
163      for (size_t k=0; k<3; k++) {      if (fs.empty())
164          int mine = 0, overlap = 0;          return false;
165          for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {      vector<bool> hasclass(7, false);
166              if (eFiles[k]->getOwnerVector()[i] == m_mpiInfo->rank)      vector<int> hasline(3, 0);
167                  mine++;      bool hasnodes=false;
168              else      bool hasrednodes=false;
169                  overlap++;      for (size_t i=0; i<fs.size(); ++i) {
170          }          switch (fs[i]) {
171          cout << "\t" << eNames[k] << ": " << eFiles[k]->getName()              case Nodes: hasnodes=true; // fall through
172              << " " << eFiles[k]->getNumElements()              case DegreesOfFreedom:
173              << " (TypeId=" << eFiles[k]->getTypeId() << ") owner="                  hasclass[0]=true;
174              << mine << " overlap=" << overlap << endl;                  break;
175          int NN = eFiles[k]->getNumNodes();              case ReducedNodes: hasrednodes=true; // fall through
176          if (full) {              case ReducedDegreesOfFreedom:
177              cout << "\t     Id   Tag Owner Color:  Nodes" << endl;                  hasclass[1]=true;
178              for (int i = 0; i < eFiles[k]->getNumElements(); i++) {                  break;
179                  cout << "\t  " << setw(5) << eFiles[k]->getIdVector()[i] << " "              case Points:
180                     << setw(5) << eFiles[k]->getTagVector()[i] << " "                  hasline[0]=1;
181                     << setw(5) << eFiles[k]->getOwnerVector()[i] << " "                  hasclass[2]=true;
182                     << setw(5) << eFiles[k]->getColorVector()[i] << ": ";                  break;
183                  for (int j = 0; j < NN; j++)              case Elements:
184                      cout << " " << setw(5) << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];                  hasclass[3]=true;
185                  cout << endl;                  hasline[1]=1;
186              }                  break;
187                case ReducedElements:
188                    hasclass[4]=true;
189                    hasline[1]=1;
190                    break;
191                case FaceElements:
192                    hasclass[5]=true;
193                    hasline[2]=1;
194                    break;
195                case ReducedFaceElements:
196                    hasclass[6]=true;
197                    hasline[2]=1;
198                    break;
199                default:
200                    return false;
201          }          }
202      }      }
203        int numLines=hasline[0]+hasline[1]+hasline[2];
204    
205        // fail if we have more than one leaf group
206      // write tags      // = there are at least two branches we can't interpolate between
207      if (m_tagMap.size() > 0) {      if (numLines > 1)
208          cout << "\tTags:" << endl;          return false;
209          TagMap::const_iterator it;      else if (numLines==1) {
210          for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {          // we have points
211              cout << "\t  " << setw(5) << it->second << " "          if (hasline[0]==1)
212                  << it->first << endl;              resultcode=Points;
213            else if (hasline[1]==1) {
214                if (hasclass[4])
215                    resultcode=ReducedElements;
216                else
217                    resultcode=Elements;
218            } else { // hasline[2]==1
219                if (hasclass[6])
220                    resultcode=ReducedFaceElements;
221                else
222                    resultcode=FaceElements;
223          }          }
224        } else { // numLines==0
225            if (hasclass[1])
226                // something from class 1
227                resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228            else // something from class 0
229                resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230      }      }
231        return true;
232  }  }
233    
234  void RipleyDomain::dump(const string& fileName) const  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235                                                  int fsType_target) const
236  {  {
237  #ifdef USE_NETCDF      if (!isValidFunctionSpaceType(fsType_target)) {
238      NcDim* ncdim = NULL;          stringstream msg;
239      int mpi_size = m_mpiInfo->size;          msg << "probeInterpolationOnDomain(): Invalid functionspace type "
240      int mpi_rank = m_mpiInfo->rank;              << fsType_target << " for " << getDescription();
241      int numDim = getDim();          throw RipleyException(msg.str());
   
 /* Incoming token indicates it's my turn to write */  
 #ifdef ESYS_MPI  
     MPI_Status status;  
     if (mpi_rank>0) {  
         int dummy;  
         MPI_Recv(&dummy, 0, MPI_INT, mpi_rank-1, 81800, m_mpiInfo->comm, &status);  
242      }      }
 #endif  
   
     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),  
                                                       mpi_size, mpi_rank);  
243    
244      // NetCDF error handler      switch (fsType_source) {
245      NcError err(NcError::verbose_nonfatal);          case Nodes:
246      // Create the file.          case DegreesOfFreedom:
247      NcFile dataFile(newFileName, NcFile::Replace);              return true;
248      const string msgPrefix("dump: netCDF operation failed - ");          case ReducedNodes:
249            case ReducedDegreesOfFreedom:
250      // check if writing was successful              return (fsType_target != Nodes &&
251      if (!dataFile.is_valid())                      fsType_target != DegreesOfFreedom);
252          throw RipleyException(msgPrefix+"Open file for output");          case Elements:
253                return (fsType_target==Elements ||
254      const size_t numTags = m_tagMap.size();                      fsType_target==ReducedElements);
255            case FaceElements:
256      if (numTags > 0)              return (fsType_target==FaceElements ||
257          if (! (ncdim = dataFile.add_dim("dim_Tags", numTags)) )                      fsType_target==ReducedFaceElements);
258              throw RipleyException(msgPrefix+"add_dim(dim_Tags)");          case ReducedElements:
259            case ReducedFaceElements:
260      // Attributes: MPI size, MPI rank, Name, order, reduced_order          case Points:
261      if (!dataFile.add_att("mpi_size", mpi_size))              return (fsType_target==fsType_source);
         throw RipleyException(msgPrefix+"add_att(mpi_size)");  
     if (!dataFile.add_att("mpi_rank", mpi_rank))  
         throw RipleyException(msgPrefix+"add_att(mpi_rank)");  
     if (!dataFile.add_att("Name", m_name.c_str()))  
         throw RipleyException(msgPrefix+"add_att(Name)");  
     if (!dataFile.add_att("numDim", numDim))  
         throw RipleyException(msgPrefix+"add_att(numDim)");  
     if (!dataFile.add_att("order", 1))  
         throw RipleyException(msgPrefix+"add_att(order)");  
     if (!dataFile.add_att("reduced_order", 1))  
         throw RipleyException(msgPrefix+"add_att(reduced_order)");  
     if (!dataFile.add_att("num_Tags", static_cast<int>(numTags)))  
         throw RipleyException(msgPrefix+"add_att(num_Tags)");  
   
     m_nodes->dumpToNetCDF(dataFile);  
     m_elements->dumpToNetCDF(dataFile, "Elements");  
     m_faceElements->dumpToNetCDF(dataFile, "FaceElements");  
     m_points->dumpToNetCDF(dataFile, "Points");  
   
     // // // // // TagMap // // // // //  
     if (numTags > 0) {  
   
         // Copy tag keys into temp array  
         vector<int> Tags_keys;  
         TagMap::const_iterator it;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             Tags_keys.push_back(it->second);  
         }  
262    
263          // Tags_keys          default: {
264          NcVar *ncVar;              stringstream msg;
265          if (! (ncVar = dataFile.add_var("Tags_keys", ncInt, ncdim)) )              msg << "probeInterpolationOnDomain(): Invalid functionspace type "
266              throw RipleyException(msgPrefix+"add_var(Tags_keys)");                  << fsType_source << " for " << getDescription();
267          if (!ncVar->put(&Tags_keys[0], numTags))              throw RipleyException(msg.str());
             throw RipleyException(msgPrefix+"put(Tags_keys)");  
   
         // Tags_names_*  
         // This is an array of strings, it should be stored as an array but  
         // instead I have hacked in one attribute per string because the NetCDF  
         // manual doesn't tell how to do an array of strings  
         int i = 0;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             stringstream name;  
             name << "Tags_name_" << i;  
             if (!dataFile.add_att(name.str().c_str(), it->first.c_str()) )  
                 throw RipleyException(msgPrefix+"add_att(Tags_names_XX)");  
             i++;  
268          }          }
269      }      }
   
     // Send token to next MPI process so he can take his turn  
 #ifdef ESYS_MPI  
     if (mpi_rank<mpi_size-1) {  
         int dummy = 0;  
         MPI_Send(&dummy, 0, MPI_INT, mpi_rank+1, 81800, m_mpiInfo->comm);  
     }  
 #endif  
   
     // netCDF file is closed by destructor of NcFile object  
 #else  
     throw RipleyException("dump: not configured with netCDF. Please contact your installation manager.");  
 #endif  /* USE_NETCDF */  
270  }  }
271    
272  string RipleyDomain::getDescription() const  void RipleyDomain::interpolateOnDomain(escript::Data& target,
273  {                                         const escript::Data& in) const
     return "RipleyMesh";  
 }  
   
 string RipleyDomain::functionSpaceTypeAsString(int functionSpaceType) const  
 {  
     FunctionSpaceNamesMapType::iterator loc;  
     loc=m_functionSpaceTypeNames.find(functionSpaceType);  
     if (loc!=m_functionSpaceTypeNames.end())  
         return loc->second;  
   
     return "Invalid function space type code";  
 }  
   
 bool RipleyDomain::isValidFunctionSpaceType(int functionSpaceType) const  
 {  
     FunctionSpaceNamesMapType::iterator loc;  
     loc=m_functionSpaceTypeNames.find(functionSpaceType);  
     return (loc!=m_functionSpaceTypeNames.end());  
 }  
   
 void RipleyDomain::setFunctionSpaceTypeNames()  
 {  
     if (m_functionSpaceTypeNames.size() > 0)  
         return;  
   
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Ripley_DegreesOfFreedom"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Ripley_ReducedDegreesOfFreedom"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Nodes,"Ripley_Nodes"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedNodes,"Ripley_Reduced_Nodes"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Elements,"Ripley_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedElements,"Ripley_Reduced_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(FaceElements,"Ripley_Face_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Ripley_Reduced_Face_Elements"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Points,"Ripley_Points"));  
 }  
   
 int RipleyDomain::getContinuousFunctionCode() const  
 {  
     return Nodes;  
 }  
   
 int RipleyDomain::getReducedContinuousFunctionCode() const  
274  {  {
275      return ReducedNodes;      const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
276  }      const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
277        if (inDomain != *this)
278            throw RipleyException("Illegal domain of interpolant");
279        if (targetDomain != *this)
280            throw RipleyException("Illegal domain of interpolation target");
281    
282  int RipleyDomain::getFunctionCode() const      stringstream msg;
283  {      msg << "interpolateOnDomain() not implemented for function space "
284      return Elements;          << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
285  }          << " -> "
286            << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
287    
288        const int inFS = in.getFunctionSpace().getTypeCode();
289        const int outFS = target.getFunctionSpace().getTypeCode();
290    
291        // simplest case: 1:1 copy
292        if (inFS==outFS) {
293            copyData(target, *const_cast<escript::Data*>(&in));
294        // not allowed: reduced->non-reduced
295        } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
296                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
297            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
298        } else if ((inFS==Elements && outFS==ReducedElements)
299                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
300            if (in.actsExpanded())
301                averageData(target, *const_cast<escript::Data*>(&in));
302            else
303                copyData(target, *const_cast<escript::Data*>(&in));
304        } else {
305            switch (inFS) {
306                case Nodes:
307                case ReducedNodes: //FIXME: reduced
308                    switch (outFS) {
309                        case DegreesOfFreedom:
310                        case ReducedDegreesOfFreedom: //FIXME: reduced
311                            if (getMPISize()==1)
312                                copyData(target, *const_cast<escript::Data*>(&in));
313                            else
314                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
315                            break;
316    
317                        case Nodes:
318                        case ReducedNodes: //FIXME: reduced
319                            copyData(target, *const_cast<escript::Data*>(&in));
320                            break;
321    
322                        case Elements:
323                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
324                            break;
325    
326                        case ReducedElements:
327                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
328                            break;
329    
330                        case FaceElements:
331                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
332                            break;
333    
334                        case ReducedFaceElements:
335                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
336                            break;
337    
338  int RipleyDomain::getReducedFunctionCode() const                      default:
339  {                          throw RipleyException(msg.str());
340      return ReducedElements;                  }
341  }                  break;
342    
343  int RipleyDomain::getFunctionOnBoundaryCode() const              case DegreesOfFreedom:
344  {              case ReducedDegreesOfFreedom: //FIXME: reduced
345      return FaceElements;                  switch (outFS) {
346  }                      case Nodes:
347                        case ReducedNodes: //FIXME: reduced
348                            if (getMPISize()==1)
349                                copyData(target, *const_cast<escript::Data*>(&in));
350                            else
351                                dofToNodes(target, *const_cast<escript::Data*>(&in));
352                            break;
353    
354                        case DegreesOfFreedom:
355                        case ReducedDegreesOfFreedom: //FIXME: reduced
356                            copyData(target, *const_cast<escript::Data*>(&in));
357                            break;
358    
359                        case Elements:
360                        case ReducedElements:
361                            if (getMPISize()==1) {
362                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
363                            } else {
364                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
365                                            escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
366                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
367                            }
368                            break;
369    
370                        case FaceElements:
371                        case ReducedFaceElements:
372                            if (getMPISize()==1) {
373                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
374                            } else {
375                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
376                                         escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
377                                interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
378                            }
379                            break;
380    
381  int RipleyDomain::getReducedFunctionOnBoundaryCode() const                      default:
382  {                          throw RipleyException(msg.str());
383      return ReducedFaceElements;                  }
384                    break;
385                default:
386                    throw RipleyException(msg.str());
387            }
388        }
389  }  }
390    
391  int RipleyDomain::getFunctionOnContactZeroCode() const  escript::Data RipleyDomain::getX() const
392  {  {
393      throw RipleyException("Ripley does not support contact elements");      return escript::continuousFunction(*this).getX();
394  }  }
395    
396  int RipleyDomain::getReducedFunctionOnContactZeroCode() const  escript::Data RipleyDomain::getNormal() const
397  {  {
398      throw RipleyException("Ripley does not support contact elements");      return escript::functionOnBoundary(*this).getNormal();
399  }  }
400    
401  int RipleyDomain::getFunctionOnContactOneCode() const  escript::Data RipleyDomain::getSize() const
402  {  {
403      throw RipleyException("Ripley does not support contact elements");      return escript::function(*this).getSize();
404  }  }
405    
406  int RipleyDomain::getReducedFunctionOnContactOneCode() const  void RipleyDomain::setToX(escript::Data& arg) const
407  {  {
408      throw RipleyException("Ripley does not support contact elements");      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
409  }              *(arg.getFunctionSpace().getDomain()));
410        if (argDomain != *this)
411            throw RipleyException("setToX: Illegal domain of data point locations");
412        if (!arg.isExpanded())
413            throw RipleyException("setToX: Expanded Data object expected");
414    
415  int RipleyDomain::getSolutionCode() const      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
416  {          assembleCoordinates(arg);
417      return DegreesOfFreedom;      } else {
418            // interpolate the result
419            escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
420            assembleCoordinates(contData);
421            interpolateOnDomain(arg, contData);
422        }
423  }  }
424    
425  int RipleyDomain::getReducedSolutionCode() const  void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426  {  {
427      return ReducedDegreesOfFreedom;      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
428  }              *(arg.getFunctionSpace().getDomain()));
429        if (argDomain != *this)
430            throw RipleyException("setToGradient: Illegal domain of gradient argument");
431        const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
432                *(grad.getFunctionSpace().getDomain()));
433        if (gradDomain != *this)
434            throw RipleyException("setToGradient: Illegal domain of gradient");
435    
436  int RipleyDomain::getDiracDeltaFunctionsCode() const      switch (grad.getFunctionSpace().getTypeCode()) {
437  {          case Elements:
438      return Points;          case ReducedElements:
439  }          case FaceElements:
440            case ReducedFaceElements:
441                break;
442            default: {
443                stringstream msg;
444                msg << "setToGradient(): not supported for "
445                    << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
446                throw RipleyException(msg.str());
447            }
448        }
449    
450  //      if (getMPISize()>1) {
451  // returns the spatial dimension of the mesh          if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
452  //              escript::Data contArg(arg, escript::continuousFunction(*this));
453  int RipleyDomain::getDim() const              assembleGradient(grad, contArg);
454  {          } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
455      return m_nodes->getNumDim();              escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
456                assembleGradient(grad, contArg);
457            } else {
458                assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459            }
460        } else {
461            assembleGradient(grad, *const_cast<escript::Data*>(&arg));
462        }
463  }  }
464    
465  //  void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
 // returns the number of data points summed across all MPI processes  
 //  
 int RipleyDomain::getNumDataPointsGlobal() const  
466  {  {
467      return m_nodes->getGlobalNumNodes();      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
468  }              *(arg.getFunctionSpace().getDomain()));
469        if (argDomain != *this)
470            throw RipleyException("setToIntegrals: Illegal domain of integration kernel");
471    
472  //      switch (arg.getFunctionSpace().getTypeCode()) {
 // returns the number of data points per sample and the number of samples  
 // needed to represent data on a parts of the mesh.  
 //  
 pair<int,int> RipleyDomain::getDataShape(int functionSpaceCode) const  
 {  
     int numDataPointsPerSample=0;  
     int numSamples=0;  
     switch (functionSpaceCode) {  
473          case Nodes:          case Nodes:
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumNodes();  
             break;  
474          case ReducedNodes:          case ReducedNodes:
475              numDataPointsPerSample=1;          case DegreesOfFreedom:
476              numSamples=m_nodes->getNumReducedNodes();          case ReducedDegreesOfFreedom:
477                {
478                    escript::Data funcArg(arg, escript::function(*this));
479                    assembleIntegrate(integrals, funcArg);
480                }
481              break;              break;
482          case Elements:          case Elements:
             numSamples=m_elements->getNumElements();  
             numDataPointsPerSample=m_elements->getNumLocalDim()+1;  
             break;  
483          case ReducedElements:          case ReducedElements:
             numSamples=m_elements->getNumElements();  
             numDataPointsPerSample=(m_elements->getNumLocalDim()==0)?0:1;  
             break;  
484          case FaceElements:          case FaceElements:
             numDataPointsPerSample=m_faceElements->getNumLocalDim()+1;  
             numSamples=m_faceElements->getNumElements();  
             break;  
485          case ReducedFaceElements:          case ReducedFaceElements:
486              numDataPointsPerSample=(m_faceElements->getNumLocalDim()==0)?0:1;              assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
             numSamples=m_faceElements->getNumElements();  
             break;  
         case Points:  
             numDataPointsPerSample=1;  
             numSamples=m_points->getNumElements();  
             break;  
         case DegreesOfFreedom:  
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumDegreesOfFreedom();  
487              break;              break;
488          case ReducedDegreesOfFreedom:          default: {
489              numDataPointsPerSample=1;              stringstream msg;
490              numSamples=m_nodes->getNumReducedDegreesOfFreedom();              msg << "setToIntegrals(): not supported for "
491              break;                  << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
492          default:              throw RipleyException(msg.str());
493              stringstream temp;          }
             temp << "Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();  
             throw RipleyException(temp.str());  
494      }      }
     return pair<int,int>(numDataPointsPerSample,numSamples);  
 }  
495    
 //  
 // adds linear PDE of second order into a given stiffness matrix and right hand side  
 //  
 void RipleyDomain::addPDEToSystem(  
         escript::AbstractSystemMatrix& mat, escript::Data& rhs,  
         const escript::Data& A, const escript::Data& B, const escript::Data& C,  
         const escript::Data& D, const escript::Data& X, const escript::Data& Y,  
         const escript::Data& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac,const escript::Data& y_dirac) const  
 {  
     if (!d_contact.isEmpty() || !y_contact.isEmpty())  
         throw RipleyException("Ripley does not support contact elements");  
   
     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);  
     if (smat==NULL)  
         throw RipleyException("Ripley only accepts its own system matrices");  
 /*  
     escriptDataC _rhs=rhs.getDataC();  
     escriptDataC _A =A.getDataC();  
     escriptDataC _B=B.getDataC();  
     escriptDataC _C=C.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     Ripley_Assemble_PDE(m_nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d_dirac, 0, &_y_dirac);  
     checkPasoError();  
 */  
 }  
   
 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  
         const escript::Data& D, const escript::Data& d,  
         const escript::Data& d_dirac, const bool useHRZ) const  
 {  
 /*  
     escriptDataC _mat=mat.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     Ripley_Assemble_LumpedSystem(m_nodes, mesh->Elements, &_mat, &_D, useHRZ);  
     checkPasoError();  
   
     Ripley_Assemble_LumpedSystem(m_nodes, mesh->FaceElements, &_mat, &_d, useHRZ);  
     checkPasoError();  
   
     Ripley_Assemble_LumpedSystem(m_nodes, mesh->FaceElements, &_mat, &_d_dirac, useHRZ);  
     checkPasoError();  
 */  
 }  
   
   
 //  
 // adds linear PDE of second order into the right hand side only  
 //  
 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,  
         const escript::Data& Y, const escript::Data& y,  
         const escript::Data& y_contact, const escript::Data& y_dirac) const  
 {  
     if (!y_contact.isEmpty())  
         throw RipleyException("Ripley does not support contact elements");  
 /*  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
   
     escriptDataC _rhs=rhs.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
     Ripley_Assemble_PDE(m_nodes, mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->Points, 0, &_rhs, 0, 0, 0, 0, 0, &_y_dirac );  
     checkPasoError();  
 */  
496  }  }
497    
498  //  bool RipleyDomain::isCellOriented(int fsType) const
 // adds PDE of second order into a transport problem  
 //  
 void RipleyDomain::addPDEToTransportProblem(  
         escript::AbstractTransportProblem& tp,  
         escript::Data& source, const escript::Data& M,  
         const escript::Data& A, const escript::Data& B, const escript::Data& C,  
         const escript::Data& D, const escript::Data& X, const escript::Data& Y,  
         const escript::Data& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac, const escript::Data& y_dirac) const  
499  {  {
500      if (!d_contact.isEmpty() || !y_contact.isEmpty())      switch(fsType) {
501          throw RipleyException("Ripley does not support contact elements");          case Nodes:
502      TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);          case ReducedNodes:
503      if (tpa==NULL)          case DegreesOfFreedom:
504          throw RipleyException("Ripley only accepts its own transport problems");          case ReducedDegreesOfFreedom:
505  /*              return false;
506      DataTypes::ShapeType shape;          case Elements:
507      source.expand();          case FaceElements:
508      escriptDataC _source=source.getDataC();          case Points:
509      escriptDataC _M=M.getDataC();          case ReducedElements:
510      escriptDataC _A=A.getDataC();          case ReducedFaceElements:
511      escriptDataC _B=B.getDataC();              return true;
512      escriptDataC _C=C.getDataC();          default:
513      escriptDataC _D=D.getDataC();              break;
514      escriptDataC _X=X.getDataC();      }
515      escriptDataC _Y=Y.getDataC();      stringstream msg;
516      escriptDataC _d=d.getDataC();      msg << "isCellOriented(): Illegal function space type " << fsType
517      escriptDataC _y=y.getDataC();          << " on " << getDescription();
518      escriptDataC _d_dirac=d_dirac.getDataC();      throw RipleyException(msg.str());
     escriptDataC _y_dirac=y_dirac.getDataC();  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();  
     Ripley_Assemble_PDE(m_nodes, mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y);  
     checkPasoError();  
   
     Ripley_Assemble_PDE(m_nodes, mesh->Points, _tp->transport_matrix, &_source, 0, 0, 0, &_d_dirac, 0, &_y_dirac);  
     checkPasoError();  
 */  
519  }  }
520    
521  //  bool RipleyDomain::canTag(int fsType) const
 // interpolates data between different function spaces  
 //  
 void RipleyDomain::interpolateOnDomain(escript::Data& target,  
                                       const escript::Data& in) const  
522  {  {
523      const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));      switch(fsType) {
     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));  
     if (inDomain != *this)  
         throw RipleyException("Illegal domain of interpolant");  
     if (targetDomain != *this)  
         throw RipleyException("Illegal domain of interpolation target");  
 /*  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _target=target.getDataC();  
     escriptDataC _in=in.getDataC();  
   
     switch (in.getFunctionSpace().getTypeCode()) {  
524          case Nodes:          case Nodes:
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                 case ReducedNodes:  
                 case DegreesOfFreedom:  
                 case ReducedDegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes, &_target, &_in);  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     Ripley_Assemble_interpolate(m_nodes, mesh->Elements, &_in, &_target);  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements, &_in, &_target);  
                     break;  
                 case Points:  
                     Ripley_Assemble_interpolate(m_nodes, mesh->Points, &_in, &_target);  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation on Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
         case ReducedNodes:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                 case ReducedNodes:  
                 case DegreesOfFreedom:  
                 case ReducedDegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     Ripley_Assemble_interpolate(m_nodes,mesh->Elements,&_in,&_target);  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in,&_target);  
                     break;  
                 case Points:  
                     Ripley_Assemble_interpolate(m_nodes,mesh->Points,&_in,&_target);  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation on Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
525          case Elements:          case Elements:
             if (target.getFunctionSpace().getTypeCode()==Elements) {  
                 Ripley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
             } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
                 Ripley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on elements possible.");  
             }  
             break;  
526          case ReducedElements:          case ReducedElements:
             if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
                 Ripley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on elements with reduced integration order possible.");  
             }  
             break;  
527          case FaceElements:          case FaceElements:
             if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
                 Ripley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
             } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
                 Ripley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on face elements possible.");  
             }  
             break;  
528          case ReducedFaceElements:          case ReducedFaceElements:
529              if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {              return true;
                 Ripley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on face elements with reduced integration order possible.");  
             }  
             break;  
         case Points:  
             if (target.getFunctionSpace().getTypeCode()==Points) {  
                 Ripley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on points possible.");  
             }  
             break;  
530          case DegreesOfFreedom:          case DegreesOfFreedom:
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case ReducedDegreesOfFreedom:  
                 case DegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     break;  
   
                 case Nodes:  
                 case ReducedNodes:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in);  
                         temp.expand();  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in2);  
                     } else {  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     }  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, continuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in,&_target);  
                     }  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, continuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements,&_in,&_target);  
                     }  
                     break;  
                 case Points:  
                     if (getMPISize()>1) {  
                         //escript::Data temp=escript::Data(in, continuousFunction(*this) );  
                         //escriptDataC _in2 = temp.getDataC();  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes,mesh->Points,&_in,&_target);  
                     }  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
531          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
532              switch (target.getFunctionSpace().getTypeCode()) {          case Points:
533                  case Nodes:          case ReducedNodes:
534                      throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to mesh nodes.");              return false;
                     break;  
                 case ReducedNodes:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in);  
                         temp.expand();  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in2);  
                     } else {  
                         Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     }  
                     break;  
                 case DegreesOfFreedom:  
                     throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
                 case ReducedDegreesOfFreedom:  
                     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);  
                     break;  
                 case Elements:  
                 case ReducedElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in,&_target);  
                     }  
                     break;  
                 case FaceElements:  
                 case ReducedFaceElements:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in,&_target);  
                     }  
                     break;  
                 case Points:  
                     if (getMPISize()>1) {  
                         escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );  
                         escriptDataC _in2 = temp.getDataC();  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Points,&_in2,&_target);  
                     } else {  
                         Ripley_Assemble_interpolate(m_nodes, mesh->Points,&_in,&_target);  
                     }  
                     break;  
                 default:  
                     stringstream temp;  
                     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
                     throw RipleyException(temp.str());  
             }  
             break;  
535          default:          default:
536              stringstream temp;              break;
             temp << "Interpolation On Domain: Ripley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
537      }      }
538  */      stringstream msg;
539      checkPasoError();      msg << "canTag(): Illegal function space type " << fsType << " on "
540            << getDescription();
541        throw RipleyException(msg.str());
542  }  }
543    
544  //  void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
 // copies the locations of sample points into arg:  
 //  
 void RipleyDomain::setToX(escript::Data& arg) const  
545  {  {
546      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));      IndexVector* target=NULL;
547      if (argDomain != *this)      dim_t num=0;
         throw RipleyException("setToX: Illegal domain of data point locations");  
548    
549      // do we need to interpolate?      switch(fsType) {
     if (arg.getFunctionSpace().getTypeCode()==Nodes) {  
         m_nodes->assembleCoordinates(arg);  
     } else {  
         escript::Data tmp_data=Vector(0.0, continuousFunction(*this), true);  
         m_nodes->assembleCoordinates(tmp_data);  
         // this is then interpolated onto arg:  
         interpolateOnDomain(arg, tmp_data);  
     }  
 }  
   
 //  
 // returns the normal vectors at the location of data points as a Data object  
 //  
 void RipleyDomain::setToNormal(escript::Data& normal) const  
 {  
 /*  
     const RipleyDomain& normalDomain=dynamic_cast<const RipleyDomain&>(*(normal.getFunctionSpace().getDomain()));  
     if (normalDomain!=*this)  
         throw RipleyException("Illegal domain of normal locations");  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _normal=normal.getDataC();  
     switch(normal.getFunctionSpace().getTypeCode()) {  
550          case Nodes:          case Nodes:
551              throw RipleyException("Ripley does not support surface normal vectors for nodes");              num=getNumNodes();
552          case ReducedNodes:              target=&m_nodeTags;
553              throw RipleyException("Ripley does not support surface normal vectors for reduced nodes");              break;
554          case Elements:          case Elements:
             throw RipleyException("Ripley does not support surface normal vectors for elements");  
555          case ReducedElements:          case ReducedElements:
556              throw RipleyException("Ripley does not support surface normal vectors for elements with reduced integration order");              num=getNumElements();
557          case FaceElements:              target=&m_elementTags;
             Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);  
558              break;              break;
559            case FaceElements:
560          case ReducedFaceElements:          case ReducedFaceElements:
561              Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);              num=getNumFaceElements();
562                target=&m_faceTags;
563              break;              break;
564          case Points:          default: {
565              throw RipleyException("Ripley does not support surface normal vectors for point elements");              stringstream msg;
566          case DegreesOfFreedom:              msg << "setTags(): not implemented for "
567              throw RipleyException("Ripley does not support surface normal vectors for degrees of freedom.");                  << functionSpaceTypeAsString(fsType);
568          case ReducedDegreesOfFreedom:              throw RipleyException(msg.str());
569              throw RipleyException("Ripley does not support surface normal vectors for reduced degrees of freedom.");          }
570          default:      }
571              stringstream temp;      if (target->size() != num) {
572              temp << "Normal Vectors: Ripley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();          target->assign(num, -1);
             throw RipleyException(temp.str());  
573      }      }
 */  
 }  
   
 //  
 // interpolates data to other domain  
 //  
 void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  
 {  
     escript::const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();  
     const RipleyDomain* targetDomain=dynamic_cast<const RipleyDomain*>(targetDomain_p.get());  
     if (targetDomain!=this)  
         throw RipleyException("Illegal domain of interpolation target");  
574    
575      throw RipleyException("Ripley does not allow interpolation across domains yet");      escript::Data& mask=*const_cast<escript::Data*>(&cMask);
576    #pragma omp parallel for
577        for (index_t i=0; i<num; i++) {
578            if (mask.getSampleDataRO(i)[0] > 0) {
579                (*target)[i]=newTag;
580            }
581        }
582        updateTagsInUse(fsType);
583  }  }
584    
585  //  int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
 // calculates the integral of a function arg  
 //  
 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  
586  {  {
587  /*      switch(fsType) {
     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));  
     if (argDomain!=*this)  
         throw RipleyException("Illegal domain of integration kernel");  
   
     double blocktimer_start = blocktimer_time();  
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _temp;  
     escript::Data temp;  
     escriptDataC _arg=arg.getDataC();  
     switch(arg.getFunctionSpace().getTypeCode()) {  
588          case Nodes:          case Nodes:
589          case ReducedNodes:              if (m_nodeTags.size() > sampleNo)
590          case DegreesOfFreedom:                  return m_nodeTags[sampleNo];
         case ReducedDegreesOfFreedom:  
             temp=escript::Data( arg, escript::function(*this) );  
             _temp=temp.getDataC();  
             Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_temp,&integrals[0]);  
591              break;              break;
592          case Elements:          case Elements:
593          case ReducedElements:          case ReducedElements:
594              Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_arg,&integrals[0]);              if (m_elementTags.size() > sampleNo)
595                    return m_elementTags[sampleNo];
596              break;              break;
597          case FaceElements:          case FaceElements:
598          case ReducedFaceElements:          case ReducedFaceElements:
599              Ripley_Assemble_integrate(m_nodes,mesh->FaceElements,&_arg,&integrals[0]);              if (m_faceTags.size() > sampleNo)
600                    return m_faceTags[sampleNo];
601              break;              break;
602          case Points:          default: {
603              throw RipleyException("Integral of data on points is not supported.");              stringstream msg;
604          default:              msg << "getTagFromSampleNo(): not implemented for "
605              stringstream temp;                  << functionSpaceTypeAsString(fsType);
606              temp << "Integrals: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();              throw RipleyException(msg.str());
607              throw RipleyException(temp.str());          }
608      }      }
609      blocktimer_increment("integrate()", blocktimer_start);      return -1;
 */  
610  }  }
611    
612  //  int RipleyDomain::getNumberOfTagsInUse(int fsType) const
 // calculates the gradient of arg  
 //  
 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const  
613  {  {
614  /*      switch(fsType) {
     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));  
     if (argDomain!=*this)  
         throw RipleyException("Illegal domain of gradient argument");  
     const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(*(grad.getFunctionSpace().getDomain()));  
     if (gradDomain!=*this)  
         throw RipleyException("Illegal domain of gradient");  
   
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC _grad=grad.getDataC();  
     escriptDataC nodeDataC;  
     escript::Data temp;  
     if (getMPISize()>1) {  
         if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {  
             temp=escript::Data(arg, continuousFunction(*this) );  
             nodeDataC = temp.getDataC();  
         } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {  
             temp=escript::Data(arg, reducedContinuousFunction(*this) );  
             nodeDataC = temp.getDataC();  
         } else {  
             nodeDataC = arg.getDataC();  
         }  
     } else {  
         nodeDataC = arg.getDataC();  
     }  
     switch(grad.getFunctionSpace().getTypeCode()) {  
615          case Nodes:          case Nodes:
616              throw RipleyException("Gradient at nodes is not supported.");              return m_nodeTagsInUse.size();
         case ReducedNodes:  
             throw RipleyException("Gradient at reduced nodes is not supported.");  
617          case Elements:          case Elements:
             Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);  
             break;  
618          case ReducedElements:          case ReducedElements:
619              Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);              return m_elementTagsInUse.size();
             break;  
620          case FaceElements:          case FaceElements:
             Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);  
             break;  
621          case ReducedFaceElements:          case ReducedFaceElements:
622              Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);              return m_faceTagsInUse.size();
623              break;          default: {
624          case Points:              stringstream msg;
625              throw RipleyException("Gradient at points is not supported");              msg << "getNumberOfTagsInUse(): not implemented for "
626          case DegreesOfFreedom:                  << functionSpaceTypeAsString(fsType);
627              throw RipleyException("Gradient at degrees of freedom is not supported");              throw RipleyException(msg.str());
628          case ReducedDegreesOfFreedom:          }
             throw RipleyException("Gradient at reduced degrees of freedom is not supported");  
         default:  
             stringstream temp;  
             temp << "Gradient: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
629      }      }
 */  
630  }  }
631    
632  //  const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
 // returns the size of elements  
 //  
 void RipleyDomain::setToSize(escript::Data& size) const  
633  {  {
634  /*      switch(fsType) {
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC tmp=size.getDataC();  
     switch(size.getFunctionSpace().getTypeCode()) {  
635          case Nodes:          case Nodes:
636              throw RipleyException("Size of nodes is not supported");              return &m_nodeTagsInUse[0];
         case ReducedNodes:  
             throw RipleyException("Size of reduced nodes is not supported");  
637          case Elements:          case Elements:
             Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);  
             break;  
638          case ReducedElements:          case ReducedElements:
639              Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);              return &m_elementTagsInUse[0];
             break;  
640          case FaceElements:          case FaceElements:
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
             break;  
641          case ReducedFaceElements:          case ReducedFaceElements:
642              Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);              return &m_faceTagsInUse[0];
643              break;          default: {
644          case Points:              stringstream msg;
645              throw RipleyException("Size of point elements is not supported");              msg << "borrowListOfTagsInUse(): not implemented for "
646          case DegreesOfFreedom:                  << functionSpaceTypeAsString(fsType);
647              throw RipleyException("Size of degrees of freedom is not supported");              throw RipleyException(msg.str());
648          case ReducedDegreesOfFreedom:          }
             throw RipleyException("Size of reduced degrees of freedom is not supported");  
         default:  
             stringstream temp;  
             temp << "Element size: Ripley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
649      }      }
 */  
650  }  }
651    
652  //  void RipleyDomain::Print_Mesh_Info(const bool full) const
 // sets the location of nodes  
 //  
 void RipleyDomain::setNewX(const escript::Data& new_x)  
653  {  {
654      const RipleyDomain& newDomain=dynamic_cast<const RipleyDomain&>(*(new_x.getFunctionSpace().getDomain()));      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
655      if (newDomain!=*this)          << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
656          throw RipleyException("Illegal domain of new node locations");      cout << "Number of dimensions: " << m_numDim << endl;
657    
658      if (new_x.getFunctionSpace() == continuousFunction(*this)) {      // write tags
659          m_nodes->setCoordinates(new_x);      if (m_tagMap.size() > 0) {
660      } else {          cout << "Tags:" << endl;
661          escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this) );          TagMap::const_iterator it;
662          m_nodes->setCoordinates(new_x_inter);          for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
663                cout << "  " << setw(5) << it->second << " "
664                    << it->first << endl;
665            }
666      }      }
667  }  }
668    
669  bool RipleyDomain::ownSample(int fs_code, index_t id) const  int RipleyDomain::getSystemMatrixTypeId(const int solver,
670            const int preconditioner, const int package, const bool symmetry) const
671  {  {
672  #ifdef ESYS_MPI      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
673      index_t myFirstNode, myLastNode, k;              package, symmetry, m_mpiInfo);
     if (fs_code == ReducedNodes) {  
         myFirstNode = m_nodes->getFirstReducedNode();  
         myLastNode = m_nodes->getLastReducedNode();  
         k = m_nodes->getIndexForGlobalReducedNode(id);  
     } else {  
         myFirstNode = m_nodes->getFirstNode();  
         myLastNode = m_nodes->getLastNode();  
         k = m_nodes->getIndexForGlobalNode(id);  
     }  
     return static_cast<bool>((myFirstNode <= k) && (k < myLastNode));  
 #endif  
     return true;  
674  }  }
675    
676    int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
677            const int package, const bool symmetry) const
678    {
679        return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
680                package, symmetry, m_mpiInfo);
681    }
682    
 //  
 // creates a SystemMatrixAdapter stiffness matrix and initializes it with zeros  
 //  
683  escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,  escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
684          const escript::FunctionSpace& row_functionspace,          const escript::FunctionSpace& row_functionspace,
685          const int column_blocksize,          const int column_blocksize,
# Line 1135  escript::ASM_ptr RipleyDomain::newSystem Line 690  escript::ASM_ptr RipleyDomain::newSystem
690      bool reduceColOrder=false;      bool reduceColOrder=false;
691      // is the domain right?      // is the domain right?
692      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
693      if (row_domain!=*this)      if (row_domain != *this)
694          throw RipleyException("Domain of row function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
695      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
696      if (col_domain!=*this)      if (col_domain != *this)
697          throw RipleyException("Domain of columnn function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
698      // is the function space type right      // is the function space type right?
699      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
700          reduceRowOrder=true;          reduceRowOrder=true;
701      } else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
702          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
703      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
704          reduceColOrder=true;          reduceColOrder=true;
705      } else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
706          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
707    
708      // generate matrix:      // generate matrix
709      Paso_SystemMatrixPattern* fsystemMatrixPattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
710      checkPasoError();      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
711      Paso_SystemMatrix* fsystemMatrix = Paso_SystemMatrix_alloc(type,              row_blocksize, column_blocksize, FALSE);
712              fsystemMatrixPattern, row_blocksize, column_blocksize, FALSE);      paso::checkPasoError();
713      checkPasoError();      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
714      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);                  row_functionspace, column_blocksize, column_functionspace));
715      SystemMatrixAdapter* sma = new SystemMatrixAdapter(fsystemMatrix,      return sma;
             row_blocksize, row_functionspace, column_blocksize, column_functionspace);  
     return escript::ASM_ptr(sma);  
716  }  }
717    
718  //  void RipleyDomain::addPDEToSystem(
719  // creates a TransportProblemAdapter          escript::AbstractSystemMatrix& mat, escript::Data& rhs,
720  //          const escript::Data& A, const escript::Data& B, const escript::Data& C,
721  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,          const escript::Data& D, const escript::Data& X, const escript::Data& Y,
722          const int blocksize, const escript::FunctionSpace& functionspace,          const escript::Data& d, const escript::Data& y,
723          const int type) const          const escript::Data& d_contact, const escript::Data& y_contact,
724            const escript::Data& d_dirac,const escript::Data& y_dirac) const
725  {  {
726      int reduceOrder=0;      if (!d_contact.isEmpty() || !y_contact.isEmpty())
727      // is the domain right?          throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
     const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));  
     if (domain!=*this)  
         throw RipleyException("Domain of function space does not match the domain of transport problem generator");  
     // is the function space type right  
     if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceOrder=1;  
     } else if (functionspace.getTypeCode()!=DegreesOfFreedom)  
         throw RipleyException("Illegal function space type for system matrix rows");  
728    
729      // generate matrix      paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
730      Paso_SystemMatrixPattern* fsystemMatrixPattern = getPattern(reduceOrder, reduceOrder);      if (!sma)
731      Paso_TransportProblem* transportProblem = Paso_TransportProblem_alloc(          throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
732              useBackwardEuler, fsystemMatrixPattern, blocksize);  
733      checkPasoError();      if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
734      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);          throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
735      escript::AbstractTransportProblem* atp=new TransportProblemAdapter(  
736              transportProblem, useBackwardEuler, blocksize, functionspace);      int fsType=UNKNOWN;
737      return escript::ATP_ptr(atp);      fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
738  }      fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
739        fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
740        fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
741        fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
742        fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
743        fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
744        fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
745        fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
746        fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
747        if (fsType==UNKNOWN)
748            return;
749    
750  // returns true if data on the function space is considered as being cell      const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
751  // centered  
752  bool RipleyDomain::isCellOriented(int functionSpaceCode) const      //TODO: more input param checks (shape, etc)
753  {  
754      switch(functionSpaceCode) {      Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
755          case Nodes:  
756          case DegreesOfFreedom:      if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
757          case ReducedDegreesOfFreedom:          throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
758              return false;  
759          case Elements:      const int numEq=S->logical_row_block_size;
760          case FaceElements:      const int numComp=S->logical_col_block_size;
761          case Points:  
762          case ReducedElements:      if (numEq != numComp)
763          case ReducedFaceElements:          throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
764              return true;      //TODO: more system matrix checks
765          default:  
766              stringstream temp;      if (numEq==1) {
767              temp << "Ripley does not know anything about function space type " << functionSpaceCode;          if (reducedOrder) {
768              throw RipleyException(temp.str());              assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
769                assemblePDEBoundarySingleReduced(S, rhs, d, y);
770            } else {
771                assemblePDESingle(S, rhs, A, B, C, D, X, Y);
772                assemblePDEBoundarySingle(S, rhs, d, y);
773            }
774        } else {
775            if (reducedOrder) {
776                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
777                assemblePDEBoundarySystemReduced(S, rhs, d, y);
778            } else {
779                assemblePDESystem(S, rhs, A, B, C, D, X, Y);
780                assemblePDEBoundarySystem(S, rhs, d, y);
781            }
782      }      }
     return false;  
783  }  }
784    
785  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
786            const escript::Data& Y, const escript::Data& y,
787            const escript::Data& y_contact, const escript::Data& y_dirac) const
788  {  {
789     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]      if (!y_contact.isEmpty())
790      class 1: DOF <-> Nodes          throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
     class 2: ReducedDOF <-> ReducedNodes  
     class 3: Points  
     class 4: Elements  
     class 5: ReducedElements  
     class 6: FaceElements  
     class 7: ReducedFaceElements  
     class 8: ContactElementZero <-> ContactElementOne  
     class 9: ReducedContactElementZero <-> ReducedContactElementOne  
   
    There is also a set of lines. Interpolation is possible down a line but not between lines.  
    class 1 and 2 belong to all lines so aren't considered.  
     line 0: class 3  
     line 1: class 4,5  
     line 2: class 6,7  
     line 3: class 8,9  
791    
792      For classes with multiple members (eg class 2) we have vars to record if      if (rhs.isEmpty()) {
793      there is at least one instance. e.g. hasnodes is true if we have at least          if (!X.isEmpty() || !Y.isEmpty())
794      one instance of Nodes.              throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
795      */          else
796      if (fs.empty())              return;
         return false;  
     vector<int> hasclass(10);  
     vector<int> hasline(4);  
     bool hasnodes=false;  
     bool hasrednodes=false;  
     for (int i=0;i<fs.size();++i) {  
         switch (fs[i]) {  
             case Nodes: hasnodes=true; // no break is deliberate  
             case DegreesOfFreedom:  
                 hasclass[1]=1;  
                 break;  
             case ReducedNodes: hasrednodes=true; // no break is deliberate  
             case ReducedDegreesOfFreedom:  
                 hasclass[2]=1;  
                 break;  
             case Points:  
                 hasline[0]=1;  
                 hasclass[3]=1;  
                 break;  
             case Elements:  
                 hasclass[4]=1;  
                 hasline[1]=1;  
                 break;  
             case ReducedElements:  
                 hasclass[5]=1;  
                 hasline[1]=1;  
                 break;  
             case FaceElements:  
                 hasclass[6]=1;  
                 hasline[2]=1;  
                 break;  
             case ReducedFaceElements:  
                 hasclass[7]=1;  
                 hasline[2]=1;  
                 break;  
             default:  
                 return false;  
         }  
797      }      }
     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];  
798    
799      // fail if we have more than one leaf group      if (rhs.getDataPointSize() == 1) {
800      // = there are at least two branches we can't interpolate between          assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
801      if (totlines>1)          assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
802          return false;      } else {
803      else if (totlines==1) {          assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
804          // we have points          assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
805          if (hasline[0]==1) {      }
806              resultcode=Points;  }
807          } else if (hasline[1]==1) {  
808              if (hasclass[5]==1) {  void RipleyDomain::setNewX(const escript::Data& arg)
809                  resultcode=ReducedElements;  {
810              } else {      throw RipleyException("setNewX(): Operation not supported");
811                  resultcode=Elements;  }
812              }  
813          } else if (hasline[2]==1) {  //protected
814              if (hasclass[7]==1) {  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
815                  resultcode=ReducedFaceElements;  {
816              } else {      const dim_t numComp = in.getDataPointSize();
817                  resultcode=FaceElements;      const dim_t dpp = in.getNumDataPointsPerSample();
818              }      out.requireWrite();
819          } else { // so we must be in line3  #pragma omp parallel for
820              throw RipleyException("Programmer Error - choosing between contact elements - we should never get here");      for (index_t i=0; i<in.getNumSamples(); i++) {
821          }          const double* src = in.getSampleDataRO(i);
822      } else { // totlines==0          double* dest = out.getSampleDataRW(i);
823          if (hasclass[2]==1) {          for (index_t c=0; c<numComp; c++) {
824              // something from class 2              double res=0.;
825              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              for (index_t q=0; q<dpp; q++)
826          } else { // something from class 1                  res+=src[c+q*numComp];
827              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              *dest++ = res/dpp;
828          }          }
829      }      }
     return true;  
830  }  }
831    
832  bool RipleyDomain::probeInterpolationOnDomain(int functionSpaceType_source,  //protected
833                                               int functionSpaceType_target) const  void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
834  {  {
835      if (functionSpaceType_target != Nodes &&      const dim_t numComp = in.getDataPointSize();
836              functionSpaceType_target != ReducedNodes &&      out.requireWrite();
837              functionSpaceType_target != ReducedDegreesOfFreedom &&  #pragma omp parallel for
838              functionSpaceType_target != DegreesOfFreedom &&      for (index_t i=0; i<in.getNumSamples(); i++) {
839              functionSpaceType_target != Elements &&          const double* src = in.getSampleDataRO(i);
840              functionSpaceType_target != ReducedElements &&          copy(src, src+numComp, out.getSampleDataRW(i));
             functionSpaceType_target != FaceElements &&  
             functionSpaceType_target != ReducedFaceElements &&  
             functionSpaceType_target != Points) {  
         stringstream temp;  
         temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_target;  
         throw RipleyException(temp.str());  
841      }      }
842    }
843    
844      switch (functionSpaceType_source) {  //protected
845    void RipleyDomain::updateTagsInUse(int fsType) const
846    {
847        IndexVector* tagsInUse=NULL;
848        const IndexVector* tags=NULL;
849        switch(fsType) {
850          case Nodes:          case Nodes:
851          case DegreesOfFreedom:              tags=&m_nodeTags;
852              return true;              tagsInUse=&m_nodeTagsInUse;
853          case ReducedNodes:              break;
         case ReducedDegreesOfFreedom:  
             return (functionSpaceType_target != Nodes &&  
                     functionSpaceType_target != DegreesOfFreedom);  
854          case Elements:          case Elements:
             return (functionSpaceType_target==Elements ||  
                     functionSpaceType_target==ReducedElements);  
855          case ReducedElements:          case ReducedElements:
856              return (functionSpaceType_target==ReducedElements);              tags=&m_elementTags;
857                tagsInUse=&m_elementTagsInUse;
858                break;
859          case FaceElements:          case FaceElements:
             return (functionSpaceType_target==FaceElements ||  
                     functionSpaceType_target==ReducedFaceElements);  
860          case ReducedFaceElements:          case ReducedFaceElements:
861              return (functionSpaceType_target==ReducedFaceElements);              tags=&m_faceTags;
862          case Points:              tagsInUse=&m_faceTagsInUse;
863              return (functionSpaceType_target==Points);              break;
   
864          default:          default:
865              stringstream temp;              return;
866              temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_source;      }
867              throw RipleyException(temp.str());  
868        // gather global unique tag values from tags into tagsInUse
869        tagsInUse->clear();
870        index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
871    
872        while (true) {
873            // find smallest value bigger than lastFoundValue
874            minFoundValue = INDEX_T_MAX;
875    #pragma omp parallel private(local_minFoundValue)
876            {
877                local_minFoundValue = minFoundValue;
878    #pragma omp for schedule(static) nowait
879                for (size_t i = 0; i < tags->size(); i++) {
880                    const index_t v = (*tags)[i];
881                    if ((v > lastFoundValue) && (v < local_minFoundValue))
882                        local_minFoundValue = v;
883                }
884    #pragma omp critical
885                {
886                    if (local_minFoundValue < minFoundValue)
887                        minFoundValue = local_minFoundValue;
888                }
889            }
890    #ifdef ESYS_MPI
891            local_minFoundValue = minFoundValue;
892            MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
893    #endif
894    
895            // if we found a new value add it to the tagsInUse vector
896            if (minFoundValue < INDEX_T_MAX) {
897                tagsInUse->push_back(minFoundValue);
898                lastFoundValue = minFoundValue;
899            } else
900                break;
901      }      }
     return false;  
902  }  }
903    
904  bool RipleyDomain::probeInterpolationACross(int functionSpaceType_source,  //protected
905          const AbstractDomain& targetDomain, int functionSpaceType_target) const  Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
906            const IndexVector& index, const dim_t M, const dim_t N) const
907  {  {
908      return false;      // paso will manage the memory
909        index_t* indexC = MEMALLOC(index.size(), index_t);
910        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
911        copy(index.begin(), index.end(), indexC);
912        copy(ptr.begin(), ptr.end(), ptrC);
913        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
914  }  }
915    
916  bool RipleyDomain::operator==(const AbstractDomain& other) const  //protected
917    Paso_Pattern* RipleyDomain::createMainPattern() const
918  {  {
919      const RipleyDomain* temp=dynamic_cast<const RipleyDomain*>(&other);      IndexVector ptr(1,0);
920      if (temp != NULL) {      IndexVector index;
921          return (m_name==temp->m_name && m_nodes==temp->m_nodes &&  
922                  m_elements==temp->m_elements &&      for (index_t i=0; i<getNumDOF(); i++) {
923                  m_faceElements==temp->m_faceElements &&          // add the DOF itself
924                  m_points==temp->m_points);          index.push_back(i);
925            const dim_t num=insertNeighbourNodes(index, i);
926            ptr.push_back(ptr.back()+num+1);
927      }      }
928      return false;  
929        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
930  }  }
931    
932  bool RipleyDomain::operator!=(const AbstractDomain& other) const  //protected
933    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
934            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
935  {  {
936      return !(operator==(other));      IndexVector ptr(1,0);
937        IndexVector index;
938        for (index_t i=0; i<getNumDOF(); i++) {
939            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
940            ptr.push_back(ptr.back()+colIndices[i].size());
941        }
942    
943        const dim_t M=ptr.size()-1;
944        *colPattern=createPasoPattern(ptr, index, M, N);
945    
946        IndexVector rowPtr(1,0);
947        IndexVector rowIndex;
948        for (dim_t id=0; id<N; id++) {
949            dim_t n=0;
950            for (dim_t i=0; i<M; i++) {
951                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
952                    if (index[j]==id) {
953                        rowIndex.push_back(i);
954                        n++;
955                        break;
956                    }
957                }
958            }
959            rowPtr.push_back(rowPtr.back()+n);
960        }
961    
962        // M and N are now reversed
963        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
964  }  }
965    
966  int RipleyDomain::getSystemMatrixTypeId(const int solver,  //protected
967          const int preconditioner, const int package, const bool symmetry) const  void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
968  {         const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
969      int out=Paso_SystemMatrix_getSystemMatrixTypeId(         dim_t num_Sol, const vector<double>& array) const
970              SystemMatrixAdapter::mapOptionToPaso(solver),  {
971              SystemMatrixAdapter::mapOptionToPaso(preconditioner),      const dim_t numMyCols = mat->pattern->mainPattern->numInput;
972              SystemMatrixAdapter::mapOptionToPaso(package),      const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
973              symmetry?1:0, m_mpiInfo);      const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
974      checkPasoError();      const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
975      return out;  
976        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
977        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
978        double* mainBlock_val = mat->mainBlock->val;
979        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
980        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
981        double* col_coupleBlock_val = mat->col_coupleBlock->val;
982        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
983        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
984        double* row_coupleBlock_val = mat->row_coupleBlock->val;
985    
986        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
987            // down columns of array
988            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
989                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
990                // only look at the matrix rows stored on this processor
991                if (i_row < numMyRows) {
992                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
993                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
994                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
995                            if (i_col < numMyCols) {
996                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
997                                    if (mainBlock_index[k] == i_col) {
998                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
999                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1000                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1001                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1002                                                mainBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1003                                            }
1004                                        }
1005                                        break;
1006                                    }
1007                                }
1008                            } else {
1009                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1010                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
1011                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1012                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1013                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1014                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1015                                                col_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1016                                            }
1017                                        }
1018                                        break;
1019                                    }
1020                                }
1021                            }
1022                        }
1023                    }
1024                } else {
1025                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1026                        // across rows of array
1027                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1028                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1029                            if (i_col < numMyCols) {
1030                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1031                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1032                                {
1033                                    if (row_coupleBlock_index[k] == i_col) {
1034                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1035                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1036                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1037                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1038                                            row_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1039                                            }
1040                                        }
1041                                        break;
1042                                    }
1043                                }
1044                            }
1045                        }
1046                    }
1047                }
1048            }
1049        }
1050  }  }
1051    
1052  int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,  //
1053          const int package, const bool symmetry) const  // the methods that follow have to be implemented by the subclasses
1054    //
1055    
1056    string RipleyDomain::getDescription() const
1057  {  {
1058      int out=Paso_TransportProblem_getTypeId(      throw RipleyException("getDescription() not implemented");
             SystemMatrixAdapter::mapOptionToPaso(solver),  
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
1059  }  }
1060    
1061  escript::Data RipleyDomain::getX() const  void RipleyDomain::write(const string& filename) const
1062  {  {
1063      return continuousFunction(*this).getX();      throw RipleyException("write() not implemented");
1064  }  }
1065    
1066  escript::Data RipleyDomain::getNormal() const  void RipleyDomain::dump(const string& filename) const
1067  {  {
1068      return functionOnBoundary(*this).getNormal();      throw RipleyException("dump() not implemented");
1069  }  }
1070    
1071  escript::Data RipleyDomain::getSize() const  const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1072  {  {
1073      return escript::function(*this).getSize();      throw RipleyException("borrowSampleReferenceIDs() not implemented");
1074  }  }
1075    
1076  const int* RipleyDomain::borrowSampleReferenceIDs(int functionSpaceType) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1077  {  {
1078      switch (functionSpaceType) {      throw RipleyException("interpolateACross() not implemented");
         case Nodes:  
             return &m_nodes->getIdVector()[0];  
         case ReducedNodes:  
             return &m_nodes->getReducedNodesId()[0];  
         case Elements:  
         case ReducedElements:  
             return &m_elements->getIdVector()[0];  
         case FaceElements:  
         case ReducedFaceElements:  
             return &m_faceElements->getIdVector()[0];  
         case Points:  
             return &m_points->getIdVector()[0];  
         case DegreesOfFreedom:  
             return &m_nodes->getDegreesOfFreedomId()[0];  
         case ReducedDegreesOfFreedom:  
             return &m_nodes->getReducedDegreesOfFreedomId()[0];  
         default:  
             stringstream msg;  
             msg << "borrowSampleReferenceIDs: Invalid function space type " << functionSpaceType << " for domain: " << getDescription();  
             throw RipleyException(msg.str());  
     }  
1079  }  }
1080    
1081  int RipleyDomain::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  bool RipleyDomain::probeInterpolationACross(int fsType_source,
1082            const escript::AbstractDomain&, int fsType_target) const
1083  {  {
1084      throw RipleyException("not implemented");      throw RipleyException("probeInterpolationACross() not implemented");
1085  }  }
1086    
1087    void RipleyDomain::setToNormal(escript::Data& normal) const
 void RipleyDomain::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  
1088  {  {
1089      switch (functionSpaceType) {      throw RipleyException("setToNormal() not implemented");
         case Nodes:  
             m_nodes->setTags(newTag, mask);  
             break;  
         case ReducedNodes:  
             throw RipleyException("ReducedNodes does not support tags");  
         case DegreesOfFreedom:  
             throw RipleyException("DegreesOfFreedom does not support tags");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("ReducedDegreesOfFreedom does not support tags");  
         case Elements:  
             m_elements->setTags(newTag, mask);  
             break;  
         case ReducedElements:  
             m_elements->setTags(newTag, mask);  
             break;  
         case FaceElements:  
             m_faceElements->setTags(newTag, mask);  
             break;  
         case ReducedFaceElements:  
             m_faceElements->setTags(newTag, mask);  
             break;  
         case Points:  
             m_points->setTags(newTag, mask);  
             break;  
         default:  
             stringstream msg;  
             msg << "Ripley does not know anything about function space type " << functionSpaceType;  
             throw RipleyException(msg.str());  
     }  
1090  }  }
1091    
1092  void RipleyDomain::setTagMap(const string& name, int tag)  void RipleyDomain::setToSize(escript::Data& size) const
1093  {  {
1094      m_tagMap[name] = tag;      throw RipleyException("setToSize() not implemented");
1095  }  }
1096    
1097  int RipleyDomain::getTag(const string& name) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
1098  {  {
1099      return m_tagMap.find(name)->second;      throw RipleyException("ownSample() not implemented");
1100  }  }
1101    
1102  bool RipleyDomain::isValidTagName(const string& name) const  void RipleyDomain::addPDEToTransportProblem(
1103            escript::AbstractTransportProblem& tp,
1104            escript::Data& source, const escript::Data& M,
1105            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1106            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1107            const escript::Data& d, const escript::Data& y,
1108            const escript::Data& d_contact, const escript::Data& y_contact,
1109            const escript::Data& d_dirac, const escript::Data& y_dirac) const
1110  {  {
1111      return (m_tagMap.find(name)!=m_tagMap.end());      throw RipleyException("addPDEToTransportProblem() not implemented");
1112  }  }
1113    
1114  string RipleyDomain::showTagNames() const  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1115            const int blocksize, const escript::FunctionSpace& functionspace,
1116            const int type) const
1117  {  {
1118      stringstream ret;      throw RipleyException("newTransportProblem() not implemented");
     TagMap::const_iterator it;  
     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
         if (it!=m_tagMap.begin()) ret << ", ";  
         ret << it->first;  
     }  
     return ret.str();  
1119  }  }
1120    
1121  int RipleyDomain::getNumberOfTagsInUse(int functionSpaceCode) const  Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
1122                                                       bool reducedColOrder) const
1123  {  {
1124      dim_t numTags=0;      throw RipleyException("getPattern() not implemented");
     switch (functionSpaceCode) {  
         case Nodes:  
             numTags=m_nodes->getNumberOfTagsInUse();  
             break;  
         case ReducedNodes:  
             throw RipleyException("ReducedNodes does not support tags");  
         case DegreesOfFreedom:  
             throw RipleyException("DegreesOfFreedom does not support tags");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("ReducedDegreesOfFreedom does not support tags");  
         case Elements:  
         case ReducedElements:  
             numTags=m_elements->getNumberOfTagsInUse();  
             break;  
         case FaceElements:  
         case ReducedFaceElements:  
             numTags=m_faceElements->getNumberOfTagsInUse();  
             break;  
         case Points:  
             numTags=m_points->getNumberOfTagsInUse();  
             break;  
         default:  
             stringstream msg;  
             msg << "Ripley does not know anything about function space type " << functionSpaceCode;  
             throw RipleyException(msg.str());  
     }  
     return numTags;  
1125  }  }
1126    
1127  const int* RipleyDomain::borrowListOfTagsInUse(int functionSpaceCode) const  dim_t RipleyDomain::getNumDataPointsGlobal() const
1128  {  {
1129      switch (functionSpaceCode) {      throw RipleyException("getNumDataPointsGlobal() not implemented");
         case Nodes:  
             return &m_nodes->getTagsInUse()[0];  
         case ReducedNodes:  
             throw RipleyException("ReducedNodes does not support tags");  
         case DegreesOfFreedom:  
             throw RipleyException("DegreesOfFreedom does not support tags");  
         case ReducedDegreesOfFreedom:  
             throw RipleyException("ReducedDegreesOfFreedom does not support tags");  
         case Elements:  
         case ReducedElements:  
             return &m_elements->getTagsInUse()[0];  
         case FaceElements:  
         case ReducedFaceElements:  
             return &m_faceElements->getTagsInUse()[0];  
         case Points:  
             return &m_points->getTagsInUse()[0];  
         default:  
             stringstream msg;  
             msg << "Ripley does not know anything about function space type " << functionSpaceCode;  
             throw RipleyException(msg.str());  
     }  
1130  }  }
1131    
1132  bool RipleyDomain::canTag(int functionSpaceCode) const  IndexVector RipleyDomain::getNumNodesPerDim() const
1133  {  {
1134      switch(functionSpaceCode) {      throw RipleyException("getNumNodesPerDim() not implemented");
         case Nodes:  
         case Elements:  
         case ReducedElements:  
         case FaceElements:  
         case ReducedFaceElements:  
         case Points:  
             return true;  
     }  
     return false;  
1135  }  }
1136    
1137  escript::AbstractDomain::StatusType RipleyDomain::getStatus() const  IndexVector RipleyDomain::getNumElementsPerDim() const
1138  {  {
1139      return m_nodes->getStatus();      throw RipleyException("getNumElementsPerDim() not implemented");
1140  }  }
1141    
1142  void RipleyDomain::prepare(bool optimize)  IndexVector RipleyDomain::getNumFacesPerBoundary() const
1143  {  {
1144      resolveNodeIds();      throw RipleyException("getNumFacesPerBoundary() not implemented");
   
     // first step is to distribute the elements according to a global  
     // distribution of DOF  
     // - create dense labeling for the DOFs  
     dim_t newGlobalNumDOFs=m_nodes->createDenseDOFLabeling();  
   
     // - create a distribution of the global DOFs and determine  
     //   the MPI_rank controlling the DOFs on this processor  
     IndexVector distribution(m_mpiInfo->size+1);  
     Esys_MPIInfo_setDistribution(m_mpiInfo, 0, newGlobalNumDOFs-1, &distribution[0]);  
     checkPasoError();  
   
     // Now the mesh is redistributed according to the mpiRankOfDOF vector.  
     // This will redistribute the Nodes and Elements including overlap and will  
     // create an element coloring but will not create any mappings (see later  
     // in this function).  
     distributeByRankOfDOF(distribution);  
   
     // At this stage we are able to start an optimization of the DOF  
     // distribution using ParMetis. On return distribution is altered and new  
     // DOF IDs have been assigned.  
     if (optimize) {  
         if (m_mpiInfo->size > 1) {  
             optimizeDOFDistribution(distribution);  
             distributeByRankOfDOF(distribution);  
         }  
   
         // optimize the local labeling of the degrees of freedom  
         optimizeDOFLabeling(distribution);  
     }  
   
     // Rearrange elements in order to try and bring them closer to memory  
     // locations of the nodes (distributed shared memory).  
     optimizeElementOrdering();  
   
 /* useful DEBUG:  
     {  
         printf("prepare: global DOF : %d\n",newGlobalNumDOFs);  
         IndexPair range = m_nodes->getGlobalIdRange();  
         printf("prepare: global node id range = %d :%d\n",range.first,range.second);  
         range = m_nodes->getIdRange();  
         printf("prepare: local node id range = %d :%d\n",range.first,range.second);  
     }  
 */  
   
     // create the global indices  
     IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);  
     markNodes(maskReducedNodes, 0);  
     IndexVector indexReducedNodes = packMask(maskReducedNodes);  
   
     IndexVector nodeDistribution(m_mpiInfo->size+1);  
     m_nodes->createDenseNodeLabeling(nodeDistribution, distribution);  
     m_nodes->createDenseReducedDOFLabeling(maskReducedNodes);  
     m_nodes->createDenseReducedNodeLabeling(maskReducedNodes);  
   
     // create the missing mappings  
     m_nodes->createNodeFileMappings(indexReducedNodes, distribution,  
                                     nodeDistribution);  
   
     updateTagsInUse();  
1145  }  }
1146    
1147  void RipleyDomain::optimizeElementOrdering()  IndexVector RipleyDomain::getNodeDistribution() const
1148  {  {
1149      m_elements->optimizeOrdering();      throw RipleyException("getNodeDistribution() not implemented");
     m_faceElements->optimizeOrdering();  
     m_points->optimizeOrdering();  
1150  }  }
1151    
1152  void RipleyDomain::updateTagsInUse()  IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1153  {  {
1154      m_nodes->updateTagsInUse();      throw RipleyException("getNumSubdivisionsPerDim() not implemented");
     m_elements->updateTagsInUse();  
     m_faceElements->updateTagsInUse();  
     m_points->updateTagsInUse();  
1155  }  }
1156    
1157  void RipleyDomain::createColoring(const IndexVector &node_localDOF_map)  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1158  {  {
1159      m_elements->createColoring(node_localDOF_map);      throw RipleyException("getFirstCoordAndSpacing() not implemented");
     m_faceElements->createColoring(node_localDOF_map);  
     m_points->createColoring(node_localDOF_map);  
1160  }  }
1161    
1162  void RipleyDomain::distributeByRankOfDOF(const IndexVector &dofDistribution)  dim_t RipleyDomain::getNumFaceElements() const
1163  {  {
1164      RankVector mpiRankOfDOF = m_nodes->getOwnerOfDOFs(dofDistribution);      throw RipleyException("getNumFaceElements() not implemented");
   
     // First the elements are redistributed according to mpiRankOfDOF.  
     // At the input the Node tables refer to the local labeling of  
     // the nodes while at the output they refer to the global labeling  
     // which is rectified in the next step  
     m_elements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_faceElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
     m_points->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());  
   
     // resolve the node IDs  
     resolveNodeIds();  
   
     // create a local labeling of the DOFs  
     IndexPair dofRange = m_nodes->getDOFRange();  
     dim_t len = dofRange.second - dofRange.first + 1;  
   
     // local mask for used nodes  
     IndexVector tmp_node_localDOF_mask(len, -1);  
     IndexVector tmp_node_localDOF_map(m_nodes->getNumNodes(), -1);  
   
 #pragma omp parallel for schedule(static)  
     for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {  
 #ifdef BOUNDS_CHECK  
         if ((m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) >= len  
                 || (m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) < 0) {  
             printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);  
             exit(1);  
         }  
 #endif  
         tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first] = n;  
     }  
   
     dim_t numDOFs = 0;  
     for (dim_t n = 0; n < len; n++) {  
         if (tmp_node_localDOF_mask[n] >= 0) {  
             tmp_node_localDOF_mask[n] = numDOFs;  
             numDOFs++;  
         }  
     }  
 #pragma omp parallel for  
     for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {  
         tmp_node_localDOF_map[n] = tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first];  
     }  
     // create element coloring  
     createColoring(tmp_node_localDOF_map);  
1165  }  }
1166    
1167  /**************************************************************  dim_t RipleyDomain::getNumElements() const
1168     Check whether there is any node which has no vertex. In case  {
1169     such a node exists, we don't use ParMetis since it requires      throw RipleyException("getNumElements() not implemented");
    that every node has at least 1 vertex (at line 129 of file  
    "xyzpart.c" in parmetis 3.1.1, variable "nvtxs" would be 0 if  
    any node has no vertex).  
  **************************************************************/  
 #ifdef USE_PARMETIS  
 static int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank,  
                         const vector<dim_t> &distribution, MPI_Comm *comm)  
 {  
     dim_t i, len;  
     int ret_val = 1;  
   
     if (rank == 0) {  
         for (i = 0; i < mpiSize; i++) {  
             len = distribution[i + 1] - distribution[i];  
             if (len == 0) {  
                 ret_val = 0;  
                 break;  
             }  
         }  
     }  
     MPI_Bcast(&ret_val, 1, MPI_INTEGER, 0, *comm);  
     if (ret_val == 0)  
         printf("INFO: Not using ParMetis since some nodes have no vertex!\n");  
     return ret_val;  
1170  }  }
 #endif  
1171    
1172  void RipleyDomain::optimizeDOFDistribution(RankVector &distribution)  dim_t RipleyDomain::getNumNodes() const
1173  {  {
1174      const Esys_MPI_rank myRank = m_mpiInfo->rank;      throw RipleyException("getNumNodes() not implemented");
1175      dim_t mpiSize = m_mpiInfo->size;  }
     dim_t dim = getDim();  
   
     // first step is to distribute the elements according to a global X of DOF  
     const index_t myFirstVertex = distribution[myRank];  
     const index_t myLastVertex = distribution[myRank + 1];  
     const dim_t myNumVertices = myLastVertex - myFirstVertex;  
     const dim_t globalNumVertices = distribution[mpiSize];  
     dim_t len = 0;  
     for (dim_t p = 0; p < mpiSize; ++p)  
         len = MAX(len, distribution[p+1] - distribution[p]);  
     index_t *partition = TMPMEMALLOC(len, index_t);  
     float *xyz = TMPMEMALLOC(myNumVertices * dim, float);  
     if (!(Esys_checkPtr(partition) || Esys_checkPtr(xyz))) {  
         // set the coordinates  
         // it is assumed that at least one node on this processor provides a  
         // coordinate  
 #pragma omp parallel for  
         for (dim_t i = 0; i < m_nodes->getNumNodes(); ++i) {  
             dim_t k = m_nodes->getGlobalDegreesOfFreedom()[i] - myFirstVertex;  
             if ((k >= 0) && (k < myNumVertices)) {  
                 for (dim_t j = 0; j < dim; ++j)  
                     xyz[k * dim + j] = (float)(m_nodes->getCoordinates()[INDEX2(j, i, dim)]);  
             }  
         }  
   
         IndexMatrix indexList(myNumVertices);  
         // ksteube CSR of DOF IDs  
         // create the adjacency structure xadj and adjncy  
         // insert contributions from element matrices into columns of indexList  
         m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 myFirstVertex, myLastVertex);  
         m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 myFirstVertex, myLastVertex);  
         m_points->insertIntoIndexMatrixNoMainDiagonal(indexList,  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 m_nodes->getGlobalDegreesOfFreedom(),  
                 myFirstVertex, myLastVertex);  
   
         // create the local matrix pattern  
         Paso_Pattern *pattern = createPatternFromIndexMatrix(  
                 indexList, 0, myNumVertices, 0, globalNumVertices, 0);  
   
 #ifdef USE_PARMETIS  
         if (mpiSize > 1 && Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(m_mpiInfo->comm)) > 0) {  
             int wgtflag = 0;  
             int numflag = 0; /* pattern->ptr is C style: starting from 0 instead of 1 */  
             int ncon = 1;  
             int edgecut;  
             int options[2];  
             float *tpwgts = TMPMEMALLOC(ncon * mpiSize, float);  
             float *ubvec = TMPMEMALLOC(ncon, float);  
             for (dim_t i = 0; i < ncon * mpiSize; i++)  
                 tpwgts[i] = 1.0 / (float)mpiSize;  
             for (dim_t i = 0; i < ncon; i++)  
                 ubvec[i] = 1.05;  
             options[0] = 3;  
             options[1] = 15;  
             ParMETIS_V3_PartGeomKway(&distribution[0], pattern->ptr,  
                     pattern->index, NULL, NULL, &wgtflag, &numflag,  
                     &dim, xyz, &ncon, &mpiSize, tpwgts, ubvec, options,  
                     &edgecut, partition, // new CPU ownership of elements  
                     &m_mpiInfo->comm);  
             TMPMEMFREE(ubvec);  
             TMPMEMFREE(tpwgts);  
         } else {  
             for (dim_t i = 0; i < myNumVertices; ++i)  
                 partition[i] = 0; // CPU 0 owns it  
         }  
 #else  
         for (dim_t i = 0; i < myNumVertices; ++i)  
             partition[i] = myRank; // CPU myRank owns it  
 #endif  
   
         Paso_Pattern_free(pattern);  
   
         // create a new distribution and labeling of the DOF  
         RankVector newDistribution(mpiSize+1, 0);  
 #pragma omp parallel  
         {  
             RankVector loc_partition_count(mpiSize, 0);  
 #pragma omp for  
             for (dim_t i = 0; i < myNumVertices; ++i)  
                 loc_partition_count[partition[i]]++;  
 #pragma omp critical  
             {  
                 for (dim_t i = 0; i < mpiSize; ++i)  
                     newDistribution[i] += loc_partition_count[i];  
             }  
         }  
1176    
1177          // recvbuf will be the concatenation of each CPU's contribution  dim_t RipleyDomain::getNumDOF() const
1178          // to newDistribution  {
1179          dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);      throw RipleyException("getNumDOF() not implemented");
1180    }
1181    
1182  #ifdef ESYS_MPI  dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1183          MPI_Allgather(&newDistribution[0], mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, m_mpiInfo->comm);  {
1184  #else      throw RipleyException("insertNeighbourNodes() not implemented");
1185          for (dim_t i = 0; i < mpiSize; ++i)  }
             recvbuf[i] = newDistribution[i];  
 #endif  
         newDistribution[0] = 0;  
         IndexVector newGlobalDOFid(len);  
         for (Esys_MPI_rank rank = 0; rank < mpiSize; rank++) {  
             int c = 0;  
             for (dim_t i = 0; i < myRank; ++i)  
                 c += recvbuf[rank + mpiSize * i];  
             for (dim_t i = 0; i < myNumVertices; ++i) {  
                 if (rank == partition[i]) {  
                     newGlobalDOFid[i] = newDistribution[rank] + c;  
                     c++;  
                 }  
             }  
             for (dim_t i = myRank + 1; i < mpiSize; ++i)  
                 c += recvbuf[rank + mpiSize * i];  
             newDistribution[rank + 1] = newDistribution[rank] + c;  
         }  
         TMPMEMFREE(recvbuf);  
1186    
1187          // now the overlap needs to be created by sending the partition  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1188          // around  {
1189          m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);      throw RipleyException("assembleCoordinates() not implemented");
1190          for (dim_t i = 0; i < mpiSize+1; ++i)  }
             distribution[i] = newDistribution[i];  
     }  
     TMPMEMFREE(partition);  
     TMPMEMFREE(xyz);  
 }  
   
 void RipleyDomain::optimizeDOFLabeling(const IndexVector &distribution)  
 {  
     Esys_MPI_rank myRank = m_mpiInfo->rank;  
     index_t myFirstVertex = distribution[myRank];  
     index_t myLastVertex = distribution[myRank + 1];  
     dim_t myNumVertices = myLastVertex - myFirstVertex;  
     dim_t len = 0;  
     for (dim_t p = 0; p < m_mpiInfo->size; ++p)  
         len = MAX(len, distribution[p + 1] - distribution[p]);  
   
     IndexMatrix indexList(myNumVertices);  
   
     // create the adjacency structure xadj and adjncy  
     // insert contributions from element matrices into columns of indexList  
     m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
             m_nodes->getGlobalDegreesOfFreedom(),  
             m_nodes->getGlobalDegreesOfFreedom(),  
             myFirstVertex, myLastVertex);  
     m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,  
             m_nodes->getGlobalDegreesOfFreedom(),  
             m_nodes->getGlobalDegreesOfFreedom(),  
             myFirstVertex, myLastVertex);  
     m_points->insertIntoIndexMatrixNoMainDiagonal(indexList,  
             m_nodes->getGlobalDegreesOfFreedom(),  
             m_nodes->getGlobalDegreesOfFreedom(),  
             myFirstVertex, myLastVertex);  
   
     // create the local matrix pattern  
     Paso_Pattern *pattern = createPatternFromIndexMatrix(indexList, 0,  
             myNumVertices, myFirstVertex, myLastVertex, -myFirstVertex);  
   
     IndexVector newGlobalDOFid(len);  
     Paso_Pattern_reduceBandwidth(pattern, &newGlobalDOFid[0]);  
     Paso_Pattern_free(pattern);  
1191    
1192      Esys_MPIInfo_noError(m_mpiInfo);  void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1193      checkPasoError();  {
1194        throw RipleyException("assembleGradient() not implemented");
1195    }
1196    
1197      /* shift new labeling to create a global id */  void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1198  #pragma omp parallel for  {
1199      for (dim_t i = 0; i < myNumVertices; ++i)      throw RipleyException("assembleIntegrate() not implemented");
1200          newGlobalDOFid[i] += myFirstVertex;  }
1201    
1202      /* distribute new labeling to other processors */  void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1203      m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);          const escript::Data& A, const escript::Data& B, const escript::Data& C,
1204  #if 0          const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1205      for (i = 0; i < m_nodes->getNumNodes(); ++i)  {
1206          printf("%d ", m_nodes->getGlobalDegreesOfFreedom()[i]);      throw RipleyException("assemblePDESingle() not implemented");
     printf("\n");  
 #endif  
1207  }  }
1208    
1209  void RipleyDomain::resolveNodeIds()  void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1210          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1211  {  {
1212      // find the minimum and maximum node ID used by elements      throw RipleyException("assemblePDEBoundarySingle() not implemented");
1213      IndexPair minMaxId = m_elements->getNodeRange();  }
     index_t minId = minMaxId.first;  
     index_t maxId = minMaxId.second;  
     minMaxId = m_faceElements->getNodeRange();  
     minId = MIN(minId, minMaxId.first);  
     maxId = MAX(maxId, minMaxId.second);  
     minMaxId = m_points->getNodeRange();  
     minId = MIN(minId, minMaxId.first);  
     maxId = MAX(maxId, minMaxId.second);  
1214    
1215  #ifdef Ripley_TRACE  void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1216      index_t globalMinId, globalMaxId;          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1217  #ifdef ESYS_MPI          const escript::Data& C, const escript::Data& D, const escript::Data& X,
1218      index_t idRange[2], globalIdRange[2];          const escript::Data& Y) const
1219      idRange[0] = -minId;  {
1220      idRange[1] = maxId;      throw RipleyException("assemblePDESingleReduced() not implemented");
1221      MPI_Allreduce(idRange, globalIdRange, 2, MPI_INT, MPI_MAX, m_mpiInfo->comm);  }
     globalMinId = -globalIdRange[0];  
     globalMaxId = globalIdRange[1];  
 #else  
     globalMinId = minId;  
     globalMaxId = maxId;  
 #endif // ESYS_MPI  
     cout << "Node id range used by elements is " << globalMinId << ":"  
         << globalMaxId << endl;  
 #endif // Ripley_TRACE  
   
     // this is only true if we have no elements at all  
     if (minId > maxId) {  
         maxId = -1;  
         minId = 0;  
     }  
   
     // allocate mappings for new local node labeling to global node labeling  
     // (newLocalToGlobalNodeLabels) and global node labeling to the new local  
     // node labeling. globalToNewLocalNodeLabels[i-minId] is the new local id  
     // of global node i  
     dim_t len = (maxId >= minId) ? maxId - minId + 1 : 0;  
     IndexVector globalToNewLocalNodeLabels(len, -1);  
   
     // mark the nodes referred by elements in globalToNewLocalNodeLabels  
     // which is currently used as a mask  
     markNodes(globalToNewLocalNodeLabels, minId);  
   
     // create a local labeling newLocalToGlobalNodeLabels of the local  
     // nodes by packing the mask globalToNewLocalNodeLabels  
     IndexVector newLocalToGlobalNodeLabels = packMask(globalToNewLocalNodeLabels);  
     const dim_t newNumNodes = newLocalToGlobalNodeLabels.size();  
   
     // invert the new labeling and shift the index newLocalToGlobalNodeLabels  
     // to global node IDs  
 #pragma omp parallel for schedule(static)  
     for (dim_t n = 0; n < newNumNodes; n++) {  
 #ifdef BOUNDS_CHECK  
         if (n >= len || n < 0) {  
             printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);  
             exit(1);  
         }  
         if (newLocalToGlobalNodeLabels[n] >= len || newLocalToGlobalNodeLabels[n] < 0) {  
             printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);  
             exit(1);  
         }  
 #endif  
         globalToNewLocalNodeLabels[newLocalToGlobalNodeLabels[n]] = n;  
         newLocalToGlobalNodeLabels[n] += minId;  
     }  
1222    
1223      // create a new node file  void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1224      m_nodes = m_nodes->gatherGlobal(newLocalToGlobalNodeLabels);        escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1225    {
1226        throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1227    }
1228    
1229      // relabel nodes of the elements  void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1230      relabelElementNodes(globalToNewLocalNodeLabels, minId);          const escript::Data& A, const escript::Data& B, const escript::Data& C,
1231            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1232    {
1233        throw RipleyException("assemblePDESystem() not implemented");
1234  }  }
1235    
1236  void RipleyDomain::relabelElementNodes(const IndexVector &newNode, index_t offset)  void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1237          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1238  {  {
1239      m_elements->relabelNodes(newNode, offset);      throw RipleyException("assemblePDEBoundarySystem() not implemented");
     m_faceElements->relabelNodes(newNode, offset);  
     m_points->relabelNodes(newNode, offset);  
1240  }  }
1241    
1242  void RipleyDomain::markNodes(IndexVector &mask, index_t offset)  void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1243            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1244            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1245            const escript::Data& Y) const
1246  {  {
1247      m_elements->markNodes(mask, offset);      throw RipleyException("assemblePDESystemReduced() not implemented");
     m_faceElements->markNodes(mask, offset);  
     m_points->markNodes(mask, offset);  
1248  }  }
1249    
1250  void RipleyDomain::createMappings(const IndexVector &dofDistribution, const IndexVector &nodeDistribution)  void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1251          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1252  {  {
1253      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);      throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1254      markNodes(maskReducedNodes, 0);  }
1255    
1256      IndexVector indexReducedNodes = packMask(maskReducedNodes);  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1257      m_nodes->createNodeFileMappings(indexReducedNodes, dofDistribution,  {
1258                                      nodeDistribution);      throw RipleyException("interpolateNodesOnElements() not implemented");
1259  }  }
1260    
1261  Paso_SystemMatrixPattern *RipleyDomain::makePattern(bool reduce_row_order, bool reduce_col_order) const  void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1262  {  {
1263      Paso_Connector *colConnector, *rowConnector;      throw RipleyException("interpolateNodesOnFaces() not implemented");
1264      NodeMapping_ptr colMap, rowMap;  }
     Paso_Distribution *colDistribution = NULL, *rowDistribution = NULL;  
1265    
1266      if (reduce_col_order) {  void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1267          colMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
1268          colDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("nodesToDOF() not implemented");
1269          colConnector = m_nodes->getReducedDegreesOfFreedomConnector();  }
     } else {  
         colMap = m_nodes->getDegreesOfFreedomMapping();  
         colDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         colConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
1270    
1271      if (reduce_row_order) {  void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1272          rowMap = m_nodes->getReducedDegreesOfFreedomMapping();  {
1273          rowDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();      throw RipleyException("dofToNodes() not implemented");
         rowConnector = m_nodes->getReducedDegreesOfFreedomConnector();  
     } else {  
         rowMap = m_nodes->getDegreesOfFreedomMapping();  
         rowDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         rowConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
   
     //  insert contributions from element matrices into columns in indexList  
     IndexMatrix indexList(rowMap->numTargets);  
     m_elements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
     m_faceElements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
     m_points->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);  
   
     // create patterns  
     Paso_Pattern *mainPattern = createPatternFromIndexMatrix(indexList, 0,  
             Paso_Distribution_getMyNumComponents(rowDistribution),  
             0, Paso_Distribution_getMyNumComponents(colDistribution), 0);  
     Paso_Pattern *colCouplePattern = createPatternFromIndexMatrix(indexList, 0,  
             Paso_Distribution_getMyNumComponents(rowDistribution),  
             Paso_Distribution_getMyNumComponents(colDistribution),  
             colMap->numTargets,  
             -Paso_Distribution_getMyNumComponents(colDistribution));  
     Paso_Pattern *rowCouplePattern = createPatternFromIndexMatrix(indexList,  
             Paso_Distribution_getMyNumComponents(rowDistribution),  
             rowMap->numTargets, 0,  
             Paso_Distribution_getMyNumComponents(colDistribution), 0);  
   
     // if everything is in order we can create the return value  
     Paso_SystemMatrixPattern *out = Paso_SystemMatrixPattern_alloc(  
             MATRIX_FORMAT_DEFAULT, rowDistribution, colDistribution,  
             mainPattern, colCouplePattern, rowCouplePattern,  
             colConnector, rowConnector);  
   
     // clean up  
     Paso_Pattern_free(mainPattern);  
     Paso_Pattern_free(colCouplePattern);  
     Paso_Pattern_free(rowCouplePattern);  
     Esys_MPIInfo_noError(m_mpiInfo);  
     return out;  
 }  
   
 Paso_SystemMatrixPattern *RipleyDomain::getPattern(bool reduce_row_order, bool reduce_col_order) const  
 {  
     return Paso_SystemMatrixPattern_getReference(makePattern(reduce_row_order, reduce_col_order));  
     /*  
     Paso_SystemMatrixPattern *out = NULL;  
   
     // make sure that the requested pattern is available  
     if (reduce_row_order) {  
         if (reduce_col_order) {  
             if (m_reducedReducedPattern == NULL)  
                 m_reducedReducedPattern = makePattern(reduce_row_order, reduce_col_order);  
         } else {  
             if (m_reducedFullPattern == NULL)  
                 m_reducedFullPattern = makePattern(reduce_row_order, reduce_col_order);  
         }  
     } else {  
         if (reduce_col_order) {  
             if (m_fullReducedPattern == NULL)  
                 m_fullReducedPattern = makePattern(reduce_row_order, reduce_col_order);  
         } else {  
             if (m_fullFullPattern == NULL)  
                 m_fullFullPattern = makePattern(reduce_row_order, reduce_col_order);  
         }  
     }  
     if (Ripley_noError()) {  
         if (reduce_row_order) {  
             if (reduce_col_order) {  
                 out = Paso_SystemMatrixPattern_getReference(m_reducedReducedPattern);  
             } else {  
                 out = Paso_SystemMatrixPattern_getReference(m_reducedFullPattern);  
             }  
         } else {  
             if (reduce_col_order) {  
                 out = Paso_SystemMatrixPattern_getReference(m_fullReducedPattern);  
             } else {  
                 out = Paso_SystemMatrixPattern_getReference(m_fullFullPattern);  
             }  
         }  
     }  
     return out;  
     */  
1274  }  }
1275    
1276  } // end of namespace ripley  } // end of namespace ripley

Legend:
Removed from v.3670  
changed lines
  Added in v.3785

  ViewVC Help
Powered by ViewVC 1.1.26