Panzer_BlockedEpetraLinearObjFactory.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//           Panzer: A partial differential equation assembly
//       engine for strongly coupled complex multiphysics systems
//                 Copyright (2011) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
// Eric C. Cyr (eccyr@sandia.gov)
// ***********************************************************************
// @HEADER
 
#ifndef   __Panzer_BlockedEpetraLinearObjFactory_hpp__
#define   __Panzer_BlockedEpetraLinearObjFactory_hpp__
 
#include <map>
 
// Epetra includes
#include "Epetra_Map.h"
#include "Epetra_CrsGraph.h"
#include "Epetra_Import.h"
#include "Epetra_Export.h"
 
#include "PanzerDiscFE_config.hpp"
#include "Panzer_GlobalIndexer.hpp"
#include "Panzer_LinearObjFactory.hpp"
#include "Panzer_EpetraLinearObjContainer.hpp"
#include "Panzer_BlockedEpetraLinearObjContainer.hpp"
#include "Panzer_BlockedDOFManager.hpp"
#include "Panzer_CloneableEvaluator.hpp"
 
#include "Panzer_GatherOrientation.hpp"
#include "Panzer_GatherSolution_BlockedEpetra.hpp"
#include "Panzer_GatherTangent_BlockedEpetra.hpp"
#include "Panzer_ScatterResidual_BlockedEpetra.hpp"
#include "Panzer_ScatterDirichletResidual_BlockedEpetra.hpp"
#include "Panzer_ScatterResidual_Epetra.hpp"
#include "Panzer_ScatterDirichletResidual_Epetra.hpp"
#include "Panzer_GatherSolution_Epetra.hpp"
#include "Panzer_GatherTangent_Epetra.hpp"
#include "Panzer_GatherOrientation.hpp"
#include "Panzer_ThyraObjFactory.hpp"
#include "Panzer_HashUtils.hpp"
 
#include "Thyra_BlockedLinearOpBase.hpp"
#include "Thyra_ProductVectorBase.hpp"
 
#include "Teuchos_RCP.hpp"
#include "Teuchos_DefaultMpiComm.hpp"
#include "Teuchos_OpaqueWrapper.hpp"
 
namespace panzer
{
 
template <typename Traits,typename LocalOrdinalT>
class BlockedEpetraLinearObjFactory : public LinearObjFactory<Traits>
                                    , public ThyraObjFactory<double> {
public:
 
   BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
                                 const Teuchos::RCP<const GlobalIndexer> & gidProvider,
                                 bool useDiscreteAdjoint=false);
 
   BlockedEpetraLinearObjFactory(const Teuchos::RCP<const Teuchos::MpiComm<int> > & comm,
                                 const Teuchos::RCP<const GlobalIndexer> & gidProvider,
                                 const Teuchos::RCP<const GlobalIndexer> & colGidProvider,
                                 bool useDiscreteAdjoint=false);
 
   virtual ~BlockedEpetraLinearObjFactory();
 
/*************** Linear object factory methods *******************/
 
   virtual void readVector(const std::string & identifier,LinearObjContainer & loc,int id) const;
 
   virtual void writeVector(const std::string & identifier,const LinearObjContainer & loc,int id) const;
 
   virtual Teuchos::RCP<LinearObjContainer> buildLinearObjContainer() const;
 
   virtual Teuchos::RCP<LinearObjContainer> buildPrimitiveLinearObjContainer() const 
   { return buildLinearObjContainer(); }
 
   virtual Teuchos::RCP<LinearObjContainer> buildGhostedLinearObjContainer() const;
 
   virtual Teuchos::RCP<LinearObjContainer> buildPrimitiveGhostedLinearObjContainer() const 
   { return buildGhostedLinearObjContainer(); }
 
   virtual void globalToGhostContainer(const LinearObjContainer & container,
                                       LinearObjContainer & ghostContainer,int) const;
   virtual void ghostToGlobalContainer(const LinearObjContainer & ghostContainer,
                                       LinearObjContainer & container,int) const;
 
   virtual void adjustForDirichletConditions(const LinearObjContainer & localBCRows,
                                             const LinearObjContainer & globalBCRows,
                                             LinearObjContainer & ghostedObjs,
                                             bool zeroVectorRows=false, bool adjustX=false) const;
 
