/[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 3546 by jfenwick, Thu Jul 28 01:32:20 2011 UTC revision 3998 by caltinay, Thu Sep 27 01:17:28 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
# Line 26  extern "C" { Line 29  extern "C" {
29    
30  using namespace std;  using namespace std;
31  using namespace escript;  using namespace escript;
32    using namespace paso;
33    
34  namespace finley {  namespace finley {
35    
# Line 777  void MeshAdapter::addPDEToSystem( Line 781  void MeshAdapter::addPDEToSystem(
781     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
782     if (smat==0)     if (smat==0)
783     {     {
784      throw FinleyAdapterException("finley only supports its own system matrices.");      throw FinleyAdapterException("finley only supports Paso system matrices.");
785     }     }
786     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
787     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
# Line 875  void MeshAdapter::addPDEToTransportProbl Line 879  void MeshAdapter::addPDEToTransportProbl
879     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
880     if (tpa==0)     if (tpa==0)
881     {     {
882      throw FinleyAdapterException("finley only supports its own transport problems.");      throw FinleyAdapterException("finley only supports Paso transport problems.");
883     }     }
884    
885    
# Line 1510  void MeshAdapter::setNewX(const escript: Line 1514  void MeshAdapter::setNewX(const escript:
1514         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1515         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1516     } else {     } else {
1517         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );         throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1518    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1519         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1520         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);*/
1521     }     }
1522     checkFinleyError();     checkFinleyError();
1523  }  }
# Line 1564  void MeshAdapter::saveDX(const string& f Line 1569  void MeshAdapter::saveDX(const string& f
1569        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1570     }     }
1571     TMPMEMFREE(names);     TMPMEMFREE(names);
   
    return;  
1572  }  }
1573    
1574  //  //
# Line 1580  void MeshAdapter::saveVTK(const string& Line 1583  void MeshAdapter::saveVTK(const string&
1583    
1584  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1585  {  {
1586        if (getMPISize() > 1) {
1587  #ifdef ESYS_MPI  #ifdef ESYS_MPI
1588      index_t myFirstNode=0, myLastNode=0, k=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1589      index_t* globalNodeIndex=0;          index_t* globalNodeIndex=0;
1590      Finley_Mesh* mesh_p=m_finleyMesh.get();          Finley_Mesh* mesh_p=m_finleyMesh.get();
1591      if (fs_code == FINLEY_REDUCED_NODES)          /*
1592      {           * this method is only used by saveDataCSV which would use the returned
1593      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1594      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES)
1595      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          {
1596      }              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1597      else              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1598      {              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1599      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);          }
1600      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          else
1601      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          */
1602      }          if (fs_code == FINLEY_NODES)
1603      k=globalNodeIndex[id];          {
1604      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );              myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1605                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1606                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1607            }
1608            else
1609            {
1610                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1611            }
1612    
1613            k=globalNodeIndex[id];
1614            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1615  #endif  #endif
1616        }
1617      return true;      return true;
1618  }  }
1619    
1620    
   
