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

Diff of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

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

revision 2641 by jfenwick, Mon Aug 31 07:41:49 2009 UTC revision 3317 by caltinay, Thu Oct 28 00:50:41 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2009 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
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
21  extern "C" {  extern "C" {
22  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
# Line 85  int MeshAdapter::getMPIRank() const Line 81  int MeshAdapter::getMPIRank() const
81  }  }
82  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
83  {  {
84  #ifdef PASO_MPI  #ifdef ESYS_MPI
85     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
86  #endif  #endif
87     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 92  bool MeshAdapter::onMasterProcessor() co
92  }  }
93    
94    
95  #ifdef PASO_MPI  #ifdef ESYS_MPI
96    MPI_Comm    MPI_Comm
97  #else  #else
98    unsigned int    unsigned int
99  #endif  #endif
100  MeshAdapter::getMPIComm() const  MeshAdapter::getMPIComm() const
101  {  {
102  #ifdef PASO_MPI  #ifdef ESYS_MPI
103      return m_finleyMesh->MPIInfo->comm;      return m_finleyMesh->MPIInfo->comm;
104  #else  #else
105      return 0;      return 0;
# Line 115  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 111  Finley_Mesh* MeshAdapter::getFinley_Mesh
111     return m_finleyMesh.get();     return m_finleyMesh.get();
112  }  }
113    
114  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
115  {  {
116     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;
117     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 129  void MeshAdapter::Print_Mesh_Info(const Line 125  void MeshAdapter::Print_Mesh_Info(const
125     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
126  }  }
127    
128  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
129  {  {
130  #ifdef USE_NETCDF  #ifdef USE_NETCDF
131     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 149  void MeshAdapter::dump(const std::string Line 145  void MeshAdapter::dump(const std::string
145     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
146     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
147     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
148  #ifdef PASO_MPI  #ifdef ESYS_MPI
149     MPI_Status status;     MPI_Status status;
150  #endif  #endif
151    
152  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
153  #ifdef PASO_MPI  #ifdef ESYS_MPI
154     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);
155  #endif  #endif
156    
157     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
158                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
159    
160     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
161     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
162     num_Tags = 0;     num_Tags = 0;
163     if (tag_map) {     while (tag_map) {
164        while (tag_map) {        num_Tags++;
165           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
166     }     }
167    
168     // NetCDF error handler     // NetCDF error handler
169     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
170     // Create the file.     // Create the file.
171     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
172       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
173     // check if writing was successful     // check if writing was successful
174     if (!dataFile.is_valid())     if (!dataFile.is_valid())
175        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
176    
177     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
178       // dim_Elements only appears if > 0)
179     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
180        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
181     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
182        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
183     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)) )
184        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
185     if (num_Elements>0)     if (num_Elements>0)
186        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
187           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
188     if (num_FaceElements>0)     if (num_FaceElements>0)
189        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
190           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
191     if (num_ContactElements>0)     if (num_ContactElements>0)
192        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
193           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
194     if (num_Points>0)     if (num_Points>0)
195        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
196           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
197     if (num_Elements>0)     if (num_Elements>0)
198        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
199           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
200     if (num_FaceElements>0)     if (num_FaceElements>0)
201        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
202           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
203     if (num_ContactElements>0)     if (num_ContactElements>0)
204        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
205           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
206     if (num_Tags>0)     if (num_Tags>0)
207        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
208           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
209    
210     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
211     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
212        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
213     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
214        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
215     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
216        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
217     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
218        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
219     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
220        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
221     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
222        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
223     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
224        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
225     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
226        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
227     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
228        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
229     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
230        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
231     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
232        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
233     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
234        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
235     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
236        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
237     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
238        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
239     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
240        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
241     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
242        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
243     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
245     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
247     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
248        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
249    
250     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
251    
252     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
253       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
254          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
255       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
256       if (! (ids->put(int_ptr, mpi_size+1)) )
257          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
258    
259       // Nodes degreesOfFreedomDistribution
260       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
261          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
262       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
263       if (! (ids->put(int_ptr, mpi_size+1)) )
264          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
265    
266       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
267       // (it treats them as NC_UNLIMITED)
268     if (numNodes>0) {     if (numNodes>0) {
269    
270        // Nodes Id        // Nodes Id
271        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
272           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
273        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
274        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
275           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
276    
277        // Nodes Tag        // Nodes Tag
278        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
279           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
280        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
281        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
282           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
283    
284        // Nodes gDOF        // Nodes gDOF
285        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
286           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
287        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
288        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
289           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
290    
291        // Nodes global node index        // Nodes global node index
292        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
293           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
294        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
295        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
296           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
297    
298        // Nodes grDof        // Nodes grDof
299        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
300           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
301        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
302        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
303           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
304    
305        // Nodes grNI        // Nodes grNI
306        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
307           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
308        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
309        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
310           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
311    
312        // Nodes Coordinates        // Nodes Coordinates
313        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
314           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
315        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
316           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);  
317    
318     }     }
319    
# Line 326  void MeshAdapter::dump(const std::string Line 323  void MeshAdapter::dump(const std::string
323    
324        // Elements_Id        // Elements_Id
325        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
326           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
327        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
328        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
329           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
330    
331        // Elements_Tag        // Elements_Tag
332        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
333           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
334        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
335        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
336           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
337    
338        // Elements_Owner        // Elements_Owner
339        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
340           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
341        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
342        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
343           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
344    
345        // Elements_Color        // Elements_Color
346        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
347           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
348        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
349        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
350           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
351    
352        // Elements_Nodes        // Elements_Nodes
353        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
354           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
355        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
356           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
357    
358     }     }
359    
# Line 366  void MeshAdapter::dump(const std::string Line 363  void MeshAdapter::dump(const std::string
363    
364        // FaceElements_Id        // FaceElements_Id
365        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
366           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
367        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
368        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
369           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
370    
371        // FaceElements_Tag        // FaceElements_Tag
372        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
373           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
374        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
375        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
376           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
377    
378        // FaceElements_Owner        // FaceElements_Owner
379        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
380           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
381        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
382        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
383           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
384    
385        // FaceElements_Color        // FaceElements_Color
386        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
387           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
388        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
389        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
390           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
391    
392        // FaceElements_Nodes        // FaceElements_Nodes
393        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
394           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
395        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
396           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
397    
398     }     }
399    
# Line 406  void MeshAdapter::dump(const std::string Line 403  void MeshAdapter::dump(const std::string
403    
404        // ContactElements_Id        // ContactElements_Id
405        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
406           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
407        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
408        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
409           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
410    
411        // ContactElements_Tag        // ContactElements_Tag
412        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
413           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
414        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
415        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
416           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
417    
418        // ContactElements_Owner        // ContactElements_Owner
419        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
420           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
421        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
422        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
423           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
424    
425        // ContactElements_Color        // ContactElements_Color
426        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
427           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
428        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
429        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
430           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
431    
432        // ContactElements_Nodes        // ContactElements_Nodes
433        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
434           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
435        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
436           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
437    
438     }     }
439    
# Line 448  void MeshAdapter::dump(const std::string Line 445  void MeshAdapter::dump(const std::string
445    
446        // Points_Id        // Points_Id
447        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
448           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
449        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
450        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
451           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
452    
453        // Points_Tag        // Points_Tag
454        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
455           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
456        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
457        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
458           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
459    
460        // Points_Owner        // Points_Owner
461        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
462           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
463        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
464        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
465           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
466    
467        // Points_Color        // Points_Color
468        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
469           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
470        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
471        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
472           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
473    
474        // Points_Nodes        // Points_Nodes
475        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
476        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
477           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
478        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
479           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
480    
481     }     }
482    
# Line 503  void MeshAdapter::dump(const std::string Line 500  void MeshAdapter::dump(const std::string
500    
501        // Tags_keys        // Tags_keys
502        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
503           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
504        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
505        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
506           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
507    
508        // Tags_names_*        // Tags_names_*
509        // 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
510        // 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
511          // manual doesn't tell how to do an array of strings
512        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
513        if (tag_map) {        if (tag_map) {
514           int i = 0;           int i = 0;
515           while (tag_map) {           while (tag_map) {
516              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
517              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
518                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
519              tag_map=tag_map->next;              tag_map=tag_map->next;
520              i++;              i++;
521           }           }
522        }        }
523    
524        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
525     }     }
526    
527  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
528  #ifdef PASO_MPI  #ifdef ESYS_MPI
529     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);
530  #endif  #endif
531    
# Line 693  pair<int,int> MeshAdapter::getDataShape( Line 690  pair<int,int> MeshAdapter::getDataShape(
690     case(Elements):     case(Elements):
691     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
692        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
693        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
694     }     }
695     break;     break;
696     case(ReducedElements):     case(ReducedElements):
697     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
698        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
699        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
700     }     }
701     break;     break;
702     case(FaceElements):     case(FaceElements):
703     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
704        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
705        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
706     }     }
707     break;     break;
708     case(ReducedFaceElements):     case(ReducedFaceElements):
709     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
710        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
711        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
712     }     }
713     break;     break;
# Line 722  pair<int,int> MeshAdapter::getDataShape( Line 719  pair<int,int> MeshAdapter::getDataShape(
719     break;     break;
720     case(ContactElementsZero):     case(ContactElementsZero):
721     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
722        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
723        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
724     }     }
725     break;     break;
726     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
727     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
728        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
729        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
730     }     }
731     break;     break;
732     case(ContactElementsOne):     case(ContactElementsOne):
733     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
734        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
735        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
736     }     }
737     break;     break;
738     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
739     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
740        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
741        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
742     }     }
743     break;     break;
# Line 769  pair<int,int> MeshAdapter::getDataShape( Line 766  pair<int,int> MeshAdapter::getDataShape(
766  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
767  //  //
768  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
769                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
770                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
771                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
772                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
773  {  {
774       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
775       if (smat==0)
776       {
777        throw FinleyAdapterException("finley only supports its own system matrices.");
778       }
779     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
780     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
781     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 788  void MeshAdapter::addPDEToSystem( Line 790  void MeshAdapter::addPDEToSystem(
790    
791     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
792    
793     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
794     checkFinleyError();     checkFinleyError();
795    
796     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
797     checkFinleyError();     checkFinleyError();
798    
799     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
800     checkFinleyError();     checkFinleyError();
801  }  }
802    
# Line 810  void  MeshAdapter::addPDEToLumpedSystem( Line 812  void  MeshAdapter::addPDEToLumpedSystem(
812     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
813    
814     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
815       checkFinleyError();
816      
817     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
   
818     checkFinleyError();     checkFinleyError();
819    
820  }  }
821    
822    
# Line 842  void MeshAdapter::addPDEToRHS( escript:: Line 846  void MeshAdapter::addPDEToRHS( escript::
846  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
847  //  //
848  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
849                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
850                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
851                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
852                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
853                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
854  {  {
855       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
856       if (tpa==0)
857       {
858        throw FinleyAdapterException("finley only supports its own transport problems.");
859       }
860    
861    
862     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
863     source.expand();     source.expand();
864     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 864  void MeshAdapter::addPDEToTransportProbl Line 875  void MeshAdapter::addPDEToTransportProbl
875     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
876    
877     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
878     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
879    
880     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 );
881     checkFinleyError();     checkFinleyError();
# Line 1041  void MeshAdapter::interpolateOnDomain(es Line 1052  void MeshAdapter::interpolateOnDomain(es
1052        case(Elements):        case(Elements):
1053        case(ReducedElements):        case(ReducedElements):
1054        if (getMPISize()>1) {        if (getMPISize()>1) {
1055           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1056           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1057           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1058        } else {        } else {
# Line 1051  void MeshAdapter::interpolateOnDomain(es Line 1062  void MeshAdapter::interpolateOnDomain(es
1062        case(FaceElements):        case(FaceElements):
1063        case(ReducedFaceElements):        case(ReducedFaceElements):
1064        if (getMPISize()>1) {        if (getMPISize()>1) {
1065           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1066           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1067           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1068        
# Line 1061  void MeshAdapter::interpolateOnDomain(es Line 1072  void MeshAdapter::interpolateOnDomain(es
1072        break;        break;
1073        case(Points):        case(Points):
1074        if (getMPISize()>1) {        if (getMPISize()>1) {
1075           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1076           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1077        } else {        } else {
1078           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
# Line 1072  void MeshAdapter::interpolateOnDomain(es Line 1083  void MeshAdapter::interpolateOnDomain(es
1083        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1084        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1085        if (getMPISize()>1) {        if (getMPISize()>1) {
1086           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1087           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1088           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1089        } else {        } else {
# Line 1110  void MeshAdapter::interpolateOnDomain(es Line 1121  void MeshAdapter::interpolateOnDomain(es
1121        case(Elements):        case(Elements):
1122        case(ReducedElements):        case(ReducedElements):
1123        if (getMPISize()>1) {        if (getMPISize()>1) {
1124           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1125           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1126           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1127        } else {        } else {
# Line 1120  void MeshAdapter::interpolateOnDomain(es Line 1131  void MeshAdapter::interpolateOnDomain(es
1131        case(FaceElements):        case(FaceElements):
1132        case(ReducedFaceElements):        case(ReducedFaceElements):
1133        if (getMPISize()>1) {        if (getMPISize()>1) {
1134           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1135           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1136           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1137        } else {        } else {
# Line 1129  void MeshAdapter::interpolateOnDomain(es Line 1140  void MeshAdapter::interpolateOnDomain(es
1140        break;        break;
1141        case(Points):        case(Points):
1142        if (getMPISize()>1) {        if (getMPISize()>1) {
1143           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1144           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1145           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1146        } else {        } else {
# Line 1141  void MeshAdapter::interpolateOnDomain(es Line 1152  void MeshAdapter::interpolateOnDomain(es
1152        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1153        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1154        if (getMPISize()>1) {        if (getMPISize()>1) {
1155           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1156           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1157           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1158        } else {        } else {
# Line 1178  void MeshAdapter::setToX(escript::Data& Line 1189  void MeshAdapter::setToX(escript::Data&
1189        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1190        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1191     } else {     } else {
1192        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1193        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1194        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1195        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1259  void MeshAdapter::interpolateACross(escr Line 1270  void MeshAdapter::interpolateACross(escr
1270  //  //
1271  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1272  //  //
1273  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1274  {  {
1275     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1276     if (argDomain!=*this)     if (argDomain!=*this)
# Line 1272  void MeshAdapter::setToIntegrals(std::ve Line 1283  void MeshAdapter::setToIntegrals(std::ve
1283     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1284     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1285     case(Nodes):     case(Nodes):
1286     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1287     _temp=temp.getDataC();     _temp=temp.getDataC();
1288     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1289     break;     break;
1290     case(ReducedNodes):     case(ReducedNodes):
1291     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1292     _temp=temp.getDataC();     _temp=temp.getDataC();
1293     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1294     break;     break;
# Line 1309  void MeshAdapter::setToIntegrals(std::ve Line 1320  void MeshAdapter::setToIntegrals(std::ve
1320     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1321     break;     break;
1322     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1323     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1324     _temp=temp.getDataC();     _temp=temp.getDataC();
1325     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1326     break;     break;
1327     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1328     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1329     _temp=temp.getDataC();     _temp=temp.getDataC();
1330     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1331     break;     break;
# Line 1346  void MeshAdapter::setToGradient(escript: Line 1357  void MeshAdapter::setToGradient(escript:
1357     escript::Data temp;     escript::Data temp;
1358     if (getMPISize()>1) {     if (getMPISize()>1) {
1359        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1360           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1361           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1362        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1363           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1364           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1365        } else {        } else {
1366           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1468  void MeshAdapter::setNewX(const escript: Line 1479  void MeshAdapter::setNewX(const escript:
1479     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1480     if (newDomain!=*this)     if (newDomain!=*this)
1481        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1482     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1483         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1484         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1485     } else {     } else {
1486         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1487         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1488         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1489     }     }
# Line 1493  void MeshAdapter::extractArgsFromDict(co Line 1504  void MeshAdapter::extractArgsFromDict(co
1504    
1505     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1506     for (int i=0; i<numData; ++i) {     for (int i=0; i<numData; ++i) {
1507        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1508        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1509        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1510           throw FinleyAdapterException("Error: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
# Line 1507  void MeshAdapter::extractArgsFromDict(co Line 1518  void MeshAdapter::extractArgsFromDict(co
1518  //  //
1519  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1520  //  //
1521  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
1522  {  {
1523     int num_data;     int num_data;
1524     char **names;     char **names;
# Line 1533  void MeshAdapter::saveDX(const std::stri Line 1544  void MeshAdapter::saveDX(const std::stri
1544  //  //
1545  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1546  //  //
1547  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1548  {  {
1549     int num_data;     int num_data;
1550     char **names;     char **names;
# Line 1554  void MeshAdapter::saveVTK(const std::str Line 1565  void MeshAdapter::saveVTK(const std::str
1565     TMPMEMFREE(names);     TMPMEMFREE(names);
1566  }  }
1567    
1568    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1569    {
1570    #ifdef ESYS_MPI
1571        index_t myFirstNode=0, myLastNode=0, k=0;
1572        index_t* globalNodeIndex=0;
1573        Finley_Mesh* mesh_p=m_finleyMesh.get();
1574        if (fs_code == FINLEY_REDUCED_NODES)
1575        {
1576        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1577        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1578        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1579        }
1580        else
1581        {
1582        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1583        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1584        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1585        }
1586        k=globalNodeIndex[id];
1587        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1588    #endif
1589        return true;
1590    }
1591    
1592    
1593    
1594  //  //
1595  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1596  //  //
1597  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1598                                                   const int row_blocksize,                                                   const int row_blocksize,
1599                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1600                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1604  SystemMatrixAdapter MeshAdapter::newSyst Line 1641  SystemMatrixAdapter MeshAdapter::newSyst
1641     }     }
1642     checkPasoError();     checkPasoError();
1643     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1644     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1645       return ASM_ptr(sma);
1646    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1647  }  }
1648    
1649  //  //
1650  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1651  //  //
1652  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1653                                                           const double theta,                                                           const bool useBackwardEuler,
1654                                                           const int blocksize,                                                           const int blocksize,
1655                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1656                                                           const int type) const                                                           const int type) const
# Line 1633  TransportProblemAdapter MeshAdapter::new Line 1672  TransportProblemAdapter MeshAdapter::new
1672    
1673     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1674     checkFinleyError();     checkFinleyError();
1675     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1676     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1677     checkPasoError();     checkPasoError();
1678     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1679     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1680       return ATP_ptr(tpa);
1681    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1682  }  }
1683    
1684  //  //
# Line 1677  bool MeshAdapter::isCellOriented(int fun Line 1718  bool MeshAdapter::isCellOriented(int fun
1718  }  }
1719    
1720  bool  bool
1721  MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1722  {  {
1723     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1724      class 1: DOF <-> Nodes      class 1: DOF <-> Nodes
# Line 1700  MeshAdapter::commonFunctionSpace(const s Line 1741  MeshAdapter::commonFunctionSpace(const s
1741     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1742     eg hasnodes is true if we have at least one instance of Nodes.     eg hasnodes is true if we have at least one instance of Nodes.
1743     */     */
1744      if (fs.size()==0)      if (fs.empty())
1745      {      {
1746      return false;          return false;
1747      }      }
1748      std::vector<int> hasclass(10);      vector<int> hasclass(10);
1749      std::vector<int> hasline(4);          vector<int> hasline(4);
1750      bool hasnodes=false;      bool hasnodes=false;
1751      bool hasrednodes=false;      bool hasrednodes=false;
1752      bool hascez=false;      bool hascez=false;
# Line 1782  MeshAdapter::commonFunctionSpace(const s Line 1823  MeshAdapter::commonFunctionSpace(const s
1823      }      }
1824      else if (hasline[2]==1)      else if (hasline[2]==1)
1825      {      {
1826          if (hasclass[7]==ReducedFaceElements)          if (hasclass[7]==1)
1827          {          {
1828          resultcode=ReducedFaceElements;          resultcode=ReducedFaceElements;
1829          }          }
# Line 2003  int MeshAdapter::getSystemMatrixTypeId(c Line 2044  int MeshAdapter::getSystemMatrixTypeId(c
2044  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
2045  {  {
2046     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2047     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);
2048     checkPasoError();     checkPasoError();
2049     return out;     return out;
2050  }  }
2051    
2052  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2053  {  {
2054     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2055  }  }
2056    
2057  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2058  {  {
2059     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2060  }  }
2061    
2062  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2063  {  {
2064     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2065  }  }
2066    
2067  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2182  void MeshAdapter::setTags(const int func Line 2223  void MeshAdapter::setTags(const int func
2223     return;     return;
2224  }  }
2225    
2226  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2227  {  {
2228     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2229     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
# Line 2190  void MeshAdapter::setTagMap(const std::s Line 2231  void MeshAdapter::setTagMap(const std::s
2231     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2232  }  }
2233    
2234  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2235  {  {
2236     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2237     int tag=0;     int tag=0;
# Line 2200  int MeshAdapter::getTag(const std::strin Line 2241  int MeshAdapter::getTag(const std::strin
2241     return tag;     return tag;
2242  }  }
2243    
2244  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2245  {  {
2246     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2247     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2248  }  }
2249    
2250  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2251  {  {
2252     stringstream temp;     stringstream temp;
2253     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2333  AbstractDomain::StatusType MeshAdapter:: Line 2374  AbstractDomain::StatusType MeshAdapter::
2374    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2375  }  }
2376    
2377    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2378    {
2379      
2380      Finley_Mesh* mesh=m_finleyMesh.get();
2381      int order =-1;
2382      switch(functionSpaceCode) {
2383       case(Nodes):
2384       case(DegreesOfFreedom):
2385              order=mesh->approximationOrder;
2386              break;
2387       case(ReducedNodes):
2388       case(ReducedDegreesOfFreedom):
2389              order=mesh->reducedApproximationOrder;
2390              break;
2391       case(Elements):
2392       case(FaceElements):
2393       case(Points):
2394       case(ContactElementsZero):
2395       case(ContactElementsOne):
2396              order=mesh->integrationOrder;
2397              break;
2398       case(ReducedElements):
2399       case(ReducedFaceElements):
2400       case(ReducedContactElementsZero):
2401       case(ReducedContactElementsOne):
2402              order=mesh->reducedIntegrationOrder;
2403              break;
2404       default:
2405          stringstream temp;
2406          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2407          throw FinleyAdapterException(temp.str());
2408      }
2409      return order;
2410    }
2411    
2412    bool MeshAdapter::supportsContactElements() const
2413    {
2414      return true;
2415    }
2416    
2417    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2418    {
2419      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2420    }
2421    
2422    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2423    {
2424      Finley_ReferenceElementSet_dealloc(m_refSet);
2425    }
2426    
2427  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2641  
changed lines
  Added in v.3317

  ViewVC Help
Powered by ViewVC 1.1.26