/[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 1059 by gross, Fri Mar 23 11:01:53 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 45  const int MeshAdapter::DegreesOfFreedom= Line 33  const int MeshAdapter::DegreesOfFreedom=
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::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=FINLEY_ELEMENTS;
36    const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
37  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
38    const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
39  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=FINLEY_POINTS;
40  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
41    const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
42  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
43    const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
44    
45  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
46  {  {
# Line 85  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 77  Finley_Mesh* MeshAdapter::getFinley_Mesh
77    
78  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
79  {  {
80    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
81    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
82    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
83    checkFinleyError();    checkFinleyError();
84      TMPMEMFREE(fName);
85  }  }
86    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
87  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
88  {  {
89    return "FinleyMesh";    return "FinleyMesh";
# Line 132  void MeshAdapter::setFunctionSpaceTypeNa Line 118  void MeshAdapter::setFunctionSpaceTypeNa
118    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
119      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
120    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
121        (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
122      m_functionSpaceTypeNames.insert
123      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
124    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
125        (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
126      m_functionSpaceTypeNames.insert
127      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
128    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
129      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
130    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
131        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
132      m_functionSpaceTypeNames.insert
133      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
134      m_functionSpaceTypeNames.insert
135        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
136  }  }
137    
138  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
# Line 150  int MeshAdapter::getFunctionCode() const Line 144  int MeshAdapter::getFunctionCode() const
144  {  {
145    return Elements;    return Elements;
146  }  }
147    int MeshAdapter::getReducedFunctionCode() const
148    {
149      return ReducedElements;
150    }
151    
152  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
153  {  {
154    return FaceElements;    return FaceElements;
155  }  }
156    int MeshAdapter::getReducedFunctionOnBoundaryCode() const
157    {
158      return ReducedFaceElements;
159    }
160    
161  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
162  {  {
163    return ContactElementsZero;    return ContactElementsZero;
164  }  }
165    int MeshAdapter::getReducedFunctionOnContactZeroCode() const
166    {
167      return ReducedContactElementsZero;
168    }
169    
170  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
171  {  {
172    return ContactElementsOne;    return ContactElementsOne;
173  }  }
174    int MeshAdapter::getReducedFunctionOnContactOneCode() const
175    {
176      return ReducedContactElementsOne;
177    }
178    
179  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
180  {  {
# Line 182  int MeshAdapter::getDiracDeltaFunctionCo Line 192  int MeshAdapter::getDiracDeltaFunctionCo
192  }  }
193    
194  //  //
 // 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;  
 }  
   
 //  
195  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
196  //  //
197  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 384  pair<int,int> MeshAdapter::getDataShape( Line 248  pair<int,int> MeshAdapter::getDataShape(
248        case(DegreesOfFreedom):        case(DegreesOfFreedom):
249             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
250               numDataPointsPerSample=1;               numDataPointsPerSample=1;
251    #ifndef PASO_MPI
252               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
253    #else
254                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
255    #endif
256             }             }
257             break;             break;
258        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
259             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
260               numDataPointsPerSample=1;               numDataPointsPerSample=1;
261    #ifndef PASO_MPI
262               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
263    #else
264                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
265    #endif
266             }             }
267             break;             break;
268        default:        default:
269          stringstream temp;          stringstream temp;
270          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
271          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
272          break;          break;
273        }        }
# Line 407  pair<int,int> MeshAdapter::getDataShape( Line 278  pair<int,int> MeshAdapter::getDataShape(
278  // 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:
279  //  //
280  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
281                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
282                       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,
283                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
284                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
285  {  {
286       escriptDataC _rhs=rhs.getDataC();
287       escriptDataC _A  =A.getDataC();
288       escriptDataC _B=B.getDataC();
289       escriptDataC _C=C.getDataC();
290       escriptDataC _D=D.getDataC();
291       escriptDataC _X=X.getDataC();
292       escriptDataC _Y=Y.getDataC();
293       escriptDataC _d=d.getDataC();
294       escriptDataC _y=y.getDataC();
295       escriptDataC _d_contact=d_contact.getDataC();
296       escriptDataC _y_contact=y_contact.getDataC();
297    
298     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
299     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),  
300                         &(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 );
301     checkFinleyError();     checkFinleyError();
302     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
303                    mat.getFinley_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
                   &(rhs.getDataC()),  
                                   &(d.getDataC()),&(y.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Mean_out);  
304     checkFinleyError();     checkFinleyError();
305     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  
306                    mat.getFinley_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
                   &(rhs.getDataC()),  
                                   &(d_contact.getDataC()),  
                   &(y_contact.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Step_out);  
307     checkFinleyError();     checkFinleyError();
308  }  }
309    
310  //  //
311  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
312  //  //
313  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  
314  {  {
315     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
316    
317     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
318     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     escriptDataC _X=X.getDataC();
319     checkFinleyError();     escriptDataC _Y=Y.getDataC();
320       escriptDataC _y=y.getDataC();
321       escriptDataC _y_contact=y_contact.getDataC();
322    
323     // 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 );
324     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     checkFinleyError();
325    
326       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
327     checkFinleyError();     checkFinleyError();
328     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
329     // 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 );
330     checkFinleyError();     checkFinleyError();
331  }  }
332    
333  //  //
334  // interpolates data between different function spaces:  // interpolates data between different function spaces:
335  //  //
336  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
337  {  {
338    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
339    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
340    if (inDomain!=*this)    if (inDomain!=*this)  
341       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
342    if (targetDomain!=*this)    if (targetDomain!=*this)
343       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
344    
345    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
346    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 489  void MeshAdapter::interpolateOnDomain(Da Line 367  void MeshAdapter::interpolateOnDomain(Da
367                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
368                 break;                 break;
369             default:             default:
370                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
371                 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();
372                   throw FinleyAdapterException(temp.str());
373                 break;                 break;
374          }          }
375          break;          break;
# Line 498  void MeshAdapter::interpolateOnDomain(Da Line 377  void MeshAdapter::interpolateOnDomain(Da
377          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
378             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
379          } else {          } else {
380             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
381          }          }
382          break;          break;
383       case(FaceElements):       case(FaceElements):
384          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
385             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
386          } else {          } else {
387             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");  
388             break;             break;
389         }         }
390       case(Points):       case(Points):
391          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
392             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
393          } else {          } else {
394             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
395             break;             break;
396          }          }
397          break;          break;
# Line 524  void MeshAdapter::interpolateOnDomain(Da Line 400  void MeshAdapter::interpolateOnDomain(Da
400          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
401             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
402          } else {          } else {
403             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.");  
404             break;             break;
405          }          }
406          break;          break;
407       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
408          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
409             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
410             case(DegreesOfFreedom):             case(DegreesOfFreedom):
411             case(Nodes):             case(Nodes):
412                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
413                break;                break;
414    #ifndef PASO_MPI
415             case(Elements):             case(Elements):
416                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
417                break;                break;
# Line 549  void MeshAdapter::interpolateOnDomain(Da Line 425  void MeshAdapter::interpolateOnDomain(Da
425             case(ContactElementsOne):             case(ContactElementsOne):
426                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
427               break;               break;
428    #else
429               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
430               case(Elements):
431               {
432                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
433                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
434                  break;
435               }
436               case(FaceElements):
437               {
438                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
439                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
440                  break;
441               }
442               case(Points):
443               {
444                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
445                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
446                  break;
447               }
448               case(ContactElementsZero):
449               case(ContactElementsOne):
450               {
451                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
452                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
453                 break;
454               }
455    #endif
456             default:             default:
457               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
458               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();
459                 throw FinleyAdapterException(temp.str());
460               break;               break;
461          }          }
462          break;          break;
# Line 574  void MeshAdapter::interpolateOnDomain(Da Line 479  void MeshAdapter::interpolateOnDomain(Da
479               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
480               break;               break;
481            case(Nodes):            case(Nodes):
482               Finley_ErrorCode=TYPE_ERROR;               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
              sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");  
