/[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 3269 by jfenwick, Wed Oct 13 03:21:50 2010 UTC revision 3675 by jfenwick, Thu Nov 17 00:53:38 2011 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
14    #include <pasowrap/PasoException.h>
15    #include <pasowrap/TransportProblemAdapter.h>
16  #include "MeshAdapter.h"  #include "MeshAdapter.h"
17  #include "escript/Data.h"  #include "escript/Data.h"
18  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
19  #ifdef USE_NETCDF  #ifdef USE_NETCDF
20  #include <netcdfcpp.h>  #include <netcdfcpp.h>
21  #endif  #endif
 #ifdef ESYS_MPI  
 #include <mpi.h>  
 #include "esysUtils/Esys_MPI.h"  
 #endif  
22  extern "C" {  extern "C" {
23  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
24  }  }
25    
26    #include <boost/python/import.hpp>
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    using namespace paso;
31    
32  namespace finley {  namespace finley {
33    
# Line 649  int MeshAdapter::getReducedSolutionCode( Line 649  int MeshAdapter::getReducedSolutionCode(
649     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
650  }  }
651    
652  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
653  {  {
654     return Points;     return Points;
655  }  }
# Line 773  void MeshAdapter::addPDEToSystem( Line 773  void MeshAdapter::addPDEToSystem(
773                                   AbstractSystemMatrix& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
774                                   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,
775                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
776                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
777                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
778  {  {
779     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
780     if (smat==0)     if (smat==0)
# Line 791  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
793     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
794     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
795       escriptDataC _d_dirac=d_dirac.getDataC();
796       escriptDataC _y_dirac=y_dirac.getDataC();
797    
798     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
799    
# Line 802  void MeshAdapter::addPDEToSystem( Line 805  void MeshAdapter::addPDEToSystem(
805    
806     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 );
807     checkFinleyError();     checkFinleyError();
808    
809        Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
810       checkFinleyError();
811  }  }
812    
813  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
814                                          escript::Data& mat,                                          escript::Data& mat,
815                                          const escript::Data& D,                                          const escript::Data& D,
816                                          const escript::Data& d) const                                          const escript::Data& d,
817                                            const escript::Data& d_dirac,
818                                            const bool useHRZ) const
819  {  {
820     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
821     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
822     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
823       escriptDataC _d_dirac=d_dirac.getDataC();
824    
825     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
826    
827     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
828     checkFinleyError();     checkFinleyError();
829        
830     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
831       checkFinleyError();
832    
833       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
834     checkFinleyError();     checkFinleyError();
835    
836  }  }
# Line 827  void  MeshAdapter::addPDEToLumpedSystem( Line 839  void  MeshAdapter::addPDEToLumpedSystem(
839  //  //
840  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
841  //  //
842  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
843  {  {
844     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
845    
# Line 836  void MeshAdapter::addPDEToRHS( escript:: Line 848  void MeshAdapter::addPDEToRHS( escript::
848     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
849     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
850     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
851       escriptDataC _y_dirac=y_dirac.getDataC();
852    
853     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 );
854     checkFinleyError();     checkFinleyError();
# Line 845  void MeshAdapter::addPDEToRHS( escript:: Line 858  void MeshAdapter::addPDEToRHS( escript::
858    
859     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 );
860     checkFinleyError();     checkFinleyError();
861    
862       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
863       checkFinleyError();
864    
865  }  }
866  //  //
867  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
# Line 854  void MeshAdapter::addPDEToTransportProbl Line 871  void MeshAdapter::addPDEToTransportProbl
871                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
872                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
873                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
874                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
875                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
876  {  {
877     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);     paso::TransportProblemAdapter* tpa=dynamic_cast<paso::TransportProblemAdapter*>(&tp);
878     if (tpa==0)     if (tpa==0)
879     {     {
880      throw FinleyAdapterException("finley only supports its own transport problems.");      throw FinleyAdapterException("finley only supports its own transport problems.");
# Line 877  void MeshAdapter::addPDEToTransportProbl Line 895  void MeshAdapter::addPDEToTransportProbl
895     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
896     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
897     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
898       escriptDataC _d_dirac=d_dirac.getDataC();
899       escriptDataC _y_dirac=y_dirac.getDataC();
900    
901    
902     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
903     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
# Line 892  void MeshAdapter::addPDEToTransportProbl Line 913  void MeshAdapter::addPDEToTransportProbl
913    
914     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 );
915     checkFinleyError();     checkFinleyError();
916    
917       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
918       checkFinleyError();
919    
920  }  }
921    
922  //  //
# Line 1076  void MeshAdapter::interpolateOnDomain(es Line 1101  void MeshAdapter::interpolateOnDomain(es
1101        break;        break;
1102        case(Points):        case(Points):
1103        if (getMPISize()>1) {        if (getMPISize()>1) {
1104           escript::Data temp=escript::Data( in,  continuousFunction(*this) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1105           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1106        } else {        } else {
1107           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1108        }        }
# Line 1550  void MeshAdapter::saveDX(const string& f Line 1575  void MeshAdapter::saveDX(const string& f
1575  //  //
1576  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
1577  {  {
1578     int num_data;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1579     char **names;      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1580     escriptDataC *data;                metadata, metadata_schema, arg);
    escriptDataC **ptr_data;  
   
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());  
    checkFinleyError();  
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1581  }  }
1582    
1583  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
# Line 2428  ReferenceElementSetWrapper::~ReferenceEl Line 2439  ReferenceElementSetWrapper::~ReferenceEl
2439    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
2440  }  }
2441    
2442    // points will be flattened
2443    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2444    {
2445          const int dim = getDim();
2446          int numPoints=points.size()/dim;
2447          int numTags=tags.size();
2448          Finley_Mesh* mesh=m_finleyMesh.get();
2449          
2450          if ( points.size() % dim != 0 )
2451          {
2452        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2453          }
2454          
2455          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2456         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2457          
2458          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2459          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2460          
2461          for (int i=0;i<numPoints;++i) {
2462           points_ptr[ i * dim     ] = points[i * dim ];
2463           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2464           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2465               tags_ptr[i]=tags[i];
2466          }
2467          
2468          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2469          checkPasoError();
2470          
2471          TMPMEMFREE(points_ptr);
2472          TMPMEMFREE(tags_ptr);
2473    }
2474    
2475    
2476    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2477    // {
2478    //       const int dim = getDim();
2479    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2480    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2481    //       Finley_Mesh* mesh=m_finleyMesh.get();
2482    //      
2483    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2484    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2485    //      
2486    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2487    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2488    //      
2489    //       for (int i=0;i<numPoints;++i) {
2490    //     int tag_id=-1;
2491    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2492    //     if  ( numComps !=   dim ) {
2493    //                stringstream temp;            
2494    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2495    //                throw FinleyAdapterException(temp.str());
2496    //     }
2497    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2498    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2499    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2500    //    
2501    //     if ( numTags > 0) {
2502    //            boost::python::extract<string> ex_str(tags[i]);
2503    //        if  ( ex_str.check() ) {
2504    //            tag_id=getTag( ex_str());
2505    //        } else {
2506    //             boost::python::extract<int> ex_int(tags[i]);
2507    //             if ( ex_int.check() ) {
2508    //                 tag_id=ex_int();
2509    //             } else {
2510    //              stringstream temp;          
2511    //                  temp << "Error - unable to extract tag for point " << i;
2512    //              throw FinleyAdapterException(temp.str());
2513    //            }
2514    //        }
2515    //     }      
2516    //            tags_ptr[i]=tag_id;
2517    //       }
2518    //      
2519    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2520    //       checkPasoError();
2521    //      
2522    //       TMPMEMFREE(points_ptr);
2523    //       TMPMEMFREE(tags_ptr);
2524    // }
2525    
2526    /*
2527    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2528    {  
2529        boost::python::list points =  boost::python::list();
2530        boost::python::list tags =  boost::python::list();
2531        points.append(point);
2532        tags.append(tag);
2533        addDiracPoints(points, tags);
2534    }
2535    */
2536    
2537    /*
2538    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2539    {
2540            boost::python::list points =   boost::python::list();
2541            boost::python::list tags =   boost::python::list();
2542            points.append(point);
2543            tags.append(tag);
2544            addDiracPoints(points, tags);
2545    }
2546    */
2547  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.3269  
changed lines
  Added in v.3675

  ViewVC Help
Powered by ViewVC 1.1.26