/[escript]/branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapter.cpp

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

revision 6078 by caltinay, Thu Mar 10 07:04:41 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 39  using esys_trilinos::const_TrilinosGraph Line 39  using esys_trilinos::const_TrilinosGraph
39    
40  using namespace std;  using namespace std;
41  namespace bp = boost::python;  namespace bp = boost::python;
42    using escript::ValueError;
43    
44  namespace dudley {  namespace dudley {
45    
46  // define the static constants  // define the static constants
47  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
48    // dudley only supports single approximation order
49  const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;
 const int MeshAdapter::ReducedDegreesOfFreedom=DUDLEY_REDUCED_DEGREES_OF_FREEDOM;  
50  const int MeshAdapter::Nodes=DUDLEY_NODES;  const int MeshAdapter::Nodes=DUDLEY_NODES;
 const int MeshAdapter::ReducedNodes=DUDLEY_REDUCED_NODES;  
51  const int MeshAdapter::Elements=DUDLEY_ELEMENTS;  const int MeshAdapter::Elements=DUDLEY_ELEMENTS;
52  const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;  const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;
53  const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;
54  const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;  const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;
55  const int MeshAdapter::Points=DUDLEY_POINTS;  const int MeshAdapter::Points=DUDLEY_POINTS;
56    
57  MeshAdapter::MeshAdapter(Dudley_Mesh* dudleyMesh)  MeshAdapter::MeshAdapter(Mesh* dudleyMesh) :
58        m_dudleyMesh(dudleyMesh)
59  {  {
60      setFunctionSpaceTypeNames();      setFunctionSpaceTypeNames();
     // need to use a null_deleter as Dudley_Mesh_free deletes the pointer  
     // for us.  
     m_dudleyMesh.reset(dudleyMesh,null_deleter());  
61  }  }
62    
 //  
63  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
64  MeshAdapter::MeshAdapter(const MeshAdapter& in) :  MeshAdapter::MeshAdapter(const MeshAdapter& in) :
65      m_dudleyMesh(in.m_dudleyMesh)      m_dudleyMesh(in.m_dudleyMesh)
# Line 72  MeshAdapter::MeshAdapter(const MeshAdapt Line 69  MeshAdapter::MeshAdapter(const MeshAdapt
69    
70  MeshAdapter::~MeshAdapter()  MeshAdapter::~MeshAdapter()
71  {  {
     // I hope the case for the pointer being zero has been taken care of  
     if (m_dudleyMesh.unique()) {  
         Dudley_Mesh_free(m_dudleyMesh.get());  
     }  
72  }  }
73    
74  escript::JMPI MeshAdapter::getMPI() const  escript::JMPI MeshAdapter::getMPI() const
75  {  {
76      return m_dudleyMesh.get()->MPIInfo;      return m_dudleyMesh->MPIInfo;
77  }  }
78    
79  int MeshAdapter::getMPISize() const  int MeshAdapter::getMPISize() const
80  {  {
81      return m_dudleyMesh.get()->MPIInfo->size;      return getMPI()->size;
82  }  }
83    
84  int MeshAdapter::getMPIRank() const  int MeshAdapter::getMPIRank() const
85  {  {
86      return m_dudleyMesh.get()->MPIInfo->rank;      return getMPI()->rank;
87  }  }
88    
89  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
90  {  {
91  #ifdef ESYS_MPI  #ifdef ESYS_MPI
92      MPI_Barrier(m_dudleyMesh.get()->MPIInfo->comm);      MPI_Barrier(getMPIComm());
93  #endif  #endif
94  }  }
95    
96  bool MeshAdapter::onMasterProcessor() const  bool MeshAdapter::onMasterProcessor() const
97  {  {
98      return m_dudleyMesh.get()->MPIInfo->rank == 0;      return getMPIRank() == 0;
99  }  }
100    
101  MPI_Comm MeshAdapter::getMPIComm() const  MPI_Comm MeshAdapter::getMPIComm() const
102  {  {
103      return m_dudleyMesh->MPIInfo->comm;      return getMPI()->comm;
104  }  }
105    
106  Dudley_Mesh* MeshAdapter::getDudley_Mesh() const  Mesh* MeshAdapter::getMesh() const
107  {  {
108      return m_dudleyMesh.get();      return m_dudleyMesh.get();
109  }  }
110    
111  void MeshAdapter::write(const string& fileName) const  void MeshAdapter::write(const string& fileName) const
112  {  {
113      char *fName = (fileName.size()+1>0) ? new char[fileName.size()+1] : (char*)NULL;      m_dudleyMesh->write(fileName);
     strcpy(fName, fileName.c_str());  
     Dudley_Mesh_write(m_dudleyMesh.get(),fName);  
     delete[] fName;  
114  }  }
115    
116  void MeshAdapter::Print_Mesh_Info(const bool full) const  void MeshAdapter::Print_Mesh_Info(bool full) const
117  {  {
118      Dudley_PrintMesh_Info(m_dudleyMesh.get(), full);      m_dudleyMesh->printInfo(full);
119  }  }
120    
121  void MeshAdapter::dump(const string& fileName) const  void MeshAdapter::dump(const string& fileName) const
122  {  {
123  #ifdef USE_NETCDF  #ifdef USE_NETCDF
124     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
125     NcVar *ids;     NcVar* ids;
126     int *int_ptr;     index_t* index_ptr;
127     Dudley_Mesh *mesh = m_dudleyMesh.get();  #ifdef ESYS_INDEXTYPE_LONG
128     Dudley_TagMap* tag_map;      NcType ncIdxType = ncLong;
129    #else
130        NcType ncIdxType = ncInt;
131    #endif
132       Mesh* mesh = m_dudleyMesh.get();
133     int num_Tags = 0;     int num_Tags = 0;
134     int mpi_size                         = mesh->MPIInfo->size;     int mpi_size                  = getMPISize();
135     int mpi_rank                         = mesh->MPIInfo->rank;     int mpi_rank                  = getMPIRank();
136     int numDim                           = mesh->Nodes->numDim;     int numDim                    = mesh->Nodes->numDim;
137     int numNodes                         = mesh->Nodes->numNodes;     dim_t numNodes                = mesh->Nodes->getNumNodes();
138     int num_Elements                     = mesh->Elements->numElements;     dim_t num_Elements            = mesh->Elements->numElements;
139     int num_FaceElements                 = mesh->FaceElements->numElements;     dim_t num_FaceElements        = mesh->FaceElements->numElements;
140     int num_Points                       = mesh->Points->numElements;     dim_t num_Points              = mesh->Points->numElements;
141     int num_Elements_numNodes            = mesh->Elements->numNodes;     int num_Elements_numNodes     = mesh->Elements->numNodes;
142     int num_FaceElements_numNodes        = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
143  #ifdef ESYS_MPI  #ifdef ESYS_MPI
144     MPI_Status status;     MPI_Status status;
145  #endif  #endif
146    
147  /* Incoming token indicates it's my turn to write */      // Incoming token indicates it's my turn to write
148  #ifdef ESYS_MPI  #ifdef ESYS_MPI
149     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);      if (mpi_rank > 0)
150            MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, getMPIComm(), &status);
151  #endif  #endif
152    
153     string newFileName(mesh->MPIInfo->appendRankToFileName(fileName));      const string newFileName(mesh->MPIInfo->appendRankToFileName(fileName));
154    
155        // Figure out how much storage is required for tags
156        num_Tags = mesh->tagMap.size();
157    
158     /* Figure out how much storage is required for tags */      // NetCDF error handler
159     tag_map = mesh->TagMap;      NcError err(NcError::verbose_nonfatal);
160     num_Tags = 0;      // Create the file
161     while (tag_map) {      NcFile dataFile(newFileName.c_str(), NcFile::Replace);
162        num_Tags++;      string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
163        tag_map=tag_map->next;      // check if writing was successful
164     }      if (!dataFile.is_valid())
165            throw DudleyException(msgPrefix + "Open file for output");
166     // NetCDF error handler  
167     NcError err(NcError::verbose_nonfatal);      // Define dimensions (num_Elements and dim_Elements are identical,
168     // Create the file.      // dim_Elements only appears if > 0)
169     NcFile dataFile(newFileName.c_str(), NcFile::Replace);      if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
170     string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");          throw DudleyException(msgPrefix+"add_dim(numNodes)");
171     // check if writing was successful      if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
172     if (!dataFile.is_valid())          throw DudleyException(msgPrefix+"add_dim(numDim)");
173        throw DudleyException(msgPrefix+"Open file for output");      if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
174            throw DudleyException(msgPrefix+"add_dim(mpi_size)");
175     // Define dimensions (num_Elements and dim_Elements are identical,      if (num_Elements > 0)
176     // dim_Elements only appears if > 0)          if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
177     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )              throw DudleyException(msgPrefix+"add_dim(dim_Elements)");
178        throw DudleyException(msgPrefix+"add_dim(numNodes)");      if (num_FaceElements > 0)
179     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )          if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
       throw DudleyException(msgPrefix+"add_dim(numDim)");  
    if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )  
       throw DudleyException(msgPrefix+"add_dim(mpi_size)");  
    if (num_Elements>0)  
       if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )  
          throw DudleyException(msgPrefix+"add_dim(dim_Elements)");  
    if (num_FaceElements>0)  
       if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )  
