Support for superelement analysis with recovery

ANDES/wavgmConstrEqn.f90

@@ -10,12 +10,12 @@
 !> @brief Weighted Average Motion constraint handling.
 
 !!==============================================================================
-!> @brief Computes constraint equation coefficients for a WAVGM element.
+!> @brief Computes the constraint equation coefficients for a WAVGM element.
 !>
 !> @param[in] iel Element index
 !> @param[in] lDof Local index of dependent DOF to compute coefficients for
 !> @param[in] nM Number of independent nodes in current element
-!> @param[in] nW SiNumber of independent nodes in current element
+!> @param[in] nW Size of the @a weight array
 !> @param[in] indC Nodal component indices (common for all nodes)
 !> @param[in] tenc Table of nodal coordinates for current element
 !> @param[in] weight Independent DOF weights for current element
@@ -26,10 +26,13 @@
 !> @param[in] ipsw Print switch
 !> @param[in] lpu File unit number for res-file output
 !>
-!> @details This subroutine pre-processes a Weighted Average Motion (WAVGM)
-!> elements (also known as RBE3 in Nastran). The coefficients that couples
-!> a dependent DOF at the reference node of a Weighted AVerage Motion element
-!> to the independent nodal DOFs of the same element are computed.
+!> @details This subroutine pre-processes a Weighted AVeraGe Motion (WAVGM)
+!> element (also known as RBE3 in Nastran). The coefficients that couples
+!> a dependent DOF at the reference node of a WAVGM element to the
+!> independent nodal DOFs of the same element are computed.
+!>
+!> @note This is a slightly modified version of the subroutine with same name
+!> from the openfedem project. See https://github.com/openfedem/fedem-solvers.
 !>
 !> @callergraph
 !>
@@ -88,6 +91,7 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
         call DAXPY (nM,1.0_dp/sumWM,weight(iFrst),1,omega(lDof),nndof)
      else if (ipsw > 0) then
         write(lpu,600) iel,lDof,lDof,sumWM,epsDiv0_p
+        write(lpu,620) weight(iFrst:iLast)
      end if
 
   else if (indC(lDof) < 0) then
@@ -125,7 +129,8 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
            call DAXPY (nM,-1.0_dp/sumWM,work(1),1,omega(j),nndof)
         else if (ipsw > 0) then
            write(lpu,600) iel,lDof,j,sumWM,tolWM
-           write(lpu,620) work(1:nM) / sumWM
+           write(lpu,620) dX(j,2:1+nM)
+           write(lpu,620) work(1:nM)
         end if
      end if
      if (indC(k) /= 0) then
@@ -152,7 +157,8 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
            call DAXPY (nM,1.0_dp/sumWM,work(1),1,omega(k),nndof)
         else if (ipsw > 0) then
            write(lpu,600) iel,lDof,k,sumWM,tolWM
-           write(lpu,620) work(1:nM) / sumWM
+           write(lpu,620) dX(j,2:1+nM)
+           write(lpu,620) work(1:nM)
         end if
      end if
 
@@ -187,7 +193,8 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
               call DAXPY (nM,dX(j,1)/sumWM,work(1),1,omega(i),nndof)
            else if (ipsw > 0) then
               write(lpu,600) iel,lDof,i,sumWM,tolWM
-              write(lpu,620) work(1:nM) * dX(j,1)/sumWM
+              write(lpu,620) dX(i,2:1+nM)
+              write(lpu,620) work(1:nM) * dX(j,1)
            end if
         else if (ipsw > 0) then
            write(lpu,610) iel,lDof,i,dX(j,1),tolX(j)
@@ -213,11 +220,13 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
               sumWM = sumSqr(dX(i,2:1+nM),dX(j,2:1+nM)) * real(nM,dp)
               call DCOPY (nM,dX(i,2),size(dX,1),work(1),1)
            end if
+           tolWM = tolX(i)*tolX(i) + tolX(j)*tolX(j)
            if (sumWM > tolWM) then
               call DAXPY (nM,-dX(j,1)/sumWM,work(1),1,omega(j),nndof)
            else if (ipsw > 0) then
               write(lpu,600) iel,lDof,j,sumWM,tolWM
-              write(lpu,620) work(1:nM) * dX(j,1)/sumWM
+              write(lpu,620) dX(i,2:1+nM)
+              write(lpu,620) work(1:nM) * dX(j,1)
            end if
         else if (ipsw > 0) then
            write(lpu,610) iel,lDof,j,dX(j,1),tolX(j)
