/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

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

revision 1116 by gross, Mon Apr 23 03:36:34 2007 UTC revision 2635 by jfenwick, Thu Aug 27 04:54:41 2009 UTC
# Line 1  Line 1 
 // $Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2  #ifdef PASO_MPI  /*******************************************************
3  #include <mpi.h>  *
4  #endif  * Copyright (c) 2003-2009 by University of Queensland
5  #include "MeshAdapter.h"  * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15    #include "MeshAdapter.h"
16  #include "escript/Data.h"  #include "escript/Data.h"
17  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
18    #ifdef USE_NETCDF
19    #include <netcdfcpp.h>
20    #endif
21    #ifdef PASO_MPI
22    #include <mpi.h>
23    #include "paso/Paso_MPI.h"
24    #endif
25    extern "C" {
26    #include "esysUtils/blocktimer.h"
27    }
28    
29  using namespace std;  using namespace std;
30  using namespace escript;  using namespace escript;
# Line 45  const int MeshAdapter::ReducedContactEle Line 50  const int MeshAdapter::ReducedContactEle
50    
51  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
52  {  {
53    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
54    //     //
55    // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer     // need to use a null_deleter as Finley_Mesh_free deletes the pointer
56    // for us.     // for us.
57    m_finleyMesh.reset(finleyMesh,null_deleter());     m_finleyMesh.reset(finleyMesh,null_deleter());
58  }  }
59    
60  //  //
# Line 57  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 62  MeshAdapter::MeshAdapter(Finley_Mesh* fi
62  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
63  m_finleyMesh(in.m_finleyMesh)  m_finleyMesh(in.m_finleyMesh)
64  {  {
65    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
66  }  }
67    
68  MeshAdapter::~MeshAdapter()  MeshAdapter::~MeshAdapter()
69  {  {
70    //     //
71    // I hope the case for the pointer being zero has been taken care of.     // I hope the case for the pointer being zero has been taken care of.
72    //  cout << "In MeshAdapter destructor." << endl;     //  cout << "In MeshAdapter destructor." << endl;
73    if (m_finleyMesh.unique()) {     if (m_finleyMesh.unique()) {
74      //   cout << "Calling dealloc." << endl;        Finley_Mesh_free(m_finleyMesh.get());
75      Finley_Mesh_dealloc(m_finleyMesh.get());     }
76      //   cout << "Finished dealloc." << endl;  }
77    }  
78    int MeshAdapter::getMPISize() const
79    {
80       return m_finleyMesh.get()->MPIInfo->size;
81    }
82    int MeshAdapter::getMPIRank() const
83    {
84       return m_finleyMesh.get()->MPIInfo->rank;
85    }
86    void MeshAdapter::MPIBarrier() const
87    {
88    #ifdef PASO_MPI
89       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
90    #endif
91       return;
92    }
93    bool MeshAdapter::onMasterProcessor() const
94    {
95       return m_finleyMesh.get()->MPIInfo->rank == 0;
96  }  }
97    
98    
99  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
100     return m_finleyMesh.get();     return m_finleyMesh.get();
101  }  }
102    
103  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
104  {  {
105    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;
106    strcpy(fName,fileName.c_str());     strcpy(fName,fileName.c_str());
107    Finley_Mesh_write(m_finleyMesh.get(),fName);     Finley_Mesh_write(m_finleyMesh.get(),fName);
108    checkFinleyError();     checkFinleyError();
109    TMPMEMFREE(fName);     TMPMEMFREE(fName);
110    }
111    
112    void MeshAdapter::Print_Mesh_Info(const bool full) const
113    {
114       Finley_PrintMesh_Info(m_finleyMesh.get(), full);
115    }
116    
117    void MeshAdapter::dump(const std::string& fileName) const
118    {
119    #ifdef USE_NETCDF
120       const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
121       NcVar *ids;
122       int *int_ptr;
123       Finley_Mesh *mesh = m_finleyMesh.get();
124       Finley_TagMap* tag_map;
125       int num_Tags = 0;
126       int mpi_size             = mesh->MPIInfo->size;
127       int mpi_rank             = mesh->MPIInfo->rank;
128       int numDim               = mesh->Nodes->numDim;
129       int numNodes             = mesh->Nodes->numNodes;
130       int num_Elements         = mesh->Elements->numElements;
131       int num_FaceElements         = mesh->FaceElements->numElements;
132       int num_ContactElements      = mesh->ContactElements->numElements;
133       int num_Points           = mesh->Points->numElements;
134       int num_Elements_numNodes        = mesh->Elements->numNodes;
135       int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
136       int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
137    #ifdef PASO_MPI
138       MPI_Status status;
139    #endif
140    
141    /* Incoming token indicates it's my turn to write */
142    #ifdef PASO_MPI
143       if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
144    #endif
145    
146       char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),
147                                                         mpi_size, mpi_rank);
148    
149       /* Figure out how much storage is required for tags */
150       tag_map = mesh->TagMap;
151       num_Tags = 0;
152       if (tag_map) {
153          while (tag_map) {
154             num_Tags++;
155             tag_map=tag_map->next;
156          }
157       }
158    
159       // NetCDF error handler
160       NcError err(NcError::verbose_nonfatal);
161       // Create the file.
162       NcFile dataFile(newFileName, NcFile::Replace);
163       // check if writing was successful
164       if (!dataFile.is_valid())
165          throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);
166    
167       // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)
168       if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
169          throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);
170       if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
171          throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);
172       if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
173          throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);
174       if (num_Elements>0)
175          if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
176             throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);
177       if (num_FaceElements>0)
178          if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
179             throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);
180       if (num_ContactElements>0)
181          if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
182             throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);
183       if (num_Points>0)
184          if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
185             throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);
186       if (num_Elements>0)
187          if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
188             throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);
189       if (num_FaceElements>0)
190          if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
191             throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);
192       if (num_ContactElements>0)
193          if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
194             throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);
195       if (num_Tags>0)
196          if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
197             throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);
198    
199       // Attributes: MPI size, MPI rank, Name, order, reduced_order
200       if (!dataFile.add_att("mpi_size", mpi_size) )
201          throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);
202       if (!dataFile.add_att("mpi_rank", mpi_rank) )
203          throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);
204       if (!dataFile.add_att("Name",mesh->Name) )
205          throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);
206       if (!dataFile.add_att("numDim",numDim) )
207          throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);
208       if (!dataFile.add_att("order",mesh->order) )
209          throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);
210       if (!dataFile.add_att("reduced_order",mesh->reduced_order) )
211          throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);
212       if (!dataFile.add_att("numNodes",numNodes) )
213          throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);
214       if (!dataFile.add_att("num_Elements",num_Elements) )
215          throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);
216       if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
217          throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);
218       if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
219          throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);
220       if (!dataFile.add_att("num_Points",num_Points) )
221          throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);
222       if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
223          throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);
224       if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
225          throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);
226       if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
227          throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);
228       if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )
229          throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);
230       if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )
231          throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);
232       if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )
233          throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);
234       if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )
235          throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);
236       if (!dataFile.add_att("num_Tags", num_Tags) )
237          throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);
238    
239       // // // // // Nodes // // // // //
240    
241       // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)
242       if (numNodes>0) {
243    
244          // Nodes Id
245          if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
246             throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);
247          int_ptr = &mesh->Nodes->Id[0];
248          if (! (ids->put(int_ptr, numNodes)) )
249             throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);
250    
251          // Nodes Tag
252          if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
253             throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);
254          int_ptr = &mesh->Nodes->Tag[0];
255          if (! (ids->put(int_ptr, numNodes)) )
256             throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);
257    
258          // Nodes gDOF
259          if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
260             throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);
261          int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
262          if (! (ids->put(int_ptr, numNodes)) )
263             throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);
264    
265          // Nodes global node index
266          if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
267             throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);
268          int_ptr = &mesh->Nodes->globalNodesIndex[0];
269          if (! (ids->put(int_ptr, numNodes)) )
270             throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);
271    
272          // Nodes grDof
273          if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
274             throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);
275          int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
276          if (! (ids->put(int_ptr, numNodes)) )
277             throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);
278    
279          // Nodes grNI
280          if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
281             throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);
282          int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
283          if (! (ids->put(int_ptr, numNodes)) )
284             throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);
285    
286          // Nodes Coordinates
287          if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
288             throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);
289          if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
290             throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);
291    
292          // Nodes degreesOfFreedomDistribution
293          if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
294             throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);
295          int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
296          if (! (ids->put(int_ptr, mpi_size+1)) )
297             throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);
298    
299          // Nodes nodeDistribution
300          if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
301             throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);
302          int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
303          if (! (ids->put(int_ptr, mpi_size+1)) )
304             throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);
305    
306       }
307    
308       // // // // // Elements // // // // //
309    
310       if (num_Elements>0) {
311    
312          // Elements_Id
313          if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
314             throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);
315          int_ptr = &mesh->Elements->Id[0];
316          if (! (ids->put(int_ptr, num_Elements)) )
317             throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);
318    
319          // Elements_Tag
320          if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
321             throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);
322          int_ptr = &mesh->Elements->Tag[0];
323          if (! (ids->put(int_ptr, num_Elements)) )
324             throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);
325    
326          // Elements_Owner
327          if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
328             throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);
329          int_ptr = &mesh->Elements->Owner[0];
330          if (! (ids->put(int_ptr, num_Elements)) )
331             throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);
332    
333          // Elements_Color
334          if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
335             throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);
336          int_ptr = &mesh->Elements->Color[0];
337          if (! (ids->put(int_ptr, num_Elements)) )
338             throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);
339    
340          // Elements_Nodes
341          if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
342             throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);
343          if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
344             throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);
345    
346       }
347    
348       // // // // // Face_Elements // // // // //
349    
350       if (num_FaceElements>0) {
351    
352          // FaceElements_Id
353          if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
354             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);
355          int_ptr = &mesh->FaceElements->Id[0];
356          if (! (ids->put(int_ptr, num_FaceElements)) )
357             throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);
358    
359          // FaceElements_Tag
360          if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
361             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);
362          int_ptr = &mesh->FaceElements->Tag[0];
363          if (! (ids->put(int_ptr, num_FaceElements)) )
364             throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);
365    
366          // FaceElements_Owner
367          if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
368             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);
369          int_ptr = &mesh->FaceElements->Owner[0];
370          if (! (ids->put(int_ptr, num_FaceElements)) )
371             throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);
372    
373          // FaceElements_Color
374          if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
375             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);
376          int_ptr = &mesh->FaceElements->Color[0];
377          if (! (ids->put(int_ptr, num_FaceElements)) )
378             throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);
379    
380          // FaceElements_Nodes
381          if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
382             throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);
383          if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
384             throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);
385    
386       }
387    
388       // // // // // Contact_Elements // // // // //
389    
390       if (num_ContactElements>0) {
391    
392          // ContactElements_Id
393          if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
394             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);
395          int_ptr = &mesh->ContactElements->Id[0];
396          if (! (ids->put(int_ptr, num_ContactElements)) )
397             throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);
398    
399          // ContactElements_Tag
400          if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
401             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);
402          int_ptr = &mesh->ContactElements->Tag[0];
403          if (! (ids->put(int_ptr, num_ContactElements)) )
404             throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);
405    
406          // ContactElements_Owner
407          if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
408             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);
409          int_ptr = &mesh->ContactElements->Owner[0];
410          if (! (ids->put(int_ptr, num_ContactElements)) )
411             throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);
412    
413          // ContactElements_Color
414          if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
415             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);
416          int_ptr = &mesh->ContactElements->Color[0];
417          if (! (ids->put(int_ptr, num_ContactElements)) )
418             throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);
419    
420          // ContactElements_Nodes
421          if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
422             throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);
423          if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
424             throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);
425    
426       }
427    
428       // // // // // Points // // // // //
429    
430       if (num_Points>0) {
431    
432          fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
433    
434          // Points_Id
435          if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
436             throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);
437          int_ptr = &mesh->Points->Id[0];
438          if (! (ids->put(int_ptr, num_Points)) )
439             throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);
440    
441          // Points_Tag
442          if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
443             throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);
444          int_ptr = &mesh->Points->Tag[0];
445          if (! (ids->put(int_ptr, num_Points)) )
446             throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);
447    
448          // Points_Owner
449          if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
450             throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);
451          int_ptr = &mesh->Points->Owner[0];
452          if (! (ids->put(int_ptr, num_Points)) )
453             throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);
454    
455          // Points_Color
456          if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
457             throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);
458          int_ptr = &mesh->Points->Color[0];
459          if (! (ids->put(int_ptr, num_Points)) )
460             throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);
461    
462          // Points_Nodes
463          // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
464          if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
465             throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);
466          if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
467             throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);
468    
469       }
470    
471       // // // // // TagMap // // // // //
472    
473       if (num_Tags>0) {
474    
475          // Temp storage to gather node IDs
476          int *Tags_keys = TMPMEMALLOC(num_Tags, int);
477          char name_temp[4096];
478    
479          /* Copy tag data into temp arrays */
480          tag_map = mesh->TagMap;
481          if (tag_map) {
482             int i = 0;
483             while (tag_map) {
484                Tags_keys[i++] = tag_map->tag_key;
485                tag_map=tag_map->next;
486             }
487          }
488    
489          // Tags_keys
490          if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
491             throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);
492          int_ptr = &Tags_keys[0];
493          if (! (ids->put(int_ptr, num_Tags)) )
494             throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);
495    
496          // Tags_names_*
497          // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string
498          // because the NetCDF manual doesn't tell how to do an array of strings
499          tag_map = mesh->TagMap;
500          if (tag_map) {
501             int i = 0;
502             while (tag_map) {
503                sprintf(name_temp, "Tags_name_%d", i);
504                if (!dataFile.add_att(name_temp, tag_map->name) )
505                   throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);
506                tag_map=tag_map->next;
507                i++;
508             }
509          }
510    
511          TMPMEMFREE(Tags_keys);
512    
513       }
514    
515    /* Send token to next MPI process so he can take his turn */
516    #ifdef PASO_MPI
517       if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
518    #endif
519    
520       // NetCDF file is closed by destructor of NcFile object
521    
522    #else
523       Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
524    #endif  /* USE_NETCDF */
525       checkFinleyError();
526  }  }
527    
528  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
529  {  {
530    return "FinleyMesh";     return "FinleyMesh";
531  }  }
532    
533  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
534  {  {
535    FunctionSpaceNamesMapType::iterator loc;     FunctionSpaceNamesMapType::iterator loc;
536    loc=m_functionSpaceTypeNames.find(functionSpaceType);     loc=m_functionSpaceTypeNames.find(functionSpaceType);
537    if (loc==m_functionSpaceTypeNames.end()) {     if (loc==m_functionSpaceTypeNames.end()) {
538      return "Invalid function space type code.";        return "Invalid function space type code.";
539    } else {     } else {
540      return loc->second;        return loc->second;
541    }     }
542  }  }
543    
544  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const  bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
545  {  {
546    FunctionSpaceNamesMapType::iterator loc;     FunctionSpaceNamesMapType::iterator loc;
547    loc=m_functionSpaceTypeNames.find(functionSpaceType);     loc=m_functionSpaceTypeNames.find(functionSpaceType);
548    return (loc!=m_functionSpaceTypeNames.end());     return (loc!=m_functionSpaceTypeNames.end());
549  }  }
550    
551  void MeshAdapter::setFunctionSpaceTypeNames()  void MeshAdapter::setFunctionSpaceTypeNames()
552  {  {
553    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
554      (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
555    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
556      (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
557    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
558      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));     (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
559    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
560      (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
561    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
562      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
563    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
564      (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
565    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
566      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
567    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
568      (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
569    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
570      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
571    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
572      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
573    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
574      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
575    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
576      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
577    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
578      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
579  }  }
580    
581  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
582  {  {
583    return Nodes;     return Nodes;
584  }  }
585  int MeshAdapter::getReducedContinuousFunctionCode() const  int MeshAdapter::getReducedContinuousFunctionCode() const
586  {  {
587    return ReducedNodes;     return ReducedNodes;
588  }  }
589    
590  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
591  {  {
592    return Elements;     return Elements;
593  }  }
594  int MeshAdapter::getReducedFunctionCode() const  int MeshAdapter::getReducedFunctionCode() const
595  {  {
596    return ReducedElements;     return ReducedElements;
597  }  }
598    
599  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
600  {  {
601    return FaceElements;     return FaceElements;
602  }  }
603  int MeshAdapter::getReducedFunctionOnBoundaryCode() const  int MeshAdapter::getReducedFunctionOnBoundaryCode() const
604  {  {
605    return ReducedFaceElements;     return ReducedFaceElements;
606  }  }
607    
608  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
609  {  {
610    return ContactElementsZero;     return ContactElementsZero;
611  }  }
612  int MeshAdapter::getReducedFunctionOnContactZeroCode() const  int MeshAdapter::getReducedFunctionOnContactZeroCode() const
613  {  {
614    return ReducedContactElementsZero;     return ReducedContactElementsZero;
615  }  }
616    
617  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
618  {  {
619    return ContactElementsOne;     return ContactElementsOne;
620  }  }
621  int MeshAdapter::getReducedFunctionOnContactOneCode() const  int MeshAdapter::getReducedFunctionOnContactOneCode() const
622  {  {
623    return ReducedContactElementsOne;     return ReducedContactElementsOne;
624  }  }
625    
626  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
627  {  {
628    return DegreesOfFreedom;     return DegreesOfFreedom;
629  }  }
630    
631  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
632  {  {
633    return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
634  }  }
635    
636  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionCode() const
637  {  {
638    return Points;     return Points;
639  }  }
640    
641  //  //
# Line 203  int MeshAdapter::getDiracDeltaFunctionCo Line 643  int MeshAdapter::getDiracDeltaFunctionCo
643  //  //
644  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
645  {  {
646    int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
647    checkFinleyError();     int numDim=Finley_Mesh_getDim(mesh);
648    return numDim;     checkFinleyError();
649       return numDim;
650    }
651    
652    //
653    // Return the number of data points summed across all MPI processes
654    //
655    int MeshAdapter::getNumDataPointsGlobal() const
656    {
657       return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
658  }  }
659    
660  //  //
# Line 218  pair<int,int> MeshAdapter::getDataShape( Line 667  pair<int,int> MeshAdapter::getDataShape(
667     int numSamples=0;     int numSamples=0;
668     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
669     switch (functionSpaceCode) {     switch (functionSpaceCode) {
670        case(Nodes):     case(Nodes):
671             numDataPointsPerSample=1;     numDataPointsPerSample=1;
672             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
673             break;     break;
674        case(ReducedNodes):     case(ReducedNodes):
675             /* TODO: add ReducedNodes */     numDataPointsPerSample=1;
676             numDataPointsPerSample=1;     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
677             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;     break;
678             throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");     case(Elements):
679             break;     if (mesh->Elements!=NULL) {
680        case(Elements):        numSamples=mesh->Elements->numElements;
681             if (mesh->Elements!=NULL) {        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
682               numSamples=mesh->Elements->numElements;     }
683               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;     break;
684             }     case(ReducedElements):
685             break;     if (mesh->Elements!=NULL) {
686        case(ReducedElements):        numSamples=mesh->Elements->numElements;
687             if (mesh->Elements!=NULL) {        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
688               numSamples=mesh->Elements->numElements;     }
689               numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;     break;
690             }     case(FaceElements):
691             break;     if (mesh->FaceElements!=NULL) {
692        case(FaceElements):        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
693             if (mesh->FaceElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
694                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;     }
695                  numSamples=mesh->FaceElements->numElements;     break;
696             }     case(ReducedFaceElements):
697             break;     if (mesh->FaceElements!=NULL) {
698        case(ReducedFaceElements):        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
699             if (mesh->FaceElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
700                  numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;     }
701                  numSamples=mesh->FaceElements->numElements;     break;
702             }     case(Points):
703             break;     if (mesh->Points!=NULL) {
704        case(Points):        numDataPointsPerSample=1;
705             if (mesh->Points!=NULL) {        numSamples=mesh->Points->numElements;
706               numDataPointsPerSample=1;     }
707               numSamples=mesh->Points->numElements;     break;
708             }     case(ContactElementsZero):
709             break;     if (mesh->ContactElements!=NULL) {
710        case(ContactElementsZero):        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
711             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
712               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
713               numSamples=mesh->ContactElements->numElements;     break;
714             }     case(ReducedContactElementsZero):
715             break;     if (mesh->ContactElements!=NULL) {
716        case(ReducedContactElementsZero):        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
717             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
718               numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;     }
719               numSamples=mesh->ContactElements->numElements;     break;
720             }     case(ContactElementsOne):
721             break;     if (mesh->ContactElements!=NULL) {
722        case(ContactElementsOne):        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
723             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
724               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
725               numSamples=mesh->ContactElements->numElements;     break;
726             }     case(ReducedContactElementsOne):
727             break;     if (mesh->ContactElements!=NULL) {
728        case(ReducedContactElementsOne):        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
729             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
730               numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;     }
731               numSamples=mesh->ContactElements->numElements;     break;
732             }     case(DegreesOfFreedom):
733             break;     if (mesh->Nodes!=NULL) {
734        case(DegreesOfFreedom):        numDataPointsPerSample=1;
735             if (mesh->Nodes!=NULL) {        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
736               numDataPointsPerSample=1;     }
737  #ifndef PASO_MPI     break;
738               numSamples=mesh->Nodes->numDegreesOfFreedom;     case(ReducedDegreesOfFreedom):
739  #else     if (mesh->Nodes!=NULL) {
740               numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;        numDataPointsPerSample=1;
741  #endif        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
742             }     }
743             break;     break;
744        case(ReducedDegreesOfFreedom):     default:
745             if (mesh->Nodes!=NULL) {        stringstream temp;
746               numDataPointsPerSample=1;        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
747  #ifndef PASO_MPI        throw FinleyAdapterException(temp.str());
748               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;        break;
749  #else     }
750               numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;     return pair<int,int>(numDataPointsPerSample,numSamples);
 #endif  
            }  
            break;  
       default:  
         stringstream temp;  
         temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();  
         throw FinleyAdapterException(temp.str());  
         break;  
       }  
       return pair<int,int>(numDataPointsPerSample,numSamples);  
751  }  }
752    
753  //  //
754  // 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:
755  //  //
756  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
757                       SystemMatrixAdapter& mat, escript::Data& rhs,                                   SystemMatrixAdapter& mat, escript::Data& rhs,
758                       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,
759                       const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
760                       const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
761  {  {
762     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
763     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
# Line 344  void MeshAdapter::addPDEToSystem( Line 783  void MeshAdapter::addPDEToSystem(
783     checkFinleyError();     checkFinleyError();
784  }  }
785    
786    void  MeshAdapter::addPDEToLumpedSystem(
787                                            escript::Data& mat,
788                                            const escript::Data& D,
789                                            const escript::Data& d) const
790    {
791       escriptDataC _mat=mat.getDataC();
792       escriptDataC _D=D.getDataC();
793       escriptDataC _d=d.getDataC();
794    
795       Finley_Mesh* mesh=m_finleyMesh.get();
796    
797       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
798       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
799    
800       checkFinleyError();
801    }
802    
803    
804  //  //
805  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
806  //  //
# Line 366  void MeshAdapter::addPDEToRHS( escript:: Line 823  void MeshAdapter::addPDEToRHS( escript::
823     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 );
824     checkFinleyError();     checkFinleyError();
825  }  }
826    //
827    // adds PDE of second order into a transport problem
828    //
829    void MeshAdapter::addPDEToTransportProblem(
830                                               TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
831                                               const escript::Data& A, const escript::Data& B, const escript::Data& C,
832                                               const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
833                                               const escript::Data& d, const escript::Data& y,
834                                               const escript::Data& d_contact,const escript::Data& y_contact) const
835    {
836       DataTypes::ShapeType shape;
837       source.expand();
838       escriptDataC _source=source.getDataC();
839       escriptDataC _M=M.getDataC();
840       escriptDataC _A=A.getDataC();
841       escriptDataC _B=B.getDataC();
842       escriptDataC _C=C.getDataC();
843       escriptDataC _D=D.getDataC();
844       escriptDataC _X=X.getDataC();
845       escriptDataC _Y=Y.getDataC();
846       escriptDataC _d=d.getDataC();
847       escriptDataC _y=y.getDataC();
848       escriptDataC _d_contact=d_contact.getDataC();
849       escriptDataC _y_contact=y_contact.getDataC();
850    
851       Finley_Mesh* mesh=m_finleyMesh.get();
852       Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();
853    
854       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
855       checkFinleyError();
856    
857       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
858       checkFinleyError();
859    
860       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
861       checkFinleyError();
862    
863       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
864       checkFinleyError();
865    }
866    
867  //  //
868  // interpolates data between different function spaces:  // interpolates data between different function spaces:
869  //  //
870  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
871  {  {
872    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
873    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
874    if (inDomain!=*this)       if (inDomain!=*this)  
875      throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
876    if (targetDomain!=*this)     if (targetDomain!=*this)
877      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
878    
879    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
880    escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
881    escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
882    switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
883       case(Nodes):     case(Nodes):
884          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
885             case(Nodes):        case(Nodes):
886             case(ReducedNodes):        case(ReducedNodes):
887             case(DegreesOfFreedom):        case(DegreesOfFreedom):
888             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
889                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
890                 break;        break;
891             case(Elements):        case(Elements):
892             case(ReducedElements):        case(ReducedElements):
893                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
894                 break;        break;
895             case(FaceElements):        case(FaceElements):
896             case(ReducedFaceElements):        case(ReducedFaceElements):
897                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
898                 break;        break;
899             case(Points):        case(Points):
900                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
901                 break;        break;
902             case(ContactElementsZero):        case(ContactElementsZero):
903             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
904                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
905                 break;        break;
906             case(ContactElementsOne):        case(ContactElementsOne):
907             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
908                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
909                 break;        break;
910             default:        default:
911                 stringstream temp;           stringstream temp;
912                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
913                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
914                 break;           break;
915          }        }
916          break;        break;
917       case(ReducedNodes):     case(ReducedNodes):
918          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
919             case(Nodes):        case(Nodes):
920             case(ReducedNodes):        case(ReducedNodes):
921             case(DegreesOfFreedom):        case(DegreesOfFreedom):
922             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
923                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
924                 break;        break;
925             case(Elements):        case(Elements):
926             case(ReducedElements):        case(ReducedElements):
927                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
928                 break;        break;
929             case(FaceElements):        case(FaceElements):
930             case(ReducedFaceElements):        case(ReducedFaceElements):
931                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
932                 break;        break;
933             case(Points):        case(Points):
934                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
935                 break;        break;
936             case(ContactElementsZero):        case(ContactElementsZero):
937             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
938                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
939                 break;        break;
940             case(ContactElementsOne):        case(ContactElementsOne):
941             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
942                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
943                 break;        break;
944             default:        default:
945                 stringstream temp;           stringstream temp;
946                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
947                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
948                 break;           break;
949          }        }
950          break;        break;
951       case(Elements):     case(Elements):
952          if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
953             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
954          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
955             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
956          } else {        } else {
957             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
958          }        }
959          break;        break;
960       case(ReducedElements):     case(ReducedElements):
961          if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
962             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
963          } else {        } else {
964             throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
965          }        }
966          break;        break;
967       case(FaceElements):     case(FaceElements):
968          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
969             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
970          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
971             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
972          } else {        } else {
973             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
974         }        }
975         break;        break;
976       case(ReducedFaceElements):     case(ReducedFaceElements):
977          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
978             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
979          } else {        } else {
980             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
981         }        }
982         break;        break;
983       case(Points):     case(Points):
984          if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
985             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
986          } else {        } else {
987             throw FinleyAdapterException("Error - No interpolation with data on points possible.");           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
988          }        }
989          break;        break;
990       case(ContactElementsZero):     case(ContactElementsZero):
991       case(ContactElementsOne):     case(ContactElementsOne):
992          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
993             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
994          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
995             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
996          } else {        } else {
997             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
998             break;        }
999          }        break;
1000          break;     case(ReducedContactElementsZero):
1001       case(ReducedContactElementsZero):     case(ReducedContactElementsOne):
1002       case(ReducedContactElementsOne):        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1003          if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1004             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);        } else {
1005          } else {           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1006             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");        }
1007             break;        break;
1008          }     case(DegreesOfFreedom):      
1009          break;        switch(target.getFunctionSpace().getTypeCode()) {
1010       case(DegreesOfFreedom):              case(ReducedDegreesOfFreedom):
1011          switch(target.getFunctionSpace().getTypeCode()) {        case(DegreesOfFreedom):
1012             case(ReducedDegreesOfFreedom):        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1013             case(DegreesOfFreedom):        break;
1014             case(Nodes):    
1015             case(ReducedNodes):        case(Nodes):
1016                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(ReducedNodes):
1017                break;        if (getMPISize()>1) {
1018  #ifndef PASO_MPI           escript::Data temp=escript::Data(in);
1019             case(Elements):           temp.expand();
1020             case(ReducedElements):           escriptDataC _in2 = temp.getDataC();
1021                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1022                break;        } else {
1023             case(FaceElements):           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1024             case(ReducedFaceElements):        }
1025                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
1026                break;        case(Elements):
1027             case(Points):        case(ReducedElements):
1028                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        if (getMPISize()>1) {
1029                break;           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1030             case(ContactElementsZero):           escriptDataC _in2 = temp.getDataC();
1031             case(ContactElementsOne):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1032             case(ReducedContactElementsZero):        } else {
1033             case(ReducedContactElementsOne):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1034                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        }
1035               break;        break;
1036  #else        case(FaceElements):
1037             /* need to copy Degrees of freedom data to nodal data so that the external values are available */        case(ReducedFaceElements):
1038             case(Elements):        if (getMPISize()>1) {
1039             case(ReducedElements):           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1040                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );           escriptDataC _in2 = temp.getDataC();
1041                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1042                break;    
1043             case(FaceElements):        } else {
1044             case(ReducedFaceElements):           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1045                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );        }
1046                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);        break;
1047                break;        case(Points):
1048             case(Points):        if (getMPISize()>1) {
1049                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1050                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);           escriptDataC _in2 = temp.getDataC();
1051                break;        } else {
1052             case(ContactElementsZero):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1053             case(ContactElementsOne):        }
1054             case(ReducedContactElementsZero):        break;
1055             case(ReducedContactElementsOne):        case(ContactElementsZero):
1056                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );        case(ContactElementsOne):
1057                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);        case(ReducedContactElementsZero):
1058               break;        case(ReducedContactElementsOne):
1059  #endif        if (getMPISize()>1) {
1060             default:           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1061               stringstream temp;           escriptDataC _in2 = temp.getDataC();
1062               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1063               throw FinleyAdapterException(temp.str());        } else {
1064               break;           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1065          }        }
1066          break;        break;
1067       case(ReducedDegreesOfFreedom):        default:
1068         switch(target.getFunctionSpace().getTypeCode()) {           stringstream temp;
1069            case(Nodes):           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1070               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           throw FinleyAdapterException(temp.str());
1071               break;           break;
1072            case(ReducedNodes):        }
1073               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        break;
1074               break;     case(ReducedDegreesOfFreedom):
1075            case(DegreesOfFreedom):        switch(target.getFunctionSpace().getTypeCode()) {
1076               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        case(Nodes):
1077               break;        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1078            case(ReducedDegreesOfFreedom):        break;
1079               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(ReducedNodes):
1080               break;        if (getMPISize()>1) {
1081            case(Elements):           escript::Data temp=escript::Data(in);
1082            case(ReducedElements):           temp.expand();
1083               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           escriptDataC _in2 = temp.getDataC();
1084               break;           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1085            case(FaceElements):        } else {
1086            case(ReducedFaceElements):           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1087               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        }
1088               break;        break;
1089            case(Points):        case(DegreesOfFreedom):
1090               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1091               break;        break;
1092            case(ContactElementsZero):        case(ReducedDegreesOfFreedom):
1093            case(ContactElementsOne):        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1094            case(ReducedContactElementsZero):        break;
1095            case(ReducedContactElementsOne):        case(Elements):
1096               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        case(ReducedElements):
1097               break;        if (getMPISize()>1) {
1098            default:           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1099               stringstream temp;           escriptDataC _in2 = temp.getDataC();
1100               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1101               throw FinleyAdapterException(temp.str());        } else {
1102               break;           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1103         }        }
1104         break;        break;
1105       default:        case(FaceElements):
1106          stringstream temp;        case(ReducedFaceElements):
1107          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        if (getMPISize()>1) {
1108          throw FinleyAdapterException(temp.str());           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1109          break;           escriptDataC _in2 = temp.getDataC();
1110    }           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1111    checkFinleyError();        } else {
1112             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1113          }
1114          break;
1115          case(Points):
1116          if (getMPISize()>1) {
1117             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1118             escriptDataC _in2 = temp.getDataC();
1119             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1120          } else {
1121             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1122          }
1123          break;
1124          case(ContactElementsZero):
1125          case(ContactElementsOne):
1126          case(ReducedContactElementsZero):
1127          case(ReducedContactElementsOne):
1128          if (getMPISize()>1) {
1129             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1130             escriptDataC _in2 = temp.getDataC();
1131             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1132          } else {
1133             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1134          }
1135          break;
1136          default:
1137             stringstream temp;
1138             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1139             throw FinleyAdapterException(temp.str());
1140             break;
1141          }
1142          break;
1143       default:
1144          stringstream temp;
1145          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1146          throw FinleyAdapterException(temp.str());
1147          break;
1148       }
1149       checkFinleyError();
1150  }  }
1151    
1152  //  //
# Line 619  void MeshAdapter::interpolateOnDomain(es Line 1154  void MeshAdapter::interpolateOnDomain(es
1154  //  //
1155  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1156  {  {
1157    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1158    if (argDomain!=*this)     if (argDomain!=*this)
1159       throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1160    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1161    // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1162    if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1163       escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1164       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1165    } else {     } else {
1166       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1167       escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1168       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1169       // this is then interpolated onto arg:        // this is then interpolated onto arg:
1170       interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1171    }     }
1172    checkFinleyError();     checkFinleyError();
1173  }  }
1174    
1175  //  //
# Line 642  void MeshAdapter::setToX(escript::Data& Line 1177  void MeshAdapter::setToX(escript::Data&
1177  //  //
1178  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1179  {  {
1180    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1181    if (normalDomain!=*this)     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1182       throw FinleyAdapterException("Error - Illegal domain of normal locations");     if (normalDomain!=*this)
1183    Finley_Mesh* mesh=m_finleyMesh.get();        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1184    escriptDataC _normal=normal.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1185    switch(normal.getFunctionSpace().getTypeCode()) {     escriptDataC _normal=normal.getDataC();
1186      case(Nodes):     switch(normal.getFunctionSpace().getTypeCode()) {
1187        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     case(Nodes):
1188        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1189      case(ReducedNodes):     break;
1190        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     case(ReducedNodes):
1191        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1192      case(Elements):     break;
1193        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     case(Elements):
1194        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1195      case(ReducedElements):     break;
1196        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     case(ReducedElements):
1197        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1198      case (FaceElements):     break;
1199        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (FaceElements):
1200        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1201      case (ReducedFaceElements):     break;
1202        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (ReducedFaceElements):
1203        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1204      case(Points):     break;
1205        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     case(Points):
1206        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1207      case (ContactElementsOne):     break;
1208      case (ContactElementsZero):     case (ContactElementsOne):
1209        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ContactElementsZero):
1210        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1211      case (ReducedContactElementsOne):     break;
1212      case (ReducedContactElementsZero):     case (ReducedContactElementsOne):
1213        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ReducedContactElementsZero):
1214        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1215      case(DegreesOfFreedom):     break;
1216        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case(DegreesOfFreedom):
1217        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1218      case(ReducedDegreesOfFreedom):     break;
1219        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(ReducedDegreesOfFreedom):
1220        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1221      default:     break;
1222       default:
1223        stringstream temp;        stringstream temp;
1224        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1225        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1226        break;        break;
1227    }     }
1228    checkFinleyError();     checkFinleyError();
1229  }  }
1230    
1231  //  //
# Line 697  void MeshAdapter::setToNormal(escript::D Line 1233  void MeshAdapter::setToNormal(escript::D
1233  //  //
1234  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1235  {  {
1236    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1237    if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1238       throw FinleyAdapterException("Error - Illegal domain of interpolation target");     if (targetDomain!=this)
1239          throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1240    
1241    throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1242  }  }
1243    
1244  //  //
# Line 709  void MeshAdapter::interpolateACross(escr Line 1246  void MeshAdapter::interpolateACross(escr
1246  //  //
1247  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1248  {  {
1249    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1250    if (argDomain!=*this)     if (argDomain!=*this)
1251       throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1252    
1253    Finley_Mesh* mesh=m_finleyMesh.get();     double blocktimer_start = blocktimer_time();
1254    escriptDataC _arg=arg.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1255    switch(arg.getFunctionSpace().getTypeCode()) {     escriptDataC _temp;
1256       case(Nodes):     escript::Data temp;
1257          /* TODO */     escriptDataC _arg=arg.getDataC();
1258          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");     switch(arg.getFunctionSpace().getTypeCode()) {
1259          break;     case(Nodes):
1260       case(ReducedNodes):     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1261          /* TODO */     _temp=temp.getDataC();
1262          throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1263          break;     break;
1264       case(Elements):     case(ReducedNodes):
1265          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1266          break;     _temp=temp.getDataC();
1267       case(ReducedElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1268          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     break;
1269          break;     case(Elements):
1270       case(FaceElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1271          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     break;
1272          break;     case(ReducedElements):
1273       case(ReducedFaceElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1274          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     break;
1275          break;     case(FaceElements):
1276       case(Points):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1277          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     break;
1278          break;     case(ReducedFaceElements):
1279       case(ContactElementsZero):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1280          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1281          break;     case(Points):
1282       case(ReducedContactElementsZero):     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1283          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1284          break;     case(ContactElementsZero):
1285       case(ContactElementsOne):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1286          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1287          break;     case(ReducedContactElementsZero):
1288       case(ReducedContactElementsOne):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1289          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1290          break;     case(ContactElementsOne):
1291       case(DegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1292          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");     break;
1293          break;     case(ReducedContactElementsOne):
1294       case(ReducedDegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1295          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");     break;
1296          break;     case(DegreesOfFreedom):
1297       default:     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1298          stringstream temp;     _temp=temp.getDataC();
1299          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1300          throw FinleyAdapterException(temp.str());     break;
1301          break;     case(ReducedDegreesOfFreedom):
1302    }     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1303    checkFinleyError();     _temp=temp.getDataC();
1304       Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1305       break;
1306       default:
1307          stringstream temp;
1308          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1309          throw FinleyAdapterException(temp.str());
1310          break;
1311       }
1312       checkFinleyError();
1313       blocktimer_increment("integrate()", blocktimer_start);
1314  }  }
1315    
1316  //  //
# Line 771  void MeshAdapter::setToIntegrals(std::ve Line 1318  void MeshAdapter::setToIntegrals(std::ve
1318  //  //
1319  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1320  {  {
1321    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1322    if (argDomain!=*this)     if (argDomain!=*this)
1323      throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1324    const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1325    if (gradDomain!=*this)     if (gradDomain!=*this)
1326       throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1327    
1328    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1329    escriptDataC _grad=grad.getDataC();     escriptDataC _grad=grad.getDataC();
1330    escriptDataC nodeDataC;     escriptDataC nodeDataC;
1331  #ifdef PASO_MPI     escript::Data temp;
1332    escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );     if (getMPISize()>1) {
1333    if( arg.getFunctionSpace().getTypeCode() != Nodes )        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1334    {           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );
1335      Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));           nodeDataC = temp.getDataC();
1336      nodeDataC = nodeTemp.getDataC();        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1337    }           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1338    else           nodeDataC = temp.getDataC();
1339      nodeDataC = arg.getDataC();        } else {
1340  #else           nodeDataC = arg.getDataC();
1341    nodeDataC = arg.getDataC();        }
1342  #endif     } else {
1343    switch(grad.getFunctionSpace().getTypeCode()) {        nodeDataC = arg.getDataC();
1344         case(Nodes):     }
1345            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     switch(grad.getFunctionSpace().getTypeCode()) {
1346            break;     case(Nodes):
1347         case(ReducedNodes):     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1348            throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");     break;
1349            break;     case(ReducedNodes):
1350         case(Elements):     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1351            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     break;
1352            break;     case(Elements):
1353         case(ReducedElements):     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1354            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     break;
1355            break;     case(ReducedElements):
1356         case(FaceElements):     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1357            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     break;
1358            break;     case(FaceElements):
1359         case(ReducedFaceElements):     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1360            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     break;
1361            break;     case(ReducedFaceElements):
1362         case(Points):     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1363            throw FinleyAdapterException("Error - Gradient at points is not supported.");     break;
1364            break;     case(Points):
1365         case(ContactElementsZero):     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1366            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1367            break;     case(ContactElementsZero):
1368         case(ReducedContactElementsZero):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1369            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1370            break;     case(ReducedContactElementsZero):
1371         case(ContactElementsOne):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1372            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1373            break;     case(ContactElementsOne):
1374         case(ReducedContactElementsOne):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1375            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1376            break;     case(ReducedContactElementsOne):
1377         case(DegreesOfFreedom):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1378            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     break;
1379            break;     case(DegreesOfFreedom):
1380         case(ReducedDegreesOfFreedom):     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1381            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     break;
1382            break;     case(ReducedDegreesOfFreedom):
1383         default:     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1384            stringstream temp;     break;
1385            temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     default:
1386            throw FinleyAdapterException(temp.str());        stringstream temp;
1387            break;        temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1388    }        throw FinleyAdapterException(temp.str());
1389    checkFinleyError();        break;
1390       }
1391       checkFinleyError();
1392  }  }
1393    
1394  //  //
# Line 847  void MeshAdapter::setToGradient(escript: Line 1396  void MeshAdapter::setToGradient(escript:
1396  //  //
1397  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1398  {  {
1399    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1400    escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1401    switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1402         case(Nodes):     case(Nodes):
1403            throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1404            break;     break;
1405         case(ReducedNodes):     case(ReducedNodes):
1406            throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1407            break;     break;
1408         case(Elements):     case(Elements):
1409            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1410            break;     break;
1411         case(ReducedElements):     case(ReducedElements):
1412            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1413            break;     break;
1414         case(FaceElements):     case(FaceElements):
1415            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1416            break;     break;
1417         case(ReducedFaceElements):     case(ReducedFaceElements):
1418            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1419            break;     break;
1420         case(Points):     case(Points):
1421            throw FinleyAdapterException("Error - Size of point elements is not supported.");     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1422            break;     break;
1423         case(ContactElementsZero):     case(ContactElementsZero):
1424         case(ContactElementsOne):     case(ContactElementsOne):
1425            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1426            break;     break;
1427         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1428         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1429            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1430            break;     break;
1431         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1432            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1433            break;     break;
1434         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1435            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1436            break;     break;
1437         default:     default:
1438            stringstream temp;        stringstream temp;
1439            temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();        temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1440            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1441            break;        break;
1442    }     }
1443    checkFinleyError();     checkFinleyError();
1444  }  }
1445    
1446  // sets the location of nodes:  //
1447    // sets the location of nodes
1448    //
1449  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1450  {  {
1451    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1452    escriptDataC tmp;     escriptDataC tmp;
1453    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1454    if (newDomain!=*this)     if (newDomain!=*this)
1455       throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1456    tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1457    Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1458    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1459       } else {
1460           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );
1461           tmp = new_x_inter.getDataC();
1462           Finley_Mesh_setCoordinates(mesh,&tmp);
1463       }
1464       checkFinleyError();
1465  }  }
1466    
1467  // saves a data array in openDX format:  //
1468    // Helper for the save* methods. Extracts optional data variable names and the
1469    // corresponding pointers from python dictionary. Caller must free arrays.
1470    //
1471    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1472    {
1473       numData = boost::python::extract<int>(arg.attr("__len__")());
1474       /* win32 refactor */
1475       names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1476       data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1477       dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1478    
1479       boost::python::list keys=arg.keys();
1480       for (int i=0; i<numData; ++i) {
1481          std::string n=boost::python::extract<std::string>(keys[i]);
1482          escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1483          if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1484             throw FinleyAdapterException("Error: Data must be defined on same Domain");
1485          data[i] = d.getDataC();
1486          dataPtr[i] = &(data[i]);
1487          names[i] = TMPMEMALLOC(n.length()+1, char);
1488          strcpy(names[i], n.c_str());
1489       }
1490    }
1491    
1492    //
1493    // saves mesh and optionally data arrays in openDX format
1494    //
1495  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1496  {  {
1497      int MAX_namelength=256;     int num_data;
1498      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1499    /* win32 refactor */     escriptDataC *data;
1500    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     escriptDataC **ptr_data;
   for(int i=0;i<num_data;i++)  
   {  
     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
   }  
1501    
1502    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1503    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1504    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     checkFinleyError();
1505    
1506      boost::python::list keys=arg.keys();     /* win32 refactor */
1507      for (int i=0;i<num_data;++i) {     TMPMEMFREE(data);
1508           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);     TMPMEMFREE(ptr_data);
1509           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)     for(int i=0; i<num_data; i++)
1510               throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");     {
1511           data[i]=d.getDataC();        TMPMEMFREE(names[i]);
1512           ptr_data[i]=&(data[i]);     }
1513           std::string n=boost::python::extract<std::string>(keys[i]);     TMPMEMFREE(names);
          c_names[i]=&(names[i][0]);  
          if (n.length()>MAX_namelength-1) {  
             strncpy(c_names[i],n.c_str(),MAX_namelength-1);  
             c_names[i][MAX_namelength-1]='\0';  
          } else {  
             strcpy(c_names[i],n.c_str());  
          }  
     }  
     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
     checkFinleyError();  
       
       /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1514    