180           throw DudleyException(msgPrefix+"add_dim(dim_FaceElements)");           throw DudleyException(msgPrefix+"add_dim(dim_FaceElements)");
181     if (num_Points>0)      if (num_Points > 0)
182        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )          if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
183           throw DudleyException(msgPrefix+"add_dim(dim_Points)");              throw DudleyException(msgPrefix+"add_dim(dim_Points)");
184     if (num_Elements>0)      if (num_Elements > 0)
185        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )          if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
186           throw DudleyException(msgPrefix+"add_dim(dim_Elements_Nodes)");              throw DudleyException(msgPrefix+"add_dim(dim_Elements_Nodes)");
187     if (num_FaceElements>0)      if (num_FaceElements > 0)
188        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )          if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
189           throw DudleyException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");              throw DudleyException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
190     if (num_Tags>0)      if (num_Tags > 0)
191        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )          if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
192           throw DudleyException(msgPrefix+"add_dim(dim_Tags)");              throw DudleyException(msgPrefix+"add_dim(dim_Tags)");
193    
194     // Attributes: MPI size, MPI rank, Name, order, reduced_order      // Attributes: MPI size, MPI rank, Name, order, reduced_order
195     if (!dataFile.add_att("mpi_size", mpi_size) )      if (!dataFile.add_att("index_size", (int)sizeof(index_t)))
196        throw DudleyException(msgPrefix+"add_att(mpi_size)");          throw DudleyException(msgPrefix+"add_att(index_size)");
197     if (!dataFile.add_att("mpi_rank", mpi_rank) )      if (!dataFile.add_att("mpi_size", mpi_size) )
198        throw DudleyException(msgPrefix+"add_att(mpi_rank)");          throw DudleyException(msgPrefix+"add_att(mpi_size)");
199     if (!dataFile.add_att("Name",mesh->Name) )      if (!dataFile.add_att("mpi_rank", mpi_rank) )
200        throw DudleyException(msgPrefix+"add_att(Name)");          throw DudleyException(msgPrefix+"add_att(mpi_rank)");
201     if (!dataFile.add_att("numDim",numDim) )      if (!dataFile.add_att("Name",mesh->m_name.c_str()) )
202        throw DudleyException(msgPrefix+"add_att(order)");          throw DudleyException(msgPrefix+"add_att(Name)");
203     if (!dataFile.add_att("order",mesh->integrationOrder) )      if (!dataFile.add_att("numDim",numDim) )
204        throw DudleyException(msgPrefix+"add_att(order)");          throw DudleyException(msgPrefix+"add_att(order)");
205     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )      if (!dataFile.add_att("order",mesh->integrationOrder) )
206        throw DudleyException(msgPrefix+"add_att(reduced_order)");          throw DudleyException(msgPrefix+"add_att(order)");
207     if (!dataFile.add_att("numNodes",numNodes) )      if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
208        throw DudleyException(msgPrefix+"add_att(numNodes)");          throw DudleyException(msgPrefix+"add_att(reduced_order)");
209     if (!dataFile.add_att("num_Elements",num_Elements) )      if (!dataFile.add_att("numNodes",numNodes) )
210        throw DudleyException(msgPrefix+"add_att(num_Elements)");          throw DudleyException(msgPrefix+"add_att(numNodes)");
211     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )      if (!dataFile.add_att("num_Elements",num_Elements) )
212        throw DudleyException(msgPrefix+"add_att(num_FaceElements)");          throw DudleyException(msgPrefix+"add_att(num_Elements)");
213     if (!dataFile.add_att("num_Points",num_Points) )      if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
214        throw DudleyException(msgPrefix+"add_att(num_Points)");          throw DudleyException(msgPrefix+"add_att(num_FaceElements)");
215     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )      if (!dataFile.add_att("num_Points",num_Points) )
216        throw DudleyException(msgPrefix+"add_att(num_Elements_numNodes)");          throw DudleyException(msgPrefix+"add_att(num_Points)");
217     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )      if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
218        throw DudleyException(msgPrefix+"add_att(num_FaceElements_numNodes)");          throw DudleyException(msgPrefix+"add_att(num_Elements_numNodes)");
219     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->etype) )      if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
220        throw DudleyException(msgPrefix+"add_att(Elements_TypeId)");          throw DudleyException(msgPrefix+"add_att(num_FaceElements_numNodes)");
221     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->etype) )      if (!dataFile.add_att("Elements_TypeId", mesh->Elements->etype) )
222        throw DudleyException(msgPrefix+"add_att(FaceElements_TypeId)");          throw DudleyException(msgPrefix+"add_att(Elements_TypeId)");
223     if (!dataFile.add_att("Points_TypeId", mesh->Points->etype) )      if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->etype) )
224        throw DudleyException(msgPrefix+"add_att(Points_TypeId)");          throw DudleyException(msgPrefix+"add_att(FaceElements_TypeId)");
225     if (!dataFile.add_att("num_Tags", num_Tags) )      if (!dataFile.add_att("Points_TypeId", mesh->Points->etype) )
226        throw DudleyException(msgPrefix+"add_att(num_Tags)");          throw DudleyException(msgPrefix+"add_att(Points_TypeId)");
227        if (!dataFile.add_att("num_Tags", num_Tags) )
228     // // // // // Nodes // // // // //          throw DudleyException(msgPrefix+"add_att(num_Tags)");
229    
230     // Nodes nodeDistribution      // // // // // Nodes // // // // //
231     if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
232        throw DudleyException(msgPrefix+"add_var(Nodes_NodeDistribution)");      // Nodes nodeDistribution
233     int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];      if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncIdxType, ncdims[2])) )
234     if (! (ids->put(int_ptr, mpi_size+1)) )          throw DudleyException(msgPrefix+"add_var(Nodes_NodeDistribution)");
235        throw DudleyException(msgPrefix+"put(Nodes_NodeDistribution)");      index_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
236        if (! (ids->put(index_ptr, mpi_size+1)) )
237     // Nodes degreesOfFreedomDistribution          throw DudleyException(msgPrefix+"put(Nodes_NodeDistribution)");
238     if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
239        throw DudleyException(msgPrefix+"add_var(Nodes_DofDistribution)");      // Nodes degreesOfFreedomDistribution
240     int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];      if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncIdxType, ncdims[2])) )
241     if (! (ids->put(int_ptr, mpi_size+1)) )          throw DudleyException(msgPrefix+"add_var(Nodes_DofDistribution)");
242        throw DudleyException(msgPrefix+"put(Nodes_DofDistribution)");      index_ptr = &mesh->Nodes->dofDistribution->first_component[0];
243        if (! (ids->put(index_ptr, mpi_size+1)) )
244     // Only write nodes if non-empty because NetCDF doesn't like empty arrays          throw DudleyException(msgPrefix+"put(Nodes_DofDistribution)");
245     // (it treats them as NC_UNLIMITED)  
246     if (numNodes>0) {      // Only write nodes if non-empty because NetCDF doesn't like empty arrays
247        // (it treats them as NC_UNLIMITED)
248        // Nodes Id      if (numNodes > 0) {
249        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )          // Nodes Id
250           throw DudleyException(msgPrefix+"add_var(Nodes_Id)");          if (! ( ids = dataFile.add_var("Nodes_Id", ncIdxType, ncdims[0])) )
251        int_ptr = &mesh->Nodes->Id[0];              throw DudleyException(msgPrefix+"add_var(Nodes_Id)");
252        if (! (ids->put(int_ptr, numNodes)) )          if (! (ids->put(mesh->Nodes->Id, numNodes)) )
253           throw DudleyException(msgPrefix+"put(Nodes_Id)");              throw DudleyException(msgPrefix+"put(Nodes_Id)");
254    
255        // Nodes Tag          // Nodes Tag
256        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )          if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
257           throw DudleyException(msgPrefix+"add_var(Nodes_Tag)");              throw DudleyException(msgPrefix+"add_var(Nodes_Tag)");
258        int_ptr = &mesh->Nodes->Tag[0];          if (! (ids->put(mesh->Nodes->Tag, numNodes)) )
259        if (! (ids->put(int_ptr, numNodes)) )              throw DudleyException(msgPrefix+"put(Nodes_Tag)");
260           throw DudleyException(msgPrefix+"put(Nodes_Tag)");  
261            // Nodes gDOF
262        // Nodes gDOF          if (! ( ids = dataFile.add_var("Nodes_gDOF", ncIdxType, ncdims[0])) )
263        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )              throw DudleyException(msgPrefix+"add_var(Nodes_gDOF)");
264           throw DudleyException(msgPrefix+"add_var(Nodes_gDOF)");          if (! (ids->put(mesh->Nodes->globalDegreesOfFreedom, numNodes)) )
265        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];              throw DudleyException(msgPrefix+"put(Nodes_gDOF)");
266        if (! (ids->put(int_ptr, numNodes)) )  
267           throw DudleyException(msgPrefix+"put(Nodes_gDOF)");          // Nodes global node index
268            if (! ( ids = dataFile.add_var("Nodes_gNI", ncIdxType, ncdims[0])) )
269        // Nodes global node index              throw DudleyException(msgPrefix+"add_var(Nodes_gNI)");
270        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )          if (! (ids->put(mesh->Nodes->globalNodesIndex, numNodes)) )
271           throw DudleyException(msgPrefix+"add_var(Nodes_gNI)");              throw DudleyException(msgPrefix+"put(Nodes_gNI)");
272        int_ptr = &mesh->Nodes->globalNodesIndex[0];  
273        if (! (ids->put(int_ptr, numNodes)) )          // Nodes Coordinates
274           throw DudleyException(msgPrefix+"put(Nodes_gNI)");          if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
275                throw DudleyException(msgPrefix+"add_var(Nodes_Coordinates)");
276        // Nodes grDof          if (! (ids->put(mesh->Nodes->Coordinates, numNodes, numDim)) )
277        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )              throw DudleyException(msgPrefix+"put(Nodes_Coordinates)");
278           throw DudleyException(msgPrefix+"add_var(Nodes_grDfI)");      }
279        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];  
280        if (! (ids->put(int_ptr, numNodes)) )      // // // // // Elements // // // // //
281           throw DudleyException(msgPrefix+"put(Nodes_grDfI)");      if (num_Elements > 0) {
282            // Elements_Id
283        // Nodes grNI          if (! ( ids = dataFile.add_var("Elements_Id", ncIdxType, ncdims[3])) )
284        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )              throw DudleyException(msgPrefix+"add_var(Elements_Id)");
285           throw DudleyException(msgPrefix+"add_var(Nodes_grNI)");          if (! (ids->put(mesh->Elements->Id, num_Elements)) )
286        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];              throw DudleyException(msgPrefix+"put(Elements_Id)");
287        if (! (ids->put(int_ptr, numNodes)) )  
288           throw DudleyException(msgPrefix+"put(Nodes_grNI)");          // Elements_Tag
289            if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
290        // Nodes Coordinates              throw DudleyException(msgPrefix+"add_var(Elements_Tag)");
291        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )          if (! (ids->put(mesh->Elements->Tag, num_Elements)) )
292           throw DudleyException(msgPrefix+"add_var(Nodes_Coordinates)");              throw DudleyException(msgPrefix+"put(Elements_Tag)");
293        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )  
294           throw DudleyException(msgPrefix+"put(Nodes_Coordinates)");          // Elements_Owner
295            if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
296     }              throw DudleyException(msgPrefix+"add_var(Elements_Owner)");
297            if (! (ids->put(mesh->Elements->Owner, num_Elements)) )
298     // // // // // Elements // // // // //              throw DudleyException(msgPrefix+"put(Elements_Owner)");
299    
300     if (num_Elements>0) {          // Elements_Color
301            if (! ( ids = dataFile.add_var("Elements_Color", ncIdxType, ncdims[3])) )
302        // Elements_Id              throw DudleyException(msgPrefix+"add_var(Elements_Color)");
303        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )          if (! (ids->put(mesh->Elements->Color, num_Elements)) )
304           throw DudleyException(msgPrefix+"add_var(Elements_Id)");              throw DudleyException(msgPrefix+"put(Elements_Color)");
305        int_ptr = &mesh->Elements->Id[0];  
306        if (! (ids->put(int_ptr, num_Elements)) )          // Elements_Nodes
307           throw DudleyException(msgPrefix+"put(Elements_Id)");          if (! ( ids = dataFile.add_var("Elements_Nodes", ncIdxType, ncdims[3], ncdims[7]) ) )
308                throw DudleyException(msgPrefix+"add_var(Elements_Nodes)");
309        // Elements_Tag          if (! (ids->put(mesh->Elements->Nodes, num_Elements, num_Elements_numNodes)) )
310        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )              throw DudleyException(msgPrefix+"put(Elements_Nodes)");
311           throw DudleyException(msgPrefix+"add_var(Elements_Tag)");      }
312        int_ptr = &mesh->Elements->Tag[0];  
313        if (! (ids->put(int_ptr, num_Elements)) )      // // // // // Face_Elements // // // // //
314           throw DudleyException(msgPrefix+"put(Elements_Tag)");      if (num_FaceElements > 0) {
315            // FaceElements_Id
316        // Elements_Owner          if (!(ids = dataFile.add_var("FaceElements_Id", ncIdxType, ncdims[4])))
317        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )              throw DudleyException(msgPrefix+"add_var(FaceElements_Id)");
318           throw DudleyException(msgPrefix+"add_var(Elements_Owner)");          if (!(ids->put(mesh->FaceElements->Id, num_FaceElements)))
319        int_ptr = &mesh->Elements->Owner[0];              throw DudleyException(msgPrefix+"put(FaceElements_Id)");
320        if (! (ids->put(int_ptr, num_Elements)) )  
321           throw DudleyException(msgPrefix+"put(Elements_Owner)");          // FaceElements_Tag
322            if (!(ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])))
323        // Elements_Color              throw DudleyException(msgPrefix+"add_var(FaceElements_Tag)");
324        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )          if (!(ids->put(mesh->FaceElements->Tag, num_FaceElements)))
325           throw DudleyException(msgPrefix+"add_var(Elements_Color)");              throw DudleyException(msgPrefix+"put(FaceElements_Tag)");
326        int_ptr = &mesh->Elements->Color[0];  
327        if (! (ids->put(int_ptr, num_Elements)) )          // FaceElements_Owner
328           throw DudleyException(msgPrefix+"put(Elements_Color)");          if (!(ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])))
329                throw DudleyException(msgPrefix+"add_var(FaceElements_Owner)");
330        // Elements_Nodes          if (!(ids->put(mesh->FaceElements->Owner, num_FaceElements)))
331        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )              throw DudleyException(msgPrefix+"put(FaceElements_Owner)");
332           throw DudleyException(msgPrefix+"add_var(Elements_Nodes)");  
333        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )          // FaceElements_Color
334           throw DudleyException(msgPrefix+"put(Elements_Nodes)");          if (!(ids = dataFile.add_var("FaceElements_Color", ncIdxType, ncdims[4])))
335                throw DudleyException(msgPrefix+"add_var(FaceElements_Color)");
336     }          if (!(ids->put(mesh->FaceElements->Color, num_FaceElements)))
337                throw DudleyException(msgPrefix+"put(FaceElements_Color)");
338     // // // // // Face_Elements // // // // //  
339            // FaceElements_Nodes
340     if (num_FaceElements>0) {          if (!(ids = dataFile.add_var("FaceElements_Nodes", ncIdxType, ncdims[4], ncdims[8])))
341                throw DudleyException(msgPrefix+"add_var(FaceElements_Nodes)");
342        // FaceElements_Id          if (!(ids->put(mesh->FaceElements->Nodes, num_FaceElements, num_FaceElements_numNodes)))
343        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )              throw DudleyException(msgPrefix+"put(FaceElements_Nodes)");
344           throw DudleyException(msgPrefix+"add_var(FaceElements_Id)");      }
345        int_ptr = &mesh->FaceElements->Id[0];  
346        if (! (ids->put(int_ptr, num_FaceElements)) )      // // // // // Points // // // // //
347           throw DudleyException(msgPrefix+"put(FaceElements_Id)");      if (num_Points > 0) {
348            // Points_Id
349        // FaceElements_Tag          if (!(ids = dataFile.add_var("Points_Id", ncIdxType, ncdims[6])))
350        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )              throw DudleyException(msgPrefix+"add_var(Points_Id)");
351           throw DudleyException(msgPrefix+"add_var(FaceElements_Tag)");          if (!(ids->put(mesh->Points->Id, num_Points)))
352        int_ptr = &mesh->FaceElements->Tag[0];              throw DudleyException(msgPrefix+"put(Points_Id)");
353        if (! (ids->put(int_ptr, num_FaceElements)) )  
354           throw DudleyException(msgPrefix+"put(FaceElements_Tag)");          // Points_Tag
355            if (!(ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])))
356        // FaceElements_Owner              throw DudleyException(msgPrefix+"add_var(Points_Tag)");
357        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )          if (!(ids->put(mesh->Points->Tag, num_Points)))
358           throw DudleyException(msgPrefix+"add_var(FaceElements_Owner)");              throw DudleyException(msgPrefix+"put(Points_Tag)");
359        int_ptr = &mesh->FaceElements->Owner[0];  
360        if (! (ids->put(int_ptr, num_FaceElements)) )          // Points_Owner
361           throw DudleyException(msgPrefix+"put(FaceElements_Owner)");          if (!(ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])))
362                throw DudleyException(msgPrefix+"add_var(Points_Owner)");
363        // FaceElements_Color          if (!(ids->put(mesh->Points->Owner, num_Points)))
364        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )              throw DudleyException(msgPrefix+"put(Points_Owner)");
365           throw DudleyException(msgPrefix+"add_var(FaceElements_Color)");  
366        int_ptr = &mesh->FaceElements->Color[0];          // Points_Color
367        if (! (ids->put(int_ptr, num_FaceElements)) )          if (!(ids = dataFile.add_var("Points_Color", ncIdxType, ncdims[6])))
368           throw DudleyException(msgPrefix+"put(FaceElements_Color)");              throw DudleyException(msgPrefix+"add_var(Points_Color)");
369            if (!(ids->put(mesh->Points->Color, num_Points)))
370        // FaceElements_Nodes              throw DudleyException(msgPrefix+"put(Points_Color)");
371        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )  
372           throw DudleyException(msgPrefix+"add_var(FaceElements_Nodes)");          // Points_Nodes
373        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )          if (!(ids = dataFile.add_var("Points_Nodes", ncIdxType, ncdims[6])))
374           throw DudleyException(msgPrefix+"put(FaceElements_Nodes)");              throw DudleyException(msgPrefix+"add_var(Points_Nodes)");
375            if (!(ids->put(mesh->Points->Nodes, num_Points)))
376     }              throw DudleyException(msgPrefix+"put(Points_Nodes)");
377        }
378     // // // // // Points // // // // //  
379        // // // // // TagMap // // // // //
380     if (num_Points>0) {      if (num_Tags > 0) {
381            // Temp storage to gather node IDs
382        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");          vector<int> Tags_keys;
383    
384        // Points_Id          // Copy tag data into temp arrays
385        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )          TagMap::const_iterator it;
386           throw DudleyException(msgPrefix+"add_var(Points_Id)");          for (it = mesh->tagMap.begin(); it != mesh->tagMap.end(); it++) {
387        int_ptr = &mesh->Points->Id[0];              Tags_keys.push_back(it->second);
388        if (! (ids->put(int_ptr, num_Points)) )          }
          throw DudleyException(msgPrefix+"put(Points_Id)");  
   
       // Points_Tag  
       if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )  
          throw DudleyException(msgPrefix+"add_var(Points_Tag)");  
       int_ptr = &mesh->Points->Tag[0];  
       if (! (ids->put(int_ptr, num_Points)) )  
          throw DudleyException(msgPrefix+"put(Points_Tag)");  
   
       // Points_Owner  
       if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )  
          throw DudleyException(msgPrefix+"add_var(Points_Owner)");  
       int_ptr = &mesh->Points->Owner[0];  
       if (! (ids->put(int_ptr, num_Points)) )  
          throw DudleyException(msgPrefix+"put(Points_Owner)");  
   
       // Points_Color  
       if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )  
          throw DudleyException(msgPrefix+"add_var(Points_Color)");  
       int_ptr = &mesh->Points->Color[0];  
       if (! (ids->put(int_ptr, num_Points)) )  
          throw DudleyException(msgPrefix+"put(Points_Color)");  
   
       // Points_Nodes  
       // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]  
       if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )  
          throw DudleyException(msgPrefix+"add_var(Points_Nodes)");  
       if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )  
          throw DudleyException(msgPrefix+"put(Points_Nodes)");  
   
    }  
   
    // // // // // TagMap // // // // //  
   
    if (num_Tags>0) {  
   
       // Temp storage to gather node IDs  
       int *Tags_keys = new int[num_Tags];  
       char name_temp[4096];  
   
       /* Copy tag data into temp arrays */  
       tag_map = mesh->TagMap;  
       if (tag_map) {  
          int i = 0;  
          while (tag_map) {  
             Tags_keys[i++] = tag_map->tag_key;  
             tag_map=tag_map->next;  
          }  
       }  
   
       // Tags_keys  
       if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )  
          throw DudleyException(msgPrefix+"add_var(Tags_keys)");  
       int_ptr = &Tags_keys[0];  
       if (! (ids->put(int_ptr, num_Tags)) )  
          throw DudleyException(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  
       tag_map = mesh->TagMap;  
       if (tag_map) {  
          int i = 0;  
          while (tag_map) {  
             sprintf(name_temp, "Tags_name_%d", i);  
             if (!dataFile.add_att(name_temp, tag_map->name) )  
                throw DudleyException(msgPrefix+"add_att(Tags_names_XX)");  
             tag_map=tag_map->next;  
             i++;  
          }  
       }  
389    
390        delete[] Tags_keys;          // Tags_keys
391     }          if (!(ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])))
392                throw DudleyException(msgPrefix+"add_var(Tags_keys)");
393            if (!(ids->put(&Tags_keys[0], num_Tags)))
394                throw DudleyException(msgPrefix+"put(Tags_keys)");
395    
396            // Tags_names_*
397            // This is an array of strings, it should be stored as an array but
398            // instead I have hacked in one attribute per string because the NetCDF
399            // manual doesn't tell how to do an array of strings
400            int i = 0;
401            for (it = mesh->tagMap.begin(); it != mesh->tagMap.end(); it++, i++) {
402                stringstream ss;
403                ss << "Tags_name_" << i;
404                const string name(ss.str());
405                if (!dataFile.add_att(name.c_str(), it->first.c_str()))
406                    throw DudleyException(msgPrefix+"add_att(Tags_names_XX)");
407            }
408        }
409    
410  /* Send token to next MPI process so he can take his turn */      // Send token to next MPI process so he can take his turn
411  #ifdef ESYS_MPI  #ifdef ESYS_MPI
412     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);      if (mpi_rank < mpi_size-1)
413            MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, getMPIComm());
414  #endif  #endif
415    
416     // NetCDF file is closed by destructor of NcFile object      // NetCDF file is closed by destructor of NcFile object
417    
418  #else  #else // USE_NETCDF
419     throw DudleyException("MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");      throw DudleyException("MeshAdapter::dump: not configured with netCDF. "
420  #endif  /* USE_NETCDF */                            "Please contact your installation manager.");
421    #endif // USE_NETCDF
422  }  }
423    
424  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
# Line 491  string MeshAdapter::getDescription() con Line 428  string MeshAdapter::getDescription() con
428    
429  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
430  {  {
431     FunctionSpaceNamesMapType::iterator loc;      FunctionSpaceNamesMapType::iterator loc;
432     loc=m_functionSpaceTypeNames.find(functionSpaceType);      loc = m_functionSpaceTypeNames.find(functionSpaceType);
433     if (loc==m_functionSpaceTypeNames.end()) {      if (loc == m_functionSpaceTypeNames.end()) {
434        return "Invalid function space type code.";          return "Invalid function space type code.";
435     } else {      } else {
436        return loc->second;          return loc->second;
437     }      }
438  }  }
439    
440  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
441  {  {
442      FunctionSpaceNamesMapType::iterator loc;      FunctionSpaceNamesMapType::iterator loc;
443      loc=m_functionSpaceTypeNames.find(functionSpaceType);      loc = m_functionSpaceTypeNames.find(functionSpaceType);
444      return (loc!=m_functionSpaceTypeNames.end());      return (loc != m_functionSpaceTypeNames.end());
445  }  }
446    
447  void MeshAdapter::setFunctionSpaceTypeNames()  void MeshAdapter::setFunctionSpaceTypeNames()
# Line 512  void MeshAdapter::setFunctionSpaceTypeNa Line 449  void MeshAdapter::setFunctionSpaceTypeNa
449      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
450                  DegreesOfFreedom,"Dudley_DegreesOfFreedom [Solution(domain)]"));                  DegreesOfFreedom,"Dudley_DegreesOfFreedom [Solution(domain)]"));
451      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
                 ReducedDegreesOfFreedom,"Dudley_ReducedDegreesOfFreedom [ReducedSolution(domain)]"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(  
452                  Nodes,"Dudley_Nodes [ContinuousFunction(domain)]"));                  Nodes,"Dudley_Nodes [ContinuousFunction(domain)]"));
453      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
                 ReducedNodes,"Dudley_Reduced_Nodes [ReducedContinuousFunction(domain)]"));  
     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(  
454                  Elements,"Dudley_Elements [Function(domain)]"));                  Elements,"Dudley_Elements [Function(domain)]"));
455      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(      m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(
456                  ReducedElements,"Dudley_Reduced_Elements [ReducedFunction(domain)]"));                  ReducedElements,"Dudley_Reduced_Elements [ReducedFunction(domain)]"));
# Line 536  int MeshAdapter::getContinuousFunctionCo Line 469  int MeshAdapter::getContinuousFunctionCo
469    
470  int MeshAdapter::getReducedContinuousFunctionCode() const  int MeshAdapter::getReducedContinuousFunctionCode() const
471  {  {
472      return ReducedNodes;      return Nodes;
473  }  }
474    
475  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
# Line 586  int MeshAdapter::getSolutionCode() const Line 519  int MeshAdapter::getSolutionCode() const
519    
520  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
521  {  {
522      return ReducedDegreesOfFreedom;      return DegreesOfFreedom;
523  }  }
524    
525  int MeshAdapter::getDiracDeltaFunctionsCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
# Line 594  int MeshAdapter::getDiracDeltaFunctionsC Line 527  int MeshAdapter::getDiracDeltaFunctionsC
527      return Points;      return Points;
528  }  }
529    
 //  
 // return the spatial dimension of the Mesh:  
 //  
