Intrepid2_OrientationToolsDefModifyBasis.hpp

Go to the documentation of this file.

// @HEADER
// ************************************************************************
//
//                           Intrepid2 Package
//                 Copyright (2007) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// 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 Kyungjoo Kim  (kyukim@sandia.gov), or
//                    Mauro Perego  (mperego@sandia.gov)
//
// ************************************************************************
// @HEADER


#ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_BASIS_HPP__
#define __INTREPID2_ORIENTATIONTOOLS_DEF_MODIFY_BASIS_HPP__

// disable clang warnings
#if defined (__clang__) && !defined (__INTEL_COMPILER)
#pragma clang system_header
#endif

#include "Intrepid2_Orientation.hpp"

namespace Intrepid2 {

  template<typename DT>
  template<typename elemOrtValueType, class ...elemOrtProperties,
           typename elemNodeValueType, class ...elemNodeProperties>
  void
  OrientationTools<DT>::
  getOrientation(      Kokkos::DynRankView<elemOrtValueType,elemOrtProperties...> elemOrts,
                 const Kokkos::DynRankView<elemNodeValueType,elemNodeProperties...> elemNodes,
                 const shards::CellTopology cellTopo,
                 bool isSide) {
    // small meta data modification and it uses shards; let's do this on host
    auto elemOrtsHost = Kokkos::create_mirror_view(elemOrts);
    auto elemNodesHost = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elemNodes);
    
    const ordinal_type numCells = elemNodes.extent(0);
    for (auto cell=0;cell<numCells;++cell) {
      const auto nodes = Kokkos::subview(elemNodesHost, cell, Kokkos::ALL());
      elemOrtsHost(cell) = Orientation::getOrientation(cellTopo, nodes, isSide);
    }

