/[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 1116 by gross, Mon Apr 23 03:36:34 2007 UTC revision 3998 by caltinay, Thu Sep 27 01:17:28 2012 UTC
# Line 1  Line 1 
 // $Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2  #ifdef PASO_MPI  /*****************************************************************************
3  #include <mpi.h>  *
4  #endif  * Copyright (c) 2003-2012 by University of Queensland
5  #include "MeshAdapter.h"  * 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 45  const int MeshAdapter::ReducedContactEle Line 52  const int MeshAdapter::ReducedContactEle
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 57  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());     }
     //   cout << "Finished dealloc." << endl;  
   }  
78  }  }
79    
80    int MeshAdapter::getMPISize() const
81    {
82       return m_finleyMesh.get()->MPIInfo->size;
83    }
84    int MeshAdapter::getMPIRank() const
85    {
86       return m_finleyMesh.get()->MPIInfo->rank;
87    }
88    void MeshAdapter::MPIBarrier() const
89    {
90    #ifdef ESYS_MPI
91       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
92    #endif
93       return;
94    }
95    bool MeshAdapter::onMasterProcessor() const
96    {
97       return m_finleyMesh.get()->MPIInfo->rank == 0;
98    }
99    
100    
101    #ifdef ESYS_MPI
102      MPI_Comm
103    #else
104      unsigned int
105    #endif
106    MeshAdapter::getMPIComm() const
107    {
108    #ifdef ESYS_MPI
109        return m_finleyMesh->MPIInfo->comm;
110    #else
111        return 0;
112    #endif
113    }
114    
115    
116  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
117     return m_finleyMesh.get();     return m_finleyMesh.get();
118  }  }
119    
120  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const string& fileName) const
121    {
122       char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
123       strcpy(fName,fileName.c_str());
124       Finley_Mesh_write(m_finleyMesh.get(),fName);
125       checkFinleyError();
126       TMPMEMFREE(fName);
127    }
128    
129    void MeshAdapter::Print_Mesh_Info(const bool full) const
130  {  {
131    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
132    strcpy(fName,fileName.c_str());  }
133    Finley_Mesh_write(m_finleyMesh.get(),fName);  
134    checkFinleyError();  void MeshAdapter::dump(const string& fileName) const
135    TMPMEMFREE(fName);  {
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(ReducedNodes,"Finley_Reduced_Nodes"));     (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
579    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
580      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));     (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
581    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
582      (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
583    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
584      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
585    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
586      (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));     (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
587    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
588      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));     (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
589    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
590      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
591    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
592      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
593    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
594      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));     (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
595    m_functionSpaceTypeNames.insert     m_functionSpaceTypeNames.insert
596      (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));     (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  int MeshAdapter::getReducedContinuousFunctionCode() const
604  {  {
605    return ReducedNodes;     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  int MeshAdapter::getReducedFunctionCode() const
613  {  {
614    return ReducedElements;     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  int MeshAdapter::getReducedFunctionOnBoundaryCode() const
622  {  {
623    return ReducedFaceElements;     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  int MeshAdapter::getReducedFunctionOnContactZeroCode() const
631  {  {
632    return ReducedContactElementsZero;     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  int MeshAdapter::getReducedFunctionOnContactOneCode() const
640  {  {
641    return ReducedContactElementsOne;     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 203  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 218  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(ReducedNodes):     case(ReducedNodes):
693             /* TODO: add ReducedNodes */     numDataPointsPerSample=1;
694             numDataPointsPerSample=1;     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
695             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;     break;
696             throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");     case(Elements):
697             break;     if (mesh->Elements!=NULL) {
698        case(Elements):        numSamples=mesh->Elements->numElements;
699             if (mesh->Elements!=NULL) {        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
700               numSamples=mesh->Elements->numElements;     }
701               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;     break;
702             }     case(ReducedElements):
703             break;     if (mesh->Elements!=NULL) {
704        case(ReducedElements):        numSamples=mesh->Elements->numElements;
705             if (mesh->Elements!=NULL) {        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
706               numSamples=mesh->Elements->numElements;     }
707               numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;     break;
708             }     case(FaceElements):
709             break;     if (mesh->FaceElements!=NULL) {
710        case(FaceElements):        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
711             if (mesh->FaceElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
712                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;     }
713                  numSamples=mesh->FaceElements->numElements;     break;
714             }     case(ReducedFaceElements):
715             break;     if (mesh->FaceElements!=NULL) {
716        case(ReducedFaceElements):        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
717             if (mesh->FaceElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
718                  numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;     }
719                  numSamples=mesh->FaceElements->numElements;     break;
720             }     case(Points):
721             break;     if (mesh->Points!=NULL) {
722        case(Points):        numDataPointsPerSample=1;
723             if (mesh->Points!=NULL) {        numSamples=mesh->Points->numElements;
724               numDataPointsPerSample=1;     }
725               numSamples=mesh->Points->numElements;     break;
726             }     case(ContactElementsZero):
727             break;     if (mesh->ContactElements!=NULL) {
728        case(ContactElementsZero):        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
729             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
730               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
731               numSamples=mesh->ContactElements->numElements;     break;
732             }     case(ReducedContactElementsZero):
733             break;     if (mesh->ContactElements!=NULL) {
734        case(ReducedContactElementsZero):        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
735             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
736               numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;     }
737               numSamples=mesh->ContactElements->numElements;     break;
738             }     case(ContactElementsOne):
739             break;     if (mesh->ContactElements!=NULL) {
740        case(ContactElementsOne):        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
741             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
742               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
743               numSamples=mesh->ContactElements->numElements;     break;
744             }     case(ReducedContactElementsOne):
745             break;     if (mesh->ContactElements!=NULL) {
746        case(ReducedContactElementsOne):        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
747             if (mesh->ContactElements!=NULL) {        numSamples=mesh->ContactElements->numElements;
748               numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;     }
749               numSamples=mesh->ContactElements->numElements;     break;
750             }     case(DegreesOfFreedom):
751             break;     if (mesh->Nodes!=NULL) {
752        case(DegreesOfFreedom):        numDataPointsPerSample=1;
753             if (mesh->Nodes!=NULL) {        numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
754               numDataPointsPerSample=1;     }
755  #ifndef PASO_MPI     break;
756               numSamples=mesh->Nodes->numDegreesOfFreedom;     case(ReducedDegreesOfFreedom):
757  #else     if (mesh->Nodes!=NULL) {
758               numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;        numDataPointsPerSample=1;
759  #endif        numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
760             }     }
761             break;     break;
762        case(ReducedDegreesOfFreedom):     default:
763             if (mesh->Nodes!=NULL) {        stringstream temp;
764               numDataPointsPerSample=1;        temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
765  #ifndef PASO_MPI        throw FinleyAdapterException(temp.str());
766               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;        break;
767  #else     }
768               numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;     return pair<int,int>(numDataPointsPerSample,numSamples);
 #endif  
            }  
            break;  
       default:  
         stringstream temp;  
         temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();  
         throw FinleyAdapterException(temp.str());  
         break;  
       }  
       return pair<int,int>(numDataPointsPerSample,numSamples);  
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    {
781       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
782       if (smat==0)
783       {
784        throw FinleyAdapterException("finley only supports Paso system matrices.");
785       }
786     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
787     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
788     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 331  void MeshAdapter::addPDEToSystem( Line 794  void MeshAdapter::addPDEToSystem(
794     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
795     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
796     escriptDataC _y_contact=y_contact.getDataC();     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();     Finley_Mesh* mesh=m_finleyMesh.get();
801    
802     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
803     checkFinleyError();     checkFinleyError();
804    
805     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
806     checkFinleyError();     checkFinleyError();
807    
808     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
809     checkFinleyError();     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();
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, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  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
845  {  {
846     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
847    
# Line 356  void MeshAdapter::addPDEToRHS( escript:: Line 850  void MeshAdapter::addPDEToRHS( escript::
850     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
851     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
852     escriptDataC _y_contact=y_contact.getDataC();     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 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
856     checkFinleyError();     checkFinleyError();
# Line 365  void MeshAdapter::addPDEToRHS( escript:: Line 860  void MeshAdapter::addPDEToRHS( escript::
860    
861     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
862     checkFinleyError();     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();
909    
910       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
911       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();
918    
919       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
920       checkFinleyError();
921    
922  }  }
923    
924  //  //
# Line 372  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    escriptDataC _target=target.getDataC();     escriptDataC _target=target.getDataC();
938    escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
939    switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
940       case(Nodes):     case(Nodes):
941          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
942             case(Nodes):        case(Nodes):
943             case(ReducedNodes):        case(ReducedNodes):
944             case(DegreesOfFreedom):        case(DegreesOfFreedom):
945             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
946                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
947                 break;        break;
948             case(Elements):        case(Elements):
949             case(ReducedElements):        case(ReducedElements):
950                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
951                 break;        break;
952             case(FaceElements):        case(FaceElements):
953             case(ReducedFaceElements):        case(ReducedFaceElements):
954                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
955                 break;        break;
956             case(Points):        case(Points):
957                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
958                 break;        break;
959             case(ContactElementsZero):        case(ContactElementsZero):
960             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
961                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
962                 break;        break;
963             case(ContactElementsOne):        case(ContactElementsOne):
964             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
965                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
966                 break;        break;
967             default:        default:
968                 stringstream temp;           stringstream temp;
969                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
970                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
971                 break;           break;
972          }        }
973          break;        break;
974       case(ReducedNodes):     case(ReducedNodes):
975          switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
976             case(Nodes):        case(Nodes):
977             case(ReducedNodes):        case(ReducedNodes):
978             case(DegreesOfFreedom):        case(DegreesOfFreedom):
979             case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
980                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
981                 break;        break;
982             case(Elements):        case(Elements):
983             case(ReducedElements):        case(ReducedElements):
984                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
985                 break;        break;
986             case(FaceElements):        case(FaceElements):
987             case(ReducedFaceElements):        case(ReducedFaceElements):
988                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
989                 break;        break;
990             case(Points):        case(Points):
991                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
992                 break;        break;
993             case(ContactElementsZero):        case(ContactElementsZero):
994             case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
995                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
996                 break;        break;
997             case(ContactElementsOne):        case(ContactElementsOne):
998             case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
999                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1000                 break;        break;
1001             default:        default:
1002                 stringstream temp;           stringstream temp;
1003                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1004                 throw FinleyAdapterException(temp.str());           throw FinleyAdapterException(temp.str());
1005                 break;           break;
1006          }        }
1007          break;        break;
1008       case(Elements):     case(Elements):
1009          if (target.getFunctionSpace().getTypeCode()==Elements) {        if (target.getFunctionSpace().getTypeCode()==Elements) {
1010             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1011          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1012             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
1013          } else {        } else {
1014             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
1015          }        }
1016          break;        break;
1017       case(ReducedElements):     case(ReducedElements):
1018          if (target.getFunctionSpace().getTypeCode()==ReducedElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
1019             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
1020          } else {        } else {
1021             throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
1022          }        }
1023          break;        break;
1024       case(FaceElements):     case(FaceElements):
1025          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
1026             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1027          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1028             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1029          } else {        } else {
1030             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1031         }        }
1032         break;        break;
1033       case(ReducedFaceElements):     case(ReducedFaceElements):
1034          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1035             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1036          } else {        } else {
1037             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1038         }        }
1039         break;        break;
1040       case(Points):     case(Points):
1041          if (target.getFunctionSpace().getTypeCode()==Points) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1042             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1043          } else {        } else {
1044             throw FinleyAdapterException("Error - No interpolation with data on points possible.");           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1045          }        }
1046          break;        break;
1047       case(ContactElementsZero):     case(ContactElementsZero):
1048       case(ContactElementsOne):     case(ContactElementsOne):
1049          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1050             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1051          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1052             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1053          } else {        } else {
1054             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1055             break;        }
1056          }        break;
1057          break;     case(ReducedContactElementsZero):
1058       case(ReducedContactElementsZero):     case(ReducedContactElementsOne):
1059       case(ReducedContactElementsOne):        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1060          if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1061             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);        } else {
1062          } else {           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1063             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");        }
1064             break;        break;
1065          }     case(DegreesOfFreedom):      
1066          break;        switch(target.getFunctionSpace().getTypeCode()) {
1067       case(DegreesOfFreedom):              case(ReducedDegreesOfFreedom):
1068          switch(target.getFunctionSpace().getTypeCode()) {        case(DegreesOfFreedom):
1069             case(ReducedDegreesOfFreedom):        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1070             case(DegreesOfFreedom):        break;
1071             case(Nodes):    
1072             case(ReducedNodes):        case(Nodes):
1073                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(ReducedNodes):
1074                break;        if (getMPISize()>1) {
1075  #ifndef PASO_MPI           escript::Data temp=escript::Data(in);
1076             case(Elements):           temp.expand();
1077             case(ReducedElements):           escriptDataC _in2 = temp.getDataC();
1078                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1079                break;        } else {
1080             case(FaceElements):           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1081             case(ReducedFaceElements):        }
1082                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
1083                break;        case(Elements):
1084             case(Points):        case(ReducedElements):
1085                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        if (getMPISize()>1) {
1086                break;           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1087             case(ContactElementsZero):           escriptDataC _in2 = temp.getDataC();
1088             case(ContactElementsOne):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1089             case(ReducedContactElementsZero):        } else {
1090             case(ReducedContactElementsOne):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1091                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        }
1092               break;        break;
1093  #else        case(FaceElements):
1094             /* need to copy Degrees of freedom data to nodal data so that the external values are available */        case(ReducedFaceElements):
1095             case(Elements):        if (getMPISize()>1) {
1096             case(ReducedElements):           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1097                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );           escriptDataC _in2 = temp.getDataC();
1098                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1099                break;    
1100             case(FaceElements):        } else {
1101             case(ReducedFaceElements):           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1102                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );        }
1103                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);        break;
1104                break;        case(Points):
1105             case(Points):        if (getMPISize()>1) {
1106                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1107                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);           //escriptDataC _in2 = temp.getDataC();
1108                break;        } else {
1109             case(ContactElementsZero):           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1110             case(ContactElementsOne):        }
1111             case(ReducedContactElementsZero):        break;
1112             case(ReducedContactElementsOne):        case(ContactElementsZero):
1113                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );        case(ContactElementsOne):
1114                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);        case(ReducedContactElementsZero):
1115               break;        case(ReducedContactElementsOne):
1116  #endif        if (getMPISize()>1) {
1117             default:           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1118               stringstream temp;           escriptDataC _in2 = temp.getDataC();
1119               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1120               throw FinleyAdapterException(temp.str());        } else {
1121               break;           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1122          }        }
1123          break;        break;
1124       case(ReducedDegreesOfFreedom):        default:
1125         switch(target.getFunctionSpace().getTypeCode()) {           stringstream temp;
1126            case(Nodes):           temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1127               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           throw FinleyAdapterException(temp.str());
1128               break;           break;
1129            case(ReducedNodes):        }
1130               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        break;
1131               break;     case(ReducedDegreesOfFreedom):
1132            case(DegreesOfFreedom):        switch(target.getFunctionSpace().getTypeCode()) {
1133               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");        case(Nodes):
1134               break;        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1135            case(ReducedDegreesOfFreedom):        break;
1136               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(ReducedNodes):
1137               break;        if (getMPISize()>1) {
1138            case(Elements):           escript::Data temp=escript::Data(in);
1139            case(ReducedElements):           temp.expand();
1140               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);           escriptDataC _in2 = temp.getDataC();
1141               break;           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1142            case(FaceElements):        } else {
1143            case(ReducedFaceElements):           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1144               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        }
1145               break;        break;
1146            case(Points):        case(DegreesOfFreedom):
1147               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1148               break;        break;
1149            case(ContactElementsZero):        case(ReducedDegreesOfFreedom):
1150            case(ContactElementsOne):        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1151            case(ReducedContactElementsZero):        break;
1152            case(ReducedContactElementsOne):        case(Elements):
1153               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        case(ReducedElements):
1154               break;        if (getMPISize()>1) {
1155            default:           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1156               stringstream temp;           escriptDataC _in2 = temp.getDataC();
1157               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1158               throw FinleyAdapterException(temp.str());        } else {
1159               break;           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1160         }        }
1161         break;        break;
1162       default:        case(FaceElements):
1163          stringstream temp;        case(ReducedFaceElements):
1164          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        if (getMPISize()>1) {
1165          throw FinleyAdapterException(temp.str());           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1166          break;           escriptDataC _in2 = temp.getDataC();
1167    }           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1168    checkFinleyError();        } 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 619  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       escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1221       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1222    } else {     } else {
1223       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1224       escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1225       Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1226       // this is then interpolated onto arg:        // this is then interpolated onto arg:
1227       interpolateOnDomain(arg,tmp_data);        interpolateOnDomain(arg,tmp_data);
1228    }     }
1229    checkFinleyError();     checkFinleyError();
1230  }  }
1231    
1232  //  //
# Line 642  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    escriptDataC _normal=normal.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1242    switch(normal.getFunctionSpace().getTypeCode()) {     escriptDataC _normal=normal.getDataC();
1243      case(Nodes):     switch(normal.getFunctionSpace().getTypeCode()) {
1244        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     case(Nodes):
1245        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1246      case(ReducedNodes):     break;
1247        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");     case(ReducedNodes):
1248        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1249      case(Elements):     break;
1250        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     case(Elements):
1251        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1252      case(ReducedElements):     break;
1253        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");     case(ReducedElements):
1254        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1255      case (FaceElements):     break;
1256        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (FaceElements):
1257        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1258      case (ReducedFaceElements):     break;
1259        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);     case (ReducedFaceElements):
1260        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1261      case(Points):     break;
1262        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     case(Points):
1263        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1264      case (ContactElementsOne):     break;
1265      case (ContactElementsZero):     case (ContactElementsOne):
1266        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ContactElementsZero):
1267        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1268      case (ReducedContactElementsOne):     break;
1269      case (ReducedContactElementsZero):     case (ReducedContactElementsOne):
1270        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);     case (ReducedContactElementsZero):
1271        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1272      case(DegreesOfFreedom):     break;
1273        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case(DegreesOfFreedom):
1274        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1275      case(ReducedDegreesOfFreedom):     break;
1276        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(ReducedDegreesOfFreedom):
1277        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1278      default:     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 697  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    escriptDataC _arg=arg.getDataC();     Finley_Mesh* mesh=m_finleyMesh.get();
1312    switch(arg.getFunctionSpace().getTypeCode()) {     escriptDataC _temp;
1313       case(Nodes):     escript::Data temp;
1314          /* TODO */     escriptDataC _arg=arg.getDataC();
1315          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");     switch(arg.getFunctionSpace().getTypeCode()) {
1316          break;     case(Nodes):
1317       case(ReducedNodes):     temp=escript::Data( arg, escript::function(*this) );
1318          /* TODO */     _temp=temp.getDataC();
1319          throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1320          break;     break;
1321       case(Elements):     case(ReducedNodes):
1322          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     temp=escript::Data( arg, escript::function(*this) );
1323          break;     _temp=temp.getDataC();
1324       case(ReducedElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1325          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);     break;
1326          break;     case(Elements):
1327       case(FaceElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1328          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     break;
1329          break;     case(ReducedElements):
1330       case(ReducedFaceElements):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1331          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);     break;
1332          break;     case(FaceElements):
1333       case(Points):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1334          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     break;
1335          break;     case(ReducedFaceElements):
1336       case(ContactElementsZero):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1337          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1338          break;     case(Points):
1339       case(ReducedContactElementsZero):     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1340          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1341          break;     case(ContactElementsZero):
1342       case(ContactElementsOne):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1343          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1344          break;     case(ReducedContactElementsZero):
1345       case(ReducedContactElementsOne):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1346          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     break;
1347          break;     case(ContactElementsOne):
1348       case(DegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1349          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");     break;
1350          break;     case(ReducedContactElementsOne):
1351       case(ReducedDegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1352          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");     break;
1353          break;     case(DegreesOfFreedom):
1354       default:     temp=escript::Data( arg, escript::function(*this) );
1355          stringstream temp;     _temp=temp.getDataC();
1356          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1357          throw FinleyAdapterException(temp.str());     break;
1358          break;     case(ReducedDegreesOfFreedom):
1359    }     temp=escript::Data( arg, escript::function(*this) );
1360    checkFinleyError();     _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 771  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    escriptDataC _grad=grad.getDataC();     escriptDataC _grad=grad.getDataC();
1387    escriptDataC nodeDataC;     escriptDataC nodeDataC;
1388  #ifdef PASO_MPI     escript::Data temp;
1389    escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );     if (getMPISize()>1) {
1390    if( arg.getFunctionSpace().getTypeCode() != Nodes )        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1391    {           temp=escript::Data( arg,  continuousFunction(*this) );
1392      Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));           nodeDataC = temp.getDataC();
1393      nodeDataC = nodeTemp.getDataC();        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1394    }           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1395    else           nodeDataC = temp.getDataC();
1396      nodeDataC = arg.getDataC();        } else {
1397  #else           nodeDataC = arg.getDataC();
1398    nodeDataC = arg.getDataC();        }
1399  #endif     } else {
1400    switch(grad.getFunctionSpace().getTypeCode()) {        nodeDataC = arg.getDataC();
1401         case(Nodes):     }
1402            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     switch(grad.getFunctionSpace().getTypeCode()) {
1403            break;     case(Nodes):
1404         case(ReducedNodes):     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1405            throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");     break;
1406            break;     case(ReducedNodes):
1407         case(Elements):     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1408            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     break;
1409            break;     case(Elements):
1410         case(ReducedElements):     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1411            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);     break;
1412            break;     case(ReducedElements):
1413         case(FaceElements):     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1414            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     break;
1415            break;     case(FaceElements):
1416         case(ReducedFaceElements):     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1417            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);     break;
1418            break;     case(ReducedFaceElements):
1419         case(Points):     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1420            throw FinleyAdapterException("Error - Gradient at points is not supported.");     break;
1421            break;     case(Points):
1422         case(ContactElementsZero):     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1423            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1424            break;     case(ContactElementsZero):
1425         case(ReducedContactElementsZero):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1426            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1427            break;     case(ReducedContactElementsZero):
1428         case(ContactElementsOne):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1429            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1430            break;     case(ContactElementsOne):
1431         case(ReducedContactElementsOne):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1432            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);     break;
1433            break;     case(ReducedContactElementsOne):
1434         case(DegreesOfFreedom):     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1435            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     break;
1436            break;     case(DegreesOfFreedom):
1437         case(ReducedDegreesOfFreedom):     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1438            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     break;
1439            break;     case(ReducedDegreesOfFreedom):
1440         default:     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1441            stringstream temp;     break;
1442            temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     default:
1443            throw FinleyAdapterException(temp.str());        stringstream temp;
1444            break;        temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1445    }        throw FinleyAdapterException(temp.str());
1446    checkFinleyError();        break;
1447       }
1448       checkFinleyError();
1449  }  }
1450    
1451  //  //
# Line 847  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(ReducedNodes):     case(ReducedNodes):
1463            throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1464            break;     break;
1465         case(Elements):     case(Elements):
1466            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1467            break;     break;
1468         case(ReducedElements):     case(ReducedElements):
1469            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1470            break;     break;
1471         case(FaceElements):     case(FaceElements):
1472            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1473            break;     break;
1474         case(ReducedFaceElements):     case(ReducedFaceElements):
1475            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1476            break;     break;
1477         case(Points):     case(Points):
1478            throw FinleyAdapterException("Error - Size of point elements is not supported.");     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1479            break;     break;
1480         case(ContactElementsZero):     case(ContactElementsZero):
1481         case(ContactElementsOne):     case(ContactElementsOne):
1482            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1483            break;     break;
1484         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1485         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1486            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1487            break;     break;
1488         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1489            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1490            break;     break;
1491         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1492            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1493            break;     break;
1494         default:     default:
1495            stringstream temp;        stringstream temp;
1496            temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();        temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1497            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1498            break;        break;
1499    }     }
1500    checkFinleyError();     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    escriptDataC tmp;     escriptDataC tmp;
1510    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1511    if (newDomain!=*this)     if (newDomain!=*this)
1512       throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1513    tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1514    Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1515    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1516  }     } else {
1517           throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1518  // saves a data array in openDX format:  /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1519  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         tmp = new_x_inter.getDataC();
1520  {         Finley_Mesh_setCoordinates(mesh,&tmp);*/
1521      int MAX_namelength=256;     }
1522      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     checkFinleyError();
1523    /* win32 refactor */  }
   char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  
   for(int i=0;i<num_data;i++)  
   {  
     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
   }  
