Select Git revision
JacobianBasis.cpp
JacobianBasis.cpp 34.94 KiB
// 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);
}