From 52030637db13eabcbc8b1c8eec43171a62281728 Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Tue, 17 Mar 2026 16:40:29 +0100
Subject: [PATCH 1/6] Added: class IFEMNastranReader in separate implementation
 file. Fixed: Spell error in console message.

---
 ASMu2DNastran.C       |  97 ++++++-------------------------------
 CMakeLists.txt        |   7 +++
 IFEMNastranReader.C   | 108 ++++++++++++++++++++++++++++++++++++++++++
 IFEMNastranReader.h   |  41 ++++++++++++++++
 Test/DroppedBox.reg   |   2 +-
 Test/Fine-dump.reg    |   2 +-
 Test/Fine-load.reg    |   2 +-
 Test/Fine-partial.reg |   2 +-
 Test/Q4-dyn.reg       |   2 +-
 Test/Q4-dynp.reg      |   2 +-
 Test/Robaat-irr.reg   |   2 +-
 Test/Robaat.reg       |   2 +-
 Test/Wheel.reg        |   2 +-
 main.C                |  15 ++++--
 14 files changed, 190 insertions(+), 96 deletions(-)
 create mode 100644 IFEMNastranReader.C
 create mode 100644 IFEMNastranReader.h

diff --git a/ASMu2DNastran.C b/ASMu2DNastran.C
index 5a70f0e..3c958b2 100644
--- a/ASMu2DNastran.C
+++ b/ASMu2DNastran.C
@@ -25,20 +25,19 @@
 
 namespace ASM {
   char useBeam = 1; //!< If nonzero, include beam elements as a separate patch
-  char skipMass = 0; //!< 1: ignore one-noded mass elements, 2: no VTF output
-  bool skipRBE2 = false; //!< If \e true, ignore all RBE2 elements
   bool replRBE3 = false; //!< If \e true, convert RBE3 elements to RBE2 elements
   bool skipLoad = false; //!< If \e true, ignore all FE surface loads
-  bool readSets = true; //!< If \e true, read pre-bulk Nastran SET definitions
   double Ktra = 0.0; //!< Translation stiffness for added springs in slave nodes
   double Krot = 0.0; //!< Rotation stiffness for added sprins in slave nodes
   IntVec fixRBE3; //!< List of RBE3 elements to be constrained
-  IntVec ignoredElms; //!< List of elements to ignore
+#ifdef HAS_FFLLIB
+  extern char skipMass;
+#endif
 }
 
 #ifdef HAS_FFLLIB
+#include "IFEMNastranReader.h"
 #include "FFlLinkHandler.H"
-#include "FFlNastranReader.H"
 #include "FFlElementBase.H"
 #include "FFlLoadBase.H"
 #include "FFlGroup.H"
