Panzer_BasisValues2.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_BASIS_VALUES2_HPP
#define PANZER_BASIS_VALUES2_HPP

#include "Teuchos_RCP.hpp"

#include "Intrepid2_Basis.hpp"
#include "Intrepid2_Orientation.hpp"

#include "Panzer_BasisIRLayout.hpp"
#include "Panzer_ArrayTraits.hpp"

namespace panzer {

  // For applying orientations with lazy evaluation
  class OrientationsInterface;

  template <typename Scalar>
  class BasisValues2 {
  public:
    typedef typename ArrayTraits<Scalar,PHX::MDField<Scalar> >::size_type size_type;
    using IntrepidBasis = Intrepid2::Basis<PHX::Device::execution_space,Scalar,Scalar>;

    typedef PHX::MDField<Scalar>                    ArrayDynamic;
    typedef PHX::MDField<Scalar,BASIS,IP>           Array_BasisIP;
    typedef PHX::MDField<Scalar,BASIS,IP,Dim>       Array_BasisIPDim;
    typedef PHX::MDField<Scalar,BASIS,Dim>          Array_BasisDim;
    typedef PHX::MDField<Scalar,Cell,BASIS,IP>      Array_CellBasisIP;
    typedef PHX::MDField<Scalar,Cell,BASIS,IP,Dim>  Array_CellBasisIPDim;
    typedef PHX::MDField<Scalar,Cell,BASIS,Dim>     Array_CellBasisDim;

    typedef PHX::MDField<const Scalar>                    ConstArrayDynamic;
    typedef PHX::MDField<const Scalar,BASIS,IP>           ConstArray_BasisIP;
    typedef PHX::MDField<const Scalar,BASIS,IP,Dim>       ConstArray_BasisIPDim;
    typedef PHX::MDField<const Scalar,BASIS,Dim>          ConstArray_BasisDim;
    typedef PHX::MDField<const Scalar,Cell,BASIS,IP>      ConstArray_CellBasisIP;
    typedef PHX::MDField<const Scalar,Cell,BASIS,IP,Dim>  ConstArray_CellBasisIPDim;
    typedef PHX::MDField<const Scalar,Cell,BASIS,Dim>     ConstArray_CellBasisDim;

    //=======================================================================================================
    // DEPRECATED Interface

    BasisValues2(const std::string & prefix="",
                 const bool allocArrays=false,
                 const bool buildWeighted=false);

    void setupArrays(const Teuchos::RCP<const panzer::BasisIRLayout>& basis,
                     bool computeDerivatives=true);

    void evaluateValues(const PHX::MDField<Scalar,IP,Dim> & cub_points,
                        const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
                        const PHX::MDField<Scalar,Cell,IP> & jac_det,
                        const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
                        const int in_num_cells = -1);

    void evaluateValues(const PHX::MDField<Scalar,IP,Dim> & cub_points,
                        const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
                        const PHX::MDField<Scalar,Cell,IP> & jac_det,
                        const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
                        const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
                        const PHX::MDField<Scalar,Cell,NODE,Dim> & node_coordinates,
                        bool use_node_coordinates=true,
                        const int in_num_cells = -1);

    void evaluateValuesCV(const PHX::MDField<Scalar,Cell,IP,Dim> & cell_cub_points,
                          const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
                          const PHX::MDField<Scalar,Cell,IP> & jac_det,
                          const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv);


    void evaluateValuesCV(const PHX::MDField<Scalar,Cell,IP,Dim> & cell_cub_points,
                          const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
                          const PHX::MDField<Scalar,Cell,IP> & jac_det,
                          const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
                          const PHX::MDField<Scalar,Cell,NODE,Dim> & node_coordinates,
                          bool use_node_coordinates=true,
                          const int in_num_cells = -1);


    void evaluateValues(const PHX::MDField<Scalar,Cell,IP,Dim> & cub_points,
                        const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac,
                        const PHX::MDField<Scalar,Cell,IP> & jac_det,
                        const PHX::MDField<Scalar,Cell,IP,Dim,Dim> & jac_inv,
                        const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
                        const PHX::MDField<Scalar,Cell,NODE,Dim> & node_coordinates,
                        bool use_node_coordinates=true,
                        const int in_num_cells = -1);


    // some evaluators use this apply orientation (this will be deprecated)
    void applyOrientations(const PHX::MDField<const Scalar,Cell,BASIS> & orientations);

    // this is used in workset factory
    void applyOrientations(const std::vector<Intrepid2::Orientation> & orientations,
                           const int in_num_cells = -1);

    void setExtendedDimensions(const std::vector<PHX::index_size_type> & ddims)
    { ddims_ = ddims; }

