Skip to content

GitLab

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

Moved debugging functions in ProblemVec to a new class ProblemVecDbg.

parent c7167424
Hide whitespace changes
Inline Side-by-side
AMDiS/bin/Makefile.am
View file @ ded81b7f
......@@ -82,6 +82,7 @@ $(SOURCE_DIR)/ProblemIterationInterface.h \
$(SOURCE_DIR)/StandardProblemIteration.h $(SOURCE_DIR)/StandardProblemIteration.cc \
$(SOURCE_DIR)/ProblemScal.h $(SOURCE_DIR)/ProblemScal.cc \
$(SOURCE_DIR)/ProblemVec.h $(SOURCE_DIR)/ProblemVec.cc \
$(SOURCE_DIR)/ProblemVecDbg.h $(SOURCE_DIR)/ProblemVecDbg.cc \
$(SOURCE_DIR)/DualTraverse.h $(SOURCE_DIR)/DualTraverse.cc \
$(SOURCE_DIR)/ElementPartition_ED.h $(SOURCE_DIR)/SurfacePartition_ED.h \
$(SOURCE_DIR)/ElementData.h $(SOURCE_DIR)/ElementData.cc \
......
AMDiS/bin/Makefile.in
View file @ ded81b7f
......@@ -109,6 +109,7 @@ am__libamdis_la_SOURCES_DIST = $(PARALLEL_DIR)/ParallelDomainBase.h \
$(SOURCE_DIR)/StandardProblemIteration.cc \
$(SOURCE_DIR)/ProblemScal.h $(SOURCE_DIR)/ProblemScal.cc \
$(SOURCE_DIR)/ProblemVec.h $(SOURCE_DIR)/ProblemVec.cc \
$(SOURCE_DIR)/ProblemVecDbg.h $(SOURCE_DIR)/ProblemVecDbg.cc \
$(SOURCE_DIR)/DualTraverse.h $(SOURCE_DIR)/DualTraverse.cc \
$(SOURCE_DIR)/ElementPartition_ED.h \
$(SOURCE_DIR)/SurfacePartition_ED.h \
......@@ -248,7 +249,8 @@ am_libamdis_la_OBJECTS = $(am__objects_2) libamdis_la-DOFIndexed.lo \
libamdis_la-AdaptBase.lo \
libamdis_la-StandardProblemIteration.lo \
libamdis_la-ProblemScal.lo libamdis_la-ProblemVec.lo \
libamdis_la-DualTraverse.lo libamdis_la-ElementData.lo \
libamdis_la-ProblemVecDbg.lo libamdis_la-DualTraverse.lo \
libamdis_la-ElementData.lo \
libamdis_la-ComponentTraverseInfo.lo libamdis_la-CreatorMap.lo \
libamdis_la-ProblemInterpolScal.lo \
libamdis_la-ProblemInterpolVec.lo libamdis_la-MacroReader.lo \
......@@ -493,6 +495,7 @@ $(SOURCE_DIR)/ProblemIterationInterface.h \
$(SOURCE_DIR)/StandardProblemIteration.h $(SOURCE_DIR)/StandardProblemIteration.cc \
$(SOURCE_DIR)/ProblemScal.h $(SOURCE_DIR)/ProblemScal.cc \
$(SOURCE_DIR)/ProblemVec.h $(SOURCE_DIR)/ProblemVec.cc \
$(SOURCE_DIR)/ProblemVecDbg.h $(SOURCE_DIR)/ProblemVecDbg.cc \
$(SOURCE_DIR)/DualTraverse.h $(SOURCE_DIR)/DualTraverse.cc \
$(SOURCE_DIR)/ElementPartition_ED.h $(SOURCE_DIR)/SurfacePartition_ED.h \
$(SOURCE_DIR)/ElementData.h $(SOURCE_DIR)/ElementData.cc \
......@@ -777,6 +780,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemNonLin.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemScal.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemVec.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-ProblemVecDbg.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Projection.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-QPsiPhi.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libamdis_la-Quadrature.Plo@am__quote@
......@@ -969,6 +973,13 @@ libamdis_la-ProblemVec.lo: $(SOURCE_DIR)/ProblemVec.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-ProblemVec.lo `test -f '$(SOURCE_DIR)/ProblemVec.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/ProblemVec.cc
libamdis_la-ProblemVecDbg.lo: $(SOURCE_DIR)/ProblemVecDbg.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-ProblemVecDbg.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-ProblemVecDbg.Tpo" -c -o libamdis_la-ProblemVecDbg.lo `test -f '$(SOURCE_DIR)/ProblemVecDbg.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/ProblemVecDbg.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-ProblemVecDbg.Tpo" "$(DEPDIR)/libamdis_la-ProblemVecDbg.Plo"; else rm -f "$(DEPDIR)/libamdis_la-ProblemVecDbg.Tpo"; exit 1; fi
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ source='$(SOURCE_DIR)/ProblemVecDbg.cc' object='libamdis_la-ProblemVecDbg.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-ProblemVecDbg.lo `test -f '$(SOURCE_DIR)/ProblemVecDbg.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/ProblemVecDbg.cc
libamdis_la-DualTraverse.lo: $(SOURCE_DIR)/DualTraverse.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-DualTraverse.lo -MD -MP -MF "$(DEPDIR)/libamdis_la-DualTraverse.Tpo" -c -o libamdis_la-DualTraverse.lo `test -f '$(SOURCE_DIR)/DualTraverse.cc' || echo '$(srcdir)/'`$(SOURCE_DIR)/DualTraverse.cc; \
@am__fastdepCXX_TRUE@ then mv -f "$(DEPDIR)/libamdis_la-DualTraverse.Tpo" "$(DEPDIR)/libamdis_la-DualTraverse.Plo"; else rm -f "$(DEPDIR)/libamdis_la-DualTraverse.Tpo"; exit 1; fi
......
AMDiS/libtool
View file @ ded81b7f
......@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
# Libtool was configured on host deimos104:
# Libtool was configured on host p2d125:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
# Libtool was configured on host deimos104:
# Libtool was configured on host p2d125:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
# Libtool was configured on host deimos104:
# Libtool was configured on host p2d125:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......
AMDiS/src/ParallelDomainBase.cc
View file @ ded81b7f
......@@ -16,7 +16,6 @@
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
#include "ElementFileWriter.h"
#include "Serializer.h"
#include "petscksp.h"
......@@ -35,13 +34,13 @@ namespace AMDiS {
return (*dof1 < *dof2);
}
ParallelDomainBase::ParallelDomainBase(std::string name,
ProblemIterationInterface *iIF,
ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF,
ProblemTimeInterface *tIF,
FiniteElemSpace *fe,
RefinementManager *refineManager)
: iterationIF(iIF),
timeIF(tIF),
name(iIF->getName()),
feSpace(fe),
mesh(fe->getMesh()),
refinementManager(refineManager),
......@@ -1605,7 +1604,7 @@ namespace AMDiS {
SerUtil::serialize(out, mapLocalToDofIndex);
SerUtil::serialize(out, isRankDof);
// SerUtil::serialize(out, vertexDof);
serialize(out, vertexDof);
SerUtil::serialize(out, &rstart);
SerUtil::serialize(out, &nRankRows);
......@@ -1649,7 +1648,7 @@ namespace AMDiS {
SerUtil::deserialize(in, mapLocalToDofIndex);
SerUtil::deserialize(in, isRankDof);
// SerUtil::deserialize(in, vertexDof);
deserialize(in, vertexDof, dofMap);
SerUtil::deserialize(in, &rstart);
SerUtil::deserialize(in, &nRankRows);
......
AMDiS/src/ParallelDomainBase.h
View file @ ded81b7f
......@@ -32,6 +32,7 @@
#include "FiniteElemSpace.h"
#include "AdaptInfo.h"
#include "InteriorBoundary.h"
#include "Serializer.h"
#include "AMDiS_fwd.h"
#include "petsc.h"
......@@ -81,8 +82,7 @@ namespace AMDiS {
typedef std::map<const DegreeOfFreedom*, DegreeOfFreedom> DofIndexMap;
public:
ParallelDomainBase(std::string name,
ProblemIterationInterface *iterationIF,
ParallelDomainBase(ProblemIterationInterface *iterationIF,
ProblemTimeInterface *timeIF,
FiniteElemSpace *feSpace,
RefinementManager *refineManager);
......@@ -305,7 +305,38 @@ namespace AMDiS {
/// Reads a \ref RankToDofContainer from an input stream.
void deserialize(std::istream &in, RankToDofContainer &data,
std::map<int, const DegreeOfFreedom*> &dofMap);
/// Writes a mapping from dof pointers to some values to an output stream.
template<typename T>
void serialize(std::ostream &out, std::map<const DegreeOfFreedom*, T> &data)
{
int mapSize = data.size();
SerUtil::serialize(out, &mapSize);
for (typename std::map<const DegreeOfFreedom*, T>::iterator it = data.begin();
it != data.end(); ++it) {
int v1 = (*(it->first));
T v2 = it->second;
SerUtil::serialize(out, &v1);
SerUtil::serialize(out, &v2);
}
}
/// Reads a mapping from dof pointer to some values from an input stream.
template<typename T>
void deserialize(std::istream &in, std::map<const DegreeOfFreedom*, T> &data,
std::map<int, const DegreeOfFreedom*> &dofMap)
{
int mapSize = 0;
SerUtil::deserialize(in, &mapSize);
for (int i = 0; i < mapSize; i++) {
int v1 = 0;
T v2;
SerUtil::deserialize(in, &v1);
SerUtil::deserialize(in, &v2);
data[dofMap[v1]] = v2;
}
}
inline void orderDOFs(const DegreeOfFreedom* dof1,
const DegreeOfFreedom* dof2,
const DegreeOfFreedom* dof3,
......
AMDiS/src/ParallelDomainScal.cc
View file @ ded81b7f
......@@ -5,11 +5,9 @@
namespace AMDiS {
ParallelDomainScal::ParallelDomainScal(std::string name,
ProblemScal *problem,
ParallelDomainScal::ParallelDomainScal(ProblemScal *problem,
ProblemInstatScal *problemInstat)
: ParallelDomainBase(name,
problem,
: ParallelDomainBase(problem,
problemInstat,
problem->getFESpace(),
problem->getRefinementManager()),
......
AMDiS/src/ParallelDomainScal.h
View file @ ded81b7f
......@@ -29,8 +29,7 @@ namespace AMDiS {
class ParallelDomainScal : public ParallelDomainBase
{
public:
ParallelDomainScal(std::string name,
ProblemScal *problem,
ParallelDomainScal(ProblemScal *problem,
ProblemInstatScal *problemInstat);
void initParallelization(AdaptInfo *adaptInfo);
......
AMDiS/src/ParallelDomainVec.cc
View file @ ded81b7f
......@@ -5,11 +5,9 @@
namespace AMDiS {
ParallelDomainVec::ParallelDomainVec(std::string name,
ProblemVec *problem,
ParallelDomainVec::ParallelDomainVec(ProblemVec *problem,
ProblemInstatVec *problemInstat)
: ParallelDomainBase(name,
problem,
: ParallelDomainBase(problem,
problemInstat,
problem->getFESpace(0),
problem->getRefinementManager(0)),
......
AMDiS/src/ParallelDomainVec.h
View file @ ded81b7f
......@@ -30,8 +30,7 @@ namespace AMDiS {
class ParallelDomainVec : public ParallelDomainBase
{
public:
ParallelDomainVec(std::string name,
ProblemVec *problem,
ParallelDomainVec(ProblemVec *problem,
ProblemInstatVec *problemInstat);
void initParallelization(AdaptInfo *adaptInfo);
......
AMDiS/src/ProblemScal.cc
View file @ ded81b7f
......@@ -485,10 +485,14 @@ namespace AMDiS {
{
fileWriters.push_back(new FileWriter(name + "->output", mesh, solution));
int writeSerialization = 0;
GET_PARAMETER(0, name + "->output->write serialization", "%d",
&writeSerialization);
GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
// Serialization is not allowed to be done by the problem, if its part of a parallel
// problem definition. Than, the parallel problem object must be serialized.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
if (writeSerialization)
fileWriters.push_back(new Serializer<ProblemScal>(this));
#endif
}
void ProblemScal::estimate(AdaptInfo *adaptInfo)
......
AMDiS/src/ProblemVec.cc
View file @ ded81b7f
......@@ -470,8 +470,13 @@ namespace AMDiS {
}
GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
// Serialization is not allowed to be done by the problem, if its part of a parallel
// problem definition. Than, the parallel problem object must be serialized.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
if (writeSerialization)
fileWriters.push_back(new Serializer<ProblemVec>(this));
#endif
}
void ProblemVec::doOtherStuff()
......@@ -1385,404 +1390,6 @@ namespace AMDiS {
delete sol;
}
}
void ProblemVec::createDofToCoordMap(DofToCoord &dofMap)
{
const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
WorldVector<double> coords;
ElementDofIterator elDofIter(componentSpaces[0], true);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(componentMeshes[0], -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
elDofIter.reset(elInfo->getElement());
int i = 0;
do {
DimVec<double> *baryCoord = basisFcts->getCoords(i);
elInfo->coordToWorld(*baryCoord, coords);
i++;
dofMap[elDofIter.getDof()] = coords;
} while (elDofIter.next());
elInfo = stack.traverseNext(elInfo);
}
}
void ProblemVec::createCoordToDofMap(CoordToDof &dofMap)
{
const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
WorldVector<double> coords;
ElementDofIterator elDofIter(componentSpaces[0], true);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(componentMeshes[0], -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
elDofIter.reset(elInfo->getElement());
int i = 0;
do {
DimVec<double> *baryCoord = basisFcts->getCoords(i);
elInfo->coordToWorld(*baryCoord, coords);
i++;
dofMap[coords] = elDofIter.getDof();
} while (elDofIter.next());
elInfo = stack.traverseNext(elInfo);
}
}
void ProblemVec::writeDbgMatrix(std::string filename)
{
FUNCNAME("ProblemVec::writeDbMatrix()");
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
typedef traits::range_generator<major, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
// Create map with all world coordinates of all DOFs.
DofToCoord dofToCoords;
createDofToCoordMap(dofToCoords);
// Start writing output file.
std::ofstream out(filename.c_str());
out.precision(15);
// First, we write the number of DOFs the file contains.
out << dofToCoords.size() << std::endl;
int dimOfWorld = Global::getGeo(WORLD);
// And now write all the DOF's coords. The first line contains the dof index, all
// the other lines contain one component of the coordinates.
for (DofToCoord::iterator it = dofToCoords.begin(); it != dofToCoords.end(); ++it) {
out << it->first << std::endl;
for (int j = 0; j < dimOfWorld; j++)
out << (it->second)[j] << std::endl;
}
// Write the matrices.
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
if (!dofmatrix)
continue;
Matrix& matrix = dofmatrix->getBaseMatrix();
int nnz = matrix.nnz();
int testNnz = 0;
// Write to file, how many entries the matrix conatins.
out << nnz << std::endl;
traits::row<Matrix>::type row(matrix);
traits::col<Matrix>::type col(matrix);
traits::const_value<Matrix>::type value(matrix);
for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) {
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) {
out << row(*icursor) << " "
<< col(*icursor) << " "
<< value(*icursor) << std::endl;
testNnz++;
}
}
// Just to check, if nnz() is correct.
TEST_EXIT(testNnz == nnz)("NNZ does not fit together!\n");
}
}
for (int i = 0; i < nComponents; i++) {
out << rhs->getDOFVector(i)->getUsedSize() << std::endl;
DOFIterator<double> it(rhs->getDOFVector(i), USED_DOFS);
int c = 0;
for (it.reset(); !it.end(); ++it) {
out << c << " " << *it << std::endl;
c++;
}
TEST_EXIT(c == rhs->getDOFVector(i)->getUsedSize())
("RHS NNZ does not fit together!\n");
}
for (int i = 0; i < nComponents; i++) {
out << solution->getDOFVector(i)->getUsedSize() << std::endl;
DOFIterator<double> it(solution->getDOFVector(i), USED_DOFS);
int c = 0;
for (it.reset(); !it.end(); ++it) {
out << c << " " << *it << std::endl;
c++;
}
TEST_EXIT(c == rhs->getDOFVector(i)->getUsedSize())
("RHS NNZ does not fit together!\n");
}
out.close();
}
void ProblemVec::readAndCompareDbgMatrix(std::vector<std::string> filenames)
{
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
typedef traits::range_generator<major, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
// Create a map from coords of all DOFs, to the DOF indices in this problem.
CoordToDof coordToDof;
createCoordToDofMap(coordToDof);
int dimOfWorld = Global::getGeo(WORLD);
std::vector<std::vector<std::map<std::pair<int, int>, double> > > nnzValues(nComponents);
std::vector<std::map<int, double> > rhsValues(nComponents);
std::vector<std::map<int, double> > solValues(nComponents);
for (int i = 0; i < nComponents; i++)
nnzValues[i].resize(nComponents);
// Stores to each dof index of this problem a map from rank indices (of each rank
// that also has this dof) to the corresponding local dof index.
std::map<DegreeOfFreedom, std::vector<DegreeOfFreedom> > dofMapHereToFiles;
for (std::vector<std::string>::iterator fileIt = filenames.begin();
fileIt != filenames.end(); ++fileIt) {
// Open file and read the number of DOFs the file contains.
std::ifstream in(fileIt->c_str());
int nReadDofs;
in >> nReadDofs;
// Is used to map the dof indices in the files to the global coordinates.
DofToCoord dofToCoord;
// Read all DOF indices and their world coordinates.
for (int i = 0; i < nReadDofs; i++) {
DegreeOfFreedom dof;
WorldVector<double> coords;
in >> dof;
for (int j = 0; j < dimOfWorld; j++)
in >> coords[j];
dofToCoord[dof] = coords;
}
std::map<DegreeOfFreedom, DegreeOfFreedom> dofMapFileToHere;
for (DofToCoord::iterator it = dofToCoord.begin(); it != dofToCoord.end(); ++it) {
DegreeOfFreedom dofIndexInFile = it->first;
WorldVector<double> &dofCoords = it->second;
if (coordToDof.find(dofCoords) == coordToDof.end()) {
std::cout << "Cannot find dof index for coords: " << std::endl;
dofCoords.print();
exit(0);
}
DegreeOfFreedom dofIndexHere = coordToDof[dofCoords];
dofMapHereToFiles[dofIndexHere].push_back(dofIndexInFile);
dofMapFileToHere[dofIndexInFile] = dofIndexHere;
}
// Now we traverse all matrices and check them.
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
if (!dofmatrix)
continue;
int readNnz;
in >> readNnz;
// Read each entry in file and check it with the corresponding matrix value
// in this problem.
for (int k = 0; k < readNnz; k++) {
DegreeOfFreedom row, col;
double value;
in >> row;
in >> col;
in >> value;
if (dofMapFileToHere.count(row) == 0) {
std::cout << "Cannot find row index for: " << row << std::endl;
exit(0);
}
if (dofMapFileToHere.count(col) == 0) {
std::cout << "Cannot find col index for: " << col << std::endl;
exit(0);
}
// Get dof indices for row and col of this problem matrix.
DegreeOfFreedom rowHere = dofMapFileToHere[row];
DegreeOfFreedom colHere = dofMapFileToHere[col];
std::pair<int, int> rowcol = make_pair(rowHere, colHere);
if (nnzValues[i][j].count(rowcol) == 0)
nnzValues[i][j][rowcol] = value;
else
nnzValues[i][j][rowcol] += value;
}
}
}
for (int i = 0; i < nComponents; i++) {
int readNnz;
in >> readNnz;
for (int k = 0; k < readNnz; k++) {
int row;
double value;
in >> row;
in >> value;
WorldVector<double> rowCoords = dofToCoord[row];
if (coordToDof.find(rowCoords) == coordToDof.end()) {
std::cout << "Cannot find row index for coords: " << std::endl;
rowCoords.print();
exit(0);
}
DegreeOfFreedom rowHere = coordToDof[rowCoords];
if (rhsValues[i].count(rowHere) == 0)
rhsValues[i][rowHere] = value;
else
rhsValues[i][rowHere] += value;
}
}
for (int i = 0; i < nComponents; i++) {
int readNnz;
in >> readNnz;
for (int k = 0; k < readNnz; k++) {
int row;
double value;
in >> row;
in >> value;
WorldVector<double> rowCoords = dofToCoord[row];
if (coordToDof.find(rowCoords) == coordToDof.end()) {
std::cout << "Cannot find row index for coords: " << std::endl;
rowCoords.print();
exit(0);
}
DegreeOfFreedom rowHere = coordToDof[rowCoords];
if (solValues[i].count(rowHere) == 0) {
solValues[i][rowHere] = value;
} else {
double diff = fabs(solValues[i][rowHere] - value);
if (diff > 1e-8) {
std::cout << "DIFFERENT values in solution vector!" << std::endl;
exit(0);
}
}
}
}
in.close();
std::cout << "File read!" << std::endl;
}
std::cout << "Now checking ..." << std::endl;
std::cout.precision(15);
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
DOFMatrix *dofmatrix = (*systemMatrix)[i][j];
if (!dofmatrix)
continue;
for (std::map<std::pair<int, int>, double>::iterator nnzIt =
nnzValues[i][j].begin(); nnzIt != nnzValues[i][j].end(); ++nnzIt) {
int row = nnzIt->first.first;
int col = nnzIt->first.second;
double value = nnzIt->second;