Panzer_IntegrationValues2.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_INTEGRATION_VALUES2_HPP
#define PANZER_INTEGRATION_VALUES2_HPP

#include "Teuchos_RCP.hpp"

#include "PanzerDiscFE_config.hpp"
#include "Panzer_IntegrationRule.hpp"
#include "Panzer_ArrayTraits.hpp"
#include "Panzer_Dimension.hpp"
#include "Phalanx_MDField.hpp"
#include "Intrepid2_Cubature.hpp"

namespace panzer {

  class SubcellConnectivity;

  template <typename Scalar>
  class IntegrationValues2 {
  public:
    typedef typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type size_type;

    typedef PHX::MDField<Scalar> ArrayDynamic;
    typedef PHX::MDField<double> DblArrayDynamic;

    typedef PHX::MDField<Scalar,IP>               Array_IP;
    typedef PHX::MDField<Scalar,IP,Dim>           Array_IPDim;
    typedef PHX::MDField<Scalar,Point>            Array_Point;
    typedef PHX::MDField<Scalar,Cell,IP>          Array_CellIP;
    typedef PHX::MDField<Scalar,Cell,IP,Dim>      Array_CellIPDim;
    typedef PHX::MDField<Scalar,Cell,IP,Dim,Dim>  Array_CellIPDimDim;
    typedef PHX::MDField<Scalar,Cell,BASIS,Dim>   Array_CellBASISDim;

    typedef PHX::MDField<const Scalar,IP>               ConstArray_IP;
    typedef PHX::MDField<const Scalar,IP,Dim>           ConstArray_IPDim;
    typedef PHX::MDField<const Scalar,Point>            ConstArray_Point;
    typedef PHX::MDField<const Scalar,Cell,IP>          ConstArray_CellIP;
    typedef PHX::MDField<const Scalar,Cell,IP,Dim>      ConstArray_CellIPDim;
    typedef PHX::MDField<const Scalar,Cell,IP,Dim,Dim>  ConstArray_CellIPDimDim;
    typedef PHX::MDField<const Scalar,Cell,BASIS,Dim>   ConstArray_CellBASISDim;

    IntegrationValues2(const std::string & pre="",
                       const bool allocArrays=false);


    // =====================================================================================================
    // Classic Interface (DEPRECATED)

    void setupArrays(const Teuchos::RCP<const panzer::IntegrationRule>& ir);

    void evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
                        const int num_cells = -1,
                        const Teuchos::RCP<const SubcellConnectivity> & face_connectivity = Teuchos::null,
                        const int num_virtual_cells = -1);

    void evaluateValues(const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
                        const PHX::MDField<Scalar,Cell,IP,Dim> & other_ip_coordinates,
                        const int num_cells = -1);

    // Reference space quantities
    mutable Array_IPDim cub_points;              // <IP,Dim>
    mutable Array_IPDim side_cub_points;         // <IP,Dim> points on face topology (dim-1)
    mutable Array_IP cub_weights;                // <IP>

    // Physical space quantities
    mutable Array_CellBASISDim node_coordinates; // <Cell,BASIS,Dim>
    mutable Array_CellIPDimDim jac;              // <Cell,IP,Dim,Dim>
    mutable Array_CellIPDimDim jac_inv;          // <Cell,IP,Dim,Dim>
    mutable Array_CellIP jac_det;                // <Cell,IP>
    mutable Array_CellIP weighted_measure;       // <Cell,IP>
    mutable Array_CellIPDim weighted_normals;    // <Cell,IP,Dim>
    mutable Array_CellIPDim surface_normals;    // <Cell,IP,Dim>
    mutable Array_CellIPDimDim surface_rotation_matrices;    // <Cell,IP,Dim,Dim>
      // this (appears) is a matrix where the first row is the "normal" direction
      // and the remaining two rows lie in the hyperplane

    // for Shakib stabilization <Cell,IP,Dim,Dim>
    mutable Array_CellIPDimDim covarient;
    mutable Array_CellIPDimDim contravarient;
    mutable Array_CellIP norm_contravarient;

    // integration points
    mutable Array_CellIPDim ip_coordinates;      // <Cell,IP,Dim>
    mutable Array_CellIPDim ref_ip_coordinates;  // <Cell,IP,Dim> for Control Volumes or Surface integrals


    Teuchos::RCP<const panzer::IntegrationRule> int_rule;

    Teuchos::RCP<Intrepid2::Cubature<PHX::Device::execution_space,double,double>> intrepid_cubature;