@@ -58,60 +57,6 @@ namespace ASM {
 namespace FFlNastran {
   extern std::string mainPath;
 }
-
-
-/*!
-  \brief Class for reading Nastran bulk data into a FFlLinkHandler object.
-*/
-
-class MyNastranReader : public FFlNastranReader
-{
-  int nPreBulk; //!< Number of lines before BEGIN BULK
-
-public:
-  //! \brief The constructor forwards to the parent class constructor.
-  MyNastranReader(FFlLinkHandler& fePart, int lCount)
-    : FFlNastranReader(&fePart,lCount), nPreBulk(lCount) {}
-
-  //! \brief Reads the FE model from the Nastran file stream.
-  bool readFE(std::istream& is, std::istream& iset)
-  {
-    PROFILE("Nastran file parser");
-
-    if (!this->resolve(this->read(is)))
-      myLink->deleteGeometry(); // Parsing failure, delete all FE data
-    else if (nWarnings+nNotes > 0)
-      IFEM::cout <<"\n  ** Parsing FE data succeeded."
-                 <<"\n     However, "<< nWarnings
-                 <<" warning(s) and "<< nNotes <<" note(s) were reported.\n"
-                 <<"     Review the messages and check the FE data file.\n"
-                 << std::endl;
-    if (!myLink->hasGeometry())
-      return false;
-
-    // Remove all solid elements (not yet supported),
-    // and (optionally) all the mass and beam elements, and ...
-    ElementsVec toBeErased;
-    ElementsCIter it;
-    for (it = myLink->elementsBegin(); it != myLink->elementsEnd(); ++it)
-      if ((*it)->getCathegory() == FFlTypeInfoSpec::SOLID_ELM ||
-          utl::findIndex(ASM::ignoredElms,(*it)->getID()) >= 0 ||
-          (ASM::skipRBE2 && (*it)->getTypeName() == "RGD") ||
-          (ASM::skipMass == 1 && (*it)->getTypeName() == "CMASS") ||
-          (!ASM::useBeam && (*it)->getCathegory() == FFlTypeInfoSpec::BEAM_ELM))
-        toBeErased.push_back(*it);
-    if (!toBeErased.empty())
-    {
-      IFEM::cout <<"  ** Erasing "<< toBeErased.size()
-                 <<" elements from the model."<< std::endl;
-      myLink->removeElements(toBeErased);
-    }
-
-    // Now parse the element set definitions, if any
-    lastComment = { 0, "" };
-    return iset ? this->processAllSets(iset,nPreBulk) : true;
-  }
-};
 #endif
 
 
@@ -155,38 +100,21 @@ bool ASMu2DNastran::read (std::istream& is)
   int nErr = 0;
   int iErr = 0;
 
-  // Fast-forward until "BEGIN BULK"
-  int lCount = 0;
-  char cline[256];
+#ifdef HAS_FFLLIB
   std::stringstream sets;
-  while (is.getline(cline,255))
-    if (!strncmp(cline,"BEGIN BULK",10))
-      break;
-    else
-    {
-      ++lCount;
-      // Copy element SET definitions to a second stream,
-      // since they have to be parsed after the FE model is loaded
-      if (ASM::readSets && (!strncmp(cline,"SET ",4) || sets.tellp() > 0))
-        sets << cline << '\n';
-    }
-
-  if (!is) return false; // No bulk data file
-
-  if (lCount > 0)
-    IFEM::cout <<"\tNastran bulk data starting at line "
-               << lCount+1 << std::endl;
+  int lCount = IFEMNastranReader::parseSets(is,sets);
+  if (lCount < 0) return false; // No bulk data file
 
-#ifdef HAS_FFLLIB
   bool ok = false;
   FFlLinkHandler fem;
-  if (MyNastranReader reader(fem,lCount); reader.readFE(is,sets))
+  if (IFEMNastranReader reader(fem,lCount);
+      reader.readFE(is,sets,!ASM::useBeam))
   {
     PROFILE("Resolve cross references");
     ok = fem.resolve();
   }
   if (ok)
-    IFEM::cout <<"\nParsing Nastran bulk data file succceeded."<< std::endl;
+    IFEM::cout <<"\nParsing Nastran bulk data file succeeded."<< std::endl;
   else
   {
     std::cerr <<"\n *** Parsing/resolving FE data failed.\n"
@@ -1001,7 +929,10 @@ void ASMu2DNastran::addBlock (int idx, ElementBlock* blk)
 ElementBlock* ASMu2DNastran::immersedGeometry (char* name) const
 {
   ElementBlock* geo = nullptr;
-  if (myMass.empty() || ASM::skipMass) return geo;
+  if (myMass.empty()) return geo;
+#ifdef HAS_FFLLIB
+  if (ASM::skipMass) return geo;
+#endif
 
   // Let the largest point mass be visualized as a sphere
   // with diameter ~5% of the total length of the structure
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7e6a40b..53de266 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -97,9 +97,16 @@ ifem_add_application(
     ShellSim
   SOURCES
     main.C
+  HEADERS
+    IFEMNastranReader.h
   LIBRARIES
     ShellEx
 )
+target_compile_definitions(ShellSim PUBLIC HAS_PROFILER=1)
+
+if(IFEM_USE_FFLLIB)
+  target_sources(ShellSim PRIVATE IFEMNastranReader.C)
+endif()
 
 # Regression tests
 enable_testing()
diff --git a/IFEMNastranReader.C b/IFEMNastranReader.C
new file mode 100644
index 0000000..09b175d
--- /dev/null
+++ b/IFEMNastranReader.C
@@ -0,0 +1,108 @@
+// $Id$
+//==============================================================================
+//!
+//! \file IFEMNastranReader.C
+//!
+//! \date Feb 27 2024
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief A Nastran Bulk Data File (BDF) parser for IFEM.
+//!
+//==============================================================================
+
+#include <cstring>
+
+#include "IFEMNastranReader.h"
+#include "FFlLinkHandler.H"
+#include "FFlElementBase.H"
+#include "FFaLib/FFaDefinitions/FFaMsg.H"
+#ifdef HAS_PROFILER
+#include "Profiler.h"
+#else
+//! \cond DO_NOT_DOCUMENT
+#define PROFILE(s)
+//! \endcond
+#endif
+
+namespace ASM {
+  char skipMass = 0; //!< 1: ignore one-noded mass elements, 2: no VTF output
+  bool skipRBE2 = false; //!< If \e true, ignore all RBE2 elements
+  bool readSets = true; //!< If \e true, read pre-bulk Nastran SET definitions
+  std::vector<int> ignoredElms; //!< List of elements to ignore
+}
+
+
+int IFEMNastranReader::parseSets (std::istream& is, std::ostream& os,
+                                  const char* fileName)
+{
+  // Fast-forward until "BEGIN BULK"
+  int lCount = 0;
+  char cline[256];
+  while (is.getline(cline,255))
+    if (!strncmp(cline,"BEGIN BULK",10))
+      break;
+    else
+    {
+      ++lCount;
+      // Copy element SET definitions to a second stream,
+      // since they have to be parsed after the FE model is loaded
+      if (ASM::readSets && (!strncmp(cline,"SET ",4) || os.tellp() > 0))
+        os << cline << '\n';
+    }
+
+  if (!is)
+  {
+    if (fileName)
+      std::cerr <<"\n *** Parsing FE data file "<< fileName
+                <<" failed."<< std::endl;
+    return -1; // No bulk data file, or second pass when parsing for more
+  }
+  else if (fileName)
+    ListUI <<"Parsing FE data file "<< fileName <<"\n";
+
+  if (lCount > 0)
+    ListUI <<"\tNastran bulk data starting at line "<< lCount+1 <<"\n";
+
+  return lCount;
+}
+
+
+bool IFEMNastranReader::readFE (std::istream& is, std::istream& iset,
+                                bool skipBeams)
+{
+  PROFILE("Nastran file parser");
+
+  if (!this->resolve(this->read(is)))
+    myLink->deleteGeometry(); // Parsing failure, delete all FE data
+  else if (nWarnings+nNotes > 0)
+    ListUI <<"\n  ** Parsing FE data succeeded."
+           <<"\n     However, "<< nWarnings
+           <<" warning(s) and "<< nNotes <<" note(s) were reported.\n"
+           <<"     Review the messages and check the FE data file.\n\n";
+  if (!myLink->hasGeometry())
+    return false;
+
+  // Remove all solid elements (not yet supported),
+  // and (optionally) all the mass and beam elements, and ...
+  ElementsVec toBeErased;
+  ElementsCIter it;
+  for (it = myLink->elementsBegin(); it != myLink->elementsEnd(); ++it)
+    if ((*it)->getCathegory() == FFlTypeInfoSpec::SOLID_ELM ||
+        std::find(ASM::ignoredElms.begin(),ASM::ignoredElms.end(),
+                  (*it)->getID()) != ASM::ignoredElms.end() ||
+        (ASM::skipRBE2 && (*it)->getTypeName() == "RGD") ||
+        (ASM::skipMass == 1 && (*it)->getTypeName() == "CMASS") ||
+        (skipBeams && (*it)->getCathegory() == FFlTypeInfoSpec::BEAM_ELM))
+        toBeErased.push_back(*it);
+  if (!toBeErased.empty())
+  {
+    ListUI <<"  ** Erasing "<< toBeErased.size()
+           <<" elements from the model.\n";
+    myLink->removeElements(toBeErased);
+  }
+
+  // Now parse the element set definitions, if any
+  lastComment = { 0, "" };
+  return iset ? this->processAllSets(iset,nPreBulk) : true;
+}
diff --git a/IFEMNastranReader.h b/IFEMNastranReader.h
new file mode 100644
index 0000000..22a8a30
--- /dev/null
+++ b/IFEMNastranReader.h
@@ -0,0 +1,41 @@
+// $Id$
+//==============================================================================
+//!
+//! \file IFEMNastranReader.h
+//!
+//! \date Feb 27 2024
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief A Nastran Bulk Data File (BDF) parser for IFEM.
+//!
+//==============================================================================
+
+#ifndef _IFEM_NASTRAN_READER_H
+#define _IFEM_NASTRAN_READER_H
+
+#include "FFlNastranReader.H"
+
+
+/*!
+  \brief Class for reading Nastran bulk data into a FFlLinkHandler object.
+*/
+
+class IFEMNastranReader : public FFlNastranReader
+{
+  int nPreBulk; //!< Number of lines before BEGIN BULK
+
+public:
+  //! \brief The constructor forwards to the parent class constructor.
+  IFEMNastranReader(FFlLinkHandler& fePart, int lCount)
+    : FFlNastranReader(&fePart,lCount), nPreBulk(lCount) {}
+
+  //! \brief Reads the FE model from the Nastran file stream.
+  bool readFE(std::istream& is, std::istream& iset, bool skipBeams = false);
+
+  //! \brief Reads through a Nastran file to extract pre-bulk SET definitions.
+  static int parseSets(std::istream& is, std::ostream& os,
+                       const char* fileName = nullptr);
+};
+
+#endif
diff --git a/Test/DroppedBox.reg b/Test/DroppedBox.reg
index 08f715c..20a2ba1 100644
--- a/Test/DroppedBox.reg
+++ b/Test/DroppedBox.reg
@@ -7,7 +7,7 @@ Parsing <discretization>
 Parsing <geometry>
   Parsing <patchfile>
 	Reading data file Q4box.nas
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          8
 Number of shell elements:       6
 Number of constraint elements:  0
diff --git a/Test/Fine-dump.reg b/Test/Fine-dump.reg
index 9bd3316..7636868 100644
--- a/Test/Fine-dump.reg
+++ b/Test/Fine-dump.reg
@@ -8,7 +8,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file FEM_Fine.nas
 	Nastran bulk data starting at line 36
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          951
 Number of shell elements:       902
 Number of constraint elements:  0
diff --git a/Test/Fine-load.reg b/Test/Fine-load.reg
index 530e6dc..fcb8a95 100644
--- a/Test/Fine-load.reg
+++ b/Test/Fine-load.reg
@@ -8,7 +8,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file FEM_Fine.nas
 	Nastran bulk data starting at line 36
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          951
 Number of shell elements:       902
 Number of constraint elements:  0
diff --git a/Test/Fine-partial.reg b/Test/Fine-partial.reg
index e260912..82882be 100644
--- a/Test/Fine-partial.reg
+++ b/Test/Fine-partial.reg
@@ -8,7 +8,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file FEM_Fine.nas
 	Nastran bulk data starting at line 36
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          951
 Number of shell elements:       902
 Number of constraint elements:  0
diff --git a/Test/Q4-dyn.reg b/Test/Q4-dyn.reg
index 4a29880..358dbdc 100644
--- a/Test/Q4-dyn.reg
+++ b/Test/Q4-dyn.reg
@@ -8,7 +8,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file CQUAD04_.nas
 	Nastran bulk data starting at line 31
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          27
 Number of shell elements:       16
 Number of constraint elements:  2
diff --git a/Test/Q4-dynp.reg b/Test/Q4-dynp.reg
index e9343d3..8821b56 100644
--- a/Test/Q4-dynp.reg
+++ b/Test/Q4-dynp.reg
@@ -8,7 +8,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file CQUAD04_.nas
 	Nastran bulk data starting at line 31
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          27
 Number of shell elements:       16
 Number of constraint elements:  2
diff --git a/Test/Robaat-irr.reg b/Test/Robaat-irr.reg
index 850d083..746eea4 100644
--- a/Test/Robaat-irr.reg
+++ b/Test/Robaat-irr.reg
@@ -9,7 +9,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file Small_Boat.dat
 	Nastran bulk data starting at line 71
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          518
 Number of shell elements:       504
 Number of constraint elements:  0
diff --git a/Test/Robaat.reg b/Test/Robaat.reg
index 7caacd5..ac5e9b3 100644
--- a/Test/Robaat.reg
+++ b/Test/Robaat.reg
@@ -9,7 +9,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file Small_Boat.dat
 	Nastran bulk data starting at line 71
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          518
 Number of shell elements:       504
 Number of constraint elements:  0
diff --git a/Test/Wheel.reg b/Test/Wheel.reg
index 94fff60..f2028ef 100644
--- a/Test/Wheel.reg
+++ b/Test/Wheel.reg
@@ -8,7 +8,7 @@ Parsing <geometry>
   Parsing <patchfile>
 	Reading data file sykkelhjul.nas
 	Nastran bulk data starting at line 2
-Parsing Nastran bulk data file succceeded.
+Parsing Nastran bulk data file succeeded.
 Total number of nodes:          31
 Number of shell elements:       20
 Number of constraint elements:  1
diff --git a/main.C b/main.C
index 5d7b227..321116a 100644
--- a/main.C
+++ b/main.C
@@ -36,10 +36,13 @@
 #include <unistd.h>
 
 namespace ASM {
-  extern char useBeam, skipMass;
-  extern bool skipRBE2, replRBE3, skipLoad, readSets;
+  extern bool replRBE3, skipLoad;
   extern double Ktra, Krot;
+#ifdef HAS_FFLLIB
+  extern bool skipRBE2, readSets;
+  extern char useBeam, skipMass;
   extern std::vector<int> ignoredElms;
+#endif
 }
 
 
@@ -138,11 +141,10 @@ int main (int argc, char** argv)
       iop = 100;
     else if (!strcmp(argv[i],"-ignoreSol"))
       iop = 200;
+#ifdef HAS_FFLLIB
     else if (!strcmp(argv[i],"-ignore"))
       while (i < argc-1 && isdigit(argv[i+1][0]))
         utl::parseIntegers(ASM::ignoredElms,argv[++i]);
-    else if (!strcmp(argv[i],"-vizRHS"))
-      vizRHS = true;
     else if (!strcmp(argv[i],"-noSets"))
       ASM::readSets = false;
     else if (!strcmp(argv[i],"-noMasses"))
@@ -153,6 +155,7 @@ int main (int argc, char** argv)
       ASM::useBeam = 2;
     else if (!strcmp(argv[i],"-noRBE2"))
       ASM::skipRBE2 = true;
+#endif
     else if (!strcmp(argv[i],"-noRBE3"))
       ASM::replRBE3 = true;
     else if (!strcmp(argv[i],"-noLoad"))
@@ -164,6 +167,8 @@ int main (int argc, char** argv)
       if (i+1 < argc && argv[i+1][0] != '-')
         ASM::Krot = atof(argv[++i]);
     }
+    else if (!strcmp(argv[i],"-vizRHS"))
+      vizRHS = true;
     else if (!strcmp(argv[i],"-fixDup"))
     {
       fixDup = true;
@@ -186,8 +191,10 @@ int main (int argc, char** argv)
     }
     else if (!strcmp(argv[i],"-keep-previous-state"))
       nstates = 2; // both current and previous state will reside in core
+#ifdef HAS_FFLLIB
     else if (!strcmp(argv[i],"-no-vtfmass"))
       ASM::skipMass = 2;
+#endif
     else if (!strncmp(argv[i],"-vtfres",6))
     {
       while (i+1 < argc && argv[i+1][0] != '-')

From 705f7153dffaf51d15d449c9327bbc2cdc7cf395 Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Wed, 18 Mar 2026 13:25:23 +0100
Subject: [PATCH 2/6] Use std::array as internal container for FFaTensor2 and
 FFaTensor3. Add method fill() for the tensor classes.

---
 FFlLib/FFaLib/FFaAlgebra/FFaTensor1.C |  1 -
 FFlLib/FFaLib/FFaAlgebra/FFaTensor1.H |  2 ++
 FFlLib/FFaLib/FFaAlgebra/FFaTensor2.C | 11 +++++------
 FFlLib/FFaLib/FFaAlgebra/FFaTensor2.H | 23 ++++++++++++-----------
 FFlLib/FFaLib/FFaAlgebra/FFaTensor3.C | 23 ++++++++++++-----------
 FFlLib/FFaLib/FFaAlgebra/FFaTensor3.H | 23 ++++++++++++-----------
 6 files changed, 43 insertions(+), 40 deletions(-)

diff --git a/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.C b/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.C
index c7d41a6..5366476 100644
--- a/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.C
+++ b/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.C
@@ -8,7 +8,6 @@
 #include "FFaLib/FFaAlgebra/FFaTensor1.H"
 #include "FFaLib/FFaAlgebra/FFaTensor2.H"
 #include "FFaLib/FFaAlgebra/FFaTensor3.H"
-#include <cctype>
 #include <cmath>
 
 
diff --git a/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.H b/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.H
index 8a54e88..a6e7deb 100644
--- a/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.H
+++ b/FFlLib/FFaLib/FFaAlgebra/FFaTensor1.H
@@ -33,6 +33,8 @@ public:
   FFaTensor1(const FFaTensor2& t);
   FFaTensor1(const FFaTensor3& t);
 
+  void fill(double d) { myT = d; }
+
   // Local operators
 
   FFaTensor1& operator= (const FFaTensor1& t);
diff --git a/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.C b/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.C
index c945907..ac48335 100644
--- a/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.C
+++ b/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.C
@@ -11,7 +11,6 @@
 #include "FFaLib/FFaAlgebra/FFaMat33.H"
 #include "FFaLib/FFaAlgebra/FFaMat34.H"
 #include "FFaLib/FFaAlgebra/FFaTensorTransforms.H"
-#include <cctype>
 
 
 FFaTensor2::FFaTensor2(const FFaTensor3& t)
@@ -64,7 +63,7 @@ FFaTensor2& FFaTensor2::operator= (const FFaTensor1& t)
 
 FFaTensor2& FFaTensor2::rotate(const double ex[2], const double ey[2])
 {
-  FFaTensorTransforms::rotate(myT,ex,ey, myT);
+  FFaTensorTransforms::rotate(myT.data(), ex,ey, myT.data());
   return *this;
 }
 
@@ -102,7 +101,7 @@ double FFaTensor2::maxShear() const
 void FFaTensor2::maxShear(FaVec3& v) const
 {
   double values[2], max[2], min[2];
-  if (FFaTensorTransforms::principalDirs(myT,values,max,min) == 0)
+  if (FFaTensorTransforms::principalDirs(myT.data(),values,max,min) == 0)
   {
     v[2] = 0.0;
     FFaTensorTransforms::maxShearDir(2,max,min,v.getPt());
@@ -166,7 +165,7 @@ void FFaTensor2::prinsipalValues(FaVec3& values, FaMat33& rotation) const
 {
   values[2] = 0.0;
   rotation.setIdentity();
-  if (FFaTensorTransforms::principalDirs(myT,
+  if (FFaTensorTransforms::principalDirs(myT.data(),
                                          values.getPt(),
                                          rotation[0].getPt(),
                                          rotation[1].getPt()))
@@ -245,7 +244,7 @@ FFaTensor3 operator* (const FFaTensor2& a, const FaMat33 & m)
   FFaTensor3 result;
   FFaTensor3 in(a);
   FFaTensorTransforms::rotate(in.getPt(),
-			      m[0].getPt(), m[1].getPt(), m[2].getPt(),
+                              m[0].getPt(), m[1].getPt(), m[2].getPt(),
                               result.getPt());
   return result;
 }
@@ -255,7 +254,7 @@ FFaTensor3 operator* (const FFaTensor2& a, const FaMat34& m)
   FFaTensor3 result;
   FFaTensor3 in(a);
   FFaTensorTransforms::rotate(in.getPt(),
-			      m[0].getPt(), m[1].getPt(), m[2].getPt(),
+                              m[0].getPt(), m[1].getPt(), m[2].getPt(),
                               result.getPt());
   return result;
 }
diff --git a/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.H b/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.H
index aa95aa2..9d08a3e 100644
--- a/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.H
+++ b/FFlLib/FFaLib/FFaAlgebra/FFaTensor2.H
@@ -8,6 +8,7 @@
 #ifndef FFA_TENSOR2_H
 #define FFA_TENSOR2_H
 
+#include <array>
 #include <iostream>
 
 class FaVec3;
@@ -25,21 +26,21 @@ class FFaTensor3;
 
 class FFaTensor2
 {
-  double myT[3];
+  std::array<double,3> myT;
 
 public:
 
   // Constructors
 
-  FFaTensor2()                    { myT[0] = myT[1] = myT[2] = 0.0; }
-  FFaTensor2(double d)            { myT[0] = myT[1] = myT[2] = d; }
+  FFaTensor2(double d = 0.0)      { myT.fill(d); }
   FFaTensor2(const float* t)      { for (int i=0; i<3; i++) myT[i] = t[i]; }
   FFaTensor2(const double* t)     { for (int i=0; i<3; i++) myT[i] = t[i]; }
   FFaTensor2(const FFaTensor2& t) { for (int i=0; i<3; i++) myT[i] = t.myT[i]; }
   FFaTensor2(const FFaTensor1& t);
   FFaTensor2(const FFaTensor3& t);
-  FFaTensor2(double t11, double t22, double t12 = 0.0)
-  { myT[0] = t11; myT[1] = t22; myT[2] = t12; }
+  FFaTensor2(double t11, double t22, double t12 = 0.0) { myT = {t11,t22,t12}; }
+
+  void fill(double d) { myT.fill(d); }
 
   // Local operators
 
@@ -73,8 +74,8 @@ public:
 
   // Access to internal array
 
-  const double* getPt() const { return myT; }
-  double*       getPt()       { return myT; }
+  const double* getPt() const { return myT.data(); }
+  double*       getPt()       { return myT.data(); }
 
   // Indexing
 
@@ -122,8 +123,8 @@ inline const double& FFaTensor2::operator[] (int i) const
 {
 #ifdef FFA_INDEXCHECK
   if (i < 0 || i > 2)
-    std::cerr <<"FFaTensor2::operator[]: index i="<< i <<" is out of range [0,2]"
-	      << std::endl;
+    std::cerr <<"FFaTensor2::operator[]: index i="<< i
+              <<" is out of range [0,2]"<< std::endl;
 #endif
   return myT[i];
 }
@@ -132,8 +133,8 @@ inline double& FFaTensor2::operator[] (int i)
 {
 #ifdef FFA_INDEXCHECK
   if (i < 0 || i > 2)
-    std::cerr <<"FFaTensor2::operator[]: index i="<< i <<" is out of range [0,2]"
-	      << std::endl;
+    std::cerr <<"FFaTensor2::operator[]: index i="<< i
+              <<" is out of range [0,2]"<< std::endl;
 #endif
   return myT[i];
 }
diff --git a/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.C b/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.C
index e0fa11b..a0e8ec6 100644
--- a/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.C
+++ b/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.C
@@ -11,7 +11,6 @@
 #include "FFaLib/FFaAlgebra/FFaMat33.H"
 #include "FFaLib/FFaAlgebra/FFaMat34.H"
 #include "FFaLib/FFaAlgebra/FFaTensorTransforms.H"
-#include <cctype>
 
 
 FFaTensor3::FFaTensor3(const FFaTensor2& t)
@@ -71,22 +70,22 @@ FFaTensor3& FFaTensor3::operator= (const FFaTensor1& t)
 
 FFaTensor3& FFaTensor3::rotate(const FaMat33& rotMx)
 {
-  FFaTensorTransforms::rotate(myT,
+  FFaTensorTransforms::rotate(myT.data(),
                               rotMx[0].getPt(),
                               rotMx[1].getPt(),
                               rotMx[2].getPt(),
-                              myT);
+                              myT.data());
   return *this;
 }
 
 
 FFaTensor3& FFaTensor3::rotate(const FaMat34& rotMx)
 {
-  FFaTensorTransforms::rotate(myT,
+  FFaTensorTransforms::rotate(myT.data(),
                               rotMx[0].getPt(),
                               rotMx[1].getPt(),
                               rotMx[2].getPt(),
-                              myT);
+                              myT.data());
   return *this;
 }
 
@@ -174,7 +173,7 @@ double FFaTensor3::maxShear() const
 void FFaTensor3::maxShear(FaVec3& v) const
 {
   double values[3], max[3], middle[3], min[3];
-  if (FFaTensorTransforms::principalDirs(myT,values,max,middle,min) == 0)
+  if (FFaTensorTransforms::principalDirs(myT.data(),values,max,middle,min) == 0)
   {
     FFaTensorTransforms::maxShearDir(3,max,min,v.getPt());
     v *= FFaTensorTransforms::maxShearValue(values[0],values[2]);
@@ -258,7 +257,7 @@ void FFaTensor3::prinsipalValues(double& max, double& middle, double& min) const
 
 void FFaTensor3::prinsipalValues(FaVec3& values, FaMat33& rotation) const
 {
-  if (FFaTensorTransforms::principalDirs(myT,
+  if (FFaTensorTransforms::principalDirs(myT.data(),
                                          values.getPt(),
                                          rotation[0].getPt(),
                                          rotation[1].getPt(),
@@ -355,16 +354,18 @@ bool operator!= (const FFaTensor3& a, const FFaTensor3& b)
 FFaTensor3 operator* (const FFaTensor3& a, const FaMat33& m)
 {
   FFaTensor3 result;
-  FFaTensorTransforms::rotate(a.myT, m[0].getPt(), m[1].getPt(), m[2].getPt(),
-                              result.myT);
+  FFaTensorTransforms::rotate(a.myT.data(),
+                              m[0].getPt(), m[1].getPt(), m[2].getPt(),
+                              result.myT.data());
   return result;
 }
 
 FFaTensor3 operator* (const FFaTensor3& a, const FaMat34& m)
 {
   FFaTensor3 result;
-  FFaTensorTransforms::rotate(a.myT, m[0].getPt(), m[1].getPt(), m[2].getPt(),
-                              result.myT);
+  FFaTensorTransforms::rotate(a.myT.data(),
+                              m[0].getPt(), m[1].getPt(), m[2].getPt(),
+                              result.myT.data());
   return result;
 }
 
diff --git a/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.H b/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.H
index eb14673..04908fe 100644
--- a/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.H
+++ b/FFlLib/FFaLib/FFaAlgebra/FFaTensor3.H
@@ -8,6 +8,7 @@
 #ifndef FFA_TENSOR3_H
 #define FFA_TENSOR3_H
 
+#include <array>
 #include <iostream>
 
 class FaVec3;
@@ -25,14 +26,13 @@ class FFaTensor2;
 
 class FFaTensor3
 {
-  double myT[6];
+  std::array<double,6> myT;
 
 public:
 
   // Constructors
 
-  FFaTensor3()                    { for (int i=0; i<6; i++) myT[i] = 0.0; }
-  FFaTensor3(double d)            { for (int i=0; i<6; i++) myT[i] = d; }
+  FFaTensor3(double d = 0.0)      { myT.fill(d); }
   FFaTensor3(const float* t)      { for (int i=0; i<6; i++) myT[i] = t[i]; }
   FFaTensor3(const double* t)     { for (int i=0; i<6; i++) myT[i] = t[i]; }
   FFaTensor3(const FFaTensor3& t) { for (int i=0; i<6; i++) myT[i] = t.myT[i]; }
@@ -42,14 +42,15 @@ public:
   FFaTensor3(double t11      , double t22      , double t33,
              double t12 = 0.0, double t13 = 0.0, double t23 = 0.0)
   {
-    myT[0] = t11; myT[1] = t22; myT[2] = t33;
-    myT[3] = t12; myT[4] = t13; myT[5] = t23;
+    myT = { t11, t22, t33, t12, t13, t23 };
   }
   FFaTensor3(const FaVec3& v1, const FaVec3& v2, const FaVec3& v3)
   {
     this->makeInertia(v1,v2,v3);
   }
 
+  void fill(double d) { myT.fill(d); }
+
   // Local operators
 
   FFaTensor3& operator= (const FFaTensor3& t);
@@ -82,8 +83,8 @@ public:
 
   // Access to internal array
 
-  const double* getPt() const { return myT; }
-  double*       getPt()       { return myT; }
+  const double* getPt() const { return myT.data(); }
+  double*       getPt()       { return myT.data(); }
 
   // Indexing
 
@@ -136,8 +137,8 @@ inline const double& FFaTensor3::operator[] (int i) const
 {
 #ifdef FFA_INDEXCHECK
   if (i < 0 || i > 5)
-    std::cerr <<"FFaTensor3::operator[]: index i="<< i <<" is out of range [0,5]"
-	      << std::endl;
+    std::cerr <<"FFaTensor3::operator[]: index i="<< i
+              <<" is out of range [0,5]"<< std::endl;
 #endif
   return myT[i];
 }
@@ -146,8 +147,8 @@ inline double& FFaTensor3::operator[] (int i)
 {
 #ifdef FFA_INDEXCHECK
   if (i < 0 || i > 5)
-    std::cerr <<"FFaTensor3::operator[]: index i="<< i <<" is out of range [0,5]"
-	      << std::endl;
+    std::cerr <<"FFaTensor3::operator[]: index i="<< i
+              <<" is out of range [0,5]"<< std::endl;
 #endif
   return myT[i];
 }

From fa9dc7a561723433a6212f62189b4e8cda34619e Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Wed, 18 Mar 2026 14:25:15 +0100
Subject: [PATCH 3/6] Replace virtual method getVolumeAndInertia() by
 getVolumeAndCoG() where the inertia calculation is optional

---
 FFlLib/FFlElementBase.C       |  15 ++--
 FFlLib/FFlElementBase.H       |   3 +-
 FFlLib/FFlFEParts/FFlBeams.C  |  42 +++++-----
 FFlLib/FFlFEParts/FFlBeams.H  |   8 +-
 FFlLib/FFlFEParts/FFlCMASS.C  |  35 ++++----
 FFlLib/FFlFEParts/FFlCMASS.H  |   4 +-
 FFlLib/FFlFEParts/FFlShells.C | 145 ++++++++++++++++++----------------
 FFlLib/FFlFEParts/FFlShells.H |  17 +---
 FFlLib/FFlFEParts/FFlSolids.C |  54 +++++++------
 FFlLib/FFlFEParts/FFlSolids.H |  24 ++----
 10 files changed, 168 insertions(+), 179 deletions(-)

diff --git a/FFlLib/FFlElementBase.C b/FFlLib/FFlElementBase.C
index 8b7fa97..55d81d6 100644
--- a/FFlLib/FFlElementBase.C
+++ b/FFlLib/FFlElementBase.C
@@ -81,20 +81,21 @@ double FFlElementBase::getMassDensity() const
 }
 
 
-bool FFlElementBase::getVolumeAndInertia(double& volume, FaVec3& cog,
-					 FFaTensor3& inertia) const
+double FFlElementBase::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inert) const
 {
-  volume  = 0.0;
-  cog     = FaVec3(0.0,0.0,0.0);
-  inertia = FFaTensor3(0.0);
-  return false;
+  cog.clear();
+  if (inert)
+    inert->fill(0.0);
+
+  return 0.0;
 }
 
 
 bool FFlElementBase::getMassProperties(double& mass, FaVec3& cog,
 				       FFaTensor3& inertia) const
 {
-  if (!this->getVolumeAndInertia(mass,cog,inertia)) return false;
+  if ((mass = this->getVolumeAndCoG(cog,&inertia)) <= 0.0)
+    return false;
 
   double rho = this->getMassDensity();
   mass    *= rho;
diff --git a/FFlLib/FFlElementBase.H b/FFlLib/FFlElementBase.H
index 05d1168..955ce8f 100644
--- a/FFlLib/FFlElementBase.H
+++ b/FFlLib/FFlElementBase.H
@@ -73,8 +73,7 @@ public:
 
   // Structural properties:
   virtual double getMassDensity() const;
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-				   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inert = NULL) const;
   bool getMassProperties(double& mass, FaVec3& cog, FFaTensor3& inertia) const;
 
   // Geometrical mapping:
diff --git a/FFlLib/FFlFEParts/FFlBeams.C b/FFlLib/FFlFEParts/FFlBeams.C
index b632d76..00ed99d 100644
--- a/FFlLib/FFlFEParts/FFlBeams.C
+++ b/FFlLib/FFlFEParts/FFlBeams.C
@@ -71,11 +71,10 @@ double FFlBEAM2::getMassDensity() const
 }
 
 
-bool FFlBEAM2::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlBEAM2::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   FFlPBEAMSECTION* psec = dynamic_cast<FFlPBEAMSECTION*>(this->getAttribute("PBEAMSECTION"));
-  if (!psec) return false; // Should not happen
+  if (!psec) return -1.0; // Should not happen
 
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
@@ -90,8 +89,9 @@ bool FFlBEAM2::getVolumeAndInertia(double& volume, FaVec3& cog,
 
   double A = psec->crossSectionArea.getValue();
   double L = (v2-v1).length();
+  double V = A*L;
   cog      = (v1+v2)*0.5;
-  volume   = A*L;
+  if (!inertia) return V;
 
   // Compute local coordinate system
   FaMat33 Telm;
@@ -109,14 +109,15 @@ bool FFlBEAM2::getVolumeAndInertia(double& volume, FaVec3& cog,
 
   // Set up the inertia tensor in local axes
   double Ixx = psec->Iy.getValue() + psec->Iz.getValue();
-  inertia[0] = L*(Ixx > 0.0 ? Ixx : psec->It.getValue());
-  inertia[1] = L*(psec->Iy.getValue() + A*L*L/12.0);
-  inertia[2] = L*(psec->Iz.getValue() + A*L*L/12.0);
-  inertia[3] = inertia[4] = inertia[5] = 0.0;
+  FFaTensor3& I = *inertia;
+  I[0] = L*(Ixx > 0.0 ? Ixx : psec->It.getValue());
+  I[1] = L*(psec->Iy.getValue() + A*L*L/12.0);
+  I[2] = L*(psec->Iz.getValue() + A*L*L/12.0);
+  I[3] = I[4] = I[5] = 0.0;
 
   // Transform to global axes
-  inertia.rotate(Telm.transpose());
-  return true;
+  inertia->rotate(Telm.transpose());
+  return V;
 }
 
 
@@ -312,11 +313,10 @@ double FFlBEAM3::getMassDensity() const
 }
 
 
-bool FFlBEAM3::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlBEAM3::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   FFlPBEAMSECTION* psec = dynamic_cast<FFlPBEAMSECTION*>(this->getAttribute("PBEAMSECTION"));
-  if (!psec) return false; // Should not happen
+  if (!psec) return -1.0; // Should not happen
 
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
@@ -333,8 +333,9 @@ bool FFlBEAM3::getVolumeAndInertia(double& volume, FaVec3& cog,
 
   double A = psec->crossSectionArea.getValue();
   double L = (v2-v1).length() + (v3-v2).length();
+  double V = A*L;
   cog      = (v1+v3)*0.25 + v2*0.5;
-  volume   = A*L;
+  if (!inertia) return V;
 
   // Compute local coordinate system
   FFlPORIENT* pori = dynamic_cast<FFlPORIENT*>(this->getAttribute("PORIENT"));
@@ -355,14 +356,15 @@ bool FFlBEAM3::getVolumeAndInertia(double& volume, FaVec3& cog,
 
   // Set up the inertia tensor in local axes
   double Ixx = psec->Iy.getValue() + psec->Iz.getValue();
-  inertia[0] = L*(Ixx > 0.0 ? Ixx : psec->It.getValue());
-  inertia[1] = L*(psec->Iy.getValue() + A*L*L/12.0);
-  inertia[2] = L*(psec->Iz.getValue() + A*L*L/12.0);
-  inertia[3] = inertia[4] = inertia[5] = 0.0;
+  FFaTensor3& I = *inertia;
+  I[0] = L*(Ixx > 0.0 ? Ixx : psec->It.getValue());
+  I[1] = L*(psec->Iy.getValue() + A*L*L/12.0);
+  I[2] = L*(psec->Iz.getValue() + A*L*L/12.0);
+  I[3] = I[4] = I[5] = 0.0;
 
   // Transform to global axes
-  inertia.rotate(Telm.transpose());
-  return true;
+  inertia->rotate(Telm.transpose());
+  return V;
 }
 
 #ifdef FF_NAMESPACE
diff --git a/FFlLib/FFlFEParts/FFlBeams.H b/FFlLib/FFlFEParts/FFlBeams.H
index 569b151..bb83110 100644
--- a/FFlLib/FFlFEParts/FFlBeams.H
+++ b/FFlLib/FFlFEParts/FFlBeams.H
@@ -25,7 +25,6 @@ class FFlBEAM2 : public FFlElementBase
 {
 public:
   FFlBEAM2(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlBEAM2() {}
 
   static void init();
 
@@ -35,8 +34,7 @@ public:
   FFL_TYPE_INFO(FFlBEAM2);
 
   virtual double getMassDensity() const;
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
   virtual int getNodalCoor(double* X, double* Y, double* Z) const;
   FaVec3 getLocalZdirection() const;
 
@@ -50,7 +48,6 @@ class FFlBEAM3 : public FFlElementBase
 {
 public:
   FFlBEAM3(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlBEAM3() {}
 
   static void init();
 
@@ -60,8 +57,7 @@ public:
   FFL_TYPE_INFO(FFlBEAM3);
 
   virtual double getMassDensity() const;
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
   //! \brief Splits the 3-noded beam element into two 2-noded elements.
   virtual bool split(Elements& newElm, FFlLinkHandler* owner, int);
diff --git a/FFlLib/FFlFEParts/FFlCMASS.C b/FFlLib/FFlFEParts/FFlCMASS.C
index c8d5e98..497312d 100644
--- a/FFlLib/FFlFEParts/FFlCMASS.C
+++ b/FFlLib/FFlFEParts/FFlCMASS.C
@@ -37,27 +37,26 @@ void FFlCMASS::init()
 }
 
 
-bool FFlCMASS::getVolumeAndInertia(double& volume, FaVec3& cog,
-				   FFaTensor3& inertia) const
+double FFlCMASS::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
-  volume = 0.0;
   cog = this->getNode(1)->getPos();
 
   FFlPMASS* pmass = dynamic_cast<FFlPMASS*>(this->getAttribute("PMASS"));
   if (!pmass)
   {
     // Property-less mass element, ignore
-    inertia = FFaTensor3(0.0);
-    return true;
+    if (inertia)
+      inertia->fill(0.0);
+    return 0.0;
   }
 
   const std::vector<double>& M = pmass->M.getValue();
   if (M.size() > 0 && M.size() <= 6)
   {
     // No inertia for this element, only mass
-    volume = M[0];
-    inertia = FFaTensor3(0.0);
-    return true;
+    if (inertia)
+      inertia->fill(0.0);
+    return M[0];
   }
   else if (M.size() != 21)
   {
@@ -82,21 +81,23 @@ bool FFlCMASS::getVolumeAndInertia(double& volume, FaVec3& cog,
   bool haveInertia = RtMR(M,R,II);
 
   // Compute the mass = volume
-  volume = (II[0][0]+II[1][1]+II[2][2]) / 3.0;
+  double V = (II[0][0]+II[1][1]+II[2][2]) / 3.0;
   if (!haveInertia)
   {
-    inertia = FFaTensor3(0.0);
-    return true;
+    if (inertia)
+      inertia->fill(0.0);
+    return V;
   }
-  else if (volume < 1.0e-16)
+  else if (V < 1.0e-16)
   {
     ListUI <<"\n  ** Warning: CMASS element "<< this->getID()
            <<" has zero mass, but non-zero inertia. This is non-physical.\n";
-    return false;
+    return -2.0;
   }
 
   // Adjust the center of gravity for possible offset
-  cog = FaVec3(II[1][5],II[2][3],II[0][4]) / volume;
+  cog = FaVec3(II[1][5],II[2][3],II[0][4]) / V;
+  if (!inertia) return V;
 
   // Compute the inertia tensor about the center of gravity
   R[1][3] += cog[2];
@@ -107,9 +108,9 @@ bool FFlCMASS::getVolumeAndInertia(double& volume, FaVec3& cog,
   R[1][5] -= cog[0];
   RtMR(M,R,II);
 
-  inertia = FFaTensor3(II[3][3],II[4][4],II[5][5],
-                       II[3][4],II[3][5],II[4][5]);
-  return true;
+  *inertia = FFaTensor3(II[3][3],II[4][4],II[5][5],
+                        II[3][4],II[3][5],II[4][5]);
+  return V;
 }
 
 
diff --git a/FFlLib/FFlFEParts/FFlCMASS.H b/FFlLib/FFlFEParts/FFlCMASS.H
index 00055ce..d77f7a1 100644
--- a/FFlLib/FFlFEParts/FFlCMASS.H
+++ b/FFlLib/FFlFEParts/FFlCMASS.H
@@ -25,7 +25,6 @@ class FFlCMASS : public FFlElementBase
 {
 public:
   FFlCMASS(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlCMASS() {}
 
   static void init();
 
@@ -35,8 +34,7 @@ public:
   FFL_TYPE_INFO(FFlCMASS);
 
   virtual double getMassDensity() const { return 1.0; }
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-				   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
 private:
   static bool RtMR(const std::vector<double>& M, double R[6][6], double II[6][6]);
diff --git a/FFlLib/FFlFEParts/FFlShells.C b/FFlLib/FFlFEParts/FFlShells.C
index 5edf027..9936b1f 100644
--- a/FFlLib/FFlFEParts/FFlShells.C
+++ b/FFlLib/FFlFEParts/FFlShells.C
@@ -109,35 +109,37 @@ bool FFlTRI3::getFaceNormals(std::vector<FaVec3>& normals, short int,
 }
 
 
-bool FFlTRI3::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                  FFaTensor3& inertia) const
+double FFlTRI3::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   double th = this->getThickness();
-  if (th <= 0.0) return false; // Should not happen
+  if (th <= 0.0) return -1.0; // Should not happen
 
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
   FaVec3 v3(this->getNode(3)->getPos());
+  cog = (v1 + v2 + v3) / 3.0;
 
   FaVec3 normal = (v2-v1) ^ (v3-v1);
   double length = normal.length();
-  volume = 0.5 * length * th;
-  cog = (v1 + v2 + v3) / 3.0;
-  if (length < 1.0e-16) return false;
+  if (length < 1.0e-16) return -2.0;
 
   // Compute inertias by expanding the shell into a solid
 
-  normal *= th/length;
-  v1 -= cog + 0.5*normal;
-  v2 -= cog + 0.5*normal;
-  v3 -= cog + 0.5*normal;
+  if (inertia)
+  {
+    normal *= th/length;
+    v1 -= cog + 0.5*normal;
+    v2 -= cog + 0.5*normal;
+    v3 -= cog + 0.5*normal;
+
+    FaVec3 v4(v1+normal);
+    FaVec3 v5(v2+normal);
+    FaVec3 v6(v3+normal);
 
-  FaVec3 v4(v1+normal);
-  FaVec3 v5(v2+normal);
-  FaVec3 v6(v3+normal);
+    FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,*inertia);
+  }
 
-  FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,inertia);
-  return true;
+  return 0.5 * length * th; // volume
 }
 
 
@@ -295,11 +297,10 @@ bool FFlTRI6::getFaceNormals(std::vector<FaVec3>& normals, short int,
 }
 
 
-bool FFlTRI6::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                  FFaTensor3& inertia) const
+double FFlTRI6::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   double th = this->getThickness();
-  if (th <= 0.0) return false; // Should not happen
+  if (th <= 0.0) return -1.0; // Should not happen
 
   //TODO,kmo: Account for curved edges and face (by numerical integration)
   FaVec3 v1(this->getNode(1)->getPos());
@@ -308,23 +309,25 @@ bool FFlTRI6::getVolumeAndInertia(double& volume, FaVec3& cog,
 
   FaVec3 normal = (v2-v1) ^ (v3-v1);
   double length = normal.length();
-  volume = 0.5 * length * th;
   cog = (v1 + v2 + v3) / 3.0;
-  if (length < 1.0e-16) return false;
+  if (length < 1.0e-16) return -2.0;
 
   // Compute inertias by expanding the shell into a solid
 
-  normal *= th/length;
-  v1 -= cog + 0.5*normal;
-  v2 -= cog + 0.5*normal;
-  v3 -= cog + 0.5*normal;
+  if (inertia)
+  {
+    normal *= th/length;
+    v1 -= cog + 0.5*normal;
+    v2 -= cog + 0.5*normal;
+    v3 -= cog + 0.5*normal;
 
-  FaVec3 v4(v1+normal);
-  FaVec3 v5(v2+normal);
-  FaVec3 v6(v3+normal);
+    FaVec3 v4(v1+normal);
+    FaVec3 v5(v2+normal);
+    FaVec3 v6(v3+normal);
 
-  FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,inertia);
-  return true;
+    FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,*inertia);
+  }
+  return 0.5 * length * th; // volume
 }
 
 
@@ -401,11 +404,10 @@ bool FFlQUAD4::getFaceNormals(std::vector<FaVec3>& normals, short int,
 }
 
 
-bool FFlQUAD4::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlQUAD4::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   double th = this->getThickness();
-  if (th <= 0.0) return false; // Should not happen
+  if (th <= 0.0) return -1.0; // Should not happen
 
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
@@ -415,28 +417,31 @@ bool FFlQUAD4::getVolumeAndInertia(double& volume, FaVec3& cog,
   double a1 = ((v2-v1) ^ (v3-v1)).length();
   double a2 = ((v3-v1) ^ (v4-v1)).length();
 
-  volume = 0.5*(a1+a2)*th;
   cog = ((v1+v3)*(a1+a2) + v2*a1 + v4*a2) / ((a1+a2)*3.0);
 
   // Compute inertias by expanding the shell into a solid
 
-  FaVec3 normal = (v3-v1) ^ (v4-v2);
-  double length = normal.length();
-  if (length < 1.0e-16) return false;
-
-  normal *= th/length;
-  v1 -= cog + 0.5*normal;
-  v2 -= cog + 0.5*normal;
-  v3 -= cog + 0.5*normal;
-  v4 -= cog + 0.5*normal;
-
-  FaVec3 v5(v1+normal);
-  FaVec3 v6(v2+normal);
-  FaVec3 v7(v3+normal);
-  FaVec3 v8(v4+normal);
+  if (inertia)
+  {
+    FaVec3 normal = (v3-v1) ^ (v4-v2);
+    double length = normal.length();
+    if (length < 1.0e-16) return -2.0;
+
+    normal *= th/length;
+    v1 -= cog + 0.5*normal;
+    v2 -= cog + 0.5*normal;
+    v3 -= cog + 0.5*normal;
+    v4 -= cog + 0.5*normal;
+
+    FaVec3 v5(v1+normal);
+    FaVec3 v6(v2+normal);
+    FaVec3 v7(v3+normal);
+    FaVec3 v8(v4+normal);
+
+    FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,*inertia);
+  }
 
-  FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,inertia);
-  return true;
+  return 0.5*(a1+a2)*th; // volume
 }
 
 
@@ -663,11 +668,10 @@ bool FFlQUAD8::getFaceNormals(std::vector<FaVec3>& normals, short int,
 }
 
 
-bool FFlQUAD8::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlQUAD8::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   double th = this->getThickness();
-  if (th <= 0.0) return false; // Should not happen
+  if (th <= 0.0) return -1.0; // Should not happen
 
   //TODO,kmo: Account for curved edges and face (by numerical integration)
   FaVec3 v1(this->getNode(1)->getPos());
@@ -678,28 +682,31 @@ bool FFlQUAD8::getVolumeAndInertia(double& volume, FaVec3& cog,
   double a1 = ((v2-v1) ^ (v3-v1)).length();
   double a2 = ((v3-v1) ^ (v4-v1)).length();
 
-  volume = 0.5*(a1+a2)*th;
   cog = ((v1+v3)*(a1+a2) + v2*a1 + v4*a2) / ((a1+a2)*3.0);
 
   // Compute inertias by expanding the shell into a solid
 
-  FaVec3 normal = (v3-v1) ^ (v4-v2);
-  double length = normal.length();
-  if (length < 1.0e-16) return false;
-
-  normal *= th/length;
-  v1 -= cog + 0.5*normal;
-  v2 -= cog + 0.5*normal;
-  v3 -= cog + 0.5*normal;
-  v4 -= cog + 0.5*normal;
-
-  FaVec3 v5(v1+normal);
-  FaVec3 v6(v2+normal);
-  FaVec3 v7(v3+normal);
-  FaVec3 v8(v4+normal);
+  if (inertia)
+  {
+    FaVec3 normal = (v3-v1) ^ (v4-v2);
+    double length = normal.length();
+    if (length < 1.0e-16) return -2.0;
+
+    normal *= th/length;
+    v1 -= cog + 0.5*normal;
+    v2 -= cog + 0.5*normal;
+    v3 -= cog + 0.5*normal;
+    v4 -= cog + 0.5*normal;
+
+    FaVec3 v5(v1+normal);
+    FaVec3 v6(v2+normal);
+    FaVec3 v7(v3+normal);
+    FaVec3 v8(v4+normal);
+
+    FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,*inertia);
+  }
 
-  FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,inertia);
-  return true;
+  return 0.5*(a1+a2)*th; // volume
 }
 
 #ifdef FF_NAMESPACE
diff --git a/FFlLib/FFlFEParts/FFlShells.H b/FFlLib/FFlFEParts/FFlShells.H
index 32d0105..500c1ac 100644
--- a/FFlLib/FFlFEParts/FFlShells.H
+++ b/FFlLib/FFlFEParts/FFlShells.H
@@ -25,7 +25,6 @@ class FFlShellElementBase : public FFlElementBase
 {
 public:
   FFlShellElementBase(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlShellElementBase() {}
 
   double getThickness() const;
 
@@ -39,7 +38,6 @@ class FFlTRI3 : public FFlShellElementBase
 {
 public:
   FFlTRI3(int ID) : FFlShellElementBase(ID) {}
-  virtual ~FFlTRI3() {}
 
   static void init();
 
@@ -52,8 +50,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
   virtual FaVec3 mapping(double xi, double eta, double) const;
   virtual bool invertMapping(const FaVec3& X, double* Xi) const;
@@ -69,7 +66,6 @@ class FFlTRI6 : public FFlShellElementBase
 {
 public:
   FFlTRI6(int ID) : FFlShellElementBase(ID) {}
-  virtual ~FFlTRI6() {}
 
   static void init();
 
@@ -85,8 +81,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
   virtual int getNodalCoor(double* X, double* Y, double* Z) const;
 
 #ifdef FT_USE_MEMPOOL
@@ -99,7 +94,6 @@ class FFlQUAD4 : public FFlShellElementBase
 {
 public:
   FFlQUAD4(int ID) : FFlShellElementBase(ID) {}
-  virtual ~FFlQUAD4() {}
 
   static void init();
 
@@ -112,8 +106,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
   virtual FaVec3 mapping(double xi, double eta, double = 0.0) const;
   virtual bool invertMapping(const FaVec3& X, double* Xi) const;
@@ -129,7 +122,6 @@ class FFlQUAD8 : public FFlShellElementBase
 {
 public:
   FFlQUAD8(int ID) : FFlShellElementBase(ID) {}
-  virtual ~FFlQUAD8() {}
 
   static void init();
 
@@ -145,8 +137,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
 #ifdef FT_USE_MEMPOOL
   FFA_MAKE_MEMPOOL;
diff --git a/FFlLib/FFlFEParts/FFlSolids.C b/FFlLib/FFlFEParts/FFlSolids.C
index 120b3da..742d7c7 100644
--- a/FFlLib/FFlFEParts/FFlSolids.C
+++ b/FFlLib/FFlFEParts/FFlSolids.C
@@ -59,19 +59,20 @@ bool FFlTET4::getFaceNormals(std::vector<FaVec3>& normals, short int face,
 }
 
 
-bool FFlTET4::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                  FFaTensor3& inertia) const
+double FFlTET4::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
   FaVec3 v3(this->getNode(3)->getPos());
   FaVec3 v4(this->getNode(4)->getPos());
 
+  double volume = 0.0;
   FFaVolume::tetVolume(v1,v2,v3,v4,volume);
   FFaVolume::tetCenter(v1,v2,v3,v4,cog);
-  FFaVolume::tetMoment(v1,v2,v3,v4,inertia);
+  if (inertia)
+    FFaVolume::tetMoment(v1,v2,v3,v4,*inertia);
 
-  return true;
+  return volume;
 }
 
 
@@ -147,8 +148,7 @@ bool FFlWEDG6::getFaceNormals(std::vector<FaVec3>& normals, short int face,
 }
 
 
-bool FFlWEDG6::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlWEDG6::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
@@ -157,11 +157,13 @@ bool FFlWEDG6::getVolumeAndInertia(double& volume, FaVec3& cog,
   FaVec3 v5(this->getNode(5)->getPos());
   FaVec3 v6(this->getNode(6)->getPos());
 
+  double volume = 0.0;
   FFaVolume::wedVolume(v1,v2,v3,v4,v5,v6,volume);
   FFaVolume::wedCenter(v1,v2,v3,v4,v5,v6,cog);
-  FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,inertia);
+  if (inertia)
+    FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,*inertia);
 
-  return true;
+  return volume;
 }
 
 
@@ -211,8 +213,7 @@ bool FFlHEX8::getFaceNormals(std::vector<FaVec3>& normals, short int face,
 }
 
 
-bool FFlHEX8::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                  FFaTensor3& inertia) const
+double FFlHEX8::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(2)->getPos());
@@ -223,11 +224,13 @@ bool FFlHEX8::getVolumeAndInertia(double& volume, FaVec3& cog,
   FaVec3 v7(this->getNode(7)->getPos());
   FaVec3 v8(this->getNode(8)->getPos());
 
+  double volume = 0.0;
   FFaVolume::hexVolume(v1,v2,v3,v4,v5,v6,v7,v8,volume);
   FFaVolume::hexCenter(v1,v2,v3,v4,v5,v6,v7,v8,cog);
-  FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,inertia);
+  if (inertia)
+    FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,*inertia);
 
-  return true;
+  return volume;
 }
 
 
@@ -266,8 +269,7 @@ bool FFlTET10::getFaceNormals(std::vector<FaVec3>& normals, short int face,
 }
 
 
-bool FFlTET10::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlTET10::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   //TODO: Account for possibly curved edges
   FaVec3 v1(this->getNode(1)->getPos());
@@ -275,11 +277,13 @@ bool FFlTET10::getVolumeAndInertia(double& volume, FaVec3& cog,
   FaVec3 v3(this->getNode(5)->getPos());
   FaVec3 v4(this->getNode(10)->getPos());
 
+  double volume = 0.0;
   FFaVolume::tetVolume(v1,v2,v3,v4,volume);
   FFaVolume::tetCenter(v1,v2,v3,v4,cog);
-  FFaVolume::tetMoment(v1,v2,v3,v4,inertia);
+  if (inertia)
+    FFaVolume::tetMoment(v1,v2,v3,v4,*inertia);
 
-  return true;
+  return volume;
 }
 
 
@@ -353,8 +357,7 @@ bool FFlWEDG15::getFaceNormals(std::vector<FaVec3>& normals, short int face,
 }
 
 
-bool FFlWEDG15::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                    FFaTensor3& inertia) const
+double FFlWEDG15::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   FaVec3 v1(this->getNode(1)->getPos());
   FaVec3 v2(this->getNode(3)->getPos());
@@ -363,11 +366,13 @@ bool FFlWEDG15::getVolumeAndInertia(double& volume, FaVec3& cog,
   FaVec3 v5(this->getNode(12)->getPos());
   FaVec3 v6(this->getNode(14)->getPos());
 
+  double volume = 0.0;
   FFaVolume::wedVolume(v1,v2,v3,v4,v5,v6,volume);
   FFaVolume::wedCenter(v1,v2,v3,v4,v5,v6,cog);
-  FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,inertia);
+  if (inertia)
+    FFaVolume::wedMoment(v1,v2,v3,v4,v5,v6,*inertia);
 
-  return true;
+  return volume;
 }
 
 
@@ -407,8 +412,7 @@ bool FFlHEX20::getFaceNormals(std::vector<FaVec3>& normals, short int face,
 }
 
 
-bool FFlHEX20::getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const
+double FFlHEX20::getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const
 {
   //TODO: Account for possibly curved egdges
   FaVec3 v1(this->getNode(1)->getPos());
@@ -420,11 +424,13 @@ bool FFlHEX20::getVolumeAndInertia(double& volume, FaVec3& cog,
   FaVec3 v7(this->getNode(17)->getPos());
   FaVec3 v8(this->getNode(19)->getPos());
 
+  double volume = 0.0;
   FFaVolume::hexVolume(v1,v2,v3,v4,v5,v6,v7,v8,volume);
   FFaVolume::hexCenter(v1,v2,v3,v4,v5,v6,v7,v8,cog);
-  FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,inertia);
+  if (inertia)
+    FFaVolume::hexMoment(v1,v2,v3,v4,v5,v6,v7,v8,*inertia);
 
-  return true;
+  return volume;
 }
 
 #ifdef FF_NAMESPACE
diff --git a/FFlLib/FFlFEParts/FFlSolids.H b/FFlLib/FFlFEParts/FFlSolids.H
index 3ad5f7f..18b0145 100644
--- a/FFlLib/FFlFEParts/FFlSolids.H
+++ b/FFlLib/FFlFEParts/FFlSolids.H
@@ -25,7 +25,6 @@ class FFlTET4 : public FFlElementBase
 {
 public:
   FFlTET4(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlTET4() {}
 
   static void init();
 
@@ -37,8 +36,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
   virtual int checkOrientation(bool fixIt = false);
 
@@ -52,7 +50,6 @@ class FFlWEDG6 : public FFlElementBase
 {
 public:
   FFlWEDG6(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlWEDG6() {}
 
   static void init();
 
@@ -64,8 +61,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
 #ifdef FT_USE_MEMPOOL
   FFA_MAKE_MEMPOOL;
@@ -77,7 +73,6 @@ class FFlHEX8 : public FFlElementBase
 {
 public:
   FFlHEX8(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlHEX8() {}
 
   static void init();
 
@@ -89,8 +84,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
 #ifdef FT_USE_MEMPOOL
   FFA_MAKE_MEMPOOL;
@@ -102,7 +96,6 @@ class FFlTET10 : public FFlElementBase
 {
 public:
   FFlTET10(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlTET10() {}
 
   static void init();
 
@@ -114,8 +107,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
   virtual int checkOrientation(bool fixIt = false);
 
@@ -129,7 +121,6 @@ class FFlWEDG15 : public FFlElementBase
 {
 public:
   FFlWEDG15(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlWEDG15() {}
 
   static void init();
 
@@ -141,8 +132,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
 #ifdef FT_USE_MEMPOOL
   FFA_MAKE_MEMPOOL;
@@ -154,7 +144,6 @@ class FFlHEX20 : public FFlElementBase
 {
 public:
   FFlHEX20(int ID) : FFlElementBase(ID) {}
-  virtual ~FFlHEX20() {}
 
   static void init();
 
@@ -166,8 +155,7 @@ public:
   virtual bool getFaceNormals(std::vector<FaVec3>& normals, short int face = 1,
                               bool switchNormal = false) const;
 
-  virtual bool getVolumeAndInertia(double& volume, FaVec3& cog,
-                                   FFaTensor3& inertia) const;
+  virtual double getVolumeAndCoG(FaVec3& cog, FFaTensor3* inertia) const;
 
 #ifdef FT_USE_MEMPOOL
   FFA_MAKE_MEMPOOL;

From fb6e1a3f64c60ea635434a95a4e788c1c7acda3e Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Wed, 18 Mar 2026 14:39:11 +0100
Subject: [PATCH 4/6] Added: Utility to extract a point cloud from a FE model

---
 CMakeLists.txt |  14 +++
 main_pcloud.C  | 254 +++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 268 insertions(+)
 create mode 100644 main_pcloud.C

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 53de266..5e0b4b0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -108,6 +108,20 @@ if(IFEM_USE_FFLLIB)
   target_sources(ShellSim PRIVATE IFEMNastranReader.C)
 endif()
 
+ifem_add_application(
+  NAME
+    PointCloud
+  CONDITIONS
+    IFEM_USE_FFLLIB
+  SOURCES
+    main_pcloud.C
+    IFEMNastranReader.C
+  HEADERS
+    IFEMNastranReader.h
+  LIBRARIES
+    FFlLib
+)
+
 # Regression tests
 enable_testing()
 
diff --git a/main_pcloud.C b/main_pcloud.C
new file mode 100644
index 0000000..58c2ee6
--- /dev/null
+++ b/main_pcloud.C
@@ -0,0 +1,254 @@
+// $Id$
+//==============================================================================
+//!
+//! \file main_pcloud.C
+//!
+//! \date Mar 17 2026
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Utility for extracting a point cloud from a FE model.
+//!
+//==============================================================================
+
+#include "IFEMNastranReader.h"
+#include "FFlLinkHandler.H"
+#include "FFlGroup.H"
+#include "FFlFEParts/FFlAllFEParts.H"
+
+#include <algorithm>
+#include <numeric>
+#include <fstream>
+#include <sstream>
+#include <cstring>
+#include <cstdlib>
+
+
+/*!
+  \brief Auxiliary class to handle initialization of heap-allocated singeltons.
+*/
+
+class FFlInit
+{
+public:
+  //! \brief The constructor initializes the finite element library.
+  FFlInit() { FFl::initAllElements(); }
+  //! \brief The destructor releases the finite element library.
+  ~FFlInit() { FFl::releaseAllElements(); }
+};
+
+
+/*!
+  \brief Reads a Nastran bulk data file into a FFlLinkHandler object.
+*/
+
+FFlLinkHandler* readNastran (const char* fileName)
+{
+  std::ifstream is(fileName);
+  std::stringstream sets;
+
+  int lCount = IFEMNastranReader::parseSets(is,sets,fileName);
+  if (lCount < 0) return nullptr; // invalid Nastran file
+
+  FFlLinkHandler* fem = new FFlLinkHandler();
+  if (IFEMNastranReader reader(*fem,lCount);
+      reader.readFE(is,sets) && fem->resolve())
+    std::cout <<"\nParsing Nastran bulk data file succeeded."<< std::endl;
+  else
+  {
+    delete fem;
+    std::cerr <<"\n *** Parsing/resolving FE data failed.\n"
+              <<"     The FE model is probably not consistent and has not been"
+              <<" resolved completely."<< std::endl;
+    return nullptr;
+  }
+
+  int nnod = fem->getNodeCount(FFlLinkHandler::FFL_FEM);
+  int nel  = fem->getElementCount(FFlTypeInfoSpec::SHELL_ELM);
+  int nBel = fem->getElementCount(FFlTypeInfoSpec::BEAM_ELM);
+  std::cout <<"\nTotal number of nodes:          "<< nnod
+            <<"\nNumber of shell elements:       "<< nel;
+  if (nBel > 0)
+    std::cout <<"\nNumber of beam elements:        "<< nBel;
+  std::cout <<"\nNumber of constraint elements:  "
+            << fem->getElementCount(FFlTypeInfoSpec::CONSTRAINT_ELM)
+            <<"\nNumber of other elements:       "
+            << fem->getElementCount(FFlTypeInfoSpec::OTHER_ELM)
+            << std::endl;
+  if (size_t allN = fem->getNodeCount(FFlLinkHandler::FFL_ALL); allN > nnod)
+    std::cout <<"\n  ** Warning: This model contains "<< allN-nnod
+              <<" node(s) without any element connections (ignored)."
+              <<"\n     Please check the FE data file.\n"<< std::endl;
+
+  return fem;
+}
+
+
+/*!
+  \brief Main program for the point cloud extraction utility.
+*/
+
+int main (int argc, char** argv)
+{
+  if (argc < 2)
+  {
+    std::cout <<"usage: "<< argv[0] <<" <bdf-file>";
+    for (char d = 'x'; d <= 'z'; d++)
+      std::cout <<" [-"<< d <<"min <"<< d <<"0>] [-"<< d <<"max <"<< d << "1>]";
+    std::cout <<" [-group <ID>] [-out <filename>] [-distancesort]"<< std::endl;
+    return 0;
+  }
+
+  FFlInit _initializer;
+  FFlLinkHandler* fem = readNastran(argv[1]);
+  if (!fem) return 1;
+
+  FaVec3 max, min;
+  fem->getExtents(max,min);
+  std::cout <<"Model extension (diameter):     "
+            << (max-min).length() << std::endl;
+
+  // Lambda function parsing the point cloud bound.
+  auto&& parseRange = [](const char* option, const char* value,
+                         double& x0, double& x1, char d)
+  {
+    static char xmin[6] = "-xmin";
+    static char xmax[6] = "-xmax";
+    xmin[1] = d;
+    xmax[1] = d;
+    if (!strcmp(option,xmin))
+    {
+      if (double xval = atof(value); xval > x0)
+        x0 = xval;
+      return true;
+    }
+    else if (!strcmp(option,xmax))
+    {
+      if (double xval = atof(value); xval < x1)
+        x1 = xval;
+      return true;
+    }
+    return false;
+  };
+
+  bool distansort = !strncmp(argv[argc-1],"-dist",5);
+  const char* out = nullptr;
+  FFlGroup* group = nullptr;
+  for (int i = 2; i+1 < argc; i += 2)
+    if (!strcmp(argv[i],"-out"))
+      out = argv[i+1];
+    else if (!strncmp(argv[i],"-dist",5))
+      distansort = true, --i;
+    else if (!strcmp(argv[i],"-group"))
+    {
+      group = fem->getGroup(atoi(argv[i+1]));
+      std::cout <<"Element groups:";
+      for (GroupCIter it = fem->groupsBegin(); it != fem->groupsEnd(); ++it)
+        std::cout << (group == it->second ? "\n*\t" : "\n\t")
+                  << it->first <<"\t"<< it->second->getName();
+      std::cout << std::endl;
+    }
+    else for (int d = 0; d < 3; d++)
+      if (parseRange(argv[i],argv[i+1],min[d],max[d],'x'+d))
+        break;
+
+  std::cout <<"\nExtracting a point cloud for the domain ["<< min
+            <<"] x ["<< max <<"]";
+  if (group)
+    std::cout <<"\nfrom the element group "
+              << group->getID() <<": "<< group->getName();
+  std::cout << std::endl;
+
+  std::vector<FaVec3> cloudPts;
+  FaVec3 X0;
+
+  // Lambda function adding an element center to the point cloud.
+  auto&& addElmCenter = [&cloudPts,&X0,&min,&max](const FFlElementBase* elm)
+  {
+    FaVec3 X;
+    if (elm->getVolumeAndCoG(X) <= 0.0)
+      return; // zero-olum element, ignore
+
+    for (int i = 0; i < 3; i++)
+      if (X[i] < min[i] || X[i] > max[i])
+        return; // point is outisde the bounding box
+
+    cloudPts.push_back(X);
+    X0 += X;
+  };
+
+  if (group)
+    for (const GroupElemRef& elm : *group)
+      addElmCenter(elm.getReference());
+  else // no group specified, use all shell elements instead
+    for (ElementsCIter e = fem->elementsBegin(); e != fem->elementsEnd(); ++e)
+      if ((*e)->getCathegory() == FFlTypeInfoSpec::SHELL_ELM)
+        addElmCenter(*e);
+
+  delete fem; // release the FE model
+  if (cloudPts.empty())
+    return 2; // no points, invalid bounding box perhaps
+
+  X0 /= cloudPts.size();
+  std::cout <<"\nFound "<< cloudPts.size() <<" points, centered at "<< X0
+            << std::endl;
+
+  // Compute the theta angle of a cylindrical coordinate system,
+  // with the global X-axis as the cylinder axis and where
+  // zero angle is for points along the positive global Z-axis.
+  // Therefore, the angles are shifted by 90 degrees.
+  std::vector<double> angles;
+  angles.reserve(cloudPts.size());
+  for (const FaVec3& X : cloudPts)
+    if (double theta = (X-X0).getAsCylCoords(VX).y(); theta < 0.5*M_PI)
+      angles.push_back(theta*180.0/M_PI + 270.0);
+    else
+      angles.push_back(theta*180.0/M_PI - 90.0);
+
+  // Index-sort the points w.r.t. the angles
+  std::vector<int> indices(cloudPts.size());
+  std::iota(indices.begin(),indices.end(),0);
+  std::sort(indices.begin(),indices.end(),[&angles](int a, int b)
+            {
+              return angles[a] < angles[b];
+            });
+
+  std::cout <<"End points: "<< cloudPts[indices.front()]
+            <<" (theta=" << angles[indices.front()]
+            <<") and "<< cloudPts[indices.back()]
+            <<" (theta=" << angles[indices.back()] <<")"<< std::endl;
+
+  if (distansort)
+  {
+    // Reorder the points using the smallest-distance-to-neighbor approach
+    std::vector<int> newOrder;
+    newOrder.reserve(indices.size());
+    newOrder.push_back(indices.front());
+    indices.erase(indices.begin());
+    while (!indices.empty())
+    {
+      size_t indxmin = 0;
+      double distmin = (max-min).length();
+      const FaVec3& Xlast = cloudPts[newOrder.back()];
+      for (size_t i = 0; i < indices.size(); i++)
+        if (double dist = (cloudPts[indices[i]]-Xlast).length(); dist < distmin)
+        {
+          distmin = dist;
+          indxmin = i;
+        }
+      newOrder.push_back(indices[indxmin]);
+      indices.erase(indices.begin()+indxmin);
+    }
+    indices.swap(newOrder);
+  }
+
+  if (out)
+  {
+    // Print out the point cloud in sorted order
+    std::ofstream fs(out);
+    for (int idx : indices)
+      fs << cloudPts[idx] <<"\n";
+  }
+
+  return 0;
+}

From 63674b5572cbbf9669e86905ed7e9f677943f77a Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Thu, 19 Mar 2026 09:43:50 +0100
Subject: [PATCH 5/6] Add optional heading argument to FFlLinkHandler::dump().
 Fix issue openfedem/fedem-solvers#49: Don't enforce 6 DOFs on the WAVGM
 reference nodes if it couples translations only. Reduce scope of some
 temporaries using new if-syntax.

---
 FFlLib/FFlLinkHandler.C | 214 +++++++++++++++++-----------------------
 FFlLib/FFlLinkHandler.H |  19 ++--
 2 files changed, 100 insertions(+), 133 deletions(-)

diff --git a/FFlLib/FFlLinkHandler.C b/FFlLib/FFlLinkHandler.C
index a550dbb..937d1b5 100644
--- a/FFlLib/FFlLinkHandler.C
+++ b/FFlLib/FFlLinkHandler.C
@@ -436,36 +436,30 @@ bool FFlLinkHandler::updateCalculationFlag(int groupID, bool status)
 
 bool FFlLinkHandler::updateCalculationFlag(FFlPartBase* part, bool status)
 {
-  FFlGroup* tmpGroup = dynamic_cast<FFlGroup*>(part);
-  if (tmpGroup) {
-    for (const GroupElemRef& elm : *tmpGroup)
+  if (FFlGroup* group = dynamic_cast<FFlGroup*>(part); group)
+  {
+    for (const GroupElemRef& elm : *group)
       elm->setUpForCalculations(status);
-    return true;
   }
-
-  FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part);
-  if (attr) {
+  else if (FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part); attr)
+  {
     for (FFlElementBase* elm : myElements)
       if (elm->hasAttribute(attr))
         elm->setUpForCalculations(status);
-    return true;
   }
-
-  FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part);
-  if (elm) {
+  else if (FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part); elm)
     elm->setUpForCalculations(status);
-    return true;
-  }
+  else
+    return false;
 
-  return false;
+  return true;
 }
 
 
 bool FFlLinkHandler::updateCalculationFlag(const std::string& type,
                                            int id, bool status)
 {
-  FFlAttributeBase* attrib = this->getAttribute(type,id);
-  if (attrib)
+  if (FFlAttributeBase* attrib = this->getAttribute(type,id); attrib)
     return this->updateCalculationFlag(attrib,status);
   else
     return false;
@@ -521,28 +515,23 @@ bool FFlLinkHandler::setVisDetail(const FFlVDetail* detail)
 
 bool FFlLinkHandler::setVisDetail(FFlPartBase* part, const FFlVDetail* detail)
 {
-  FFlGroup* tmpGroup = dynamic_cast<FFlGroup*>(part);
-  if (tmpGroup) {
-    for (const GroupElemRef& elm : *tmpGroup)
+  if (FFlGroup* group = dynamic_cast<FFlGroup*>(part); group)
+  {
+    for (const GroupElemRef& elm : *group)
       elm->setDetail(detail);
-    return true;
   }
-
-  FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part);
-  if (attr) {
+  else if (FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part); attr)
+  {
     for (FFlElementBase* elm : myElements)
       if (elm->hasAttribute(attr))
         elm->setDetail(detail);
-    return true;
   }
-
-  FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part);
-  if (elm) {
+  else if (FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part); elm)
     elm->setDetail(detail);
-    return true;
-  }
+  else
+    return false;
 
-  return false;
+  return true;
 }
 
 
@@ -551,25 +540,14 @@ bool FFlLinkHandler::setVisDetail(const std::vector<FFlPartBase*>& parts,
 {
   std::vector<FFlAttributeBase*> attribParts;
 
-  for (FFlPartBase* part : parts)
-  {
-    FFlGroup* tmpGroup = dynamic_cast<FFlGroup*>(part);
-    if (tmpGroup)
-      for (const GroupElemRef& elm : *tmpGroup)
+  for (FFlPartBase* pt : parts)
+    if (FFlGroup* group = dynamic_cast<FFlGroup*>(pt); group)
+      for (const GroupElemRef& elm : *group)
         elm->setDetail(detail);
-    else
-    {
-      FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part);
-      if (attr)
-        attribParts.push_back(attr);
-      else
-      {
-        FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part);
-        if (elm)
-          elm->setDetail(detail);
-      }
-    }
-  }
+    else if (FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(pt); attr)
+      attribParts.push_back(attr);
+    else if (FFlElementBase* elm = dynamic_cast<FFlElementBase*>(pt); elm)
+      elm->setDetail(detail);
 
   if (!attribParts.empty())
     for (FFlElementBase* elm : myElements)
@@ -584,28 +562,23 @@ bool FFlLinkHandler::setVisDetail(const std::vector<FFlPartBase*>& parts,
 bool FFlLinkHandler::setVisAppearance(FFlPartBase* part,
                                       const FFlVAppearance* app)
 {
-  FFlGroup* tmpGroup = dynamic_cast<FFlGroup*>(part);
-  if (tmpGroup) {
-    for (const GroupElemRef& elm : *tmpGroup)
+  if (FFlGroup* group = dynamic_cast<FFlGroup*>(part); group)
+  {
+    for (const GroupElemRef& elm : *group)
       elm->setAppearance(app);
-    return true;
   }
-
-  FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part);
-  if (attr) {
+  else if (FFlAttributeBase* attr = dynamic_cast<FFlAttributeBase*>(part); attr)
+  {
     for (FFlElementBase* elm : myElements)
       if (elm->hasAttribute(attr))
         elm->setAppearance(app);
-    return true;
   }
-
-  FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part);
-  if (elm) {
+  else if (FFlElementBase* elm = dynamic_cast<FFlElementBase*>(part); elm)
     elm->setAppearance(app);
-    return true;
-  }
+  else
+    return false;
 
-  return false;
+  return true;
 }
 #endif
 
@@ -701,9 +674,9 @@ void FFlLinkHandler::addVisual(FFlVisualBase* visual, bool sortOnInsert)
 
 void FFlLinkHandler::setRunningIdxOnAppearances()
 {
-  FFlVAppearance* vapp = NULL; int idx = 0;
+  int idx = 0;
   for (FFlVisualBase* vis : myVisuals)
-    if ((vapp = dynamic_cast<FFlVAppearance*>(vis)))
+    if (FFlVAppearance* vapp = dynamic_cast<FFlVAppearance*>(vis); vapp)
       vapp->runningIdx = idx++;
 }
 #endif
@@ -769,14 +742,13 @@ FFlVAppearance* FFlLinkHandler::getAppearance(int ID) const
                                                         myVisuals.end(), ID,
                                                         FFlFEPartBaseLess());
 
-  FFlVAppearance* vapp = NULL;
   while (ep.first != ep.second)
-    if ((vapp = dynamic_cast<FFlVAppearance*>(*ep.first)))
-      break;
+    if (FFlVAppearance* vapp = dynamic_cast<FFlVAppearance*>(*ep.first); vapp)
+      return vapp;
     else
       ep.first++;
 
-  return vapp;
+  return NULL;
 }
 
 
@@ -788,14 +760,13 @@ FFlVDetail* FFlLinkHandler::getDetail(int ID) const
                                                         myVisuals.end(), ID,
                                                         FFlFEPartBaseLess());
 
-  FFlVDetail* vdet = NULL;
   while (ep.first != ep.second)
-    if ((vdet = dynamic_cast<FFlVDetail*>(*ep.first)))
-      break;
+    if (FFlVDetail* vdet = dynamic_cast<FFlVDetail*>(*ep.first); vdet)
+      return vdet;
     else
       ep.first++;
 
-  return vdet;
+  return NULL;
 }
 #endif
 
@@ -1217,7 +1188,7 @@ bool FFlLinkHandler::addGroup(FFlGroup* group, bool silence)
 {
   if (!group) return false;
 
-  if (myGroupMap.insert({group->getID(),group}).second)
+  if (myGroupMap.emplace(group->getID(),group).second)
   {
     isResolved = false;
     return true;
@@ -1237,7 +1208,7 @@ bool FFlLinkHandler::addAttribute(FFlAttributeBase* attr, bool silence,
 {
   if (!attr) return false;
 
-  if (myAttributes[name].insert({attr->getID(),attr}).second)
+  if (myAttributes[name].emplace(attr->getID(),attr).second)
   {
     isResolved = false;
     return true;
@@ -1334,13 +1305,12 @@ bool FFlLinkHandler::removeAttribute(const std::string& typeName,
 
 FFlVDetail* FFlLinkHandler::getPredefDetail(int detailType)
 {
-  FFlVDetail* vdet = NULL;
   for (FFlVisualBase* vis : myVisuals)
-    if ((vdet = dynamic_cast<FFlVDetail*>(vis)))
+    if (FFlVDetail* vdet = dynamic_cast<FFlVDetail*>(vis); vdet)
       if (vdet->detail.getValue() == detailType)
 	return vdet;
 
-  vdet = new FFlVDetail(this->getNewVisualID());
+  FFlVDetail* vdet = new FFlVDetail(this->getNewVisualID());
   vdet->detail.setValue(detailType);
   this->addVisual(vdet,true);
 
@@ -1371,10 +1341,9 @@ void FFlLinkHandler::getAllInternalCoordSys(std::vector<FaMat34>& mxes) const
   AttributeTypeCIter mapit = myAttributes.find("PCOORDSYS");
   if (mapit == myAttributes.end()) return;
 
-  FFlPCOORDSYS* lcs;
   mxes.reserve(mapit->second.size());
   for (const AttributeMap::value_type& attr : mapit->second)
-    if ((lcs = dynamic_cast<FFlPCOORDSYS*>(attr.second)))
+    if (FFlPCOORDSYS* lcs = dynamic_cast<FFlPCOORDSYS*>(attr.second); lcs)
     {
       mxes.push_back(FaMat34());
       mxes.back().makeCS_Z_XZ(lcs->Origo.getValue(),
@@ -1410,14 +1379,13 @@ FFlNode* FFlLinkHandler::findFreeNodeAtPoint(const FaVec3& point,
   FFlNode* closest = NULL;
   double closestdist = DBL_MAX;
   double sqrTol = tol > sqrt(DBL_MAX) ? DBL_MAX : tol*tol;
-  double sqrdist, xd, yd, zd;
 
   for (FFlNode* node : myNodes)
     if (node->hasDOFs(dofFilter) || node->isExternal() || node->isRefNode())
-      if ((xd = fabs(node->getPos().x() - point.x())) < tol)
-        if ((yd = fabs(node->getPos().y() - point.y())) < tol)
-          if ((zd = fabs(node->getPos().z() - point.z())) < tol)
-            if ((sqrdist = xd*xd+yd*yd+zd*zd) < sqrTol)
+      if (double xd = fabs(node->getPos().x() - point.x()); xd < tol)
+        if (double yd = fabs(node->getPos().y() - point.y()); yd < tol)
+          if (double zd = fabs(node->getPos().z() - point.z()); zd < tol)
+            if (double sqrdist = xd*xd+yd*yd+zd*zd; sqrdist < sqrTol)
             {
               if (closestdist >= sqrTol)
               {
@@ -1617,13 +1585,12 @@ FFlNode* FFlLinkHandler::findClosestNode(const FaVec3& point) const
 {
   FFlNode* closest = NULL;
   double closestdist = DBL_MAX;
-  double sqrdist;
 
   for (FFlNode* node : myNodes)
-    if ((sqrdist = (point - node->getPos()).sqrLength()) <= closestdist)
+    if (double d2 = (point - node->getPos()).sqrLength(); d2 <= closestdist)
     {
       closest     = node;
-      closestdist = sqrdist;
+      closestdist = d2;
     }
 
   return closest;
@@ -1641,19 +1608,18 @@ FFlElementBase* FFlLinkHandler::findClosestElement(const FaVec3& point,
 {
   FFlElementBase* closest = NULL;
   double closestdist = DBL_MAX;
-  double sqrdist;
 
   // OTHER_ELM = the total number of element cathegories
   std::vector<bool> typeOk(FFlTypeInfoSpec::OTHER_ELM+1,wantedTypes.empty());
   for (FFlTypeInfoSpec::Cathegory type : wantedTypes)
     typeOk[type] = true;
 
-  for (FFlElementBase* elm : myElements)
-    if (typeOk[elm->getCathegory()])
-      if ((sqrdist = (point - elm->getNodeCenter()).sqrLength()) <= closestdist)
+  for (FFlElementBase* e : myElements)
+    if (typeOk[e->getCathegory()])
+      if (double d2 = (point - e->getNodeCenter()).sqrLength(); d2 <= closestdist)
       {
-        closest     = elm;
-        closestdist = sqrdist;
+        closest     = e;
+        closestdist = d2;
       }
 
   return closest;
@@ -1670,13 +1636,12 @@ FFlElementBase* FFlLinkHandler::findClosestElement(const FaVec3& point,
 {
   FFlElementBase* closest = NULL;
   double closestdist = DBL_MAX;
-  double sqrdist;
 
-  for (const GroupElemRef& elm : group)
-    if ((sqrdist = (point - elm->getNodeCenter()).sqrLength()) <= closestdist)
+  for (const GroupElemRef& e : group)
+    if (double d2 = (point - e->getNodeCenter()).sqrLength(); d2 <= closestdist)
     {
-      closest     = elm.getReference();
-      closestdist = sqrdist;
+      closest     = e.getReference();
+      closestdist = d2;
     }
 
   return closest;
@@ -1820,15 +1785,10 @@ void FFlLinkHandler::findWindowedNodes(std::map<int,FaVec3>& nodes,
   int nVertices = myVertices.size();
   for (int idx : indices)
     if (idx < nVertices && visited.insert(idx).second)
-    {
-      const FFlVertex* vrtx = static_cast<const FFlVertex*>(myVertices[idx]);
-      if (vrtx && vrtx->getNode())
-      {
-        FaVec3 globalPos = lCS * (*vrtx);
-        if (isInsideWindow(globalPos))
-          nodes[vrtx->getNode()->getID()] = globalPos;
-      }
-    }
+      if (const FFlVertex* vx = static_cast<const FFlVertex*>(myVertices[idx]);
+          vx && vx->getNode())
+        if (FaVec3 globalPos = lCS * (*vx); isInsideWindow(globalPos))
+          nodes[vx->getNode()->getID()] = globalPos;
 }
 
 #endif
@@ -1945,10 +1905,9 @@ bool FFlLinkHandler::resolve(bool subdivParabolic, bool fromSESAM)
                  << attr.second->getID() <<" failed\n";
 
   // Remove loose nodes from WAVGM and RGD elements
-  int nNod = 0;
   ElementsVec toBeErased;
   for (FFlElementBase* elm : myElements)
-    if ((nNod = elm->getNodeCount()) < 1)
+    if (int nNod = elm->getNodeCount(); nNod < 1)
       toBeErased.push_back(elm);
     else if (elm->getTypeName() == "WAVGM")
     {
@@ -1983,7 +1942,15 @@ bool FFlLinkHandler::resolve(bool subdivParabolic, bool fromSESAM)
         continue;
       }
 
-      (*refn)->pushDOFs(6); // Ensure the reference node is not flagged DOF-less
+      FFlAttributeBase* pAtt = elm->getAttribute("PWAVGM");
+
+      // Ensure the reference node is not flagged DOF-less.
+      // Issue fedem-solvers#49: Don't enforce rotational DOFs on the reference
+      // node if the constraint element couples translational DOFs only.
+      if (pAtt && abs(static_cast<FFlPWAVGM*>(pAtt)->refC.getValue()) < 322)
+        (*refn)->pushDOFs(3);
+      else
+        (*refn)->pushDOFs(6);
 
       if (looseNodes.empty()) continue;
 
@@ -1992,18 +1959,16 @@ bool FFlLinkHandler::resolve(bool subdivParabolic, bool fromSESAM)
              <<" loose nodes from WAVGM element "<< elm->getID();
       if (!static_cast<FFlWAVGM*>(elm)->removeMasterNodes(looseNodes))
         ++nError;
-      else
+      else if (pAtt)
       {
-        FFlAttributeBase* newAtt = elm->getAttribute("PWAVGM");
-        if (newAtt)
-        {
-          int ID = 1;
-          while (this->getAttribute("PWAVGM",ID)) ID++;
-          newAtt = static_cast<FFlPWAVGM*>(newAtt)->removeWeights(looseNodes,nNod);
-          newAtt->setID(ID);
-          elm->clearAttribute("PWAVGM"); // Replace the old attribute with the new one
-          elm->setAttribute(this->getAttribute("PWAVGM",this->addUniqueAttribute(newAtt)));
-        }
+        int ID = 1;
+        while (this->getAttribute("PWAVGM",ID)) ID++;
+        pAtt = static_cast<FFlPWAVGM*>(pAtt)->removeWeights(looseNodes,nNod);
+        pAtt->setID(ID);
+        // Replace the old attribute with the new one
+        elm->clearAttribute("PWAVGM");
+        ID = this->addUniqueAttribute(pAtt);
+        elm->setAttribute(this->getAttribute("PWAVGM",ID));
       }
     }
     else if (elm->getTypeName() == "RGD")
@@ -2250,13 +2215,13 @@ void FFlLinkHandler::countElements() const
 }
 
 
-void FFlLinkHandler::dump() const
+void FFlLinkHandler::dump(const std::string& heading) const
 {
   if (myNumElements.empty())
     this->countElements();
 
-  std::cout <<"\nFFlLinkHandler::dump()"
-            <<"\n   Elements:        "<< myElements.size();
+  std::cout <<"\nFFlLinkHandler::dump("<< heading
+            <<")\n   Elements:        "<< myElements.size();
   for (const ElmTypeCount::value_type& ec : myNumElements)
     std::cout <<"\n\t"<< ec.first <<"\t"<< ec.second;
 
@@ -2419,7 +2384,6 @@ int FFlLinkHandler::createConnector(const FFaCompoundGeometry& compound,
   cItems.clear();
 
   FFlNode* refNode = new FFlNode(this->getNewNodeID(),nodePos);
-  refNode->pushDOFs(6);
   this->addNode(refNode,true);
   cItems.addNode(refNode->getID());
   ListUI <<"  -> Creating FE node "<< refNode->getID() <<"\n";
diff --git a/FFlLib/FFlLinkHandler.H b/FFlLib/FFlLinkHandler.H
index f4cfb17..0674ff1 100644
--- a/FFlLib/FFlLinkHandler.H
+++ b/FFlLib/FFlLinkHandler.H
@@ -96,10 +96,11 @@ public:
   bool addGroup(FFlGroup* group, bool silence = false);
   void addLoad(FFlLoadBase* load, bool sortOnInsert = false);
   bool addAttribute(FFlAttributeBase* attr, bool silence = false);
-  bool addAttribute(FFlAttributeBase* attr, bool silence, const std::string& name);
+  bool addAttribute(FFlAttributeBase* attr, bool silence,
+                    const std::string& name);
   int  addUniqueAttribute(FFlAttributeBase* attr, bool silence = false);
   int  addUniqueAttributeCS(FFlAttributeBase*& attr);
-  bool removeAttribute(const std::string& typeName, int ID, bool silence = false);
+  bool removeAttribute(const std::string& name, int ID, bool silence = false);
 #ifdef FT_USE_VISUALS
   void addVisual(FFlVisualBase* visual, bool sortOnInsert = false);
   void setRunningIdxOnAppearances();
@@ -138,10 +139,12 @@ public:
   bool getLoads(int ID, LoadsVec& loads) const;
   void getLoadCases(std::set<int>& IDs) const;
 
-  // External-to-internal node and (finite) element number mapping
+  //! \brief External-to-internal node number mapping.
   int getIntNodeID(int extID) const;
+  //! \brief External-to-internal (finite) element number mapping.
   int getIntElementID(int extID) const;
 
+  //! \brief Removes the given elements from the myElements container.
   void removeElements(const ElementsVec& elms);
 
   int getNewElmID() const;
@@ -197,12 +200,10 @@ public:
   const FFlrVxToElmMap& getVxToElementMapping();
 
   using WindowTester = bool(*)(const FaVec3&);
-
   void findWindowedNodes(std::map<int,FaVec3>& nodes,
                          const std::vector<int>& indices,
                          const FaMat34& lCS, bool lFirst,
                          WindowTester isInsideWindow) const;
-
 #endif
 
   void getAllInternalCoordSys(std::vector<FaMat34>& mxes) const;
@@ -274,13 +275,15 @@ public:
   void initiateCalculationFlag(bool status = false);
   bool updateCalculationFlag(int groupId, bool status = true);
   bool updateCalculationFlag(FFlPartBase* part, bool status = true);
-  bool updateCalculationFlag(const std::string& attType, int id, bool status = true);
+  bool updateCalculationFlag(const std::string& attType, int id,
+                             bool status = true);
 
 #ifdef FT_USE_VISUALS
   // Visual settings
   bool setVisDetail(const FFlVDetail* det); // on all elements
   bool setVisDetail(FFlPartBase* part, const FFlVDetail* det);
-  bool setVisDetail(const std::vector<FFlPartBase*>& parts, const FFlVDetail* det);
+  bool setVisDetail(const std::vector<FFlPartBase*>& parts,
+                    const FFlVDetail* det);
   bool setVisAppearance(FFlPartBase* partSpec, const FFlVAppearance* app);
   void updateGroupVisibilityStatus();
 
@@ -297,7 +300,7 @@ public:
   const Strings& getOP2files() const { return myOP2files; }
   size_t getNumberOfGenDofs() const { return nGenDofs; }
 
-  void dump() const;
+  void dump(const std::string& heading = "") const;
 
 #ifdef FT_USE_CONNECTORS
   int createConnector(const FFaCompoundGeometry& compound,

From afcf74232a9c150ef2dbe643f585bdb1c92d3e9b Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Mon, 23 Mar 2026 10:38:51 +0100
Subject: [PATCH 6/6] Avoid doxy warning when referring to member of base class

---
 ASMu2DNastran.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ASMu2DNastran.h b/ASMu2DNastran.h
index ce4aaab..eb13381 100644
--- a/ASMu2DNastran.h
+++ b/ASMu2DNastran.h
@@ -220,7 +220,7 @@ class ASMuBeam : public ASMu1DLag
   //! \brief Empty destructor.
   virtual ~ASMuBeam() {}
 
-  //! \brief Sets \ref nGauss to 1 due to explicit matrices.
+  //! \brief Sets \a nGauss to 1 due to explicit matrices.
   virtual void setGauss(int) { nGauss = 1; }
 
   //! \brief Retrieves the properties for element with index \a id.