    PureBasis::EElementSpace getElementSpace() const;

    mutable Array_BasisIP     basis_ref_scalar;           // <BASIS,IP>
    mutable Array_CellBasisIP basis_scalar;               // <Cell,BASIS,IP>

    mutable Array_BasisIPDim     basis_ref_vector;           // <BASIS,IP,Dim>
    mutable Array_CellBasisIPDim basis_vector;               // <Cell,BASIS,IP,Dim>

    mutable Array_BasisIPDim     grad_basis_ref;             // <BASIS,IP,Dim>
    mutable Array_CellBasisIPDim grad_basis;                 // <Cell,BASIS,IP,Dim>

    mutable Array_BasisIP     curl_basis_ref_scalar;         // <BASIS,IP> - 2D!
    mutable Array_CellBasisIP curl_basis_scalar;             // <Cell,BASIS,IP> - 2D!

    mutable Array_BasisIPDim     curl_basis_ref_vector;      // <BASIS,IP,Dim> - 3D!
    mutable Array_CellBasisIPDim curl_basis_vector;          // <Cell,BASIS,IP,Dim> - 3D!

    mutable Array_BasisIP     div_basis_ref;           // <BASIS,IP>
    mutable Array_CellBasisIP div_basis;               // <Cell,BASIS,IP>

    mutable Array_CellBasisIP weighted_basis_scalar;                  // <Cell,BASIS,IP>
    mutable Array_CellBasisIPDim weighted_basis_vector;               // <Cell,BASIS,IP,Dim>
    mutable Array_CellBasisIPDim weighted_grad_basis;                 // <Cell,BASIS,IP,Dim>
    mutable Array_CellBasisIP weighted_curl_basis_scalar;             // <Cell,BASIS,IP>
    mutable Array_CellBasisIPDim weighted_curl_basis_vector;          // <Cell,BASIS,IP,Dim>
    mutable Array_CellBasisIP weighted_div_basis;                     // <Cell,BASIS,IP>

    mutable Array_BasisDim basis_coordinates_ref;     // <BASIS,Dim>

    mutable Array_CellBasisDim basis_coordinates;     // <Cell,BASIS,Dim>

    Teuchos::RCP<const panzer::BasisIRLayout> basis_layout;

    Teuchos::RCP<IntrepidBasis> intrepid_basis;

    bool compute_derivatives;
    bool build_weighted;
    bool alloc_arrays;
    std::string prefix;
    std::vector<PHX::index_size_type> ddims_;

    bool orientationsApplied() const
    {return orientations_applied_;}

    void evaluateBasisCoordinates(const PHX::MDField<Scalar,Cell,NODE,Dim> & node_coordinates,
                                  const int in_num_cells = -1);

  private:

    bool references_evaluated;

    // Have orientations been applied
    bool orientations_applied_;

    //=======================================================================================================
    // New Interface

  protected:

    // Reset internal data structures for lazy evaluation
    void resetArrays();

    // Number of total cells (owned + ghost + virtual)
    int num_cells_;

    // Number of cells to evaluate (HACK for backward compatability)
    int num_evaluate_cells_;

    // True if a uniform point space is used
    bool is_uniform_;

    // Only valid for uniform reference space
    PHX::MDField<const Scalar,IP,Dim>           cubature_points_uniform_ref_;

    // For non-uniform reference space
    PHX::MDField<const Scalar,Cell,IP,Dim>      cubature_points_ref_;

    // Geometry objects that represent cubature/point values
    PHX::MDField<const Scalar,Cell,IP,Dim,Dim>  cubature_jacobian_;
    PHX::MDField<const Scalar,Cell,IP>          cubature_jacobian_determinant_;
    PHX::MDField<const Scalar,Cell,IP,Dim,Dim>  cubature_jacobian_inverse_;
    PHX::MDField<const Scalar,Cell,IP>          cubature_weights_;

    // Coordinates of the mesh nodes
    PHX::MDField<const Scalar,Cell,NODE,Dim> cell_node_coordinates_;

    // Cell topology from the mesh
    Teuchos::RCP<const shards::CellTopology> cell_topology_;

    // Number of cells to apply orientations to (required in situations where virtual cells exist)
    int num_orientations_cells_;

    // Orientations object
    std::vector<Intrepid2::Orientation>   orientations_;

