Commit 5c5731cd authored by Radu Popescu's avatar Radu Popescu
Browse files

Merge branch 'static_graph'

parents 3e658f51 948fa42a
......@@ -472,6 +472,21 @@ public:
DataType* const* const localValues,
Int format = Epetra_FECrsMatrix::COLUMN_MAJOR );
//! Add a set of values to the corresponding set of coefficient in the closed matrix
/*!
@param numRows Number of rows into the list given in "localValues"
@param numColumns Number of columns into the list given in "localValues"
@param rowIndices List of row indices
@param columnIndices List of column indices
@param localValues 2D array containing the coefficient related to "rowIndices" and "columnIndices"
@param format Format of the matrix (Epetra_FECrsMatrix::COLUMN_MAJOR or Epetra_FECrsMatrix::ROW_MAJOR)
*/
void sumIntoCoefficients ( Int const numRows, Int const numColumns,
std::vector<Int> const& rowIndices,
std::vector<Int> const& columnIndices,
DataType* const* const localValues,
Int format = Epetra_FECrsMatrix::COLUMN_MAJOR );
//! Add a value at a coefficient of the matrix
/*!
@param row Row index of the value to be added
......@@ -495,6 +510,12 @@ public:
//! @name Get Methods
//@{
//! Return the fill-complete status of the Epetra_FECrsMatrix
bool filled() const
{
return M_epetraCrs->Filled();
}
//! Return the shared_pointer of the Epetra_FECrsMatrix
matrix_ptrtype& matrixPtr()
{
......@@ -1513,6 +1534,23 @@ addToCoefficients ( Int const numRows, Int const numColumns,
}
template <typename DataType>
void MatrixEpetra<DataType>::
sumIntoCoefficients ( Int const numRows, Int const numColumns,
std::vector<Int> const& rowIndices, std::vector<Int> const& columnIndices,
DataType* const* const localValues,
Int format )
{
Int ierr = M_epetraCrs->SumIntoGlobalValues ( numRows, &rowIndices[0], numColumns,
&columnIndices[0], localValues, format );
std::stringstream errorMessage;
errorMessage << " error in matrix insertion [addToCoefficients] " << ierr
<< " when inserting in (" << rowIndices[0] << ", " << columnIndices[0] << ")" << std::endl;
ASSERT ( ierr >= 0, errorMessage.str() );
}
// ===================================================
// Get Methods
// ===================================================
......
......@@ -104,6 +104,11 @@ public:
std::vector<Int> const& blockRowIndices, std::vector<Int> const& blockColumnIndices,
DataType* const* const localValues,
Int format = Epetra_FECrsMatrix::COLUMN_MAJOR ) const;
//! Function to assemble an elemental matrix in a block (for closed matrices)
void sumIntoCoefficients ( UInt const numRows, UInt const numColumns,
std::vector<Int> const& blockRowIndices, std::vector<Int> const& blockColumnIndices,
DataType* const* const localValues,
Int format = Epetra_FECrsMatrix::COLUMN_MAJOR ) const;
//@}
//! @name Set Methods
......@@ -161,6 +166,12 @@ public:
return M_lastColumnIndex;
}
//! Return the fill-complete status of the inner Epetra_FECrsMatrix
bool filled() const
{
return M_matrix->matrixPtr()->Filled();
}
//! Return the pointer of the full matrix
matrix_Type* matrixPtr() const
{
......@@ -266,6 +277,31 @@ addToCoefficients ( UInt const numRows, UInt const numColumns,
localValues, format);
}
template<typename DataType>
void
MatrixEpetraStructuredView<DataType>::
sumIntoCoefficients ( UInt const numRows, UInt const numColumns,
std::vector<Int> const& blockRowIndices,
std::vector<Int> const& blockColumnIndices,
DataType* const* const localValues,
Int format) const
{
std::vector<Int> rowIndices (blockRowIndices);
std::vector<Int> columnIndices (blockColumnIndices);
for (UInt i (0); i < numRows; ++i)
{
rowIndices[i] += M_firstRowIndex;
}
for (UInt i (0); i < numColumns; ++i)
{
columnIndices[i] += M_firstColumnIndex;
}
M_matrix->sumIntoCoefficients (numRows, numColumns,
rowIndices, columnIndices,
localValues, format);
}
// ===================================================
// Set Methods
......
......@@ -120,6 +120,23 @@ public:
M_rawData, Epetra_FECrsMatrix::ROW_MAJOR);
}
//! Assembly procedure for a matrix or a block of a matrix
/*!
This method puts the values stored in this elemental matrix into the global
matrix passed as argument, using the positions given in the global indices
stored.
The method is used when the global matrix is closed
*/
// Method defined in class to allow compiler optimization
// as this class is used repeatedly during the assembly
template <typename MatrixType>
void pushToClosedGlobal (MatrixType& mat)
{
mat.sumIntoCoefficients ( M_nbRow, M_nbColumn,
rowIndices(), columnIndices(),
M_rawData, Epetra_FECrsMatrix::ROW_MAJOR);
}
//! Assembly procedure for a matrix or a block of a matrix passed in a shared_ptr
/*!
This method puts the values stored in this elemental matrix into the global
......@@ -138,6 +155,25 @@ public:
M_rawData, Epetra_FECrsMatrix::ROW_MAJOR);
}
//! Assembly procedure for a matrix or a block of a matrix passed in a shared_ptr
/*!
This method puts the values stored in this elemental matrix into the global
matrix passed as argument, using the positions given in the global indices
stored.
The method is used when the global matrix is closed
This is a partial specialization of the other assembly procedure of this class.
*/
// Method defined in class to allow compiler optimization
// as this class is used repeatedly during the assembly
template <typename MatrixType>
void pushToClosedGlobal (boost::shared_ptr<MatrixType> mat)
{
mat->sumIntoCoefficients ( M_nbRow, M_nbColumn,
rowIndices(), columnIndices(),
M_rawData, Epetra_FECrsMatrix::ROW_MAJOR);
}
//! Ouput method for the sizes and the stored values
void showMe ( std::ostream& out = std::cout ) const;
......
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory UNiversity
This file is part of the LifeV library
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, see <http://www.gnu.org/licenses/>
*******************************************************************************
*/
//@HEADER
/*!
* @file
@brief This file contains the definition of the buildGraph function.
This function is used to precompute the graph of a finite element
matrix, allowing the matrix to be build in closed, optimized form,
which makes the assembly procedure more efficient
@date 03/2013
@author Radu Popescu <radu.popescu@epfl.ch>
*/
#ifndef BUILD_GRAPH_HPP
#define BUILD_GRAPH_HPP
#include <lifev/core/LifeV.hpp>
#include <lifev/eta/expression/RequestLoopElement.hpp>
#include <lifev/core/fem/QuadratureRule.hpp>
#include <lifev/eta/expression/GraphElement.hpp>
#include <boost/shared_ptr.hpp>
namespace LifeV
{
/*!
\namespace ExpressionAssembly
Namespace for the assembly via expressions
*/
namespace ExpressionAssembly
{
//! Function to precompute the matrix graph
/*!
@author Radu Popescu <radu.popescu@epfl.ch>
This is a helper function to precompute the Crs graph used to build
a FECrsMatrix in closed optimized form
*/
template < typename MeshType,
typename TestSpaceType,
typename SolutionSpaceType,
typename ExpressionType >
GraphElement < MeshType,
TestSpaceType,
SolutionSpaceType,
ExpressionType >
buildGraph ( const RequestLoopElement<MeshType>& request,
const QuadratureRule& quadrature,
const boost::shared_ptr<TestSpaceType>& testSpace,
const boost::shared_ptr<SolutionSpaceType>& solutionSpace,
const ExpressionType& expression,
const UInt offsetUp = 0,
const UInt offsetLeft = 0);
template < typename MeshType,
typename TestSpaceType,
typename SolutionSpaceType,
typename ExpressionType >
GraphElement < MeshType,
TestSpaceType,
SolutionSpaceType,
ExpressionType >
buildGraph ( const RequestLoopElement<MeshType>& request,
const QuadratureRule& quadrature,
const boost::shared_ptr<TestSpaceType>& testSpace,
const boost::shared_ptr<SolutionSpaceType>& solutionSpace,
const ExpressionType& expression,
const UInt offsetUp,
const UInt offsetLeft)
{
return GraphElement < MeshType,
TestSpaceType,
SolutionSpaceType,
ExpressionType >
(request.mesh(), quadrature, testSpace, solutionSpace, expression,
offsetUp, offsetLeft );
}
} // Namespace ExpressionAssembly
} // Namespace LifeV
#endif // BUILD_GRAPH_HPP
SET(expression_HEADERS
expression/BuildGraph.hpp
expression/EvaluationAddition.hpp
expression/EvaluationDivI.hpp
expression/EvaluationDivision.hpp
......@@ -41,6 +42,7 @@ SET(expression_HEADERS
expression/ExpressionSubstraction.hpp
expression/ExpressionToEvaluation.hpp
expression/ExpressionVector.hpp
expression/GraphElement.hpp
expression/Integrate.hpp
expression/IntegrateMatrixElement.hpp
expression/IntegrateValueElement.hpp
......
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory UNiversity
This file is part of the LifeV library
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, see <http://www.gnu.org/licenses/>
*******************************************************************************
*/
//@HEADER
/*!
* @file
@brief This file contains the definition of the GraphElement class.
@date 06/2011
@author Radu Popescu <radu.popescu@epfl.ch>
*/
#ifndef GRAPH_ELEMENT_HPP
#define GRAPH_ELEMENT_HPP
#include <lifev/core/LifeV.hpp>
#include <lifev/core/fem/QuadratureRule.hpp>
#include <lifev/eta/fem/ETCurrentFE.hpp>
#include <lifev/eta/fem/MeshGeometricMap.hpp>
#include <lifev/eta/expression/ExpressionToEvaluation.hpp>
#include <lifev/eta/array/ETMatrixElemental.hpp>
#include <boost/shared_ptr.hpp>
namespace LifeV
{
namespace ExpressionAssembly
{
//! The class to actually perform the loop over the elements to precompute a graph
/*!
@author Radu Popescu <radu.popescu@epfl.ch>
This class is used to store the data required for building the graph of a
matrix
*/
template < typename MeshType,
typename TestSpaceType,
typename SolutionSpaceType,
typename ExpressionType >
class GraphElement
{
public:
//! @name Public Types
//@{
//! Type of the Evaluation
typedef typename ExpressionToEvaluation < ExpressionType,
TestSpaceType::S_fieldDim,
SolutionSpaceType::S_fieldDim,
3 >::evaluation_Type evaluation_Type;
//@}
//! @name Constructors, destructor
//@{
//! Full data constructor
GraphElement (const boost::shared_ptr<MeshType>& mesh,
const QuadratureRule& quadrature,
const boost::shared_ptr<TestSpaceType>& testSpace,
const boost::shared_ptr<SolutionSpaceType>& solutionSpace,
const ExpressionType& expression,
const UInt offsetUp = 0,
const UInt offsetLeft = 0);
//! Copy constructor
GraphElement (const GraphElement < MeshType,
TestSpaceType,
SolutionSpaceType,
ExpressionType > & integrator);
//! Destructor
~GraphElement();
//@}
//! @name Operators
//@{
//! Operator wrapping the addTo method (for shared_ptr)
template <typename GraphType>
inline void operator>> (boost::shared_ptr<GraphType> graph)
{
addTo (graph);
}
//@}
//! @name Methods
//@{
//! Ouput method
void check (std::ostream& out = std::cout);
//! Method that builds the graph
/*!
The loop over the elements is located right
in this method.
*/
template <typename GraphType>
void addTo (GraphType& graph);
//! Method that performs the assembly
/*!
The loop over the elements is located right
in this method.
Specialized for the case where the matrix is passed as a shared_ptr
*/
template <typename GraphType>
inline void addTo (boost::shared_ptr<GraphType> graph)
{
ASSERT (graph != 0, " Cannot assemble with an empty graph");
addTo (*graph);
}
//@}
private:
//! @name Private Methods
//@{
//! No empty constructor
GraphElement();
//@}
// Pointer on the mesh
boost::shared_ptr<MeshType> M_mesh;
// Quadrature to be used
QuadratureRule M_quadrature;
// Shared pointer on the Spaces
boost::shared_ptr<TestSpaceType> M_testSpace;
boost::shared_ptr<SolutionSpaceType> M_solutionSpace;
// Offsets
UInt M_offsetUp;
UInt M_offsetLeft;
};
// ===================================================
// IMPLEMENTATION
// ===================================================
// ===================================================
// Constructors & Destructor
// ===================================================
template < typename MeshType, typename TestSpaceType, typename SolutionSpaceType, typename ExpressionType>
GraphElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>::
GraphElement (const boost::shared_ptr<MeshType>& mesh,
const QuadratureRule& quadrature,
const boost::shared_ptr<TestSpaceType>& testSpace,
const boost::shared_ptr<SolutionSpaceType>& solutionSpace,
const ExpressionType& /*expression*/,
const UInt offsetUp,
const UInt offsetLeft)
: M_mesh (mesh),
M_quadrature (quadrature),
M_testSpace (testSpace),
M_solutionSpace (solutionSpace),
M_offsetUp (offsetUp),
M_offsetLeft (offsetLeft)
{
}
template < typename MeshType, typename TestSpaceType, typename SolutionSpaceType, typename ExpressionType>
GraphElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>::
GraphElement (const GraphElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>& integrator)
: M_mesh (integrator.M_mesh),
M_quadrature (integrator.M_quadrature),
M_testSpace (integrator.M_testSpace),
M_solutionSpace (integrator.M_solutionSpace),
M_offsetUp (integrator.M_offsetUp),
M_offsetLeft (integrator.M_offsetLeft)
{
}
template < typename MeshType, typename TestSpaceType, typename SolutionSpaceType, typename ExpressionType>
GraphElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>::
~GraphElement()
{
}
// ===================================================
// Methods
// ===================================================
template < typename MeshType, typename TestSpaceType, typename SolutionSpaceType, typename ExpressionType>
void
GraphElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>::
check (std::ostream& out)
{
}
template < typename MeshType, typename TestSpaceType, typename SolutionSpaceType, typename ExpressionType>
template <typename GraphType>
void
GraphElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>::
addTo (GraphType& graph)
{
UInt nbElements (M_mesh->numElements() );
UInt nbTestDof (M_testSpace->refFE().nbDof() );
UInt nbSolutionDof (M_solutionSpace->refFE().nbDof() );
ETMatrixElemental elementalMatrix (TestSpaceType::S_fieldDim * M_testSpace->refFE().nbDof(),
SolutionSpaceType::S_fieldDim * M_solutionSpace->refFE().nbDof() );
for (UInt iElement = 0; iElement < nbElements; ++iElement)
{
// Zeros out the matrix
elementalMatrix.zero();
// Loop on the blocks
for (UInt iblock (0); iblock < TestSpaceType::S_fieldDim; ++iblock)
{
for (UInt jblock (0); jblock < SolutionSpaceType::S_fieldDim; ++jblock)
{
// Set the row global indices in the local matrix
for (UInt i (0); i < nbTestDof; ++i)
{
elementalMatrix.setRowIndex
(i + iblock * nbTestDof,
M_testSpace->dof().localToGlobalMap (iElement, i) + iblock * M_testSpace->dof().numTotalDof() + M_offsetUp);
}
// Set the column global indices in the local matrix
for (UInt j (0); j < nbSolutionDof; ++j)
{
elementalMatrix.setColumnIndex
(j + jblock * nbSolutionDof,
M_solutionSpace->dof().localToGlobalMap (iElement, j) + jblock * M_solutionSpace->dof().numTotalDof() + M_offsetLeft);
}
}
}
const std::vector<Int>& rowIdx = elementalMatrix.rowIndices();
const std::vector<Int>& colIdx = elementalMatrix.columnIndices();
graph.InsertGlobalIndices (rowIdx.size(), &rowIdx[0],
colIdx.size(), &colIdx[0]);
}
}
} // Namespace ExpressionAssembly
} // Namespace LifeV
#endif
......@@ -74,10 +74,21 @@ integrate ( const RequestLoopElement<MeshType>& request,
const QuadratureRule& quadrature,
const boost::shared_ptr<TestSpaceType>& testSpace,
const boost::shared_ptr<SolutionSpaceType>& solutionSpace,
const ExpressionType& expression)
const ExpressionType& expression,
const UInt offsetUp = 0,
const UInt offsetLeft = 0);
template < typename MeshType, typename TestSpaceType, typename SolutionSpaceType, typename ExpressionType>
IntegrateMatrixElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>
integrate ( const RequestLoopElement<MeshType>& request,
const QuadratureRule& quadrature,
const boost::shared_ptr<TestSpaceType>& testSpace,
const boost::shared_ptr<SolutionSpaceType>& solutionSpace,
const ExpressionType& expression,
const UInt offsetUp,
const UInt offsetLeft)
{
return IntegrateMatrixElement<MeshType, TestSpaceType, SolutionSpaceType, ExpressionType>
(request.mesh(), quadrature, testSpace, solutionSpace, expression);
(request.mesh(), quadrature, testSpace, solutionSpace, expression, offsetUp, offsetLeft);
}
//! Integrate function for vectorial expressions
......@@ -93,10 +104,18 @@ IntegrateVectorElement<MeshType, TestSpaceType, ExpressionType>
integrate ( const RequestLoopElement<MeshType>& request,
const QuadratureRule& quadrature,