Panzer_BCStrategy_Dirichlet_DefaultImpl_impl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//           Panzer: A partial differential equation assembly
//       engine for strongly coupled complex multiphysics systems
//
// Copyright 2011 NTESS and the Panzer contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef PANZER_BCSTRATEGY_DIRICHLET_DEFAULT_IMPL_IMPL_HPP
#define PANZER_BCSTRATEGY_DIRICHLET_DEFAULT_IMPL_IMPL_HPP

#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_Assert.hpp"

#include "Phalanx_DataLayout_MDALayout.hpp"
#include "Phalanx_FieldManager.hpp"

#include "Phalanx_MDField.hpp"
#include "Phalanx_DataLayout.hpp"
#include "Phalanx_DataLayout_MDALayout.hpp"

// Evaluators
#include "Panzer_PhysicsBlock.hpp"
#include "Panzer_PureBasis.hpp"
#include "Panzer_Dirichlet_Residual.hpp"
#include "Panzer_Dirichlet_Residual_EdgeBasis.hpp"
#include "Panzer_Dirichlet_Residual_FaceBasis.hpp"

#ifdef PANZER_HAVE_EPETRA_STACK
#include "Panzer_GatherSolution_Epetra.hpp"
#include "Panzer_ScatterDirichletResidual_Epetra.hpp"
#endif

#include "Panzer_GatherBasisCoordinates.hpp"
#include "Panzer_BasisValues_Evaluator.hpp"
#include "Panzer_PointValues_Evaluator.hpp"
#include "Panzer_DOF.hpp"
#include "Panzer_DOF_PointValues.hpp"

// ***********************************************************************
template <typename EvalT>
panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
BCStrategy_Dirichlet_DefaultImpl(const panzer::BC& bc,
         const Teuchos::RCP<panzer::GlobalData>& global_data,
         const bool in_check_apply_bc) :
  panzer::BCStrategy<EvalT>(bc),
  panzer::GlobalDataAcceptorDefaultImpl(global_data),
  check_apply_bc(in_check_apply_bc),
  descriptor_map_built(false)
{

}

// ***********************************************************************
template <typename EvalT>
panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
~BCStrategy_Dirichlet_DefaultImpl()
{

}

// ***********************************************************************
template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
buildAndRegisterGatherScatterEvaluators(PHX::FieldManager<panzer::Traits>& fm,
                      const panzer::PhysicsBlock& pb,
                const panzer::LinearObjFactory<panzer::Traits> & lof,
          const Teuchos::ParameterList& user_data) const
{
  // for deprecated interface support
  buildDescriptorMapFromVectors();

  buildAndRegisterGatherAndOrientationEvaluators(fm,pb,lof,user_data);
  buildAndRegisterScatterEvaluators(fm,pb,lof,user_data);
}

// ***********************************************************************

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
buildAndRegisterScatterEvaluators(PHX::FieldManager<panzer::Traits>& fm,
                                  const panzer::PhysicsBlock& pb,
                const LinearObjFactory<panzer::Traits> & lof,
                const Teuchos::ParameterList& /* user_data */) const
{
  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using std::vector;
  using std::map;
  using std::string;
  using std::pair;

  // for deprecated interface support
  buildDescriptorMapFromVectors();

  // Scatter
  // for (map<string,string>::const_iterator res_to_dof = residual_to_dof_names_map.begin();
  //      res_to_dof != residual_to_dof_names_map.end(); ++res_to_dof) {
  for(DescriptorIterator itr=m_provided_dofs_desc.begin();
      itr!=m_provided_dofs_desc.end(); ++itr) {

    const DOFDescriptor & desc = itr->second;

    // there is no residual to scatter
    if(!desc.residualName.first)
      continue;

    // std::string residualName = res_to_dof->first;
    // std::string dofName = res_to_dof->second;
    std::string dofName      = desc.dofName;
    std::string residualName = desc.residualName.second;

    ParameterList p("Scatter: "+residualName + " to " + dofName);

    // Set name
    string scatter_field_name = "Dummy Scatter: " + this->m_bc.identifier() + residualName;
    p.set("Scatter Name", scatter_field_name);

    // Set basis
    const vector<pair<string,RCP<panzer::PureBasis> > >& dofBasisPair = pb.getProvidedDOFs();
    RCP<panzer::PureBasis> basis;
    for (vector<pair<string,RCP<panzer::PureBasis> > >::const_iterator it =
     dofBasisPair.begin(); it != dofBasisPair.end(); ++it) {
      if (it->first == dofName)
  basis = it->second;
    }

    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(basis), std::runtime_error,
           "Error the name \"" << dofName
           << "\" is not a valid DOF for the boundary condition:\n"
           << this->m_bc << "\n");

    p.set("Basis", basis);

    RCP<vector<string> > residual_names = rcp(new vector<string>);
    residual_names->push_back(residualName);
    p.set("Dependent Names", residual_names);

    RCP<map<string,string> > names_map = rcp(new map<string,string>);
    names_map->insert(std::make_pair(residualName,dofName));
    p.set("Dependent Map", names_map);

    TEUCHOS_TEST_FOR_EXCEPTION(!pb.cellData().isSide(), std::logic_error,
           "Error - physics block is not a side set!");

    p.set<int>("Side Subcell Dimension",
         pb.getBaseCellTopology().getDimension() - 1);
    p.set<int>("Local Side ID", pb.cellData().side());

    p.set("Check Apply BC",check_apply_bc);

    RCP< PHX::Evaluator<panzer::Traits> > op = lof.buildScatterDirichlet<EvalT>(p);

    this->template registerEvaluator<EvalT>(fm, op);

    // Require variables
    {
      using panzer::Dummy;
      PHX::Tag<typename EvalT::ScalarT> tag(scatter_field_name,
              rcp(new PHX::MDALayout<Dummy>(0)));
      fm.template requireField<EvalT>(tag);
    }

  }
}