530  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
531  {  {
532      Dudley_Mesh* mesh=m_dudleyMesh.get();      return m_dudleyMesh->getDim();
     int numDim=Dudley_Mesh_getDim(mesh);  
     return numDim;  
533  }  }
534    
535  //  //
# Line 609  int MeshAdapter::getDim() const Line 537  int MeshAdapter::getDim() const
537  //  //
538  int MeshAdapter::getNumDataPointsGlobal() const  int MeshAdapter::getNumDataPointsGlobal() const
539  {  {
540      return Dudley_NodeFile_getGlobalNumNodes(m_dudleyMesh.get()->Nodes);      return m_dudleyMesh->Nodes->getGlobalNumNodes();
541  }  }
542    
543  //  //
544  // return the number of data points per sample and the number of samples  // return the number of data points per sample and the number of samples
545  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
546  //  //
547  pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const  pair<int,dim_t> MeshAdapter::getDataShape(int functionSpaceCode) const
548  {  {
549     int numDataPointsPerSample=0;      int numDataPointsPerSample = 0;
550     int numSamples=0;      dim_t numSamples = 0;
551     Dudley_Mesh* mesh=m_dudleyMesh.get();      Mesh* mesh = getMesh();
552     switch (functionSpaceCode) {      switch (functionSpaceCode) {
553     case(Nodes):          case Nodes:
554     numDataPointsPerSample=1;              numDataPointsPerSample = 1;
555     numSamples=Dudley_NodeFile_getNumNodes(mesh->Nodes);              numSamples = mesh->Nodes->getNumNodes();
556     break;          break;
557     case(ReducedNodes):          case Elements:
558     numDataPointsPerSample=1;              if (mesh->Elements) {
559     numSamples=Dudley_NodeFile_getNumReducedNodes(mesh->Nodes);                  numSamples = mesh->Elements->numElements;
560     break;                  numDataPointsPerSample = mesh->Elements->numLocalDim+1;
561     case(Elements):              }
562     if (mesh->Elements!=NULL) {          break;
563        numSamples=mesh->Elements->numElements;          case ReducedElements:
564        numDataPointsPerSample=mesh->Elements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;              if (mesh->Elements) {
565     }                  numSamples = mesh->Elements->numElements;
566     break;                  numDataPointsPerSample =(mesh->Elements->numLocalDim==0)?0:1;
567     case(ReducedElements):              }
568     if (mesh->Elements!=NULL) {          break;
569        numSamples=mesh->Elements->numElements;          case FaceElements:
570        numDataPointsPerSample=(mesh->Elements->numLocalDim==0)?0:1;              if (mesh->FaceElements) {
571     }                  numDataPointsPerSample = mesh->FaceElements->numLocalDim+1;
572     break;                  numSamples = mesh->FaceElements->numElements;
573     case(FaceElements):              }
574     if (mesh->FaceElements!=NULL) {          break;
575        numDataPointsPerSample=mesh->FaceElements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;          case ReducedFaceElements:
576        numSamples=mesh->FaceElements->numElements;              if (mesh->FaceElements) {
577     }                  numDataPointsPerSample = (mesh->FaceElements->numLocalDim==0)?0:1;
578     break;                  numSamples = mesh->FaceElements->numElements;
579     case(ReducedFaceElements):              }
580     if (mesh->FaceElements!=NULL) {          break;
581        numDataPointsPerSample=(mesh->FaceElements->numLocalDim==0)?0:1/*referenceElementSet->referenceElementReducedQuadrature->BasisFunctions->numQuadNodes*/;          case Points:
582        numSamples=mesh->FaceElements->numElements;              if (mesh->Points) {
583     }                  numDataPointsPerSample = 1;
584     break;                  numSamples = mesh->Points->numElements;
585     case(Points):              }
586     if (mesh->Points!=NULL) {          break;
587        numDataPointsPerSample=1;          case DegreesOfFreedom:
588        numSamples=mesh->Points->numElements;              if (mesh->Nodes) {
589     }                  numDataPointsPerSample = 1;
590     break;                  numSamples = mesh->Nodes->getNumDegreesOfFreedom();
591     case(DegreesOfFreedom):              }
592     if (mesh->Nodes!=NULL) {          break;
593        numDataPointsPerSample=1;          default:
594        numSamples=Dudley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);              stringstream ss;
595     }              ss << "Invalid function space type: " << functionSpaceCode
596     break;                  << " for domain " << getDescription();
597     case(ReducedDegreesOfFreedom):              throw DudleyException(ss.str());
598     if (mesh->Nodes!=NULL) {              break;
599        numDataPointsPerSample=1;      }
600        numSamples=Dudley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);      return pair<int,dim_t>(numDataPointsPerSample,numSamples);
    }  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();  
       throw DudleyException(temp.str());  
       break;  
    }  
    return pair<int,int>(numDataPointsPerSample,numSamples);  
