/[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 3269 by jfenwick, Wed Oct 13 03:21:50 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 754  pair<int,int> MeshAdapter::getDataShape( Line 770  pair<int,int> MeshAdapter::getDataShape(
770  // 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:
771  //  //
772  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
773                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
774                                   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,
775                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
776                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
777  {  {
778       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
779       if (smat==0)
780       {
781        throw FinleyAdapterException("finley only supports its own system matrices.");
782       }
783     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
784     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
785     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 773  void MeshAdapter::addPDEToSystem( Line 794  void MeshAdapter::addPDEToSystem(
794    
795     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
796    
797     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 );
798     checkFinleyError();     checkFinleyError();
799    
800     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 );
801     checkFinleyError();     checkFinleyError();
802    
803     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 );
804     checkFinleyError();     checkFinleyError();
805  }  }
806    
# Line 795  void  MeshAdapter::addPDEToLumpedSystem( Line 816  void  MeshAdapter::addPDEToLumpedSystem(
816     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
817    
818     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
819       checkFinleyError();
820      
821     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
   
822     checkFinleyError();     checkFinleyError();
823    
824  }  }
825    
826    
# Line 827  void MeshAdapter::addPDEToRHS( escript:: Line 850  void MeshAdapter::addPDEToRHS( escript::
850  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
851  //  //
852  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
853                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
854                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
855                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
856                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
857                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
858  {  {
859       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
860       if (tpa==0)
861       {
862        throw FinleyAdapterException("finley only supports its own transport problems.");
863       }
864    
865    
866     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
867     source.expand();     source.expand();
868     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 849  void MeshAdapter::addPDEToTransportProbl Line 879  void MeshAdapter::addPDEToTransportProbl
879     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
880    
881     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
882     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
883    
884     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 );
885     checkFinleyError();     checkFinleyError();
# Line 881  void MeshAdapter::interpolateOnDomain(es Line 911  void MeshAdapter::interpolateOnDomain(es
911     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
912     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
913     case(Nodes):     case(Nodes):
914     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
915     case(Nodes):        case(Nodes):
916     case(ReducedNodes):        case(ReducedNodes):
917     case(DegreesOfFreedom):        case(DegreesOfFreedom):
918     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
919     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());  
920        break;        break;
921     }        case(Elements):
922     break;        case(ReducedElements):
923     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
924     switch(target.getFunctionSpace().getTypeCode()) {        break;
925     case(Nodes):        case(FaceElements):
926     case(ReducedNodes):        case(ReducedFaceElements):
927     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
928     case(ReducedDegreesOfFreedom):        break;
929     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
930     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
931     case(Elements):        break;
932     case(ReducedElements):        case(ContactElementsZero):
933     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
934     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
935     case(FaceElements):        break;
936     case(ReducedFaceElements):        case(ContactElementsOne):
937     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
938     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
939     case(Points):        break;
940     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
941     break;           stringstream temp;
942     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
943     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
944     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
945     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());  
946        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):  
947     case(ReducedNodes):     case(ReducedNodes):
948     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
949        escript::Data temp=escript::Data(in);        case(Nodes):
950        temp.expand();        case(ReducedNodes):
951        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
952        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
953        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
954     }        break;
955     break;        case(Elements):
956     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 {  
957        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
958     }        break;
959     break;        case(FaceElements):
960     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 {  
961        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
962     }        break;
963     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
964        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
965     }        break;
966     break;        case(ContactElementsZero):
967     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 {  
968        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());  
969        break;        break;
970     }        case(ContactElementsOne):
971     break;        case(ReducedContactElementsOne):
972     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
973     switch(target.getFunctionSpace().getTypeCode()) {        break;
974     case(Nodes):        default:
975     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
976     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
977     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
978     if (getMPISize()>1) {           break;
979        escript::Data temp=escript::Data(in);        }
980        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;  
981     case(Elements):     case(Elements):
982          if (target.getFunctionSpace().getTypeCode()==Elements) {
983             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
984          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
985             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
986          } else {
987             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
988          }
989          break;
990     case(ReducedElements):     case(ReducedElements):
991     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
992        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
993        escriptDataC _in2 = temp.getDataC();        } else {
994        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
995     } else {        }
996        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
997     case(FaceElements):     case(FaceElements):
998          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
999             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1000          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1001             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1002          } else {
1003             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1004          }
1005          break;
1006     case(ReducedFaceElements):     case(ReducedFaceElements):
1007     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1008        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1009        escriptDataC _in2 = temp.getDataC();        } else {
1010        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1011     } else {        }
1012        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1013     case(Points):     case(Points):
1014     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1015        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1016        escriptDataC _in2 = temp.getDataC();        } else {
1017        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1018     } else {        }
1019        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1020     case(ContactElementsZero):     case(ContactElementsZero):
1021     case(ContactElementsOne):     case(ContactElementsOne):
1022          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1023             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1024          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1025             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1026          } else {
1027             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1028          }
1029          break;
1030     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1031     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1032     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1033        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1034        escriptDataC _in2 = temp.getDataC();        } else {
1035        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1036     } else {        }
1037        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1038     }     case(DegreesOfFreedom):      
1039     break;        switch(target.getFunctionSpace().getTypeCode()) {
1040     default:        case(ReducedDegreesOfFreedom):
1041        stringstream temp;        case(DegreesOfFreedom):
1042        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1043        throw FinleyAdapterException(temp.str());        break;
1044      
1045          case(Nodes):
1046          case(ReducedNodes):
1047          if (getMPISize()>1) {
1048             escript::Data temp=escript::Data(in);
1049             temp.expand();
1050             escriptDataC _in2 = temp.getDataC();
1051             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1052          } else {
1053             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1054          }
1055          break;
1056          case(Elements):
1057          case(ReducedElements):
1058          if (getMPISize()>1) {
1059             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1060             escriptDataC _in2 = temp.getDataC();
1061             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1062          } else {
1063             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1064          }
1065          break;
1066          case(FaceElements):
1067          case(ReducedFaceElements):
1068          if (getMPISize()>1) {
1069             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1070             escriptDataC _in2 = temp.getDataC();
1071             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1072      
1073          } else {
1074             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1075          }
1076          break;
1077          case(Points):
1078          if (getMPISize()>1) {
1079             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1080             escriptDataC _in2 = temp.getDataC();
1081          } else {
1082             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1083          }
1084          break;
1085          case(ContactElementsZero):
1086          case(ContactElementsOne):
1087          case(ReducedContactElementsZero):
1088          case(ReducedContactElementsOne):
1089          if (getMPISize()>1) {
1090             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1091             escriptDataC _in2 = temp.getDataC();
1092             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1093          } else {
1094             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1095          }
1096          break;
1097          default:
1098             stringstream temp;
1099             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1100             throw FinleyAdapterException(temp.str());
1101             break;
1102          }
1103          break;
1104       case(ReducedDegreesOfFreedom):
1105          switch(target.getFunctionSpace().getTypeCode()) {
1106          case(Nodes):
1107          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1108          break;
1109          case(ReducedNodes):
1110          if (getMPISize()>1) {
1111             escript::Data temp=escript::Data(in);
1112             temp.expand();
1113             escriptDataC _in2 = temp.getDataC();
1114             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1115          } else {
1116             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1117          }
1118          break;
1119          case(DegreesOfFreedom):
1120          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1121          break;
1122          case(ReducedDegreesOfFreedom):
1123          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1124          break;
1125          case(Elements):
1126          case(ReducedElements):
1127          if (getMPISize()>1) {
1128             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1129             escriptDataC _in2 = temp.getDataC();
1130             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1131          } else {
1132             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1133          }
1134          break;
1135          case(FaceElements):
1136          case(ReducedFaceElements):
1137          if (getMPISize()>1) {
1138             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1139             escriptDataC _in2 = temp.getDataC();
1140             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1141          } else {
1142             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1143          }
1144          break;
1145          case(Points):
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->Points,&_in2,&_target);
1150          } else {
1151             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1152          }
1153          break;
1154          case(ContactElementsZero):
1155          case(ContactElementsOne):
1156          case(ReducedContactElementsZero):
1157          case(ReducedContactElementsOne):
1158          if (getMPISize()>1) {
1159             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1160             escriptDataC _in2 = temp.getDataC();
1161             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1162          } else {
1163             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1164          }
1165          break;
1166          default:
1167             stringstream temp;
1168             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1169             throw FinleyAdapterException(temp.str());
1170             break;
1171          }
1172        break;        break;
    }  
    break;  