1524    
1525    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  //
1526    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  // Helper for the save* methods. Extracts optional data variable names and the
1527    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  // corresponding pointers from python dictionary. Caller must free arrays.
1528    //
1529      boost::python::list keys=arg.keys();  void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1530      for (int i=0;i<num_data;++i) {  {
1531           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);     numData = boost::python::extract<int>(arg.attr("__len__")());
1532           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)     /* win32 refactor */
1533               throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1534           data[i]=d.getDataC();     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1535           ptr_data[i]=&(data[i]);     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1536           std::string n=boost::python::extract<std::string>(keys[i]);  
1537           c_names[i]=&(names[i][0]);     boost::python::list keys=arg.keys();
1538           if (n.length()>MAX_namelength-1) {     for (int i=0; i<numData; ++i) {
1539              strncpy(c_names[i],n.c_str(),MAX_namelength-1);        string n=boost::python::extract<string>(keys[i]);
1540              c_names[i][MAX_namelength-1]='\0';        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1541           } else {        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1542              strcpy(c_names[i],n.c_str());           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1543           }        data[i] = d.getDataC();
1544      }        dataPtr[i] = &(data[i]);
1545      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);        names[i] = TMPMEMALLOC(n.length()+1, char);
1546      checkFinleyError();        strcpy(names[i], n.c_str());
1547           }
1548        /* win32 refactor */  }
1549    TMPMEMFREE(c_names);  
1550    TMPMEMFREE(data);  //
1551    TMPMEMFREE(ptr_data);  // saves mesh and optionally data arrays in openDX format
1552    for(int i=0;i<num_data;i++)  //
1553    {  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1554      TMPMEMFREE(names[i]);  {
1555    }     int num_data;
1556    TMPMEMFREE(names);     char **names;
1557       escriptDataC *data;
1558       escriptDataC **ptr_data;
1559    
1560      return;     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1561       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1562       checkFinleyError();
1563    
1564       /* win32 refactor */
1565       TMPMEMFREE(data);
1566       TMPMEMFREE(ptr_data);
1567       for(int i=0; i<num_data; i++)
1568       {
1569          TMPMEMFREE(names[i]);
1570       }
1571       TMPMEMFREE(names);
1572  }  }
1573    
1574  // saves a data array in openVTK format:  //
1575  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1576  {  //
1577      int MAX_namelength=256;  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1578      const int num_data=boost::python::extract<int>(arg.attr("__len__")());  {
1579    /* win32 refactor */      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1580    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1581    for(int i=0;i<num_data;i++)                metadata, metadata_schema, arg);
1582    {  }
1583      names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
1584    }  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1585    {
1586        if (getMPISize() > 1) {
1587    #ifdef ESYS_MPI
1588            index_t myFirstNode=0, myLastNode=0, k=0;
1589            index_t* globalNodeIndex=0;
1590            Finley_Mesh* mesh_p=m_finleyMesh.get();
1591            /*
1592             * this method is only used by saveDataCSV which would use the returned
1593             * values for reduced nodes wrongly so this case is disabled for now
1594            if (fs_code == FINLEY_REDUCED_NODES)
1595            {
1596                myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1597                myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1598                globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1599            }
1600            else
1601            */
1602            if (fs_code == FINLEY_NODES)
1603            {
1604                myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1605                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1606                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1607            }
1608            else
1609            {
1610                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1611            }
1612    
1613    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;          k=globalNodeIndex[id];
1614    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;          return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
   escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
   
     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());  
          }  
     }  
 #ifndef PASO_MPI      
     Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
 #else  
     Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
