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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26