/[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 793 by dhawcroft, Sat Jul 29 19:40:22 2006 UTC revision 3355 by caltinay, Tue Nov 16 06:35:06 2010 UTC
# Line 1  Line 1 
 // $Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2  #include "MeshAdapter.h"  /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15    #include "MeshAdapter.h"
16  #include "escript/Data.h"  #include "escript/Data.h"
17  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
18    #ifdef USE_NETCDF
19    #include <netcdfcpp.h>
20    #endif
21    extern "C" {
22    #include "esysUtils/blocktimer.h"
23    }
24    
25    #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::getDiracDeltaFunctionCode() const
651  {  {
652    return Points;     return Points;
653  }  }
654    
655  //  //
# Line 165  int MeshAdapter::getDiracDeltaFunctionCo Line 657  int MeshAdapter::getDiracDeltaFunctionCo
657  //  //
658  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
659  {  {
660    int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
661    checkFinleyError();     int numDim=Finley_Mesh_getDim(mesh);
662    return numDim;     checkFinleyError();
663       return numDim;
664    }
665    
666    //
667    // Return the number of data points summed across all MPI processes
668    //
669    int MeshAdapter::getNumDataPointsGlobal() const
670    {
671       return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
672  }  }
673    
674  //  //
# Line 180  pair<int,int> MeshAdapter::getDataShape( Line 681  pair<int,int> MeshAdapter::getDataShape(
681     int numSamples=0;     int numSamples=0;
682     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
683     switch (functionSpaceCode) {     switch (functionSpaceCode) {
684        case(Nodes):     case(Nodes):
685             numDataPointsPerSample=1;     numDataPointsPerSample=1;
686             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;     numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
687             break;     break;
688        case(Elements):     case(ReducedNodes):
689             if (mesh->Elements!=NULL) {     numDataPointsPerSample=1;
690               numSamples=mesh->Elements->numElements;     numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
691               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;     break;
692             }     case(Elements):
693             break;     if (mesh->Elements!=NULL) {
694        case(FaceElements):        numSamples=mesh->Elements->numElements;
695             if (mesh->FaceElements!=NULL) {        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
696                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;     }
697                  numSamples=mesh->FaceElements->numElements;     break;
698             }     case(ReducedElements):
699             break;     if (mesh->Elements!=NULL) {
700        case(Points):        numSamples=mesh->Elements->numElements;
701             if (mesh->Points!=NULL) {        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
702               numDataPointsPerSample=1;     }
703               numSamples=mesh->Points->numElements;     break;
704             }     case(FaceElements):
705             break;     if (mesh->FaceElements!=NULL) {
706        case(ContactElementsZero):        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
707             if (mesh->ContactElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
708               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
709               numSamples=mesh->ContactElements->numElements;     break;
710             }     case(ReducedFaceElements):
711             break;     if (mesh->FaceElements!=NULL) {
712        case(ContactElementsOne):        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
713             if (mesh->ContactElements!=NULL) {        numSamples=mesh->FaceElements->numElements;
714               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;     }
715               numSamples=mesh->ContactElements->numElements;     break;
716             }     case(Points):
717             break;     if (mesh->Points!=NULL) {
718        case(DegreesOfFreedom):        numDataPointsPerSample=1;
719             if (mesh->Nodes!=NULL) {        numSamples=mesh->Points->numElements;
720               numDataPointsPerSample=1;     }
721  #ifndef PASO_MPI     break;
722               numSamples=mesh->Nodes->numDegreesOfFreedom;     case(ContactElementsZero):
723  #else     if (mesh->ContactElements!=NULL) {
724               numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
725  #endif        numSamples=mesh->ContactElements->numElements;
726             }     }
727             break;     break;
728        case(ReducedDegreesOfFreedom):     case(ReducedContactElementsZero):
729             if (mesh->Nodes!=NULL) {     if (mesh->ContactElements!=NULL) {
730               numDataPointsPerSample=1;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
731  #ifndef PASO_MPI        numSamples=mesh->ContactElements->numElements;
732               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;     }
733  #else     break;
734               numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;     case(ContactElementsOne):
735  #endif     if (mesh->ContactElements!=NULL) {
736             }        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
737             break;        numSamples=mesh->ContactElements->numElements;
738        default:     }
739          stringstream temp;     break;
740          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();     case(ReducedContactElementsOne):
741          throw FinleyAdapterException(temp.str());     if (mesh->ContactElements!=NULL) {
742          break;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
743        }        numSamples=mesh->ContactElements->numElements;
744        return pair<int,int>(numDataPointsPerSample,numSamples);     }
745       break;
746       case(DegreesOfFreedom):
747       if (mesh->Nodes!=NULL) {
748          numDataPointsPerSample=1;
749          numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
750       }
751       break;
752       case(ReducedDegreesOfFreedom):
753       if (mesh->Nodes!=NULL) {
754          numDataPointsPerSample=1;
755          numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
756       }
757       break;
758       default:
759          stringstream temp;
760          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
761          throw FinleyAdapterException(temp.str());
762          break;
763       }
764       return pair<int,int>(numDataPointsPerSample,numSamples);
765  }  }
766    
767  //  //
768  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
769  //  //
770  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
771                       SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                       const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
773                       const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                       const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact) const
775  {  {
776    escriptDataC _rhs=rhs.getDataC();     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
777    escriptDataC _A  =A.getDataC();     if (smat==0)
778    escriptDataC _B=B.getDataC();     {
779    escriptDataC _C=C.getDataC();      throw FinleyAdapterException("finley only supports its own system matrices.");
780    escriptDataC _D=D.getDataC();     }
781    escriptDataC _X=X.getDataC();     escriptDataC _rhs=rhs.getDataC();
782    escriptDataC _Y=Y.getDataC();     escriptDataC _A  =A.getDataC();
783    escriptDataC _d=d.getDataC();     escriptDataC _B=B.getDataC();
784    escriptDataC _y=y.getDataC();     escriptDataC _C=C.getDataC();
785    escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _D=D.getDataC();
786    escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _X=X.getDataC();
787       escriptDataC _Y=Y.getDataC();
788       escriptDataC _d=d.getDataC();
789       escriptDataC _y=y.getDataC();
790       escriptDataC _d_contact=d_contact.getDataC();
791       escriptDataC _y_contact=y_contact.getDataC();
792    
793     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 );  
794    
795       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
796       checkFinleyError();
797    
798       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_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->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
802     checkFinleyError();     checkFinleyError();
803    }
804    
805    void  MeshAdapter::addPDEToLumpedSystem(
806                                            escript::Data& mat,
807                                            const escript::Data& D,
808                                            const escript::Data& d) const
809    {
810       escriptDataC _mat=mat.getDataC();
811       escriptDataC _D=D.getDataC();
812       escriptDataC _d=d.getDataC();
813    
814       Finley_Mesh* mesh=m_finleyMesh.get();
815    
816     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , &_d_contact, &_y_contact ,             Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
817     checkFinleyError();     checkFinleyError();
818      
819       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
820       checkFinleyError();
821    
822  }  }
823    
824    
825  //  //
826  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
827  //  //
# Line 283  void MeshAdapter::addPDEToRHS( escript:: Line 829  void MeshAdapter::addPDEToRHS( escript::
829  {  {
830     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
831    
832     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
833     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     escriptDataC _X=X.getDataC();
834       escriptDataC _Y=Y.getDataC();
835       escriptDataC _y=y.getDataC();
836       escriptDataC _y_contact=y_contact.getDataC();
837    
838       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
839       checkFinleyError();
840    
841       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
842       checkFinleyError();
843    
844       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
845       checkFinleyError();
846    }
847    //
848    // adds PDE of second order into a transport problem
849    //
850    void MeshAdapter::addPDEToTransportProblem(
851                                               AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
852                                               const escript::Data& A, const escript::Data& B, const escript::Data& C,
853                                               const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
854                                               const escript::Data& d, const escript::Data& y,
855                                               const escript::Data& d_contact,const escript::Data& y_contact) const
856    {
857       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
858       if (tpa==0)
859       {
860        throw FinleyAdapterException("finley only supports its own transport problems.");
861       }
862    
863    
864       DataTypes::ShapeType shape;
865       source.expand();
866       escriptDataC _source=source.getDataC();
867       escriptDataC _M=M.getDataC();
868       escriptDataC _A=A.getDataC();
869       escriptDataC _B=B.getDataC();
870       escriptDataC _C=C.getDataC();
871       escriptDataC _D=D.getDataC();
872       escriptDataC _X=X.getDataC();
873       escriptDataC _Y=Y.getDataC();
874       escriptDataC _d=d.getDataC();
875       escriptDataC _y=y.getDataC();
876       escriptDataC _d_contact=d_contact.getDataC();
877       escriptDataC _y_contact=y_contact.getDataC();
878    
879       Finley_Mesh* mesh=m_finleyMesh.get();
880       Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
881    
882       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
883     checkFinleyError();     checkFinleyError();
884    
885     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
886     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     checkFinleyError();
887    
888       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
889     checkFinleyError();     checkFinleyError();
890     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
891     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
892     checkFinleyError();     checkFinleyError();
893  }  }
894    
# Line 301  void MeshAdapter::addPDEToRHS( escript:: Line 897  void MeshAdapter::addPDEToRHS( escript::
897  //  //
898  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
899  {  {
900    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
901    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
902    if (inDomain!=*this)       if (inDomain!=*this)  
903      throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
904    if (targetDomain!=*this)     if (targetDomain!=*this)
905      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");        throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
906    
907    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
908    switch(in.getFunctionSpace().getTypeCode()) {     escriptDataC _target=target.getDataC();
909       case(Nodes):     escriptDataC _in=in.getDataC();
910          switch(target.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
911             case(Nodes):     case(Nodes):
912             case(ReducedDegreesOfFreedom):        switch(target.getFunctionSpace().getTypeCode()) {
913             case(DegreesOfFreedom):        case(Nodes):
914                 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));        case(ReducedNodes):
915                 break;        case(DegreesOfFreedom):
916             case(Elements):        case(ReducedDegreesOfFreedom):
917                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
918                 break;        break;
919             case(FaceElements):        case(Elements):
920                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));        case(ReducedElements):
921                 break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
922             case(Points):        break;
923                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));        case(FaceElements):
924                 break;        case(ReducedFaceElements):
925             case(ContactElementsZero):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
926                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        break;
927                 break;        case(Points):
928             case(ContactElementsOne):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
929                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        break;
930                 break;        case(ContactElementsZero):
931             default:        case(ReducedContactElementsZero):
932                 stringstream temp;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
933                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        break;
934                 throw FinleyAdapterException(temp.str());        case(ContactElementsOne):
935                 break;        case(ReducedContactElementsOne):
936          }        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
937          break;        break;
938       case(Elements):        default:
939          if (target.getFunctionSpace().getTypeCode()==Elements) {           stringstream temp;
940             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();
941          } else {           throw FinleyAdapterException(temp.str());
942             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");           break;
943          }        }
944          break;        break;
945       case(FaceElements):     case(ReducedNodes):
946          if (target.getFunctionSpace().getTypeCode()==FaceElements) {        switch(target.getFunctionSpace().getTypeCode()) {
947             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));        case(Nodes):
948          } else {        case(ReducedNodes):
949             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");        case(DegreesOfFreedom):
950             break;        case(ReducedDegreesOfFreedom):
951         }        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
952       case(Points):        break;
953          if (target.getFunctionSpace().getTypeCode()==Points) {        case(Elements):
954             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));        case(ReducedElements):
955          } else {        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
956             throw FinleyAdapterException("Error - No interpolation with data on points possible.");        break;
957             break;        case(FaceElements):
958          }        case(ReducedFaceElements):
959          break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
960       case(ContactElementsZero):        break;
961       case(ContactElementsOne):        case(Points):
962          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
963             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));        break;
964          } else {        case(ContactElementsZero):
965             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");        case(ReducedContactElementsZero):
966             break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
967          }        break;
968          break;        case(ContactElementsOne):
969       case(DegreesOfFreedom):              case(ReducedContactElementsOne):
970          switch(target.getFunctionSpace().getTypeCode()) {        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
971             case(ReducedDegreesOfFreedom):        break;
972             case(DegreesOfFreedom):        default:
973             case(Nodes):           stringstream temp;
974                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();
975                break;           throw FinleyAdapterException(temp.str());
976  #ifndef PASO_MPI           break;
977             case(Elements):        }
978                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));        break;
979                break;     case(Elements):
980             case(FaceElements):        if (target.getFunctionSpace().getTypeCode()==Elements) {
981                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
982                break;        } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
983             case(Points):           Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
984                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));        } else {
985                break;           throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
986             case(ContactElementsZero):        }
987             case(ContactElementsOne):        break;
988                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));     case(ReducedElements):
989               break;        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
990  #else           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
991             /* need to copy Degrees of freedom data to nodal data so that the external values are available */        } else {
992             case(Elements):           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
993             {        }
994                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );        break;
995                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));     case(FaceElements):
996                break;        if (target.getFunctionSpace().getTypeCode()==FaceElements) {
997             }           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
998             case(FaceElements):        } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
999             {           Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
1000                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );        } else {
1001                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));           throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
1002                break;        }
1003             }        break;
1004             case(Points):     case(ReducedFaceElements):
1005             {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
1006                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
1007                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));        } else {
1008                break;           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
1009             }        }
1010             case(ContactElementsZero):        break;
1011             case(ContactElementsOne):     case(Points):
1012             {        if (target.getFunctionSpace().getTypeCode()==Points) {
1013                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1014                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));        } else {
1015               break;           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1016             }        }
1017  #endif        break;
1018             default:     case(ContactElementsZero):
1019               stringstream temp;     case(ContactElementsOne):
1020               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) {
1021               throw FinleyAdapterException(temp.str());           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1022               break;        } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1023          }           Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1024          break;        } else {
1025       case(ReducedDegreesOfFreedom):           throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1026         switch(target.getFunctionSpace().getTypeCode()) {        }
1027            case(ReducedDegreesOfFreedom):        break;
1028               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));     case(ReducedContactElementsZero):
1029               break;     case(ReducedContactElementsOne):
1030            case(Elements):        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1031               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1032               break;        } else {
1033            case(FaceElements):           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1034               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));        }
1035               break;        break;
1036            case(Points):     case(DegreesOfFreedom):      
1037               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));        switch(target.getFunctionSpace().getTypeCode()) {
1038               break;        case(ReducedDegreesOfFreedom):
1039            case(ContactElementsZero):        case(DegreesOfFreedom):
1040            case(ContactElementsOne):        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1041               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));        break;
1042               break;    
1043            case(Nodes):        case(Nodes):
1044               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");        case(ReducedNodes):
1045               break;        if (getMPISize()>1) {
1046            case(DegreesOfFreedom):           escript::Data temp=escript::Data(in);
1047               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");           temp.expand();
1048               break;           escriptDataC _in2 = temp.getDataC();
1049            default:           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1050               stringstream temp;        } else {
1051               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();           Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1052               throw FinleyAdapterException(temp.str());        }
1053               break;        break;
1054         }        case(Elements):
1055         break;        case(ReducedElements):
1056       default:        if (getMPISize()>1) {
1057          stringstream temp;           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1058          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();           escriptDataC _in2 = temp.getDataC();
1059          throw FinleyAdapterException(temp.str());           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1060          break;        } else {
1061    }           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1062    checkFinleyError();        }
1063          break;
1064          case(FaceElements):
1065          case(ReducedFaceElements):
1066          if (getMPISize()>1) {
1067             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1068             escriptDataC _in2 = temp.getDataC();
1069             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1070      
1071          } else {
1072             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1073          }
1074          break;
1075          case(Points):
1076          if (getMPISize()>1) {
1077             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1078             escriptDataC _in2 = temp.getDataC();
1079          } else {
1080             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1081          }
1082          break;
1083          case(ContactElementsZero):
1084          case(ContactElementsOne):
1085          case(ReducedContactElementsZero):
1086          case(ReducedContactElementsOne):
1087          if (getMPISize()>1) {
1088             escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1089             escriptDataC _in2 = temp.getDataC();
1090             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1091          } else {
1092             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1093          }
1094          break;
1095          default:
1096             stringstream temp;
1097             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1098             throw FinleyAdapterException(temp.str());
1099             break;
1100          }
1101          break;
1102       case(ReducedDegreesOfFreedom):
1103          switch(target.getFunctionSpace().getTypeCode()) {
1104          case(Nodes):
1105          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1106          break;
1107          case(ReducedNodes):
1108          if (getMPISize()>1) {
1109             escript::Data temp=escript::Data(in);
1110             temp.expand();
1111             escriptDataC _in2 = temp.getDataC();
1112             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1113          } else {
1114             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1115          }
1116          break;
1117          case(DegreesOfFreedom):
1118          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1119          break;
1120          case(ReducedDegreesOfFreedom):
1121          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1122          break;
1123          case(Elements):
1124          case(ReducedElements):
1125          if (getMPISize()>1) {
1126             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1127             escriptDataC _in2 = temp.getDataC();
1128             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1129          } else {
1130             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1131          }
1132          break;
1133          case(FaceElements):
1134          case(ReducedFaceElements):
1135          if (getMPISize()>1) {
1136             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1137             escriptDataC _in2 = temp.getDataC();
1138             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1139          } else {
1140             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1141          }
1142          break;
1143          case(Points):
1144          if (getMPISize()>1) {
1145             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1146             escriptDataC _in2 = temp.getDataC();
1147             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1148          } else {
1149             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1150          }
1151          break;
1152          case(ContactElementsZero):
1153          case(ContactElementsOne):
1154          case(ReducedContactElementsZero):
1155          case(ReducedContactElementsOne):
1156          if (getMPISize()>1) {
1157             escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1158             escriptDataC _in2 = temp.getDataC();
1159             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1160          } else {
1161             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1162          }
1163          break;
1164          default:
1165             stringstream temp;
1166             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1167             throw FinleyAdapterException(temp.str());
1168             break;
1169          }
1170          break;
1171       default:
1172          stringstream temp;
1173          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1174          throw FinleyAdapterException(temp.str());
1175          break;
1176       }
1177       checkFinleyError();
1178  }  }
1179    
1180  //  //
# Line 471  void MeshAdapter::interpolateOnDomain(es Line 1182  void MeshAdapter::interpolateOnDomain(es
1182  //  //
1183  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1184  {  {
1185    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1186    if (argDomain!=*this)     if (argDomain!=*this)
1187       throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1188    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1189    // in case of values node coordinates we can do the job directly:     // in case of values node coordinates we can do the job directly:
1190    if (arg.getFunctionSpace().getTypeCode()==Nodes) {     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1191       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));        escriptDataC _arg=arg.getDataC();
1192    } else {        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1193       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);     } else {
1194       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1195       // this is then interpolated onto arg:        escriptDataC _tmp_data=tmp_data.getDataC();
1196       interpolateOnDomain(arg,tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1197    }        // this is then interpolated onto arg:
1198    checkFinleyError();        interpolateOnDomain(arg,tmp_data);
1199       }
1200       checkFinleyError();
1201  }  }
1202    
1203  //  //
# Line 492  void MeshAdapter::setToX(escript::Data& Line 1205  void MeshAdapter::setToX(escript::Data&
1205  //  //
1206  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1207  {  {
1208    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1209    if (normalDomain!=*this)     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1210       throw FinleyAdapterException("Error - Illegal domain of normal locations");     if (normalDomain!=*this)
1211    Finley_Mesh* mesh=m_finleyMesh.get();        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1212    switch(normal.getFunctionSpace().getTypeCode()) {     Finley_Mesh* mesh=m_finleyMesh.get();
1213      case(Nodes):     escriptDataC _normal=normal.getDataC();
1214        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");     switch(normal.getFunctionSpace().getTypeCode()) {
1215        break;     case(Nodes):
1216      case(Elements):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1217        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");     break;
1218        break;     case(ReducedNodes):
1219      case (FaceElements):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1220        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));     break;
1221        break;     case(Elements):
1222      case(Points):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1223        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");     break;
1224        break;     case(ReducedElements):
1225      case (ContactElementsOne):     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1226      case (ContactElementsZero):     break;
1227        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));     case (FaceElements):
1228        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1229      case(DegreesOfFreedom):     break;
1230        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");     case (ReducedFaceElements):
1231        break;     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1232      case(ReducedDegreesOfFreedom):     break;
1233        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");     case(Points):
1234        break;     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1235      default:     break;
1236       case (ContactElementsOne):
1237       case (ContactElementsZero):
1238       Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1239       break;
1240       case (ReducedContactElementsOne):
1241       case (ReducedContactElementsZero):
1242       Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1243       break;
1244       case(DegreesOfFreedom):
1245       throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1246       break;
1247       case(ReducedDegreesOfFreedom):
1248       throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1249       break;
1250       default:
1251        stringstream temp;        stringstream temp;
1252        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();
1253        throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
1254        break;        break;
1255    }     }
1256    checkFinleyError();     checkFinleyError();
1257  }  }
1258    
1259  //  //
# Line 533  void MeshAdapter::setToNormal(escript::D Line 1261  void MeshAdapter::setToNormal(escript::D
1261  //  //
1262  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1263  {  {
1264    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1265    if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1266       throw FinleyAdapterException("Error - Illegal domain of interpolation target");     if (targetDomain!=this)
1267          throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1268    
1269    throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1270  }  }
1271    
1272  //  //
1273  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
1274  //  //
1275  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1276  {  {
1277    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1278    if (argDomain!=*this)     if (argDomain!=*this)
1279       throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1280    
1281    Finley_Mesh* mesh=m_finleyMesh.get();     double blocktimer_start = blocktimer_time();
1282    switch(arg.getFunctionSpace().getTypeCode()) {     Finley_Mesh* mesh=m_finleyMesh.get();
1283       case(Nodes):     escriptDataC _temp;
1284          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");     escript::Data temp;
1285          break;     escriptDataC _arg=arg.getDataC();
1286       case(Elements):     switch(arg.getFunctionSpace().getTypeCode()) {
1287          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);     case(Nodes):
1288          break;     temp=escript::Data( arg, escript::function(*this) );
1289       case(FaceElements):     _temp=temp.getDataC();
1290          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1291          break;     break;
1292       case(Points):     case(ReducedNodes):
1293          throw FinleyAdapterException("Error - Integral of data on points is not supported.");     temp=escript::Data( arg, escript::function(*this) );
1294          break;     _temp=temp.getDataC();
1295       case(ContactElementsZero):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1296          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);     break;
1297          break;     case(Elements):
1298       case(ContactElementsOne):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1299          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);     break;
1300          break;     case(ReducedElements):
1301       case(DegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1302          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");     break;
1303          break;     case(FaceElements):
1304       case(ReducedDegreesOfFreedom):     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1305          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");     break;
1306          break;     case(ReducedFaceElements):
1307       default:     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1308          stringstream temp;     break;
1309          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     case(Points):
1310          throw FinleyAdapterException(temp.str());     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1311          break;     break;
1312    }     case(ContactElementsZero):
1313    checkFinleyError();     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1314       break;
1315       case(ReducedContactElementsZero):
1316       Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1317       break;
1318       case(ContactElementsOne):
1319       Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1320       break;
1321       case(ReducedContactElementsOne):
1322       Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1323       break;
1324       case(DegreesOfFreedom):
1325       temp=escript::Data( arg, escript::function(*this) );
1326       _temp=temp.getDataC();
1327       Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1328       break;
1329       case(ReducedDegreesOfFreedom):
1330       temp=escript::Data( arg, escript::function(*this) );
1331       _temp=temp.getDataC();
1332       Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1333       break;
1334       default:
1335          stringstream temp;
1336          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1337          throw FinleyAdapterException(temp.str());
1338          break;
1339       }
1340       checkFinleyError();
1341       blocktimer_increment("integrate()", blocktimer_start);
1342  }  }
1343    
1344  //  //
# Line 589  void MeshAdapter::setToIntegrals(std::ve Line 1346  void MeshAdapter::setToIntegrals(std::ve
1346  //  //
1347  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1348  {  {
1349    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1350    if (argDomain!=*this)     if (argDomain!=*this)
1351      throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1352    const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1353    if (gradDomain!=*this)     if (gradDomain!=*this)
1354       throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1355    
1356    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1357    escriptDataC nodeDataC;     escriptDataC _grad=grad.getDataC();
1358  #ifdef PASO_MPI     escriptDataC nodeDataC;
1359    escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );     escript::Data temp;
1360    if( arg.getFunctionSpace().getTypeCode() != Nodes )     if (getMPISize()>1) {
1361    {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1362      Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));           temp=escript::Data( arg,  continuousFunction(*this) );
1363      nodeDataC = nodeTemp.getDataC();           nodeDataC = temp.getDataC();
1364    }        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1365    else           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1366      nodeDataC = arg.getDataC();           nodeDataC = temp.getDataC();
1367  #else        } else {
1368    nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
1369  #endif        }
1370    switch(grad.getFunctionSpace().getTypeCode()) {     } else {
1371         case(Nodes):        nodeDataC = arg.getDataC();
1372            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");     }
1373            break;     switch(grad.getFunctionSpace().getTypeCode()) {
1374         case(Elements):     case(Nodes):
1375            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1376            break;     break;
1377         case(FaceElements):     case(ReducedNodes):
1378            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1379            break;     break;
1380         case(Points):     case(Elements):
1381            throw FinleyAdapterException("Error - Gradient at points is not supported.");     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1382            break;     break;
1383         case(ContactElementsZero):     case(ReducedElements):
1384            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1385            break;     break;
1386         case(ContactElementsOne):     case(FaceElements):
1387            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1388            break;     break;
1389         case(DegreesOfFreedom):     case(ReducedFaceElements):
1390            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1391            break;     break;
1392         case(ReducedDegreesOfFreedom):     case(Points):
1393            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1394            break;     break;
1395         default:     case(ContactElementsZero):
1396            stringstream temp;     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1397            temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();     break;
1398            throw FinleyAdapterException(temp.str());     case(ReducedContactElementsZero):
1399            break;     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1400    }     break;
1401    checkFinleyError();     case(ContactElementsOne):
1402       Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1403       break;
1404       case(ReducedContactElementsOne):
1405       Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1406       break;
1407       case(DegreesOfFreedom):
1408       throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1409       break;
1410       case(ReducedDegreesOfFreedom):
1411       throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1412       break;
1413       default:
1414          stringstream temp;
1415          temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1416          throw FinleyAdapterException(temp.str());
1417          break;
1418       }
1419       checkFinleyError();
1420  }  }
1421    
1422  //  //
# Line 649  void MeshAdapter::setToGradient(escript: Line 1424  void MeshAdapter::setToGradient(escript:
1424  //  //
1425  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1426  {  {
1427    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1428    escriptDataC tmp=size.getDataC();     escriptDataC tmp=size.getDataC();
1429    switch(size.getFunctionSpace().getTypeCode()) {     switch(size.getFunctionSpace().getTypeCode()) {
1430         case(Nodes):     case(Nodes):
1431            throw FinleyAdapterException("Error - Size of nodes is not supported.");     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1432            break;     break;
1433         case(Elements):     case(ReducedNodes):
1434            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1435            break;     break;
1436         case(FaceElements):     case(Elements):
1437            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1438            break;     break;
1439         case(Points):     case(ReducedElements):
1440            throw FinleyAdapterException("Error - Size of point elements is not supported.");     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1441            break;     break;
1442         case(ContactElementsZero):     case(FaceElements):
1443         case(ContactElementsOne):     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1444            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);     break;
1445            break;     case(ReducedFaceElements):
1446         case(DegreesOfFreedom):     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1447            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");     break;
1448            break;     case(Points):
1449         case(ReducedDegreesOfFreedom):     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1450            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");     break;
1451            break;     case(ContactElementsZero):
1452         default:     case(ContactElementsOne):
1453            stringstream temp;     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1454            temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();     break;
1455            throw FinleyAdapterException(temp.str());     case(ReducedContactElementsZero):
1456            break;     case(ReducedContactElementsOne):
1457    }     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1458    checkFinleyError();     break;
1459       case(DegreesOfFreedom):
1460       throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1461       break;
1462       case(ReducedDegreesOfFreedom):
1463       throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1464       break;
1465       default:
1466          stringstream temp;
1467          temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1468          throw FinleyAdapterException(temp.str());
1469          break;
1470       }
1471       checkFinleyError();
1472  }  }
1473    
1474  // sets the location of nodes:  //
1475    // sets the location of nodes
1476    //
1477  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1478  {  {
1479    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1480    escriptDataC tmp;     escriptDataC tmp;
1481    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1482    if (newDomain!=*this)     if (newDomain!=*this)
1483       throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1484    tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1485    Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1486    checkFinleyError();         Finley_Mesh_setCoordinates(mesh,&tmp);
1487  }     } else {
1488           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1489  // saves a data array in openDX format:         tmp = new_x_inter.getDataC();
1490  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const         Finley_Mesh_setCoordinates(mesh,&tmp);
1491  {     }
1492      int MAX_namelength=256;     checkFinleyError();
1493      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;  
   }  
