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

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

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

revision 532 by gross, Wed Feb 15 09:45:53 2006 UTC revision 1326 by ksteube, Mon Oct 1 08:10:41 2007 UTC
# Line 1  Line 1 
 // $Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2  #include "MeshAdapter.h"  /* $Id$ */
3    
4    /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        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 "Data.h"  #include "MeshAdapter.h"
17  #include "DataFactory.h"  #include "escript/Data.h"
18    #include "escript/DataFactory.h"
19    extern "C" {
20    #include "escript/blocktimer.h"
21    }
22    #include <vector>
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
# Line 29  MeshAdapter::FunctionSpaceNamesMapType M Line 32  MeshAdapter::FunctionSpaceNamesMapType M
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  }  }
# Line 58  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
97    {
98      Finley_PrintMesh_Info(m_finleyMesh.get());
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 115  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  {  {
# Line 167  int MeshAdapter::getDiracDeltaFunctionCo Line 221  int MeshAdapter::getDiracDeltaFunctionCo
221  }  }
222    
223  //  //
 // returns a pointer to the tag list of samples of functionSpaceType  
 //  
 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,  
                  int* numTags) const  
 {  
   *tagList=NULL;  
   *numTags=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *tagList=mesh->Nodes->Tag;  
       *numTags=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *tagList=mesh->Elements->Tag;  
       *numTags=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *tagList=mesh->FaceElements->Tag;  
       *numTags=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *tagList=mesh->Points->Tag;  
       *numTags=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*tagList==NULL) {  
     stringstream temp;  
     temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
 }  
   
 //  
 // returns a pointer to the reference no list of samples of functionSpaceType  
 //  
 void MeshAdapter::getReferenceNoList(int functionSpaceType, int** referenceNoList,  
                  int* numReferenceNo) const  
 {  
   *referenceNoList=NULL;  
   *numReferenceNo=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=mesh->Nodes->Id;  
       *numReferenceNo=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *referenceNoList=mesh->Elements->Id;  
       *numReferenceNo=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *referenceNoList=mesh->FaceElements->Id;  
       *numReferenceNo=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *referenceNoList=mesh->Points->Id;  
       *numReferenceNo=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*referenceNoList==NULL) {  
     stringstream temp;  
     temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
 }  
   
 //  
224  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
225  //  //
226  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 332  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 340  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 358  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:
# Line 394  void MeshAdapter::addPDEToSystem( Line 332  void MeshAdapter::addPDEToSystem(
332                       const escript::Data& d, const escript::Data& y,                       const escript::Data& d, const escript::Data& y,
333                       const escript::Data& d_contact,const escript::Data& y_contact) const                       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.getPaso_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();     checkFinleyError();
351     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
352                    mat.getPaso_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
353                    &(rhs.getDataC()),     checkFinleyError();
354                                    &(d.getDataC()),&(y.getDataC()),  
355                                    Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
356     checkFinleyError();     checkFinleyError();
357     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  }
358                    mat.getPaso_SystemMatrix(),  
359                    &(rhs.getDataC()),  void  MeshAdapter::addPDEToLumpedSystem(
360                                    &(d_contact.getDataC()),                       escript::Data& mat,
361                    &(y_contact.getDataC()),                       const escript::Data& D,
362                                    Finley_Assemble_handelShapeMissMatch_Step_out);                       const escript::Data& d) const
363    {
364       escriptDataC _mat=mat.getDataC();
365       escriptDataC _D=D.getDataC();
366       escriptDataC _d=d.getDataC();
367    
368       Finley_Mesh* mesh=m_finleyMesh.get();
369    
370       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
371       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
372    
373     checkFinleyError();     checkFinleyError();
374  }  }
375    
376    
377  //  //
378  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
379  //  //
380  void MeshAdapter::addPDEToRHS( escript::Data& rhs,  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
                      const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  
381  {  {
382     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
383    
384     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
385     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     escriptDataC _X=X.getDataC();
386     checkFinleyError();     escriptDataC _Y=Y.getDataC();
387       escriptDataC _y=y.getDataC();
388       escriptDataC _y_contact=y_contact.getDataC();
389    
390     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
391     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     checkFinleyError();
392    
393       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
394     checkFinleyError();     checkFinleyError();
395     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
396     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
397     checkFinleyError();     checkFinleyError();
398  }  }
399    
# Line 441  void MeshAdapter::interpolateOnDomain(es Line 404  void MeshAdapter::interpolateOnDomain(es
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                 stringstream temp;                 stringstream temp;
# Line 479  void MeshAdapter::interpolateOnDomain(es Line 483  void MeshAdapter::interpolateOnDomain(es
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             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
491          }          }
492          break;          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;
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             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
            break;  