    mutable bool basis_ref_scalar_evaluated_;
    mutable bool basis_scalar_evaluated_;
    mutable bool basis_ref_vector_evaluated_;
    mutable bool basis_vector_evaluated_;
    mutable bool grad_basis_ref_evaluated_;
    mutable bool grad_basis_evaluated_;
    mutable bool curl_basis_ref_scalar_evaluated_;
    mutable bool curl_basis_scalar_evaluated_;
    mutable bool curl_basis_ref_vector_evaluated_;
    mutable bool curl_basis_vector_evaluated_;
    mutable bool div_basis_ref_evaluated_;
    mutable bool div_basis_evaluated_;
    mutable bool weighted_basis_scalar_evaluated_;
    mutable bool weighted_basis_vector_evaluated_;
    mutable bool weighted_grad_basis_evaluated_;
    mutable bool weighted_curl_basis_scalar_evaluated_;
    mutable bool weighted_curl_basis_vector_evaluated_;
    mutable bool weighted_div_basis_evaluated_;
    mutable bool basis_coordinates_ref_evaluated_;
    mutable bool basis_coordinates_evaluated_;

  public:

    void
    setup(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
          PHX::MDField<const Scalar, Cell, IP, Dim>         reference_points,
          PHX::MDField<const Scalar, Cell, IP, Dim, Dim>    point_jacobian,
          PHX::MDField<const Scalar, Cell, IP>              point_jacobian_determinant,
          PHX::MDField<const Scalar, Cell, IP, Dim, Dim>    point_jacobian_inverse,
          const int                                         num_evaluated_cells = -1);

    void
    setupUniform(const Teuchos::RCP<const panzer::BasisIRLayout> & basis,
                 PHX::MDField<const Scalar, IP, Dim>               reference_points,
                 PHX::MDField<const Scalar, Cell, IP, Dim, Dim>    point_jacobian,
                 PHX::MDField<const Scalar, Cell, IP>              point_jacobian_determinant,
                 PHX::MDField<const Scalar, Cell, IP, Dim, Dim>    point_jacobian_inverse,
                 const int                                         num_evaluated_cells = -1);

    void
    setOrientations(const std::vector<Intrepid2::Orientation> & orientations,
                    const int num_orientations_cells = -1);

    void
    setWeightedMeasure(PHX::MDField<const Scalar, Cell, IP> weighted_measure);

    void
    setCellVertexCoordinates(PHX::MDField<Scalar,Cell,NODE,Dim> vertex_coordinates);

    void
    setCellNodeCoordinates(PHX::MDField<Scalar,Cell,NODE,Dim> node_coordinates);

    bool
    hasUniformReferenceSpace() const
    {return is_uniform_;}

    panzer::BasisDescriptor getBasisDescriptor() const;

    const std::vector<PHX::index_size_type> &
    getExtendedDimensions() const
    {return ddims_;}

    //=======================================================================================================
    // Reference space evaluations

    ConstArray_BasisDim
    getBasisCoordinatesRef(const bool cache = true,
                           const bool force = false) const;

    ConstArray_BasisIP
    getBasisValuesRef(const bool cache = true,
                      const bool force = false) const;

    ConstArray_BasisIPDim
    getVectorBasisValuesRef(const bool cache = true,
                            const bool force = false) const;

    ConstArray_BasisIPDim
    getGradBasisValuesRef(const bool cache = true,
                          const bool force = false) const;

    ConstArray_BasisIP
    getCurl2DVectorBasisRef(const bool cache = true,
                            const bool force = false) const;

    ConstArray_BasisIPDim
    getCurlVectorBasisRef(const bool cache = true,
                          const bool force = false) const;

    ConstArray_BasisIP
    getDivVectorBasisRef(const bool cache = true,
                         const bool force = false) const;

    //=======================================================================================================
    // Mesh evaluations

    ConstArray_CellBasisDim
    getBasisCoordinates(const bool cache = true,
                        const bool force = false) const;

    ConstArray_CellBasisIP
    getBasisValues(const bool weighted,
                   const bool cache = true,
                   const bool force = false) const;

    ConstArray_CellBasisIPDim
    getVectorBasisValues(const bool weighted,
                         const bool cache = true,
                         const bool force = false) const;

    ConstArray_CellBasisIPDim
    getGradBasisValues(const bool weighted,
                       const bool cache = true,
                       const bool force = false) const;

    ConstArray_CellBasisIP
    getCurl2DVectorBasis(const bool weighted,
                         const bool cache = true,
                         const bool force = false) const;

    ConstArray_CellBasisIPDim
    getCurlVectorBasis(const bool weighted,
                       const bool cache = true,
                       const bool force = false) const;

    ConstArray_CellBasisIP
    getDivVectorBasis(const bool weighted,
                      const bool cache = true,
                      const bool force = false) const;

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

  };

} // namespace panzer

#endif