1494    
1495    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  //
1496    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  // Helper for the save* methods. Extracts optional data variable names and the
1497    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  // corresponding pointers from python dictionary. Caller must free arrays.
1498    //
1499      boost::python::list keys=arg.keys();  void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1500      for (int i=0;i<num_data;++i) {  {
1501           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);     numData = boost::python::extract<int>(arg.attr("__len__")());
1502           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)     /* win32 refactor */
1503               throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1504           data[i]=d.getDataC();     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1505           ptr_data[i]=&(data[i]);     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1506           std::string n=boost::python::extract<std::string>(keys[i]);  
1507           c_names[i]=&(names[i][0]);     boost::python::list keys=arg.keys();
1508           if (n.length()>MAX_namelength-1) {     for (int i=0; i<numData; ++i) {
1509              strncpy(c_names[i],n.c_str(),MAX_namelength-1);        string n=boost::python::extract<string>(keys[i]);
1510              c_names[i][MAX_namelength-1]='\0';        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1511           } else {        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1512              strcpy(c_names[i],n.c_str());           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1513           }        data[i] = d.getDataC();
1514      }        dataPtr[i] = &(data[i]);
1515      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);        names[i] = TMPMEMALLOC(n.length()+1, char);
1516      checkFinleyError();        strcpy(names[i], n.c_str());
1517           }
1518        /* win32 refactor */  }
1519    TMPMEMFREE(c_names);  
1520    TMPMEMFREE(data);  //
1521    TMPMEMFREE(ptr_data);  // saves mesh and optionally data arrays in openDX format
1522    for(int i=0;i<num_data;i++)  //
1523    {  void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1524      TMPMEMFREE(names[i]);  {
1525    }     int num_data;
1526    TMPMEMFREE(names);     char **names;
1527       escriptDataC *data;
1528       escriptDataC **ptr_data;
1529    
1530       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1531       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1532       checkFinleyError();
1533    
1534       /* win32 refactor */
1535       TMPMEMFREE(data);
1536       TMPMEMFREE(ptr_data);
1537       for(int i=0; i<num_data; i++)
1538       {
1539          TMPMEMFREE(names[i]);
1540       }
1541       TMPMEMFREE(names);
1542    
1543      return;     return;
1544  }  }
1545    
1546  // saves a data array in openVTK format:  //
1547  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1548    //
1549    void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1550  {  {
1551      int MAX_namelength=256;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1552      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1553    /* win32 refactor */                metadata, metadata_schema, arg);
1554    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;  
   }  
