/[escript]/trunk/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

revision 2315 by gross, Wed Mar 18 00:38:48 2009 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21  #ifdef PASO_MPI  #ifdef ESYS_MPI
22  #include <mpi.h>  #include <mpi.h>
23  #include "paso/Paso_MPI.h"  #include "esysUtils/Esys_MPI.h"
24  #endif  #endif
25  extern "C" {  extern "C" {
26  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
# Line 85  int MeshAdapter::getMPIRank() const Line 85  int MeshAdapter::getMPIRank() const
85  }  }
86  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
87  {  {
88  #ifdef PASO_MPI  #ifdef ESYS_MPI
89     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
90  #endif  #endif
91     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 96  bool MeshAdapter::onMasterProcessor() co
96  }  }
97    
98    
99    #ifdef ESYS_MPI
100      MPI_Comm
101    #else
102      unsigned int
103    #endif
104    MeshAdapter::getMPIComm() const
105    {
106    #ifdef ESYS_MPI
107        return m_finleyMesh->MPIInfo->comm;
108    #else
109        return 0;
110    #endif
111    }
112    
113    
114  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
115     return m_finleyMesh.get();     return m_finleyMesh.get();
116  }  }
117    
118  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
119  {  {
120     char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;     char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
121     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 109  void MeshAdapter::write(const std::strin Line 124  void MeshAdapter::write(const std::strin
124     TMPMEMFREE(fName);     TMPMEMFREE(fName);
125  }  }
126    
127  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
128  {  {
129     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
130  }  }
131    
132  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
133  {  {
134  #ifdef USE_NETCDF  #ifdef USE_NETCDF
135     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};
# Line 134  void MeshAdapter::dump(const std::string Line 149  void MeshAdapter::dump(const std::string
149     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
150     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
151     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
152  #ifdef PASO_MPI  #ifdef ESYS_MPI
153     MPI_Status status;     MPI_Status status;
154  #endif  #endif
155    
156  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
157  #ifdef PASO_MPI  #ifdef ESYS_MPI
158     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
159  #endif  #endif
160    
161     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
162                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
163    
164     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
165     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
166     num_Tags = 0;     num_Tags = 0;
167     if (tag_map) {     while (tag_map) {
168        while (tag_map) {        num_Tags++;
169           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
170     }     }
171    
172     // NetCDF error handler     // NetCDF error handler
173     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
174     // Create the file.     // Create the file.
175     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
176       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
177     // check if writing was successful     // check if writing was successful
178     if (!dataFile.is_valid())     if (!dataFile.is_valid())
179        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
180    
181     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
182       // dim_Elements only appears if > 0)
183     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
184        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
185     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
186        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
187     if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )     if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
188        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
189     if (num_Elements>0)     if (num_Elements>0)
190        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
191           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
192     if (num_FaceElements>0)     if (num_FaceElements>0)
193        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
194           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
195     if (num_ContactElements>0)     if (num_ContactElements>0)
196        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
197           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
198     if (num_Points>0)     if (num_Points>0)
199        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
200           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
201     if (num_Elements>0)     if (num_Elements>0)
202        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
203           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
204     if (num_FaceElements>0)     if (num_FaceElements>0)
205        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
206           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
207     if (num_ContactElements>0)     if (num_ContactElements>0)
208        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
209           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
210     if (num_Tags>0)     if (num_Tags>0)
211        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
212           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
213    
214     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
215     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
216        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
217     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
218        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
219     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
220        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
221     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
222        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
223     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
224        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
225     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
226        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
227     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
228        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
229     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
230        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
231     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
232        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
233     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
234        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
235     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
236        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
237     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
238        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
239     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
240        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
241     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
242        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
243     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
245     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
247     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
249     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
250        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
251     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
252        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
253    
254     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
255    
256     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
257       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
258          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
259       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
260       if (! (ids->put(int_ptr, mpi_size+1)) )
261          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
262    
263       // Nodes degreesOfFreedomDistribution
264       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
265          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
266       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
267       if (! (ids->put(int_ptr, mpi_size+1)) )
268          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
269    
270       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
271       // (it treats them as NC_UNLIMITED)
272     if (numNodes>0) {     if (numNodes>0) {
273    
274        // Nodes Id        // Nodes Id
275        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
276           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
277        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
278        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
279           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
280    
281        // Nodes Tag        // Nodes Tag
282        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
283           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
284        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
285        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
286           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
287    
288        // Nodes gDOF        // Nodes gDOF
289        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
290           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
291        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
292        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
293           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
294    
295        // Nodes global node index        // Nodes global node index
296        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
297           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
298        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
299        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
300           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
301    
302        // Nodes grDof        // Nodes grDof
303        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
304           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
305        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
306        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
307           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
308    
309        // Nodes grNI        // Nodes grNI
310        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
311           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
312        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
313        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
314           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
315    
316        // Nodes Coordinates        // Nodes Coordinates
317        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
318           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
319        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
320           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Coordinates)");
   
       // Nodes degreesOfFreedomDistribution  
       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);  
   
       // Nodes nodeDistribution  
       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);  
321    
322     }     }
323    
# Line 311  void MeshAdapter::dump(const std::string Line 327  void MeshAdapter::dump(const std::string
327    
328        // Elements_Id        // Elements_Id
329        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
330           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
331        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
332        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
333           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
334    
335        // Elements_Tag        // Elements_Tag
336        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
337           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
338        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
339        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
340           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
341    
342        // Elements_Owner        // Elements_Owner
343        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
344           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
345        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
346        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
347           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
348    
349        // Elements_Color        // Elements_Color
350        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
351           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
352        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
353        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
354           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
355    
356        // Elements_Nodes        // Elements_Nodes
357        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
358           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
359        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
360           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
361    
362     }     }
363    
# Line 351  void MeshAdapter::dump(const std::string Line 367  void MeshAdapter::dump(const std::string
367    
368        // FaceElements_Id        // FaceElements_Id
369        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
370           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
371        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
372        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
373           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
374    
375        // FaceElements_Tag        // FaceElements_Tag
376        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
377           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
378        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
379        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
380           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
381    
382        // FaceElements_Owner        // FaceElements_Owner
383        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
384           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
385        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
386        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
387           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
388    
389        // FaceElements_Color        // FaceElements_Color
390        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
391           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
392        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
393        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
394           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
395    
396        // FaceElements_Nodes        // FaceElements_Nodes
397        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
398           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
399        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
400           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
401    
402     }     }
403    
# Line 391  void MeshAdapter::dump(const std::string Line 407  void MeshAdapter::dump(const std::string
407    
408        // ContactElements_Id        // ContactElements_Id
409        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
410           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
411        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
412        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
413           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
414    
415        // ContactElements_Tag        // ContactElements_Tag
416        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
417           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
418        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
419        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
420           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
421    
422        // ContactElements_Owner        // ContactElements_Owner
423        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
424           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
425        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
426        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
427           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
428    
429        // ContactElements_Color        // ContactElements_Color
430        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
431           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
432        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
433        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
434           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
435    
436        // ContactElements_Nodes        // ContactElements_Nodes
437        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
438           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
439        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
440           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
441    
442     }     }
443    
# Line 433  void MeshAdapter::dump(const std::string Line 449  void MeshAdapter::dump(const std::string
449    
450        // Points_Id        // Points_Id
451        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
452           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
453        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
454        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
455           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
456    
457        // Points_Tag        // Points_Tag
458        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
459           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
460        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
461        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
462           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
463    
464        // Points_Owner        // Points_Owner
465        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
466           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
467        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
468        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
469           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
470    
471        // Points_Color        // Points_Color
472        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
473           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
474        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
475        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
476           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
477    
478        // Points_Nodes        // Points_Nodes
479        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
480        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
481           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
482        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
483           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
484    
485     }     }
486    
# Line 488  void MeshAdapter::dump(const std::string Line 504  void MeshAdapter::dump(const std::string
504    
505        // Tags_keys        // Tags_keys
506        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
507           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
508        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
509        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
510           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
511    
512        // Tags_names_*        // Tags_names_*
513        // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string        // This is an array of strings, it should be stored as an array but
514        // because the NetCDF manual doesn't tell how to do an array of strings        // instead I have hacked in one attribute per string because the NetCDF
515          // manual doesn't tell how to do an array of strings
516        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
517        if (tag_map) {        if (tag_map) {
518           int i = 0;           int i = 0;
519           while (tag_map) {           while (tag_map) {
520              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
521              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
522                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
523              tag_map=tag_map->next;              tag_map=tag_map->next;
524              i++;              i++;
525           }           }
526        }        }
527    
528        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
529     }     }
530    
531  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
532  #ifdef PASO_MPI  #ifdef ESYS_MPI
533     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) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
534  #endif  #endif
535    
# Line 678  pair<int,int> MeshAdapter::getDataShape( Line 694  pair<int,int> MeshAdapter::getDataShape(
694     case(Elements):     case(Elements):
695     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
696        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
697        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
698     }     }
699     break;     break;
700     case(ReducedElements):     case(ReducedElements):
701     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
702        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
703        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
704     }     }
705     break;     break;
706     case(FaceElements):     case(FaceElements):
707     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
708        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
709        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
710     }     }
711     break;     break;
712     case(ReducedFaceElements):     case(ReducedFaceElements):
713     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
714        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
715        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
716     }     }
717     break;     break;
# Line 707  pair<int,int> MeshAdapter::getDataShape( Line 723  pair<int,int> MeshAdapter::getDataShape(
723     break;     break;
724     case(ContactElementsZero):     case(ContactElementsZero):
725     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
726        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
727        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
728     }     }
729     break;     break;
730     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
731     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
732        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
733        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
734     }     }
735     break;     break;
736     case(ContactElementsOne):     case(ContactElementsOne):
737     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
738        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
739        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
740     }     }
741     break;     break;
742     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
743     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
744        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
745        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
746     }     }
747     break;     break;
# Line 795  void  MeshAdapter::addPDEToLumpedSystem( Line 811  void  MeshAdapter::addPDEToLumpedSystem(
811     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
812    
813     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
814       checkFinleyError();
815      
816     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
   
