Separate utility for extracting an ordered point cloud from a Nastran file

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

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.

CMakeLists.txt

@@ -97,9 +97,30 @@ 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()
+
+ifem_add_application(
+  NAME
+    PointCloud
+  CONDITIONS
+    IFEM_USE_FFLLIB
+  SOURCES
+    main_pcloud.C
+    IFEMNastranReader.C
+  HEADERS
+    IFEMNastranReader.h
+  LIBRARIES
+    FFlLib
+)
 
 # Regression tests
 enable_testing()

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>
 
 

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

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

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];
 }

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