1615  #endif  #endif
1616        }
1617        return true;
1618    }
1619    
 checkFinleyError();  
   /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1620    
1621      return;  //
1622    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1623    //
1624    ASM_ptr MeshAdapter::newSystemMatrix(
1625                                                     const int row_blocksize,
1626                                                     const escript::FunctionSpace& row_functionspace,
1627                                                     const int column_blocksize,
1628                                                     const escript::FunctionSpace& column_functionspace,
1629                                                     const int type) const
1630    {
1631       int reduceRowOrder=0;
1632       int reduceColOrder=0;
1633       // is the domain right?
1634       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1635       if (row_domain!=*this)
1636          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1637       const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1638       if (col_domain!=*this)
1639          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1640       // is the function space type right
1641       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1642          reduceRowOrder=0;
1643       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1644          reduceRowOrder=1;
1645       } else {
1646          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1647       }
1648       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1649          reduceColOrder=0;
1650       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1651          reduceColOrder=1;
1652       } else {
1653          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1654       }
1655       // generate matrix:
1656    
1657       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1658       checkFinleyError();
1659       Paso_SystemMatrix* fsystemMatrix;
1660       int trilinos = 0;
1661       if (trilinos) {
1662    #ifdef TRILINOS
1663          /* Allocation Epetra_VrbMatrix here */
1664    #endif
1665       }
1666       else {
1667          fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1668       }
1669       checkPasoError();
1670       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1671       SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1672       return ASM_ptr(sma);
1673    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1674  }  }
1675                                                                                                                                                                            
1676                                                                                                                                                                            //
1677  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a TransportProblemAdapter
1678  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  //
1679                        const int row_blocksize,  ATP_ptr MeshAdapter::newTransportProblem(
1680                        const escript::FunctionSpace& row_functionspace,                                                           const int blocksize,
1681                        const int column_blocksize,                                                           const escript::FunctionSpace& functionspace,
1682                        const escript::FunctionSpace& column_functionspace,                                                           const int type) const
1683                        const int type) const  {
1684  {     int reduceOrder=0;
1685      int reduceRowOrder=0;     // is the domain right?
1686      int reduceColOrder=0;     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1687      // is the domain right?     if (domain!=*this)
1688      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1689      if (row_domain!=*this)     // is the function space type right
1690            throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1691      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());        reduceOrder=0;
1692      if (col_domain!=*this)     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1693            throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        reduceOrder=1;
1694      // is the function space type right     } else {
1695      if (row_functionspace.getTypeCode()==DegreesOfFreedom) {        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1696          reduceRowOrder=0;     }
1697      } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     // generate matrix:
1698          reduceRowOrder=1;  
1699      } else {     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1700          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");     checkFinleyError();
1701      }     Paso_TransportProblem* transportProblem;
1702      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1703          reduceColOrder=0;     checkPasoError();
1704      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1705          reduceColOrder=1;     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1706      } else {     return ATP_ptr(tpa);
1707          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");  //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
     }  
     // generate matrix:  
       
     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);  
     checkFinleyError();  
     Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);  
     checkPasoError();  
     Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);  
     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);  
