From 6b322da4b64033715329e28ecd8dd03c9e3698e3 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 17 Mar 2022 07:09:18 +0100 Subject: [PATCH] Changed: Generalize ReactionsOnly class for internal force calculation. Added: SIMbase::assembleForces() and IntegrandBase::setSecondaryInt(). Added: virtual method for calculation/print of interface force resultants. --- src/ASM/IntegrandBase.h | 2 ++ src/ASM/ReactionsOnly.C | 33 ++++++++++++++++++--------------- src/ASM/ReactionsOnly.h | 23 +++++++++++++---------- src/ASM/TimeDomain.h | 3 +++ src/SIM/SIMbase.C | 24 ++++++++++++++++++++++++ src/SIM/SIMbase.h | 24 ++++++++++++++++++------ src/SIM/SIMoutput.h | 2 ++ 7 files changed, 80 insertions(+), 31 deletions(-) diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index ba3f9021..1b81a49f 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -84,6 +84,8 @@ public: virtual void initResultPoints(double, bool = false) {} //! \brief Initializes the global node number mapping for current patch. virtual void initNodeMap(const std::vector&) {} + //! \brief Assigns a secondary integral to be computed (for reaction forces). + virtual void setSecondaryInt(GlobalIntegral* = nullptr) {} //! \brief Returns the system quantity to be integrated by \a *this. virtual GlobalIntegral& getGlobalInt(GlobalIntegral* gq) const; diff --git a/src/ASM/ReactionsOnly.C b/src/ASM/ReactionsOnly.C index 3a2ade3c..0772be66 100644 --- a/src/ASM/ReactionsOnly.C +++ b/src/ASM/ReactionsOnly.C @@ -18,30 +18,34 @@ #include "ProcessAdm.h" -ReactionsOnly::ReactionsOnly (Vector& rf, const SAM* sam, - const ProcessAdm& adm, Vector* sf) +ReactionsOnly::ReactionsOnly (const SAM* sam, const ProcessAdm& adm, + Vector* rf, Vector* sf) : mySam(sam), myAdm(adm), R(rf), S(sf) { - mySam->initForAssembly(b,&R); + mySam->initForAssembly(b,R); } void ReactionsOnly::initialize (bool) { b.init(); - R.fill(0.0); + if (R) + R->fill(0.0); } bool ReactionsOnly::finalize (bool) { - if (myAdm.dd.isPartitioned()) - myAdm.allReduceAsSum(R); + if (R) + { + if (myAdm.dd.isPartitioned()) + myAdm.allReduceAsSum(*R); - R *= -1.0; + (*R) *= -1.0; #if SP_DEBUG > 2 - std::cout <<"\nReaction forces:"<< R; + std::cout <<"\nReaction forces:"<< *R; #endif + } if (!S) return true; @@ -60,7 +64,7 @@ bool ReactionsOnly::finalize (bool) bool ReactionsOnly::assemble (const LocalIntegral* elmObj, int elmId) { const ElmMats* elMat = dynamic_cast(elmObj); - if (elMat && mySam->assembleSystem(b,elMat->getRHSVector(),elmId,&R)) + if (elMat && mySam->assembleSystem(b,elMat->getRHSVector(),elmId,R)) return true; std::cerr <<" *** ReactionsOnly::assemble: Failure for element "<< elmId @@ -69,23 +73,22 @@ bool ReactionsOnly::assemble (const LocalIntegral* elmObj, int elmId) } -/*! - Returns \e true if the patch \a pidx have any dirichlet boundary conditions. -*/ - bool ReactionsOnly::haveContributions (size_t pidx, const PropertyVec& pvec) const { - auto&& isConstrained = [this,pidx](const Property& p) + // Lambda function checking if a property have force contribution on a patch. + auto&& contributed = [this,pidx](const Property& p) { if (p.patch != pidx) return false; else if (S) return p.pcode == Property::OTHER; + else if (!R) + return false; return (p.pcode >= Property::DIRICHLET && p.pcode <= Property::DIRICHLET_ANASOL); }; - return std::find_if(pvec.begin(),pvec.end(),isConstrained) != pvec.end(); + return std::find_if(pvec.begin(),pvec.end(),contributed) != pvec.end(); } diff --git a/src/ASM/ReactionsOnly.h b/src/ASM/ReactionsOnly.h index 06f1c1c6..f087cd38 100644 --- a/src/ASM/ReactionsOnly.h +++ b/src/ASM/ReactionsOnly.h @@ -7,7 +7,7 @@ //! //! \author Knut Morten Okstad / SINTEF //! -//! \brief Global integral sub-class for reaction force integration. +//! \brief Global integral class for reaction- and interface force calculation. //! //============================================================================== @@ -19,7 +19,7 @@ /*! - \brief Class for assembly of reaction forces only. + \brief Class for assembly of reaction- and interface forces. \details This class can be used for linear problems which normally consists of a single assembly loop. For non-linear problems, the reaction forces for @@ -27,18 +27,21 @@ the solution vector of the previous iteration is used. In linear problems we therefore need a separate assembly loop where only the reaction forces are calculated. This class is provided to facilitate such calculations. + + The class can also be used to calculate interface force resultants, + for boundaries not associated with dirichlet conditions. */ class ReactionsOnly : public GlobalIntegral { public: //! \brief The constructor initializes the data members. - //! \param[in] rf The reaction force vector to be assembled //! \param[in] sam Data for FE assembly management //! \param[in] adm Parallell processing administrator + //! \param[in] rf Reaction force vector to be assembled //! \param[in] sf Internal force vector to be assembled - ReactionsOnly(Vector& rf, const SAM* sam, const ProcessAdm& adm, - Vector* sf = nullptr); + ReactionsOnly(const SAM* sam, const ProcessAdm& adm, + Vector* rf = nullptr, Vector* sf = nullptr); //! \brief Empty destructor. virtual ~ReactionsOnly() {} @@ -49,10 +52,10 @@ public: //! \brief Adds a LocalIntegral object into a corresponding global object. //! \param[in] elmObj The local integral object to add into \a *this. - //! \param[in] elmId Global number of the element associated with elmObj + //! \param[in] elmId Global number of the element associated with \a elmObj virtual bool assemble(const LocalIntegral* elmObj, int elmId); - //! \brief Returns \e false if no contributions from the specified patch. + //! \brief Returns \e true if the patch \a pidx have any force contributions. //! \param[in] pidx 1-based patch index //! \param[in] pvec Physical property mapping virtual bool haveContributions(size_t pidx, @@ -62,9 +65,9 @@ private: const SAM* mySam; //!< Data for FE assembly management const ProcessAdm& myAdm; //!< Parallel processing administrator - Vector& R; //!< Nodal reaction forces - StdVector b; //!< Dummy right-hand-side vector - Vector* S; //!< Internal force vector + StdVector b; //!< Internal right-hand-side vector used in the assembly process + Vector* R; //!< Nodal reaction forces + Vector* S; //!< Nodal internal forces }; #endif diff --git a/src/ASM/TimeDomain.h b/src/ASM/TimeDomain.h index 3eb0db99..bd6f1d05 100644 --- a/src/ASM/TimeDomain.h +++ b/src/ASM/TimeDomain.h @@ -31,6 +31,9 @@ struct TimeDomain //! \brief Default constructor. explicit TimeDomain(int i = 0, bool f = true) : it(i), first(f) { t = dt = dtn = CFL = 0.0; } + //! \brief Constructor for linear problems with fixed (non-zero) time. + explicit TimeDomain(double t0) : t(t0), it(0), first(true) + { dt = dtn = CFL = 0.0; } }; #endif diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index e08726fb..67281e7f 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -22,6 +22,7 @@ #endif #include "IntegrandBase.h" #include "AlgEqSystem.h" +#include "ReactionsOnly.h" #include "SystemMatrix.h" #include "LinSolParams.h" #include "EigSolver.h" @@ -2479,3 +2480,26 @@ void SIMbase::dumpEqSys () utl::zero_print_tol = old_tol; } } + + +bool SIMbase::assembleForces (const Vector& solution, double t0, + Vector* R, Vector* S) +{ + if (!myProblem) return false; + + // Assign secondary integral to the integrand + myProblem->setSecondaryInt(new ReactionsOnly(mySam,adm,R,S)); + + // Temporarily nullify myEqSys such that the secondary integral will be used + AlgEqSystem* tmpEqSys = myEqSys; + myEqSys = nullptr; + + // Integrate the reaction and/or interface forces + bool ok = this->setMode(SIM::RHS_ONLY) && this->assembleSystem(t0,{solution}); + + // Release the secondary integral and restore the system matrix pointer + myProblem->setSecondaryInt(); + myEqSys = tmpEqSys; + + return ok; +} diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index 5352f4f6..db776bb5 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -278,11 +278,12 @@ public: bool newLHSmatrix = true, bool poorConvg = false); //! \brief Administers assembly of the linear equation system. + //! \param[in] t0 Time for evaluation of time-dependent property functions //! \param[in] pSol Primary solution vectors in DOF-order //! //! \details Use this version for linear/stationary problems only. - bool assembleSystem(const Vectors& pSol = Vectors()) - { return this->assembleSystem(TimeDomain(),pSol); } + bool assembleSystem(double t0 = 0.0, const Vectors& pSol = Vectors()) + { return this->assembleSystem(TimeDomain(t0),pSol); } //! \brief Extracts the assembled load vector for inspection/visualization. //! \param[out] loadVec Global load vector in DOF-order @@ -464,20 +465,31 @@ public: //! \param[in] ncv Number of Arnoldi vectors (see ARPack documentation) //! \param[in] shift Eigenvalue shift //! \param[out] solution Computed eigenvalues and associated eigenvectors - //! \param[in] iA Index of system matrix \b A in \a myEqSys->A - //! \param[in] iB Index of system matrix \b B in \a myEqSys->A + //! \param[in] iA Index of system matrix \b A in \ref myEqSys + //! \param[in] iB Index of system matrix \b B in \ref myEqSys bool systemModes(std::vector& solution, int nev, int ncv, int iop, double shift, size_t iA = 0, size_t iB = 1); //! \brief Performs a generalized eigenvalue analysis of the assembled system. //! \param[out] solution Computed eigenvalues and associated eigenvectors - //! \param[in] iA Index of system matrix \b A in \a myEqSys->A - //! \param[in] iB Index of system matrix \b B in \a myEqSys->A + //! \param[in] iA Index of system matrix \b A in \ref myEqSys + //! \param[in] iB Index of system matrix \b B in \ref myEqSys bool systemModes(std::vector& solution, size_t iA = 0, size_t iB = 1) { return this->systemModes(solution,opt.nev,opt.ncv,opt.eig,opt.shift,iA,iB); } + //! \brief Returns whether reaction forces are to be computed or not. + virtual bool haveBoundaryReactions(bool = false) const { return false; } + + //! \brief Assembles reaction and interface forces for specified boundaries. + //! \param[in] solution Current primary solution vector + //! \param[in] t0 Current time (for time-dependent loads) + //! \param[out] R Nodal reaction force container + //! \param[out] S Nodal interface force container + bool assembleForces(const Vector& solution, double t0, + Vector* R, Vector* S = nullptr); + // Post-processing methods // ======================= diff --git a/src/SIM/SIMoutput.h b/src/SIM/SIMoutput.h index 9f2f3878..a2598095 100644 --- a/src/SIM/SIMoutput.h +++ b/src/SIM/SIMoutput.h @@ -297,6 +297,8 @@ public: virtual double getEffectivityIndex(const Vectors&, size_t, size_t) const = 0; //! \brief Prints integrated solution norms to the log stream. virtual void printNorms(const Vectors&, size_t = 36) const = 0; + //! \brief Prints out interface force resultants to log stream. + virtual void printIFforces(const Vector&, RealArray&) {} //! \brief Writes out the additional functions to VTF-file. virtual bool writeAddFuncs(int iStep, int& nBlock, int idBlock, double time);