817     checkFinleyError();     checkFinleyError();
818    
819  }  }
820    
821    
# Line 849  void MeshAdapter::addPDEToTransportProbl Line 867  void MeshAdapter::addPDEToTransportProbl
867     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
868    
869     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
870     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();
871    
872     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
873     checkFinleyError();     checkFinleyError();
# Line 881  void MeshAdapter::interpolateOnDomain(es Line 899  void MeshAdapter::interpolateOnDomain(es
899     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
900     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
901     case(Nodes):     case(Nodes):
902     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
903     case(Nodes):        case(Nodes):
904     case(ReducedNodes):        case(ReducedNodes):
905     case(DegreesOfFreedom):        case(DegreesOfFreedom):
906     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
907     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
    break;  
    case(Elements):  
    case(ReducedElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);  
    break;  
    case(FaceElements):  
    case(ReducedFaceElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);  
    break;  
    case(Points):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);  
    break;  
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
908        break;        break;
909     }        case(Elements):
910     break;        case(ReducedElements):
911     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
912     switch(target.getFunctionSpace().getTypeCode()) {        break;
913     case(Nodes):        case(FaceElements):
914     case(ReducedNodes):        case(ReducedFaceElements):
915     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
916     case(ReducedDegreesOfFreedom):        break;
917     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
918     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
919     case(Elements):        break;
920     case(ReducedElements):        case(ContactElementsZero):
921     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
922     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
923     case(FaceElements):        break;
924     case(ReducedFaceElements):        case(ContactElementsOne):
925     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
926     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
927     case(Points):        break;
928     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
929     break;           stringstream temp;
930     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
931     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
932     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
933     break;        }
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
934        break;        break;
    }  
    break;  
    case(Elements):  
    if (target.getFunctionSpace().getTypeCode()==Elements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements possible.");  
    }  
    break;  
    case(ReducedElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");  
    }  
    break;  
    case(FaceElements):  
    if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");  
    }  
    break;  
    case(ReducedFaceElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");  
    }  
    break;  
    case(Points):  
    if (target.getFunctionSpace().getTypeCode()==Points) {  
       Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
    }  
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");  
    }  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");  
    }  
    break;  
    case(DegreesOfFreedom):        
    switch(target.getFunctionSpace().getTypeCode()) {  
    case(ReducedDegreesOfFreedom):  
    case(DegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
   
    case(Nodes):  
935     case(ReducedNodes):     case(ReducedNodes):
936     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
937        escript::Data temp=escript::Data(in);        case(Nodes):
938        temp.expand();        case(ReducedNodes):
939        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
940        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
941        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
942     }        break;
943     break;        case(Elements):
944     case(Elements):        case(ReducedElements):
    case(ReducedElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
    } else {  
945        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
946     }        break;
947     break;        case(FaceElements):
948     case(FaceElements):        case(ReducedFaceElements):
    case(ReducedFaceElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
   
    } else {  
949        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
950     }        break;
951     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
952        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
953     }        break;
954     break;        case(ContactElementsZero):
955     case(ContactElementsZero):        case(ReducedContactElementsZero):
    case(ContactElementsOne):  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
    } else {  
956        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
    }  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
