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

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

Legend:
Removed from v.682  
changed lines
  Added in v.4255

  ViewVC Help
Powered by ViewVC 1.1.26