/[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 3746 by caltinay, Thu Dec 15 00:02:22 2011 UTC
# 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 Nodes:
65            case ReducedNodes:
66            case Elements:
67            case ReducedElements:
68            case FaceElements:
69            case ReducedFaceElements:
70            case Points:
71                return true;
72            default:
73                break;
74        }
75        return false;
76  }  }
77    
78  #ifdef ESYS_MPI  string RipleyDomain::functionSpaceTypeAsString(int fsType) const
 MPI_Comm  
 #else  
 unsigned int  
 #endif  
 RipleyDomain::getMPIComm() const  
79  {  {
80  #ifdef ESYS_MPI      switch (fsType) {
81      return m_mpiInfo->comm;          case Nodes: return "Ripley_Nodes";
82  #else          case ReducedNodes: return "Ripley_Reduced_Nodes";
83      return 0;          case Elements: return "Ripley_Elements";
84  #endif          case ReducedElements: return "Ripley_Reduced_Elements";
85            case FaceElements: return "Ripley_Face_Elements";
86            case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";
87            case Points: return "Ripley_Points";
88            default:
89                break;
90        }
91        return "Invalid function space type code";
92  }  }
93    
94  void RipleyDomain::write(const string& filename) const  pair<int,int> RipleyDomain::getDataShape(int fsType) const
95  {  {
96      if (m_mpiInfo->size > 1)      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
97          throw RipleyException("write: only single processor runs are supported.");      switch (fsType) {
98            case Nodes:
99      // open file          case ReducedNodes: //FIXME: reduced
100      ofstream f(filename.c_str());              return pair<int,int>(1, getNumNodes());
101      if (!f.is_open()) {          case Elements:
102          stringstream msg;              return pair<int,int>(ptsPerSample, getNumElements());
103          msg << "write: Opening file " << filename << " for writing failed.";          case FaceElements:
104          throw RipleyException(msg.str());              return pair<int,int>(ptsPerSample/2, getNumFaceElements());
105            case ReducedElements:
106                return pair<int,int>(1, getNumElements());
107            case ReducedFaceElements:
108                return pair<int,int>(1, getNumFaceElements());
109            case Points:
110                return pair<int,int>(1, 0); //FIXME: dirac
111            default:
112                break;
113      }      }
114    
115      // write header      stringstream msg;
116      f << m_name << endl;      msg << "getDataShape(): Unsupported function space type " << fsType
117            << " for " << getDescription();
118      f.precision(15);      throw RipleyException(msg.str());
119      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;  
         }  
     }  
120    
121      // write tags  string RipleyDomain::showTagNames() const
122      if (m_tagMap.size() > 0) {  {
123          f << "Tags" << endl;      stringstream ret;
124          TagMap::const_iterator it;      TagMap::const_iterator it;
125          for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {      for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
126              f << it->first << " " << it->second << endl;          if (it!=m_tagMap.begin()) ret << ", ";
127          }          ret << it->first;
128      }      }
129      f.close();      return ret.str();
 #ifdef Ripley_TRACE  
     cout << "mesh " << m_name << " has been written to file " << filename << endl;  
 #endif  
