/[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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
24  extern "C" {  extern "C" {
25  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
26  }  }
27    
28    #include <boost/python/import.hpp>
29    
30  using namespace std;  using namespace std;
31  using namespace escript;  using namespace escript;
32    using namespace paso;
33    
34  namespace finley {  namespace finley {
35    
# Line 85  int MeshAdapter::getMPIRank() const Line 87  int MeshAdapter::getMPIRank() const
87  }  }
88  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
89  {  {
90  #ifdef PASO_MPI  #ifdef ESYS_MPI
91     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
92  #endif  #endif
93     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 98  bool MeshAdapter::onMasterProcessor() co
98  }  }
99    
100    
101    #ifdef ESYS_MPI
102      MPI_Comm
103    #else
104      unsigned int
105    #endif
106    MeshAdapter::getMPIComm() const
107    {
108    #ifdef ESYS_MPI
109        return m_finleyMesh->MPIInfo->comm;
110    #else
111        return 0;
112    #endif
113    }
114    
115    
116  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
117     return m_finleyMesh.get();     return m_finleyMesh.get();
118  }  }
119    
120  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
121  {  {
122     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;
123     strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
# Line 114  void MeshAdapter::Print_Mesh_Info(const Line 131  void MeshAdapter::Print_Mesh_Info(const
131     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
132  }  }
133    
134  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const string& fileName) const
135  {  {
136  #ifdef USE_NETCDF  #ifdef USE_NETCDF
137     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
# Line 134  void MeshAdapter::dump(const std::string Line 151  void MeshAdapter::dump(const std::string
151     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
152     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
153     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
154  #ifdef PASO_MPI  #ifdef ESYS_MPI
155     MPI_Status status;     MPI_Status status;
156  #endif  #endif
157    
158  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
159  #ifdef PASO_MPI  #ifdef ESYS_MPI
160     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
161  #endif  #endif
162    
163     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
164                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
165    
166     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
167     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
168     num_Tags = 0;     num_Tags = 0;
169     if (tag_map) {     while (tag_map) {
170        while (tag_map) {        num_Tags++;
171           num_Tags++;        tag_map=tag_map->next;
          tag_map=tag_map->next;  
       }  
172     }     }
173    
174     // NetCDF error handler     // NetCDF error handler
175     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
176     // Create the file.     // Create the file.
177     NcFile dataFile(newFileName, NcFile::Replace);     NcFile dataFile(newFileName, NcFile::Replace);
178       string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
179     // check if writing was successful     // check if writing was successful
180     if (!dataFile.is_valid())     if (!dataFile.is_valid())
181        throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);        throw DataException(msgPrefix+"Open file for output");
182    
183     // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)     // Define dimensions (num_Elements and dim_Elements are identical,
184       // dim_Elements only appears if > 0)
185     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )     if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
186        throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numNodes)");
187     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )     if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
188        throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(numDim)");
189     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)) )
190        throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_dim(mpi_size)");
191     if (num_Elements>0)     if (num_Elements>0)
192        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )        if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
193           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements)");
194     if (num_FaceElements>0)     if (num_FaceElements>0)
195        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )        if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
196           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
197     if (num_ContactElements>0)     if (num_ContactElements>0)
198        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )        if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
199           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements)");
200     if (num_Points>0)     if (num_Points>0)
201        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )        if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
202           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Points)");
203     if (num_Elements>0)     if (num_Elements>0)
204        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )        if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
205           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
206     if (num_FaceElements>0)     if (num_FaceElements>0)
207        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )        if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
208           throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
209     if (num_ContactElements>0)     if (num_ContactElements>0)
210        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )        if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
211           throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_ContactElements_numNodes)");
212     if (num_Tags>0)     if (num_Tags>0)
213        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )        if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
214           throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_dim(dim_Tags)");
215    
216     // Attributes: MPI size, MPI rank, Name, order, reduced_order     // Attributes: MPI size, MPI rank, Name, order, reduced_order
217     if (!dataFile.add_att("mpi_size", mpi_size) )     if (!dataFile.add_att("mpi_size", mpi_size) )
218        throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_size)");
219     if (!dataFile.add_att("mpi_rank", mpi_rank) )     if (!dataFile.add_att("mpi_rank", mpi_rank) )
220        throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(mpi_rank)");
221     if (!dataFile.add_att("Name",mesh->Name) )     if (!dataFile.add_att("Name",mesh->Name) )
222        throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Name)");
223     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
224        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
225     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
226        throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(order)");
227     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
228        throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(reduced_order)");
229     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
230        throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(numNodes)");
231     if (!dataFile.add_att("num_Elements",num_Elements) )     if (!dataFile.add_att("num_Elements",num_Elements) )
232        throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements)");
233     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )     if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
234        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements)");
235     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )     if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
236        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements)");
237     if (!dataFile.add_att("num_Points",num_Points) )     if (!dataFile.add_att("num_Points",num_Points) )
238        throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Points)");
239     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )     if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
240        throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
241     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )     if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
242        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
243     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
244        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_ContactElements_numNodes)");
245     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Elements_TypeId)");
247     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
249     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
250        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(ContactElements_TypeId)");
251     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
252        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(Points_TypeId)");
253     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
254        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException(msgPrefix+"add_att(num_Tags)");
255    
256     // // // // // Nodes // // // // //     // // // // // Nodes // // // // //
257    
258     // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)     // Nodes nodeDistribution
259       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
260          throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
261       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
262       if (! (ids->put(int_ptr, mpi_size+1)) )
263          throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
264    
265       // Nodes degreesOfFreedomDistribution
266       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
267          throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
268       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
269       if (! (ids->put(int_ptr, mpi_size+1)) )
270          throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
271    
272       // Only write nodes if non-empty because NetCDF doesn't like empty arrays
273       // (it treats them as NC_UNLIMITED)
274     if (numNodes>0) {     if (numNodes>0) {
275    
276        // Nodes Id        // Nodes Id
277        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
278           throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Id)");
279        int_ptr = &mesh->Nodes->Id[0];        int_ptr = &mesh->Nodes->Id[0];
280        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
281           throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Id)");
282    
283        // Nodes Tag        // Nodes Tag
284        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
285           throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Tag)");
286        int_ptr = &mesh->Nodes->Tag[0];        int_ptr = &mesh->Nodes->Tag[0];
287        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
288           throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Tag)");
289    
290        // Nodes gDOF        // Nodes gDOF
291        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
292           throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
293        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];        int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
294        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
295           throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gDOF)");
296    
297        // Nodes global node index        // Nodes global node index
298        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
299           throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_gNI)");
300        int_ptr = &mesh->Nodes->globalNodesIndex[0];        int_ptr = &mesh->Nodes->globalNodesIndex[0];
301        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
302           throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_gNI)");
303    
304        // Nodes grDof        // Nodes grDof
305        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
306           throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
307        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];        int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
308        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
309           throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grDfI)");
310    
311        // Nodes grNI        // Nodes grNI
312        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )        if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
313           throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_grNI)");
314        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];        int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
315        if (! (ids->put(int_ptr, numNodes)) )        if (! (ids->put(int_ptr, numNodes)) )
316           throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_grNI)");
317    
318        // Nodes Coordinates        // Nodes Coordinates
319        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )        if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
320           throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
321        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )        if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
322           throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Nodes_Coordinates)");
   
       // Nodes degreesOfFreedomDistribution  
       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);  
   
       // Nodes nodeDistribution  
       if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )  
          throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);  
       int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];  
       if (! (ids->put(int_ptr, mpi_size+1)) )  
          throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);  
