/[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

trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp revision 149 by jgs, Thu Sep 1 03:31:39 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 1204 by gross, Sat Jun 23 11:43:12 2007 UTC
# Line 12  Line 12 
12   *                                                                            *   *                                                                            *
13   ******************************************************************************   ******************************************************************************
14  */  */
15  extern "C" {  
16  #include "finley/finleyC/Finley.h"  #ifdef PASO_MPI
17  #include "finley/finleyC/Assemble.h"  #include <mpi.h>
18  #include "finley/finleyC/Mesh.h"  #endif
19  #include "finley/finleyC/Finley.h"  #include "MeshAdapter.h"
20  #include "finley/finleyC/System.h"  
21  }  #include "escript/Data.h"
22  #include "finley/CPPAdapter/SystemMatrixAdapter.h"  #include "escript/DataFactory.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>  
 #include <vector>  
 #include <sstream>  
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
# Line 44  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  {  {
# Line 85  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 78  Finley_Mesh* MeshAdapter::getFinley_Mesh
78    
79  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
80  {  {
81    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
82    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
83    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
84    checkFinleyError();    checkFinleyError();
85      TMPMEMFREE(fName);
86  }  }
87    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
88  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
89  {  {
90    return "FinleyMesh";    return string("FinleyMesh");
91  }  }
92    
93  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const  string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
# Line 110  string MeshAdapter::functionSpaceTypeAsS Line 97  string MeshAdapter::functionSpaceTypeAsS
97    if (loc==m_functionSpaceTypeNames.end()) {    if (loc==m_functionSpaceTypeNames.end()) {
98      return "Invalid function space type code.";      return "Invalid function space type code.";
99    } else {    } else {
100      return loc->second;      return string(loc->second);
101    }    }
102  }  }
103    
# Line 130  void MeshAdapter::setFunctionSpaceTypeNa Line 117  void MeshAdapter::setFunctionSpaceTypeNa
117    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
118      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
119    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
120        (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
121      m_functionSpaceTypeNames.insert
122      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
123    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
124        (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
125      m_functionSpaceTypeNames.insert
126      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
127    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
128        (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
129      m_functionSpaceTypeNames.insert
130      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
131    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
132      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
133    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
134        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
135      m_functionSpaceTypeNames.insert
136      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
137      m_functionSpaceTypeNames.insert
138        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
139  }  }
140    
141  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
142  {  {
143    return Nodes;    return Nodes;
144  }  }
145    int MeshAdapter::getReducedContinuousFunctionCode() const
146    {
147      return ReducedNodes;
148    }
149    
150  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
151  {  {
152    return Elements;    return Elements;
153  }  }
154    int MeshAdapter::getReducedFunctionCode() const
155    {
156      return ReducedElements;
157    }
158    
159  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
160  {  {
161    return FaceElements;    return FaceElements;
162  }  }
163    int MeshAdapter::getReducedFunctionOnBoundaryCode() const
164    {
165      return ReducedFaceElements;
166    }
167    
168  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
169  {  {
170    return ContactElementsZero;    return ContactElementsZero;
171  }  }
172    int MeshAdapter::getReducedFunctionOnContactZeroCode() const
173    {
174      return ReducedContactElementsZero;
175    }
176    
177  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
178  {  {
179    return ContactElementsOne;    return ContactElementsOne;
180  }  }
181    int MeshAdapter::getReducedFunctionOnContactOneCode() const
182    {
183      return ReducedContactElementsOne;
184    }
185    
186  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
187  {  {
# Line 182  int MeshAdapter::getDiracDeltaFunctionCo Line 199  int MeshAdapter::getDiracDeltaFunctionCo
199  }  }
200    
201  //  //
 // 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;  
 }  
   
 //  