130  }  }
131    
132  void RipleyDomain::Print_Mesh_Info(const bool full) const  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
133  {  {
134      cout << "Ripley_PrintMesh_Info running on CPU " <<     /*
135          m_mpiInfo->rank << " of " << m_mpiInfo->size << endl;      The idea is to use equivalence classes (i.e. types which can be
136      cout << "\tMesh name '" << m_name << "'" << endl;      interpolated back and forth):
137        class 0: Nodes
138      // write nodes      class 1: ReducedNodes
139      int numDim = getDim();      class 2: Points
140      cout << "\tNodes: " << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;      class 3: Elements
141      if (full) {      class 4: ReducedElements
142          cout << "\t     Id   Tag  gDOF   gNI grDfI  grNI:  Coordinates" << endl;      class 5: FaceElements
143          cout.precision(15);      class 6: ReducedFaceElements
144          cout.setf(ios::scientific, ios::floatfield);  
145          for (int i = 0; i < m_nodes->getNumNodes(); i++) {      There is also a set of lines. Interpolation is possible down a line but not
146              cout << "\t  " << setw(5) << m_nodes->getIdVector()[i] << " "      between lines.
147                 << setw(5) << m_nodes->getTagVector()[i] << " "      class 0 and 1 belong to all lines so aren't considered.
148                 << setw(5) << m_nodes->getGlobalDegreesOfFreedom()[i] << " "      line 0: class 2
149                 << setw(5) << m_nodes->getGlobalNodesIndex()[i] << " "      line 1: class 3,4
150                 << setw(5) << m_nodes->getGlobalReducedDOFIndex()[i] << " "      line 2: class 5,6
151                 << setw(5) << m_nodes->getGlobalReducedNodesIndex()[i] << ": ";      */
152              for (int j = 0; j < numDim; j++)      if (fs.empty())
153                  cout << " " << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];          return false;
154              cout << endl;      vector<bool> hasclass(7, false);
155        vector<int> hasline(3, 0);
156        for (size_t i=0; i<fs.size(); ++i) {
157            switch (fs[i]) {
158                case Nodes:
159                    hasclass[0]=true;
160                    break;
161                case ReducedNodes:
162                    hasclass[1]=true;
163                    break;
164                case Points:
165                    hasline[0]=1;
166                    hasclass[2]=true;
167                    break;
168                case Elements:
169                    hasline[1]=1;
170                    break;
171                case ReducedElements:
172                    hasclass[4]=true;
173                    hasline[1]=1;
174                    break;
175                case FaceElements:
176                    hasline[2]=1;
177                    break;
178                case ReducedFaceElements:
179                    hasclass[6]=true;
180                    hasline[2]=1;
181                    break;
182                default:
183                    return false;
184          }          }
185      }      }
186        int numLines=hasline[0]+hasline[1]+hasline[2];
187    
188      // write all element types      // fail if we have more than one leaf group
189      const char *eNames[3] = { "Elements", "Face elements", "Points" };      // = there are at least two branches we can't interpolate between
190      const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };      if (numLines > 1)
191            return false;
192      for (size_t k=0; k<3; k++) {      else if (numLines==1) {
193          int mine = 0, overlap = 0;          if (hasline[0]==1)
194          for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {              resultcode=Points;
195              if (eFiles[k]->getOwnerVector()[i] == m_mpiInfo->rank)          else if (hasline[1]==1) {
196                  mine++;              if (hasclass[4])
197                    resultcode=ReducedElements;
198              else              else
199                  overlap++;                  resultcode=Elements;
200          }          } else { // hasline[2]==1
201          cout << "\t" << eNames[k] << ": " << eFiles[k]->getName()              if (hasclass[6])
202              << " " << eFiles[k]->getNumElements()                  resultcode=ReducedFaceElements;
203              << " (TypeId=" << eFiles[k]->getTypeId() << ") owner="              else
204              << mine << " overlap=" << overlap << endl;                  resultcode=FaceElements;
         int NN = eFiles[k]->getNumNodes();  
         if (full) {  
             cout << "\t     Id   Tag Owner Color:  Nodes" << endl;  
             for (int i = 0; i < eFiles[k]->getNumElements(); i++) {  
                 cout << "\t  " << setw(5) << eFiles[k]->getIdVector()[i] << " "  
                    << setw(5) << eFiles[k]->getTagVector()[i] << " "  
                    << setw(5) << eFiles[k]->getOwnerVector()[i] << " "  
                    << setw(5) << eFiles[k]->getColorVector()[i] << ": ";  
                 for (int j = 0; j < NN; j++)  
                     cout << " " << setw(5) << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];  
                 cout << endl;  
             }  
         }  
     }  
   
   
     // write tags  
     if (m_tagMap.size() > 0) {  
         cout << "\tTags:" << endl;  
         TagMap::const_iterator it;  
         for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {  
             cout << "\t  " << setw(5) << it->second << " "  
                 << it->first << endl;  
205          }          }
206        } else { // numLines==0
207            if (hasclass[1])
208                resultcode=ReducedNodes;
209            else
210                resultcode=Nodes;
211      }      }
212        return true;
213  }  }
214    
215  void RipleyDomain::dump(const string& fileName) const  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
216                                                  int fsType_target) const
217  {  {
218  #ifdef USE_NETCDF      if (fsType_target != Nodes &&
219      NcDim* ncdim = NULL;              fsType_target != ReducedNodes &&
220      int mpi_size = m_mpiInfo->size;              fsType_target != Elements &&
221      int mpi_rank = m_mpiInfo->rank;              fsType_target != ReducedElements &&
222      int numDim = getDim();              fsType_target != FaceElements &&
223                fsType_target != ReducedFaceElements &&
224  /* Incoming token indicates it's my turn to write */              fsType_target != Points) {
225  #ifdef ESYS_MPI          stringstream msg;
226      MPI_Status status;          msg << "probeInterpolationOnDomain(): Invalid functionspace type "
227      if (mpi_rank>0) {              << fsType_target << " for " << getDescription();
228          int dummy;          throw RipleyException(msg.str());
         MPI_Recv(&dummy, 0, MPI_INT, mpi_rank-1, 81800, m_mpiInfo->comm, &status);  
229      }      }
 #endif  
230    
231      char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),      switch (fsType_source) {
232                                                        mpi_size, mpi_rank);          case Nodes:
233                return true;
234      // NetCDF error handler          case ReducedNodes:
235      NcError err(NcError::verbose_nonfatal);              return (fsType_target != Nodes);
236      // Create the file.          case Elements:
237      NcFile dataFile(newFileName, NcFile::Replace);              return (fsType_target==Elements ||
238      const string msgPrefix("dump: netCDF operation failed - ");                      fsType_target==ReducedElements);
239            case ReducedElements:
240      // check if writing was successful              return (fsType_target==ReducedElements);
241      if (!dataFile.is_valid())          case FaceElements:
242          throw RipleyException(msgPrefix+"Open file for output");              return (fsType_target==FaceElements ||
243                        fsType_target==ReducedFaceElements);
244      const size_t numTags = m_tagMap.size();          case ReducedFaceElements:
245                return (fsType_target==ReducedFaceElements);
246      if (numTags > 0)          case Points:
247          if (! (ncdim = dataFile.add_dim("dim_Tags", numTags)) )              return (fsType_target==Points);
             throw RipleyException(msgPrefix+"add_dim(dim_Tags)");  
   
     // Attributes: MPI size, MPI rank, Name, order, reduced_order  
     if (!dataFile.add_att("mpi_size", mpi_size))  
         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);  
         }  
248    
249          // Tags_keys          default: {
250          NcVar *ncVar;              stringstream msg;
251          if (! (ncVar = dataFile.add_var("Tags_keys", ncInt, ncdim)) )              msg << "probeInterpolationOnDomain(): Invalid functionspace type "
252              throw RipleyException(msgPrefix+"add_var(Tags_keys)");                  << fsType_source << " for " << getDescription();
253          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++;  
254          }          }
255      }      }
   
     // 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 */  
 }  
   
 string RipleyDomain::getDescription() 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  
 {  
     return ReducedNodes;  
 }  
   
 int RipleyDomain::getFunctionCode() const  
 {  
     return Elements;  
 }  
   
 int RipleyDomain::getReducedFunctionCode() const  
 {  
     return ReducedElements;  
 }  
   
 int RipleyDomain::getFunctionOnBoundaryCode() const  
 {  
     return FaceElements;  
 }  
   
 int RipleyDomain::getReducedFunctionOnBoundaryCode() const  
 {  
     return ReducedFaceElements;  
256  }  }
257    
258  int RipleyDomain::getFunctionOnContactZeroCode() const  void RipleyDomain::interpolateOnDomain(escript::Data& target,
259  {                                         const escript::Data& in) const
     throw RipleyException("Ripley does not support contact elements");  
 }  
   
 int RipleyDomain::getReducedFunctionOnContactZeroCode() const  
