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

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

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

trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 1339 by ksteube, Wed Nov 7 01:53:12 2007 UTC
# Line 1  Line 1 
1  /*  
2   ******************************************************************************  /* $Id$ */
3   *                                                                            *  
4   *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  /*******************************************************
5   *                                                                            *   *
6   * This software is the property of ACcESS. No part of this code              *   *           Copyright 2003-2007 by ACceSS MNRF
7   * may be copied in any form or by any means without the expressed written    *   *       Copyright 2007 by University of Queensland
8   * consent of ACcESS.  Copying, use or modification of this software          *   *
9   * by any unauthorised person is illegal unless that person has a software    *   *                http://esscc.uq.edu.au
10   * license agreement with ACcESS.                                             *   *        Primary Business: Queensland, Australia
11   *                                                                            *   *  Licensed under the Open Software License version 3.0
12   ******************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16    #include "MeshAdapter.h"
17    #include "escript/Data.h"
18    #include "escript/DataFactory.h"
19  extern "C" {  extern "C" {
20  #include "finley/finleyC/Finley.h"  #include "escript/blocktimer.h"
21  #include "finley/finleyC/Assemble.h"  }
 #include "finley/finleyC/Mesh.h"  
 #include "finley/finleyC/Finley.h"  
 #include "finley/finleyC/System.h"  
 }  
 #include "finley/CPPAdapter/SystemMatrixAdapter.h"  
 #include "finley/CPPAdapter/MeshAdapter.h"  
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/FinleyAdapterException.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/Data.h"  
 #include "escript/Data/DataArrayView.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/DataFactory.h"  
 #include <iostream>  
22  #include <vector>  #include <vector>
 #include <sstream>  
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
26    
27  namespace finley {  namespace finley {
28    
 struct null_deleter  
 {  
   void operator()(void const *ptr) const  
   {  
   }  
 };  
   
29  //  //
30  // define the statics  // define the static constants
31  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
32  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
33  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
34  const int MeshAdapter::Nodes=FINLEY_NODES;  const int MeshAdapter::Nodes=FINLEY_NODES;
35    const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
36  const int MeshAdapter::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=FINLEY_ELEMENTS;
37    const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
38  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
39    const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
40  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=FINLEY_POINTS;
41  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
42    const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
43  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
44    const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
45    
46  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
47  {  {
48    setFunctionSpaceTypeNames();    setFunctionSpaceTypeNames();
49    //    //
50    // 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
51    // for us.    // for us.
52    m_finleyMesh.reset(finleyMesh,null_deleter());    m_finleyMesh.reset(finleyMesh,null_deleter());
53  }  }
54    
55  //  //
56  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
57  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
# Line 77  MeshAdapter::~MeshAdapter() Line 66  MeshAdapter::~MeshAdapter()
66    // I hope the case for the pointer being zero has been taken care of.    // I hope the case for the pointer being zero has been taken care of.
67    //  cout << "In MeshAdapter destructor." << endl;    //  cout << "In MeshAdapter destructor." << endl;
68    if (m_finleyMesh.unique()) {    if (m_finleyMesh.unique()) {
69      //   cout << "Calling dealloc." << endl;      Finley_Mesh_free(m_finleyMesh.get());
     Finley_Mesh_dealloc(m_finleyMesh.get());  
     //   cout << "Finished dealloc." << endl;  
70    }    }
71  }  }
72    
73    int MeshAdapter::getMPISize() const
74    {
75       return m_finleyMesh.get()->MPIInfo->size;
76    }
77    int MeshAdapter::getMPIRank() const
78    {
79       return m_finleyMesh.get()->MPIInfo->rank;
80    }
81    
82    
83  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
84     return m_finleyMesh.get();     return m_finleyMesh.get();
85  }  }
86    
87  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
88  {  {
89    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
90    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
91    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
92    checkFinleyError();    checkFinleyError();
93      TMPMEMFREE(fName);
94    }
95    
96    void MeshAdapter::Print_Mesh_Info(const bool full=false) const
97    {
98      Finley_PrintMesh_Info(m_finleyMesh.get(), full);
99  }  }
100    
101  // void MeshAdapter::getTagList(int functionSpaceType,  void MeshAdapter::dump(const std::string& fileName) const
102  //                  int* numTags) const  {
103  // {    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
104  //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);    strcpy(fName,fileName.c_str());
105  //   return;    Finley_Mesh_dump(m_finleyMesh.get(),fName);
106  // }    checkFinleyError();
107      TMPMEMFREE(fName);
108    }
109    
110  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
111  {  {
# Line 134  void MeshAdapter::setFunctionSpaceTypeNa Line 139  void MeshAdapter::setFunctionSpaceTypeNa
139    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
140      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
141    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
142        (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
143      m_functionSpaceTypeNames.insert
144      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
145    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
146        (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
147      m_functionSpaceTypeNames.insert
148      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
149    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
150        (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
151      m_functionSpaceTypeNames.insert
152      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
153    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
154      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
155    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
156        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
157      m_functionSpaceTypeNames.insert
158      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
159      m_functionSpaceTypeNames.insert
160        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
161  }  }
162    
163  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
164  {  {
165    return Nodes;    return Nodes;
166  }  }
167    int MeshAdapter::getReducedContinuousFunctionCode() const
168    {
169      return ReducedNodes;
170    }
171    
172  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
173  {  {
174    return Elements;    return Elements;
175  }  }
176    int MeshAdapter::getReducedFunctionCode() const
177    {
178      return ReducedElements;
179    }
180    
181  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
182  {  {
183    return FaceElements;    return FaceElements;
184  }  }
185    int MeshAdapter::getReducedFunctionOnBoundaryCode() const
186    {
187      return ReducedFaceElements;
188    }
189    
190  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
191  {  {
192    return ContactElementsZero;    return ContactElementsZero;
193  }  }
194    int MeshAdapter::getReducedFunctionOnContactZeroCode() const
195    {
196      return ReducedContactElementsZero;
197    }
198    
199  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
200  {  {
201    return ContactElementsOne;    return ContactElementsOne;
202  }  }
203    int MeshAdapter::getReducedFunctionOnContactOneCode() const
204    {
205      return ReducedContactElementsOne;
206    }
207    
208  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
209  {  {
210    return DegreesOfFreedom;    return DegreesOfFreedom;
211  }  }
212    
213  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
214  {  {
215    return ReducedDegreesOfFreedom;    return ReducedDegreesOfFreedom;
216  }  }
217    
218  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionCode() const
219  {  {
220    return Points;    return Points;
221  }  }
222    
 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  
 {  
   //  
   // It is assumed the sampleNo has been checked  
   // before calling this function.  
   int* tagList;  
   int numTags;  
   getTagList(functionSpaceType, &tagList, &numTags);  
   return tagList[sampleNo];  
 }  
   
 //  
 // returns a pointer to the tag list of samples of functionSpaceType  
 //  
 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,  
                  int* numTags) const  
 {  
   **tagList=0;  
   *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;  
   }  
   return;  
 }  
   
223  //  //
224  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
225  //  //
# Line 267  int MeshAdapter::getDim() const Line 229  int MeshAdapter::getDim() const
229    checkFinleyError();    checkFinleyError();
230    return numDim;    return numDim;
231  }  }
232    
233  //  //
234  // return the number of data points per sample and the number of samples  // return the number of data points per sample and the number of samples
235  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
# Line 279  pair<int,int> MeshAdapter::getDataShape( Line 242  pair<int,int> MeshAdapter::getDataShape(
242     switch (functionSpaceCode) {     switch (functionSpaceCode) {
243        case(Nodes):        case(Nodes):
244             numDataPointsPerSample=1;             numDataPointsPerSample=1;
245             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;             numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
246               break;
247          case(ReducedNodes):
248               numDataPointsPerSample=1;
249               numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
250             break;             break;
251        case(Elements):        case(Elements):
252             if (mesh->Elements!=NULL) {             if (mesh->Elements!=NULL) {
# Line 287  pair<int,int> MeshAdapter::getDataShape( Line 254  pair<int,int> MeshAdapter::getDataShape(
254               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
255             }             }
256             break;             break;
257          case(ReducedElements):
258               if (mesh->Elements!=NULL) {
259                 numSamples=mesh->Elements->numElements;
260                 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
261               }
262               break;
263        case(FaceElements):        case(FaceElements):
264             if (mesh->FaceElements!=NULL) {             if (mesh->FaceElements!=NULL) {
265                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
266                  numSamples=mesh->FaceElements->numElements;                  numSamples=mesh->FaceElements->numElements;
267             }             }
268             break;             break;
269          case(ReducedFaceElements):
270               if (mesh->FaceElements!=NULL) {
271                    numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
272                    numSamples=mesh->FaceElements->numElements;
273               }
274               break;
275        case(Points):        case(Points):
276             if (mesh->Points!=NULL) {             if (mesh->Points!=NULL) {
277               numDataPointsPerSample=1;               numDataPointsPerSample=1;
# Line 305  pair<int,int> MeshAdapter::getDataShape( Line 284  pair<int,int> MeshAdapter::getDataShape(
284               numSamples=mesh->ContactElements->numElements;               numSamples=mesh->ContactElements->numElements;
285             }             }
286             break;             break;
287          case(ReducedContactElementsZero):
288               if (mesh->ContactElements!=NULL) {
289                 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
290                 numSamples=mesh->ContactElements->numElements;
291               }
292               break;
293        case(ContactElementsOne):        case(ContactElementsOne):
294             if (mesh->ContactElements!=NULL) {             if (mesh->ContactElements!=NULL) {
295               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
296               numSamples=mesh->ContactElements->numElements;               numSamples=mesh->ContactElements->numElements;
297             }             }
298             break;             break;
299          case(ReducedContactElementsOne):
300               if (mesh->ContactElements!=NULL) {
301                 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
302                 numSamples=mesh->ContactElements->numElements;
303               }
304               break;
305        case(DegreesOfFreedom):        case(DegreesOfFreedom):
306             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
307               numDataPointsPerSample=1;               numDataPointsPerSample=1;
308               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
309             }             }
310             break;             break;
311        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
312             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
313               numDataPointsPerSample=1;               numDataPointsPerSample=1;
314               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
315             }             }
316             break;             break;
317        default:        default:
318          stringstream temp;          stringstream temp;
319          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
320          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
321          break;          break;
322        }        }
323        return pair<int,int>(numDataPointsPerSample,numSamples);        return pair<int,int>(numDataPointsPerSample,numSamples);
324  }  }
325    
326  //  //
327  // 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:
328  //  //
329  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
330                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
331                       const Data& A, const Data& B, const Data& C,const  Data& D,const  Data& X,const  Data& Y) const                       const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
332  {                       const escript::Data& d, const escript::Data& y,
333                         const escript::Data& d_contact,const escript::Data& y_contact) const
334    {
335       escriptDataC _rhs=rhs.getDataC();
336       escriptDataC _A  =A.getDataC();
337       escriptDataC _B=B.getDataC();
338       escriptDataC _C=C.getDataC();
339       escriptDataC _D=D.getDataC();
340       escriptDataC _X=X.getDataC();
341       escriptDataC _Y=Y.getDataC();
342       escriptDataC _d=d.getDataC();
343       escriptDataC _y=y.getDataC();
344       escriptDataC _d_contact=d_contact.getDataC();
345       escriptDataC _y_contact=y_contact.getDataC();
346    
347     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
348     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),  
349                         &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
350       checkFinleyError();
351    
352       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
353       checkFinleyError();
354    
355       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
356     checkFinleyError();     checkFinleyError();
357  }  }
358  //  
359  // adds Robin boundary conditions as natural boundary condition into a given stiffness matrix and right hand side:  void  MeshAdapter::addPDEToLumpedSystem(
360  //                       escript::Data& mat,
361  void MeshAdapter::addRobinConditionsToSystem(                       const escript::Data& D,
362                       SystemMatrixAdapter& mat, Data& rhs,                       const escript::Data& d) const
363                       const Data& d, const Data& y) const  {
364  {     escriptDataC _mat=mat.getDataC();
365       escriptDataC _D=D.getDataC();
366       escriptDataC _d=d.getDataC();
367    
368     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
369     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
370                    mat.getFinley_SystemMatrix(),     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
371                    &(rhs.getDataC()),     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
372                                    &(d.getDataC()),&(y.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Mean_out);  
373     checkFinleyError();     checkFinleyError();
374  }  }
375    
376    
377  //  //
378  // adds contact conditions as natural boundary condition into a given stiffness matrix and right hand side:  // adds linear PDE of second order into the right hand side only
379  //  //
380  void MeshAdapter::addContactToSystem(  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
                      SystemMatrixAdapter& mat, Data& rhs,  
                      const Data& d_contact,const Data& y_contact) const  
381  {  {
382     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
383     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
384                    mat.getFinley_SystemMatrix(),     escriptDataC _rhs=rhs.getDataC();
385                    &(rhs.getDataC()),     escriptDataC _X=X.getDataC();
386                                    &(d_contact.getDataC()),     escriptDataC _Y=Y.getDataC();
387                    &(y_contact.getDataC()),     escriptDataC _y=y.getDataC();
388                                    Finley_Assemble_handelShapeMissMatch_Step_out);     escriptDataC _y_contact=y_contact.getDataC();
389    
390       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
391       checkFinleyError();
392    
393       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
394       checkFinleyError();
395    
396       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
397     checkFinleyError();     checkFinleyError();
398  }  }
399    
400  //  //
401  // interpolates data between different function spaces:  // interpolates data between different function spaces:
402  //  //
403  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
404  {  {
405    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
406    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
407    if (inDomain!=*this)    if (inDomain!=*this)  
408       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
409    if (targetDomain!=*this)    if (targetDomain!=*this)
410       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
411    
412    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
413      escriptDataC _target=target.getDataC();
414      escriptDataC _in=in.getDataC();
415    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
416       case(Nodes):       case(Nodes):
417          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
418             case(Nodes):             case(Nodes):
419               case(ReducedNodes):
420               case(DegreesOfFreedom):
421             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
422                   Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
423                   break;
424               case(Elements):
425               case(ReducedElements):
426                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
427                   break;
428               case(FaceElements):
429               case(ReducedFaceElements):
430                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
431                   break;
432               case(Points):
433                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
434                   break;
435               case(ContactElementsZero):
436               case(ReducedContactElementsZero):
437                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
438                   break;
439               case(ContactElementsOne):
440               case(ReducedContactElementsOne):
441                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
442                   break;
443               default:
444                   stringstream temp;
445                   temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
446                   throw FinleyAdapterException(temp.str());
447                   break;
448            }
449            break;
450         case(ReducedNodes):
451            switch(target.getFunctionSpace().getTypeCode()) {
452               case(Nodes):
453               case(ReducedNodes):
454             case(DegreesOfFreedom):             case(DegreesOfFreedom):
455                 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));             case(ReducedDegreesOfFreedom):
456                   Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
457                 break;                 break;
458             case(Elements):             case(Elements):
459                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));             case(ReducedElements):
460                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
461                 break;                 break;
462             case(FaceElements):             case(FaceElements):
463                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedFaceElements):
464                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
465                 break;                 break;
466             case(Points):             case(Points):
467                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
468                 break;                 break;
469             case(ContactElementsZero):             case(ContactElementsZero):
470                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsZero):
471                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
472                 break;                 break;
473             case(ContactElementsOne):             case(ContactElementsOne):
474                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsOne):
475                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
476                 break;                 break;
477             default:             default:
478                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
479                 sprintf(Finley_ErrorMsg,"Interpolation on Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());                 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
480                   throw FinleyAdapterException(temp.str());
481                 break;                 break;
482          }          }
483          break;          break;
484       case(Elements):       case(Elements):
485          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
486             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
487            } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
488               Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
489          } else {          } else {
490             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
491             sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");          }
492            break;
493         case(ReducedElements):
494            if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
495               Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
496            } else {
497               throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
498          }          }
499          break;          break;
500       case(FaceElements):       case(FaceElements):
501          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
502             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
503            } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
504               Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
505          } else {          } else {
506             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
507             sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");         }
508             break;         break;
509         case(ReducedFaceElements):
510            if (target.getFunctionSpace().getTypeCode()==FaceElements) {
511               Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
512            } else {
513               throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
514         }         }
515           break;
516       case(Points):       case(Points):
517          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
518             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
519          } else {          } else {
520             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
            break;  
521          }          }
522          break;          break;
523       case(ContactElementsZero):       case(ContactElementsZero):
524       case(ContactElementsOne):       case(ContactElementsOne):
525          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
526             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
527            } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
528               Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
529          } else {          } else {
530             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on contact elements possible.");  
531             break;             break;
532          }          }
533          break;          break;
534       case(DegreesOfFreedom):       case(ReducedContactElementsZero):
535         case(ReducedContactElementsOne):
536            if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
537               Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
538            } else {
539               throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
540               break;
541            }
542            break;
543         case(DegreesOfFreedom):      
544          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
545             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
546             case(DegreesOfFreedom):             case(DegreesOfFreedom):
547                  Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
548                  break;
549    
550             case(Nodes):             case(Nodes):
551                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));             case(ReducedNodes):
552                  if (getMPISize()>1) {
553                      escript::Data temp=escript::Data(in);
554                      temp.expand();
555                      escriptDataC _in2 = temp.getDataC();
556                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
557                  } else {
558                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
559                  }
560                break;                break;
561             case(Elements):             case(Elements):
562                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));             case(ReducedElements):
563                  if (getMPISize()>1) {
564                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
565                     escriptDataC _in2 = temp.getDataC();
566                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
567                  } else {
568                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
569                  }
570                break;                break;
571             case(FaceElements):             case(FaceElements):
572                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedFaceElements):
573                  if (getMPISize()>1) {
574                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
575                     escriptDataC _in2 = temp.getDataC();
576                     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
577    
578                  } else {
579                     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
580                  }
581                break;                break;
582             case(Points):             case(Points):
583                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                if (getMPISize()>1) {
584                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
585                     escriptDataC _in2 = temp.getDataC();
586                  } else {
587                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
588                  }
589                break;                break;
590             case(ContactElementsZero):             case(ContactElementsZero):
591             case(ContactElementsOne):             case(ContactElementsOne):
592                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsZero):
593               case(ReducedContactElementsOne):
594                  if (getMPISize()>1) {
595                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
596                     escriptDataC _in2 = temp.getDataC();
597                     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
598                  } else {
599                     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
600                  }
601               break;               break;
602             default:             default:
603               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
604               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
605                 throw FinleyAdapterException(temp.str());
606               break;               break;
607          }          }
608          break;          break;
609       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
610         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
           case(ReducedDegreesOfFreedom):  
611            case(Nodes):            case(Nodes):
612               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
613                 break;
614              case(ReducedNodes):
615                  if (getMPISize()>1) {
616                      escript::Data temp=escript::Data(in);
617                      temp.expand();
618                      escriptDataC _in2 = temp.getDataC();
619                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
620                  } else {
621                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
622                  }
623                  break;
624              case(DegreesOfFreedom):
625                 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
626                 break;
627              case(ReducedDegreesOfFreedom):
628                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
629               break;               break;
630            case(Elements):            case(Elements):
631               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));            case(ReducedElements):
632                  if (getMPISize()>1) {
633                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
634                     escriptDataC _in2 = temp.getDataC();
635                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
636                  } else {
637                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
638                 }
639               break;               break;
640            case(FaceElements):            case(FaceElements):
641               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));            case(ReducedFaceElements):
642                  if (getMPISize()>1) {
643                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
644                     escriptDataC _in2 = temp.getDataC();
645                     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
646                  } else {
647                     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
648                  }
649               break;               break;
650            case(Points):            case(Points):
651               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                if (getMPISize()>1) {
652                     escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
653                     escriptDataC _in2 = temp.getDataC();
654                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
655                  } else {
656                     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
657                 }
658               break;               break;
659            case(ContactElementsZero):            case(ContactElementsZero):
660            case(ContactElementsOne):            case(ContactElementsOne):
661               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));            case(ReducedContactElementsZero):
662               break;            case(ReducedContactElementsOne):
663            case(DegreesOfFreedom):                if (getMPISize()>1) {
664               Finley_ErrorCode=TYPE_ERROR;                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
665               sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");                   escriptDataC _in2 = temp.getDataC();
666                     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
667                  } else {
668                     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
669                 }
670               break;               break;
671            default:            default:
672               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
673               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
674                 throw FinleyAdapterException(temp.str());
675               break;               break;
676         }         }
677         break;         break;
678       default:       default:
679          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
680          sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",in.getFunctionSpace().getTypeCode());          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
681            throw FinleyAdapterException(temp.str());
682          break;          break;
683    }    }
684    checkFinleyError();    checkFinleyError();
# Line 518  void MeshAdapter::interpolateOnDomain(Da Line 687  void MeshAdapter::interpolateOnDomain(Da
687  //  //
688  // copies the locations of sample points into x:  // copies the locations of sample points into x:
689  //  //
690  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
691  {  {
692    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
693    if (argDomain!=*this)    if (argDomain!=*this)
# Line 526  void MeshAdapter::setToX(Data& arg) cons Line 695  void MeshAdapter::setToX(Data& arg) cons
695    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
696    // in case of values node coordinates we can do the job directly:    // in case of values node coordinates we can do the job directly:
697    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
698       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       escriptDataC _arg=arg.getDataC();
699         Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
700    } else {    } else {
701       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
702       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       escriptDataC _tmp_data=tmp_data.getDataC();
703         Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
704       // this is then interpolated onto arg:       // this is then interpolated onto arg:
705       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
706    }    }
707    checkFinleyError();    checkFinleyError();
708  }  }
709    
710  //  //
711  // 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:
712  //  //
713  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
714  {  {
715    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
716    if (normalDomain!=*this)    if (normalDomain!=*this)
717       throw FinleyAdapterException("Error - Illegal domain of normal locations");       throw FinleyAdapterException("Error - Illegal domain of normal locations");
718    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
719      escriptDataC _normal=normal.getDataC();
720    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
721      case(Nodes):      case(Nodes):
722        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
723        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");        break;
724        case(ReducedNodes):
725          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
726        break;        break;
727      case(Elements):      case(Elements):
728        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
729        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");        break;
730        case(ReducedElements):
731          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
732        break;        break;
733      case (FaceElements):      case (FaceElements):
734        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
735          break;
736        case (ReducedFaceElements):
737          Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
738        break;        break;
739      case(Points):      case(Points):
740        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
       sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for point elements");  
741        break;        break;
742      case (ContactElementsOne):      case (ContactElementsOne):
743      case (ContactElementsZero):      case (ContactElementsZero):
744        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
745          break;
746        case (ReducedContactElementsOne):
747        case (ReducedContactElementsZero):
748          Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
749        break;        break;
750      case(DegreesOfFreedom):      case(DegreesOfFreedom):
751        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
       sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for degrees of freedom.");  
752        break;        break;
753      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
754        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
       sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for reduced degrees of freedom.");  
755        break;        break;
756      default:      default:
757        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
758        sprintf(Finley_ErrorMsg,"Normal Vectors: Finley does not know anything about function space type %d",normal.getFunctionSpace().getTypeCode());        temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
759          throw FinleyAdapterException(temp.str());
760        break;        break;
761    }    }
762    checkFinleyError();    checkFinleyError();
763  }  }
764    
765  //  //
766  // interpolates data to other domain:  // interpolates data to other domain:
767  //  //
768  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
769  {  {
770    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
771    if (targetDomain!=*this)    if (targetDomain!=*this)
772       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
773    
774    Finley_ErrorCode=SYSTEM_ERROR;    throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
   sprintf(Finley_ErrorMsg,"Finley does not allow interpolation across domains yet.");  
   checkFinleyError();  
775  }  }
776    
777  //  //
778  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
779  //  //
780  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
781  {  {
782    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
783    if (argDomain!=*this)    if (argDomain!=*this)
784       throw FinleyAdapterException("Error - Illegal domain of integration kernel");       throw FinleyAdapterException("Error - Illegal domain of integration kernel");
785    
786      double blocktimer_start = blocktimer_time();
787    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
788      escriptDataC _temp;
789      escript::Data temp;
790      escriptDataC _arg=arg.getDataC();
791    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
792       case(Nodes):       case(Nodes):
793          Finley_ErrorCode=TYPE_ERROR;          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
794          sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");          _temp=temp.getDataC();
795            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
796            break;
797         case(ReducedNodes):
798            temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
799            _temp=temp.getDataC();
800            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
801          break;          break;
802       case(Elements):       case(Elements):
803          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
804            break;
805         case(ReducedElements):
806            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
807          break;          break;
808       case(FaceElements):       case(FaceElements):
809          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
810            break;
811         case(ReducedFaceElements):
812            Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
813          break;          break;
814       case(Points):       case(Points):
815          Finley_ErrorCode=TYPE_ERROR;          throw FinleyAdapterException("Error - Integral of data on points is not supported.");
         sprintf(Finley_ErrorMsg,"Integral of data on points is not supported.");  
816          break;          break;
817       case(ContactElementsZero):       case(ContactElementsZero):
818          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
819            break;
820         case(ReducedContactElementsZero):
821            Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
822          break;          break;
823       case(ContactElementsOne):       case(ContactElementsOne):
824          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
825            break;
826         case(ReducedContactElementsOne):
827            Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
828          break;          break;
829       case(DegreesOfFreedom):       case(DegreesOfFreedom):
830          Finley_ErrorCode=TYPE_ERROR;          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
831          sprintf(Finley_ErrorMsg,"Integral of data on degrees of freedom is not supported.");          _temp=temp.getDataC();
832            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
833          break;          break;
834       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
835          Finley_ErrorCode=TYPE_ERROR;          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
836          sprintf(Finley_ErrorMsg,"Integral of data on reduced degrees of freedom is not supported.");          _temp=temp.getDataC();
837            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
838          break;          break;
839       default:       default:
840          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
841          sprintf(Finley_ErrorMsg,"Integrals: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());          temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
842            throw FinleyAdapterException(temp.str());
843          break;          break;
844    }    }
845    checkFinleyError();    checkFinleyError();
846      blocktimer_increment("integrate()", blocktimer_start);
847  }  }
848    
849  //  //
850  // calculates the gradient of arg:  // calculates the gradient of arg:
851  //  //
852  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
853  {  {
854    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
855    if (argDomain!=*this)    if (argDomain!=*this)
# Line 651  void MeshAdapter::setToGradient(Data& gr Line 859  void MeshAdapter::setToGradient(Data& gr
859       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
860    
861    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
862      escriptDataC _grad=grad.getDataC();
863      escriptDataC nodeDataC;
864      escript::Data temp;
865      if (getMPISize()>1) {
866          if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
867            temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );
868            nodeDataC = temp.getDataC();
869          } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
870            temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );
871            nodeDataC = temp.getDataC();
872          } else {
873            nodeDataC = arg.getDataC();
874          }
875      } else {
876         nodeDataC = arg.getDataC();
877      }
878    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
879         case(Nodes):         case(Nodes):
880            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
881            sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");            break;
882           case(ReducedNodes):
883              throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
884            break;            break;
885         case(Elements):         case(Elements):
886            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
887              break;
888           case(ReducedElements):
889              Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
890            break;            break;
891         case(FaceElements):         case(FaceElements):
892            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
893              break;
894           case(ReducedFaceElements):
895              Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
896            break;            break;
897         case(Points):         case(Points):
898            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
899            break;            break;
900         case(ContactElementsZero):         case(ContactElementsZero):
901            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
902              break;
903           case(ReducedContactElementsZero):
904              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
905            break;            break;
906         case(ContactElementsOne):         case(ContactElementsOne):
907            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
908              break;
909           case(ReducedContactElementsOne):
910              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
911            break;            break;
912         case(DegreesOfFreedom):         case(DegreesOfFreedom):
913            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at degrees of freedom is not supported.");  
914            break;            break;
915         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
916            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at reduced degrees of freedom is not supported.");  
917            break;            break;
918         default:         default:
919            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
920            sprintf(Finley_ErrorMsg,"Gradient: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());            temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
921              throw FinleyAdapterException(temp.str());
922            break;            break;
923    }    }
924    checkFinleyError();    checkFinleyError();
925  }  }
926    
927  //  //
928  // returns the size of elements:  // returns the size of elements:
929  //  //
930  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
931  {  {
932    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
933    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
934    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
935         case(Nodes):         case(Nodes):
936            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
937            sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");            break;
938           case(ReducedNodes):
939              throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
940            break;            break;
941         case(Elements):         case(Elements):
942            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
943            break;            break;
944           case(ReducedElements):
945              Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
946              break;
947         case(FaceElements):         case(FaceElements):
948            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
949            break;            break;
950           case(ReducedFaceElements):
951              Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
952              break;
953         case(Points):         case(Points):
954            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
955            break;            break;
956         case(ContactElementsZero):         case(ContactElementsZero):
957         case(ContactElementsOne):         case(ContactElementsOne):
958            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
959            break;            break;
960           case(ReducedContactElementsZero):
961           case(ReducedContactElementsOne):
962              Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
963              break;
964         case(DegreesOfFreedom):         case(DegreesOfFreedom):
965            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
           sprintf(Finley_ErrorMsg,"Size of degrees of freedom is not supported.");  
966            break;            break;
967         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
968            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
           sprintf(Finley_ErrorMsg,"Size of reduced degrees of freedom is not supported.");  
969            break;            break;
970         default:         default:
971            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
972            sprintf(Finley_ErrorMsg,"Element size: Finley does not know anything about function space type %d",size.getFunctionSpace().getTypeCode());            temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
973              throw FinleyAdapterException(temp.str());
974            break;            break;
975    }    }
976    checkFinleyError();    checkFinleyError();
977  }  }
978    
979  // sets the location of nodes:  // sets the location of nodes:
980  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
981  {  {
982    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
983      escriptDataC tmp;
984    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
985    if (newDomain!=*this)    if (newDomain!=*this)
986       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
987    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
988      Finley_Mesh_setCoordinates(mesh,&tmp);
989    checkFinleyError();    checkFinleyError();
990  }  }
991    
992  // saves a data array in openDX format:  // saves a data array in openDX format:
993  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
994  {  {
995    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
996        const int num_data=boost::python::extract<int>(arg.attr("__len__")());
997      /* win32 refactor */
998      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
999      for(int i=0;i<num_data;i++)
1000      {
1001        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1002      }
1003    
1004      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1005      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1006      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1007    
1008      boost::python::list keys=arg.keys();
1009      for (int i=0;i<num_data;++i) {
1010             std::string n=boost::python::extract<std::string>(keys[i]);
1011             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1012             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1013                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
1014             data[i]=d.getDataC();
1015             ptr_data[i]=&(data[i]);
1016             c_names[i]=&(names[i][0]);
1017             if (n.length()>MAX_namelength-1) {
1018                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1019                c_names[i][MAX_namelength-1]='\0';
1020             } else {
1021                strcpy(c_names[i],n.c_str());
1022             }
1023        }
1024        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1025        checkFinleyError();
1026        
1027          /* win32 refactor */
1028      TMPMEMFREE(c_names);
1029      TMPMEMFREE(data);
1030      TMPMEMFREE(ptr_data);
1031      for(int i=0;i<num_data;i++)
1032      {
1033        TMPMEMFREE(names[i]);
1034      }
1035      TMPMEMFREE(names);
1036    
1037        return;
1038    }
1039    
1040    // saves a data array in openVTK format:
1041    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1042    {
1043        int MAX_namelength=256;
1044        const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1045      /* win32 refactor */
1046      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1047      for(int i=0;i<num_data;i++)
1048      {
1049        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1050      }
1051    
1052      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1053      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1054      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1055    
1056        boost::python::list keys=arg.keys();
1057        for (int i=0;i<num_data;++i) {
1058             std::string n=boost::python::extract<std::string>(keys[i]);
1059             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1060             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1061                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
1062             data[i]=d.getDataC();
1063             ptr_data[i]=&(data[i]);
1064             c_names[i]=&(names[i][0]);
1065             if (n.length()>MAX_namelength-1) {
1066                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1067                c_names[i][MAX_namelength-1]='\0';
1068             } else {
1069                strcpy(c_names[i],n.c_str());
1070             }
1071        }
1072        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1073    
1074    checkFinleyError();    checkFinleyError();
1075      /* win32 refactor */
1076      TMPMEMFREE(c_names);
1077      TMPMEMFREE(data);
1078      TMPMEMFREE(ptr_data);
1079      for(int i=0;i<num_data;i++)
1080      {
1081        TMPMEMFREE(names[i]);
1082      }
1083      TMPMEMFREE(names);
1084    
1085        return;
1086  }  }
1087                                                                                                                                                                            
1088                                                                                                                                                                            
1089  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1090  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1091                        const int row_blocksize,                        const int row_blocksize,
1092                        const escript::FunctionSpace& row_functionspace,                        const escript::FunctionSpace& row_functionspace,
1093                        const int column_blocksize,                        const int column_blocksize,
1094                        const escript::FunctionSpace& column_functionspace,                        const escript::FunctionSpace& column_functionspace,
1095                        const int type,                        const int type) const
                       const bool sym) const  
