/[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 102 by jgs, Wed Dec 15 07:08:39 2004 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 798 by gross, Fri Aug 4 01:05:36 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    
 struct null_deleter  
 {  
   void operator()(void const *ptr) const  
   {  
   }  
 };  
   
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 64  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 90  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 150  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 171  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  }  }
162    
 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  
 {  
   //  
   // It is assumed the sampleNo has been checked  
   // before calling this function.  
   int* tagList;  
   int numTags;  
   getTagList(functionSpaceType, &tagList, &numTags);  
   return tagList[sampleNo];  
 }  
   
 //  
 // 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;  
   }  
   return;  
 }  
   
163  //  //
164  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
165  //  //
# Line 268  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 315  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()),  
269                         &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
270     checkFinleyError();     checkFinleyError();
271     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
272                    mat.getFinley_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
                   &(rhs.getDataC()),  
                                   &(d.getDataC()),&(y.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Mean_out);  
273     checkFinleyError();     checkFinleyError();
274     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
275                    mat.getFinley_SystemMatrix(),     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
                   &(rhs.getDataC()),  
                                   &(d_contact.getDataC()),  
                   &(y_contact.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Step_out);  
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     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));  
286       escriptDataC _rhs=rhs.getDataC();
287       escriptDataC _X=X.getDataC();
288       escriptDataC _Y=Y.getDataC();
289       escriptDataC _y=y.getDataC();
290       escriptDataC _y_contact=y_contact.getDataC();
291    
292       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
293     checkFinleyError();     checkFinleyError();
294     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),  
295                                    Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
 // cout << "Calling :addPDEToRHS." << endl;  
296     checkFinleyError();     checkFinleyError();
297     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),  
298                                    Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
 // cout << "Calling :addPDEToRHS." << endl;  
299     checkFinleyError();     checkFinleyError();
300  }  }
301    
302  //  //
303  // interpolates data between different function spaces:  // interpolates data between different function spaces:
304  //  //
305  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
306  {  {
307    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
308    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
309    if (inDomain!=*this)    if (inDomain!=*this)  
310       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
311    if (targetDomain!=*this)    if (targetDomain!=*this)
312       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
313    
314    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
315    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 415  void MeshAdapter::interpolateOnDomain(Da Line 336  void MeshAdapter::interpolateOnDomain(Da
336                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
337                 break;                 break;
338             default:             default:
339                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
340                 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();
341                   throw FinleyAdapterException(temp.str());
342                 break;                 break;
343          }          }
344          break;          break;
# Line 424  void MeshAdapter::interpolateOnDomain(Da Line 346  void MeshAdapter::interpolateOnDomain(Da
346          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
347             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
348          } else {          } else {
349             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
350          }          }
351          break;          break;
352       case(FaceElements):       case(FaceElements):
353          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
354             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
355          } else {          } else {
356             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.");  
357             break;             break;
358         }         }
359       case(Points):       case(Points):
360          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
361             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
362          } else {          } else {
363             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
364             break;             break;
365          }          }
366          break;          break;
# Line 450  void MeshAdapter::interpolateOnDomain(Da Line 369  void MeshAdapter::interpolateOnDomain(Da
369          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
370             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
371          } else {          } else {
372             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.");  
373             break;             break;
374          }          }
375          break;          break;
376       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
377          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
378             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
379             case(DegreesOfFreedom):             case(DegreesOfFreedom):
380             case(Nodes):             case(Nodes):
381                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
382                break;                break;
383    #ifndef PASO_MPI
384             case(Elements):             case(Elements):
385                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
386                break;                break;
# Line 475  void MeshAdapter::interpolateOnDomain(Da Line 394  void MeshAdapter::interpolateOnDomain(Da
394             case(ContactElementsOne):             case(ContactElementsOne):
395                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
396               break;               break;
397    #else
398               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
399               case(Elements):
400               {
401                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
402                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
403                  break;
404               }
405               case(FaceElements):
406               {
407                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
408                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
409                  break;
410               }
411               case(Points):
412               {
413                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
414                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
415                  break;
416               }
417               case(ContactElementsZero):
418               case(ContactElementsOne):
419               {
420                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
421                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
422                 break;
423               }
424    #endif
425             default:             default:
426               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
427               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();
428                 throw FinleyAdapterException(temp.str());
429               break;               break;
430          }          }
431          break;          break;
432       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
433         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
434            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
           case(Nodes):  
