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

Legend:
Removed from v.1044  
changed lines
  Added in v.2635

  ViewVC Help
Powered by ViewVC 1.1.26