// ***********************************************************************

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
buildAndRegisterGatherAndOrientationEvaluators(PHX::FieldManager<panzer::Traits>& fm,
                                         const panzer::PhysicsBlock& pb,
                 const LinearObjFactory<panzer::Traits> & lof,
                 const Teuchos::ParameterList& /* user_data */) const
{
  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using std::vector;
  using std::map;
  using std::string;
  using std::pair;

  // for deprecated interface support
  buildDescriptorMapFromVectors();

  // volume cell data object (used for handling vector valued fields)
  panzer::CellData cellData(pb.cellData().numCells(),pb.cellData().getCellTopology());

  // **************************
  // Coordinates for basis functions (no integration points needed)
  // **************************
  {
    const std::map<std::string,Teuchos::RCP<panzer::PureBasis> > & bases = pb.getBases();
    for (std::map<std::string,Teuchos::RCP<panzer::PureBasis> >::const_iterator it=bases.begin();
         it!=bases.end();it++) {

       Teuchos::RCP<panzer::PureBasis> basis = it->second;

       // add basis coordinates no matter what
       {
         RCP< PHX::Evaluator<panzer::Traits> > basis_op
            = rcp(new panzer::GatherBasisCoordinates<EvalT,panzer::Traits>(*it->second));
         this->template registerEvaluator<EvalT>(fm, basis_op);
       }

       // add point values and basis values
       if(basis->isVectorBasis()) {
         RCP<const panzer::PointRule> pointRule = rcp(new panzer::PointRule(basis->name()+":BasisPoints",basis->cardinality(),cellData));

         {
           RCP< PHX::Evaluator<panzer::Traits> > eval
             = rcp(new panzer::PointValues_Evaluator<EvalT,panzer::Traits>(pointRule,basis));

           this->template registerEvaluator<EvalT>(fm, eval);
         }
         {
           // note basis values are not constructed!
           RCP< PHX::Evaluator<panzer::Traits> > eval
             = rcp(new panzer::BasisValues_Evaluator<EvalT,panzer::Traits>(pointRule,basis,false));

           this->template registerEvaluator<EvalT>(fm, eval);
         }
       }
    }
  }

  // Gather
  for(DescriptorIterator itr=m_provided_dofs_desc.begin();
      itr!=m_provided_dofs_desc.end(); ++itr) {

    // get the dofName from the descriptor
    std::string dofName = itr->second.dofName;
    std::string fieldDof = !itr->second.timeDerivative.first
                              ? itr->second.dofName : itr->second.timeDerivative.second;

    const vector<pair<string,RCP<panzer::PureBasis> > >& dofBasisPair = pb.getProvidedDOFs();
    RCP<panzer::PureBasis> basis;
    for (vector<pair<string,RCP<panzer::PureBasis> > >::const_iterator it =
     dofBasisPair.begin(); it != dofBasisPair.end(); ++it) {
      if (it->first == dofName)
  basis = it->second;
    }

    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(basis), std::runtime_error,
           "Error the name \"" << dofName
           << "\" is not a valid DOF for the boundary condition:\n"
           << this->m_bc << "\n");

    {
      ParameterList p("BC Gather");
      RCP<vector<string> > gather_field_names_vec = rcp(new vector<string>);
      RCP<vector<string> > gather_names_vec = rcp(new vector<string>);
      gather_field_names_vec->push_back(fieldDof);
      gather_names_vec->push_back(dofName);

      p.set("DOF Names", gather_field_names_vec);
      p.set("Indexer Names", gather_names_vec);
      p.set("Basis", basis);
      p.set("Use Time Derivative Solution Vector",itr->second.timeDerivative.first);

      RCP< PHX::Evaluator<panzer::Traits> > op = lof.buildGather<EvalT>(p);
      this->template registerEvaluator<EvalT>(fm, op);
    }

    if(basis->requiresOrientations())  {
      ParameterList p("Gather Orientation");
      RCP<vector<string> > gather_field_names_vec = rcp(new vector<string>);
      RCP<vector<string> > gather_names_vec = rcp(new vector<string>);
      gather_field_names_vec->push_back(fieldDof);
      gather_names_vec->push_back(dofName);

      p.set("DOF Names", gather_field_names_vec);
      p.set("Indexer Names", gather_names_vec);
      p.set("Basis", basis);

      RCP< PHX::Evaluator<panzer::Traits> > op = lof.buildGatherOrientation<EvalT>(p);

      this->template registerEvaluator<EvalT>(fm, op);
    }

    // evaluator a vector basis at the basis points
    if(basis->isVectorBasis()) {
      RCP<const panzer::PointRule> pointRule = rcp(new panzer::PointRule(basis->name()+":BasisPoints",basis->cardinality(),cellData));

      ParameterList p;
      p.set("Name",fieldDof);
      p.set("Basis",basis.getConst());
      p.set("Point Rule",pointRule);

      RCP< PHX::Evaluator<panzer::Traits> > eval
             = rcp(new panzer::DOF_PointValues<EvalT,panzer::Traits>(p));
      this->template registerEvaluator<EvalT>(fm, eval);
    }

  }

  // Dirichlet Residual: residual = dof_value - target_value
  for(DescriptorIterator itr=m_provided_dofs_desc.begin();
      itr!=m_provided_dofs_desc.end(); ++itr) {

    const DOFDescriptor & desc = itr->second;

    // there is no residual to scatter
    if(!desc.residualName.first)
      continue;

    // std::string dofName = desc.dofName;
    std::string dofName = !itr->second.timeDerivative.first
                              ? itr->second.dofName : itr->second.timeDerivative.second;
    std::string residualName = desc.residualName.second;
    std::string targetName = desc.targetName.second;

    const vector<pair<string,RCP<panzer::PureBasis> > >& dofBasisPair = pb.getProvidedDOFs();
    RCP<panzer::PureBasis> basis;
    for (vector<pair<string,RCP<panzer::PureBasis> > >::const_iterator it =
     dofBasisPair.begin(); it != dofBasisPair.end(); ++it) {
      if (it->first == itr->second.dofName)
  basis = it->second;
    }

    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(basis), std::runtime_error,
           "Error the name \"" << itr->second.dofName
           << "\" is not a valid DOF for the boundary condition:\n"
           << this->m_bc << "\n");

    if(basis->isScalarBasis() || desc.coefficientResidual) {
      ParameterList p("Dirichlet Residual: "+residualName + " to " + dofName);
      p.set("Residual Name", residualName);
      p.set("DOF Name", dofName);
      p.set("Value Name", targetName);
      p.set("Data Layout", basis->functional);

      RCP< PHX::Evaluator<panzer::Traits> > op =
        rcp(new panzer::DirichletResidual<EvalT,panzer::Traits>(p));

      this->template registerEvaluator<EvalT>(fm, op);
    }
    // This assumes that dofs on faces are all named "<dof>_face"
    else if(basis->isVectorBasis()&&basis->supportsDiv()) {
      RCP<const panzer::PointRule> pointRule = rcp(new panzer::PointRule(basis->name()+":BasisPoints",basis->cardinality(),cellData));

      ParameterList p;
      p.set("Residual Name", residualName);
      p.set("DOF Name", dofName);
      p.set("Value Name", targetName);
      p.set("Basis", basis.getConst());
      p.set("Point Rule", pointRule);

      RCP< PHX::Evaluator<panzer::Traits> > op =
        rcp(new panzer::DirichletResidual_FaceBasis<EvalT,panzer::Traits>(p));

      this->template registerEvaluator<EvalT>(fm, op);
    }
    else if(basis->isVectorBasis()) {
      RCP<const panzer::PointRule> pointRule = rcp(new panzer::PointRule(basis->name()+":BasisPoints",basis->cardinality(),cellData));

      ParameterList p;
      p.set("Residual Name", residualName);
      p.set("DOF Name", dofName);
      p.set("Value Name", targetName);
      p.set("Basis", basis.getConst());
      p.set("Point Rule", pointRule);

      RCP< PHX::Evaluator<panzer::Traits> > op =
        rcp(new panzer::DirichletResidual_EdgeBasis<EvalT,panzer::Traits>(p));

      this->template registerEvaluator<EvalT>(fm, op);
    }

  }
}