260  {  {
261      throw RipleyException("Ripley does not support contact elements");      const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
262  }      const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
263        if (inDomain != *this)
264            throw RipleyException("Illegal domain of interpolant");
265        if (targetDomain != *this)
266            throw RipleyException("Illegal domain of interpolation target");
267    
268  int RipleyDomain::getFunctionOnContactOneCode() const      stringstream msg;
269  {      msg << "interpolateOnDomain() not implemented for function space "
270      throw RipleyException("Ripley does not support contact elements");          << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
271  }          << " -> "
272            << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
273    
274        const int inFS = in.getFunctionSpace().getTypeCode();
275        const int outFS = target.getFunctionSpace().getTypeCode();
276    
277        // simplest case: 1:1 copy
278        if (inFS==outFS) {
279            copyData(target, *const_cast<escript::Data*>(&in));
280        // not allowed: reduced->non-reduced
281        } else if (inFS==ReducedNodes && outFS==Nodes) {
282            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
283        } else if ((inFS==Elements && outFS==ReducedElements)
284                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
285            averageData(target, *const_cast<escript::Data*>(&in));
286        } else {
287            switch (inFS) {
288                case Nodes:
289                case ReducedNodes:
290                    switch (outFS) {
291                        case Nodes:
292                        case ReducedNodes: //FIXME: reduced
293                            copyData(target, *const_cast<escript::Data*>(&in));
294                            break;
295    
296                        case Elements:
297                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
298                            break;
299    
300                        case ReducedElements:
301                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
302                            break;
303    
304                        case FaceElements:
305                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
306                            break;
307    
308                        case ReducedFaceElements:
309                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
310                            break;
311    
312  int RipleyDomain::getReducedFunctionOnContactOneCode() const                      default:
313  {                          throw RipleyException(msg.str());
314      throw RipleyException("Ripley does not support contact elements");                  }
315                    break;
316                default:
317                    throw RipleyException(msg.str());
318            }
319        }
320  }  }
321    
322  int RipleyDomain::getSolutionCode() const  escript::Data RipleyDomain::getX() const
323  {  {
324      return DegreesOfFreedom;      return escript::continuousFunction(*this).getX();
325  }  }
326    
327  int RipleyDomain::getReducedSolutionCode() const  escript::Data RipleyDomain::getNormal() const
328  {  {
329      return ReducedDegreesOfFreedom;      return escript::functionOnBoundary(*this).getNormal();
330  }  }
331    
332  int RipleyDomain::getDiracDeltaFunctionsCode() const  escript::Data RipleyDomain::getSize() const
333  {  {
334      return Points;      return escript::function(*this).getSize();
335  }  }
336    
337  //  void RipleyDomain::setToX(escript::Data& arg) const
 // returns the spatial dimension of the mesh  
 //  
 int RipleyDomain::getDim() const  
338  {  {
339      return m_nodes->getNumDim();      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
340  }              *(arg.getFunctionSpace().getDomain()));
341        if (argDomain != *this)
342            throw RipleyException("setToX: Illegal domain of data point locations");
343        if (!arg.isExpanded())
344            throw RipleyException("setToX: Expanded Data object expected");
345    
346  //      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
347  // returns the number of data points summed across all MPI processes          assembleCoordinates(arg);
348  //      } else {
349  int RipleyDomain::getNumDataPointsGlobal() const          // interpolate the result
350  {          escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
351      return m_nodes->getGlobalNumNodes();          assembleCoordinates(contData);
352            interpolateOnDomain(arg, contData);
353        }
354  }  }
355    
356  //  bool RipleyDomain::isCellOriented(int fsType) const
 // 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  
357  {  {
358      int numDataPointsPerSample=0;      switch(fsType) {
     int numSamples=0;  
     switch (functionSpaceCode) {  
359          case Nodes:          case Nodes:
360              numDataPointsPerSample=1;              return false;
             numSamples=m_nodes->getNumNodes();  
             break;  
         case ReducedNodes:  
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumReducedNodes();  
             break;  
361          case Elements:          case Elements:
             numSamples=m_elements->getNumElements();  
             numDataPointsPerSample=m_elements->getNumLocalDim()+1;  
             break;  
         case ReducedElements:  
             numSamples=m_elements->getNumElements();  
             numDataPointsPerSample=(m_elements->getNumLocalDim()==0)?0:1;  
             break;  
362          case FaceElements:          case FaceElements:
             numDataPointsPerSample=m_faceElements->getNumLocalDim()+1;  
             numSamples=m_faceElements->getNumElements();  
             break;  
         case ReducedFaceElements:  
             numDataPointsPerSample=(m_faceElements->getNumLocalDim()==0)?0:1;  
             numSamples=m_faceElements->getNumElements();  
             break;  
363          case Points:          case Points:
364              numDataPointsPerSample=1;          case ReducedElements:
365              numSamples=m_points->getNumElements();          case ReducedFaceElements:
366              break;              return true;
         case DegreesOfFreedom:  
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumDegreesOfFreedom();  
             break;  
         case ReducedDegreesOfFreedom:  
             numDataPointsPerSample=1;  
             numSamples=m_nodes->getNumReducedDegreesOfFreedom();  
             break;  
367          default:          default:
368              stringstream temp;              break;
             temp << "Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();  
             throw RipleyException(temp.str());  
369      }      }
370      return pair<int,int>(numDataPointsPerSample,numSamples);      stringstream msg;
371  }      msg << "isCellOriented(): Illegal function space type " << fsType
372            << " on " << getDescription();
373  //      throw RipleyException(msg.str());
 // 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();  
 */  
 }  
   
 //  
 // 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  
 {  
     if (!d_contact.isEmpty() || !y_contact.isEmpty())  
         throw RipleyException("Ripley does not support contact elements");  
     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);  
     if (tpa==NULL)  
         throw RipleyException("Ripley only accepts its own transport problems");  
 /*  
     DataTypes::ShapeType shape;  
     source.expand();  
     escriptDataC _source=source.getDataC();  
     escriptDataC _M=M.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();  
     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();  
 */  
374  }  }
375    
376  //  bool RipleyDomain::canTag(int fsType) const
 // interpolates data between different function spaces  
 //  
 void RipleyDomain::interpolateOnDomain(escript::Data& target,  
                                       const escript::Data& in) const  