601  }  }
602    
603  //  //
# Line 702  void MeshAdapter::addPDEToSystem( Line 621  void MeshAdapter::addPDEToSystem(
621      }      }
622  #endif  #endif
623    
624      Dudley_Mesh* mesh = m_dudleyMesh.get();      Mesh* mesh = m_dudleyMesh.get();
625      Assemble_PDE(mesh->Nodes, mesh->Elements, mat.getPtr(), rhs,      Assemble_PDE(mesh->Nodes, mesh->Elements, mat.getPtr(), rhs,
626                   A, B, C, D, X, Y);                   A, B, C, D, X, Y);
627      Assemble_PDE(mesh->Nodes, mesh->FaceElements, mat.getPtr(), rhs,      Assemble_PDE(mesh->Nodes, mesh->FaceElements, mat.getPtr(), rhs,
# Line 725  void MeshAdapter::addPDEToLumpedSystem(e Line 644  void MeshAdapter::addPDEToLumpedSystem(e
644                                         const escript::Data& d_dirac,                                         const escript::Data& d_dirac,
645                                         bool useHRZ) const                                         bool useHRZ) const
646  {  {
647      Dudley_Mesh* mesh = m_dudleyMesh.get();      Mesh* mesh = m_dudleyMesh.get();
648      Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);      Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);
649      Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, mat, d, useHRZ);      Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, mat, d, useHRZ);
650      Assemble_LumpedSystem(mesh->Nodes, mesh->Points, mat, d_dirac, useHRZ);      Assemble_LumpedSystem(mesh->Nodes, mesh->Points, mat, d_dirac, useHRZ);
# Line 741  void MeshAdapter::addPDEToRHS(escript::D Line 660  void MeshAdapter::addPDEToRHS(escript::D
660      if (!y_contact.isEmpty())      if (!y_contact.isEmpty())
661          throw DudleyException("Dudley does not support y_contact");          throw DudleyException("Dudley does not support y_contact");
662    
663      Dudley_Mesh* mesh=m_dudleyMesh.get();      Mesh* mesh=m_dudleyMesh.get();
664    
665      Assemble_PDE(mesh->Nodes, mesh->Elements, escript::ASM_ptr(), rhs,      Assemble_PDE(mesh->Nodes, mesh->Elements, escript::ASM_ptr(), rhs,
666                   escript::Data(), escript::Data(), escript::Data(),                   escript::Data(), escript::Data(), escript::Data(),
# Line 778  void MeshAdapter::addPDEToTransportProbl Line 697  void MeshAdapter::addPDEToTransportProbl
697    
698      source.expand();      source.expand();
699    
700      Dudley_Mesh* mesh=m_dudleyMesh.get();      Mesh* mesh = m_dudleyMesh.get();
701    
702      Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowMassMatrix(), source,      Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowMassMatrix(), source,
703                   escript::Data(), escript::Data(), escript::Data(), M,                   escript::Data(), escript::Data(), escript::Data(), M,
# Line 799  void MeshAdapter::addPDEToTransportProbl Line 718  void MeshAdapter::addPDEToTransportProbl
718  //  //
719  // interpolates data between different function spaces:  // interpolates data between different function spaces:
720  //  //
721  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,
722                                          const escript::Data& in) const
723  {  {
724     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));      if (*in.getFunctionSpace().getDomain() != *this)  
725     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));          throw DudleyException("Illegal domain of interpolant.");
726     if (inDomain!=*this)        if (*target.getFunctionSpace().getDomain() != *this)
727        throw DudleyException("Illegal domain of interpolant.");          throw DudleyException("Illegal domain of interpolation target.");
728     if (targetDomain!=*this)  
729        throw DudleyException("Illegal domain of interpolation target.");      Mesh* mesh = m_dudleyMesh.get();
730        switch (in.getFunctionSpace().getTypeCode()) {
731     Dudley_Mesh* mesh=m_dudleyMesh.get();          case Nodes:
732     switch(in.getFunctionSpace().getTypeCode()) {              switch (target.getFunctionSpace().getTypeCode()) {
733     case(Nodes):                  case Nodes:
734        switch(target.getFunctionSpace().getTypeCode()) {                  case DegreesOfFreedom:
735        case(Nodes):                      Assemble_CopyNodalData(mesh->Nodes, target, in);
736        case(ReducedNodes):                  break;
737        case(DegreesOfFreedom):                  case Elements:
738        case(ReducedDegreesOfFreedom):                  case ReducedElements:
739        Assemble_CopyNodalData(mesh->Nodes,&target,&in);                      Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
740        break;                  break;
741        case(Elements):                  case FaceElements:
742        case(ReducedElements):                  case ReducedFaceElements:
743        Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);                      Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
744        break;                  break;
745        case(FaceElements):                  case Points:
746        case(ReducedFaceElements):                      Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
747        Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);                  break;
748        break;                  default:
749        case(Points):                      stringstream ss;
750        Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);                      ss << "interpolateOnDomain: Dudley does not know anything "
751        break;                            "about function space type "
752        default:                            << target.getFunctionSpace().getTypeCode();
753           stringstream temp;                      throw DudleyException(ss.str());
754           temp << "Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                      break;
755           throw DudleyException(temp.str());              }
756           break;          break;
757        }          case Elements:
758        break;              if (target.getFunctionSpace().getTypeCode() == Elements) {
759     case(ReducedNodes):                  Assemble_CopyElementData(mesh->Elements, target, in);
760        switch(target.getFunctionSpace().getTypeCode()) {              } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
761        case(Nodes):                  Assemble_AverageElementData(mesh->Elements, target, in);
762        case(ReducedNodes):              } else {
763        case(DegreesOfFreedom):                  throw DudleyException("No interpolation with data on elements possible.");
764        case(ReducedDegreesOfFreedom):              }
765        Assemble_CopyNodalData(mesh->Nodes,&target,&in);              break;
766        break;          case ReducedElements:
767        case(Elements):              if (target.getFunctionSpace().getTypeCode() == ReducedElements) {
768        case(ReducedElements):                  Assemble_CopyElementData(mesh->Elements, target, in);
769        Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);              } else {
770        break;                  throw DudleyException("No interpolation with data on elements "
771        case(FaceElements):                                     "with reduced integration order possible.");
772        case(ReducedFaceElements):              }
773        Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);              break;
774        break;          case FaceElements:
775        case(Points):              if (target.getFunctionSpace().getTypeCode() == FaceElements) {
776        Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);                  Assemble_CopyElementData(mesh->FaceElements, target, in);
777        break;              } else if (target.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
778        default:                  Assemble_AverageElementData(mesh->FaceElements, target, in);
779           stringstream temp;              } else {
780           temp << "Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                  throw DudleyException("No interpolation with data on face elements possible.");
781           throw DudleyException(temp.str());              }
782           break;              break;
783        }          case ReducedFaceElements:
784        break;              if (target.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
785     case(Elements):                  Assemble_CopyElementData(mesh->FaceElements, target, in);
786        if (target.getFunctionSpace().getTypeCode()==Elements) {              } else {
787           Assemble_CopyElementData(mesh->Elements,&target,&in);                  throw DudleyException("No interpolation with data on face "
788        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {                            "elements with reduced integration order possible.");
789           Assemble_AverageElementData(mesh->Elements,&target,&in);              }
790        } else {              break;
791           throw DudleyException("No interpolation with data on elements possible.");          case Points:
792        }              if (target.getFunctionSpace().getTypeCode() == Points) {
793        break;                  Assemble_CopyElementData(mesh->Points, target, in);
794     case(ReducedElements):              } else {
795        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {                  throw DudleyException("No interpolation with data on points possible.");
796           Assemble_CopyElementData(mesh->Elements,&target,&in);              }
797        } else {              break;
798           throw DudleyException("No interpolation with data on elements with reduced integration order possible.");          case DegreesOfFreedom:
799        }              switch (target.getFunctionSpace().getTypeCode()) {
800        break;                  case DegreesOfFreedom:
801     case(FaceElements):                  Assemble_CopyNodalData(mesh->Nodes, target, in);
802        if (target.getFunctionSpace().getTypeCode()==FaceElements) {                  break;
803           Assemble_CopyElementData(mesh->FaceElements,&target,&in);              
804        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {                  case Nodes:
805           Assemble_AverageElementData(mesh->FaceElements,&target,&in);                  if (getMPISize() > 1) {
806        } else {                      escript::Data temp = escript::Data(in);
807           throw DudleyException("No interpolation with data on face elements possible.");                      temp.expand();
808        }                      Assemble_CopyNodalData(mesh->Nodes, target, temp);
809        break;                  } else {
810     case(ReducedFaceElements):                      Assemble_CopyNodalData(mesh->Nodes, target, in);
811        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {                  }
812           Assemble_CopyElementData(mesh->FaceElements,&target,&in);                  break;
813        } else {                  case Elements:
814           throw DudleyException("No interpolation with data on face elements with reduced integration order possible.");                  case ReducedElements:
815        }                  if (getMPISize() > 1) {
816        break;                      escript::Data temp = escript::Data(in, continuousFunction(*this) );
817     case(Points):                      Assemble_interpolate(mesh->Nodes, mesh->Elements, temp, target);
818        if (target.getFunctionSpace().getTypeCode()==Points) {                  } else {
819           Assemble_CopyElementData(mesh->Points,&target,&in);                      Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
820        } else {                  }
821           throw DudleyException("No interpolation with data on points possible.");                  break;
822        }                  case FaceElements:
823        break;                  case ReducedFaceElements:
824     case(DegreesOfFreedom):                        if (getMPISize() > 1) {
825        switch(target.getFunctionSpace().getTypeCode()) {                      escript::Data temp = escript::Data(in, continuousFunction(*this) );
826        case(ReducedDegreesOfFreedom):                      Assemble_interpolate(mesh->Nodes, mesh->FaceElements, temp, target);
827        case(DegreesOfFreedom):                  } else {
828        Assemble_CopyNodalData(mesh->Nodes,&target,&in);                      Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
829        break;                  }
830                      break;
831        case(Nodes):                  case Points:
832        case(ReducedNodes):                  if (getMPISize() > 1) {
833        if (getMPISize()>1) {                      //escript::Data temp=escript::Data(in, continuousFunction(*this) );
834           escript::Data temp=escript::Data(in);                  } else {
835           temp.expand();                      Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
836           Assemble_CopyNodalData(mesh->Nodes,&target,&temp);                  }
837        } else {                  break;
838           Assemble_CopyNodalData(mesh->Nodes,&target,&in);                  default:
839        }                      stringstream ss;
840        break;                      ss << "interpolateOnDomain: Dudley does not know anything "
841        case(Elements):                            "about function space type "
842        case(ReducedElements):                         << target.getFunctionSpace().getTypeCode();
843        if (getMPISize()>1) {                      throw DudleyException(ss.str());
844           escript::Data temp=escript::Data( in,  continuousFunction(*this) );                      break;
845           Assemble_interpolate(mesh->Nodes,mesh->Elements,&temp,&target);              }
846        } else {              break;
847           Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);         default:
848        }            stringstream ss;
849        break;            ss << "interpolateOnDomain: Dudley does not know anything about "
850        case(FaceElements):                  "function space type " << in.getFunctionSpace().getTypeCode();
851        case(ReducedFaceElements):            throw DudleyException(ss.str());
852        if (getMPISize()>1) {            break;
853           escript::Data temp=escript::Data( in,  continuousFunction(*this) );      }
          Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&temp,&target);  
     
       } else {  
          Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);  
       }  
       break;  
       case(Points):  
       if (getMPISize()>1) {  
          //escript::Data temp=escript::Data( in,  continuousFunction(*this) );  
          //escriptDataC _in2 = temp.getDataC();  
       } else {  
          Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);  
       }  
       break;  
       default:  
          stringstream temp;  
          temp << "Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
          throw DudleyException(temp.str());  
          break;  
       }  
       break;  
    case(ReducedDegreesOfFreedom):  
       switch(target.getFunctionSpace().getTypeCode()) {  
       case(Nodes):  
       throw DudleyException("Dudley 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();  
          Assemble_CopyNodalData(mesh->Nodes,&target,&temp);  
       } else {  
          Assemble_CopyNodalData(mesh->Nodes,&target,&in);  
       }  
       break;  
       case(DegreesOfFreedom):  
       throw DudleyException("Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
       break;  
       case(ReducedDegreesOfFreedom):  
       Assemble_CopyNodalData(mesh->Nodes,&target,&in);  
       break;  
       case(Elements):  
       case(ReducedElements):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );  
          Assemble_interpolate(mesh->Nodes,mesh->Elements,&temp,&target);  
       } else {  
          Assemble_interpolate(mesh->Nodes,mesh->Elements,&in,&target);  
       }  
       break;  
       case(FaceElements):  
       case(ReducedFaceElements):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );  
          Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&temp,&target);  
       } else {  
          Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&in,&target);  
       }  
       break;  
       case(Points):  
       if (getMPISize()>1) {  
          escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );  
          Assemble_interpolate(mesh->Nodes,mesh->Points,&temp,&target);  
       } else {  
          Assemble_interpolate(mesh->Nodes,mesh->Points,&in,&target);  
       }  
       break;  
       default:  
          stringstream temp;  
          temp << "Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
          throw DudleyException(temp.str());  
          break;  
       }  
       break;  
    default:  
       stringstream temp;  
       temp << "Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();  
       throw DudleyException(temp.str());  
       break;  
    }  