// ***********************************************************************

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
buildDescriptorMapFromVectors() const
{
  if(descriptor_map_built)
    return;

  // build the reverse lookup (this assumes that the map is invertible, though it should be!)
  std::map<std::string,std::string> dof_names_to_residual_map;
  for(std::map<std::string,std::string>::const_iterator itr=residual_to_dof_names_map.begin();
      itr!=residual_to_dof_names_map.end();++itr) {
    dof_names_to_residual_map[itr->second] = itr->first;
  }

  for(std::size_t i=0;i<required_dof_names.size();i++) {
    std::string dof_name      = required_dof_names[i];

    // add the DOF right away
    const_cast<BCStrategy_Dirichlet_DefaultImpl<EvalT> *>(this)->addDOF(dof_name);

    // add target information if its required
    if(dof_names_to_residual_map.find(dof_name)!=dof_names_to_residual_map.end()) {
      std::string residual_name = dof_names_to_residual_map[dof_name];
      std::string target_name   = residual_to_target_field_map.find(residual_name)->second;

      // add the descriptor from the data structures: Note the const_cast is neccessary
      // only because this is protecting deprecated code
      const_cast<BCStrategy_Dirichlet_DefaultImpl<EvalT> *>(this)->
      addTarget(target_name,
                dof_name,
                residual_name);
    }
  }

  descriptor_map_built = true;
}