377  {  {
378      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()) {  
379          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;  
380          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;  
381          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;  
382          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;  
383          case ReducedFaceElements:          case ReducedFaceElements:
384              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;  
385          case Points:          case Points:
386              if (target.getFunctionSpace().getTypeCode()==Points) {              return false;
                 Ripley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
             } else {  
                 throw RipleyException("No interpolation with data on points possible.");  
             }  
             break;  
         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;  
         case ReducedDegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                     throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to mesh nodes.");  
                     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;  
387          default:          default:
388              stringstream temp;              break;
             temp << "Interpolation On Domain: Ripley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();  
             throw RipleyException(temp.str());  
389      }      }
390  */      stringstream msg;
391      checkPasoError();      msg << "canTag(): Illegal function space type " << fsType << " on "
392            << getDescription();
393        throw RipleyException(msg.str());
394  }  }
395    
396  //  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  
397  {  {
398      const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));      IndexVector* target=NULL;
399      if (argDomain != *this)      dim_t num=0;
         throw RipleyException("setToX: Illegal domain of data point locations");  
400    
401      // 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()) {  
402          case Nodes:          case Nodes:
403              throw RipleyException("Ripley does not support surface normal vectors for nodes");              num=getNumNodes();
404          case ReducedNodes:              target=&m_nodeTags;
405              throw RipleyException("Ripley does not support surface normal vectors for reduced nodes");              break;
406          case Elements:          case Elements:
             throw RipleyException("Ripley does not support surface normal vectors for elements");  
407          case ReducedElements:          case ReducedElements:
408              throw RipleyException("Ripley does not support surface normal vectors for elements with reduced integration order");              num=getNumElements();
409          case FaceElements:              target=&m_elementTags;
             Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);  
410              break;              break;
411            case FaceElements:
412          case ReducedFaceElements:          case ReducedFaceElements:
413              Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);              num=getNumFaceElements();
414                target=&m_faceTags;
415              break;              break;
416          case Points:          default: {
417              throw RipleyException("Ripley does not support surface normal vectors for point elements");              stringstream msg;
418          case DegreesOfFreedom:              msg << "setTags(): not implemented for "
419              throw RipleyException("Ripley does not support surface normal vectors for degrees of freedom.");                  << functionSpaceTypeAsString(fsType);
420          case ReducedDegreesOfFreedom:              throw RipleyException(msg.str());
421              throw RipleyException("Ripley does not support surface normal vectors for reduced degrees of freedom.");          }
422          default:      }
423              stringstream temp;      if (target->size() != num) {
424              temp << "Normal Vectors: Ripley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();          target->assign(num, -1);
             throw RipleyException(temp.str());  
425      }      }
 */  
 }  
   
 //  
 // 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");  
426    
427      throw RipleyException("Ripley does not allow interpolation across domains yet");      escript::Data& mask=*const_cast<escript::Data*>(&cMask);
428    #pragma omp parallel for
429        for (index_t i=0; i<num; i++) {
430            if (mask.getSampleDataRO(i)[0] > 0) {
431                (*target)[i]=newTag;
432            }
433        }
434        updateTagsInUse(fsType);
435  }  }
436    
437  //  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  
438  {  {
439  /*      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()) {  
440          case Nodes:          case Nodes:
441          case ReducedNodes:              if (m_nodeTags.size() > sampleNo)
442          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]);  
443              break;              break;
444          case Elements:          case Elements:
445          case ReducedElements:          case ReducedElements:
446              Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_arg,&integrals[0]);              if (m_elementTags.size() > sampleNo)
447                    return m_elementTags[sampleNo];
448              break;              break;
449          case FaceElements:          case FaceElements:
450          case ReducedFaceElements:          case ReducedFaceElements:
451              Ripley_Assemble_integrate(m_nodes,mesh->FaceElements,&_arg,&integrals[0]);              if (m_faceTags.size() > sampleNo)
452                    return m_faceTags[sampleNo];
453              break;              break;
454          case Points:          default: {
455              throw RipleyException("Integral of data on points is not supported.");              stringstream msg;
456          default:              msg << "getTagFromSampleNo(): not implemented for "
457              stringstream temp;                  << functionSpaceTypeAsString(fsType);
458              temp << "Integrals: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();              throw RipleyException(msg.str());
459              throw RipleyException(temp.str());          }
460      }      }
461      blocktimer_increment("integrate()", blocktimer_start);      return -1;
 */  
462  }  }
463    
464  //  int RipleyDomain::getNumberOfTagsInUse(int fsType) const
 // calculates the gradient of arg  
 //  
 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const  
465  {  {
466  /*      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()) {  
467          case Nodes:          case Nodes:
468              throw RipleyException("Gradient at nodes is not supported.");              return m_nodeTagsInUse.size();
         case ReducedNodes:  
             throw RipleyException("Gradient at reduced nodes is not supported.");  
469          case Elements:          case Elements:
             Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);  
             break;  
470          case ReducedElements:          case ReducedElements:
471              Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);              return m_elementTagsInUse.size();
             break;  
472          case FaceElements:          case FaceElements:
             Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);  
             break;  
473          case ReducedFaceElements:          case ReducedFaceElements:
474              Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);              return m_faceTagsInUse.size();
475              break;          default: {
476          case Points:              stringstream msg;
477              throw RipleyException("Gradient at points is not supported");              msg << "getNumberOfTagsInUse(): not implemented for "
478          case DegreesOfFreedom:                  << functionSpaceTypeAsString(fsType);
479              throw RipleyException("Gradient at degrees of freedom is not supported");              throw RipleyException(msg.str());
480          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());  
481      }      }
 */  
482  }  }
483    
484  //  const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
 // returns the size of elements  
 //  
 void RipleyDomain::setToSize(escript::Data& size) const  