854  }  }
855    
856  //  //
# Line 1018  void MeshAdapter::interpolateOnDomain(es Line 858  void MeshAdapter::interpolateOnDomain(es
858  //  //
859  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
860  {  {
861     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      if (*arg.getFunctionSpace().getDomain() != *this)
862     if (argDomain!=*this)          throw DudleyException("setToX: Illegal domain of data point locations");
863        throw DudleyException("Illegal domain of data point locations");  
864     Dudley_Mesh* mesh=m_dudleyMesh.get();      Mesh* mesh = m_dudleyMesh.get();
865     // in case of values node coordinates we can do the job directly:      // in case of appropriate function space we can do the job directly:
866     if (arg.getFunctionSpace().getTypeCode()==Nodes) {      if (arg.getFunctionSpace().getTypeCode() == Nodes) {
867        Assemble_NodeCoordinates(mesh->Nodes,&arg);          Assemble_NodeCoordinates(mesh->Nodes, arg);
868     } else {      } else {
869        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);          escript::Data tmp_data = Vector(0., continuousFunction(*this), true);
870        Assemble_NodeCoordinates(mesh->Nodes,&tmp_data);          Assemble_NodeCoordinates(mesh->Nodes, tmp_data);
871        // this is then interpolated onto arg:          // this is then interpolated onto arg:
872        interpolateOnDomain(arg,tmp_data);          interpolateOnDomain(arg, tmp_data);
873     }      }
874  }  }
875    
876  //  //
# Line 1038  void MeshAdapter::setToX(escript::Data& Line 878  void MeshAdapter::setToX(escript::Data&
878  //  //
879  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
880  {  {
881  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/      if (*normal.getFunctionSpace().getDomain() != *this)
882     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));          throw ValueError("setToNormal: Illegal domain of normal locations");
883     if (normalDomain!=*this)  
884        throw DudleyException("Illegal domain of normal locations");      Mesh* mesh=m_dudleyMesh.get();
885     Dudley_Mesh* mesh=m_dudleyMesh.get();      if (normal.getFunctionSpace().getTypeCode() == FaceElements ||
886     switch(normal.getFunctionSpace().getTypeCode()) {              normal.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
887     case(Nodes):          Assemble_getNormal(mesh->Nodes, mesh->FaceElements, normal);
888     throw DudleyException("Dudley does not support surface normal vectors for nodes");      } else {
889     break;          stringstream ss;
890     case(ReducedNodes):          ss << "setToNormal: Illegal function space type "
891     throw DudleyException("Dudley does not support surface normal vectors for reduced nodes");             << normal.getFunctionSpace().getTypeCode();
892     break;          throw ValueError(ss.str());
893     case(Elements):      }
    throw DudleyException("Dudley does not support surface normal vectors for elements");  
    break;  
    case(ReducedElements):  
    throw DudleyException("Dudley does not support surface normal vectors for elements with reduced integration order");  
    break;  
    case (FaceElements):  
    Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&normal);  
    break;  
    case (ReducedFaceElements):  
    Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&normal);  
    break;  
    case(Points):  
    throw DudleyException("Dudley does not support surface normal vectors for point elements");  
    break;  
    case(DegreesOfFreedom):  
    throw DudleyException("Dudley does not support surface normal vectors for degrees of freedom.");  
    break;  
    case(ReducedDegreesOfFreedom):  
    throw DudleyException("Dudley does not support surface normal vectors for reduced degrees of freedom.");  
    break;  
    default:  
       stringstream temp;  
       temp << "Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();  
       throw DudleyException(temp.str());  
       break;  
    }  
894  }  }
895    
896  //  //
897  // interpolates data to other domain  // interpolates data to other domain
898  //  //
899  void MeshAdapter::interpolateAcross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateAcross(escript::Data& target,
900                                        const escript::Data& source) const
901  {  {
902      throw escript::NotImplementedError("Dudley does not allow interpolation "      throw escript::NotImplementedError("Dudley does not allow interpolation "
903                                         "across domains.");                                         "across domains.");
# Line 1091  void MeshAdapter::interpolateAcross(escr Line 906  void MeshAdapter::interpolateAcross(escr
906  //  //
907  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
908  //  //
909  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,
910                                     const escript::Data& arg) const
911  {  {
912     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      if (*arg.getFunctionSpace().getDomain() != *this)
913     if (argDomain!=*this)          throw ValueError("setToIntegrals: Illegal domain of integration kernel");
914        throw DudleyException("Illegal domain of integration kernel");  
915        Mesh* mesh = m_dudleyMesh.get();
916     Dudley_Mesh* mesh=m_dudleyMesh.get();      switch (arg.getFunctionSpace().getTypeCode()) {
917     escript::Data temp;          case Nodes: // fall through
918     switch(arg.getFunctionSpace().getTypeCode()) {          case DegreesOfFreedom:
919     case(Nodes):          {
920     temp=escript::Data( arg, escript::function(*this) );              escript::Data temp(arg, escript::function(*this));
921     Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);              Assemble_integrate(mesh->Nodes, mesh->Elements, temp, integrals);
922     break;          }
923     case(ReducedNodes):          break;
924     temp=escript::Data( arg, escript::function(*this) );          case Elements: // fall through
925     Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);          case ReducedElements:
926     break;              Assemble_integrate(mesh->Nodes,mesh->Elements, arg, integrals);
927     case(Elements):          break;
928     Assemble_integrate(mesh->Nodes,mesh->Elements,&arg,&integrals[0]);          case FaceElements: // fall through
929     break;          case ReducedFaceElements:
930     case(ReducedElements):              Assemble_integrate(mesh->Nodes,mesh->FaceElements, arg, integrals);
931     Assemble_integrate(mesh->Nodes,mesh->Elements,&arg,&integrals[0]);          break;
932     break;          case Points:
933     case(FaceElements):              throw ValueError("Integral of data on points is not supported.");
934     Assemble_integrate(mesh->Nodes,mesh->FaceElements,&arg,&integrals[0]);          break;
935     break;          default:
936     case(ReducedFaceElements):              stringstream ss;
937     Assemble_integrate(mesh->Nodes,mesh->FaceElements,&arg,&integrals[0]);              ss << "setToIntegrals: Dudley does not know anything about "
938     break;                  "function space type " << arg.getFunctionSpace().getTypeCode();
939     case(Points):              throw DudleyException(ss.str());
940     throw DudleyException("Integral of data on points is not supported.");      }
    break;  
    case(DegreesOfFreedom):  
    temp=escript::Data( arg, escript::function(*this) );  
    Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);  
    break;  
    case(ReducedDegreesOfFreedom):  
    temp=escript::Data( arg, escript::function(*this) );  
    Assemble_integrate(mesh->Nodes,mesh->Elements,&temp,&integrals[0]);  
    break;  
    default:  
       stringstream temp;  
       temp << "Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();  
       throw DudleyException(temp.str());  
       break;  
    }  
941  }  }
942    
943  //  //
944  // calculates the gradient of arg:  // calculates the gradient of arg:
945  //  //
946  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad, const escript::Data& arg) const
947  {  {
948     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      if (*arg.getFunctionSpace().getDomain() != *this)
949     if (argDomain!=*this)          throw ValueError("setToGradient: Illegal domain of gradient argument");
950        throw DudleyException("Illegal domain of gradient argument");      if (*grad.getFunctionSpace().getDomain() != *this)
951     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));          throw ValueError("setToGradient: Illegal domain of gradient");
952     if (gradDomain!=*this)  
953        throw DudleyException("Illegal domain of gradient");      Mesh* mesh = m_dudleyMesh.get();
954        escript::Data nodeData;
955     Dudley_Mesh* mesh=m_dudleyMesh.get();      if (getMPISize() > 1) {
956     const escript::Data* nodeData=0;          if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
957     escript::Data temp;              nodeData = escript::Data(arg, continuousFunction(*this));
958     if (getMPISize()>1) {          } else {
959        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {              nodeData = arg;
960           temp=escript::Data( arg,  continuousFunction(*this) );          }
961           nodeData = &temp;      } else {
962        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {          nodeData = arg;
963           temp=escript::Data( arg,  reducedContinuousFunction(*this) );      }
964           nodeData = &temp;      switch (grad.getFunctionSpace().getTypeCode()) {
965        } else {          case Nodes:
966           nodeData = &arg;              throw DudleyException("Gradient at nodes is not supported.");
967        }          break;
968     } else {          case Elements:
969        nodeData = &arg;              Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
970     }          break;
971     switch(grad.getFunctionSpace().getTypeCode()) {          case ReducedElements:
972     case(Nodes):              Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
973     throw DudleyException("Gradient at nodes is not supported.");          break;
974     break;          case FaceElements:
975     case(ReducedNodes):              Assemble_gradient(mesh->Nodes,mesh->FaceElements, grad, nodeData);
976     throw DudleyException("Gradient at reduced nodes is not supported.");          break;
977     break;          case ReducedFaceElements:
978     case(Elements):              Assemble_gradient(mesh->Nodes, mesh->FaceElements, grad, nodeData);
979     Assemble_gradient(mesh->Nodes,mesh->Elements,&grad, nodeData);          break;
980     break;          case Points:
981     case(ReducedElements):              throw DudleyException("Gradient at points is not supported.");
982     Assemble_gradient(mesh->Nodes,mesh->Elements,&grad, nodeData);          break;
983     break;          case DegreesOfFreedom:
984     case(FaceElements):              throw DudleyException("Gradient at degrees of freedom is not supported.");
985     Assemble_gradient(mesh->Nodes,mesh->FaceElements,&grad, nodeData);          break;
986     break;          default:
987     case(ReducedFaceElements):              stringstream ss;
988     Assemble_gradient(mesh->Nodes,mesh->FaceElements,&grad, nodeData);              ss << "Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
989     break;              throw DudleyException(ss.str());
990     case(Points):      }
    throw DudleyException("Gradient at points is not supported.");  
    break;  
    case(DegreesOfFreedom):  
    throw DudleyException("Gradient at degrees of freedom is not supported.");  
    break;  
    case(ReducedDegreesOfFreedom):  
    throw DudleyException("Gradient at reduced degrees of freedom is not supported.");  
    break;  
    default:  
       stringstream temp;  
       temp << "Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();  
       throw DudleyException(temp.str());  
       break;  
    }  
