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

Legend:
Removed from v.969  
changed lines
  Added in v.3344

  ViewVC Help
Powered by ViewVC 1.1.26