/[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

revision 532 by gross, Wed Feb 15 09:45:53 2006 UTC revision 1059 by gross, Fri Mar 23 11:01:53 2007 UTC
# Line 13  Line 13 
13   ******************************************************************************   ******************************************************************************
14  */  */
15    
16    #ifdef PASO_MPI
17    #include <mpi.h>
18    #endif
19  #include "MeshAdapter.h"  #include "MeshAdapter.h"
20    
21  #include "Data.h"  #include "escript/Data.h"
22  #include "DataFactory.h"  #include "escript/DataFactory.h"
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
# Line 30  const int MeshAdapter::DegreesOfFreedom= Line 33  const int MeshAdapter::DegreesOfFreedom=
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::Elements=FINLEY_ELEMENTS;  const int MeshAdapter::Elements=FINLEY_ELEMENTS;
36    const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
37  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;  const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
38    const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
39  const int MeshAdapter::Points=FINLEY_POINTS;  const int MeshAdapter::Points=FINLEY_POINTS;
40  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;  const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
41    const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
42  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
43    const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
44    
45  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
46  {  {
# Line 70  Finley_Mesh* MeshAdapter::getFinley_Mesh Line 77  Finley_Mesh* MeshAdapter::getFinley_Mesh
77    
78  void MeshAdapter::write(const std::string& fileName) const  void MeshAdapter::write(const std::string& fileName) const
79  {  {
80    char fName[fileName.size()+1];    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
81    strcpy(fName,fileName.c_str());    strcpy(fName,fileName.c_str());
82    Finley_Mesh_write(m_finleyMesh.get(),fName);    Finley_Mesh_write(m_finleyMesh.get(),fName);
83    checkFinleyError();    checkFinleyError();
84      TMPMEMFREE(fName);
85  }  }
86    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
87  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
88  {  {
89    return "FinleyMesh";    return "FinleyMesh";
# Line 117  void MeshAdapter::setFunctionSpaceTypeNa Line 118  void MeshAdapter::setFunctionSpaceTypeNa
118    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
119      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));      (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
120    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
121        (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
122      m_functionSpaceTypeNames.insert
123      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));      (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
124    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
125        (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
126      m_functionSpaceTypeNames.insert
127      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));      (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
128    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
129      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));      (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
130    m_functionSpaceTypeNames.insert    m_functionSpaceTypeNames.insert
131        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
132      m_functionSpaceTypeNames.insert
133      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));      (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
134      m_functionSpaceTypeNames.insert
135        (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
136  }  }
137    
138  int MeshAdapter::getContinuousFunctionCode() const  int MeshAdapter::getContinuousFunctionCode() const
# Line 135  int MeshAdapter::getFunctionCode() const Line 144  int MeshAdapter::getFunctionCode() const
144  {  {
145    return Elements;    return Elements;
146  }  }
147    int MeshAdapter::getReducedFunctionCode() const
148    {
149      return ReducedElements;
150    }
151    
152  int MeshAdapter::getFunctionOnBoundaryCode() const  int MeshAdapter::getFunctionOnBoundaryCode() const
153  {  {
154    return FaceElements;    return FaceElements;
155  }  }
156    int MeshAdapter::getReducedFunctionOnBoundaryCode() const
157    {
158      return ReducedFaceElements;
159    }
160    
161  int MeshAdapter::getFunctionOnContactZeroCode() const  int MeshAdapter::getFunctionOnContactZeroCode() const
162  {  {
163    return ContactElementsZero;    return ContactElementsZero;
164  }  }
165    int MeshAdapter::getReducedFunctionOnContactZeroCode() const
166    {
167      return ReducedContactElementsZero;
168    }
169    
170  int MeshAdapter::getFunctionOnContactOneCode() const  int MeshAdapter::getFunctionOnContactOneCode() const
171  {  {
172    return ContactElementsOne;    return ContactElementsOne;
173  }  }
174    int MeshAdapter::getReducedFunctionOnContactOneCode() const
175    {
176      return ReducedContactElementsOne;
177    }
178    
179  int MeshAdapter::getSolutionCode() const  int MeshAdapter::getSolutionCode() const
180  {  {
# Line 167  int MeshAdapter::getDiracDeltaFunctionCo Line 192  int MeshAdapter::getDiracDeltaFunctionCo
192  }  }
193    
194  //  //
 // 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;  
 }  
   
 //  