485  {  {
486  /*      switch(fsType) {
     Ripley_Mesh* mesh=m_ripleyMesh.get();  
     escriptDataC tmp=size.getDataC();  
     switch(size.getFunctionSpace().getTypeCode()) {  
487          case Nodes:          case Nodes:
488              throw RipleyException("Size of nodes is not supported");              return &m_nodeTagsInUse[0];
         case ReducedNodes:  
             throw RipleyException("Size of reduced nodes is not supported");  
489          case Elements:          case Elements:
             Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);  
             break;  
490          case ReducedElements:          case ReducedElements:
491              Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);              return &m_elementTagsInUse[0];
             break;  
492          case FaceElements:          case FaceElements:
             Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);  
             break;  
493          case ReducedFaceElements:          case ReducedFaceElements:
494              Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);              return &m_faceTagsInUse[0];
495              break;          default: {
496          case Points:              stringstream msg;
497              throw RipleyException("Size of point elements is not supported");              msg << "borrowListOfTagsInUse(): not implemented for "
498          case DegreesOfFreedom:                  << functionSpaceTypeAsString(fsType);
499              throw RipleyException("Size of degrees of freedom is not supported");              throw RipleyException(msg.str());
500          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());  
501      }      }
 */  
502  }  }
503    
504  //  void RipleyDomain::Print_Mesh_Info(const bool full) const
 // sets the location of nodes  
 //  
 void RipleyDomain::setNewX(const escript::Data& new_x)  
505  {  {
506      const RipleyDomain& newDomain=dynamic_cast<const RipleyDomain&>(*(new_x.getFunctionSpace().getDomain()));      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
507      if (newDomain!=*this)          << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
508          throw RipleyException("Illegal domain of new node locations");      cout << "Number of dimensions: " << m_numDim << endl;
509    
510      if (new_x.getFunctionSpace() == continuousFunction(*this)) {      // write tags
511          m_nodes->setCoordinates(new_x);      if (m_tagMap.size() > 0) {
512      } else {          cout << "Tags:" << endl;
513          escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this) );          TagMap::const_iterator it;
514          m_nodes->setCoordinates(new_x_inter);          for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
515                cout << "  " << setw(5) << it->second << " "
516                    << it->first << endl;
517            }
518      }      }
519  }  }
520    
521  bool RipleyDomain::ownSample(int fs_code, index_t id) const  int RipleyDomain::getSystemMatrixTypeId(const int solver,
522            const int preconditioner, const int package, const bool symmetry) const
523  {  {
524  #ifdef ESYS_MPI      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
525      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;  
526  }  }
527    
528    int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
529            const int package, const bool symmetry) const
530    {
531        return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
532                package, symmetry, m_mpiInfo);
533    }
534    
 //  
 // creates a SystemMatrixAdapter stiffness matrix and initializes it with zeros  
 //  
535  escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,  escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
536          const escript::FunctionSpace& row_functionspace,          const escript::FunctionSpace& row_functionspace,
537          const int column_blocksize,          const int column_blocksize,
# Line 1136  escript::ASM_ptr RipleyDomain::newSystem Line 543  escript::ASM_ptr RipleyDomain::newSystem
543      // is the domain right?      // is the domain right?
544      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
545      if (row_domain!=*this)      if (row_domain!=*this)
546          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");
547      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
548      if (col_domain!=*this)      if (col_domain!=*this)
549          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");
550      // is the function space type right      // is the function space type right?
551      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      if (row_functionspace.getTypeCode()==ReducedNodes)
552          reduceRowOrder=true;          reduceRowOrder=true;
553      } else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=Nodes)
554          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
555      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      if (column_functionspace.getTypeCode()==ReducedNodes)
556          reduceColOrder=true;          reduceColOrder=true;
557      } else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=Nodes)
558          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
   
     // generate matrix:  
     Paso_SystemMatrixPattern* fsystemMatrixPattern=getPattern(reduceRowOrder, reduceColOrder);  
     checkPasoError();  
     Paso_SystemMatrix* fsystemMatrix = Paso_SystemMatrix_alloc(type,  
             fsystemMatrixPattern, row_blocksize, column_blocksize, FALSE);  
     checkPasoError();  
     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);  
     SystemMatrixAdapter* sma = new SystemMatrixAdapter(fsystemMatrix,  
             row_blocksize, row_functionspace, column_blocksize, column_functionspace);  
     return escript::ASM_ptr(sma);  
 }  
   
 //  
 // creates a TransportProblemAdapter  
 //  
 escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,  
         const int blocksize, const escript::FunctionSpace& functionspace,  
         const int type) const  
 {  
     int reduceOrder=0;  
     // is the domain right?  
     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");  
559    
560      // generate matrix      // generate matrix
561      Paso_SystemMatrixPattern* fsystemMatrixPattern = getPattern(reduceOrder, reduceOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
562      Paso_TransportProblem* transportProblem = Paso_TransportProblem_alloc(      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
563              useBackwardEuler, fsystemMatrixPattern, blocksize);              row_blocksize, column_blocksize, FALSE);
564      checkPasoError();      paso::checkPasoError();
565      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      Paso_SystemMatrixPattern_free(pattern);
566      escript::AbstractTransportProblem* atp=new TransportProblemAdapter(      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
567              transportProblem, useBackwardEuler, blocksize, functionspace);                  row_functionspace, column_blocksize, column_functionspace));
568      return escript::ATP_ptr(atp);      return sma;
569  }  }
570    
571  // returns true if data on the function space is considered as being cell  void RipleyDomain::setNewX(const escript::Data& arg)
 // centered  
 bool RipleyDomain::isCellOriented(int functionSpaceCode) const  
