/[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 3675 by jfenwick, Thu Nov 17 00:53:38 2011 UTC revision 4286 by caltinay, Thu Mar 7 04:28:11 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2013 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>  #include <pasowrap/PasoException.h>
17  #include <pasowrap/TransportProblemAdapter.h>  #include <pasowrap/TransportProblemAdapter.h>
# Line 779  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 874  void MeshAdapter::addPDEToTransportProbl Line 876  void MeshAdapter::addPDEToTransportProbl
876                                             const escript::Data& d_contact,const escript::Data& y_contact,                                             const escript::Data& d_contact,const escript::Data& y_contact,
877                                             const escript::Data& d_dirac, const escript::Data& y_dirac) const                                             const escript::Data& d_dirac, const escript::Data& y_dirac) const
878  {  {
879     paso::TransportProblemAdapter* tpa=dynamic_cast<paso::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 1512  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  }  }
1524    
 //  
 // Helper for the save* methods. Extracts optional data variable names and the  
 // corresponding pointers from python dictionary. Caller must free arrays.  
 //  
 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const  
 {  
    numData = boost::python::extract<int>(arg.attr("__len__")());  
    /* win32 refactor */  
    names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;  
    data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;  
    dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;  
   
    boost::python::list keys=arg.keys();  
    for (int i=0; i<numData; ++i) {  
       string n=boost::python::extract<string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)  
          throw FinleyAdapterException("Error: Data must be defined on same Domain");  
       data[i] = d.getDataC();  
       dataPtr[i] = &(data[i]);  
       names[i] = TMPMEMALLOC(n.length()+1, char);  
       strcpy(names[i], n.c_str());  
    }  
 }  
   
 //  
 // saves mesh and optionally data arrays in openDX format  
 //  
 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const  
 {  
    int num_data;  
    char **names;  
    escriptDataC *data;  
    escriptDataC **ptr_data;  
   
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);  
    checkFinleyError();  
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
   
    return;  
 }  
   
 //  
 // saves mesh and optionally data arrays in VTK format  
 //  
 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const  
 {  
     boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");  
     pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),  
               metadata, metadata_schema, arg);  
 }  
   
1525  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1526  {  {
1527        if (getMPISize() > 1) {
1528  #ifdef ESYS_MPI  #ifdef ESYS_MPI
1529      index_t myFirstNode=0, myLastNode=0, k=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1530      index_t* globalNodeIndex=0;          index_t* globalNodeIndex=0;
1531      Finley_Mesh* mesh_p=m_finleyMesh.get();          Finley_Mesh* mesh_p=m_finleyMesh.get();
1532      if (fs_code == FINLEY_REDUCED_NODES)          /*
1533      {           * this method is only used by saveDataCSV which would use the returned
1534      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1535      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES)
1536      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          {
1537      }              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1538      else              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1539      {              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1540      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);          }
1541      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          else
1542      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          */
1543      }          if (fs_code == FINLEY_NODES)
1544      k=globalNodeIndex[id];          {
1545      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );              myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1546                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1547                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1548            }
1549            else
1550            {
1551                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1552            }
1553    
1554            k=globalNodeIndex[id];
1555            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1556  #endif  #endif
1557        }
1558      return true;      return true;
1559  }  }
1560    
1561    
   
1562  //  //
1563  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1564  //  //
# Line 1624  ASM_ptr MeshAdapter::newSystemMatrix( Line 1577  ASM_ptr MeshAdapter::newSystemMatrix(
1577        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1578     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1579     if (col_domain!=*this)     if (col_domain!=*this)
1580        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1581     // is the function space type right     // is the function space type right
1582     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1583        reduceRowOrder=0;        reduceRowOrder=0;
# Line 1665  ASM_ptr MeshAdapter::newSystemMatrix( Line 1618  ASM_ptr MeshAdapter::newSystemMatrix(
1618  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1619  //  //
1620  ATP_ptr MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const bool useBackwardEuler,  
1621                                                           const int blocksize,                                                           const int blocksize,
1622                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1623                                                           const int type) const                                                           const int type) const
# Line 1688  ATP_ptr MeshAdapter::newTransportProblem Line 1640  ATP_ptr MeshAdapter::newTransportProblem
1640     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1641     checkFinleyError();     checkFinleyError();
1642     Paso_TransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1643     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1644     checkPasoError();     checkPasoError();
1645     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1646     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1647     return ATP_ptr(tpa);     return ATP_ptr(tpa);
1648  //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);  //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1649  }  }
1650    
1651  //  //
# Line 2029  bool MeshAdapter::probeInterpolationOnDo Line 1981  bool MeshAdapter::probeInterpolationOnDo
1981     return false;     return false;
1982  }  }
1983    
1984    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1985    {
1986        if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1987        {
1988            return 1;
1989        }
1990        else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1991        {
1992            return -1;
1993        }
1994        return 0;
1995    }  
1996      
1997  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1998  {  {
1999     return false;     return false;
# Line 2052  bool MeshAdapter::operator!=(const Abstr Line 2017  bool MeshAdapter::operator!=(const Abstr
2017  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
2018  {  {
2019     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2020     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2021     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2022  }  }
2023    
2024  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
2025  {  {
2026     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2027     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2028     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2029  }  }
2030    
2031  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
# Line 2242  void MeshAdapter::setTagMap(const string Line 2206  void MeshAdapter::setTagMap(const string
2206  {  {
2207     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2208     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2209     checkPasoError();     checkFinleyError();
2210     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2211  }  }
2212    
# Line 2251  int MeshAdapter::getTag(const string& na Line 2215  int MeshAdapter::getTag(const string& na
2215     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2216     int tag=0;     int tag=0;
2217     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2218     checkPasoError();     checkFinleyError();
2219     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2220     return tag;     return tag;
2221  }  }
# Line 2466  void MeshAdapter:: addDiracPoints(const Line 2430  void MeshAdapter:: addDiracPoints(const
2430        }        }
2431                
2432        Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);        Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2433        checkPasoError();        checkFinleyError();
2434                
2435        TMPMEMFREE(points_ptr);        TMPMEMFREE(points_ptr);
2436        TMPMEMFREE(tags_ptr);        TMPMEMFREE(tags_ptr);

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

  ViewVC Help
Powered by ViewVC 1.1.26