    Kokkos::deep_copy(elemOrts, elemOrtsHost);
  }

  template<typename ortViewType,
           typename OutputViewType,
           typename inputViewType,
           typename o2tViewType,
           typename t2oViewType,
           typename dataViewType>
  struct F_modifyBasisByOrientation {
    ortViewType orts;
    OutputViewType output;
    inputViewType input;
    o2tViewType ordinalToTag;
    t2oViewType tagToOrdinal;

    const dataViewType matData;
    const ordinal_type cellDim, numVerts, numEdges, numFaces, numPoints, dimBasis;

    F_modifyBasisByOrientation(ortViewType orts_,
                               OutputViewType output_,
                               inputViewType input_,
                               o2tViewType ordinalToTag_,
                               t2oViewType tagToOrdinal_,
                               const dataViewType matData_,
                               const ordinal_type cellDim_,
                               const ordinal_type numVerts_,
                               const ordinal_type numEdges_,
                               const ordinal_type numFaces_,
                               const ordinal_type numPoints_,
                               const ordinal_type dimBasis_)
    : orts(orts_),
      output(output_),
      input(input_),
      ordinalToTag(ordinalToTag_),
      tagToOrdinal(tagToOrdinal_),
      matData(matData_),
      cellDim(cellDim_),
      numVerts(numVerts_),
      numEdges(numEdges_),
      numFaces(numFaces_),
      numPoints(numPoints_),
      dimBasis(dimBasis_)
    {}

    KOKKOS_INLINE_FUNCTION
    void operator()(const ordinal_type cell) const {
      typedef typename inputViewType::non_const_value_type input_value_type;

      auto out = Kokkos::subview(output, cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());
      auto in  = (input.rank() == output.rank()) ?
                 Kokkos::subview(input,  cell, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL())
               : Kokkos::subview(input,        Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL());

      // edge transformation
      ordinal_type existEdgeDofs = 0;
      if (numEdges > 0) {
        ordinal_type ortEdges[12];
        orts(cell).getEdgeOrientation(ortEdges, numEdges);

        // apply coeff matrix
        for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId) {
          const ordinal_type ordEdge = (1 < tagToOrdinal.extent(0) ? (static_cast<size_type>(edgeId) < tagToOrdinal.extent(1) ? tagToOrdinal(1, edgeId, 0) : -1) : -1);

          if (ordEdge != -1) {
            existEdgeDofs = 1;
            const ordinal_type ndofEdge = ordinalToTag(ordEdge, 3);
            const auto mat = Kokkos::subview(matData,
                                             edgeId, ortEdges[edgeId],
                                             Kokkos::ALL(), Kokkos::ALL());

            for (ordinal_type j=0;j<numPoints;++j)
              for (ordinal_type i=0;i<ndofEdge;++i) {
                const ordinal_type ii = tagToOrdinal(1, edgeId, i);

                for (ordinal_type k=0;k<dimBasis;++k) {
                  input_value_type temp = 0.0;
                  for (ordinal_type l=0;l<ndofEdge;++l) {
                    const ordinal_type ll = tagToOrdinal(1, edgeId, l);
                    temp += mat(i,l)*in(ll, j, k);
                  }
                  out(ii, j, k) = temp;
                }
              }
          }
        }
      }

      // face transformation
      if (numFaces > 0) {
        ordinal_type ortFaces[12];
        orts(cell).getFaceOrientation(ortFaces, numFaces);

        // apply coeff matrix
        for (ordinal_type faceId=0;faceId<numFaces;++faceId) {
          const ordinal_type ordFace = (2 < tagToOrdinal.extent(0) ? (static_cast<size_type>(faceId) < tagToOrdinal.extent(1) ? tagToOrdinal(2, faceId, 0) : -1) : -1);

          if (ordFace != -1) {
            const ordinal_type ndofFace = ordinalToTag(ordFace, 3);
            const auto mat = Kokkos::subview(matData,
                                             numEdges*existEdgeDofs+faceId, ortFaces[faceId],
                                             Kokkos::ALL(), Kokkos::ALL());

            for (ordinal_type j=0;j<numPoints;++j)
              for (ordinal_type i=0;i<ndofFace;++i) {
                const ordinal_type ii = tagToOrdinal(2, faceId, i);

                for (ordinal_type k=0;k<dimBasis;++k) {
                  input_value_type temp = 0.0;
                  for (ordinal_type l=0;l<ndofFace;++l) {
                    const ordinal_type ll = tagToOrdinal(2, faceId, l);
                    temp += mat(i,l)*in(ll, j, k);
                  }
                  out(ii, j, k) = temp;
                }
              }
          }
        }
      }

      //side orientations
      ordinal_type faceOrt(0), edgeOrt(0);
      if(cellDim == 2) orts(cell).getFaceOrientation(&faceOrt, 1);
      if (faceOrt != 0) {
        const ordinal_type ordFace = (2 < tagToOrdinal.extent(0) ? (static_cast<size_type>(0) < tagToOrdinal.extent(1) ? tagToOrdinal(2, 0, 0) : -1) : -1);

        if (ordFace != -1) {
          const ordinal_type ndofFace = ordinalToTag(ordFace, 3);
          const auto mat = Kokkos::subview(matData,
                                           numEdges*existEdgeDofs, faceOrt,
                                           Kokkos::ALL(), Kokkos::ALL());

          for (ordinal_type j=0;j<numPoints;++j)
            for (ordinal_type i=0;i<ndofFace;++i) {
              const ordinal_type ii = tagToOrdinal(2, 0, i);

              for (ordinal_type k=0;k<dimBasis;++k) {
                input_value_type temp = 0.0;
                for (ordinal_type l=0;l<ndofFace;++l) {
                  const ordinal_type ll = tagToOrdinal(2, 0, l);
                  temp += mat(i,l)*in(ll, j, k);
                }
                out(ii, j, k) = temp;
              }
            }
        }
      }

      if(cellDim == 1) orts(cell).getEdgeOrientation(&edgeOrt, 1);
      if (edgeOrt != 0) {
        const ordinal_type ordEdge = (1 < tagToOrdinal.extent(0) ? (static_cast<size_type>(0) < tagToOrdinal.extent(1) ? tagToOrdinal(1, 0, 0) : -1) : -1);

        if (ordEdge != -1) {
          const ordinal_type ndofEdge = ordinalToTag(ordEdge, 3);
          const auto mat = Kokkos::subview(matData,
                                           0, edgeOrt,
                                           Kokkos::ALL(), Kokkos::ALL());

          for (ordinal_type j=0;j<numPoints;++j)
            for (ordinal_type i=0;i<ndofEdge;++i) {
              const ordinal_type ii = tagToOrdinal(1, 0, i);

              for (ordinal_type k=0;k<dimBasis;++k) {
                input_value_type temp = 0.0;
                for (ordinal_type l=0;l<ndofEdge;++l) {
                  const ordinal_type ll = tagToOrdinal(1, 0, l);
                  temp += mat(i,l)*in(ll, j, k);
                }
                out(ii, j, k) = temp;
              }
            }
        }
      }
    }
  };

  template<typename DT>
  template<typename outputValueType, class ...outputProperties,
           typename inputValueType,  class ...inputProperties,
           typename OrientationViewType,
           typename BasisType>
  void
  OrientationTools<DT>::
  modifyBasisByOrientation(      Kokkos::DynRankView<outputValueType,outputProperties...> output,
                           const Kokkos::DynRankView<inputValueType, inputProperties...>  input,
                           const OrientationViewType orts,
                           const BasisType* basis ) {
#ifdef HAVE_INTREPID2_DEBUG
    {
      if (input.rank() == output.rank())
      {
        for (size_type i=0;i<input.rank();++i)
          INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i), std::invalid_argument,
                                        ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input and output dimension does not match.");
      }
      else if (input.rank() == output.rank() - 1)
      {
        for (size_type i=0;i<input.rank();++i)
          INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(i+1), std::invalid_argument,
                                       ">>> ERROR (OrientationTools::modifyBasisByOrientation): Input dimensions must match output dimensions exactly, or else match all but the first dimension (in the case that input does not have a 'cell' dimension).");
      }
      else
      {
        INTREPID2_TEST_FOR_EXCEPTION(true, std::invalid_argument,
                                     ">>> ERROR (OrientationTools::modifyBasisByOrientation): input and output ranks must either match, or input rank must be one less than that of output.")
      }

      INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(output.extent(1)) != basis->getCardinality(), std::invalid_argument,
                                    ">>> ERROR (OrientationTools::modifyBasisByOrientation): Field dimension of input/output does not match to basis cardinality.");
    }