1515      return;     return;
1516  }  }
1517    
1518  // saves a data array in openVTK format:  //
1519  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1520    //
1521    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const
1522  {  {
1523      int MAX_namelength=256;     int num_data;
1524      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1525    /* win32 refactor */     escriptDataC *data;
1526    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     escriptDataC **ptr_data;
   for(int i=0;i<num_data;i++)  
   {  
     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
   }  
1527    
1528    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1529    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1530    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     checkFinleyError();
   
     boost::python::list keys=arg.keys();  
     for (int i=0;i<num_data;++i) {  
          escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
          if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  
              throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
          data[i]=d.getDataC();  
          ptr_data[i]=&(data[i]);  
          std::string n=boost::python::extract<std::string>(keys[i]);  
          c_names[i]=&(names[i][0]);  
          if (n.length()>MAX_namelength-1) {  
             strncpy(c_names[i],n.c_str(),MAX_namelength-1);  
             c_names[i][MAX_namelength-1]='\0';  
          } else {  
             strcpy(c_names[i],n.c_str());  
          }  
     }  
 #ifndef PASO_MPI      
     Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
 #else  
     Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
 #endif  
   
 checkFinleyError();  
   /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1531    
1532      return;     /* win32 refactor */
1533       TMPMEMFREE(data);
1534       TMPMEMFREE(ptr_data);
1535       for(int i=0; i<num_data; i++)
1536       {
1537          TMPMEMFREE(names[i]);
1538       }
1539       TMPMEMFREE(names);
1540  }  }
1541                                                                                                                                                                            
1542                                                                                                                                                                            //
1543  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1544    //
1545  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1546                        const int row_blocksize,                                                   const int row_blocksize,
1547                        const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1548                        const int column_blocksize,                                                   const int column_blocksize,
1549                        const escript::FunctionSpace& column_functionspace,                                                   const escript::FunctionSpace& column_functionspace,
1550                        const int type) const                                                   const int type) const
1551  {  {
1552      int reduceRowOrder=0;     int reduceRowOrder=0;
1553      int reduceColOrder=0;     int reduceColOrder=0;
1554      // is the domain right?     // is the domain right?
1555      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1556      if (row_domain!=*this)     if (row_domain!=*this)
1557            throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1558      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1559      if (col_domain!=*this)     if (col_domain!=*this)
1560            throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1561      // is the function space type right     // is the function space type right
1562      if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1563          reduceRowOrder=0;        reduceRowOrder=0;
1564      } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1565          reduceRowOrder=1;        reduceRowOrder=1;
1566      } else {     } else {
1567          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1568      }     }
1569      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1570          reduceColOrder=0;        reduceColOrder=0;
1571      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1572          reduceColOrder=1;        reduceColOrder=1;
1573      } else {     } else {
1574          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");        throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1575      }     }
1576      // generate matrix:     // generate matrix:
1577        
1578      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1579      checkFinleyError();     checkFinleyError();
1580      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);     Paso_SystemMatrix* fsystemMatrix;
1581      checkPasoError();     int trilinos = 0;
1582      Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);     if (trilinos) {
1583      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);  #ifdef TRILINOS
1584          /* Allocation Epetra_VrbMatrix here */
1585    #endif
1586       }
1587       else {
1588          fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1589       }
1590       checkPasoError();
1591       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1592       return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1593    }
1594    
1595    //
1596    // creates a TransportProblemAdapter
1597    //
1598    TransportProblemAdapter MeshAdapter::newTransportProblem(
1599                                                             const double theta,
1600                                                             const int blocksize,
1601                                                             const escript::FunctionSpace& functionspace,
1602                                                             const int type) const
1603    {
1604       int reduceOrder=0;
1605       // is the domain right?
1606       const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1607       if (domain!=*this)
1608          throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1609       // is the function space type right
1610       if (functionspace.getTypeCode()==DegreesOfFreedom) {
1611          reduceOrder=0;
1612       } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1613          reduceOrder=1;
1614       } else {
1615          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1616       }
1617       // generate matrix:
1618    
1619       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1620       checkFinleyError();
1621       Paso_FCTransportProblem* transportProblem;
1622       transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1623       checkPasoError();
1624       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1625       return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
1626  }  }
1627    
1628  //  //
# Line 1059  SystemMatrixAdapter MeshAdapter::newSyst Line 1634  SystemMatrixAdapter MeshAdapter::newSyst
1634  // returns true if data at the atom_type is considered as being cell centered:  // returns true if data at the atom_type is considered as being cell centered:
1635  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1636  {  {
1637    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1638         case(Nodes):     case(Nodes):
1639         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1640         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1641            return false;     return false;
1642            break;     break;
1643         case(Elements):     case(Elements):
1644         case(FaceElements):     case(FaceElements):
1645         case(Points):     case(Points):
1646         case(ContactElementsZero):     case(ContactElementsZero):
1647         case(ContactElementsOne):     case(ContactElementsOne):
1648         case(ReducedElements):     case(ReducedElements):
1649         case(ReducedFaceElements):     case(ReducedFaceElements):
1650         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1651         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1652            return true;     return true;
1653            break;     break;
1654         default:     default:
1655            stringstream temp;        stringstream temp;
1656            temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1657            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1658            break;        break;
1659    }     }
1660    checkFinleyError();     checkFinleyError();
1661    return false;     return false;
1662    }
1663    
1664    bool
1665    MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1666    {
1667       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1668        class 1: DOF <-> Nodes
1669        class 2: ReducedDOF <-> ReducedNodes
1670        class 3: Points
1671        class 4: Elements
1672        class 5: ReducedElements
1673        class 6: FaceElements
1674        class 7: ReducedFaceElements
1675        class 8: ContactElementZero <-> ContactElementOne
1676        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1677    
1678       There is also a set of lines. Interpolation is possible down a line but not between lines.
1679       class 1 and 2 belong to all lines so aren't considered.
1680        line 0: class 3
1681        line 1: class 4,5
1682        line 2: class 6,7
1683        line 3: class 8,9
1684    
1685       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1686       eg hasnodes is true if we have at least one instance of Nodes.
1687       */
1688        if (fs.size()==0)
1689        {
1690        return false;
1691        }
1692        std::vector<int> hasclass(10);
1693        std::vector<int> hasline(4);    
1694        bool hasnodes=false;
1695        bool hasrednodes=false;
1696        bool hascez=false;
1697        bool hasrcez=false;
1698        for (int i=0;i<fs.size();++i)
1699        {
1700        switch(fs[i])
1701        {
1702        case(Nodes):   hasnodes=true;   // no break is deliberate
1703        case(DegreesOfFreedom):
1704            hasclass[1]=1;
1705            break;
1706        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1707        case(ReducedDegreesOfFreedom):
1708            hasclass[2]=1;
1709            break;
1710        case(Points):
1711            hasline[0]=1;
1712            hasclass[3]=1;
1713            break;
1714        case(Elements):
1715            hasclass[4]=1;
1716            hasline[1]=1;
1717            break;
1718        case(ReducedElements):
1719            hasclass[5]=1;
1720            hasline[1]=1;
1721            break;
1722        case(FaceElements):
1723            hasclass[6]=1;
1724            hasline[2]=1;
1725            break;
1726        case(ReducedFaceElements):
1727            hasclass[7]=1;
1728            hasline[2]=1;
1729            break;
1730        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1731        case(ContactElementsOne):
1732            hasclass[8]=1;
1733            hasline[3]=1;
1734            break;
1735        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1736        case(ReducedContactElementsOne):
1737            hasclass[9]=1;
1738            hasline[3]=1;
1739            break;
1740        default:
1741            return false;
1742        }
1743        }
1744        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1745        // fail if we have more than one leaf group
1746    
1747        if (totlines>1)
1748        {
1749        return false;   // there are at least two branches we can't interpolate between
1750        }
1751        else if (totlines==1)
1752        {
1753        if (hasline[0]==1)      // we have points
1754        {
1755            resultcode=Points;
1756        }
1757        else if (hasline[1]==1)
1758        {
1759            if (hasclass[5]==1)
1760            {
1761            resultcode=ReducedElements;
1762            }
1763            else
1764            {
1765            resultcode=Elements;
1766            }
1767        }
1768        else if (hasline[2]==1)
1769        {
1770            if (hasclass[7]==ReducedFaceElements)
1771            {
1772            resultcode=ReducedFaceElements;
1773            }
1774            else
1775            {
1776            resultcode=FaceElements;
1777            }
1778        }
1779        else    // so we must be in line3
1780        {
1781            if (hasclass[9]==1)
1782            {
1783            // need something from class 9
1784            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1785            }
1786            else
1787            {
1788            // something from class 8
1789            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1790            }
1791        }
1792        }
1793        else    // totlines==0
1794        {
1795        if (hasclass[2]==1)
1796        {
1797            // something from class 2
1798            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1799        }
1800        else
1801        {   // something from class 1
1802            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1803        }
1804        }
1805        return true;
1806  }  }
1807    
1808  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1809  {  {
1810    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1811       case(Nodes):     case(Nodes):
1812          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1813             case(Nodes):      case(Nodes):
1814             case(ReducedNodes):      case(ReducedNodes):
1815             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1816             case(DegreesOfFreedom):      case(DegreesOfFreedom):
1817             case(Elements):      case(Elements):
1818             case(ReducedElements):      case(ReducedElements):
1819             case(FaceElements):      case(FaceElements):
1820             case(ReducedFaceElements):      case(ReducedFaceElements):
1821             case(Points):      case(Points):
1822             case(ContactElementsZero):      case(ContactElementsZero):
1823             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1824             case(ContactElementsOne):      case(ContactElementsOne):
1825             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1826                 return true;      return true;
1827             default:      default:
1828                 stringstream temp;            stringstream temp;
1829                 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;
1830                 throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1831          }     }
1832          break;     break;
1833       case(ReducedNodes):     case(ReducedNodes):
1834          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1835             case(ReducedNodes):      case(ReducedNodes):
1836             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1837             case(Elements):      case(Elements):
1838             case(ReducedElements):      case(ReducedElements):
1839             case(FaceElements):      case(FaceElements):
1840             case(ReducedFaceElements):      case(ReducedFaceElements):
1841             case(Points):      case(Points):
1842             case(ContactElementsZero):      case(ContactElementsZero):
1843             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1844             case(ContactElementsOne):      case(ContactElementsOne):
1845             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1846                 return true;      return true;
1847            case(Nodes):      case(Nodes):
1848            case(DegreesOfFreedom):      case(DegreesOfFreedom):
1849               return false;      return false;
1850             default:      default:
1851                 stringstream temp;          stringstream temp;
1852                 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;
1853                 throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1854          }     }
1855          break;     break;
1856       case(Elements):     case(Elements):
1857          if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1858             return true;        return true;
1859          } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1860             return true;        return true;
         } else {  
            return false;  
         }  
      case(ReducedElements):  
         if (functionSpaceType_target==ReducedElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(FaceElements):  
         if (functionSpaceType_target==FaceElements) {  
            return true;  
         } else if (functionSpaceType_target==ReducedFaceElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ReducedFaceElements):  
         if (functionSpaceType_target==ReducedFaceElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(Points):  
         if (functionSpaceType_target==Points) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ContactElementsZero):  
      case(ContactElementsOne):  
         if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
            return true;  
         } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ReducedContactElementsZero):  
      case(ReducedContactElementsOne):  
         if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
            return true;  
