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

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.472  
changed lines
  Added in v.3344

  ViewVC Help
Powered by ViewVC 1.1.26