1621  //  //
1622  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1623  //  //
# Line 1663  ASM_ptr MeshAdapter::newSystemMatrix( Line 1677  ASM_ptr MeshAdapter::newSystemMatrix(
1677  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1678  //  //
1679  ATP_ptr MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const bool useBackwardEuler,  
1680                                                           const int blocksize,                                                           const int blocksize,
1681                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1682                                                           const int type) const                                                           const int type) const
# Line 1686  ATP_ptr MeshAdapter::newTransportProblem Line 1699  ATP_ptr MeshAdapter::newTransportProblem
1699     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1700     checkFinleyError();     checkFinleyError();
1701     Paso_TransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1702     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1703     checkPasoError();     checkPasoError();
1704     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1705     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1706     return ATP_ptr(tpa);     return ATP_ptr(tpa);
1707  //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);  //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1708  }  }
1709    
1710  //  //
# Line 2050  bool MeshAdapter::operator!=(const Abstr Line 2063  bool MeshAdapter::operator!=(const Abstr
2063  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2064  {  {
2065     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2066     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2067     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2068  }  }
2069    
2070  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2071  {  {
2072     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2073     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2074     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2075  }  }
2076    
2077  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
# Line 2240  void MeshAdapter::setTagMap(const string Line 2252  void MeshAdapter::setTagMap(const string
2252  {  {
2253     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2254     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2255     checkPasoError();     checkFinleyError();
2256     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2257  }  }
2258    
# Line 2249  int MeshAdapter::getTag(const string& na Line 2261  int MeshAdapter::getTag(const string& na
2261     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2262     int tag=0;     int tag=0;
2263     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2264     checkPasoError();     checkFinleyError();
2265     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2266     return tag;     return tag;
2267  }  }
# Line 2437  ReferenceElementSetWrapper::~ReferenceEl Line 2449  ReferenceElementSetWrapper::~ReferenceEl
2449    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
2450  }  }
2451    
2452  void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const  // points will be flattened
2453    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2454  {  {
2455        const int dim = getDim();        const int dim = getDim();
2456        int numPoints=boost::python::extract<int>(points.attr("__len__")());        int numPoints=points.size()/dim;
2457        int numTags=boost::python::extract<int>(tags.attr("__len__")());        int numTags=tags.size();
2458        Finley_Mesh* mesh=m_finleyMesh.get();        Finley_Mesh* mesh=m_finleyMesh.get();
2459                
2460          if ( points.size() % dim != 0 )
2461          {
2462        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2463          }
2464          
2465        if  ( (numTags > 0) && ( numPoints !=  numTags ) )        if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2466       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2467                
# Line 2451  void MeshAdapter:: addDiracPoints(const Line 2469  void MeshAdapter:: addDiracPoints(const
2469        int*    tags_ptr= TMPMEMALLOC(numPoints, int);        int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2470                
2471        for (int i=0;i<numPoints;++i) {        for (int i=0;i<numPoints;++i) {
2472         int tag_id=-1;         points_ptr[ i * dim     ] = points[i * dim ];
2473         int numComps=boost::python::extract<int>(points[i].attr("__len__")());         if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2474         if  ( numComps !=   dim ) {         if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2475                 stringstream temp;                       tags_ptr[i]=tags[i];
                temp << "Error - illegal number of components " << numComps << " for point " << i;                
                throw FinleyAdapterException(temp.str());  
        }  
        points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);  
        if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);  
        if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);  
         
        if ( numTags > 0) {  
               boost::python::extract<string> ex_str(tags[i]);  
           if  ( ex_str.check() ) {  
               tag_id=getTag( ex_str());  
           } else {  
                boost::python::extract<int> ex_int(tags[i]);  
                if ( ex_int.check() ) {  
                    tag_id=ex_int();  
                } else {  
                 stringstream temp;            
                     temp << "Error - unable to extract tag for point " << i;  
                 throw FinleyAdapterException(temp.str());  
               }  
           }  
        }        
            tags_ptr[i]=tag_id;  
2476        }        }
2477                
2478        Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);        Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2479        checkPasoError();        checkFinleyError();
2480                
2481        TMPMEMFREE(points_ptr);        TMPMEMFREE(points_ptr);
2482        TMPMEMFREE(tags_ptr);        TMPMEMFREE(tags_ptr);
2483  }  }
2484    
2485    
2486    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2487    // {
2488    //       const int dim = getDim();
2489    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2490    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2491    //       Finley_Mesh* mesh=m_finleyMesh.get();
2492    //      
2493    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2494    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2495    //      
2496    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2497    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2498    //      
2499    //       for (int i=0;i<numPoints;++i) {
2500    //     int tag_id=-1;
2501    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2502    //     if  ( numComps !=   dim ) {
2503    //                stringstream temp;            
2504    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2505    //                throw FinleyAdapterException(temp.str());
2506    //     }
2507    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2508    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2509    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2510    //    
2511    //     if ( numTags > 0) {
2512    //            boost::python::extract<string> ex_str(tags[i]);
2513    //        if  ( ex_str.check() ) {
2514    //            tag_id=getTag( ex_str());
2515    //        } else {
2516    //             boost::python::extract<int> ex_int(tags[i]);
2517    //             if ( ex_int.check() ) {
2518    //                 tag_id=ex_int();
2519    //             } else {
2520    //              stringstream temp;          
2521    //                  temp << "Error - unable to extract tag for point " << i;
2522    //              throw FinleyAdapterException(temp.str());
2523    //            }
2524    //        }
2525    //     }      
2526    //            tags_ptr[i]=tag_id;
2527    //       }
2528    //      
2529    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2530    //       checkPasoError();
2531    //      
2532    //       TMPMEMFREE(points_ptr);
2533    //       TMPMEMFREE(tags_ptr);
2534    // }
2535    
2536    /*
2537  void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const  void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2538  {    {  
2539      boost::python::list points =  boost::python::list();      boost::python::list points =  boost::python::list();
# Line 2495  void MeshAdapter:: addDiracPoint( const Line 2542  void MeshAdapter:: addDiracPoint( const
2542      tags.append(tag);      tags.append(tag);
2543      addDiracPoints(points, tags);      addDiracPoints(points, tags);
2544  }  }
2545    */
2546    
2547    /*
2548  void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const  void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2549  {  {
2550          boost::python::list points =   boost::python::list();          boost::python::list points =   boost::python::list();
# Line 2502  void MeshAdapter:: addDiracPointWithTagN Line 2552  void MeshAdapter:: addDiracPointWithTagN
2552          points.append(point);          points.append(point);
2553          tags.append(tag);          tags.append(tag);
2554          addDiracPoints(points, tags);          addDiracPoints(points, tags);
2555  }    }
2556    */
2557  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.3546  
changed lines
  Added in v.3998

  ViewVC Help
Powered by ViewVC 1.1.26