1555    
1556    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1557    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  {
1558    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  #ifdef ESYS_MPI
1559        index_t myFirstNode=0, myLastNode=0, k=0;
1560      boost::python::list keys=arg.keys();      index_t* globalNodeIndex=0;
1561      for (int i=0;i<num_data;++i) {      Finley_Mesh* mesh_p=m_finleyMesh.get();
1562           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);      if (fs_code == FINLEY_REDUCED_NODES)
1563           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)      {
1564               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1565           data[i]=d.getDataC();      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1566           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());  
          }  
1567      }      }
1568  #ifndef PASO_MPI          else
1569      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      {
1570  #else      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1571      Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1572        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1573        }
1574        k=globalNodeIndex[id];
1575        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1576  #endif  #endif
1577        return true;
1578    }
1579    
 checkFinleyError();  
   /* win32 refactor */  
   TMPMEMFREE(c_names);  
   TMPMEMFREE(data);  
   TMPMEMFREE(ptr_data);  
   for(int i=0;i<num_data;i++)  
   {  
     TMPMEMFREE(names[i]);  
   }  
   TMPMEMFREE(names);  
1580    
1581      return;  
1582    //
1583    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1584    //
1585    ASM_ptr MeshAdapter::newSystemMatrix(
1586                                                     const int row_blocksize,
1587                                                     const escript::FunctionSpace& row_functionspace,
1588                                                     const int column_blocksize,
1589                                                     const escript::FunctionSpace& column_functionspace,
1590                                                     const int type) const
1591    {
1592       int reduceRowOrder=0;
1593       int reduceColOrder=0;
1594       // is the domain right?
1595       const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1596       if (row_domain!=*this)
1597          throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1598       const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1599       if (col_domain!=*this)
1600          throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1601       // is the function space type right
1602       if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1603          reduceRowOrder=0;
1604       } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1605          reduceRowOrder=1;
1606       } else {
1607          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1608       }
1609       if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1610          reduceColOrder=0;
1611       } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1612          reduceColOrder=1;
1613       } else {
1614          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1615       }
1616       // generate matrix:
1617    
1618       Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1619       checkFinleyError();
1620       Paso_SystemMatrix* fsystemMatrix;
1621       int trilinos = 0;
1622       if (trilinos) {
1623    #ifdef TRILINOS
1624          /* Allocation Epetra_VrbMatrix here */
1625    #endif
1626       }
1627       else {
1628          fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1629       }
1630       checkPasoError();
1631       Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1632       SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1633       return ASM_ptr(sma);
1634    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1635  }  }
1636                                                                                                                                                                            
1637                                                                                                                                                                            //
1638  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a TransportProblemAdapter
1639  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  //
1640                        const int row_blocksize,  ATP_ptr MeshAdapter::newTransportProblem(
1641                        const escript::FunctionSpace& row_functionspace,                                                           const bool useBackwardEuler,
1642                        const int column_blocksize,                                                           const int blocksize,
1643                        const escript::FunctionSpace& column_functionspace,                                                           const escript::FunctionSpace& functionspace,
1644                        const int type) const                                                           const int type) const
1645  {  {
1646      int reduceRowOrder=0;     int reduceOrder=0;
1647      int reduceColOrder=0;     // is the domain right?
1648      // is the domain right?     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1649      const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     if (domain!=*this)
1650      if (row_domain!=*this)        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1651            throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");     // is the function space type right
1652      const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1653      if (col_domain!=*this)        reduceOrder=0;
1654            throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1655      // is the function space type right        reduceOrder=1;
1656      if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     } else {
1657          reduceRowOrder=0;        throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1658      } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     }
1659          reduceRowOrder=1;     // generate matrix:
1660      } else {  
1661          throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1662      }     checkFinleyError();
1663      if (column_functionspace.getTypeCode()==DegreesOfFreedom) {     Paso_TransportProblem* transportProblem;
1664          reduceColOrder=0;     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1665      } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {     checkPasoError();
1666          reduceColOrder=1;     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1667      } else {     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1668          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");     return ATP_ptr(tpa);
1669      }  //   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);  
1670  }  }
1671    
1672  //  //
# Line 848  SystemMatrixAdapter MeshAdapter::newSyst Line 1678  SystemMatrixAdapter MeshAdapter::newSyst
1678  // 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:
1679  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1680  {  {
1681    switch(functionSpaceCode) {     switch(functionSpaceCode) {
1682         case(Nodes):     case(Nodes):
1683         case(DegreesOfFreedom):     case(DegreesOfFreedom):
1684         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1685            return false;     return false;
1686            break;     break;
1687         case(Elements):     case(Elements):
1688         case(FaceElements):     case(FaceElements):
1689         case(Points):     case(Points):
1690         case(ContactElementsZero):     case(ContactElementsZero):
1691         case(ContactElementsOne):     case(ContactElementsOne):
1692            return true;     case(ReducedElements):
1693            break;     case(ReducedFaceElements):
1694         default:     case(ReducedContactElementsZero):
1695            stringstream temp;     case(ReducedContactElementsOne):
1696            temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;     return true;
1697            throw FinleyAdapterException(temp.str());     break;
1698            break;     default:
1699    }        stringstream temp;
1700    checkFinleyError();        temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1701    return false;        throw FinleyAdapterException(temp.str());
1702          break;
1703       }
1704       checkFinleyError();
1705       return false;
1706    }
1707    
1708    bool
1709    MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1710    {
1711       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1712        class 1: DOF <-> Nodes
1713        class 2: ReducedDOF <-> ReducedNodes
1714        class 3: Points
1715        class 4: Elements
1716        class 5: ReducedElements
1717        class 6: FaceElements
1718        class 7: ReducedFaceElements
1719        class 8: ContactElementZero <-> ContactElementOne
1720        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1721    
1722       There is also a set of lines. Interpolation is possible down a line but not between lines.
1723       class 1 and 2 belong to all lines so aren't considered.
1724        line 0: class 3
1725        line 1: class 4,5
1726        line 2: class 6,7
1727        line 3: class 8,9
1728    
1729       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1730       eg hasnodes is true if we have at least one instance of Nodes.
1731       */
1732        if (fs.empty())
1733        {
1734            return false;
1735        }
1736        vector<int> hasclass(10);
1737        vector<int> hasline(4);
1738        bool hasnodes=false;
1739        bool hasrednodes=false;
1740        bool hascez=false;
1741        bool hasrcez=false;
1742        for (int i=0;i<fs.size();++i)
1743        {
1744        switch(fs[i])
1745        {
1746        case(Nodes):   hasnodes=true;   // no break is deliberate
1747        case(DegreesOfFreedom):
1748            hasclass[1]=1;
1749            break;
1750        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1751        case(ReducedDegreesOfFreedom):
1752            hasclass[2]=1;
1753            break;
1754        case(Points):
1755            hasline[0]=1;
1756            hasclass[3]=1;
1757            break;
1758        case(Elements):
1759            hasclass[4]=1;
1760            hasline[1]=1;
1761            break;
1762        case(ReducedElements):
1763            hasclass[5]=1;
1764            hasline[1]=1;
1765            break;
1766        case(FaceElements):
1767            hasclass[6]=1;
1768            hasline[2]=1;
1769            break;
1770        case(ReducedFaceElements):
1771            hasclass[7]=1;
1772            hasline[2]=1;
1773            break;
1774        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1775        case(ContactElementsOne):
1776            hasclass[8]=1;
1777            hasline[3]=1;
1778            break;
1779        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1780        case(ReducedContactElementsOne):
1781            hasclass[9]=1;
1782            hasline[3]=1;
1783            break;
1784        default:
1785            return false;
1786        }
1787        }
1788        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1789        // fail if we have more than one leaf group
1790    
1791        if (totlines>1)
1792        {
1793        return false;   // there are at least two branches we can't interpolate between
1794        }
1795        else if (totlines==1)
1796        {
1797        if (hasline[0]==1)      // we have points
1798        {
1799            resultcode=Points;
1800        }
1801        else if (hasline[1]==1)
1802        {
1803            if (hasclass[5]==1)
1804            {
1805            resultcode=ReducedElements;
1806            }
1807            else
1808            {
1809            resultcode=Elements;
1810            }
1811        }
1812        else if (hasline[2]==1)
1813        {
1814            if (hasclass[7]==1)
1815            {
1816            resultcode=ReducedFaceElements;
1817            }
1818            else
1819            {
1820            resultcode=FaceElements;
1821            }
1822        }
1823        else    // so we must be in line3
1824        {
1825            if (hasclass[9]==1)
1826            {
1827            // need something from class 9
1828            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1829            }
1830            else
1831            {
1832            // something from class 8
1833            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1834            }
1835        }
1836        }
1837        else    // totlines==0
1838        {
1839        if (hasclass[2]==1)
1840        {
1841            // something from class 2
1842            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1843        }
1844        else
1845        {   // something from class 1
1846            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1847        }
1848        }
1849        return true;
1850  }  }
1851    
1852  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1853  {  {
1854    switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1855       case(Nodes):     case(Nodes):
1856          switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1857             case(Nodes):      case(Nodes):
1858             case(ReducedDegreesOfFreedom):      case(ReducedNodes):
1859             case(DegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1860             case(Elements):      case(DegreesOfFreedom):
1861             case(FaceElements):      case(Elements):
1862             case(Points):      case(ReducedElements):
1863             case(ContactElementsZero):      case(FaceElements):
1864             case(ContactElementsOne):      case(ReducedFaceElements):
1865                 return true;      case(Points):
1866             default:      case(ContactElementsZero):
1867                 stringstream temp;      case(ReducedContactElementsZero):
1868                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;      case(ContactElementsOne):
1869                 throw FinleyAdapterException(temp.str());      case(ReducedContactElementsOne):
1870          }      return true;
1871          break;      default:
1872       case(Elements):            stringstream temp;
1873          if (functionSpaceType_target==Elements) {            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1874             return true;            throw FinleyAdapterException(temp.str());
1875          } else {     }
1876             return false;     break;
1877          }     case(ReducedNodes):
1878       case(FaceElements):      switch(functionSpaceType_target) {
1879          if (functionSpaceType_target==FaceElements) {      case(ReducedNodes):
1880             return true;      case(ReducedDegreesOfFreedom):
1881          } else {      case(Elements):
1882             return false;      case(ReducedElements):
1883          }      case(FaceElements):
1884       case(Points):      case(ReducedFaceElements):
1885          if (functionSpaceType_target==Points) {      case(Points):
1886             return true;      case(ContactElementsZero):
1887          } else {      case(ReducedContactElementsZero):
1888             return false;      case(ContactElementsOne):
1889          }      case(ReducedContactElementsOne):
1890       case(ContactElementsZero):      return true;
1891       case(ContactElementsOne):      case(Nodes):
1892          if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      case(DegreesOfFreedom):
1893             return true;      return false;
1894        default:
1895            stringstream temp;
1896            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1897            throw FinleyAdapterException(temp.str());
1898       }
1899       break;
1900       case(Elements):
1901        if (functionSpaceType_target==Elements) {
1902          return true;
1903        } else if (functionSpaceType_target==ReducedElements) {
1904          return true;
1905          } else {          } else {
1906             return false;            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());  
1907          }          }
1908          break;     case(ReducedElements):
1909       case(ReducedDegreesOfFreedom):      if (functionSpaceType_target==ReducedElements) {
1910         switch(functionSpaceType_target) {        return true;
1911            case(ReducedDegreesOfFreedom):      } else {
1912            case(Elements):            return false;
1913            case(FaceElements):      }
1914            case(Points):     case(FaceElements):
1915            case(ContactElementsZero):      if (functionSpaceType_target==FaceElements) {
1916            case(ContactElementsOne):              return true;
1917                return true;      } else if (functionSpaceType_target==ReducedFaceElements) {
1918            case(Nodes):              return true;
1919            case(DegreesOfFreedom):      } else {
1920               return false;              return false;
1921            default:      }
1922               stringstream temp;     case(ReducedFaceElements):
1923               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;      if (functionSpaceType_target==ReducedFaceElements) {
1924               throw FinleyAdapterException(temp.str());              return true;
1925         }      } else {
1926         break;          return false;
1927       default:      }
1928          stringstream temp;     case(Points):
1929          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;      if (functionSpaceType_target==Points) {
1930          throw FinleyAdapterException(temp.str());              return true;
1931          break;      } else {
1932    }              return false;
1933    checkFinleyError();      }
1934    return false;     case(ContactElementsZero):
1935       case(ContactElementsOne):
1936        if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1937                return true;
1938        } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1939                return true;
1940        } else {
1941                return false;
1942        }
1943       case(ReducedContactElementsZero):
1944       case(ReducedContactElementsOne):
1945        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1946                return true;
1947        } else {
1948                return false;
1949        }
1950       case(DegreesOfFreedom):
1951        switch(functionSpaceType_target) {
1952        case(ReducedDegreesOfFreedom):
1953        case(DegreesOfFreedom):
1954        case(Nodes):
1955        case(ReducedNodes):
1956        case(Elements):
1957        case(ReducedElements):
1958        case(Points):
1959        case(FaceElements):
1960        case(ReducedFaceElements):
1961        case(ContactElementsZero):
1962        case(ReducedContactElementsZero):
1963        case(ContactElementsOne):
1964        case(ReducedContactElementsOne):
1965        return true;
1966        default:
1967            stringstream temp;
1968            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1969            throw FinleyAdapterException(temp.str());
1970        }
1971        break;
1972       case(ReducedDegreesOfFreedom):
1973       switch(functionSpaceType_target) {
1974        case(ReducedDegreesOfFreedom):
1975        case(ReducedNodes):
1976        case(Elements):
1977        case(ReducedElements):
1978        case(FaceElements):
1979        case(ReducedFaceElements):
1980        case(Points):
1981        case(ContactElementsZero):
1982        case(ReducedContactElementsZero):
1983        case(ContactElementsOne):
1984        case(ReducedContactElementsOne):
1985        return true;
1986        case(Nodes):
1987        case(DegreesOfFreedom):
1988        return false;
1989        default:
1990            stringstream temp;
1991            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1992            throw FinleyAdapterException(temp.str());
1993        }
1994        break;
1995       default:
1996          stringstream temp;
1997          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1998          throw FinleyAdapterException(temp.str());
1999          break;
2000       }
2001       checkFinleyError();
2002       return false;
2003  }  }
2004    
2005  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
2006  {  {
2007      return false;     return false;
2008  }  }
2009    
2010  bool MeshAdapter::operator==(const AbstractDomain& other) const  bool MeshAdapter::operator==(const AbstractDomain& other) const
2011  {  {
2012    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);     const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2013    if (temp!=0) {     if (temp!=0) {
2014      return (m_finleyMesh==temp->m_finleyMesh);        return (m_finleyMesh==temp->m_finleyMesh);
2015    } else {     } else {
2016      return false;        return false;
2017    }     }
2018  }  }
2019    
2020  bool MeshAdapter::operator!=(const AbstractDomain& other) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
2021  {  {
2022    return !(operator==(other));     return !(operator==(other));
2023  }  }
2024    
2025  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
2026    {
2027       Finley_Mesh* mesh=m_finleyMesh.get();
2028       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2029       checkPasoError();
2030       return out;
2031    }
2032    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2033  {  {
2034     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2035       int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2036     checkPasoError();     checkPasoError();
2037     return out;     return out;
2038  }  }
2039    
2040  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2041  {  {
2042    return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2043  }  }
2044    
2045  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2046  {  {
2047    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2048  }  }
2049    
2050  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2051  {  {
2052    return function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2053  }  }
2054    
2055  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2056  {  {
2057    int out=0;     int *out = NULL;
2058    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2059    switch (functionSpaceType) {     switch (functionSpaceType) {
2060    case(Nodes):     case(Nodes):
2061      out=mesh->Nodes->Tag[sampleNo];     out=mesh->Nodes->Id;
2062      break;     break;
2063    case(Elements):     case(ReducedNodes):
2064      out=mesh->Elements->Tag[sampleNo];     out=mesh->Nodes->reducedNodesId;
2065      break;     break;
2066    case(FaceElements):     case(Elements):
2067      out=mesh->FaceElements->Tag[sampleNo];     out=mesh->Elements->Id;
2068      break;     break;
2069    case(Points):     case(ReducedElements):
2070      out=mesh->Points->Tag[sampleNo];     out=mesh->Elements->Id;
2071      break;     break;
2072    case(ContactElementsZero):     case(FaceElements):
2073      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->FaceElements->Id;
2074      break;     break;
2075    case(ContactElementsOne):     case(ReducedFaceElements):
2076      out=mesh->ContactElements->Tag[sampleNo];     out=mesh->FaceElements->Id;
2077      break;     break;
2078    case(DegreesOfFreedom):     case(Points):
2079      break;     out=mesh->Points->Id;
2080    case(ReducedDegreesOfFreedom):     break;
2081      break;     case(ContactElementsZero):
2082    default:     out=mesh->ContactElements->Id;
2083      stringstream temp;     break;
2084      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();     case(ReducedContactElementsZero):
2085      throw FinleyAdapterException(temp.str());     out=mesh->ContactElements->Id;
2086      break;     break;
2087    }     case(ContactElementsOne):
2088    return out;     out=mesh->ContactElements->Id;
2089       break;
2090       case(ReducedContactElementsOne):
2091       out=mesh->ContactElements->Id;
2092       break;
2093       case(DegreesOfFreedom):
2094       out=mesh->Nodes->degreesOfFreedomId;
2095       break;
2096       case(ReducedDegreesOfFreedom):
2097       out=mesh->Nodes->reducedDegreesOfFreedomId;
2098       break;
2099       default:
2100          stringstream temp;
2101          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2102          throw FinleyAdapterException(temp.str());
2103          break;
2104       }
2105       return out;
2106  }  }
2107  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2108  {  {
2109    int out=0,i;     int out=0;
2110    Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2111    switch (functionSpaceType) {     switch (functionSpaceType) {
2112    case(Nodes):     case(Nodes):
2113      if (mesh->Nodes!=NULL) {     out=mesh->Nodes->Tag[sampleNo];
2114        out=mesh->Nodes->Id[sampleNo];     break;
2115       case(ReducedNodes):
2116       throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2117       break;
2118       case(Elements):
2119       out=mesh->Elements->Tag[sampleNo];
2120       break;
2121       case(ReducedElements):
2122       out=mesh->Elements->Tag[sampleNo];
2123       break;
2124       case(FaceElements):
2125       out=mesh->FaceElements->Tag[sampleNo];
2126       break;
2127       case(ReducedFaceElements):
2128       out=mesh->FaceElements->Tag[sampleNo];
2129       break;
2130       case(Points):
2131       out=mesh->Points->Tag[sampleNo];
2132       break;
2133       case(ContactElementsZero):
2134       out=mesh->ContactElements->Tag[sampleNo];
2135       break;
2136       case(ReducedContactElementsZero):
2137       out=mesh->ContactElements->Tag[sampleNo];
2138       break;
2139       case(ContactElementsOne):
2140       out=mesh->ContactElements->Tag[sampleNo];
2141       break;
2142       case(ReducedContactElementsOne):
2143       out=mesh->ContactElements->Tag[sampleNo];
2144       break;
2145       case(DegreesOfFreedom):
2146       throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2147       break;
2148       case(ReducedDegreesOfFreedom):
2149       throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2150       break;
2151       default:
2152          stringstream temp;
2153          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2154          throw FinleyAdapterException(temp.str());
2155        break;        break;
2156      }     }
2157    case(Elements):     return out;
     out=mesh->Elements->Id[sampleNo];  
     break;  
   case(FaceElements):  
     out=mesh->FaceElements->Id[sampleNo];  
     break;  
   case(Points):  
     out=mesh->Points->Id[sampleNo];  
     break;  
   case(ContactElementsZero):  
     out=mesh->ContactElements->Id[sampleNo];  
     break;  
   case(ContactElementsOne):  
     out=mesh->ContactElements->Id[sampleNo];  
     break;  
   case(DegreesOfFreedom):  
     for (i=0;i<mesh->Nodes->numNodes; ++i) {  
        if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {  
           out=mesh->Nodes->Id[i];  
           break;  
        }  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     for (i=0;i<mesh->Nodes->numNodes; ++i) {  
        if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {  
           out=mesh->Nodes->Id[i];  
           break;  
        }  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   return out;  
2158  }  }
2159    
2160    
2161  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
2162  {  {
2163       Finley_Mesh* mesh=m_finleyMesh.get();
2164       escriptDataC tmp=mask.getDataC();
2165       switch(functionSpaceType) {
2166       case(Nodes):
2167       Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
2168       break;
2169       case(ReducedNodes):
2170       throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2171       break;
2172       case(DegreesOfFreedom):
2173       throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2174       break;
2175       case(ReducedDegreesOfFreedom):
2176       throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2177       break;
2178       case(Elements):
2179       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2180       break;
2181       case(ReducedElements):
2182       Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
2183       break;
2184       case(FaceElements):
2185       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2186       break;
2187       case(ReducedFaceElements):
2188       Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
2189       break;
2190       case(Points):
2191       Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
2192       break;
2193       case(ContactElementsZero):
2194       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2195       break;
2196       case(ReducedContactElementsZero):
2197       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2198       break;
2199       case(ContactElementsOne):
2200       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2201       break;
2202       case(ReducedContactElementsOne):
2203       Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
2204       break;
2205       default:
2206          stringstream temp;
2207          temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
2208          throw FinleyAdapterException(temp.str());
2209       }
2210       checkFinleyError();
2211       return;
2212    }
2213    
2214    void MeshAdapter::setTagMap(const string& name,  int tag)
2215    {
2216       Finley_Mesh* mesh=m_finleyMesh.get();
2217       Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2218       checkPasoError();
2219       // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2220    }
2221    
2222    int MeshAdapter::getTag(const string& name) const
2223    {
2224       Finley_Mesh* mesh=m_finleyMesh.get();
2225       int tag=0;
2226       tag=Finley_Mesh_getTag(mesh, name.c_str());
2227       checkPasoError();
2228       // throwStandardException("MeshAdapter::getTag is not implemented.");
2229       return tag;
2230    }
2231    
2232    bool MeshAdapter::isValidTagName(const string& name) const
2233    {
2234       Finley_Mesh* mesh=m_finleyMesh.get();
2235       return Finley_Mesh_isValidTagName(mesh,name.c_str());
2236    }
2237    
2238    string MeshAdapter::showTagNames() const
2239    {
2240       stringstream temp;
2241       Finley_Mesh* mesh=m_finleyMesh.get();
2242       Finley_TagMap* tag_map=mesh->TagMap;
2243       while (tag_map) {
2244          temp << tag_map->name;
2245          tag_map=tag_map->next;
2246          if (tag_map) temp << ", ";
2247       }
2248       return temp.str();
2249    }
2250    
2251    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2252    {
2253    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2254    escriptDataC tmp=mask.getDataC();    dim_t numTags=0;
2255    switch(functionSpaceType) {    switch(functionSpaceCode) {
2256         case(Nodes):     case(Nodes):
2257            Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);            numTags=mesh->Nodes->numTagsInUse;
2258            break;            break;
2259         case(DegreesOfFreedom):     case(ReducedNodes):
2260              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2261              break;
2262       case(DegreesOfFreedom):
2263            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2264            break;            break;
2265         case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
2266            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2267            break;            break;
2268         case(Elements):     case(Elements):
2269            Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);     case(ReducedElements):
2270              numTags=mesh->Elements->numTagsInUse;
2271              break;
2272       case(FaceElements):
2273       case(ReducedFaceElements):
2274              numTags=mesh->FaceElements->numTagsInUse;
2275              break;
2276       case(Points):
2277              numTags=mesh->Points->numTagsInUse;
2278              break;
2279       case(ContactElementsZero):
2280       case(ReducedContactElementsZero):
2281       case(ContactElementsOne):
2282       case(ReducedContactElementsOne):
2283              numTags=mesh->ContactElements->numTagsInUse;
2284              break;
2285       default:
2286          stringstream temp;
2287          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2288          throw FinleyAdapterException(temp.str());
2289      }
2290      return numTags;
2291    }
2292    
2293    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2294    {
2295      Finley_Mesh* mesh=m_finleyMesh.get();
2296      index_t* tags=NULL;
2297      switch(functionSpaceCode) {
2298       case(Nodes):
2299              tags=mesh->Nodes->tagsInUse;
2300              break;
2301       case(ReducedNodes):
2302              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2303            break;            break;
2304         case(FaceElements):     case(DegreesOfFreedom):
2305            Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);            throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2306            break;            break;
2307         case(Points):     case(ReducedDegreesOfFreedom):
2308            Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);            throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2309            break;            break;
2310         case(ContactElementsZero):     case(Elements):
2311            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);     case(ReducedElements):
2312              tags=mesh->Elements->tagsInUse;
2313              break;
2314       case(FaceElements):
2315       case(ReducedFaceElements):
2316              tags=mesh->FaceElements->tagsInUse;
2317              break;
2318       case(Points):
2319              tags=mesh->Points->tagsInUse;
2320              break;
2321       case(ContactElementsZero):
2322       case(ReducedContactElementsZero):
2323       case(ContactElementsOne):
2324       case(ReducedContactElementsOne):
2325              tags=mesh->ContactElements->tagsInUse;
2326            break;            break;
2327         case(ContactElementsOne):     default:
2328            Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);        stringstream temp;
2329          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2330          throw FinleyAdapterException(temp.str());
2331      }
2332      return tags;
2333    }
2334    
2335    
2336    bool MeshAdapter::canTag(int functionSpaceCode) const
2337    {
2338      switch(functionSpaceCode) {
2339       case(Nodes):
2340       case(Elements):
2341       case(ReducedElements):
2342       case(FaceElements):
2343       case(ReducedFaceElements):
2344       case(Points):
2345       case(ContactElementsZero):
2346       case(ReducedContactElementsZero):
2347       case(ContactElementsOne):
2348       case(ReducedContactElementsOne):
2349              return true;
2350       case(ReducedNodes):
2351       case(DegreesOfFreedom):
2352       case(ReducedDegreesOfFreedom):
2353          return false;
2354       default:
2355        return false;
2356      }
2357    }
2358    
2359    AbstractDomain::StatusType MeshAdapter::getStatus() const
2360    {
2361      Finley_Mesh* mesh=m_finleyMesh.get();
2362      return Finley_Mesh_getStatus(mesh);
2363    }
2364    
2365    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2366    {
2367      
2368      Finley_Mesh* mesh=m_finleyMesh.get();
2369      int order =-1;
2370      switch(functionSpaceCode) {
2371       case(Nodes):
2372       case(DegreesOfFreedom):
2373              order=mesh->approximationOrder;
2374              break;
2375       case(ReducedNodes):
2376       case(ReducedDegreesOfFreedom):
2377              order=mesh->reducedApproximationOrder;
2378              break;
2379       case(Elements):
2380       case(FaceElements):
2381       case(Points):
2382       case(ContactElementsZero):
2383       case(ContactElementsOne):
2384              order=mesh->integrationOrder;
2385              break;
2386       case(ReducedElements):
2387       case(ReducedFaceElements):
2388       case(ReducedContactElementsZero):
2389       case(ReducedContactElementsOne):
2390              order=mesh->reducedIntegrationOrder;
2391            break;            break;
2392         default:     default:
2393            stringstream temp;        stringstream temp;
2394            temp << "Error - Finley does not know anything about function space type " << functionSpaceType;        temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2395            throw FinleyAdapterException(temp.str());        throw FinleyAdapterException(temp.str());
2396    }    }
2397    checkFinleyError();    return order;
2398    return;  }
2399    
2400    bool MeshAdapter::supportsContactElements() const
2401    {
2402      return true;
2403  }  }
2404    
2405    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2406    {
2407      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2408    }
2409    
2410    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2411    {
2412      Finley_ReferenceElementSet_dealloc(m_refSet);
2413    }
2414    
2415  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.793  
changed lines
  Added in v.3355

  ViewVC Help
Powered by ViewVC 1.1.26