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

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

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

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

Legend:
Removed from v.682  
changed lines
  Added in v.3317

  ViewVC Help
Powered by ViewVC 1.1.26