/[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 121 by jgs, Fri May 6 04:26:16 2005 UTC trunk/finley/src/CPPAdapter/MeshAdapter.cpp revision 817 by ksteube, Sat Aug 26 03:08:52 2006 UTC
# Line 12  Line 12 
12   *                                                                            *   *                                                                            *
13   ******************************************************************************   ******************************************************************************
14  */  */
15  extern "C" {  
16  #include "finley/finleyC/Finley.h"  #ifdef PASO_MPI
17  #include "finley/finleyC/Assemble.h"  #include <mpi.h>
18  #include "finley/finleyC/Mesh.h"  #endif
19  #include "finley/finleyC/Finley.h"  #include "MeshAdapter.h"
20  #include "finley/finleyC/System.h"  
21  }  #include "escript/Data.h"
22  #include "finley/CPPAdapter/SystemMatrixAdapter.h"  #include "escript/DataFactory.h"
 #include "finley/CPPAdapter/MeshAdapter.h"  
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/FinleyAdapterException.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/Data.h"  
 #include "escript/Data/DataArrayView.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/DataFactory.h"  
 #include <iostream>  
 #include <vector>  
 #include <sstream>  
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
26    
27  namespace finley {  namespace finley {
28    
 struct null_deleter  
 {  
   void operator()(void const *ptr) const  
   {  
   }  
 };  
   
29  //  //
30  // define the statics  // define the static constants
31  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;  MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
32  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;  const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
33  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;  const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
# Line 64  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 46  MeshAdapter::MeshAdapter(Finley_Mesh* fi
46    // for us.    // for us.
47    m_finleyMesh.reset(finleyMesh,null_deleter());    m_finleyMesh.reset(finleyMesh,null_deleter());
48  }  }
49    
50  //  //
51  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
52  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
# Line 90  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 73  Finley_Mesh* MeshAdapter::getFinley_Mesh
73    
74  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
75  {  {
76    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
77    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
78    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
79    checkFinleyError();    checkFinleyError();
80      TMPMEMFREE(fName);
81  }  }
82    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
83  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
84  {  {
85    return "FinleyMesh";    return "FinleyMesh";
# Line 150  int MeshAdapter::getContinuousFunctionCo Line 127  int MeshAdapter::getContinuousFunctionCo
127  {  {
128    return Nodes;    return Nodes;
129  }  }
130    
131  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
132  {  {
133    return Elements;    return Elements;
134  }  }
135    
136  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
137  {  {
138    return FaceElements;    return FaceElements;
139  }  }
140    
141  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
142  {  {
143    return ContactElementsZero;    return ContactElementsZero;
144  }  }
145    
146  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
147  {  {
148    return ContactElementsOne;    return ContactElementsOne;
# Line 171  int MeshAdapter::getSolutionCode() const Line 152  int MeshAdapter::getSolutionCode() const
152  {  {
153    return DegreesOfFreedom;    return DegreesOfFreedom;
154  }  }
155    
156  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
157  {  {
158    return ReducedDegreesOfFreedom;    return ReducedDegreesOfFreedom;
159  }  }
160    
161  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionCode() const
162  {  {
163    return Points;    return Points;
164  }  }
 //  
 // 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;  
 }  
165    
166  //  //
167  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
# Line 333  int MeshAdapter::getDim() const Line 172  int MeshAdapter::getDim() const
172    checkFinleyError();    checkFinleyError();
173    return numDim;    return numDim;
174  }  }
175    
176  //  //
177  // 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
178  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
# Line 380  pair<int,int> MeshAdapter::getDataShape( Line 220  pair<int,int> MeshAdapter::getDataShape(
220        case(DegreesOfFreedom):        case(DegreesOfFreedom):
221             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
222               numDataPointsPerSample=1;               numDataPointsPerSample=1;
223    #ifndef PASO_MPI
224               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
225    #else
226                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
227    #endif
228             }             }
229             break;             break;
230        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
231             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
232               numDataPointsPerSample=1;               numDataPointsPerSample=1;
233    #ifndef PASO_MPI
234               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
235    #else
236                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
237    #endif
238             }             }
239             break;             break;
240        default:        default:
241          stringstream temp;          stringstream temp;
242          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
243          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
244          break;          break;
245        }        }
246        return pair<int,int>(numDataPointsPerSample,numSamples);        return pair<int,int>(numDataPointsPerSample,numSamples);
247  }  }
248    
249  //  //
250  // 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:
251  //  //
252  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
253                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
254                       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,
255                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
256                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
257  {  {
258       escriptDataC _rhs=rhs.getDataC();
259       escriptDataC _A  =A.getDataC();
260       escriptDataC _B=B.getDataC();
261       escriptDataC _C=C.getDataC();
262       escriptDataC _D=D.getDataC();
263       escriptDataC _X=X.getDataC();
264       escriptDataC _Y=Y.getDataC();
265       escriptDataC _d=d.getDataC();
266       escriptDataC _y=y.getDataC();
267       escriptDataC _d_contact=d_contact.getDataC();
268       escriptDataC _y_contact=y_contact.getDataC();
269    
270     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
271     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),  
272                         &(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 );
273     checkFinleyError();     checkFinleyError();
274     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
275                    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);  
276     checkFinleyError();     checkFinleyError();
277     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
278                    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);  
279     checkFinleyError();     checkFinleyError();
280  }  }
281    
282  //  //
283  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
284  //  //
285  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  
286  {  {
287     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
288     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));  
289       escriptDataC _rhs=rhs.getDataC();
290       escriptDataC _X=X.getDataC();
291       escriptDataC _Y=Y.getDataC();
292       escriptDataC _y=y.getDataC();
293       escriptDataC _y_contact=y_contact.getDataC();
294    
295       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
296     checkFinleyError();     checkFinleyError();
297     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),  
298                                    Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
 // cout << "Calling :addPDEToRHS." << endl;  