991  }  }
992    
993  //  //
994  // returns the size of elements:  // returns the size of elements
995  //  //
996  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
997  {  {
998     Dudley_Mesh* mesh=m_dudleyMesh.get();      Mesh* mesh=m_dudleyMesh.get();
999     switch(size.getFunctionSpace().getTypeCode()) {      switch (size.getFunctionSpace().getTypeCode()) {
1000     case(Nodes):          case Nodes:
1001     throw DudleyException("Size of nodes is not supported.");              throw DudleyException("Size of nodes is not supported.");
1002     break;          break;
1003     case(ReducedNodes):          case Elements:
1004     throw DudleyException("Size of reduced nodes is not supported.");              Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1005     break;          break;
1006     case(Elements):          case ReducedElements:
1007     Assemble_getSize(mesh->Nodes,mesh->Elements,&size);              Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1008     break;          break;
1009     case(ReducedElements):          case FaceElements:
1010     Assemble_getSize(mesh->Nodes,mesh->Elements,&size);              Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1011     break;          break;
1012     case(FaceElements):          case ReducedFaceElements:
1013     Assemble_getSize(mesh->Nodes,mesh->FaceElements,&size);              Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1014     break;          break;
1015     case(ReducedFaceElements):          case Points:
1016     Assemble_getSize(mesh->Nodes,mesh->FaceElements,&size);              throw DudleyException("Size of point elements is not supported.");
1017     break;          break;
1018     case(Points):          case DegreesOfFreedom:
1019     throw DudleyException("Size of point elements is not supported.");              throw DudleyException("Size of degrees of freedom is not supported.");
1020     break;          break;
1021     case(DegreesOfFreedom):          default:
1022     throw DudleyException("Size of degrees of freedom is not supported.");              stringstream ss;
1023     break;              ss << "setToSize: Dudley does not know anything about function "
1024     case(ReducedDegreesOfFreedom):                    "space type " << size.getFunctionSpace().getTypeCode();
1025     throw DudleyException("Size of reduced degrees of freedom is not supported.");              throw ValueError(ss.str());
1026     break;      }
    default:  
       stringstream temp;  
       temp << "Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();  
       throw DudleyException(temp.str());  
       break;  
    }  
1027  }  }
1028    
1029  //  //
1030  // sets the location of nodes  // sets the location of nodes
1031  //  //
1032  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& newX)
1033  {  {
1034     Dudley_Mesh* mesh=m_dudleyMesh.get();      if (*newX.getFunctionSpace().getDomain() != *this)
1035     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));          throw DudleyException("Illegal domain of new point locations");
    if (newDomain!=*this)  
       throw DudleyException("Illegal domain of new point locations");  
    if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {  
        Dudley_Mesh_setCoordinates(mesh,&new_x);  
    } else {  
        throw DudleyException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");        
1036    
1037     }      if (newX.getFunctionSpace() == continuousFunction(*this)) {
1038            m_dudleyMesh->setCoordinates(newX);
1039        } else {
1040            throw DudleyException("As of escript version 3.3 - setNewX only "
1041                    "accepts ContinuousFunction arguments. Please interpolate.");
1042        }
1043  }  }
1044    
1045  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1046  {  {
     if (getMPISize()>1) {  
1047  #ifdef ESYS_MPI  #ifdef ESYS_MPI
1048          index_t myFirstNode=0, myLastNode=0, k=0;      if (getMPISize() > 1) {
1049          index_t* globalNodeIndex=0;          Mesh* mesh = m_dudleyMesh.get();
1050          Dudley_Mesh* mesh_p=m_dudleyMesh.get();          if (fs_code == DUDLEY_NODES) {
1051          if (fs_code == DUDLEY_REDUCED_NODES)              const index_t myFirstNode = mesh->Nodes->getFirstNode();
1052          {              const index_t myLastNode = mesh->Nodes->getLastNode();
1053              myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);              const index_t k = mesh->Nodes->borrowGlobalNodesIndex()[id];
1054              myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);              return (myFirstNode <= k && k < myLastNode);
1055              globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          } else {
1056          }              throw ValueError("ownSample: unsupported function space type");
         else if (fs_code == DUDLEY_NODES)  
         {  
             myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);  
             myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);  
             globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);  
         }  
         else  
         {  
             throw DudleyException("unsupported function space type for ownSample()");  
1057          }          }
         k=globalNodeIndex[id];  
         return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );  
 #endif  
1058      }      }
1059    #endif
1060      return true;      return true;
1061  }  }
1062    
# Line 1296  bool MeshAdapter::ownSample(int fs_code, Line 1064  bool MeshAdapter::ownSample(int fs_code,
1064  const_TrilinosGraph_ptr MeshAdapter::getTrilinosGraph() const  const_TrilinosGraph_ptr MeshAdapter::getTrilinosGraph() const
1065  {  {
1066      if (m_graph.is_null()) {      if (m_graph.is_null()) {
1067          m_graph = createTrilinosGraph(m_dudleyMesh.get());          m_graph = m_dudleyMesh->createTrilinosGraph();
1068      }      }
1069      return m_graph;      return m_graph;
1070  }  }
1071  #endif  #endif
1072    
1073  //  //
1074  // creates a stiffness matrix an initializes it with zeros  // creates a stiffness matrix and initializes it with zeros
1075  //  //
1076  escript::ASM_ptr MeshAdapter::newSystemMatrix(int row_blocksize,  escript::ASM_ptr MeshAdapter::newSystemMatrix(int row_blocksize,
1077                              const escript::FunctionSpace& row_functionspace,                              const escript::FunctionSpace& row_functionspace,
# Line 1312  escript::ASM_ptr MeshAdapter::newSystemM Line 1080  escript::ASM_ptr MeshAdapter::newSystemM
1080                              int type) const                              int type) const
1081  {  {
1082      // is the domain right?      // is the domain right?
1083      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));      if (*row_functionspace.getDomain() != *this)
1084      if (row_domain!=*this)          throw ValueError("domain of row function space does not match the domain of matrix generator.");
1085          throw DudleyException("domain of row function space does not match the domain of matrix generator.");      if (*column_functionspace.getDomain() != *this)
     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));  
     if (col_domain!=*this)  
1086          throw DudleyException("domain of column function space does not match the domain of matrix generator.");          throw DudleyException("domain of column function space does not match the domain of matrix generator.");
1087    
1088      int reduceRowOrder=0;      // is the function space type right?
1089      int reduceColOrder=0;      if (row_functionspace.getTypeCode() != DegreesOfFreedom) {
     // is the function space type right  
     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {  
         reduceRowOrder=0;  
     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceRowOrder=1;  
     } else {  
1090          throw DudleyException("illegal function space type for system matrix rows.");          throw DudleyException("illegal function space type for system matrix rows.");
1091      }      }
1092      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      if (column_functionspace.getTypeCode() != DegreesOfFreedom) {
         reduceColOrder=0;  
     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceColOrder=1;  
     } else {  
1093          throw DudleyException("illegal function space type for system matrix columns.");          throw DudleyException("illegal function space type for system matrix columns.");
1094      }      }
1095    
# Line 1350  escript::ASM_ptr MeshAdapter::newSystemM Line 1106  escript::ASM_ptr MeshAdapter::newSystemM
1106                  "used.");                  "used.");
1107  #endif  #endif
1108      } else if (type & (int)SMT_PASO) {      } else if (type & (int)SMT_PASO) {
1109          paso::SystemMatrixPattern_ptr pattern(Dudley_getPattern(          paso::SystemMatrixPattern_ptr pattern(getMesh()->getPasoPattern());
                     getDudley_Mesh(), reduceRowOrder, reduceColOrder));  
1110          paso::SystemMatrix_ptr sm(new paso::SystemMatrix(type, pattern,          paso::SystemMatrix_ptr sm(new paso::SystemMatrix(type, pattern,
1111                    row_blocksize, column_blocksize, false, row_functionspace,                    row_blocksize, column_blocksize, false, row_functionspace,
1112                    column_functionspace));                    column_functionspace));
# Line 1368  escript::ATP_ptr MeshAdapter::newTranspo Line 1123  escript::ATP_ptr MeshAdapter::newTranspo
1123                                               const escript::FunctionSpace& fs,                                               const escript::FunctionSpace& fs,
1124                                               int type) const                                               int type) const
1125  {  {
     int reduceOrder=0;  
1126      // is the domain right?      // is the domain right?
1127      const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(fs.getDomain()));      if (*fs.getDomain() != *this)
     if (domain!=*this)  
1128          throw DudleyException("domain of function space does not match the domain of transport problem generator.");          throw DudleyException("domain of function space does not match the domain of transport problem generator.");
1129      // is the function space type right      // is the function space type right
1130      if (fs.getTypeCode()==DegreesOfFreedom) {      if (fs.getTypeCode() != DegreesOfFreedom) {
         reduceOrder=0;  
     } else if (fs.getTypeCode()==ReducedDegreesOfFreedom) {  
         reduceOrder=1;  
     } else {  
1131          throw DudleyException("illegal function space type for system matrix rows.");          throw DudleyException("illegal function space type for system matrix rows.");
1132      }      }
1133    
1134      // generate matrix      // generate matrix
1135      paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(      paso::SystemMatrixPattern_ptr pattern(getMesh()->getPasoPattern());
                 getDudley_Mesh(),reduceOrder,reduceOrder));  
1136      paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(      paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(
1137                                                fsystemMatrixPattern, blocksize,                                                pattern, blocksize, fs));
                                               fs));  