195  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
196  //  //
197  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 367  pair<int,int> MeshAdapter::getDataShape( Line 248  pair<int,int> MeshAdapter::getDataShape(
248        case(DegreesOfFreedom):        case(DegreesOfFreedom):
249             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
250               numDataPointsPerSample=1;               numDataPointsPerSample=1;
251    #ifndef PASO_MPI
252               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
253    #else
254                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
255    #endif
256             }             }
257             break;             break;
258        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
259             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
260               numDataPointsPerSample=1;               numDataPointsPerSample=1;
261    #ifndef PASO_MPI
262               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
263    #else
264                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
265    #endif
266             }             }
267             break;             break;
268        default:        default:
# Line 394  void MeshAdapter::addPDEToSystem( Line 283  void MeshAdapter::addPDEToSystem(
283                       const escript::Data& d, const escript::Data& y,                       const escript::Data& d, const escript::Data& y,
284                       const escript::Data& d_contact,const escript::Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
285  {  {
286       escriptDataC _rhs=rhs.getDataC();
287       escriptDataC _A  =A.getDataC();
288       escriptDataC _B=B.getDataC();
289       escriptDataC _C=C.getDataC();
290       escriptDataC _D=D.getDataC();
291       escriptDataC _X=X.getDataC();
292       escriptDataC _Y=Y.getDataC();
293       escriptDataC _d=d.getDataC();
294       escriptDataC _y=y.getDataC();
295       escriptDataC _d_contact=d_contact.getDataC();
296       escriptDataC _y_contact=y_contact.getDataC();
297    
298     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
299     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),  
300                         &(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 );
301     checkFinleyError();     checkFinleyError();
302     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
303                    mat.getPaso_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);  
304     checkFinleyError();     checkFinleyError();
305     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  
306                    mat.getPaso_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);  
307     checkFinleyError();     checkFinleyError();
308  }  }
309    
310  //  //
311  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
312  //  //
313  void MeshAdapter::addPDEToRHS( escript::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  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  
314  {  {
315     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
316    
317     // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));     escriptDataC _rhs=rhs.getDataC();
318     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));     escriptDataC _X=X.getDataC();
319     checkFinleyError();     escriptDataC _Y=Y.getDataC();
320       escriptDataC _y=y.getDataC();
321       escriptDataC _y_contact=y_contact.getDataC();
322    
323     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
324     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);     checkFinleyError();
325    
326       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
327     checkFinleyError();     checkFinleyError();
328     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);  
329     // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
330     checkFinleyError();     checkFinleyError();
331  }  }
332    
# Line 441  void MeshAdapter::interpolateOnDomain(es Line 337  void MeshAdapter::interpolateOnDomain(es
337  {  {
338    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
339    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
340    if (inDomain!=*this)    if (inDomain!=*this)  
341       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
342    if (targetDomain!=*this)    if (targetDomain!=*this)
343       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
344    
345    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
346    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 508  void MeshAdapter::interpolateOnDomain(es Line 404  void MeshAdapter::interpolateOnDomain(es
404             break;             break;
405          }          }
406          break;          break;
407       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
408          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
409             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
410             case(DegreesOfFreedom):             case(DegreesOfFreedom):
411             case(Nodes):             case(Nodes):
412                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
413                break;                break;
414    #ifndef PASO_MPI
415             case(Elements):             case(Elements):
416                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
417                break;                break;
# Line 528  void MeshAdapter::interpolateOnDomain(es Line 425  void MeshAdapter::interpolateOnDomain(es
425             case(ContactElementsOne):             case(ContactElementsOne):
426                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
427               break;               break;
428    #else
429               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
430               case(Elements):
431               {
432                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
433                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
434                  break;
435               }
436               case(FaceElements):
437               {
438                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
439                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
440                  break;
441               }
442               case(Points):
443               {
444                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
445                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
446                  break;
447               }
448               case(ContactElementsZero):
449               case(ContactElementsOne):
450               {
451                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
452                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
453                 break;
454               }
455    #endif
456             default:             default:
457               stringstream temp;               stringstream temp;
458               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
# Line 706  void MeshAdapter::setToGradient(escript: Line 631  void MeshAdapter::setToGradient(escript:
631       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
632    
633    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
634      escriptDataC nodeDataC;
635    #ifdef PASO_MPI
636      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
637      if( arg.getFunctionSpace().getTypeCode() != Nodes )
638      {
639        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
640        nodeDataC = nodeTemp.getDataC();
641      }
642      else
643        nodeDataC = arg.getDataC();
644    #else
645      nodeDataC = arg.getDataC();
646    #endif
647    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
648         case(Nodes):         case(Nodes):
649            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
650            break;            break;
651         case(Elements):         case(Elements):
652            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
653            break;            break;
654         case(FaceElements):         case(FaceElements):
655            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
656            break;            break;
657         case(Points):         case(Points):
658            throw FinleyAdapterException("Error - Gradient at points is not supported.");            throw FinleyAdapterException("Error - Gradient at points is not supported.");
659            break;            break;
660         case(ContactElementsZero):         case(ContactElementsZero):
661            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
662            break;            break;
663         case(ContactElementsOne):         case(ContactElementsOne):
664            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
665            break;            break;
666         case(DegreesOfFreedom):         case(DegreesOfFreedom):
667            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
# Line 783  void MeshAdapter::setToSize(escript::Dat Line 721  void MeshAdapter::setToSize(escript::Dat
721  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
722  {  {
723    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
724      escriptDataC tmp;
725    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
726    if (newDomain!=*this)    if (newDomain!=*this)
727       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
728    Finley_Mesh_setCoordinates(mesh,&(new_x.getDataC()));    tmp = new_x.getDataC();
729      Finley_Mesh_setCoordinates(mesh,&tmp);
730    checkFinleyError();    checkFinleyError();
731  }  }
732    
# Line 795  void MeshAdapter::saveDX(const std::stri Line 735  void MeshAdapter::saveDX(const std::stri
735  {  {
736      int MAX_namelength=256;      int MAX_namelength=256;
737      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
738      char names[num_data][MAX_namelength];    /* win32 refactor */
739      char* c_names[num_data];    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
740      escriptDataC data[num_data];    for(int i=0;i<num_data;i++)
741      escriptDataC* ptr_data[num_data];    {
742        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
743      }
744    
745      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
746      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
747      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
748    
749      boost::python::list keys=arg.keys();      boost::python::list keys=arg.keys();
750      for (int i=0;i<num_data;++i) {      for (int i=0;i<num_data;++i) {
751           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
752           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
753               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
754           data[i]=d.getDataC();           data[i]=d.getDataC();
755           ptr_data[i]=&(data[i]);           ptr_data[i]=&(data[i]);
756           std::string n=boost::python::extract<std::string>(keys[i]);           std::string n=boost::python::extract<std::string>(keys[i]);
# Line 818  void MeshAdapter::saveDX(const std::stri Line 764  void MeshAdapter::saveDX(const std::stri
764      }      }
765      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
766      checkFinleyError();      checkFinleyError();
767        
768          /* win32 refactor */
769      TMPMEMFREE(c_names);
770      TMPMEMFREE(data);
771      TMPMEMFREE(ptr_data);
772      for(int i=0;i<num_data;i++)
773      {
774        TMPMEMFREE(names[i]);
775      }
776      TMPMEMFREE(names);
777    
778      return;      return;
779  }  }
780    
# Line 826  void MeshAdapter::saveVTK(const std::str Line 783  void MeshAdapter::saveVTK(const std::str
783  {  {
784      int MAX_namelength=256;      int MAX_namelength=256;
785      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
786      char names[num_data][MAX_namelength];    /* win32 refactor */
787      char* c_names[num_data];    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
788      escriptDataC data[num_data];    for(int i=0;i<num_data;i++)
789      escriptDataC* ptr_data[num_data];    {
790        names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
791      }
792    
793      char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
794      escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
795      escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
796    
797      boost::python::list keys=arg.keys();      boost::python::list keys=arg.keys();
798      for (int i=0;i<num_data;++i) {      for (int i=0;i<num_data;++i) {
# Line 847  void MeshAdapter::saveVTK(const std::str Line 810  void MeshAdapter::saveVTK(const std::str
810              strcpy(c_names[i],n.c_str());              strcpy(c_names[i],n.c_str());
811           }           }
812      }      }
813    #ifndef PASO_MPI    
814      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
815      checkFinleyError();  #else
816        Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
817    #endif
818    
819    checkFinleyError();
820      /* win32 refactor */
821      TMPMEMFREE(c_names);
822      TMPMEMFREE(data);
823      TMPMEMFREE(ptr_data);
824      for(int i=0;i<num_data;i++)
825      {
826        TMPMEMFREE(names[i]);
827      }
828      TMPMEMFREE(names);
829    
830      return;      return;
831  }  }
832                                                                                                                                                                                                                                                                                                                                                    
# Line 1061  escript::Data MeshAdapter::getSize() con Line 1039  escript::Data MeshAdapter::getSize() con
1039    
1040  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1041  {  {
1042    int* tagList;    int out=0;
1043    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
1044    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
1045    return tagList[sampleNo];    case(Nodes):
1046        out=mesh->Nodes->Tag[sampleNo];
1047        break;
1048      case(Elements):
1049        out=mesh->Elements->Tag[sampleNo];
1050        break;
1051      case(FaceElements):
1052        out=mesh->FaceElements->Tag[sampleNo];
1053        break;
1054      case(Points):
1055        out=mesh->Points->Tag[sampleNo];
1056        break;
1057      case(ContactElementsZero):
1058        out=mesh->ContactElements->Tag[sampleNo];
1059        break;
1060      case(ContactElementsOne):
1061        out=mesh->ContactElements->Tag[sampleNo];
1062        break;
1063      case(DegreesOfFreedom):
1064        break;
1065      case(ReducedDegreesOfFreedom):
1066        break;
1067      default:
1068        stringstream temp;
1069        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1070        throw FinleyAdapterException(temp.str());
1071        break;
1072      }
1073      return out;
1074    }
1075    int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1076    {
1077      int *out=0,i;
1078      Finley_Mesh* mesh=m_finleyMesh.get();
1079      switch (functionSpaceType) {
1080      case(Nodes):
1081        if (mesh->Nodes!=NULL) {
1082          out=mesh->Nodes->Id;
1083          break;
1084        }
1085      case(Elements):
1086        out=mesh->Elements->Id;
1087        break;
1088      case(FaceElements):
1089        out=mesh->FaceElements->Id;
1090        break;
1091      case(Points):
1092        out=mesh->Points->Id;
1093        break;
1094      case(ContactElementsZero):
1095        out=mesh->ContactElements->Id;
1096        break;
1097      case(ContactElementsOne):
1098        out=mesh->ContactElements->Id;
1099        break;
1100      case(DegreesOfFreedom):
1101        out=mesh->Nodes->degreeOfFreedomId;
1102        break;
1103      case(ReducedDegreesOfFreedom):
1104        out=mesh->Nodes->reducedDegreeOfFreedomId;
1105        break;
1106      default:
1107        stringstream temp;
1108        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1109        throw FinleyAdapterException(temp.str());
1110        break;
1111      }
1112      return out;
1113    }
1114    
1115    void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1116    {
1117      Finley_Mesh* mesh=m_finleyMesh.get();
1118      escriptDataC tmp=mask.getDataC();
1119      switch(functionSpaceType) {
1120           case(Nodes):
1121              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1122              break;
1123           case(DegreesOfFreedom):
1124              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1125              break;
1126           case(ReducedDegreesOfFreedom):
1127              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1128              break;
1129           case(Elements):
1130              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1131              break;
1132           case(FaceElements):
1133              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1134              break;
1135           case(Points):
1136              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1137              break;
1138           case(ContactElementsZero):
1139              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1140              break;
1141           case(ContactElementsOne):
1142              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1143              break;
1144           default:
1145              stringstream temp;
1146              temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1147              throw FinleyAdapterException(temp.str());
1148      }
1149      checkFinleyError();
1150      return;
1151    }
1152    
1153    void MeshAdapter::setTagMap(const std::string& name,  int tag)
1154    {
1155      Finley_Mesh* mesh=m_finleyMesh.get();
1156      Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1157      checkPasoError();
1158      // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1159    }
1160    
1161    int MeshAdapter::getTag(const std::string& name) const
1162    {
1163      Finley_Mesh* mesh=m_finleyMesh.get();
1164      int tag=0;
1165      tag=Finley_Mesh_getTag(mesh, name.c_str());
1166      checkPasoError();
1167      // throwStandardException("MeshAdapter::getTag is not implemented.");
1168      return tag;
1169  }  }
1170    
1171  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  bool MeshAdapter::isValidTagName(const std::string& name) const
1172  {  {
1173    int* referenceNoList;    Finley_Mesh* mesh=m_finleyMesh.get();
1174    int numReferenceNo;    return Finley_Mesh_isValidTagName(mesh,name.c_str());
1175    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);  }
1176    return referenceNoList[sampleNo];  
1177    std::string MeshAdapter::showTagNames() const
1178    {
1179      stringstream temp;
1180      Finley_Mesh* mesh=m_finleyMesh.get();
1181      Finley_TagMap* tag_map=mesh->TagMap;
1182      while (tag_map) {
1183         temp << tag_map->name;
1184         tag_map=tag_map->next;
1185         if (tag_map) temp << ", ";
1186      }
1187      return temp.str();
1188  }  }
1189    
1190  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.532  
changed lines
  Added in v.1059

  ViewVC Help
Powered by ViewVC 1.1.26