299     checkFinleyError();     checkFinleyError();
300     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),  
301                                    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;  
302     checkFinleyError();     checkFinleyError();
303  }  }
304    
305  //  //
306  // interpolates data between different function spaces:  // interpolates data between different function spaces:
307  //  //
308  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
309  {  {
310    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
311    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
312    if (inDomain!=*this)    if (inDomain!=*this)  
313       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
314    if (targetDomain!=*this)    if (targetDomain!=*this)
315       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
316    
317    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
318    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 480  void MeshAdapter::interpolateOnDomain(Da Line 339  void MeshAdapter::interpolateOnDomain(Da
339                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
340                 break;                 break;
341             default:             default:
342                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
343                 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();
344                   throw FinleyAdapterException(temp.str());
345                 break;                 break;
346          }          }
347          break;          break;
# Line 489  void MeshAdapter::interpolateOnDomain(Da Line 349  void MeshAdapter::interpolateOnDomain(Da
349          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
350             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
351          } else {          } else {
352             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
353          }          }
354          break;          break;
355       case(FaceElements):       case(FaceElements):
356          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
357             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
358          } else {          } else {
359             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.");  
360             break;             break;
361         }         }
362       case(Points):       case(Points):
363          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
364             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
365          } else {          } else {
366             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
367             break;             break;
368          }          }
369          break;          break;
# Line 515  void MeshAdapter::interpolateOnDomain(Da Line 372  void MeshAdapter::interpolateOnDomain(Da
372          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
373             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
374          } else {          } else {
375             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.");  
376             break;             break;
377          }          }
378          break;          break;
379       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
380          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
381             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
382             case(DegreesOfFreedom):             case(DegreesOfFreedom):
383             case(Nodes):             case(Nodes):
384                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
385                break;                break;
386    #ifndef PASO_MPI
387             case(Elements):             case(Elements):
388                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
389                break;                break;
# Line 540  void MeshAdapter::interpolateOnDomain(Da Line 397  void MeshAdapter::interpolateOnDomain(Da
397             case(ContactElementsOne):             case(ContactElementsOne):
398                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
399               break;               break;
400    #else
401               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
402               case(Elements):
403               {
404                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
405                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
406                  break;
407               }
408               case(FaceElements):
409               {
410                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
411                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
412                  break;
413               }
414               case(Points):
415               {
416                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
417                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
418                  break;
419               }
420               case(ContactElementsZero):
421               case(ContactElementsOne):
422               {
423                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
424                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
425                 break;
426               }
427    #endif
428             default:             default:
429               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
430               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();
431                 throw FinleyAdapterException(temp.str());
432               break;               break;
433          }          }
434          break;          break;
435       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
436         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
437            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
           case(Nodes):  