323    
324     }     }
325    
# Line 311  void MeshAdapter::dump(const std::string Line 329  void MeshAdapter::dump(const std::string
329    
330        // Elements_Id        // Elements_Id
331        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
332           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Id)");
333        int_ptr = &mesh->Elements->Id[0];        int_ptr = &mesh->Elements->Id[0];
334        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
335           throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Id)");
336    
337        // Elements_Tag        // Elements_Tag
338        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
339           throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Tag)");
340        int_ptr = &mesh->Elements->Tag[0];        int_ptr = &mesh->Elements->Tag[0];
341        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
342           throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Tag)");
343    
344        // Elements_Owner        // Elements_Owner
345        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
346           throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Owner)");
347        int_ptr = &mesh->Elements->Owner[0];        int_ptr = &mesh->Elements->Owner[0];
348        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
349           throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Owner)");
350    
351        // Elements_Color        // Elements_Color
352        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
353           throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Color)");
354        int_ptr = &mesh->Elements->Color[0];        int_ptr = &mesh->Elements->Color[0];
355        if (! (ids->put(int_ptr, num_Elements)) )        if (! (ids->put(int_ptr, num_Elements)) )
356           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Color)");
357    
358        // Elements_Nodes        // Elements_Nodes
359        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
360           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Elements_Nodes)");
361        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
362           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Elements_Nodes)");
363    
364     }     }
365    
# Line 351  void MeshAdapter::dump(const std::string Line 369  void MeshAdapter::dump(const std::string
369    
370        // FaceElements_Id        // FaceElements_Id
371        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
372           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Id)");
373        int_ptr = &mesh->FaceElements->Id[0];        int_ptr = &mesh->FaceElements->Id[0];
374        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
375           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Id)");
376    
377        // FaceElements_Tag        // FaceElements_Tag
378        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
379           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
380        int_ptr = &mesh->FaceElements->Tag[0];        int_ptr = &mesh->FaceElements->Tag[0];
381        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
382           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Tag)");
383    
384        // FaceElements_Owner        // FaceElements_Owner
385        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
386           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
387        int_ptr = &mesh->FaceElements->Owner[0];        int_ptr = &mesh->FaceElements->Owner[0];
388        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
389           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Owner)");
390    
391        // FaceElements_Color        // FaceElements_Color
392        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
393           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Color)");
394        int_ptr = &mesh->FaceElements->Color[0];        int_ptr = &mesh->FaceElements->Color[0];
395        if (! (ids->put(int_ptr, num_FaceElements)) )        if (! (ids->put(int_ptr, num_FaceElements)) )
396           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Color)");
397    
398        // FaceElements_Nodes        // FaceElements_Nodes
399        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
400           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
401        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
402           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(FaceElements_Nodes)");
403    
404     }     }
405    
# Line 391  void MeshAdapter::dump(const std::string Line 409  void MeshAdapter::dump(const std::string
409    
410        // ContactElements_Id        // ContactElements_Id
411        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
412           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Id)");
413        int_ptr = &mesh->ContactElements->Id[0];        int_ptr = &mesh->ContactElements->Id[0];
414        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
415           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Id)");
416    
417        // ContactElements_Tag        // ContactElements_Tag
418        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
419           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Tag)");
420        int_ptr = &mesh->ContactElements->Tag[0];        int_ptr = &mesh->ContactElements->Tag[0];
421        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
422           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Tag)");
423    
424        // ContactElements_Owner        // ContactElements_Owner
425        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
426           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Owner)");
427        int_ptr = &mesh->ContactElements->Owner[0];        int_ptr = &mesh->ContactElements->Owner[0];
428        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
429           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Owner)");
430    
431        // ContactElements_Color        // ContactElements_Color
432        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
433           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Color)");
434        int_ptr = &mesh->ContactElements->Color[0];        int_ptr = &mesh->ContactElements->Color[0];
435        if (! (ids->put(int_ptr, num_ContactElements)) )        if (! (ids->put(int_ptr, num_ContactElements)) )
436           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Color)");
437    
438        // ContactElements_Nodes        // ContactElements_Nodes
439        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
440           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(ContactElements_Nodes)");
441        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
442           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(ContactElements_Nodes)");
443    
444     }     }
445    
# Line 433  void MeshAdapter::dump(const std::string Line 451  void MeshAdapter::dump(const std::string
451    
452        // Points_Id        // Points_Id
453        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
454           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Id)");
455        int_ptr = &mesh->Points->Id[0];        int_ptr = &mesh->Points->Id[0];
456        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
457           throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Id)");
458    
459        // Points_Tag        // Points_Tag
460        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
461           throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Tag)");
462        int_ptr = &mesh->Points->Tag[0];        int_ptr = &mesh->Points->Tag[0];
463        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
464           throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Tag)");
465    
466        // Points_Owner        // Points_Owner
467        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
468           throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Owner)");
469        int_ptr = &mesh->Points->Owner[0];        int_ptr = &mesh->Points->Owner[0];
470        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
471           throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Owner)");
472    
473        // Points_Color        // Points_Color
474        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
475           throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Color)");
476        int_ptr = &mesh->Points->Color[0];        int_ptr = &mesh->Points->Color[0];
477        if (! (ids->put(int_ptr, num_Points)) )        if (! (ids->put(int_ptr, num_Points)) )
478           throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Color)");
479    
480        // Points_Nodes        // Points_Nodes
481        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
482        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
483           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Points_Nodes)");
484        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
485           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Points_Nodes)");
486    
487     }     }
488    
# Line 488  void MeshAdapter::dump(const std::string Line 506  void MeshAdapter::dump(const std::string
506    
507        // Tags_keys        // Tags_keys
508        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )        if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
509           throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);           throw DataException(msgPrefix+"add_var(Tags_keys)");
510        int_ptr = &Tags_keys[0];        int_ptr = &Tags_keys[0];
511        if (! (ids->put(int_ptr, num_Tags)) )        if (! (ids->put(int_ptr, num_Tags)) )
512           throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);           throw DataException(msgPrefix+"put(Tags_keys)");
513    
514        // Tags_names_*        // Tags_names_*
515        // 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
516        // 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
517          // manual doesn't tell how to do an array of strings
518        tag_map = mesh->TagMap;        tag_map = mesh->TagMap;
519        if (tag_map) {        if (tag_map) {
520           int i = 0;           int i = 0;
521           while (tag_map) {           while (tag_map) {
522              sprintf(name_temp, "Tags_name_%d", i);              sprintf(name_temp, "Tags_name_%d", i);
523              if (!dataFile.add_att(name_temp, tag_map->name) )              if (!dataFile.add_att(name_temp, tag_map->name) )
524                 throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);                 throw DataException(msgPrefix+"add_att(Tags_names_XX)");
525              tag_map=tag_map->next;              tag_map=tag_map->next;
526              i++;              i++;
527           }           }
528        }        }
529    
530        TMPMEMFREE(Tags_keys);        TMPMEMFREE(Tags_keys);
   
