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

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

  ViewVC Help
Powered by ViewVC 1.1.26