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

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 3736 by jfenwick, Mon Dec 12 01:26:51 2011 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ******************************************************************************  *
4   *                                                                            *  * Copyright (c) 2003-2010 by University of Queensland
5   *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  * Earth Systems Science Computational Center (ESSCC)
6   *                                                                            *  * http://www.uq.edu.au/esscc
7   * This software is the property of ACcESS. No part of this code              *  *
8   * may be copied in any form or by any means without the expressed written    *  * Primary Business: Queensland, Australia
9   * consent of ACcESS.  Copying, use or modification of this software          *  * Licensed under the Open Software License version 3.0
10   * by any unauthorised person is illegal unless that person has a software    *  * http://www.opensource.org/licenses/osl-3.0.php
11   * license agreement with ACcESS.                                             *  *
12   *                                                                            *  *******************************************************/
13   ******************************************************************************  
14  */  #include <pasowrap/PasoException.h>
15    #include <pasowrap/TransportProblemAdapter.h>
16    #include "MeshAdapter.h"
17    #include "escript/Data.h"
18    #include "escript/DataFactory.h"
19    #ifdef USE_NETCDF
20    #include <netcdfcpp.h>
21    #endif
22  extern "C" {  extern "C" {
23  #include "finley/finleyC/Finley.h"  #include "esysUtils/blocktimer.h"
24  #include "finley/finleyC/Assemble.h"  }
25  #include "finley/finleyC/Mesh.h"  
26  #include "finley/finleyC/Finley.h"  #include <boost/python/import.hpp>
 #include "finley/finleyC/System.h"  
 }  
 #include "finley/CPPAdapter/SystemMatrixAdapter.h"  
 #include "finley/CPPAdapter/MeshAdapter.h"  
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/FinleyAdapterException.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/Data.h"  
 #include "escript/Data/DataArrayView.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/DataFactory.h"  
 #include <iostream>  
 #include <vector>  
 #include <sstream>  
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    
34  //  //
35  // define the statics  // define the static constants
36  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
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  //  //
61  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
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
638    {
639       return ReducedContactElementsOne;
640  }  }
641    
642  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
643  {  {
644    return DegreesOfFreedom;     return DegreesOfFreedom;
645  }  }
646    
647  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
648  {  {
649    return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
 }  
 int MeshAdapter::getDiracDeltaFunctionCode() const  
 {  
   return Points;  
650  }  }
651  //  
652  // returns a pointer to the tag list of samples of functionSpaceType  int MeshAdapter::getDiracDeltaFunctionsCode() const
 //  
 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  //  //