531     }     }
532    
533  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
534  #ifdef PASO_MPI  #ifdef ESYS_MPI
535     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
536  #endif  #endif
537    
# Line 633  int MeshAdapter::getReducedSolutionCode( Line 651  int MeshAdapter::getReducedSolutionCode(
651     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
652  }  }
653    
654  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
655  {  {
656     return Points;     return Points;
657  }  }
# Line 678  pair<int,int> MeshAdapter::getDataShape( Line 696  pair<int,int> MeshAdapter::getDataShape(
696     case(Elements):     case(Elements):
697     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
698        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
699        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
700     }     }
701     break;     break;
702     case(ReducedElements):     case(ReducedElements):
703     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
704        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
705        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
706     }     }
707     break;     break;
708     case(FaceElements):     case(FaceElements):
709     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
710        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
711        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
712     }     }
713     break;     break;
714     case(ReducedFaceElements):     case(ReducedFaceElements):
715     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
716        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
717        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
718     }     }
719     break;     break;
# Line 707  pair<int,int> MeshAdapter::getDataShape( Line 725  pair<int,int> MeshAdapter::getDataShape(
725     break;     break;
726     case(ContactElementsZero):     case(ContactElementsZero):
727     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
728        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
729        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
730     }     }
731     break;     break;
732     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
733     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
734        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
735        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
736     }     }
737     break;     break;
738     case(ContactElementsOne):     case(ContactElementsOne):
739     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
740        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
741        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
742     }     }
743     break;     break;
744     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
745     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
746        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
747        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
748     }     }
749     break;     break;
# Line 754  pair<int,int> MeshAdapter::getDataShape( Line 772  pair<int,int> MeshAdapter::getDataShape(
772  // 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:
773  //  //
774  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
775                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
776                                   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,
777                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
778                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
779                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
780  {  {
781       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
782       if (smat==0)
783       {
784        throw FinleyAdapterException("finley only supports Paso system matrices.");
785       }
786     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
787     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
788     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 770  void MeshAdapter::addPDEToSystem( Line 794  void MeshAdapter::addPDEToSystem(
794     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
795     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
796     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
797       escriptDataC _d_dirac=d_dirac.getDataC();
798       escriptDataC _y_dirac=y_dirac.getDataC();
799    
800     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
801    
802     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 );
803       checkFinleyError();
804    
805       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
806     checkFinleyError();     checkFinleyError();
807    
808     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
809     checkFinleyError();     checkFinleyError();
810    
811     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->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
812     checkFinleyError();     checkFinleyError();
813  }  }
814    
815  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
816                                          escript::Data& mat,                                          escript::Data& mat,
817                                          const escript::Data& D,                                          const escript::Data& D,
818                                          const escript::Data& d) const                                          const escript::Data& d,
819                                            const escript::Data& d_dirac,
820                                            const bool useHRZ) const
821  {  {
822     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
823     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
824     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
825       escriptDataC _d_dirac=d_dirac.getDataC();
826    
827     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
828    
829     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
830     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
831      
832       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
833       checkFinleyError();
834    
835       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
836     checkFinleyError();     checkFinleyError();
837    
838  }  }
839    
840    
841  //  //
842  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
843  //  //
844  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact,  const escript::Data& y_dirac) const
845  {  {
846     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
847    
# Line 813  void MeshAdapter::addPDEToRHS( escript:: Line 850  void MeshAdapter::addPDEToRHS( escript::
850     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
851     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
852     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
853       escriptDataC _y_dirac=y_dirac.getDataC();
854    
855     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
856     checkFinleyError();     checkFinleyError();
# Line 822  void MeshAdapter::addPDEToRHS( escript:: Line 860  void MeshAdapter::addPDEToRHS( escript::
860    
861     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
862     checkFinleyError();     checkFinleyError();
863    
864       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
865       checkFinleyError();
866    
867  }  }
868  //  //
869  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
870  //  //
871  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
872                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
873                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
874                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
875                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
876                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
877                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
878  {  {
879       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
880       if (tpa==0)
881       {
882        throw FinleyAdapterException("finley only supports Paso transport problems.");
883       }
884    
885    
886     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
887     source.expand();     source.expand();
888     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 847  void MeshAdapter::addPDEToTransportProbl Line 897  void MeshAdapter::addPDEToTransportProbl
897     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
898     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
899     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
900       escriptDataC _d_dirac=d_dirac.getDataC();
901       escriptDataC _y_dirac=y_dirac.getDataC();
902    
903    
904     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
905     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
906    
907     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 );
908     checkFinleyError();     checkFinleyError();
# Line 862  void MeshAdapter::addPDEToTransportProbl Line 915  void MeshAdapter::addPDEToTransportProbl
915    
916     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
917     checkFinleyError();     checkFinleyError();
918    
919       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
920       checkFinleyError();
921    
922  }  }
923    
924  //  //
# Line 1026  void MeshAdapter::interpolateOnDomain(es Line 1083  void MeshAdapter::interpolateOnDomain(es
1083        case(Elements):        case(Elements):
1084        case(ReducedElements):        case(ReducedElements):
1085        if (getMPISize()>1) {        if (getMPISize()>1) {
1086           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1087           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1088           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1089        } else {        } else {
# Line 1036  void MeshAdapter::interpolateOnDomain(es Line 1093  void MeshAdapter::interpolateOnDomain(es
1093        case(FaceElements):        case(FaceElements):
1094        case(ReducedFaceElements):        case(ReducedFaceElements):
1095        if (getMPISize()>1) {        if (getMPISize()>1) {
1096           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1097           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1098           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1099        
# Line 1046  void MeshAdapter::interpolateOnDomain(es Line 1103  void MeshAdapter::interpolateOnDomain(es
1103        break;        break;
1104        case(Points):        case(Points):
1105        if (getMPISize()>1) {        if (getMPISize()>1) {
1106           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1107           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1108        } else {        } else {
1109           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1110        }        }
# Line 1057  void MeshAdapter::interpolateOnDomain(es Line 1114  void MeshAdapter::interpolateOnDomain(es
1114        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1115        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1116        if (getMPISize()>1) {        if (getMPISize()>1) {
1117           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1118           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1119           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1120        } else {        } else {
# Line 1095  void MeshAdapter::interpolateOnDomain(es Line 1152  void MeshAdapter::interpolateOnDomain(es
1152        case(Elements):        case(Elements):
1153        case(ReducedElements):        case(ReducedElements):
1154        if (getMPISize()>1) {        if (getMPISize()>1) {
1155           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1156           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1157           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1158        } else {        } else {
# Line 1105  void MeshAdapter::interpolateOnDomain(es Line 1162  void MeshAdapter::interpolateOnDomain(es
1162        case(FaceElements):        case(FaceElements):
1163        case(ReducedFaceElements):        case(ReducedFaceElements):
1164        if (getMPISize()>1) {        if (getMPISize()>1) {
1165           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1166           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1167           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1168        } else {        } else {
# Line 1114  void MeshAdapter::interpolateOnDomain(es Line 1171  void MeshAdapter::interpolateOnDomain(es
1171        break;        break;
1172        case(Points):        case(Points):
1173        if (getMPISize()>1) {        if (getMPISize()>1) {
1174           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1175           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1176           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1177        } else {        } else {
# Line 1126  void MeshAdapter::interpolateOnDomain(es Line 1183  void MeshAdapter::interpolateOnDomain(es
1183        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1184        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1185        if (getMPISize()>1) {        if (getMPISize()>1) {
1186           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1187           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1188           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1189        } else {        } else {
# Line 1163  void MeshAdapter::setToX(escript::Data& Line 1220  void MeshAdapter::setToX(escript::Data&
1220        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1221        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1222     } else {     } else {
1223        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1224        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1225        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1226        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1244  void MeshAdapter::interpolateACross(escr Line 1301  void MeshAdapter::interpolateACross(escr
1301  //  //
1302  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1303  //  //
1304  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1305  {  {
1306     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1307     if (argDomain!=*this)     if (argDomain!=*this)
# Line 1257  void MeshAdapter::setToIntegrals(std::ve Line 1314  void MeshAdapter::setToIntegrals(std::ve
1314     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1315     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1316     case(Nodes):     case(Nodes):
1317     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1318     _temp=temp.getDataC();     _temp=temp.getDataC();
1319     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1320     break;     break;
1321     case(ReducedNodes):     case(ReducedNodes):
1322     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1323     _temp=temp.getDataC();     _temp=temp.getDataC();
1324     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1325     break;     break;
# Line 1294  void MeshAdapter::setToIntegrals(std::ve Line 1351  void MeshAdapter::setToIntegrals(std::ve
1351     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1352     break;     break;
1353     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1354     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1355     _temp=temp.getDataC();     _temp=temp.getDataC();
1356     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1357     break;     break;
1358     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1359     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1360     _temp=temp.getDataC();     _temp=temp.getDataC();
1361     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1362     break;     break;
# Line 1331  void MeshAdapter::setToGradient(escript: Line 1388  void MeshAdapter::setToGradient(escript:
1388     escript::Data temp;     escript::Data temp;
1389     if (getMPISize()>1) {     if (getMPISize()>1) {
1390        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1391           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1392           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1393        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1394           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1395           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1396        } else {        } else {
1397           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1453  void MeshAdapter::setNewX(const escript: Line 1510  void MeshAdapter::setNewX(const escript:
1510     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1511     if (newDomain!=*this)     if (newDomain!=*this)
1512        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1513     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1514         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1515         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1516     } else {     } else {
1517         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1518    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1519         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1520         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);*/
1521     }     }
1522     checkFinleyError();     checkFinleyError();
1523  }  }
# Line 1478  void MeshAdapter::extractArgsFromDict(co Line 1536  void MeshAdapter::extractArgsFromDict(co
1536    
1537     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1538     for (int i=0; i<numData; ++i) {     for (int i=0; i<numData; ++i) {
1539        std::string n=boost::python::extract<std::string>(keys[i]);        string n=boost::python::extract<string>(keys[i]);
1540        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1541        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1542           throw FinleyAdapterException("Error: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
# Line 1492  void MeshAdapter::extractArgsFromDict(co Line 1550  void MeshAdapter::extractArgsFromDict(co
1550  //  //
1551  // saves mesh and optionally data arrays in openDX format  // saves mesh and optionally data arrays in openDX format
1552  //  //
1553  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1554  {  {
1555     int num_data;     int num_data;
1556     char **names;     char **names;
# Line 1518  void MeshAdapter::saveDX(const std::stri Line 1576  void MeshAdapter::saveDX(const std::stri
1576  //  //
1577  // saves mesh and optionally data arrays in VTK format  // saves mesh and optionally data arrays in VTK format
1578  //  //
1579  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1580  {  {
1581     int num_data;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1582     char **names;      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1583     escriptDataC *data;                metadata, metadata_schema, arg);
1584     escriptDataC **ptr_data;  }
1585    
1586    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1587    {
1588    #ifdef ESYS_MPI
1589        index_t myFirstNode=0, myLastNode=0, k=0;
1590        index_t* globalNodeIndex=0;
1591        Finley_Mesh* mesh_p=m_finleyMesh.get();
1592        if (fs_code == FINLEY_REDUCED_NODES)
1593        {
1594        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1595        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1596        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1597        }
1598        else
1599        {
1600        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1601        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1602        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1603        }
1604        k=globalNodeIndex[id];
1605        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1606    #endif
1607        return true;
1608    }
1609    
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());  
    checkFinleyError();  
1610    
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
 }  
1611    
1612  //  //
1613  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1614  //  //
1615  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1616                                                   const int row_blocksize,                                                   const int row_blocksize,
1617                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1618                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1585  SystemMatrixAdapter MeshAdapter::newSyst Line 1655  SystemMatrixAdapter MeshAdapter::newSyst
1655  #endif  #endif
1656     }     }
1657     else {     else {
1658        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1659     }     }
1660     checkPasoError();     checkPasoError();
1661     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1662     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1663       return ASM_ptr(sma);
1664    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1665  }  }
1666    
1667  //  //
1668  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1669  //  //
1670  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const double theta,  
1671                                                           const int blocksize,                                                           const int blocksize,
1672                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1673                                                           const int type) const                                                           const int type) const
# Line 1618  TransportProblemAdapter MeshAdapter::new Line 1689  TransportProblemAdapter MeshAdapter::new
1689    
1690     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1691     checkFinleyError();     checkFinleyError();
1692     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1693     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1694     checkPasoError();     checkPasoError();
1695     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1696     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1697       return ATP_ptr(tpa);
1698    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1699  }  }
1700    
1701  //  //
# Line 1661  bool MeshAdapter::isCellOriented(int fun Line 1734  bool MeshAdapter::isCellOriented(int fun
1734     return false;     return false;
1735  }  }
1736    
1737    bool
1738    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1739    {
1740       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1741        class 1: DOF <-> Nodes
1742        class 2: ReducedDOF <-> ReducedNodes
1743        class 3: Points
1744        class 4: Elements
1745        class 5: ReducedElements
1746        class 6: FaceElements
1747        class 7: ReducedFaceElements
1748        class 8: ContactElementZero <-> ContactElementOne
1749        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1750    
1751       There is also a set of lines. Interpolation is possible down a line but not between lines.
1752       class 1 and 2 belong to all lines so aren't considered.
1753        line 0: class 3
1754        line 1: class 4,5
1755        line 2: class 6,7
1756        line 3: class 8,9
1757    
1758       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1759       eg hasnodes is true if we have at least one instance of Nodes.
1760       */
1761        if (fs.empty())
1762        {
1763            return false;
1764        }
1765        vector<int> hasclass(10);
1766        vector<int> hasline(4);
1767        bool hasnodes=false;
1768        bool hasrednodes=false;
1769        bool hascez=false;
1770        bool hasrcez=false;
1771        for (int i=0;i<fs.size();++i)
1772        {
1773        switch(fs[i])
1774        {
1775        case(Nodes):   hasnodes=true;   // no break is deliberate
1776        case(DegreesOfFreedom):
1777            hasclass[1]=1;
1778            break;
1779        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1780        case(ReducedDegreesOfFreedom):
1781            hasclass[2]=1;
1782            break;
1783        case(Points):
1784            hasline[0]=1;
1785            hasclass[3]=1;
1786            break;
1787        case(Elements):
1788            hasclass[4]=1;
1789            hasline[1]=1;
1790            break;
1791        case(ReducedElements):
1792            hasclass[5]=1;
1793            hasline[1]=1;
1794            break;
1795        case(FaceElements):
1796            hasclass[6]=1;
1797            hasline[2]=1;
1798            break;
1799        case(ReducedFaceElements):
1800            hasclass[7]=1;
1801            hasline[2]=1;
1802            break;
1803        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1804        case(ContactElementsOne):
1805            hasclass[8]=1;
1806            hasline[3]=1;
1807            break;
1808        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1809        case(ReducedContactElementsOne):
1810            hasclass[9]=1;
1811            hasline[3]=1;
1812            break;
1813        default:
1814            return false;
1815        }
1816        }
1817        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1818        // fail if we have more than one leaf group
1819    
1820        if (totlines>1)
1821        {
1822        return false;   // there are at least two branches we can't interpolate between
1823        }
1824        else if (totlines==1)
1825        {
1826        if (hasline[0]==1)      // we have points
1827        {
1828            resultcode=Points;
1829        }
1830        else if (hasline[1]==1)
1831        {
1832            if (hasclass[5]==1)
1833            {
1834            resultcode=ReducedElements;
1835            }
1836            else
1837            {
1838            resultcode=Elements;
1839            }
1840        }
1841        else if (hasline[2]==1)
1842        {
1843            if (hasclass[7]==1)
1844            {
1845            resultcode=ReducedFaceElements;
1846            }
1847            else
1848            {
1849            resultcode=FaceElements;
1850            }
1851        }
1852        else    // so we must be in line3
1853        {
1854            if (hasclass[9]==1)
1855            {
1856            // need something from class 9
1857            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1858            }
1859            else
1860            {
1861            // something from class 8
1862            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1863            }
1864        }
1865        }
1866        else    // totlines==0
1867        {
1868        if (hasclass[2]==1)
1869        {
1870            // something from class 2
1871            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1872        }
1873        else
1874        {   // something from class 1
1875            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1876        }
1877        }
1878        return true;
1879    }
1880    
1881  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1882  {  {
1883     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1884     case(Nodes):     case(Nodes):
1885     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1886     case(Nodes):      case(Nodes):
1887     case(ReducedNodes):      case(ReducedNodes):
1888     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1889     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1890     case(Elements):      case(Elements):
1891     case(ReducedElements):      case(ReducedElements):
1892     case(FaceElements):      case(FaceElements):
1893     case(ReducedFaceElements):      case(ReducedFaceElements):
1894     case(Points):      case(Points):
1895     case(ContactElementsZero):      case(ContactElementsZero):
1896     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1897     case(ContactElementsOne):      case(ContactElementsOne):
1898     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1899     return true;      return true;
1900     default:      default:
1901        stringstream temp;            stringstream temp;
1902        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;
1903        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1904     }     }
1905     break;     break;
1906     case(ReducedNodes):     case(ReducedNodes):
1907     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1908     case(ReducedNodes):      case(ReducedNodes):
1909     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1910     case(Elements):      case(Elements):
1911     case(ReducedElements):      case(ReducedElements):
1912     case(FaceElements):      case(FaceElements):
1913     case(ReducedFaceElements):      case(ReducedFaceElements):
1914     case(Points):      case(Points):
1915     case(ContactElementsZero):      case(ContactElementsZero):
1916     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1917     case(ContactElementsOne):      case(ContactElementsOne):
1918     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1919     return true;      return true;
1920     case(Nodes):      case(Nodes):
1921     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1922     return false;      return false;
1923     default:      default:
1924        stringstream temp;          stringstream temp;
1925        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;
1926        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1927     }     }
1928     break;     break;
1929     case(Elements):     case(Elements):
1930     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1931        return true;        return true;
1932     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1933        return true;        return true;
1934     } else {          } else {
1935        return false;            return false;
1936     }          }
1937     case(ReducedElements):     case(ReducedElements):
1938     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1939        return true;        return true;
1940     } else {      } else {
1941        return false;            return false;
1942     }      }
1943     case(FaceElements):     case(FaceElements):
1944     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1945        return true;              return true;
1946     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1947        return true;              return true;
1948     } else {      } else {
1949        return false;              return false;
1950     }      }
1951     case(ReducedFaceElements):     case(ReducedFaceElements):
1952     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1953        return true;              return true;
1954     } else {      } else {
1955        return false;          return false;
1956     }      }
1957     case(Points):     case(Points):
1958     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1959        return true;              return true;
1960     } else {      } else {
1961        return false;              return false;
1962     }      }
1963     case(ContactElementsZero):     case(ContactElementsZero):
1964     case(ContactElementsOne):     case(ContactElementsOne):
1965     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1966        return true;              return true;
1967     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1968        return true;              return true;
1969     } else {      } else {
1970        return false;              return false;
1971     }      }
1972     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1973     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1974     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1975        return true;              return true;
1976     } else {      } else {
1977        return false;              return false;
1978     }      }
1979     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1980     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1981     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1982     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1983     case(Nodes):      case(Nodes):
1984     case(ReducedNodes):      case(ReducedNodes):
1985     case(Elements):      case(Elements):
1986     case(ReducedElements):      case(ReducedElements):
1987     case(Points):      case(Points):
1988     case(FaceElements):      case(FaceElements):
1989     case(ReducedFaceElements):      case(ReducedFaceElements):
1990     case(ContactElementsZero):      case(ContactElementsZero):
1991     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1992     case(ContactElementsOne):      case(ContactElementsOne):
1993     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1994     return true;      return true;
1995     default:      default:
1996        stringstream temp;          stringstream temp;
1997        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;
1998        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1999     }      }
2000     break;      break;
2001     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2002     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
2003     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
2004     case(ReducedNodes):      case(ReducedNodes):
2005     case(Elements):      case(Elements):
2006     case(ReducedElements):      case(ReducedElements):
2007     case(FaceElements):      case(FaceElements):
2008     case(ReducedFaceElements):      case(ReducedFaceElements):
2009     case(Points):      case(Points):
2010     case(ContactElementsZero):      case(ContactElementsZero):
2011     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
2012     case(ContactElementsOne):      case(ContactElementsOne):
2013     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
2014     return true;      return true;
2015     case(Nodes):      case(Nodes):
2016     case(DegreesOfFreedom):      case(DegreesOfFreedom):
2017     return false;      return false;
2018     default:      default:
2019        stringstream temp;          stringstream temp;
2020        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;
2021        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
2022     }      }
2023     break;      break;
2024     default:     default:
2025        stringstream temp;        stringstream temp;
2026        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 1837  bool MeshAdapter::operator!=(const Abstr Line 2054  bool MeshAdapter::operator!=(const Abstr
2054  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2055  {  {
2056     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2057     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2058     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2059  }  }
2060    
2061  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2062  {  {
2063     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2064     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2065     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2066  }  }
2067    
2068  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2069  {  {
2070     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2071  }  }
2072    
2073  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2074  {  {
2075     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2076  }  }
2077    
2078  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2079  {  {
2080     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2081  }  }
2082    
2083  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2023  void MeshAdapter::setTags(const int func Line 2239  void MeshAdapter::setTags(const int func
2239     return;     return;
2240  }  }
2241    
2242  void MeshAdapter::setTagMap(const std::string& name,  int tag)  void MeshAdapter::setTagMap(const string& name,  int tag)
2243  {  {
2244     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2245     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2246     checkPasoError();     checkFinleyError();
2247     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2248  }  }
2249    
2250  int MeshAdapter::getTag(const std::string& name) const  int MeshAdapter::getTag(const string& name) const
2251  {  {
2252     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2253     int tag=0;     int tag=0;
2254     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2255     checkPasoError();     checkFinleyError();
2256     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2257     return tag;     return tag;
2258  }  }
2259    
2260  bool MeshAdapter::isValidTagName(const std::string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2261  {  {
2262     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2263     return Finley_Mesh_isValidTagName(mesh,name.c_str());     return Finley_Mesh_isValidTagName(mesh,name.c_str());
2264  }  }
2265    
2266  std::string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2267  {  {
2268     stringstream temp;     stringstream temp;
2269     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2174  AbstractDomain::StatusType MeshAdapter:: Line 2390  AbstractDomain::StatusType MeshAdapter::
2390    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2391  }  }
2392    
2393    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2394    {
2395      
2396      Finley_Mesh* mesh=m_finleyMesh.get();
2397      int order =-1;
2398      switch(functionSpaceCode) {
2399       case(Nodes):
2400       case(DegreesOfFreedom):
2401              order=mesh->approximationOrder;
2402              break;
2403       case(ReducedNodes):
2404       case(ReducedDegreesOfFreedom):
2405              order=mesh->reducedApproximationOrder;
2406              break;
2407       case(Elements):
2408       case(FaceElements):
2409       case(Points):
2410       case(ContactElementsZero):
2411       case(ContactElementsOne):
2412              order=mesh->integrationOrder;
2413              break;
2414       case(ReducedElements):
2415       case(ReducedFaceElements):
2416       case(ReducedContactElementsZero):
2417       case(ReducedContactElementsOne):
2418              order=mesh->reducedIntegrationOrder;
2419              break;
2420       default:
2421          stringstream temp;
2422          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2423          throw FinleyAdapterException(temp.str());
2424      }
2425      return order;
2426    }
2427    
2428    bool MeshAdapter::supportsContactElements() const
2429    {
2430      return true;
2431    }
2432    
2433    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2434    {
2435      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2436    }
2437    
2438    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2439    {
2440      Finley_ReferenceElementSet_dealloc(m_refSet);
2441    }
2442    
2443    // points will be flattened
2444    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2445    {
2446          const int dim = getDim();
2447          int numPoints=points.size()/dim;
2448          int numTags=tags.size();
2449          Finley_Mesh* mesh=m_finleyMesh.get();
2450          
2451          if ( points.size() % dim != 0 )
2452          {
2453        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2454          }
2455          
2456          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2457         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2458          
2459          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2460          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2461          
2462          for (int i=0;i<numPoints;++i) {
2463           points_ptr[ i * dim     ] = points[i * dim ];
2464           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2465           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2466               tags_ptr[i]=tags[i];
2467          }
2468          
2469          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2470          checkFinleyError();
2471          
2472          TMPMEMFREE(points_ptr);
2473          TMPMEMFREE(tags_ptr);
2474    }
2475    
2476    
2477    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2478    // {
2479    //       const int dim = getDim();
2480    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2481    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2482    //       Finley_Mesh* mesh=m_finleyMesh.get();
2483    //      
2484    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2485    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2486    //      
2487    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2488    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2489    //      
2490    //       for (int i=0;i<numPoints;++i) {
2491    //     int tag_id=-1;
2492    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2493    //     if  ( numComps !=   dim ) {
2494    //                stringstream temp;            
2495    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2496    //                throw FinleyAdapterException(temp.str());
2497    //     }
2498    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2499    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2500    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2501    //    
2502    //     if ( numTags > 0) {
2503    //            boost::python::extract<string> ex_str(tags[i]);
2504    //        if  ( ex_str.check() ) {
2505    //            tag_id=getTag( ex_str());
2506    //        } else {
2507    //             boost::python::extract<int> ex_int(tags[i]);
2508    //             if ( ex_int.check() ) {
2509    //                 tag_id=ex_int();
2510    //             } else {
2511    //              stringstream temp;          
2512    //                  temp << "Error - unable to extract tag for point " << i;
2513    //              throw FinleyAdapterException(temp.str());
2514    //            }
2515    //        }
2516    //     }      
2517    //            tags_ptr[i]=tag_id;
2518    //       }
2519    //      
2520    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2521    //       checkPasoError();
2522    //      
2523    //       TMPMEMFREE(points_ptr);
2524    //       TMPMEMFREE(tags_ptr);
2525    // }
2526    
2527    /*
2528    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2529    {  
2530        boost::python::list points =  boost::python::list();
2531        boost::python::list tags =  boost::python::list();
2532        points.append(point);
2533        tags.append(tag);
2534        addDiracPoints(points, tags);
2535    }
2536    */
2537    
2538    /*
2539    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2540    {
2541            boost::python::list points =   boost::python::list();
2542            boost::python::list tags =   boost::python::list();
2543            points.append(point);
2544            tags.append(tag);
2545            addDiracPoints(points, tags);
2546    }
2547    */
2548  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2548  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26