   virtual void applyDirichletBCs(const LinearObjContainer & counter,
                                  LinearObjContainer & result) const;
 
   virtual Teuchos::RCP<ReadOnlyVector_GlobalEvaluationData> buildReadOnlyDomainContainer() const;
 
   virtual Teuchos::RCP<WriteVector_GlobalEvaluationData> buildWriteDomainContainer() const;
 
   virtual Teuchos::MpiComm<int> getComm() const;
 
   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator> buildScatter() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new ScatterResidual_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0],
                                                                                      colDOFManagerContainer_->getFieldDOFManagers()[0],
                                                                                      useDiscreteAdjoint_));
 
     return Teuchos::rcp(new ScatterResidual_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers(),
                                                                                           colDOFManagerContainer_->getFieldDOFManagers(),
                                                                                           useDiscreteAdjoint_)); 
   }
 
   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGather() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new GatherSolution_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0]));
     return Teuchos::rcp(new GatherSolution_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers())); 
   }
 
   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGatherTangent() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new GatherTangent_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0]));
     return Teuchos::rcp(new GatherTangent_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers())); 
   }
 
   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGatherDomain() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new GatherSolution_Epetra<EvalT,Traits,LocalOrdinalT,int>(colDOFManagerContainer_->getFieldDOFManagers()[0]));
     return Teuchos::rcp(new GatherSolution_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(colDOFManagerContainer_->getFieldDOFManagers())); 
   }
 
   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator > buildGatherOrientation() const
  { return Teuchos::rcp(new GatherOrientation<EvalT,Traits,LocalOrdinalT,panzer::GlobalOrdinal>(rowDOFManagerContainer_->getFieldDOFManagers())); }
 
   template <typename EvalT>
   Teuchos::RCP<panzer::CloneableEvaluator> buildScatterDirichlet() const
   { 
     if(!colDOFManagerContainer_->containsBlockedDOFManager() &&
        !rowDOFManagerContainer_->containsBlockedDOFManager())
       return Teuchos::rcp(new ScatterDirichletResidual_Epetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers()[0],
                                                                                               colDOFManagerContainer_->getFieldDOFManagers()[0])); 
     return Teuchos::rcp(new ScatterDirichletResidual_BlockedEpetra<EvalT,Traits,LocalOrdinalT,int>(rowDOFManagerContainer_->getFieldDOFManagers(),
                                                                                                    colDOFManagerContainer_->getFieldDOFManagers())); 
   }
 
/*************** Generic helper functions for container setup *******************/
   
   void initializeContainer(int,LinearObjContainer & loc) const;
 
   void initializeGhostedContainer(int,LinearObjContainer & loc) const;
 
/*************** Thyra based methods *******************/
 
   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getThyraDomainSpace() const;
 
   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getThyraRangeSpace() const;
 
   Teuchos::RCP<Thyra::VectorBase<double> > getThyraDomainVector() const;
 
   Teuchos::RCP<Thyra::VectorBase<double> > getThyraRangeVector() const;
 
   Teuchos::RCP<Thyra::LinearOpBase<double> > getThyraMatrix() const;
 
   // and now the ghosted versions
 
   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getGhostedThyraDomainSpace() const;
 
   Teuchos::RCP<const Thyra::VectorSpaceBase<double>>
   getGhostedThyraDomainSpace2() const;
 
   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > getGhostedThyraRangeSpace() const;
 
   Teuchos::RCP<Thyra::VectorBase<double> > getGhostedThyraDomainVector() const;
 
   Teuchos::RCP<Thyra::VectorBase<double> > getGhostedThyraRangeVector() const;
 
   Teuchos::RCP<Thyra::LinearOpBase<double> > getGhostedThyraMatrix() const;
 
