dgcartesianmortar.h

Go to the documentation of this file.

// Copyright (C) 2024 The m-AIA AUTHORS
//
// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
//
// SPDX-License-Identifier: LGPL-3.0-only
 
#ifndef DGMORTAR_H_
#define DGMORTAR_H_
 
#include "INCLUDE/maiatypes.h"
#include "UTIL/tensor.h"
#include "dgcartesianinterpolation.h"
 
#include <string>
 
namespace maia {
namespace dg {
namespace mortar {
 
 
const MBool forward = true;
const MBool reverse = false;
const MInt lower = 0;
const MInt upper = 1;
 
 
inline void calcMortarProjectionMatrixP(const MInt sourcePolyDeg,
                                        const MFloat* const sourceNodes,
                                        const MFloat* const sourceBaryWeights,
                                        const MInt targetPolyDeg,
                                        const MFloat* const targetNodes,
                                        MFloat* const matrix) {
  const MInt sourceNoNodes = sourcePolyDeg + 1;
  const MInt targetNoNodes = targetPolyDeg + 1;
 
  // Matrix is Lagrange polynomials of source polynomial degree evaluated at
  // nodes of target polynomial degree
  for(MInt i = 0; i < targetNoNodes; i++) {
    interpolation::calcLagrangeInterpolatingPolynomials(
        targetNodes[i], sourcePolyDeg, sourceNodes, sourceBaryWeights, &matrix[i * sourceNoNodes]);
  }
}
 
 
inline void calcMortarProjectionMatrixHForward(const MInt polyDeg,
                                               const MFloat* const nodes,
                                               const MFloat* const wBary,
                                               const MInt position,
                                               MFloat* const matrix) {
  const MInt noNodes = polyDeg + 1;
 
  // Forward matrix is the Lagrange polynomials evaluated at other nodes
  const MFloat shift = (position == dg::mortar::lower) ? -0.5 : 0.5;
  for(MInt i = 0; i < noNodes; i++) {
    maia::dg::interpolation::calcLagrangeInterpolatingPolynomials(
        0.5 * nodes[i] + shift, polyDeg, nodes, wBary, &matrix[i * noNodes]);
  }
}
 
 
inline void calcMortarProjectionMatrixHReverse(const MInt polyDeg,
                                               const MFloat* const nodes,
                                               const MFloat* const wBary,
                                               const MInt position,
                                               MFloat* const matrix) {
  const MInt noNodes = polyDeg + 1;
 
  // Calculate nodes and quadrature weights for projection (always Gauss)
  MFloatVector projectionNodes;
  MFloatVector projectionWInt;
  projectionNodes.resize(noNodes);
  projectionWInt.resize(noNodes);
  maia::dg::interpolation::calcLegendreGaussNodesAndWeights(polyDeg, &projectionNodes[0], &projectionWInt[0]);
 
  // Calculate barycentric weights for projection nodes
  MFloatVector projectionWBary;
  projectionWBary.resize(noNodes);
  maia::dg::interpolation::calcBarycentricWeights(polyDeg, &projectionNodes[0], &projectionWBary[0]);
 
  // Calculate projection matrix from mortar element to coarse element in projection nodes
  const MFloat shift = (position == dg::mortar::lower) ? -0.5 : 0.5;
  MFloatTensor mortarToCoarseElement(noNodes, noNodes);
  for(MInt i = 0; i < noNodes; i++) {
    maia::dg::interpolation::calcLagrangeInterpolatingPolynomials(0.5 * projectionNodes[i] + shift,
                                                                  polyDeg,
                                                                  &projectionNodes[0],
                                                                  &projectionWBary[0],
                                                                  &mortarToCoarseElement(i, 0));
  }
 
  // Transpose and multiply
  for(MInt i = 0; i < noNodes; i++) {
    for(MInt j = 0; j < i + 1; j++) {
      const MFloat temp = mortarToCoarseElement(i, j);
      mortarToCoarseElement(i, j) = 0.5 * mortarToCoarseElement(j, i) * projectionWInt[j] / projectionWInt[i];
      mortarToCoarseElement(j, i) = 0.5 * temp * projectionWInt[i] / projectionWInt[j];
    }
  }
 
 
  // Transformation matrix from grid to projection nodes (this is the identity when using Gauss
  // nodes anyways)
  MFloatTensor p2dg(noNodes, noNodes);
  calcMortarProjectionMatrixP(polyDeg, &projectionNodes[0], &projectionWBary[0], polyDeg, nodes, &p2dg[0]);
 
  // Transformation matrix from projection to grid nodes (this is the identity when using Gauss
  // nodes anyways)
  MFloatTensor dg2p(noNodes, noNodes);
  calcMortarProjectionMatrixP(polyDeg, nodes, wBary, polyDeg, &projectionNodes[0], &dg2p[0]);
 
  // The total projection matrix m is a product of the transformation matrices to(dg2p) and
  // from(p2dg) quadrature nodes and the projection matrix in quadrature nodes, in this order: m =
  // p2dg * mortarToCoarseElement * dg2p (where "*" denotes a matrix-matrix product)
  MFloatTensor m(matrix, noNodes, noNodes);
  m.set(0.0);
  for(MInt q = 0; q < noNodes; q++) {
    for(MInt j = 0; j < noNodes; j++) {
      for(MInt i = 0; i < noNodes; i++) {
        for(MInt n = 0; n < noNodes; n++) {
          m(q, j) += p2dg(q, i) * mortarToCoarseElement(i, n) * dg2p(n, j);
        }
      }
    }
  }
}
 
} // namespace mortar
} // namespace dg
} // namespace maia
 
#endif // define DGMORTAR_H_