/[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 472 by jgs, Fri Jan 27 01:50:59 2006 UTC revision 682 by robwdcock, Mon Mar 27 02:43:09 2006 UTC
# Line 15  Line 15 
15    
16  #include "MeshAdapter.h"  #include "MeshAdapter.h"
17    
18    #include "escript/Data.h"
19    #include "escript/DataFactory.h"
20    
21  using namespace std;  using namespace std;
22  using namespace escript;  using namespace escript;
23    
# Line 73  void MeshAdapter::write(const std::strin Line 76  void MeshAdapter::write(const std::strin
76    checkFinleyError();    checkFinleyError();
77  }  }
78    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
79  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
80  {  {
81    return "FinleyMesh";    return "FinleyMesh";
# Line 164  int MeshAdapter::getDiracDeltaFunctionCo Line 160  int MeshAdapter::getDiracDeltaFunctionCo
160  }  }
161    
162  //  //
 // 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;  
 }  
   
 //  
163  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
164  //  //
165  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 386  pair<int,int> MeshAdapter::getDataShape( Line 238  pair<int,int> MeshAdapter::getDataShape(
238  // 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:
239  //  //
240  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
241                       SystemMatrixAdapter& mat, Data& rhs,                       SystemMatrixAdapter& mat, escript::Data& rhs,
242                       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,
243                       const Data& d, const Data& y,                       const escript::Data& d, const escript::Data& y,
244                       const Data& d_contact,const Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
245  {  {
246     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
247     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),
# Line 413  void MeshAdapter::addPDEToSystem( Line 265  void MeshAdapter::addPDEToSystem(
265  //  //
266  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
267  //  //
268  void MeshAdapter::addPDEToRHS( Data& rhs,  void MeshAdapter::addPDEToRHS( escript::Data& rhs,
269                       const  Data& X,const  Data& Y, const Data& y, const Data& y_contact) const                       const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
270  {  {
271     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
272    
# Line 434  void MeshAdapter::addPDEToRHS( Data& rhs Line 286  void MeshAdapter::addPDEToRHS( Data& rhs
286  //  //
287  // interpolates data between different function spaces:  // interpolates data between different function spaces:
288  //  //
289  void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
290  {  {
291    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
292    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
# Line 575  void MeshAdapter::interpolateOnDomain(Da Line 427  void MeshAdapter::interpolateOnDomain(Da
427  //  //
428  // copies the locations of sample points into x:  // copies the locations of sample points into x:
429  //  //
430  void MeshAdapter::setToX(Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
431  {  {
432    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
433    if (argDomain!=*this)    if (argDomain!=*this)
# Line 585  void MeshAdapter::setToX(Data& arg) cons Line 437  void MeshAdapter::setToX(Data& arg) cons
437    if (arg.getFunctionSpace().getTypeCode()==Nodes) {    if (arg.getFunctionSpace().getTypeCode()==Nodes) {
438       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
439    } else {    } else {
440       Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);       escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
441       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));       Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
442       // this is then interpolated onto arg:       // this is then interpolated onto arg:
443       interpolateOnDomain(arg,tmp_data);       interpolateOnDomain(arg,tmp_data);
# Line 596  void MeshAdapter::setToX(Data& arg) cons Line 448  void MeshAdapter::setToX(Data& arg) cons
448  //  //
449  // 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:
450  //  //
451  void MeshAdapter::setToNormal(Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
452  {  {
453    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());    const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
454    if (normalDomain!=*this)    if (normalDomain!=*this)
# Line 637  void MeshAdapter::setToNormal(Data& norm Line 489  void MeshAdapter::setToNormal(Data& norm
489  //  //
490  // interpolates data to other domain:  // interpolates data to other domain:
491  //  //
492  void MeshAdapter::interpolateACross(Data& target,const Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
493  {  {
494    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
495    if (targetDomain!=*this)    if (targetDomain!=*this)
# Line 649  void MeshAdapter::interpolateACross(Data Line 501  void MeshAdapter::interpolateACross(Data
501  //  //
502  // calculates the integral of a function defined of arg:  // calculates the integral of a function defined of arg:
503  //  //
504  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
505  {  {
506    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
507    if (argDomain!=*this)    if (argDomain!=*this)
# Line 693  void MeshAdapter::setToIntegrals(std::ve Line 545  void MeshAdapter::setToIntegrals(std::ve
545  //  //
546  // calculates the gradient of arg:  // calculates the gradient of arg:
547  //  //
548  void MeshAdapter::setToGradient(Data& grad,const Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
549  {  {
550    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());    const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
551    if (argDomain!=*this)    if (argDomain!=*this)
# Line 740  void MeshAdapter::setToGradient(Data& gr Line 592  void MeshAdapter::setToGradient(Data& gr
592  //  //
593  // returns the size of elements:  // returns the size of elements:
594  //  //
595  void MeshAdapter::setToSize(Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
596  {  {
597    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
598    escriptDataC tmp=size.getDataC();    escriptDataC tmp=size.getDataC();
# Line 777  void MeshAdapter::setToSize(Data& size) Line 629  void MeshAdapter::setToSize(Data& size)
629  }  }
630    
631  // sets the location of nodes:  // sets the location of nodes:
632  void MeshAdapter::setNewX(const Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
633  {  {
634    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
635    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
636    if (newDomain!=*this)    if (newDomain!=*this)
637       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
638    Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));    Finley_Mesh_setCoordinates(mesh,&(new_x.getDataC()));
639    checkFinleyError();    checkFinleyError();
640  }  }
641    
# Line 799  void MeshAdapter::saveDX(const std::stri Line 651  void MeshAdapter::saveDX(const std::stri
651    
652      boost::python::list keys=arg.keys();      boost::python::list keys=arg.keys();
653      for (int i=0;i<num_data;++i) {      for (int i=0;i<num_data;++i) {
654           Data& d=boost::python::extract<Data&>(arg[keys[i]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
655           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
656               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
657           data[i]=d.getDataC();           data[i]=d.getDataC();
# Line 830  void MeshAdapter::saveVTK(const std::str Line 682  void MeshAdapter::saveVTK(const std::str
682    
683      boost::python::list keys=arg.keys();      boost::python::list keys=arg.keys();
684      for (int i=0;i<num_data;++i) {      for (int i=0;i<num_data;++i) {
685           Data& d=boost::python::extract<Data&>(arg[keys[i]]);           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
686           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
687               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");
688           data[i]=d.getDataC();           data[i]=d.getDataC();
# Line 1041  int MeshAdapter::getSystemMatrixTypeId(c Line 893  int MeshAdapter::getSystemMatrixTypeId(c
893     return out;     return out;
894  }  }
895    
896  Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
897  {  {
898    return continuousFunction(asAbstractContinuousDomain()).getX();    return continuousFunction(asAbstractContinuousDomain()).getX();
899  }  }
900    
901  Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
902  {  {
903    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();    return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
904  }  }
905    
906  Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
907  {  {
908    return function(asAbstractContinuousDomain()).getSize();    return function(asAbstractContinuousDomain()).getSize();
909  }  }
910    
911  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
912  {  {
913    int* tagList;    int out=0;
914    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
915    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
916    return tagList[sampleNo];    case(Nodes):
917        out=mesh->Nodes->Tag[sampleNo];
918        break;
919      case(Elements):
920        out=mesh->Elements->Tag[sampleNo];
921        break;
922      case(FaceElements):
923        out=mesh->FaceElements->Tag[sampleNo];
924        break;
925      case(Points):
926        out=mesh->Points->Tag[sampleNo];
927        break;
928      case(ContactElementsZero):
929        out=mesh->ContactElements->Tag[sampleNo];
930        break;
931      case(ContactElementsOne):
932        out=mesh->ContactElements->Tag[sampleNo];
933        break;
934      case(DegreesOfFreedom):
935        break;
936      case(ReducedDegreesOfFreedom):
937        break;
938      default:
939        stringstream temp;
940        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
941        throw FinleyAdapterException(temp.str());
942        break;
943      }
944      return out;
945  }  }
   
946  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
947  {  {
948    int* referenceNoList;    int out=0,i;
949    int numReferenceNo;    Finley_Mesh* mesh=m_finleyMesh.get();
950    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    switch (functionSpaceType) {
951    return referenceNoList[sampleNo];    case(Nodes):
952        if (mesh->Nodes!=NULL) {
953          out=mesh->Nodes->Id[sampleNo];
954          break;
955        }
956      case(Elements):
957        out=mesh->Elements->Id[sampleNo];
958        break;
959      case(FaceElements):
960        out=mesh->FaceElements->Id[sampleNo];
961        break;
962      case(Points):
963        out=mesh->Points->Id[sampleNo];
964        break;
965      case(ContactElementsZero):
966        out=mesh->ContactElements->Id[sampleNo];
967        break;
968      case(ContactElementsOne):
969        out=mesh->ContactElements->Id[sampleNo];
970        break;
971      case(DegreesOfFreedom):
972        for (i=0;i<mesh->Nodes->numNodes; ++i) {
973           if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
974              out=mesh->Nodes->Id[i];
975              break;
976           }
977        }
978        break;
979      case(ReducedDegreesOfFreedom):
980        for (i=0;i<mesh->Nodes->numNodes; ++i) {
981           if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
982              out=mesh->Nodes->Id[i];
983              break;
984           }
985        }
986        break;
987      default:
988        stringstream temp;
989        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
990        throw FinleyAdapterException(temp.str());
991        break;
992      }
993      return out;
994  }  }
995    
996  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.472  
changed lines
  Added in v.682

  ViewVC Help
Powered by ViewVC 1.1.26