/[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 1755 by ksteube, Mon Sep 8 01:34:40 2008 UTC revision 3355 by caltinay, Tue Nov 16 06:35:06 2010 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2010 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  #include "MeshAdapter.h"  #include "MeshAdapter.h"
16  #include "escript/Data.h"  #include "escript/Data.h"
# Line 20  Line 19 
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21  extern "C" {  extern "C" {
22  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
 #include <vector>  
24    
25  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  #include <boost/python/import.hpp>
26    
27  using namespace std;  using namespace std;
28  using namespace escript;  using namespace escript;
# Line 83  int MeshAdapter::getMPIRank() const Line 81  int MeshAdapter::getMPIRank() const
81  {  {
82     return m_finleyMesh.get()->MPIInfo->rank;     return m_finleyMesh.get()->MPIInfo->rank;
83  }  }
84    void MeshAdapter::MPIBarrier() const
85    {
86    #ifdef ESYS_MPI
87       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88    #endif
89       return;
90    }
91    bool MeshAdapter::onMasterProcessor() const
92    {
93       return m_finleyMesh.get()->MPIInfo->rank == 0;
94    }
95    
96    
97    #ifdef ESYS_MPI
98      MPI_Comm
99    #else
100      unsigned int
101    #endif
102    MeshAdapter::getMPIComm() const
103    {
104    #ifdef ESYS_MPI
105        return m_finleyMesh->MPIInfo->comm;
106    #else
107        return 0;
108    #endif
109    }
110    
111    
112  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
113     return m_finleyMesh.get();     return m_finleyMesh.get();
114  }  }
115    
116  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
117  {  {
118     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;
119     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 98  void MeshAdapter::write(const std::strin Line 122  void MeshAdapter::write(const std::strin
122     TMPMEMFREE(fName);     TMPMEMFREE(fName);
123  }  }
124    
125  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
126  {  {
127     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
128  }  }
129    
130  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
131  {  {
132  #ifdef USE_NETCDF  #ifdef USE_NETCDF
133     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
134     NcVar *ids;     NcVar *ids;
135     int *int_ptr;     int *int_ptr;
136     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
# Line 123  void MeshAdapter::dump(const std::string Line 147  void MeshAdapter::dump(const std::string
147     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
148     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
149     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150    #ifdef ESYS_MPI
151       MPI_Status status;
152    #endif
153    
154    /* Incoming token indicates it's my turn to write */
155    #ifdef ESYS_MPI
156       if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
157    #endif
158    
159     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
161    
162     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
163     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
164     if (tag_map) {     num_Tags = 0;
165        while (tag_map) {     while (tag_map) {
166           num_Tags++;        num_Tags++;
167           tag_map=tag_map->next;        tag_map=tag_map->next;
       }  
168     }     }
169    
170     // NetCDF error handler     // NetCDF error handler
171     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
172     // Create the file.     // Create the file.
173     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
174       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
175     // check if writing was successful     // check if writing was successful
176     if (!dataFile.is_valid())     if (!dataFile.is_valid())
177        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
178    
179     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
180       // dim_Elements only appears if > 0)
181     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
182        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
183     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
184        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
185     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)) )
186        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
187     if (num_Elements>0)     if (num_Elements>0)
188        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
189           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
190     if (num_FaceElements>0)     if (num_FaceElements>0)
191        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
192           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
193     if (num_ContactElements>0)     if (num_ContactElements>0)
194        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
195           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
196     if (num_Points>0)     if (num_Points>0)
197        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
198           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
199     if (num_Elements>0)     if (num_Elements>0)
200        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
201           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
202     if (num_FaceElements>0)     if (num_FaceElements>0)
203        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
204           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
205     if (num_ContactElements>0)     if (num_ContactElements>0)
206        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
207           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
208     if (num_Tags>0)     if (num_Tags>0)
209        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
210           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
211    
212     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
213     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
214        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
215     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
216        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
217     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
218        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
219     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
220        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
221     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
222        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
223     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
225     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
226        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
227     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
228        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
229     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
230        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
231     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
232        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
233     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
234        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
235     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
236        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
237     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
238        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
239     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
240        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
241     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
242        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
243     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
245     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
247     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
249     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
250        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
251    
252     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
253    
254     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
255       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
256          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
257       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
258       if (! (ids->put(int_ptr, mpi_size+1)) )
259          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
260    
261       // Nodes degreesOfFreedomDistribution
262       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
263          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
264       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
265       if (! (ids->put(int_ptr, mpi_size+1)) )
266          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
267    
268       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
269       // (it treats them as NC_UNLIMITED)
270     if (numNodes>0) {     if (numNodes>0) {
271    
272        // Nodes Id        // Nodes Id
273        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
274           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
275        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
276        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
277           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
278    
279        // Nodes Tag        // Nodes Tag
280        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
281           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
282        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
283        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
284           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
285    
286        // Nodes gDOF        // Nodes gDOF
287        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
288           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
289        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
290        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
291           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
292    
293        // Nodes global node index        // Nodes global node index
294        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
295           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
296        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
297        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
298           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
299    
300        // Nodes grDof        // Nodes grDof
301        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
302           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
303        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
304        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
305           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
306    
307        // Nodes grNI        // Nodes grNI
308        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
309           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
310        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
311        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
312           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
313    
314        // Nodes Coordinates        // Nodes Coordinates
315        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
316           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
317        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
318           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);  
 { int i; printf("ksteube wrote Nodes_DofDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", mesh->Nodes->degreesOfFreedomDistribution->first_component[i]); printf("\n"); }  
   
       // 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);  
 { int i; printf("ksteube wrote Nodes_NodeDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", mesh->Nodes->nodesDistribution->first_component[i]); printf("\n"); }  
319    
320     }     }
321    
# Line 291  void MeshAdapter::dump(const std::string Line 323  void MeshAdapter::dump(const std::string
323    
324     if (num_Elements>0) {     if (num_Elements>0) {
325    
       // Temp storage to gather node IDs  
       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);  
   
326        // Elements_Id        // Elements_Id
327        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
328           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
329        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
330        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
331           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
332    
333        // Elements_Tag        // Elements_Tag
334        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
335           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
336        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
337        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
338           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
339    
340        // Elements_Owner        // Elements_Owner
341        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
342           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
343        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
344        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
345           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
346    
347        // Elements_Color        // Elements_Color
348        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
349           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
350        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
351        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
352           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
353    
354        // Elements_Nodes        // Elements_Nodes
       for (int i=0; i<num_Elements; i++)  
          for (int j=0; j<num_Elements_numNodes; j++)  
             Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];  
355        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
356           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
357        if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
358           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
   
       TMPMEMFREE(Elements_Nodes);  
359    
360     }     }
361    
# Line 339  void MeshAdapter::dump(const std::string Line 363  void MeshAdapter::dump(const std::string
363    
364     if (num_FaceElements>0) {     if (num_FaceElements>0) {
365    
       // Temp storage to gather node IDs  
       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);  
   
366        // FaceElements_Id        // FaceElements_Id
367        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
368           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
369        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
370        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
371           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
372    
373        // FaceElements_Tag        // FaceElements_Tag
374        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
375           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
376        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
377        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
378           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
379    
380        // FaceElements_Owner        // FaceElements_Owner
381        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
382           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
383        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
384        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
385           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
386    
387        // FaceElements_Color        // FaceElements_Color
388        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
389           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
390        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
391        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
392           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
393    
394        // FaceElements_Nodes        // FaceElements_Nodes
       for (int i=0; i<num_FaceElements; i++)  
          for (int j=0; j<num_FaceElements_numNodes; j++)  
             FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];  
