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

Legend:
Removed from v.102  
changed lines
  Added in v.1464

  ViewVC Help
Powered by ViewVC 1.1.26