1708  }  }
1709    
1710  //  //
# Line 1059  SystemMatrixAdapter MeshAdapter::newSyst Line 1716  SystemMatrixAdapter MeshAdapter::newSyst
1716  // 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:
1717  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1718  {  {
1719    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1720         case(Nodes):     case(Nodes):
1721         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1722         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1723            return false;     return false;
1724            break;     break;
1725         case(Elements):     case(Elements):
1726         case(FaceElements):     case(FaceElements):
1727         case(Points):     case(Points):
1728         case(ContactElementsZero):     case(ContactElementsZero):
1729         case(ContactElementsOne):     case(ContactElementsOne):
1730         case(ReducedElements):     case(ReducedElements):
1731         case(ReducedFaceElements):     case(ReducedFaceElements):
1732         case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1733         case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1734            return true;     return true;
1735            break;     break;
1736         default:     default:
1737            stringstream temp;        stringstream temp;
1738            temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1739            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1740            break;        break;
1741    }     }
1742    checkFinleyError();     checkFinleyError();
1743    return false;     return false;
1744    }
1745    
1746    bool
1747    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1748    {
1749       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1750        class 1: DOF <-> Nodes
1751        class 2: ReducedDOF <-> ReducedNodes
1752        class 3: Points
1753        class 4: Elements
1754        class 5: ReducedElements
1755        class 6: FaceElements
1756        class 7: ReducedFaceElements
1757        class 8: ContactElementZero <-> ContactElementOne
1758        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1759    
1760       There is also a set of lines. Interpolation is possible down a line but not between lines.
1761       class 1 and 2 belong to all lines so aren't considered.
1762        line 0: class 3
1763        line 1: class 4,5
1764        line 2: class 6,7
1765        line 3: class 8,9
1766    
1767       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1768       eg hasnodes is true if we have at least one instance of Nodes.
1769       */
1770        if (fs.empty())
1771        {
1772            return false;
1773        }
1774        vector<int> hasclass(10);
1775        vector<int> hasline(4);
1776        bool hasnodes=false;
1777        bool hasrednodes=false;
1778        bool hascez=false;
1779        bool hasrcez=false;
1780        for (int i=0;i<fs.size();++i)
1781        {
1782        switch(fs[i])
1783        {
1784        case(Nodes):   hasnodes=true;   // no break is deliberate
1785        case(DegreesOfFreedom):
1786            hasclass[1]=1;
1787            break;
1788        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1789        case(ReducedDegreesOfFreedom):
1790            hasclass[2]=1;
1791            break;
1792        case(Points):
1793            hasline[0]=1;
1794            hasclass[3]=1;
1795            break;
1796        case(Elements):
1797            hasclass[4]=1;
1798            hasline[1]=1;
1799            break;
1800        case(ReducedElements):
1801            hasclass[5]=1;
1802            hasline[1]=1;
1803            break;
1804        case(FaceElements):
1805            hasclass[6]=1;
1806            hasline[2]=1;
1807            break;
1808        case(ReducedFaceElements):
1809            hasclass[7]=1;
1810            hasline[2]=1;
1811            break;
1812        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1813        case(ContactElementsOne):
1814            hasclass[8]=1;
1815            hasline[3]=1;
1816            break;
1817        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1818        case(ReducedContactElementsOne):
1819            hasclass[9]=1;
1820            hasline[3]=1;
1821            break;
1822        default:
1823            return false;
1824        }
1825        }
1826        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1827        // fail if we have more than one leaf group
1828    
1829        if (totlines>1)
1830        {
1831        return false;   // there are at least two branches we can't interpolate between
1832        }
1833        else if (totlines==1)
1834        {
1835        if (hasline[0]==1)      // we have points
1836        {
1837            resultcode=Points;
1838        }
1839        else if (hasline[1]==1)
1840        {
1841            if (hasclass[5]==1)
1842            {
1843            resultcode=ReducedElements;
1844            }
1845            else
1846            {
1847            resultcode=Elements;
1848            }
1849        }
1850        else if (hasline[2]==1)
1851        {
1852            if (hasclass[7]==1)
1853            {
1854            resultcode=ReducedFaceElements;
1855            }
1856            else
1857            {
1858            resultcode=FaceElements;
1859            }
1860        }
1861        else    // so we must be in line3
1862        {
1863            if (hasclass[9]==1)
1864            {
1865            // need something from class 9
1866            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1867            }
1868            else
1869            {
1870            // something from class 8
1871            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1872            }
1873        }
1874        }
1875        else    // totlines==0
1876        {
1877        if (hasclass[2]==1)
1878        {
1879            // something from class 2
1880            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1881        }
1882        else
1883        {   // something from class 1
1884            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1885        }
1886        }
1887        return true;
1888  }  }
1889    
1890  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1891  {  {
1892    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1893       case(Nodes):     case(Nodes):
1894          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1895             case(Nodes):      case(Nodes):
1896             case(ReducedNodes):      case(ReducedNodes):
1897             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1898             case(DegreesOfFreedom):      case(DegreesOfFreedom):
1899             case(Elements):      case(Elements):
1900             case(ReducedElements):      case(ReducedElements):
1901             case(FaceElements):      case(FaceElements):
1902             case(ReducedFaceElements):      case(ReducedFaceElements):
1903             case(Points):      case(Points):
1904             case(ContactElementsZero):      case(ContactElementsZero):
1905             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1906             case(ContactElementsOne):      case(ContactElementsOne):
1907             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1908                 return true;      return true;
1909             default:      default:
1910                 stringstream temp;            stringstream temp;
1911                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1912                 throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1913          }     }
1914          break;     break;
1915       case(ReducedNodes):     case(ReducedNodes):
1916          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1917             case(ReducedNodes):      case(ReducedNodes):
1918             case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1919             case(Elements):      case(Elements):
1920             case(ReducedElements):      case(ReducedElements):
1921             case(FaceElements):      case(FaceElements):
1922             case(ReducedFaceElements):      case(ReducedFaceElements):
1923             case(Points):      case(Points):
1924             case(ContactElementsZero):      case(ContactElementsZero):
1925             case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1926             case(ContactElementsOne):      case(ContactElementsOne):
1927             case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1928                 return true;      return true;
1929            case(Nodes):      case(Nodes):
1930            case(DegreesOfFreedom):      case(DegreesOfFreedom):
1931               return false;      return false;
1932             default:      default:
1933                 stringstream temp;          stringstream temp;
1934                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1935                 throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1936          }     }
1937          break;     break;
1938       case(Elements):     case(Elements):
1939          if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1940             return true;        return true;
1941          } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1942             return true;        return true;
         } else {  
            return false;  
         }  
      case(ReducedElements):  
         if (functionSpaceType_target==ReducedElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(FaceElements):  
         if (functionSpaceType_target==FaceElements) {  
            return true;  
         } else if (functionSpaceType_target==ReducedFaceElements) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ReducedFaceElements):  
         if (functionSpaceType_target==ReducedFaceElements) {  
            return true;  
1943          } else {          } else {
1944             return false;            return false;
         }  
      case(Points):  
         if (functionSpaceType_target==Points) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ContactElementsZero):  
      case(ContactElementsOne):  
         if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
            return true;  
         } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
            return true;  
         } else {  
            return false;  
         }  
      case(ReducedContactElementsZero):  
      case(ReducedContactElementsOne):  
         if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
            return true;  
         } else {  
            return false;  
         }  
      case(DegreesOfFreedom):  
         switch(functionSpaceType_target) {  
            case(ReducedDegreesOfFreedom):  
            case(DegreesOfFreedom):  
            case(Nodes):  
            case(ReducedNodes):  
            case(Elements):  
            case(ReducedElements):  
            case(Points):  
            case(FaceElements):  
            case(ReducedFaceElements):  
            case(ContactElementsZero):  
            case(ReducedContactElementsZero):  
            case(ContactElementsOne):  
            case(ReducedContactElementsOne):  
               return true;  
            default:  
              stringstream temp;  
              temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
              throw FinleyAdapterException(temp.str());  
1945          }          }
1946          break;     case(ReducedElements):
1947       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1948         switch(functionSpaceType_target) {        return true;
1949            case(ReducedDegreesOfFreedom):      } else {
1950            case(ReducedNodes):            return false;
1951            case(Elements):      }
1952            case(ReducedElements):     case(FaceElements):
1953            case(FaceElements):      if (functionSpaceType_target==FaceElements) {
1954            case(ReducedFaceElements):              return true;
1955            case(Points):      } else if (functionSpaceType_target==ReducedFaceElements) {
1956            case(ContactElementsZero):              return true;
1957            case(ReducedContactElementsZero):      } else {
1958            case(ContactElementsOne):              return false;
1959            case(ReducedContactElementsOne):      }
1960                return true;     case(ReducedFaceElements):
1961            case(Nodes):      if (functionSpaceType_target==ReducedFaceElements) {
1962            case(DegreesOfFreedom):              return true;
1963               return false;      } else {
1964            default:          return false;
1965               stringstream temp;      }
1966               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;     case(Points):
1967               throw FinleyAdapterException(temp.str());      if (functionSpaceType_target==Points) {
1968         }              return true;
1969         break;      } else {
1970       default:              return false;
1971          stringstream temp;      }
1972          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;     case(ContactElementsZero):
1973          throw FinleyAdapterException(temp.str());     case(ContactElementsOne):
1974          break;      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1975    }              return true;
1976    checkFinleyError();      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1977    return false;              return true;
1978        } else {
1979                return false;
1980        }
1981       case(ReducedContactElementsZero):
1982       case(ReducedContactElementsOne):
1983        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1984                return true;
1985        } else {
1986                return false;
1987        }
1988       case(DegreesOfFreedom):
1989        switch(functionSpaceType_target) {
1990        case(ReducedDegreesOfFreedom):
1991        case(DegreesOfFreedom):
1992        case(Nodes):
1993        case(ReducedNodes):
1994        case(Elements):
1995        case(ReducedElements):
1996        case(Points):
1997        case(FaceElements):
1998        case(ReducedFaceElements):
1999        case(ContactElementsZero):
2000        case(ReducedContactElementsZero):
2001        case(ContactElementsOne):
2002        case(ReducedContactElementsOne):
2003        return true;
2004        default:
2005            stringstream temp;
2006            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2007            throw FinleyAdapterException(temp.str());
2008        }
2009        break;
2010       case(ReducedDegreesOfFreedom):
2011       switch(functionSpaceType_target) {
2012        case(ReducedDegreesOfFreedom):
2013        case(ReducedNodes):
2014        case(Elements):
2015        case(ReducedElements):
2016        case(FaceElements):
2017        case(ReducedFaceElements):
2018        case(Points):
2019        case(ContactElementsZero):
2020        case(ReducedContactElementsZero):
2021        case(ContactElementsOne):
2022        case(ReducedContactElementsOne):
2023        return true;
2024        case(Nodes):
2025        case(DegreesOfFreedom):
2026        return false;
2027        default:
2028            stringstream temp;
2029            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
2030            throw FinleyAdapterException(temp.str());
2031        }
2032        break;
2033       default:
2034          stringstream temp;
2035          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2036          throw FinleyAdapterException(temp.str());
2037          break;
2038       }
2039       checkFinleyError();
2040       return false;
2041  }  }
2042    
2043  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
2044  {  {
2045      return false;     return false;
2046  }  }
2047    
2048  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2049  {  {
2050    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2051    if (temp!=0) {     if (temp!=0) {
2052      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2053    } else {     } else {
2054      return false;        return false;
2055    }     }
2056  }  }
2057    
2058  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2059  {  {
2060    return !(operator==(other));     return !(operator==(other));
2061  }  }
2062    
2063  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
2064  {  {
2065     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2066     checkPasoError();     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2067     return out;             package, symmetry, mesh->MPIInfo);
2068    }
2069    
2070    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2071    {
2072       Finley_Mesh* mesh=m_finleyMesh.get();
2073       return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2074               package, symmetry, mesh->MPIInfo);
2075  }  }
2076    
2077  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2078  {  {
2079    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2080  }  }
2081    
2082  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2083  {  {
2084    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2085  }  }
2086    
2087  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2088  {  {
2089    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2090  }  }
2091    
2092  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2093  {  {
2094    int *out=0,i;     int *out = NULL;
2095    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2096    switch (functionSpaceType) {     switch (functionSpaceType) {
2097      case(Nodes):     case(Nodes):
2098        if (mesh->Nodes!=NULL) {     out=mesh->Nodes->Id;
2099          out=mesh->Nodes->Id;     break;
2100          break;     case(ReducedNodes):
2101        }     out=mesh->Nodes->reducedNodesId;
2102      case(ReducedNodes):     break;
2103        throw FinleyAdapterException("Error -  ReducedNodes not supported yet.");     case(Elements):
2104        break;     out=mesh->Elements->Id;
2105      case(Elements):     break;
2106        out=mesh->Elements->Id;     case(ReducedElements):
2107        break;     out=mesh->Elements->Id;
2108      case(ReducedElements):     break;
2109        out=mesh->Elements->Id;     case(FaceElements):
2110        break;     out=mesh->FaceElements->Id;
2111      case(FaceElements):     break;
2112        out=mesh->FaceElements->Id;     case(ReducedFaceElements):
2113        break;     out=mesh->FaceElements->Id;
2114      case(ReducedFaceElements):     break;
2115        out=mesh->FaceElements->Id;     case(Points):
2116        break;     out=mesh->Points->Id;
2117      case(Points):     break;
2118        out=mesh->Points->Id;     case(ContactElementsZero):
2119        break;     out=mesh->ContactElements->Id;
2120      case(ContactElementsZero):     break;
2121        out=mesh->ContactElements->Id;     case(ReducedContactElementsZero):
2122        break;     out=mesh->ContactElements->Id;
2123      case(ReducedContactElementsZero):     break;
2124        out=mesh->ContactElements->Id;     case(ContactElementsOne):
2125        break;     out=mesh->ContactElements->Id;
2126      case(ContactElementsOne):     break;
2127        out=mesh->ContactElements->Id;     case(ReducedContactElementsOne):
2128        break;     out=mesh->ContactElements->Id;
2129      case(ReducedContactElementsOne):     break;
2130        out=mesh->ContactElements->Id;     case(DegreesOfFreedom):
2131        break;     out=mesh->Nodes->degreesOfFreedomId;
2132      case(DegreesOfFreedom):     break;
2133        out=mesh->Nodes->degreeOfFreedomId;     case(ReducedDegreesOfFreedom):
2134        break;     out=mesh->Nodes->reducedDegreesOfFreedomId;
2135      case(ReducedDegreesOfFreedom):     break;
2136        out=mesh->Nodes->reducedDegreeOfFreedomId;     default:
       break;  
     default:  
2137        stringstream temp;        stringstream temp;
2138        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2139        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2140        break;        break;
2141    }     }
2142    return out;     return out;
2143  }  }
2144  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2145  {  {
2146    int out=0;     int out=0;
2147    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2148    switch (functionSpaceType) {     switch (functionSpaceType) {
2149    case(Nodes):     case(Nodes):
2150      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Tag[sampleNo];
2151      break;     break;
2152    case(ReducedNodes):     case(ReducedNodes):
2153      throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2154      break;     break;
2155    case(Elements):     case(Elements):
2156      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2157      break;     break;
2158    case(ReducedElements):     case(ReducedElements):
2159      out=mesh->Elements->Tag[sampleNo];     out=mesh->Elements->Tag[sampleNo];
2160      break;     break;
2161    case(FaceElements):     case(FaceElements):
2162      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2163      break;     break;
2164    case(ReducedFaceElements):     case(ReducedFaceElements):
2165      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->FaceElements->Tag[sampleNo];
2166      break;     break;
2167    case(Points):     case(Points):
2168      out=mesh->Points->Tag[sampleNo];     out=mesh->Points->Tag[sampleNo];
2169      break;     break;
2170    case(ContactElementsZero):     case(ContactElementsZero):
2171      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2172      break;     break;
2173    case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
2174      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2175      break;     break;
2176    case(ContactElementsOne):     case(ContactElementsOne):
2177      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2178      break;     break;
2179    case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
2180      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->ContactElements->Tag[sampleNo];
2181      break;     break;
2182    case(DegreesOfFreedom):     case(DegreesOfFreedom):
2183      throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2184      break;     break;
2185    case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2186      throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2187      break;     break;
2188    default:     default:
2189      stringstream temp;        stringstream temp;
2190      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2191      throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2192      break;        break;
2193    }     }
2194    return out;     return out;
2195  }  }
2196    
2197    
2198  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2199  {  {
2200       Finley_Mesh* mesh=m_finleyMesh.get();
2201       escriptDataC tmp=mask.getDataC();
2202       switch(functionSpaceType) {
2203       case(Nodes):
2204       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2205       break;
2206       case(ReducedNodes):
2207       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2208       break;
2209       case(DegreesOfFreedom):
2210       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2211       break;
2212       case(ReducedDegreesOfFreedom):
2213       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2214       break;
2215       case(Elements):
2216       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2217       break;
2218       case(ReducedElements):
2219       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2220       break;
2221       case(FaceElements):
2222       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2223       break;
2224       case(ReducedFaceElements):
2225       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2226       break;
2227       case(Points):
2228       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2229       break;
2230       case(ContactElementsZero):
2231       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2232       break;
2233       case(ReducedContactElementsZero):
2234       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2235       break;
2236       case(ContactElementsOne):
2237       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2238       break;
2239       case(ReducedContactElementsOne):
2240       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2241       break;
2242       default:
2243          stringstream temp;
2244          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2245          throw FinleyAdapterException(temp.str());
2246       }
2247       checkFinleyError();
2248       return;
2249    }
2250    
2251    void MeshAdapter::setTagMap(const string& name,  int tag)
2252    {
2253       Finley_Mesh* mesh=m_finleyMesh.get();
2254       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2255       checkFinleyError();
2256       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2257    }
2258    
2259    int MeshAdapter::getTag(const string& name) const
2260    {
2261       Finley_Mesh* mesh=m_finleyMesh.get();
2262       int tag=0;
2263       tag=Finley_Mesh_getTag(mesh, name.c_str());
2264       checkFinleyError();
2265       // throwStandardException("MeshAdapter::getTag is not implemented.");
2266       return tag;
2267    }
2268    
2269    bool MeshAdapter::isValidTagName(const string& name) const
2270    {
2271       Finley_Mesh* mesh=m_finleyMesh.get();
2272       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2273    }
2274    
2275    string MeshAdapter::showTagNames() const
2276    {
2277       stringstream temp;
2278       Finley_Mesh* mesh=m_finleyMesh.get();
2279       Finley_TagMap* tag_map=mesh->TagMap;
2280       while (tag_map) {
2281          temp << tag_map->name;
2282          tag_map=tag_map->next;
2283          if (tag_map) temp << ", ";
2284       }
2285       return temp.str();
2286    }
2287    
2288    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2289    {
2290    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2291    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2292    switch(functionSpaceType) {    switch(functionSpaceCode) {
2293         case(Nodes):     case(Nodes):
2294            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2295            break;            break;
2296         case(ReducedNodes):     case(ReducedNodes):
2297            throw FinleyAdapterException("Error - ReducedNodes does not support tags");            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2298            break;            break;
2299         case(DegreesOfFreedom):     case(DegreesOfFreedom):
2300            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2301            break;            break;
2302         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2303            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2304            break;            break;
2305         case(Elements):     case(Elements):
2306            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2307              numTags=mesh->Elements->numTagsInUse;
2308              break;
2309       case(FaceElements):
2310       case(ReducedFaceElements):
2311              numTags=mesh->FaceElements->numTagsInUse;
2312              break;
2313       case(Points):
2314              numTags=mesh->Points->numTagsInUse;
2315              break;
2316       case(ContactElementsZero):
2317       case(ReducedContactElementsZero):
2318       case(ContactElementsOne):
2319       case(ReducedContactElementsOne):
2320              numTags=mesh->ContactElements->numTagsInUse;
2321            break;            break;
2322         case(ReducedElements):     default:
2323            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);        stringstream temp;
2324            break;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2325         case(FaceElements):        throw FinleyAdapterException(temp.str());
2326            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);    }
2327            break;    return numTags;
2328         case(ReducedFaceElements):  }
2329            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);  
2330            break;  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2331         case(Points):  {
2332            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);    Finley_Mesh* mesh=m_finleyMesh.get();
2333      index_t* tags=NULL;
2334      switch(functionSpaceCode) {
2335       case(Nodes):
2336              tags=mesh->Nodes->tagsInUse;
2337            break;            break;
2338         case(ContactElementsZero):     case(ReducedNodes):
2339            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2340            break;            break;
2341         case(ReducedContactElementsZero):     case(DegreesOfFreedom):
2342            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2343            break;            break;
2344         case(ContactElementsOne):     case(ReducedDegreesOfFreedom):
2345            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2346            break;            break;
2347         case(ReducedContactElementsOne):     case(Elements):
2348            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2349              tags=mesh->Elements->tagsInUse;
2350              break;
2351       case(FaceElements):
2352       case(ReducedFaceElements):
2353              tags=mesh->FaceElements->tagsInUse;
2354              break;
2355       case(Points):
2356              tags=mesh->Points->tagsInUse;
2357              break;
2358       case(ContactElementsZero):
2359       case(ReducedContactElementsZero):
2360       case(ContactElementsOne):
2361       case(ReducedContactElementsOne):
2362              tags=mesh->ContactElements->tagsInUse;
2363            break;            break;
2364         default:     default:
2365            stringstream temp;        stringstream temp;
2366            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2367            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2368    }    }
2369    checkFinleyError();    return tags;
   return;  