1861          } else {          } else {
1862             return false;            return false;
         }  
      case(DegreesOfFreedom):  
         switch(functionSpaceType_target) {  
            case(ReducedDegreesOfFreedom):  
            case(DegreesOfFreedom):  
            case(Nodes):  
            case(ReducedNodes):  
            case(Elements):  
            case(ReducedElements):  
            case(Points):  
            case(FaceElements):  
            case(ReducedFaceElements):  
            case(ContactElementsZero):  
            case(ReducedContactElementsZero):  
            case(ContactElementsOne):  
            case(ReducedContactElementsOne):  
               return true;  
            default:  
              stringstream temp;  
              temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
              throw FinleyAdapterException(temp.str());  
1863          }          }
1864          break;     case(ReducedElements):
1865       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1866         switch(functionSpaceType_target) {        return true;
1867            case(ReducedDegreesOfFreedom):      } else {
1868            case(ReducedNodes):            return false;
1869            case(Elements):      }
1870            case(ReducedElements):     case(FaceElements):
1871            case(FaceElements):      if (functionSpaceType_target==FaceElements) {
1872            case(ReducedFaceElements):              return true;
1873            case(Points):      } else if (functionSpaceType_target==ReducedFaceElements) {
1874            case(ContactElementsZero):              return true;
1875            case(ReducedContactElementsZero):      } else {
1876            case(ContactElementsOne):              return false;
1877            case(ReducedContactElementsOne):      }
1878                return true;     case(ReducedFaceElements):
1879            case(Nodes):      if (functionSpaceType_target==ReducedFaceElements) {
1880            case(DegreesOfFreedom):              return true;
1881               return false;      } else {
1882            default:          return false;
1883               stringstream temp;      }
1884               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;     case(Points):
1885               throw FinleyAdapterException(temp.str());      if (functionSpaceType_target==Points) {
1886         }              return true;
1887         break;      } else {
1888       default:              return false;
1889          stringstream temp;      }
1890          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;     case(ContactElementsZero):
1891          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
1892          break;      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1893    }              return true;
1894    checkFinleyError();      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1895    return false;              return true;
1896        } else {
1897                return false;
1898        }
1899       case(ReducedContactElementsZero):
1900       case(ReducedContactElementsOne):
1901        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1902                return true;
1903        } else {
1904                return false;
1905        }
1906       case(DegreesOfFreedom):
1907        switch(functionSpaceType_target) {
1908        case(ReducedDegreesOfFreedom):
1909        case(DegreesOfFreedom):
1910        case(Nodes):
1911        case(ReducedNodes):
1912        case(Elements):
1913        case(ReducedElements):
1914        case(Points):
1915        case(FaceElements):
1916        case(ReducedFaceElements):
1917        case(ContactElementsZero):
1918        case(ReducedContactElementsZero):
1919        case(ContactElementsOne):
1920        case(ReducedContactElementsOne):
1921        return true;
1922        default:
1923            stringstream temp;
1924            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1925            throw FinleyAdapterException(temp.str());
1926        }
1927        break;
1928       case(ReducedDegreesOfFreedom):
1929       switch(functionSpaceType_target) {
1930        case(ReducedDegreesOfFreedom):
1931        case(ReducedNodes):
1932        case(Elements):
1933        case(ReducedElements):
1934        case(FaceElements):
1935        case(ReducedFaceElements):
1936        case(Points):
1937        case(ContactElementsZero):
1938        case(ReducedContactElementsZero):
1939        case(ContactElementsOne):
1940        case(ReducedContactElementsOne):
1941        return true;
1942        case(Nodes):
1943        case(DegreesOfFreedom):
1944        return false;
1945        default:
1946            stringstream temp;
1947            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1948            throw FinleyAdapterException(temp.str());
1949        }
1950        break;
1951       default:
1952          stringstream temp;
1953          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1954          throw FinleyAdapterException(temp.str());
1955          break;
1956       }
1957       checkFinleyError();
1958       return false;
1959  }  }
1960    
1961  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1962  {  {
1963      return false;     return false;
1964  }  }
1965    
1966  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
1967  {  {
1968    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1969    if (temp!=0) {     if (temp!=0) {
1970      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
1971    } else {     } else {
1972      return false;        return false;
1973    }     }
1974  }  }
1975    
1976  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
1977  {  {
1978    return !(operator==(other));     return !(operator==(other));
1979  }  }
1980    
1981  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1982  {  {
1983     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
1984       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1985       checkPasoError();
1986       return out;
1987    }
1988    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1989    {
1990       Finley_Mesh* mesh=m_finleyMesh.get();
1991       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1992     checkPasoError();     checkPasoError();
1993     return out;     return out;
1994  }  }
1995    
1996  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
1997  {  {
1998    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(asAbstractContinuousDomain()).getX();
1999  }  }
2000    
2001  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2002  {  {
2003    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
2004  }  }
2005    
2006  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2007  {  {
2008    return function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
2009  }  }
2010    
2011  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2012  {  {
2013    int *out=0,i;     int *out = NULL;
2014    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2015    switch (functionSpaceType) {     switch (functionSpaceType) {
2016      case(Nodes):     case(Nodes):
2017        if (mesh->Nodes!=NULL) {     out=mesh->Nodes->Id;
2018          out=mesh->Nodes->Id;     break;
2019          break;     case(ReducedNodes):
2020        }     out=mesh->Nodes->reducedNodesId;
2021      case(ReducedNodes):     break;
2022        throw FinleyAdapterException("Error -  ReducedNodes not supported yet.");     case(Elements):
2023        break;     out=mesh->Elements->Id;
2024      case(Elements):     break;
2025        out=mesh->Elements->Id;     case(ReducedElements):
2026        break;     out=mesh->Elements->Id;
2027      case(ReducedElements):     break;
2028        out=mesh->Elements->Id;     case(FaceElements):
2029        break;     out=mesh->FaceElements->Id;
2030      case(FaceElements):     break;
2031        out=mesh->FaceElements->Id;     case(ReducedFaceElements):
2032        break;     out=mesh->FaceElements->Id;
2033      case(ReducedFaceElements):     break;
2034        out=mesh->FaceElements->Id;     case(Points):
2035        break;     out=mesh->Points->Id;
2036      case(Points):     break;
2037        out=mesh->Points->Id;     case(ContactElementsZero):
2038        break;     out=mesh->ContactElements->Id;
2039      case(ContactElementsZero):     break;
2040        out=mesh->ContactElements->Id;     case(ReducedContactElementsZero):
2041        break;     out=mesh->ContactElements->Id;
2042      case(ReducedContactElementsZero):     break;
2043        out=mesh->ContactElements->Id;     case(ContactElementsOne):
2044        break;     out=mesh->ContactElements->Id;
2045      case(ContactElementsOne):     break;
2046        out=mesh->ContactElements->Id;     case(ReducedContactElementsOne):
2047        break;     out=mesh->ContactElements->Id;
2048      case(ReducedContactElementsOne):     break;
2049        out=mesh->ContactElements->Id;     case(DegreesOfFreedom):
2050        break;     out=mesh->Nodes->degreesOfFreedomId;
2051      case(DegreesOfFreedom):     break;
2052        out=mesh->Nodes->degreeOfFreedomId;     case(ReducedDegreesOfFreedom):
2053        break;     out=mesh->Nodes->reducedDegreesOfFreedomId;
2054      case(ReducedDegreesOfFreedom):     break;
2055        out=mesh->Nodes->reducedDegreeOfFreedomId;     default:
       break;  
     default:  
2056        stringstream temp;        stringstream temp;
2057        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2058        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2059        break;        break;
2060    }     }
2061    return out;     return out;
2062  }  }
2063  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2064  {  {
2065    int out=0;     int out=0;
2066    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2067    switch (functionSpaceType) {     switch (functionSpaceType) {
2068    case(Nodes):     case(Nodes):
2069      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
2070      break;     break;
2071    case(ReducedNodes):     case(ReducedNodes):
2072      throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2073      break;     break;
2074    case(Elements):     case(Elements):
2075      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2076      break;     break;
2077    case(ReducedElements):     case(ReducedElements):
2078      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2079      break;     break;
2080    case(FaceElements):     case(FaceElements):
2081      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2082      break;     break;
2083    case(ReducedFaceElements):     case(ReducedFaceElements):
2084      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2085      break;     break;
2086    case(Points):     case(Points):
2087      out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
2088      break;     break;
2089    case(ContactElementsZero):     case(ContactElementsZero):
2090      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2091      break;     break;
2092    case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2093      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2094      break;     break;
2095    case(ContactElementsOne):     case(ContactElementsOne):
2096      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2097      break;     break;
2098    case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2099      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2100      break;     break;
2101    case(DegreesOfFreedom):     case(DegreesOfFreedom):
2102      throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2103      break;     break;
2104    case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2105      throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2106      break;     break;
2107    default:     default:
2108      stringstream temp;        stringstream temp;
2109      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2110      throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2111      break;        break;
2112    }     }
2113    return out;     return out;
2114  }  }
2115    
2116    
2117  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2118  {  {
2119       Finley_Mesh* mesh=m_finleyMesh.get();
2120       escriptDataC tmp=mask.getDataC();
2121       switch(functionSpaceType) {
2122       case(Nodes):
2123       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2124       break;
2125       case(ReducedNodes):
2126       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2127       break;
2128       case(DegreesOfFreedom):
2129       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2130       break;
2131       case(ReducedDegreesOfFreedom):
2132       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2133       break;
2134       case(Elements):
2135       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2136       break;
2137       case(ReducedElements):
2138       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2139       break;
2140       case(FaceElements):
2141       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2142       break;
2143       case(ReducedFaceElements):
2144       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2145       break;
2146       case(Points):
2147       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2148       break;
2149       case(ContactElementsZero):
2150       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2151       break;
2152       case(ReducedContactElementsZero):
2153       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2154       break;
2155       case(ContactElementsOne):
2156       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2157       break;
2158       case(ReducedContactElementsOne):
2159       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2160       break;
2161       default:
2162          stringstream temp;
2163          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2164          throw FinleyAdapterException(temp.str());
2165       }
2166       checkFinleyError();
2167       return;
2168    }
2169    
2170    void MeshAdapter::setTagMap(const std::string& name,  int tag)
2171    {
2172       Finley_Mesh* mesh=m_finleyMesh.get();
2173       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2174       checkPasoError();
2175       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2176    }
2177    
2178    int MeshAdapter::getTag(const std::string& name) const
2179    {
2180       Finley_Mesh* mesh=m_finleyMesh.get();
2181       int tag=0;
2182       tag=Finley_Mesh_getTag(mesh, name.c_str());
2183       checkPasoError();
2184       // throwStandardException("MeshAdapter::getTag is not implemented.");
2185       return tag;
2186    }
2187    
2188    bool MeshAdapter::isValidTagName(const std::string& name) const
2189    {
2190       Finley_Mesh* mesh=m_finleyMesh.get();
2191       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2192    }
2193    
2194    std::string MeshAdapter::showTagNames() const
2195    {
2196       stringstream temp;
2197       Finley_Mesh* mesh=m_finleyMesh.get();
2198       Finley_TagMap* tag_map=mesh->TagMap;
2199       while (tag_map) {
2200          temp << tag_map->name;
2201          tag_map=tag_map->next;
2202          if (tag_map) temp << ", ";
2203       }
2204       return temp.str();
2205    }
2206    
2207    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2208    {
2209    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2210    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2211    switch(functionSpaceType) {    switch(functionSpaceCode) {
2212         case(Nodes):     case(Nodes):
2213            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2214            break;            break;
2215         case(ReducedNodes):     case(ReducedNodes):
2216            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2217            break;            break;
2218         case(DegreesOfFreedom):     case(DegreesOfFreedom):
2219            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2220            break;            break;
2221         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2222            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2223            break;            break;
2224         case(Elements):     case(Elements):
2225            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2226            break;            numTags=mesh->Elements->numTagsInUse;
2227         case(ReducedElements):            break;
2228            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(FaceElements):
2229       case(ReducedFaceElements):
2230              numTags=mesh->FaceElements->numTagsInUse;
2231              break;
2232       case(Points):
2233              numTags=mesh->Points->numTagsInUse;
2234              break;
2235       case(ContactElementsZero):
2236       case(ReducedContactElementsZero):
2237       case(ContactElementsOne):
2238       case(ReducedContactElementsOne):
2239              numTags=mesh->ContactElements->numTagsInUse;
2240            break;            break;
2241         case(FaceElements):     default:
2242            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);        stringstream temp;
2243            break;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2244         case(ReducedFaceElements):        throw FinleyAdapterException(temp.str());
2245            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);    }
2246            break;    return numTags;
2247         case(Points):  }
2248            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);  
2249    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2250    {
2251      Finley_Mesh* mesh=m_finleyMesh.get();
2252      index_t* tags=NULL;
2253      switch(functionSpaceCode) {
2254       case(Nodes):
2255              tags=mesh->Nodes->tagsInUse;
2256            break;            break;
2257         case(ContactElementsZero):     case(ReducedNodes):
2258            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2259            break;            break;
2260         case(ReducedContactElementsZero):     case(DegreesOfFreedom):
2261            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2262            break;            break;
2263         case(ContactElementsOne):     case(ReducedDegreesOfFreedom):
2264            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2265            break;            break;
2266         case(ReducedContactElementsOne):     case(Elements):
2267            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2268              tags=mesh->Elements->tagsInUse;
2269              break;
2270       case(FaceElements):
2271       case(ReducedFaceElements):
2272              tags=mesh->FaceElements->tagsInUse;
2273              break;
2274       case(Points):
2275              tags=mesh->Points->tagsInUse;
2276              break;
2277       case(ContactElementsZero):
2278       case(ReducedContactElementsZero):
2279       case(ContactElementsOne):
2280       case(ReducedContactElementsOne):
2281              tags=mesh->ContactElements->tagsInUse;
2282            break;            break;
2283         default:     default:
2284            stringstream temp;        stringstream temp;
2285            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2286            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2287    }    }
2288    checkFinleyError();    return tags;
   return;  
