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

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

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

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

Legend:
Removed from v.757  
changed lines
  Added in v.3911

  ViewVC Help
Powered by ViewVC 1.1.26