435               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
436               break;               break;
437            case(Elements):            case(Elements):
# Line 500  void MeshAdapter::interpolateOnDomain(Da Line 447  void MeshAdapter::interpolateOnDomain(Da
447            case(ContactElementsOne):            case(ContactElementsOne):
448               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
449               break;               break;
450              case(Nodes):
451                 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
452                 break;
453            case(DegreesOfFreedom):            case(DegreesOfFreedom):
454               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");  
455               break;               break;
456            default:            default:
457               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
458               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
459                 throw FinleyAdapterException(temp.str());
460               break;               break;
461         }         }
462         break;         break;
463       default:       default:
464          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
465          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();
466            throw FinleyAdapterException(temp.str());
467          break;          break;
468    }    }
469    checkFinleyError();    checkFinleyError();
# Line 521  void MeshAdapter::interpolateOnDomain(Da Line 472  void MeshAdapter::interpolateOnDomain(Da
472  //  //
473  // copies the locations of sample points into x:  // copies the locations of sample points into x:
474  //  //
475  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
476  {  {
477    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
478    if (argDomain!=*this)    if (argDomain!=*this)
# Line 531  void MeshAdapter::setToX(Data& arg) cons Line 482  void MeshAdapter::setToX(Data& arg) cons
482    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
483       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
484    } else {    } else {
485       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
486       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
487       // this is then interpolated onto arg:       // this is then interpolated onto arg:
488       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
489    }    }
490    checkFinleyError();    checkFinleyError();
491  }  }
492    
493  //  //
494  // 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:
495  //  //
496  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
497  {  {
498    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
499    if (normalDomain!=*this)    if (normalDomain!=*this)
# Line 549  void MeshAdapter::setToNormal(Data& norm Line 501  void MeshAdapter::setToNormal(Data& norm
501    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
502    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
503      case(Nodes):      case(Nodes):
504        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");  
505        break;        break;
506      case(Elements):      case(Elements):
507        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");  
508        break;        break;
509      case (FaceElements):      case (FaceElements):
510        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
511        break;        break;
512      case(Points):      case(Points):
513        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");  
514        break;        break;
515      case (ContactElementsOne):      case (ContactElementsOne):
516      case (ContactElementsZero):      case (ContactElementsZero):
517        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
518        break;        break;
519      case(DegreesOfFreedom):      case(DegreesOfFreedom):
520        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.");  
521        break;        break;
522      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
523        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.");  
524        break;        break;
525      default:      default:
526        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
527        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();
528          throw FinleyAdapterException(temp.str());
529        break;        break;
530    }    }
531    checkFinleyError();    checkFinleyError();
532  }  }
533    
534  //  //
535  // interpolates data to other domain:  // interpolates data to other domain:
536  //  //
537  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
538  {  {
539    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
540    if (targetDomain!=*this)    if (targetDomain!=*this)
541       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
542    
543    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();  
544  }  }
545    
546  //  //
547  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
548  //  //
549  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
550  {  {
551    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
552    if (argDomain!=*this)    if (argDomain!=*this)
# Line 607  void MeshAdapter::setToIntegrals(std::ve Line 555  void MeshAdapter::setToIntegrals(std::ve
555    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
556    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
557       case(Nodes):       case(Nodes):
558          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.");  
559          break;          break;
560       case(Elements):       case(Elements):
561          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
# Line 617  void MeshAdapter::setToIntegrals(std::ve Line 564  void MeshAdapter::setToIntegrals(std::ve
564          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
565          break;          break;
566       case(Points):       case(Points):
567          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.");  
568          break;          break;
569       case(ContactElementsZero):       case(ContactElementsZero):
570          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
# Line 627  void MeshAdapter::setToIntegrals(std::ve Line 573  void MeshAdapter::setToIntegrals(std::ve
573          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
574          break;          break;
575       case(DegreesOfFreedom):       case(DegreesOfFreedom):
576          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.");  
577          break;          break;
578       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
579          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.");  
580          break;          break;
581       default:       default:
582          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
583          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();
584            throw FinleyAdapterException(temp.str());
585          break;          break;
586    }    }
587    checkFinleyError();    checkFinleyError();
588  }  }
589    
590  //  //
591  // calculates the gradient of arg:  // calculates the gradient of arg:
592  //  //
593  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
594  {  {
595    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
596    if (argDomain!=*this)    if (argDomain!=*this)
# Line 654  void MeshAdapter::setToGradient(Data& gr Line 600  void MeshAdapter::setToGradient(Data& gr
600       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
601    
602    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
603      escriptDataC nodeDataC;
604    #ifdef PASO_MPI
605      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
606      if( arg.getFunctionSpace().getTypeCode() != Nodes )
607      {
608        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
609        nodeDataC = nodeTemp.getDataC();
610      }
611      else
612        nodeDataC = arg.getDataC();
613    #else
614      nodeDataC = arg.getDataC();
615    #endif
616    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
617         case(Nodes):         case(Nodes):
618            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");  
619            break;            break;
620         case(Elements):         case(Elements):
621            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
622            break;            break;
623         case(FaceElements):         case(FaceElements):
624            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
625            break;            break;
626         case(Points):         case(Points):
627            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
628            break;            break;
629         case(ContactElementsZero):         case(ContactElementsZero):
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(ContactElementsOne):         case(ContactElementsOne):
633            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
634            break;            break;
635         case(DegreesOfFreedom):         case(DegreesOfFreedom):
636            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.");  
637            break;            break;
638         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
639            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.");  
640            break;            break;
641         default:         default:
642            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
643            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();
644              throw FinleyAdapterException(temp.str());
645            break;            break;
646    }    }
647    checkFinleyError();    checkFinleyError();
648  }  }
649    
650  //  //
651  // returns the size of elements:  // returns the size of elements:
652  //  //
653  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
654  {  {
655    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
656    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
657    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
658         case(Nodes):         case(Nodes):
659            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");  
660            break;            break;
661         case(Elements):         case(Elements):
662            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
# Line 709  void MeshAdapter::setToSize(Data& size) Line 665  void MeshAdapter::setToSize(Data& size)
665            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
666            break;            break;
667         case(Points):         case(Points):
668            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
669            break;            break;
670         case(ContactElementsZero):         case(ContactElementsZero):
671         case(ContactElementsOne):         case(ContactElementsOne):
672            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
673            break;            break;
674         case(DegreesOfFreedom):         case(DegreesOfFreedom):
675            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.");  
676            break;            break;
677         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
678            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.");  
679            break;            break;
680         default:         default:
681            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
682            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();
683              throw FinleyAdapterException(temp.str());
684            break;            break;
685    }    }
686    checkFinleyError();    checkFinleyError();
687  }  }
688    
689  // sets the location of nodes:  // sets the location of nodes:
690  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
691  {  {
692    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
693      escriptDataC tmp;
694    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
695    if (newDomain!=*this)    if (newDomain!=*this)
696       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
697    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
698      Finley_Mesh_setCoordinates(mesh,&tmp);
699    checkFinleyError();    checkFinleyError();
700  }  }
701    
702  // saves a data array in openDX format:  // saves a data array in openDX format:
703  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
704  {  {
705    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
706    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
707      /* win32 refactor */
708      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
709      for(int i=0;i<num_data;i++)
710      {
711        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
712      }
713    
714      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
715      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
716      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
717    
718        boost::python::list keys=arg.keys();
719        for (int i=0;i<num_data;++i) {
720             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
721             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
722                 throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
723             data[i]=d.getDataC();
724             ptr_data[i]=&(data[i]);
725             std::string n=boost::python::extract<std::string>(keys[i]);
726             c_names[i]=&(names[i][0]);
727             if (n.length()>MAX_namelength-1) {
728                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
729                c_names[i][MAX_namelength-1]='\0';
730             } else {
731                strcpy(c_names[i],n.c_str());
732             }
733        }
734        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
735        checkFinleyError();
736        
737          /* win32 refactor */
738      TMPMEMFREE(c_names);
739      TMPMEMFREE(data);
740      TMPMEMFREE(ptr_data);
741      for(int i=0;i<num_data;i++)
742      {
743        TMPMEMFREE(names[i]);
744      }
745      TMPMEMFREE(names);
746    
747        return;
748  }  }
749    
750    // saves a data array in openVTK format:
751    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
752    {
753        int MAX_namelength=256;
754        const int num_data=boost::python::extract<int>(arg.attr("__len__")());
755      /* win32 refactor */
756      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
757      for(int i=0;i<num_data;i++)
758      {
759        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
760      }
761    
762      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
763      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
764      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
765    
766        boost::python::list keys=arg.keys();
767        for (int i=0;i<num_data;++i) {
768             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
769             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
770                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
771             data[i]=d.getDataC();
772             ptr_data[i]=&(data[i]);
773             std::string n=boost::python::extract<std::string>(keys[i]);
774             c_names[i]=&(names[i][0]);
775             if (n.length()>MAX_namelength-1) {
776                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
777                c_names[i][MAX_namelength-1]='\0';
778             } else {
779                strcpy(c_names[i],n.c_str());
780             }
781        }
782    #ifndef PASO_MPI    
783        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
784    #else
785        Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
786    #endif
787    
788    checkFinleyError();
789      /* win32 refactor */
790      TMPMEMFREE(c_names);
791      TMPMEMFREE(data);
792      TMPMEMFREE(ptr_data);
793      for(int i=0;i<num_data;i++)
794      {
795        TMPMEMFREE(names[i]);
796      }
797      TMPMEMFREE(names);
798    
799        return;
800    }
801                                                                                                                                                                            
802                                                                                                                                                                            
803  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
804  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
805                        const int row_blocksize,                        const int row_blocksize,
# Line 781  SystemMatrixAdapter MeshAdapter::newSyst Line 834  SystemMatrixAdapter MeshAdapter::newSyst
834      }      }
835      // generate matrix:      // generate matrix:
836            
837      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);  
838      checkFinleyError();      checkFinleyError();
839        Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
840        checkPasoError();
841        Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
842      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
843  }  }
844    
845  //  //
846  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
847  // TODO:  // TODO:
848  //  //
849    
850  //  //
851  // 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:
852  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
# Line 809  bool MeshAdapter::isCellOriented(int fun Line 865  bool MeshAdapter::isCellOriented(int fun
865            return true;            return true;
866            break;            break;
867         default:         default:
868            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
869            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;
870              throw FinleyAdapterException(temp.str());
871            break;            break;
872    }    }
873    checkFinleyError();    checkFinleyError();
874    return false;    return false;
875  }  }
876    
877  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
878  {  {
879    switch(functionSpaceType_source) {    switch(functionSpaceType_source) {
# Line 831  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(Elements):       case(Elements):
# Line 872  bool MeshAdapter::probeInterpolationOnDo Line 931  bool MeshAdapter::probeInterpolationOnDo
931             case(ContactElementsOne):             case(ContactElementsOne):
932                return true;                return true;
933             default:             default:
934               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
935               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;
936                 throw FinleyAdapterException(temp.str());
937          }          }
938          break;          break;
939       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
940         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
941            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
           case(Nodes):  
942            case(Elements):            case(Elements):
943            case(FaceElements):            case(FaceElements):
944            case(Points):            case(Points):
945            case(ContactElementsZero):            case(ContactElementsZero):
946            case(ContactElementsOne):            case(ContactElementsOne):
947                return true;                return true;
948              case(Nodes):
949            case(DegreesOfFreedom):            case(DegreesOfFreedom):
950               return false;               return false;
951            default:            default:
952               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
953               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;
954                 throw FinleyAdapterException(temp.str());
955         }         }
956         break;         break;
957       default:       default:
958          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
959          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;
960            throw FinleyAdapterException(temp.str());
961          break;          break;
962    }    }
963    checkFinleyError();    checkFinleyError();
964    return false;    return false;
965  }  }
966    
967  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
968  {  {
969      return false;      return false;
970  }  }
971  bool MeshAdapter::operator==(const MeshAdapter& other) const  
972    bool MeshAdapter::operator==(const AbstractDomain& other) const
973    {
974      const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
975      if (temp!=0) {
976        return (m_finleyMesh==temp->m_finleyMesh);
977      } else {
978        return false;
979      }
980    }
981    
982    bool MeshAdapter::operator!=(const AbstractDomain& other) const
983  {  {
984    return (m_finleyMesh==other.m_finleyMesh);    return !(operator==(other));
985  }  }
986    
987  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
988  {  {
989     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
990     checkFinleyError();     checkPasoError();
991     return out;     return out;
992  }  }
993  Data MeshAdapter::getX() const  
994    escript::Data MeshAdapter::getX() const
995  {  {
996    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
997  }  }
998  Data MeshAdapter::getNormal() const  
999    escript::Data MeshAdapter::getNormal() const
1000  {  {
1001    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1002  }  }
1003  Data MeshAdapter::getSize() const  
1004    escript::Data MeshAdapter::getSize() const
1005  {  {
1006    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
1007  }  }
1008    
1009    int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1010    {
1011      int out=0;
1012      Finley_Mesh* mesh=m_finleyMesh.get();
1013      switch (functionSpaceType) {
1014      case(Nodes):
1015        out=mesh->Nodes->Tag[sampleNo];
1016        break;
1017      case(Elements):
1018        out=mesh->Elements->Tag[sampleNo];
1019        break;
1020      case(FaceElements):
1021        out=mesh->FaceElements->Tag[sampleNo];
1022        break;
1023      case(Points):
1024        out=mesh->Points->Tag[sampleNo];
1025        break;
1026      case(ContactElementsZero):
1027        out=mesh->ContactElements->Tag[sampleNo];
1028        break;
1029      case(ContactElementsOne):
1030        out=mesh->ContactElements->Tag[sampleNo];
1031        break;
1032      case(DegreesOfFreedom):
1033        break;
1034      case(ReducedDegreesOfFreedom):
1035        break;
1036      default:
1037        stringstream temp;
1038        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1039        throw FinleyAdapterException(temp.str());
1040        break;
1041      }
1042      return out;
1043    }
1044    int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1045    {
1046      int out=0,i;
1047      Finley_Mesh* mesh=m_finleyMesh.get();
1048      switch (functionSpaceType) {
1049      case(Nodes):
1050        if (mesh->Nodes!=NULL) {
1051          out=mesh->Nodes->Id[sampleNo];
1052          break;
1053        }
1054      case(Elements):
1055        out=mesh->Elements->Id[sampleNo];
1056        break;
1057      case(FaceElements):
1058        out=mesh->FaceElements->Id[sampleNo];
1059        break;
1060      case(Points):
1061        out=mesh->Points->Id[sampleNo];
1062        break;
1063      case(ContactElementsZero):
1064        out=mesh->ContactElements->Id[sampleNo];
1065        break;
1066      case(ContactElementsOne):
1067        out=mesh->ContactElements->Id[sampleNo];
1068        break;
1069      case(DegreesOfFreedom):
1070        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1071           if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
1072              out=mesh->Nodes->Id[i];
1073              break;
1074           }
1075        }
1076        break;
1077      case(ReducedDegreesOfFreedom):
1078        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1079           if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
1080              out=mesh->Nodes->Id[i];
1081              break;
1082           }
1083        }
1084        break;
1085      default:
1086        stringstream temp;
1087        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1088        throw FinleyAdapterException(temp.str());
1089        break;
1090      }
1091      return out;
1092    }
1093    
1094  bool MeshAdapter::operator!=(const MeshAdapter& other) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1095  {  {
1096    return !operator==(other);    Finley_Mesh* mesh=m_finleyMesh.get();
1097      escriptDataC tmp=mask.getDataC();
1098      switch(functionSpaceType) {
1099           case(Nodes):
1100              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1101              break;
1102           case(DegreesOfFreedom):
1103              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1104              break;
1105           case(ReducedDegreesOfFreedom):
1106              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1107              break;
1108           case(Elements):
1109              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1110              break;
1111           case(FaceElements):
1112              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1113              break;
1114           case(Points):
1115              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1116              break;
1117           case(ContactElementsZero):
1118              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1119              break;
1120           case(ContactElementsOne):
1121              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1122              break;
1123           default:
1124              stringstream temp;
1125              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1126              throw FinleyAdapterException(temp.str());
1127      }
1128      checkFinleyError();
1129      return;
1130  }  }
1131  // bool MeshAdapter::operator==(const AbstractDomain& other) const  
1132  // {  
   // try {  
     // const MeshAdapter& ref = dynamic_cast<const MeshAdapter&>(other);  
     // return (operator==(ref));  
   // }  
   // catch (bad_cast) {  
     // return false;  
   // }  
 // }  
1133  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.102  
changed lines
  Added in v.798

  ViewVC Help
Powered by ViewVC 1.1.26