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

trunk/finley/src/finley/CPPAdapter/MeshAdapter.cpp revision 201 by jgs, Wed Nov 23 04:10:21 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 1  Line 1 
 // $Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
 extern "C" {  
 #include "finley/finleyC/Finley.h"  
 #include "finley/finleyC/Assemble.h"  
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/Finley.h"  
 #include "paso/SystemMatrix.h"  
 }  
 #include "finley/CPPAdapter/SystemMatrixAdapter.h"  
 #include "finley/CPPAdapter/MeshAdapter.h"  
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/FinleyAdapterException.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/Data.h"  
 #include "escript/Data/DataArrayView.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/DataFactory.h"  
 #include <boost/python/extract.hpp>  
1    
2  #include <iostream>  /*******************************************************
3  #include <sstream>  *
4    * Copyright (c) 2003-2010 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    
15    #include "MeshAdapter.h"
16    #include "escript/Data.h"
17    #include "escript/DataFactory.h"
18    #ifdef USE_NETCDF
19    #include <netcdfcpp.h>
20    #endif
21    #ifdef ESYS_MPI
22    #include <mpi.h>
23    #include "esysUtils/Esys_MPI.h"
24    #endif
25    extern "C" {
26    #include "esysUtils/blocktimer.h"
27    }
28    
29  using namespace std;  using namespace std;
30  using namespace escript;  using namespace escript;
# Line 44  MeshAdapter::FunctionSpaceNamesMapType M Line 37  MeshAdapter::FunctionSpaceNamesMapType M
37  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
38  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
39  const int MeshAdapter::Nodes=FINLEY_NODES;  const int MeshAdapter::Nodes=FINLEY_NODES;
40    const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
41  const int MeshAdapter::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=FINLEY_ELEMENTS;
42    const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
43  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
44    const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
45  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=FINLEY_POINTS;
46  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
47    const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
48  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
49    const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
50    
51  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
52  {  {
53    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
54    //     //
55    // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer     // need to use a null_deleter as Finley_Mesh_free deletes the pointer
56    // for us.     // for us.
57    m_finleyMesh.reset(finleyMesh,null_deleter());     m_finleyMesh.reset(finleyMesh,null_deleter());
58  }  }
59    
60  //  //
# Line 64  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 62  MeshAdapter::MeshAdapter(Finley_Mesh* fi
62  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
63  m_finleyMesh(in.m_finleyMesh)  m_finleyMesh(in.m_finleyMesh)
64  {  {
65    setFunctionSpaceTypeNames();     setFunctionSpaceTypeNames();
66  }  }
67    
68  MeshAdapter::~MeshAdapter()  MeshAdapter::~MeshAdapter()
69  {  {
70    //     //
71    // I hope the case for the pointer being zero has been taken care of.     // I hope the case for the pointer being zero has been taken care of.
72    //  cout << "In MeshAdapter destructor." << endl;     //  cout << "In MeshAdapter destructor." << endl;
73    if (m_finleyMesh.unique()) {     if (m_finleyMesh.unique()) {
74      //   cout << "Calling dealloc." << endl;        Finley_Mesh_free(m_finleyMesh.get());
75      Finley_Mesh_dealloc(m_finleyMesh.get());     }
     //   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];     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
   strcpy(fName,fileName.c_str());  
   Finley_Mesh_write(m_finleyMesh.get(),fName);  
   checkFinleyError();  
130  }  }
131    
132  // void MeshAdapter::getTagList(int functionSpaceType,  void MeshAdapter::dump(const string& fileName) const
133  //                  int* numTags) const  {
134  // {  #ifdef USE_NETCDF
135  //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
136  //   return;     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::getDiracDeltaFunctionCode() const
653  {  {
654    return Points;     return Points;
 }  
   
 //  
 // returns a pointer to the tag list of samples of functionSpaceType  
 //  
 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,  
                  int* numTags) const  
 {  
   *tagList=NULL;  
   *numTags=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *tagList=mesh->Nodes->Tag;  
       *numTags=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *tagList=mesh->Elements->Tag;  
       *numTags=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *tagList=mesh->FaceElements->Tag;  
       *numTags=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *tagList=mesh->Points->Tag;  
       *numTags=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*tagList==NULL) {  
     stringstream temp;  
     temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
655  }  }
656    
657  //  //
658  // returns a pointer to the reference no list of samples of functionSpaceType  // return the spatial dimension of the Mesh:
659  //  //
660  void MeshAdapter::getReferenceNoList(int functionSpaceType, int** referenceNoList,  int MeshAdapter::getDim() const
                  int* numReferenceNo) const  
