Panzer_GeometricAggFieldPattern.cpp

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
 
#include "Panzer_GeometricAggFieldPattern.hpp"
 
using Teuchos::RCP;
 
namespace panzer {
 
GeometricAggFieldPattern::GeometricAggFieldPattern()
   : patternBuilt_(false), dimension_(0)
{}
 
  GeometricAggFieldPattern::GeometricAggFieldPattern(std::vector<std::pair<FieldType,Teuchos::RCP<const FieldPattern>>> & patterns) 
   : patternBuilt_(false), dimension_(0)
{ 
   buildPattern(patterns); 
}
 
GeometricAggFieldPattern::GeometricAggFieldPattern(const FieldType& fieldType,
                                                   const Teuchos::RCP<const FieldPattern> & pattern) 
   : patternBuilt_(false), dimension_(0)
{ 
  buildPattern(fieldType,pattern);
}
 
void GeometricAggFieldPattern::buildPattern(const std::vector<std::pair<FieldType,Teuchos::RCP<const FieldPattern>>> & patterns)
{
   std::size_t numPat = patterns.size();
 
   // must be at least one field to do something
   if(numPat<1) {
      bool no_patterns_to_construct = true;
      TEUCHOS_TEST_FOR_EXCEPTION(no_patterns_to_construct,std::logic_error,
                         "GeometricAggFieldPattern::buildPattern requires at least one field pattern");
      return;
   }
 
   bool sameGeometry=true;
   for(std::size_t i=1;i<patterns.size();i++)
     sameGeometry &= patterns[0].second->sameGeometry(*(patterns[i].second));       
   TEUCHOS_TEST_FOR_EXCEPTION(not sameGeometry,std::logic_error,
             "GeometricAggFieldPattern::buildPattern(): Patterns must "
             "have the same geometry!");
 
   // copy cell topology
   cellTopo_ = patterns[0].second->getCellTopology();
 
   // grab the dimension
   dimension_ = patterns[0].second->getDimension();
   patternData_.resize(dimension_+1);
 
   // build space for subcells
   std::vector<int> subcellCount(dimension_+1);
   for(std::size_t d=0;d<dimension_+1;d++) {
      subcellCount[d] = patterns[0].second->getSubcellCount(d);
      patternData_[d].resize(subcellCount[d]);
   }
 
   // Build geometric pattern: increment it logically over all the
   // subcells. Start with CG fields first. Then all DG fields. This
   // is done so that we can use individual field pattern offsets when
   // mapping DOFs to subcells.
   int counter = 0;
   for(std::size_t d=0;d<dimension_+1;d++) {
      for(int s=0;s<subcellCount[d];s++) {
         std::vector<int> & current = patternData_[d][s];
         for(std::size_t p=0;p<patterns.size();p++) {
           if ( (patterns[p].first) == FieldType::CG) {
              RCP<const FieldPattern> field = (patterns[p]).second;
              // if dofs exist, we have a geometric entity 
              const std::size_t num = ( (field->getSubcellIndices(d,s).size() > 0) ? 1 : 0 );
              if(current.size()<num) { 
                for(int i=num-current.size();i>0;i--,counter++) 
                  current.push_back(counter);
              }
            }
         } 
      }
   }
 
   // Add DG fields. These fields are considered internal "cell"
   // fields for DOF numbering.
   const int cellDim = dimension_;
   for(std::size_t d=0;d<dimension_+1;d++) {
      for(int s=0;s<subcellCount[d];s++) {
         std::vector<int> & current = patternData_[cellDim][0];
         for(std::size_t p=0;p<patterns.size();p++) {
           if ( (patterns[p].first) == FieldType::DG) {
              RCP<const FieldPattern> field = (patterns[p]).second;
              // if dofs exist, we have a geometric entity
              const std::size_t num = ( (field->getSubcellIndices(d,s).size() > 0) ? 1 : 0 );
              if(current.size()<num) { 
                for(int i=num-current.size();i>0;i--,counter++) 
                  current.push_back(counter);
              }
            }
         } 
      }
   }
   
   // record that the pattern has been built
   patternBuilt_ = true;
}
 
void GeometricAggFieldPattern::buildPattern(const FieldType& fieldType,
                                            const Teuchos::RCP<const FieldPattern> & pattern)
{
  std::vector<std::pair<FieldType,Teuchos::RCP<const FieldPattern>>> patterns;
  patterns.push_back(std::make_pair(fieldType,pattern));
  buildPattern(patterns);
}
 
int GeometricAggFieldPattern::getSubcellCount(int dim) const
{
   if(patternBuilt_) return patternData_[dim].size();
 
   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
            "GeometricAggFieldPattern::getSubcellCount() cannot be called before "
            "GeometricAggFieldPattern::buildPattern()");
}
 
const std::vector<int> & GeometricAggFieldPattern::getSubcellIndices(int dim,int cellIndex) const
{
   if(patternBuilt_) return patternData_[dim][cellIndex];
 
   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
            "GeometricAggFieldPattern::getSubcellIndices() cannot be called before "
            "GeometricAggFieldPattern::buildPattern()");
}
 
int GeometricAggFieldPattern::getDimension() const
{
   if(patternBuilt_) return dimension_;
 
   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
            "GeometricAggFieldPattern::getDimension() cannot be called before "
            "GeometricAggFieldPattern::buildPattern()");
}
 
shards::CellTopology GeometricAggFieldPattern::getCellTopology() const
{
   if(patternBuilt_) return cellTopo_;
 
   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
            "GeometricAggFieldPattern::getCellTopology() cannot be called before "
            "GeometricAggFieldPattern::buildPattern()");
}
 
}