diff --git a/AssemblerLib/GlobalSetup.h b/AssemblerLib/GlobalSetup.h new file mode 100644 index 0000000000000000000000000000000000000000..282ecb5915bc4d457797b9c133e713c9cc236269 --- /dev/null +++ b/AssemblerLib/GlobalSetup.h @@ -0,0 +1,52 @@ +/** + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_GLOBALSETUP_H_ +#define ASSEMBLERLIB_GLOBALSETUP_H_ + +#include <functional> + +#include "MeshComponentMap.h" + +namespace AssemblerLib +{ + +/// The GlobalSetup collects vector and matrix builder and corresponding global +/// loop executor. +template <typename VectorMatrixBuilder, typename Executor> +struct GlobalSetup +{ + typedef typename VectorMatrixBuilder::VectorType VectorType; + typedef typename VectorMatrixBuilder::MatrixType MatrixType; + + template <typename... Args> + static + VectorType* createVector(Args&& ... args) + { + return VectorMatrixBuilder::createVector(std::forward<Args>(args)...); + } + + template <typename... Args> + static + MatrixType* createMatrix(Args&& ... args) + { + return VectorMatrixBuilder::createMatrix(std::forward<Args>(args)...); + } + + template <typename... Args> + static + void execute(Args&& ... args) + { + return Executor::execute(std::forward<Args>(args)...); + } +}; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_GLOBALSETUP_H_ diff --git a/AssemblerLib/LocalToGlobalIndexMap.h b/AssemblerLib/LocalToGlobalIndexMap.h new file mode 100644 index 0000000000000000000000000000000000000000..378ea13a1f8f8fde86b8cf4a50ec7109878ffc7e --- /dev/null +++ b/AssemblerLib/LocalToGlobalIndexMap.h @@ -0,0 +1,76 @@ +/** + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_LOCALTOGLOBALINDEXMAP_H_ +#define ASSEMBLERLIB_LOCALTOGLOBALINDEXMAP_H_ + +#include <vector> + +#include "MathLib/LinAlg/RowColumnIndices.h" + +namespace AssemblerLib +{ + +/// Row and column indices in global linear algebra objects for each mesh item. +/// +/// The row and column indices in global matrix and rhs vector for example, +/// required for addition of local contributions from every mesh item (like node +/// or cell) to global objects. +/// +/// The number of rows should be equal to the number of mesh items and the +/// number of columns should be equal to the number of the components on that +/// mesh item. +class LocalToGlobalIndexMap +{ +public: + typedef typename MathLib::RowColumnIndices<std::size_t> RowColumnIndices; + typedef RowColumnIndices::LineIndex LineIndex; + +public: + LocalToGlobalIndexMap( + std::vector<LineIndex> const& rows, + std::vector<LineIndex> const & columns) + : _rows(rows), _columns(columns) + { + assert(_rows.size() == _columns.size()); + assert(!_rows.empty()); + } + + explicit LocalToGlobalIndexMap(std::vector<LineIndex> const& rows) + : _rows(rows), _columns(rows) + { } + + std::size_t size() const + { + return _rows.size(); + } + + RowColumnIndices operator[](std::size_t const mesh_item_id) const + { + return RowColumnIndices(_rows[mesh_item_id], _columns[mesh_item_id]); + } + + LineIndex rowIndices(std::size_t const mesh_item_id) const + { + return _rows[mesh_item_id]; + } + + LineIndex columnIndices(std::size_t const mesh_item_id) const + { + return _columns[mesh_item_id]; + } + +private: + std::vector<LineIndex> const& _rows; + std::vector<LineIndex> const& _columns; +}; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_LOCALTOGLOBALINDEXMAP_H_ diff --git a/AssemblerLib/SerialDenseSetup.h b/AssemblerLib/SerialDenseSetup.h new file mode 100644 index 0000000000000000000000000000000000000000..fbac129d4c4eb8f43340e7e80114a231403e1dbd --- /dev/null +++ b/AssemblerLib/SerialDenseSetup.h @@ -0,0 +1,31 @@ +/** + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_SERIALDENSESETUP_H_ +#define ASSEMBLERLIB_SERIALDENSESETUP_H_ + + +#include "AssemblerLib/GlobalSetup.h" + +#include "AssemblerLib/SerialDenseVectorMatrixBuilder.h" +#include "AssemblerLib/SerialExecutor.h" + +namespace AssemblerLib +{ + +/// Using GlobalDenseMatrix and DenseVector for global entities and serial +/// global assembly loop. +typedef GlobalSetup< + AssemblerLib::SerialDenseVectorMatrixBuilder, + AssemblerLib::SerialExecutor> + SerialDenseSetup; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_SERIALDENSESETUP_H_ diff --git a/AssemblerLib/SerialDenseVectorMatrixBuilder.h b/AssemblerLib/SerialDenseVectorMatrixBuilder.h new file mode 100644 index 0000000000000000000000000000000000000000..db37bfaccd725acc0d7564d51b249c5fbd564846 --- /dev/null +++ b/AssemblerLib/SerialDenseVectorMatrixBuilder.h @@ -0,0 +1,36 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_SERIALDENSEVECTORMATRIXBUILDER_H_ +#define ASSEMBLERLIB_SERIALDENSEVECTORMATRIXBUILDER_H_ + +#include "AssemblerLib/SerialVectorMatrixBuilder.h" + +#include "MathLib/LinAlg/Dense/DenseVector.h" +#include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h" + +namespace AssemblerLib +{ + +/// Serial vector/matrix builder using the GlobalDenseMatrix and the DenseVector +/// implementations. +/// +/// \attention Using of GlobalDenseMatrix requires substantial amounts of +/// memory and this implementation is meant to be used for test purpose only. +typedef SerialVectorMatrixBuilder< + MathLib::GlobalDenseMatrix<double>, + MathLib::DenseVector<double> + > SerialDenseVectorMatrixBuilder; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_SERIALVECTOR_MATRIXBUILDER_H_ diff --git a/AssemblerLib/SerialExecutor.h b/AssemblerLib/SerialExecutor.h new file mode 100644 index 0000000000000000000000000000000000000000..edfac6e973e5803dbd08d997e8af2c320d107af5 --- /dev/null +++ b/AssemblerLib/SerialExecutor.h @@ -0,0 +1,43 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_SERIALEXECUTOR_H_H +#define ASSEMBLERLIB_SERIALEXECUTOR_H_H + +namespace AssemblerLib +{ + +struct SerialExecutor +{ + /// Executes a \c f for each element from the input container. + /// Return values of the function call are ignored. + /// + /// \tparam F \c f type. + /// \tparam C input container type. + /// + /// \param f a function that accepts a pointer to container's elements and + /// an index as arguments. + /// \param c a container supporting access over operator[]. + template <typename F, typename C> + static + void + execute(F const& f, C const& c) + { + for (std::size_t i = 0; i < c.size(); i++) + f(c[i], i); + }; + +}; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_SERIALEXECUTOR_H_H diff --git a/AssemblerLib/SerialLisVectorMatrixBuilder.h b/AssemblerLib/SerialLisVectorMatrixBuilder.h new file mode 100644 index 0000000000000000000000000000000000000000..43d1a1625d2fffd4dbff35ce9531d809f3023d39 --- /dev/null +++ b/AssemblerLib/SerialLisVectorMatrixBuilder.h @@ -0,0 +1,32 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_SERIALLISVECTORMATRIXBUILDER_H_ +#define ASSEMBLERLIB_SERIALLISVECTORMATRIXBUILDER_H_ + +#include "AssemblerLib/SerialVectorMatrixBuilder.h" + +#include "MathLib/LinAlg/Lis/LisMatrix.h" +#include "MathLib/LinAlg/Lis/LisVector.h" + +namespace AssemblerLib +{ + +/// Serial vector/matrix builder using LIS vectors and matrices. +typedef SerialVectorMatrixBuilder< + MathLib::LisMatrix, + MathLib::LisVector + > SerialLisVectorMatrixBuilder; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_SERIALLISVECTORMATRIXBUILDER_H_ diff --git a/AssemblerLib/SerialVectorMatrixBuilder.h b/AssemblerLib/SerialVectorMatrixBuilder.h new file mode 100644 index 0000000000000000000000000000000000000000..bf88cafcde8bdbeb07120b5b26d6dc4a5e06e00b --- /dev/null +++ b/AssemblerLib/SerialVectorMatrixBuilder.h @@ -0,0 +1,47 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-16 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_SERIALVECTORMATRIXBUILDER_H_ +#define ASSEMBLERLIB_SERIALVECTORMATRIXBUILDER_H_ + +#include "AssemblerLib/MeshComponentMap.h" + +namespace AssemblerLib +{ + +template <typename MatrixType_, typename VectorType_> +class SerialVectorMatrixBuilder +{ +public: + typedef VectorType_ VectorType; + typedef MatrixType_ MatrixType; + +public: + static + VectorType* createVector(const MeshComponentMap &dist_layout) + { + VectorType* vec = new VectorType(dist_layout.size()); + return vec; + } + + static + MatrixType* createMatrix(const MeshComponentMap &dist_layout) + { + MatrixType* mat = new MatrixType(dist_layout.size()); + return mat; + } + +}; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_SERIALVECTORMATRIXBUILDER_H_ diff --git a/AssemblerLib/VectorMatrixAssembler.h b/AssemblerLib/VectorMatrixAssembler.h new file mode 100644 index 0000000000000000000000000000000000000000..04ae0cd90d046c1c06932169d71d98304c211757 --- /dev/null +++ b/AssemblerLib/VectorMatrixAssembler.h @@ -0,0 +1,75 @@ +/** + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef ASSEMBLERLIB_VECTORMATRIXASSEMBLER_H_ +#define ASSEMBLERLIB_VECTORMATRIXASSEMBLER_H_ + +#include "LocalToGlobalIndexMap.h" + +namespace AssemblerLib +{ + +/// Adds result of local assembler into a global vector and a global matrix. +/// The VectorMatrixAssembler executes the local assembler for a given mesh item +/// and adds the local vector and matrix entries into the global vector and +/// the global matrix. The indices in global objects are provided by +/// the LocalToGlobalIndexMap in the construction. +template< + typename GLOBAL_MATRIX_, + typename GLOBAL_VECTOR_, + typename MESH_ITEM_, + typename ASSEMBLER_, + typename LOCAL_MATRIX_, + typename LOCAL_VECTOR_> +class VectorMatrixAssembler +{ +public: + typedef GLOBAL_MATRIX_ GLOBAL_MATRIX; + typedef GLOBAL_VECTOR_ GLOBAL_VECTOR; + typedef LOCAL_MATRIX_ LOCAL_MATRIX; + typedef LOCAL_VECTOR_ LOCAL_VECTOR; + +public: + VectorMatrixAssembler( + GLOBAL_MATRIX_ &A, + GLOBAL_VECTOR_ &rhs, + ASSEMBLER_ &local_assembler, + LocalToGlobalIndexMap const& data_pos) + : _A(A), _rhs(rhs), _local_assembler(local_assembler), _data_pos(data_pos) {} + + ~VectorMatrixAssembler() {} + + /// Executes local assembler for the given mesh item and adds the result + /// into the global matrix and vector. + /// The positions in the global matrix/vector are taken from + /// the LocalToGlobalIndexMap provided in the constructor at index \c id. + /// \attention The index \c id is not necesserily the mesh item's id. + void operator()(const MESH_ITEM_* item, std::size_t id) const + { + assert(_data_pos.size() > id); + + LocalToGlobalIndexMap::RowColumnIndices const& indices = _data_pos[id]; + LOCAL_MATRIX_ local_A(indices.rows.size(), indices.columns.size()); + LOCAL_VECTOR_ local_rhs(indices.rows.size()); + + _local_assembler(*item, local_A, local_rhs); + _A.add(indices, local_A); + _rhs.add(indices.rows, local_rhs); + } + +protected: + GLOBAL_MATRIX_ &_A; + GLOBAL_VECTOR_ &_rhs; + ASSEMBLER_ &_local_assembler; + LocalToGlobalIndexMap const& _data_pos; +}; + +} // namespace AssemblerLib + +#endif // ASSEMBLERLIB_VECTORMATRIXASSEMBLER_H_ diff --git a/MathLib/LinAlg/Dense/GlobalDenseMatrix.h b/MathLib/LinAlg/Dense/GlobalDenseMatrix.h index 7272c9888b8bc4f3806d52277ffbecbabad83c2e..1215d880bf254537e739b2ab2e0ef4a42c86c1ef 100644 --- a/MathLib/LinAlg/Dense/GlobalDenseMatrix.h +++ b/MathLib/LinAlg/Dense/GlobalDenseMatrix.h @@ -33,9 +33,7 @@ public: public: /// Dense square matrix constructor. - GlobalDenseMatrix (IDX_TYPE rows) - : GlobalDenseMatrix(rows, rows) - { }; + GlobalDenseMatrix (IDX_TYPE rows); /// Dense rectangular matrix constructor. GlobalDenseMatrix (IDX_TYPE rows, IDX_TYPE cols); diff --git a/MathLib/LinAlg/Dense/GlobalDenseMatrix.tpp b/MathLib/LinAlg/Dense/GlobalDenseMatrix.tpp index 2fbf20a3c86d24991c03ce3e2c5fab08d0ccef7e..d7248dedcaeba7841e7d2542d034b468cbd6ba29 100644 --- a/MathLib/LinAlg/Dense/GlobalDenseMatrix.tpp +++ b/MathLib/LinAlg/Dense/GlobalDenseMatrix.tpp @@ -14,6 +14,11 @@ namespace MathLib { +template<typename FP_TYPE, typename IDX_TYPE> +GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::GlobalDenseMatrix(IDX_TYPE rows) : + DenseMatrix<FP_TYPE, IDX_TYPE>(rows, rows) +{} + template<typename FP_TYPE, typename IDX_TYPE> GlobalDenseMatrix<FP_TYPE, IDX_TYPE>::GlobalDenseMatrix(IDX_TYPE rows, IDX_TYPE cols) : DenseMatrix<FP_TYPE, IDX_TYPE>(rows, cols) diff --git a/Tests/AssemblerLib/SteadyDiffusion2DExample1.h b/Tests/AssemblerLib/SteadyDiffusion2DExample1.h new file mode 100644 index 0000000000000000000000000000000000000000..c63f59816c674c2ef889b45d32c29fe5bbdf4fe6 --- /dev/null +++ b/Tests/AssemblerLib/SteadyDiffusion2DExample1.h @@ -0,0 +1,103 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-18 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef STEADYDIFFUSION2DEXAMPLE1_H_ +#define STEADYDIFFUSION2DEXAMPLE1_H_ + +#include <cmath> +#include <vector> + +#include "MathLib/LinAlg/Dense/DenseMatrix.h" +#include "MathLib/LinAlg/Dense/DenseVector.h" + +#include "MeshLib/Elements/Element.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/MeshGenerators/MeshGenerator.h" +#include "MeshLib/Node.h" + +struct SteadyDiffusion2DExample1 +{ + class LocalAssembler + { + public: + LocalAssembler() + : _m(4,4) + { + _m(0,0) = 4.0; _m(0,1) = -1.0; _m(0,2) = -2.0; _m(0,3) = -1.0; + _m(1,1) = 4.0; _m(1,2) = -1.0; _m(1,3) = -2.0; + _m(2,2) = 4.0; _m(2,3) = -1.0; + _m(3,3) = 4.0; + + // copy upper triangle to lower to make symmetric + for (std::size_t i = 0; i < 4; i++) + for (std::size_t j = 0; j < i; j++) + _m(i,j) = _m(j,i); + + //_m *= 1.e-11/6.0; + for (std::size_t i = 0; i < 4; i++) + for (std::size_t j = 0; j < 4; j++) + _m(i,j) *= 1.e-11 / 6.0; + } + + void operator()(const MeshLib::Element & /*e*/, + MathLib::DenseMatrix<double> &localA, + MathLib::DenseVector<double> & /*rhs*/) + { + for (std::size_t i = 0; i < 4; i++) + for (std::size_t j = 0; j < 4; j++) + localA(i,j) = _m(i,j); + } + + private: + MathLib::DenseMatrix<double> _m; + }; + + SteadyDiffusion2DExample1() + { + msh = MeshLib::MeshGenerator::generateRegularQuadMesh(2.0, mesh_subdivs); + for (auto* node : msh->getNodes()) + vec_nodeIDs.push_back(node->getID()); + vec_DirichletBC_id.resize(2 * mesh_stride); + for (std::size_t i = 0; i < mesh_stride; i++) + { + // left side + vec_DirichletBC_id[i] = i * mesh_stride; + + // right side + vec_DirichletBC_id[mesh_stride + i] = i * mesh_stride + mesh_subdivs; + } + + vec_DirichletBC_value.resize(2 * mesh_stride); + std::fill_n(vec_DirichletBC_value.begin(), mesh_stride, 0); + std::fill_n(vec_DirichletBC_value.begin() + mesh_stride, mesh_stride, 1); + exact_solutions.resize(dim_eqs); + for (std::size_t i = 0; i < mesh_stride; i++) + for (std::size_t j = 0; j < mesh_stride; j++) + exact_solutions[i*mesh_stride + j] = j * 1./mesh_subdivs; + } + + ~SteadyDiffusion2DExample1() + { + delete msh; + } + + static const std::size_t mesh_subdivs = 10; + static const std::size_t mesh_stride = mesh_subdivs + 1; + static const std::size_t dim_eqs = mesh_stride * mesh_stride; + MeshLib::Mesh* msh; + std::vector<std::size_t> vec_DirichletBC_id; + std::vector<double> vec_DirichletBC_value; + std::vector<double> exact_solutions; + std::vector<std::size_t> vec_nodeIDs; +}; + +#endif //STEADYDIFFUSION2DEXAMPLE1_H_ diff --git a/Tests/AssemblerLib/TestSerialExecutor.cpp b/Tests/AssemblerLib/TestSerialExecutor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..433398c7d363b694994d27b4dd729188682b5d27 --- /dev/null +++ b/Tests/AssemblerLib/TestSerialExecutor.cpp @@ -0,0 +1,71 @@ +/** + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include <gtest/gtest.h> + +#include <algorithm> +#include <vector> +#include "AssemblerLib/SerialExecutor.h" + +template <typename Container> +class AssemblerLibSerialExecutor : public ::testing::Test +{ + public: + AssemblerLibSerialExecutor() + { + reference.resize(size); + std::iota(reference.begin(), reference.end(), 0); + + std::copy(reference.begin(), reference.end(), + std::back_inserter(container)); + } + + void + subtractFromReference(int const value, std::size_t const index) + { + reference[index] -= value; + } + + bool + referenceIsZero() const + { + return std::all_of(reference.begin(), reference.end(), + [](int const reference_value) + { + return reference_value == 0; + }); + } + + static std::size_t const size = 100; + std::vector<int> reference; + Container container; +}; + +TYPED_TEST_CASE_P(AssemblerLibSerialExecutor); + +using namespace std::placeholders; + +TYPED_TEST_P(AssemblerLibSerialExecutor, ContainerArgument) +{ + AssemblerLib::SerialExecutor::execute( + std::bind(&TestFixture::subtractFromReference, this, _1, _2), + this->container); + ASSERT_TRUE(this->referenceIsZero()); +} + +REGISTER_TYPED_TEST_CASE_P(AssemblerLibSerialExecutor, + ContainerArgument); + + +typedef ::testing::Types < + std::vector<int> + > TestTypes; + +INSTANTIATE_TYPED_TEST_CASE_P(templated, AssemblerLibSerialExecutor, + TestTypes); diff --git a/Tests/AssemblerLib/TestSerialLinearSolver.cpp b/Tests/AssemblerLib/TestSerialLinearSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3c71c2a00c4f691eb24e1227ad555b511d31e281 --- /dev/null +++ b/Tests/AssemblerLib/TestSerialLinearSolver.cpp @@ -0,0 +1,140 @@ +/** + * \author Norihiro Watanabe + * \date 2013-04-18 + * + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include <cmath> +#include <memory> +#include <vector> + +#include <gtest/gtest.h> + +#include "AssemblerLib/VectorMatrixAssembler.h" +#include "AssemblerLib/MeshComponentMap.h" +#include "AssemblerLib/SerialDenseSetup.h" + + +#include "MathLib/LinAlg/Dense/DenseTools.h" +#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h" +#include "MathLib/LinAlg/FinalizeMatrixAssembly.h" +#include "MathLib/MathTools.h" + +#include "MeshLib/Elements/Element.h" +#include "MeshLib/Elements/Quad.h" +#include "MeshLib/Location.h" +#include "MeshLib/MeshSubsets.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/Node.h" + +#include "../TestTools.h" +#include "SteadyDiffusion2DExample1.h" + +TEST(AssemblerLibSerialLinearSolver, Steady2DdiffusionQuadElem) +{ + // example + SteadyDiffusion2DExample1 ex1; + + //-------------------------------------------------------------------------- + // Choose implementation type + //-------------------------------------------------------------------------- + typedef AssemblerLib::SerialDenseSetup GlobalSetup; + const GlobalSetup globalSetup; + + //-------------------------------------------------------------------------- + // Prepare mesh items where data are assigned + //-------------------------------------------------------------------------- + const MeshLib::MeshSubset mesh_items_all_nodes(*ex1.msh, + ex1.msh->getNodes()); + + //-------------------------------------------------------------------------- + // Allocate a coefficient matrix, RHS and solution vectors + //-------------------------------------------------------------------------- + // define a mesh item composition in a vector + std::vector<MeshLib::MeshSubsets*> vec_comp_dis; + vec_comp_dis.push_back( + new MeshLib::MeshSubsets(&mesh_items_all_nodes)); + AssemblerLib::MeshComponentMap vec1_composition( + vec_comp_dis, AssemblerLib::ComponentOrder::BY_COMPONENT); + + // allocate a vector and matrix + typedef GlobalSetup::VectorType GlobalVector; + typedef GlobalSetup::MatrixType GlobalMatrix; + std::unique_ptr<GlobalMatrix> A(globalSetup.createMatrix(vec1_composition)); + A->setZero(); + std::unique_ptr<GlobalVector> rhs(globalSetup.createVector(vec1_composition)); + std::unique_ptr<GlobalVector> x(globalSetup.createVector(vec1_composition)); + + //-------------------------------------------------------------------------- + // Construct a linear system + //-------------------------------------------------------------------------- + // create a mapping table from element nodes to entries in the linear system + auto &all_eles = ex1.msh->getElements(); + std::vector<std::vector<std::size_t> > map_ele_nodes2vec_entries; + map_ele_nodes2vec_entries.reserve(all_eles.size()); + for (auto e = all_eles.cbegin(); e != all_eles.cend(); ++e) + { + std::size_t const nnodes = (*e)->getNNodes(); + std::size_t const mesh_id = ex1.msh->getID(); + std::vector<MeshLib::Location> vec_items; + vec_items.reserve(nnodes); + for (std::size_t j = 0; j < nnodes; j++) + vec_items.emplace_back( + mesh_id, + MeshLib::MeshItemType::Node, + (*e)->getNode(j)->getID()); + + map_ele_nodes2vec_entries.push_back( + vec1_composition.getGlobalIndices + <AssemblerLib::ComponentOrder::BY_COMPONENT>(vec_items)); + } + + // Local and global assemblers. + typedef SteadyDiffusion2DExample1::LocalAssembler LocalAssembler; + LocalAssembler local_assembler; + + typedef AssemblerLib::VectorMatrixAssembler< + GlobalMatrix, GlobalVector, + MeshLib::Element, LocalAssembler, + MathLib::DenseMatrix<double>, + MathLib::DenseVector<double> + > GlobalAssembler; + + GlobalAssembler assembler(*A.get(), *rhs.get(), local_assembler, + AssemblerLib::LocalToGlobalIndexMap(map_ele_nodes2vec_entries)); + + // Call global assembler for each mesh element. + globalSetup.execute(assembler, ex1.msh->getElements()); + + //std::cout << "A=\n"; + //A->write(std::cout); + //std::cout << "rhs=\n"; + //rhs->write(std::cout); + + // apply Dirichlet BC + MathLib::applyKnownSolution(*A, *rhs, ex1.vec_DirichletBC_id, + ex1.vec_DirichletBC_value); + //std::cout << "A=\n"; + //A->write(std::cout); + //std::cout << "rhs=\n"; + //rhs->write(std::cout); + + MathLib::finalizeMatrixAssembly(*A); + //-------------------------------------------------------------------------- + // solve x=A^-1 rhs + //-------------------------------------------------------------------------- + MathLib::GaussAlgorithm<GlobalMatrix, GlobalVector> ls(*A); + ls.solve(*rhs, *x); + + double* px = &(*x)[0]; + ASSERT_DOUBLE_ARRAY_EQ(&ex1.exact_solutions[0], px, ex1.dim_eqs, 1.e-5); + + std::remove_if(vec_comp_dis.begin(), vec_comp_dis.end(), + [](MeshLib::MeshSubsets * p) { delete p; return true; }); +} diff --git a/Tests/AssemblerLib/TestSerialVectorMatrixBuilder.cpp b/Tests/AssemblerLib/TestSerialVectorMatrixBuilder.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7476391194c631c61d86539ca171af364617d66d --- /dev/null +++ b/Tests/AssemblerLib/TestSerialVectorMatrixBuilder.cpp @@ -0,0 +1,106 @@ +/** + * \copyright + * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include <gtest/gtest.h> + +#include "AssemblerLib/MeshComponentMap.h" + +#include "MeshLib/MeshGenerators/MeshGenerator.h" +#include "MeshLib/MeshSubsets.h" + +#include "AssemblerLib/SerialDenseVectorMatrixBuilder.h" + +#ifdef USE_LIS +#include "AssemblerLib/SerialLisVectorMatrixBuilder.h" +#endif // USE_LIS + +template <typename SerialBuilder> +class AssemblerLibSerialVectorMatrixBuilder : public ::testing::Test +{ + public: + typedef MeshLib::MeshItemType MeshItemType; + typedef MeshLib::Location Location; + typedef AssemblerLib::MeshComponentMap MeshComponentMap; + + typedef typename SerialBuilder::VectorType VectorType; + typedef typename SerialBuilder::MatrixType MatrixType; + + public: + AssemblerLibSerialVectorMatrixBuilder() + : mesh(nullptr), nodesSubset(nullptr), cmap(nullptr) + { + mesh = MeshLib::MeshGenerator::generateLineMesh(1.0, mesh_size); + nodesSubset = new MeshLib::MeshSubset(*mesh, mesh->getNodes()); + + // Add two components both based on the same nodesSubset. + components.emplace_back(new MeshLib::MeshSubsets(nodesSubset)); + components.emplace_back(new MeshLib::MeshSubsets(nodesSubset)); + + cmap = new MeshComponentMap(components, + AssemblerLib::ComponentOrder::BY_COMPONENT); + } + + ~AssemblerLibSerialVectorMatrixBuilder() + { + delete cmap; + std::remove_if(components.begin(), components.end(), + [](MeshLib::MeshSubsets* p) { delete p; return true; }); + delete nodesSubset; + delete mesh; + } + + static std::size_t const mesh_size = 9; + MeshLib::Mesh const* mesh; + MeshLib::MeshSubset const* nodesSubset; + + std::vector<MeshLib::MeshSubsets*> components; + MeshComponentMap const* cmap; +}; + +TYPED_TEST_CASE_P(AssemblerLibSerialVectorMatrixBuilder); + +TYPED_TEST_P(AssemblerLibSerialVectorMatrixBuilder, createVector) +{ + typedef typename TestFixture::VectorType V; + typedef TypeParam SerialBuilder; + V* v = SerialBuilder::createVector(*this->cmap); + + ASSERT_TRUE(v != nullptr); + ASSERT_EQ(this->cmap->size(), v->size()); + + delete v; +} + +TYPED_TEST_P(AssemblerLibSerialVectorMatrixBuilder, createMatrix) +{ + typedef typename TestFixture::MatrixType M; + typedef TypeParam SerialBuilder; + M* v = SerialBuilder::createMatrix(*this->cmap); + + ASSERT_TRUE(v != nullptr); + ASSERT_EQ(this->cmap->size(), v->getNRows()); + ASSERT_EQ(this->cmap->size(), v->getNCols()); + + delete v; +} + +REGISTER_TYPED_TEST_CASE_P(AssemblerLibSerialVectorMatrixBuilder, + createVector, createMatrix); + + +typedef ::testing::Types + < AssemblerLib::SerialDenseVectorMatrixBuilder +#ifdef USE_LIS + , AssemblerLib::SerialLisVectorMatrixBuilder +#endif // USE_LIS + > TestTypes; + +INSTANTIATE_TYPED_TEST_CASE_P(templated, AssemblerLibSerialVectorMatrixBuilder, + TestTypes); +