202  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
203  //  //
204  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 351  pair<int,int> MeshAdapter::getDataShape( Line 222  pair<int,int> MeshAdapter::getDataShape(
222             numDataPointsPerSample=1;             numDataPointsPerSample=1;
223             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
224             break;             break;
225          case(ReducedNodes):
226               /* TODO: add ReducedNodes */
227               numDataPointsPerSample=1;
228               if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
229               throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");
230               break;
231        case(Elements):        case(Elements):
232             if (mesh->Elements!=NULL) {             if (mesh->Elements!=NULL) {
233               numSamples=mesh->Elements->numElements;               numSamples=mesh->Elements->numElements;
234               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
235             }             }
236             break;             break;
237          case(ReducedElements):
238               if (mesh->Elements!=NULL) {
239                 numSamples=mesh->Elements->numElements;
240                 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
241               }
242               break;
243        case(FaceElements):        case(FaceElements):
244             if (mesh->FaceElements!=NULL) {             if (mesh->FaceElements!=NULL) {
245                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
246                  numSamples=mesh->FaceElements->numElements;                  numSamples=mesh->FaceElements->numElements;
247             }             }
248             break;             break;
249          case(ReducedFaceElements):
250               if (mesh->FaceElements!=NULL) {
251                    numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
252                    numSamples=mesh->FaceElements->numElements;
253               }
254               break;
255        case(Points):        case(Points):
256             if (mesh->Points!=NULL) {             if (mesh->Points!=NULL) {
257               numDataPointsPerSample=1;               numDataPointsPerSample=1;
# Line 375  pair<int,int> MeshAdapter::getDataShape( Line 264  pair<int,int> MeshAdapter::getDataShape(
264               numSamples=mesh->ContactElements->numElements;               numSamples=mesh->ContactElements->numElements;
265             }             }
266             break;             break;
267          case(ReducedContactElementsZero):
268               if (mesh->ContactElements!=NULL) {
269                 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
270                 numSamples=mesh->ContactElements->numElements;
271               }
272               break;
273        case(ContactElementsOne):        case(ContactElementsOne):
274             if (mesh->ContactElements!=NULL) {             if (mesh->ContactElements!=NULL) {
275               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
276               numSamples=mesh->ContactElements->numElements;               numSamples=mesh->ContactElements->numElements;
277             }             }
278             break;             break;
279          case(ReducedContactElementsOne):
280               if (mesh->ContactElements!=NULL) {
281                 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
282                 numSamples=mesh->ContactElements->numElements;
283               }
284               break;
285        case(DegreesOfFreedom):        case(DegreesOfFreedom):
286             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
287               numDataPointsPerSample=1;               numDataPointsPerSample=1;
288    #ifndef PASO_MPI
289               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
290    #else
291                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
292    #endif
293             }             }
294             break;             break;
295        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
296             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
297               numDataPointsPerSample=1;               numDataPointsPerSample=1;
298    #ifndef PASO_MPI
299               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
300    #else
301                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
302    #endif
303             }             }
304             break;             break;
305        default:        default:
306          stringstream temp;          stringstream temp;
307          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
308          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
309          break;          break;
310        }        }
# Line 407  pair<int,int> MeshAdapter::getDataShape( Line 315  pair<int,int> MeshAdapter::getDataShape(
315  // 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:
316  //  //
317  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
318                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
319                       const Data& A, const Data& B, const Data& C,const  Data& D,const  Data& X,const  Data& Y,                       const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
320                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
321                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
322  {  {
323       escriptDataC _rhs=rhs.getDataC();
324       escriptDataC _A  =A.getDataC();
325       escriptDataC _B=B.getDataC();
326       escriptDataC _C=C.getDataC();
327       escriptDataC _D=D.getDataC();
328       escriptDataC _X=X.getDataC();
329       escriptDataC _Y=Y.getDataC();
330       escriptDataC _d=d.getDataC();
331       escriptDataC _y=y.getDataC();
332       escriptDataC _d_contact=d_contact.getDataC();
333       escriptDataC _y_contact=y_contact.getDataC();
334    
335     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
336     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),  
337                         &(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 );
338     checkFinleyError();     checkFinleyError();
339     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
340                    mat.getFinley_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
341                    &(rhs.getDataC()),     checkFinleyError();
342                                    &(d.getDataC()),&(y.getDataC()),  
343                                    Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
344     checkFinleyError();     checkFinleyError();
345     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  }
346                    mat.getFinley_SystemMatrix(),  
347                    &(rhs.getDataC()),  void  MeshAdapter::addPDEToLumpedSystem(
348                                    &(d_contact.getDataC()),                       escript::Data& mat,
349                    &(y_contact.getDataC()),                       const escript::Data& D,
350                                    Finley_Assemble_handelShapeMissMatch_Step_out);                       const escript::Data& d) const
351    {
352       escriptDataC _mat=mat.getDataC();
353       escriptDataC _D=D.getDataC();
354       escriptDataC _d=d.getDataC();
355    
356       Finley_Mesh* mesh=m_finleyMesh.get();
357      
358       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
359       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
360    
361     checkFinleyError();     checkFinleyError();
362  }  }
363    
364    
365  //  //
366  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
367  //  //
368  void MeshAdapter::addPDEToRHS( Data& rhs,  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
                      const  Data& X,const  Data& Y, const Data& y, const Data& y_contact) const  
