/[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 1064 by gross, Tue Mar 27 06:21:02 2007 UTC
# Line 12  Line 12 
12   *                                                                            *   *                                                                            *
13   ******************************************************************************   ******************************************************************************
14  */  */
15  extern "C" {  
16  #include "finley/finleyC/Finley.h"  #ifdef PASO_MPI
17  #include "finley/finleyC/Assemble.h"  #include <mpi.h>
18  #include "finley/finleyC/Mesh.h"  #endif
19  #include "finley/finleyC/Finley.h"  #include "MeshAdapter.h"
20  #include "finley/finleyC/System.h"  
21  }  #include "escript/Data.h"
22  #include "finley/CPPAdapter/SystemMatrixAdapter.h"  #include "escript/DataFactory.h"
 #include "finley/CPPAdapter/MeshAdapter.h"  
 #include "finley/CPPAdapter/FinleyError.h"  
 #include "finley/CPPAdapter/FinleyAdapterException.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/Data.h"  
 #include "escript/Data/DataArrayView.h"  
 #include "escript/Data/FunctionSpace.h"  
 #include "escript/Data/DataFactory.h"  
 #include <iostream>  
 #include <vector>  
 #include <sstream>  
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
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;
34  const int MeshAdapter::Nodes=FINLEY_NODES;  const int MeshAdapter::Nodes=FINLEY_NODES;
35    const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
36  const int MeshAdapter::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=FINLEY_ELEMENTS;
37    const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
38  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
39    const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
40  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=FINLEY_POINTS;
41  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
42    const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
43  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
44    const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
45    
46  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
47  {  {
# Line 64  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 51  MeshAdapter::MeshAdapter(Finley_Mesh* fi
51    // for us.    // for us.
52    m_finleyMesh.reset(finleyMesh,null_deleter());    m_finleyMesh.reset(finleyMesh,null_deleter());
53  }  }
54    
55  //  //
56  // The copy constructor should just increment the use count  // The copy constructor should just increment the use count
57  MeshAdapter::MeshAdapter(const MeshAdapter& in):  MeshAdapter::MeshAdapter(const MeshAdapter& in):
# Line 90  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 78  Finley_Mesh* MeshAdapter::getFinley_Mesh
78    
79  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
80  {  {
81    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
82    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
83    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
84    checkFinleyError();    checkFinleyError();
85      TMPMEMFREE(fName);
86  }  }
87    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
88  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
89  {  {
90    return "FinleyMesh";    return "FinleyMesh";
# Line 135  void MeshAdapter::setFunctionSpaceTypeNa Line 117  void MeshAdapter::setFunctionSpaceTypeNa
117    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
118      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));      (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
119    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
120        (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
121      m_functionSpaceTypeNames.insert
122      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
123    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
124        (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
125      m_functionSpaceTypeNames.insert
126      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
127    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
128        (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
129      m_functionSpaceTypeNames.insert
130      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
131    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
132      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
133    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
134        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
135      m_functionSpaceTypeNames.insert
136      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
137      m_functionSpaceTypeNames.insert
138        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
139  }  }
140    
141  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
142  {  {
143    return Nodes;    return Nodes;
144  }  }
145    int MeshAdapter::getReducedContinuousFunctionCode() const
146    {
147      return ReducedNodes;
148    }
149    
150  int MeshAdapter::getFunctionCode() const  int MeshAdapter::getFunctionCode() const
151  {  {
152    return Elements;    return Elements;
153  }  }
154    int MeshAdapter::getReducedFunctionCode() const
155    {
156      return ReducedElements;
157    }
158    
159  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
160  {  {
161    return FaceElements;    return FaceElements;
162  }  }
163    int MeshAdapter::getReducedFunctionOnBoundaryCode() const
164    {
165      return ReducedFaceElements;
166    }
167    
168  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
169  {  {
170    return ContactElementsZero;    return ContactElementsZero;
171  }  }
172    int MeshAdapter::getReducedFunctionOnContactZeroCode() const
173    {
174      return ReducedContactElementsZero;
175    }
176    
177  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
178  {  {
179    return ContactElementsOne;    return ContactElementsOne;
180  }  }
181    int MeshAdapter::getReducedFunctionOnContactOneCode() const
182    {
183      return ReducedContactElementsOne;
184    }
185    
186  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
187  {  {
188    return DegreesOfFreedom;    return DegreesOfFreedom;
189  }  }
190    
191  int MeshAdapter::getReducedSolutionCode() const  int MeshAdapter::getReducedSolutionCode() const
192  {  {
193    return ReducedDegreesOfFreedom;    return ReducedDegreesOfFreedom;
194  }  }
195    
196  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionCode() const
197  {  {
198    return Points;    return Points;
199  }  }
200    
 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;  
 }  
   
201  //  //
202  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
203  //  //
# Line 268  int MeshAdapter::getDim() const Line 207  int MeshAdapter::getDim() const
207    checkFinleyError();    checkFinleyError();
208    return numDim;    return numDim;
209  }  }
210    
211  //  //
212  // 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
213  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
# Line 282  pair<int,int> MeshAdapter::getDataShape( Line 222  pair<int,int> MeshAdapter::getDataShape(
222             numDataPointsPerSample=1;             numDataPointsPerSample=1;
223             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
224             break;             break;
225          case(ReducedNodes):
226               /* TODO: add ReducedNodes */
227               numDataPointsPerSample=1;
228               if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
229               throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");
230               break;
231        case(Elements):        case(Elements):
232             if (mesh->Elements!=NULL) {             if (mesh->Elements!=NULL) {
233               numSamples=mesh->Elements->numElements;               numSamples=mesh->Elements->numElements;
234               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;               numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
235             }             }
236             break;             break;
237          case(ReducedElements):
238               if (mesh->Elements!=NULL) {
239                 numSamples=mesh->Elements->numElements;
240                 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
241               }
242               break;
243        case(FaceElements):        case(FaceElements):
244             if (mesh->FaceElements!=NULL) {             if (mesh->FaceElements!=NULL) {
245                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;                  numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
246                  numSamples=mesh->FaceElements->numElements;                  numSamples=mesh->FaceElements->numElements;
247             }             }
248             break;             break;
249          case(ReducedFaceElements):
250               if (mesh->FaceElements!=NULL) {
251                    numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
252                    numSamples=mesh->FaceElements->numElements;
253               }
254               break;
255        case(Points):        case(Points):
256             if (mesh->Points!=NULL) {             if (mesh->Points!=NULL) {
257               numDataPointsPerSample=1;               numDataPointsPerSample=1;
# Line 306  pair<int,int> MeshAdapter::getDataShape( Line 264  pair<int,int> MeshAdapter::getDataShape(
264               numSamples=mesh->ContactElements->numElements;               numSamples=mesh->ContactElements->numElements;
265             }             }
266             break;             break;
267          case(ReducedContactElementsZero):
268               if (mesh->ContactElements!=NULL) {
269                 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
270                 numSamples=mesh->ContactElements->numElements;
271               }
272               break;
273        case(ContactElementsOne):        case(ContactElementsOne):
274             if (mesh->ContactElements!=NULL) {             if (mesh->ContactElements!=NULL) {
275               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;               numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
276               numSamples=mesh->ContactElements->numElements;               numSamples=mesh->ContactElements->numElements;
277             }             }
278             break;             break;
279          case(ReducedContactElementsOne):
280               if (mesh->ContactElements!=NULL) {
281                 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
282                 numSamples=mesh->ContactElements->numElements;
283               }
284               break;
285        case(DegreesOfFreedom):        case(DegreesOfFreedom):
286             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
287               numDataPointsPerSample=1;               numDataPointsPerSample=1;
288    #ifndef PASO_MPI
289               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
290    #else
291                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
292    #endif
293             }             }
294             break;             break;
295        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
296             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
297               numDataPointsPerSample=1;               numDataPointsPerSample=1;
298    #ifndef PASO_MPI
299               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
300    #else
301                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
302    #endif
303             }             }
304             break;             break;
305        default:        default:
306          stringstream temp;          stringstream temp;
307          temp << "Error - Invalid function space type: "          temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
          << functionSpaceCode << " for domain: " << getDescription();  