395        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
396           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
397        if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
398           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
   
       TMPMEMFREE(FaceElements_Nodes);  
399    
400     }     }
401    
# Line 387  void MeshAdapter::dump(const std::string Line 403  void MeshAdapter::dump(const std::string
403    
404     if (num_ContactElements>0) {     if (num_ContactElements>0) {
405    
       // Temp storage to gather node IDs  
       int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);  
   
406        // ContactElements_Id        // ContactElements_Id
407        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
408           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
409        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
410        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
411           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
412    
413        // ContactElements_Tag        // ContactElements_Tag
414        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
415           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
416        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
417        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
418           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
419    
420        // ContactElements_Owner        // ContactElements_Owner
421        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
422           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
423        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
424        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
425           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
426    
427        // ContactElements_Color        // ContactElements_Color
428        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
429           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
430        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
431        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
432           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
433    
434        // ContactElements_Nodes        // ContactElements_Nodes
       for (int i=0; i<num_ContactElements; i++)  
          for (int j=0; j<num_ContactElements_numNodes; j++)  
             ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];  
435        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
436           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
437        if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
438           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
   
       TMPMEMFREE(ContactElements_Nodes);  
439    
440     }     }
441    
# Line 437  void MeshAdapter::dump(const std::string Line 445  void MeshAdapter::dump(const std::string
445    
446        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
447    
       // Temp storage to gather node IDs  
       int *Points_Nodes = TMPMEMALLOC(num_Points,int);  
   
448        // Points_Id        // Points_Id
449        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
450           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
451        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
452        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
453           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
454    
455        // Points_Tag        // Points_Tag
456        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
457           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
458        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
459        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
460           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
461    
462        // Points_Owner        // Points_Owner
463        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
464           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
465        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
466        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
467           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
468    
469        // Points_Color        // Points_Color
470        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
471           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
472        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
473        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
474           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
475    
476        // Points_Nodes        // Points_Nodes
477        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
       for (int i=0; i<num_Points; i++)  
          Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];  