957        break;        break;
958     }        case(ContactElementsOne):
959     break;        case(ReducedContactElementsOne):
960     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
961     switch(target.getFunctionSpace().getTypeCode()) {        break;
962     case(Nodes):        default:
963     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
964     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
965     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
966     if (getMPISize()>1) {           break;
967        escript::Data temp=escript::Data(in);        }
968        temp.expand();        break;
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
    } else {  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    }  
    break;  
    case(DegreesOfFreedom):  
    throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
    break;  
    case(ReducedDegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
969     case(Elements):     case(Elements):
970          if (target.getFunctionSpace().getTypeCode()==Elements) {
971             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
972          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
973             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
974          } else {
975             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
976          }
977          break;
978     case(ReducedElements):     case(ReducedElements):
979     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
980        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
981        escriptDataC _in2 = temp.getDataC();        } else {
982        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
983     } else {        }
984        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
985     case(FaceElements):     case(FaceElements):
986          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
987             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
988          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
989             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
990          } else {
991             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
992          }
993          break;
994     case(ReducedFaceElements):     case(ReducedFaceElements):
995     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
996        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
997        escriptDataC _in2 = temp.getDataC();        } else {
998        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
999     } else {        }
1000        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1001     case(Points):     case(Points):
1002     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1003        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1004        escriptDataC _in2 = temp.getDataC();        } else {
1005        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1006     } else {        }
1007        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1008     case(ContactElementsZero):     case(ContactElementsZero):
1009     case(ContactElementsOne):     case(ContactElementsOne):
1010          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1011             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1012          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1013             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1014          } else {
1015             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1016          }
1017          break;
1018     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1019     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1020     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1021        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1022        escriptDataC _in2 = temp.getDataC();        } else {
1023        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1024     } else {        }
1025        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1026     }     case(DegreesOfFreedom):      
1027     break;        switch(target.getFunctionSpace().getTypeCode()) {
1028     default:        case(ReducedDegreesOfFreedom):
1029        stringstream temp;        case(DegreesOfFreedom):
1030        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1031        throw FinleyAdapterException(temp.str());        break;
1032      
1033          case(Nodes):
1034          case(ReducedNodes):
1035          if (getMPISize()>1) {
1036             escript::Data temp=escript::Data(in);
1037             temp.expand();
1038             escriptDataC _in2 = temp.getDataC();
1039             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1040          } else {
1041             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1042          }
1043          break;
1044          case(Elements):
1045          case(ReducedElements):
1046          if (getMPISize()>1) {
1047             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1048             escriptDataC _in2 = temp.getDataC();
1049             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1050          } else {
1051             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1052          }
1053          break;
1054          case(FaceElements):
1055          case(ReducedFaceElements):
1056          if (getMPISize()>1) {
1057             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1058             escriptDataC _in2 = temp.getDataC();
1059             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1060      
1061          } else {
1062             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1063          }
1064          break;
1065          case(Points):
1066          if (getMPISize()>1) {
1067             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1068             escriptDataC _in2 = temp.getDataC();
1069          } else {
1070             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1071          }
1072          break;
1073          case(ContactElementsZero):
1074          case(ContactElementsOne):
1075          case(ReducedContactElementsZero):
1076          case(ReducedContactElementsOne):
1077          if (getMPISize()>1) {
1078             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1079             escriptDataC _in2 = temp.getDataC();
1080             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1081          } else {
1082             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1083          }
1084          break;
1085          default:
1086             stringstream temp;
1087             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1088             throw FinleyAdapterException(temp.str());
1089             break;
1090          }
1091          break;
1092       case(ReducedDegreesOfFreedom):
1093          switch(target.getFunctionSpace().getTypeCode()) {
1094          case(Nodes):
1095          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1096          break;
1097          case(ReducedNodes):
1098          if (getMPISize()>1) {
1099             escript::Data temp=escript::Data(in);
1100             temp.expand();
1101             escriptDataC _in2 = temp.getDataC();
1102             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1103          } else {
1104             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1105          }
1106          break;
1107          case(DegreesOfFreedom):
1108          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1109          break;
1110          case(ReducedDegreesOfFreedom):
1111          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1112          break;
1113          case(Elements):
1114          case(ReducedElements):
1115          if (getMPISize()>1) {
1116             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1117             escriptDataC _in2 = temp.getDataC();
1118             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1119          } else {
1120             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1121          }
1122          break;
1123          case(FaceElements):
1124          case(ReducedFaceElements):
1125          if (getMPISize()>1) {
1126             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1127             escriptDataC _in2 = temp.getDataC();
1128             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1129          } else {
1130             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1131          }
1132          break;
1133          case(Points):
1134          if (getMPISize()>1) {
1135             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1136             escriptDataC _in2 = temp.getDataC();
1137             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1138          } else {
1139             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1140          }
1141          break;
1142          case(ContactElementsZero):
1143          case(ContactElementsOne):
1144          case(ReducedContactElementsZero):
1145          case(ReducedContactElementsOne):
1146          if (getMPISize()>1) {
1147             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1148             escriptDataC _in2 = temp.getDataC();
1149             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1150          } else {
1151             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1152          }
1153          break;
1154          default:
1155             stringstream temp;
1156             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1157             throw FinleyAdapterException(temp.str());
1158             break;
1159          }
1160        break;        break;
    }  
    break;  