677  // return the number of data points per sample and the number of samples  // return the number of data points per sample and the number of samples
678  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
# Line 335  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: "     }
735           << functionSpaceCode << " for domain: " << getDescription();     break;
736          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
737          break;     if (mesh->ContactElements!=NULL) {
738        }        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
739        return pair<int,int>(numDataPointsPerSample,numSamples);        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.getFinley_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.getFinley_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.getFinley_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();     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();
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();
855    
856       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
857       checkFinleyError();
858    
859       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
860       checkFinleyError();
861    
862       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
863       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();     checkFinleyError();
907    
908     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
909     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     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();     checkFinleyError();
916     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
917     // 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, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
918     checkFinleyError();     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                 Finley_ErrorCode=TYPE_ERROR;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
960                 sprintf(Finley_ErrorMsg,"Interpolation on Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());        break;
961                 break;        case(ContactElementsOne):
962          }        case(ReducedContactElementsOne):
963          break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
964       case(Elements):        break;
965          if (target.getFunctionSpace().getTypeCode()==Elements) {        default:
966             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));           stringstream temp;
967          } else {           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
968             Finley_ErrorCode=TYPE_ERROR;           throw FinleyAdapterException(temp.str());
969             sprintf(Finley_ErrorMsg,"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             Finley_ErrorCode=TYPE_ERROR;        case(DegreesOfFreedom):
977             sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");        case(ReducedDegreesOfFreedom):
978             break;        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
979         }        break;
980       case(Points):        case(Elements):
981          if (target.getFunctionSpace().getTypeCode()==Points) {        case(ReducedElements):
982             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
983          } else {        break;
984             Finley_ErrorCode=TYPE_ERROR;        case(FaceElements):
985             sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");        case(ReducedFaceElements):
986             break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
987          }        break;
988          break;        case(Points):
989       case(ContactElementsZero):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
990       case(ContactElementsOne):        break;
991          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        case(ContactElementsZero):
992             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));        case(ReducedContactElementsZero):
993          } else {        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
994             Finley_ErrorCode=TYPE_ERROR;        break;
995             sprintf(Finley_ErrorMsg,"No interpolation with data on contact elements possible.");        case(ContactElementsOne):
996             break;        case(ReducedContactElementsOne):
997          }        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
998          break;        break;
999       case(DegreesOfFreedom):        default:
1000          switch(target.getFunctionSpace().getTypeCode()) {           stringstream temp;
1001             case(ReducedDegreesOfFreedom):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1002             case(DegreesOfFreedom):           throw FinleyAdapterException(temp.str());
1003             case(Nodes):           break;
1004                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));        }
1005                break;        break;
1006             case(Elements):     case(Elements):
1007                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));        if (target.getFunctionSpace().getTypeCode()==Elements) {
1008                break;           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1009             case(FaceElements):        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1010                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1011                break;        } else {
1012             case(Points):           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1013                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));        }
1014                break;        break;
1015             case(ContactElementsZero):     case(ReducedElements):
1016             case(ContactElementsOne):        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1017                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1018               break;        } else {
1019             default:           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
1020               Finley_ErrorCode=TYPE_ERROR;        }
1021               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());        break;
1022               break;     case(FaceElements):
1023          }        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1024          break;           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1025       case(ReducedDegreesOfFreedom):        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1026         switch(target.getFunctionSpace().getTypeCode()) {           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1027            case(ReducedDegreesOfFreedom):        } else {
1028               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1029               break;        }
1030            case(Elements):        break;
1031               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));     case(ReducedFaceElements):
1032               break;        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1033            case(FaceElements):           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1034               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));        } else {
1035               break;           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1036            case(Points):        }
1037               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));        break;
1038               break;     case(Points):
1039            case(ContactElementsZero):        if (target.getFunctionSpace().getTypeCode()==Points) {
1040            case(ContactElementsOne):           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1041               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        } else {
1042               break;           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1043            case(Nodes):        }
1044               Finley_ErrorCode=TYPE_ERROR;        break;
1045               sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");     case(ContactElementsZero):
1046               break;     case(ContactElementsOne):
1047            case(DegreesOfFreedom):        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1048               Finley_ErrorCode=TYPE_ERROR;           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1049               sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1050               break;           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1051            default:        } else {
1052               Finley_ErrorCode=TYPE_ERROR;           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1053               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());        }
1054               break;        break;
1055         }     case(ReducedContactElementsZero):
1056         break;     case(ReducedContactElementsOne):
1057       default:        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1058          Finley_ErrorCode=TYPE_ERROR;           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1059          sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",in.getFunctionSpace().getTypeCode());        } else {
1060          break;           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1061    }        }
1062    checkFinleyError();        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        Finley_ErrorCode=TYPE_ERROR;     switch(normal.getFunctionSpace().getTypeCode()) {
1242        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");     case(Nodes):
1243        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1244      case(Elements):     break;
1245        Finley_ErrorCode=TYPE_ERROR;     case(ReducedNodes):
1246        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1247        break;     break;
1248      case (FaceElements):     case(Elements):
1249        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1250        break;     break;
1251      case(Points):     case(ReducedElements):
1252        Finley_ErrorCode=TYPE_ERROR;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1253        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for point elements");     break;
1254        break;     case (FaceElements):
1255      case (ContactElementsOne):     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1256      case (ContactElementsZero):     break;
1257        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));     case (ReducedFaceElements):
1258        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1259      case(DegreesOfFreedom):     break;
1260        Finley_ErrorCode=TYPE_ERROR;     case(Points):
1261        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for degrees of freedom.");     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1262        break;     break;
1263      case(ReducedDegreesOfFreedom):     case (ContactElementsOne):
1264        Finley_ErrorCode=TYPE_ERROR;     case (ContactElementsZero):
1265        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for reduced degrees of freedom.");     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1266        break;     break;
1267      default:     case (ReducedContactElementsOne):
1268        Finley_ErrorCode=TYPE_ERROR;     case (ReducedContactElementsZero):
1269        sprintf(Finley_ErrorMsg,"Normal Vectors: Finley does not know anything about function space type %d",normal.getFunctionSpace().getTypeCode());     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
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;
1279          temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1280          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    Finley_ErrorCode=SYSTEM_ERROR;  
1296    sprintf(Finley_ErrorMsg,"Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
   checkFinleyError();  
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          Finley_ErrorCode=TYPE_ERROR;     escript::Data temp;
1312          sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");     escriptDataC _arg=arg.getDataC();
1313          break;     switch(arg.getFunctionSpace().getTypeCode()) {
1314       case(Elements):     case(Nodes):
1315          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);     temp=escript::Data( arg, escript::function(*this) );
1316          break;     _temp=temp.getDataC();
1317       case(FaceElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1318          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);     break;
1319          break;     case(ReducedNodes):
1320       case(Points):     temp=escript::Data( arg, escript::function(*this) );
1321          Finley_ErrorCode=TYPE_ERROR;     _temp=temp.getDataC();
1322          sprintf(Finley_ErrorMsg,"Integral of data on points is not supported.");     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1323          break;     break;
1324       case(ContactElementsZero):     case(Elements):
1325          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1326          break;     break;
1327       case(ContactElementsOne):     case(ReducedElements):
1328          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1329          break;     break;
1330       case(DegreesOfFreedom):     case(FaceElements):
1331          Finley_ErrorCode=TYPE_ERROR;     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1332          sprintf(Finley_ErrorMsg,"Integral of data on degrees of freedom is not supported.");     break;
1333          break;     case(ReducedFaceElements):
1334       case(ReducedDegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1335          Finley_ErrorCode=TYPE_ERROR;     break;
1336          sprintf(Finley_ErrorMsg,"Integral of data on reduced degrees of freedom is not supported.");     case(Points):
1337          break;     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1338       default:     break;
1339          Finley_ErrorCode=TYPE_ERROR;     case(ContactElementsZero):
1340          sprintf(Finley_ErrorMsg,"Integrals: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1341          break;     break;
1342    }     case(ReducedContactElementsZero):
1343    checkFinleyError();     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            Finley_ErrorCode=TYPE_ERROR;     escript::Data temp;
1387            sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");     if (getMPISize()>1) {
1388            break;        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1389         case(Elements):           temp=escript::Data( arg,  continuousFunction(*this) );
1390            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));           nodeDataC = temp.getDataC();
1391            break;        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1392         case(FaceElements):           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1393            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));           nodeDataC = temp.getDataC();
1394            break;        } else {
1395         case(Points):           nodeDataC = arg.getDataC();
1396            Finley_ErrorCode=TYPE_ERROR;        }
1397            sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");     } else {
1398            break;        nodeDataC = arg.getDataC();
1399         case(ContactElementsZero):     }
1400            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));     switch(grad.getFunctionSpace().getTypeCode()) {
1401            break;     case(Nodes):
1402         case(ContactElementsOne):     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1403            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));     break;
1404            break;     case(ReducedNodes):
1405         case(DegreesOfFreedom):     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1406            Finley_ErrorCode=TYPE_ERROR;     break;
1407            sprintf(Finley_ErrorMsg,"Gradient at degrees of freedom is not supported.");     case(Elements):
1408            break;     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1409         case(ReducedDegreesOfFreedom):     break;
1410            Finley_ErrorCode=TYPE_ERROR;     case(ReducedElements):
1411            sprintf(Finley_ErrorMsg,"Gradient at reduced degrees of freedom is not supported.");     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1412            break;     break;
1413         default:     case(FaceElements):
1414            Finley_ErrorCode=TYPE_ERROR;     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1415            sprintf(Finley_ErrorMsg,"Gradient: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());     break;
1416            break;     case(ReducedFaceElements):
1417    }     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1418    checkFinleyError();     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            Finley_ErrorCode=TYPE_ERROR;     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1459            sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");     break;
1460            break;     case(ReducedNodes):
1461         case(Elements):     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1462            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     break;
1463            break;     case(Elements):
1464         case(FaceElements):     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1465            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     break;
1466            break;     case(ReducedElements):
1467         case(Points):     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1468            Finley_ErrorCode=TYPE_ERROR;     break;
1469            sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");     case(FaceElements):
1470            break;     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1471         case(ContactElementsZero):     break;
1472         case(ContactElementsOne):     case(ReducedFaceElements):
1473            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1474            break;     break;
1475         case(DegreesOfFreedom):     case(Points):
1476            Finley_ErrorCode=TYPE_ERROR;     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1477            sprintf(Finley_ErrorMsg,"Size of degrees of freedom is not supported.");     break;
1478            break;     case(ContactElementsZero):
1479         case(ReducedDegreesOfFreedom):     case(ContactElementsOne):
1480            Finley_ErrorCode=TYPE_ERROR;     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1481            sprintf(Finley_ErrorMsg,"Size of reduced degrees of freedom is not supported.");     break;
1482            break;     case(ReducedContactElementsZero):
1483         default:     case(ReducedContactElementsOne):
1484            Finley_ErrorCode=TYPE_ERROR;     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1485            sprintf(Finley_ErrorMsg,"Element size: Finley does not know anything about function space type %d",size.getFunctionSpace().getTypeCode());     break;
1486            break;     case(DegreesOfFreedom):
1487    }     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1488    checkFinleyError();     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  // sets the location of nodes:  
1501  void MeshAdapter::setNewX(const Data& new_x)  //
1502    // 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  // saves a data array in openDX format:         Finley_Mesh_setCoordinates(mesh,&tmp);
1514  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const     } else {
1515  {         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1516    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));         tmp = new_x_inter.getDataC();
1517    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1518  }     }
1519  // saves a data array in openVTK format:     checkFinleyError();
1520  void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const  }
1521  {  
1522    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));  //
1523    checkFinleyError();  // Helper for the save* methods. Extracts optional data variable names and the
1524  }  // corresponding pointers from python dictionary. Caller must free arrays.
1525  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  //
1526  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1527                        const int row_blocksize,  {
1528                        const escript::FunctionSpace& row_functionspace,     numData = boost::python::extract<int>(arg.attr("__len__")());
1529                        const int column_blocksize,     /* win32 refactor */
1530                        const escript::FunctionSpace& column_functionspace,     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1531                        const int type) const     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1532  {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1533      int reduceRowOrder=0;  
1534      int reduceColOrder=0;     boost::python::list keys=arg.keys();
1535      // is the domain right?     for (int i=0; i<numData; ++i) {
1536      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());        string n=boost::python::extract<string>(keys[i]);
1537      if (row_domain!=*this)        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1538            throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1539      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1540      if (col_domain!=*this)        data[i] = d.getDataC();
1541            throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        dataPtr[i] = &(data[i]);
1542      // is the function space type right        names[i] = TMPMEMALLOC(n.length()+1, char);
1543      if (row_functionspace.getTypeCode()==DegreesOfFreedom) {        strcpy(names[i], n.c_str());
1544          reduceRowOrder=0;     }
1545      } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  }
1546          reduceRowOrder=1;  
1547      } else {  //
1548          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");  // saves mesh and optionally data arrays in openDX format
1549    //
1550    void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1551    {
1552       int num_data;
1553       char **names;
1554       escriptDataC *data;
1555       escriptDataC **ptr_data;
1556    
1557       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1558       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1559       checkFinleyError();
1560    
1561       /* win32 refactor */
1562       TMPMEMFREE(data);
1563       TMPMEMFREE(ptr_data);
1564       for(int i=0; i<num_data; i++)
1565       {
1566          TMPMEMFREE(names[i]);
1567       }
1568       TMPMEMFREE(names);
1569    
1570       return;
1571    }
1572    
1573    //
1574    // saves mesh and optionally data arrays in VTK format
1575    //
1576    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1577    {
1578        boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1579        pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1580                  metadata, metadata_schema, arg);
1581    }
1582    
1583    bool MeshAdapter::ownSample(int fs_code, index_t id) const
1584    {
1585    #ifdef ESYS_MPI
1586        index_t myFirstNode=0, myLastNode=0, k=0;
1587        index_t* globalNodeIndex=0;
1588        Finley_Mesh* mesh_p=m_finleyMesh.get();
1589        if (fs_code == FINLEY_REDUCED_NODES)
1590        {
1591        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1592        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1593        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1594      }      }
1595      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      else
1596          reduceColOrder=0;      {
1597      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1598          reduceColOrder=1;      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1599      } else {      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
         throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");  
1600      }      }
1601      // generate matrix:      k=globalNodeIndex[id];
1602            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1603      Finley_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);  #endif
1604      checkFinleyError();      return true;
1605      Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);  }
1606      checkFinleyError();  
1607      Finley_SystemMatrixPattern_dealloc(fsystemMatrixPattern);  
1608      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);  
1609    //
1610    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1611    //
1612    ASM_ptr MeshAdapter::newSystemMatrix(
1613                                                     const int row_blocksize,
1614                                                     const escript::FunctionSpace& row_functionspace,
1615                                                     const int column_blocksize,
1616                                                     const escript::FunctionSpace& column_functionspace,
1617                                                     const int type) const
1618    {
1619       int reduceRowOrder=0;
1620       int reduceColOrder=0;
1621       // is the domain right?
1622       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1623       if (row_domain!=*this)
1624          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1625       const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1626       if (col_domain!=*this)
1627          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1628       // is the function space type right
1629       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1630          reduceRowOrder=0;
1631       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1632          reduceRowOrder=1;
1633       } else {
1634          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1635       }
1636       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1637          reduceColOrder=0;
1638       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1639          reduceColOrder=1;
1640       } else {
1641          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1642       }
1643       // generate matrix:
1644    
1645       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1646       checkFinleyError();
1647       Paso_SystemMatrix* fsystemMatrix;
1648       int trilinos = 0;
1649       if (trilinos) {
1650    #ifdef TRILINOS
1651          /* Allocation Epetra_VrbMatrix here */
1652    #endif
1653       }
1654       else {
1655          fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1656       }
1657       checkPasoError();
1658       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1659       SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1660       return ASM_ptr(sma);
1661    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1662    }
1663    
1664    //
1665    // creates a TransportProblemAdapter
1666    //
1667    ATP_ptr MeshAdapter::newTransportProblem(
1668                                                             const bool useBackwardEuler,
1669                                                             const int blocksize,
1670                                                             const escript::FunctionSpace& functionspace,
1671                                                             const int type) const
1672    {
1673       int reduceOrder=0;
1674       // is the domain right?
1675       const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1676       if (domain!=*this)
1677          throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1678       // is the function space type right
1679       if (functionspace.getTypeCode()==DegreesOfFreedom) {
1680          reduceOrder=0;
1681       } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1682          reduceOrder=1;
1683       } else {
1684          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1685       }
1686       // generate matrix:
1687    
1688       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1689       checkFinleyError();
1690       Paso_TransportProblem* transportProblem;
1691       transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1692       checkPasoError();
1693       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1694       TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1695       return ATP_ptr(tpa);
1696    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1697  }  }
1698    
1699  //  //
1700  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
1701  // TODO:  // TODO:
1702  //  //
1703    
1704  //  //
1705  // 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:
1706  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1707  {  {
1708    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1709         case(Nodes):     case(Nodes):
1710         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1711         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1712            return false;     return false;
1713            break;     break;
1714         case(Elements):     case(Elements):
1715         case(FaceElements):     case(FaceElements):
1716         case(Points):     case(Points):
1717         case(ContactElementsZero):     case(ContactElementsZero):
1718         case(ContactElementsOne):     case(ContactElementsOne):
1719            return true;     case(ReducedElements):
1720            break;     case(ReducedFaceElements):
1721         default:     case(ReducedContactElementsZero):
1722            Finley_ErrorCode=TYPE_ERROR;     case(ReducedContactElementsOne):
1723            sprintf(Finley_ErrorMsg,"Cell: Finley does not know anything about function space type %d",functionSpaceCode);     return true;
1724            break;     break;
1725    }     default:
1726    checkFinleyError();        stringstream temp;
1727    return false;        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1728          throw FinleyAdapterException(temp.str());
1729          break;
1730       }
1731       checkFinleyError();
1732       return false;
1733    }
1734    
1735    bool
1736    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1737    {
1738       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1739        class 1: DOF <-> Nodes
1740        class 2: ReducedDOF <-> ReducedNodes
1741        class 3: Points
1742        class 4: Elements
1743        class 5: ReducedElements
1744        class 6: FaceElements
1745        class 7: ReducedFaceElements
1746        class 8: ContactElementZero <-> ContactElementOne
1747        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1748    
1749       There is also a set of lines. Interpolation is possible down a line but not between lines.
1750       class 1 and 2 belong to all lines so aren't considered.
1751        line 0: class 3
1752        line 1: class 4,5
1753        line 2: class 6,7
1754        line 3: class 8,9
1755    
1756       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1757       eg hasnodes is true if we have at least one instance of Nodes.
1758       */
1759        if (fs.empty())
1760        {
1761            return false;
1762        }
1763        vector<int> hasclass(10);
1764        vector<int> hasline(4);
1765        bool hasnodes=false;
1766        bool hasrednodes=false;
1767        bool hascez=false;
1768        bool hasrcez=false;
1769        for (int i=0;i<fs.size();++i)
1770        {
1771        switch(fs[i])
1772        {
1773        case(Nodes):   hasnodes=true;   // no break is deliberate
1774        case(DegreesOfFreedom):
1775            hasclass[1]=1;
1776            break;
1777        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1778        case(ReducedDegreesOfFreedom):
1779            hasclass[2]=1;
1780            break;
1781        case(Points):
1782            hasline[0]=1;
1783            hasclass[3]=1;
1784            break;
1785        case(Elements):
1786            hasclass[4]=1;
1787            hasline[1]=1;
1788            break;
1789        case(ReducedElements):
1790            hasclass[5]=1;
1791            hasline[1]=1;
1792            break;
1793        case(FaceElements):
1794            hasclass[6]=1;
1795            hasline[2]=1;
1796            break;
1797        case(ReducedFaceElements):
1798            hasclass[7]=1;
1799            hasline[2]=1;
1800            break;
1801        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1802        case(ContactElementsOne):
1803            hasclass[8]=1;
1804            hasline[3]=1;
1805            break;
1806        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1807        case(ReducedContactElementsOne):
1808            hasclass[9]=1;
1809            hasline[3]=1;
1810            break;
1811        default:
1812            return false;
1813        }
1814        }
1815        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1816        // fail if we have more than one leaf group
1817    
1818        if (totlines>1)
1819        {
1820        return false;   // there are at least two branches we can't interpolate between
1821        }
1822        else if (totlines==1)
1823        {
1824        if (hasline[0]==1)      // we have points
1825        {
1826            resultcode=Points;
1827        }
1828        else if (hasline[1]==1)
1829        {
1830            if (hasclass[5]==1)
1831            {
1832            resultcode=ReducedElements;
1833            }
1834            else
1835            {
1836            resultcode=Elements;
1837            }
1838        }
1839        else if (hasline[2]==1)
1840        {
1841            if (hasclass[7]==1)
1842            {
1843            resultcode=ReducedFaceElements;
1844            }
1845            else
1846            {
1847            resultcode=FaceElements;
1848            }
1849        }
1850        else    // so we must be in line3
1851        {
1852            if (hasclass[9]==1)
1853            {
1854            // need something from class 9
1855            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1856            }
1857            else
1858            {
1859            // something from class 8
1860            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1861            }
1862        }
1863        }
1864        else    // totlines==0
1865        {
1866        if (hasclass[2]==1)
1867        {
1868            // something from class 2
1869            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1870        }
1871        else
1872        {   // something from class 1
1873            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1874        }
1875        }
1876        return true;
1877  }  }
1878    
1879  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1880  {  {
1881    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1882       case(Nodes):     case(Nodes):
1883          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1884             case(Nodes):      case(Nodes):
1885             case(ReducedDegreesOfFreedom):      case(ReducedNodes):
1886             case(DegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1887             case(Elements):      case(DegreesOfFreedom):
1888             case(FaceElements):      case(Elements):
1889             case(Points):      case(ReducedElements):
1890             case(ContactElementsZero):      case(FaceElements):
1891             case(ContactElementsOne):      case(ReducedFaceElements):
1892                 return true;      case(Points):
1893             default:      case(ContactElementsZero):
1894                 Finley_ErrorCode=TYPE_ERROR;      case(ReducedContactElementsZero):
1895                 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);      case(ContactElementsOne):
1896          }      case(ReducedContactElementsOne):
1897          break;      return true;
1898       case(Elements):      default:
1899          if (functionSpaceType_target==Elements) {            stringstream temp;
1900             return true;            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1901              throw FinleyAdapterException(temp.str());
1902       }
1903       break;
1904       case(ReducedNodes):
1905        switch(functionSpaceType_target) {
1906        case(ReducedNodes):
1907        case(ReducedDegreesOfFreedom):
1908        case(Elements):
1909        case(ReducedElements):
1910        case(FaceElements):
1911        case(ReducedFaceElements):
1912        case(Points):
1913        case(ContactElementsZero):
1914        case(ReducedContactElementsZero):
1915        case(ContactElementsOne):
1916        case(ReducedContactElementsOne):
1917        return true;
1918        case(Nodes):
1919        case(DegreesOfFreedom):
1920        return false;
1921        default:
1922            stringstream temp;
1923            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1924            throw FinleyAdapterException(temp.str());
1925       }
1926       break;
1927       case(Elements):
1928        if (functionSpaceType_target==Elements) {
1929          return true;
1930        } else if (functionSpaceType_target==ReducedElements) {
1931          return true;
1932          } else {          } else {
1933             return false;            return false;
         }  
      case(FaceElements):  
         if (functionSpaceType_target==FaceElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(Points):  
         if (functionSpaceType_target==Points) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ContactElementsZero):  
      case(ContactElementsOne):  
         if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
            return true;  
         } else {  
            return false;  
         }  
      case(DegreesOfFreedom):  
         switch(functionSpaceType_target) {  
            case(ReducedDegreesOfFreedom):  
            case(DegreesOfFreedom):  
            case(Nodes):  
            case(Elements):  
            case(FaceElements):  
            case(Points):  
            case(ContactElementsZero):  
            case(ContactElementsOne):  
               return true;  
            default:  
              Finley_ErrorCode=TYPE_ERROR;  
              sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);  
