/[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 147 by jgs, Fri Aug 12 01:45:47 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 757 by woo409, Mon Jun 26 13:12:56 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 82  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 70  Finley_Mesh* MeshAdapter::getFinley_Mesh
70    
71  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
72  {  {
73    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
74    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
75    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
76    checkFinleyError();    checkFinleyError();
77      TMPMEMFREE(fName);
78  }  }
79    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
80  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
81  {  {
82    return "FinleyMesh";    return "FinleyMesh";
# Line 142  int MeshAdapter::getContinuousFunctionCo Line 124  int MeshAdapter::getContinuousFunctionCo
124  {  {
125    return Nodes;    return Nodes;
126  }  }
127    
128  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
129  {  {
130    return Elements;    return Elements;
131  }  }
132    
133  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
134  {  {
135    return FaceElements;    return FaceElements;
136  }  }
137    
138  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
139  {  {
140    return ContactElementsZero;    return ContactElementsZero;
141  }  }
142    
143  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
144  {  {
145    return ContactElementsOne;    return ContactElementsOne;
# Line 163  int MeshAdapter::getSolutionCode() const Line 149  int MeshAdapter::getSolutionCode() const
149  {  {
150    return DegreesOfFreedom;    return DegreesOfFreedom;
151  }  }
152    
153  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
154  {  {
155    return ReducedDegreesOfFreedom;    return ReducedDegreesOfFreedom;
156  }  }
157    
158  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionCode() const
159  {  {
160    return Points;    return Points;
161  }  }
 //  
 // 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;  
 }  
162    
163  //  //
164  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
# Line 325  int MeshAdapter::getDim() const Line 169  int MeshAdapter::getDim() const
169    checkFinleyError();    checkFinleyError();
170    return numDim;    return numDim;
171  }  }
172    
173  //  //
174  // 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
175  // 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 217  pair<int,int> MeshAdapter::getDataShape(
217        case(DegreesOfFreedom):        case(DegreesOfFreedom):
218             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
219               numDataPointsPerSample=1;               numDataPointsPerSample=1;
220    #ifndef PASO_MPI
221               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
222    #else
223                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
224    #endif
225             }             }
226             break;             break;
227        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
228             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
229               numDataPointsPerSample=1;               numDataPointsPerSample=1;
230    #ifndef PASO_MPI
231               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
232    #else
233                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
234    #endif
235             }             }
236             break;             break;
237        default:        default:
238          stringstream temp;          stringstream temp;
239          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
240          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
241          break;          break;
242        }        }
243        return pair<int,int>(numDataPointsPerSample,numSamples);        return pair<int,int>(numDataPointsPerSample,numSamples);
244  }  }
245    
246  //  //
247  // 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:
248  //  //
249  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
250                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
251                       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,
252                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
253                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
254  {  {
255      escriptDataC _rhs=rhs.getDataC();
256      escriptDataC _A  =A.getDataC();
257      escriptDataC _B=B.getDataC();
258      escriptDataC _C=C.getDataC();
259      escriptDataC _D=D.getDataC();
260      escriptDataC _X=X.getDataC();
261      escriptDataC _Y=Y.getDataC();
262      escriptDataC _d=d.getDataC();
263      escriptDataC _y=y.getDataC();
264      escriptDataC _d_contact=d_contact.getDataC();
265      escriptDataC _y_contact=y_contact.getDataC();
266    
267     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
268     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 );
269                         &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));  
270     checkFinleyError();     checkFinleyError();
271     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
272                    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);  
273     checkFinleyError();     checkFinleyError();
274     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  
275                    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);  
276     checkFinleyError();     checkFinleyError();
277  }  }
278    
279  //  //
280  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
281  //  //
282  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  
283  {  {
284     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
285    
# Line 437  void MeshAdapter::addPDEToRHS( Data& rhs Line 295  void MeshAdapter::addPDEToRHS( Data& rhs
295     // 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);
296     checkFinleyError();     checkFinleyError();
297  }  }
298    
299  //  //
300  // interpolates data between different function spaces:  // interpolates data between different function spaces:
301  //  //
302  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
303  {  {
304    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
305    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
306    if (inDomain!=*this)    if (inDomain!=*this)  
307       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
308    if (targetDomain!=*this)    if (targetDomain!=*this)
309       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
310    
311    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
312    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 474  void MeshAdapter::interpolateOnDomain(Da Line 333  void MeshAdapter::interpolateOnDomain(Da
333                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
334                 break;                 break;
335             default:             default:
336                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
337                 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();
338                   throw FinleyAdapterException(temp.str());
339                 break;                 break;
340          }          }
341          break;          break;
# Line 483  void MeshAdapter::interpolateOnDomain(Da Line 343  void MeshAdapter::interpolateOnDomain(Da
343          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
344             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
345          } else {          } else {
346             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
347          }          }
348          break;          break;
349       case(FaceElements):       case(FaceElements):
350          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
351             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
352          } else {          } else {
353             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.");  
354             break;             break;
355         }         }
356       case(Points):       case(Points):
357          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
358             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
359          } else {          } else {
360             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
361             break;             break;
362          }          }
363          break;          break;
# Line 509  void MeshAdapter::interpolateOnDomain(Da Line 366  void MeshAdapter::interpolateOnDomain(Da
366          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
367             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
368          } else {          } else {
369             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.");  
370             break;             break;
371          }          }
372          break;          break;
373       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
374          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
375             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
376             case(DegreesOfFreedom):             case(DegreesOfFreedom):
377             case(Nodes):             case(Nodes):
378                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
379                break;                break;
380    #ifndef PASO_MPI
381             case(Elements):             case(Elements):
382                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
383                break;                break;
# Line 534  void MeshAdapter::interpolateOnDomain(Da Line 391  void MeshAdapter::interpolateOnDomain(Da
391             case(ContactElementsOne):             case(ContactElementsOne):
392                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
393               break;               break;
394    #else
395               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
396               case(Elements):
397               {
398                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
399                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
400                  break;
401               }
402               case(FaceElements):
403               {
404                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
405                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
406                  break;
407               }
408               case(Points):
409               {
410                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
411                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
412                  break;
413               }
414               case(ContactElementsZero):
415               case(ContactElementsOne):
416               {
417                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
418                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
419                 break;
420               }
421    #endif
422             default:             default:
423               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
424               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();
425                 throw FinleyAdapterException(temp.str());
426               break;               break;
427          }          }
428          break;          break;
# Line 559  void MeshAdapter::interpolateOnDomain(Da Line 445  void MeshAdapter::interpolateOnDomain(Da
445               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
446               break;               break;
447            case(Nodes):            case(Nodes):
448               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.");  
449               break;               break;
450            case(DegreesOfFreedom):            case(DegreesOfFreedom):
451               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");  
452               break;               break;
453            default:            default:
454               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
455               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();
456                 throw FinleyAdapterException(temp.str());
457               break;               break;
458         }         }
459         break;         break;
460       default:       default:
461          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
462          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();
463            throw FinleyAdapterException(temp.str());
464          break;          break;
465    }    }
466    checkFinleyError();    checkFinleyError();
# Line 583  void MeshAdapter::interpolateOnDomain(Da Line 469  void MeshAdapter::interpolateOnDomain(Da
469  //  //
470  // copies the locations of sample points into x:  // copies the locations of sample points into x:
471  //  //
472  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
473  {  {
474    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
475    if (argDomain!=*this)    if (argDomain!=*this)
# Line 593  void MeshAdapter::setToX(Data& arg) cons Line 479  void MeshAdapter::setToX(Data& arg) cons
479    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
480       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
481    } else {    } else {
482       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
483       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
484       // this is then interpolated onto arg:       // this is then interpolated onto arg:
485       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
486    }    }
487    checkFinleyError();    checkFinleyError();
488  }  }
489    
490  //  //
491  // 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:
492  //  //
493  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
494  {  {
495    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
496    if (normalDomain!=*this)    if (normalDomain!=*this)
# Line 611  void MeshAdapter::setToNormal(Data& norm Line 498  void MeshAdapter::setToNormal(Data& norm
498    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
499    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
500      case(Nodes):      case(Nodes):
501        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");  
502        break;        break;
503      case(Elements):      case(Elements):
504        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");  
505        break;        break;
506      case (FaceElements):      case (FaceElements):
507        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
508        break;        break;
509      case(Points):      case(Points):
510        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");  
511        break;        break;
512      case (ContactElementsOne):      case (ContactElementsOne):
513      case (ContactElementsZero):      case (ContactElementsZero):
514        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
515        break;        break;
516      case(DegreesOfFreedom):      case(DegreesOfFreedom):
517        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.");  
518        break;        break;
519      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
520        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.");  
521        break;        break;
522      default:      default:
523        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
524        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();
525          throw FinleyAdapterException(temp.str());
526        break;        break;
527    }    }
528    checkFinleyError();    checkFinleyError();
529  }  }
530    
531  //  //
532  // interpolates data to other domain:  // interpolates data to other domain:
533  //  //
534  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
535  {  {
536    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
537    if (targetDomain!=*this)    if (targetDomain!=*this)
538       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
539    
540    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();  
541  }  }
542    
543  //  //
544  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
545  //  //
546  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
547  {  {
548    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
549    if (argDomain!=*this)    if (argDomain!=*this)
# Line 669  void MeshAdapter::setToIntegrals(std::ve Line 552  void MeshAdapter::setToIntegrals(std::ve
552    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
553    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
554       case(Nodes):       case(Nodes):
555          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.");  
556          break;          break;
557       case(Elements):       case(Elements):
558          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 561  void MeshAdapter::setToIntegrals(std::ve
561          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
562          break;          break;
563       case(Points):       case(Points):
564          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.");  
565          break;          break;
566       case(ContactElementsZero):       case(ContactElementsZero):
567          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 570  void MeshAdapter::setToIntegrals(std::ve
570          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
571          break;          break;
572       case(DegreesOfFreedom):       case(DegreesOfFreedom):
573          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.");  
574          break;          break;
575       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
576          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.");  
577          break;          break;
578       default:       default:
579          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
580          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();
581            throw FinleyAdapterException(temp.str());
582          break;          break;
583    }    }
584    checkFinleyError();    checkFinleyError();
585  }  }
586    
587  //  //
588  // calculates the gradient of arg:  // calculates the gradient of arg:
589  //  //
590  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
591  {  {
592    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
593    if (argDomain!=*this)    if (argDomain!=*this)
# Line 716  void MeshAdapter::setToGradient(Data& gr Line 597  void MeshAdapter::setToGradient(Data& gr
597       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
598    
599    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
600      escriptDataC nodeDataC;
601    #ifdef PASO_MPI
602      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
603      if( arg.getFunctionSpace().getTypeCode() != Nodes )
604      {
605        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
606        nodeDataC = nodeTemp.getDataC();
607      }
608      else
609        nodeDataC = arg.getDataC();
610    #else
611      nodeDataC = arg.getDataC();
612    #endif
613    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
614         case(Nodes):         case(Nodes):
615            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");  
616            break;            break;
617         case(Elements):         case(Elements):
618            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
619            break;            break;
620         case(FaceElements):         case(FaceElements):
621            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
622            break;            break;
623         case(Points):         case(Points):
624            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
625            break;            break;
626         case(ContactElementsZero):         case(ContactElementsZero):
627            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
628            break;            break;
629         case(ContactElementsOne):         case(ContactElementsOne):
630            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
631            break;            break;
632         case(DegreesOfFreedom):         case(DegreesOfFreedom):
633            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.");  
634            break;            break;
635         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
636            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.");  
637            break;            break;
638         default:         default:
639            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
640            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();
641              throw FinleyAdapterException(temp.str());
642            break;            break;
643    }    }
644    checkFinleyError();    checkFinleyError();
645  }  }
646    
647  //  //
648  // returns the size of elements:  // returns the size of elements:
649  //  //
650  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
651  {  {
652    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
653    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
654    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
655         case(Nodes):         case(Nodes):
656            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");  
657            break;            break;
658         case(Elements):         case(Elements):
659            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
# Line 771  void MeshAdapter::setToSize(Data& size) Line 662  void MeshAdapter::setToSize(Data& size)
662            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
663            break;            break;
664         case(Points):         case(Points):
665            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
666            break;            break;
667         case(ContactElementsZero):         case(ContactElementsZero):
668         case(ContactElementsOne):         case(ContactElementsOne):
669            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
670            break;            break;
671         case(DegreesOfFreedom):         case(DegreesOfFreedom):
672            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.");  
673            break;            break;
674         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
675            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.");  
676            break;            break;
677         default:         default:
678            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
679            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();
680              throw FinleyAdapterException(temp.str());
681            break;            break;
682    }    }
683    checkFinleyError();    checkFinleyError();
684  }  }
685    
686  // sets the location of nodes:  // sets the location of nodes:
687  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
688  {  {
689    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
690      escriptDataC tmp;
691    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
692    if (newDomain!=*this)    if (newDomain!=*this)
693       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
694    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
695      Finley_Mesh_setCoordinates(mesh,&tmp);
696    checkFinleyError();    checkFinleyError();
697  }  }
698    
699  // saves a data array in openDX format:  // saves a data array in openDX format:
700  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
701  {  {
702    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
703    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
704      /* win32 refactor */
705      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
706      for(int i=0;i<num_data;i++)
707      {
708        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
709      }
710    
711      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
712      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
713      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
714    
715        boost::python::list keys=arg.keys();
716        for (int i=0;i<num_data;++i) {
717             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
718             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
719                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
720             data[i]=d.getDataC();
721             ptr_data[i]=&(data[i]);
722             std::string n=boost::python::extract<std::string>(keys[i]);
723             c_names[i]=&(names[i][0]);
724             if (n.length()>MAX_namelength-1) {
725                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
726                c_names[i][MAX_namelength-1]='\0';
727             } else {
728                strcpy(c_names[i],n.c_str());
729             }
730        }
731        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
732        checkFinleyError();
733        
734          /* win32 refactor */
735      TMPMEMFREE(c_names);
736      TMPMEMFREE(data);
737      TMPMEMFREE(ptr_data);
738      for(int i=0;i<num_data;i++)
739      {
740        TMPMEMFREE(names[i]);
741      }
742      TMPMEMFREE(names);
743    
744        return;
745  }  }
746    
747  // saves a data array in openVTK format:  // saves a data array in openVTK format:
748  void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
749  {  {
750    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
751    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
752      /* win32 refactor */
753      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
754      for(int i=0;i<num_data;i++)
755      {
756        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
757      }
758    
759      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
760      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
761      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
762    
763        boost::python::list keys=arg.keys();
764        for (int i=0;i<num_data;++i) {
765             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
766             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
767                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
768             data[i]=d.getDataC();
769             ptr_data[i]=&(data[i]);
770             std::string n=boost::python::extract<std::string>(keys[i]);
771             c_names[i]=&(names[i][0]);
772             if (n.length()>MAX_namelength-1) {
773                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
774                c_names[i][MAX_namelength-1]='\0';
775             } else {
776                strcpy(c_names[i],n.c_str());
777             }
778        }
779        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
780        checkFinleyError();
781      /* win32 refactor */
782      TMPMEMFREE(c_names);
783      TMPMEMFREE(data);
784      TMPMEMFREE(ptr_data);
785      for(int i=0;i<num_data;i++)
786      {
787        TMPMEMFREE(names[i]);
788      }
789      TMPMEMFREE(names);
790    
791        return;
792  }  }
793                                                                                                                                                                            
794                                                                                                                                                                            
795  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
796  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
797                        const int row_blocksize,                        const int row_blocksize,
# Line 849  SystemMatrixAdapter MeshAdapter::newSyst Line 826  SystemMatrixAdapter MeshAdapter::newSyst
826      }      }
827      // generate matrix:      // generate matrix:
828            
829      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);  
830      checkFinleyError();      checkFinleyError();
831        Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
832        checkPasoError();
833        Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
834      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
835  }  }
836    
837  //  //
838  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
839  // TODO:  // TODO:
840  //  //
841    
842  //  //
843  // 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:
844  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
# Line 877  bool MeshAdapter::isCellOriented(int fun Line 857  bool MeshAdapter::isCellOriented(int fun
857            return true;            return true;
858            break;            break;
859         default:         default:
860            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
861            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;
862              throw FinleyAdapterException(temp.str());
863            break;            break;
864    }    }
865    checkFinleyError();    checkFinleyError();
866    return false;    return false;
867  }  }
868    
869  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
870  {  {
871    switch(functionSpaceType_source) {    switch(functionSpaceType_source) {
# Line 899  bool MeshAdapter::probeInterpolationOnDo Line 881  bool MeshAdapter::probeInterpolationOnDo
881             case(ContactElementsOne):             case(ContactElementsOne):
882                 return true;                 return true;
883             default:             default:
884                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
885                 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;
886                   throw FinleyAdapterException(temp.str());
887          }          }
888          break;          break;
889       case(Elements):       case(Elements):
# Line 940  bool MeshAdapter::probeInterpolationOnDo Line 923  bool MeshAdapter::probeInterpolationOnDo
923             case(ContactElementsOne):             case(ContactElementsOne):
924                return true;                return true;
925             default:             default:
926               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
927               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;
928                 throw FinleyAdapterException(temp.str());
929          }          }
930          break;          break;
931       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
# Line 957  bool MeshAdapter::probeInterpolationOnDo Line 941  bool MeshAdapter::probeInterpolationOnDo
941            case(DegreesOfFreedom):            case(DegreesOfFreedom):
942               return false;               return false;
943            default:            default:
944               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
945               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;
946                 throw FinleyAdapterException(temp.str());
947         }         }
948         break;         break;
949       default:       default:
950          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
951          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;
952            throw FinleyAdapterException(temp.str());
953          break;          break;
954    }    }
955    checkFinleyError();    checkFinleyError();
956    return false;    return false;
957  }  }
958    
959  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
960  {  {
961      return false;      return false;
# Line 989  bool MeshAdapter::operator!=(const Abstr Line 976  bool MeshAdapter::operator!=(const Abstr
976    return !(operator==(other));    return !(operator==(other));
977  }  }
978    
979  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
980  {  {
981     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
982     checkFinleyError();     checkPasoError();
983     return out;     return out;
984  }  }
985  Data MeshAdapter::getX() const  
986    escript::Data MeshAdapter::getX() const
987  {  {
988    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
989  }  }
990  Data MeshAdapter::getNormal() const  
991    escript::Data MeshAdapter::getNormal() const
992  {  {
993    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
994  }  }
995  Data MeshAdapter::getSize() const  
996    escript::Data MeshAdapter::getSize() const
997  {  {
998    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
999  }  }
1000    
1001  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1002  {  {
1003    int* tagList;    int out=0;
1004    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1005    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1006    return tagList[sampleNo];    case(Nodes):
1007        out=mesh->Nodes->Tag[sampleNo];
1008        break;
1009      case(Elements):
1010        out=mesh->Elements->Tag[sampleNo];
1011        break;
1012      case(FaceElements):
1013        out=mesh->FaceElements->Tag[sampleNo];
1014        break;
1015      case(Points):
1016        out=mesh->Points->Tag[sampleNo];
1017        break;
1018      case(ContactElementsZero):
1019        out=mesh->ContactElements->Tag[sampleNo];
1020        break;
1021      case(ContactElementsOne):
1022        out=mesh->ContactElements->Tag[sampleNo];
1023        break;
1024      case(DegreesOfFreedom):
1025        break;
1026      case(ReducedDegreesOfFreedom):
1027        break;
1028      default:
1029        stringstream temp;
1030        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1031        throw FinleyAdapterException(temp.str());
1032        break;
1033      }
1034      return out;
1035  }  }
   