// ***********************************************************************

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
addDOF(const std::string & dofName)
{
  DescriptorIterator itr = m_provided_dofs_desc.find(dofName);
  TEUCHOS_ASSERT(itr==m_provided_dofs_desc.end());

  DOFDescriptor & desc = m_provided_dofs_desc[dofName];
  desc.dofName = dofName;
}

// ***********************************************************************

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
addTarget(const std::string & targetName,
          const std::string & dofName,
          const std::string & residualName)
{
  typename std::map<std::string,DOFDescriptor>::iterator itr = m_provided_dofs_desc.find(dofName);
  TEUCHOS_ASSERT(itr!=m_provided_dofs_desc.end());

  DOFDescriptor & desc = itr->second;
  desc.dofName = dofName;
  desc.coefficientResidual = false;
  desc.targetName = std::make_pair(true,targetName);
  desc.residualName = (residualName=="") ? std::make_pair(true,"RESIDUAL_"+dofName)
                                         : std::make_pair(true,residualName);
}

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
addCoefficientTarget(const std::string & targetName,
                     const std::string & dofName,
                     const std::string & residualName)
{
  typename std::map<std::string,DOFDescriptor>::iterator itr = m_provided_dofs_desc.find(dofName);
  TEUCHOS_ASSERT(itr!=m_provided_dofs_desc.end());

  DOFDescriptor & desc = itr->second;
  desc.dofName = dofName;
  desc.coefficientResidual = true;
  desc.targetName = std::make_pair(true,targetName);
  desc.residualName = (residualName=="") ? std::make_pair(true,"RESIDUAL_"+dofName)
                                         : std::make_pair(true,residualName);
}

// ***********************************************************************

template <typename EvalT>
void panzer::BCStrategy_Dirichlet_DefaultImpl<EvalT>::
addDotTarget(const std::string & targetName,
             const std::string & dofName,
             const std::string & dotName,
             const std::string & residualName)
{
  typename std::map<std::string,DOFDescriptor>::iterator itr = m_provided_dofs_desc.find(dofName);
  TEUCHOS_ASSERT(itr!=m_provided_dofs_desc.end());

  DOFDescriptor & desc = itr->second;
  desc.dofName = dofName;
  desc.targetName = std::make_pair(true,targetName);
  desc.timeDerivative = (dotName=="") ? std::make_pair(true,"DXDT_"+dofName)
                                      : std::make_pair(true,dotName);
  desc.residualName = (residualName=="") ? std::make_pair(true,"RESIDUAL_"+dofName)
                                         : std::make_pair(true,residualName);
}

// ***********************************************************************

#endif