308          throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
309          break;          break;
310        }        }
311        return pair<int,int>(numDataPointsPerSample,numSamples);        return pair<int,int>(numDataPointsPerSample,numSamples);
312  }  }
313    
314  //  //
315  // 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:
316  //  //
317  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
318                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
319                       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,
320                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
321                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
322  {  {
323       escriptDataC _rhs=rhs.getDataC();
324       escriptDataC _A  =A.getDataC();
325       escriptDataC _B=B.getDataC();
326       escriptDataC _C=C.getDataC();
327       escriptDataC _D=D.getDataC();
328       escriptDataC _X=X.getDataC();
329       escriptDataC _Y=Y.getDataC();
330       escriptDataC _d=d.getDataC();
331       escriptDataC _y=y.getDataC();
332       escriptDataC _d_contact=d_contact.getDataC();
333       escriptDataC _y_contact=y_contact.getDataC();
334    
335     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
336     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),  
337                         &(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 );
338     checkFinleyError();     checkFinleyError();
339     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
340                    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);  
341     checkFinleyError();     checkFinleyError();
342     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
343                    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);  
344     checkFinleyError();     checkFinleyError();
345  }  }
346    
347  //  //
348  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
349  //  //
350  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  
351  {  {
352     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
353     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));  
354       escriptDataC _rhs=rhs.getDataC();
355       escriptDataC _X=X.getDataC();
356       escriptDataC _Y=Y.getDataC();
357       escriptDataC _y=y.getDataC();
358       escriptDataC _y_contact=y_contact.getDataC();
359    
360       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
361     checkFinleyError();     checkFinleyError();
362     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),  
363                                    Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
 // cout << "Calling :addPDEToRHS." << endl;  