369  {  {
370     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
371    
372     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
373     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     escriptDataC _X=X.getDataC();
374     checkFinleyError();     escriptDataC _Y=Y.getDataC();
375       escriptDataC _y=y.getDataC();
376       escriptDataC _y_contact=y_contact.getDataC();
377    
378     // 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 );
379     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     checkFinleyError();
380    
381       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
382     checkFinleyError();     checkFinleyError();
383     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
384     // 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 );
385     checkFinleyError();     checkFinleyError();
386  }  }
387    
388  //  //
389  // interpolates data between different function spaces:  // interpolates data between different function spaces:
390  //  //
391  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
392  {  {
393    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
394    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
395    if (inDomain!=*this)    if (inDomain!=*this)  
396       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
397    if (targetDomain!=*this)    if (targetDomain!=*this)
398       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
399    
400    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
401      escriptDataC _target=target.getDataC();
402      escriptDataC _in=in.getDataC();
403    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
404       case(Nodes):       case(Nodes):
405          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
406             case(Nodes):             case(Nodes):
407               case(ReducedNodes):
408               case(DegreesOfFreedom):
409             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
410                   Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
411                   break;
412               case(Elements):
413               case(ReducedElements):
414                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
415                   break;
416               case(FaceElements):
417               case(ReducedFaceElements):
418                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
419                   break;
420               case(Points):
421                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
422                   break;
423               case(ContactElementsZero):
424               case(ReducedContactElementsZero):
425                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
426                   break;
427               case(ContactElementsOne):
428               case(ReducedContactElementsOne):
429                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
430                   break;
431               default:
432                   stringstream temp;
433                   temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
434                   throw FinleyAdapterException(temp.str());
435                   break;
436            }
437            break;
438         case(ReducedNodes):
439            switch(target.getFunctionSpace().getTypeCode()) {
440               case(Nodes):
441               case(ReducedNodes):
442             case(DegreesOfFreedom):             case(DegreesOfFreedom):
443                 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));             case(ReducedDegreesOfFreedom):
444                   Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
445                 break;                 break;
446             case(Elements):             case(Elements):
447                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));             case(ReducedElements):
448                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
449                 break;                 break;
450             case(FaceElements):             case(FaceElements):
451                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedFaceElements):
452                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
453                 break;                 break;
454             case(Points):             case(Points):
455                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
456                 break;                 break;
457             case(ContactElementsZero):             case(ContactElementsZero):
458                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsZero):
459                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
460                 break;                 break;
461             case(ContactElementsOne):             case(ContactElementsOne):
462                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsOne):
463                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
464                 break;                 break;
465             default:             default:
466                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
467                 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();
468                   throw FinleyAdapterException(temp.str());
469                 break;                 break;
470          }          }
471          break;          break;
472       case(Elements):       case(Elements):
473          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
474             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
475            } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
476               Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
477            } else {
478               throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
479            }
480            break;
481         case(ReducedElements):
482            if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
483               Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
484          } else {          } else {
485             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
486          }          }
487          break;          break;
488       case(FaceElements):       case(FaceElements):
489          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
490             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
491            } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
492               Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
493          } else {          } else {
494             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
495             sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");         }
496             break;         break;
497         case(ReducedFaceElements):
498            if (target.getFunctionSpace().getTypeCode()==FaceElements) {
499               Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
500            } else {
501               throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
502         }         }
503           break;
504       case(Points):       case(Points):
505          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
506             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
507          } else {          } else {
508             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;  
509          }          }
510          break;          break;
511       case(ContactElementsZero):       case(ContactElementsZero):
512       case(ContactElementsOne):       case(ContactElementsOne):
513          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
514             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
515            } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
516               Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
517          } else {          } else {
518             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.");  
519             break;             break;
520          }          }
521          break;          break;
522       case(DegreesOfFreedom):       case(ReducedContactElementsZero):
523         case(ReducedContactElementsOne):
524            if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
525               Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
526            } else {
527               throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
528               break;
529            }
530            break;
531         case(DegreesOfFreedom):      
532          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
533             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
534             case(DegreesOfFreedom):             case(DegreesOfFreedom):
535             case(Nodes):             case(Nodes):
536                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));             case(ReducedNodes):
537                  Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
538                  break;
539    #ifndef PASO_MPI
540               case(Elements):
541               case(ReducedElements):
542                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
543                  break;
544               case(FaceElements):
545               case(ReducedFaceElements):
546                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
547                break;                break;
548               case(Points):
549                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
550                  break;
551               case(ContactElementsZero):
552               case(ContactElementsOne):
553               case(ReducedContactElementsZero):
554               case(ReducedContactElementsOne):
555                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
556                 break;
557    #else
558               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
559             case(Elements):             case(Elements):
560                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));             case(ReducedElements):
561                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
562                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);
563                break;                break;
564             case(FaceElements):             case(FaceElements):
565                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedFaceElements):
566                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
567                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);
568                break;                break;
569             case(Points):             case(Points):
570                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
571                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);
572                break;                break;
573             case(ContactElementsZero):             case(ContactElementsZero):
574             case(ContactElementsOne):             case(ContactElementsOne):
575                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsZero):
576               case(ReducedContactElementsOne):
577                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
578                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);
579               break;               break;
580    #endif
581             default:             default:
582               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
583               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();
584                 throw FinleyAdapterException(temp.str());
585               break;               break;
586          }          }
587          break;          break;
588       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
589         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
590              case(Nodes):
591                 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
592                 break;
593              case(ReducedNodes):
594                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
595                 break;
596              case(DegreesOfFreedom):
597                 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
598                 break;
599            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
600               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
601               break;               break;
602            case(Elements):            case(Elements):
603               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));            case(ReducedElements):
604                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
605               break;               break;
606            case(FaceElements):            case(FaceElements):
607               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));            case(ReducedFaceElements):
608                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
609               break;               break;
610            case(Points):            case(Points):
611               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
612               break;               break;
613            case(ContactElementsZero):            case(ContactElementsZero):
614            case(ContactElementsOne):            case(ContactElementsOne):
615               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));            case(ReducedContactElementsZero):
616               break;            case(ReducedContactElementsOne):
617            case(Nodes):               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
              Finley_ErrorCode=TYPE_ERROR;  
              sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");  
              break;  
           case(DegreesOfFreedom):  
              Finley_ErrorCode=TYPE_ERROR;  
              sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
