/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.102  
changed lines
  Added in v.3993

  ViewVC Help
Powered by ViewVC 1.1.26