1138      return transportProblem;      return transportProblem;
1139  }  }
1140    
# Line 1396  escript::ATP_ptr MeshAdapter::newTranspo Line 1143  escript::ATP_ptr MeshAdapter::newTranspo
1143  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1144  {  {
1145      switch (functionSpaceCode) {      switch (functionSpaceCode) {
1146          case(Nodes):          case Nodes:
1147          case(DegreesOfFreedom):          case DegreesOfFreedom:
         case(ReducedDegreesOfFreedom):  
1148              return false;              return false;
1149              break;              break;
1150          case(Elements):          case Elements:
1151          case(FaceElements):          case FaceElements:
1152          case(Points):          case Points:
1153          case(ReducedElements):          case ReducedElements:
1154          case(ReducedFaceElements):          case ReducedFaceElements:
1155              return true;              return true;
1156              break;              break;
1157          default:          default:
1158              stringstream temp;              stringstream ss;
1159              temp << "Cell: Dudley does not know anything about function space type " << functionSpaceCode;              ss << "isCellOriented: Dudley does not know anything about "
1160              throw DudleyException(temp.str());                    "function space type " << functionSpaceCode;
1161              break;              throw ValueError(ss.str());
1162      }      }
1163      return false;      return false;
1164  }  }
# Line 1420  bool MeshAdapter::isCellOriented(int fun Line 1166  bool MeshAdapter::isCellOriented(int fun
1166  bool  bool
1167  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1168  {  {
    /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]  
         class 1: DOF <-> Nodes  
         class 2: ReducedDOF <-> ReducedNodes  
         class 3: Points  
         class 4: Elements  
         class 5: ReducedElements  
         class 6: FaceElements  
         class 7: ReducedFaceElements  
         class 8: ContactElementZero <-> ContactElementOne  
         class 9: ReducedContactElementZero <-> ReducedContactElementOne  
   
    There is also a set of lines. Interpolation is possible down a line but not between lines.  
    class 1 and 2 belong to all lines so aren't considered.  
         line 0: class 3  
         line 1: class 4,5  
         line 2: class 6,7  
         line 3: class 8,9  
   
    For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.  
    eg hasnodes is true if we have at least one instance of Nodes.  
    */  
1169      if (fs.empty())      if (fs.empty())
     {  
1170          return false;          return false;
1171      }      // The idea is to use equivalence classes, i.e. types which can be
1172      vector<int> hasclass(10);      // interpolated back and forth
1173      vector<int> hasline(4);          //    class 1: DOF <-> Nodes
1174      bool hasnodes=false;      //    class 3: Points
1175      bool hasrednodes=false;      //    class 4: Elements
1176      for (int i=0;i<fs.size();++i)      //    class 5: ReducedElements
1177      {      //    class 6: FaceElements
1178          switch(fs[i]) {      //    class 7: ReducedFaceElements
1179              case(Nodes):  
1180                  hasnodes=true;   // no break is deliberate      // There is also a set of lines. Interpolation is possible down a line but
1181              case(DegreesOfFreedom):      // not between lines.
1182                  hasclass[1]=1;      // class 1 and 2 belong to all lines so aren't considered.
1183        //    line 0: class 3
1184        //    line 1: class 4,5
1185        //    line 2: class 6,7
1186    
1187        // For classes with multiple members (class 1) we have vars to record
1188        // if there is at least one instance -> hasnodes is true if we have at
1189        // least one instance of Nodes.
1190        vector<int> hasclass(8);
1191        vector<int> hasline(3);
1192        bool hasnodes = false;
1193        for (int i = 0; i < fs.size(); ++i) {
1194            switch (fs[i]) {
1195                case Nodes:
1196                    hasnodes = true; // fall through
1197                case DegreesOfFreedom:
1198                    hasclass[1] = 1;
1199                  break;                  break;
1200              case(ReducedNodes):              case Points:
1201                  hasrednodes=true; // no break is deliberate                  hasline[0] = 1;
1202              case(ReducedDegreesOfFreedom):                  hasclass[3] = 1;
                 hasclass[2]=1;  
1203                  break;                  break;
1204              case(Points):              case Elements:
1205                  hasline[0]=1;                  hasclass[4] = 1;
1206                  hasclass[3]=1;                  hasline[1] = 1;
1207                  break;                  break;
1208              case(Elements):              case ReducedElements:
1209                  hasclass[4]=1;                  hasclass[5] = 1;
1210                  hasline[1]=1;                  hasline[1] = 1;
1211                  break;                  break;
1212              case(ReducedElements):              case FaceElements:
1213                  hasclass[5]=1;                  hasclass[6] = 1;
1214                  hasline[1]=1;                  hasline[2] = 1;
1215                  break;                  break;
1216              case(FaceElements):              case ReducedFaceElements:
1217                  hasclass[6]=1;                  hasclass[7] = 1;
1218                  hasline[2]=1;                  hasline[2] = 1;
                 break;  
             case(ReducedFaceElements):  
                 hasclass[7]=1;  
                 hasline[2]=1;  
1219                  break;                  break;
1220              default:              default:
1221                  return false;                  return false;
1222          }          }
1223      }      }
1224      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int totlines = hasline[0]+hasline[1]+hasline[2];
1225      // fail if we have more than one leaf group      // fail if we have more than one leaf group
1226        if (totlines > 1)
1227            // there are at least two branches we can't interpolate between
1228            return false;
1229    
1230      if (totlines>1)      if (totlines == 1) {
1231      {          if (hasline[0] == 1) // we have points
1232          return false;   // there are at least two branches we can't interpolate between              resultcode = Points;
1233      }          else if (hasline[1] == 1) {
     else if (totlines==1)  
     {  
         if (hasline[0]==1)              // we have points  
         {  
             resultcode=Points;  
         }  
         else if (hasline[1]==1)  
         {  
1234              if (hasclass[5]==1)              if (hasclass[5]==1)
             {  
1235                  resultcode=ReducedElements;                  resultcode=ReducedElements;
             }  
1236              else              else
             {  
1237                  resultcode=Elements;                  resultcode=Elements;
1238              }          } else if (hasline[2]==1) {
         }  
         else if (hasline[2]==1)  
         {  
1239              if (hasclass[7]==1)              if (hasclass[7]==1)
             {  
1240                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
             }  
1241              else              else
             {  
1242                  resultcode=FaceElements;                  resultcode=FaceElements;
             }  
         }  
         else    // so we must be in line3  
         {  
   
             throw DudleyException("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);  
1243          }          }
1244        } else { // totlines==0
1245            // something from class 1
1246            resultcode = (hasnodes ? Nodes : DegreesOfFreedom);
1247      }      }
1248      return true;      return true;
1249  }  }
1250    
1251  signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,
1252                                                 int functionSpaceType_target) const
1253    {
1254        switch(functionSpaceType_source) {
1255            case Nodes:
1256                switch (functionSpaceType_target) {
1257                    case Nodes:
1258                    case DegreesOfFreedom:
1259                    case Elements:
1260                    case ReducedElements:
1261                    case FaceElements:
1262                    case ReducedFaceElements:
1263                    case Points:
1264                        return true;
1265                    default:
1266                        stringstream ss;
1267                        ss << "Interpolation On Domain: Dudley does not know "
1268                            "anything about function space type "
1269                            << functionSpaceType_target;
1270                        throw ValueError(ss.str());
1271                }
1272            break;
1273            case Elements:
1274                return (functionSpaceType_target == Elements ||
1275                        functionSpaceType_target == ReducedElements);
1276            case ReducedElements:
1277                return (functionSpaceType_target == ReducedElements);
1278            case FaceElements:
1279                return (functionSpaceType_target == FaceElements ||
1280                        functionSpaceType_target == ReducedFaceElements);
1281            case ReducedFaceElements:
1282                return (functionSpaceType_target == ReducedFaceElements);
1283            case Points:
1284                return (functionSpaceType_target == Points);
1285            case DegreesOfFreedom:
1286                switch (functionSpaceType_target) {
1287                    case DegreesOfFreedom:
1288                    case Nodes:
1289                    case Elements:
1290                    case ReducedElements:
1291                    case Points:
1292                    case FaceElements:
1293                    case ReducedFaceElements:
1294                        return true;
1295                    default:
1296                        stringstream ss;
1297                        ss << "Interpolation On Domain: Dudley does not know "
1298                              "anything about function space type "
1299                           << functionSpaceType_target;
1300                        throw DudleyException(ss.str());
1301                }
1302                break;
1303            default:
1304                stringstream ss;
1305                ss << "Interpolation On Domain: Dudley does not know anything "
1306                      "about function space type " << functionSpaceType_source;
1307                throw DudleyException(ss.str());
1308        }
1309        return false;
1310    }
1311    
1312    signed char MeshAdapter::preferredInterpolationOnDomain(
1313            int functionSpaceType_source,int functionSpaceType_target) const
1314  {  {
1315      if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))      if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1316          return 1;          return 1;
# Line 1552  signed char MeshAdapter::preferredInterp Line 1319  signed char MeshAdapter::preferredInterp
1319      return 0;      return 0;
1320  }  }
1321    
   
   
 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  
 {  
    switch(functionSpaceType_source) {  
    case(Nodes):  
         switch(functionSpaceType_target) {  
         case(Nodes):  
         case(ReducedNodes):  
         case(ReducedDegreesOfFreedom):  
         case(DegreesOfFreedom):  
         case(Elements):  
         case(ReducedElements):  
         case(FaceElements):  
         case(ReducedFaceElements):  
         case(Points):  
         return true;  
         default:  
               stringstream temp;  
               temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;  
               throw DudleyException(temp.str());  
    }  
    break;  
    case(ReducedNodes):  
         switch(functionSpaceType_target) {  
         case(ReducedNodes):  
         case(ReducedDegreesOfFreedom):  
         case(Elements):  
         case(ReducedElements):  
         case(FaceElements):  
         case(ReducedFaceElements):  
         case(Points):  
         return true;  
         case(Nodes):  
         case(DegreesOfFreedom):  
         return false;  
         default:  
                 stringstream temp;  
                 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;  
                 throw DudleyException(temp.str());  
    }  
    break;  
    case(Elements):  
         if (functionSpaceType_target==Elements) {  
           return true;  
         } else if (functionSpaceType_target==ReducedElements) {  
           return true;  
         } else {  
           return false;  
         }  
    case(ReducedElements):  
         if (functionSpaceType_target==ReducedElements) {  
           return true;  
         } else {  
           return false;  
         }  
    case(FaceElements):  
         if (functionSpaceType_target==FaceElements) {  
                 return true;  
         } else if (functionSpaceType_target==ReducedFaceElements) {  
                 return true;  
         } else {  
                 return false;  
         }  
    case(ReducedFaceElements):  
         if (functionSpaceType_target==ReducedFaceElements) {  
                 return true;  
         } else {  
                 return false;  
         }  
    case(Points):  
         if (functionSpaceType_target==Points) {  
                 return true;  
         } else {  
                 return false;  
         }  
    case(DegreesOfFreedom):  
         switch(functionSpaceType_target) {  
         case(ReducedDegreesOfFreedom):  
         case(DegreesOfFreedom):  
         case(Nodes):  
         case(ReducedNodes):  
         case(Elements):  
         case(ReducedElements):  
         case(Points):  
         case(FaceElements):  
         case(ReducedFaceElements):  
         return true;  
         default:  
                 stringstream temp;  
                 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;  
                 throw DudleyException(temp.str());  
         }  
         break;  
    case(ReducedDegreesOfFreedom):  
    switch(functionSpaceType_target) {  
         case(ReducedDegreesOfFreedom):  
         case(ReducedNodes):  
         case(Elements):  
         case(ReducedElements):  
         case(FaceElements):  
         case(ReducedFaceElements):  
         case(Points):  
         return true;  
         case(Nodes):  
         case(DegreesOfFreedom):  
         return false;  
         default:  
                 stringstream temp;  
                 temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;  
                 throw DudleyException(temp.str());  
         }  
         break;  
    default:  
       stringstream temp;  
       temp << "Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;  
       throw DudleyException(temp.str());  
       break;  
    }  
    return false;  
 }  
   