2370  }  }
2371    
2372  void MeshAdapter::setTagMap(const std::string& name,  int tag)  
2373    bool MeshAdapter::canTag(int functionSpaceCode) const
2374  {  {
2375    Finley_Mesh* mesh=m_finleyMesh.get();    switch(functionSpaceCode) {
2376    Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     case(Nodes):
2377    checkPasoError();     case(Elements):
2378    // throwStandardException("MeshAdapter::set TagMap is not implemented.");     case(ReducedElements):
2379       case(FaceElements):
2380       case(ReducedFaceElements):
2381       case(Points):
2382       case(ContactElementsZero):
2383       case(ReducedContactElementsZero):
2384       case(ContactElementsOne):
2385       case(ReducedContactElementsOne):
2386              return true;
2387       case(ReducedNodes):
2388       case(DegreesOfFreedom):
2389       case(ReducedDegreesOfFreedom):
2390          return false;
2391       default:
2392        return false;
2393      }
2394  }  }
2395    
2396  int MeshAdapter::getTag(const std::string& name) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2397  {  {
2398    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2399    int tag=0;    return Finley_Mesh_getStatus(mesh);
   tag=Finley_Mesh_getTag(mesh, name.c_str());  
   checkPasoError();  
   // throwStandardException("MeshAdapter::getTag is not implemented.");  
   return tag;  
2400  }  }
2401    
2402  bool MeshAdapter::isValidTagName(const std::string& name) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2403  {  {
2404      
2405    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2406    return Finley_Mesh_isValidTagName(mesh,name.c_str());    int order =-1;
2407      switch(functionSpaceCode) {
2408       case(Nodes):
2409       case(DegreesOfFreedom):
2410              order=mesh->approximationOrder;
2411              break;
2412       case(ReducedNodes):
2413       case(ReducedDegreesOfFreedom):
2414              order=mesh->reducedApproximationOrder;
2415              break;
2416       case(Elements):
2417       case(FaceElements):
2418       case(Points):
2419       case(ContactElementsZero):
2420       case(ContactElementsOne):
2421              order=mesh->integrationOrder;
2422              break;
2423       case(ReducedElements):
2424       case(ReducedFaceElements):
2425       case(ReducedContactElementsZero):
2426       case(ReducedContactElementsOne):
2427              order=mesh->reducedIntegrationOrder;
2428              break;
2429       default:
2430          stringstream temp;
2431          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2432          throw FinleyAdapterException(temp.str());
2433      }
2434      return order;
2435  }  }
2436    
2437  std::string MeshAdapter::showTagNames() const  bool MeshAdapter::supportsContactElements() const
2438  {  {
2439    stringstream temp;    return true;
2440    Finley_Mesh* mesh=m_finleyMesh.get();  }
2441    Finley_TagMap* tag_map=mesh->TagMap;  
2442    while (tag_map) {  ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2443       temp << tag_map->name;  {
2444       tag_map=tag_map->next;    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2445       if (tag_map) temp << ", ";  }
2446    }  
2447    return temp.str();  ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2448    {
2449      Finley_ReferenceElementSet_dealloc(m_refSet);
2450    }
2451    
2452    // points will be flattened
2453    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2454    {
2455          const int dim = getDim();
2456          int numPoints=points.size()/dim;
2457          int numTags=tags.size();
2458          Finley_Mesh* mesh=m_finleyMesh.get();
2459          
2460          if ( points.size() % dim != 0 )
2461          {
2462        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2463          }
2464          
2465          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2466         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2467          
2468          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2469          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2470          
2471          for (int i=0;i<numPoints;++i) {
2472           points_ptr[ i * dim     ] = points[i * dim ];
2473           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2474           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2475               tags_ptr[i]=tags[i];
2476          }
2477          
2478          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2479          checkFinleyError();
2480          
2481          TMPMEMFREE(points_ptr);
2482          TMPMEMFREE(tags_ptr);
2483    }
2484    
2485    
2486    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2487    // {
2488    //       const int dim = getDim();
2489    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2490    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2491    //       Finley_Mesh* mesh=m_finleyMesh.get();
2492    //      
2493    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2494    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2495    //      
2496    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2497    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2498    //      
2499    //       for (int i=0;i<numPoints;++i) {
2500    //     int tag_id=-1;
2501    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2502    //     if  ( numComps !=   dim ) {
2503    //                stringstream temp;            
2504    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2505    //                throw FinleyAdapterException(temp.str());
2506    //     }
2507    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2508    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2509    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2510    //    
2511    //     if ( numTags > 0) {
2512    //            boost::python::extract<string> ex_str(tags[i]);
2513    //        if  ( ex_str.check() ) {
2514    //            tag_id=getTag( ex_str());
2515    //        } else {
2516    //             boost::python::extract<int> ex_int(tags[i]);
2517    //             if ( ex_int.check() ) {
2518    //                 tag_id=ex_int();
2519    //             } else {
2520    //              stringstream temp;          
2521    //                  temp << "Error - unable to extract tag for point " << i;
2522    //              throw FinleyAdapterException(temp.str());
2523    //            }
2524    //        }
2525    //     }      
2526    //            tags_ptr[i]=tag_id;
2527    //       }
2528    //      
2529    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2530    //       checkPasoError();
2531    //      
2532    //       TMPMEMFREE(points_ptr);
2533    //       TMPMEMFREE(tags_ptr);
2534    // }
2535    
2536    /*
2537    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2538    {  
2539        boost::python::list points =  boost::python::list();
2540        boost::python::list tags =  boost::python::list();
2541        points.append(point);
2542        tags.append(tag);
2543        addDiracPoints(points, tags);
2544  }  }
2545    */
2546    
2547    /*
2548    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2549    {
2550            boost::python::list points =   boost::python::list();
2551            boost::python::list tags =   boost::python::list();
2552            points.append(point);
2553            tags.append(tag);
2554            addDiracPoints(points, tags);
2555    }
2556    */
2557  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26