From ea41ff1f65fd3d86f83a5ddca5b73ed0f816faa9 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Mon, 18 Apr 2016 20:31:41 +0200 Subject: [PATCH 1/7] Fixed: Missing initialization of basis member in struct Property. Added: IntegrandBase::getSolutions(). Changed: Made IntegrandBase::registerVector public (needed by PoroElasticity). --- src/ASM/IntegrandBase.h | 26 +++++++++++++++----------- src/SIM/Property.h | 8 ++++---- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/ASM/IntegrandBase.h b/src/ASM/IntegrandBase.h index 65d89cd38..de52ddb7f 100644 --- a/src/ASM/IntegrandBase.h +++ b/src/ASM/IntegrandBase.h @@ -59,7 +59,7 @@ class IntegrandBase : public Integrand //! \brief Defines the solution mode before the element assembly is started. virtual void setMode(SIM::SolutionMode mode) { m_mode = mode; } - //! \brief Obtain current solution mode + //! \brief Returns current solution mode. SIM::SolutionMode getMode() const { return m_mode; } //! \brief Initializes an integration parameter for the integrand. virtual void setIntegrationPrm(unsigned short int, double) {} @@ -114,7 +114,7 @@ class IntegrandBase : public Integrand const FiniteElement& fe, const Vec3& X0, size_t nPt, LocalIntegral& elmInt); //! \brief Initializes current element for numerical integration (mixed). - //! \param[in] MNPC Nodal point correspondance for the bases + //! \param[in] MNPC Matrix of nodal point correspondance for current element //! \param[in] elem_sizes Size of each basis on the element //! \param[in] basis_sizes Size of each basis on the patch //! \param elmInt Local integral for element @@ -129,7 +129,7 @@ class IntegrandBase : public Integrand virtual bool initElementBou(const std::vector& MNPC, LocalIntegral& elmInt); //! \brief Initializes current element for boundary integration (mixed). - //! \param[in] MNPC Nodal point correspondance for the bases + //! \param[in] MNPC Matrix of nodal point correspondance for current element //! \param[in] elem_sizes Size of each basis on the element //! \param[in] basis_sizes Size of each basis on the patch //! \param elmInt Local integral for element @@ -153,7 +153,7 @@ class IntegrandBase : public Integrand //! \param[in] X Cartesian coordinates of current point //! \param[in] MNPC Nodal point correspondance for the basis function values virtual bool evalSol(Vector& s, const FiniteElement& fe, const Vec3& X, - const std::vector& MNPC) const; + const std::vector& MNPC) const; //! \brief Evaluates the secondary solution at a result point (mixed problem). //! \param[out] s The solution field values at current point @@ -163,7 +163,7 @@ class IntegrandBase : public Integrand //! \param[in] elem_sizes Size of each basis on the element //! \param[in] basis_sizes Size of each basis on the patch virtual bool evalSol(Vector& s, const MxFiniteElement& fe, const Vec3& X, - const std::vector& MNPC, + const std::vector& MNPC, const std::vector& elem_sizes, const std::vector& basis_sizes) const; @@ -172,21 +172,21 @@ class IntegrandBase : public Integrand //! \param[in] asol The analytical solution field (tensor field) //! \param[in] X Cartesian coordinates of current point virtual bool evalSol(Vector& s, - const TensorFunc& asol, const Vec3& X) const; + const TensorFunc& asol, const Vec3& X) const; //! \brief Evaluates the analytical secondary solution at a result point. //! \param[out] s The solution field values at current point //! \param[in] asol The analytical solution field (symmetric tensor field) //! \param[in] X Cartesian coordinates of current point virtual bool evalSol(Vector& s, - const STensorFunc& asol, const Vec3& X) const; + const STensorFunc& asol, const Vec3& X) const; //! \brief Evaluates the analytical secondary solution at a result point. //! \param[out] s The solution field values at current point //! \param[in] asol The analytical solution field (vector field) //! \param[in] X Cartesian coordinates of current point virtual bool evalSol(Vector& s, - const VecFunc& asol, const Vec3& X) const; + const VecFunc& asol, const Vec3& X) const; //! \brief Returns an evaluated principal direction vector field for plotting. virtual bool getPrincipalDir(Matrix&, size_t, size_t) const { return false; } @@ -234,6 +234,8 @@ class IntegrandBase : public Integrand //! \brief Accesses the primary solution vector of current patch. Vector& getSolution(size_t n = 0) { return primsol[n]; } + //! \brief Accesses the primary solution vectors of current patch. + Vectors& getSolutions() { return primsol; } //! \brief Resets the primary solution vectors. void resetSolution(); @@ -255,7 +257,6 @@ class IntegrandBase : public Integrand return LinAlg::GENERAL_MATRIX; } -protected: //! \brief Registers a vector to inject a named field into. //! \param[in] name Name of field //! \param[in] vec Vector to inject field into @@ -297,6 +298,7 @@ class NormBase : public Integrand //! \brief Sets a vector of LocalIntegrals to be used during norm integration. void setLocalIntegrals(LintegralVec* elementNorms) { lints = elementNorms; } + using Integrand::getLocalIntegral; //! \brief Returns a local integral container for the element \a iEl. virtual LocalIntegral* getLocalIntegral(size_t, size_t iEl, bool) const; @@ -398,6 +400,7 @@ class ForceBase : public Integrand //! \brief Initializes the integrand with the number of integration points. virtual void initIntegration(size_t, size_t) {} + using Integrand::getLocalIntegral; //! \brief Returns a local integral container for the element \a iEl. virtual LocalIntegral* getLocalIntegral(size_t, size_t iEl, bool = false) const; @@ -412,7 +415,8 @@ class ForceBase : public Integrand { return false; } //! \brief Dummy implementation (only boundary integration is relevant). - virtual bool initElement(const std::vector&, const std::vector&, + virtual bool initElement(const std::vector&, + const std::vector&, const std::vector&, LocalIntegral&) { return false; } @@ -421,7 +425,7 @@ class ForceBase : public Integrand virtual bool initElementBou(const std::vector& MNPC, LocalIntegral& elmInt); //! \brief Initializes current element for boundary integration (mixed). - virtual bool initElementBou(const std::vector& MNPC1, + virtual bool initElementBou(const std::vector& MNPC, const std::vector& elem_sizes, const std::vector& basis_sizes, LocalIntegral& elmInt); diff --git a/src/SIM/Property.h b/src/SIM/Property.h index 43366ee5e..7254fc45b 100644 --- a/src/SIM/Property.h +++ b/src/SIM/Property.h @@ -42,18 +42,18 @@ struct Property }; Type pcode; //!< Physical property code - int pindx; //!< Physical property index (1-based) + int pindx; //!< Physical property index (0-based) size_t patch; //!< Patch index [1,nPatch] char lindx; //!< Local entity index (1-based) which is assigned the property char ldim; //!< Local entity dimension flag [0,3] char basis; //!< Which basis the property is defined on //! \brief Default constructor. - Property() : pcode(UNDEFINED), pindx(0), patch(0), lindx(0), ldim(0) {} + Property() : pcode(UNDEFINED), pindx(0), patch(0) { lindx=ldim = basis = 0; } //! \brief Constructor creating an initialized property instance. - Property(Type t, int px, size_t p, char ld, char lx = 0) : - pcode(t), pindx(px), patch(p), lindx(lx), ldim(ld) {} + Property(Type t, int px, size_t p, char ld, char lx = 0, char b = 0) : + pcode(t), pindx(px), patch(p), lindx(lx), ldim(ld), basis(b) {} }; typedef std::vector PropertyVec; //!< Vector of properties From 6aac59a0d9b3c5e6eaa11773bb37bdc2a113cd95 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 19 Apr 2016 13:58:47 +0200 Subject: [PATCH 2/7] Added: operator=(const std::vector&), begin(), end() and expandRows(int) for utl::matrix, and a unit test for the latter method. --- src/LinAlg/Test/TestMatrix.C | 57 ++++++++++++++++++++++++------------ src/LinAlg/matrix.h | 54 +++++++++++++++++++++++++++++++++- 2 files changed, 91 insertions(+), 20 deletions(-) diff --git a/src/LinAlg/Test/TestMatrix.C b/src/LinAlg/Test/TestMatrix.C index b4fc38722..5e78e863c 100644 --- a/src/LinAlg/Test/TestMatrix.C +++ b/src/LinAlg/Test/TestMatrix.C @@ -17,25 +17,44 @@ TEST(TestMatrix, AddBlock) { - utl::matrix a(3,3); - int data_a[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9}; - a.fill(data_a, 9); - - utl::matrix b(2,2); - int data_b[9] = {1, 2, 3, 4}; - b.fill(data_b, 4); - - a.addBlock(b, 2, 2, 2, false); - - ASSERT_EQ(a(2,2), 7); - ASSERT_EQ(a(3,2), 10); - ASSERT_EQ(a(2,3), 14); - ASSERT_EQ(a(3,3), 17); + utl::matrix a(3,3),b(2,2); + std::iota(a.begin(),a.end(),1); + std::iota(b.begin(),b.end(),1); + + a.addBlock(b, 2, 2, 2, false); + ASSERT_EQ(a(2,2), 7); + ASSERT_EQ(a(3,2), 10); + ASSERT_EQ(a(2,3), 14); + ASSERT_EQ(a(3,3), 17); + + a.addBlock(b, 1, 1, 1, true); + ASSERT_EQ(a(1,1), 2); + ASSERT_EQ(a(2,1), 5); + ASSERT_EQ(a(1,2), 6); + ASSERT_EQ(a(2,2), 11); +} - a.addBlock(b, 1, 1, 1, true); - ASSERT_EQ(a(1,1), 2); - ASSERT_EQ(a(2,1), 5); - ASSERT_EQ(a(1,2), 6); - ASSERT_EQ(a(2,2), 11); +TEST(TestMatrix, AddRows) +{ + utl::matrix a(3,5); + std::iota(a.begin(),a.end(),1); + std::cout <<"A:"<< a; + + a.expandRows(1); + std::cout <<"B:"<< a; + int fasit = 1; + for (size_t j = 1; j <= 5; j++) + { + for (size_t i = 1; i <= 3; i++, fasit++) + ASSERT_EQ(a(i,j), fasit); + ASSERT_EQ(a(4,j), 0); + } + + a.expandRows(-2); + std::cout <<"C:"<< a; + fasit = 1; + for (size_t j = 1; j <= 5; j++, fasit++) + for (size_t i = 1; i <= 2; i++, fasit++) + ASSERT_EQ(a(i,j), fasit); } diff --git a/src/LinAlg/matrix.h b/src/LinAlg/matrix.h index e2770e67e..307d79c47 100644 --- a/src/LinAlg/matrix.h +++ b/src/LinAlg/matrix.h @@ -292,6 +292,42 @@ namespace utl //! General utility classes and functions. elem.std::template vector::resize(r*c,T(0)); } + //! \brief Increase or decrease the number of rows in the matrix. + matrix& expandRows(int incRows) + { + int newRows = nrow + incRows; + if (newRows < 1 || ncol < 1) + { + // The matrix is empty + nrow = 0; + elem.clear(); + } + else if (incRows < 0) + { + // The matrix size is reduced + T* newMat = this->ptr() + newRows; + for (size_t c = 1; c < ncol; c++, newMat += newRows) + memcpy(newMat,this->ptr(c),newRows*sizeof(T)); + nrow = newRows; + elem.std::template vector::resize(nrow*ncol); + } + else if (incRows > 0) + { + // The matrix size is increased + size_t oldRows = nrow; + T* oldMat = this->ptr(ncol-1); + nrow = newRows; + elem.std::template vector::resize(nrow*ncol,T(0)); + for (size_t c = ncol-1; c > 0; c--, oldMat -= oldRows) + { + memcpy(this->ptr(c),oldMat,oldRows*sizeof(T)); + for (size_t r = nrow-1; r >= oldRows; r--) + elem[r+nrow*(c-1)] = T(0); + } + } + return *this; + } + //! \brief Query number of matrix rows. size_t rows() const { return nrow; } //! \brief Query number of matrix columns. @@ -308,6 +344,21 @@ namespace utl //! General utility classes and functions. //! \brief Type casting to a one-dimensional vector. operator const vector&() const { return elem; } + //! \brief Iterator to the start of the matrix elements. + typename std::vector::iterator begin() { return elem.begin(); } + //! \brief Iterator to the end of the matrix elements. + typename std::vector::iterator end() { return elem.end(); } + + //! \brief Overloaded assignment operator. + matrix& operator=(const std::vector& X) + { + // Do not use vector::operator= because we don't want to alter size + size_t nval = X.size() < elem.size() ? X.size() : elem.size(); + std::copy(X.begin(),X.begin()+nval,elem.begin()); + std::fill(elem.begin()+nval,elem.end(),T(0)); + return *this; + } + //! \brief Index-1 based element access. //! \details Assuming column-wise storage as in Fortran. T& operator()(size_t r, size_t c) @@ -389,7 +440,8 @@ namespace utl //! General utility classes and functions. } //! \brief Add a scalar multiple of another matrix to a block of the matrix. - void addBlock(const matrix& block, const T& s, size_t row, size_t col, bool transposed = false) + void addBlock(const matrix& block, const T& s, size_t row, size_t col, + bool transposed = false) { size_t nr = transposed ? block.cols() : block.rows(); size_t nc = transposed ? block.rows() : block.cols(); From 2a512315b5083a63c66dc1d9b99e9ac0c20c1ed6 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 24 Apr 2016 08:08:16 +0200 Subject: [PATCH 3/7] Changed: The addTo argument of the matrix-vector multiplication method is now a char (instead of bool) to also support subtraction. --- src/LinAlg/matrix.h | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/LinAlg/matrix.h b/src/LinAlg/matrix.h index 307d79c47..fb9b5fde7 100644 --- a/src/LinAlg/matrix.h +++ b/src/LinAlg/matrix.h @@ -620,9 +620,11 @@ namespace utl //! General utility classes and functions. -# \f$ {\bf Y} = {\bf A}^T {\bf X} \f$ -# \f$ {\bf Y} = {\bf Y} + {\bf A} {\bf X} \f$ -# \f$ {\bf Y} = {\bf Y} + {\bf A}^T {\bf X} \f$ + -# \f$ {\bf Y} = {\bf Y} - {\bf A} {\bf X} \f$ + -# \f$ {\bf Y} = {\bf Y} - {\bf A}^T {\bf X} \f$ */ bool multiply(const std::vector& X, std::vector& Y, - bool transA = false, bool addTo = false) const; + bool transA = false, char addTo = 0) const; //! \brief Outer product between two vectors. bool outer_product(const std::vector& X, const std::vector& Y, @@ -1032,14 +1034,14 @@ namespace utl //! General utility classes and functions. template<> inline bool matrix::multiply(const std::vector& X, std::vector& Y, - bool transA, bool addTo) const + bool transA, char addTo) const { if (!this->compatible(X,transA)) return false; if (!addTo) Y.resize(transA ? ncol : nrow); cblas_sgemv(CblasColMajor, transA ? CblasTrans : CblasNoTrans, - nrow, ncol, 1.0f, + nrow, ncol, addTo < 0 ? -1.0f : 1.0f, this->ptr(), nrow, &X.front(), 1, addTo ? 1.0f : 0.0f, &Y.front(), 1); @@ -1050,14 +1052,14 @@ namespace utl //! General utility classes and functions. template<> inline bool matrix::multiply(const std::vector& X, std::vector& Y, - bool transA, bool addTo) const + bool transA, char addTo) const { if (!this->compatible(X,transA)) return false; if (!addTo) Y.resize(transA ? ncol : nrow); cblas_dgemv(CblasColMajor, transA ? CblasTrans : CblasNoTrans, - nrow, ncol, 1.0, + nrow, ncol, addTo < 0 ? -1.0 : 1.0, this->ptr(), nrow, &X.front(), 1, addTo ? 1.0 : 0.0, &Y.front(), 1); @@ -1233,7 +1235,7 @@ namespace utl //! General utility classes and functions. //============================================================================ template inline - vector& vector::operator*=(const T& c) + vector& vector::operator*=(T c) { for (size_t i = 0; i < this->size(); i++) std::vector::operator[](i) *= c; @@ -1319,7 +1321,7 @@ namespace utl //! General utility classes and functions. } template inline - matrix& matrix::multiply(const T& c) + matrix& matrix::multiply(T c) { for (size_t i = 0; i < elem.size(); i++) elem[i] *= c; @@ -1336,9 +1338,9 @@ namespace utl //! General utility classes and functions. template inline bool matrix::multiply(const std::vector& X, std::vector& Y, - bool transA, bool addTo) const + bool transA, char addTo) const { - if (!this->compatible(X,Y,transA)) return false; + if (!this->compatible(X,transA)) return false; if (!addTo) { Y.clear(); @@ -1348,9 +1350,9 @@ namespace utl //! General utility classes and functions. for (size_t i = 0; i < Y.size(); i++) for (size_t j = 0; j < X.size(); j++) if (transA) - Y[i] += THIS(j+1,i+1) * X[j]; + Y[i] += THIS(j+1,i+1) * (addTo > 0 ? X[j] : -X[j]); else - Y[i] += THIS(i+1,j+1) * X[j]; + Y[i] += THIS(i+1,j+1) * (addTo > 0 ? X[j] : -X[j]); return true; } From e6bec44721eb3ada5cc866290118158c3dd14fc1 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 10 May 2016 14:54:56 +0200 Subject: [PATCH 4/7] Added: A new utl::getAttribute method returning a Vec3 value --- src/Utility/Test/TestUtilities.C | 107 +++++++++++++++------- src/Utility/Test/refdata/getattribute.xml | 3 + src/Utility/Utilities.C | 39 ++++++-- src/Utility/Utilities.h | 16 +++- 4 files changed, 125 insertions(+), 40 deletions(-) diff --git a/src/Utility/Test/TestUtilities.C b/src/Utility/Test/TestUtilities.C index 1cceb3120..c427e3c5b 100644 --- a/src/Utility/Test/TestUtilities.C +++ b/src/Utility/Test/TestUtilities.C @@ -11,6 +11,7 @@ //============================================================================== #include "Utilities.h" +#include "Vec3.h" #include "tinyxml.h" #include "gtest/gtest.h" @@ -27,21 +28,20 @@ TEST(TestUtilities, ParseIntegers) ASSERT_EQ(values.size(), 1U); ASSERT_EQ(values2.size(), 5U); ASSERT_EQ(values[0], 1); - for (int i=0;i<5;++i) + for (int i = 0; i < 5; i++) ASSERT_EQ(values2[i], i+1); } TEST(TestUtilities, ParseKnots) { - const char* simple = "0 0 0.5 1.0 1.0"; - char* meh = strdup(simple); - strtok(meh, " "); - std::vector xi; + std::string simple("0 0 0.5 0.9 1.0"); + strtok(const_cast(simple.c_str())," "); + std::vector xi; utl::parseKnots(xi); ASSERT_EQ(xi.size(), 4U); ASSERT_FLOAT_EQ(xi[0], 0.0); ASSERT_FLOAT_EQ(xi[1], 0.5); - ASSERT_FLOAT_EQ(xi[2], 1.0); + ASSERT_FLOAT_EQ(xi[2], 0.9); ASSERT_FLOAT_EQ(xi[3], 1.0); } @@ -51,10 +51,8 @@ TEST(TestUtilities, ReadLine) str << "Blah blah\n" << "# Commented line\n" << "Blah bluh\n"; - char* tmp = utl::readLine(str); - ASSERT_STREQ(tmp, "Blah blah"); - tmp = utl::readLine(str); - ASSERT_STREQ(tmp, "Blah bluh"); + ASSERT_STREQ(utl::readLine(str), "Blah blah"); + ASSERT_STREQ(utl::readLine(str), "Blah bluh"); } TEST(TestUtilities, IgnoreComments) @@ -75,39 +73,86 @@ TEST(TestUtilities, GetAttribute) { TiXmlDocument doc; doc.LoadFile("src/Utility/Test/refdata/getattribute.xml"); + ASSERT_TRUE(doc.RootElement()); - if (!doc.RootElement()) - ASSERT_TRUE(false); - - const std::vector truebool = {"booltrue", "boolone", - "boolon", "boolyes"}; - const std::vector falsebool = {"boolfalse", "boolzero", - "booloff", "boolno"}; + const char* bTrue[4] = { "booltrue" , "boolone" , "boolon" , "boolyes" }; + const char* bFalse[4] = { "boolfalse", "boolzero", "booloff", "boolno" }; - for (size_t i=0;iFirstChildElement(truebool[i]); - bool b=false; - ASSERT_TRUE(elem != NULL); + int i; + for (i = 0; i < 4; i++) { + const TiXmlElement* elem = doc.RootElement()->FirstChildElement(bTrue[i]); + ASSERT_TRUE(elem != nullptr); + bool b = false; ASSERT_TRUE(utl::getAttribute(elem, "bar", b)); ASSERT_TRUE(b); } - for (size_t i=0;iFirstChildElement(falsebool[i]); - bool b=false; + + for (i = 0; i < 4; i++) { + const TiXmlElement* elem = doc.RootElement()->FirstChildElement(bFalse[i]); + ASSERT_TRUE(elem != nullptr); + bool b = false; ASSERT_TRUE(utl::getAttribute(elem, "bar", b)); ASSERT_FALSE(b); } + const TiXmlElement* elem = doc.RootElement()->FirstChildElement("intval"); - int val; - ASSERT_TRUE(elem != NULL); - ASSERT_TRUE(utl::getAttribute(elem, "bar", val)); - ASSERT_EQ(val, 1); - size_t val2; + ASSERT_TRUE(elem != nullptr); + int val1 = 0; + ASSERT_TRUE(utl::getAttribute(elem, "bar", val1)); + ASSERT_EQ(val1, 1); + size_t val2 = 0; ASSERT_TRUE(utl::getAttribute(elem, "bar", val2)); ASSERT_EQ(val2, 1U); - Real val3; + elem = doc.RootElement()->FirstChildElement("realval"); - ASSERT_TRUE(elem != NULL); + ASSERT_TRUE(elem != nullptr); + double val3 = 0.0; ASSERT_TRUE(utl::getAttribute(elem, "bar", val3)); ASSERT_FLOAT_EQ(val3, 2.01); + + elem = doc.RootElement()->FirstChildElement("vecval1"); + ASSERT_TRUE(elem != nullptr); + Vec3 val4; + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 0.0); + ASSERT_FLOAT_EQ(val4.z, 0.0); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4, 2)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 1.2); + ASSERT_FLOAT_EQ(val4.z, 0.0); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4, 3)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 1.2); + ASSERT_FLOAT_EQ(val4.z, 1.2); + + elem = doc.RootElement()->FirstChildElement("vecval2"); + ASSERT_TRUE(elem != nullptr); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 3.4); + ASSERT_FLOAT_EQ(val4.z, 0.0); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4, 2)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 3.4); + ASSERT_FLOAT_EQ(val4.z, 0.0); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4, 3)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 3.4); + ASSERT_FLOAT_EQ(val4.z, 3.4); + + elem = doc.RootElement()->FirstChildElement("vecval3"); + ASSERT_TRUE(elem != nullptr); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 3.4); + ASSERT_FLOAT_EQ(val4.z, 5.6); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4, 2)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 3.4); + ASSERT_FLOAT_EQ(val4.z, 0.0); + ASSERT_TRUE(utl::getAttribute(elem, "bar", val4, 3)); + ASSERT_FLOAT_EQ(val4.x, 1.2); + ASSERT_FLOAT_EQ(val4.y, 3.4); + ASSERT_FLOAT_EQ(val4.z, 5.6); } diff --git a/src/Utility/Test/refdata/getattribute.xml b/src/Utility/Test/refdata/getattribute.xml index c1a63be46..a5b1a0343 100644 --- a/src/Utility/Test/refdata/getattribute.xml +++ b/src/Utility/Test/refdata/getattribute.xml @@ -9,4 +9,7 @@ + + + diff --git a/src/Utility/Utilities.C b/src/Utility/Utilities.C index 36e1ba5bb..48a418979 100644 --- a/src/Utility/Utilities.C +++ b/src/Utility/Utilities.C @@ -12,6 +12,7 @@ //============================================================================== #include "Utilities.h" +#include "Vec3.h" #include "tinyxml.h" #include #include @@ -167,6 +168,33 @@ bool utl::getAttribute (const TiXmlElement* xml, const char* att, Real& val) } +/*! + If fewer than \a ncomp components specified, use the last value for the rest. + If \a ncomp is zero, use the value zero for the missing components, if any. +*/ + +bool utl::getAttribute (const TiXmlElement* xml, const char* att, + Vec3& val, int ncomp) +{ + if (!xml || !xml->Attribute(att)) + return false; + + std::string value(xml->Attribute(att)); + char* cval = const_cast(value.c_str()); + val.x = atof(strtok(cval," ")); + + for (int i = 1; i < 3; i++) + if (i >= ncomp && ncomp > 0) + val[i] = Real(0); + else if ((cval = strtok(nullptr," "))) + val[i] = atof(cval); + else + val[i] = ncomp > 0 ? val[i-1] : Real(0); + + return true; +} + + bool utl::getAttribute (const TiXmlElement* xml, const char* att, std::string& val, bool toLower) { @@ -402,15 +430,14 @@ void utl::merge (std::vector& a1, const std::vector& a2, } -void utl::interleave(const std::vector& v1, const std::vector& v2, - std::vector& out, size_t n1, size_t n2) +void utl::interleave (const std::vector& v1, const std::vector& v2, + std::vector& out, size_t n1, size_t n2) { out.resize(v1.size()+v2.size()); - auto it_v1 = v1.begin(); - auto it_v2 = v2.begin(); - for (auto it_out = out.begin(); it_out != out.end(); - it_out += n1+n2, it_v1 += n1, it_v2 += n2) { + std::vector::iterator it_out = out.begin(); + std::vector::const_iterator it_v1 = v1.begin(), it_v2 = v2.begin(); + for (; it_out != out.end(); it_out += n1+n2, it_v1 += n1, it_v2 += n2) { std::copy(it_v1, it_v1+n1, it_out); std::copy(it_v2, it_v2+n2, it_out+n1); } diff --git a/src/Utility/Utilities.h b/src/Utility/Utilities.h index 3f646fb04..9841ccc1b 100644 --- a/src/Utility/Utilities.h +++ b/src/Utility/Utilities.h @@ -22,6 +22,7 @@ #include #include +class Vec3; class TiXmlElement; class TiXmlNode; @@ -83,6 +84,15 @@ namespace utl //! \return \e true if the attribute \a att is found in \a xml, //! otherwise \e false bool getAttribute(const TiXmlElement* xml, const char* att, Real& val); + //! \brief Extracts a vector attribute value from the specified XML-element. + //! \param[in] xml Pointer to XML-element to extract from + //! \param[in] att The attribute tag + //! \param[out] val The attribute value + //! \param[in] ncomp Maximum number of components to read + //! \return \e true if the attribute \a att is found in \a xml, + //! otherwise \e false + bool getAttribute(const TiXmlElement* xml, const char* att, Vec3& val, + int ncomp = 0); //! \brief Extracts a string attribute value from the specified XML-element. //! \param[in] xml Pointer to XML-element to extract from //! \param[in] att The attribute tag @@ -152,12 +162,12 @@ namespace utl const std::vector& in, std::vector& out, size_t offset_in = 0, int shift_idx = 0); - //! \brief Interleave vectors + //! \brief Merges two arrays of nodal values into one array by interleaving. //! \param[in] v1 The first array //! \param[in] v2 The second array //! \param[out] out The output array - //! \param[in] n1 Number of entries per node in first vector - //! \param[in] n2 Number of entries per node in second vector + //! \param[in] n1 Number of entries per node in first array + //! \param[in] n2 Number of entries per node in second array void interleave(const std::vector& v1, const std::vector& v2, std::vector& out, size_t n1 = 1, size_t n2 = 1); From c5860d0f3a910d58ea4f5074014c88fab364e33d Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Mon, 16 May 2016 12:28:05 +0200 Subject: [PATCH 5/7] Added: Unit test for the Newmark time integrator --- src/SIM/SIMbase.C | 6 +- src/SIM/SIMbase.h | 4 +- src/SIM/Test/TestInitialConditions.C | 19 +-- src/SIM/Test/TestNewmark.C | 201 +++++++++++++++++++++++++++ 4 files changed, 212 insertions(+), 18 deletions(-) create mode 100644 src/SIM/Test/TestNewmark.C diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 148d28992..e6ba82c9c 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -1537,10 +1537,8 @@ bool SIMbase::updateDirichlet (double time, const Vector* prevSol) #endif return false; - if (mySam) - return static_cast(mySam)->updateConstraintEqs(myModel,prevSol); - else - return true; + SAMpatch* pSam = dynamic_cast(mySam); + return pSam ? pSam->updateConstraintEqs(myModel,prevSol) : true; } diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index 1e227d25d..07abb492f 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -289,8 +289,8 @@ class SIMbase : public SIMinput, public SIMdependency //! \param[in] pSol Previous primary solution vectors in DOF-order //! \param[in] newLHSmatrix If \e false, only integrate the RHS vector //! \param[in] poorConvg If \e true, the nonlinear driver is converging poorly - bool assembleSystem(const TimeDomain& time, const Vectors& pSol, - bool newLHSmatrix = true, bool poorConvg = false); + virtual bool assembleSystem(const TimeDomain& time, const Vectors& pSol, + bool newLHSmatrix = true, bool poorConvg = false); //! \brief Administers assembly of the linear equation system. //! \param[in] pSol Primary solution vectors in DOF-order diff --git a/src/SIM/Test/TestInitialConditions.C b/src/SIM/Test/TestInitialConditions.C index 7564e6701..fb9823286 100644 --- a/src/SIM/Test/TestInitialConditions.C +++ b/src/SIM/Test/TestInitialConditions.C @@ -16,22 +16,17 @@ TEST(TestInitialConditions, Parse) { - SIM2D sim({4}, false); + SIM2D sim(4); EXPECT_TRUE(sim.read("src/SIM/Test/input.xinp")); // Recognize both comp and component attributes and correct priority // Boundary conditions - RealFunc* func; - const Vec3 X; - for (int i = 1; i < 5; i++) { - func = sim.getSclFunc(i); - ASSERT_FLOAT_EQ((float) i, (*func)(X)); - } + for (int i = 1; i < 5; i++) + ASSERT_FLOAT_EQ((float)i,(*sim.getSclFunc(i))(Vec3())); // Initial conditions - const SIMdependency::InitialCondMap& ic = sim.getICs(); - SIMdependency::InitialCondMap::const_iterator foo = ic.find("nofile"); - std::vector bar = foo->second; - for (std::vector::iterator info = bar.begin(); + ASSERT_TRUE(sim.getICs().begin() != sim.getICs().end()); + const std::vector& bar = sim.getICs().begin()->second; + for (std::vector::const_iterator info = bar.begin(); info != bar.end(); info++) switch (info->component) { case 1: @@ -47,6 +42,6 @@ TEST(TestInitialConditions, Parse) ASSERT_EQ(info->function, "4"); break; default: - ASSERT_TRUE( false ); + ASSERT_TRUE(false); } } diff --git a/src/SIM/Test/TestNewmark.C b/src/SIM/Test/TestNewmark.C new file mode 100644 index 000000000..6635fad79 --- /dev/null +++ b/src/SIM/Test/TestNewmark.C @@ -0,0 +1,201 @@ +// $Id$ +//============================================================================== +//! +//! \file TestNewmark.C +//! +//! \date May 14 2016 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Tests the Newmark time integrator. +//! +//============================================================================== + +#include "SAM.h" +#include "IntegrandBase.h" +#include "SIMoutput.h" + +#include "GenAlphaSIM.h" +#include "NewmarkMats.h" +#include "AlgEqSystem.h" +#include "TimeStep.h" + +#include "gtest/gtest.h" + + +// SAM class representing a single-DOF system. +class SAM1DOF : public SAM +{ +public: + SAM1DOF() + { + nmmnpc = nel = nnod = ndof = neq = 1; + mmnpc = new int[1]; mmnpc[0] = 1; + mpmnpc = new int[2]; std::iota(mpmnpc,mpmnpc+2,1); + madof = new int[2]; std::iota(madof,madof+2,1); + msc = new int[1]; msc[0] = 1; + EXPECT_TRUE(this->initSystemEquations()); + } + virtual ~SAM1DOF() {} +}; + + +// Dummy integrand class, for time integration scheme testing. +class Problem : public IntegrandBase +{ +public: + Problem() : IntegrandBase(1) {} + virtual ~Problem() {} + virtual void setIntegrationPrm(unsigned short int i, double p) { prm[i] = p; } + virtual double getIntegrationPrm(unsigned short int i) const + { return m_mode == SIM::DYNAMIC ? prm[i] : 0.0; } + const double* getIntPrm() { return prm; } +private: + double prm[5]; +}; + + +// Simulator wrapper implementing dummy versions of the pure virtual methods. +class SIMwrapper : public SIMoutput +{ +public: + SIMwrapper(IntegrandBase* p) : SIMoutput(p) {} + virtual ~SIMwrapper() {} + virtual unsigned short int getNoParamDim() const { return 0; } + virtual bool addConstraint(int,int,int,int,int,int&,char) + { return false; } + virtual ASMbase* readPatch(std::istream&,int,const CharVec&) const + { return nullptr; } + virtual bool readPatches(std::istream&,SIMdependency::PatchVec&,const char*) + { return false; } + virtual ASMbase* createDefaultGeometry(const TiXmlElement*) const + { return nullptr; } +}; + + +// Simulator class for a single-DOF oscillator. +class SIM1DOF : public SIMwrapper +{ +public: + SIM1DOF(IntegrandBase* p) : SIMwrapper(p) { mySam = new SAM1DOF(); } + virtual ~SIM1DOF() {} + virtual bool assembleSystem(const TimeDomain& time, + const Vectors& prevSol, + bool newLHSmatrix, bool) + { + const double M = 10.0; // Mass of the oscillator + const double K = 1000.0; // Stiffness of the oscillator + const double F = 1.0; // External load (constant) + + myEqSys->initialize(newLHSmatrix); + + bool ok; + if (myProblem->getMode() == SIM::MASS_ONLY) { + ElmMats elm; + elm.resize(1,1); elm.redim(1); + elm.A[0].fill(M); // Mass Matrix + ok = myEqSys->assemble(&elm,1); + } + else { + const double* intPrm = static_cast(myProblem)->getIntPrm(); + NewmarkMats elm(intPrm[0],intPrm[1],intPrm[2],intPrm[3]); + elm.resize(3,1); elm.redim(1); elm.vec.resize(3); + elm.setStepSize(time.dt,0); + elm.A[1].fill(M); // Mass matrix + elm.A[2].fill(K); // Stiffness matrix + elm.b[0] = -K*prevSol.front(); // Elastic forces + for (int i = 0; i < 3; i++) elm.vec[i] = prevSol[i]; + ok = myEqSys->assemble(&elm,1); + } + + // Add in the external load + ok &= mySam->assembleSystem(*myEqSys->getVector(),&F,1); + + return ok && myEqSys->finalize(newLHSmatrix); + } +}; + + +// Newmark time integrator with numerical damping (alpha_H = -0.1). +class Newmark : public NewmarkSIM +{ +public: + Newmark(SIMbase& sim) : NewmarkSIM(sim) + { + beta = 0.3025; gamma = 0.6; + this->initPrm(); + this->initSol(3); + } + virtual ~Newmark() {} +}; + + +// Generalized-alpha time integrator with numerical damping (alpha_H = -0.1). +class GenAlpha : public GenAlphaSIM +{ +public: + GenAlpha(SIMbase& sim) : GenAlphaSIM(sim) + { + this->initPrm(); + this->initSol(3); + } + virtual ~GenAlpha() {} +}; + + +void runSingleDof (SIMbase& model, NewmarkSIM& solver, double rtol = 0.5e-11) +{ + TimeStep tp; + tp.time.dt = 0.01; + tp.stopTime = 0.65; + + ASSERT_TRUE(model.initSystem(0)); + ASSERT_TRUE(solver.initAcc()); + EXPECT_FLOAT_EQ(solver.getAcceleration().front(),0.1); + + // at t=0.1 at t=0.25 at t=0.5 + double u[3] = { 0.000457484252515, 0.00178698471292, 0.000732016593476 }; + double v[3] = { 0.008368445734720, 0.00592764975245,-0.009365075630580 }; + double a[3] = { 0.054251574748500,-0.07869847129160, 0.026798340652400 }; + + while (solver.advanceStep(tp)) + { + ASSERT_TRUE(solver.solveStep(tp) == SIM::CONVERGED); + double dis = solver.getSolution().front(); + double vel = solver.getVelocity().front(); + double acc = solver.getAcceleration().front(); + + // Check the response at three randomly selected time steps + if (tp.step == 10) { + EXPECT_NEAR(dis,u[0], (dis+u[0])*rtol); + EXPECT_NEAR(vel,v[0], (vel+v[0])*rtol); + EXPECT_NEAR(acc,a[0], (acc+a[0])*rtol); + } + else if (tp.step == 25) { + EXPECT_NEAR(dis,u[1], (dis+u[1])*rtol); + EXPECT_NEAR(vel,v[1], (vel+v[1])*rtol); + EXPECT_NEAR(acc,a[1],-(acc+a[1])*rtol); + } + else if (tp.step == 50) { + EXPECT_NEAR(dis,u[2], (dis+u[2])*rtol); + EXPECT_NEAR(vel,v[2],-(vel+v[2])*rtol); + EXPECT_NEAR(acc,a[2], (acc+a[2])*rtol); + } + } +} + + +TEST(TestNewmark, SingleDOF) +{ + SIM1DOF simulator(new Problem()); + Newmark integrator(simulator); + runSingleDof(simulator,integrator); +} + + +TEST(TestGenAlpha, SingleDOF) +{ + SIM1DOF simulator(new Problem()); + GenAlpha integrator(simulator); + runSingleDof(simulator,integrator,0.02); +} From 9a8795f4c67b82cf056bfe514f272c14599a2423 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 20 May 2016 22:44:59 +0200 Subject: [PATCH 6/7] Added: Option to use incremental displacements as primary unknowns instead of acceleration. This must be used when inhomogeneous Dirichlet conditions are present, combined with the constant displacment predictor. --- src/ASM/NewmarkMats.C | 8 +- src/ASM/NewmarkMats.h | 1 + src/SIM/GenAlphaSIM.C | 10 ++- src/SIM/NewmarkSIM.C | 16 ++-- src/SIM/NewmarkSIM.h | 1 + src/SIM/Test/TestNewmark.C | 163 +++++++++++++++++++++++++++++++++++-- 6 files changed, 180 insertions(+), 19 deletions(-) diff --git a/src/ASM/NewmarkMats.C b/src/ASM/NewmarkMats.C index dd211cf36..da1264751 100644 --- a/src/ASM/NewmarkMats.C +++ b/src/ASM/NewmarkMats.C @@ -23,7 +23,7 @@ NewmarkMats::NewmarkMats (double a1, double a2, double b, double c, if (generalizedAlpha) { - alpha_m = b; + alpha_m = fabs(b); alpha_f = c; double alpha = alpha_f - alpha_m; beta = 0.25*(1.0-alpha)*(1.0-alpha); @@ -32,9 +32,11 @@ NewmarkMats::NewmarkMats (double a1, double a2, double b, double c, else { alpha_m = alpha_f = 1.0; - beta = b; + beta = fabs(b); gamma = c; } + + slvDisp = b < 0.0; // Displacement increments are used as primary unknowns } @@ -45,7 +47,7 @@ const Matrix& NewmarkMats::getNewtonMatrix () const N = A[1]; N.multiply(alpha_m + alpha_f*alpha1*gamma*h); N.add(A[2],alpha_f*(alpha2*gamma + beta*h)*h); - + if (slvDisp) N.multiply(1.0/(beta*h*h)); #if SP_DEBUG > 2 std::cout <<"\nElement mass matrix"<< A[1]; std::cout <<"Element stiffness matrix"<< A[2]; diff --git a/src/ASM/NewmarkMats.h b/src/ASM/NewmarkMats.h index 620cb70aa..c43e3a9dc 100644 --- a/src/ASM/NewmarkMats.h +++ b/src/ASM/NewmarkMats.h @@ -60,6 +60,7 @@ class NewmarkMats : public ElmMats private: double alpha_m; //!< Generalized-alpha parameter double alpha_f; //!< Generalized-alpha parameter + bool slvDisp; //!< If \e true, solve for displacement increments }; #endif diff --git a/src/SIM/GenAlphaSIM.C b/src/SIM/GenAlphaSIM.C index c2a32d246..6e7992d4a 100644 --- a/src/SIM/GenAlphaSIM.C +++ b/src/SIM/GenAlphaSIM.C @@ -61,7 +61,7 @@ void GenAlphaSIM::initPrm () { model.setIntegrationPrm(0,alpha1); model.setIntegrationPrm(1,fabs(alpha2)); - model.setIntegrationPrm(2,alpha_m); + model.setIntegrationPrm(2,solveDisp ? -alpha_m : alpha_m); model.setIntegrationPrm(3,alpha_f); model.setIntegrationPrm(4,2.0); } @@ -142,9 +142,11 @@ bool GenAlphaSIM::correctStep (TimeStep& param, bool converged) size_t ipV = solution.size() - 2; // Update current displacement, velocity and acceleration solutions - tempSol[0].add(linsol,beta*param.time.dt*param.time.dt); - tempSol[1].add(linsol,gamma*param.time.dt); - tempSol[2].add(linsol,1.0); + double bdt = beta*param.time.dt; + double bdt2 = bdt*param.time.dt; + tempSol[0].add(linsol, solveDisp ? 1.0 : bdt2); + tempSol[1].add(linsol, solveDisp ? gamma/bdt : gamma*param.time.dt); + tempSol[2].add(linsol, solveDisp ? 1.0/bdt2 : 1.0); if (converged) { diff --git a/src/SIM/NewmarkSIM.C b/src/SIM/NewmarkSIM.C index cff04d24d..ffb843971 100644 --- a/src/SIM/NewmarkSIM.C +++ b/src/SIM/NewmarkSIM.C @@ -29,6 +29,7 @@ NewmarkSIM::NewmarkSIM (SIMbase& sim) : MultiStepSIM(sim) alpha1 = 0.0; alpha2 = 0.0; + solveDisp = false; // default use acceleration as primary variables predictor = 'a'; // default predictor (zero acceleration) cNorm = 1; // default convergence check, force residual @@ -92,6 +93,8 @@ bool NewmarkSIM::parse (const TiXmlElement* elem) } else if ((value = utl::getValue(child,"rotation"))) rotUpd = tolower(value[0]); + else if (!strncasecmp(child->Value(),"solve_dis",9)) + solveDisp = true; // no need for value here return true; } @@ -108,7 +111,8 @@ void NewmarkSIM::printProblem () const case 'v': IFEM::cout <<"\n- using constant velocity predictor"; break; case 'a': IFEM::cout <<"\n- using zero acceleration predictor"; break; } - + if (solveDisp) + IFEM::cout <<"\n- using displacement increments as primary unknowns"; if (alpha1 > 0.0) IFEM::cout <<"\nMass-proportional damping (alpha1): "<< alpha1; if (alpha2 != 0.0) @@ -122,7 +126,7 @@ void NewmarkSIM::initPrm () { model.setIntegrationPrm(0,alpha1); model.setIntegrationPrm(1,fabs(alpha2)); - model.setIntegrationPrm(2,beta); + model.setIntegrationPrm(2,solveDisp ? -beta : beta); model.setIntegrationPrm(3,gamma); } @@ -272,13 +276,13 @@ bool NewmarkSIM::correctStep (TimeStep& param, bool) size_t iV = solution.size() - 2; // Corrected displacement - solution[iD].add(linsol,beta*dt*dt); + solution[iD].add(linsol, solveDisp ? 1.0 : beta*dt*dt); // Corrected velocity - solution[iV].add(linsol,gamma*dt); + solution[iV].add(linsol, solveDisp ? gamma/(beta*dt) : gamma*dt); // Corrected acceleration - solution[iA].add(linsol,1.0); + solution[iA].add(linsol, solveDisp ? 1.0/(beta*dt*dt) : 1.0); #if SP_DEBUG > 1 std::cout <<"Corrected displacement:"<< solution[iD]; @@ -289,7 +293,7 @@ bool NewmarkSIM::correctStep (TimeStep& param, bool) if (rotUpd == 't') model.updateRotations(solution[iD]); else if (rotUpd) - model.updateRotations(linsol,beta*dt*dt); + model.updateRotations(linsol, solveDisp ? 1.0 : beta*dt*dt); return model.updateConfiguration(solution[iD]); } diff --git a/src/SIM/NewmarkSIM.h b/src/SIM/NewmarkSIM.h index 98793b652..bbf8ee122 100644 --- a/src/SIM/NewmarkSIM.h +++ b/src/SIM/NewmarkSIM.h @@ -97,6 +97,7 @@ class NewmarkSIM : public MultiStepSIM double gamma; //!< Newmark time integration parameter // Solution algorithm parameters + bool solveDisp; //!< If \e true, use incremental displacements as unknowns char predictor; //!< Predictor type flag int maxit; //!< Maximum number of iterations in a time step int saveIts; //!< Time step for which iteration result should be saved diff --git a/src/SIM/Test/TestNewmark.C b/src/SIM/Test/TestNewmark.C index 6635fad79..afa593db1 100644 --- a/src/SIM/Test/TestNewmark.C +++ b/src/SIM/Test/TestNewmark.C @@ -40,6 +40,38 @@ public: }; +// SAM class representing a 2-DOF system where the 2nd DOF is prescribed. +class SAM2DOF : public SAM +{ +public: + SAM2DOF() + { + nmmnpc = nnod = ndof = 2; + nmmceq = nceq = nel = neq = 1; + mmnpc = new int[2]; std::iota(mmnpc,mmnpc+2,1); + mpmnpc = new int[2]; mpmnpc[0] = 1; mpmnpc[1] = 3; + madof = new int[3]; std::iota(madof,madof+3,1); + msc = new int[2]; msc[0] = 1; msc[1] = 0; + mmceq = new int[1]; mmceq[0] = 2; + mpmceq = new int[2]; std::iota(mpmceq,mpmceq+2,1); + ttcc = new double[1]; ttcc[0] = 0.0; + EXPECT_TRUE(this->initSystemEquations()); + } + virtual ~SAM2DOF() {} + + bool updateContrainEqs(double time, const Vector* prevSol) + { + if (!prevSol) + ttcc[0] = 0.0; + else if (prevSol->size() >= 2) + ttcc[0] = sin(10.0*time) - (*prevSol)(2); + else + ttcc[0] = sin(10.0*time); + return true; + } +}; + + // Dummy integrand class, for time integration scheme testing. class Problem : public IntegrandBase { @@ -116,13 +148,61 @@ public: }; +// Simulator class for a two-DOF oscillator with prescribed motion. +class SIM2DOF : public SIMwrapper +{ +public: + SIM2DOF(IntegrandBase* p) : SIMwrapper(p) { mySam = new SAM2DOF(); } + virtual ~SIM2DOF() {} + virtual bool assembleSystem(const TimeDomain& time, + const Vectors& prevSol, + bool newLHSmatrix, bool) + { + const double M = 1.0; // Mass of the oscillator + const double K = 1000.0; // Stiffness of the oscillator + + myEqSys->initialize(newLHSmatrix); + + bool ok; + if (myProblem->getMode() == SIM::MASS_ONLY) { + ElmMats elm; + elm.resize(1,1); elm.redim(2); + elm.A[0](1,1) = elm.A[0](2,2) = M; // Mass matrix + ok = myEqSys->assemble(&elm,1); + } + else { + const double* intPrm = static_cast(myProblem)->getIntPrm(); + NewmarkMats elm(intPrm[0],intPrm[1],intPrm[2],intPrm[3]); + elm.resize(3,1); elm.redim(2); elm.vec.resize(3); + elm.setStepSize(time.dt,0); + elm.A[1](1,1) = elm.A[1](2,2) = M; // Mass Matrix + elm.A[2](1,1) = elm.A[2](2,2) = K; // Stiffness matrix + elm.A[2](2,1) = elm.A[2](1,2) = -K; + elm.b[0] = elm.A[2]*prevSol.front(); // Elastic forces + elm.b[0] *= -1.0; + for (int i = 0; i < 3; i++) elm.vec[i] = prevSol[i]; + ok = myEqSys->assemble(&elm,1); + } + + return ok && myEqSys->finalize(newLHSmatrix); + } + + virtual bool updateDirichlet (double time, const Vector* prevSol) + { + return static_cast(mySam)->updateContrainEqs(time,prevSol); + } +}; + + // Newmark time integrator with numerical damping (alpha_H = -0.1). class Newmark : public NewmarkSIM { public: - Newmark(SIMbase& sim) : NewmarkSIM(sim) + Newmark(SIMbase& sim, bool useDispl) : NewmarkSIM(sim) { beta = 0.3025; gamma = 0.6; + solveDisp = useDispl; + predictor = useDispl ? 'd' : 'a'; this->initPrm(); this->initSol(3); } @@ -134,8 +214,10 @@ public: class GenAlpha : public GenAlphaSIM { public: - GenAlpha(SIMbase& sim) : GenAlphaSIM(sim) + GenAlpha(SIMbase& sim, bool useDispl) : GenAlphaSIM(sim) { + solveDisp = useDispl; + predictor = useDispl ? 'd' : 'a'; this->initPrm(); this->initSol(3); } @@ -185,17 +267,86 @@ void runSingleDof (SIMbase& model, NewmarkSIM& solver, double rtol = 0.5e-11) } -TEST(TestNewmark, SingleDOF) +void runPrescribed (SIMbase& model, NewmarkSIM& solver, double rtol = 0.5e-11) +{ + TimeStep tp; + tp.time.dt = 0.01; + tp.stopTime = 0.65; + + ASSERT_TRUE(model.initSystem(0)); + + // at t=0.1 at t=0.25 at t=0.5 + double u[3] = {0.9312639435267, 0.3547915343361, -1.075289543029 }; + double v[3] = { 16.61378563136, -8.930847244282, 11.85822968077 }; + double a[3] = {-89.79295871881, 243.6806097678, 116.3652683660 }; + + while (solver.advanceStep(tp)) + { + ASSERT_TRUE(solver.solveStep(tp) == SIM::CONVERGED); + double dis = solver.getSolution().front(); + double vel = solver.getVelocity().front(); + double acc = solver.getAcceleration().front(); + + // Check the response at three randomly selected time steps + if (tp.step == 10) { + EXPECT_NEAR(dis,u[0], (dis+u[0])*rtol); + EXPECT_NEAR(vel,v[0], (vel+v[0])*rtol); + EXPECT_NEAR(acc,a[0],-(acc+a[0])*rtol); + } + else if (tp.step == 25) { + EXPECT_NEAR(dis,u[1], (dis+u[1])*rtol); + EXPECT_NEAR(vel,v[1],-(vel+v[1])*rtol); + EXPECT_NEAR(acc,a[1], (acc+a[1])*rtol); + } + else if (tp.step == 50) { + EXPECT_NEAR(dis,u[2],-(dis+u[2])*rtol); + EXPECT_NEAR(vel,v[2], (vel+v[2])*rtol); + EXPECT_NEAR(acc,a[2], (acc+a[2])*rtol); + } + } +} + + +TEST(TestNewmark, SingleDOFa) { SIM1DOF simulator(new Problem()); - Newmark integrator(simulator); + Newmark integrator(simulator,false); runSingleDof(simulator,integrator); } +TEST(TestGenAlpha, SingleDOFa) +{ + SIM1DOF simulator(new Problem()); + GenAlpha integrator(simulator,false); + runSingleDof(simulator,integrator,0.02); +} -TEST(TestGenAlpha, SingleDOF) +TEST(TestNewmark, SingleDOFu) { SIM1DOF simulator(new Problem()); - GenAlpha integrator(simulator); + Newmark integrator(simulator,true); + runSingleDof(simulator,integrator); +} + +TEST(TestNewmark, Prescribed) +{ + SIM2DOF simulator(new Problem()); + Newmark integrator(simulator,true); + runPrescribed(simulator,integrator); +} + +/* does not work, yet +TEST(TestGenAlpha, SingleDOFu) +{ + SIM1DOF simulator(new Problem()); + GenAlpha integrator(simulator,true); runSingleDof(simulator,integrator,0.02); } + +TEST(TestGenAlpha, Prescribed) +{ + SIM2DOF simulator(new Problem()); + GenAlpha integrator(simulator,true); + runPrescribed(simulator,integrator,0.02); +} +*/ From c450f8fe067c05c739f1cd0bf888f5d5e1ed0fcb Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 24 May 2016 10:53:17 +0200 Subject: [PATCH 7/7] Fixed: Crash when trying to access empty nodeInd array when using Lagrange --- src/ASM/ASMs2D.C | 8 ++++++-- src/ASM/ASMs3D.C | 9 ++++++--- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 9691ac4d5..ee9e7c06d 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -578,8 +578,12 @@ bool ASMs2D::assignNodeNumbers (BlockNodes& nodes, int basis) #if SP_DEBUG > 2 if (basis > 0) std::cout <<"\nBasis "<< basis <<":"; for (int i = inod-n1*n2; i < inod; i++) - std::cout <<"\nNode "<< i+1 <<"\t: "<< nodeInd[i].I <<" "<< nodeInd[i].J - <<"\tglobal no. "<< MLGN[i]; + { + std::cout <<"\nNode "<< i+1 <<"\t: "; + if (!nodeInd.empty()) + std::cout << nodeInd[i].I <<" "<< nodeInd[i].J; + std::cout <<"\tglobal no. "<< MLGN[i]; + } std::cout << std::endl; #endif return true; diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index c7368275d..a57401ce4 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -604,9 +604,12 @@ bool ASMs3D::assignNodeNumbers (BlockNodes& nodes, int basis) #if SP_DEBUG > 2 if (basis > 0) std::cout <<"\nBasis "<< basis <<":"; for (int i = inod-n1*n2*n3; i < inod; i++) - std::cout <<"\nNode "<< i+1 <<"\t: " - << nodeInd[i].I <<" "<< nodeInd[i].J <<" "<< nodeInd[i].K - <<"\tglobal no. "<< MLGN[i]; + { + std::cout <<"\nNode "<< i+1 <<"\t: "; + if (!nodeInd.empty()) + std::cout << nodeInd[i].I <<" "<< nodeInd[i].J <<" "<< nodeInd[i].K; + std::cout <<"\tglobal no. "<< MLGN[i]; + } std::cout << std::endl; #endif return true;