507         }         }
508           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             throw FinleyAdapterException("Error - No interpolation with data on points possible.");             throw FinleyAdapterException("Error - 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             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");             throw FinleyAdapterException("Error - 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               stringstream temp;               stringstream temp;
# Line 537  void MeshAdapter::interpolateOnDomain(es Line 608  void MeshAdapter::interpolateOnDomain(es
608          break;          break;
609       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
610         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
611              case(Nodes):
612                 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):            case(ReducedDegreesOfFreedom):
628               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               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(Nodes):                if (getMPISize()>1) {
664               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");                   escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
665               break;                   escriptDataC _in2 = temp.getDataC();
666            case(DegreesOfFreedom):                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
667               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");                } else {
668                     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
669                 }
670               break;               break;
671            default:            default:
672               stringstream temp;               stringstream temp;
# Line 586  void MeshAdapter::setToX(escript::Data& Line 695  void MeshAdapter::setToX(escript::Data&
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       escript::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    }    }
# Line 605  void MeshAdapter::setToNormal(escript::D Line 716  void MeshAdapter::setToNormal(escript::D
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        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
723        break;        break;
724        case(ReducedNodes):
725          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
726          break;
727      case(Elements):      case(Elements):
728        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
729        break;        break;
730        case(ReducedElements):
731          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
732          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        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");        throw FinleyAdapterException("Error - 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        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
# Line 658  void MeshAdapter::setToIntegrals(std::ve Line 783  void MeshAdapter::setToIntegrals(std::ve
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          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
794            _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          throw FinleyAdapterException("Error - Integral of data on points is not supported.");          throw FinleyAdapterException("Error - 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          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
831            _temp=temp.getDataC();
832            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
833          break;          break;
834       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
835          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
836            _temp=temp.getDataC();
837            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
838          break;          break;
839       default:       default:
840          stringstream temp;          stringstream temp;
# Line 691  void MeshAdapter::setToIntegrals(std::ve Line 843  void MeshAdapter::setToIntegrals(std::ve
843          break;          break;
844    }    }
845    checkFinleyError();    checkFinleyError();
846      blocktimer_increment("integrate()", blocktimer_start);
847  }  }
848    
849  //  //
# Line 706  void MeshAdapter::setToGradient(escript: Line 859  void MeshAdapter::setToGradient(escript:
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            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
881            break;            break;
882           case(ReducedNodes):
883              throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
884              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            throw FinleyAdapterException("Error - Gradient at points is not supported.");            throw FinleyAdapterException("Error - 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            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
# Line 751  void MeshAdapter::setToSize(escript::Dat Line 935  void MeshAdapter::setToSize(escript::Dat
935         case(Nodes):         case(Nodes):
936            throw FinleyAdapterException("Error - Size of nodes is not supported.");            throw FinleyAdapterException("Error - Size of nodes is not supported.");
937            break;            break;
938           case(ReducedNodes):
939              throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
940              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            throw FinleyAdapterException("Error - Size of point elements is not supported.");            throw FinleyAdapterException("Error - Size of point elements is not supported.");
955            break;            break;
# Line 764  void MeshAdapter::setToSize(escript::Dat Line 957  void MeshAdapter::setToSize(escript::Dat
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            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");            throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
966            break;            break;
# Line 783  void MeshAdapter::setToSize(escript::Dat Line 980  void MeshAdapter::setToSize(escript::Dat
980  void MeshAdapter::setNewX(const escript::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_Mesh_setCoordinates(mesh,&(new_x.getDataC()));    tmp = new_x.getDataC();
988      Finley_Mesh_setCoordinates(mesh,&tmp);
989    checkFinleyError();    checkFinleyError();
990  }  }
991    
# Line 795  void MeshAdapter::saveDX(const std::stri Line 994  void MeshAdapter::saveDX(const std::stri
994  {  {
995      int MAX_namelength=256;      int MAX_namelength=256;
996      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
997      char names[num_data][MAX_namelength];    /* win32 refactor */
998      char* c_names[num_data];    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
999      escriptDataC data[num_data];    for(int i=0;i<num_data;i++)
1000      escriptDataC* ptr_data[num_data];    {
1001        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1002      }
1003    
1004      boost::python::list keys=arg.keys();    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1005      for (int i=0;i<num_data;++i) {    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]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1012           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1013               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
1014           data[i]=d.getDataC();           data[i]=d.getDataC();
1015           ptr_data[i]=&(data[i]);           ptr_data[i]=&(data[i]);
          std::string n=boost::python::extract<std::string>(keys[i]);  
1016           c_names[i]=&(names[i][0]);           c_names[i]=&(names[i][0]);
1017           if (n.length()>MAX_namelength-1) {           if (n.length()>MAX_namelength-1) {
1018              strncpy(c_names[i],n.c_str(),MAX_namelength-1);              strncpy(c_names[i],n.c_str(),MAX_namelength-1);
# Line 818  void MeshAdapter::saveDX(const std::stri Line 1023  void MeshAdapter::saveDX(const std::stri
1023      }      }
1024      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1025      checkFinleyError();      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;      return;
1038  }  }
1039    
# Line 826  void MeshAdapter::saveVTK(const std::str Line 1042  void MeshAdapter::saveVTK(const std::str
1042  {  {
1043      int MAX_namelength=256;      int MAX_namelength=256;
1044      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1045      char names[num_data][MAX_namelength];    /* win32 refactor */
1046      char* c_names[num_data];    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1047      escriptDataC data[num_data];    for(int i=0;i<num_data;i++)
1048      escriptDataC* ptr_data[num_data];    {
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();      boost::python::list keys=arg.keys();
1057      for (int i=0;i<num_data;++i) {      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]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1060           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1061               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
1062           data[i]=d.getDataC();           data[i]=d.getDataC();
1063           ptr_data[i]=&(data[i]);           ptr_data[i]=&(data[i]);
          std::string n=boost::python::extract<std::string>(keys[i]);  
1064           c_names[i]=&(names[i][0]);           c_names[i]=&(names[i][0]);
1065           if (n.length()>MAX_namelength-1) {           if (n.length()>MAX_namelength-1) {
1066              strncpy(c_names[i],n.c_str(),MAX_namelength-1);              strncpy(c_names[i],n.c_str(),MAX_namelength-1);
# Line 848  void MeshAdapter::saveVTK(const std::str Line 1070  void MeshAdapter::saveVTK(const std::str
1070           }           }
1071      }      }
1072      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1073      checkFinleyError();  
1074      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;      return;
1086  }  }
1087                                                                                                                                                                                                                                                                                                                                                    
# Line 889  SystemMatrixAdapter MeshAdapter::newSyst Line 1122  SystemMatrixAdapter MeshAdapter::newSyst
1122            
1123      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1124      checkFinleyError();      checkFinleyError();
1125      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);      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();      checkPasoError();
1136      Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);      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    
# Line 915  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:
# Line 933  bool MeshAdapter::probeInterpolationOnDo Line 1179  bool MeshAdapter::probeInterpolationOnDo
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):
1188               case(ReducedFaceElements):
1189               case(Points):
1190               case(ContactElementsZero):
1191               case(ReducedContactElementsZero):
1192               case(ContactElementsOne):
1193               case(ReducedContactElementsOne):
1194                   return true;
1195               default:
1196                   stringstream temp;
1197                   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):             case(FaceElements):
1208               case(ReducedFaceElements):
1209             case(Points):             case(Points):
1210             case(ContactElementsZero):             case(ContactElementsZero):
1211               case(ReducedContactElementsZero):
1212             case(ContactElementsOne):             case(ContactElementsOne):
1213               case(ReducedContactElementsOne):
1214                 return true;                 return true;
1215              case(Nodes):
1216              case(DegreesOfFreedom):
1217                 return false;
1218             default:             default:
1219                 stringstream temp;                 stringstream temp;
1220                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;                 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
# Line 950  bool MeshAdapter::probeInterpolationOnDo Line 1224  bool MeshAdapter::probeInterpolationOnDo
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 969  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 977  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               stringstream temp;               stringstream temp;
# Line 992  bool MeshAdapter::probeInterpolationOnDo Line 1296  bool MeshAdapter::probeInterpolationOnDo
1296       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1297         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
1298            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
1299              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):            case(Nodes):
1311            case(DegreesOfFreedom):            case(DegreesOfFreedom):
# Line 1059  escript::Data MeshAdapter::getSize() con Line 1368  escript::Data MeshAdapter::getSize() con
1368    return function(asAbstractContinuousDomain()).getSize();    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  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1424  {  {
1425    int* tagList;    int out=0;
1426    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1427    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1428    return tagList[sampleNo];    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  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  
1477  {  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1478    int* referenceNoList;  {
1479    int numReferenceNo;    Finley_Mesh* mesh=m_finleyMesh.get();
1480    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    escriptDataC tmp=mask.getDataC();
1481    return referenceNoList[sampleNo];    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    std::string MeshAdapter::showTagNames() const
1555    {
1556      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    
1567  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.532  
changed lines
  Added in v.1326

  ViewVC Help
Powered by ViewVC 1.1.26