1934          }          }
1935          break;     case(ReducedElements):
1936       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1937         switch(functionSpaceType_target) {        return true;
1938            case(ReducedDegreesOfFreedom):      } else {
1939            case(Elements):            return false;
1940            case(FaceElements):      }
1941            case(Points):     case(FaceElements):
1942            case(ContactElementsZero):      if (functionSpaceType_target==FaceElements) {
1943            case(ContactElementsOne):              return true;
1944                return true;      } else if (functionSpaceType_target==ReducedFaceElements) {
1945            case(Nodes):              return true;
1946            case(DegreesOfFreedom):      } else {
1947               return false;              return false;
1948            default:      }
1949               Finley_ErrorCode=TYPE_ERROR;     case(ReducedFaceElements):
1950               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);      if (functionSpaceType_target==ReducedFaceElements) {
1951         }              return true;
1952         break;      } else {
1953       default:          return false;
1954          Finley_ErrorCode=TYPE_ERROR;      }
1955          sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_source);     case(Points):
1956          break;      if (functionSpaceType_target==Points) {
1957    }              return true;
1958    checkFinleyError();      } else {
1959    return false;              return false;
1960        }
1961       case(ContactElementsZero):
1962       case(ContactElementsOne):
1963        if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1964                return true;
1965        } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1966                return true;
1967        } else {
1968                return false;
1969        }
1970       case(ReducedContactElementsZero):
1971       case(ReducedContactElementsOne):
1972        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1973                return true;
1974        } else {
1975                return false;
1976        }
1977       case(DegreesOfFreedom):
1978        switch(functionSpaceType_target) {
1979        case(ReducedDegreesOfFreedom):
1980        case(DegreesOfFreedom):
1981        case(Nodes):
1982        case(ReducedNodes):
1983        case(Elements):
1984        case(ReducedElements):
1985        case(Points):
1986        case(FaceElements):
1987        case(ReducedFaceElements):
1988        case(ContactElementsZero):
1989        case(ReducedContactElementsZero):
1990        case(ContactElementsOne):
1991        case(ReducedContactElementsOne):
1992        return true;
1993        default:
1994            stringstream temp;
1995            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1996            throw FinleyAdapterException(temp.str());
1997        }
1998        break;
1999       case(ReducedDegreesOfFreedom):
2000       switch(functionSpaceType_target) {
2001        case(ReducedDegreesOfFreedom):
2002        case(ReducedNodes):
2003        case(Elements):
2004        case(ReducedElements):
2005        case(FaceElements):
2006        case(ReducedFaceElements):
2007        case(Points):
2008        case(ContactElementsZero):
2009        case(ReducedContactElementsZero):
2010        case(ContactElementsOne):
2011        case(ReducedContactElementsOne):
2012        return true;
2013        case(Nodes):
2014        case(DegreesOfFreedom):
2015        return false;
2016        default:
2017            stringstream temp;
2018            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2019            throw FinleyAdapterException(temp.str());
2020        }
2021        break;
2022       default:
2023          stringstream temp;
2024          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2025          throw FinleyAdapterException(temp.str());
2026          break;
2027       }
2028       checkFinleyError();
2029       return false;
2030  }  }
2031    
2032  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
2033  {  {
2034      return false;     return false;
2035  }  }
2036    
2037  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2038  {  {
2039    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2040    if (temp!=0) {     if (temp!=0) {
2041      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2042    } else {     } else {
2043      return false;        return false;
2044    }     }
2045  }  }
2046    
2047  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2048  {  {
2049    return !(operator==(other));     return !(operator==(other));
2050  }  }
2051    
2052  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2053  {  {
2054     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2055     checkFinleyError();     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2056     return out;             package, symmetry, mesh->MPIInfo);
2057  }  }
2058  Data MeshAdapter::getX() const  
2059    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2060  {  {
2061    return continuousFunction(asAbstractContinuousDomain()).getX();     Finley_Mesh* mesh=m_finleyMesh.get();
2062       return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2063               package, symmetry, mesh->MPIInfo);
2064  }  }
2065  Data MeshAdapter::getNormal() const  
2066    escript::Data MeshAdapter::getX() const
2067    {
2068       return continuousFunction(*this).getX();
2069    }
2070    
2071    escript::Data MeshAdapter::getNormal() const
2072  {  {
2073    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2074  }  }
2075  Data MeshAdapter::getSize() const  
2076    escript::Data MeshAdapter::getSize() const
2077  {  {
2078    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2079  }  }
2080    
2081    const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2082    {
2083       int *out = NULL;
2084       Finley_Mesh* mesh=m_finleyMesh.get();
2085       switch (functionSpaceType) {
2086       case(Nodes):
2087       out=mesh->Nodes->Id;
2088       break;
2089       case(ReducedNodes):
2090       out=mesh->Nodes->reducedNodesId;
2091       break;
2092       case(Elements):
2093       out=mesh->Elements->Id;
2094       break;
2095       case(ReducedElements):
2096       out=mesh->Elements->Id;
2097       break;
2098       case(FaceElements):
2099       out=mesh->FaceElements->Id;
2100       break;
2101       case(ReducedFaceElements):
2102       out=mesh->FaceElements->Id;
2103       break;
2104       case(Points):
2105       out=mesh->Points->Id;
2106       break;
2107       case(ContactElementsZero):
2108       out=mesh->ContactElements->Id;
2109       break;
2110       case(ReducedContactElementsZero):
2111       out=mesh->ContactElements->Id;
2112       break;
2113       case(ContactElementsOne):
2114       out=mesh->ContactElements->Id;
2115       break;
2116       case(ReducedContactElementsOne):
2117       out=mesh->ContactElements->Id;
2118       break;
2119       case(DegreesOfFreedom):
2120       out=mesh->Nodes->degreesOfFreedomId;
2121       break;
2122       case(ReducedDegreesOfFreedom):
2123       out=mesh->Nodes->reducedDegreesOfFreedomId;
2124       break;
2125       default:
2126          stringstream temp;
2127          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2128          throw FinleyAdapterException(temp.str());
2129          break;
2130       }
2131       return out;
2132    }
2133  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2134  {  {
2135    int* tagList;     int out=0;
2136    int numTags;     Finley_Mesh* mesh=m_finleyMesh.get();
2137    getTagList(functionSpaceType, &tagList, &numTags);     switch (functionSpaceType) {
2138    return tagList[sampleNo];     case(Nodes):
2139       out=mesh->Nodes->Tag[sampleNo];
2140       break;
2141       case(ReducedNodes):
2142       throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2143       break;
2144       case(Elements):
2145       out=mesh->Elements->Tag[sampleNo];
2146       break;
2147       case(ReducedElements):
2148       out=mesh->Elements->Tag[sampleNo];
2149       break;
2150       case(FaceElements):
2151       out=mesh->FaceElements->Tag[sampleNo];
2152       break;
2153       case(ReducedFaceElements):
2154       out=mesh->FaceElements->Tag[sampleNo];
2155       break;
2156       case(Points):
2157       out=mesh->Points->Tag[sampleNo];
2158       break;
2159       case(ContactElementsZero):
2160       out=mesh->ContactElements->Tag[sampleNo];
2161       break;
2162       case(ReducedContactElementsZero):
2163       out=mesh->ContactElements->Tag[sampleNo];
2164       break;
2165       case(ContactElementsOne):
2166       out=mesh->ContactElements->Tag[sampleNo];
2167       break;
2168       case(ReducedContactElementsOne):
2169       out=mesh->ContactElements->Tag[sampleNo];
2170       break;
2171       case(DegreesOfFreedom):
2172       throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2173       break;
2174       case(ReducedDegreesOfFreedom):
2175       throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2176       break;
2177       default:
2178          stringstream temp;
2179          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2180          throw FinleyAdapterException(temp.str());
2181          break;
2182       }
2183       return out;
2184    }
2185    
2186    
2187    void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2188    {
2189       Finley_Mesh* mesh=m_finleyMesh.get();
2190       escriptDataC tmp=mask.getDataC();
2191       switch(functionSpaceType) {
2192       case(Nodes):
2193       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2194       break;
2195       case(ReducedNodes):
2196       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2197       break;
2198       case(DegreesOfFreedom):
2199       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2200       break;
2201       case(ReducedDegreesOfFreedom):
2202       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2203       break;
2204       case(Elements):
2205       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2206       break;
2207       case(ReducedElements):
2208       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2209       break;
2210       case(FaceElements):
2211       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2212       break;
2213       case(ReducedFaceElements):
2214       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2215       break;
2216       case(Points):
2217       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2218       break;
2219       case(ContactElementsZero):
2220       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2221       break;
2222       case(ReducedContactElementsZero):
2223       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2224       break;
2225       case(ContactElementsOne):
2226       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2227       break;
2228       case(ReducedContactElementsOne):
2229       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2230       break;
2231       default:
2232          stringstream temp;
2233          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2234          throw FinleyAdapterException(temp.str());
2235       }
2236       checkFinleyError();
2237       return;
2238    }
2239    
2240    void MeshAdapter::setTagMap(const string& name,  int tag)
2241    {
2242       Finley_Mesh* mesh=m_finleyMesh.get();
2243       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2244       checkFinleyError();
2245       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2246    }
2247    
2248    int MeshAdapter::getTag(const string& name) const
2249    {
2250       Finley_Mesh* mesh=m_finleyMesh.get();
2251       int tag=0;
2252       tag=Finley_Mesh_getTag(mesh, name.c_str());
2253       checkFinleyError();
2254       // throwStandardException("MeshAdapter::getTag is not implemented.");
2255       return tag;
2256    }
2257    
2258    bool MeshAdapter::isValidTagName(const string& name) const
2259    {
2260       Finley_Mesh* mesh=m_finleyMesh.get();
2261       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2262    }
2263    
2264    string MeshAdapter::showTagNames() const
2265    {
2266       stringstream temp;
2267       Finley_Mesh* mesh=m_finleyMesh.get();
2268       Finley_TagMap* tag_map=mesh->TagMap;
2269       while (tag_map) {
2270          temp << tag_map->name;
2271          tag_map=tag_map->next;
2272          if (tag_map) temp << ", ";
2273       }
2274       return temp.str();
2275    }
2276    
2277    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2278    {
2279      Finley_Mesh* mesh=m_finleyMesh.get();
2280      dim_t numTags=0;
2281      switch(functionSpaceCode) {
2282       case(Nodes):
2283              numTags=mesh->Nodes->numTagsInUse;
2284              break;
2285       case(ReducedNodes):
2286              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2287              break;
2288       case(DegreesOfFreedom):
2289              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2290              break;
2291       case(ReducedDegreesOfFreedom):
2292              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2293              break;
2294       case(Elements):
2295       case(ReducedElements):
2296              numTags=mesh->Elements->numTagsInUse;
2297              break;
2298       case(FaceElements):
2299       case(ReducedFaceElements):
2300              numTags=mesh->FaceElements->numTagsInUse;
2301              break;
2302       case(Points):
2303              numTags=mesh->Points->numTagsInUse;
2304              break;
2305       case(ContactElementsZero):
2306       case(ReducedContactElementsZero):
2307       case(ContactElementsOne):
2308       case(ReducedContactElementsOne):
2309              numTags=mesh->ContactElements->numTagsInUse;
2310              break;
2311       default:
2312          stringstream temp;
2313          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2314          throw FinleyAdapterException(temp.str());
2315      }
2316      return numTags;
2317    }
2318    
2319    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2320    {
2321      Finley_Mesh* mesh=m_finleyMesh.get();
2322      index_t* tags=NULL;
2323      switch(functionSpaceCode) {
2324       case(Nodes):
2325              tags=mesh->Nodes->tagsInUse;
2326              break;
2327       case(ReducedNodes):
2328              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2329              break;
2330       case(DegreesOfFreedom):
2331              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2332              break;
2333       case(ReducedDegreesOfFreedom):
2334              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2335              break;
2336       case(Elements):
2337       case(ReducedElements):
2338              tags=mesh->Elements->tagsInUse;
2339              break;
2340       case(FaceElements):
2341       case(ReducedFaceElements):
2342              tags=mesh->FaceElements->tagsInUse;
2343              break;
2344       case(Points):
2345              tags=mesh->Points->tagsInUse;
2346              break;
2347       case(ContactElementsZero):
2348       case(ReducedContactElementsZero):
2349       case(ContactElementsOne):
2350       case(ReducedContactElementsOne):
2351              tags=mesh->ContactElements->tagsInUse;
2352              break;
2353       default:
2354          stringstream temp;
2355          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2356          throw FinleyAdapterException(temp.str());
2357      }
2358      return tags;
2359    }
2360    
2361    
2362    bool MeshAdapter::canTag(int functionSpaceCode) const
2363    {
2364      switch(functionSpaceCode) {
2365       case(Nodes):
2366       case(Elements):
2367       case(ReducedElements):
2368       case(FaceElements):
2369       case(ReducedFaceElements):
2370       case(Points):
2371       case(ContactElementsZero):
2372       case(ReducedContactElementsZero):
2373       case(ContactElementsOne):
2374       case(ReducedContactElementsOne):
2375              return true;
2376       case(ReducedNodes):
2377       case(DegreesOfFreedom):
2378       case(ReducedDegreesOfFreedom):
2379          return false;
2380       default:
2381        return false;
2382      }
2383    }
2384    
2385    AbstractDomain::StatusType MeshAdapter::getStatus() const
2386    {
2387      Finley_Mesh* mesh=m_finleyMesh.get();
2388      return Finley_Mesh_getStatus(mesh);
2389    }
2390    
2391    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2392    {
2393      
2394      Finley_Mesh* mesh=m_finleyMesh.get();
2395      int order =-1;
2396      switch(functionSpaceCode) {
2397       case(Nodes):
2398       case(DegreesOfFreedom):
2399              order=mesh->approximationOrder;
2400              break;
2401       case(ReducedNodes):
2402       case(ReducedDegreesOfFreedom):
2403              order=mesh->reducedApproximationOrder;
2404              break;
2405       case(Elements):
2406       case(FaceElements):
2407       case(Points):
2408       case(ContactElementsZero):
2409       case(ContactElementsOne):
2410              order=mesh->integrationOrder;
2411              break;
2412       case(ReducedElements):
2413       case(ReducedFaceElements):
2414       case(ReducedContactElementsZero):
2415       case(ReducedContactElementsOne):
2416              order=mesh->reducedIntegrationOrder;
2417              break;
2418       default:
2419          stringstream temp;
2420          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2421          throw FinleyAdapterException(temp.str());
2422      }
2423      return order;
2424    }
2425    
2426    bool MeshAdapter::supportsContactElements() const
2427    {
2428      return true;
2429    }
2430    
2431    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2432    {
2433      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2434    }
2435    
2436    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2437    {
2438      Finley_ReferenceElementSet_dealloc(m_refSet);
2439  }  }
2440    
2441  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  // points will be flattened
2442  {  void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2443    int* referenceNoList;  {
2444    int numReferenceNo;        const int dim = getDim();
2445    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);        int numPoints=points.size()/dim;
2446    return referenceNoList[sampleNo];        int numTags=tags.size();
2447          Finley_Mesh* mesh=m_finleyMesh.get();
2448          
2449          if ( points.size() % dim != 0 )
2450          {
2451        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2452          }
2453          
2454          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2455         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2456          
2457          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2458          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2459          
2460          for (int i=0;i<numPoints;++i) {
2461           points_ptr[ i * dim     ] = points[i * dim ];
2462           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2463           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2464               tags_ptr[i]=tags[i];
2465          }
2466          
2467          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2468          checkFinleyError();
2469          
2470          TMPMEMFREE(points_ptr);
2471          TMPMEMFREE(tags_ptr);
2472  }  }
2473    
2474    
2475    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2476    // {
2477    //       const int dim = getDim();
2478    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2479    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2480    //       Finley_Mesh* mesh=m_finleyMesh.get();
2481    //      
2482    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2483    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2484    //      
2485    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2486    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2487    //      
2488    //       for (int i=0;i<numPoints;++i) {
2489    //     int tag_id=-1;
2490    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2491    //     if  ( numComps !=   dim ) {
2492    //                stringstream temp;            
2493    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2494    //                throw FinleyAdapterException(temp.str());
2495    //     }
2496    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2497    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2498    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2499    //    
2500    //     if ( numTags > 0) {
2501    //            boost::python::extract<string> ex_str(tags[i]);
2502    //        if  ( ex_str.check() ) {
2503    //            tag_id=getTag( ex_str());
2504    //        } else {
2505    //             boost::python::extract<int> ex_int(tags[i]);
2506    //             if ( ex_int.check() ) {
2507    //                 tag_id=ex_int();
2508    //             } else {
2509    //              stringstream temp;          
2510    //                  temp << "Error - unable to extract tag for point " << i;
2511    //              throw FinleyAdapterException(temp.str());
2512    //            }
2513    //        }
2514    //     }      
2515    //            tags_ptr[i]=tag_id;
2516    //       }
2517    //      
2518    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2519    //       checkPasoError();
2520    //      
2521    //       TMPMEMFREE(points_ptr);
2522    //       TMPMEMFREE(tags_ptr);
2523    // }
2524    
2525    /*
2526    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2527    {  
2528        boost::python::list points =  boost::python::list();
2529        boost::python::list tags =  boost::python::list();
2530        points.append(point);
2531        tags.append(tag);
2532        addDiracPoints(points, tags);
2533    }
2534    */
2535    
2536    /*
2537    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2538    {
2539            boost::python::list points =   boost::python::list();
2540            boost::python::list tags =   boost::python::list();
2541            points.append(point);
2542            tags.append(tag);
2543            addDiracPoints(points, tags);
2544    }
2545    */
2546  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.148  
changed lines
  Added in v.3736

  ViewVC Help
Powered by ViewVC 1.1.26