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

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

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

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

Legend:
Removed from v.757  
changed lines
  Added in v.3546

  ViewVC Help
Powered by ViewVC 1.1.26