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) 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

#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()");
}

}