1036  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1037  {  {
1038    int* referenceNoList;    int out=0,i;
1039    int numReferenceNo;    Finley_Mesh* mesh=m_finleyMesh.get();
1040    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    switch (functionSpaceType) {
1041    return referenceNoList[sampleNo];    case(Nodes):
1042        if (mesh->Nodes!=NULL) {
1043          out=mesh->Nodes->Id[sampleNo];
1044          break;
1045        }
1046      case(Elements):
1047        out=mesh->Elements->Id[sampleNo];
1048        break;
1049      case(FaceElements):
1050        out=mesh->FaceElements->Id[sampleNo];
1051        break;
1052      case(Points):
1053        out=mesh->Points->Id[sampleNo];
1054        break;
1055      case(ContactElementsZero):
1056        out=mesh->ContactElements->Id[sampleNo];
1057        break;
1058      case(ContactElementsOne):
1059        out=mesh->ContactElements->Id[sampleNo];
1060        break;
1061      case(DegreesOfFreedom):
1062        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1063           if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
1064              out=mesh->Nodes->Id[i];
1065              break;
1066           }
1067        }
1068        break;
1069      case(ReducedDegreesOfFreedom):
1070        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1071           if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
1072              out=mesh->Nodes->Id[i];
1073              break;
1074           }
1075        }
1076        break;
1077      default:
1078        stringstream temp;
1079        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1080        throw FinleyAdapterException(temp.str());
1081        break;
1082      }
1083      return out;
1084  }  }
1085    
1086  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.147  
changed lines
  Added in v.757

  ViewVC Help
Powered by ViewVC 1.1.26