1322  bool MeshAdapter::probeInterpolationAcross(int functionSpaceType_source,  bool MeshAdapter::probeInterpolationAcross(int functionSpaceType_source,
1323          const AbstractDomain& targetDomain, int functionSpaceType_target) const          const AbstractDomain& targetDomain, int functionSpaceType_target) const
1324  {  {
# Line 1682  bool MeshAdapter::probeInterpolationAcro Line 1327  bool MeshAdapter::probeInterpolationAcro
1327    
1328  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
1329  {  {
1330      const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);      const MeshAdapter* temp = dynamic_cast<const MeshAdapter*>(&other);
1331      if (temp!=0) {      if (temp) {
1332          return (m_dudleyMesh==temp->m_dudleyMesh);          return (m_dudleyMesh == temp->m_dudleyMesh);
     } else {  
         return false;  
1333      }      }
1334        return false;
1335  }  }
1336    
1337  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
# Line 1704  int MeshAdapter::getSystemMatrixTypeId(c Line 1348  int MeshAdapter::getSystemMatrixTypeId(c
1348  #ifdef USE_TRILINOS  #ifdef USE_TRILINOS
1349          return (int)SMT_TRILINOS;          return (int)SMT_TRILINOS;
1350  #else  #else
1351          throw DudleyAdapterException("Trilinos requested but not built with Trilinos.");                throw DudleyException("Trilinos requested but not built with Trilinos.");      
1352  #endif  #endif
1353      }      }
1354      return (int)SMT_PASO | paso::SystemMatrix::getSystemMatrixTypeId(      return (int)SMT_PASO | paso::SystemMatrix::getSystemMatrixTypeId(
# Line 1715  int MeshAdapter::getSystemMatrixTypeId(c Line 1359  int MeshAdapter::getSystemMatrixTypeId(c
1359  int MeshAdapter::getTransportTypeId(int solver, int preconditioner,  int MeshAdapter::getTransportTypeId(int solver, int preconditioner,
1360                                      int package, bool symmetry) const                                      int package, bool symmetry) const
1361  {  {
     Dudley_Mesh* mesh=m_dudleyMesh.get();  
1362      return paso::TransportProblem::getTypeId(solver, preconditioner, package,      return paso::TransportProblem::getTypeId(solver, preconditioner, package,
1363                                               symmetry, mesh->MPIInfo);                                               symmetry, getMPI());
1364  }  }
1365    
1366  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
# Line 1735  escript::Data MeshAdapter::getSize() con Line 1378  escript::Data MeshAdapter::getSize() con
1378      return escript::function(*this).getSize();      return escript::function(*this).getSize();
1379  }  }
1380    
1381  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const index_t* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1382  {  {
1383     int *out = NULL;      index_t* out = NULL;
1384     Dudley_Mesh* mesh=m_dudleyMesh.get();      switch (functionSpaceType) {
1385     switch (functionSpaceType) {          case Nodes:
1386     case(Nodes):              out = getMesh()->Nodes->Id;
1387     out=mesh->Nodes->Id;          break;
1388     break;          case Elements:
1389     case(ReducedNodes):              out = getMesh()->Elements->Id;
1390     out=mesh->Nodes->reducedNodesId;          break;
1391     break;          case ReducedElements:
1392     case(Elements):              out = getMesh()->Elements->Id;
1393     out=mesh->Elements->Id;          break;
1394     break;          case FaceElements:
1395     case(ReducedElements):              out = getMesh()->FaceElements->Id;
1396     out=mesh->Elements->Id;          break;
1397     break;          case ReducedFaceElements:
1398     case(FaceElements):              out = getMesh()->FaceElements->Id;
1399     out=mesh->FaceElements->Id;          break;
1400     break;          case Points:
1401     case(ReducedFaceElements):              out = getMesh()->Points->Id;
1402     out=mesh->FaceElements->Id;          break;
1403     break;          case DegreesOfFreedom:
1404     case(Points):              out = getMesh()->Nodes->degreesOfFreedomId;
1405     out=mesh->Points->Id;          break;
1406     break;          default:
1407     case(DegreesOfFreedom):              stringstream ss;
1408     out=mesh->Nodes->degreesOfFreedomId;              ss << "Invalid function space type: " << functionSpaceType
1409     break;                 << " for domain: " << getDescription();
1410     case(ReducedDegreesOfFreedom):              throw ValueError(ss.str());
1411     out=mesh->Nodes->reducedDegreesOfFreedomId;      }
1412     break;      return out;
1413     default:  }
1414        stringstream temp;  
1415        temp << "Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, index_t sampleNo) const
1416        throw DudleyException(temp.str());  {
1417        break;      int out = 0;
1418     }      switch (functionSpaceType) {
1419     return out;          case Nodes:
1420  }              out = getMesh()->Nodes->Tag[sampleNo];
1421  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const          break;
1422  {          case Elements:
1423     int out=0;              out = getMesh()->Elements->Tag[sampleNo];
1424     Dudley_Mesh* mesh=m_dudleyMesh.get();          break;
1425     switch (functionSpaceType) {          case ReducedElements:
1426     case(Nodes):              out = getMesh()->Elements->Tag[sampleNo];
1427     out=mesh->Nodes->Tag[sampleNo];          break;
1428     break;          case FaceElements:
1429     case(ReducedNodes):              out = getMesh()->FaceElements->Tag[sampleNo];
1430     throw DudleyException("ReducedNodes does not support tags.");          break;
1431     break;          case ReducedFaceElements:
1432     case(Elements):              out = getMesh()->FaceElements->Tag[sampleNo];
1433     out=mesh->Elements->Tag[sampleNo];          break;
1434     break;          case Points:
1435     case(ReducedElements):              out = getMesh()->Points->Tag[sampleNo];
1436     out=mesh->Elements->Tag[sampleNo];          break;
1437     break;          case DegreesOfFreedom:
1438     case(FaceElements):              throw DudleyException("DegreesOfFreedom does not support tags.");
1439     out=mesh->FaceElements->Tag[sampleNo];          break;
1440     break;          default:
1441     case(ReducedFaceElements):              stringstream ss;
1442     out=mesh->FaceElements->Tag[sampleNo];              ss << "Invalid function space type: " << functionSpaceType
1443     break;                 << " for domain: " << getDescription();
1444     case(Points):              throw DudleyException(ss.str());
1445     out=mesh->Points->Tag[sampleNo];      }
1446     break;      return out;
    case(DegreesOfFreedom):  
    throw DudleyException("DegreesOfFreedom does not support tags.");  
    break;  
    case(ReducedDegreesOfFreedom):  
    throw DudleyException("ReducedDegreesOfFreedom does not support tags.");  
    break;  
    default:  
       stringstream temp;  
       temp << "Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
       throw DudleyException(temp.str());  
       break;  
    }  
    return out;  
1447  }  }
1448    
1449    
1450  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(int functionSpaceType, int newTag, const escript::Data& mask) const
1451  {  {
     Dudley_Mesh* mesh=m_dudleyMesh.get();  
1452      switch (functionSpaceType) {      switch (functionSpaceType) {
1453          case Nodes:          case Nodes:
1454              Dudley_NodeFile_setTags(mesh->Nodes,newTag,&mask);              getMesh()->Nodes->setTags(newTag, mask);
             break;  
         case ReducedNodes:  
             throw DudleyException("ReducedNodes does not support tags");  
1455              break;              break;
1456          case DegreesOfFreedom:          case DegreesOfFreedom:
1457              throw DudleyException("DegreesOfFreedom does not support tags");              throw DudleyException("DegreesOfFreedom does not support tags");
1458              break;              break;
1459          case ReducedDegreesOfFreedom:          case Elements: // fall through
             throw DudleyException("ReducedDegreesOfFreedom does not support tags");  
             break;  
         case Elements:  
             Dudley_ElementFile_setTags(mesh->Elements,newTag,&mask);  
             break;  
1460          case ReducedElements:          case ReducedElements:
1461              Dudley_ElementFile_setTags(mesh->Elements,newTag,&mask);              getMesh()->Elements->setTags(newTag, mask);
1462              break;              break;
1463          case FaceElements:          case FaceElements:
             Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&mask);  
             break;  
1464          case ReducedFaceElements:          case ReducedFaceElements:
1465              Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&mask);              getMesh()->FaceElements->setTags(newTag, mask);
1466              break;              break;
1467          case Points:          case Points:
1468              Dudley_ElementFile_setTags(mesh->Points,newTag,&mask);              getMesh()->Points->setTags(newTag, mask);
1469              break;              break;
1470          default:          default:
1471              stringstream temp;              stringstream ss;
1472              temp << "Dudley does not know anything about function space type " << functionSpaceType;              ss << "Dudley does not know anything about function space type "
1473              throw DudleyException(temp.str());                 << functionSpaceType;
1474                throw ValueError(ss.str());
1475      }      }
1476  }  }
1477    
1478  void MeshAdapter::setTagMap(const string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
1479  {  {
1480      Dudley_Mesh* mesh=m_dudleyMesh.get();      getMesh()->addTagMap(name, tag);
     Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);  
1481  }  }
1482    
1483  int MeshAdapter::getTag(const string& name) const  int MeshAdapter::getTag(const string& name) const
1484  {  {
1485      Dudley_Mesh* mesh=m_dudleyMesh.get();      return getMesh()->getTag(name);
     int tag=0;  
     tag=Dudley_Mesh_getTag(mesh, name.c_str());  
     return tag;  
1486  }  }
1487    
1488  bool MeshAdapter::isValidTagName(const string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
1489  {  {
1490      Dudley_Mesh* mesh=m_dudleyMesh.get();      return getMesh()->isValidTagName(name);
     return Dudley_Mesh_isValidTagName(mesh,name.c_str());  
1491  }  }
1492    
1493  string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
1494  {  {
1495      stringstream temp;      stringstream ss;
1496      Dudley_Mesh* mesh=m_dudleyMesh.get();      TagMap::const_iterator it = getMesh()->tagMap.begin();
1497      Dudley_TagMap* tag_map=mesh->TagMap;      while (it != getMesh()->tagMap.end()) {
1498      while (tag_map) {          ss << it->first;
1499          temp << tag_map->name;          ++it;
1500          tag_map=tag_map->next;          if (it != getMesh()->tagMap.end())
1501          if (tag_map) temp << ", ";              ss << ", ";
1502      }      }
1503      return temp.str();      return ss.str();
1504  }  }
1505    
1506  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1507  {  {
1508    Dudley_Mesh* mesh=m_dudleyMesh.get();      switch (functionSpaceCode) {
1509    dim_t numTags=0;          case Nodes:
1510    switch(functionSpaceCode) {              return getMesh()->Nodes->tagsInUse.size();
1511     case(Nodes):          case DegreesOfFreedom:
1512            numTags=mesh->Nodes->numTagsInUse;              throw ValueError("DegreesOfFreedom does not support tags");
1513            break;          case Elements: // fall through
1514     case(ReducedNodes):          case ReducedElements:
1515            throw DudleyException("ReducedNodes does not support tags");              return getMesh()->Elements->tagsInUse.size();
1516            break;          case FaceElements: // fall through
1517     case(DegreesOfFreedom):          case ReducedFaceElements:
1518            throw DudleyException("DegreesOfFreedom does not support tags");              return getMesh()->FaceElements->tagsInUse.size();
1519            break;          case Points:
1520     case(ReducedDegreesOfFreedom):              return getMesh()->Points->tagsInUse.size();
1521            throw DudleyException("ReducedDegreesOfFreedom does not support tags");          default:
1522            break;              stringstream ss;
1523     case(Elements):              ss << "Dudley does not know anything about function space type "
1524     case(ReducedElements):                 << functionSpaceCode;
1525            numTags=mesh->Elements->numTagsInUse;              throw ValueError(ss.str());
1526            break;      }
1527     case(FaceElements):      return 0;
    case(ReducedFaceElements):  
           numTags=mesh->FaceElements->numTagsInUse;  
           break;  
    case(Points):  
           numTags=mesh->Points->numTagsInUse;  
           break;  
    default:  
       stringstream temp;  
       temp << "Dudley does not know anything about function space type " << functionSpaceCode;  
       throw DudleyException(temp.str());  
   }  
   return numTags;  
1528  }  }
1529    
1530  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
1531  {  {
     Dudley_Mesh* mesh=m_dudleyMesh.get();  
     index_t* tags=NULL;  
1532      switch (functionSpaceCode) {      switch (functionSpaceCode) {
1533          case Nodes:          case Nodes:
1534              tags=mesh->Nodes->tagsInUse;              if (getMesh()->Nodes->tagsInUse.empty())
1535              break;                  return NULL;
1536          case ReducedNodes:              else
1537              throw DudleyException("ReducedNodes does not support tags");                  return &getMesh()->Nodes->tagsInUse[0];
             break;  
1538          case DegreesOfFreedom:          case DegreesOfFreedom:
1539              throw DudleyException("DegreesOfFreedom does not support tags");              throw DudleyException("DegreesOfFreedom does not support tags");
1540              break;          case Elements: // fall through
         case ReducedDegreesOfFreedom:  
             throw DudleyException("ReducedDegreesOfFreedom does not support tags");  
             break;  
         case Elements:  
1541          case ReducedElements:          case ReducedElements:
1542              tags=mesh->Elements->tagsInUse;              if (getMesh()->Elements->tagsInUse.empty())
1543              break;                  return NULL;
1544          case FaceElements:              else
1545                    return &getMesh()->Elements->tagsInUse[0];
1546            case FaceElements: // fall through
1547          case ReducedFaceElements:          case ReducedFaceElements:
1548              tags=mesh->FaceElements->tagsInUse;              if (getMesh()->FaceElements->tagsInUse.empty())
1549              break;                  return NULL;
1550                else
1551                    return &getMesh()->FaceElements->tagsInUse[0];
1552          case Points:          case Points:
1553              tags=mesh->Points->tagsInUse;              if (getMesh()->Points->tagsInUse.empty())
1554              break;                  return NULL;
1555                else
1556                    return &getMesh()->Points->tagsInUse[0];
1557          default:          default:
1558              stringstream temp;              stringstream ss;
1559              temp << "Dudley does not know anything about function space type " << functionSpaceCode;              ss << "Dudley does not know anything about function space type "
1560              throw DudleyException(temp.str());                 << functionSpaceCode;
1561                throw DudleyException(ss.str());
1562      }      }
1563      return tags;      return NULL;
1564  }  }
1565    
1566    
# Line 1978  bool MeshAdapter::canTag(int functionSpa Line 1581  bool MeshAdapter::canTag(int functionSpa
1581    
1582  MeshAdapter::StatusType MeshAdapter::getStatus() const  MeshAdapter::StatusType MeshAdapter::getStatus() const
1583  {  {
1584      Dudley_Mesh* mesh=m_dudleyMesh.get();      return getMesh()->getStatus();
     return Dudley_Mesh_getStatus(mesh);  
1585  }  }
1586    
1587  int MeshAdapter::getApproximationOrder(int functionSpaceCode) const  int MeshAdapter::getApproximationOrder(int functionSpaceCode) const
1588  {  {
1589          switch (functionSpaceCode) {
1590    Dudley_Mesh* mesh=m_dudleyMesh.get();          case Nodes:
1591    int order =-1;          case DegreesOfFreedom:
1592    switch(functionSpaceCode) {              return getMesh()->approximationOrder;
1593     case(Nodes):          case Elements:
1594     case(DegreesOfFreedom):          case FaceElements:
1595            order=mesh->approximationOrder;          case Points:
1596            break;              return getMesh()->integrationOrder;
1597     case(ReducedNodes):          case ReducedElements:
1598     case(ReducedDegreesOfFreedom):          case ReducedFaceElements:
1599            order=mesh->reducedApproximationOrder;              return getMesh()->reducedIntegrationOrder;
1600            break;          default:
1601     case(Elements):              stringstream ss;
1602     case(FaceElements):              ss << "Dudley does not know anything about function space type "
1603     case(Points):                 << functionSpaceCode;
1604            order=mesh->integrationOrder;              throw ValueError(ss.str());
1605            break;      }
1606     case(ReducedElements):      return 0;
    case(ReducedFaceElements):  
           order=mesh->reducedIntegrationOrder;  
           break;  
    default:  
       stringstream temp;  
       temp << "Dudley does not know anything about function space type " << functionSpaceCode;  
       throw DudleyException(temp.str());  
   }  
   return order;  
1607  }  }
1608    
1609  bool MeshAdapter::supportsContactElements() const  bool MeshAdapter::supportsContactElements() const
# Line 2018  bool MeshAdapter::supportsContactElement Line 1611  bool MeshAdapter::supportsContactElement
1611      return false;      return false;
1612  }  }
1613    
1614  escript::Data MeshAdapter::randomFill(const escript::DataTypes::ShapeType& shape,  escript::Data MeshAdapter::randomFill(
1615                                  const escript::FunctionSpace& what, long seed,          const escript::DataTypes::ShapeType& shape,
1616                                  const bp::tuple& filter) const          const escript::FunctionSpace& what, long seed,
1617            const bp::tuple& filter) const
1618  {  {
1619      escript::Data towipe(0, shape, what, true);      escript::Data towipe(0, shape, what, true);
1620      // since we just made this object, no sharing is possible and we don't      // since we just made this object, no sharing is possible and we don't
1621      // need to check for exlusive write      // need to check for exclusive write
1622      escript::DataTypes::RealVectorType& dv=towipe.getExpandedVectorReference();      escript::DataTypes::RealVectorType& dv(towipe.getExpandedVectorReference());
1623      const size_t dvsize=dv.size();      escript::randomFillArray(seed, &dv[0], dv.size());
     escript::randomFillArray(seed, &(dv[0]), dvsize);  
1624      return towipe;            return towipe;      
1625  }  }
1626    

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26