/[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

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

Legend:
Removed from v.155  
changed lines
  Added in v.4346

  ViewVC Help
Powered by ViewVC 1.1.26