Commit 25f1eb4e authored by Tom Fischer's avatar Tom Fischer

Merge branch 'ConstexprIf' into 'master'

Constexpr if instead of enable_if.

See merge request ogs/ogs!3469
parents 0ed92c36 df2b56d9
......@@ -76,26 +76,18 @@ inline void computeMappingMatrices(
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline typename std::enable_if<T_SHAPE_FUNC::DIM != 0>::type
computeMappingMatrices(
const T_MESH_ELEMENT& /*ele*/,
const double* natural_pt,
const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
T_SHAPE_MATRICES& shapemat,
FieldType<ShapeMatrixType::DNDR> /*unused*/)
{
double* const dNdr = shapemat.dNdr.data();
T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline typename std::enable_if<T_SHAPE_FUNC::DIM == 0>::type
computeMappingMatrices(
inline void computeMappingMatrices(
const T_MESH_ELEMENT& /*ele*/,
const double* /*natural_pt*/,
[[maybe_unused]] const double* natural_pt,
const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
T_SHAPE_MATRICES& /*shapemat*/,
[[maybe_unused]] T_SHAPE_MATRICES& shapemat,
FieldType<ShapeMatrixType::DNDR> /*unused*/)
{
if constexpr (T_SHAPE_FUNC::DIM != 0)
{
double* const dNdr = shapemat.dNdr.data();
T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
}
}
static void checkJacobianDeterminant(const double detJ,
......@@ -136,45 +128,48 @@ static void checkJacobianDeterminant(const double detJ,
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline typename std::enable_if<T_SHAPE_FUNC::DIM != 0>::type
computeMappingMatrices(
inline void computeMappingMatrices(
const T_MESH_ELEMENT& ele,
const double* natural_pt,
[[maybe_unused]] const double* natural_pt,
const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
T_SHAPE_MATRICES& shapemat,
FieldType<ShapeMatrixType::DNDR_J> /*unused*/)
{
computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
(ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR>());
auto const dim = T_MESH_ELEMENT::dimension;
auto const nnodes = T_MESH_ELEMENT::n_all_nodes;
//jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
for (auto k = decltype(nnodes){0}; k<nnodes; k++) {
const MathLib::Point3d& mapped_pt = ele_local_coord.getMappedCoordinates(k);
// outer product of dNdr and mapped_pt for a particular node
for (auto i_r = decltype(dim){0}; i_r<dim; i_r++) {
for (auto j_x = decltype(dim){0}; j_x<dim; j_x++) {
shapemat.J(i_r,j_x) += shapemat.dNdr(i_r,k) * mapped_pt[j_x];
if constexpr (T_SHAPE_FUNC::DIM != 0)
{
computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
ele,
natural_pt,
ele_local_coord,
shapemat,
FieldType<ShapeMatrixType::DNDR>());
auto const dim = T_MESH_ELEMENT::dimension;
auto const nnodes = T_MESH_ELEMENT::n_all_nodes;
// jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
for (auto k = decltype(nnodes){0}; k < nnodes; k++)
{
const MathLib::Point3d& mapped_pt =
ele_local_coord.getMappedCoordinates(k);
// outer product of dNdr and mapped_pt for a particular node
for (auto i_r = decltype(dim){0}; i_r < dim; i_r++)
{
for (auto j_x = decltype(dim){0}; j_x < dim; j_x++)
{
shapemat.J(i_r, j_x) +=
shapemat.dNdr(i_r, k) * mapped_pt[j_x];
}
}
}
}
shapemat.detJ = shapemat.J.determinant();
checkJacobianDeterminant(shapemat.detJ, ele);
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline typename std::enable_if<T_SHAPE_FUNC::DIM == 0>::type
computeMappingMatrices(
const T_MESH_ELEMENT& /*ele*/,
const double* /*natural_pt*/,
const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
T_SHAPE_MATRICES& shapemat,
FieldType<ShapeMatrixType::DNDR_J> /*unused*/)
{
shapemat.detJ = 1.0;
shapemat.detJ = shapemat.J.determinant();
checkJacobianDeterminant(shapemat.detJ, ele);
}
else
{
shapemat.detJ = 1.0;
}
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
......@@ -192,49 +187,47 @@ inline void computeMappingMatrices(
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline typename std::enable_if<T_SHAPE_FUNC::DIM != 0>::type
computeMappingMatrices(
inline void computeMappingMatrices(
const T_MESH_ELEMENT& ele,
const double* natural_pt,
const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
T_SHAPE_MATRICES& shapemat,
FieldType<ShapeMatrixType::DNDX> /*unused*/)
{
computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
(ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
checkJacobianDeterminant(shapemat.detJ, ele);
//J^-1, dshape/dx
shapemat.invJ.noalias() = shapemat.J.inverse();
auto const nnodes(shapemat.dNdr.cols());
auto const ele_dim(shapemat.dNdr.rows());
assert(shapemat.dNdr.rows()==ele.getDimension());
const unsigned global_dim = ele_local_coord.getGlobalDimension();
if (global_dim==ele_dim) {
shapemat.dNdx.topLeftCorner(ele_dim, nnodes).noalias() = shapemat.invJ * shapemat.dNdr;
} else {
auto const& matR = ele_local_coord.getRotationMatrixToGlobal(); // 3 x 3
auto invJ_dNdr = shapemat.invJ * shapemat.dNdr;
auto dshape_global = matR.topLeftCorner(3u, ele_dim) * invJ_dNdr; //3 x nnodes
shapemat.dNdx = dshape_global.topLeftCorner(global_dim, nnodes);;
computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
ele,
natural_pt,
ele_local_coord,
shapemat,
FieldType<ShapeMatrixType::DNDR_J>());
if constexpr (T_SHAPE_FUNC::DIM != 0)
{
checkJacobianDeterminant(shapemat.detJ, ele);
// J^-1, dshape/dx
shapemat.invJ.noalias() = shapemat.J.inverse();
auto const nnodes(shapemat.dNdr.cols());
auto const ele_dim(shapemat.dNdr.rows());
assert(shapemat.dNdr.rows() == ele.getDimension());
const unsigned global_dim = ele_local_coord.getGlobalDimension();
if (global_dim == ele_dim)
{
shapemat.dNdx.topLeftCorner(ele_dim, nnodes).noalias() =
shapemat.invJ * shapemat.dNdr;
}
else
{
auto const& matR =
ele_local_coord.getRotationMatrixToGlobal(); // 3 x 3
auto invJ_dNdr = shapemat.invJ * shapemat.dNdr;
auto dshape_global =
matR.topLeftCorner(3u, ele_dim) * invJ_dNdr; // 3 x nnodes
shapemat.dNdx = dshape_global.topLeftCorner(global_dim, nnodes);
}
}
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline typename std::enable_if<T_SHAPE_FUNC::DIM == 0>::type
computeMappingMatrices(
const T_MESH_ELEMENT& ele,
const double* natural_pt,
const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
T_SHAPE_MATRICES& shapemat,
FieldType<ShapeMatrixType::DNDX> /*unused*/)
{
computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
(ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline void computeMappingMatrices(
const T_MESH_ELEMENT& ele,
......
......@@ -64,152 +64,116 @@ public:
//! per node.
using LaplaceMatrix = Matrix<Dim*NodalDOF, Dim*NodalDOF>;
//! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF != 0,
decltype(std::declval<const Mat>().template block<Dim, Dim>(0u, 0u))
>::type
blockDimDim(Mat const& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
{
assert(nrows==Dim && ncols==Dim);
(void) nrows; (void) ncols;
return mat.template block<Dim, Dim>(top, left);
}
//! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF == 0,
decltype(std::declval<const Mat>().block(0u, 0u, 0u, 0u))
>::type
blockDimDim(Mat const& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
}
//! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF != 0,
decltype(std::declval<Mat>().template block<Dim, Dim>(0u, 0u))
>::type
blockDimDim(Mat& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
//! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c
//! left.
template <typename Mat>
static auto blockDimDim(Mat const& mat, unsigned top, unsigned left,
unsigned nrows, unsigned ncols)
{
assert(nrows==Dim && ncols==Dim);
(void) nrows; (void) ncols;
return mat.template block<Dim, Dim>(top, left);
}
//! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF == 0,
decltype(std::declval<Mat>().block(0u, 0u, 0u, 0u))
>::type
blockDimDim(Mat& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
if constexpr (NodalDOF != 0)
{
assert(nrows == Dim && ncols == Dim);
(void)nrows;
(void)ncols;
return mat.template block<Dim, Dim>(top, left);
}
else
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
}
}
//! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF != 0,
decltype(std::declval<const Mat>().template block<NNodes, NNodes>(0u, 0u))
>::type
blockShpShp(Mat const& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
{
assert(nrows==NNodes && ncols==NNodes);
(void) nrows; (void) ncols;
return mat.template block<NNodes, NNodes>(top, left);
}
//! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF == 0,
decltype(std::declval<const Mat>().block(0u, 0u, 0u, 0u))
>::type
blockShpShp(Mat const& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
//! Get a block \c Dim x \c Dim whose upper left corner is at \c top and \c
//! left.
template <typename Mat>
static auto blockDimDim(Mat& mat, unsigned top, unsigned left,
unsigned nrows, unsigned ncols)
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
}
//! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF != 0,
decltype(std::declval<Mat>().template block<NNodes, NNodes>(0u, 0u))
>::type
blockShpShp(Mat& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
{
assert(nrows==NNodes && ncols==NNodes);
(void) nrows; (void) ncols;
return mat.template block<NNodes, NNodes>(top, left);
}
//! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top and \c left.
template<typename Mat>
static
typename std::enable_if<NodalDOF == 0,
decltype(std::declval<Mat>().block(0u, 0u, 0u, 0u))
>::type
blockShpShp(Mat& mat,
unsigned top, unsigned left, unsigned nrows, unsigned ncols)
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
if constexpr (NodalDOF != 0)
{
assert(nrows == Dim && ncols == Dim);
(void)nrows;
(void)ncols;
return mat.template block<Dim, Dim>(top, left);
}
else
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
}
}
//! Get a block \c NNodes x 1 starting at the \c top'th row.
template<typename Vec>
static
typename std::enable_if<NodalDOF != 0,
decltype(std::declval<const Vec>().template block<NNodes, 1>(0u, 0u))
>::type
blockShp(Vec const& vec, unsigned top, unsigned nrows)
//! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top
//! and \c left.
template <typename Mat>
static auto blockShpShp(Mat const& mat, unsigned top, unsigned left,
unsigned nrows, unsigned ncols)
{
assert(nrows==NNodes);
(void) nrows;
return vec.template block<NNodes, 1>(top, 0);
if constexpr (NodalDOF != 0)
{
assert(nrows == NNodes && ncols == NNodes);
(void)nrows;
(void)ncols;
return mat.template block<NNodes, NNodes>(top, left);
}
else
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
}
}
//! Get a block \c NNodes x 1 starting at the \c top'th row.
template<typename Vec>
static
typename std::enable_if<NodalDOF == 0,
decltype(std::declval<const Vec>().block(0u, 0u, 0u, 0u))
>::type
blockShp(Vec const& vec, unsigned top, unsigned nrows)
//! Get a block \c NNodes x \c NNodes whose upper left corner is at \c top
//! and \c left.
template <typename Mat>
static auto blockShpShp(Mat& mat, unsigned top, unsigned left,
unsigned nrows, unsigned ncols)
{
return vec.block(top, 0, nrows, 1);
if constexpr (NodalDOF != 0)
{
assert(nrows == NNodes && ncols == NNodes);
(void)nrows;
(void)ncols;
return mat.template block<NNodes, NNodes>(top, left);
}
else
{
assert(nrows == ncols);
return mat.block(top, left, nrows, ncols);
}
}
//! Get a block \c NNodes x 1 starting at the \c top'th row.
template<typename Vec>
static
typename std::enable_if<NodalDOF != 0,
decltype(std::declval<Vec>().template block<NNodes, 1>(0u, 0u))
>::type
blockShp(Vec& vec, unsigned top, unsigned nrows)
template <typename Vec>
static auto blockShp(Vec const& vec, unsigned top, unsigned nrows)
{
assert(nrows==NNodes);
(void) nrows;
return vec.template block<NNodes, 1>(top, 0);
if constexpr (NodalDOF != 0)
{
assert(nrows == NNodes);
(void)nrows;
return vec.template block<NNodes, 1>(top, 0);
}
else
{
return vec.block(top, 0, nrows, 1);
}
}
//! Get a block \c NNodes x 1 starting at the \c top'th row.
template<typename Vec>
static
typename std::enable_if<NodalDOF == 0,
decltype(std::declval<Vec>().block(0u, 0u, 0u, 0u))
>::type
blockShp(Vec& vec, unsigned top, unsigned nrows)
template <typename Vec>
static auto blockShp(Vec& vec, unsigned top, unsigned nrows)
{
return vec.block(top, 0, nrows, 1);
if constexpr (NodalDOF != 0)
{
assert(nrows == NNodes);
(void)nrows;
return vec.template block<NNodes, 1>(top, 0);
}
else
{
return vec.block(top, 0, nrows, 1);
}
}
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment