Commit cd4895db authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

New structure for PETSc solver objects.

parent aee5a801
......@@ -36,7 +36,8 @@ else
endif
PARMETIS_LIB = -L$(PARMETIS_DIR) -lparmetis -lmetis
ZOLTAN_LIB = -L$(AMDIS_DIR)/lib/zoltan_build/lib -lzoltan
#ZOLTAN_LIB = -L$(AMDIS_DIR)/lib/zoltan_build/lib -lzoltan
ZOLTAN_LIB =
LIBS += $(AMDIS_LIB) $(PNG_LIB)
LIBS += -lboost_iostreams -lboost_filesystem -lboost_system -lboost_date_time
......
......@@ -399,7 +399,7 @@ namespace AMDiS {
void serialize(std::ostream &out);
/// Reads a matrix from an input stream.
void deserialize(::std::istream &in);
void deserialize(std::istream &in);
///
int memsize();
......
......@@ -52,7 +52,7 @@ namespace AMDiS {
public:
/// The constructor reads needed parameters and sets solvers \ref name.
ITL_OEMSolver(::std::string name) : OEMSolver(name) {}
ITL_OEMSolver(std::string name) : OEMSolver(name) {}
~ITL_OEMSolver() {}
......@@ -102,7 +102,7 @@ namespace AMDiS {
typedef DOFMatrix::value_type value_type;
public:
/// The constructor reads needed parameters and sets solvers \ref name.
ITL_OEMSolver_para(::std::string name)
ITL_OEMSolver_para(std::string name)
: OEMSolver(name), ell(OEMSolver::max_iter)
{
GET_PARAMETER(0, name + "->ell", "%d", &ell);
......
......@@ -20,8 +20,8 @@ if USE_PARALLEL_DOMAIN_AMDIS
parallel/ParallelDebug.cc \
parallel/ParallelProblemStatBase.cc \
parallel/ParMetisPartitioner.cc \
parallel/PetscProblemStat.cc \
parallel/PetscSolver.cc \
parallel/PetscSolverSchur.cc \
parallel/StdMpi.cc
libamdis_la_CXXFLAGS += -DHAVE_PARALLEL_DOMAIN_AMDIS=1
AMDIS_INCLUDES += -I$(PARMETIS_DIR)
......@@ -252,8 +252,8 @@ parallel/MpiHelper.h \
parallel/ParallelDebug.h \
parallel/ParallelProblemStatBase.h \
parallel/ParMetisPartitioner.h \
parallel/PetscProblemStat.h \
parallel/PetscSolver.h \
parallel/PetscSolverSchur.h \
parallel/StdMpi.h \
parallel/ZoltanPartitioner.h \
reinit/BoundaryElementDist.h \
......
......@@ -46,8 +46,8 @@ host_triplet = @host@
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/ParallelDebug.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/ParallelProblemStatBase.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/ParMetisPartitioner.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/PetscProblemStat.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/PetscSolver.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/PetscSolverSchur.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ parallel/StdMpi.cc
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@am__append_2 = -DHAVE_PARALLEL_DOMAIN_AMDIS=1
......@@ -91,8 +91,8 @@ am__libamdis_la_SOURCES_DIST = parallel/ElementObjectData.cc \
parallel/MeshDistributor.cc parallel/MeshManipulation.cc \
parallel/MeshPartitioner.cc parallel/MpiHelper.cc \
parallel/ParallelDebug.cc parallel/ParallelProblemStatBase.cc \
parallel/ParMetisPartitioner.cc parallel/PetscSolver.cc \
parallel/PetscSolverSchur.cc parallel/StdMpi.cc \
parallel/ParMetisPartitioner.cc parallel/PetscProblemStat.cc \
parallel/PetscSolver.cc parallel/StdMpi.cc \
parallel/ZoltanPartitioner.cc AdaptBase.cc AdaptInfo.cc \
AdaptInstationary.cc AdaptStationary.cc Assembler.cc \
BasisFunction.cc Boundary.cc BoundaryManager.cc Cholesky.cc \
......@@ -137,8 +137,8 @@ am__libamdis_la_SOURCES_DIST = parallel/ElementObjectData.cc \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelDebug.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParallelProblemStatBase.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-ParMetisPartitioner.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-PetscProblemStat.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-PetscSolver.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-PetscSolverSchur.lo \
@USE_PARALLEL_DOMAIN_AMDIS_TRUE@ libamdis_la-StdMpi.lo
@USE_ZOLTAN_TRUE@am__objects_2 = libamdis_la-ZoltanPartitioner.lo
am__objects_3 = $(am__objects_1) $(am__objects_2)
......@@ -574,8 +574,8 @@ parallel/MpiHelper.h \
parallel/ParallelDebug.h \
parallel/ParallelProblemStatBase.h \
parallel/ParMetisPartitioner.h \
parallel/PetscProblemStat.h \
parallel/PetscSolver.h \
parallel/PetscSolverSchur.h \
parallel/StdMpi.h \
parallel/ZoltanPartitioner.h \
reinit/BoundaryElementDist.h \
......@@ -879,8 +879,8 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Parameters.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Parametric.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PeriodicBC.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PetscProblemStat.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PetscSolver.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PetscSolverSchur.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PngWriter.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-PovrayWriter.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemImplicit.Plo@am__quote@
......@@ -1029,6 +1029,13 @@ libamdis_la-ParMetisPartitioner.lo: parallel/ParMetisPartitioner.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-ParMetisPartitioner.lo `test -f 'parallel/ParMetisPartitioner.cc' || echo '$(srcdir)/'`parallel/ParMetisPartitioner.cc
libamdis_la-PetscProblemStat.lo: parallel/PetscProblemStat.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-PetscProblemStat.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-PetscProblemStat.Tpo" -c -o libamdis_la-PetscProblemStat.lo `test -f 'parallel/PetscProblemStat.cc' || echo '$(srcdir)/'`parallel/PetscProblemStat.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-PetscProblemStat.Tpo" "$(DEPDIR)/libamdis_la-PetscProblemStat.Plo"; else rm -f "$(DEPDIR)/libamdis_la-PetscProblemStat.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='parallel/PetscProblemStat.cc' object='libamdis_la-PetscProblemStat.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-PetscProblemStat.lo `test -f 'parallel/PetscProblemStat.cc' || echo '$(srcdir)/'`parallel/PetscProblemStat.cc
libamdis_la-PetscSolver.lo: parallel/PetscSolver.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-PetscSolver.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-PetscSolver.Tpo" -c -o libamdis_la-PetscSolver.lo `test -f 'parallel/PetscSolver.cc' || echo '$(srcdir)/'`parallel/PetscSolver.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-PetscSolver.Tpo" "$(DEPDIR)/libamdis_la-PetscSolver.Plo"; else rm -f "$(DEPDIR)/libamdis_la-PetscSolver.Tpo"; exit 1; fi
......@@ -1036,13 +1043,6 @@ libamdis_la-PetscSolver.lo: parallel/PetscSolver.cc
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-PetscSolver.lo `test -f 'parallel/PetscSolver.cc' || echo '$(srcdir)/'`parallel/PetscSolver.cc
libamdis_la-PetscSolverSchur.lo: parallel/PetscSolverSchur.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-PetscSolverSchur.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-PetscSolverSchur.Tpo" -c -o libamdis_la-PetscSolverSchur.lo `test -f 'parallel/PetscSolverSchur.cc' || echo '$(srcdir)/'`parallel/PetscSolverSchur.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-PetscSolverSchur.Tpo" "$(DEPDIR)/libamdis_la-PetscSolverSchur.Plo"; else rm -f "$(DEPDIR)/libamdis_la-PetscSolverSchur.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='parallel/PetscSolverSchur.cc' object='libamdis_la-PetscSolverSchur.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCXX_FALSE@ $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -c -o libamdis_la-PetscSolverSchur.lo `test -f 'parallel/PetscSolverSchur.cc' || echo '$(srcdir)/'`parallel/PetscSolverSchur.cc
libamdis_la-StdMpi.lo: parallel/StdMpi.cc
@am__fastdepCXX_TRUE@ if $(LIBTOOL) --tag=CXX --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libamdis_la_CXXFLAGS) $(CXXFLAGS) -MT libamdis_la-StdMpi.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-StdMpi.Tpo" -c -o libamdis_la-StdMpi.lo `test -f 'parallel/StdMpi.cc' || echo '$(srcdir)/'`parallel/StdMpi.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-StdMpi.Tpo" "$(DEPDIR)/libamdis_la-StdMpi.Plo"; else rm -f "$(DEPDIR)/libamdis_la-StdMpi.Tpo"; exit 1; fi
......
......@@ -39,6 +39,8 @@
namespace AMDiS {
using namespace std;
struct OperatorPos
{
int row, col;
......@@ -51,14 +53,14 @@ namespace AMDiS {
{
protected:
// Defines a mapping type from dof indices to world coordinates.
typedef std::map<DegreeOfFreedom, WorldVector<double> > DofToCoord;
typedef map<DegreeOfFreedom, WorldVector<double> > DofToCoord;
// Defines a mapping type from dof indices to world coordinates.
typedef std::map<WorldVector<double>, DegreeOfFreedom> CoordToDof;
typedef map<WorldVector<double>, DegreeOfFreedom> CoordToDof;
public:
/// Constructor
ProblemVec(std::string nameStr,
ProblemVec(string nameStr,
ProblemIterationInterface *problemIteration = NULL)
: StandardProblemIteration(this),
name(nameStr),
......@@ -227,7 +229,7 @@ namespace AMDiS {
void writeFiles(AdaptInfo &adaptInfo, bool force);
/// Interpolates fct to \ref solution.
void interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct);
void interpolInitialSolution(vector<AbstractFunction<double, WorldVector<double> >*> *fct);
/// Adds an operator to \ref A.
void addMatrixOperator(Operator *op, int i, int j,
......@@ -346,7 +348,7 @@ namespace AMDiS {
}
/// Returns \ref meshes
inline std::vector<Mesh*> getMeshes()
inline vector<Mesh*> getMeshes()
{
return meshes;
}
......@@ -361,19 +363,19 @@ namespace AMDiS {
}
/// Returns \ref feSpaces.
inline std::vector<FiniteElemSpace*> getFeSpaces()
inline vector<FiniteElemSpace*> getFeSpaces()
{
return feSpaces;
}
/// Returns \ref componentSpaces;
inline std::vector<FiniteElemSpace*> getComponentFeSpaces()
inline vector<FiniteElemSpace*> getComponentFeSpaces()
{
return componentSpaces;
}
/// Returns \ref estimator.
inline std::vector<Estimator*> getEstimators()
inline vector<Estimator*> getEstimators()
{
return estimator;
}
......@@ -409,13 +411,13 @@ namespace AMDiS {
}
/// Returns \ref marker.
inline std::vector<Marker*> getMarkers()
inline vector<Marker*> getMarkers()
{
return marker;
}
/// Returns the name of the problem
inline virtual std::string getName()
inline virtual string getName()
{
return name;
}
......@@ -445,7 +447,7 @@ namespace AMDiS {
*/
/// Sets \ref estimator.
inline void setEstimator(std::vector<Estimator*> est)
inline void setEstimator(vector<Estimator*> est)
{
estimator = est;
}
......@@ -511,17 +513,17 @@ namespace AMDiS {
* Outputs the mesh of the given component, but the values are taken from the
* residual error estimator.
*/
void writeResidualMesh(int comp, AdaptInfo *adaptInfo, std::string name);
void writeResidualMesh(int comp, AdaptInfo *adaptInfo, string name);
/// Function that implements the serialization procedure.
virtual void serialize(std::ostream &out);
virtual void serialize(ostream &out);
/// Function that implements the deserialization procedure.
virtual void deserialize(std::istream &in);
virtual void deserialize(istream &in);
/// Returns \ref fileWriters.
std::vector<FileWriterInterface*>& getFileWriterList()
vector<FileWriterInterface*>& getFileWriterList()
{
return fileWriters;
}
......@@ -536,7 +538,7 @@ namespace AMDiS {
protected:
/// Name of this problem.
std::string name;
string name;
/// Number of problem components
int nComponents;
......@@ -549,16 +551,16 @@ namespace AMDiS {
int nMeshes;
/// FE spaces of this problem.
std::vector<FiniteElemSpace*> feSpaces;
vector<FiniteElemSpace*> feSpaces;
/// Meshes of this problem.
std::vector<Mesh*> meshes;
vector<Mesh*> meshes;
/// Pointer to the fe spaces for the different problem components
std::vector<FiniteElemSpace*> componentSpaces;
vector<FiniteElemSpace*> componentSpaces;
/// Pointer to the meshes for the different problem components
std::vector<Mesh*> componentMeshes;
vector<Mesh*> componentMeshes;
/** \brief
* Stores information about which meshes must be traversed to assemble the
......@@ -568,10 +570,10 @@ namespace AMDiS {
ComponentTraverseInfo traverseInfo;
/// Responsible for element marking.
std::vector<Marker*> marker;
vector<Marker*> marker;
/// Estimator of this problem. Used in \ref estimate().
std::vector<Estimator*> estimator;
vector<Estimator*> estimator;
/// Linear solver of this problem. Used in \ref solve().
OEMSolver *solver;
......@@ -595,7 +597,7 @@ namespace AMDiS {
* be assembled only once. All other matrices will be assembled at every
* time step.
*/
std::vector< std::vector< bool > > assembleMatrixOnlyOnce;
vector<vector<bool> > assembleMatrixOnlyOnce;
/** \brief
* If [i][j] of this field is set to true, the corresponding DOFMatrix of
......@@ -603,13 +605,13 @@ namespace AMDiS {
* to determine, if assembling of a matrix can be ommitted, if it is set
* to be assembled only once.
*/
std::vector< std::vector< bool > > assembledMatrix;
vector<vector<bool> > assembledMatrix;
/// Determines whether domain boundaries should be considered at assembling.
bool useGetBound;
/// Writes the meshes and solution after the adaption loop.
std::vector<FileWriterInterface*> fileWriters;
vector<FileWriterInterface*> fileWriters;
/** \brief
* All actions of mesh refinement are performed by refinementManager.
......@@ -636,7 +638,7 @@ namespace AMDiS {
* the problem. This may be used to compute the real error of the computed
* solution.
*/
std::vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
vector<AbstractFunction<double, WorldVector<double> >*> exactSolutionFcts;
/** \brief
* If true, the error is not estimated but computed from the exact solution
......@@ -654,9 +656,9 @@ namespace AMDiS {
/// If true, AMDiS prints information about the assembling procedure to the screen.
bool writeAsmInfo;
std::map<Operator*, std::vector<OperatorPos> > operators;
map<Operator*, vector<OperatorPos> > operators;
std::map<Operator*, Flag> opFlags;
map<Operator*, Flag> opFlags;
};
}
......
......@@ -60,9 +60,9 @@ namespace AMDiS {
}
MeshDistributor::MeshDistributor(string str)
MeshDistributor::MeshDistributor()
: probStat(0),
name(str),
name("parallel"),
feSpace(NULL),
mesh(NULL),
refineManager(NULL),
......
......@@ -72,13 +72,13 @@ namespace AMDiS {
typedef map<const DegreeOfFreedom*, DegreeOfFreedom> DofIndexMap;
/// Mapps a boundar type, i.e., a boundary identifier index, to a periodic
/// dof mapping.
/// DOF mapping.
typedef map<BoundaryType, DofMapping> PeriodicDofMap;
typedef vector<MeshStructure> MeshCodeVec;
public:
MeshDistributor(string str);
MeshDistributor();
virtual ~MeshDistributor() {}
......@@ -106,7 +106,7 @@ namespace AMDiS {
* This function checks if the mesh has changed on at least on rank. In this case,
* the interior boundaries are adapted on all ranks such that they fit together on
* all ranks. Furthermore the function \ref updateLocalGlobalNumbering() is called
* to update the dof numberings and mappings on all rank due to the new mesh
* to update the DOF numberings and mappings on all rank due to the new mesh
* structure.
*
* \param[in] tryRepartition If this parameter is true, repartitioning may be
......@@ -160,7 +160,7 @@ namespace AMDiS {
return mapLocalGlobalDofs;
}
/// Maps a local dof to its global index.
/// Maps a local DOF to its global index.
inline DegreeOfFreedom mapLocalToGlobal(DegreeOfFreedom dof)
{
return mapLocalGlobalDofs[dof];
......@@ -168,7 +168,7 @@ namespace AMDiS {
DegreeOfFreedom mapGlobalToLocal(DegreeOfFreedom dof);
/// Maps a local dof to its local index.
/// Maps a local DOF to its local index.
inline DegreeOfFreedom mapLocalToDofIndex(DegreeOfFreedom dof)
{
return mapLocalDofIndex[dof];
......
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.
#include <vector>
#include <set>
#include "parallel/PetscProblemStat.h"
#include "parallel/StdMpi.h"
#include "parallel/ParallelDebug.h"
#include "parallel/PetscSolver.h"
#include "DOFVector.h"
#include "Debug.h"
#include "SystemVector.h"
#include "petscksp.h"
namespace AMDiS {
using namespace std;
void PetscProblemStat::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{
FUNCNAME("PetscProblemStat::solve()");
TEST_EXIT(meshDistributor)("Should not happen!\n");
double wtime = MPI::Wtime();
fillPetscMatrix(systemMatrix, rhs);
solvePetscMatrix(*solution, adaptInfo);
INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
MPI::Wtime() - wtime);
}
void PetscProblemStat::setDofMatrix(DOFMatrix* mat, int dispMult,
int dispAddRow, int dispAddCol)
{
FUNCNAME("PetscProblemStat::setDofMatrix()");
TEST_EXIT(mat)("No DOFMatrix!\n");
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits= mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
traits::col<Matrix>::type col(mat->getBaseMatrix());
traits::const_value<Matrix>::type value(mat->getBaseMatrix());
typedef traits::range_generator<row, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
vector<int> cols;
vector<double> values;
cols.reserve(300);
values.reserve(300);
vector<int> globalCols;
// === Traverse all rows of the dof matrix and insert row wise the values ===
// === to the PETSc matrix. ===
for (cursor_type cursor = begin<row>(mat->getBaseMatrix()),
cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {
// Global index of the current row dof.
int globalRowDof = meshDistributor->mapLocalToGlobal(*cursor);
// Test if the current row dof is a periodic dof.
bool periodicRow = meshDistributor->isPeriodicDof(globalRowDof);
if (!periodicRow) {
// === Row DOF index is not periodic. ===
// Calculate PETSc row index.
int rowIndex = globalRowDof * dispMult + dispAddRow;
cols.clear();
values.clear();
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
// Global index of the current column index.
int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
// Test if the current col dof is a periodic dof.
bool periodicCol = meshDistributor->isPeriodicDof(globalColDof);
// Calculate PETSc col index.
int colIndex = globalColDof * dispMult + dispAddCol;
// Ignore all zero entries, expect it is a diagonal entry.
if (value(*icursor) == 0.0 && rowIndex != colIndex)
continue;
if (!periodicCol) {
// Calculate the exact position of the column index in the PETSc matrix.
cols.push_back(globalColDof * dispMult + dispAddCol);
values.push_back(value(*icursor));
} else {
// === Row index is not periodic, but column index is. ===
// Create set of all periodic associations of the column index.
std::set<int> perAsc;
std::set<int>& perColAsc =
meshDistributor->getPerDofAssociations(globalColDof);
for (std::set<int>::iterator it = perColAsc.begin();
it != perColAsc.end(); ++it)
if (*it >= -3)
perAsc.insert(*it);
// Scale value to the number of periodic associations of the column index.
double scaledValue =
value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));
// === Create set of all matrix column indices due to the periodic ===
// === associations of the column DOF index. ===
vector<int> newCols;
// First, add the original matrix index.
newCols.push_back(globalColDof);
// And add all periodic matrix indices.
for (std::set<int>::iterator it = perAsc.begin();
it != perAsc.end(); ++it) {
int nCols = static_cast<int>(newCols.size());
for (int i = 0; i < nCols; i++) {
TEST_EXIT_DBG(meshDistributor->isPeriodicDof(newCols[i], *it))
("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
*it, newCols[i]);
newCols.push_back(meshDistributor->getPeriodicMapping(newCols[i], *it));
}
}
for (unsigned int i = 0; i < newCols.size(); i++) {
cols.push_back(newCols[i] * dispMult + dispAddCol);
values.push_back(scaledValue);
}
}
}
MatSetValues(petscMatrix, 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
} else {
// === Row DOF index is periodic. ===
// Because this row is periodic, we will have to add the entries of this
// matrix row to multiple rows. The following maps store to each row an
// array of column indices and values of the entries that must be added to
// the PETSc matrix.
map<int, vector<int> > colsMap;
map<int, vector<double> > valsMap;
// Traverse all column entries.
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
// Global index of the current column index.
int globalColDof = meshDistributor->mapLocalToGlobal(col(*icursor));
// Ignore all zero entries, expect it is a diagonal entry.
if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
continue;
// === Add all periodic associations of both, the row and the column ===
// === indices to the set perAsc. ===
std::set<int> perAsc;
if (meshDistributor->isPeriodicDof(globalColDof)) {
std::set<int>& perColAsc =
meshDistributor->getPerDofAssociations(globalColDof);
for (std::set<int>::iterator it = perColAsc.begin();
it != perColAsc.end(); ++it)
if (*it >= -3)
perAsc.insert(*it);
}
std::set<int>& perRowAsc =
meshDistributor->getPerDofAssociations(globalRowDof);
for (std::set<int>::iterator it = perRowAsc.begin();
it != perRowAsc.end(); ++it)
if (*it >= -3)
perAsc.insert(*it);
// Scale the value with respect to the number of periodic associations.
double scaledValue =
value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));
// === Create all matrix entries with respect to the periodic ===
// === associations of the row and column indices. ===
vector<pair<int, int> > entry;
// First, add the original entry.
entry.push_back(make_pair(globalRowDof, globalColDof));
// Then, traverse the periodic associations of the row and column indices
// and create the corresponding entries.
for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
int nEntry = static_cast<int>(entry.size());
for (int i = 0; i < nEntry; i++) {
int perRowDof = 0;
if (meshDistributor->getPeriodicMapping()[*it].count(entry[i].first))
perRowDof = meshDistributor->getPeriodicMapping(entry[i].first, *it);