@@ -237,6 +246,7 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
               call DAXPY (nM,-dX(j,1)/sumWM,weight(iFrst),1,omega(3+k),nndof)
            else if (ipsw > 0) then
               write(lpu,600) iel,lDof,3+k,sumWM,epsDiv0_p
+              write(lpu,620) weight(iFrst:iLast)
            end if
         else if (ipsw > 0) then
            write(lpu,610) iel,lDof,3+k,dX(j,1),tolX(j)
@@ -272,7 +282,8 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
               call DAXPY (nM,dX(k,1)/sumWM,work(1),1,omega(i),nndof)
            else if (ipsw > 0) then
               write(lpu,600) iel,lDof,i,sumWM,tolWM
-              write(lpu,620) work(1:nM) * dX(k,1)/sumWM
+              write(lpu,620) dX(k,2:1+nM)
+              write(lpu,620) work(1:nM) * dX(k,1)
            end if
         else if (ipsw > 0) then
            write(lpu,610) iel,lDof,i,dX(k,1),tolX(k)
@@ -303,7 +314,8 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
               call DAXPY (nM,-dX(k,1)/sumWM,work(1),1,omega(k),nndof)
            else if (ipsw > 0) then
               write(lpu,600) iel,lDof,k,sumWM,tolWM
-              write(lpu,620) work(1:nM) * dX(k,1)/sumWM
+              write(lpu,620) dX(k,2:1+nM)
+              write(lpu,620) work(1:nM) * dX(k,1)
            end if
         else if (ipsw > 0) then
            write(lpu,610) iel,lDof,k,dX(k,1),tolX(k)
@@ -323,6 +335,7 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
               call DAXPY (nM,dX(k,1)/sumWM,weight(iFrst),1,omega(3+j),nndof)
            else if (ipsw > 0) then
               write(lpu,600) iel,lDof,3+j,sumWM,epsDiv0_p
+              write(lpu,620) weight(iFrst:iLast)
            end if
         else if (ipsw > 0) then
            write(lpu,610) iel,lDof,3+j,dX(k,1),tolX(k)
@@ -349,14 +362,14 @@ subroutine wavgmConstrEqn (iel,lDof,nM,nW,indC,tenc,weight, &
 610 format('  ** Warning: WAVGM element',I10, &
          & ': Ignored coupling by dependent DOF',I2,' to independent DOFs',I2 &
          / 14X, 'due to small eccentricity',1PE13.5,'  tolerance =',E12.5 )
-620 format(14X,1P10E13.5/(14X,1P10E13.5))
+620 format(12X,1P10E13.5/(12X,1P10E13.5))
 690 format(' *** Error: Indices out of range:',3I6,/12X,'For WAVGM element',I10)
 
 contains
 
   !> @brief Computes the nodal coordinates relative to the centre of gravity.
-  !> @details If the weights (W) are present, a weighted CoG is used.
-  !> This subroutine also recomputes the geometric tolerances (tolX).
+  !> @details If the weights @a W are present, a weighted CoG is used.
+  !> This subroutine also recomputes the geometric tolerances, @a tolX.
   subroutine computePointCoords (W)
     real(dp), intent(in), optional :: W(:)
     integer  :: n

ASMu2DNastran.C

@@ -160,7 +160,11 @@ bool ASMu2DNastran::read (std::istream& is)
         sets << cline << '\n';
     }
 
-  if (!is) return false; // No bulk data file
+  if (!is)
+  {
+    std::cerr <<"\n *** Invalid Nastran bulk data file."<< std::endl;
+    return false;
+  }
 
 #ifdef HAS_FFLLIB
 
@@ -170,9 +174,11 @@ bool ASMu2DNastran::read (std::istream& is)
   {
     std::cerr <<"\n *** Parsing/resolving FE data failed.\n"
               <<"     The FE model is probably not consistent and has not been"
-              <<" resolved completely.\n";
+              <<" resolved completely."<< std::endl;
     return false;
   }
