DGUtils.h

Go to the documentation of this file.

// This file is part of the ExaHyPE2 project. For conditions of distribution and
// use, please see the copyright notice at www.peano-framework.org
#pragma once
 
 
#include <vector>
 
#include "tarch/la/la.h"
#include "tarch/la/Vector.h"
#include "tarch/multicore/multicore.h"
#include "tarch/Assertions.h"
 
#include "peano4/grid/GridControlEvent.h"
#include "peano4/datamanagement/CellMarker.h"
#include "peano4/utils/Globals.h"
#include "peano4/utils/Loop.h"
 
#include "Functors.h"
 
#include "exahype2/fv/PatchUtils.h"
 
 
 
namespace exahype2 {
  namespace dg {
    template<typename T>
    tarch::la::Vector<Dimensions,double>  getQuadraturePoint(
      const tarch::la::Vector<Dimensions,double>&  cellCentre,
      const tarch::la::Vector<Dimensions,double>&  cellSize,
      const tarch::la::Vector<Dimensions,int>&     index,
      int                                 polynomialOrder,
      const T* __restrict__          quadraturePoints
    );
 
    template<typename T>
    void  getQuadraturePointInFace(
      const tarch::la::Vector<Dimensions,double>&  faceCentre,
      const tarch::la::Vector<Dimensions,double>&  cellSize,
      const tarch::la::Vector<Dimensions,int>&     index,
      int                                          polynomialOrder,
      int                                          faceOrientation,
      const T* __restrict__          quadraturePoints,
      tarch::la::Vector<Dimensions,double>*       output
    );
 
    double  getQuadratureWeight(
      const tarch::la::Vector<3,double>&  cellSize,
      const tarch::la::Vector<3,int>&     index,
      const double* __restrict__          quadratureWeights
    );
 
 
    int getNodesPerCell(int nodesPerAxis);
 
 
    /* returns a vector with the powers of nodesPerAxis up to the Dimensions - 1
    ** e.g. for dims = 2, nodesPerAxis = 2 will return (1, 2)
    **      for dims = 3, nodesPerAxis = 4 will return (1, 4, 16)
    **
    ** This is used to iterate over the nodes of a cell in a given dimension,
    ** since the cells are all in the same array one after another advancing by
    ** n*strides[i] will move you to the cell which is n steps in direction i
    ** from the current cell.
    ** i.e. node+1 would be the "right" neighbour, node+nodesPerAxis would be the
    ** "upper" neighbour, and so on.
    */
    tarch::la::Vector<Dimensions,int> getStrides(int nodesPerAxis);
 
 
    /*for a given linear node number, returns the index of that node
    **in cartesian coordinate system.
    **e.g. for dims = 2, nodesPerAxis = 4, there will be 16 total nodes
    **  so node 0 is [0,0], 1 is [1,0] and so on until 15 which is [3,3]
    **     for dims = 3, nodesPerAxis = 4, there will be 16*4 = 64 total nodes
    **     ranging from 0 which is [0,0,0] to 63 which is [3,3,3]
    */
    tarch::la::Vector<Dimensions,int> getIndex(
        int                                node,
        tarch::la::Vector<Dimensions,int>  strides
    );
 
    /*
     * for a given node index within the cell, returns the index of the neighbouring
     * node on the surface in a given direction with a given orientation.
     *
     * Ordering on the surface is the same as for the previous functions, but only
     * considering one "side" on each face
     *
     * So with the same notation as the previous function, if given a 2*2 cell
     * in 2 dimensions, there are a total of 4 faces (left,right,down,up), each
     * with 2 nodes.
     * If given index which corresponds to [0,1],
     *   direction 0, orientation 0 would return node 1 on the left face
     *   direction 1, orientation 1 would return node 0 on the upper face
     *       6 7
     *   1 [ o o ] 3
     *   0 [ o o ] 2
     *       4 5
     */
    int cellIndexToHullIndex(
        const tarch::la::Vector<Dimensions,int>&  indexCell,
        const int                                 direction,
        const int                                 orientation,
        const int                                 nodesPerAxis
        );
 
    /*
     * Computes the gradient of Q in a given direction for a given node
     *
     * Mostly used for the non-conservative product.
     *
     * Since the solution is represented by the sum of basis functions, its
     * derivative at a given node is the sum of the derivatives of the basis
     * functions evaluated at that node.
     * The derivative operator contains the value of the derivatives evaluated
     * at the quadrature nodes, therefore it only needs to be multiplied by
     * their coefficients.
     *
     *
     * u(x) = \sum_i \phi_i(x) * u_i => du(x)/dx = d\phi_i(x_i)/dx * u_i
     *
     * @param invDx: the inverse of the size of the cell in any one dimension, this
     * assumes that the cells are cubic
     *
     * @param derivativeOperator: array containing the value of the derivatives of
     * the n base functions of the dg polynomial evaluated at the n quadrature points
     * (in 1 dimension)
     *
     * @param scalarIndex: index of the node at which the gradient should be evaluated
     *
     * @param gradQ: where the gradient values should be placed, has size unknowns*Dimensions
     * since there is a gradient in each direction and the gradient is only nonzero for the
     * unknowns
     *
     * @param strideQ: total number of variables over which the gradient should be computed
     */
    void computeGradient(
        const double* __restrict__ const QCell,
        const double* __restrict__ const derivativeOperator,
        const double                     invDx,
        const int                        nodesPerAxis,
        const int                        strideQ,
        const int                        scalarIndex,
        double* __restrict__             gradQ
        );
 
    /*
     * Substracts contents of one cell array from another.
     */
    void subtractCell(
        double* __restrict__        QOut,
        const double* __restrict__  Qsubstract,
        const int                   order,
        const int                   unknowns,
        const int                   auxiliaryVariables
        );
 
 
    /*
     * Plot patch
     *
     * usually used for debugging
     */
 
    std::string plotCell(
        const double* __restrict__ Q,
        const int order,
        const int unknowns,
        const int auxiliaryVariables
    );
 
    /*
     * Plot Face
     *
     * Usually used for debugging. We plot all unknowns, i.e. both the unknowns
     * form the left and the right adjacent side. The order is linear, i.e. you
     * have to take the normal into account when you want to find out how the
     * unknowns are spatially arranged.
     *
     * See the discussion around ExaHyPE's lexicographic ordering.
     */
    std::string plotFace(
      const double* __restrict__ Q,
      const int order,
      const int unknowns,
      const int auxiliaryVariables,
      int       normal,
      int       numberOfQuantitiesProjectedOntoFace
    );
 
    template<typename QStoreType>
    QStoreType evaluatePolynomial(
      const peano4::datamanagement::CellMarker&    marker,
      int                                          order,
      const double* __restrict__                   QuadratureNodes1d,
      int                                          unknownsPerDoF,
      const QStoreType* __restrict__               Q,
      const tarch::la::Vector<Dimensions,double>&  x,
      int                                          unknown
    );
 
    void copyOneSideOfFaceProjection(
      int    unknownsPlusAuxiliaryVariables,
      int    order,
      int    numberOfProjectedQuantities,
      int    normal,
      int    isRightFaceHalf,
      const double* __restrict__   srcQ,
      double* __restrict__         destQ
    );
  }
}
 
 
#include "DGUtils.cpph"