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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26