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();
}