618               break;               break;
619            default:            default:
620               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
621               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();
622                 throw FinleyAdapterException(temp.str());
623               break;               break;
624         }         }
625         break;         break;
626       default:       default:
627          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
628          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();
629            throw FinleyAdapterException(temp.str());
630          break;          break;
631    }    }
632    checkFinleyError();    checkFinleyError();
# Line 598  void MeshAdapter::interpolateOnDomain(Da Line 635  void MeshAdapter::interpolateOnDomain(Da
635  //  //
636  // copies the locations of sample points into x:  // copies the locations of sample points into x:
637  //  //
638  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
639  {  {
640    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
641    if (argDomain!=*this)    if (argDomain!=*this)
# Line 606  void MeshAdapter::setToX(Data& arg) cons Line 643  void MeshAdapter::setToX(Data& arg) cons
643    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
644    // in case of values node coordinates we can do the job directly:    // in case of values node coordinates we can do the job directly:
645    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
646       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       escriptDataC _arg=arg.getDataC();
647         Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
648    } else {    } else {
649       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
650       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       escriptDataC _tmp_data=tmp_data.getDataC();
651         Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
652       // this is then interpolated onto arg:       // this is then interpolated onto arg:
653       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
654    }    }
# Line 619  void MeshAdapter::setToX(Data& arg) cons Line 658  void MeshAdapter::setToX(Data& arg) cons
658  //  //
659  // 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:
660  //  //
661  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
662  {  {
663    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
664    if (normalDomain!=*this)    if (normalDomain!=*this)
665       throw FinleyAdapterException("Error - Illegal domain of normal locations");       throw FinleyAdapterException("Error - Illegal domain of normal locations");
666    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
667      escriptDataC _normal=normal.getDataC();
668    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
669      case(Nodes):      case(Nodes):
670        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
671        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");        break;
672        case(ReducedNodes):
673          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
674        break;        break;
675      case(Elements):      case(Elements):
676        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
677        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");        break;
678        case(ReducedElements):
679          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
680        break;        break;
681      case (FaceElements):      case (FaceElements):
682        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
683          break;
684        case (ReducedFaceElements):
685          Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
686        break;        break;
687      case(Points):      case(Points):
688        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");  
689        break;        break;
690      case (ContactElementsOne):      case (ContactElementsOne):
691      case (ContactElementsZero):      case (ContactElementsZero):
692        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
693          break;
694        case (ReducedContactElementsOne):
695        case (ReducedContactElementsZero):
696          Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
697        break;        break;
698      case(DegreesOfFreedom):      case(DegreesOfFreedom):
699        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.");  
700        break;        break;
701      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
702        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.");  
703        break;        break;
704      default:      default:
705        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
706        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();
707          throw FinleyAdapterException(temp.str());
708        break;        break;
709    }    }
710    checkFinleyError();    checkFinleyError();
# Line 664  void MeshAdapter::setToNormal(Data& norm Line 713  void MeshAdapter::setToNormal(Data& norm
713  //  //
714  // interpolates data to other domain:  // interpolates data to other domain:
715  //  //
716  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
717  {  {
718    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
719    if (targetDomain!=*this)    if (targetDomain!=*this)
720       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
721    
722    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();  
723  }  }
724    
725  //  //
726  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
727  //  //
728  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
729  {  {
730    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
731    if (argDomain!=*this)    if (argDomain!=*this)
732       throw FinleyAdapterException("Error - Illegal domain of integration kernel");       throw FinleyAdapterException("Error - Illegal domain of integration kernel");
733    
734    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
735      escriptDataC _arg=arg.getDataC();
736    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
737       case(Nodes):       case(Nodes):
738          Finley_ErrorCode=TYPE_ERROR;          /* TODO */
739          sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
740            break;
741         case(ReducedNodes):
742            /* TODO */
743            throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");
744          break;          break;
745       case(Elements):       case(Elements):
746          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
747            break;
748         case(ReducedElements):
749            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
750          break;          break;
751       case(FaceElements):       case(FaceElements):
752          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
753            break;
754         case(ReducedFaceElements):
755            Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
756          break;          break;
757       case(Points):       case(Points):
758          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.");  
759          break;          break;
760       case(ContactElementsZero):       case(ContactElementsZero):
761          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
762            break;
763         case(ReducedContactElementsZero):
764            Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
765          break;          break;
766       case(ContactElementsOne):       case(ContactElementsOne):
767          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
768            break;
769         case(ReducedContactElementsOne):
770            Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
771          break;          break;
772       case(DegreesOfFreedom):       case(DegreesOfFreedom):
773          Finley_ErrorCode=TYPE_ERROR;          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
         sprintf(Finley_ErrorMsg,"Integral of data on degrees of freedom is not supported.");  
774          break;          break;
775       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
776          Finley_ErrorCode=TYPE_ERROR;          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
         sprintf(Finley_ErrorMsg,"Integral of data on reduced degrees of freedom is not supported.");  
777          break;          break;
778       default:       default:
779          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
780          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();
781            throw FinleyAdapterException(temp.str());
782          break;          break;
783    }    }
784    checkFinleyError();    checkFinleyError();
# Line 725  void MeshAdapter::setToIntegrals(std::ve Line 787  void MeshAdapter::setToIntegrals(std::ve
787  //  //
788  // calculates the gradient of arg:  // calculates the gradient of arg:
789  //  //
790  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
791  {  {
792    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
793    if (argDomain!=*this)    if (argDomain!=*this)
# Line 735  void MeshAdapter::setToGradient(Data& gr Line 797  void MeshAdapter::setToGradient(Data& gr
797       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
798    
799    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
800      escriptDataC _grad=grad.getDataC();
801      escriptDataC nodeDataC;
802    #ifdef PASO_MPI
803      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
804      if( arg.getFunctionSpace().getTypeCode() != Nodes )
805      {
806        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
807        nodeDataC = nodeTemp.getDataC();
808      }
809      else
810        nodeDataC = arg.getDataC();
811    #else
812      nodeDataC = arg.getDataC();
813    #endif
814    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
815         case(Nodes):         case(Nodes):
816            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
817            sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");            break;
818           case(ReducedNodes):
819              throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
820            break;            break;
821         case(Elements):         case(Elements):
822            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
823              break;
824           case(ReducedElements):
825              Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
826            break;            break;
827         case(FaceElements):         case(FaceElements):
828            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
829              break;
830           case(ReducedFaceElements):
831              Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
832            break;            break;
833         case(Points):         case(Points):
834            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
835            break;            break;
836         case(ContactElementsZero):         case(ContactElementsZero):
837            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
838              break;
839           case(ReducedContactElementsZero):
840              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
841            break;            break;
842         case(ContactElementsOne):         case(ContactElementsOne):
843            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
844              break;
845           case(ReducedContactElementsOne):
846              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
847            break;            break;
848         case(DegreesOfFreedom):         case(DegreesOfFreedom):
849            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.");  
850            break;            break;
851         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
852            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.");  
853            break;            break;
854         default:         default:
855            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
856            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();
857              throw FinleyAdapterException(temp.str());
858            break;            break;
859    }    }
860    checkFinleyError();    checkFinleyError();
# Line 775  void MeshAdapter::setToGradient(Data& gr Line 863  void MeshAdapter::setToGradient(Data& gr
863  //  //
864  // returns the size of elements:  // returns the size of elements:
865  //  //
866  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
867  {  {
868    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
869    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
870    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
871         case(Nodes):         case(Nodes):
872            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
873            sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");            break;
874           case(ReducedNodes):
875              throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
876            break;            break;
877         case(Elements):         case(Elements):
878            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
879            break;            break;
880           case(ReducedElements):
881              Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
882              break;
883         case(FaceElements):         case(FaceElements):
884            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
885            break;            break;
886           case(ReducedFaceElements):
887              Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
888              break;
889         case(Points):         case(Points):
890            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
891            break;            break;
892         case(ContactElementsZero):         case(ContactElementsZero):
893         case(ContactElementsOne):         case(ContactElementsOne):
894            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
895            break;            break;
896           case(ReducedContactElementsZero):
897           case(ReducedContactElementsOne):
898              Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
899              break;
900         case(DegreesOfFreedom):         case(DegreesOfFreedom):
901            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.");  
902            break;            break;
903         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
904            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.");  
905            break;            break;
906         default:         default:
907            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
908            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();
909              throw FinleyAdapterException(temp.str());
910            break;            break;
911    }    }
912    checkFinleyError();    checkFinleyError();
913  }  }
914    
915  // sets the location of nodes:  // sets the location of nodes:
916  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
917  {  {
918    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
919      escriptDataC tmp;
920    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
921    if (newDomain!=*this)    if (newDomain!=*this)
922       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
923    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
924      Finley_Mesh_setCoordinates(mesh,&tmp);
925    checkFinleyError();    checkFinleyError();
926  }  }
927    
928  // saves a data array in openDX format:  // saves a data array in openDX format:
929  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
930  {  {
931    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
932    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
933      /* win32 refactor */
934      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
935      for(int i=0;i<num_data;i++)
936      {
937        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
938      }
939    
940      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
941      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
942      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
943    
944        boost::python::list keys=arg.keys();
945        for (int i=0;i<num_data;++i) {
946             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
947             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
948                 throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
949             data[i]=d.getDataC();
950             ptr_data[i]=&(data[i]);
951             std::string n=boost::python::extract<std::string>(keys[i]);
952             c_names[i]=&(names[i][0]);
953             if (n.length()>MAX_namelength-1) {
954                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
955                c_names[i][MAX_namelength-1]='\0';
956             } else {
957                strcpy(c_names[i],n.c_str());
958             }
959        }
960        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
961        checkFinleyError();
962        
963          /* win32 refactor */
964      TMPMEMFREE(c_names);
965      TMPMEMFREE(data);
966      TMPMEMFREE(ptr_data);
967      for(int i=0;i<num_data;i++)
968      {
969        TMPMEMFREE(names[i]);
970      }
971      TMPMEMFREE(names);
972    
973        return;
974  }  }
975    
976  // saves a data array in openVTK format:  // saves a data array in openVTK format:
977  void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
978  {  {
979    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
980    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
981  }    /* win32 refactor */
982      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
983      for(int i=0;i<num_data;i++)
984      {
985        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
986      }
987    
988      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
989      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
990      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
991    
992        boost::python::list keys=arg.keys();
993        for (int i=0;i<num_data;++i) {
994             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
995             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
996                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
997             data[i]=d.getDataC();
998             ptr_data[i]=&(data[i]);
999             std::string n=boost::python::extract<std::string>(keys[i]);
1000             c_names[i]=&(names[i][0]);
1001             if (n.length()>MAX_namelength-1) {
1002                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1003                c_names[i][MAX_namelength-1]='\0';
1004             } else {
1005                strcpy(c_names[i],n.c_str());
1006             }
1007        }
1008    #ifndef PASO_MPI    
1009        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1010    #else
1011        Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1012    #endif
1013    
1014    checkFinleyError();
1015      /* win32 refactor */
1016      TMPMEMFREE(c_names);
1017      TMPMEMFREE(data);
1018      TMPMEMFREE(ptr_data);
1019      for(int i=0;i<num_data;i++)
1020      {
1021        TMPMEMFREE(names[i]);
1022      }
1023      TMPMEMFREE(names);
1024    
1025        return;
1026    }
1027                                                                                                                                                                            
1028                                                                                                                                                                            
1029  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1030  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1031                        const int row_blocksize,                        const int row_blocksize,
# Line 873  SystemMatrixAdapter MeshAdapter::newSyst Line 1060  SystemMatrixAdapter MeshAdapter::newSyst
1060      }      }
1061      // generate matrix:      // generate matrix:
1062            
1063      Finley_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
     checkFinleyError();  
     Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);  