483               break;               break;
484            case(DegreesOfFreedom):            case(DegreesOfFreedom):
485               Finley_ErrorCode=TYPE_ERROR;               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
              sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
486               break;               break;
487            default:            default:
488               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
489               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();
490                 throw FinleyAdapterException(temp.str());
491               break;               break;
492         }         }
493         break;         break;
494       default:       default:
495          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
496          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();
497            throw FinleyAdapterException(temp.str());
498          break;          break;
499    }    }
500    checkFinleyError();    checkFinleyError();
# Line 598  void MeshAdapter::interpolateOnDomain(Da Line 503  void MeshAdapter::interpolateOnDomain(Da
503  //  //
504  // copies the locations of sample points into x:  // copies the locations of sample points into x:
505  //  //
506  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
507  {  {
508    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
509    if (argDomain!=*this)    if (argDomain!=*this)
# Line 608  void MeshAdapter::setToX(Data& arg) cons Line 513  void MeshAdapter::setToX(Data& arg) cons
513    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
514       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
515    } else {    } else {
516       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
517       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
518       // this is then interpolated onto arg:       // this is then interpolated onto arg:
519       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
# Line 619  void MeshAdapter::setToX(Data& arg) cons Line 524  void MeshAdapter::setToX(Data& arg) cons
524  //  //
525  // 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:
526  //  //
527  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
528  {  {
529    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
530    if (normalDomain!=*this)    if (normalDomain!=*this)
# Line 627  void MeshAdapter::setToNormal(Data& norm Line 532  void MeshAdapter::setToNormal(Data& norm
532    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
533    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
534      case(Nodes):      case(Nodes):
535        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
       sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");  
536        break;        break;
537      case(Elements):      case(Elements):
538        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
       sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");  
539        break;        break;
540      case (FaceElements):      case (FaceElements):
541        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
542        break;        break;
543      case(Points):      case(Points):
544        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");  
545        break;        break;
546      case (ContactElementsOne):      case (ContactElementsOne):
547      case (ContactElementsZero):      case (ContactElementsZero):
548        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
549        break;        break;
550      case(DegreesOfFreedom):      case(DegreesOfFreedom):
551        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.");  
552        break;        break;
553      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
554        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.");  
555        break;        break;
556      default:      default:
557        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
558        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();
559          throw FinleyAdapterException(temp.str());
560        break;        break;
561    }    }
562    checkFinleyError();    checkFinleyError();
# Line 664  void MeshAdapter::setToNormal(Data& norm Line 565  void MeshAdapter::setToNormal(Data& norm
565  //  //
566  // interpolates data to other domain:  // interpolates data to other domain:
567  //  //
568  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
569  {  {
570    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
571    if (targetDomain!=*this)    if (targetDomain!=*this)
572       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
573    
574    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();  
575  }  }
576    
577  //  //
578  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
579  //  //
580  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
581  {  {
582    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
583    if (argDomain!=*this)    if (argDomain!=*this)
# Line 687  void MeshAdapter::setToIntegrals(std::ve Line 586  void MeshAdapter::setToIntegrals(std::ve
586    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
587    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
588       case(Nodes):       case(Nodes):
589          Finley_ErrorCode=TYPE_ERROR;          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
         sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");  
590          break;          break;
591       case(Elements):       case(Elements):
592          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
# Line 697  void MeshAdapter::setToIntegrals(std::ve Line 595  void MeshAdapter::setToIntegrals(std::ve
595          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
596          break;          break;
597       case(Points):       case(Points):
598          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.");  
599          break;          break;
600       case(ContactElementsZero):       case(ContactElementsZero):
601          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
# Line 707  void MeshAdapter::setToIntegrals(std::ve Line 604  void MeshAdapter::setToIntegrals(std::ve
604          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
605          break;          break;
606       case(DegreesOfFreedom):       case(DegreesOfFreedom):
607          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.");  
608          break;          break;
609       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
610          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.");  
611          break;          break;
612       default:       default:
613          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
614          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();
615            throw FinleyAdapterException(temp.str());
616          break;          break;
617    }    }
618    checkFinleyError();    checkFinleyError();
# Line 725  void MeshAdapter::setToIntegrals(std::ve Line 621  void MeshAdapter::setToIntegrals(std::ve
621  //  //
622  // calculates the gradient of arg:  // calculates the gradient of arg:
623  //  //
624  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
625  {  {
626    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
627    if (argDomain!=*this)    if (argDomain!=*this)
# Line 735  void MeshAdapter::setToGradient(Data& gr Line 631  void MeshAdapter::setToGradient(Data& gr
631       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
632    
633    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
634      escriptDataC nodeDataC;
635    #ifdef PASO_MPI
636      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
637      if( arg.getFunctionSpace().getTypeCode() != Nodes )
638      {
639        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
640        nodeDataC = nodeTemp.getDataC();
641      }
642      else
643        nodeDataC = arg.getDataC();
644    #else
645      nodeDataC = arg.getDataC();
646    #endif
647    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
648         case(Nodes):         case(Nodes):
649            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");  
650            break;            break;
651         case(Elements):         case(Elements):
652            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
653            break;            break;
654         case(FaceElements):         case(FaceElements):
655            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
656            break;            break;
657         case(Points):         case(Points):
658            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
659            break;            break;
660         case(ContactElementsZero):         case(ContactElementsZero):
661            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
662            break;            break;
663         case(ContactElementsOne):         case(ContactElementsOne):
664            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
665            break;            break;
666         case(DegreesOfFreedom):         case(DegreesOfFreedom):
667            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.");  
668            break;            break;
669         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
670            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.");  
671            break;            break;
672         default:         default:
673            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
674            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();
675              throw FinleyAdapterException(temp.str());
676            break;            break;
677    }    }
678    checkFinleyError();    checkFinleyError();
# Line 775  void MeshAdapter::setToGradient(Data& gr Line 681  void MeshAdapter::setToGradient(Data& gr
681  //  //
682  // returns the size of elements:  // returns the size of elements:
683  //  //
684  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
685  {  {
686    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
687    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
688    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
689         case(Nodes):         case(Nodes):
690            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");  
691            break;            break;
692         case(Elements):         case(Elements):
693            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
# Line 791  void MeshAdapter::setToSize(Data& size) Line 696  void MeshAdapter::setToSize(Data& size)
696            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
697            break;            break;
698         case(Points):         case(Points):
699            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
700            break;            break;
701         case(ContactElementsZero):         case(ContactElementsZero):
702         case(ContactElementsOne):         case(ContactElementsOne):
703            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
704            break;            break;
705         case(DegreesOfFreedom):         case(DegreesOfFreedom):
706            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.");  
707            break;            break;
708         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
709            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.");  
710            break;            break;
711         default:         default:
712            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
713            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();
714              throw FinleyAdapterException(temp.str());
715            break;            break;
716    }    }
717    checkFinleyError();    checkFinleyError();
718  }  }
719    
720  // sets the location of nodes:  // sets the location of nodes:
721  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
722  {  {
723    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
724      escriptDataC tmp;
725    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
726    if (newDomain!=*this)    if (newDomain!=*this)
727       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
728    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
729      Finley_Mesh_setCoordinates(mesh,&tmp);
730    checkFinleyError();    checkFinleyError();
731  }  }
732    
733  // saves a data array in openDX format:  // saves a data array in openDX format:
734  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
735  {  {
736    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
737    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
738      /* win32 refactor */
739      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
740      for(int i=0;i<num_data;i++)
741      {
742        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
743      }
744    
745      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
746      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
747      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
748    
749        boost::python::list keys=arg.keys();
750        for (int i=0;i<num_data;++i) {
751             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
752             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
753                 throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
754             data[i]=d.getDataC();
755             ptr_data[i]=&(data[i]);
756             std::string n=boost::python::extract<std::string>(keys[i]);
757             c_names[i]=&(names[i][0]);
758             if (n.length()>MAX_namelength-1) {
759                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
760                c_names[i][MAX_namelength-1]='\0';
761             } else {
762                strcpy(c_names[i],n.c_str());
763             }
764        }
765        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
766        checkFinleyError();
767        
768          /* win32 refactor */
769      TMPMEMFREE(c_names);
770      TMPMEMFREE(data);
771      TMPMEMFREE(ptr_data);
772      for(int i=0;i<num_data;i++)
773      {
774        TMPMEMFREE(names[i]);
775      }
776      TMPMEMFREE(names);
777    
778        return;
779  }  }
780    
781  // saves a data array in openVTK format:  // saves a data array in openVTK format:
782  void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
783  {  {
784    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
785    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
786  }    /* win32 refactor */
787      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
788      for(int i=0;i<num_data;i++)
789      {
790        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
791      }
792    
793      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
794      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
795      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
796    
797        boost::python::list keys=arg.keys();
798        for (int i=0;i<num_data;++i) {
799             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
800             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
801                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
802             data[i]=d.getDataC();
803             ptr_data[i]=&(data[i]);
804             std::string n=boost::python::extract<std::string>(keys[i]);
805             c_names[i]=&(names[i][0]);
806             if (n.length()>MAX_namelength-1) {
807                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
808                c_names[i][MAX_namelength-1]='\0';
809             } else {
810                strcpy(c_names[i],n.c_str());
811             }
812        }
813    #ifndef PASO_MPI    
814        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
815    #else
816        Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
817    #endif
818    
819    checkFinleyError();
820      /* win32 refactor */
821      TMPMEMFREE(c_names);
822      TMPMEMFREE(data);
823      TMPMEMFREE(ptr_data);
824      for(int i=0;i<num_data;i++)
825      {
826        TMPMEMFREE(names[i]);
827      }
828      TMPMEMFREE(names);
829    
830        return;
831    }
832                                                                                                                                                                            
833                                                                                                                                                                            
834  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
835  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
836                        const int row_blocksize,                        const int row_blocksize,
# Line 873  SystemMatrixAdapter MeshAdapter::newSyst Line 865  SystemMatrixAdapter MeshAdapter::newSyst
865      }      }
866      // generate matrix:      // generate matrix:
867            
868      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);  
869      checkFinleyError();      checkFinleyError();
870      Finley_SystemMatrixPattern_dealloc(fsystemMatrixPattern);      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
871        checkPasoError();
872        Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
873      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
874  }  }
875    
# Line 904  bool MeshAdapter::isCellOriented(int fun Line 896  bool MeshAdapter::isCellOriented(int fun
896            return true;            return true;
897            break;            break;
898         default:         default:
899            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
900            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;
901              throw FinleyAdapterException(temp.str());
902            break;            break;
903    }    }
904    checkFinleyError();    checkFinleyError();
# Line 927  bool MeshAdapter::probeInterpolationOnDo Line 920  bool MeshAdapter::probeInterpolationOnDo
920             case(ContactElementsOne):             case(ContactElementsOne):
921                 return true;                 return true;
922             default:             default:
923                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
924                 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;
925                   throw FinleyAdapterException(temp.str());
926          }          }
927          break;          break;
928       case(Elements):       case(Elements):
# Line 968  bool MeshAdapter::probeInterpolationOnDo Line 962  bool MeshAdapter::probeInterpolationOnDo
962             case(ContactElementsOne):             case(ContactElementsOne):
963                return true;                return true;
964             default:             default:
965               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
966               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;
967                 throw FinleyAdapterException(temp.str());
968          }          }
969          break;          break;
970       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
# Line 985  bool MeshAdapter::probeInterpolationOnDo Line 980  bool MeshAdapter::probeInterpolationOnDo
980            case(DegreesOfFreedom):            case(DegreesOfFreedom):
981               return false;               return false;
982            default:            default:
983               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
984               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;
985                 throw FinleyAdapterException(temp.str());
986         }         }
987         break;         break;
988       default:       default:
989          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
990          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;
991            throw FinleyAdapterException(temp.str());
992          break;          break;
993    }    }
994    checkFinleyError();    checkFinleyError();
# Line 1018  bool MeshAdapter::operator!=(const Abstr Line 1015  bool MeshAdapter::operator!=(const Abstr
1015    return !(operator==(other));    return !(operator==(other));
1016  }  }
1017    
1018  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1019  {  {
1020     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1021     checkFinleyError();     checkPasoError();
1022     return out;     return out;
1023  }  }
1024    
1025  Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
1026  {  {
1027    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
1028  }  }
1029    
1030  Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
1031  {  {
1032    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1033  }  }
1034    
1035  Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
1036  {  {
1037    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
1038  }  }
1039    
1040  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1041  {  {
1042    int* tagList;    int out=0;
1043    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1044    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1045    return tagList[sampleNo];    case(Nodes):
1046        out=mesh->Nodes->Tag[sampleNo];
1047        break;
1048      case(Elements):
1049        out=mesh->Elements->Tag[sampleNo];
1050        break;
1051      case(FaceElements):
1052        out=mesh->FaceElements->Tag[sampleNo];
1053        break;
1054      case(Points):
1055        out=mesh->Points->Tag[sampleNo];
1056        break;
1057      case(ContactElementsZero):
1058        out=mesh->ContactElements->Tag[sampleNo];
1059        break;
1060      case(ContactElementsOne):
1061        out=mesh->ContactElements->Tag[sampleNo];
1062        break;
1063      case(DegreesOfFreedom):
1064        break;
1065      case(ReducedDegreesOfFreedom):
1066        break;
1067      default:
1068        stringstream temp;
1069        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1070        throw FinleyAdapterException(temp.str());
1071        break;
1072      }
1073      return out;
1074    }
1075    int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1076    {
1077      int *out=0,i;
1078      Finley_Mesh* mesh=m_finleyMesh.get();
1079      switch (functionSpaceType) {
1080      case(Nodes):
1081        if (mesh->Nodes!=NULL) {
1082          out=mesh->Nodes->Id;
1083          break;
1084        }
1085      case(Elements):
1086        out=mesh->Elements->Id;
1087        break;
1088      case(FaceElements):
1089        out=mesh->FaceElements->Id;
1090        break;
1091      case(Points):
1092        out=mesh->Points->Id;
1093        break;
1094      case(ContactElementsZero):
1095        out=mesh->ContactElements->Id;
1096        break;
1097      case(ContactElementsOne):
1098        out=mesh->ContactElements->Id;
1099        break;
1100      case(DegreesOfFreedom):
1101        out=mesh->Nodes->degreeOfFreedomId;
1102        break;
1103      case(ReducedDegreesOfFreedom):
1104        out=mesh->Nodes->reducedDegreeOfFreedomId;
1105        break;
1106      default:
1107        stringstream temp;
1108        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1109        throw FinleyAdapterException(temp.str());
1110        break;
1111      }
1112      return out;
1113  }  }
1114    
1115  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1116  {  {
1117    int* referenceNoList;    Finley_Mesh* mesh=m_finleyMesh.get();
1118    int numReferenceNo;    escriptDataC tmp=mask.getDataC();
1119    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    switch(functionSpaceType) {
1120    return referenceNoList[sampleNo];         case(Nodes):
1121              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1122              break;
1123           case(DegreesOfFreedom):
1124              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1125              break;
1126           case(ReducedDegreesOfFreedom):
1127              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1128              break;
1129           case(Elements):
1130              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1131              break;
1132           case(FaceElements):
1133              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1134              break;
1135           case(Points):
1136              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1137              break;
1138           case(ContactElementsZero):
1139              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1140              break;
1141           case(ContactElementsOne):
1142              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1143              break;
1144           default:
1145              stringstream temp;
1146              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1147              throw FinleyAdapterException(temp.str());
1148      }
1149      checkFinleyError();
1150      return;
1151    }
1152    
1153    void MeshAdapter::setTagMap(const std::string& name,  int tag)
1154    {
1155      Finley_Mesh* mesh=m_finleyMesh.get();
1156      Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1157      checkPasoError();
1158      // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1159    }
1160    
1161    int MeshAdapter::getTag(const std::string& name) const
1162    {
1163      Finley_Mesh* mesh=m_finleyMesh.get();
1164      int tag=0;
1165      tag=Finley_Mesh_getTag(mesh, name.c_str());
1166      checkPasoError();
1167      // throwStandardException("MeshAdapter::getTag is not implemented.");
1168      return tag;
1169    }
1170    
1171    bool MeshAdapter::isValidTagName(const std::string& name) const
1172    {
1173      Finley_Mesh* mesh=m_finleyMesh.get();
1174      return Finley_Mesh_isValidTagName(mesh,name.c_str());
1175    }
1176    
1177    std::string MeshAdapter::showTagNames() const
1178    {
1179      stringstream temp;
1180      Finley_Mesh* mesh=m_finleyMesh.get();
1181      Finley_TagMap* tag_map=mesh->TagMap;
1182      while (tag_map) {
1183         temp << tag_map->name;
1184         tag_map=tag_map->next;
1185         if (tag_map) temp << ", ";
1186      }
1187      return temp.str();
1188  }  }
1189    
1190  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26