438               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
439               break;               break;
440            case(Elements):            case(Elements):
# Line 565  void MeshAdapter::interpolateOnDomain(Da Line 450  void MeshAdapter::interpolateOnDomain(Da
450            case(ContactElementsOne):            case(ContactElementsOne):
451               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
452               break;               break;
453              case(Nodes):
454                 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
455                 break;
456            case(DegreesOfFreedom):            case(DegreesOfFreedom):
457               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");  
458               break;               break;
459            default:            default:
460               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
461               sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
462                 throw FinleyAdapterException(temp.str());
463               break;               break;
464         }         }
465         break;         break;
466       default:       default:
467          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
468          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();
469            throw FinleyAdapterException(temp.str());
470          break;          break;
471    }    }
472    checkFinleyError();    checkFinleyError();
# Line 586  void MeshAdapter::interpolateOnDomain(Da Line 475  void MeshAdapter::interpolateOnDomain(Da
475  //  //
476  // copies the locations of sample points into x:  // copies the locations of sample points into x:
477  //  //
478  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
479  {  {
480    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
481    if (argDomain!=*this)    if (argDomain!=*this)
# Line 596  void MeshAdapter::setToX(Data& arg) cons Line 485  void MeshAdapter::setToX(Data& arg) cons
485    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
486       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
487    } else {    } else {
488       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
489       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
490       // this is then interpolated onto arg:       // this is then interpolated onto arg:
491       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
492    }    }
493    checkFinleyError();    checkFinleyError();
494  }  }
495    
496  //  //
497  // 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:
498  //  //
499  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
500  {  {
501    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
502    if (normalDomain!=*this)    if (normalDomain!=*this)
# Line 614  void MeshAdapter::setToNormal(Data& norm Line 504  void MeshAdapter::setToNormal(Data& norm
504    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
505    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
506      case(Nodes):      case(Nodes):
507        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");  
508        break;        break;
509      case(Elements):      case(Elements):
510        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");  
511        break;        break;
512      case (FaceElements):      case (FaceElements):
513        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
514        break;        break;
515      case(Points):      case(Points):
516        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");  
517        break;        break;
518      case (ContactElementsOne):      case (ContactElementsOne):
519      case (ContactElementsZero):      case (ContactElementsZero):
520        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
521        break;        break;
522      case(DegreesOfFreedom):      case(DegreesOfFreedom):
523        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.");  
524        break;        break;
525      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
526        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.");  
527        break;        break;
528      default:      default:
529        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
530        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();
531          throw FinleyAdapterException(temp.str());
532        break;        break;
533    }    }
534    checkFinleyError();    checkFinleyError();
535  }  }
536    
537  //  //
538  // interpolates data to other domain:  // interpolates data to other domain:
539  //  //
540  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
541  {  {
542    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
543    if (targetDomain!=*this)    if (targetDomain!=*this)
544       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
545    
546    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();  
547  }  }
548    
549  //  //
550  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
551  //  //
552  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
553  {  {
554    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
555    if (argDomain!=*this)    if (argDomain!=*this)
# Line 672  void MeshAdapter::setToIntegrals(std::ve Line 558  void MeshAdapter::setToIntegrals(std::ve
558    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
559    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
560       case(Nodes):       case(Nodes):
561          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.");  
562          break;          break;
563       case(Elements):       case(Elements):
564          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
# Line 682  void MeshAdapter::setToIntegrals(std::ve Line 567  void MeshAdapter::setToIntegrals(std::ve
567          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
568          break;          break;
569       case(Points):       case(Points):
570          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.");  
571          break;          break;
572       case(ContactElementsZero):       case(ContactElementsZero):
573          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
# Line 692  void MeshAdapter::setToIntegrals(std::ve Line 576  void MeshAdapter::setToIntegrals(std::ve
576          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
577          break;          break;
578       case(DegreesOfFreedom):       case(DegreesOfFreedom):
579          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.");  
580          break;          break;
581       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
582          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.");  
583          break;          break;
584       default:       default:
585          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
586          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();
587            throw FinleyAdapterException(temp.str());
588          break;          break;
589    }    }
590    checkFinleyError();    checkFinleyError();
591  }  }
592    
593  //  //
594  // calculates the gradient of arg:  // calculates the gradient of arg:
595  //  //
596  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
597  {  {
598    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
599    if (argDomain!=*this)    if (argDomain!=*this)
# Line 719  void MeshAdapter::setToGradient(Data& gr Line 603  void MeshAdapter::setToGradient(Data& gr
603       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
604    
605    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
606      escriptDataC nodeDataC;
607    #ifdef PASO_MPI
608      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
609      if( arg.getFunctionSpace().getTypeCode() != Nodes )
610      {
611        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
612        nodeDataC = nodeTemp.getDataC();
613      }
614      else
615        nodeDataC = arg.getDataC();
616    #else
617      nodeDataC = arg.getDataC();
618    #endif
619    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
620         case(Nodes):         case(Nodes):
621            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");  
622            break;            break;
623         case(Elements):         case(Elements):
624            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
625            break;            break;
626         case(FaceElements):         case(FaceElements):
627            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
628            break;            break;
629         case(Points):         case(Points):
630            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
631            break;            break;
632         case(ContactElementsZero):         case(ContactElementsZero):
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(ContactElementsOne):         case(ContactElementsOne):
636            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
637            break;            break;
638         case(DegreesOfFreedom):         case(DegreesOfFreedom):
639            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.");  
640            break;            break;
641         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
642            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.");  
643            break;            break;
644         default:         default:
645            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
646            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();
647              throw FinleyAdapterException(temp.str());
648            break;            break;
649    }    }
650    checkFinleyError();    checkFinleyError();
651  }  }
652    
653  //  //
654  // returns the size of elements:  // returns the size of elements:
655  //  //
656  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
657  {  {
658    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
659    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
660    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
661         case(Nodes):         case(Nodes):
662            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
           sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");  
663            break;            break;
664         case(Elements):         case(Elements):
665            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
# Line 774  void MeshAdapter::setToSize(Data& size) Line 668  void MeshAdapter::setToSize(Data& size)
668            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
669            break;            break;
670         case(Points):         case(Points):
671            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
672            break;            break;
673         case(ContactElementsZero):         case(ContactElementsZero):
674         case(ContactElementsOne):         case(ContactElementsOne):
675            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
676            break;            break;
677         case(DegreesOfFreedom):         case(DegreesOfFreedom):
678            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.");  
679            break;            break;
680         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
681            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.");  
682            break;            break;
683         default:         default:
684            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
685            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();
686              throw FinleyAdapterException(temp.str());
687            break;            break;
688    }    }
689    checkFinleyError();    checkFinleyError();
690  }  }
691    
692  // sets the location of nodes:  // sets the location of nodes:
693  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
694  {  {
695    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
696      escriptDataC tmp;
697    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
698    if (newDomain!=*this)    if (newDomain!=*this)
699       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
700    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
701      Finley_Mesh_setCoordinates(mesh,&tmp);
702    checkFinleyError();    checkFinleyError();
703  }  }
704    
705  // saves a data array in openDX format:  // saves a data array in openDX format:
706  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
707  {  {
708    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
709    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
710      /* win32 refactor */
711      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
712      for(int i=0;i<num_data;i++)
713      {
714        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
715      }
716    
717      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
718      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
719      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
720    
721        boost::python::list keys=arg.keys();
722        for (int i=0;i<num_data;++i) {
723             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
724             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
725                 throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
726             data[i]=d.getDataC();
727             ptr_data[i]=&(data[i]);
728             std::string n=boost::python::extract<std::string>(keys[i]);
729             c_names[i]=&(names[i][0]);
730             if (n.length()>MAX_namelength-1) {
731                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
732                c_names[i][MAX_namelength-1]='\0';
733             } else {
734                strcpy(c_names[i],n.c_str());
735             }
736        }
737        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
738        checkFinleyError();
739        
740          /* win32 refactor */
741      TMPMEMFREE(c_names);
742      TMPMEMFREE(data);
743      TMPMEMFREE(ptr_data);
744      for(int i=0;i<num_data;i++)
745      {
746        TMPMEMFREE(names[i]);
747      }
748      TMPMEMFREE(names);
749    
750        return;
751  }  }
752    
753  // saves a data array in openVTK format:  // saves a data array in openVTK format:
754  void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
755  {  {
756    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
757    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
758      /* win32 refactor */
759      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
760      for(int i=0;i<num_data;i++)
761      {
762        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
763      }
764    
765      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
766      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
767      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
768    
769        boost::python::list keys=arg.keys();
770        for (int i=0;i<num_data;++i) {
771             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
772             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
773                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
774             data[i]=d.getDataC();
775             ptr_data[i]=&(data[i]);
776             std::string n=boost::python::extract<std::string>(keys[i]);
777             c_names[i]=&(names[i][0]);
778             if (n.length()>MAX_namelength-1) {
779                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
780                c_names[i][MAX_namelength-1]='\0';
781             } else {
782                strcpy(c_names[i],n.c_str());
783             }
784        }
785    #ifndef PASO_MPI    
786        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
787    #else
788        Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
789    #endif
790    
791    checkFinleyError();
792      /* win32 refactor */
793      TMPMEMFREE(c_names);
794      TMPMEMFREE(data);
795      TMPMEMFREE(ptr_data);
796      for(int i=0;i<num_data;i++)
797      {
798        TMPMEMFREE(names[i]);
799      }
800      TMPMEMFREE(names);
801    
802        return;
803  }  }
804                                                                                                                                                                            
805                                                                                                                                                                            
806  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
807  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
808                        const int row_blocksize,                        const int row_blocksize,
# Line 852  SystemMatrixAdapter MeshAdapter::newSyst Line 837  SystemMatrixAdapter MeshAdapter::newSyst
837      }      }
838      // generate matrix:      // generate matrix:
839            
840      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);  
841      checkFinleyError();      checkFinleyError();
842        Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
843        checkPasoError();
844        Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
845      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
846  }  }
847    
848  //  //
849  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
850  // TODO:  // TODO:
851  //  //
852    
853  //  //
854  // 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:
855  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
# Line 880  bool MeshAdapter::isCellOriented(int fun Line 868  bool MeshAdapter::isCellOriented(int fun
868            return true;            return true;
869            break;            break;
870         default:         default:
871            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
872            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;
873              throw FinleyAdapterException(temp.str());
874            break;            break;
875    }    }
876    checkFinleyError();    checkFinleyError();
877    return false;    return false;
878  }  }
879    
880  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
881  {  {
882    switch(functionSpaceType_source) {    switch(functionSpaceType_source) {
# Line 902  bool MeshAdapter::probeInterpolationOnDo Line 892  bool MeshAdapter::probeInterpolationOnDo
892             case(ContactElementsOne):             case(ContactElementsOne):
893                 return true;                 return true;
894             default:             default:
895                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
896                 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;
897                   throw FinleyAdapterException(temp.str());
898          }          }
899          break;          break;
900       case(Elements):       case(Elements):
# Line 943  bool MeshAdapter::probeInterpolationOnDo Line 934  bool MeshAdapter::probeInterpolationOnDo
934             case(ContactElementsOne):             case(ContactElementsOne):
935                return true;                return true;
936             default:             default:
937               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
938               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;
939                 throw FinleyAdapterException(temp.str());
940          }          }
941          break;          break;
942       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
943         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
944            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
           case(Nodes):  