/*************** Epetra based methods *******************/
 
   virtual const Teuchos::RCP<Epetra_Map> getMap(int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map> getColMap(int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map> getGhostedMap(int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map>
   getGhostedMap2(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map> getGhostedColMap(int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map>
   getGhostedColMap2(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_CrsGraph> getGraph(int i,int j) const;
 
   virtual const Teuchos::RCP<Epetra_CrsGraph> getGhostedGraph(int i,int j) const;
 
   virtual const Teuchos::RCP<Epetra_Import> getGhostedImport(int i) const;
 
   virtual const Teuchos::RCP<Epetra_Import>
   getGhostedImport2(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Import> getGhostedColImport(int i) const;
 
   virtual const Teuchos::RCP<Epetra_Import>
   getGhostedColImport2(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Export> getGhostedExport(int j) const;
 
   virtual const Teuchos::RCP<Epetra_Export>
   getGhostedExport2(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Export> getGhostedColExport(int j) const;
 
   virtual const Teuchos::RCP<Epetra_Export>
   getGhostedColExport2(
     int i) const;
 
   virtual const Teuchos::RCP<const Epetra_Comm> getEpetraComm() const;
 
   Teuchos::RCP<Epetra_CrsMatrix> getEpetraMatrix(int i,int j) const;
   Teuchos::RCP<Epetra_CrsMatrix> getGhostedEpetraMatrix(int i,int j) const;
 
   int getBlockRowCount() const;
 
   int getBlockColCount() const;
 
   Teuchos::RCP<const panzer::BlockedDOFManager> getGlobalIndexer() const
   { return rowDOFManagerContainer_->getBlockedIndexer(); }
 
   Teuchos::RCP<const panzer::GlobalIndexer> getRangeGlobalIndexer() const
   { return rowDOFManagerContainer_->getGlobalIndexer(); }
 
   Teuchos::RCP<const panzer::GlobalIndexer> getDomainGlobalIndexer() const
   { return colDOFManagerContainer_->getGlobalIndexer(); }
 
   const std::vector<Teuchos::RCP<const GlobalIndexer> > & getRangeGlobalIndexers() const
   { return rowDOFManagerContainer_->getFieldDOFManagers(); }
 
   const std::vector<Teuchos::RCP<const GlobalIndexer> > & getDomainGlobalIndexers() const
   { return colDOFManagerContainer_->getFieldDOFManagers(); }
 
   void addExcludedPair(int rowBlock,int colBlock);
 
   void addExcludedPairs(const std::vector<std::pair<int,int> > & exPairs);
 
protected:
 
   void initializeContainer_internal(int mem,ThyraObjContainer<double> & loc) const;
 
   void initializeGhostedContainer_internal(int mem,ThyraObjContainer<double> & loc) const;
 
/*************** Utility class for handling blocked and nonblocked DOF managers *******************/
 
  class DOFManagerContainer {
  public:
    DOFManagerContainer() {} 
    DOFManagerContainer(const Teuchos::RCP<const GlobalIndexer> & ugi) 
    { setGlobalIndexer(ugi); } 
 
    void setGlobalIndexer(const Teuchos::RCP<const GlobalIndexer> & ugi)
    {
      using Teuchos::RCP;
      using Teuchos::rcp_dynamic_cast;
 
      auto blockedDOFManager = rcp_dynamic_cast<const BlockedDOFManager>(ugi);
      auto flatDOFManager    = rcp_dynamic_cast<const GlobalIndexer>(ugi);
 
      if(blockedDOFManager!=Teuchos::null) {
        // set BlockedDOFManager
        blockedDOFManager_ = blockedDOFManager;
 
        // get all GID providers
        auto dofManagers =  blockedDOFManager_->getFieldDOFManagers();
        for(auto itr=dofManagers.begin();itr!=dofManagers.end();++itr)
          gidProviders_.push_back(*itr);
      }
      else if(flatDOFManager!=Teuchos::null) {
        // for absolute clarity, nullify the blockedDOFManager_
        blockedDOFManager_ = Teuchos::null;
 
        // you have only a single GID provider
        gidProviders_.push_back(flatDOFManager);
      }
      else {
        TEUCHOS_ASSERT(false);
      }
    }
 
    int getFieldBlocks() const
    { return Teuchos::as<int>(gidProviders_.size()); }
 
    bool containsBlockedDOFManager() const
    { return blockedDOFManager_ !=Teuchos::null; }
 
    Teuchos::RCP<const BlockedDOFManager> getBlockedIndexer() const
    {
      TEUCHOS_ASSERT(containsBlockedDOFManager());
      return blockedDOFManager_;
    }
 
    Teuchos::RCP<const GlobalIndexer> getGlobalIndexer() const
    {
      if(blockedDOFManager_!=Teuchos::null) 
        return blockedDOFManager_;
 
      TEUCHOS_ASSERT(gidProviders_.size()==1);
      return gidProviders_[0];
    }
 
    const std::vector<Teuchos::RCP<const GlobalIndexer> > & getFieldDOFManagers() const
    { return gidProviders_; }
 
  private:
    Teuchos::RCP<const BlockedDOFManager> blockedDOFManager_;
    std::vector<Teuchos::RCP<const GlobalIndexer> > gidProviders_;
  };
 
/*************** Generic methods/members *******************/
 
   // Get the global indexer associated with a particular block
   Teuchos::RCP<const GlobalIndexer> getGlobalIndexer(int i) const;
 
   Teuchos::RCP<const GlobalIndexer> getColGlobalIndexer(int i) const;
 
   void makeRoomForBlocks(std::size_t blockCnt,std::size_t colBlockCnt=0);
 
   Teuchos::RCP<const DOFManagerContainer> rowDOFManagerContainer_;
   Teuchos::RCP<const DOFManagerContainer> colDOFManagerContainer_;
 
   bool useColGidProviders_;
 
   // which block entries are ignored
   std::unordered_set<std::pair<int,int>,panzer::pair_hash> excludedPairs_;
  
/*************** Thyra based methods/members *******************/
 
   void ghostToGlobalThyraVector(const Teuchos::RCP<const Thyra::VectorBase<double> > & in,
                                 const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const;
   void ghostToGlobalThyraMatrix(const Thyra::LinearOpBase<double> & in,Thyra::LinearOpBase<double> & out) const;
   void globalToGhostThyraVector(const Teuchos::RCP<const Thyra::VectorBase<double> > & in,
                                 const Teuchos::RCP<Thyra::VectorBase<double> > & out,bool col) const;
 
   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > rangeSpace_;
   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > domainSpace_;
 
   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > ghostedRangeSpace_;
   mutable Teuchos::RCP<const Thyra::VectorSpaceBase<double> > ghostedDomainSpace_;
 
/*************** Epetra based methods/members *******************/
 
   void adjustForDirichletConditions(const Epetra_Vector & local_bcs,
                                     const Epetra_Vector & global_bcs,
                                     const Teuchos::Ptr<Epetra_Vector> & f,
                                     const Teuchos::Ptr<Epetra_CrsMatrix> & A,
                                     bool zeroVectorRows) const;
 
   void ghostToGlobalEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const;
   void globalToGhostEpetraVector(int i,const Epetra_Vector & in,Epetra_Vector & out,bool col) const;
   void ghostToGlobalEpetraMatrix(int blockRow,const Epetra_CrsMatrix & in,Epetra_CrsMatrix & out) const;
 
   // get the map from the matrix
 
   virtual const Teuchos::RCP<Epetra_Map>
   buildMap(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map>
   buildGhostedMap(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map>
   buildGhostedMap2(
     int i) const;
 
   // get the map from the matrix
 
   virtual const Teuchos::RCP<Epetra_Map>
   buildColMap(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map>
   buildColGhostedMap(
     int i) const;
 
   virtual const Teuchos::RCP<Epetra_Map>
   buildColGhostedMap2(
     int i) const;
 
   // get the graph of the crs matrix
   virtual const Teuchos::RCP<Epetra_CrsGraph> buildGraph(int i,int j) const;
   virtual const Teuchos::RCP<Epetra_CrsGraph> buildGhostedGraph(int i,int j,bool optimizeStorage) const;
   virtual const Teuchos::RCP<Epetra_CrsGraph> buildFilteredGhostedGraph(int i,int j) const;
 
   // storage for Epetra graphs and maps
   Teuchos::RCP<const Epetra_Comm> eComm_;
   Teuchos::RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm_;
   Teuchos::RCP<Teuchos::MpiComm<int> > tComm_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Map>> maps_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Map>> ghostedMaps_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Map>> ghostedMaps2_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Import>> importers_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Import>> importers2_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Export>> exporters_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Map>> colMaps_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Map>> colGhostedMaps_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Map>> colGhostedMaps2_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Import>> colImporters_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Import>> colImporters2_;
 
   mutable std::vector<Teuchos::RCP<Epetra_Export>> colExporters_;
 
   mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> graphs_ ;
   mutable std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsGraph>,panzer::pair_hash> ghostedGraphs_;
 
   bool useDiscreteAdjoint_;
};
 
} // end of namespace panzer
 
#endif // __Panzer_BlockedEpetraLinearObjFactory_hpp__