1064      checkFinleyError();      checkFinleyError();
1065      Finley_SystemMatrixPattern_dealloc(fsystemMatrixPattern);      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1066        checkPasoError();
1067        Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
1068      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1069  }  }
1070    
# Line 901  bool MeshAdapter::isCellOriented(int fun Line 1088  bool MeshAdapter::isCellOriented(int fun
1088         case(Points):         case(Points):
1089         case(ContactElementsZero):         case(ContactElementsZero):
1090         case(ContactElementsOne):         case(ContactElementsOne):
1091           case(ReducedElements):
1092           case(ReducedFaceElements):
1093           case(ReducedContactElementsZero):
1094           case(ReducedContactElementsOne):
1095            return true;            return true;
1096            break;            break;
1097         default:         default:
1098            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
1099            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;
1100              throw FinleyAdapterException(temp.str());
1101            break;            break;
1102    }    }
1103    checkFinleyError();    checkFinleyError();
# Line 918  bool MeshAdapter::probeInterpolationOnDo Line 1110  bool MeshAdapter::probeInterpolationOnDo
1110       case(Nodes):       case(Nodes):
1111          switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1112             case(Nodes):             case(Nodes):
1113               case(ReducedNodes):
1114             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
1115             case(DegreesOfFreedom):             case(DegreesOfFreedom):
1116             case(Elements):             case(Elements):
1117               case(ReducedElements):
1118             case(FaceElements):             case(FaceElements):
1119               case(ReducedFaceElements):
1120             case(Points):             case(Points):
1121             case(ContactElementsZero):             case(ContactElementsZero):
1122               case(ReducedContactElementsZero):
1123             case(ContactElementsOne):             case(ContactElementsOne):
1124               case(ReducedContactElementsOne):
1125                 return true;                 return true;
1126             default:             default:
1127                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
1128                 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;
1129                   throw FinleyAdapterException(temp.str());
1130            }
1131            break;
1132         case(ReducedNodes):
1133            switch(functionSpaceType_target) {
1134               case(ReducedNodes):
1135               case(ReducedDegreesOfFreedom):
1136               case(Elements):
1137               case(ReducedElements):
1138               case(FaceElements):
1139               case(ReducedFaceElements):
1140               case(Points):
1141               case(ContactElementsZero):
1142               case(ReducedContactElementsZero):
1143               case(ContactElementsOne):
1144               case(ReducedContactElementsOne):
1145                   return true;
1146              case(Nodes):
1147              case(DegreesOfFreedom):
1148                 return false;
1149               default:
1150                   stringstream temp;
1151                   temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1152                   throw FinleyAdapterException(temp.str());
1153          }          }
1154          break;          break;
1155       case(Elements):       case(Elements):
1156          if (functionSpaceType_target==Elements) {          if (functionSpaceType_target==Elements) {
1157             return true;             return true;
1158            } else if (functionSpaceType_target==ReducedElements) {
1159               return true;
1160            } else {
1161               return false;
1162            }
1163         case(ReducedElements):
1164            if (functionSpaceType_target==ReducedElements) {
1165               return true;
1166          } else {          } else {
1167             return false;             return false;
1168          }          }
1169       case(FaceElements):       case(FaceElements):
1170          if (functionSpaceType_target==FaceElements) {          if (functionSpaceType_target==FaceElements) {
1171             return true;             return true;
1172            } else if (functionSpaceType_target==ReducedFaceElements) {
1173               return true;
1174            } else {
1175               return false;
1176            }
1177         case(ReducedFaceElements):
1178            if (functionSpaceType_target==ReducedFaceElements) {
1179               return true;
1180          } else {          } else {
1181             return false;             return false;
1182          }          }
# Line 953  bool MeshAdapter::probeInterpolationOnDo Line 1190  bool MeshAdapter::probeInterpolationOnDo
1190       case(ContactElementsOne):       case(ContactElementsOne):
1191          if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {          if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1192             return true;             return true;
1193            } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1194               return true;
1195            } else {
1196               return false;
1197            }
1198         case(ReducedContactElementsZero):
1199         case(ReducedContactElementsOne):
1200            if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1201               return true;
1202          } else {          } else {
1203             return false;             return false;
1204          }          }
# Line 961  bool MeshAdapter::probeInterpolationOnDo Line 1207  bool MeshAdapter::probeInterpolationOnDo
1207             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
1208             case(DegreesOfFreedom):             case(DegreesOfFreedom):
1209             case(Nodes):             case(Nodes):
1210               case(ReducedNodes):
1211             case(Elements):             case(Elements):
1212             case(FaceElements):             case(ReducedElements):
1213             case(Points):             case(Points):
1214               case(FaceElements):
1215               case(ReducedFaceElements):
1216             case(ContactElementsZero):             case(ContactElementsZero):
1217               case(ReducedContactElementsZero):
1218             case(ContactElementsOne):             case(ContactElementsOne):
1219               case(ReducedContactElementsOne):
1220                return true;                return true;
1221             default:             default:
1222               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
1223               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;
1224                 throw FinleyAdapterException(temp.str());
1225          }          }
1226          break;          break;
1227       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1228         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
1229            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
1230              case(ReducedNodes):
1231            case(Elements):            case(Elements):
1232              case(ReducedElements):
1233            case(FaceElements):            case(FaceElements):
1234              case(ReducedFaceElements):
1235            case(Points):            case(Points):
1236            case(ContactElementsZero):            case(ContactElementsZero):
1237              case(ReducedContactElementsZero):
1238            case(ContactElementsOne):            case(ContactElementsOne):
1239              case(ReducedContactElementsOne):
1240                return true;                return true;
1241            case(Nodes):            case(Nodes):
1242            case(DegreesOfFreedom):            case(DegreesOfFreedom):
1243               return false;               return false;
1244            default:            default:
1245               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
1246               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;
1247                 throw FinleyAdapterException(temp.str());
1248         }         }
1249         break;         break;
1250       default:       default:
1251          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
1252          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;
1253            throw FinleyAdapterException(temp.str());
1254          break;          break;
1255    }    }
1256    checkFinleyError();    checkFinleyError();
# Line 1018  bool MeshAdapter::operator!=(const Abstr Line 1277  bool MeshAdapter::operator!=(const Abstr
1277    return !(operator==(other));    return !(operator==(other));
1278  }  }
1279    
1280  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1281  {  {
1282     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1283     checkFinleyError();     checkPasoError();
1284     return out;     return out;
1285  }  }
1286    
1287  Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
1288  {  {
1289    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
1290  }  }
1291    
1292  Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
1293  {  {
1294    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1295  }  }
1296    
1297  Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
1298  {  {
1299    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
1300  }  }
1301    
1302    int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1303    {
1304      int *out=0,i;
1305      Finley_Mesh* mesh=m_finleyMesh.get();
1306      switch (functionSpaceType) {
1307        case(Nodes):
1308          if (mesh->Nodes!=NULL) {
1309            out=mesh->Nodes->Id;
1310            break;
1311          }
1312        case(ReducedNodes):
1313          throw FinleyAdapterException("Error -  ReducedNodes not supported yet.");
1314          break;
1315        case(Elements):
1316          out=mesh->Elements->Id;
1317          break;
1318        case(ReducedElements):
1319          out=mesh->Elements->Id;
1320          break;
1321        case(FaceElements):
1322          out=mesh->FaceElements->Id;
1323          break;
1324        case(ReducedFaceElements):
1325          out=mesh->FaceElements->Id;
1326          break;
1327        case(Points):
1328          out=mesh->Points->Id;
1329          break;
1330        case(ContactElementsZero):
1331          out=mesh->ContactElements->Id;
1332          break;
1333        case(ReducedContactElementsZero):
1334          out=mesh->ContactElements->Id;
1335          break;
1336        case(ContactElementsOne):
1337          out=mesh->ContactElements->Id;
1338          break;
1339        case(ReducedContactElementsOne):
1340          out=mesh->ContactElements->Id;
1341          break;
1342        case(DegreesOfFreedom):
1343          out=mesh->Nodes->degreeOfFreedomId;
1344          break;
1345        case(ReducedDegreesOfFreedom):
1346          out=mesh->Nodes->reducedDegreeOfFreedomId;
1347          break;
1348        default:
1349          stringstream temp;
1350          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1351          throw FinleyAdapterException(temp.str());
1352          break;
1353      }
1354      return out;
1355    }
1356  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1357  {  {
1358    int* tagList;    int out=0;
1359    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1360    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1361    return tagList[sampleNo];    case(Nodes):
1362        out=mesh->Nodes->Tag[sampleNo];
1363        break;
1364      case(ReducedNodes):
1365        throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1366        break;
1367      case(Elements):
1368        out=mesh->Elements->Tag[sampleNo];
1369        break;
1370      case(ReducedElements):
1371        out=mesh->Elements->Tag[sampleNo];
1372        break;
1373      case(FaceElements):
1374        out=mesh->FaceElements->Tag[sampleNo];
1375        break;
1376      case(ReducedFaceElements):
1377        out=mesh->FaceElements->Tag[sampleNo];
1378        break;
1379      case(Points):
1380        out=mesh->Points->Tag[sampleNo];
1381        break;
1382      case(ContactElementsZero):
1383        out=mesh->ContactElements->Tag[sampleNo];
1384        break;
1385      case(ReducedContactElementsZero):
1386        out=mesh->ContactElements->Tag[sampleNo];
1387        break;
1388      case(ContactElementsOne):
1389        out=mesh->ContactElements->Tag[sampleNo];
1390        break;
1391      case(ReducedContactElementsOne):
1392        out=mesh->ContactElements->Tag[sampleNo];
1393        break;
1394      case(DegreesOfFreedom):
1395        throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1396        break;
1397      case(ReducedDegreesOfFreedom):
1398        throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1399        break;
1400      default:
1401        stringstream temp;
1402        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1403        throw FinleyAdapterException(temp.str());
1404        break;
1405      }
1406      return out;
1407    }
1408    
1409    
1410    void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1411    {
1412      Finley_Mesh* mesh=m_finleyMesh.get();
1413      escriptDataC tmp=mask.getDataC();
1414      switch(functionSpaceType) {
1415           case(Nodes):
1416              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1417              break;
1418           case(ReducedNodes):
1419              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1420              break;
1421           case(DegreesOfFreedom):
1422              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1423              break;
1424           case(ReducedDegreesOfFreedom):
1425              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1426              break;
1427           case(Elements):
1428              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1429              break;
1430           case(ReducedElements):
1431              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1432              break;
1433           case(FaceElements):
1434              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1435              break;
1436           case(ReducedFaceElements):
1437              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1438              break;
1439           case(Points):
1440              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1441              break;
1442           case(ContactElementsZero):
1443              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1444              break;
1445           case(ReducedContactElementsZero):
1446              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1447              break;
1448           case(ContactElementsOne):
1449              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1450              break;
1451           case(ReducedContactElementsOne):
1452              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1453              break;
1454           default:
1455              stringstream temp;
1456              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1457              throw FinleyAdapterException(temp.str());
1458      }
1459      checkFinleyError();
1460      return;
1461    }
1462    
1463    void MeshAdapter::setTagMap(const std::string& name,  int tag)
1464    {
1465      Finley_Mesh* mesh=m_finleyMesh.get();
1466      Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1467      checkPasoError();
1468      // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1469  }  }
1470    
1471  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTag(const std::string& name) const
1472  {  {
1473    int* referenceNoList;    Finley_Mesh* mesh=m_finleyMesh.get();
1474    int numReferenceNo;    int tag=0;
1475    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    tag=Finley_Mesh_getTag(mesh, name.c_str());
1476    return referenceNoList[sampleNo];    checkPasoError();
1477      // throwStandardException("MeshAdapter::getTag is not implemented.");
1478      return tag;
1479    }
1480    
1481    bool MeshAdapter::isValidTagName(const std::string& name) const
1482    {
1483      Finley_Mesh* mesh=m_finleyMesh.get();
1484      return Finley_Mesh_isValidTagName(mesh,name.c_str());
1485    }
1486    
1487    std::string MeshAdapter::showTagNames() const
1488    {
1489      stringstream temp;
1490      Finley_Mesh* mesh=m_finleyMesh.get();
1491      Finley_TagMap* tag_map=mesh->TagMap;
1492      while (tag_map) {
1493         temp << tag_map->name;
1494         tag_map=tag_map->next;
1495         if (tag_map) temp << ", ";
1496      }
1497      return string(temp.str());
1498  }  }
1499    
1500  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.149  
changed lines
  Added in v.1204

  ViewVC Help
Powered by ViewVC 1.1.26