construct_system_matrix_elemental_process.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Riccardo Rossi
//
 
 
 
#if !defined(KRATOS_CONSTRUCT_SYSTEM_MATRIX_ELEMENTAL_PROCESS_H_INCLUDED )
#define  KRATOS_CONSTRUCT_SYSTEM_MATRIX_ELEMENTAL_PROCESS_H_INCLUDED
 
 
 
// System includes
#include <string>
#include <iostream>
 
 
// External includes
 
 
// Project includes
#include "includes/define.h"
#include "processes/process.h"
#include "containers/pointer_vector_set.h"
 
 
namespace Kratos
{
 
 
 
 
 
 
 
template<class TSystemSpaceType, class TDirichletSpaceType = TSystemSpaceType>
class ConstructSystemMatrixElementalProcess : public Process
{
public:
 
    KRATOS_CLASS_POINTER_DEFINITION(ConstructSystemMatrixElementalProcess);
 
    typedef typename TSystemSpaceType::MatrixType SystemMatrixType;
 
    typedef typename TDirichletSpaceType::MatrixType DirichletMatrixType;
 
    typedef Element ElementType;
 
    typedef PointerVectorSet<ElementType, IndexedObject> ElementsContainerType;
 
 
    ConstructSystemMatrixElementalProcess(ElementsContainerType & rElements, SystemMatrixType& rSystemMatrix, DirichletMatrixType& rDirichletMatrix)
        : mElements(rElements), mSystemMatrix(rSystemMatrix), mDirichletMatrix(rDirichletMatrix) {}
 
    virtual ~ConstructSystemMatrixElementalProcess() {}
 
 
 
 
 
    virtual void Execute()
    {
        std::size_t equation_size = TSystemSpaceType::Size1(mSystemMatrix);
        std::vector<std::vector<std::size_t> > indices(equation_size);
        std::vector<std::vector<std::size_t> > dirichlet_indices(TSystemSpaceType::Size1(mDirichletMatrix));
 
        ProcessInfo process_info;
        for(ElementsContainerType::iterator i_element = mElements.begin() ; i_element != mElements.end() ; i_element++)
        {
            Element::EquationIdVectorType ids(3,0);
 
            i_element->EquationIdVector(ids, process_info);
            for(std::size_t i = 0 ; i < 3 ; i++)
                if(ids[i] < equation_size)
                    for(std::size_t j = 0 ; j < 3 ; j++)
                        if(ids[j] < equation_size)
                            indices[ids[i]].push_back(ids[j]);
                        else
                            dirichlet_indices[ids[i]].push_back(ids[j]- equation_size);
 
        }
        for(std::size_t i = 0 ; i < indices.size() ; i++)
        {
            std::vector<std::size_t>& row_indices = indices[i];
            std::vector<std::size_t>& dirichlet_row_indices = dirichlet_indices[i];
            std::sort(row_indices.begin(), row_indices.end());
            std::unique(row_indices.begin(), row_indices.end());
            std::sort(dirichlet_row_indices.begin(), dirichlet_row_indices.end());
            std::unique(dirichlet_row_indices.begin(), dirichlet_row_indices.end());
            for(std::size_t j = 0 ; j < row_indices.size() ; j++)
                mSystemMatrix()(i,row_indices[j]) = 0.00;
            for(std::size_t j = 0 ; j < dirichlet_row_indices.size() ; j++)
                mDirichletMatrix()(i,dirichlet_row_indices[j]) = 0.00;
        }
        for(std::size_t i = 0 ; i < indices.size() ; i++) // To make put at least one element per each row. Pooyan.
            mDirichletMatrix()(i,0) = 0.00;
    }
 
 
 
 
 
 
 
    virtual std::string Info() const
    {
        return "ConstructSystemMatrixElementalProcess";
    }
 
    virtual void PrintInfo(std::ostream& rOStream) const
    {
        rOStream << "ConstructSystemMatrixElementalProcess";
    }
 
    virtual void PrintData(std::ostream& rOStream) const
    {
        rOStream << "    System matrix non zeros : " << mSystemMatrix().non_zeros() << std::endl;
        rOStream << "    Dirichlet matrix non zeros : " << mDirichletMatrix().non_zeros() << std::endl;
    }
 
 
 
 
 
protected:
 
 
 
    ElementsContainerType& mElements;
    SystemMatrixType& mSystemMatrix;
    DirichletMatrixType& mDirichletMatrix;
 
 
 
 
 
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
 
 
 
 
 
    ConstructSystemMatrixElementalProcess& operator=(ConstructSystemMatrixElementalProcess const& rOther);
 
    ConstructSystemMatrixElementalProcess(ConstructSystemMatrixElementalProcess const& rOther);
 
 
 
}; // Class ConstructSystemMatrixElementalProcess
 
 
 
 
 
 
 
template<class TSystemSpaceType, class TDirichletSpaceType>
inline std::istream& operator >> (std::istream& rIStream,
                                  ConstructSystemMatrixElementalProcess<TSystemSpaceType, TDirichletSpaceType>& rThis);
 
template<class TSystemSpaceType, class TDirichletSpaceType>
inline std::ostream& operator << (std::ostream& rOStream,
                                  const ConstructSystemMatrixElementalProcess<TSystemSpaceType, TDirichletSpaceType>& rThis)
{
    rThis.PrintInfo(rOStream);
    rOStream << std::endl;
    rThis.PrintData(rOStream);
 
    return rOStream;
}
 
 
}  // namespace Kratos.
 
#endif // KRATOS_CONSTRUCT_SYSTEM_MATRIX_ELEMENTAL_PROCESS_H_INCLUDED  defined