    // =====================================================================================================
    // Lazy evaluation interface

    void
    setup(const Teuchos::RCP<const panzer::IntegrationRule>& ir,
          const PHX::MDField<Scalar,Cell,NODE,Dim> & node_coordinates,
          const int num_cells = -1);

    void
    setupPermutations(const Teuchos::RCP<const SubcellConnectivity> & face_connectivity,
                      const int num_virtual_cells);

    void
    setupPermutations(const PHX::MDField<Scalar,Cell,IP,Dim> & other_ip_coordinates);

    ConstArray_IPDim
    getUniformCubaturePointsRef(const bool cache = true,
                                const bool force = false,
                                const bool apply_permutation = true) const;

    ConstArray_IPDim
    getUniformSideCubaturePointsRef(const bool cache = true,
                                    const bool force = false,
                                    const bool apply_permutation = true) const;

    ConstArray_IP
    getUniformCubatureWeightsRef(const bool cache = true,
                                 const bool force = false,
                                 const bool apply_permutation = true) const;


    ConstArray_CellBASISDim
    getNodeCoordinates() const;

    ConstArray_CellIPDimDim
    getJacobian(const bool cache = true,
                const bool force = false) const;


    ConstArray_CellIPDimDim
    getJacobianInverse(const bool cache = true,
                       const bool force = false) const;

    ConstArray_CellIP
    getJacobianDeterminant(const bool cache = true,
                           const bool force = false) const;

    ConstArray_CellIP
    getWeightedMeasure(const bool cache = true,
                       const bool force = false) const;

    ConstArray_CellIPDim
    getWeightedNormals(const bool cache = true,
                       const bool force = false) const;

    ConstArray_CellIPDim
    getSurfaceNormals(const bool cache = true,
                      const bool force = false) const;

    ConstArray_CellIPDimDim
    getSurfaceRotationMatrices(const bool cache = true,
                               const bool force = false) const;

    ConstArray_CellIPDimDim
    getCovarientMatrix(const bool cache = true,
                       const bool force = false) const;

    ConstArray_CellIPDimDim
    getContravarientMatrix(const bool cache = true,
                           const bool force = false) const;

    ConstArray_CellIP
    getNormContravarientMatrix(const bool cache = true,
                               const bool force = false) const;

    ConstArray_CellIPDim
    getCubaturePoints(const bool cache = true,
                      const bool force = false) const;

    ConstArray_CellIPDim
    getCubaturePointsRef(const bool cache = true,
                         const bool force = false) const;

    // =====================================================================================================

  protected:

    // Reset all the lazy evaluation arrays
    void
    resetArrays();

    // Number of cells in mesh
    int num_cells_;

    // Number of cells in mesh to evaluate
    int num_evaluate_cells_;

    // Number of virtual cells in the mesh - used for surface evaluations
    int num_virtual_cells_;

    // Permutations (used to re-orient arrays similar to orientations in BasisValues2)
    bool requires_permutation_;

    // Array contains the mapping from uniform reference space to permuted space
    PHX::MDField<const int,Cell,IP> permutations_;

    // TODO: There is a way around this, but it will require some work
    // Subcell connectivity is required for surface evaluations (normals and rotation matrices)
    Teuchos::RCP<const SubcellConnectivity> side_connectivity_;

    // Lazy evaluation checks
    mutable bool cub_points_evaluated_;
    mutable bool side_cub_points_evaluated_;
    mutable bool cub_weights_evaluated_;
    mutable bool node_coordinates_evaluated_;
    mutable bool jac_evaluated_;
    mutable bool jac_inv_evaluated_;
    mutable bool jac_det_evaluated_;
    mutable bool weighted_measure_evaluated_;
    mutable bool weighted_normals_evaluated_;
    mutable bool surface_normals_evaluated_;
    mutable bool surface_rotation_matrices_evaluated_;
    mutable bool covarient_evaluated_;
    mutable bool contravarient_evaluated_;
    mutable bool norm_contravarient_evaluated_;
    mutable bool ip_coordinates_evaluated_;
    mutable bool ref_ip_coordinates_evaluated_;

    // Backward compatibility call that evaluates all internal values for CV, surface, side, or volume integration schemes
    void
    evaluateEverything();

  private:

    bool alloc_arrays_;
    std::string prefix_;
    std::vector<PHX::index_size_type> ddims_;

  };

} // namespace panzer

#endif