1173     default:     default:
1174        stringstream temp;        stringstream temp;
1175        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 1193  void MeshAdapter::setToX(escript::Data&
1193        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1194        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1195     } else {     } else {
1196        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1197        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1198        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1199        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1244  void MeshAdapter::interpolateACross(escr Line 1274  void MeshAdapter::interpolateACross(escr
1274  //  //
1275  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1276  //  //
1277  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1278  {  {
1279     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1280     if (argDomain!=*this)     if (argDomain!=*this)
# Line 1257  void MeshAdapter::setToIntegrals(std::ve Line 1287  void MeshAdapter::setToIntegrals(std::ve
1287     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1288     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1289     case(Nodes):     case(Nodes):
1290     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1291     _temp=temp.getDataC();     _temp=temp.getDataC();
1292     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1293     break;     break;
1294     case(ReducedNodes):     case(ReducedNodes):
1295     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1296     _temp=temp.getDataC();     _temp=temp.getDataC();
1297     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1298     break;     break;
# Line 1294  void MeshAdapter::setToIntegrals(std::ve Line 1324  void MeshAdapter::setToIntegrals(std::ve
1324     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1325     break;     break;
1326     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1327     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1328     _temp=temp.getDataC();     _temp=temp.getDataC();
1329     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1330     break;     break;
1331     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1332     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1333     _temp=temp.getDataC();     _temp=temp.getDataC();
1334     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1335     break;     break;
# Line 1331  void MeshAdapter::setToGradient(escript: Line 1361  void MeshAdapter::setToGradient(escript:
1361     escript::Data temp;     escript::Data temp;
1362     if (getMPISize()>1) {     if (getMPISize()>1) {
1363        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1364           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1365           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1366        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1367           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1368           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1369        } else {        } else {
1370           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1453  void MeshAdapter::setNewX(const escript: Line 1483  void MeshAdapter::setNewX(const escript:
1483     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1484     if (newDomain!=*this)     if (newDomain!=*this)
1485        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1486     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1487     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1488           Finley_Mesh_setCoordinates(mesh,&tmp);
1489       } else {
1490           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1491           tmp = new_x_inter.getDataC();
1492           Finley_Mesh_setCoordinates(mesh,&tmp);
1493       }
1494     checkFinleyError();     checkFinleyError();
1495  }  }
1496    
# Line 1472  void MeshAdapter::extractArgsFromDict(co Line 1508  void MeshAdapter::extractArgsFromDict(co
1508    
1509     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1510     for (int i=0; i<numData; ++i) {     for (int i=0; i<numData; ++i) {
1511        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1512        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1513        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1514           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 1522  void MeshAdapter::extractArgsFromDict(co
1522  //  //
1523  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1524  //  //
1525  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
1526  {  {
1527     int num_data;     int num_data;
1528     char **names;     char **names;
# Line 1512  void MeshAdapter::saveDX(const std::stri Line 1548  void MeshAdapter::saveDX(const std::stri
1548  //  //
1549  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1550  //  //
1551  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
1552  {  {
1553     int num_data;     int num_data;
1554     char **names;     char **names;
# Line 1520  void MeshAdapter::saveVTK(const std::str Line 1556  void MeshAdapter::saveVTK(const std::str
1556     escriptDataC **ptr_data;     escriptDataC **ptr_data;
1557    
1558     extractArgsFromDict(arg, num_data, names, data, ptr_data);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1559     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());
1560     checkFinleyError();     checkFinleyError();
1561    
1562     /* win32 refactor */     /* win32 refactor */
# Line 1533  void MeshAdapter::saveVTK(const std::str Line 1569  void MeshAdapter::saveVTK(const std::str
1569     TMPMEMFREE(names);     TMPMEMFREE(names);
1570  }  }
1571    
1572    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1573    {
1574    #ifdef ESYS_MPI
1575        index_t myFirstNode=0, myLastNode=0, k=0;
1576        index_t* globalNodeIndex=0;
1577        Finley_Mesh* mesh_p=m_finleyMesh.get();
1578        if (fs_code == FINLEY_REDUCED_NODES)
1579        {
1580        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1581        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1582        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1583        }
1584        else
1585        {
1586        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1587        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1588        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1589        }
1590        k=globalNodeIndex[id];
1591        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1592    #endif
1593        return true;
1594    }
1595    
1596    
1597    
1598  //  //
1599  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1600  //  //
1601  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1602                                                   const int row_blocksize,                                                   const int row_blocksize,
1603                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1604                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1579  SystemMatrixAdapter MeshAdapter::newSyst Line 1641  SystemMatrixAdapter MeshAdapter::newSyst
1641  #endif  #endif
1642     }     }
1643     else {     else {
1644        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1645     }     }
1646     checkPasoError();     checkPasoError();
1647     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1648     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1649       return ASM_ptr(sma);
1650    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1651  }  }
1652    
1653  //  //
1654  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1655  //  //
1656  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1657                                                           const double theta,                                                           const bool useBackwardEuler,
1658                                                           const int blocksize,                                                           const int blocksize,
1659                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1660                                                           const int type) const                                                           const int type) const
# Line 1612  TransportProblemAdapter MeshAdapter::new Line 1676  TransportProblemAdapter MeshAdapter::new
1676    
1677     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1678     checkFinleyError();     checkFinleyError();
1679     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1680     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1681     checkPasoError();     checkPasoError();
1682     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1683     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1684       return ATP_ptr(tpa);
1685    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1686  }  }
1687    
1688  //  //
# Line 1655  bool MeshAdapter::isCellOriented(int fun Line 1721  bool MeshAdapter::isCellOriented(int fun
1721     return false;     return false;
1722  }  }
1723    
1724    bool
1725    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1726    {
1727       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1728        class 1: DOF <-> Nodes
1729        class 2: ReducedDOF <-> ReducedNodes
1730        class 3: Points
1731        class 4: Elements
1732        class 5: ReducedElements
1733        class 6: FaceElements
1734        class 7: ReducedFaceElements
1735        class 8: ContactElementZero <-> ContactElementOne
1736        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1737    
1738       There is also a set of lines. Interpolation is possible down a line but not between lines.
1739       class 1 and 2 belong to all lines so aren't considered.
1740        line 0: class 3
1741        line 1: class 4,5
1742        line 2: class 6,7
1743        line 3: class 8,9
1744    
1745       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1746       eg hasnodes is true if we have at least one instance of Nodes.
1747       */
1748        if (fs.empty())
1749        {
1750            return false;
1751        }
1752        vector<int> hasclass(10);
1753        vector<int> hasline(4);
1754        bool hasnodes=false;
1755        bool hasrednodes=false;
1756        bool hascez=false;
1757        bool hasrcez=false;
1758        for (int i=0;i<fs.size();++i)
1759        {
1760        switch(fs[i])
1761        {
1762        case(Nodes):   hasnodes=true;   // no break is deliberate
1763        case(DegreesOfFreedom):
1764            hasclass[1]=1;
1765            break;
1766        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1767        case(ReducedDegreesOfFreedom):
1768            hasclass[2]=1;
1769            break;
1770        case(Points):
1771            hasline[0]=1;
1772            hasclass[3]=1;
1773            break;
1774        case(Elements):
1775            hasclass[4]=1;
1776            hasline[1]=1;
1777            break;
1778        case(ReducedElements):
1779            hasclass[5]=1;
1780            hasline[1]=1;
1781            break;
1782        case(FaceElements):
1783            hasclass[6]=1;
1784            hasline[2]=1;
1785            break;
1786        case(ReducedFaceElements):
1787            hasclass[7]=1;
1788            hasline[2]=1;
1789            break;
1790        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1791        case(ContactElementsOne):
1792            hasclass[8]=1;
1793            hasline[3]=1;
1794            break;
1795        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1796        case(ReducedContactElementsOne):
1797            hasclass[9]=1;
1798            hasline[3]=1;
1799            break;
1800        default:
1801            return false;
1802        }
1803        }
1804        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1805        // fail if we have more than one leaf group
1806    
1807        if (totlines>1)
1808        {
1809        return false;   // there are at least two branches we can't interpolate between
1810        }
1811        else if (totlines==1)
1812        {
1813        if (hasline[0]==1)      // we have points
1814        {
1815            resultcode=Points;
1816        }
1817        else if (hasline[1]==1)
1818        {
1819            if (hasclass[5]==1)
1820            {
1821            resultcode=ReducedElements;
1822            }
1823            else
1824            {
1825            resultcode=Elements;
1826            }
1827        }
1828        else if (hasline[2]==1)
1829        {
1830            if (hasclass[7]==1)
1831            {
1832            resultcode=ReducedFaceElements;
1833            }
1834            else
1835            {
1836            resultcode=FaceElements;
1837            }
1838        }
1839        else    // so we must be in line3
1840        {
1841            if (hasclass[9]==1)
1842            {
1843            // need something from class 9
1844            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1845            }
1846            else
1847            {
1848            // something from class 8
1849            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1850            }
1851        }
1852        }
1853        else    // totlines==0
1854        {
1855        if (hasclass[2]==1)
1856        {
1857            // something from class 2
1858            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1859        }
1860        else
1861        {   // something from class 1
1862            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1863        }
1864        }
1865        return true;
1866    }
1867    
1868  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1869  {  {
1870     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1871     case(Nodes):     case(Nodes):
1872     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1873     case(Nodes):      case(Nodes):
1874     case(ReducedNodes):      case(ReducedNodes):
1875     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1876     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1877     case(Elements):      case(Elements):
1878     case(ReducedElements):      case(ReducedElements):
1879     case(FaceElements):      case(FaceElements):
1880     case(ReducedFaceElements):      case(ReducedFaceElements):
1881     case(Points):      case(Points):
1882     case(ContactElementsZero):      case(ContactElementsZero):
1883     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1884     case(ContactElementsOne):      case(ContactElementsOne):
1885     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1886     return true;      return true;
1887     default:      default:
1888        stringstream temp;            stringstream temp;
1889        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;
1890        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1891     }     }
1892     break;     break;
1893     case(ReducedNodes):     case(ReducedNodes):
1894     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1895     case(ReducedNodes):      case(ReducedNodes):
1896     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1897     case(Elements):      case(Elements):
1898     case(ReducedElements):      case(ReducedElements):
1899     case(FaceElements):      case(FaceElements):
1900     case(ReducedFaceElements):      case(ReducedFaceElements):
1901     case(Points):      case(Points):
1902     case(ContactElementsZero):      case(ContactElementsZero):
1903     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1904     case(ContactElementsOne):      case(ContactElementsOne):
1905     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1906     return true;      return true;
1907     case(Nodes):      case(Nodes):
1908     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1909     return false;      return false;
1910     default:      default:
1911        stringstream temp;          stringstream temp;
1912        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;
1913        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1914     }     }
1915     break;     break;
1916     case(Elements):     case(Elements):
1917     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1918        return true;        return true;
1919     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1920        return true;        return true;
1921     } else {          } else {
1922        return false;            return false;
1923     }          }
1924     case(ReducedElements):     case(ReducedElements):
1925     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1926        return true;        return true;
1927     } else {      } else {
1928        return false;            return false;
1929     }      }
1930     case(FaceElements):     case(FaceElements):
1931     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1932        return true;              return true;
1933     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1934        return true;              return true;
1935     } else {      } else {
1936        return false;              return false;
1937     }      }
1938     case(ReducedFaceElements):     case(ReducedFaceElements):
1939     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1940        return true;              return true;
1941     } else {      } else {
1942        return false;          return false;
1943     }      }
1944     case(Points):     case(Points):
1945     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1946        return true;              return true;
1947     } else {      } else {
1948        return false;              return false;
1949     }      }
1950     case(ContactElementsZero):     case(ContactElementsZero):
1951     case(ContactElementsOne):     case(ContactElementsOne):
1952     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1953        return true;              return true;
1954     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1955        return true;              return true;
1956     } else {      } else {
1957        return false;              return false;
1958     }      }
1959     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1960     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1961     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1962        return true;              return true;
1963     } else {      } else {
1964        return false;              return false;
1965     }      }
1966     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1967     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1968     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1969     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1970     case(Nodes):      case(Nodes):
1971     case(ReducedNodes):      case(ReducedNodes):
1972     case(Elements):      case(Elements):
1973     case(ReducedElements):      case(ReducedElements):
1974     case(Points):      case(Points):
1975     case(FaceElements):      case(FaceElements):
1976     case(ReducedFaceElements):      case(ReducedFaceElements):
1977     case(ContactElementsZero):      case(ContactElementsZero):
1978     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1979     case(ContactElementsOne):      case(ContactElementsOne):
1980     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1981     return true;      return true;
1982     default:      default:
1983        stringstream temp;          stringstream temp;
1984        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;
1985        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1986     }      }
1987     break;      break;
1988     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1989     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1990     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1991     case(ReducedNodes):      case(ReducedNodes):
1992     case(Elements):      case(Elements):
1993     case(ReducedElements):      case(ReducedElements):
1994     case(FaceElements):      case(FaceElements):
1995     case(ReducedFaceElements):      case(ReducedFaceElements):
1996     case(Points):      case(Points):
1997     case(ContactElementsZero):      case(ContactElementsZero):
1998     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1999     case(ContactElementsOne):      case(ContactElementsOne):
2000     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
2001     return true;      return true;
2002     case(Nodes):      case(Nodes):
2003     case(DegreesOfFreedom):      case(DegreesOfFreedom):
2004     return false;      return false;
2005     default:      default:
2006        stringstream temp;          stringstream temp;
2007        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;
2008        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
2009     }      }
2010     break;      break;
2011     default:     default:
2012        stringstream temp;        stringstream temp;
2013        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 2048  int MeshAdapter::getSystemMatrixTypeId(c
2048  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
2049  {  {
2050     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2051     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);
2052     checkPasoError();     checkPasoError();
2053     return out;     return out;
2054  }  }
2055    
2056  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2057  {  {
2058     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2059  }  }
2060    
2061  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2062  {  {
2063     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2064  }  }
2065    
2066  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2067  {  {
2068     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2069  }  }
2070    
2071  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2072  {  {
2073     int *out = NULL;     int *out = NULL;
2074     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2017  void MeshAdapter::setTags(const int func Line 2227  void MeshAdapter::setTags(const int func
2227     return;     return;
2228  }  }
2229    
2230  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2231  {  {
2232     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2233     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 2235  void MeshAdapter::setTagMap(const std::s
2235     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2236  }  }
2237    
2238  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2239  {  {
2240     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2241     int tag=0;     int tag=0;
# Line 2035  int MeshAdapter::getTag(const std::strin Line 2245  int MeshAdapter::getTag(const std::strin
2245     return tag;     return tag;
2246  }  }
2247    
2248  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2249  {  {
2250     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2251     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2252  }  }
2253    
2254  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2255  {  {
2256     stringstream temp;     stringstream temp;
2257     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2095  int MeshAdapter::getNumberOfTagsInUse(in Line 2305  int MeshAdapter::getNumberOfTagsInUse(in
2305    }    }
2306    return numTags;    return numTags;
2307  }  }
2308  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2309    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2310  {  {
2311    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2312    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2161  bool MeshAdapter::canTag(int functionSpa Line 2372  bool MeshAdapter::canTag(int functionSpa
2372    }    }
2373  }  }
2374    
2375    AbstractDomain::StatusType MeshAdapter::getStatus() const
2376    {
2377      Finley_Mesh* mesh=m_finleyMesh.get();
2378      return Finley_Mesh_getStatus(mesh);
2379    }
2380    
2381    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2382    {
2383      
2384      Finley_Mesh* mesh=m_finleyMesh.get();
2385      int order =-1;
2386      switch(functionSpaceCode) {
2387       case(Nodes):
2388       case(DegreesOfFreedom):
2389              order=mesh->approximationOrder;
2390              break;
2391       case(ReducedNodes):
2392       case(ReducedDegreesOfFreedom):
2393              order=mesh->reducedApproximationOrder;
2394              break;
2395       case(Elements):
2396       case(FaceElements):
2397       case(Points):
2398       case(ContactElementsZero):
2399       case(ContactElementsOne):
2400              order=mesh->integrationOrder;
2401              break;
2402       case(ReducedElements):
2403       case(ReducedFaceElements):
2404       case(ReducedContactElementsZero):
2405       case(ReducedContactElementsOne):
2406              order=mesh->reducedIntegrationOrder;
2407              break;
2408       default:
2409          stringstream temp;
2410          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2411          throw FinleyAdapterException(temp.str());
2412      }
2413      return order;
2414    }
2415    
2416    bool MeshAdapter::supportsContactElements() const
2417    {
2418      return true;
2419    }
2420    
2421    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2422    {
2423      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2424    }
2425    
2426    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2427    {
2428      Finley_ReferenceElementSet_dealloc(m_refSet);
2429    }
2430    
2431  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26