364     checkFinleyError();     checkFinleyError();
365     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),  
366                                    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;  
367     checkFinleyError();     checkFinleyError();
368  }  }
369    
370  //  //
371  // interpolates data between different function spaces:  // interpolates data between different function spaces:
372  //  //
373  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
374  {  {
375    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
376    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
377    if (inDomain!=*this)    if (inDomain!=*this)  
378       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
379    if (targetDomain!=*this)    if (targetDomain!=*this)
380       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
381    
382    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
383      escriptDataC _target=target.getDataC();
384      escriptDataC _in=in.getDataC();
385    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
386       case(Nodes):       case(Nodes):
387          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
388             case(Nodes):             case(Nodes):
389               case(ReducedNodes):
390               case(DegreesOfFreedom):
391             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
392                   Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
393                   break;
394               case(Elements):
395               case(ReducedElements):
396                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
397                   break;
398               case(FaceElements):
399               case(ReducedFaceElements):
400                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
401                   break;
402               case(Points):
403                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
404                   break;
405               case(ContactElementsZero):
406               case(ReducedContactElementsZero):
407                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
408                   break;
409               case(ContactElementsOne):
410               case(ReducedContactElementsOne):
411                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
412                   break;
413               default:
414                   stringstream temp;
415                   temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
416                   throw FinleyAdapterException(temp.str());
417                   break;
418            }
419            break;
420         case(ReducedNodes):
421            switch(target.getFunctionSpace().getTypeCode()) {
422               case(Nodes):
423               case(ReducedNodes):
424             case(DegreesOfFreedom):             case(DegreesOfFreedom):
425                 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));             case(ReducedDegreesOfFreedom):
426                   Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
427                 break;                 break;
428             case(Elements):             case(Elements):
429                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));             case(ReducedElements):
430                   Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
431                 break;                 break;
432             case(FaceElements):             case(FaceElements):
433                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedFaceElements):
434                   Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
435                 break;                 break;
436             case(Points):             case(Points):
437                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
438                 break;                 break;
439             case(ContactElementsZero):             case(ContactElementsZero):
440                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsZero):
441                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
442                 break;                 break;
443             case(ContactElementsOne):             case(ContactElementsOne):
444                 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsOne):
445                   Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
446                 break;                 break;
447             default:             default:
448                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
449                 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();
450                   throw FinleyAdapterException(temp.str());
451                 break;                 break;
452          }          }
453          break;          break;
454       case(Elements):       case(Elements):
455          if (target.getFunctionSpace().getTypeCode()==Elements) {          if (target.getFunctionSpace().getTypeCode()==Elements) {
456             Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
457            } else {
458               throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
459            }
460            break;
461         case(ReducedElements):
462            if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
463               Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
464          } else {          } else {
465             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");  
466          }          }
467          break;          break;
468       case(FaceElements):       case(FaceElements):
469          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
470             Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
471          } else {          } else {
472             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
473             sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");         }
474             break;         break;
475         case(ReducedFaceElements):
476            if (target.getFunctionSpace().getTypeCode()==FaceElements) {
477               Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
478            } else {
479               throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
480         }         }
481           break;
482       case(Points):       case(Points):
483          if (target.getFunctionSpace().getTypeCode()==Points) {          if (target.getFunctionSpace().getTypeCode()==Points) {
484             Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
485          } else {          } else {
486             Finley_ErrorCode=TYPE_ERROR;             throw FinleyAdapterException("Error - No interpolation with data on points possible.");
            sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");  
            break;  
487          }          }
488          break;          break;
489       case(ContactElementsZero):       case(ContactElementsZero):
490       case(ContactElementsOne):       case(ContactElementsOne):
491          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
492             Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
493          } else {          } else {
494             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.");  
495             break;             break;
496          }          }
497          break;          break;
498       case(DegreesOfFreedom):       case(ReducedContactElementsZero):
499         case(ReducedContactElementsOne):
500            if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
501               Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
502            } else {
503               throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
504               break;
505            }
506            break;
507         case(DegreesOfFreedom):      
508          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
509             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
510             case(DegreesOfFreedom):             case(DegreesOfFreedom):
511             case(Nodes):             case(Nodes):
512                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));             case(ReducedNodes):
513                  Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
514                  break;
515    #ifndef PASO_MPI
516               case(Elements):
517               case(ReducedElements):
518                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
519                  break;
520               case(FaceElements):
521               case(ReducedFaceElements):
522                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
523                break;                break;
524               case(Points):
525                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
526                  break;
527               case(ContactElementsZero):
528               case(ContactElementsOne):
529               case(ReducedContactElementsZero):
530               case(ReducedContactElementsOne):
531                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
532                 break;
533    #else
534               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
535             case(Elements):             case(Elements):
536                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));             case(ReducedElements):
537                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
538                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);
539                break;                break;
540             case(FaceElements):             case(FaceElements):
541                Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedFaceElements):
542                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
543                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);
544                break;                break;
545             case(Points):             case(Points):
546                Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));                escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
547                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);
548                break;                break;
549             case(ContactElementsZero):             case(ContactElementsZero):
550             case(ContactElementsOne):             case(ContactElementsOne):
551                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));             case(ReducedContactElementsZero):
552               case(ReducedContactElementsOne):
553                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
554                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);
555               break;               break;
556    #endif
557             default:             default:
558               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
559               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();
560                 throw FinleyAdapterException(temp.str());
561               break;               break;
562          }          }
563          break;          break;
564       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
565         switch(target.getFunctionSpace().getTypeCode()) {         switch(target.getFunctionSpace().getTypeCode()) {
           case(ReducedDegreesOfFreedom):  
566            case(Nodes):            case(Nodes):
567               Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));               throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
568                 break;
569              case(ReducedNodes):
570                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
571                 break;
572              case(DegreesOfFreedom):
573                 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
574                 break;
575              case(ReducedDegreesOfFreedom):
576                 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
577               break;               break;
578            case(Elements):            case(Elements):
579               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));            case(ReducedElements):
580                 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
581               break;               break;
582            case(FaceElements):            case(FaceElements):
583               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));            case(ReducedFaceElements):
584                 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
585               break;               break;
586            case(Points):            case(Points):
587               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
588               break;               break;
589            case(ContactElementsZero):            case(ContactElementsZero):
590            case(ContactElementsOne):            case(ContactElementsOne):
591               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));            case(ReducedContactElementsZero):
592               break;            case(ReducedContactElementsOne):
593            case(DegreesOfFreedom):               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
              Finley_ErrorCode=TYPE_ERROR;  
              sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