661  {  {
662    *referenceNoList=NULL;     Finley_Mesh* mesh=m_finleyMesh.get();
663    *numReferenceNo=0;     int numDim=Finley_Mesh_getDim(mesh);
664    Finley_Mesh* mesh=m_finleyMesh.get();     checkFinleyError();
665    switch (functionSpaceType) {     return numDim;
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=mesh->Nodes->Id;  
       *numReferenceNo=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *referenceNoList=mesh->Elements->Id;  
       *numReferenceNo=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *referenceNoList=mesh->FaceElements->Id;  
       *numReferenceNo=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *referenceNoList=mesh->Points->Id;  
       *numReferenceNo=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*referenceNoList==NULL) {  
     stringstream temp;  
     temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
666  }  }
667    
668  //  //
669  // return the spatial dimension of the Mesh:  // Return the number of data points summed across all MPI processes
670  //  //
671  int MeshAdapter::getDim() const  int MeshAdapter::getNumDataPointsGlobal() const
672  {  {
673    int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
   checkFinleyError();  
   return numDim;  
674  }  }
675    
676  //  //
# Line 345  pair<int,int> MeshAdapter::getDataShape( Line 683  pair<int,int> MeshAdapter::getDataShape(
683     int numSamples=0;     int numSamples=0;
684     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
685     switch (functionSpaceCode) {     switch (functionSpaceCode) {
686        case(Nodes):     case(Nodes):
687             numDataPointsPerSample=1;     numDataPointsPerSample=1;
688             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
689             break;     break;
690        case(Elements):     case(ReducedNodes):
691             if (mesh->Elements!=NULL) {     numDataPointsPerSample=1;
692               numSamples=mesh->Elements->numElements;     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
693               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;     break;
694             }     case(Elements):
695             break;     if (mesh->Elements!=NULL) {
696        case(FaceElements):        numSamples=mesh->Elements->numElements;
697             if (mesh->FaceElements!=NULL) {        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
698                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;     }
699                  numSamples=mesh->FaceElements->numElements;     break;
700             }     case(ReducedElements):
701             break;     if (mesh->Elements!=NULL) {
702        case(Points):        numSamples=mesh->Elements->numElements;
703             if (mesh->Points!=NULL) {        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
704               numDataPointsPerSample=1;     }
705               numSamples=mesh->Points->numElements;     break;
706             }     case(FaceElements):
707             break;     if (mesh->FaceElements!=NULL) {
708        case(ContactElementsZero):        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
709             if (mesh->ContactElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
710               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
711               numSamples=mesh->ContactElements->numElements;     break;
712             }     case(ReducedFaceElements):
713             break;     if (mesh->FaceElements!=NULL) {
714        case(ContactElementsOne):        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
715             if (mesh->ContactElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
716               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
717               numSamples=mesh->ContactElements->numElements;     break;
718             }     case(Points):
719             break;     if (mesh->Points!=NULL) {
720        case(DegreesOfFreedom):        numDataPointsPerSample=1;
721             if (mesh->Nodes!=NULL) {        numSamples=mesh->Points->numElements;
722               numDataPointsPerSample=1;     }
723               numSamples=mesh->Nodes->numDegreesOfFreedom;     break;
724             }     case(ContactElementsZero):
725             break;     if (mesh->ContactElements!=NULL) {
726        case(ReducedDegreesOfFreedom):        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
727             if (mesh->Nodes!=NULL) {        numSamples=mesh->ContactElements->numElements;
728               numDataPointsPerSample=1;     }
729               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;     break;
730             }     case(ReducedContactElementsZero):
731             break;     if (mesh->ContactElements!=NULL) {
732        default:        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
733          stringstream temp;        numSamples=mesh->ContactElements->numElements;
734          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();     }
735          throw FinleyAdapterException(temp.str());     break;
736          break;     case(ContactElementsOne):
737        }     if (mesh->ContactElements!=NULL) {
738        return pair<int,int>(numDataPointsPerSample,numSamples);        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
739          numSamples=mesh->ContactElements->numElements;
740       }
741       break;
742       case(ReducedContactElementsOne):
743       if (mesh->ContactElements!=NULL) {
744          numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
745          numSamples=mesh->ContactElements->numElements;
746       }
747       break;
748       case(DegreesOfFreedom):
749       if (mesh->Nodes!=NULL) {
750          numDataPointsPerSample=1;
751          numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
752       }
753       break;
754       case(ReducedDegreesOfFreedom):
755       if (mesh->Nodes!=NULL) {
756          numDataPointsPerSample=1;
757          numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
758       }
759       break;
760       default:
761          stringstream temp;
762          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
763          throw FinleyAdapterException(temp.str());
764          break;
765       }
766       return pair<int,int>(numDataPointsPerSample,numSamples);
767  }  }
768    
769  //  //
770  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
771  //  //
772  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
773                       SystemMatrixAdapter& mat, Data& rhs,                                   SystemMatrixAdapter& mat, escript::Data& rhs,
774                       const Data& A, const Data& B, const Data& C,const  Data& D,const  Data& X,const  Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
775                       const Data& d, const Data& y,                                   const escript::Data& d, const escript::Data& y,
776                       const Data& d_contact,const Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
777  {  {
778     Finley_Mesh* mesh=m_finleyMesh.get();     escriptDataC _rhs=rhs.getDataC();
779     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),     escriptDataC _A  =A.getDataC();
780                         &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _B=B.getDataC();
781     checkFinleyError();     escriptDataC _C=C.getDataC();
782     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,     escriptDataC _D=D.getDataC();
783                    mat.getPaso_SystemMatrix(),     escriptDataC _X=X.getDataC();
784                    &(rhs.getDataC()),     escriptDataC _Y=Y.getDataC();
785                                    &(d.getDataC()),&(y.getDataC()),     escriptDataC _d=d.getDataC();
786                                    Finley_Assemble_handelShapeMissMatch_Mean_out);     escriptDataC _y=y.getDataC();
787     checkFinleyError();     escriptDataC _d_contact=d_contact.getDataC();
788     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,     escriptDataC _y_contact=y_contact.getDataC();
789                    mat.getPaso_SystemMatrix(),  
790                    &(rhs.getDataC()),     Finley_Mesh* mesh=m_finleyMesh.get();
791                                    &(d_contact.getDataC()),  
792                    &(y_contact.getDataC()),     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
793                                    Finley_Assemble_handelShapeMissMatch_Step_out);     checkFinleyError();
794    
795       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
796       checkFinleyError();
797    
798       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
799     checkFinleyError();     checkFinleyError();
800  }  }
801    
802    void  MeshAdapter::addPDEToLumpedSystem(
803                                            escript::Data& mat,
804                                            const escript::Data& D,
805                                            const escript::Data& d) const
806    {
807       escriptDataC _mat=mat.getDataC();
808       escriptDataC _D=D.getDataC();
809       escriptDataC _d=d.getDataC();
810    
811       Finley_Mesh* mesh=m_finleyMesh.get();
812    
813       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
814       checkFinleyError();
815      
816       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
817       checkFinleyError();
818    
819    }
820    
821    
822  //  //
823  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
824  //  //
825  void MeshAdapter::addPDEToRHS( Data& rhs,  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
                      const  Data& X,const  Data& Y, const Data& y, const Data& y_contact) const  
