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

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

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

revision 532 by gross, Wed Feb 15 09:45:53 2006 UTC 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"  #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" {
23    #include "escript/blocktimer.h"
24    }
25    #include <vector>
26    
27  #include "Data.h"  #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256
 #include "DataFactory.h"  
28    
29  using namespace std;  using namespace std;
30  using namespace escript;  using namespace escript;
# Line 29  MeshAdapter::FunctionSpaceNamesMapType M Line 37  MeshAdapter::FunctionSpaceNamesMapType M
37  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
38  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
39  const int MeshAdapter::Nodes=FINLEY_NODES;  const int MeshAdapter::Nodes=FINLEY_NODES;
40    const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
41  const int MeshAdapter::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=FINLEY_ELEMENTS;
42    const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
43  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
44    const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
45  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=FINLEY_POINTS;
46  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
47    const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
48  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
49    const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
50    
51  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
52  {  {
53    setFunctionSpaceTypeNames();    setFunctionSpaceTypeNames();
54    //    //
55    // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer    // need to use a null_deleter as Finley_Mesh_free deletes the pointer
56    // for us.    // for us.
57    m_finleyMesh.reset(finleyMesh,null_deleter());    m_finleyMesh.reset(finleyMesh,null_deleter());
58  }  }
# Line 58  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::Print_Mesh_Info(const bool full=false) const
102    {
103      Finley_PrintMesh_Info(m_finleyMesh.get(), full);
104  }  }
105    
106  // void MeshAdapter::getTagList(int functionSpaceType,  void MeshAdapter::dump(const std::string& fileName) const
107  //                  int* numTags) const  {
108  // {  #ifdef USE_NETCDF
109  //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);     const NcDim* ncdims[12];
110  //   return;     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 115  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  {  {
# Line 167  int MeshAdapter::getDiracDeltaFunctionCo Line 636  int MeshAdapter::getDiracDeltaFunctionCo
636  }  }
637    
638  //  //
 // returns a pointer to the tag list of samples of functionSpaceType  
 //  
 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,  
                  int* numTags) const  
 {  
   *tagList=NULL;  
   *numTags=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *tagList=mesh->Nodes->Tag;  
       *numTags=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *tagList=mesh->Elements->Tag;  
       *numTags=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *tagList=mesh->FaceElements->Tag;  
       *numTags=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *tagList=mesh->Points->Tag;  
       *numTags=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*tagList==NULL) {  
     stringstream temp;  
     temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
 }  
   
 //  
 // returns a pointer to the reference no list of samples of functionSpaceType  
 //  
 void MeshAdapter::getReferenceNoList(int functionSpaceType, int** referenceNoList,  
                  int* numReferenceNo) const  
 {  
   *referenceNoList=NULL;  
   *numReferenceNo=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=mesh->Nodes->Id;  
       *numReferenceNo=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *referenceNoList=mesh->Elements->Id;  
       *numReferenceNo=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *referenceNoList=mesh->FaceElements->Id;  
       *numReferenceNo=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *referenceNoList=mesh->Points->Id;  
       *numReferenceNo=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*referenceNoList==NULL) {  
     stringstream temp;  
     temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
 }  
   
 //  
639  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
640  //  //
641  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 332  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 340  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 358  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:
# Line 394  void MeshAdapter::addPDEToSystem( Line 747  void MeshAdapter::addPDEToSystem(
747                       const escript::Data& d, const escript::Data& y,                       const escript::Data& d, const escript::Data& y,
748                       const escript::Data& d_contact,const escript::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.getPaso_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();     checkFinleyError();
766     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
767                    mat.getPaso_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
                   &(rhs.getDataC()),  
                                   &(d.getDataC()),&(y.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Mean_out);  