#endif

    const shards::CellTopology cellTopo = basis->getBaseCellTopology();
    const ordinal_type  cellDim = cellTopo.getDimension();

    //Initialize output with values from input
    if(input.rank() == output.rank())
      Kokkos::deep_copy(output, input);
    else
      RealSpaceTools<DT>::clone(output, input);

    if ((cellDim < 3) || basis->requireOrientation()) {
      auto ordinalToTag = Kokkos::create_mirror_view_and_copy(typename DT::memory_space(), basis->getAllDofTags());
      auto tagToOrdinal = Kokkos::create_mirror_view_and_copy(typename DT::memory_space(), basis->getAllDofOrdinal());

      const ordinal_type
        numCells  = output.extent(0),
        //numBasis  = output.extent(1),
        numPoints = output.extent(2),
        dimBasis  = output.extent(3); //returns 1 when output.rank() < 4;

      const CoeffMatrixDataViewType matData = createCoeffMatrix(basis);

      ordinal_type numVerts(0), numEdges(0), numFaces(0);

      if (basis->requireOrientation()) {
        numVerts = cellTopo.getVertexCount()*ordinal_type(basis->getDofCount(0, 0) > 0);
        numEdges = cellTopo.getEdgeCount()*ordinal_type(basis->getDofCount(1, 0) > 0);
        numFaces = cellTopo.getFaceCount()*ordinal_type(basis->getDofCount(2, 0) > 0);
      }

      const Kokkos::RangePolicy<typename DT::execution_space> policy(0, numCells);
      typedef F_modifyBasisByOrientation
        <decltype(orts),
         decltype(output),decltype(input),
         decltype(ordinalToTag),decltype(tagToOrdinal),
         decltype(matData)> FunctorType;
      Kokkos::parallel_for
        (policy,
         FunctorType(orts,
                     output, input,
                     ordinalToTag, tagToOrdinal,
                     matData,
                     cellDim, numVerts, numEdges, numFaces,
                     numPoints, dimBasis));
    }
  }
}

#endif