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

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

  ViewVC Help
Powered by ViewVC 1.1.26