+  else
+    while (is.ignore()); // To avoid trying to read another patch
 
   for (int eId : fixRBE3)
   {
@@ -298,7 +304,7 @@ bool ASMu2DNastran::read (std::istream& is)
     else if ((*e)->getTypeName() == "RGD" && mnpc.size() > 1)
     {
 #if INT_DEBUG > 5
-      std::cout <<"Rigid element "<< eid <<": master = "<< MLGN[mnpc.front()]
+      std::cout <<"\nRigid element "<< eid <<": master = "<< MLGN[mnpc.front()]
                 <<" slaves =";
       for (size_t i = 1; i < mnpc.size(); i++) std::cout <<" "<< MLGN[mnpc[i]];
       std::cout << std::endl;
@@ -311,10 +317,12 @@ bool ASMu2DNastran::read (std::istream& is)
     else if ((*e)->getTypeName() == "WAVGM" && mnpc.size() > 1)
     {
 #if INT_DEBUG > 5
-      std::cout <<"Weighted average motion element "<< eid
+      std::cout <<"\nWeighted average motion element "<< eid
                 <<": reference node = "<< MLGN[mnpc.front()]
-                <<"\n\tmasters =";
-      for (size_t i = 1; i < mnpc.size(); i++) std::cout <<" "<< MLGN[mnpc[i]];
+                <<"\n     masters =";
+      for (size_t i = 1; i < mnpc.size(); i++)
+        std::cout << (i > 1 && i%10 == 1 ? "\n               " : " ")
+                  << MLGN[mnpc[i]];
       std::cout << std::endl;
 #endif
       spiders.push_back(mnpc);
@@ -804,12 +812,14 @@ void ASMu2DNastran::addFlexibleCoupling (int eId, int lDof, const int* indC,
       for (int mDof = 1; mDof <= 6; mDof++, omega++)
         if (fabs(*omega) > Zero)
           cons->addMaster(MLGN[mnpc[iM]],mDof,*omega);
-#if INT_DEBUG > 1
-        else
+        else if (ips > 0)
           std::cout <<"  ** Ignoring small coupling coefficient "<< *omega
                     <<" to local dof "<< mDof
                     <<" of master node "<< MLGN[mnpc[iM]]
                     <<" in RBE3 element "<< eId << std::endl;
+
+#if INT_DEBUG > 1
+  std::cout <<"Added constraint: "<< *cons << std::endl;
 #endif
 #endif
 }

CMakeLists.txt

@@ -77,12 +77,21 @@ include(IFEMTesting)
 if(IFEM_USE_ANDES AND TARGET FFlLib)
   file(GLOB SHELL_TESTFILES RELATIVE ${PROJECT_SOURCE_DIR}/Test
                                      ${PROJECT_SOURCE_DIR}/Test/*.reg)
+  file(GLOB SHELL_VTF_FILES RELATIVE ${PROJECT_SOURCE_DIR}/Test
+                                     ${PROJECT_SOURCE_DIR}/Test/*.vreg)
+  if(NOT IFEM_USE_FMX)
+    list(REMOVE_ITEM SHELL_TESTFILES Q4-sup.reg Q4andBeam-sup.reg)
+    list(REMOVE_ITEM SHELL_VTF_FILES Q4-sup.vreg)
+  endif()
 endif()
 
 if(NOT MPI_FOUND OR IFEM_SERIAL_TESTS_IN_PARALLEL)
   foreach(TESTFILE ${SHELL_TESTFILES})
     ifem_add_test(${TESTFILE} ShellSim)
   endforeach()
+  foreach(TESTFILE ${SHELL_VTF_FILES})
+    ifem_add_vtf_test(${TESTFILE} ShellSim)
+  endforeach()
 endif()
 
 list(APPEND TEST_APPS ShellSim)

HasPointLoads.C

@@ -0,0 +1,105 @@
+// $Id$
+//==============================================================================
+//!
+//! \file HasPointLoads.C
+//!
+//! \date Jun 12 2025
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Class representing a nodal point load container.
+//!
+//==============================================================================
+
+#include "HasPointLoads.h"
+#include "AlgEqSystem.h"
+#include "SystemMatrix.h"
+#include "SAM.h"
+#include "Functions.h"
+#include "Utilities.h"
+#include "IFEM.h"
+#include "tinyxml2.h"
+
+
+HasPointLoads::~HasPointLoads ()
+{
+  for (PointLoad& load : myLoads)
+    delete load.func;
+}
+
+
+bool HasPointLoads::parseLoad (const tinyxml2::XMLElement* elem)
+{
+  int inod = 0, ldof = 0;
+  utl::getAttribute(elem,"node",inod);
+  utl::getAttribute(elem,"dof",ldof);
+  if (inod <= 0 || ldof <= 0 || !elem->FirstChild())
+    return false;
+
+  IFEM::cout <<"\tNode "<< inod <<" dof "<< ldof <<" Load: ";
+  std::string type("constant");
+  utl::getAttribute(elem,"type",type);
+  if (type == "file")
+  {
+    IntSet dofs = utl::getDigits(ldof);
+    char* funcv = strdup(elem->FirstChild()->Value());
+    char* fname = strtok(funcv," ");
+    char* colmn = nullptr;
+    IFEM::cout <<"\""<< fname <<"\"";
+    for (int ldof : dofs)
+      if ((colmn = strtok(nullptr," ")))
+      {
+	IFEM::cout <<" ("<< ldof <<","<< colmn <<")";
+	myLoads.emplace_back(inod, ldof, new LinearFunc(fname,atoi(colmn)));
+      }
+    free(funcv);
+    IFEM::cout << std::endl;
+  }
+  else if (ldof <= 6)
+  {
+    ScalarFunc* f = nullptr;
+    if (type == "constant")
+    {
+      f = new ConstantFunc(atof(elem->FirstChild()->Value()));
+      IFEM::cout << (*f)(0.0) << std::endl;
+    }
+    else
+      f = utl::parseTimeFunc(elem->FirstChild()->Value(),type);
+
+    myLoads.emplace_back(inod,ldof,f);
+  }
+  else
+    IFEM::cout <<" (invalid, ignored)."<< std::endl;
+
+  return true;
+}
+
+
+bool HasPointLoads::renumberLoadedNodes (const std::map<int,int>& nodeMap)
+{
+  bool ok = true;
+  for (PointLoad& load : myLoads)
+    if (load.inod > 0)
+      ok &= utl::renumber(load.inod,nodeMap,true);
+
+  return ok;
+}
+
+
+bool HasPointLoads::assembleLoads (AlgEqSystem& eqSys, double t) const
+{
+  bool ok = true;
+  const size_t nrhs = eqSys.getNoRHS();
+  SystemVector* R = nrhs > 0 ? eqSys.getVector(nrhs-1) : nullptr;
+  if (!R) return ok; // No external force vector
+
+  for (const PointLoad& load : myLoads)
+  {
+    double P = (*load.func)(t);
+    int ldof = load.ldof;
+    eqSys.addScalar(P,ldof-1);
+    ok &= eqSys.getSAM().assembleSystem(*R,P,{load.inod,ldof});
+  }
+
+  return ok;
+}

HasPointLoads.h

@@ -0,0 +1,61 @@
+// $Id$
+//==============================================================================
+//!
+//! \file HasPointLoads.h
+//!
+//! \date Jun 12 2025
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Class representing a nodal point load container.
+//!
+//==============================================================================
+
+#ifndef _HAS_POINT_LOADS_H
+#define _HAS_POINT_LOADS_H
+
+#include <vector>
+#include <map>
+
+class ScalarFunc;
+class AlgEqSystem;
+namespace tinyxml2 { class XMLElement; }
+
+
+/*!
+  \brief Class representing a nodal point load container.
+*/
+
+class HasPointLoads
+{
+protected:
+  //! \brief Struct defining a nodal point load.
+  struct PointLoad
+  {
+    int         inod; //!< Node index
+    int         ldof; //!< Local DOF number
+    ScalarFunc* func; //!< The load magnitude as function of time
+    //! \brief Default constructor.
+    explicit PointLoad(int n = 0, int d = 0, ScalarFunc* f = nullptr)
+      : inod(n), ldof(d), func(f) {}
+  };
+
+  //! \brief Empty default constructor.
+  HasPointLoads() {}
+  //! \brief The destructor deletes the load functions.
+  virtual ~HasPointLoads();
+
+  //! \brief Parses a point load definition from an XML-element
+  bool parseLoad(const tinyxml2::XMLElement* elem);
+
+  //! \brief Renumbers the global node numbers of the point loads.
+  bool renumberLoadedNodes(const std::map<int,int>& nodeMap);
+
+  //! \brief Assemble right-hand-side vector contributions from the point loads.
+  bool assembleLoads(AlgEqSystem& eqSys, double t) const;
+
+private:
+  std::vector<PointLoad> myLoads; //!< Nodal point load container
+};
+
+#endif