From dfbde55f40c7b7036b6d5e0b17aab53961aa59f3 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Mon, 25 Aug 2025 16:02:53 +0200 Subject: [PATCH 1/7] Added: Use mkstemps to create the temporary input file --- main.C | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/main.C b/main.C index 0824725..9a7c0a5 100644 --- a/main.C +++ b/main.C @@ -29,9 +29,10 @@ #endif #include #include -#include -#include #include +#include +#include +#include namespace ASM { extern double Ktra, Krot; extern bool skipVTFmass; } @@ -204,13 +205,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 From 2156bd53d85a162f203fed115071240553768c88 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sat, 24 May 2025 08:52:27 +0200 Subject: [PATCH 2/7] Add class for superelement analysis --- SIMShellSup.C | 301 ++++++++++++++++++++++++++++++++++++++++++++++++++ SIMShellSup.h | 79 +++++++++++++ main.C | 77 +++++++------ 3 files changed, 422 insertions(+), 35 deletions(-) create mode 100644 SIMShellSup.C create mode 100644 SIMShellSup.h diff --git a/SIMShellSup.C b/SIMShellSup.C new file mode 100644 index 0000000..047a8e4 --- /dev/null +++ b/SIMShellSup.C @@ -0,0 +1,301 @@ +// $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 (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::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; +} + + +IntegrandBase* SIMShellSup::getMyProblem () const +{ + for (const std::pair& sub : mySubSim) + if (sub.second.sim) + return const_cast(sub.second.sim->getProblem()); + + return myProblem; +} + + +/*! + \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 + + 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], this->getMyProblem()->getField1Name(i).c_str(), + idBlock++, iStep); + + return ok ? idBlock : -4; +} diff --git a/SIMShellSup.h b/SIMShellSup.h new file mode 100644 index 0000000..e727124 --- /dev/null +++ b/SIMShellSup.h @@ -0,0 +1,79 @@ +// $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" + + +/*! + \brief Solution driver for linear elastic superelement FEM analysis. +*/ + +class SIMShellSup : public SIMsupel +{ +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 preprocess the FE substructures. + virtual bool preprocessB(); + + 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 Returns a pointer to the problem-specific data object. + IntegrandBase* getMyProblem() 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/main.C b/main.C index 9a7c0a5..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" @@ -79,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) @@ -259,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; } @@ -284,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(); @@ -296,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. @@ -360,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; @@ -471,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; @@ -480,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); From 2da16d34050bfadfda0657dcfc0f8bc84777480d Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 3 Jun 2025 16:31:55 +0200 Subject: [PATCH 3/7] Read through to EOF after successfully parsing the Nastran bulk data, to avoid that the read() method is invoked a second time on the same file. Ignore empty patches, i.e., a valid Nastran file but with no elements. Cosmetic changes in debug print when processing the RBE2/RBE3 elements. --- ASMu2DNastran.C | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) 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 } From 9347d1627c6acaad5b60031feff9a0d0db33123d Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 3 Jun 2025 16:55:06 +0200 Subject: [PATCH 4/7] Use a dummy/empty integrand if a negative value is passed as first argument to the SIMAndesShell constructor --- SIMAndesShell.C | 16 ++++++++++++++-- SIMAndesShell.h | 4 ++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/SIMAndesShell.C b/SIMAndesShell.C index 71d9e65..c9891f0 100644 --- a/SIMAndesShell.C +++ b/SIMAndesShell.C @@ -35,7 +35,7 @@ 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; @@ -51,8 +51,20 @@ SIMAndesShell::~SIMAndesShell () ElasticBase* SIMAndesShell::getIntegrand () { + // Dummy empty integrand class + class NoProblem : public IntegrandBase + { + public: + NoProblem() : IntegrandBase(0) {} + }; + if (!myProblem) - myProblem = new AndesShell(nss,modal,useBeams); + { + if (nss >= 0) + myProblem = new AndesShell(nss,modal,useBeams); + else + myProblem = new NoProblem(); + } return dynamic_cast(myProblem); } diff --git a/SIMAndesShell.h b/SIMAndesShell.h index 7f8e795..a6c073d 100644 --- a/SIMAndesShell.h +++ b/SIMAndesShell.h @@ -30,7 +30,7 @@ class SIMAndesShell : public SIMElasticity //! \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); + explicit SIMAndesShell(short int n = 0, bool m = false); //! \brief The destructor deletes the nodal point load functions. virtual ~SIMAndesShell(); @@ -128,7 +128,7 @@ class SIMAndesShell : public SIMElasticity 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 }; From 0b8969da0950b101d1ce1a68feeb4a9ccebcbcb2 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 3 Jun 2025 17:19:22 +0200 Subject: [PATCH 5/7] Added: Regression test for using FEDEM superelement matrices --- CMakeLists.txt | 9 ++++ Test/CQUAD04.dat | 20 ++++++++ Test/Q4-sup.reg | 79 ++++++++++++++++++++++++++++++++ Test/Q4-sup.vreg | 13 ++++++ Test/Q4-sup.xinp | 35 ++++++++++++++ Test/superelms/CQUAD04__B.fmx | Bin 0 -> 12142 bytes Test/superelms/CQUAD04__G.fmx | Bin 0 -> 334 bytes Test/superelms/CQUAD04__M.fmx | Bin 0 -> 1198 bytes Test/superelms/CQUAD04__S.fmx | Bin 0 -> 1198 bytes Test/superelms/CQUAD04__SAM.fsm | Bin 0 -> 3674 bytes 10 files changed, 156 insertions(+) create mode 100644 Test/CQUAD04.dat create mode 100644 Test/Q4-sup.reg create mode 100644 Test/Q4-sup.vreg create mode 100644 Test/Q4-sup.xinp create mode 100644 Test/superelms/CQUAD04__B.fmx create mode 100644 Test/superelms/CQUAD04__G.fmx create mode 100644 Test/superelms/CQUAD04__M.fmx create mode 100644 Test/superelms/CQUAD04__S.fmx create mode 100644 Test/superelms/CQUAD04__SAM.fsm 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/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..86f6ef1 --- /dev/null +++ b/Test/Q4-sup.reg @@ -0,0 +1,79 @@ +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 1236: (fixed) +Parsing + Parsing + Parsing +Parsing input file succeeded. +Equation solver: 0 +Number of Gauss points: 1 +Lagrangian basis functions are used +11. Shell superelement +Resolving Dirichlet boundary conditions + Constraining P1 in direction(s) 1236 + >>> SAM model summary <<< +Number of elements 1 +Number of nodes 2 +Number of dofs 12 +Number of unknowns 4 +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. +Solving the equation system ... + Condition number: 10.7647 + >>> Solution summary <<< +L2-norm : 7.41869e-06 +Max x-displacement : 1.8172e-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..4c03096 --- /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: 'primary solution 1' + Scalar: 'primary solution 2' + Scalar: 'primary solution 3' + Scalar: 'primary solution 4' + Scalar: 'primary solution 5' + Scalar: 'primary solution 6' diff --git a/Test/Q4-sup.xinp b/Test/Q4-sup.xinp new file mode 100644 index 0000000..143018f --- /dev/null +++ b/Test/Q4-sup.xinp @@ -0,0 +1,35 @@ + + + + + + + CQUAD04_.nas + + + + + + + CQUAD04.dat + + + 1 2 + + + + + + + + + + 1 + + + + + + + + diff --git a/Test/superelms/CQUAD04__B.fmx b/Test/superelms/CQUAD04__B.fmx new file mode 100644 index 0000000000000000000000000000000000000000..f77ec776105f98cfe7ec32b46d2c8a42d2cd982a GIT binary patch literal 12142 zcmai)d0fs}8^(Vql+jeSEHlDTBt0mZX=Oyp(j+9k$ufkHZANMcO~{tD$(pfbor%;~ z+Qd-yC50@BY!%f+hLN}5In#BV+w*?j=Z`0M=X399${vXp5&z;oos*lE7Tetfy`gj_z6vf}T za#X!%6{xKVe{~Dvv0p>@d*|90?YS%TRo&9N>B%+vH=_2$8NVxs{(}$d<398%b&|SU zoj8eg_t0l6hWz{v^?kx$<2?4;50rrWK!rALr?*vAtIP4gHeZS0d*%`V)G| z{CcroXe9UNIpJ56{^SzgF*{(0@6yxISNjrsED3%>Nq>I2Z1c^OHQC@@>m4{eE;b#0 zlB`aA@!dn=e;(n@IgkC@34biwV`t8Eq0e3?*EIO~>1C=#^0P{u*UqGHOjc4M^u`lE zbu_Y{HiUO1er^#zI!AI3Ha;u#nV)eP@bd@p<9X_3c!PiT!Oxj)nXbXUzoDLU>#jy; z-4pm9Z=QF^@}J#UXZnqILw{=Ohnbd!kHPy9{W7kP^`D6T={kq{IeEN*eyi@XpHjll z8d&g^LzLi|pRgnF^OpFb`oUZuejIPh*UR*(pM0k6{*Y%{1Yflb8os~#Svl%8ZlMd_ zPELeQ4B_`Ql=+Q*jK}_=OIA|cBK&hR#;4Wvzc_B&chl^JzDr?vi^ai0U*&olh}SEz+2>Ua zD#U%6nkxJ0Wh?9Jn#lY{jchx8J_t-=&vaF)KA!LGxTRrarvOq&O34aiBID6 zzw~&7IxnhY@98^IaJ@~JUG!=bEAD%}p?8%ZvjwkS@16&`(03@3^(!9B`j^%-eVlje zs^L3Linx!%?#q1aTA9CXmhk-jdcmg;{M3)47{kwTqVJP0>$@%H*L&=~=bskJHeC0} z*Ogvx8$X55vUa0}H{ZPr^=G#kk2-+xwYxJOdTO7dEyW= zhPsOI$tJY}181Lvz7grqnhvH1njW}}^Oo}d!0&a^A1}g>x03s*UPu1BeE4ZHrmkU~ zEkd7N=akxbT}?KnjP>mz`ZLdN$owhAx$v`Q)Pm90b3)*!*E+^yeJtT^erG(^*?Fgm z;3pz_?4pt)q0iR6KA=Bah@U=VhB;^HiTfmbv8BHl^_6U_HoD zPhaRWKbtGDe&OAoLsJd5;d(BPc^|FmnTC2>cvq|Q$1_l?`_tECVt=E-R>RwUohq#> zZQtX1=e=k8Y3$tUr;cM~eVc_$AM0Mk&(0xA$=@c$LjPx(=epVO)2eBU^?5x%;#J^j z71j-_WIr2;ehBv?&f^Ct=Rzg~{4gs6J{=qurNj=(1UIqI@B^N&4xw%kD)XN^%KYB` zjE7!l|KweZ6`?<%G~>L%DWT77jMT+_AMkzipwM3F*gve>vvx-|T|w>JBCX^6MqA-0 zjqtMwf60dN*uQ*E;NSs9Lcf=rm1TjR(9bb8?YVGk1NhmjYa7@6bR70Ky#LBG(Muau4=v{2~|G{&gL-$~n&wuk)F(S4k`P zCt>}`t9ikd8*iXazr4lTqJF4TjlRzUnRk3G^GmvIRe9_`<63Xt!gWG_jjqwmQw@c_ zs+*}>7ySus@yC$L_NTD_l9|$H$J|8J6(2hHo47-K&b~%?Z^D-~v{ZTQ4;*oN)0nH` z{nNjL{++R-6R@uO33JnfALVO-sa^C0_gSv3tZzW{^9C?|=$Q=(u-VY{80vMFE0?`7EWr7(lV>$~ z5hD6^E_Z8Hu-j4aUd};RW4Av5uliBK{(_&Up0fS`qW>g-={x!L+0gfK&`Ic=_UqWD z)tPA2CfSi^TYP&A^(exZa6XNlU;PwL_?L|mbu%6dea6|Whfl+1Yu%prKMnmRzpTGv zRVy8JNdC@EOOLIE{X(i!%rIFzux}ixZw1?SZDoN zxEg+TT_`%`WRwH_Q1|iy-+%Q7cOo_Jw^P!Tz9d4uuaHVoF5VLnajbMXHhRC{LAxuh3My+X6N?~xOx!$XqNzs z@!{g{)>6-L5pN=n%6^n^_%Zl;@QN*!f~WcsTp#;=h(5c%>L-gAX2)D|6X4@c`*jaJ z^pa6uSsR&Xf9@%CsyL52r48dz7Z9GEGxe5z&3Woq=OhI`kByFKW~W~e`ul&0sebIm z`K_jWn>DVR@Z-pY#B~R5Vm!`Y$$8Yn2tR+)(k=5R7s8L`+AxDjt20pdCH&+4XBG_B zi-w;oZRP%SIDDgH!`urv@9=MoN3Cln_vahJtH*^k$-((*ZMZ^wZno~~U{Q5J=(k*} zSh%nF6@Ds5EZsiRQuL=X)2ghgRW8mCTF-dYFMJt~x^0%stLMjdDvZ!ZopM zs7B~J4sNbf;TaD<4Y)q)vX%jhX9bAQ;qyHhkMoak|ETp_$-L@E`QahXSIhJX;`du` z-}d!xOc458A{#XQbLT?nO-!v*;rYAxz2zS29Wuc?3G1rf@-xYjzOwfX>c4t^7af0G z+_x{5F#R-jpZe%`g6X4kR+0E|L5wig-V=L;S_(JA+$I zO7_+lydT%c`D0ULKO6tme?;=+l_XEL<9RZCJSO~Ml2`X3yseh}+5g?4Wahv6sbCe$ z=b=yY^Bj`5+lR1x{(tfnB+vK$KYZt#p7WPmB*Tw=5XcaAH_-{{mUy`2}6TY6dJX{g;lwB`aJ`ep6lK)!} zKj(;kDB)d*pFF}RYMU2C9YuerpHn0cA3*Zz(wi)=hW{Db@^I&P=%=Qzd>;HsqOV8t z_Ee%jjO6(RB+t(#yuY@2LH>65IY<1s6a5r2kGMqiYwo*Hk|&Sjd2(9K`ZaGke;(Ce z&Z`RvU(i?1pZ!Vx{F3L-&{xkfHEruqbKWwNpU3jN9sK_M(SyBo1Yb`0GlU;j8gJIg z-WK{t2)}{n^U%LH^^TvLuLgV+nMWk@c?9^gByZ0le0{=O_{;r_BJ+z&Jy<>ueYG$C z`sVOcN%HV?o^ONi)K$*I{YV}@m++ZC%6WJm$-^h{JRJI}pG?m#@H2zt|7Upq58l~F z)^D97&llQcusk39;_GHE^z$mmc{HEzQ`UFj$jZ0i^KZ%Xh)HA~ zVcA2z-V)NE5rnrlmHVmsDouF`Kfzk&7h$u-++`ujx7|s;-TSeehv$;~yf@)%5C1d| z_kR0-^6*0q(f{+r56%B?Zj|*SiJxFHFE|n=`&awp9QG3Z3D-6+&|JZ~4#}&dNnRba zKuaDTdma0==HbQSzVat|dmPd48zJZazmh!PoaFhj+Vb*?KWE82 zVke(R;J#P;sXlh{{Mm_Diqf6r$@>0sp1g?X&*1kFKF(gwtKaba8N7#<{5gGnExeB+ zc-{^^n&j=Zn#+0n6T%-P{A61>&;On0^Ux0_`FyKWSw`lI#q0cx&o98ok@>|6KEFt- z>0c*4zW{%a%rDOJ`32V9`TPQPp|<%&ybk)aljqyu^~~iw+@I&+X*GUowPg7@_$M7? ze#qx656Aijo`<8>)slzjA8i0X5BR(Qyd}y1BhRsUL0XNU2tF?WA42$`9&BEK_2+zE zfI3UdyudJ64}N;`JRE#7$*Zk-9**_k>>?rME?Yt7Z?(KuBRw{56JnmUpqNZ){s2e zm*mfp9pt>)f$-`%Q46lXw^}Ff7N6sLko>u}8n=Ix&s&ka-BlI%=>Oq#Bl&zP$@Bk1 zc*VDFi)97*@RM;daJGN3_*^)c@NQ&&v7FB@oETR9beGH{O3D1>zn~CeOps|5&Yg_z6j0$*)9hO!Dxl zBoB{Imh<+#L_dJc3v~Fr0Dk;Pp3kn=G!f^s`lG)?|0eSSv!(L9Ac4;d@c!KLul|1l DpQ;)s literal 0 HcmV?d00001 diff --git a/Test/superelms/CQUAD04__G.fmx b/Test/superelms/CQUAD04__G.fmx new file mode 100644 index 0000000000000000000000000000000000000000..22b0c9dac7c826b87410aab825ff55b82b9fe060 GIT binary patch literal 334 zcmY#(b9Hg`RY)&NEXypZR7lG&N={WMOHD4xFDh10FcD$^f{7cuT&)fD46M04%0(`m z=?QiKiNe8?u)Q$0%2Hj{#PgaCFupq$^$njR;OaFX>JPxw!)T`pRD0?}OuqCz(L_}z z&E(Q`MMVS=}NWesh0SG2;>~gg>)HAT=`f=!1#}5Zz2aqTfG;A#M z*$d-X1W&*4>6p9&jQ>Auj-Jy-V>lluzTjcDH(dSeZ-*_uZxtV)_z6Uv-!&8*D6KA4lsEjozeW=9Hbt_y_qG_mUS6k z$m$rx6V>7F1FBc(koJJ{i@r#mGfT2WQLnLH748n8dN4qDZ&d2_KEvc6_Aq%6?GP>l zH}^Kgy$4|OxAGsqoTd9|Ka3Bg6WR`m?nMt@kUX+`1%ppC%GUpc+w%qzz78ODFd+35 z;lIDvsvmv%fnpC#8A|wq1d-j_{F~86!0i&;9+3Q$w{|GuaHK{W&S%+jy3$qgfIXJ< zfD*pw?){0G-X5f=`ohD_29n+mfb2&G=;;AHd_iiF-CKc{-V}1c>CFMj%nXJVRv`bN zq+e|50Xcku>QUT_CA~H1vS`A?31%O9`o)qS7?8sk-M!l}(;LTUZE;&{`3*V!qNfKe q;fvy43AFSEb3c%V9&Al22x#hkKW~16OWh&fUnoZDjLoQo4K#m4n2Iq;z;h@ zu#{ylS8E4c-(N`h9)vgwO6G5b+pBqT>bKon%i!u!)FFj0SRBc{24PKR4~xT4^(|P7 z8V(7)aDGxaSp7k4=>a)>(cSBC0yDieFg}I5-vN@|Q2hsVEP8rC4_~l8*Y{Je)x4K>$7dVoeW7;fwCx69+KU zTf?`mec19Fa{5J24_Lz&$-T?a(i_C}5E7Q(;9|h^3uD951B?%%;o%Dv1i2R`{~#2S F-T-iP@ZbOd literal 0 HcmV?d00001 diff --git a/Test/superelms/CQUAD04__SAM.fsm b/Test/superelms/CQUAD04__SAM.fsm new file mode 100644 index 0000000000000000000000000000000000000000..504cbe455316b3fc9a198e79d33ea117ec7dcb38 GIT binary patch literal 3674 zcmeI!)p8U;6o6q|g1bWk1h-%Tf`ueNfZ*=I-Q9vqxC(ATpMXoQcpco`!v9m-RbgQY zwn*-ns;_#sXIiGa&#CP&Ei*59S^l>Cdx?7l&k=CnEMp9pXwSv{RaR9Z^ zJci=OZ1W% zcjwv|P0$q0;5s`0mT-T}bt7z6uqJo__Ny*F-+*HC9^+#0{$`^>RJ zou#pCOo#Y3u0x%r@$Fng%p2>W&eCr%_Cp>p{zILmdBA+&z0dsKP-kgwFh_)ZVV($e zmgWoBC*%-whjrn5U5AiUe9!j(!gZGYhV^vK{C?NL?+@SQoE+0Rxqs&7_DDr%xR2f8 zUiLvhxS#28PwmgWaqV1x*WUiac&@8^<9Lqm`uUD?vz==Y+Bw(I&hub9_uO{wyX}f` z5!Z1E*Kh-uaTPal1-Ibu+wh$t*zOEMyR+PfcIUXa-FevV0&Ms5^HEYdLugE0gd7z%Su514Oyp*PGsY3Pf7 z=nr$yK%}A*I-?8BSKTlV^RWO6u?UNik0n@&Wmt|CScz3wjWt+{b$H6&Kfz->!b3d3 zF`gaA37o_!6yi0{Ug0HP;5nY5BF`$JGD`1bb$-`CO(elwXPz^!np4fC=FkRch^X~X rbRVPEEz$Ll+Q&rCLoVMo4&zb!d?tF%%WB`tD&C{=Ow{>F^qhPL{_Y`V literal 0 HcmV?d00001 From 8b7773c5b64b704080cae620efa146a6ca50cc4b Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 11 Jun 2025 16:40:16 +0200 Subject: [PATCH 6/7] Add base class HasPointLoads as the nodal point load container, used in both SIMAndesShell and SIMShellSup. Override method printProblem() in SIMAndesShell to print stuff only if the integrand member is an AndesShell instance. Add another test for superelement analysis with beam elements. --- HasPointLoads.C | 105 +++++++++++++++++++++++++++++++ HasPointLoads.h | 61 ++++++++++++++++++ SIMAndesShell.C | 57 ++++------------- SIMAndesShell.h | 25 +++----- SIMShellSup.C | 34 ++++++---- SIMShellSup.h | 16 +++-- Test/Q4-sup.reg | 21 +++++-- Test/Q4-sup.vreg | 12 ++-- Test/Q4-sup.xinp | 8 ++- Test/Q4andBeam-sup.reg | 91 +++++++++++++++++++++++++++ Test/Q4andBeam-sup.xinp | 38 +++++++++++ Test/Q4andBeam.dat | 21 +++++++ Test/superelms/Q4andBeam_B.fmx | Bin 0 -> 3646 bytes Test/superelms/Q4andBeam_G.fmx | Bin 0 -> 478 bytes Test/superelms/Q4andBeam_M.fmx | Bin 0 -> 2638 bytes Test/superelms/Q4andBeam_S.fmx | Bin 0 -> 2638 bytes Test/superelms/Q4andBeam_SAM.fsm | Bin 0 -> 1322 bytes 17 files changed, 398 insertions(+), 91 deletions(-) create mode 100644 HasPointLoads.C create mode 100644 HasPointLoads.h create mode 100644 Test/Q4andBeam-sup.reg create mode 100644 Test/Q4andBeam-sup.xinp create mode 100644 Test/Q4andBeam.dat create mode 100644 Test/superelms/Q4andBeam_B.fmx create mode 100644 Test/superelms/Q4andBeam_G.fmx create mode 100644 Test/superelms/Q4andBeam_M.fmx create mode 100644 Test/superelms/Q4andBeam_S.fmx create mode 100644 Test/superelms/Q4andBeam_SAM.fsm 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 c9891f0..7468f83 100644 --- a/SIMAndesShell.C +++ b/SIMAndesShell.C @@ -42,13 +42,6 @@ SIMAndesShell::SIMAndesShell (short int n, bool m) : nss(n), modal(m) } -SIMAndesShell::~SIMAndesShell () -{ - for (PointLoad& load : myLoads) - delete load.p; -} - - ElasticBase* SIMAndesShell::getIntegrand () { // Dummy empty integrand class @@ -70,6 +63,15 @@ ElasticBase* SIMAndesShell::getIntegrand () } +bool SIMAndesShell::printProblem () const +{ + if (dynamic_cast(myProblem)) + return this->SIMElasticity::printProblem(); + + return true; +} + + bool SIMAndesShell::parse (const tinyxml2::XMLElement* elem) { if (!strcasecmp(elem->Value(),"geometry")) @@ -165,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) { @@ -312,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); } @@ -338,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 a6c073d..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(short int n = 0, bool m = false); - //! \brief The destructor deletes the nodal point load functions. - virtual ~SIMAndesShell(); + + //! \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,19 +115,7 @@ 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 diff --git a/SIMShellSup.C b/SIMShellSup.C index 047a8e4..f6eed99 100644 --- a/SIMShellSup.C +++ b/SIMShellSup.C @@ -93,6 +93,12 @@ bool SIMShellSup::parse (const tinyxml2::XMLElement* elem) 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; @@ -140,6 +146,20 @@ bool SIMShellSup::preprocessB () } +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 { @@ -194,16 +214,6 @@ ElementBlock* SIMShellSup::tesselatePatch (size_t pidx) const } -IntegrandBase* SIMShellSup::getMyProblem () const -{ - for (const std::pair& sub : mySubSim) - if (sub.second.sim) - return const_cast(sub.second.sim->getProblem()); - - return myProblem; -} - - /*! \brief Static helper to write out scalar fields to VTF-file. */ @@ -292,10 +302,10 @@ int SIMShellSup::writeGlvS1 (const Vector& psol, int iStep, int& nBlock, // 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], this->getMyProblem()->getField1Name(i).c_str(), - idBlock++, iStep); + ok = vtf->writeSblk(sID[i], fieldNames[i], idBlock++, iStep); return ok ? idBlock : -4; } diff --git a/SIMShellSup.h b/SIMShellSup.h index e727124..9dc6210 100644 --- a/SIMShellSup.h +++ b/SIMShellSup.h @@ -15,13 +15,14 @@ #define _SIM_SHELL_SUP_H #include "SIMsupel.h" +#include "HasPointLoads.h" /*! \brief Solution driver for linear elastic superelement FEM analysis. */ -class SIMShellSup : public SIMsupel +class SIMShellSup : public SIMsupel, private HasPointLoads { public: //! \brief The constructor creates a dummy integrand with gravity only. @@ -36,10 +37,20 @@ class SIMShellSup : public SIMsupel //! \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 @@ -51,9 +62,6 @@ class SIMShellSup : public SIMsupel //! \brief Tesselates the superelement associated with specified patch. virtual ElementBlock* tesselatePatch(size_t pidx) const; - //! \brief Returns a pointer to the problem-specific data object. - IntegrandBase* getMyProblem() 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 diff --git a/Test/Q4-sup.reg b/Test/Q4-sup.reg index 86f6ef1..556bf55 100644 --- a/Test/Q4-sup.reg +++ b/Test/Q4-sup.reg @@ -39,7 +39,10 @@ Parsing Parsing Topology sets: support (1,1,4D) Parsing - Dirichlet code 1236: (fixed) + Dirichlet code 123456: (fixed) +Parsing + Parsing + Node 2 dof 1 Load: (expression) 1000\*sin(t) Parsing Parsing Parsing @@ -48,13 +51,15 @@ 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) 1236 + Constraining P1 in direction(s) 123456 >>> SAM model summary <<< Number of elements 1 Number of nodes 2 Number of dofs 12 -Number of unknowns 4 +Number of unknowns 6 Setting up for recovery of superelement "Q4" Resolving Dirichlet boundary conditions >>> SAM model summary <<< @@ -67,11 +72,15 @@ 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: 10.7647 + Condition number: 1791.19 >>> Solution summary <<< -L2-norm : 7.41869e-06 -Max x-displacement : 1.8172e-05 +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 diff --git a/Test/Q4-sup.vreg b/Test/Q4-sup.vreg index 4c03096..f192a9e 100644 --- a/Test/Q4-sup.vreg +++ b/Test/Q4-sup.vreg @@ -5,9 +5,9 @@ Step 1 27 nodes 16 elements Displacement: 'Displacement' - Scalar: 'primary solution 1' - Scalar: 'primary solution 2' - Scalar: 'primary solution 3' - Scalar: 'primary solution 4' - Scalar: 'primary solution 5' - Scalar: 'primary solution 6' + 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 index 143018f..8c401cf 100644 --- a/Test/Q4-sup.xinp +++ b/Test/Q4-sup.xinp @@ -14,15 +14,19 @@ CQUAD04.dat - 1 2 + 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/Q4andBeam_B.fmx b/Test/superelms/Q4andBeam_B.fmx new file mode 100644 index 0000000000000000000000000000000000000000..0b94cd813efcac5ef666242382cb94c9aa4f65cb GIT binary patch literal 3646 zcmYk8cQ_V&`^IHeMA9-+JPM)g^bl8s2yxr1Y-Qv&vPX7^+$ejiyO8Yt)McftLS&W9 ztc=Xkqu$=*cf5|@`QP`y>vNpvbzWa?B}D~Abxtc=S9{J|=C@sJJvslM#F;7nsC&-i ztF*{jVQISmzW%H3mL{6W^lnHyImxIu@C}x;&WL8(zJjfip1s1(S#V?sry9Oipxjr> zJ(79jIK@shnt3(`?C)sVgq@c`&5M)k%}hOz0=!@9->iaX-a0X)`gTm22p*{{oyRKo zzG##BX&CC4Q+ICs0&b;zuQw|bVZ21;s$jvgBSvgouTW z4Q7u4D6^E9W3bPJ&HmlAS*HH5>x^Z?G3j`e7fa6N6?_i<2js^kMw79;u~ggN>m96! z&-q=*O@yrjgN3xk5%6VkwLU2?4bL9mZ6c#RjI!_a3g{XeAfms}xTSFb+MW9M@wd7` z=BmRza;*nE6b*EKGn>>vVa<>9;DZs8`p%aw(tU4F>j6#;+o>C6B5QypJ{=1K0 z4&^h>R$aeIc;Gve*Ik;AI~+t>*N|jaRX8;pZV+Ia#$GM%U_4XW9R}(HU#_eDLT}_+ep?&9E>&VJ!-`4Tjab z=Mk1mLwza)x>5V}6g3ND1FHB33`h&80fnEX0+$yb%2Nnw*AEE+2eSPyO|woHIXizC_@~bdJ zr`$%i|1OSno^aQ??+kn-$%r4E-k4z(WN&bo2>kDw$kc}hVK&1$Nkpj<8Ghu{UR_PV z<`{?+vbO>Aj7TxLJOYO9H%~Zw^cGTg^U^0KeS!tHr?-Wv%P?cOqU_;p6nx-dpRP#P z1A(``=M~s>vG_Tg2EVU4)}}o>`fG4t$B(Yp$F^MLuAJ05_*a@bUSs}D4(|Ai_U~{VQ0aR|W<}QDaKvC-iV#v!( zJbpJp@x9>( zWPFEXz)D2Gu-gx4__(k4#PopvH^uwwqG{m0x^hq^%o`%*#_t;+c0|*Fug-JiAv`dO zyL`^nqD@EI_lsYKft_G|O=U0z`gylcIl<5M1S4v}d-#`!g@lyB5R!bkkXchgtD+W_U9ZpuZ*p!b?v$Xd)9j;LDyxX#>&b+ok% zaB2zOf3B}i->!rYryr)K<88<(_%*b6Kns>~`@OmNbRp!y2sed$Eoy&Qe@K$w1T$9A z;CZ7mC?pt&S?JrrM3T^j6eR-Aght=3UC@Cjv7VAoYs(<^gl%5fVH5K2WYI1JxPdfF zL`mu?yB&@-5utziHJMS*{oDrmgCy=M6BSLKX~z~NHkp$300^4MRokyo2i>8DgC&vj z*c?I~wh=Fb31U}1iB2e?W_Lk>&{iYfRXFg|slW@e_cKj-yOe=5IZ#Qj<~l|hGNkdd z@}gD3XigG=4>cbXIeJwdV?yDgL(iZ7iAjGZbULFNE~T5A%70VDZ$lcQcEAdIjCQ}j zA#@T=r?N@I^N3F=DF+n=yRqi4Ymqk?8sMy;kGT50GXAfxbYA9dGRg>mDeVgynm!wB zZ`jWFyu5(soWikV)16o{ZF!|2uMfqNhb`O`ng)wRRM zc-wIHAB~u7j()J~$p8D&_giRJml)gsMg>~ZbT1c55peSEy~sQF#_&VneI;J2df5pYWKw;^k1^64S$}G zdrN5g!i9(6aInAyYLyxYR};UUZ|Q*U3(o|q>hmGpVZeDyBoJJ-t&EK>=Ao=^NRsh^ z5vZGX%W3#h4lbke_Kv6Hz+`ZO3J322^J+LD{`M_UQ+{0^IBL_y=d#Lm5`n5@wF!hi4RPOd0?de7xuXa~9e zu%U&h2AmIMJK%5MhTjsk3n^-Ip@4X1v+-mn{F>%AGc%b%yEI*fpRuoiPWK9xGUW>x z?G4eejDCYsJ*8>6^u<^b%W_74w``NnbOivVw@EMfH9c97m(4>>j;R!$lL@%NYYfea5YTxlgKPh;R z;D$0xoU<~2i&kIAKXjGEx;WQKW5d?@li?3j1stH&rRxuiGML_D6Wl;JLW2{Y6d z5+aDE*!J~P>6l3qI95MMi2oe|D?QjV)8CYP1kl4Vi&)okP7kil)A;)V475MUQ?Q=aQ9M{&D?r}3JW_>Xp!AVs7cJc^?`L%*Sb@eJ0&L?RjT z63(_3(ud+tzj*iG-R0;d`7*PD&2q;NMmMjdFK>0gtjbfmz}^A13jMg7`Ru^g;Np3y zSH2KN)5R2>VvFtq%hW8hH!)#)rtr3$5o#8Gv$)cC3-3R&F3*p0#y8DYNi$}?FfPlN zSn$XR^a_PC4eN>MPH0`5Y`BVq*Q_1IBmxGX_Wf*HZHz-_ucoZ{+X6vC)^2~JIaK#v zd66Jvfo1d44lys)V7*VHBj*$W%}yqDv^>&-cGcGzRP2A|>X!eZyUPff>Hc-y9b&P= zk(7o{R-)OMV_N%iNy)zBO|yZ!`;ycvpzWWAADOV%Uo zUn_7RhZa!H7{D~Y@ysF|r;T62* z-a!d*Ic|O|T=C#k+W*MU#R+&kdsq6kvyom#tepOy685HD_?=B91ug`0FNthnypXH6 zLST;vSxH8J$~-5O3f`i7tC5Xw$Xb+Km(T9xO4i`T?Y0Lwm~^m*fA?4{26=e^!uSf{e`R0j1MvCdJKS@9|l`6RL-3gfHXNn{b zWn$PA`-j(znb39OU&ao3JEVF-isHC#g)Ey#ZRPD_VNxXg&b|ga{Pg$b=!5B*_-pl? JAD?MF{0E}I)PMj0 literal 0 HcmV?d00001 diff --git a/Test/superelms/Q4andBeam_G.fmx b/Test/superelms/Q4andBeam_G.fmx new file mode 100644 index 0000000000000000000000000000000000000000..8878185c626c202da2148c5f00a37ce3e49470aa GIT binary patch literal 478 zcmV<40U`b)MnyzLO(17-VRmVBc_3zQa${v6c4cF9Z*p@WAT$yH0002d>UZorF)lDW z3c5@p;~6KCKrm09Bl*j1K*q05dBe8Ezx^dM@1YXaKPGlAR}J#;zkHcNrQ>sCz$^yE zl>PRRK>5T@h0Muqz%XtdwpqxsKk|b1iKxlYznQsi3h@&1zpyAU)~K&%Kp=AC;P`=X zK%mP}eDY+mzl!jGJT0ZRKcD!g^lG2Szh`|$^w**aKnd`A1P;zbz%kR*!mwSOKpso= zrs~U&KnUvj8mNYGK(OZw7~YdIK=jQnzK%joz($#9RDAY0z}{lcexKx8(6 zBNj#lz}M^ffYx8XK#~+?>XlZKKtmuH zHme+Qz=}W2-rob1z>*v9XV3;wKseby*89hAz)I-Gg%Rpez)t)w+@kJnKso;J`x_(? UK+a#xP7y6uKs0+5a`?DXKvj3{(*OVf literal 0 HcmV?d00001 diff --git a/Test/superelms/Q4andBeam_M.fmx b/Test/superelms/Q4andBeam_M.fmx new file mode 100644 index 0000000000000000000000000000000000000000..b4b9b1e415f9f0fae67963c8d819b7ade84217f7 GIT binary patch literal 2638 zcmYk8d05W*8pgZI(ViA!2rW|Th%iZ|-}APYw2&lYo0hjY)-u{??`=xE7CFvXMjTVF zoEFzO>3tB9G#W2D>Xa-IiL{7f&ir$Co!=k7zwhgQ?)!6pja;0aoLB4b+!hh>=PzSw531Lv7_Bex#sTl zSY)K)J}sQT#$CDR_Ddk=IoH^=K0CsrTa1JMElZ;3k`8uX8^$xe(|)T`wP5+!U5aK6 z{v4Y<@^qJM4PC80u|a3VlVU7(tWWogqsv~~hXW(6$<$lV(ONQ#Uu1m##$vMxJ^$y8 zwp#`k3>_*(3Vr^yNo1#NE;H&4psWi2p@Qo96dfI_Yg*^brZ=7BKYQ4+qE%Jfq8$l* zXW6@HC4bvXUtbA@7v+6fLHs`XVEaqyC4C?6wlwBywFVvcti?p}>OIqXG->trSJF8L zw^5%Yc*d}uCT)ErF)^u><6{>dhZouZAzX?`K9+G@jo}A$bOh%do-?nI%3}^!N}mMN zbj`K0ZE30Of6du!BwdlO+GOm^Td;>84CW7s`o5!I&p9Z$?B2`1pZx9&)%)^ji$y_n zZwzauDvw+4SYv0O|6yB3s3 z-4LkZUjFoji-OoQphwlu_JQy^J|IKL^>ErtXhov_ z6Ax{a9M~%)Rqu~7eew%3C@GiBA;m=6*r_%k-P_sR`tpVrM_v8gl^i{fX-u{vrF$9I z>Xs%dj!)+mexY}BtDXtkIpt?uW3LEPpP$dZc|1a>&WTp5YEmWi1vy1Na2NOi9a(R$ zmpyGdMUfTRIsT7+v0-?Dh1$B$>;^m(fg)+7k&LDSrA&Ge zL{A)ZzJBnYLjGTt9>{T3VA=BBOB=6^39_Ccd;KNDQt$)4MPHCp-am^0Esp?E6Xtj)5QOE&76-A|JR5{D6+)vG#^nA6!Xp|AiOtZ-$X? zW}r({!hBBD&eS>i!HJF)o({i8mSp+5C1|~sh-W3w_@!fNr&Ja#iWl3Mv+JtXup@e= z6FddSzz_5meL+r<58MTQK!=m^h`MKYCPnl%myYVYvzkin{lpM49j)kEzQV(ov|1k@ z>0A@U5-pj!&woEj z_1M$T^sU5b+=%kCH@Tf@ljpP2+BR{YAqkViil+T=R)`uM_o>m-q?3ERwO5X`O0h5C zI(Q0>fgk8C`huJyAGizrfR3`DqGFNV7FzMw^S3L{s*uYfangZ75zEPCdby-GNeyRe zd`$E(6&&J@9W<`vR_mr{Ea)f!?Ao$SLxHyTA|V zkS6}@{>EuQSRDM9q@;TtLQCW_ng3@N@czK=V{fr@*cWgeJO#(V5A+s&K~9km+y#C> zN4fUqlu7cs)KORPuj1b|Cf-4Ke_;2qx7a!C3%CxRf@9zZdW*gwr^pBH0zaT5#{Q|^ Jl9%7m{{UP@>zx1q literal 0 HcmV?d00001 diff --git a/Test/superelms/Q4andBeam_S.fmx b/Test/superelms/Q4andBeam_S.fmx new file mode 100644 index 0000000000000000000000000000000000000000..eae960451c216148e10b5c0a0db3bfd3a740e487 GIT binary patch literal 2638 zcmX|@i$4{1AIDEh7ah+NHI24*!a24Raw(S-=PUQ-lEPd|E2d3ySy>kHxMdz@{W}?gaS?&_mm=>m;#;Mcg;T7w~okK}xl62fF!&M)`- zn~^N<*7z#3K;jv>_encanTFx1&gP<n z`8^@*t21|NX6vezI*ult3ivRZYcH>CRZ$7}h`fW^f(!bVvEDakG4V@`=xlJM&uD|3 zisY(R`FCpBB6p`5fjPH1F26FbwAq@~ZF32}`%IIvMpu=Ymsn8XiH@;NgL3vtC;aDy zf1Y3_*4jb|rzh6;kp~6&<1)Gx|neekEj@L#eGLtC}qmlzu=*vaSD6cXlI9 z$on39rikZ!xnEu%oW$I|n|AK_=UVpYf|iC^<2Gh^ngiGwP!u{iT%^7P1%g!r7KNAPRwgOZNx3lkaUk9 zNehu$vdrbecb~bMG28`yD0LiNx17u>Nbui&QOKrZ zZ6l`5jy+Wg_4N=l(S3oZyeE@B+5L|oo&0cetiSI2hhjC$DLWUPI$%i+-Itm|UHO#z za>VzoDAwu-0hk(tSdHREXot`ouxC{J%j^gJ-CDqNXTXG{e=l=q$?^21ldtpM+;Ao-``)|c*ikk-tYm?YC9nLgvkz8#V+xDszf_;4bh3I(mPp zo7C{?B~GP-*Bu*2mfNs%VdSu!^$xbLsJ%EVm8D$1t}#9-;+!b_G{tl!{VZ={A*2WnYj}HHAg49xpaA!w@t(96k1{I#~ zIkWlm2UMKC2&7w$U1gKiDr8z4G7z!Gg%stbwJuj3m0Sl;!7=ayy+vP;Q{)48fgjM( z^Kw$n)BP4*)MG1mNOM?nV!2l7jF4P@Om_HB6Q6iz47fY|o=L|%-69UBZ)VcVTXo;b zUPIUya2-4a$G{Kt7JWfZkq_Jjen3aSp8l`OBk$6~!o?S>Jx)>kg~7IQFH3TjB`@=z zx1Ob+R4gzNAV_pY1(= z^Tw1d>zvL#zu_ve*1BU8l2T(P&hh?wK|#85Z?SXO7jPXs1;@Y-^cH+krmGLWPjMD8j*u?_p!IwIqVC#4xWNz;0Jn(z96T_ z2krttpkv~8LVx#WM=(}+_{nZq%E-v%ewVPW0J3E`{!R{Y#pQ@Eoxtcv88-_tB!!o#KAOY7f^TQ+#H$^ZF znFQCDQn7jCI0LK~+}me@rJprde*bGtYY)LV%z}5xK@se|0`9R7_xFSM^6tK;4nB)M z<6iMj_}KFkpRVyAai<=*Mw@xIruWVxFbZll0WnAGsf+rAM(X1`qE2bO)GRb~_F`Q* zg>&#sYkKAwOoFp_CNZNU;uTnhJZ!=-m>ZW;aeSxFw%r~)4h!GRKJ!>_`McQB^sV)z zPCTtA`|4lUle~Txy|3}gI&44=@?ig4P=IYHf|{#C{D;)TZ*~smVF4CFJ=I8^)h`Qb ztB%gjdBvP!4$f!?c3}_p;Q-Y89eaaUc!3OS&x!9;tJB%t-l1EKI^DsZ|Dx9q7i>FA literal 0 HcmV?d00001 From 979935e561ed5298f8533a0743bdeb74baf5da3a Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 18 Jun 2025 09:48:40 +0200 Subject: [PATCH 7/7] Avoid NaN in debug print when zero denominator --- ANDES/wavgmConstrEqn.f90 | 43 ++++++++++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 15 deletions(-) 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