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

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

  ViewVC Help
Powered by ViewVC 1.1.26