572  {  {
573      switch(functionSpaceCode) {      throw RipleyException("setNewX(): Operation not supported");
         case Nodes:  
         case DegreesOfFreedom:  
         case ReducedDegreesOfFreedom:  
             return false;  
         case Elements:  
         case FaceElements:  
         case Points:  
         case ReducedElements:  
         case ReducedFaceElements:  
             return true;  
         default:  
             stringstream temp;  
             temp << "Ripley does not know anything about function space type " << functionSpaceCode;  
             throw RipleyException(temp.str());  
     }  
     return false;  
574  }  }
575    
576  bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  //protected
577    void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
578  {  {
579     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]      const dim_t numComp = in.getDataPointSize();
580      class 1: DOF <-> Nodes      const dim_t dpp = in.getNumDataPointsPerSample();
581      class 2: ReducedDOF <-> ReducedNodes      out.requireWrite();
582      class 3: Points  #pragma omp parallel for
583      class 4: Elements      for (index_t i=0; i<in.getNumSamples(); i++) {
584      class 5: ReducedElements          const double* src = in.getSampleDataRO(i);
585      class 6: FaceElements          double* dest = out.getSampleDataRW(i);
586      class 7: ReducedFaceElements          for (index_t c=0; c<numComp; c++) {
587      class 8: ContactElementZero <-> ContactElementOne              double res=0.;
588      class 9: ReducedContactElementZero <-> ReducedContactElementOne              for (index_t q=0; q<dpp; q++)
589                    res+=src[c+q*numComp];
590     There is also a set of lines. Interpolation is possible down a line but not between lines.              *dest++ = res/dpp;
    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  
   
     For classes with multiple members (eg class 2) we have vars to record if  
     there is at least one instance. e.g. hasnodes is true if we have at least  
     one instance of Nodes.  
     */  
     if (fs.empty())  
         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;  
         }  
     }  
     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];  
   
     // fail if we have more than one leaf group  
     // = there are at least two branches we can't interpolate between  
     if (totlines>1)  
         return false;  
     else if (totlines==1) {  
         // we have points  
         if (hasline[0]==1) {  
             resultcode=Points;  
         } else if (hasline[1]==1) {  
             if (hasclass[5]==1) {  
                 resultcode=ReducedElements;  
             } else {  
                 resultcode=Elements;  
             }  
         } else if (hasline[2]==1) {  
             if (hasclass[7]==1) {  
                 resultcode=ReducedFaceElements;  
             } else {  
                 resultcode=FaceElements;  
             }  
         } else { // so we must be in line3  
             throw RipleyException("Programmer Error - choosing between contact elements - we should never get here");  
         }  
     } else { // totlines==0  
         if (hasclass[2]==1) {  
             // something from class 2  
             resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);  
         } else { // something from class 1  
             resultcode=(hasnodes ? Nodes : DegreesOfFreedom);  
591          }          }
592      }      }
     return true;  
593  }  }
594    
595  bool RipleyDomain::probeInterpolationOnDomain(int functionSpaceType_source,  //protected
596                                               int functionSpaceType_target) const  void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
597  {  {
598      if (functionSpaceType_target != Nodes &&      const dim_t numComp = in.getDataPointSize();
599              functionSpaceType_target != ReducedNodes &&      out.requireWrite();
600              functionSpaceType_target != ReducedDegreesOfFreedom &&  #pragma omp parallel for
601              functionSpaceType_target != DegreesOfFreedom &&      for (index_t i=0; i<in.getNumSamples(); i++) {
602              functionSpaceType_target != Elements &&          const double* src = in.getSampleDataRO(i);
603              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());  
604      }      }
605    }
606    
607      switch (functionSpaceType_source) {  //protected
608    void RipleyDomain::updateTagsInUse(int fsType) const
609    {
610        IndexVector* tagsInUse=NULL;
611        const IndexVector* tags=NULL;
612        switch(fsType) {
613          case Nodes:          case Nodes:
614          case DegreesOfFreedom:              tags=&m_nodeTags;
615              return true;              tagsInUse=&m_nodeTagsInUse;
616          case ReducedNodes:              break;
         case ReducedDegreesOfFreedom:  
             return (functionSpaceType_target != Nodes &&  
                     functionSpaceType_target != DegreesOfFreedom);  
617          case Elements:          case Elements:
             return (functionSpaceType_target==Elements ||  
                     functionSpaceType_target==ReducedElements);  
618          case ReducedElements:          case ReducedElements:
619              return (functionSpaceType_target==ReducedElements);              tags=&m_elementTags;
620                tagsInUse=&m_elementTagsInUse;
621                break;
622          case FaceElements:          case FaceElements:
             return (functionSpaceType_target==FaceElements ||  
                     functionSpaceType_target==ReducedFaceElements);  
623          case ReducedFaceElements:          case ReducedFaceElements:
624              return (functionSpaceType_target==ReducedFaceElements);              tags=&m_faceTags;
625          case Points:              tagsInUse=&m_faceTagsInUse;
626              return (functionSpaceType_target==Points);              break;
   
627          default:          default:
628              stringstream temp;              return;
             temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_source;  
             throw RipleyException(temp.str());  
629      }      }
     return false;  
 }  
630    
631  bool RipleyDomain::probeInterpolationACross(int functionSpaceType_source,      // gather global unique tag values from tags into tagsInUse
632          const AbstractDomain& targetDomain, int functionSpaceType_target) const      tagsInUse->clear();
633  {      index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
     return false;  
 }  
634    
635  bool RipleyDomain::operator==(const AbstractDomain& other) const      while (true) {
636  {          // find smallest value bigger than lastFoundValue
637      const RipleyDomain* temp=dynamic_cast<const RipleyDomain*>(&other);          minFoundValue = INDEX_T_MAX;
638      if (temp != NULL) {  #pragma omp parallel private(local_minFoundValue)
639          return (m_name==temp->m_name && m_nodes==temp->m_nodes &&          {
640                  m_elements==temp->m_elements &&              local_minFoundValue = minFoundValue;
641                  m_faceElements==temp->m_faceElements &&  #pragma omp for schedule(static) nowait
642                  m_points==temp->m_points);              for (size_t i = 0; i < tags->size(); i++) {
643                    const index_t v = (*tags)[i];
644                    if ((v > lastFoundValue) && (v < local_minFoundValue))
645                        local_minFoundValue = v;
646                }
647    #pragma omp critical
648                {
649                    if (local_minFoundValue < minFoundValue)
650                        minFoundValue = local_minFoundValue;
651                }
652            }
653    #ifdef ESYS_MPI
654            local_minFoundValue = minFoundValue;
655            MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
656    #endif
657    
658            // if we found a new value add it to the tagsInUse vector
659            if (minFoundValue < INDEX_T_MAX) {
660                tagsInUse->push_back(minFoundValue);
661                lastFoundValue = minFoundValue;
662            } else
663                break;
664      }      }
     return false;  
665  }  }
666    
667  bool RipleyDomain::operator!=(const AbstractDomain& other) const  //
668    // the methods that follow have to be implemented by the subclasses
669    //
670    
671    string RipleyDomain::getDescription() const
672  {  {
673      return !(operator==(other));      throw RipleyException("getDescription() not implemented");
674  }  }
675    
676  int RipleyDomain::getSystemMatrixTypeId(const int solver,  void RipleyDomain::write(const string& filename) const
         const int preconditioner, const int package, const bool symmetry) const  
