From 5438cf1b0535c791aa2eb99736d09134b14ca31b Mon Sep 17 00:00:00 2001 From: timovanopstal Date: Wed, 6 Apr 2016 11:22:23 +0200 Subject: [PATCH 1/6] Fixed: Apply homogeneous Dirichlet on initial conditions --- src/LinAlg/SAM.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/LinAlg/SAM.C b/src/LinAlg/SAM.C index 4a752ff56..7b175ff1a 100644 --- a/src/LinAlg/SAM.C +++ b/src/LinAlg/SAM.C @@ -810,6 +810,8 @@ bool SAM::applyDirichlet (Vector& dofVec) const int ip = mpmceq[iceq-1]; dofVec[idof] = ttcc[ip-1]; } + else if (iceq == 0) + dofVec[idof] = 0.0; } return true; From 7b725ad8a98eb40cff3acbc0098efb2450bc4ecf Mon Sep 17 00:00:00 2001 From: timovanopstal Date: Wed, 6 Apr 2016 13:15:48 +0200 Subject: [PATCH 2/6] Added: Support both comp and component attributes in XML. Added: Unit test for initial conditions. --- cmake/Scripts/IFEMTesting.cmake | 2 +- src/SIM/SIMbase.C | 25 +++++++------ src/SIM/SIMbase.h | 3 ++ src/SIM/Test/TestInitialConditions.C | 52 ++++++++++++++++++++++++++++ src/SIM/Test/input.xinp | 34 ++++++++++++++++++ 5 files changed, 104 insertions(+), 12 deletions(-) create mode 100644 src/SIM/Test/TestInitialConditions.C create mode 100644 src/SIM/Test/input.xinp diff --git a/cmake/Scripts/IFEMTesting.cmake b/cmake/Scripts/IFEMTesting.cmake index 4d9373485..be4b3a893 100644 --- a/cmake/Scripts/IFEMTesting.cmake +++ b/cmake/Scripts/IFEMTesting.cmake @@ -41,7 +41,7 @@ macro(IFEM_add_test_app path workdir name) endmacro() macro(IFEM_add_unittests IFEM_PATH) - IFEM_add_test_app("${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C" + IFEM_add_test_app("${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C;${IFEM_PATH}/src/SIM/Test/*.C" ${IFEM_PATH} IFEM ${IFEM_LIBRARIES} ${IFEM_DEPLIBS}) diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 7b5aa9665..db63f05b8 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -343,8 +343,10 @@ bool SIMbase::parseBCTag (const TiXmlElement* elem) else if (!strcasecmp(elem->Value(),"dirichlet") && !ignoreDirichlet) { int comp = 0; - int basis=1; + int basis = 1; std::string set, type, axes; + // long and short form supported, short prioritized + utl::getAttribute(elem,"component",comp); utl::getAttribute(elem,"comp",comp); utl::getAttribute(elem,"set",set); utl::getAttribute(elem,"type",type,true); @@ -441,38 +443,39 @@ bool SIMbase::parse (const TiXmlElement* elem) { ICInfo info(field); std::string type("file"); - utl::getAttribute(elem,"type", type); + utl::getAttribute(elem,"type",type); if (type == "file") { utl::getAttribute(elem,"file",file); utl::getAttribute(elem,"file_field",info.file_field); utl::getAttribute(elem,"file_level",info.file_level); utl::getAttribute(elem,"geo_level",info.geo_level); - } else { // a function + } + else { // function + // long and short from supported, short prioritized utl::getAttribute(elem,"component",info.component); + utl::getAttribute(elem,"comp",info.component); info.file_field = type; info.function = utl::getValue(elem,"initialcondition"); file = "nofile"; } utl::getAttribute(elem,"level",info.sim_level); - int basis=1; + int basis = 1; utl::getAttribute(elem,"basis",basis); info.basis = basis; IFEM::cout <<"\tInitial condition"; if (info.component > -1) - IFEM::cout << " function: " << info.function; + IFEM::cout <<" function: "<< info.function; else - IFEM::cout << " file: "<< file; + IFEM::cout <<" file: "<< file; IFEM::cout <<"\n\tField name: \""<< info.sim_field; - if (info.component == -1) - IFEM::cout <<"\" (on file \""<< info.file_field <<"\")\n"; + IFEM::cout <<"\" (on file \""<< info.file_field <<"\")"; else - IFEM::cout << " (component " << info.component << " basis " << basis << ")\n"; - - IFEM::cout << "\tTime level: "<< info.sim_level; + IFEM::cout <<" (component "<< info.component <<" basis "<< basis <<")"; + IFEM::cout <<"\n\tTime level: "<< info.sim_level; if (info.component == -1) IFEM::cout <<" (on file "<< info.file_level <<" with basis "<< info.geo_level <<")"; diff --git a/src/SIM/SIMbase.h b/src/SIM/SIMbase.h index f4690e81e..24e1e7e81 100644 --- a/src/SIM/SIMbase.h +++ b/src/SIM/SIMbase.h @@ -671,6 +671,9 @@ class SIMbase : public SIMinput, public SIMdependency static bool ignoreDirichlet; //!< Set to \e true for free vibration analysis static bool preserveNOrder; //!< Set to \e true to preserve node ordering + //! For testing purposes + RealFunc* getSclFunc(int i) { return (myScalars[i]); } + protected: //! \brief Scalar field container typedef std::map SclFuncMap; diff --git a/src/SIM/Test/TestInitialConditions.C b/src/SIM/Test/TestInitialConditions.C new file mode 100644 index 000000000..7564e6701 --- /dev/null +++ b/src/SIM/Test/TestInitialConditions.C @@ -0,0 +1,52 @@ +//============================================================================== +//! +//! \file TestInitialConditions.C +//! +//! \date Apr 6 2016 +//! +//! \author Timo van Opstal / NTNU +//! +//! \brief Tests initial condition handling. +//! +//============================================================================== + +#include "SIM2D.h" + +#include "gtest/gtest.h" + +TEST(TestInitialConditions, Parse) +{ + SIM2D sim({4}, false); + 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)); + } + // 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(); + info != bar.end(); info++) + switch (info->component) { + case 1: + ASSERT_EQ(info->function, "1"); + break; + case 2: + ASSERT_EQ(info->function, "2"); + break; + case 3: + ASSERT_EQ(info->function, "3"); + break; + case 4: + ASSERT_EQ(info->function, "4"); + break; + default: + ASSERT_TRUE( false ); + } +} diff --git a/src/SIM/Test/input.xinp b/src/SIM/Test/input.xinp new file mode 100644 index 000000000..e8537915a --- /dev/null +++ b/src/SIM/Test/input.xinp @@ -0,0 +1,34 @@ + + + + + + + + 1 + + + 2 + + + 3 + + + 4 + + + + + + 1 + 2 + 3 + 4 + + + 1 + 2 + 3 + 4 + + From 779ff4ee4a2485587326274722f2fadcf9993967 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 14 Apr 2016 09:32:58 +0200 Subject: [PATCH 3/6] Added: enum value MixedType::NONE. Fixed: Doxygen comments. --- src/ASM/ASMmxBase.h | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/ASM/ASMmxBase.h b/src/ASM/ASMmxBase.h index f1b5c5ad0..3717092b5 100644 --- a/src/ASM/ASMmxBase.h +++ b/src/ASM/ASMmxBase.h @@ -31,8 +31,7 @@ class ASMmxBase { protected: //! \brief The constructor sets the number of field variables. - //! \param[in] n_f1 Number of nodal variables in field 1 - //! \param[in] n_f2 Number of nodal variables in field 2 + //! \param[in] n_f Number of nodal variables in each field ASMmxBase(const std::vector& n_f); //! \brief Initializes the patch level MADOF array. @@ -48,13 +47,13 @@ class ASMmxBase int basis = 0) const; //! \brief Injects nodal results for this patch into a global vector. - //! \param[in] globVec Global solution vector in DOF-order - //! \param[out] nodeVec Nodal result vector for this patch - //! \param[in] basis Which basis to extract the nodal values for + //! \param[out] globVec Global solution vector in DOF-order + //! \param[in] nodeVec Nodal result vector for this patch + //! \param[in] basis Which basis to inject the nodal values for void injectNodeVecMx(Vector& globVec, const Vector& nodeVec, int basis = 0) const; - //! \brief Extract the primary solution field at the specified nodes. + //! \brief Extracts the primary solution field at the specified nodes. //! \param[out] sField Solution field //! \param[in] locSol Solution vector local to current patch //! \param[in] nodes 1-based local node numbers to extract solution for @@ -62,10 +61,9 @@ class ASMmxBase const std::vector& nodes) const; public: - static char geoBasis; //!< The basis representing the geometry - - //! \brief Mixed formulation type + //! \brief Enum defining available mixed formulation types enum MixedType { + NONE = 0, //!< No fixed formulation FULL_CONT_RAISE_BASIS1, //!< Full continuity, raise order and use as basis 1 REDUCED_CONT_RAISE_BASIS1, //!< Reduced continuity, raise order and use as basis 1 FULL_CONT_RAISE_BASIS2, //!< Full continuity, raise order and use as basis 2 @@ -74,27 +72,29 @@ class ASMmxBase }; static MixedType Type; //!< Type of mixed formulation used + static char geoBasis; //!< 1-based index of basis representing the geometry +protected: typedef std::vector> SurfaceVec; //!< Convenience type typedef std::vector> VolumeVec; //!< Convenience type - //! \brief Establish mixed bases - //! \param[in] surf The base basis to use. - //! \param[in] type The type of bases to establish. - //! \return Vector with bases. + //! \brief Establish mixed bases in 2D. + //! \param[in] surf The base basis to use + //! \param[in] type The type of bases to establish + //! \return Vector with bases static SurfaceVec establishBases(Go::SplineSurface* surf, MixedType type); - //! \brief Establish mixed bases - //! \param[in] vol The base basis to use. - //! \param[in] type The type of bases to establish. - //! \return Vector with bases. + //! \brief Establish mixed bases in 3D. + //! \param[in] vol The base basis to use + //! \param[in] type The type of bases to establish + //! \return Vector with bases static VolumeVec establishBases(Go::SplineVolume* vol, MixedType type); private: std::vector MADOF; //!< Matrix of accumulated DOFs for this patch protected: - std::vector nb; //!< Number of basis functions in each basis + std::vector nb; //!< Number of basis functions in each basis std::vector nfx; //!< Number of fields on each basis }; From a032e1b8f2779a60d10a02ccdbd10ba8c65ba322 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 14 Apr 2016 19:02:48 +0200 Subject: [PATCH 4/6] Changed: Look one more level up for in-tree / submodule builds. The BeamSim and Nonlinear apps are now assumed to be sub-folders of the BeamEx and FiniteDeformation folders, respectively. --- Apps/CMakeLists.txt | 4 ++-- cmake/Modules/FindIFEM.cmake | 11 ++++++++--- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/Apps/CMakeLists.txt b/Apps/CMakeLists.txt index 8f362af92..1d12fc01b 100644 --- a/Apps/CMakeLists.txt +++ b/Apps/CMakeLists.txt @@ -35,11 +35,11 @@ if(EXISTS ${PROJECT_SOURCE_DIR}/Elasticity) add_subdirectory(Elasticity/Linear) if(EXISTS ${PROJECT_SOURCE_DIR}/Elasticity/BeamEx) add_subdirectory(Elasticity/BeamEx) - add_subdirectory(Elasticity/BeamSim) + add_subdirectory(Elasticity/BeamEx/BeamSim) endif() if(EXISTS ${PROJECT_SOURCE_DIR}/Elasticity/FiniteDeformation) add_subdirectory(Elasticity/FiniteDeformation) - add_subdirectory(Elasticity/Nonlinear) + add_subdirectory(Elasticity/FiniteDeformation/Nonlinear) endif() if(EXISTS ${PROJECT_SOURCE_DIR}/ThermoElasticity) add_subdirectory(ThermoElasticity) diff --git a/cmake/Modules/FindIFEM.cmake b/cmake/Modules/FindIFEM.cmake index 441408b50..bef5e735a 100644 --- a/cmake/Modules/FindIFEM.cmake +++ b/cmake/Modules/FindIFEM.cmake @@ -27,14 +27,18 @@ if(NOT IFEM_AS_SUBMODULE) NAMES ${IFEM_BUILD_TYPE}/IFEM.h PATHS ${PROJECT_SOURCE_DIR}/.. ${PROJECT_SOURCE_DIR}/../.. - ${PROJECT_SOURCE_DIR}/../../.. NO_DEFAULT_PATHS) + ${PROJECT_SOURCE_DIR}/../../.. + ${PROJECT_SOURCE_DIR}/../../../.. + NO_DEFAULT_PATHS) set(IFEM_H_PATH ${IFEM_PATH}/${IFEM_BUILD_TYPE}) if(NOT IFEM_PATH) find_path(IFEM_PATH NAMES IFEM.h PATHS ${PROJECT_SOURCE_DIR}/.. ${PROJECT_SOURCE_DIR}/../.. - ${PROJECT_SOURCE_DIR}/../../.. NO_DEFAULT_PATHS) + ${PROJECT_SOURCE_DIR}/../../.. + ${PROJECT_SOURCE_DIR}/../../../.. + NO_DEFAULT_PATHS) set(IFEM_H_PATH ${IFEM_PATH}) endif(NOT IFEM_PATH) endif() @@ -48,7 +52,8 @@ if(IFEM_AS_SUBMODULE) find_path(IFEM_PATH NAMES src/IFEM.h.in PATHS ${PROJECT_SOURCE_DIR}/.. ${PROJECT_SOURCE_DIR}/../.. - ${PROJECT_SOURCE_DIR}/../../.. NO_DEFAULT_PATHS) + ${PROJECT_SOURCE_DIR}/../../.. + ${PROJECT_SOURCE_DIR}/../../../.. NO_DEFAULT_PATHS) if(NOT IFEM_PATH) message(FATAL_ERROR "IFEM cannot be located, and we want it as a submodule") endif() From 348eb5832ae50b7bd8c6462728b3c9ce414cfb58 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Mon, 18 Apr 2016 20:31:41 +0200 Subject: [PATCH 5/6] 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 a711f41c1..738ee9691 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; } @@ -237,6 +237,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(); @@ -258,7 +260,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 @@ -300,6 +301,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; @@ -401,6 +403,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; @@ -415,7 +418,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; } @@ -424,7 +428,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 39067169d74c8c100cd29a2a27ca657c7f8e7cb9 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 19 Apr 2016 13:58:47 +0200 Subject: [PATCH 6/6] Added: operator=(const std::vector&), begin(), end() and expandCols(int) for utl::matrix. --- 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..f478f711d 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.expandCols(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.expandCols(-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..6849792b6 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& expandCols(int incRows) + { + int newRows = nrow + incRows; + if (newRows < 1) + { + // The matrix size 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; c > 1; c--, oldMat -= oldRows) + { + memcpy(this->ptr(c-1),oldMat,oldRows*sizeof(T)); + for (size_t r = nrow; r > oldRows; r--) + elem[r-1+nrow*(c-2)] = 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();