1096  {  {
1097      int reduceRowOrder=0;      int reduceRowOrder=0;
1098      int reduceColOrder=0;      int reduceColOrder=0;
# Line 778  SystemMatrixAdapter MeshAdapter::newSyst Line 1119  SystemMatrixAdapter MeshAdapter::newSyst
1119          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");          throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1120      }      }
1121      // generate matrix:      // generate matrix:
1122      Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(getFinley_Mesh(),type,sym?1:0,      
1123                                                                   row_blocksize,reduceRowOrder,      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
                                                                  column_blocksize,reduceColOrder);  
1124      checkFinleyError();      checkFinleyError();
1125        Paso_SystemMatrix* fsystemMatrix;
1126        int trilinos = 0;
1127        if (trilinos) {
1128    #ifdef TRILINOS
1129          /* Allocation Epetra_VrbMatrix here */
1130    #endif
1131        }
1132        else {
1133          fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1134        }
1135        checkPasoError();
1136        Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1137      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1138  }  }
1139    
1140  //  //
1141  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
1142  // TODO:  // TODO:
1143  //  //
1144    
1145  //  //
1146  // 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:
1147  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
# Line 803  bool MeshAdapter::isCellOriented(int fun Line 1157  bool MeshAdapter::isCellOriented(int fun
1157         case(Points):         case(Points):
1158         case(ContactElementsZero):         case(ContactElementsZero):
1159         case(ContactElementsOne):         case(ContactElementsOne):
1160           case(ReducedElements):
1161           case(ReducedFaceElements):
1162           case(ReducedContactElementsZero):
1163           case(ReducedContactElementsOne):
1164            return true;            return true;
1165            break;            break;
1166         default:         default:
1167            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
1168            sprintf(Finley_ErrorMsg,"Cell: Finley does not know anything about function space type %d",functionSpaceCode);            temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1169              throw FinleyAdapterException(temp.str());
1170            break;            break;
1171    }    }
1172    checkFinleyError();    checkFinleyError();
1173    return false;    return false;
1174  }  }
1175    
1176  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1177  {  {
1178    switch(functionSpaceType_source) {    switch(functionSpaceType_source) {
1179       case(Nodes):       case(Nodes):
1180          switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1181             case(Nodes):             case(Nodes):
1182               case(ReducedNodes):
1183             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
1184             case(DegreesOfFreedom):             case(DegreesOfFreedom):
1185             case(Elements):             case(Elements):
1186               case(ReducedElements):
1187             case(FaceElements):             case(FaceElements):
1188               case(ReducedFaceElements):
1189             case(Points):             case(Points):
1190             case(ContactElementsZero):             case(ContactElementsZero):
1191               case(ReducedContactElementsZero):
1192             case(ContactElementsOne):             case(ContactElementsOne):
1193               case(ReducedContactElementsOne):
1194                 return true;                 return true;
1195             default:             default:
1196                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
1197                 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1198                   throw FinleyAdapterException(temp.str());
1199            }
1200            break;
1201         case(ReducedNodes):
1202            switch(functionSpaceType_target) {
1203               case(ReducedNodes):
1204               case(ReducedDegreesOfFreedom):
1205               case(Elements):
1206               case(ReducedElements):
1207               case(FaceElements):
1208               case(ReducedFaceElements):
1209               case(Points):
1210               case(ContactElementsZero):
1211               case(ReducedContactElementsZero):
1212               case(ContactElementsOne):
1213               case(ReducedContactElementsOne):
1214                   return true;
1215              case(Nodes):
1216              case(DegreesOfFreedom):
1217                 return false;
1218               default:
1219                   stringstream temp;
1220                   temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1221                   throw FinleyAdapterException(temp.str());
1222          }          }
1223          break;          break;
1224       case(Elements):       case(Elements):
1225          if (functionSpaceType_target==Elements) {          if (functionSpaceType_target==Elements) {
1226             return true;             return true;
1227            } else if (functionSpaceType_target==ReducedElements) {
1228               return true;
1229            } else {
1230               return false;
1231            }
1232         case(ReducedElements):
1233            if (functionSpaceType_target==ReducedElements) {
1234               return true;
1235          } else {          } else {
1236             return false;             return false;
1237          }          }
1238       case(FaceElements):       case(FaceElements):
1239          if (functionSpaceType_target==FaceElements) {          if (functionSpaceType_target==FaceElements) {
1240             return true;             return true;
1241            } else if (functionSpaceType_target==ReducedFaceElements) {
1242               return true;
1243            } else {
1244               return false;
1245            }
1246         case(ReducedFaceElements):
1247            if (functionSpaceType_target==ReducedFaceElements) {
1248               return true;
1249          } else {          } else {
1250             return false;             return false;
1251          }          }
# Line 854  bool MeshAdapter::probeInterpolationOnDo Line 1259  bool MeshAdapter::probeInterpolationOnDo
1259       case(ContactElementsOne):       case(ContactElementsOne):
1260          if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {          if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1261             return true;             return true;
1262            } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1263               return true;
1264            } else {
1265               return false;
1266            }
1267         case(ReducedContactElementsZero):
1268         case(ReducedContactElementsOne):
1269            if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1270               return true;
1271          } else {          } else {
1272             return false;             return false;
1273          }          }
# Line 862  bool MeshAdapter::probeInterpolationOnDo Line 1276  bool MeshAdapter::probeInterpolationOnDo
1276             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
1277             case(DegreesOfFreedom):             case(DegreesOfFreedom):
1278             case(Nodes):             case(Nodes):
1279               case(ReducedNodes):
1280             case(Elements):             case(Elements):
1281             case(FaceElements):             case(ReducedElements):
1282             case(Points):             case(Points):
1283               case(FaceElements):
1284               case(ReducedFaceElements):
1285             case(ContactElementsZero):             case(ContactElementsZero):
1286               case(ReducedContactElementsZero):
1287             case(ContactElementsOne):             case(ContactElementsOne):
1288               case(ReducedContactElementsOne):
1289                return true;                return true;
1290             default:             default:
1291               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
1292               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1293                 throw FinleyAdapterException(temp.str());
1294          }          }
1295          break;          break;
1296       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1297         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
1298            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
1299            case(Nodes):            case(ReducedNodes):
1300            case(Elements):            case(Elements):
1301              case(ReducedElements):
1302            case(FaceElements):            case(FaceElements):
1303              case(ReducedFaceElements):
1304            case(Points):            case(Points):
1305            case(ContactElementsZero):            case(ContactElementsZero):
1306              case(ReducedContactElementsZero):
1307            case(ContactElementsOne):            case(ContactElementsOne):
1308              case(ReducedContactElementsOne):
1309                return true;                return true;
1310              case(Nodes):
1311            case(DegreesOfFreedom):            case(DegreesOfFreedom):
1312               return false;               return false;
1313            default:            default:
1314               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
1315               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1316                 throw FinleyAdapterException(temp.str());
1317         }         }
1318         break;         break;
1319       default:       default:
1320          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
1321          sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_source);          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1322            throw FinleyAdapterException(temp.str());
1323          break;          break;
1324    }    }
1325    checkFinleyError();    checkFinleyError();
1326    return false;    return false;
1327  }  }
1328    
1329  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
1330  {  {
1331      return false;      return false;
1332  }  }
1333  bool MeshAdapter::operator==(const MeshAdapter& other) const  
1334    bool MeshAdapter::operator==(const AbstractDomain& other) const
1335    {
1336      const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1337      if (temp!=0) {
1338        return (m_finleyMesh==temp->m_finleyMesh);
1339      } else {
1340        return false;
1341      }
1342    }
1343    
1344    bool MeshAdapter::operator!=(const AbstractDomain& other) const
1345    {
1346      return !(operator==(other));
1347    }
1348    
1349    int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1350    {
1351       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1352       checkPasoError();
1353       return out;
1354    }
1355    
1356    escript::Data MeshAdapter::getX() const
1357    {
1358      return continuousFunction(asAbstractContinuousDomain()).getX();
1359    }
1360    
1361    escript::Data MeshAdapter::getNormal() const
1362    {
1363      return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1364    }
1365    
1366    escript::Data MeshAdapter::getSize() const
1367    {
1368      return function(asAbstractContinuousDomain()).getSize();
1369    }
1370    
1371    int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1372    {
1373      int *out=0,i;
1374      Finley_Mesh* mesh=m_finleyMesh.get();
1375      switch (functionSpaceType) {
1376        case(Nodes):
1377          out=mesh->Nodes->Id;
1378          break;
1379        case(ReducedNodes):
1380          out=mesh->Nodes->reducedNodesId;
1381          break;
1382        case(Elements):
1383          out=mesh->Elements->Id;
1384          break;
1385        case(ReducedElements):
1386          out=mesh->Elements->Id;
1387          break;
1388        case(FaceElements):
1389          out=mesh->FaceElements->Id;
1390          break;
1391        case(ReducedFaceElements):
1392          out=mesh->FaceElements->Id;
1393          break;
1394        case(Points):
1395          out=mesh->Points->Id;
1396          break;
1397        case(ContactElementsZero):
1398          out=mesh->ContactElements->Id;
1399          break;
1400        case(ReducedContactElementsZero):
1401          out=mesh->ContactElements->Id;
1402          break;
1403        case(ContactElementsOne):
1404          out=mesh->ContactElements->Id;
1405          break;
1406        case(ReducedContactElementsOne):
1407          out=mesh->ContactElements->Id;
1408          break;
1409        case(DegreesOfFreedom):
1410          out=mesh->Nodes->degreesOfFreedomId;
1411          break;
1412        case(ReducedDegreesOfFreedom):
1413          out=mesh->Nodes->reducedDegreesOfFreedomId;
1414          break;
1415        default:
1416          stringstream temp;
1417          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1418          throw FinleyAdapterException(temp.str());
1419          break;
1420      }
1421      return out;
1422    }
1423    int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1424  {  {
1425    return (m_finleyMesh==other.m_finleyMesh);    int out=0;
1426      Finley_Mesh* mesh=m_finleyMesh.get();
1427      switch (functionSpaceType) {
1428      case(Nodes):
1429        out=mesh->Nodes->Tag[sampleNo];
1430        break;
1431      case(ReducedNodes):
1432        throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1433        break;
1434      case(Elements):
1435        out=mesh->Elements->Tag[sampleNo];
1436        break;
1437      case(ReducedElements):
1438        out=mesh->Elements->Tag[sampleNo];
1439        break;
1440      case(FaceElements):
1441        out=mesh->FaceElements->Tag[sampleNo];
1442        break;
1443      case(ReducedFaceElements):
1444        out=mesh->FaceElements->Tag[sampleNo];
1445        break;
1446      case(Points):
1447        out=mesh->Points->Tag[sampleNo];
1448        break;
1449      case(ContactElementsZero):
1450        out=mesh->ContactElements->Tag[sampleNo];
1451        break;
1452      case(ReducedContactElementsZero):
1453        out=mesh->ContactElements->Tag[sampleNo];
1454        break;
1455      case(ContactElementsOne):
1456        out=mesh->ContactElements->Tag[sampleNo];
1457        break;
1458      case(ReducedContactElementsOne):
1459        out=mesh->ContactElements->Tag[sampleNo];
1460        break;
1461      case(DegreesOfFreedom):
1462        throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1463        break;
1464      case(ReducedDegreesOfFreedom):
1465        throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1466        break;
1467      default:
1468        stringstream temp;
1469        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1470        throw FinleyAdapterException(temp.str());
1471        break;
1472      }
1473      return out;
1474    }
1475    
1476    
1477    void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1478    {
1479      Finley_Mesh* mesh=m_finleyMesh.get();
1480      escriptDataC tmp=mask.getDataC();
1481      switch(functionSpaceType) {
1482           case(Nodes):
1483              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1484              break;
1485           case(ReducedNodes):
1486              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1487              break;
1488           case(DegreesOfFreedom):
1489              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1490              break;
1491           case(ReducedDegreesOfFreedom):
1492              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1493              break;
1494           case(Elements):
1495              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1496              break;
1497           case(ReducedElements):
1498              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1499              break;
1500           case(FaceElements):
1501              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1502              break;
1503           case(ReducedFaceElements):
1504              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1505              break;
1506           case(Points):
1507              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1508              break;
1509           case(ContactElementsZero):
1510              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1511              break;
1512           case(ReducedContactElementsZero):
1513              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1514              break;
1515           case(ContactElementsOne):
1516              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1517              break;
1518           case(ReducedContactElementsOne):
1519              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1520              break;
1521           default:
1522              stringstream temp;
1523              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1524              throw FinleyAdapterException(temp.str());
1525      }
1526      checkFinleyError();
1527      return;
1528    }
1529    
1530    void MeshAdapter::setTagMap(const std::string& name,  int tag)
1531    {
1532      Finley_Mesh* mesh=m_finleyMesh.get();
1533      Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1534      checkPasoError();
1535      // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1536    }
1537    
1538    int MeshAdapter::getTag(const std::string& name) const
1539    {
1540      Finley_Mesh* mesh=m_finleyMesh.get();
1541      int tag=0;
1542      tag=Finley_Mesh_getTag(mesh, name.c_str());
1543      checkPasoError();
1544      // throwStandardException("MeshAdapter::getTag is not implemented.");
1545      return tag;
1546    }
1547    
1548    bool MeshAdapter::isValidTagName(const std::string& name) const
1549    {
1550      Finley_Mesh* mesh=m_finleyMesh.get();
1551      return Finley_Mesh_isValidTagName(mesh,name.c_str());
1552  }  }
1553    
1554  bool MeshAdapter::operator!=(const MeshAdapter& other) const  std::string MeshAdapter::showTagNames() const
1555  {  {
1556    return !operator==(other);    stringstream temp;
1557      Finley_Mesh* mesh=m_finleyMesh.get();
1558      Finley_TagMap* tag_map=mesh->TagMap;
1559      while (tag_map) {
1560         temp << tag_map->name;
1561         tag_map=tag_map->next;
1562         if (tag_map) temp << ", ";
1563      }
1564      return temp.str();
1565  }  }
1566  // bool MeshAdapter::operator==(const AbstractDomain& other) const  
 // {  
   // try {  
     // const MeshAdapter& ref = dynamic_cast<const MeshAdapter&>(other);  
     // return (operator==(ref));  
   // }  
   // catch (bad_cast) {  
     // return false;  
   // }  
 // }  
1567  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.82  
changed lines
  Added in v.1339

  ViewVC Help
Powered by ViewVC 1.1.26