diff --git a/ANDES/wavgmConstrEqn.f90 b/ANDES/wavgmConstrEqn.f90 index 4e265ee..f5afae3 100644 --- a/ANDES/wavgmConstrEqn.f90 +++ b/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 diff --git a/ASMu2DNastran.C b/ASMu2DNastran.C index 89ee342..909a141 100644 --- a/ASMu2DNastran.C +++ b/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 } diff --git a/CMakeLists.txt b/CMakeLists.txt index 4c7c5ad..4db2b3a 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/HasPointLoads.C b/HasPointLoads.C new file mode 100644 index 0000000..c60481f --- /dev/null +++ b/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& 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; +} diff --git a/HasPointLoads.h b/HasPointLoads.h new file mode 100644 index 0000000..7692d27 --- /dev/null +++ b/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 +#include + +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& nodeMap); + + //! \brief Assemble right-hand-side vector contributions from the point loads. + bool assembleLoads(AlgEqSystem& eqSys, double t) const; + +private: + std::vector myLoads; //!< Nodal point load container +}; + +#endif diff --git a/SIMAndesShell.C b/SIMAndesShell.C index 71d9e65..7468f83 100644 --- a/SIMAndesShell.C +++ b/SIMAndesShell.C @@ -35,26 +35,40 @@ bool SIMAndesShell::replRBE3 = false; namespace Elastic { extern double time; } -SIMAndesShell::SIMAndesShell (unsigned short int n, bool m) : nss(n), modal(m) +SIMAndesShell::SIMAndesShell (short int n, bool m) : nss(n), modal(m) { nsd = 3; nf.front() = 6; } -SIMAndesShell::~SIMAndesShell () +ElasticBase* SIMAndesShell::getIntegrand () { - for (PointLoad& load : myLoads) - delete load.p; + // Dummy empty integrand class + class NoProblem : public IntegrandBase + { + public: + NoProblem() : IntegrandBase(0) {} + }; + + if (!myProblem) + { + if (nss >= 0) + myProblem = new AndesShell(nss,modal,useBeams); + else + myProblem = new NoProblem(); + } + + return dynamic_cast(myProblem); } -ElasticBase* SIMAndesShell::getIntegrand () +bool SIMAndesShell::printProblem () const { - if (!myProblem) - myProblem = new AndesShell(nss,modal,useBeams); + if (dynamic_cast(myProblem)) + return this->SIMElasticity::printProblem(); - return dynamic_cast(myProblem); + return true; } @@ -153,27 +167,7 @@ bool SIMAndesShell::parse (const tinyxml2::XMLElement* elem) { IFEM::cout <<" Parsing "<< std::endl; - int inod = 0, ldof = 0; - ScalarFunc* f = nullptr; - utl::getAttribute(child,"node",inod); - utl::getAttribute(child,"dof",ldof); - - if (inod > 0 && ldof > 0 && ldof <= 6) - { - std::string type("constant"); - utl::getAttribute(child,"type",type); - - IFEM::cout <<"\tNode "<< inod <<" dof "<< ldof <<" Load: "; - if (type == "constant") - { - f = new ConstantFunc(atof(child->FirstChild()->Value())); - IFEM::cout << (*f)(0.0) << std::endl; - } - else - f = utl::parseTimeFunc(child->FirstChild()->Value(),type); - - myLoads.emplace_back(inod,ldof,f); - } + this->parseLoad(child); } else if (!strcasecmp(child->Value(),"material") && elInt) { @@ -300,11 +294,7 @@ bool SIMAndesShell::renumberNodes (const std::map& nodeMap) if (spr.inod > 0) ok &= utl::renumber(spr.inod,nodeMap,true); - for (PointLoad& load : myLoads) - if (load.inod > 0) - ok &= utl::renumber(load.inod,nodeMap,true); - - return ok; + return ok && this->renumberLoadedNodes(nodeMap); } @@ -326,18 +316,7 @@ bool SIMAndesShell::assembleDiscreteTerms (const IntegrandBase* itg, for (const DOFspring& spr : mySprings) ok &= K->add(spr.coeff,mySam->getEquation(spr.inod,spr.ldof)); - const size_t nrhs = myEqSys->getNoRHS(); - SystemVector* R = nrhs > 0 ? myEqSys->getVector(nrhs-1) : nullptr; - if (R) // Assemble external nodal point loads - for (const PointLoad& load : myLoads) - { - double P = (*load.p)(time.t); - int ldof = load.ldof; - myEqSys->addScalar(P,ldof-1); - ok &= mySam->assembleSystem(*R,P,{load.inod,ldof}); - } - - return ok; + return ok && this->assembleLoads(*myEqSys,time.t); } diff --git a/SIMAndesShell.h b/SIMAndesShell.h index 7f8e795..a92dcf1 100644 --- a/SIMAndesShell.h +++ b/SIMAndesShell.h @@ -14,8 +14,9 @@ #ifndef _SIM_ANDES_SHELL_H #define _SIM_ANDES_SHELL_H -#include "SIMElasticity.h" #include "SIM2D.h" +#include "SIMElasticity.h" +#include "HasPointLoads.h" class DataExporter; @@ -24,15 +25,19 @@ class DataExporter; \brief Driver class for FE analysis using the ANDES shell elements. */ -class SIMAndesShell : public SIMElasticity +class SIMAndesShell : public SIMElasticity, private HasPointLoads { public: //! \brief Default constructor. //! \param[in] n Number of consecutive solutions in core (0 = linear analysis) //! \param[in] m If \e true, a modal linear dynamics simulation is performed - explicit SIMAndesShell(unsigned short int n = 0, bool m = false); - //! \brief The destructor deletes the nodal point load functions. - virtual ~SIMAndesShell(); + explicit SIMAndesShell(short int n = 0, bool m = false); + + //! \brief Prints out problem-specific data to the log stream. + virtual bool printProblem() const; + + //! \brief Returns the name of this simulator. + virtual std::string getName() const { return "AndesShell"; } //! \brief Enables the element matrix cache for reaction force calculation. void initForSingleStep() { if (!myRFset.empty()) this->initLHSbuffers(); } @@ -110,25 +115,13 @@ class SIMAndesShell : public SIMElasticity : inod(n), ldof(d), coeff(c) {} }; - //! \brief Struct defining a nodal point load. - struct PointLoad - { - int inod; //!< Node index - int ldof; //!< Local DOF number - ScalarFunc* p; //!< Load magnitude - //! \brief Default constructor. - explicit PointLoad(int n = 0, int d = 0, ScalarFunc* f = nullptr) - : inod(n), ldof(d), p(f) {} - }; - std::vector mySprings; //!< Global DOF springs - std::vector myLoads; //!< Nodal point loads RealArray myReact; //!< Nodal reaction forces std::string myRFset; //!< Node set for calculation of reaction forces std::string myPath; //!< Relative path of the patch file - unsigned short int nss; //!< Number of consequtive solution states in core + short int nss; //!< Number of consequtive solution states in core bool modal; //!< Modal dynamics simulation flag }; diff --git a/SIMShellSup.C b/SIMShellSup.C new file mode 100644 index 0000000..f6eed99 --- /dev/null +++ b/SIMShellSup.C @@ -0,0 +1,311 @@ +// $Id$ +//============================================================================== +//! +//! \file SIMShellSup.C +//! +//! \date May 24 2025 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Solution driver for linear elastic superelement FEM analysis. +//! +//============================================================================== + +#include "SIMShellSup.h" +#include "SIMAndesShell.h" +#include "SAM.h" +#include "ASMsupel.h" +#include "HasGravityBase.h" +#include "ElementBlock.h" +#include "Utilities.h" +#include "Tensor.h" +#include "VTF.h" +#include "IFEM.h" +#include "tinyxml2.h" + + +SIMShellSup::SIMShellSup (const char* hd, bool fd) : SIMsupel(hd,-6), fixDup(fd) +{ + // Integrand class for superelement analysis with gravity loads + class MyProblem : public HasGravityBase + { + public: + MyProblem() : HasGravityBase(3) { npv = 6; } + void printLog() const { IFEM::cout <<"Linear superelement"<< std::endl; } + }; + + myProblem = new MyProblem(); +} + + +SIMShellSup::~SIMShellSup() +{ + for (std::pair& sub : mySubSim) + { + delete sub.second.sim; + delete sub.second.blk; + } +} + + +bool SIMShellSup::parse (const tinyxml2::XMLElement* elem) +{ + if (strcasecmp(elem->Value(),"elasticity")) + return this->SIMsupel::parse(elem); + + std::string supId; + SIMgeneric* supEl = nullptr; + if (utl::getAttribute(elem,"supId",supId)) + { + // This block defines the FE model of a superelement. + // Create a separate SIM-object for it here for the recovery operation, but + // with no integrand since recovery is based on an externally reduced model. + IFEM::cout <<" Parsing FE model for superelement \""<< supId <<"\"\n"; + supEl = new SIMAndesShell(-1); + } + + const tinyxml2::XMLElement* child = elem->FirstChildElement(); + for (; child; child = child->NextSiblingElement()) + if (!strcasecmp(child->Value(),"supernodes") && supEl) + { + std::string setName; + IFEM::cout <<" Parsing for \""<< supId <<"\".\n"; + std::vector& superNodes = mySubSim[supId].superNodes; + if (utl::getAttribute(child,"set",setName)) + superNodes = supEl->getNodeSet(setName); + if (superNodes.empty()) + IFEM::cout <<" ** Empty or non-existing node set \""<< setName <<"\"."; + else + { + IFEM::cout <<"\t"<< superNodes.front(); + for (size_t i = 1; i < 10 && i < superNodes.size(); i++) + IFEM::cout <<" "<< superNodes[i]; + if (superNodes.size() > 10) + IFEM::cout <<" ... ("<< superNodes.size() <<" nodes)."; + } + IFEM::cout << std::endl; + } + else if (!strcasecmp(child->Value(),"gravity") && !myModel.empty()) + IFEM::cout <<" ** The superelement definition will be ignored" + <<" unless you specify before the" + <<" block of the superelement."<< std::endl; + + else if (myProblem->parse(child)) + continue; // parsed + + else if (!strcasecmp(child->Value(),"nodeload") && child->FirstChild()) + { + IFEM::cout <<" Parsing "<< std::endl; + + this->parseLoad(child); + } + else if (supEl && !supEl->parse(child)) // parse the superelement + { + delete supEl; + std::cerr <<" *** Failure."<< std::endl; + return false; + } + + if (supEl) + { + mySubSim[supId].sim = supEl; + IFEM::cout <<" FE model \""<< supId <<"\" loaded."<< std::endl; + } + + return true; +} + + +bool SIMShellSup::preprocessB () +{ + // Set up links to the underlying FE model for each superelement + for (SuperElm& sup : mySups) + { + std::map::const_iterator sit = mySubSim.find(sup.id); + if (sit != mySubSim.end()) sup.sim = sit->second.sim; + } + + // Preprocess the superelement FE models for the recovery process. + // We need to maintain the original nodal ordering (only eliminating holes + // in the sequence), otherwise the recovery will fail if the patch contains + // beam elements which are placed in a separate patch. + bool preserved = SIMbase::preserveNOrder; + bool ok = SIMbase::preserveNOrder = true; + for (std::pair& sub : mySubSim) + if (sub.second.sim) + { + IFEM::cout <<"\nSetting up for recovery of superelement \"" + << sub.first <<"\""<< std::endl; + ok &= (sub.second.sim->preprocess({},fixDup) && + sub.second.sim->setMode(SIM::RECOVERY)); + } + + IFEM::cout <<"\nSuperelement recovery setup done."<< std::endl; + SIMbase::preserveNOrder = preserved; + return ok; +} + + +bool SIMShellSup::renumberNodes (const std::map& nodeMap) +{ + return (this->SIMsupel::renumberNodes(nodeMap) && + this->renumberLoadedNodes(nodeMap)); +} + + +bool SIMShellSup::assembleDiscreteTerms (const IntegrandBase* itg, + const TimeDomain& time) +{ + return this->assembleLoads(*myEqSys,time.t); +} + + +bool SIMShellSup::recoverInternalDOFs (const ASMbase* pch, SuperElm& sup, + const Vector& supSol) const +{ + const ASMsupel* supch = dynamic_cast(pch); + if (supch && supch->hasRecovery()) + // Recover the internal displacement state for this superelement + return supch->recoverInternals(supSol,sup.sol); + + return this->SIMsupel::recoverInternalDOFs(pch,sup,supSol); +} + + +ElementBlock* SIMShellSup::tesselatePatch (size_t pidx) const +{ + if (pidx >= mySups.size() || !mySups[pidx].sim) + return this->SIMsupel::tesselatePatch(pidx); // use simplified tesselation + + FEmodel& sub = const_cast(this)->mySubSim[mySups[pidx].id]; + + // Create an ElementBlock containing the shell patches in this superelement. + // Normally, there should only be one such patch, but... + if (sub.sim && !sub.blk) + for (ASMbase* pch : sub.sim->getFEModel()) + if (pch && !pch->empty() && pch->getNoParamDim() > 1) // skip beam patches + { + bool ok = true; + if (!sub.blk) + { + sub.blk = new ElementBlock(4); + ok = pch->tesselate(*sub.blk,nullptr); + } + else + { + ElementBlock newblk(4); + if ((ok = pch->tesselate(newblk,nullptr))) + sub.blk->merge(newblk,false); + } + if (!ok) + { + delete sub.blk; + return sub.blk = nullptr; + } + } + + if (!sub.blk) + return nullptr; + + // Make a copy of sub.blk and apply the superelement transformation to it + ElementBlock* supblk = new ElementBlock(*sub.blk); + supblk->transform(mySups[pidx].MVP); + return supblk; +} + + +/*! + \brief Static helper to write out scalar fields to VTF-file. +*/ + +static bool writeFields (const Matrix& field, int geomID, + int& nBlock, std::vector& sID, VTF* vtf) +{ + for (size_t j = 1; j <= field.rows(); j++) + if (!vtf->writeNres(field.getRow(j), ++nBlock, geomID)) + return false; + else if (j <= sID.size()) + sID[j-1].push_back(nBlock); + else + sID.push_back({nBlock}); + + return true; +} + + +int SIMShellSup::writeGlvS1 (const Vector& psol, int iStep, int& nBlock, + double, const char*, int idBlock, int, bool) +{ + if (adm.dd.isPartitioned() && adm.getProcId() != 0) + return 0; + else if (psol.empty()) + return 0; + + VTF* vtf = this->getVTF(); + if (!vtf) return -99; + + // Recover internal displacements for all superelements + if (!this->recoverInternalDOFs(psol)) + return -9; + + IntVec vID; + std::vector sID; + sID.reserve(this->getNoFields()); + + int geomID = this->getStartGeo(); + for (size_t pidx = 0; pidx < myModel.size(); pidx++) + { + Matrix field, subfield; + + if (mySups[pidx].sim) + { + for (ASMbase* pch : mySups[pidx].sim->getFEModel()) + if (pch && !pch->empty() && pch->getNoParamDim() > 1) + { + // Extract displacement vector for this sub-patch + Vector pchvec; + pch->extractNodalVec(mySups[pidx].sol, pchvec, + mySups[pidx].sim->getSAM()->getMADOF()); +#if INT_DEBUG > 2 + std::cout <<"\nSolution vector for sub-patch "<< pch->idx+1 << pchvec; +#endif + + // Evaluate the internal displacement field on this sub-patch + if (!pch->evalSolution(subfield,pchvec,(const int*)nullptr)) + return -5; + + if (field.empty()) + field = subfield; + else + field.augmentCols(subfield); + } + } + else // Evaluate displacement field on supernodes only + if (!myModel[pidx]->evalSolution(field, mySups[pidx].sol, opt.nViz)) + return -1; + + if (msgLevel > 1) + IFEM::cout <<"Writing primary solution for patch " + << myModel[pidx]->idx+1 <<" ("<< field.rows() + <<","<< field.cols() <<")"<< std::endl; + + // Output as vector field + if (!vtf->writeVres(field, ++nBlock, ++geomID, this->getNoSpaceDim())) + return -2; + else + vID.push_back(nBlock); + + // Output as scalar fields + if (!writeFields(field, geomID, nBlock, sID, vtf)) + return -3; + } + + // Write result block identifications + + const char* fieldNames[6] = { "u_x", "u_y", "u_z", "r_x", "r_y", "r_z" }; + bool ok = vID.empty() || vtf->writeDblk(vID,"Displacement",idBlock,iStep); + for (size_t i = 0; i < sID.size() && !sID[i].empty() && ok; i++) + ok = vtf->writeSblk(sID[i], fieldNames[i], idBlock++, iStep); + + return ok ? idBlock : -4; +} diff --git a/SIMShellSup.h b/SIMShellSup.h new file mode 100644 index 0000000..9dc6210 --- /dev/null +++ b/SIMShellSup.h @@ -0,0 +1,87 @@ +// $Id$ +//============================================================================== +//! +//! \file SIMShellSup.h +//! +//! \date May 24 2025 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Solution driver for linear elastic superelement FEM analysis. +//! +//============================================================================== + +#ifndef _SIM_SHELL_SUP_H +#define _SIM_SHELL_SUP_H + +#include "SIMsupel.h" +#include "HasPointLoads.h" + + +/*! + \brief Solution driver for linear elastic superelement FEM analysis. +*/ + +class SIMShellSup : public SIMsupel, private HasPointLoads +{ +public: + //! \brief The constructor creates a dummy integrand with gravity only. + //! \param[in] hd Sub-simulator heading + //! \param[in] fd If \e true, merge duplicated FE nodes on patch interfaces + SIMShellSup(const char* hd, bool fd); + //! \brief The destructor deletes the FE substructure data. + virtual ~SIMShellSup(); + +protected: + using SIMsupel::parse; + //! \brief Parses a data section from an XML element. + virtual bool parse(const tinyxml2::XMLElement* elem); + + //! \brief Performs some pre-processing tasks on the FE model. + //! \details This method is overridden to print out the problem definition. + virtual void preprocessA() { this->printProblem(); } + //! \brief Performs some pre-processing tasks on the FE model. + //! \details This method is overridden to preprocess the FE substructures. + virtual bool preprocessB(); + + //! \brief Renumbers the global node numbers of the springs and point loads. + virtual bool renumberNodes(const std::map& nodeMap); + + //! \brief Assembles the DOF springs and nodal point loads, if any. + virtual bool assembleDiscreteTerms(const IntegrandBase* itg, + const TimeDomain& time); + + using SIMsupel::recoverInternalDOFs; + //! \brief Recovers the internal DOFs from static condensation. + //! \param[in] pch The patch associated with the superelement \a sup + //! \param sup The superelement to recover internal DOFs for + //! \param[in] supSol Supernode solution values + virtual bool recoverInternalDOFs(const ASMbase* pch, SuperElm& sup, + const Vector& supSol) const; + + //! \brief Tesselates the superelement associated with specified patch. + virtual ElementBlock* tesselatePatch(size_t pidx) const; + + //! \brief Writes primary solution for a given time step to the VTF-file. + //! \param[in] psol Primary solution vector + //! \param[in] iStep Time step identifier + //! \param nBlock Running result block counter + //! \param[in] idBlock Starting value of result block numbering + virtual int writeGlvS1(const Vector& psol, int iStep, int& nBlock, + double, const char*, int idBlock, int, bool); + +private: + //! \brief Struct representing a FE substructure. + struct FEmodel + { + SIMgeneric* sim = nullptr; //!< Underlying FE model of the substructure + ElementBlock* blk = nullptr; //!< Tesselated FE model for visualization + std::vector superNodes; //!< List of external supernode numbers + }; + + std::map mySubSim; //!< FE substructure container + + bool fixDup; //!< If \e true, merge duplicated FE nodes on patch interfaces +}; + +#endif diff --git a/Test/CQUAD04.dat b/Test/CQUAD04.dat new file mode 100644 index 0000000..5d71c6b --- /dev/null +++ b/Test/CQUAD04.dat @@ -0,0 +1,20 @@ +# Test of superelement from FEDEM in IFEM-ShellEx. +# ################################################ +# Notice the file suffixes _S.fmx, _M.fmx, etc., are _not_ specified. +# The fmx-matrices are produced by executing the FEDEM FE model Reducer as follows: +# $ fedem_reducer -linkfile CQUAD04_.nas +# +# Stiffness matrix: +Kx 12 12 superelms/CQUAD04_ +# Mass matrix: +Mx 12 12 superelms/CQUAD04_ +# Gravity force vectors: +gx 12 3 superelms/CQUAD04_ +# Displacement recovery matrix: +Bx 126 12 superelms/CQUAD04_ +# SAM data for substructure recovery: +Sx superelms/CQUAD04_ +# Global coordinates of the supernodes: +G 2 +0.0 0.0 -1.0 +0.0 0.0 0.0 diff --git a/Test/Q4-sup.reg b/Test/Q4-sup.reg new file mode 100644 index 0000000..556bf55 --- /dev/null +++ b/Test/Q4-sup.reg @@ -0,0 +1,88 @@ +Q4-sup.xinp -dense + +Input file: Q4-sup.xinp +Evaluation time for property functions: 1 +1. Shell superelement +Parsing input file Q4-sup.xinp +Parsing +Parsing + Parsing FE model for superelement "Q4" + Parsing + Reading data file CQUAD04_.nas +Total number of nodes: 27 +Number of shell elements: 16 +Number of constraint elements: 2 +Number of other elements: 0 +Model extension (diameter): 1.0198 +Pre-defined node sets: 1 + "ASET" 2 nodes + Reading patch 1 + Parsing + Parsing + Gravitation vector: 0 -9.81 0 + FE model "Q4" loaded. +Parsing + Parsing + Reading data file CQUAD04.dat + Test of superelement from FEDEM in IFEM-ShellEx. + ################################################ + Notice the file suffixes _S.fmx, _M.fmx, etc., are _not_ specified. + The fmx-matrices are produced by executing the FEDEM FE model Reducer as follows: + $ fedem_reducer -linkfile CQUAD04_.nas + Stiffness matrix: + Mass matrix: + Gravity force vectors: + Displacement recovery matrix: + SAM data for substructure recovery: + Global coordinates of the supernodes: + Reading patch 1 + Parsing + Topology sets: support (1,1,4D) + Parsing + Dirichlet code 123456: (fixed) +Parsing + Parsing + Node 2 dof 1 Load: (expression) 1000\*sin(t) +Parsing + Parsing + Parsing +Parsing input file succeeded. +Equation solver: 0 +Number of Gauss points: 1 +Lagrangian basis functions are used +11. Shell superelement +Problem definition: +Linear superelement +Resolving Dirichlet boundary conditions + Constraining P1 in direction(s) 123456 + >>> SAM model summary <<< +Number of elements 1 +Number of nodes 2 +Number of dofs 12 +Number of unknowns 6 +Setting up for recovery of superelement "Q4" +Resolving Dirichlet boundary conditions + >>> SAM model summary <<< +Number of elements 16 +Number of nodes 27 +Number of dofs 162 +Number of constraints 24 +Number of unknowns 138 +Superelement recovery setup done. +Processing integrand associated with code 0 +Assembling interior matrix terms for P1 +Done. +Total external load: Sum(Fex) = 841.471 -1534.28 0 +Solving the equation system ... + Condition number: 1791.19 + >>> Solution summary <<< +L2-norm : 2.84878e-05 +Max X-displacement : 2.08317e-05 +Max Y-displacement : 5.46372e-05 +Max x-displacement : 7.34201e-05 +Max y-displacement : 3.04787e-05 + \* Reading stiffness matrix file: superelms/CQUAD04__S.fmx 144 + \* Reading mass matrix file: superelms/CQUAD04__M.fmx 144 + \* Reading gravity force vectors file: superelms/CQUAD04__G.fmx 36 + \* Reading disk matrix file: superelms/CQUAD04__B.fmx 1512 + \* Reading SAM data file: superelms/CQUAD04__SAM.fsm diff --git a/Test/Q4-sup.vreg b/Test/Q4-sup.vreg new file mode 100644 index 0000000..f192a9e --- /dev/null +++ b/Test/Q4-sup.vreg @@ -0,0 +1,13 @@ +Q4-sup.xinp -dense + +Step 1 + Element block 1 + 27 nodes + 16 elements + Displacement: 'Displacement' + Scalar: 'u_x' + Scalar: 'u_y' + Scalar: 'u_z' + Scalar: 'r_x' + Scalar: 'r_y' + Scalar: 'r_z' diff --git a/Test/Q4-sup.xinp b/Test/Q4-sup.xinp new file mode 100644 index 0000000..8c401cf --- /dev/null +++ b/Test/Q4-sup.xinp @@ -0,0 +1,39 @@ + + + + + + + CQUAD04_.nas + + + + + + + CQUAD04.dat + + + 1 + + + + + + + + + + 1000*sin(t) + + + + 1 + + + + + + + + diff --git a/Test/Q4andBeam-sup.reg b/Test/Q4andBeam-sup.reg new file mode 100644 index 0000000..5766f26 --- /dev/null +++ b/Test/Q4andBeam-sup.reg @@ -0,0 +1,91 @@ +Q4andBeam-sup.xinp -dense + +Input file: Q4andBeam-sup.xinp +Evaluation time for property functions: 1 +1. Shell superelement +Parsing input file Q4andBeam-sup.xinp +Parsing +Parsing + Parsing FE model for superelement "Q4andBeam" + Parsing + Reading data file Q4andBeam.nas +Total number of nodes: 8 +Number of shell elements: 2 +Number of beam elements: 2 +Number of constraint elements: 1 +Number of other elements: 1 +Model extension (diameter): 3.20156 + Created beam patch 1 with 3 nodes +Pre-defined node sets: 1 + "ASET" 2 nodes + Created shell patch 2 with 3 elements + Parsing + FE model "Q4andBeam" loaded. +Parsing + Parsing + Reading data file Q4andBeam.dat + Test of superelement from FEDEM in IFEM-ShellEx. + ################################################ + Notice the file suffix, _S.fmx, _M.fmx, etc., is _not_ specified. + The fmx-matrices are produced by executing the FEDEM FE model Reducer as follows: + $ fedem_reducer -linkfile Q4andBeam.nas -extNodes 7 + Stiffness matrix: + Mass matrix: + Gravity force vectors: + Displacement recovery matrix: + SAM data for substructure recovery: + Global coordinates of the supernodes: + Reading patch 1 + Parsing + Topology sets: support (1,1,4D) + Parsing + Dirichlet code 123456: (fixed) +Parsing + Parsing + Node 3 dof 3 Load: 1e+06 +Parsing + Parsing + Parsing +Parsing input file succeeded. +Equation solver: 0 +Number of Gauss points: 1 +Lagrangian basis functions are used +11. Shell superelement +Problem definition: +Linear superelement +Resolving Dirichlet boundary conditions + Constraining P1 in direction(s) 123456 + >>> SAM model summary <<< +Number of elements 1 +Number of nodes 3 +Number of dofs 18 +Number of unknowns 6 +Setting up for recovery of superelement "Q4andBeam" +Resolving Dirichlet boundary conditions + >>> SAM model summary <<< +Number of elements 5 +Number of nodes 8 +Number of dofs 48 +Detected 1 elements out of order, reordering... +Number of constraints 6 +Number of unknowns 42 +Superelement recovery setup done. +Processing integrand associated with code 0 +Assembling interior matrix terms for P1 +Done. +Total external load: Sum(Fex) = 0 0 1e+06 +Solving the equation system ... + Condition number: 3067.94 + >>> Solution summary <<< +L2-norm : 0.167182 +Max X-displacement : 0.010265 +Max Y-displacement : 0.00232858 +Max Z-displacement : 0.449137 +Max x-displacement : 0.0273557 +Max y-displacement : 0.547576 +Max z-displacement : 0.025921 + \* Reading stiffness matrix file: superelms/Q4andBeam_S.fmx 324 + \* Reading mass matrix file: superelms/Q4andBeam_M.fmx 324 + \* Reading gravity force vectors file: superelms/Q4andBeam_G.fmx 54 + \* Reading disk matrix file: superelms/Q4andBeam_B.fmx 450 + \* Reading SAM data file: superelms/Q4andBeam_SAM.fsm diff --git a/Test/Q4andBeam-sup.xinp b/Test/Q4andBeam-sup.xinp new file mode 100644 index 0000000..2f3e706 --- /dev/null +++ b/Test/Q4andBeam-sup.xinp @@ -0,0 +1,38 @@ + + + + + + + Q4andBeam.nas + + + + + + Q4andBeam.dat + + + 1 2 + + + + + + + + + + 1.0e6 + + + + 1 + + + + + + + + diff --git a/Test/Q4andBeam.dat b/Test/Q4andBeam.dat new file mode 100644 index 0000000..dd652c3 --- /dev/null +++ b/Test/Q4andBeam.dat @@ -0,0 +1,21 @@ +# Test of superelement from FEDEM in IFEM-ShellEx. +# ################################################ +# Notice the file suffix, _S.fmx, _M.fmx, etc., is _not_ specified. +# The fmx-matrices are produced by executing the FEDEM FE model Reducer as follows: +# $ fedem_reducer -linkfile Q4andBeam.nas -extNodes 7 +# +# Stiffness matrix: +Kx 18 18 superelms/Q4andBeam +# Mass matrix: +Mx 18 18 superelms/Q4andBeam +# Gravity force vectors: +gx 18 3 superelms/Q4andBeam +# Displacement recovery matrix: +Bx 25 18 superelms/Q4andBeam +# SAM data for substructure recovery: +Sx superelms/Q4andBeam +# Global coordinates of the supernodes: +G 3 +0.0 0.0 0.0 +0.0 1.0 0.0 +3.5 0.5 0.0 diff --git a/Test/superelms/CQUAD04__B.fmx b/Test/superelms/CQUAD04__B.fmx new file mode 100644 index 0000000..f77ec77 Binary files /dev/null and b/Test/superelms/CQUAD04__B.fmx differ diff --git a/Test/superelms/CQUAD04__G.fmx b/Test/superelms/CQUAD04__G.fmx new file mode 100644 index 0000000..22b0c9d Binary files /dev/null and b/Test/superelms/CQUAD04__G.fmx differ diff --git a/Test/superelms/CQUAD04__M.fmx b/Test/superelms/CQUAD04__M.fmx new file mode 100644 index 0000000..0ef954a Binary files /dev/null and b/Test/superelms/CQUAD04__M.fmx differ diff --git a/Test/superelms/CQUAD04__S.fmx b/Test/superelms/CQUAD04__S.fmx new file mode 100644 index 0000000..0185807 Binary files /dev/null and b/Test/superelms/CQUAD04__S.fmx differ diff --git a/Test/superelms/CQUAD04__SAM.fsm b/Test/superelms/CQUAD04__SAM.fsm new file mode 100644 index 0000000..504cbe4 Binary files /dev/null and b/Test/superelms/CQUAD04__SAM.fsm differ diff --git a/Test/superelms/Q4andBeam_B.fmx b/Test/superelms/Q4andBeam_B.fmx new file mode 100644 index 0000000..0b94cd8 Binary files /dev/null and b/Test/superelms/Q4andBeam_B.fmx differ diff --git a/Test/superelms/Q4andBeam_G.fmx b/Test/superelms/Q4andBeam_G.fmx new file mode 100644 index 0000000..8878185 Binary files /dev/null and b/Test/superelms/Q4andBeam_G.fmx differ diff --git a/Test/superelms/Q4andBeam_M.fmx b/Test/superelms/Q4andBeam_M.fmx new file mode 100644 index 0000000..b4b9b1e Binary files /dev/null and b/Test/superelms/Q4andBeam_M.fmx differ diff --git a/Test/superelms/Q4andBeam_S.fmx b/Test/superelms/Q4andBeam_S.fmx new file mode 100644 index 0000000..eae9604 Binary files /dev/null and b/Test/superelms/Q4andBeam_S.fmx differ diff --git a/Test/superelms/Q4andBeam_SAM.fsm b/Test/superelms/Q4andBeam_SAM.fsm new file mode 100644 index 0000000..22af201 Binary files /dev/null and b/Test/superelms/Q4andBeam_SAM.fsm differ diff --git a/main.C b/main.C index 0824725..fedd74f 100644 --- a/main.C +++ b/main.C @@ -13,6 +13,7 @@ #include "IFEM.h" #include "SIMenums.h" +#include "SIMShellSup.h" #include "SIMShellModal.h" #include "SIMAndesSplit.h" #include "ElasticityUtils.h" @@ -29,9 +30,10 @@ #endif #include #include -#include -#include #include +#include +#include +#include namespace ASM { extern double Ktra, Krot; extern bool skipVTFmass; } @@ -78,6 +80,7 @@ namespace ASM { extern double Ktra, Krot; extern bool skipVTFmass; } state-dependent property functions \arg -no-vtfmass : Skip the sphere geometries for 1-noded mass elements \arg -Kbush : Spring stiffness(es) connecting RBE3 reference nodes to ground + \arg -sup : Perform superelement analysis (linear static) */ int main (int argc, char** argv) @@ -204,13 +207,18 @@ int main (int argc, char** argv) if (bdfile) { // Create a temporary input file referring to the Nastran bulk data file - static const char* tmpname = "/tmp/tmp.xinp"; + static char tmpname[24] = "/tmp/tmp_XXXXXX.xinp"; + int fd = mkstemps(tmpname,5); infile = const_cast(tmpname); - std::ofstream xinp(infile); - xinp <<"\n" - <<" "<< bdfile <<"\n" - <<"\n"; - + std::string xinp = "\n" + " " + std::string(bdfile) + + "\n\n"; + if (write(fd,xinp.c_str(),xinp.size()) < 0) + { + std::cerr <<" *** Failed to write temporary file "<< tmpname << std::endl; + return false; + } + close(fd); if (IFEM::getOptions().format >= 0) { // Set default vtf file name @@ -253,7 +261,7 @@ int main (int argc, char** argv) "[-no-vtfmass]]", "[-hdf5 [] [-dumpNodeMap]]","[-fixDup []]", "[-refsol ]","[-noBeams]","[-noEccs]","[-noSets]", - "[-noRBE3]","[-split]"}); + "[-noRBE3]","[-split]","[-sup]"}); delete prof; return 0; } @@ -278,8 +286,8 @@ int main (int argc, char** argv) IFEM::cout <<"\nSpecified boundary conditions are ignored"; if (fixDup) IFEM::cout <<"\nCo-located nodes will be merged," - <<" using comparison tolerance "<< Vec3::comparisonTolerance - << std::endl; + <<" using comparison tolerance "<< Vec3::comparisonTolerance; + IFEM::cout << std::endl; #ifdef HAS_FFLLIB FFl::initAllElements(); @@ -290,13 +298,16 @@ int main (int argc, char** argv) // Create the simulation model std::vector modes; DataExporter* writer = nullptr; - SIMAndesShell* model = nullptr; + SIMAndesShell* shell = nullptr; + SIMoutput* model = nullptr; if (modal) - model = new SIMShellModal(modes); + model = shell = new SIMShellModal(modes); else if (splitM) - model = new SIMAndesSplit(nstates); + model = shell = new SIMAndesSplit(nstates); + else if (args.dim == 4) + model = new SIMShellSup("Shell superelement",fixDup); else - model = new SIMAndesShell(nstates); + model = shell = new SIMAndesShell(nstates); // Lambda function for cleaning the heap-allocated objects on termination. // To ensure that their destructors are invoked also on simulation failure. @@ -354,8 +365,8 @@ int main (int argc, char** argv) } // Open HDF5 result database, if requested - if (model->opt.dumpHDF5(infile)) - writer = model->getHDF5writer(displ.front(),dumpNodeMap); + if (shell && shell->opt.dumpHDF5(infile)) + writer = shell->getHDF5writer(displ.front(),dumpNodeMap); Vector load; RealArray Fex; @@ -465,7 +476,7 @@ int main (int argc, char** argv) if (model->opt.format >= 0) { - int geoBlk = 0, nBlock = 0; + int geoBlk = 0, nBlock = 0, idBlock = 1; size_t iStep = 1, nStep = 0; double time = 0.0; Vector data; @@ -474,31 +485,33 @@ int main (int argc, char** argv) if (!model->writeGlvG(geoBlk,infile)) return terminate(12); - // Write sensor locations, if any - if (nodalR != 'm' && !model->writeGlvLoc(locfiles,nodalR,geoBlk)) - return terminate(12); - - // Write surface pressures, if any - if (!model->writeGlvT(iStep,geoBlk,nBlock)) - return terminate(13); - - // Write Dirichlet boundary conditions - if (!model->writeGlvBC(nBlock)) - return terminate(13); - - // Write global node and element numbers as scalar fields - int idBlock = 7; - if (!model->writeGlvNo(nBlock,idBlock,18)) - return terminate(14); - - // Write shell thickness as scalar field - model->getShellThicknesses(data); - if (!model->writeGlvE(data,iStep,nBlock,"Shell thickness",idBlock++,true)) - return terminate(15); - - // Write load vector to VTF-file - if (vizRHS && !model->writeGlvV(load,"Load vector",iStep,nBlock,1)) - return terminate(20); + if (shell) + { + // Write sensor locations, if any + if (nodalR != 'm' && !shell->writeGlvLoc(locfiles,nodalR,geoBlk)) + return terminate(12); + + // Write surface pressures, if any + if (!shell->writeGlvT(iStep,geoBlk,nBlock)) + return terminate(13); + + // Write Dirichlet boundary conditions + if (!shell->writeGlvBC(nBlock)) + return terminate(13); + + idBlock = 7; // Write global node and element numbers as scalar fields + if (!shell->writeGlvNo(nBlock,idBlock,18)) + return terminate(14); + + // Write shell thickness as scalar field + shell->getShellThicknesses(data); + if (!shell->writeGlvE(data,iStep,nBlock,"Shell thickness",idBlock++,true)) + return terminate(15); + + // Write load vector to VTF-file + if (vizRHS && !shell->writeGlvV(load,"Load vector",iStep,nBlock,1)) + return terminate(20); + } // Write solution fields to VTF-file model->setMode(SIM::RECOVERY);