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

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

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

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

Legend:
Removed from v.798  
changed lines
  Added in v.3269

  ViewVC Help
Powered by ViewVC 1.1.26