1161     default:     default:
1162        stringstream temp;        stringstream temp;
1163        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
# Line 1163  void MeshAdapter::setToX(escript::Data& Line 1181  void MeshAdapter::setToX(escript::Data&
1181        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1182        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1183     } else {     } else {
1184        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1185        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1186        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1187        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1244  void MeshAdapter::interpolateACross(escr Line 1262  void MeshAdapter::interpolateACross(escr
1262  //  //
1263  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1264  //  //
1265  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1266  {  {
1267     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1268     if (argDomain!=*this)     if (argDomain!=*this)
# Line 1257  void MeshAdapter::setToIntegrals(std::ve Line 1275  void MeshAdapter::setToIntegrals(std::ve
1275     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1276     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1277     case(Nodes):     case(Nodes):
1278     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1279     _temp=temp.getDataC();     _temp=temp.getDataC();
1280     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1281     break;     break;
1282     case(ReducedNodes):     case(ReducedNodes):
1283     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1284     _temp=temp.getDataC();     _temp=temp.getDataC();
1285     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1286     break;     break;
# Line 1294  void MeshAdapter::setToIntegrals(std::ve Line 1312  void MeshAdapter::setToIntegrals(std::ve
1312     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1313     break;     break;
1314     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1315     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1316     _temp=temp.getDataC();     _temp=temp.getDataC();
1317     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1318     break;     break;
1319     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1320     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1321     _temp=temp.getDataC();     _temp=temp.getDataC();
1322     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1323     break;     break;
# Line 1331  void MeshAdapter::setToGradient(escript: Line 1349  void MeshAdapter::setToGradient(escript:
1349     escript::Data temp;     escript::Data temp;
1350     if (getMPISize()>1) {     if (getMPISize()>1) {
1351        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1352           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1353           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1354        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1355           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1356           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1357        } else {        } else {
1358           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1453  void MeshAdapter::setNewX(const escript: Line 1471  void MeshAdapter::setNewX(const escript:
1471     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1472     if (newDomain!=*this)     if (newDomain!=*this)
1473        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1474     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1475     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1476           Finley_Mesh_setCoordinates(mesh,&tmp);
1477       } else {
1478           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1479           tmp = new_x_inter.getDataC();
1480           Finley_Mesh_setCoordinates(mesh,&tmp);
1481       }
1482     checkFinleyError();     checkFinleyError();
1483  }  }
1484    
# Line 1472  void MeshAdapter::extractArgsFromDict(co Line 1496  void MeshAdapter::extractArgsFromDict(co
1496    
1497     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1498     for (int i=0; i<numData; ++i) {     for (int i=0; i<numData; ++i) {
1499        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1500        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1501        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1502           throw FinleyAdapterException("Error: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
# Line 1486  void MeshAdapter::extractArgsFromDict(co Line 1510  void MeshAdapter::extractArgsFromDict(co
1510  //  //
1511  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1512  //  //
1513  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1514  {  {
1515     int num_data;     int num_data;
1516     char **names;     char **names;
# Line 1512  void MeshAdapter::saveDX(const std::stri Line 1536  void MeshAdapter::saveDX(const std::stri
1536  //  //
1537  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1538  //  //
1539  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1540  {  {
1541     int num_data;     int num_data;
1542     char **names;     char **names;
# Line 1520  void MeshAdapter::saveVTK(const std::str Line 1544  void MeshAdapter::saveVTK(const std::str
1544     escriptDataC **ptr_data;     escriptDataC **ptr_data;
1545    
1546     extractArgsFromDict(arg, num_data, names, data, ptr_data);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1547     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1548     checkFinleyError();     checkFinleyError();
1549    
1550     /* win32 refactor */     /* win32 refactor */
# Line 1533  void MeshAdapter::saveVTK(const std::str Line 1557  void MeshAdapter::saveVTK(const std::str
1557     TMPMEMFREE(names);     TMPMEMFREE(names);
1558  }  }
1559    
1560    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1561    {
1562    #ifdef ESYS_MPI
1563        index_t myFirstNode=0, myLastNode=0, k=0;
1564        index_t* globalNodeIndex=0;
1565        Finley_Mesh* mesh_p=m_finleyMesh.get();
1566        if (fs_code == FINLEY_REDUCED_NODES)
1567        {
1568        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1569        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1570        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1571        }
1572        else
1573        {
1574        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1575        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1576        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1577        }
1578        k=globalNodeIndex[id];
1579        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1580    #endif
1581        return true;
1582    }
1583    
1584    
1585    
1586  //  //
1587  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1588  //  //
1589  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1590                                                   const int row_blocksize,                                                   const int row_blocksize,
1591                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1592                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1579  SystemMatrixAdapter MeshAdapter::newSyst Line 1629  SystemMatrixAdapter MeshAdapter::newSyst
1629  #endif  #endif
1630     }     }
1631     else {     else {
1632        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1633     }     }
1634     checkPasoError();     checkPasoError();
1635     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1636     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1637       return ASM_ptr(sma);
1638    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1639  }  }
1640    
1641  //  //
1642  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1643  //  //
1644  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1645                                                           const double theta,                                                           const bool useBackwardEuler,
1646                                                           const int blocksize,                                                           const int blocksize,
1647                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1648                                                           const int type) const                                                           const int type) const
# Line 1612  TransportProblemAdapter MeshAdapter::new Line 1664  TransportProblemAdapter MeshAdapter::new
1664    
1665     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1666     checkFinleyError();     checkFinleyError();
1667     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1668     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1669     checkPasoError();     checkPasoError();
1670     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1671     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1672       return ATP_ptr(tpa);
1673    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1674  }  }
1675    
1676  //  //
# Line 1655  bool MeshAdapter::isCellOriented(int fun Line 1709  bool MeshAdapter::isCellOriented(int fun
1709     return false;     return false;
1710  }  }
1711    
1712    bool
1713    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1714    {
1715       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1716        class 1: DOF <-> Nodes
1717        class 2: ReducedDOF <-> ReducedNodes
1718        class 3: Points
1719        class 4: Elements
1720        class 5: ReducedElements
1721        class 6: FaceElements
1722        class 7: ReducedFaceElements
1723        class 8: ContactElementZero <-> ContactElementOne
1724        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1725    
1726       There is also a set of lines. Interpolation is possible down a line but not between lines.
1727       class 1 and 2 belong to all lines so aren't considered.
1728        line 0: class 3
1729        line 1: class 4,5
1730        line 2: class 6,7
1731        line 3: class 8,9
1732    
1733       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1734       eg hasnodes is true if we have at least one instance of Nodes.
1735       */
1736        if (fs.empty())
1737        {
1738            return false;
1739        }
1740        vector<int> hasclass(10);
1741        vector<int> hasline(4);
1742        bool hasnodes=false;
1743        bool hasrednodes=false;
1744        bool hascez=false;
1745        bool hasrcez=false;
1746        for (int i=0;i<fs.size();++i)
1747        {
1748        switch(fs[i])
1749        {
1750        case(Nodes):   hasnodes=true;   // no break is deliberate
1751        case(DegreesOfFreedom):
1752            hasclass[1]=1;
1753            break;
1754        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1755        case(ReducedDegreesOfFreedom):
1756            hasclass[2]=1;
1757            break;
1758        case(Points):
1759            hasline[0]=1;
1760            hasclass[3]=1;
1761            break;
1762        case(Elements):
1763            hasclass[4]=1;
1764            hasline[1]=1;
1765            break;
1766        case(ReducedElements):
1767            hasclass[5]=1;
1768            hasline[1]=1;
1769            break;
1770        case(FaceElements):
1771            hasclass[6]=1;
1772            hasline[2]=1;
1773            break;
1774        case(ReducedFaceElements):
1775            hasclass[7]=1;
1776            hasline[2]=1;
1777            break;
1778        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1779        case(ContactElementsOne):
1780            hasclass[8]=1;
1781            hasline[3]=1;
1782            break;
1783        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1784        case(ReducedContactElementsOne):
1785            hasclass[9]=1;
1786            hasline[3]=1;
1787            break;
1788        default:
1789            return false;
1790        }
1791        }
1792        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1793        // fail if we have more than one leaf group
1794    
1795        if (totlines>1)
1796        {
1797        return false;   // there are at least two branches we can't interpolate between
1798        }
1799        else if (totlines==1)
1800        {
1801        if (hasline[0]==1)      // we have points
1802        {
1803            resultcode=Points;
1804        }
1805        else if (hasline[1]==1)
1806        {
1807            if (hasclass[5]==1)
1808            {
1809            resultcode=ReducedElements;
1810            }
1811            else
1812            {
1813            resultcode=Elements;
1814            }
1815        }
1816        else if (hasline[2]==1)
1817        {
1818            if (hasclass[7]==1)
1819            {
1820            resultcode=ReducedFaceElements;
1821            }
1822            else
1823            {
1824            resultcode=FaceElements;
1825            }
1826        }
1827        else    // so we must be in line3
1828        {
1829            if (hasclass[9]==1)
1830            {
1831            // need something from class 9
1832            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1833            }
1834            else
1835            {
1836            // something from class 8
1837            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1838            }
1839        }
1840        }
1841        else    // totlines==0
1842        {
1843        if (hasclass[2]==1)
1844        {
1845            // something from class 2
1846            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1847        }
1848        else
1849        {   // something from class 1
1850            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1851        }
1852        }
1853        return true;
1854    }
1855    
1856  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1857  {  {
1858     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1859     case(Nodes):     case(Nodes):
1860     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1861     case(Nodes):      case(Nodes):
1862     case(ReducedNodes):      case(ReducedNodes):
1863     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1864     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1865     case(Elements):      case(Elements):
1866     case(ReducedElements):      case(ReducedElements):
1867     case(FaceElements):      case(FaceElements):
1868     case(ReducedFaceElements):      case(ReducedFaceElements):
1869     case(Points):      case(Points):
1870     case(ContactElementsZero):      case(ContactElementsZero):
1871     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1872     case(ContactElementsOne):      case(ContactElementsOne):
1873     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1874     return true;      return true;
1875     default:      default:
1876        stringstream temp;            stringstream temp;
1877        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1878        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1879     }     }
1880     break;     break;
1881     case(ReducedNodes):     case(ReducedNodes):
1882     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1883     case(ReducedNodes):      case(ReducedNodes):
1884     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1885     case(Elements):      case(Elements):
1886     case(ReducedElements):      case(ReducedElements):
1887     case(FaceElements):      case(FaceElements):
1888     case(ReducedFaceElements):      case(ReducedFaceElements):
1889     case(Points):      case(Points):
1890     case(ContactElementsZero):      case(ContactElementsZero):
1891     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1892     case(ContactElementsOne):      case(ContactElementsOne):
1893     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1894     return true;      return true;
1895     case(Nodes):      case(Nodes):
1896     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1897     return false;      return false;
1898     default:      default:
1899        stringstream temp;          stringstream temp;
1900        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1901        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1902     }     }
1903     break;     break;
1904     case(Elements):     case(Elements):
1905     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1906        return true;        return true;
1907     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1908        return true;        return true;
1909     } else {          } else {
1910        return false;            return false;
1911     }          }
1912     case(ReducedElements):     case(ReducedElements):
1913     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1914        return true;        return true;
1915     } else {      } else {
1916        return false;            return false;
1917     }      }
1918     case(FaceElements):     case(FaceElements):
1919     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1920        return true;              return true;
1921     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1922        return true;              return true;
1923     } else {      } else {
1924        return false;              return false;
1925     }      }
1926     case(ReducedFaceElements):     case(ReducedFaceElements):
1927     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1928        return true;              return true;
1929     } else {      } else {
1930        return false;          return false;
1931     }      }
1932     case(Points):     case(Points):
1933     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1934        return true;              return true;
1935     } else {      } else {
1936        return false;              return false;
1937     }      }
1938     case(ContactElementsZero):     case(ContactElementsZero):
1939     case(ContactElementsOne):     case(ContactElementsOne):
1940     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1941        return true;              return true;
1942     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1943        return true;              return true;
1944     } else {      } else {
1945        return false;              return false;
1946     }      }
1947     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1948     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1949     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1950        return true;              return true;
1951     } else {      } else {
1952        return false;              return false;
1953     }      }
1954     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1955     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1956     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1957     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1958     case(Nodes):      case(Nodes):
1959     case(ReducedNodes):      case(ReducedNodes):
1960     case(Elements):      case(Elements):
1961     case(ReducedElements):      case(ReducedElements):
1962     case(Points):      case(Points):
1963     case(FaceElements):      case(FaceElements):
1964     case(ReducedFaceElements):      case(ReducedFaceElements):
1965     case(ContactElementsZero):      case(ContactElementsZero):
1966     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1967     case(ContactElementsOne):      case(ContactElementsOne):
1968     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1969     return true;      return true;
1970     default:      default:
1971        stringstream temp;          stringstream temp;
1972        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1973        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1974     }      }
1975     break;      break;
1976     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1977     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1978     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1979     case(ReducedNodes):      case(ReducedNodes):
1980     case(Elements):      case(Elements):
1981     case(ReducedElements):      case(ReducedElements):
1982     case(FaceElements):      case(FaceElements):
1983     case(ReducedFaceElements):      case(ReducedFaceElements):
1984     case(Points):      case(Points):
1985     case(ContactElementsZero):      case(ContactElementsZero):
1986     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1987     case(ContactElementsOne):      case(ContactElementsOne):
1988     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1989     return true;      return true;
1990     case(Nodes):      case(Nodes):
1991     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1992     return false;      return false;
1993     default:      default:
1994        stringstream temp;          stringstream temp;
1995        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1996        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1997     }      }
1998     break;      break;
1999     default:     default:
2000        stringstream temp;        stringstream temp;
2001        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
# Line 1838  int MeshAdapter::getSystemMatrixTypeId(c Line 2036  int MeshAdapter::getSystemMatrixTypeId(c
2036  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2037  {  {
2038     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2039     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2040     checkPasoError();     checkPasoError();
2041     return out;     return out;
2042  }  }
2043    
2044  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2045  {  {
2046     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2047  }  }
2048    
2049  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2050  {  {
2051     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2052  }  }
2053    
2054  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2055  {  {
2056     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2057  }  }
2058    
2059  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2060  {  {
2061     int *out = NULL;     int *out = NULL;
2062     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2017  void MeshAdapter::setTags(const int func Line 2215  void MeshAdapter::setTags(const int func
2215     return;     return;
2216  }  }
2217    
2218  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2219  {  {
2220     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2221     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
# Line 2025  void MeshAdapter::setTagMap(const std::s Line 2223  void MeshAdapter::setTagMap(const std::s
2223     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2224  }  }
2225    
2226  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2227  {  {
2228     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2229     int tag=0;     int tag=0;
# Line 2035  int MeshAdapter::getTag(const std::strin Line 2233  int MeshAdapter::getTag(const std::strin
2233     return tag;     return tag;
2234  }  }
2235    
2236  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2237  {  {
2238     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2239     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2240  }  }
2241    
2242  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2243  {  {
2244     stringstream temp;     stringstream temp;
2245     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2095  int MeshAdapter::getNumberOfTagsInUse(in Line 2293  int MeshAdapter::getNumberOfTagsInUse(in
2293    }    }
2294    return numTags;    return numTags;
2295  }  }
2296  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2297    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2298  {  {
2299    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2300    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2161  bool MeshAdapter::canTag(int functionSpa Line 2360  bool MeshAdapter::canTag(int functionSpa
2360    }    }
2361  }  }
2362    
2363    AbstractDomain::StatusType MeshAdapter::getStatus() const
2364    {
2365      Finley_Mesh* mesh=m_finleyMesh.get();
2366      return Finley_Mesh_getStatus(mesh);
2367    }
2368    
2369    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2370    {
2371      
2372      Finley_Mesh* mesh=m_finleyMesh.get();
2373      int order =-1;
2374      switch(functionSpaceCode) {
2375       case(Nodes):
2376       case(DegreesOfFreedom):
2377              order=mesh->approximationOrder;
2378              break;
2379       case(ReducedNodes):
2380       case(ReducedDegreesOfFreedom):
2381              order=mesh->reducedApproximationOrder;
2382              break;
2383       case(Elements):
2384       case(FaceElements):
2385       case(Points):
2386       case(ContactElementsZero):
2387       case(ContactElementsOne):
2388              order=mesh->integrationOrder;
2389              break;
2390       case(ReducedElements):
2391       case(ReducedFaceElements):
2392       case(ReducedContactElementsZero):
2393       case(ReducedContactElementsOne):
2394              order=mesh->reducedIntegrationOrder;
2395              break;
2396       default:
2397          stringstream temp;
2398          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2399          throw FinleyAdapterException(temp.str());
2400      }
2401      return order;
2402    }
2403    
2404    bool MeshAdapter::supportsContactElements() const
2405    {
2406      return true;
2407    }
2408    
2409    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2410    {
2411      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2412    }
2413    
2414    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2415    {
2416      Finley_ReferenceElementSet_dealloc(m_refSet);
2417    }
2418    
2419  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2315  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26