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

Legend:
Removed from v.149  
changed lines
  Added in v.2635

  ViewVC Help
Powered by ViewVC 1.1.26