STYLE: Put local matrix variables in FEM on the stack, remove new/delete

Note that the type of these local matrices, `MatrixType = vnl_matrix<Float>`, which is just empty, when default-constructed. Following C++ Core Guidelines, February 15, 2024, "Prefer scoped objects, don’t heap-allocate unnecessarily", https://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#r5-prefer-scoped-objects-dont-heap-allocate-unnecessarily
InsightSoftwareConsortium · Apr 4, 2024 · 2e1cfba · 2e1cfba
1 parent b3bd01a
commit 2e1cfba
Show file tree

Hide file tree

Showing 4 changed files with 30 additions and 62 deletions.
diff --git a/Modules/Numerics/FEM/src/itkFEMElement2DC0LinearTriangular.cxx b/Modules/Numerics/FEM/src/itkFEMElement2DC0LinearTriangular.cxx
@@ -165,38 +165,33 @@ Element2DC0LinearTriangular::GetLocalFromGlobalCoordinates(const VectorType & gl
 Element2DC0LinearTriangular::Float
 Element2DC0LinearTriangular::JacobianDeterminant(const VectorType & pt, const MatrixType * pJ) const
 {
-
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   Float det = (((*pJ)[1][0] - (*pJ)[0][0]) * ((*pJ)[2][1] - (*pJ)[0][1])) -
               (((*pJ)[0][1] - (*pJ)[1][1]) * ((*pJ)[0][0] - (*pJ)[2][0]));
 
-  delete pJlocal;
-
   return det;
 }
 
 void
 Element2DC0LinearTriangular::JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType * pJ) const
 {
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   // Note that inverse of Jacobian is not quadratic matrix
@@ -210,8 +205,6 @@ Element2DC0LinearTriangular::JacobianInverse(const VectorType & pt, MatrixType &
   invJ[1][0] = idet * ((*pJ)[2][0] - (*pJ)[1][0]);
   invJ[1][1] = idet * ((*pJ)[0][0] - (*pJ)[2][0]);
   invJ[1][2] = idet * ((*pJ)[1][0] - (*pJ)[0][0]);
-
-  delete pJlocal;
 }
 
 void

diff --git a/Modules/Numerics/FEM/src/itkFEMElement2DC0QuadraticTriangular.cxx b/Modules/Numerics/FEM/src/itkFEMElement2DC0QuadraticTriangular.cxx
@@ -115,37 +115,33 @@ Element2DC0QuadraticTriangular::JacobianDeterminant(const VectorType & pt, const
 {
   //  return Superclass::JacobianDeterminant( pt, pJ );
 
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   Float det = (((*pJ)[1][0] - (*pJ)[0][0]) * ((*pJ)[2][1] - (*pJ)[0][1])) -
               (((*pJ)[0][1] - (*pJ)[1][1]) * ((*pJ)[0][0] - (*pJ)[2][0]));
 
-  delete pJlocal;
-
   return det;
 }
 
 void
 Element2DC0QuadraticTriangular::JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType * pJ) const
 {
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   // Note that inverse of Jacobian is not quadratic matrix
@@ -159,8 +155,6 @@ Element2DC0QuadraticTriangular::JacobianInverse(const VectorType & pt, MatrixTyp
   invJ[1][0] = idet * ((*pJ)[2][0] - (*pJ)[1][0]);
   invJ[1][1] = idet * ((*pJ)[0][0] - (*pJ)[2][0]);
   invJ[1][2] = idet * ((*pJ)[1][0] - (*pJ)[0][0]);
-
-  delete pJlocal;
 }
 
 bool

diff --git a/Modules/Numerics/FEM/src/itkFEMElement3DC0LinearTriangular.cxx b/Modules/Numerics/FEM/src/itkFEMElement3DC0LinearTriangular.cxx
@@ -258,15 +258,14 @@ Element3DC0LinearTriangular::JacobianDeterminant(const VectorType & /*HACK pt*/,
 void
 Element3DC0LinearTriangular::JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType * pJ) const
 {
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   //  invJ=vnl_svd_inverse<Float>(*pJ);
@@ -294,8 +293,6 @@ std::cout << (invJ) << std::endl;
 
 std::cout << " invJ2 " << std::endl;
 std::cout << (invJ2) << std::endl;*/
-
-  delete pJlocal;
 }
 
 void

diff --git a/Modules/Numerics/FEM/src/itkFEMElementBase.cxx b/Modules/Numerics/FEM/src/itkFEMElementBase.cxx
@@ -270,15 +270,14 @@ Element::InterpolateSolutionN(const VectorType & pt,
 void
 Element::Jacobian(const VectorType & pt, MatrixType & J, const MatrixType * pshapeD) const
 {
-  MatrixType * pshapeDlocal = nullptr;
+  MatrixType shapeDlocal{};
 
   // If derivatives of shape functions were not provided, we
   // need to compute them here
   if (pshapeD == nullptr)
   {
-    pshapeDlocal = new MatrixType();
-    this->ShapeFunctionDerivatives(pt, *pshapeDlocal);
-    pshapeD = pshapeDlocal;
+    this->ShapeFunctionDerivatives(pt, shapeDlocal);
+    pshapeD = &shapeDlocal;
   }
 
   const unsigned int Nn = pshapeD->columns();
@@ -292,52 +291,42 @@ Element::Jacobian(const VectorType & pt, MatrixType & J, const MatrixType * psha
   }
 
   J = (*pshapeD) * coords;
-
-  // Destroy local copy of derivatives of shape functions, if
-  // they were computed.
-  delete pshapeDlocal;
 }
 
 Element::Float
 Element::JacobianDeterminant(const VectorType & pt, const MatrixType * pJ) const
 {
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   //  Float det=vnl_svd<Float>(*pJ).determinant_magnitude();
   Float det = vnl_qr<Float>(*pJ).determinant();
 
-  delete pJlocal;
-
   return det;
 }
 
 void
 Element::JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType * pJ) const
 {
-  MatrixType * pJlocal = nullptr;
+  MatrixType jlocal{};
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal);
+    pJ = &jlocal;
   }
 
   //  invJ=vnl_svd_inverse<Float>(*pJ);
   invJ = vnl_qr<Float>(*pJ).inverse();
-
-  delete pJlocal;
 }
 
 void
@@ -346,34 +335,29 @@ Element::ShapeFunctionGlobalDerivatives(const VectorType & pt,
                                         const MatrixType * pJ,
                                         const MatrixType * pshapeD) const
 {
-  MatrixType * pshapeDlocal = nullptr;
-  MatrixType * pJlocal = nullptr;
+  MatrixType shapeDlocal{};
+  MatrixType jlocal{};
 
   // If derivatives of shape functions were not provided, we
   // need to compute them here
   if (pshapeD == nullptr)
   {
-    pshapeDlocal = new MatrixType();
-    this->ShapeFunctionDerivatives(pt, *pshapeDlocal);
-    pshapeD = pshapeDlocal;
+    this->ShapeFunctionDerivatives(pt, shapeDlocal);
+    pshapeD = &shapeDlocal;
   }
 
   // If Jacobian was not provided, we
   // need to compute it here
   if (pJ == nullptr)
   {
-    pJlocal = new MatrixType();
-    this->Jacobian(pt, *pJlocal, pshapeD);
-    pJ = pJlocal;
+    this->Jacobian(pt, jlocal, pshapeD);
+    pJ = &jlocal;
   }
 
   MatrixType invJ;
   this->JacobianInverse(pt, invJ, pJ);
 
   shapeDgl = invJ * (*pshapeD);
-
-  delete pJlocal;
-  delete pshapeDlocal;
 }
 
 Element::VectorType