768     checkFinleyError();     checkFinleyError();
769     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  
770                    mat.getPaso_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
                   &(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( escript::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  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  
796  {  {
797     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
798    
799     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
800     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     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();
807    
808       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
809       checkFinleyError();
810    
811       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
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_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Mesh* mesh=m_finleyMesh.get();
840     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     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();     checkFinleyError();
850     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
851     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
852     checkFinleyError();     checkFinleyError();
853  }  }
854    
# Line 441  void MeshAdapter::interpolateOnDomain(es Line 859  void MeshAdapter::interpolateOnDomain(es
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                 stringstream temp;                 stringstream temp;
# Line 479  void MeshAdapter::interpolateOnDomain(es Line 938  void MeshAdapter::interpolateOnDomain(es
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 {          } else {
945             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
946          }          }
947          break;          break;
948         case(ReducedElements):
949            if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
950               Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
951            } else {
952               throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
953            }
954            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             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
962             break;          }
963            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             throw FinleyAdapterException("Error - No interpolation with data on points possible.");             throw FinleyAdapterException("Error - 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             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");             throw FinleyAdapterException("Error - 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               stringstream temp;               stringstream temp;
# Line 537  void MeshAdapter::interpolateOnDomain(es Line 1061  void MeshAdapter::interpolateOnDomain(es
1061          break;          break;
1062       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1063         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
1064              case(Nodes):
1065                 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):            case(ReducedDegreesOfFreedom):
1081               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               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(Nodes):                if (getMPISize()>1) {
1117               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");                   escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1118               break;                   escriptDataC _in2 = temp.getDataC();
1119            case(DegreesOfFreedom):                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1120               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");                } else {
1121                     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1122                 }
1123               break;               break;
1124            default:            default:
1125               stringstream temp;               stringstream temp;
# Line 586  void MeshAdapter::setToX(escript::Data& Line 1148  void MeshAdapter::setToX(escript::Data&
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       escript::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    }    }
# Line 605  void MeshAdapter::setToNormal(escript::D Line 1169  void MeshAdapter::setToNormal(escript::D
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        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1176        break;        break;
1177        case(ReducedNodes):
1178          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1179          break;
1180      case(Elements):      case(Elements):
1181        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1182        break;        break;
1183        case(ReducedElements):
1184          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1185          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        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");        throw FinleyAdapterException("Error - 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        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
# Line 658  void MeshAdapter::setToIntegrals(std::ve Line 1236  void MeshAdapter::setToIntegrals(std::ve
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          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1247            _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          throw FinleyAdapterException("Error - Integral of data on points is not supported.");          throw FinleyAdapterException("Error - 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          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1284            _temp=temp.getDataC();
1285            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1286          break;          break;
1287       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1288          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1289            _temp=temp.getDataC();
1290            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1291          break;          break;
1292       default:       default:
1293          stringstream temp;          stringstream temp;
# Line 691  void MeshAdapter::setToIntegrals(std::ve Line 1296  void MeshAdapter::setToIntegrals(std::ve
1296          break;          break;
1297    }    }
1298    checkFinleyError();    checkFinleyError();
1299      blocktimer_increment("integrate()", blocktimer_start);
1300  }  }
1301    
1302  //  //
# Line 706  void MeshAdapter::setToGradient(escript: Line 1312  void MeshAdapter::setToGradient(escript:
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            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1334            break;            break;
1335           case(ReducedNodes):
1336              throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1337              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            throw FinleyAdapterException("Error - Gradient at points is not supported.");            throw FinleyAdapterException("Error - 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            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
# Line 751  void MeshAdapter::setToSize(escript::Dat Line 1388  void MeshAdapter::setToSize(escript::Dat
1388         case(Nodes):         case(Nodes):
1389            throw FinleyAdapterException("Error - Size of nodes is not supported.");            throw FinleyAdapterException("Error - Size of nodes is not supported.");
1390            break;            break;
1391           case(ReducedNodes):
1392              throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1393              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            throw FinleyAdapterException("Error - Size of point elements is not supported.");            throw FinleyAdapterException("Error - Size of point elements is not supported.");
1408            break;            break;
# Line 764  void MeshAdapter::setToSize(escript::Dat Line 1410  void MeshAdapter::setToSize(escript::Dat
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            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1419            break;            break;
# Line 783  void MeshAdapter::setToSize(escript::Dat Line 1433  void MeshAdapter::setToSize(escript::Dat
1433  void MeshAdapter::setNewX(const escript::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_Mesh_setCoordinates(mesh,&(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 boost::python::dict& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1447  {  {
1448      int MAX_namelength=256;      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__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1450      char names[num_data][MAX_namelength];    /* win32 refactor */
1451      char* c_names[num_data];    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1452      escriptDataC data[num_data];    for(int i=0;i<num_data;i++)
1453      escriptDataC* ptr_data[num_data];    {
1454        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1455      }
1456    
1457      boost::python::list keys=arg.keys();    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1458      for (int i=0;i<num_data;++i) {    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]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1465           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1466               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
1467           data[i]=d.getDataC();           data[i]=d.getDataC();
1468           ptr_data[i]=&(data[i]);           ptr_data[i]=&(data[i]);
          std::string n=boost::python::extract<std::string>(keys[i]);  
1469           c_names[i]=&(names[i][0]);           c_names[i]=&(names[i][0]);
1470           if (n.length()>MAX_namelength-1) {           if (n.length()>MAX_namelength-1) {
1471              strncpy(c_names[i],n.c_str(),MAX_namelength-1);              strncpy(c_names[i],n.c_str(),MAX_namelength-1);
# Line 818  void MeshAdapter::saveDX(const std::stri Line 1476  void MeshAdapter::saveDX(const std::stri
1476      }      }
1477      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1478      checkFinleyError();      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;      return;
1491  }  }
1492    
1493  // saves a data array in openVTK format:  // saves a data array in openVTK format:
1494  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1495  {  {
1496      int MAX_namelength=256;      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__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1498      char names[num_data][MAX_namelength];    /* win32 refactor */
1499      char* c_names[num_data];    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1500      escriptDataC data[num_data];    for(int i=0;i<num_data;i++)
1501      escriptDataC* ptr_data[num_data];    {
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();      boost::python::list keys=arg.keys();
1510      for (int i=0;i<num_data;++i) {      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]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1513           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1514               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
1515           data[i]=d.getDataC();           data[i]=d.getDataC();
1516           ptr_data[i]=&(data[i]);           ptr_data[i]=&(data[i]);
          std::string n=boost::python::extract<std::string>(keys[i]);  
1517           c_names[i]=&(names[i][0]);           c_names[i]=&(names[i][0]);
1518           if (n.length()>MAX_namelength-1) {           if (n.length()>MAX_namelength-1) {
1519              strncpy(c_names[i],n.c_str(),MAX_namelength-1);              strncpy(c_names[i],n.c_str(),MAX_namelength-1);
# Line 848  void MeshAdapter::saveVTK(const std::str Line 1523  void MeshAdapter::saveVTK(const std::str
1523           }           }
1524      }      }
1525      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1526      checkFinleyError();  
1527      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;      return;
1539  }  }
1540                                                                                                                                                                                                                                                                                                                                                    
# Line 889  SystemMatrixAdapter MeshAdapter::newSyst Line 1575  SystemMatrixAdapter MeshAdapter::newSyst
1575            
1576      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1577      checkFinleyError();      checkFinleyError();
1578      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);      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();      checkPasoError();
1589      Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);      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
# Line 915  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:
# Line 933  bool MeshAdapter::probeInterpolationOnDo Line 1662  bool MeshAdapter::probeInterpolationOnDo
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):             case(FaceElements):
1671               case(ReducedFaceElements):
1672             case(Points):             case(Points):
1673             case(ContactElementsZero):             case(ContactElementsZero):
1674               case(ReducedContactElementsZero):
1675             case(ContactElementsOne):             case(ContactElementsOne):
1676               case(ReducedContactElementsOne):
1677                 return true;                 return true;
1678             default:             default:
1679                 stringstream temp;                 stringstream temp;
# Line 947  bool MeshAdapter::probeInterpolationOnDo Line 1681  bool MeshAdapter::probeInterpolationOnDo
1681                 throw FinleyAdapterException(temp.str());                 throw FinleyAdapterException(temp.str());
1682          }          }
1683          break;          break;
1684         case(ReducedNodes):
1685            switch(functionSpaceType_target) {
1686               case(ReducedNodes):
1687               case(ReducedDegreesOfFreedom):
1688               case(Elements):
1689               case(ReducedElements):
1690               case(FaceElements):
1691               case(ReducedFaceElements):
1692               case(Points):
1693               case(ContactElementsZero):
1694               case(ReducedContactElementsZero):
1695               case(ContactElementsOne):
1696               case(ReducedContactElementsOne):
1697                   return true;
1698              case(Nodes):
1699              case(DegreesOfFreedom):
1700                 return false;
1701               default:
1702                   stringstream temp;
1703                   temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1704                   throw FinleyAdapterException(temp.str());
1705            }
1706            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 969  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 977  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               stringstream temp;               stringstream temp;
# Line 992  bool MeshAdapter::probeInterpolationOnDo Line 1779  bool MeshAdapter::probeInterpolationOnDo
1779       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1780         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
1781            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
1782              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):            case(Nodes):
1794            case(DegreesOfFreedom):            case(DegreesOfFreedom):
# Line 1059  escript::Data MeshAdapter::getSize() con Line 1851  escript::Data MeshAdapter::getSize() con
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  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1907  {  {
1908    int* tagList;    int out=0;
1909    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1910    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1911    return tagList[sampleNo];    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    void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1961    {
1962      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    
2013  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  void MeshAdapter::setTagMap(const std::string& name,  int tag)
2014  {  {
2015    int* referenceNoList;    Finley_Mesh* mesh=m_finleyMesh.get();
2016    int numReferenceNo;    Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2017    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    checkPasoError();
2018    return referenceNoList[sampleNo];    // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2019    }
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.532  
changed lines
  Added in v.1464

  ViewVC Help
Powered by ViewVC 1.1.26