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

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

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

trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 751 by bcumming, Mon Jun 26 01:46:34 2006 UTC
# Line 12  Line 12 
12   *                                                                            *   *                                                                            *
13   ******************************************************************************   ******************************************************************************
14  */  */
15  extern "C" {  
16  #include "finley/finleyC/Finley.h"  #include "MeshAdapter.h"
17  #include "finley/finleyC/Assemble.h"  
18  #include "finley/finleyC/Mesh.h"  #include "escript/Data.h"
19  #include "finley/finleyC/Finley.h"  #include "escript/DataFactory.h"
 #include "finley/finleyC/System.h"  
 }  
 #include "finley/CPPAdapter/SystemMatrixAdapter.h"  
 #include "finley/CPPAdapter/MeshAdapter.h"  
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/FinleyAdapterException.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/Data.h"  
 #include "escript/Data/DataArrayView.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/DataFactory.h"  
 #include <iostream>  
 #include <vector>  
 #include <sstream>  
20    
21  using namespace std;  using namespace std;
22  using namespace escript;  using namespace escript;
23    
24  namespace finley {  namespace finley {
25    
26  //  //
27  // define the statics  // define the static constants
28  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
29  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
30  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
# Line 56  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 43  MeshAdapter::MeshAdapter(Finley_Mesh* fi
43    // for us.    // for us.
44    m_finleyMesh.reset(finleyMesh,null_deleter());    m_finleyMesh.reset(finleyMesh,null_deleter());
45  }  }
46    
47  //  //
48  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
49  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
# Line 88  void MeshAdapter::write(const std::strin Line 76  void MeshAdapter::write(const std::strin
76    checkFinleyError();    checkFinleyError();
77  }  }
78    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
79  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
80  {  {
81    return "FinleyMesh";    return "FinleyMesh";
# Line 142  int MeshAdapter::getContinuousFunctionCo Line 123  int MeshAdapter::getContinuousFunctionCo
123  {  {
124    return Nodes;    return Nodes;
125  }  }
126    
127  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
128  {  {
129    return Elements;    return Elements;
130  }  }
131    
132  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
133  {  {
134    return FaceElements;    return FaceElements;
135  }  }
136    
137  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
138  {  {
139    return ContactElementsZero;    return ContactElementsZero;
140  }  }
141    
142  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
143  {  {
144    return ContactElementsOne;    return ContactElementsOne;
# Line 163  int MeshAdapter::getSolutionCode() const Line 148  int MeshAdapter::getSolutionCode() const
148  {  {
149    return DegreesOfFreedom;    return DegreesOfFreedom;
150  }  }
151    
152  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
153  {  {
154    return ReducedDegreesOfFreedom;    return ReducedDegreesOfFreedom;
155  }  }
156    
157  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionCode() const
158  {  {
159    return Points;    return Points;
160  }  }
 //  
 // 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;  
 }  
161    
162  //  //
163  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
# Line 325  int MeshAdapter::getDim() const Line 168  int MeshAdapter::getDim() const
168    checkFinleyError();    checkFinleyError();
169    return numDim;    return numDim;
170  }  }
171    
172  //  //
173  // return the number of data points per sample and the number of samples  // return the number of data points per sample and the number of samples
174  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
# Line 372  pair<int,int> MeshAdapter::getDataShape( Line 216  pair<int,int> MeshAdapter::getDataShape(
216        case(DegreesOfFreedom):        case(DegreesOfFreedom):
217             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
218               numDataPointsPerSample=1;               numDataPointsPerSample=1;
219    #ifndef PASO_MPI
220               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
221    #else
222                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
223    #endif
224             }             }
225             break;             break;
226        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
227             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
228               numDataPointsPerSample=1;               numDataPointsPerSample=1;
229    #ifndef PASO_MPI
230               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
231    #else
232                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
233    #endif
234             }             }
235             break;             break;
236        default:        default:
237          stringstream temp;          stringstream temp;
238          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
239          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
240          break;          break;
241        }        }
242        return pair<int,int>(numDataPointsPerSample,numSamples);        return pair<int,int>(numDataPointsPerSample,numSamples);
243  }  }
244    
245  //  //
246  // 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:
247  //  //
248  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
249                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
250                       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,
251                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
252                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
253  {  {
254      escriptDataC _rhs=rhs.getDataC();
255      escriptDataC _A  =A.getDataC();
256      escriptDataC _B=B.getDataC();
257      escriptDataC _C=C.getDataC();
258      escriptDataC _D=D.getDataC();
259      escriptDataC _X=X.getDataC();
260      escriptDataC _Y=Y.getDataC();
261      escriptDataC _d=d.getDataC();
262      escriptDataC _y=y.getDataC();
263      escriptDataC _d_contact=d_contact.getDataC();
264      escriptDataC _y_contact=y_contact.getDataC();
265    
266     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
267     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
268                         &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));  
269     checkFinleyError();     checkFinleyError();
270     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
271                    mat.getFinley_SystemMatrix(),     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, &_d, &_y, Finley_Assemble_handelShapeMissMatch_Mean_out);
                   &(rhs.getDataC()),  
                                   &(d.getDataC()),&(y.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Mean_out);  
272     checkFinleyError();     checkFinleyError();
273     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  
274                    mat.getFinley_SystemMatrix(),     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , &_d_contact, &_y_contact ,             Finley_Assemble_handelShapeMissMatch_Step_out);
                   &(rhs.getDataC()),  
                                   &(d_contact.getDataC()),  
                   &(y_contact.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Step_out);  
275     checkFinleyError();     checkFinleyError();
276  }  }
277    
278  //  //
279  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
280  //  //
281  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  
282  {  {
283     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
284    
# Line 437  void MeshAdapter::addPDEToRHS( Data& rhs Line 294  void MeshAdapter::addPDEToRHS( Data& rhs
294     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
295     checkFinleyError();     checkFinleyError();
296  }  }
297    
298  //  //
299  // interpolates data between different function spaces:  // interpolates data between different function spaces:
300  //  //
301  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
302  {  {
303    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
304    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
305    if (inDomain!=*this)    if (inDomain!=*this)  
306       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
307    if (targetDomain!=*this)    if (targetDomain!=*this)
308       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
309    
310    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
311    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 474  void MeshAdapter::interpolateOnDomain(Da Line 332  void MeshAdapter::interpolateOnDomain(Da
332                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
333                 break;                 break;
334             default:             default:
335                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
336                 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();
337                   throw FinleyAdapterException(temp.str());
338                 break;                 break;
339          }          }
340          break;          break;
# Line 483  void MeshAdapter::interpolateOnDomain(Da Line 342  void MeshAdapter::interpolateOnDomain(Da
342          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
343             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
344          } else {          } else {
345             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
346          }          }
347          break;          break;
348       case(FaceElements):       case(FaceElements):
349          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
350             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
351          } else {          } else {
352             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.");  
353             break;             break;
354         }         }
355       case(Points):       case(Points):
356          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
357             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
358          } else {          } else {
359             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
360             break;             break;
361          }          }
362          break;          break;
# Line 509  void MeshAdapter::interpolateOnDomain(Da Line 365  void MeshAdapter::interpolateOnDomain(Da
365          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
366             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
367          } else {          } else {
368             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.");  
369             break;             break;
370          }          }
371          break;          break;
372       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
373          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
374             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
375             case(DegreesOfFreedom):             case(DegreesOfFreedom):
376             case(Nodes):             case(Nodes):
377                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
378                break;                break;
379    #ifndef PASO_MPI
380             case(Elements):             case(Elements):
381                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
382                break;                break;
# Line 534  void MeshAdapter::interpolateOnDomain(Da Line 390  void MeshAdapter::interpolateOnDomain(Da
390             case(ContactElementsOne):             case(ContactElementsOne):
391                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
392               break;               break;
393    #else
394               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
395               case(Elements):
396               {
397                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
398                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
399                  break;
400               }
401               case(FaceElements):
402               {
403                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
404                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
405                  break;
406               }
407               case(Points):
408               {
409                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
410                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
411                  break;
412               }
413               case(ContactElementsZero):
414               case(ContactElementsOne):
415               {
416                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
417                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
418                 break;
419               }
420    #endif
421             default:             default:
422               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
423               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();
424                 throw FinleyAdapterException(temp.str());
425               break;               break;
426          }          }
427          break;          break;
# Line 559  void MeshAdapter::interpolateOnDomain(Da Line 444  void MeshAdapter::interpolateOnDomain(Da
444               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
445               break;               break;
446            case(Nodes):            case(Nodes):
447               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.");  
448               break;               break;
449            case(DegreesOfFreedom):            case(DegreesOfFreedom):
450               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");  
451               break;               break;
452            default:            default:
453               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
454               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();
455                 throw FinleyAdapterException(temp.str());
456               break;               break;
457         }         }
458         break;         break;
459       default:       default:
460          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
461          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();
462            throw FinleyAdapterException(temp.str());
463          break;          break;
464    }    }
465    checkFinleyError();    checkFinleyError();
# Line 583  void MeshAdapter::interpolateOnDomain(Da Line 468  void MeshAdapter::interpolateOnDomain(Da
468  //  //
469  // copies the locations of sample points into x:  // copies the locations of sample points into x:
470  //  //
471  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
472  {  {
473    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
474    if (argDomain!=*this)    if (argDomain!=*this)
# Line 593  void MeshAdapter::setToX(Data& arg) cons Line 478  void MeshAdapter::setToX(Data& arg) cons
478    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
479       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
480    } else {    } else {
481       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
482       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
483       // this is then interpolated onto arg:       // this is then interpolated onto arg:
484       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
485    }    }
486    checkFinleyError();    checkFinleyError();
487  }  }
488    
489  //  //
490  // 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:
491  //  //
492  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
493  {  {
494    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
495    if (normalDomain!=*this)    if (normalDomain!=*this)
# Line 611  void MeshAdapter::setToNormal(Data& norm Line 497  void MeshAdapter::setToNormal(Data& norm
497    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
498    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
499      case(Nodes):      case(Nodes):
500        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");  
501        break;        break;
502      case(Elements):      case(Elements):
503        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");  
504        break;        break;
505      case (FaceElements):      case (FaceElements):
506        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
507        break;        break;
508      case(Points):      case(Points):
509        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");  
510        break;        break;
511      case (ContactElementsOne):      case (ContactElementsOne):
512      case (ContactElementsZero):      case (ContactElementsZero):
513        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
514        break;        break;
515      case(DegreesOfFreedom):      case(DegreesOfFreedom):
516        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.");  
517        break;        break;
518      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
519        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.");  
520        break;        break;
521      default:      default:
522        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
523        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();
524          throw FinleyAdapterException(temp.str());
525        break;        break;
526    }    }
527    checkFinleyError();    checkFinleyError();
528  }  }
529    
530  //  //
531  // interpolates data to other domain:  // interpolates data to other domain:
532  //  //
533  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
534  {  {
535    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
536    if (targetDomain!=*this)    if (targetDomain!=*this)
537       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
538    
539    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();  
540  }  }
541    
542  //  //
543  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
544  //  //
545  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
546  {  {
547    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
548    if (argDomain!=*this)    if (argDomain!=*this)
# Line 669  void MeshAdapter::setToIntegrals(std::ve Line 551  void MeshAdapter::setToIntegrals(std::ve
551    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
552    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
553       case(Nodes):       case(Nodes):
554          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.");  
555          break;          break;
556       case(Elements):       case(Elements):
557          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
# Line 679  void MeshAdapter::setToIntegrals(std::ve Line 560  void MeshAdapter::setToIntegrals(std::ve
560          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
561          break;          break;
562       case(Points):       case(Points):
563          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.");  
564          break;          break;
565       case(ContactElementsZero):       case(ContactElementsZero):
566          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
# Line 689  void MeshAdapter::setToIntegrals(std::ve Line 569  void MeshAdapter::setToIntegrals(std::ve
569          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
570          break;          break;
571       case(DegreesOfFreedom):       case(DegreesOfFreedom):
572          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.");  
573          break;          break;
574       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
575          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.");  
576          break;          break;
577       default:       default:
578          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
579          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();
580            throw FinleyAdapterException(temp.str());
581          break;          break;
582    }    }
583    checkFinleyError();    checkFinleyError();
584  }  }
585    
586  //  //
587  // calculates the gradient of arg:  // calculates the gradient of arg:
588  //  //
589  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
590  {  {
591    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
592    if (argDomain!=*this)    if (argDomain!=*this)
# Line 716  void MeshAdapter::setToGradient(Data& gr Line 596  void MeshAdapter::setToGradient(Data& gr
596       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
597    
598    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
599      escriptDataC nodeDataC;
600    #ifdef PASO_MPI
601      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
602      if( arg.getFunctionSpace().getTypeCode() != Nodes )
603      {
604        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
605        nodeDataC = nodeTemp.getDataC();
606      }
607      else
608        nodeDataC = arg.getDataC();
609    #else
610      nodeDataC = arg.getDataC();
611    #endif
612    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
613         case(Nodes):         case(Nodes):
614            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");  
615            break;            break;
616         case(Elements):         case(Elements):
617            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
618            break;            break;
619         case(FaceElements):         case(FaceElements):
620            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
621            break;            break;
622         case(Points):         case(Points):
623            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
624            break;            break;
625         case(ContactElementsZero):         case(ContactElementsZero):
626            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
627            break;            break;
628         case(ContactElementsOne):         case(ContactElementsOne):
629            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
630            break;            break;
631         case(DegreesOfFreedom):         case(DegreesOfFreedom):
632            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.");  
633            break;            break;
634         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
635            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.");  
636            break;            break;
637         default:         default:
638            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
639            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();
640              throw FinleyAdapterException(temp.str());
641            break;            break;
642    }    }
643    checkFinleyError();    checkFinleyError();
644  }  }
645    
646  //  //
647  // returns the size of elements:  // returns the size of elements:
648  //  //
649  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
650  {  {
651    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
652    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
653    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
654         case(Nodes):         case(Nodes):
655            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");  
656            break;            break;
657         case(Elements):         case(Elements):
658            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
# Line 771  void MeshAdapter::setToSize(Data& size) Line 661  void MeshAdapter::setToSize(Data& size)
661            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
662            break;            break;
663         case(Points):         case(Points):
664            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
665            break;            break;
666         case(ContactElementsZero):         case(ContactElementsZero):
667         case(ContactElementsOne):         case(ContactElementsOne):
668            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
669            break;            break;
670         case(DegreesOfFreedom):         case(DegreesOfFreedom):
671            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.");  
672            break;            break;
673         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
674            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.");  
675            break;            break;
676         default:         default:
677            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
678            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();
679              throw FinleyAdapterException(temp.str());
680            break;            break;
681    }    }
682    checkFinleyError();    checkFinleyError();
683  }  }
684    
685  // sets the location of nodes:  // sets the location of nodes:
686  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
687  {  {
688    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
689      escriptDataC tmp;
690    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
691    if (newDomain!=*this)    if (newDomain!=*this)
692       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
693    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
694      Finley_Mesh_setCoordinates(mesh,&tmp);
695    checkFinleyError();    checkFinleyError();
696  }  }
697    
698  // saves a data array in openDX format:  // saves a data array in openDX format:
699  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
700  {  {
701    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
702    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
703        char names[num_data][MAX_namelength];
704        char* c_names[num_data];
705        escriptDataC data[num_data];
706        escriptDataC* ptr_data[num_data];
707    
708        boost::python::list keys=arg.keys();
709        for (int i=0;i<num_data;++i) {
710             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
711             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
712                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
713             data[i]=d.getDataC();
714             ptr_data[i]=&(data[i]);
715             std::string n=boost::python::extract<std::string>(keys[i]);
716             c_names[i]=&(names[i][0]);
717             if (n.length()>MAX_namelength-1) {
718                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
719                c_names[i][MAX_namelength-1]='\0';
720             } else {
721                strcpy(c_names[i],n.c_str());
722             }
723        }
724        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
725        checkFinleyError();
726        return;
727  }  }
728    
729  // saves a data array in openVTK format:  // saves a data array in openVTK format:
730  void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
731  {  {
732    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
733    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
734        char names[num_data][MAX_namelength];
735        char* c_names[num_data];
736        escriptDataC data[num_data];
737        escriptDataC* ptr_data[num_data];
738    
739        boost::python::list keys=arg.keys();
740        for (int i=0;i<num_data;++i) {
741             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
742             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
743                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
744             data[i]=d.getDataC();
745             ptr_data[i]=&(data[i]);
746             std::string n=boost::python::extract<std::string>(keys[i]);
747             c_names[i]=&(names[i][0]);
748             if (n.length()>MAX_namelength-1) {
749                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
750                c_names[i][MAX_namelength-1]='\0';
751             } else {
752                strcpy(c_names[i],n.c_str());
753             }
754        }
755        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
756        checkFinleyError();
757        return;
758  }  }
759                                                                                                                                                                            
760                                                                                                                                                                            
761  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
762  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
763                        const int row_blocksize,                        const int row_blocksize,
# Line 849  SystemMatrixAdapter MeshAdapter::newSyst Line 792  SystemMatrixAdapter MeshAdapter::newSyst
792      }      }
793      // generate matrix:      // generate matrix:
794            
795      Finley_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
796      checkFinleyError();      checkFinleyError();
797      Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);      Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
798      checkFinleyError();      checkPasoError();
799      Finley_SystemMatrixPattern_dealloc(fsystemMatrixPattern);      Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
800      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
801  }  }
802    
803  //  //
804  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
805  // TODO:  // TODO:
806  //  //
807    
808  //  //
809  // returns true if data at the atom_type is considered as being cell centered:  // returns true if data at the atom_type is considered as being cell centered:
810  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
# Line 878  bool MeshAdapter::isCellOriented(int fun Line 823  bool MeshAdapter::isCellOriented(int fun
823            return true;            return true;
824            break;            break;
825         default:         default:
826            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
827            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;
828              throw FinleyAdapterException(temp.str());
829            break;            break;
830    }    }
831    checkFinleyError();    checkFinleyError();
832    return false;    return false;
833  }  }
834    
835  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
836  {  {
837    switch(functionSpaceType_source) {    switch(functionSpaceType_source) {
# Line 900  bool MeshAdapter::probeInterpolationOnDo Line 847  bool MeshAdapter::probeInterpolationOnDo
847             case(ContactElementsOne):             case(ContactElementsOne):
848                 return true;                 return true;
849             default:             default:
850                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
851                 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;
852                   throw FinleyAdapterException(temp.str());
853          }          }
854          break;          break;
855       case(Elements):       case(Elements):
# Line 941  bool MeshAdapter::probeInterpolationOnDo Line 889  bool MeshAdapter::probeInterpolationOnDo
889             case(ContactElementsOne):             case(ContactElementsOne):
890                return true;                return true;
891             default:             default:
892               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
893               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;
894                 throw FinleyAdapterException(temp.str());
895          }          }
896          break;          break;
897       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
# Line 958  bool MeshAdapter::probeInterpolationOnDo Line 907  bool MeshAdapter::probeInterpolationOnDo
907            case(DegreesOfFreedom):            case(DegreesOfFreedom):
908               return false;               return false;
909            default:            default:
910               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
911               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;
912                 throw FinleyAdapterException(temp.str());
913         }         }
914         break;         break;
915       default:       default:
916          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
917          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;
918            throw FinleyAdapterException(temp.str());
919          break;          break;
920    }    }
921    checkFinleyError();    checkFinleyError();
922    return false;    return false;
923  }  }
924    
925  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
926  {  {
927      return false;      return false;
# Line 990  bool MeshAdapter::operator!=(const Abstr Line 942  bool MeshAdapter::operator!=(const Abstr
942    return !(operator==(other));    return !(operator==(other));
943  }  }
944    
945  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
946  {  {
947     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
948     checkFinleyError();     checkPasoError();
949     return out;     return out;
950  }  }
951  Data MeshAdapter::getX() const  
952    escript::Data MeshAdapter::getX() const
953  {  {
954    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
955  }  }
956  Data MeshAdapter::getNormal() const  
957    escript::Data MeshAdapter::getNormal() const
958  {  {
959    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
960  }  }
961  Data MeshAdapter::getSize() const  
962    escript::Data MeshAdapter::getSize() const
963  {  {
964    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
965  }  }
966    
967  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
968  {  {
969    int* tagList;    int out=0;
970    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
971    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
972    return tagList[sampleNo];    case(Nodes):
973        out=mesh->Nodes->Tag[sampleNo];
974        break;
975      case(Elements):
976        out=mesh->Elements->Tag[sampleNo];
977        break;
978      case(FaceElements):
979        out=mesh->FaceElements->Tag[sampleNo];
980        break;
981      case(Points):
982        out=mesh->Points->Tag[sampleNo];
983        break;
984      case(ContactElementsZero):
985        out=mesh->ContactElements->Tag[sampleNo];
986        break;
987      case(ContactElementsOne):
988        out=mesh->ContactElements->Tag[sampleNo];
989        break;
990      case(DegreesOfFreedom):
991        break;
992      case(ReducedDegreesOfFreedom):
993        break;
994      default:
995        stringstream temp;
996        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
997        throw FinleyAdapterException(temp.str());
998        break;
999      }
1000      return out;
1001  }  }
   
1002  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1003  {  {
1004    int* referenceNoList;    int out=0,i;
1005    int numReferenceNo;    Finley_Mesh* mesh=m_finleyMesh.get();
1006    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    switch (functionSpaceType) {
1007    return referenceNoList[sampleNo];    case(Nodes):
1008        if (mesh->Nodes!=NULL) {
1009          out=mesh->Nodes->Id[sampleNo];
1010          break;
1011        }
1012      case(Elements):
1013        out=mesh->Elements->Id[sampleNo];
1014        break;
1015      case(FaceElements):
1016        out=mesh->FaceElements->Id[sampleNo];
1017        break;
1018      case(Points):
1019        out=mesh->Points->Id[sampleNo];
1020        break;
1021      case(ContactElementsZero):
1022        out=mesh->ContactElements->Id[sampleNo];
1023        break;
1024      case(ContactElementsOne):
1025        out=mesh->ContactElements->Id[sampleNo];
1026        break;
1027      case(DegreesOfFreedom):
1028        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1029           if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
1030              out=mesh->Nodes->Id[i];
1031              break;
1032           }
1033        }
1034        break;
1035      case(ReducedDegreesOfFreedom):
1036        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1037           if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
1038              out=mesh->Nodes->Id[i];
1039              break;
1040           }
1041        }
1042        break;
1043      default:
1044        stringstream temp;
1045        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1046        throw FinleyAdapterException(temp.str());
1047        break;
1048      }
1049      return out;
1050  }  }
1051    
1052  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.148  
changed lines
  Added in v.751

  ViewVC Help
Powered by ViewVC 1.1.26