677  {  {
678      int out=Paso_SystemMatrix_getSystemMatrixTypeId(      throw RipleyException("write() not implemented");
             SystemMatrixAdapter::mapOptionToPaso(solver),  
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
679  }  }
680    
681  int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,  void RipleyDomain::dump(const string& filename) const
         const int package, const bool symmetry) const  
682  {  {
683      int out=Paso_TransportProblem_getTypeId(      throw RipleyException("dump() not implemented");
             SystemMatrixAdapter::mapOptionToPaso(solver),  
             SystemMatrixAdapter::mapOptionToPaso(preconditioner),  
             SystemMatrixAdapter::mapOptionToPaso(package),  
             symmetry?1:0, m_mpiInfo);  
     checkPasoError();  
     return out;  
684  }  }
685    
686  escript::Data RipleyDomain::getX() const  const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
687  {  {
688      return continuousFunction(*this).getX();      throw RipleyException("borrowSampleReferenceIDs() not implemented");
689  }  }
690    
691  escript::Data RipleyDomain::getNormal() const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
692  {  {
693      return functionOnBoundary(*this).getNormal();      throw RipleyException("interpolateACross() not implemented");
694  }  }
695    
696  escript::Data RipleyDomain::getSize() const  bool RipleyDomain::probeInterpolationACross(int fsType_source,
697            const escript::AbstractDomain&, int fsType_target) const
698  {  {
699      return escript::function(*this).getSize();      throw RipleyException("probeInterpolationACross() not implemented");
700  }  }
701    
702  const int* RipleyDomain::borrowSampleReferenceIDs(int functionSpaceType) const  void RipleyDomain::setToNormal(escript::Data& normal) const
703  {  {
704      switch (functionSpaceType) {      throw RipleyException("setToNormal() 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());  
     }  
705  }  }
706    
707  int RipleyDomain::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  void RipleyDomain::setToSize(escript::Data& size) const
708  {  {
709      throw RipleyException("not implemented");      throw RipleyException("setToSize() not implemented");
710  }  }
711    
712    void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
 void RipleyDomain::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  
713  {  {
714      switch (functionSpaceType) {      throw RipleyException("setToGradient() 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());  
     }  
715  }  }
716    
717  void RipleyDomain::setTagMap(const string& name, int tag)  void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
718  {  {
719      m_tagMap[name] = tag;      throw RipleyException("setToIntegrals() not implemented");
720  }  }
721    
722  int RipleyDomain::getTag(const string& name) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
723  {  {
724      return m_tagMap.find(name)->second;      throw RipleyException("ownSample() not implemented");
725  }  }
726    
727  bool RipleyDomain::isValidTagName(const string& name) const  void RipleyDomain::addPDEToSystem(
728            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
729            const escript::Data& A, const escript::Data& B, const escript::Data& C,
730            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
731            const escript::Data& d, const escript::Data& y,
732            const escript::Data& d_contact, const escript::Data& y_contact,
733            const escript::Data& d_dirac,const escript::Data& y_dirac) const
734  {  {
735      return (m_tagMap.find(name)!=m_tagMap.end());      throw RipleyException("addPDEToSystem() not implemented");
736  }  }
737    
738  string RipleyDomain::showTagNames() const  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
739            const escript::Data& D, const escript::Data& d,
740            const escript::Data& d_dirac, const bool useHRZ) const
741  {  {
742      stringstream ret;      throw RipleyException("addPDEToLumpedSystem() 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();  
743  }  }
744    
745  int RipleyDomain::getNumberOfTagsInUse(int functionSpaceCode) const  void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
746            const escript::Data& Y, const escript::Data& y,
747            const escript::Data& y_contact, const escript::Data& y_dirac) const
748  {  {
749      dim_t numTags=0;      throw RipleyException("addPDEToRHS() 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;  
750  }  }
751    
752  const int* RipleyDomain::borrowListOfTagsInUse(int functionSpaceCode) const  void RipleyDomain::addPDEToTransportProblem(
753            escript::AbstractTransportProblem& tp,
754            escript::Data& source, const escript::Data& M,
755            const escript::Data& A, const escript::Data& B, const escript::Data& C,
756            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
757            const escript::Data& d, const escript::Data& y,
758            const escript::Data& d_contact, const escript::Data& y_contact,
759            const escript::Data& d_dirac, const escript::Data& y_dirac) const
760  {  {
761      switch (functionSpaceCode) {      throw RipleyException("addPDEToTransportProblem() 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());  
     }  
762  }  }
763    
764  bool RipleyDomain::canTag(int functionSpaceCode) const  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
765            const int blocksize, const escript::FunctionSpace& functionspace,
766            const int type) const
767  {  {
768      switch(functionSpaceCode) {      throw RipleyException("newTransportProblem() not implemented");
         case Nodes:  
         case Elements:  
         case ReducedElements:  
         case FaceElements:  
         case ReducedFaceElements:  
         case Points:  
             return true;  
     }  
     return false;  
769  }  }
770    
771  escript::AbstractDomain::StatusType RipleyDomain::getStatus() const  Paso_SystemMatrixPattern* RipleyDomain::getPattern(bool reducedRowOrder,
772                                                       bool reducedColOrder) const
773  {  {
774      return m_nodes->getStatus();      throw RipleyException("getPattern() not implemented");
775  }  }
776    
777  void RipleyDomain::prepare(bool optimize)  dim_t RipleyDomain::getNumDataPointsGlobal() const
778  {  {
779      resolveNodeIds();      throw RipleyException("getNumDataPointsGlobal() 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();  
780  }  }
781    
782  void RipleyDomain::optimizeElementOrdering()  IndexVector RipleyDomain::getNumNodesPerDim() const
783  {  {
784      m_elements->optimizeOrdering();      throw RipleyException("getNumNodesPerDim() not implemented");
     m_faceElements->optimizeOrdering();  
     m_points->optimizeOrdering();  
785  }  }
786    
787  void RipleyDomain::updateTagsInUse()  IndexVector RipleyDomain::getNumElementsPerDim() const
788  {  {
789      m_nodes->updateTagsInUse();      throw RipleyException("getNumElementsPerDim() not implemented");
     m_elements->updateTagsInUse();  
     m_faceElements->updateTagsInUse();  
     m_points->updateTagsInUse();  
790  }  }
791    
792  void RipleyDomain::createColoring(const IndexVector &node_localDOF_map)  IndexVector RipleyDomain::getNumFacesPerBoundary() const
793  {  {
794      m_elements->createColoring(node_localDOF_map);      throw RipleyException("getNumFacesPerBoundary() not implemented");
     m_faceElements->createColoring(node_localDOF_map);  
     m_points->createColoring(node_localDOF_map);  
795  }  }
796    
797  void RipleyDomain::distributeByRankOfDOF(const IndexVector &dofDistribution)  IndexVector RipleyDomain::getNodeDistribution() const
798  {  {
799      RankVector mpiRankOfDOF = m_nodes->getOwnerOfDOFs(dofDistribution);      throw RipleyException("getNodeDistribution() 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);  
800  }  }
801    
802  /**************************************************************  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
803     Check whether there is any node which has no vertex. In case  {
804     such a node exists, we don't use ParMetis since it requires      throw RipleyException("getFirstCoordAndSpacing() 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;  
805  }  }
 #endif  
806    
807  void RipleyDomain::optimizeDOFDistribution(RankVector &distribution)  dim_t RipleyDomain::getNumFaceElements() const
808  {  {
809      const Esys_MPI_rank myRank = m_mpiInfo->rank;      throw RipleyException("getNumFaceElements() not implemented");
     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];  
             }  
         }  
   
         // recvbuf will be the concatenation of each CPU's contribution  
         // to newDistribution  
         dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);  
   
 #ifdef ESYS_MPI  
         MPI_Allgather(&newDistribution[0], mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, m_mpiInfo->comm);  
 #else  
         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);  
   
         // now the overlap needs to be created by sending the partition  
         // around  
         m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);  
         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);  
   
     Esys_MPIInfo_noError(m_mpiInfo);  
     checkPasoError();  
   
     /* shift new labeling to create a global id */  
 #pragma omp parallel for  
     for (dim_t i = 0; i < myNumVertices; ++i)  
         newGlobalDOFid[i] += myFirstVertex;  
   
     /* distribute new labeling to other processors */  
     m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);  
 #if 0  
     for (i = 0; i < m_nodes->getNumNodes(); ++i)  
         printf("%d ", m_nodes->getGlobalDegreesOfFreedom()[i]);  
     printf("\n");  
 #endif  
810  }  }
811    
812  void RipleyDomain::resolveNodeIds()  dim_t RipleyDomain::getNumElements() const
813  {  {
814      // find the minimum and maximum node ID used by elements      throw RipleyException("getNumElements() not implemented");
     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);  
   
 #ifdef Ripley_TRACE  
     index_t globalMinId, globalMaxId;  
 #ifdef ESYS_MPI  
     index_t idRange[2], globalIdRange[2];  
     idRange[0] = -minId;  
     idRange[1] = maxId;  
     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;  
     }  
   
     // create a new node file  
     m_nodes = m_nodes->gatherGlobal(newLocalToGlobalNodeLabels);  
   
     // relabel nodes of the elements  
     relabelElementNodes(globalToNewLocalNodeLabels, minId);  
815  }  }
816    
817  void RipleyDomain::relabelElementNodes(const IndexVector &newNode, index_t offset)  dim_t RipleyDomain::getNumNodes() const
818  {  {
819      m_elements->relabelNodes(newNode, offset);      throw RipleyException("getNumNodes() not implemented");
     m_faceElements->relabelNodes(newNode, offset);  
     m_points->relabelNodes(newNode, offset);  
820  }  }
821    
822  void RipleyDomain::markNodes(IndexVector &mask, index_t offset)  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
823  {  {
824      m_elements->markNodes(mask, offset);      throw RipleyException("assembleCoordinates() not implemented");
     m_faceElements->markNodes(mask, offset);  
     m_points->markNodes(mask, offset);  
825  }  }
826    
827  void RipleyDomain::createMappings(const IndexVector &dofDistribution, const IndexVector &nodeDistribution)  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
828  {  {
829      IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);      throw RipleyException("interpolateNodesOnElements() not implemented");
     markNodes(maskReducedNodes, 0);  
   
     IndexVector indexReducedNodes = packMask(maskReducedNodes);  
     m_nodes->createNodeFileMappings(indexReducedNodes, dofDistribution,  
                                     nodeDistribution);  
830  }  }
831    
832  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
833  {  {
834      Paso_Connector *colConnector, *rowConnector;      throw RipleyException("interpolateNodesOnFaces() not implemented");
     NodeMapping_ptr colMap, rowMap;  
     Paso_Distribution *colDistribution = NULL, *rowDistribution = NULL;  
   
     if (reduce_col_order) {  
         colMap = m_nodes->getReducedDegreesOfFreedomMapping();  
         colDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();  
         colConnector = m_nodes->getReducedDegreesOfFreedomConnector();  
     } else {  
         colMap = m_nodes->getDegreesOfFreedomMapping();  
         colDistribution = m_nodes->getDegreesOfFreedomDistribution();  
         colConnector = m_nodes->getDegreesOfFreedomConnector();  
     }  
   
     if (reduce_row_order) {  
         rowMap = m_nodes->getReducedDegreesOfFreedomMapping();  
         rowDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();  
         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;  
     */  
835  }  }
836    
837    
838  } // end of namespace ripley  } // end of namespace ripley
839    

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

  ViewVC Help
Powered by ViewVC 1.1.26