826  {  {
827     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
828    
829     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
830     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     escriptDataC _X=X.getDataC();
831       escriptDataC _Y=Y.getDataC();
832       escriptDataC _y=y.getDataC();
833       escriptDataC _y_contact=y_contact.getDataC();
834    
835       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
836       checkFinleyError();
837    
838       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
839       checkFinleyError();
840    
841       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
842       checkFinleyError();
843    }
844    //
845    // adds PDE of second order into a transport problem
846    //
847    void MeshAdapter::addPDEToTransportProblem(
848                                               TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
849                                               const escript::Data& A, const escript::Data& B, const escript::Data& C,
850                                               const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
851                                               const escript::Data& d, const escript::Data& y,
852                                               const escript::Data& d_contact,const escript::Data& y_contact) const
853    {
854       DataTypes::ShapeType shape;
855       source.expand();
856       escriptDataC _source=source.getDataC();
857       escriptDataC _M=M.getDataC();
858       escriptDataC _A=A.getDataC();
859       escriptDataC _B=B.getDataC();
860       escriptDataC _C=C.getDataC();
861       escriptDataC _D=D.getDataC();
862       escriptDataC _X=X.getDataC();
863       escriptDataC _Y=Y.getDataC();
864       escriptDataC _d=d.getDataC();
865       escriptDataC _y=y.getDataC();
866       escriptDataC _d_contact=d_contact.getDataC();
867       escriptDataC _y_contact=y_contact.getDataC();
868    
869       Finley_Mesh* mesh=m_finleyMesh.get();
870       Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();
871    
872       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
873     checkFinleyError();     checkFinleyError();
874    
875     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
876     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     checkFinleyError();
877    
878       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
879     checkFinleyError();     checkFinleyError();
880     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
881     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
882     checkFinleyError();     checkFinleyError();
883  }  }
884    
885  //  //
886  // interpolates data between different function spaces:  // interpolates data between different function spaces:
887  //  //
888  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
889  {  {
890    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
891    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
892    if (inDomain!=*this)     if (inDomain!=*this)  
893       throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
894    if (targetDomain!=*this)     if (targetDomain!=*this)
895       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
896    
897    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
898    switch(in.getFunctionSpace().getTypeCode()) {     escriptDataC _target=target.getDataC();
899       case(Nodes):     escriptDataC _in=in.getDataC();
900          switch(target.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
901             case(Nodes):     case(Nodes):
902             case(ReducedDegreesOfFreedom):        switch(target.getFunctionSpace().getTypeCode()) {
903             case(DegreesOfFreedom):        case(Nodes):
904                 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));        case(ReducedNodes):
905                 break;        case(DegreesOfFreedom):
906             case(Elements):        case(ReducedDegreesOfFreedom):
907                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
908                 break;        break;
909             case(FaceElements):        case(Elements):
910                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));        case(ReducedElements):
911                 break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
912             case(Points):        break;
913                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));        case(FaceElements):
914                 break;        case(ReducedFaceElements):
915             case(ContactElementsZero):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
916                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        break;
917                 break;        case(Points):
918             case(ContactElementsOne):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
919                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        break;
920                 break;        case(ContactElementsZero):
921             default:        case(ReducedContactElementsZero):
922                 stringstream temp;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
923                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        break;
924                 throw FinleyAdapterException(temp.str());        case(ContactElementsOne):
925                 break;        case(ReducedContactElementsOne):
926          }        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
927          break;        break;
928       case(Elements):        default:
929          if (target.getFunctionSpace().getTypeCode()==Elements) {           stringstream temp;
930             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();
931          } else {           throw FinleyAdapterException(temp.str());
932             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           break;
933          }        }
934          break;        break;
935       case(FaceElements):     case(ReducedNodes):
936          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        switch(target.getFunctionSpace().getTypeCode()) {
937             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));        case(Nodes):
938          } else {        case(ReducedNodes):
939             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");        case(DegreesOfFreedom):
940             break;        case(ReducedDegreesOfFreedom):
941         }        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
942       case(Points):        break;
943          if (target.getFunctionSpace().getTypeCode()==Points) {        case(Elements):
944             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));        case(ReducedElements):
945          } else {        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
946             throw FinleyAdapterException("Error - No interpolation with data on points possible.");        break;
947             break;        case(FaceElements):
948          }        case(ReducedFaceElements):
949          break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
950       case(ContactElementsZero):        break;
951       case(ContactElementsOne):        case(Points):
952          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
953             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));        break;
954          } else {        case(ContactElementsZero):
955             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");        case(ReducedContactElementsZero):
956             break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
957          }        break;
958          break;        case(ContactElementsOne):
959       case(DegreesOfFreedom):        case(ReducedContactElementsOne):
960          switch(target.getFunctionSpace().getTypeCode()) {        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
961             case(ReducedDegreesOfFreedom):        break;
962             case(DegreesOfFreedom):        default:
963             case(Nodes):           stringstream temp;
964                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();
965                break;           throw FinleyAdapterException(temp.str());
966             case(Elements):           break;
967                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));        }
968                break;        break;
969             case(FaceElements):     case(Elements):
970                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));        if (target.getFunctionSpace().getTypeCode()==Elements) {
971                break;           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
972             case(Points):        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
973                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
974                break;        } else {
975             case(ContactElementsZero):           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
976             case(ContactElementsOne):        }
977                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        break;
978               break;     case(ReducedElements):
979             default:        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
980               stringstream temp;           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
981               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        } else {
982               throw FinleyAdapterException(temp.str());           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
983               break;        }
984          }        break;
985          break;     case(FaceElements):
986       case(ReducedDegreesOfFreedom):        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
987         switch(target.getFunctionSpace().getTypeCode()) {           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
988            case(ReducedDegreesOfFreedom):        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
989               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
990               break;        } else {
991            case(Elements):           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
992               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));        }
993               break;        break;
994            case(FaceElements):     case(ReducedFaceElements):
995               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
996               break;           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
997            case(Points):        } else {
998               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
999               break;        }
1000            case(ContactElementsZero):        break;
1001            case(ContactElementsOne):     case(Points):
1002               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        if (target.getFunctionSpace().getTypeCode()==Points) {
1003               break;           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1004            case(Nodes):        } else {
1005               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1006               break;        }
1007            case(DegreesOfFreedom):        break;
1008               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");     case(ContactElementsZero):
1009               break;     case(ContactElementsOne):
1010            default:        if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1011               stringstream temp;           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1012               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1013               throw FinleyAdapterException(temp.str());           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1014               break;        } else {
1015         }           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1016         break;        }
1017       default:        break;
1018          stringstream temp;     case(ReducedContactElementsZero):
1019          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();     case(ReducedContactElementsOne):
1020          throw FinleyAdapterException(temp.str());        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1021          break;           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1022    }        } else {
1023    checkFinleyError();           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1024          }
1025          break;
1026       case(DegreesOfFreedom):      
1027          switch(target.getFunctionSpace().getTypeCode()) {
1028          case(ReducedDegreesOfFreedom):
1029          case(DegreesOfFreedom):
1030          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1031          break;
1032      
1033          case(Nodes):
1034          case(ReducedNodes):
1035          if (getMPISize()>1) {
1036             escript::Data temp=escript::Data(in);
1037             temp.expand();
1038             escriptDataC _in2 = temp.getDataC();
1039             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1040          } else {
1041             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1042          }
1043          break;
1044          case(Elements):
1045          case(ReducedElements):
1046          if (getMPISize()>1) {
1047             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1048             escriptDataC _in2 = temp.getDataC();
1049             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1050          } else {
1051             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1052          }
1053          break;
1054          case(FaceElements):
1055          case(ReducedFaceElements):
1056          if (getMPISize()>1) {
1057             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1058             escriptDataC _in2 = temp.getDataC();
1059             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1060      
1061          } else {
1062             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1063          }
1064          break;
1065          case(Points):
1066          if (getMPISize()>1) {
1067             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1068             escriptDataC _in2 = temp.getDataC();
1069          } else {
1070             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1071          }
1072          break;
1073          case(ContactElementsZero):
1074          case(ContactElementsOne):
1075          case(ReducedContactElementsZero):
1076          case(ReducedContactElementsOne):
1077          if (getMPISize()>1) {
1078             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1079             escriptDataC _in2 = temp.getDataC();
1080             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1081          } else {
1082             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1083          }
1084          break;
1085          default:
1086             stringstream temp;
1087             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1088             throw FinleyAdapterException(temp.str());
1089             break;
1090          }
1091          break;
1092       case(ReducedDegreesOfFreedom):
1093          switch(target.getFunctionSpace().getTypeCode()) {
1094          case(Nodes):
1095          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1096          break;
1097          case(ReducedNodes):
1098          if (getMPISize()>1) {
1099             escript::Data temp=escript::Data(in);
1100             temp.expand();
1101             escriptDataC _in2 = temp.getDataC();
1102             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1103          } else {
1104             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1105          }
1106          break;
1107          case(DegreesOfFreedom):
1108          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1109          break;
1110          case(ReducedDegreesOfFreedom):
1111          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1112          break;
1113          case(Elements):
1114          case(ReducedElements):
1115          if (getMPISize()>1) {
1116             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1117             escriptDataC _in2 = temp.getDataC();
1118             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1119          } else {
1120             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1121          }
1122          break;
1123          case(FaceElements):
1124          case(ReducedFaceElements):
1125          if (getMPISize()>1) {
1126             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1127             escriptDataC _in2 = temp.getDataC();
1128             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1129          } else {
1130             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1131          }
1132          break;
1133          case(Points):
1134          if (getMPISize()>1) {
1135             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1136             escriptDataC _in2 = temp.getDataC();
1137             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1138          } else {
1139             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1140          }
1141          break;
1142          case(ContactElementsZero):
1143          case(ContactElementsOne):
1144          case(ReducedContactElementsZero):
1145          case(ReducedContactElementsOne):
1146          if (getMPISize()>1) {
1147             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1148             escriptDataC _in2 = temp.getDataC();
1149             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1150          } else {
1151             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1152          }
1153          break;
1154          default:
1155             stringstream temp;
1156             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1157             throw FinleyAdapterException(temp.str());
1158             break;
1159          }
1160          break;
1161       default:
1162          stringstream temp;
1163          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1164          throw FinleyAdapterException(temp.str());
1165          break;
1166       }
1167       checkFinleyError();
1168  }  }
1169    
1170  //  //
1171  // copies the locations of sample points into x:  // copies the locations of sample points into x:
1172  //  //
1173  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1174  {  {
1175    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1176    if (argDomain!=*this)     if (argDomain!=*this)
1177       throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1178    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1179    // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1180    if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1181       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));        escriptDataC _arg=arg.getDataC();
1182    } else {        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1183       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);     } else {
1184       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1185       // this is then interpolated onto arg:        escriptDataC _tmp_data=tmp_data.getDataC();
1186       interpolateOnDomain(arg,tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1187    }        // this is then interpolated onto arg:
1188    checkFinleyError();        interpolateOnDomain(arg,tmp_data);
1189       }
1190       checkFinleyError();
1191  }  }
1192    
1193  //  //
1194  // return the normal vectors at the location of data points as a Data object:  // return the normal vectors at the location of data points as a Data object:
1195  //  //
1196  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1197  {  {
1198    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1199    if (normalDomain!=*this)     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1200       throw FinleyAdapterException("Error - Illegal domain of normal locations");     if (normalDomain!=*this)
1201    Finley_Mesh* mesh=m_finleyMesh.get();        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1202    switch(normal.getFunctionSpace().getTypeCode()) {     Finley_Mesh* mesh=m_finleyMesh.get();
1203      case(Nodes):     escriptDataC _normal=normal.getDataC();
1204        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     switch(normal.getFunctionSpace().getTypeCode()) {
1205        break;     case(Nodes):
1206      case(Elements):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1207        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     break;
1208        break;     case(ReducedNodes):
1209      case (FaceElements):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1210        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));     break;
1211        break;     case(Elements):
1212      case(Points):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1213        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     break;
1214        break;     case(ReducedElements):
1215      case (ContactElementsOne):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1216      case (ContactElementsZero):     break;
1217        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));     case (FaceElements):
1218        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1219      case(DegreesOfFreedom):     break;
1220        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case (ReducedFaceElements):
1221        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1222      case(ReducedDegreesOfFreedom):     break;
1223        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(Points):
1224        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1225      default:     break;
1226       case (ContactElementsOne):
1227       case (ContactElementsZero):
1228       Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1229       break;
1230       case (ReducedContactElementsOne):
1231       case (ReducedContactElementsZero):
1232       Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1233       break;
1234       case(DegreesOfFreedom):
1235       throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1236       break;
1237       case(ReducedDegreesOfFreedom):
1238       throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1239       break;
1240       default:
1241        stringstream temp;        stringstream temp;
1242        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();
1243        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1244        break;        break;
1245    }     }
1246    checkFinleyError();     checkFinleyError();
1247  }  }
1248    
1249  //  //
1250  // interpolates data to other domain:  // interpolates data to other domain:
1251  //  //
1252  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1253  {  {
1254    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1255    if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1256       throw FinleyAdapterException("Error - Illegal domain of interpolation target");     if (targetDomain!=this)
1257          throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1258    
1259    throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1260  }  }
1261    
1262  //  //
1263  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1264  //  //
1265  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1266  {  {
1267    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1268    if (argDomain!=*this)     if (argDomain!=*this)
1269       throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1270    
1271    Finley_Mesh* mesh=m_finleyMesh.get();     double blocktimer_start = blocktimer_time();
1272    switch(arg.getFunctionSpace().getTypeCode()) {     Finley_Mesh* mesh=m_finleyMesh.get();
1273       case(Nodes):     escriptDataC _temp;
1274          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");     escript::Data temp;
1275          break;     escriptDataC _arg=arg.getDataC();
1276       case(Elements):     switch(arg.getFunctionSpace().getTypeCode()) {
1277          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);     case(Nodes):
1278          break;     temp=escript::Data( arg, escript::function(*this) );
1279       case(FaceElements):     _temp=temp.getDataC();
1280          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1281          break;     break;
1282       case(Points):     case(ReducedNodes):
1283          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     temp=escript::Data( arg, escript::function(*this) );
1284          break;     _temp=temp.getDataC();
1285       case(ContactElementsZero):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1286          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);     break;
1287          break;     case(Elements):
1288       case(ContactElementsOne):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1289          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);     break;
1290          break;     case(ReducedElements):
1291       case(DegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1292          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");     break;
1293          break;     case(FaceElements):
1294       case(ReducedDegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1295          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");     break;
1296          break;     case(ReducedFaceElements):
1297       default:     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1298          stringstream temp;     break;
1299          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     case(Points):
1300          throw FinleyAdapterException(temp.str());     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1301          break;     break;
1302    }     case(ContactElementsZero):
1303    checkFinleyError();     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1304       break;
1305       case(ReducedContactElementsZero):
1306       Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1307       break;
1308       case(ContactElementsOne):
1309       Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1310       break;
1311       case(ReducedContactElementsOne):
1312       Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1313       break;
1314       case(DegreesOfFreedom):
1315       temp=escript::Data( arg, escript::function(*this) );
1316       _temp=temp.getDataC();
1317       Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1318       break;
1319       case(ReducedDegreesOfFreedom):
1320       temp=escript::Data( arg, escript::function(*this) );
1321       _temp=temp.getDataC();
1322       Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1323       break;
1324       default:
1325          stringstream temp;
1326          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1327          throw FinleyAdapterException(temp.str());
1328          break;
1329       }
1330       checkFinleyError();
1331       blocktimer_increment("integrate()", blocktimer_start);
1332  }  }
1333    
1334  //  //
1335  // calculates the gradient of arg:  // calculates the gradient of arg:
1336  //  //
1337  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1338  {  {
1339    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1340    if (argDomain!=*this)     if (argDomain!=*this)
1341      throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1342    const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1343    if (gradDomain!=*this)     if (gradDomain!=*this)
1344       throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1345    
1346    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1347    switch(grad.getFunctionSpace().getTypeCode()) {     escriptDataC _grad=grad.getDataC();
1348         case(Nodes):     escriptDataC nodeDataC;
1349            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     escript::Data temp;
1350            break;     if (getMPISize()>1) {
1351         case(Elements):        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1352            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));           temp=escript::Data( arg,  continuousFunction(*this) );
1353            break;           nodeDataC = temp.getDataC();
1354         case(FaceElements):        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1355            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1356            break;           nodeDataC = temp.getDataC();
1357         case(Points):        } else {
1358            throw FinleyAdapterException("Error - Gradient at points is not supported.");           nodeDataC = arg.getDataC();
1359            break;        }
1360         case(ContactElementsZero):     } else {
1361            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));        nodeDataC = arg.getDataC();
1362            break;     }
1363         case(ContactElementsOne):     switch(grad.getFunctionSpace().getTypeCode()) {
1364            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));     case(Nodes):
1365            break;     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1366         case(DegreesOfFreedom):     break;
1367            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     case(ReducedNodes):
1368            break;     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1369         case(ReducedDegreesOfFreedom):     break;
1370            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     case(Elements):
1371            break;     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1372         default:     break;
1373            stringstream temp;     case(ReducedElements):
1374            temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1375            throw FinleyAdapterException(temp.str());     break;
1376            break;     case(FaceElements):
1377    }     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1378    checkFinleyError();     break;
1379       case(ReducedFaceElements):
1380       Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1381       break;
1382       case(Points):
1383       throw FinleyAdapterException("Error - Gradient at points is not supported.");
1384       break;
1385       case(ContactElementsZero):
1386       Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1387       break;
1388       case(ReducedContactElementsZero):
1389       Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1390       break;
1391       case(ContactElementsOne):
1392       Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1393       break;
1394       case(ReducedContactElementsOne):
1395       Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1396       break;
1397       case(DegreesOfFreedom):
1398       throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1399       break;
1400       case(ReducedDegreesOfFreedom):
1401       throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1402       break;
1403       default:
1404          stringstream temp;
1405          temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1406          throw FinleyAdapterException(temp.str());
1407          break;
1408       }
1409       checkFinleyError();
1410  }  }
1411    
1412  //  //
1413  // returns the size of elements:  // returns the size of elements:
1414  //  //
1415  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1416  {  {
1417    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1418    escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1419    switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1420         case(Nodes):     case(Nodes):
1421            throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1422            break;     break;
1423         case(Elements):     case(ReducedNodes):
1424            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1425            break;     break;
1426         case(FaceElements):     case(Elements):
1427            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1428            break;     break;
1429         case(Points):     case(ReducedElements):
1430            throw FinleyAdapterException("Error - Size of point elements is not supported.");     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1431            break;     break;
1432         case(ContactElementsZero):     case(FaceElements):
1433         case(ContactElementsOne):     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1434            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     break;
1435            break;     case(ReducedFaceElements):
1436         case(DegreesOfFreedom):     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1437            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     break;
1438            break;     case(Points):
1439         case(ReducedDegreesOfFreedom):     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1440            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     break;
1441            break;     case(ContactElementsZero):
1442         default:     case(ContactElementsOne):
1443            stringstream temp;     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1444            temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();     break;
1445            throw FinleyAdapterException(temp.str());     case(ReducedContactElementsZero):
1446            break;     case(ReducedContactElementsOne):
1447    }     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1448    checkFinleyError();     break;
1449       case(DegreesOfFreedom):
1450       throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1451       break;
1452       case(ReducedDegreesOfFreedom):
1453       throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1454       break;
1455       default:
1456          stringstream temp;
1457          temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1458          throw FinleyAdapterException(temp.str());
1459          break;
1460       }
1461       checkFinleyError();
1462  }  }
1463    
1464  // sets the location of nodes:  //
1465  void MeshAdapter::setNewX(const Data& new_x)  // sets the location of nodes
1466    //
1467    void MeshAdapter::setNewX(const escript::Data& new_x)
1468  {  {
1469    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1470    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     escriptDataC tmp;
1471    if (newDomain!=*this)     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1472       throw FinleyAdapterException("Error - Illegal domain of new point locations");     if (newDomain!=*this)
1473    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1474    checkFinleyError();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1475  }         tmp = new_x.getDataC();
1476           Finley_Mesh_setCoordinates(mesh,&tmp);
1477  // saves a data array in openDX format:     } else {
1478  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1479  {         tmp = new_x_inter.getDataC();
1480      int MAX_namelength=256;         Finley_Mesh_setCoordinates(mesh,&tmp);
1481      const int num_data=boost::python::extract<int>(arg.attr("__len__")());     }
1482      char names[num_data][MAX_namelength];     checkFinleyError();
1483      char* c_names[num_data];  }
1484      escriptDataC data[num_data];  
1485      escriptDataC* ptr_data[num_data];  //
1486    // Helper for the save* methods. Extracts optional data variable names and the
1487      boost::python::list keys=arg.keys();  // corresponding pointers from python dictionary. Caller must free arrays.
1488      for (int i=0;i<num_data;++i) {  //
1489           Data& d=boost::python::extract<Data&>(arg[keys[i]]);  void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1490           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  {
1491               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");     numData = boost::python::extract<int>(arg.attr("__len__")());
1492           data[i]=d.getDataC();     /* win32 refactor */
1493           ptr_data[i]=&(data[i]);     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1494           std::string n=boost::python::extract<std::string>(keys[i]);     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1495           c_names[i]=&(names[i][0]);     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1496           if (n.length()>MAX_namelength-1) {  
1497              strncpy(c_names[i],n.c_str(),MAX_namelength-1);     boost::python::list keys=arg.keys();
1498              c_names[i][MAX_namelength-1]='\0';     for (int i=0; i<numData; ++i) {
1499           } else {        string n=boost::python::extract<string>(keys[i]);
1500              strcpy(c_names[i],n.c_str());        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1501           }        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1502      }           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1503      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);        data[i] = d.getDataC();
1504      checkFinleyError();        dataPtr[i] = &(data[i]);
1505      return;        names[i] = TMPMEMALLOC(n.length()+1, char);
1506  }        strcpy(names[i], n.c_str());
1507       }
1508  // saves a data array in openVTK format:  }
1509  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  
1510  {  //
1511      int MAX_namelength=256;  // saves mesh and optionally data arrays in openDX format
1512      const int num_data=boost::python::extract<int>(arg.attr("__len__")());  //
1513      char names[num_data][MAX_namelength];  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1514      char* c_names[num_data];  {
1515      escriptDataC data[num_data];     int num_data;
1516      escriptDataC* ptr_data[num_data];     char **names;
1517       escriptDataC *data;
1518      boost::python::list keys=arg.keys();     escriptDataC **ptr_data;
1519      for (int i=0;i<num_data;++i) {  
1520           Data& d=boost::python::extract<Data&>(arg[keys[i]]);     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1521           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1522               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");     checkFinleyError();
1523           data[i]=d.getDataC();  
1524           ptr_data[i]=&(data[i]);     /* win32 refactor */
1525           std::string n=boost::python::extract<std::string>(keys[i]);     TMPMEMFREE(data);
1526           c_names[i]=&(names[i][0]);     TMPMEMFREE(ptr_data);
1527           if (n.length()>MAX_namelength-1) {     for(int i=0; i<num_data; i++)
1528              strncpy(c_names[i],n.c_str(),MAX_namelength-1);     {
1529              c_names[i][MAX_namelength-1]='\0';        TMPMEMFREE(names[i]);
1530           } else {     }
1531              strcpy(c_names[i],n.c_str());     TMPMEMFREE(names);
1532           }  
1533      }     return;
1534      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  }
1535      checkFinleyError();  
1536      return;  //
1537  }  // saves mesh and optionally data arrays in VTK format
1538                                                                                                                                                                            //
1539                                                                                                                                                                            void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1540  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  {
1541  SystemMatrixAdapter MeshAdapter::newSystemMatrix(     int num_data;
1542                        const int row_blocksize,     char **names;
1543                        const escript::FunctionSpace& row_functionspace,     escriptDataC *data;
1544                        const int column_blocksize,     escriptDataC **ptr_data;
1545                        const escript::FunctionSpace& column_functionspace,  
1546                        const int type) const     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1547  {     Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1548      int reduceRowOrder=0;     checkFinleyError();
1549      int reduceColOrder=0;  
1550      // is the domain right?     /* win32 refactor */
1551      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     TMPMEMFREE(data);
1552      if (row_domain!=*this)     TMPMEMFREE(ptr_data);
1553            throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");     for(int i=0; i<num_data; i++)
1554      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     {
1555      if (col_domain!=*this)        TMPMEMFREE(names[i]);
1556            throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");     }
1557      // is the function space type right     TMPMEMFREE(names);
1558      if (row_functionspace.getTypeCode()==DegreesOfFreedom) {  }
1559          reduceRowOrder=0;  
1560      } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1561          reduceRowOrder=1;  {
1562      } else {  #ifdef ESYS_MPI
1563          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");      index_t myFirstNode=0, myLastNode=0, k=0;
1564        index_t* globalNodeIndex=0;
1565        Finley_Mesh* mesh_p=m_finleyMesh.get();
1566        if (fs_code == FINLEY_REDUCED_NODES)
1567        {
1568        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1569        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1570        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1571      }      }
1572      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {      else
1573          reduceColOrder=0;      {
1574      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1575          reduceColOrder=1;      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1576      } else {      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
         throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");  
1577      }      }
1578      // generate matrix:      k=globalNodeIndex[id];
1579            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1580      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);  #endif
1581      checkFinleyError();      return true;
1582      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);  }
1583      checkPasoError();  
1584      Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);  
1585      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);  
1586    //
1587    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1588    //
1589    ASM_ptr MeshAdapter::newSystemMatrix(
1590                                                     const int row_blocksize,
1591                                                     const escript::FunctionSpace& row_functionspace,
1592                                                     const int column_blocksize,
1593                                                     const escript::FunctionSpace& column_functionspace,
1594                                                     const int type) const
1595    {
1596       int reduceRowOrder=0;
1597       int reduceColOrder=0;
1598       // is the domain right?
1599       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1600       if (row_domain!=*this)
1601          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1602       const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1603       if (col_domain!=*this)
1604          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1605       // is the function space type right
1606       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1607          reduceRowOrder=0;
1608       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1609          reduceRowOrder=1;
1610       } else {
1611          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1612       }
1613       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1614          reduceColOrder=0;
1615       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1616          reduceColOrder=1;
1617       } else {
1618          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1619       }
1620       // generate matrix:
1621    
1622       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1623       checkFinleyError();
1624       Paso_SystemMatrix* fsystemMatrix;
1625       int trilinos = 0;
1626       if (trilinos) {
1627    #ifdef TRILINOS
1628          /* Allocation Epetra_VrbMatrix here */
1629    #endif
1630       }
1631       else {
1632          fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1633       }
1634       checkPasoError();
1635       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1636       SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1637       return ASM_ptr(sma);
1638    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1639    }
1640    
1641    //
1642    // creates a TransportProblemAdapter
1643    //
1644    ATP_ptr MeshAdapter::newTransportProblem(
1645                                                             const bool useBackwardEuler,
1646                                                             const int blocksize,
1647                                                             const escript::FunctionSpace& functionspace,
1648                                                             const int type) const
1649    {
1650       int reduceOrder=0;
1651       // is the domain right?
1652       const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1653       if (domain!=*this)
1654          throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1655       // is the function space type right
1656       if (functionspace.getTypeCode()==DegreesOfFreedom) {
1657          reduceOrder=0;
1658       } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1659          reduceOrder=1;
1660       } else {
1661          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1662       }
1663       // generate matrix:
1664    
1665       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1666       checkFinleyError();
1667       Paso_TransportProblem* transportProblem;
1668       transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1669       checkPasoError();
1670       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1671       TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1672       return ATP_ptr(tpa);
1673    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1674  }  }
1675    
1676  //  //
# Line 919  SystemMatrixAdapter MeshAdapter::newSyst Line 1682  SystemMatrixAdapter MeshAdapter::newSyst
1682  // 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:
1683  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1684  {  {
1685    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1686         case(Nodes):     case(Nodes):
1687         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1688         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1689            return false;     return false;
1690            break;     break;
1691         case(Elements):     case(Elements):
1692         case(FaceElements):     case(FaceElements):
1693         case(Points):     case(Points):
1694         case(ContactElementsZero):     case(ContactElementsZero):
1695         case(ContactElementsOne):     case(ContactElementsOne):
1696            return true;     case(ReducedElements):
1697            break;     case(ReducedFaceElements):
1698         default:     case(ReducedContactElementsZero):
1699            stringstream temp;     case(ReducedContactElementsOne):
1700            temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;     return true;
1701            throw FinleyAdapterException(temp.str());     break;
1702            break;     default:
1703    }        stringstream temp;
1704    checkFinleyError();        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1705    return false;        throw FinleyAdapterException(temp.str());
1706          break;
1707       }
1708       checkFinleyError();
1709       return false;
1710    }
1711    
1712    bool
1713    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1714    {
1715       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1716        class 1: DOF <-> Nodes
1717        class 2: ReducedDOF <-> ReducedNodes
1718        class 3: Points
1719        class 4: Elements
1720        class 5: ReducedElements
1721        class 6: FaceElements
1722        class 7: ReducedFaceElements
1723        class 8: ContactElementZero <-> ContactElementOne
1724        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1725    
1726       There is also a set of lines. Interpolation is possible down a line but not between lines.
1727       class 1 and 2 belong to all lines so aren't considered.
1728        line 0: class 3
1729        line 1: class 4,5
1730        line 2: class 6,7
1731        line 3: class 8,9
1732    
1733       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1734       eg hasnodes is true if we have at least one instance of Nodes.
1735       */
1736        if (fs.empty())
1737        {
1738            return false;
1739        }
1740        vector<int> hasclass(10);
1741        vector<int> hasline(4);
1742        bool hasnodes=false;
1743        bool hasrednodes=false;
1744        bool hascez=false;
1745        bool hasrcez=false;
1746        for (int i=0;i<fs.size();++i)
1747        {
1748        switch(fs[i])
1749        {
1750        case(Nodes):   hasnodes=true;   // no break is deliberate
1751        case(DegreesOfFreedom):
1752            hasclass[1]=1;
1753            break;
1754        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1755        case(ReducedDegreesOfFreedom):
1756            hasclass[2]=1;
1757            break;
1758        case(Points):
1759            hasline[0]=1;
1760            hasclass[3]=1;
1761            break;
1762        case(Elements):
1763            hasclass[4]=1;
1764            hasline[1]=1;
1765            break;
1766        case(ReducedElements):
1767            hasclass[5]=1;
1768            hasline[1]=1;
1769            break;
1770        case(FaceElements):
1771            hasclass[6]=1;
1772            hasline[2]=1;
1773            break;
1774        case(ReducedFaceElements):
1775            hasclass[7]=1;
1776            hasline[2]=1;
1777            break;
1778        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1779        case(ContactElementsOne):
1780            hasclass[8]=1;
1781            hasline[3]=1;
1782            break;
1783        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1784        case(ReducedContactElementsOne):
1785            hasclass[9]=1;
1786            hasline[3]=1;
1787            break;
1788        default:
1789            return false;
1790        }
1791        }
1792        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1793        // fail if we have more than one leaf group
1794    
1795        if (totlines>1)
1796        {
1797        return false;   // there are at least two branches we can't interpolate between
1798        }
1799        else if (totlines==1)
1800        {
1801        if (hasline[0]==1)      // we have points
1802        {
1803            resultcode=Points;
1804        }
1805        else if (hasline[1]==1)
1806        {
1807            if (hasclass[5]==1)
1808            {
1809            resultcode=ReducedElements;
1810            }
1811            else
1812            {
1813            resultcode=Elements;
1814            }
1815        }
1816        else if (hasline[2]==1)
1817        {
1818            if (hasclass[7]==1)
1819            {
1820            resultcode=ReducedFaceElements;
1821            }
1822            else
1823            {
1824            resultcode=FaceElements;
1825            }
1826        }
1827        else    // so we must be in line3
1828        {
1829            if (hasclass[9]==1)
1830            {
1831            // need something from class 9
1832            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1833            }
1834            else
1835            {
1836            // something from class 8
1837            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1838            }
1839        }
1840        }
1841        else    // totlines==0
1842        {
1843        if (hasclass[2]==1)
1844        {
1845            // something from class 2
1846            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1847        }
1848        else
1849        {   // something from class 1
1850            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1851        }
1852        }
1853        return true;
1854  }  }
1855    
1856  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1857  {  {
1858    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1859       case(Nodes):     case(Nodes):
1860          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1861             case(Nodes):      case(Nodes):
1862             case(ReducedDegreesOfFreedom):      case(ReducedNodes):
1863             case(DegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1864             case(Elements):      case(DegreesOfFreedom):
1865             case(FaceElements):      case(Elements):
1866             case(Points):      case(ReducedElements):
1867             case(ContactElementsZero):      case(FaceElements):
1868             case(ContactElementsOne):      case(ReducedFaceElements):
1869                 return true;      case(Points):
1870             default:      case(ContactElementsZero):
1871                 stringstream temp;      case(ReducedContactElementsZero):
1872                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;      case(ContactElementsOne):
1873                 throw FinleyAdapterException(temp.str());      case(ReducedContactElementsOne):
1874          }      return true;
1875          break;      default:
1876       case(Elements):            stringstream temp;
1877          if (functionSpaceType_target==Elements) {            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1878             return true;            throw FinleyAdapterException(temp.str());
1879          } else {     }
1880             return false;     break;
1881          }     case(ReducedNodes):
1882       case(FaceElements):      switch(functionSpaceType_target) {
1883          if (functionSpaceType_target==FaceElements) {      case(ReducedNodes):
1884             return true;      case(ReducedDegreesOfFreedom):
1885          } else {      case(Elements):
1886             return false;      case(ReducedElements):
1887          }      case(FaceElements):
1888       case(Points):      case(ReducedFaceElements):
1889          if (functionSpaceType_target==Points) {      case(Points):
1890             return true;      case(ContactElementsZero):
1891          } else {      case(ReducedContactElementsZero):
1892             return false;      case(ContactElementsOne):
1893          }      case(ReducedContactElementsOne):
1894       case(ContactElementsZero):      return true;
1895       case(ContactElementsOne):      case(Nodes):
1896          if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      case(DegreesOfFreedom):
1897             return true;      return false;
1898        default:
1899            stringstream temp;
1900            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1901            throw FinleyAdapterException(temp.str());
1902       }
1903       break;
1904       case(Elements):
1905        if (functionSpaceType_target==Elements) {
1906          return true;
1907        } else if (functionSpaceType_target==ReducedElements) {
1908          return true;
1909          } else {          } else {
1910             return false;            return false;
         }  
      case(DegreesOfFreedom):  
         switch(functionSpaceType_target) {  
            case(ReducedDegreesOfFreedom):  
            case(DegreesOfFreedom):  
            case(Nodes):  
            case(Elements):  
            case(FaceElements):  
            case(Points):  
            case(ContactElementsZero):  
            case(ContactElementsOne):  
               return true;  
            default:  
              stringstream temp;  
              temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
              throw FinleyAdapterException(temp.str());  
1911          }          }
1912          break;     case(ReducedElements):
1913       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1914         switch(functionSpaceType_target) {        return true;
1915            case(ReducedDegreesOfFreedom):      } else {
1916            case(Elements):            return false;
1917            case(FaceElements):      }
1918            case(Points):     case(FaceElements):
1919            case(ContactElementsZero):      if (functionSpaceType_target==FaceElements) {
1920            case(ContactElementsOne):              return true;
1921                return true;      } else if (functionSpaceType_target==ReducedFaceElements) {
1922            case(Nodes):              return true;
1923            case(DegreesOfFreedom):      } else {
1924               return false;              return false;
1925            default:      }
1926               stringstream temp;     case(ReducedFaceElements):
1927               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;      if (functionSpaceType_target==ReducedFaceElements) {
1928               throw FinleyAdapterException(temp.str());              return true;
1929         }      } else {
1930         break;          return false;
1931       default:      }
1932          stringstream temp;     case(Points):
1933          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;      if (functionSpaceType_target==Points) {
1934          throw FinleyAdapterException(temp.str());              return true;
1935          break;      } else {
1936    }              return false;
1937    checkFinleyError();      }
1938    return false;     case(ContactElementsZero):
1939       case(ContactElementsOne):
1940        if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1941                return true;
1942        } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1943                return true;
1944        } else {
1945                return false;
1946        }
1947       case(ReducedContactElementsZero):
1948       case(ReducedContactElementsOne):
1949        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1950                return true;
1951        } else {
1952                return false;
1953        }
1954       case(DegreesOfFreedom):
1955        switch(functionSpaceType_target) {
1956        case(ReducedDegreesOfFreedom):
1957        case(DegreesOfFreedom):
1958        case(Nodes):
1959        case(ReducedNodes):
1960        case(Elements):
1961        case(ReducedElements):
1962        case(Points):
1963        case(FaceElements):
1964        case(ReducedFaceElements):
1965        case(ContactElementsZero):
1966        case(ReducedContactElementsZero):
1967        case(ContactElementsOne):
1968        case(ReducedContactElementsOne):
1969        return true;
1970        default:
1971            stringstream temp;
1972            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1973            throw FinleyAdapterException(temp.str());
1974        }
1975        break;
1976       case(ReducedDegreesOfFreedom):
1977       switch(functionSpaceType_target) {
1978        case(ReducedDegreesOfFreedom):
1979        case(ReducedNodes):
1980        case(Elements):
1981        case(ReducedElements):
1982        case(FaceElements):
1983        case(ReducedFaceElements):
1984        case(Points):
1985        case(ContactElementsZero):
1986        case(ReducedContactElementsZero):
1987        case(ContactElementsOne):
1988        case(ReducedContactElementsOne):
1989        return true;
1990        case(Nodes):
1991        case(DegreesOfFreedom):
1992        return false;
1993        default:
1994            stringstream temp;
1995            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1996            throw FinleyAdapterException(temp.str());
1997        }
1998        break;
1999       default:
2000          stringstream temp;
2001          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
2002          throw FinleyAdapterException(temp.str());
2003          break;
2004       }
2005       checkFinleyError();
2006       return false;
2007  }  }
2008    
2009  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
2010  {  {
2011      return false;     return false;
2012  }  }
2013    
2014  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2015  {  {
2016    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2017    if (temp!=0) {     if (temp!=0) {
2018      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2019    } else {     } else {
2020      return false;        return false;
2021    }     }
2022  }  }
2023    
2024  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2025  {  {
2026    return !(operator==(other));     return !(operator==(other));
2027  }  }
2028    
2029  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
2030  {  {
2031     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2032       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2033       checkPasoError();
2034       return out;
2035    }
2036    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2037    {
2038       Finley_Mesh* mesh=m_finleyMesh.get();
2039       int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2040     checkPasoError();     checkPasoError();
2041     return out;     return out;
2042  }  }
2043    
2044  Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2045  {  {
2046    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2047  }  }
2048    
2049  Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2050  {  {
2051    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2052  }  }
2053    
2054  Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2055  {  {
2056    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2057  }  }
2058    
2059    const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2060    {
2061       int *out = NULL;
2062       Finley_Mesh* mesh=m_finleyMesh.get();
2063       switch (functionSpaceType) {
2064       case(Nodes):
2065       out=mesh->Nodes->Id;
2066       break;
2067       case(ReducedNodes):
2068       out=mesh->Nodes->reducedNodesId;
2069       break;
2070       case(Elements):
2071       out=mesh->Elements->Id;
2072       break;
2073       case(ReducedElements):
2074       out=mesh->Elements->Id;
2075       break;
2076       case(FaceElements):
2077       out=mesh->FaceElements->Id;
2078       break;
2079       case(ReducedFaceElements):
2080       out=mesh->FaceElements->Id;
2081       break;
2082       case(Points):
2083       out=mesh->Points->Id;
2084       break;
2085       case(ContactElementsZero):
2086       out=mesh->ContactElements->Id;
2087       break;
2088       case(ReducedContactElementsZero):
2089       out=mesh->ContactElements->Id;
2090       break;
2091       case(ContactElementsOne):
2092       out=mesh->ContactElements->Id;
2093       break;
2094       case(ReducedContactElementsOne):
2095       out=mesh->ContactElements->Id;
2096       break;
2097       case(DegreesOfFreedom):
2098       out=mesh->Nodes->degreesOfFreedomId;
2099       break;
2100       case(ReducedDegreesOfFreedom):
2101       out=mesh->Nodes->reducedDegreesOfFreedomId;
2102       break;
2103       default:
2104          stringstream temp;
2105          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2106          throw FinleyAdapterException(temp.str());
2107          break;
2108       }
2109       return out;
2110    }
2111  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2112  {  {
2113    int* tagList;     int out=0;
2114    int numTags;     Finley_Mesh* mesh=m_finleyMesh.get();
2115    getTagList(functionSpaceType, &tagList, &numTags);     switch (functionSpaceType) {
2116    return tagList[sampleNo];     case(Nodes):
2117       out=mesh->Nodes->Tag[sampleNo];
2118       break;
2119       case(ReducedNodes):
2120       throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2121       break;
2122       case(Elements):
2123       out=mesh->Elements->Tag[sampleNo];
2124       break;
2125       case(ReducedElements):
2126       out=mesh->Elements->Tag[sampleNo];
2127       break;
2128       case(FaceElements):
2129       out=mesh->FaceElements->Tag[sampleNo];
2130       break;
2131       case(ReducedFaceElements):
2132       out=mesh->FaceElements->Tag[sampleNo];
2133       break;
2134       case(Points):
2135       out=mesh->Points->Tag[sampleNo];
2136       break;
2137       case(ContactElementsZero):
2138       out=mesh->ContactElements->Tag[sampleNo];
2139       break;
2140       case(ReducedContactElementsZero):
2141       out=mesh->ContactElements->Tag[sampleNo];
2142       break;
2143       case(ContactElementsOne):
2144       out=mesh->ContactElements->Tag[sampleNo];
2145       break;
2146       case(ReducedContactElementsOne):
2147       out=mesh->ContactElements->Tag[sampleNo];
2148       break;
2149       case(DegreesOfFreedom):
2150       throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2151       break;
2152       case(ReducedDegreesOfFreedom):
2153       throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2154       break;
2155       default:
2156          stringstream temp;
2157          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2158          throw FinleyAdapterException(temp.str());
2159          break;
2160       }
2161       return out;
2162    }
2163    
2164    
2165    void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2166    {
2167       Finley_Mesh* mesh=m_finleyMesh.get();
2168       escriptDataC tmp=mask.getDataC();
2169       switch(functionSpaceType) {
2170       case(Nodes):
2171       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2172       break;
2173       case(ReducedNodes):
2174       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2175       break;
2176       case(DegreesOfFreedom):
2177       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2178       break;
2179       case(ReducedDegreesOfFreedom):
2180       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2181       break;
2182       case(Elements):
2183       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2184       break;
2185       case(ReducedElements):
2186       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2187       break;
2188       case(FaceElements):
2189       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2190       break;
2191       case(ReducedFaceElements):
2192       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2193       break;
2194       case(Points):
2195       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2196       break;
2197       case(ContactElementsZero):
2198       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2199       break;
2200       case(ReducedContactElementsZero):
2201       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2202       break;
2203       case(ContactElementsOne):
2204       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2205       break;
2206       case(ReducedContactElementsOne):
2207       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2208       break;
2209       default:
2210          stringstream temp;
2211          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2212          throw FinleyAdapterException(temp.str());
2213       }
2214       checkFinleyError();
2215       return;
2216    }
2217    
2218    void MeshAdapter::setTagMap(const string& name,  int tag)
2219    {
2220       Finley_Mesh* mesh=m_finleyMesh.get();
2221       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2222       checkPasoError();
2223       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2224    }
2225    
2226    int MeshAdapter::getTag(const string& name) const
2227    {
2228       Finley_Mesh* mesh=m_finleyMesh.get();
2229       int tag=0;
2230       tag=Finley_Mesh_getTag(mesh, name.c_str());
2231       checkPasoError();
2232       // throwStandardException("MeshAdapter::getTag is not implemented.");
2233       return tag;
2234    }
2235    
2236    bool MeshAdapter::isValidTagName(const string& name) const
2237    {
2238       Finley_Mesh* mesh=m_finleyMesh.get();
2239       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2240    }
2241    
2242    string MeshAdapter::showTagNames() const
2243    {
2244       stringstream temp;
2245       Finley_Mesh* mesh=m_finleyMesh.get();
2246       Finley_TagMap* tag_map=mesh->TagMap;
2247       while (tag_map) {
2248          temp << tag_map->name;
2249          tag_map=tag_map->next;
2250          if (tag_map) temp << ", ";
2251       }
2252       return temp.str();
2253    }
2254    
2255    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2256    {
2257      Finley_Mesh* mesh=m_finleyMesh.get();
2258      dim_t numTags=0;
2259      switch(functionSpaceCode) {
2260       case(Nodes):
2261              numTags=mesh->Nodes->numTagsInUse;
2262              break;
2263       case(ReducedNodes):
2264              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2265              break;
2266       case(DegreesOfFreedom):
2267              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2268              break;
2269       case(ReducedDegreesOfFreedom):
2270              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2271              break;
2272       case(Elements):
2273       case(ReducedElements):
2274              numTags=mesh->Elements->numTagsInUse;
2275              break;
2276       case(FaceElements):
2277       case(ReducedFaceElements):
2278              numTags=mesh->FaceElements->numTagsInUse;
2279              break;
2280       case(Points):
2281              numTags=mesh->Points->numTagsInUse;
2282              break;
2283       case(ContactElementsZero):
2284       case(ReducedContactElementsZero):
2285       case(ContactElementsOne):
2286       case(ReducedContactElementsOne):
2287              numTags=mesh->ContactElements->numTagsInUse;
2288              break;
2289       default:
2290          stringstream temp;
2291          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2292          throw FinleyAdapterException(temp.str());
2293      }
2294      return numTags;
2295    }
2296    
2297    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2298    {
2299      Finley_Mesh* mesh=m_finleyMesh.get();
2300      index_t* tags=NULL;
2301      switch(functionSpaceCode) {
2302       case(Nodes):
2303              tags=mesh->Nodes->tagsInUse;
2304              break;
2305       case(ReducedNodes):
2306              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2307              break;
2308       case(DegreesOfFreedom):
2309              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2310              break;
2311       case(ReducedDegreesOfFreedom):
2312              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2313              break;
2314       case(Elements):
2315       case(ReducedElements):
2316              tags=mesh->Elements->tagsInUse;
2317              break;
2318       case(FaceElements):
2319       case(ReducedFaceElements):
2320              tags=mesh->FaceElements->tagsInUse;
2321              break;
2322       case(Points):
2323              tags=mesh->Points->tagsInUse;
2324              break;
2325       case(ContactElementsZero):
2326       case(ReducedContactElementsZero):
2327       case(ContactElementsOne):
2328       case(ReducedContactElementsOne):
2329              tags=mesh->ContactElements->tagsInUse;
2330              break;
2331       default:
2332          stringstream temp;
2333          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2334          throw FinleyAdapterException(temp.str());
2335      }
2336      return tags;
2337    }
2338    
2339    
2340    bool MeshAdapter::canTag(int functionSpaceCode) const
2341    {
2342      switch(functionSpaceCode) {
2343       case(Nodes):
2344       case(Elements):
2345       case(ReducedElements):
2346       case(FaceElements):
2347       case(ReducedFaceElements):
2348       case(Points):
2349       case(ContactElementsZero):
2350       case(ReducedContactElementsZero):
2351       case(ContactElementsOne):
2352       case(ReducedContactElementsOne):
2353              return true;
2354       case(ReducedNodes):
2355       case(DegreesOfFreedom):
2356       case(ReducedDegreesOfFreedom):
2357          return false;
2358       default:
2359        return false;
2360      }
2361  }  }
2362    
2363  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  AbstractDomain::StatusType MeshAdapter::getStatus() const
2364  {  {
2365    int* referenceNoList;    Finley_Mesh* mesh=m_finleyMesh.get();
2366    int numReferenceNo;    return Finley_Mesh_getStatus(mesh);
2367    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);  }
2368    return referenceNoList[sampleNo];  
2369    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2370    {
2371      
2372      Finley_Mesh* mesh=m_finleyMesh.get();
2373      int order =-1;
2374      switch(functionSpaceCode) {
2375       case(Nodes):
2376       case(DegreesOfFreedom):
2377              order=mesh->approximationOrder;
2378              break;
2379       case(ReducedNodes):
2380       case(ReducedDegreesOfFreedom):
2381              order=mesh->reducedApproximationOrder;
2382              break;
2383       case(Elements):
2384       case(FaceElements):
2385       case(Points):
2386       case(ContactElementsZero):
2387       case(ContactElementsOne):
2388              order=mesh->integrationOrder;
2389              break;
2390       case(ReducedElements):
2391       case(ReducedFaceElements):
2392       case(ReducedContactElementsZero):
2393       case(ReducedContactElementsOne):
2394              order=mesh->reducedIntegrationOrder;
2395              break;
2396       default:
2397          stringstream temp;
2398          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2399          throw FinleyAdapterException(temp.str());
2400      }
2401      return order;
2402    }
2403    
2404    bool MeshAdapter::supportsContactElements() const
2405    {
2406      return true;
2407    }
2408    
2409    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2410    {
2411      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2412    }
2413    
2414    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2415    {
2416      Finley_ReferenceElementSet_dealloc(m_refSet);
2417  }  }
2418    
2419  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.201  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26