594               break;               break;
595            default:            default:
596               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
597               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();
598                 throw FinleyAdapterException(temp.str());
599               break;               break;
600         }         }
601         break;         break;
602       default:       default:
603          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
604          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();
605            throw FinleyAdapterException(temp.str());
606          break;          break;
607    }    }
608    checkFinleyError();    checkFinleyError();
# Line 521  void MeshAdapter::interpolateOnDomain(Da Line 611  void MeshAdapter::interpolateOnDomain(Da
611  //  //
612  // copies the locations of sample points into x:  // copies the locations of sample points into x:
613  //  //
614  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
615  {  {
616    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
617    if (argDomain!=*this)    if (argDomain!=*this)
# Line 529  void MeshAdapter::setToX(Data& arg) cons Line 619  void MeshAdapter::setToX(Data& arg) cons
619    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
620    // in case of values node coordinates we can do the job directly:    // in case of values node coordinates we can do the job directly:
621    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
622       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       escriptDataC _arg=arg.getDataC();
623         Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
624    } else {    } else {
625       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
626       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       escriptDataC _tmp_data=tmp_data.getDataC();
627         Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
628       // this is then interpolated onto arg:       // this is then interpolated onto arg:
629       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
630    }    }
631    checkFinleyError();    checkFinleyError();
632  }  }
633    
634  //  //
635  // 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:
636  //  //
637  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
638  {  {
639    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
640    if (normalDomain!=*this)    if (normalDomain!=*this)
641       throw FinleyAdapterException("Error - Illegal domain of normal locations");       throw FinleyAdapterException("Error - Illegal domain of normal locations");
642    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
643      escriptDataC _normal=normal.getDataC();
644    switch(normal.getFunctionSpace().getTypeCode()) {    switch(normal.getFunctionSpace().getTypeCode()) {
645      case(Nodes):      case(Nodes):
646        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
647        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");        break;
648        case(ReducedNodes):
649          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
650        break;        break;
651      case(Elements):      case(Elements):
652        Finley_ErrorCode=TYPE_ERROR;        throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
653        sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");        break;
654        case(ReducedElements):
655          throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
656        break;        break;
657      case (FaceElements):      case (FaceElements):
658        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
659          break;
660        case (ReducedFaceElements):
661          Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
662        break;        break;
663      case(Points):      case(Points):
664        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");  
665        break;        break;
666      case (ContactElementsOne):      case (ContactElementsOne):
667      case (ContactElementsZero):      case (ContactElementsZero):
668        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));        Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
669          break;
670        case (ReducedContactElementsOne):
671        case (ReducedContactElementsZero):
672          Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
673        break;        break;
674      case(DegreesOfFreedom):      case(DegreesOfFreedom):
675        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.");  
676        break;        break;
677      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
678        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.");  
679        break;        break;
680      default:      default:
681        Finley_ErrorCode=TYPE_ERROR;        stringstream temp;
682        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();
683          throw FinleyAdapterException(temp.str());
684        break;        break;
685    }    }
686    checkFinleyError();    checkFinleyError();
687  }  }
688    
689  //  //
690  // interpolates data to other domain:  // interpolates data to other domain:
691  //  //
692  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
693  {  {
694    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
695    if (targetDomain!=*this)    if (targetDomain!=*this)
696       throw FinleyAdapterException("Error - Illegal domain of interpolation target");       throw FinleyAdapterException("Error - Illegal domain of interpolation target");
697    
698    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();  
699  }  }
700    
701  //  //
702  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
703  //  //
704  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
705  {  {
706    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
707    if (argDomain!=*this)    if (argDomain!=*this)
708       throw FinleyAdapterException("Error - Illegal domain of integration kernel");       throw FinleyAdapterException("Error - Illegal domain of integration kernel");
709    
710    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
711      escriptDataC _arg=arg.getDataC();
712    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
713       case(Nodes):       case(Nodes):
714          Finley_ErrorCode=TYPE_ERROR;          /* TODO */
715          sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
716            break;
717         case(ReducedNodes):
718            /* TODO */
719            throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");
720          break;          break;
721       case(Elements):       case(Elements):
722          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
723            break;
724         case(ReducedElements):
725            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
726          break;          break;
727       case(FaceElements):       case(FaceElements):
728          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
729            break;
730         case(ReducedFaceElements):
731            Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
732          break;          break;
733       case(Points):       case(Points):
734          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.");  
735          break;          break;
736       case(ContactElementsZero):       case(ContactElementsZero):
737          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
738            break;
739         case(ReducedContactElementsZero):
740            Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
741          break;          break;
742       case(ContactElementsOne):       case(ContactElementsOne):
743          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
744            break;
745         case(ReducedContactElementsOne):
746            Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
747          break;          break;
748       case(DegreesOfFreedom):       case(DegreesOfFreedom):
749          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.");  
750          break;          break;
751       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
752          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.");  
753          break;          break;
754       default:       default:
755          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
756          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();
757            throw FinleyAdapterException(temp.str());
758          break;          break;
759    }    }
760    checkFinleyError();    checkFinleyError();
761  }  }
762    
763  //  //
764  // calculates the gradient of arg:  // calculates the gradient of arg:
765  //  //
766  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
767  {  {
768    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
769    if (argDomain!=*this)    if (argDomain!=*this)
# Line 654  void MeshAdapter::setToGradient(Data& gr Line 773  void MeshAdapter::setToGradient(Data& gr
773       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
774    
775    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
776      escriptDataC _grad=grad.getDataC();
777      escriptDataC nodeDataC;
778    #ifdef PASO_MPI
779      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
780      if( arg.getFunctionSpace().getTypeCode() != Nodes )
781      {
782        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
783        nodeDataC = nodeTemp.getDataC();
784      }
785      else
786        nodeDataC = arg.getDataC();
787    #else
788      nodeDataC = arg.getDataC();
789    #endif
790    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
791         case(Nodes):         case(Nodes):
792            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
793            sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");            break;
794           case(ReducedNodes):
795              throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
796            break;            break;
797         case(Elements):         case(Elements):
798            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
799              break;
800           case(ReducedElements):
801              Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
802            break;            break;
803         case(FaceElements):         case(FaceElements):
804            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
805              break;
806           case(ReducedFaceElements):
807              Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
808            break;            break;
809         case(Points):         case(Points):
810            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Gradient at points is not supported.");
           sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");  
811            break;            break;
812         case(ContactElementsZero):         case(ContactElementsZero):
813            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
814              break;
815           case(ReducedContactElementsZero):
816              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
817            break;            break;
818         case(ContactElementsOne):         case(ContactElementsOne):
819            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
820              break;
821           case(ReducedContactElementsOne):
822              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
823            break;            break;
824         case(DegreesOfFreedom):         case(DegreesOfFreedom):
825            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.");  
826            break;            break;
827         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
828            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.");  
829            break;            break;
830         default:         default:
831            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
832            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();
833              throw FinleyAdapterException(temp.str());
834            break;            break;
835    }    }
836    checkFinleyError();    checkFinleyError();
837  }  }
838    
839  //  //
840  // returns the size of elements:  // returns the size of elements:
841  //  //
842  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
843  {  {
844    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
845    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
846    switch(size.getFunctionSpace().getTypeCode()) {    switch(size.getFunctionSpace().getTypeCode()) {
847         case(Nodes):         case(Nodes):
848            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of nodes is not supported.");
849            sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");            break;
850           case(ReducedNodes):
851              throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
852            break;            break;
853         case(Elements):         case(Elements):
854            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
855            break;            break;
856           case(ReducedElements):
857              Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
858              break;
859         case(FaceElements):         case(FaceElements):
860            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
861            break;            break;
862           case(ReducedFaceElements):
863              Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
864              break;
865         case(Points):         case(Points):
866            Finley_ErrorCode=TYPE_ERROR;            throw FinleyAdapterException("Error - Size of point elements is not supported.");
           sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");  
867            break;            break;
868         case(ContactElementsZero):         case(ContactElementsZero):
869         case(ContactElementsOne):         case(ContactElementsOne):
870            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);            Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
871            break;            break;
872           case(ReducedContactElementsZero):
873           case(ReducedContactElementsOne):
874              Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
875              break;
876         case(DegreesOfFreedom):         case(DegreesOfFreedom):
877            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.");  
878            break;            break;
879         case(ReducedDegreesOfFreedom):         case(ReducedDegreesOfFreedom):
880            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.");  
881            break;            break;
882         default:         default:
883            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
884            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();
885              throw FinleyAdapterException(temp.str());
886            break;            break;
887    }    }
888    checkFinleyError();    checkFinleyError();
889  }  }
890    
891  // sets the location of nodes:  // sets the location of nodes:
892  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
893  {  {
894    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
895      escriptDataC tmp;
896    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
897    if (newDomain!=*this)    if (newDomain!=*this)
898       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
899    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    tmp = new_x.getDataC();
900      Finley_Mesh_setCoordinates(mesh,&tmp);
901    checkFinleyError();    checkFinleyError();
902  }  }
903    
904  // saves a data array in openDX format:  // saves a data array in openDX format:
905  void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
906  {  {
907    Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));      int MAX_namelength=256;
908    checkFinleyError();      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
909      /* win32 refactor */
910      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
911      for(int i=0;i<num_data;i++)
912      {
913        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
914      }
915    
916      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
917      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
918      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
919    
920        boost::python::list keys=arg.keys();
921        for (int i=0;i<num_data;++i) {
922             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
923             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
924                 throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
925             data[i]=d.getDataC();
926             ptr_data[i]=&(data[i]);
927             std::string n=boost::python::extract<std::string>(keys[i]);
928             c_names[i]=&(names[i][0]);
929             if (n.length()>MAX_namelength-1) {
930                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
931                c_names[i][MAX_namelength-1]='\0';
932             } else {
933                strcpy(c_names[i],n.c_str());
934             }
935        }
936        Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
937        checkFinleyError();
938        
939          /* win32 refactor */
940      TMPMEMFREE(c_names);
941      TMPMEMFREE(data);
942      TMPMEMFREE(ptr_data);
943      for(int i=0;i<num_data;i++)
944      {
945        TMPMEMFREE(names[i]);
946      }
947      TMPMEMFREE(names);
948    
949        return;
950  }  }
951    
952    // saves a data array in openVTK format:
953    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
954    {
955        int MAX_namelength=256;
956        const int num_data=boost::python::extract<int>(arg.attr("__len__")());
957      /* win32 refactor */
958      char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
959      for(int i=0;i<num_data;i++)
960      {
961        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
962      }
963    
964      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
965      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
966      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
967    
968        boost::python::list keys=arg.keys();
969        for (int i=0;i<num_data;++i) {
970             escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
971             if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
972                 throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
973             data[i]=d.getDataC();
974             ptr_data[i]=&(data[i]);
975             std::string n=boost::python::extract<std::string>(keys[i]);
976             c_names[i]=&(names[i][0]);
977             if (n.length()>MAX_namelength-1) {
978                strncpy(c_names[i],n.c_str(),MAX_namelength-1);
979                c_names[i][MAX_namelength-1]='\0';
980             } else {
981                strcpy(c_names[i],n.c_str());
982             }
983        }
984    #ifndef PASO_MPI    
985        Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
986    #else
987        Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
988    #endif
989    
990    checkFinleyError();
991      /* win32 refactor */
992      TMPMEMFREE(c_names);
993      TMPMEMFREE(data);
994      TMPMEMFREE(ptr_data);
995      for(int i=0;i<num_data;i++)
996      {
997        TMPMEMFREE(names[i]);
998      }
999      TMPMEMFREE(names);
1000    
1001        return;
1002    }
1003                                                                                                                                                                            
1004                                                                                                                                                                            
1005  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1006  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1007                        const int row_blocksize,                        const int row_blocksize,
# Line 781  SystemMatrixAdapter MeshAdapter::newSyst Line 1036  SystemMatrixAdapter MeshAdapter::newSyst
1036      }      }
1037      // generate matrix:      // generate matrix:
1038            
1039      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);  
1040      checkFinleyError();      checkFinleyError();
1041        Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1042        checkPasoError();
1043        Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
1044      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1045  }  }
1046    
1047  //  //
1048  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
1049  // TODO:  // TODO:
1050  //  //
1051    
1052  //  //
1053  // 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:
1054  bool MeshAdapter::isCellOriented(int functionSpaceCode) const  bool MeshAdapter::isCellOriented(int functionSpaceCode) const
# Line 806  bool MeshAdapter::isCellOriented(int fun Line 1064  bool MeshAdapter::isCellOriented(int fun
1064         case(Points):         case(Points):
1065         case(ContactElementsZero):         case(ContactElementsZero):
1066         case(ContactElementsOne):         case(ContactElementsOne):
1067           case(ReducedElements):
1068           case(ReducedFaceElements):
1069           case(ReducedContactElementsZero):
1070           case(ReducedContactElementsOne):
1071            return true;            return true;
1072            break;            break;
1073         default:         default:
1074            Finley_ErrorCode=TYPE_ERROR;            stringstream temp;
1075            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;
1076              throw FinleyAdapterException(temp.str());
1077            break;            break;
1078    }    }
1079    checkFinleyError();    checkFinleyError();
1080    return false;    return false;
1081  }  }
1082    
1083  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1084  {  {
1085    switch(functionSpaceType_source) {    switch(functionSpaceType_source) {
1086       case(Nodes):       case(Nodes):
1087          switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1088             case(Nodes):             case(Nodes):
1089               case(ReducedNodes):
1090             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
1091             case(DegreesOfFreedom):             case(DegreesOfFreedom):
1092             case(Elements):             case(Elements):
1093               case(ReducedElements):
1094               case(FaceElements):
1095               case(ReducedFaceElements):
1096               case(Points):
1097               case(ContactElementsZero):
1098               case(ReducedContactElementsZero):
1099               case(ContactElementsOne):
1100               case(ReducedContactElementsOne):
1101                   return true;
1102               default:
1103                   stringstream temp;
1104                   temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1105                   throw FinleyAdapterException(temp.str());
1106            }
1107            break;
1108         case(ReducedNodes):
1109            switch(functionSpaceType_target) {
1110               case(ReducedNodes):
1111               case(ReducedDegreesOfFreedom):
1112               case(Elements):
1113               case(ReducedElements):
1114             case(FaceElements):             case(FaceElements):
1115               case(ReducedFaceElements):
1116             case(Points):             case(Points):
1117             case(ContactElementsZero):             case(ContactElementsZero):
1118               case(ReducedContactElementsZero):
1119             case(ContactElementsOne):             case(ContactElementsOne):
1120               case(ReducedContactElementsOne):
1121                 return true;                 return true;
1122              case(Nodes):
1123              case(DegreesOfFreedom):
1124                 return false;
1125             default:             default:
1126                 Finley_ErrorCode=TYPE_ERROR;                 stringstream temp;
1127                 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;
1128                   throw FinleyAdapterException(temp.str());
1129          }          }
1130          break;          break;
1131       case(Elements):       case(Elements):
# Line 841  bool MeshAdapter::probeInterpolationOnDo Line 1134  bool MeshAdapter::probeInterpolationOnDo
1134          } else {          } else {
1135             return false;             return false;
1136          }          }
1137         case(ReducedElements):
1138            if (functionSpaceType_target==ReducedElements) {
1139               return true;
1140            } else {
1141               return false;
1142            }
1143       case(FaceElements):       case(FaceElements):
1144          if (functionSpaceType_target==FaceElements) {          if (functionSpaceType_target==FaceElements) {
1145             return true;             return true;
1146          } else {          } else {
1147             return false;             return false;
1148          }          }
1149         case(ReducedFaceElements):
1150            if (functionSpaceType_target==ReducedFaceElements) {
1151               return true;
1152            } else {
1153               return false;
1154            }
1155       case(Points):       case(Points):
1156          if (functionSpaceType_target==Points) {          if (functionSpaceType_target==Points) {
1157             return true;             return true;
# Line 860  bool MeshAdapter::probeInterpolationOnDo Line 1165  bool MeshAdapter::probeInterpolationOnDo
1165          } else {          } else {
1166             return false;             return false;
1167          }          }
1168         case(ReducedContactElementsZero):
1169         case(ReducedContactElementsOne):
1170            if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1171               return true;
1172            } else {
1173               return false;
1174            }
1175       case(DegreesOfFreedom):       case(DegreesOfFreedom):
1176          switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1177             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
1178             case(DegreesOfFreedom):             case(DegreesOfFreedom):
1179             case(Nodes):             case(Nodes):
1180               case(ReducedNodes):
1181             case(Elements):             case(Elements):
1182             case(FaceElements):             case(ReducedElements):
1183             case(Points):             case(Points):
1184               case(FaceElements):
1185               case(ReducedFaceElements):
1186             case(ContactElementsZero):             case(ContactElementsZero):
1187               case(ReducedContactElementsZero):
1188             case(ContactElementsOne):             case(ContactElementsOne):
1189               case(ReducedContactElementsOne):
1190                return true;                return true;
1191             default:             default:
1192               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
1193               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;
1194                 throw FinleyAdapterException(temp.str());
1195          }          }
1196          break;          break;
1197       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
1198         switch(functionSpaceType_target) {         switch(functionSpaceType_target) {
1199            case(ReducedDegreesOfFreedom):            case(ReducedDegreesOfFreedom):
1200            case(Nodes):            case(ReducedNodes):
1201            case(Elements):            case(Elements):
1202              case(ReducedElements):
1203            case(FaceElements):            case(FaceElements):
1204              case(ReducedFaceElements):
1205            case(Points):            case(Points):
1206            case(ContactElementsZero):            case(ContactElementsZero):
1207              case(ReducedContactElementsZero):
1208            case(ContactElementsOne):            case(ContactElementsOne):
1209              case(ReducedContactElementsOne):
1210                return true;                return true;
1211              case(Nodes):
1212            case(DegreesOfFreedom):            case(DegreesOfFreedom):
1213               return false;               return false;
1214            default:            default:
1215               Finley_ErrorCode=TYPE_ERROR;               stringstream temp;
1216               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;
1217                 throw FinleyAdapterException(temp.str());
1218         }         }
1219         break;         break;
1220       default:       default:
1221          Finley_ErrorCode=TYPE_ERROR;          stringstream temp;
1222          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;
1223            throw FinleyAdapterException(temp.str());
1224          break;          break;
1225    }    }
1226    checkFinleyError();    checkFinleyError();
1227    return false;    return false;
1228  }  }
1229    
1230  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
1231  {  {
1232      return false;      return false;
1233  }  }
1234  bool MeshAdapter::operator==(const MeshAdapter& other) const  
1235    bool MeshAdapter::operator==(const AbstractDomain& other) const
1236  {  {
1237    return (m_finleyMesh==other.m_finleyMesh);    const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1238      if (temp!=0) {
1239        return (m_finleyMesh==temp->m_finleyMesh);
1240      } else {
1241        return false;
1242      }
1243  }  }
1244    
1245  int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const  bool MeshAdapter::operator!=(const AbstractDomain& other) const
1246  {  {
1247     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);    return !(operator==(other));
1248     checkFinleyError();  }
1249    
1250    int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1251    {
1252       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1253       checkPasoError();
1254     return out;     return out;
1255  }  }
1256  Data MeshAdapter::getX() const  
1257    escript::Data MeshAdapter::getX() const
1258  {  {
1259    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
1260  }  }
1261  Data MeshAdapter::getNormal() const  
1262    escript::Data MeshAdapter::getNormal() const
1263  {  {
1264    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1265  }  }
1266  Data MeshAdapter::getSize() const  
1267    escript::Data MeshAdapter::getSize() const
1268  {  {
1269    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
1270  }  }
1271    
1272    int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1273    {
1274      int *out=0,i;
1275      Finley_Mesh* mesh=m_finleyMesh.get();
1276      switch (functionSpaceType) {
1277        case(Nodes):
1278          if (mesh->Nodes!=NULL) {
1279            out=mesh->Nodes->Id;
1280            break;
1281          }
1282        case(ReducedNodes):
1283          throw FinleyAdapterException("Error -  ReducedNodes not supported yet.");
1284          break;
1285        case(Elements):
1286          out=mesh->FaceElements->Id;
1287          break;
1288        case(ReducedElements):
1289          out=mesh->Elements->Id;
1290          break;
1291        case(FaceElements):
1292          out=mesh->FaceElements->Id;
1293          break;
1294        case(ReducedFaceElements):
1295          out=mesh->FaceElements->Id;
1296          break;
1297        case(Points):
1298          out=mesh->Points->Id;
1299          break;
1300        case(ContactElementsZero):
1301          out=mesh->ContactElements->Id;
1302          break;
1303        case(ReducedContactElementsZero):
1304          out=mesh->ContactElements->Id;
1305          break;
1306        case(ContactElementsOne):
1307          out=mesh->ContactElements->Id;
1308          break;
1309        case(ReducedContactElementsOne):
1310          out=mesh->ContactElements->Id;
1311          break;
1312        case(DegreesOfFreedom):
1313          out=mesh->Nodes->degreeOfFreedomId;
1314          break;
1315        case(ReducedDegreesOfFreedom):
1316          out=mesh->Nodes->reducedDegreeOfFreedomId;
1317          break;
1318        default:
1319          stringstream temp;
1320          temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1321          throw FinleyAdapterException(temp.str());
1322          break;
1323      }
1324      return out;
1325    }
1326    int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1327    {
1328      int out=0;
1329      Finley_Mesh* mesh=m_finleyMesh.get();
1330      switch (functionSpaceType) {
1331      case(Nodes):
1332        out=mesh->Nodes->Tag[sampleNo];
1333        break;
1334      case(ReducedNodes):
1335        throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1336        break;
1337      case(Elements):
1338        out=mesh->Elements->Tag[sampleNo];
1339        break;
1340      case(ReducedElements):
1341        out=mesh->Elements->Tag[sampleNo];
1342        break;
1343      case(FaceElements):
1344        out=mesh->FaceElements->Tag[sampleNo];
1345        break;
1346      case(ReducedFaceElements):
1347        out=mesh->FaceElements->Tag[sampleNo];
1348        break;
1349      case(Points):
1350        out=mesh->Points->Tag[sampleNo];
1351        break;
1352      case(ContactElementsZero):
1353        out=mesh->ContactElements->Tag[sampleNo];
1354        break;
1355      case(ReducedContactElementsZero):
1356        out=mesh->ContactElements->Tag[sampleNo];
1357        break;
1358      case(ContactElementsOne):
1359        out=mesh->ContactElements->Tag[sampleNo];
1360        break;
1361      case(ReducedContactElementsOne):
1362        out=mesh->ContactElements->Tag[sampleNo];
1363        break;
1364      case(DegreesOfFreedom):
1365        throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1366        break;
1367      case(ReducedDegreesOfFreedom):
1368        throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1369        break;
1370      default:
1371        stringstream temp;
1372        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1373        throw FinleyAdapterException(temp.str());
1374        break;
1375      }
1376      return out;
1377    }
1378    
1379    
1380  bool MeshAdapter::operator!=(const MeshAdapter& other) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1381  {  {
1382    return !operator==(other);    Finley_Mesh* mesh=m_finleyMesh.get();
1383      escriptDataC tmp=mask.getDataC();
1384      switch(functionSpaceType) {
1385           case(Nodes):
1386              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1387              break;
1388           case(ReducedNodes):
1389              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1390              break;
1391           case(DegreesOfFreedom):
1392              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1393              break;
1394           case(ReducedDegreesOfFreedom):
1395              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1396              break;
1397           case(Elements):
1398              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1399              break;
1400           case(ReducedElements):
1401              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1402              break;
1403           case(FaceElements):
1404              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1405              break;
1406           case(ReducedFaceElements):
1407              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1408              break;
1409           case(Points):
1410              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1411              break;
1412           case(ContactElementsZero):
1413              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1414              break;
1415           case(ReducedContactElementsZero):
1416              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1417              break;
1418           case(ContactElementsOne):
1419              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1420              break;
1421           case(ReducedContactElementsOne):
1422              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1423              break;
1424           default:
1425              stringstream temp;
1426              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1427              throw FinleyAdapterException(temp.str());
1428      }
1429      checkFinleyError();
1430      return;
1431    }
1432    
1433    void MeshAdapter::setTagMap(const std::string& name,  int tag)
1434    {
1435      Finley_Mesh* mesh=m_finleyMesh.get();
1436      Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1437      checkPasoError();
1438      // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1439    }
1440    
1441    int MeshAdapter::getTag(const std::string& name) const
1442    {
1443      Finley_Mesh* mesh=m_finleyMesh.get();
1444      int tag=0;
1445      tag=Finley_Mesh_getTag(mesh, name.c_str());
1446      checkPasoError();
1447      // throwStandardException("MeshAdapter::getTag is not implemented.");
1448      return tag;
1449    }
1450    
1451    bool MeshAdapter::isValidTagName(const std::string& name) const
1452    {
1453      Finley_Mesh* mesh=m_finleyMesh.get();
1454      return Finley_Mesh_isValidTagName(mesh,name.c_str());
1455  }  }
1456  // bool MeshAdapter::operator==(const AbstractDomain& other) const  
1457  // {  std::string MeshAdapter::showTagNames() const
1458    // try {  {
1459      // const MeshAdapter& ref = dynamic_cast<const MeshAdapter&>(other);    stringstream temp;
1460      // return (operator==(ref));    Finley_Mesh* mesh=m_finleyMesh.get();
1461    // }    Finley_TagMap* tag_map=mesh->TagMap;
1462    // catch (bad_cast) {    while (tag_map) {
1463      // return false;       temp << tag_map->name;
1464    // }       tag_map=tag_map->next;
1465  // }       if (tag_map) temp << ", ";
1466      }
1467      return temp.str();
1468    }
1469    
1470  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26