diff --git a/NumLib/Fem/InitShapeMatrices.h b/NumLib/Fem/InitShapeMatrices.h
index ab5f8feba2104ac73ce95a0ed223326554073e69..fb2781741bccc188dba1bd10fe54b0fbf32aaa53 100644
--- a/NumLib/Fem/InitShapeMatrices.h
+++ b/NumLib/Fem/InitShapeMatrices.h
@@ -18,11 +18,11 @@ namespace NumLib
{
template <typename ShapeFunction, typename ShapeMatricesType, int GlobalDim,
ShapeMatrixType SelectedShapeMatrixType = ShapeMatrixType::ALL,
- typename IntegrationMethod>
+ typename PointContainer>
std::vector<typename ShapeMatricesType::ShapeMatrices,
Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
-initShapeMatrices(MeshLib::Element const& e, bool is_axially_symmetric,
- IntegrationMethod const& integration_method)
+computeShapeMatrices(MeshLib::Element const& e, bool const is_axially_symmetric,
+ PointContainer const& points)
{
std::vector<
typename ShapeMatricesType::ShapeMatrices,
@@ -32,22 +32,41 @@ initShapeMatrices(MeshLib::Element const& e, bool is_axially_symmetric,
auto const fe =
createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);
- unsigned const n_integration_points =
- integration_method.getNumberOfPoints();
-
- shape_matrices.reserve(n_integration_points);
- for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ shape_matrices.reserve(points.size());
+ for (auto const& p : points)
{
shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
ShapeFunction::NPOINTS);
- fe.computeShapeFunctions(
- integration_method.getWeightedPoint(ip).getCoords(),
- shape_matrices[ip], GlobalDim, is_axially_symmetric);
+ fe.template computeShapeFunctions<SelectedShapeMatrixType>(
+ p.getCoords(), shape_matrices.back(), GlobalDim,
+ is_axially_symmetric);
}
return shape_matrices;
}
+template <typename ShapeFunction, typename ShapeMatricesType, int GlobalDim,
+ ShapeMatrixType SelectedShapeMatrixType = ShapeMatrixType::ALL,
+ typename IntegrationMethod>
+std::vector<typename ShapeMatricesType::ShapeMatrices,
+ Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
+initShapeMatrices(MeshLib::Element const& e, bool const is_axially_symmetric,
+ IntegrationMethod const& integration_method)
+{
+ int const n_integration_points = integration_method.getNumberOfPoints();
+
+ std::vector<typename IntegrationMethod::WeightedPoint> points;
+ points.reserve(n_integration_points);
+ for (int ip = 0; ip < n_integration_points; ++ip)
+ {
+ points.push_back(integration_method.getWeightedPoint(ip));
+ }
+
+ return computeShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim,
+ SelectedShapeMatrixType>(
+ e, is_axially_symmetric, points);
+}
+
template <typename ShapeFunction, typename ShapeMatricesType>
double interpolateXCoordinate(
MeshLib::Element const& e,
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendrePrism.h b/NumLib/Fem/Integration/IntegrationGaussLegendrePrism.h
index e87e50baead3275754c299f7e6b65b84492e4dca..9dc77292531b5259f61c7201dcc31076ac447328 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendrePrism.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendrePrism.h
@@ -21,9 +21,9 @@ namespace NumLib
*/
class IntegrationGaussLegendrePrism
{
+public:
using WeightedPoint = MathLib::TemplateWeightedPoint<double, double, 3>;
-public:
/**
* Construct this object with the given integration order
*
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h b/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h
index f02fe46e7d853b8d6e9f60c497a2e6f9b1d2d2a3..4e3aee5e1f0791ae3b421bd627e2366eff01986e 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h
@@ -20,9 +20,9 @@ namespace NumLib
*/
class IntegrationGaussLegendrePyramid
{
+public:
using WeightedPoint = MathLib::TemplateWeightedPoint<double, double, 3>;
-public:
/**
* Construct this object with the given integration order
*
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h b/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h
index ab9c7e26bad2a4859c04ee43fff14ade5404b960..1f3f55991b9244624ab1f16a0679ee827578b789 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h
@@ -27,10 +27,10 @@ namespace NumLib
template <unsigned N_DIM>
class IntegrationGaussLegendreRegular
{
+public:
using WeightedPoint =
typename MathLib::TemplateWeightedPoint<double, double, N_DIM>;
-public:
/// Create IntegrationGaussLegendreRegular of the given Gauss-Legendre
/// integration order.
///
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h b/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h
index 2fd50e2457bdda409b7c7e6ff278efc882133bd4..a797e4e7cd70f87dbe88c174112ec18ff4c1ff84 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h
@@ -20,9 +20,9 @@ namespace NumLib
*/
class IntegrationGaussLegendreTet
{
+public:
using WeightedPoint = MathLib::TemplateWeightedPoint<double, double, 3>;
-public:
/**
* Construct this object with the given integration order
*
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h b/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h
index f472b1430b7ad87d6d5135e15af02cfcd898e075..3ab76037dceb181b5662ff61b78a5317874c997a 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h
@@ -35,9 +35,9 @@ namespace NumLib
*/
class IntegrationGaussLegendreTri
{
+public:
using WeightedPoint = MathLib::TemplateWeightedPoint<double, double, 2>;
-public:
/**
* Construct this object with the given integration order
*
diff --git a/NumLib/Fem/Integration/IntegrationPoint.h b/NumLib/Fem/Integration/IntegrationPoint.h
index fdb72ee33752047b609a03bf8c55ec209a523214..a6cd0f08fbb0f3677ac57de6fcaec99dabc9a882 100644
--- a/NumLib/Fem/Integration/IntegrationPoint.h
+++ b/NumLib/Fem/Integration/IntegrationPoint.h
@@ -20,9 +20,9 @@ namespace NumLib
/// It is only needed to satisfy the common integration rule concepts.
class IntegrationPoint
{
+public:
using WeightedPoint = MathLib::TemplateWeightedPoint<double, double, 1>;
-public:
/// IntegrationPoint constructor for given order.
explicit IntegrationPoint(unsigned /* order */)
{
diff --git a/NumLib/Function/Interpolation.h b/NumLib/Function/Interpolation.h
index 06408b1939bd3d39c3aeef0949efdd71b3fe9eb4..0d6bb7b6ae661a1e313062e8cce6a358e07bf904 100644
--- a/NumLib/Function/Interpolation.h
+++ b/NumLib/Function/Interpolation.h
@@ -15,6 +15,7 @@
#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
namespace NumLib
{
@@ -108,11 +109,7 @@ void interpolateToHigherOrderNodes(
using SF = LowerOrderShapeFunction;
using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
- using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
- auto const fe =
- NumLib::createIsoparametricFiniteElement<SF, ShapeMatricesType>(
- element);
int const number_base_nodes = element.getNumberOfBaseNodes();
int const number_all_nodes = element.getNumberOfNodes();
@@ -123,21 +120,32 @@ void interpolateToHigherOrderNodes(
interpolated_values_global_vector[global_index] = node_values[n];
}
- // Shape matrices storage reused in the interpolation loop.
- ShapeMatrices shape_matrices{SF::DIM, GlobalDim, SF::NPOINTS};
-
+ //
// Interpolate values for higher order nodes.
- for (int n = number_base_nodes; n < number_all_nodes; ++n)
+ //
+
+ int const number_higher_order_nodes = number_all_nodes - number_base_nodes;
+ std::vector<MathLib::Point3d> higher_order_nodes;
+ higher_order_nodes.reserve(number_higher_order_nodes);
+ for (int n = 0; n < number_higher_order_nodes; ++n)
{
- // Evaluated at higher order nodes' coordinates.
- fe.template computeShapeFunctions<ShapeMatrixType::N>(
- NaturalCoordinates<HigherOrderMeshElementType>::coordinates[n]
- .data(),
- shape_matrices, GlobalDim, is_axially_symmetric);
+ higher_order_nodes.emplace_back(
+ NaturalCoordinates<HigherOrderMeshElementType>::coordinates
+ [number_base_nodes + n]);
+ }
- std::size_t const global_index = element.getNodeIndex(n);
+ // Shape matrices evaluated at higher order nodes' coordinates.
+ auto const shape_matrices =
+ computeShapeMatrices<SF, ShapeMatricesType, GlobalDim,
+ ShapeMatrixType::N>(element, is_axially_symmetric,
+ higher_order_nodes);
+
+ for (int n = 0; n < number_higher_order_nodes; ++n)
+ {
+ std::size_t const global_index =
+ element.getNodeIndex(number_base_nodes + n);
interpolated_values_global_vector[global_index] =
- shape_matrices.N * node_values;
+ shape_matrices[n].N * node_values;
}
}
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 2c3b5f5d5162e0dc9f9f3f8f4489271502dcdf8b..235e8d0f147e63fe9538428492eddeed434d9ddb 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -854,21 +854,14 @@ public:
auto const local_C = Eigen::Map<const NodalVectorType>(
&local_x[first_concentration_index], concentration_size);
- // eval shape matrices at given point
- auto const shape_matrices = [&]() {
- auto const fe = NumLib::createIsoparametricFiniteElement<
- ShapeFunction, ShapeMatricesType>(_element);
-
- typename ShapeMatricesType::ShapeMatrices sm(
- ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
-
- // Note: Axial symmetry is set to false here, because we only need
- // dNdx here, which is not affected by axial symmetry.
- fe.template computeShapeFunctions<NumLib::ShapeMatrixType::DNDX>(
- pnt_local_coords.getCoords(), sm, GlobalDim, false);
-
- return sm;
- }();
+ // Eval shape matrices at given point
+ // Note: Axial symmetry is set to false here, because we only need dNdx
+ // here, which is not affected by axial symmetry.
+ auto const shape_matrices =
+ NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
+ GlobalDim>(
+ _element, false /*is_axially_symmetric*/,
+ std::array{pnt_local_coords})[0];
ParameterLib::SpatialPosition pos;
pos.setElementID(_element.getID());
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 701d448e400a11492301f1a26807013def616d8b..643a114e91772a6fc6e8a40dbd86d06198671961 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -95,17 +95,14 @@ public:
double const t,
std::vector<double> const& local_x) const override
{
- // eval dNdx and invJ at given point
- auto const fe = NumLib::createIsoparametricFiniteElement<
- ShapeFunction, ShapeMatricesType>(_element);
-
- typename ShapeMatricesType::ShapeMatrices shape_matrices(
- ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
-
+ // Eval shape matrices at given point
// Note: Axial symmetry is set to false here, because we only need dNdx
// here, which is not affected by axial symmetry.
- fe.computeShapeFunctions(pnt_local_coords.getCoords(), shape_matrices,
- GlobalDim, false);
+ auto const shape_matrices =
+ NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
+ GlobalDim>(
+ _element, false /*is_axially_symmetric*/,
+ std::array{pnt_local_coords})[0];
ParameterLib::SpatialPosition pos;
pos.setElementID(this->_element.getID());
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 5b4922097f54afd8b57e1aaa7cd15fa4c07471c9..56e483bea270400306dfe5fc1ba3ad5a837e0324 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -76,18 +76,13 @@ LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::getFlux(
// extension of getFlux interface
double const dt = std::numeric_limits<double>::quiet_NaN();
- // eval dNdx and invJ at p
- auto const fe =
- NumLib::createIsoparametricFiniteElement<ShapeFunction,
- ShapeMatricesType>(_element);
-
- typename ShapeMatricesType::ShapeMatrices shape_matrices(
- ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
-
// Note: Axial symmetry is set to false here, because we only need dNdx
// here, which is not affected by axial symmetry.
- fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
- GlobalDim, false);
+ auto const shape_matrices =
+ NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
+ GlobalDim>(_element,
+ false /*is_axially_symmetric*/,
+ std::array{p_local_coords})[0];
// create pos object to access the correct media property
ParameterLib::SpatialPosition pos;
diff --git a/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h b/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h
index 2bcbabbc3b88d13422764267b7118cba6d1b5e2f..4cb5664408e09044fd62e110099a3d2df368942f 100644
--- a/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h
+++ b/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h
@@ -136,17 +136,14 @@ public:
// extension of getFlux interface
double const dt = std::numeric_limits<double>::quiet_NaN();
- // eval dNdx and invJ at p
- auto const fe = NumLib::createIsoparametricFiniteElement<
- ShapeFunction, ShapeMatricesType>(_element);
-
- typename ShapeMatricesType::ShapeMatrices shape_matrices(
- ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
-
+ // Eval shape matrices at given point
// Note: Axial symmetry is set to false here, because we only need dNdx
// here, which is not affected by axial symmetry.
- fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
- GlobalDim, false);
+ auto const shape_matrices =
+ NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
+ GlobalDim>(
+ _element, false /*is_axially_symmetric*/,
+ std::array{p_local_coords})[0];
// fetch hydraulic conductivity
ParameterLib::SpatialPosition pos;
diff --git a/Tests/NumLib/NaturalNodeCoordinates.cpp b/Tests/NumLib/NaturalNodeCoordinates.cpp
index 6166382d3470cf54d62b90fd5aeb922e714a4f11..93d44ad9310ce578c6cfb922ffb761d9f7543442 100644
--- a/Tests/NumLib/NaturalNodeCoordinates.cpp
+++ b/Tests/NumLib/NaturalNodeCoordinates.cpp
@@ -7,15 +7,16 @@
*
*/
+#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
+
#include <gtest/gtest.h>
#include <vector>
#include "MeshLib/ElementCoordinatesMappingLocal.h"
-
-#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -25,8 +26,6 @@ template <typename ShapeFunction, int GlobalDim>
bool test(MeshLib::Element const& element)
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
- using FemType =
- NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
// Test if the cast is possible.
bool const dynamic_cast_failed =
@@ -37,29 +36,31 @@ bool test(MeshLib::Element const& element)
{
return false;
}
- FemType fe{
- static_cast<typename ShapeFunction::MeshElement const&>(element)};
-
- bool result = true;
// For each node evaluate the shape functions at natural coordinates and
// test if only the corresponding shape function has value 1 and all others
// must return 0.
int const number_nodes = element.getNumberOfNodes();
+ std::vector<MathLib::Point3d> points;
+ points.reserve(number_nodes);
for (int n = 0; n < number_nodes; ++n)
{
// Evaluate shape matrices at natural coordinates.
- typename ShapeMatricesType::ShapeMatrices shape_matrices{
- ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS};
// Compute only N, because for pyramid the detJ becomes zero at the
// apex, and we only use N in the following test anyway.
- fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N>(
+ points.emplace_back(
NumLib::NaturalCoordinates<
- typename ShapeFunction::MeshElement>::coordinates[n]
- .data(),
- shape_matrices, GlobalDim, false /* axial symmetry */);
+ typename ShapeFunction::MeshElement>::coordinates[n]);
+ }
+ auto const shape_matrices =
+ NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
+ GlobalDim, NumLib::ShapeMatrixType::N>(
+ element, false /*is_axially_symmetric*/, points);
- auto const& N = shape_matrices.N;
+ bool result = true;
+ for (int n = 0; n < number_nodes; ++n)
+ {
+ auto const& N = shape_matrices[n].N;
for (int p = 0; p < static_cast<int>(ShapeFunction::NPOINTS); ++p)
{
if (p == n)
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 395f0b369508a3ae99e9bf43797a18baa08baf3b..d7b935057960c2aff0060a4f31686e5d37067578 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -265,22 +265,22 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
TYPED_TEST(NumLibFemIsoTest, CheckGaussLegendreIntegrationLevel)
{
// Refer to typedefs in the fixture
- using FeType = typename TestFixture::FeType;
using NodalMatrix = typename TestFixture::NodalMatrix;
- using ShapeMatricesType = typename TestFixture::ShapeMatricesType;
- // create a finite element object with gauss quadrature level 2
- FeType fe(*this->mesh_element);
+ auto shape_matrices =
+ NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
+ typename TestFixture::ShapeMatrixTypes,
+ TestFixture::global_dim>(
+ *this->mesh_element, false /*is_axially_symmetric*/,
+ this->integration_method);
// evaluate a mass matrix
NodalMatrix M(this->e_nnodes, this->e_nnodes);
M.setZero();
- ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
ASSERT_EQ(TestFixture::n_sample_pt_order2, this->integration_method.getNumberOfPoints());
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
- shape.setZero();
+ auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
}
//std::cout << "M=\n" << M << std::endl;
@@ -290,10 +290,17 @@ TYPED_TEST(NumLibFemIsoTest, CheckGaussLegendreIntegrationLevel)
this->integration_method.setIntegrationOrder(3);
M *= .0;
ASSERT_EQ(TestFixture::n_sample_pt_order3, this->integration_method.getNumberOfPoints());
+
+ shape_matrices =
+ NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
+ typename TestFixture::ShapeMatrixTypes,
+ TestFixture::global_dim>(
+ *this->mesh_element, false /*is_axially_symmetric*/,
+ this->integration_method);
+
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
- shape.setZero();
+ auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
M.noalias() += shape.N * shape.N.transpose() * shape.detJ * wp.getWeight();
}
//std::cout << "M=\n" << M << std::endl;
diff --git a/Tests/NumLib/TestGradShapeFunction.cpp b/Tests/NumLib/TestGradShapeFunction.cpp
index d8457ecc3febb9d77c326c5c3d182931a4ef6925..daaf40b3f50f164f284399740875e9ff7f20c9b1 100644
--- a/Tests/NumLib/TestGradShapeFunction.cpp
+++ b/Tests/NumLib/TestGradShapeFunction.cpp
@@ -14,36 +14,33 @@
#include <gtest/gtest.h>
-#include <vector>
-#include <cmath>
-
#include <Eigen/Eigen>
+#include <cmath>
+#include <vector>
-#include "MeshLib/ElementCoordinatesMappingLocal.h"
-#include "MeshLib/CoordinateSystem.h"
-
-#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
-#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
-
+#include "FeTestData/TestFeHEX20.h"
+#include "FeTestData/TestFeHEX8.h"
#include "FeTestData/TestFeLINE2.h"
#include "FeTestData/TestFeLINE2Y.h"
#include "FeTestData/TestFeLINE3.h"
-#include "FeTestData/TestFeTRI3.h"
-#include "FeTestData/TestFeTRI6.h"
+#include "FeTestData/TestFePRISM15.h"
+#include "FeTestData/TestFePRISM6.h"
+#include "FeTestData/TestFePYRA13.h"
+#include "FeTestData/TestFePYRA5.h"
#include "FeTestData/TestFeQUAD4.h"
#include "FeTestData/TestFeQUAD8.h"
#include "FeTestData/TestFeQUAD9.h"
-#include "FeTestData/TestFeHEX8.h"
-#include "FeTestData/TestFeHEX20.h"
-#include "FeTestData/TestFeTET4.h"
#include "FeTestData/TestFeTET10.h"
-#include "FeTestData/TestFePRISM6.h"
-#include "FeTestData/TestFePRISM15.h"
-#include "FeTestData/TestFePYRA5.h"
-#include "FeTestData/TestFePYRA13.h"
-
+#include "FeTestData/TestFeTET4.h"
+#include "FeTestData/TestFeTRI3.h"
+#include "FeTestData/TestFeTRI6.h"
+#include "MeshLib/CoordinateSystem.h"
+#include "MeshLib/ElementCoordinatesMappingLocal.h"
+#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "Tests/TestTools.h"
using namespace NumLib;
@@ -110,11 +107,6 @@ public:
// Finite element type
template <typename X>
using ShapeMatrixPolicy = typename T::template ShapeMatrixPolicy<X>;
- using FeType =
- typename TestFeType::template FeType<ShapeMatrixPolicy>::type;
-
- // Shape matrix data type
- using ShapeMatricesType = typename ShapeMatrixTypes::ShapeMatrices;
using MeshElementType = typename TestFeType::MeshElementType;
static const unsigned dim = TestFeType::dim;
@@ -182,23 +174,21 @@ TYPED_TEST_CASE(GradShapeFunctionTest, TestTypes);
TYPED_TEST(GradShapeFunctionTest,
CheckGradShapeFunctionByComputingElementVolume)
{
- // Refer to typedefs in the fixture
- using FeType = typename TestFixture::FeType;
- using ShapeMatricesType = typename TestFixture::ShapeMatricesType;
-
- // create a finite element object
- FeType fe(*this->mesh_element);
+ auto const shape_matrices =
+ NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
+ typename TestFixture::ShapeMatrixTypes,
+ TestFixture::global_dim,
+ ShapeMatrixType::N_J>(
+ *this->mesh_element, false /*is_axially_symmetric*/,
+ this->integration_method);
// Compute element volume numerically as V_e = int {1}dA_e
double computed_element_volume = 0.;
- ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
for (std::size_t i = 0; i < this->integration_method.getNumberOfPoints();
i++)
{
- shape.setZero();
+ auto const& shape = shape_matrices[i];
auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::N_J>(
- wp.getCoords(), shape, this->global_dim, false);
computed_element_volume += shape.detJ * wp.getWeight();
}