478        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
479           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
480        if (! (ids->put(&(Points_Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
481           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
   
       TMPMEMFREE(Points_Nodes);  
482    
483     }     }
484    
# Line 501  void MeshAdapter::dump(const std::string Line 502  void MeshAdapter::dump(const std::string
502    
503        // Tags_keys        // Tags_keys
504        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
505           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
506        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
507        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
508           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
509    
510        // Tags_names_*        // Tags_names_*
511        // 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
512        // 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
513          // manual doesn't tell how to do an array of strings
514        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
515        if (tag_map) {        if (tag_map) {
516           int i = 0;           int i = 0;
517           while (tag_map) {           while (tag_map) {
518              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
519              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
520                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
521              tag_map=tag_map->next;              tag_map=tag_map->next;
522              i++;              i++;
523           }           }
524        }        }
525    
526        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
527     }     }
528    
529    /* Send token to next MPI process so he can take his turn */
530    #ifdef ESYS_MPI
531       if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
532    #endif
533    
534     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
535    
536  #else  #else
537     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
538  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
# Line 651  int MeshAdapter::getDiracDeltaFunctionCo Line 657  int MeshAdapter::getDiracDeltaFunctionCo
657  //  //
658  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
659  {  {
660     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
661       int numDim=Finley_Mesh_getDim(mesh);
662     checkFinleyError();     checkFinleyError();
663     return numDim;     return numDim;
664  }  }
# Line 685  pair<int,int> MeshAdapter::getDataShape( Line 692  pair<int,int> MeshAdapter::getDataShape(
692     case(Elements):     case(Elements):
693     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
694        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
695        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696     }     }
697     break;     break;
698     case(ReducedElements):     case(ReducedElements):
699     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
700        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
701        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702     }     }
703     break;     break;
704     case(FaceElements):     case(FaceElements):
705     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
706        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
708     }     }
709     break;     break;
710     case(ReducedFaceElements):     case(ReducedFaceElements):
711     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
712        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
714     }     }
715     break;     break;
# Line 714  pair<int,int> MeshAdapter::getDataShape( Line 721  pair<int,int> MeshAdapter::getDataShape(
721     break;     break;
722     case(ContactElementsZero):     case(ContactElementsZero):
723     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
724        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
726     }     }
727     break;     break;
728     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
729     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
730        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
732     }     }
733     break;     break;
734     case(ContactElementsOne):     case(ContactElementsOne):
735     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
736        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
738     }     }
739     break;     break;
740     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
741     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
742        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
744     }     }
745     break;     break;
# Line 761  pair<int,int> MeshAdapter::getDataShape( Line 768  pair<int,int> MeshAdapter::getDataShape(
768  // 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:
769  //  //
770  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
771                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                                   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,
773                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
775  {  {
776       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
777       if (smat==0)
778       {
779        throw FinleyAdapterException("finley only supports its own system matrices.");
780       }
781     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
782     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
783     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 780  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792    
793     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
794    
795     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 );
796     checkFinleyError();     checkFinleyError();
797    
798     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 );
799     checkFinleyError();     checkFinleyError();
800    
801     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 );
802     checkFinleyError();     checkFinleyError();
803  }  }
804    
# Line 802  void  MeshAdapter::addPDEToLumpedSystem( Line 814  void  MeshAdapter::addPDEToLumpedSystem(
814     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
815    
816     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
817       checkFinleyError();
818      
819     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
   
820     checkFinleyError();     checkFinleyError();
821    
822  }  }
823    
824    
# Line 834  void MeshAdapter::addPDEToRHS( escript:: Line 848  void MeshAdapter::addPDEToRHS( escript::
848  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
849  //  //
850  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
851                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
852                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
853                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
854                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
855                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
856  {  {
857     DataArrayView::ShapeType shape;     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
858       if (tpa==0)
859       {
860        throw FinleyAdapterException("finley only supports its own transport problems.");
861       }
862    
863    
864       DataTypes::ShapeType shape;
865     source.expand();     source.expand();
866     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
867     escriptDataC _M=M.getDataC();     escriptDataC _M=M.getDataC();
# Line 856  void MeshAdapter::addPDEToTransportProbl Line 877  void MeshAdapter::addPDEToTransportProbl
877     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
878    
879     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
880     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
881    
882     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 );
883     checkFinleyError();     checkFinleyError();
# Line 876  void MeshAdapter::addPDEToTransportProbl Line 897  void MeshAdapter::addPDEToTransportProbl
897  //  //
898  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
899  {  {
900     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
901     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
902     if (inDomain!=*this)       if (inDomain!=*this)  
903        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
904     if (targetDomain!=*this)     if (targetDomain!=*this)
# Line 888  void MeshAdapter::interpolateOnDomain(es Line 909  void MeshAdapter::interpolateOnDomain(es
909     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
910     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
911     case(Nodes):     case(Nodes):
912     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
913     case(Nodes):        case(Nodes):
914     case(ReducedNodes):        case(ReducedNodes):
915     case(DegreesOfFreedom):        case(DegreesOfFreedom):
916     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
917     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());  
918        break;        break;
919     }        case(Elements):
920     break;        case(ReducedElements):
921     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
922     switch(target.getFunctionSpace().getTypeCode()) {        break;
923     case(Nodes):        case(FaceElements):
924     case(ReducedNodes):        case(ReducedFaceElements):
925     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
926     case(ReducedDegreesOfFreedom):        break;
927     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
928     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
929     case(Elements):        break;
930     case(ReducedElements):        case(ContactElementsZero):
931     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
932     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
933     case(FaceElements):        break;
934     case(ReducedFaceElements):        case(ContactElementsOne):
935     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
936     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
937     case(Points):        break;
938     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
939     break;           stringstream temp;
940     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
941     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
942     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
943     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());  
944        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):  
945     case(ReducedNodes):     case(ReducedNodes):
946     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
947        escript::Data temp=escript::Data(in);        case(Nodes):
948        temp.expand();        case(ReducedNodes):
949        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
950        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
951        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
952     }        break;
953     break;        case(Elements):
954     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 {  
955        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
956     }        break;
957     break;        case(FaceElements):
958     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 {  
959        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
960     }        break;
961     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
962        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
963     }        break;
964     break;        case(ContactElementsZero):
965     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 {  
966        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());  
967        break;        break;
968     }        case(ContactElementsOne):
969     break;        case(ReducedContactElementsOne):
970     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
971     switch(target.getFunctionSpace().getTypeCode()) {        break;
972     case(Nodes):        default:
973     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
974     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
975     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
976     if (getMPISize()>1) {           break;
977        escript::Data temp=escript::Data(in);        }
978        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;  
979     case(Elements):     case(Elements):
980          if (target.getFunctionSpace().getTypeCode()==Elements) {
981             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
982          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
983             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
984          } else {
985             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
986          }
987          break;
988     case(ReducedElements):     case(ReducedElements):
989     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
990        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
991        escriptDataC _in2 = temp.getDataC();        } else {
992        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
993     } else {        }
994        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
995     case(FaceElements):     case(FaceElements):
996          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
997             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
998          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
999             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1000          } else {
1001             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1002          }
1003          break;
1004     case(ReducedFaceElements):     case(ReducedFaceElements):
1005     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1006        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1007        escriptDataC _in2 = temp.getDataC();        } else {
1008        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1009     } else {        }
1010        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
1011     case(Points):     case(Points):
1012     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1013        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1014        escriptDataC _in2 = temp.getDataC();        } else {
1015        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1016     } else {        }
1017        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1018     case(ContactElementsZero):     case(ContactElementsZero):
1019     case(ContactElementsOne):     case(ContactElementsOne):
1020          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1021             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1022          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1023             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1024          } else {
1025             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1026          }
1027          break;
1028     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1029     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1030     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1031        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1032        escriptDataC _in2 = temp.getDataC();        } else {
1033        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1034     } else {        }
1035        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1036     }     case(DegreesOfFreedom):      
1037     break;        switch(target.getFunctionSpace().getTypeCode()) {
1038     default:        case(ReducedDegreesOfFreedom):
1039        stringstream temp;        case(DegreesOfFreedom):
1040        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1041        throw FinleyAdapterException(temp.str());        break;
1042      
1043          case(Nodes):
1044          case(ReducedNodes):
1045          if (getMPISize()>1) {
1046             escript::Data temp=escript::Data(in);
1047             temp.expand();
1048             escriptDataC _in2 = temp.getDataC();
1049             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1050          } else {
1051             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1052          }
1053          break;
1054          case(Elements):
1055          case(ReducedElements):
1056          if (getMPISize()>1) {
1057             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1058             escriptDataC _in2 = temp.getDataC();
1059             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1060          } else {
1061             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1062          }
1063          break;
1064          case(FaceElements):
1065          case(ReducedFaceElements):
1066          if (getMPISize()>1) {
1067             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1068             escriptDataC _in2 = temp.getDataC();
1069             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1070      
1071          } else {
1072             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1073          }
1074          break;
1075          case(Points):
1076          if (getMPISize()>1) {
1077             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1078             escriptDataC _in2 = temp.getDataC();
1079          } else {
1080             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1081          }
1082          break;
1083          case(ContactElementsZero):
1084          case(ContactElementsOne):
1085          case(ReducedContactElementsZero):
1086          case(ReducedContactElementsOne):
1087          if (getMPISize()>1) {
1088             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1089             escriptDataC _in2 = temp.getDataC();
1090             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1091          } else {
1092             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1093          }
1094          break;
1095          default:
1096             stringstream temp;
1097             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1098             throw FinleyAdapterException(temp.str());
1099             break;
1100          }
1101          break;
1102       case(ReducedDegreesOfFreedom):
1103          switch(target.getFunctionSpace().getTypeCode()) {
1104          case(Nodes):
1105          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1106          break;
1107          case(ReducedNodes):
1108          if (getMPISize()>1) {
1109             escript::Data temp=escript::Data(in);
1110             temp.expand();
1111             escriptDataC _in2 = temp.getDataC();
1112             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1113          } else {
1114             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1115          }
1116          break;
1117          case(DegreesOfFreedom):
1118          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1119          break;
1120          case(ReducedDegreesOfFreedom):
1121          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1122          break;
1123          case(Elements):
1124          case(ReducedElements):
1125          if (getMPISize()>1) {
1126             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1127             escriptDataC _in2 = temp.getDataC();
1128             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1129          } else {
1130             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1131          }
1132          break;
1133          case(FaceElements):
1134          case(ReducedFaceElements):
1135          if (getMPISize()>1) {
1136             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1137             escriptDataC _in2 = temp.getDataC();
1138             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1139          } else {
1140             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1141          }
1142          break;
1143          case(Points):
1144          if (getMPISize()>1) {
1145             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1146             escriptDataC _in2 = temp.getDataC();
1147             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1148          } else {
1149             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1150          }
1151          break;
1152          case(ContactElementsZero):
1153          case(ContactElementsOne):
1154          case(ReducedContactElementsZero):
1155          case(ReducedContactElementsOne):
1156          if (getMPISize()>1) {
1157             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1158             escriptDataC _in2 = temp.getDataC();
1159             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1160          } else {
1161             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1162          }
1163          break;
1164          default:
1165             stringstream temp;
1166             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1167             throw FinleyAdapterException(temp.str());
1168             break;
1169          }
1170        break;        break;
    }  
    break;  
