JacobianBasis.cpp

// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@onelab.info>.

#include "JacobianBasis.h"
#include "pointsGenerators.h"
#include "nodalBasis.h"
#include "bezierBasis.h"
#include "BasisFactory.h"
#include "Numeric.h"
#include <cmath>

namespace {

template<class T>
void calcMapFromIdealElement(int type, T &gSMatX, T &gSMatY, T &gSMatZ)
{
  // 2D scaling
  switch(type) {
    case TYPE_QUA:                                             // Quad, hex, pyramid -> square with side of length 1
    case TYPE_HEX:
    case TYPE_PYR: {
      gSMatX.scale(2.);
      gSMatY.scale(2.);
      break;
    }
    default: {                                                // Tri, tet, prism: equilateral tri with side of length 1
      static const double cTri[2] = {-1./std::sqrt(3.), 2./std::sqrt(3.)};
      gSMatY.scale(cTri[1]);
      gSMatY.axpy(gSMatX, cTri[0]);
      break;
    }
  }

  // 3D scaling
  switch(type) {
    case TYPE_HEX:                                            // Hex, prism -> side of length 1 in z
    case TYPE_PRI: {
      gSMatZ.scale(2.);
      break;
    }
    case TYPE_PYR: {                                          // Pyramid -> height sqrt(2.)/2
      static const double cPyr = sqrt(2.);
      gSMatZ.scale(cPyr);
      break;
    }
    case TYPE_TET:                                            // Tet: take into account (x, y) scaling to obtain regular tet
    {
      static const double cTet[3] = {-3./2/std::sqrt(6.),
                                     -1./2/std::sqrt(2.),
                                     std::sqrt(1.5)};
      gSMatZ.scale(cTet[2]);
      gSMatZ.axpy(gSMatX, cTet[0]);
      gSMatZ.axpy(gSMatY, cTet[1]);
      break;
    }
  }
}

// Compute the determinant of a 3x3 matrix
inline double calcDet3D(double M11, double M12, double M13,
                        double M21, double M22, double M23,
                        double M31, double M32, double M33)
{
  return M11 * (M22*M33 - M23*M32)
       - M12 * (M21*M33 - M23*M31)
       + M13 * (M21*M32 - M22*M31);
}