/[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 3344 by caltinay, Thu Nov 11 23:26:52 2010 UTC revision 3558 by jfenwick, Fri Aug 26 06:09:02 2011 UTC
# Line 23  extern "C" { Line 23  extern "C" {
23  }  }
24    
25  #include <boost/python/import.hpp>  #include <boost/python/import.hpp>
 #include <boost/python/tuple.hpp>  
26    
27  using namespace std;  using namespace std;
28  using namespace escript;  using namespace escript;
# Line 648  int MeshAdapter::getReducedSolutionCode( Line 647  int MeshAdapter::getReducedSolutionCode(
647     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
648  }  }
649    
650  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
651  {  {
652     return Points;     return Points;
653  }  }
# Line 772  void MeshAdapter::addPDEToSystem( Line 771  void MeshAdapter::addPDEToSystem(
771                                   AbstractSystemMatrix& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::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,
773                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
775                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
776  {  {
777     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
778     if (smat==0)     if (smat==0)
# Line 790  void MeshAdapter::addPDEToSystem( Line 790  void MeshAdapter::addPDEToSystem(
790     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
791     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
792     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
793       escriptDataC _d_dirac=d_dirac.getDataC();
794       escriptDataC _y_dirac=y_dirac.getDataC();
795    
796     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
797    
# Line 801  void MeshAdapter::addPDEToSystem( Line 803  void MeshAdapter::addPDEToSystem(
803    
804     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
805     checkFinleyError();     checkFinleyError();
806    
807        Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
808       checkFinleyError();
809  }  }
810    
811  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
812                                          escript::Data& mat,                                          escript::Data& mat,
813                                          const escript::Data& D,                                          const escript::Data& D,
814                                          const escript::Data& d) const                                          const escript::Data& d,
815                                            const escript::Data& d_dirac,
816                                            const bool useHRZ) const
817  {  {
818     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
819     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
820     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
821       escriptDataC _d_dirac=d_dirac.getDataC();
822    
823     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
824    
825     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
826     checkFinleyError();     checkFinleyError();
827        
828     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
829       checkFinleyError();
830    
831       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
832     checkFinleyError();     checkFinleyError();
833    
834  }  }
# Line 826  void  MeshAdapter::addPDEToLumpedSystem( Line 837  void  MeshAdapter::addPDEToLumpedSystem(
837  //  //
838  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
839  //  //
840  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact,  const escript::Data& y_dirac) const
841  {  {
842     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
843    
# Line 835  void MeshAdapter::addPDEToRHS( escript:: Line 846  void MeshAdapter::addPDEToRHS( escript::
846     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
847     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
848     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
849       escriptDataC _y_dirac=y_dirac.getDataC();
850    
851     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
852     checkFinleyError();     checkFinleyError();
# Line 844  void MeshAdapter::addPDEToRHS( escript:: Line 856  void MeshAdapter::addPDEToRHS( escript::
856    
857     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
858     checkFinleyError();     checkFinleyError();
859    
860       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
861       checkFinleyError();
862    
863  }  }
864  //  //
865  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
# Line 853  void MeshAdapter::addPDEToTransportProbl Line 869  void MeshAdapter::addPDEToTransportProbl
869                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
870                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
871                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
872                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
873                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
874  {  {
875     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
876     if (tpa==0)     if (tpa==0)
# Line 876  void MeshAdapter::addPDEToTransportProbl Line 893  void MeshAdapter::addPDEToTransportProbl
893     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
894     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
895     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
896       escriptDataC _d_dirac=d_dirac.getDataC();
897       escriptDataC _y_dirac=y_dirac.getDataC();
898    
899    
900     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
901     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
# Line 891  void MeshAdapter::addPDEToTransportProbl Line 911  void MeshAdapter::addPDEToTransportProbl
911    
912     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
913     checkFinleyError();     checkFinleyError();
914    
915       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
916       checkFinleyError();
917    
918  }  }
919    
920  //  //
# Line 1075  void MeshAdapter::interpolateOnDomain(es Line 1099  void MeshAdapter::interpolateOnDomain(es
1099        break;        break;
1100        case(Points):        case(Points):
1101        if (getMPISize()>1) {        if (getMPISize()>1) {
1102           escript::Data temp=escript::Data( in,  continuousFunction(*this) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1103           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1104        } else {        } else {
1105           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1106        }        }
# Line 1549  void MeshAdapter::saveDX(const string& f Line 1573  void MeshAdapter::saveDX(const string& f
1573  //  //
1574  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1575  {  {
1576      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("saveVTK");      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1577      pySaveVTK(*boost::python::make_tuple(filename,      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1578                 const_cast<MeshAdapter*>(this)->getPtr(),                metadata, metadata_schema, arg);
                metadata, metadata_schema), **arg);  
1579  }  }
1580    
1581  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
# Line 2414  ReferenceElementSetWrapper::~ReferenceEl Line 2437  ReferenceElementSetWrapper::~ReferenceEl
2437    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
2438  }  }
2439    
2440    // points will be flattened
2441    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2442    {
2443          const int dim = getDim();
2444          int numPoints=points.size()/dim;
2445          int numTags=tags.size();
2446          Finley_Mesh* mesh=m_finleyMesh.get();
2447          
2448          if ( points.size() % dim != 0 )
2449          {
2450        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2451          }
2452          
2453          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2454         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2455          
2456          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2457          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2458          
2459          for (int i=0;i<numPoints;++i) {
2460           points_ptr[ i * dim     ] = points[i * dim ];
2461           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2462           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2463               tags_ptr[i]=tags[i];
2464          }
2465          
2466          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2467          checkPasoError();
2468          
2469          TMPMEMFREE(points_ptr);
2470          TMPMEMFREE(tags_ptr);
2471    }
2472    
2473    
2474    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2475    // {
2476    //       const int dim = getDim();
2477    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2478    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2479    //       Finley_Mesh* mesh=m_finleyMesh.get();
2480    //      
2481    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2482    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2483    //      
2484    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2485    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2486    //      
2487    //       for (int i=0;i<numPoints;++i) {
2488    //     int tag_id=-1;
2489    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2490    //     if  ( numComps !=   dim ) {
2491    //                stringstream temp;            
2492    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2493    //                throw FinleyAdapterException(temp.str());
2494    //     }
2495    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2496    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2497    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2498    //    
2499    //     if ( numTags > 0) {
2500    //            boost::python::extract<string> ex_str(tags[i]);
2501    //        if  ( ex_str.check() ) {
2502    //            tag_id=getTag( ex_str());
2503    //        } else {
2504    //             boost::python::extract<int> ex_int(tags[i]);
2505    //             if ( ex_int.check() ) {
2506    //                 tag_id=ex_int();
2507    //             } else {
2508    //              stringstream temp;          
2509    //                  temp << "Error - unable to extract tag for point " << i;
2510    //              throw FinleyAdapterException(temp.str());
2511    //            }
2512    //        }
2513    //     }      
2514    //            tags_ptr[i]=tag_id;
2515    //       }
2516    //      
2517    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2518    //       checkPasoError();
2519    //      
2520    //       TMPMEMFREE(points_ptr);
2521    //       TMPMEMFREE(tags_ptr);
2522    // }
2523    
2524    /*
2525    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2526    {  
2527        boost::python::list points =  boost::python::list();
2528        boost::python::list tags =  boost::python::list();
2529        points.append(point);
2530        tags.append(tag);
2531        addDiracPoints(points, tags);
2532    }
2533    */
2534    
2535    /*
2536    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2537    {
2538            boost::python::list points =   boost::python::list();
2539            boost::python::list tags =   boost::python::list();
2540            points.append(point);
2541            tags.append(tag);
2542            addDiracPoints(points, tags);
2543    }
2544    */
2545  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.3344  
changed lines
  Added in v.3558

  ViewVC Help
Powered by ViewVC 1.1.26