1171     default:     default:
1172        stringstream temp;        stringstream temp;
1173        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 1161  void MeshAdapter::interpolateOnDomain(es Line 1182  void MeshAdapter::interpolateOnDomain(es
1182  //  //
1183  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1184  {  {
1185     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1186     if (argDomain!=*this)     if (argDomain!=*this)
1187        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1188     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1170  void MeshAdapter::setToX(escript::Data& Line 1191  void MeshAdapter::setToX(escript::Data&
1191        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1192        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1193     } else {     } else {
1194        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1195        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1196        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1197        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1184  void MeshAdapter::setToX(escript::Data& Line 1205  void MeshAdapter::setToX(escript::Data&
1205  //  //
1206  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1207  {  {
1208     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1209       const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1210     if (normalDomain!=*this)     if (normalDomain!=*this)
1211        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1212     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1239  void MeshAdapter::setToNormal(escript::D Line 1261  void MeshAdapter::setToNormal(escript::D
1261  //  //
1262  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1263  {  {
1264     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1265     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1266       if (targetDomain!=this)
1267        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1268    
1269     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1249  void MeshAdapter::interpolateACross(escr Line 1272  void MeshAdapter::interpolateACross(escr
1272  //  //
1273  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1274  //  //
1275  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1276  {  {
1277     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1278     if (argDomain!=*this)     if (argDomain!=*this)
1279        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1280    
# Line 1262  void MeshAdapter::setToIntegrals(std::ve Line 1285  void MeshAdapter::setToIntegrals(std::ve
1285     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1286     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1287     case(Nodes):     case(Nodes):
1288     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1289     _temp=temp.getDataC();     _temp=temp.getDataC();
1290     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1291     break;     break;
1292     case(ReducedNodes):     case(ReducedNodes):
1293     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1294     _temp=temp.getDataC();     _temp=temp.getDataC();
1295     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1296     break;     break;
# Line 1299  void MeshAdapter::setToIntegrals(std::ve Line 1322  void MeshAdapter::setToIntegrals(std::ve
1322     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1323     break;     break;
1324     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1325     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1326     _temp=temp.getDataC();     _temp=temp.getDataC();
1327     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1328     break;     break;
1329     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1330     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1331     _temp=temp.getDataC();     _temp=temp.getDataC();
1332     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1333     break;     break;
# Line 1323  void MeshAdapter::setToIntegrals(std::ve Line 1346  void MeshAdapter::setToIntegrals(std::ve
1346  //  //
1347  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1348  {  {
1349     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1350     if (argDomain!=*this)     if (argDomain!=*this)
1351        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1352     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1353     if (gradDomain!=*this)     if (gradDomain!=*this)
1354        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1355    
# Line 1336  void MeshAdapter::setToGradient(escript: Line 1359  void MeshAdapter::setToGradient(escript:
1359     escript::Data temp;     escript::Data temp;
1360     if (getMPISize()>1) {     if (getMPISize()>1) {
1361        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1362           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1363           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1364        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1365           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1366           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1367        } else {        } else {
1368           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1448  void MeshAdapter::setToSize(escript::Dat Line 1471  void MeshAdapter::setToSize(escript::Dat
1471     checkFinleyError();     checkFinleyError();
1472  }  }
1473    
1474  // sets the location of nodes:  //
1475    // sets the location of nodes
1476    //
1477  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1478  {  {
1479     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1480     escriptDataC tmp;     escriptDataC tmp;
1481     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1482     if (newDomain!=*this)     if (newDomain!=*this)
1483        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1484     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1485     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1486           Finley_Mesh_setCoordinates(mesh,&tmp);
1487       } else {
1488           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1489           tmp = new_x_inter.getDataC();
1490           Finley_Mesh_setCoordinates(mesh,&tmp);
1491       }
1492     checkFinleyError();     checkFinleyError();
1493  }  }
1494    
1495  // saves a data array in openDX format:  //
1496  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  // Helper for the save* methods. Extracts optional data variable names and the
1497    // corresponding pointers from python dictionary. Caller must free arrays.
1498    //
1499    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1500  {  {
1501     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     numData = boost::python::extract<int>(arg.attr("__len__")());
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
1502     /* win32 refactor */     /* win32 refactor */
1503     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1504     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1505     {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1506    
1507     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1508     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1509        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1510        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1511        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1512           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1513        data[i]=d.getDataC();        data[i] = d.getDataC();
1514        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1515        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1516           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1517     }     }
1518     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1519    
1520    //
1521    // saves mesh and optionally data arrays in openDX format
1522    //
1523    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1524    {
1525       int num_data;
1526       char **names;
1527       escriptDataC *data;
1528       escriptDataC **ptr_data;
1529    
1530       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1531       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1532     checkFinleyError();     checkFinleyError();
1533    
1534     /* win32 refactor */     /* win32 refactor */
1535     TMPMEMFREE(data);     TMPMEMFREE(data);
1536     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1537     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1538     {     {
1539        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1540     }     }
# Line 1506  void MeshAdapter::saveDX(const std::stri Line 1543  void MeshAdapter::saveDX(const std::stri
1543     return;     return;
1544  }  }
1545    
1546  // saves a data array in openVTK format:  //
1547  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1548    //
1549    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1550  {  {
1551     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1552     const int num_data=boost::python::extract<int>(arg.attr("__len__")());      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1553     /* win32 refactor */                metadata, metadata_schema, arg);
1554     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  }
1555     for(int i=0;i<num_data;i++)  
1556     {  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1557        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  {
1558     }  #ifdef ESYS_MPI
1559        index_t myFirstNode=0, myLastNode=0, k=0;
1560     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;      index_t* globalNodeIndex=0;
1561     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;      Finley_Mesh* mesh_p=m_finleyMesh.get();
1562        if (fs_code == FINLEY_REDUCED_NODES)
1563        {
1564        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1565        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1566        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1567        }
1568        else
1569        {
1570        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1571        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1572        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1573        }
1574        k=globalNodeIndex[id];
1575        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1576    #endif
1577        return true;
1578    }
1579    
    boost::python::list keys=arg.keys();  
    for (int i=0;i<num_data;++i) {  
       std::string n=boost::python::extract<std::string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  
          throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
       data[i]=d.getDataC();  
       ptr_data[i]=&(data[i]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
    }  
    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1580    
    checkFinleyError();  
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0;i<num_data;i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1581    
1582     return;  //
1583  }  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1584                                                                                                                                                                        //
1585                                                                                                                                                                        ASM_ptr MeshAdapter::newSystemMatrix(
 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
 SystemMatrixAdapter MeshAdapter::newSystemMatrix(  
1586                                                   const int row_blocksize,                                                   const int row_blocksize,
1587                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1588                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1563  SystemMatrixAdapter MeshAdapter::newSyst Line 1592  SystemMatrixAdapter MeshAdapter::newSyst
1592     int reduceRowOrder=0;     int reduceRowOrder=0;
1593     int reduceColOrder=0;     int reduceColOrder=0;
1594     // is the domain right?     // is the domain right?
1595     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1596     if (row_domain!=*this)     if (row_domain!=*this)
1597        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1598     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1599     if (col_domain!=*this)     if (col_domain!=*this)
1600        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1601     // is the function space type right     // is the function space type right
# Line 1596  SystemMatrixAdapter MeshAdapter::newSyst Line 1625  SystemMatrixAdapter MeshAdapter::newSyst
1625  #endif  #endif
1626     }     }
1627     else {     else {
1628        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1629     }     }
1630     checkPasoError();     checkPasoError();
1631     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1632     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1633       return ASM_ptr(sma);
1634    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1635  }  }
1636    
1637    //
1638  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1639  TransportProblemAdapter MeshAdapter::newTransportProblem(  //
1640                                                           const double theta,  ATP_ptr MeshAdapter::newTransportProblem(
1641                                                             const bool useBackwardEuler,
1642                                                           const int blocksize,                                                           const int blocksize,
1643                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1644                                                           const int type) const                                                           const int type) const
1645  {  {
1646     int reduceOrder=0;     int reduceOrder=0;
1647     // is the domain right?     // is the domain right?
1648     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1649     if (domain!=*this)     if (domain!=*this)
1650        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1651     // is the function space type right     // is the function space type right
# Line 1626  TransportProblemAdapter MeshAdapter::new Line 1660  TransportProblemAdapter MeshAdapter::new
1660    
1661     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1662     checkFinleyError();     checkFinleyError();
1663     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1664     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1665     checkPasoError();     checkPasoError();
1666     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1667     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1668       return ATP_ptr(tpa);
1669    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1670  }  }
1671    
1672  //  //
# Line 1669  bool MeshAdapter::isCellOriented(int fun Line 1705  bool MeshAdapter::isCellOriented(int fun
1705     return false;     return false;
1706  }  }
1707    
1708    bool
1709    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1710    {
1711       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1712        class 1: DOF <-> Nodes
1713        class 2: ReducedDOF <-> ReducedNodes
1714        class 3: Points
1715        class 4: Elements
1716        class 5: ReducedElements
1717        class 6: FaceElements
1718        class 7: ReducedFaceElements
1719        class 8: ContactElementZero <-> ContactElementOne
1720        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1721    
1722       There is also a set of lines. Interpolation is possible down a line but not between lines.
1723       class 1 and 2 belong to all lines so aren't considered.
1724        line 0: class 3
1725        line 1: class 4,5
1726        line 2: class 6,7
1727        line 3: class 8,9
1728    
1729       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1730       eg hasnodes is true if we have at least one instance of Nodes.
1731       */
1732        if (fs.empty())
1733        {
1734            return false;
1735        }
1736        vector<int> hasclass(10);
1737        vector<int> hasline(4);
1738        bool hasnodes=false;
1739        bool hasrednodes=false;
1740        bool hascez=false;
1741        bool hasrcez=false;
1742        for (int i=0;i<fs.size();++i)
1743        {
1744        switch(fs[i])
1745        {
1746        case(Nodes):   hasnodes=true;   // no break is deliberate
1747        case(DegreesOfFreedom):
1748            hasclass[1]=1;
1749            break;
1750        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1751        case(ReducedDegreesOfFreedom):
1752            hasclass[2]=1;
1753            break;
1754        case(Points):
1755            hasline[0]=1;
1756            hasclass[3]=1;
1757            break;
1758        case(Elements):
1759            hasclass[4]=1;
1760            hasline[1]=1;
1761            break;
1762        case(ReducedElements):
1763            hasclass[5]=1;
1764            hasline[1]=1;
1765            break;
1766        case(FaceElements):
1767            hasclass[6]=1;
1768            hasline[2]=1;
1769            break;
1770        case(ReducedFaceElements):
1771            hasclass[7]=1;
1772            hasline[2]=1;
1773            break;
1774        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1775        case(ContactElementsOne):
1776            hasclass[8]=1;
1777            hasline[3]=1;
1778            break;
1779        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1780        case(ReducedContactElementsOne):
1781            hasclass[9]=1;
1782            hasline[3]=1;
1783            break;
1784        default:
1785            return false;
1786        }
1787        }
1788        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1789        // fail if we have more than one leaf group
1790    
1791        if (totlines>1)
1792        {
1793        return false;   // there are at least two branches we can't interpolate between
1794        }
1795        else if (totlines==1)
1796        {
1797        if (hasline[0]==1)      // we have points
1798        {
1799            resultcode=Points;
1800        }
1801        else if (hasline[1]==1)
1802        {
1803            if (hasclass[5]==1)
1804            {
1805            resultcode=ReducedElements;
1806            }
1807            else
1808            {
1809            resultcode=Elements;
1810            }
1811        }
1812        else if (hasline[2]==1)
1813        {
1814            if (hasclass[7]==1)
1815            {
1816            resultcode=ReducedFaceElements;
1817            }
1818            else
1819            {
1820            resultcode=FaceElements;
1821            }
1822        }
1823        else    // so we must be in line3
1824        {
1825            if (hasclass[9]==1)
1826            {
1827            // need something from class 9
1828            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1829            }
1830            else
1831            {
1832            // something from class 8
1833            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1834            }
1835        }
1836        }
1837        else    // totlines==0
1838        {
1839        if (hasclass[2]==1)
1840        {
1841            // something from class 2
1842            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1843        }
1844        else
1845        {   // something from class 1
1846            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1847        }
1848        }
1849        return true;
1850    }
1851    
1852  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1853  {  {
1854     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1855     case(Nodes):     case(Nodes):
1856     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1857     case(Nodes):      case(Nodes):
1858     case(ReducedNodes):      case(ReducedNodes):
1859     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1860     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1861     case(Elements):      case(Elements):
1862     case(ReducedElements):      case(ReducedElements):
1863     case(FaceElements):      case(FaceElements):
1864     case(ReducedFaceElements):      case(ReducedFaceElements):
1865     case(Points):      case(Points):
1866     case(ContactElementsZero):      case(ContactElementsZero):
1867     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1868     case(ContactElementsOne):      case(ContactElementsOne):
1869     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1870     return true;      return true;
1871     default:      default:
1872        stringstream temp;            stringstream temp;
1873        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;
1874        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1875     }     }
1876     break;     break;
1877     case(ReducedNodes):     case(ReducedNodes):
1878     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1879     case(ReducedNodes):      case(ReducedNodes):
1880     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1881     case(Elements):      case(Elements):
1882     case(ReducedElements):      case(ReducedElements):
1883     case(FaceElements):      case(FaceElements):
1884     case(ReducedFaceElements):      case(ReducedFaceElements):
1885     case(Points):      case(Points):
1886     case(ContactElementsZero):      case(ContactElementsZero):
1887     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1888     case(ContactElementsOne):      case(ContactElementsOne):
1889     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1890     return true;      return true;
1891     case(Nodes):      case(Nodes):
1892     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1893     return false;      return false;
1894     default:      default:
1895        stringstream temp;          stringstream temp;
1896        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;
1897        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1898     }     }
1899     break;     break;
1900     case(Elements):     case(Elements):
1901     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1902        return true;        return true;
1903     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1904        return true;        return true;
1905     } else {          } else {
1906        return false;            return false;
1907     }          }
1908     case(ReducedElements):     case(ReducedElements):
1909     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1910        return true;        return true;
1911     } else {      } else {
1912        return false;            return false;
1913     }      }
1914     case(FaceElements):     case(FaceElements):
1915     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1916        return true;              return true;
1917     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1918        return true;              return true;
1919     } else {      } else {
1920        return false;              return false;
1921     }      }
1922     case(ReducedFaceElements):     case(ReducedFaceElements):
1923     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1924        return true;              return true;
1925     } else {      } else {
1926        return false;          return false;
1927     }      }
1928     case(Points):     case(Points):
1929     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1930        return true;              return true;
1931     } else {      } else {
1932        return false;              return false;
1933     }      }
1934     case(ContactElementsZero):     case(ContactElementsZero):
1935     case(ContactElementsOne):     case(ContactElementsOne):
1936     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1937        return true;              return true;
1938     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1939        return true;              return true;
1940     } else {      } else {
1941        return false;              return false;
1942     }      }
1943     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1944     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1945     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1946        return true;              return true;
1947     } else {      } else {
1948        return false;              return false;
1949     }      }
1950     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1951     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1952     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1953     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1954     case(Nodes):      case(Nodes):
1955     case(ReducedNodes):      case(ReducedNodes):
1956     case(Elements):      case(Elements):
1957     case(ReducedElements):      case(ReducedElements):
1958     case(Points):      case(Points):
1959     case(FaceElements):      case(FaceElements):
1960     case(ReducedFaceElements):      case(ReducedFaceElements):
1961     case(ContactElementsZero):      case(ContactElementsZero):
1962     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1963     case(ContactElementsOne):      case(ContactElementsOne):
1964     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1965     return true;      return true;
1966     default:      default:
1967        stringstream temp;          stringstream temp;
1968        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;
1969        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1970     }      }
1971     break;      break;
1972     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1973     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1974     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1975     case(ReducedNodes):      case(ReducedNodes):
1976     case(Elements):      case(Elements):
1977     case(ReducedElements):      case(ReducedElements):
1978     case(FaceElements):      case(FaceElements):
1979     case(ReducedFaceElements):      case(ReducedFaceElements):
1980     case(Points):      case(Points):
1981     case(ContactElementsZero):      case(ContactElementsZero):
1982     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1983     case(ContactElementsOne):      case(ContactElementsOne):
1984     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1985     return true;      return true;
1986     case(Nodes):      case(Nodes):
1987     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1988     return false;      return false;
1989     default:      default:
1990        stringstream temp;          stringstream temp;
1991        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;
1992        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1993     }      }
1994     break;      break;
1995     default:     default:
1996        stringstream temp;        stringstream temp;
1997        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 1842  bool MeshAdapter::operator!=(const Abstr Line 2022  bool MeshAdapter::operator!=(const Abstr
2022     return !(operator==(other));     return !(operator==(other));
2023  }  }
2024    
2025  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2026    {
2027       Finley_Mesh* mesh=m_finleyMesh.get();
2028       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2029       checkPasoError();
2030       return out;
2031    }
2032    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2033  {  {
2034     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2035       int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2036     checkPasoError();     checkPasoError();
2037     return out;     return out;
2038  }  }
2039    
2040  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2041  {  {
2042     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2043  }  }
2044    
2045  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2046  {  {
2047     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2048  }  }
2049    
2050  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2051  {  {
2052     return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2053  }  }
2054    
2055  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2056  {  {
2057     int *out = NULL;     int *out = NULL;
2058     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2023  void MeshAdapter::setTags(const int func Line 2211  void MeshAdapter::setTags(const int func
2211     return;     return;
2212  }  }
2213    
2214  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2215  {  {
2216     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2217     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
# Line 2031  void MeshAdapter::setTagMap(const std::s Line 2219  void MeshAdapter::setTagMap(const std::s
2219     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2220  }  }
2221    
2222  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2223  {  {
2224     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2225     int tag=0;     int tag=0;
# Line 2041  int MeshAdapter::getTag(const std::strin Line 2229  int MeshAdapter::getTag(const std::strin
2229     return tag;     return tag;
2230  }  }
2231    
2232  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2233  {  {
2234     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2235     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2236  }  }
2237    
2238  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2239  {  {
2240     stringstream temp;     stringstream temp;
2241     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2101  int MeshAdapter::getNumberOfTagsInUse(in Line 2289  int MeshAdapter::getNumberOfTagsInUse(in
2289    }    }
2290    return numTags;    return numTags;
2291  }  }
2292  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2293    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2294  {  {
2295    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2296    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2144  int* MeshAdapter::borrowListOfTagsInUse( Line 2333  int* MeshAdapter::borrowListOfTagsInUse(
2333  }  }
2334    
2335    
2336    bool MeshAdapter::canTag(int functionSpaceCode) const
2337    {
2338      switch(functionSpaceCode) {
2339       case(Nodes):
2340       case(Elements):
2341       case(ReducedElements):
2342       case(FaceElements):
2343       case(ReducedFaceElements):
2344       case(Points):
2345       case(ContactElementsZero):
2346       case(ReducedContactElementsZero):
2347       case(ContactElementsOne):
2348       case(ReducedContactElementsOne):
2349              return true;
2350       case(ReducedNodes):
2351       case(DegreesOfFreedom):
2352       case(ReducedDegreesOfFreedom):
2353          return false;
2354       default:
2355        return false;
2356      }
2357    }
2358    
2359    AbstractDomain::StatusType MeshAdapter::getStatus() const
2360    {
2361      Finley_Mesh* mesh=m_finleyMesh.get();
2362      return Finley_Mesh_getStatus(mesh);
2363    }
2364    
2365    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2366    {
2367      
2368      Finley_Mesh* mesh=m_finleyMesh.get();
2369      int order =-1;
2370      switch(functionSpaceCode) {
2371       case(Nodes):
2372       case(DegreesOfFreedom):
2373              order=mesh->approximationOrder;
2374              break;
2375       case(ReducedNodes):
2376       case(ReducedDegreesOfFreedom):
2377              order=mesh->reducedApproximationOrder;
2378              break;
2379       case(Elements):
2380       case(FaceElements):
2381       case(Points):
2382       case(ContactElementsZero):
2383       case(ContactElementsOne):
2384              order=mesh->integrationOrder;
2385              break;
2386       case(ReducedElements):
2387       case(ReducedFaceElements):
2388       case(ReducedContactElementsZero):
2389       case(ReducedContactElementsOne):
2390              order=mesh->reducedIntegrationOrder;
2391              break;
2392       default:
2393          stringstream temp;
2394          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2395          throw FinleyAdapterException(temp.str());
2396      }
2397      return order;
2398    }
2399    
2400    bool MeshAdapter::supportsContactElements() const
2401    {
2402      return true;
2403    }
2404    
2405    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2406    {
2407      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2408    }
2409    
2410    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2411    {
2412      Finley_ReferenceElementSet_dealloc(m_refSet);
2413    }
2414    
2415  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1755  
changed lines
  Added in v.3355

  ViewVC Help
Powered by ViewVC 1.1.26