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

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.1059  
changed lines
  Added in v.2842

  ViewVC Help
Powered by ViewVC 1.1.26