945            case(Elements):            case(Elements):
946            case(FaceElements):            case(FaceElements):
947            case(Points):            case(Points):
948            case(ContactElementsZero):            case(ContactElementsZero):
949            case(ContactElementsOne):            case(ContactElementsOne):
950                return true;                return true;
951              case(Nodes):
952            case(DegreesOfFreedom):            case(DegreesOfFreedom):
953               return false;               return false;
954            default:            default:
955               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
956               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;
957                 throw FinleyAdapterException(temp.str());
958         }         }
959         break;         break;
960       default:       default:
961          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
962          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;
963            throw FinleyAdapterException(temp.str());
964          break;          break;
965    }    }
966    checkFinleyError();    checkFinleyError();
967    return false;    return false;
968  }  }
969    
970  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
971  {  {
972      return false;      return false;
# Line 992  bool MeshAdapter::operator!=(const Abstr Line 987  bool MeshAdapter::operator!=(const Abstr
987    return !(operator==(other));    return !(operator==(other));
988  }  }
989    
990  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
991  {  {
992     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
993     checkFinleyError();     checkPasoError();
994     return out;     return out;
995  }  }
996  Data MeshAdapter::getX() const  
997    escript::Data MeshAdapter::getX() const
998  {  {
999    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
1000  }  }
1001  Data MeshAdapter::getNormal() const  
1002    escript::Data MeshAdapter::getNormal() const
1003  {  {
1004    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1005  }  }
1006  Data MeshAdapter::getSize() const  
1007    escript::Data MeshAdapter::getSize() const
1008  {  {
1009    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
1010  }  }
1011    
1012  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1013  {  {
1014    int* tagList;    int out=0;
1015    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1016    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1017    return tagList[sampleNo];    case(Nodes):
1018        out=mesh->Nodes->Tag[sampleNo];
1019        break;
1020      case(Elements):
1021        out=mesh->Elements->Tag[sampleNo];
1022        break;
1023      case(FaceElements):
1024        out=mesh->FaceElements->Tag[sampleNo];
1025        break;
1026      case(Points):
1027        out=mesh->Points->Tag[sampleNo];
1028        break;
1029      case(ContactElementsZero):
1030        out=mesh->ContactElements->Tag[sampleNo];
1031        break;
1032      case(ContactElementsOne):
1033        out=mesh->ContactElements->Tag[sampleNo];
1034        break;
1035      case(DegreesOfFreedom):
1036        break;
1037      case(ReducedDegreesOfFreedom):
1038        break;
1039      default:
1040        stringstream temp;
1041        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1042        throw FinleyAdapterException(temp.str());
1043        break;
1044      }
1045      return out;
1046    }
1047    int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1048    {
1049      int out=0,i;
1050      Finley_Mesh* mesh=m_finleyMesh.get();
1051      switch (functionSpaceType) {
1052      case(Nodes):
1053        if (mesh->Nodes!=NULL) {
1054          out=mesh->Nodes->Id[sampleNo];
1055          break;
1056        }
1057      case(Elements):
1058        out=mesh->Elements->Id[sampleNo];
1059        break;
1060      case(FaceElements):
1061        out=mesh->FaceElements->Id[sampleNo];
1062        break;
1063      case(Points):
1064        out=mesh->Points->Id[sampleNo];
1065        break;
1066      case(ContactElementsZero):
1067        out=mesh->ContactElements->Id[sampleNo];
1068        break;
1069      case(ContactElementsOne):
1070        out=mesh->ContactElements->Id[sampleNo];
1071        break;
1072      case(DegreesOfFreedom):
1073        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1074           if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
1075              out=mesh->Nodes->Id[i];
1076              break;
1077           }
1078        }
1079        break;
1080      case(ReducedDegreesOfFreedom):
1081        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1082           if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
1083              out=mesh->Nodes->Id[i];
1084              break;
1085           }
1086        }
1087        break;
1088      default:
1089        stringstream temp;
1090        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1091        throw FinleyAdapterException(temp.str());
1092        break;
1093      }
1094      return out;
1095  }  }
1096    
1097  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1098  {  {
1099    int* referenceNoList;    Finley_Mesh* mesh=m_finleyMesh.get();
1100    int numReferenceNo;    escriptDataC tmp=mask.getDataC();
1101    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    switch(functionSpaceType) {
1102    return referenceNoList[sampleNo];         case(Nodes):
1103              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1104              break;
1105           case(DegreesOfFreedom):
1106              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1107              break;
1108           case(ReducedDegreesOfFreedom):
1109              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1110              break;
1111           case(Elements):
1112              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1113              break;
1114           case(FaceElements):
1115              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1116              break;
1117           case(Points):
1118              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1119              break;
1120           case(ContactElementsZero):
1121              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1122              break;
1123           case(ContactElementsOne):
1124              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1125              break;
1126           default:
1127              stringstream temp;
1128              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1129              throw FinleyAdapterException(temp.str());
1130      }
1131      checkFinleyError();
1132      return;
1133  }  }
1134    
1135    
1136  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.121  
changed lines
  Added in v.817

  ViewVC Help
Powered by ViewVC 1.1.26