/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.1116  
changed lines
  Added in v.3242

  ViewVC Help
Powered by ViewVC 1.1.26