2289  }  }
2290    
 void MeshAdapter::setTagMap(const std::string& name,  int tag)  
 {  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   Finley_Mesh_addTagMap(mesh, name.c_str(),tag);  
   checkPasoError();  
   // throwStandardException("MeshAdapter::set TagMap is not implemented.");  
 }  
2291    
2292  int MeshAdapter::getTag(const std::string& name) const  bool MeshAdapter::canTag(int functionSpaceCode) const
2293  {  {
2294    Finley_Mesh* mesh=m_finleyMesh.get();    switch(functionSpaceCode) {
2295    int tag=0;     case(Nodes):
2296    tag=Finley_Mesh_getTag(mesh, name.c_str());     case(Elements):
2297    checkPasoError();     case(ReducedElements):
2298    // throwStandardException("MeshAdapter::getTag is not implemented.");     case(FaceElements):
2299    return tag;     case(ReducedFaceElements):
2300       case(Points):
2301       case(ContactElementsZero):
2302       case(ReducedContactElementsZero):
2303       case(ContactElementsOne):
2304       case(ReducedContactElementsOne):
2305              return true;
2306       case(ReducedNodes):
2307       case(DegreesOfFreedom):
2308       case(ReducedDegreesOfFreedom):
2309          return false;
2310       default:
2311        return false;
2312      }
2313  }  }
2314    
2315  bool MeshAdapter::isValidTagName(const std::string& name) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2316  {  {
2317    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2318    return Finley_Mesh_isValidTagName(mesh,name.c_str());    return Finley_Mesh_getStatus(mesh);
2319  }  }
2320    
2321  std::string MeshAdapter::showTagNames() const  
 {  
   stringstream temp;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   Finley_TagMap* tag_map=mesh->TagMap;  
   while (tag_map) {  
      temp << tag_map->name;  
      tag_map=tag_map->next;  
      if (tag_map) temp << ", ";  
   }  
   return temp.str();  
 }  
2322    
2323  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1116  
changed lines
  Added in v.2635

  ViewVC Help
Powered by ViewVC 1.1.26