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

Fixed bug, added some documentation.

parent 6cfe12be
...@@ -56,10 +56,9 @@ namespace AMDiS { ...@@ -56,10 +56,9 @@ namespace AMDiS {
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) if (el != lastMatEl || !operat->isOptimized()) {
initElement(elInfo); initElement(elInfo);
if (el != lastMatEl || !operat->isOptimized()) {
if (rememberElMat) if (rememberElMat)
set_to_zero(elementMatrix); set_to_zero(elementMatrix);
......
...@@ -211,7 +211,7 @@ namespace AMDiS { ...@@ -211,7 +211,7 @@ namespace AMDiS {
using namespace mtl; using namespace mtl;
#if 0 #if 0
std::cout << "----- PRINT MAT --------" << std::endl; std::cout << "----- PRINT MAT " << rowElInfo->getElement()->getIndex() << "--------" << std::endl;
std::cout << elMat << std::endl; std::cout << elMat << std::endl;
std::cout << "rows: "; std::cout << "rows: ";
for (int i = 0; i < rowIndices.size(); i++) for (int i = 0; i < rowIndices.size(); i++)
......
...@@ -560,6 +560,8 @@ namespace AMDiS { ...@@ -560,6 +560,8 @@ namespace AMDiS {
adaptInfo->setMaxSolverIterations(solver->getMaxIterations()); adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
adaptInfo->setSolverTolerance(solver->getTolerance()); adaptInfo->setSolverTolerance(solver->getTolerance());
adaptInfo->setSolverResidual(solver->getResidual()); adaptInfo->setSolverResidual(solver->getResidual());
debug::writeMatlabVector(*solution, "sol.dat");
} }
...@@ -828,6 +830,9 @@ namespace AMDiS { ...@@ -828,6 +830,9 @@ namespace AMDiS {
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first, clock())); TIME_USED(first, clock()));
#endif #endif
debug::writeMatlabMatrix(*systemMatrix, "matrix.dat");
debug::writeMatlabVector(*rhs, "rhs.dat");
} }
......
...@@ -28,10 +28,12 @@ namespace AMDiS { ...@@ -28,10 +28,12 @@ namespace AMDiS {
TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh()) TEST_EXIT_DBG(elInfo->getMesh() == vec->getFeSpace()->getMesh())
("Vector and FE space do not fit together!\n"); ("Vector and FE space do not fit together!\n");
Quadrature *localQuad = quad ? quad : quadrature; Quadrature *localQuad = quad ? quad : quadrature;
if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) { if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) {
vecAtQPs = valuesAtQPs[vec]->values; vecAtQPs = valuesAtQPs[vec]->values;
return; return;
} }
......
...@@ -74,12 +74,12 @@ namespace AMDiS { ...@@ -74,12 +74,12 @@ namespace AMDiS {
// create new assembler // create new assembler
if (!optimized) { if (!optimized) {
newAssembler = new StandardZOA(op, assembler, quad); newAssembler = new StandardZOA(op, assembler, quad);
} else { } else {
if (pwConst) if (pwConst)
newAssembler = new PrecalcZOA(op, assembler, quad); newAssembler = new PrecalcZOA(op, assembler, quad);
else else
newAssembler = new FastQuadZOA(op, assembler, quad); newAssembler = new FastQuadZOA(op, assembler, quad);
} }
subAssemblers->push_back(newAssembler); subAssemblers->push_back(newAssembler);
...@@ -103,8 +103,8 @@ namespace AMDiS { ...@@ -103,8 +103,8 @@ namespace AMDiS {
ElementVector c(nPoints, 0.0); ElementVector c(nPoints, 0.0);
std::vector<double> phival(nCol); std::vector<double> phival(nCol);
std::vector<OperatorTerm*>::iterator termIt; for (std::vector<OperatorTerm*>::iterator termIt = terms.begin();
for (termIt = terms.begin(); termIt != terms.end(); ++termIt) termIt != terms.end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c); (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
if (symmetric) { if (symmetric) {
......
...@@ -44,7 +44,7 @@ namespace AMDiS { ...@@ -44,7 +44,7 @@ namespace AMDiS {
} }
void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) void VecAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector& C)
{ {
if (f) { if (f) {
for (int iq = 0; iq < nPoints; iq++) for (int iq = 0; iq < nPoints; iq++)
......
...@@ -571,7 +571,7 @@ namespace AMDiS { ...@@ -571,7 +571,7 @@ namespace AMDiS {
("No edge %d/%d in elObjDB!\n", it->first, it->second); ("No edge %d/%d in elObjDB!\n", it->first, it->second);
std::set<int> el0, el1; std::set<int> el0, el1;
std::map<int, int> edgeNoInEl; map<int, int> edgeNoInEl;
for (unsigned int i = 0; i < edgeEls.size(); i++) { for (unsigned int i = 0; i < edgeEls.size(); i++) {
if (rankMacroEls.count(edgeEls[i].elIndex)) { if (rankMacroEls.count(edgeEls[i].elIndex)) {
...@@ -1197,7 +1197,7 @@ namespace AMDiS { ...@@ -1197,7 +1197,7 @@ namespace AMDiS {
StdMpi<MeshCodeVec> stdMpi(mpiComm, true); StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
stdMpi.send(sendCodes); stdMpi.send(sendCodes);
for (std::map<int, std::vector<int> >::iterator it = partitioner->getRecvElements().begin(); for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin();
it != partitioner->getRecvElements().end(); ++it) it != partitioner->getRecvElements().end(); ++it)
stdMpi.recv(it->first); stdMpi.recv(it->first);
stdMpi.startCommunication(); stdMpi.startCommunication();
...@@ -1208,7 +1208,7 @@ namespace AMDiS { ...@@ -1208,7 +1208,7 @@ namespace AMDiS {
StdMpi<vector<vector<double> > > stdMpi2(mpiComm, true); StdMpi<vector<vector<double> > > stdMpi2(mpiComm, true);
stdMpi2.send(sendValues); stdMpi2.send(sendValues);
for (std::map<int, std::vector<int> >::iterator it = partitioner->getRecvElements().begin(); for (map<int, vector<int> >::iterator it = partitioner->getRecvElements().begin();
it != partitioner->getRecvElements().end(); ++it) it != partitioner->getRecvElements().end(); ++it)
stdMpi2.recv(it->first); stdMpi2.recv(it->first);
stdMpi2.startCommunication(); stdMpi2.startCommunication();
...@@ -2044,7 +2044,7 @@ namespace AMDiS { ...@@ -2044,7 +2044,7 @@ namespace AMDiS {
stdMpi2.startCommunication(); stdMpi2.startCommunication();
for (std::map<int, PeriodicDofMap>::iterator it = stdMpi2.getRecvData().begin(); for (map<int, PeriodicDofMap>::iterator it = stdMpi2.getRecvData().begin();
it != stdMpi2.getRecvData().end(); ++it) { it != stdMpi2.getRecvData().end(); ++it) {
for (PeriodicDofMap::iterator perIt = it->second.begin(); for (PeriodicDofMap::iterator perIt = it->second.begin();
perIt != it->second.end(); ++perIt) { perIt != it->second.end(); ++perIt) {
...@@ -2070,9 +2070,11 @@ namespace AMDiS { ...@@ -2070,9 +2070,11 @@ namespace AMDiS {
allMacroElements.push_back(*it); allMacroElements.push_back(*it);
macroElementNeighbours[(*it)->getIndex()].resize(mesh->getGeo(NEIGH)); int elIndex = (*it)->getIndex();
macroElementNeighbours[elIndex].resize(mesh->getGeo(NEIGH));
for (int i = 0; i < mesh->getGeo(NEIGH); i++) for (int i = 0; i < mesh->getGeo(NEIGH); i++)
macroElementNeighbours[(*it)->getIndex()][i] = macroElementNeighbours[elIndex][i] =
(*it)->getNeighbour(i) ? (*it)->getNeighbour(i)->getIndex() : -1; (*it)->getNeighbour(i) ? (*it)->getNeighbour(i)->getIndex() : -1;
} }
} }
......
...@@ -62,8 +62,9 @@ namespace AMDiS { ...@@ -62,8 +62,9 @@ namespace AMDiS {
void exitParallelization(); void exitParallelization();
/// Adds a DOFVector to the set of \ref interchangeVecs. Thus, this vector will /// Adds a DOFVector to the set of \ref interchangeVecs. Thus, this vector
/// be automatically interchanged between ranks when mesh is repartitioned. /// will be automatically interchanged between ranks when mesh is
/// repartitioned.
void addInterchangeVector(DOFVector<double> *vec) void addInterchangeVector(DOFVector<double> *vec)
{ {
interchangeVectors.push_back(vec); interchangeVectors.push_back(vec);
...@@ -77,28 +78,30 @@ namespace AMDiS { ...@@ -77,28 +78,30 @@ namespace AMDiS {
} }
/** \brief /** \brief
* This function checks if the mesh has changed on at least on rank. In this case, * This function checks if the mesh has changed on at least on rank. In
* the interior boundaries are adapted on all ranks such that they fit together on * this case, the interior boundaries are adapted on all ranks such that
* all ranks. Furthermore the function \ref updateLocalGlobalNumbering() is called * they fit together on all ranks. Furthermore the function
* to update the DOF numberings and mappings on all rank due to the new mesh * \ref updateLocalGlobalNumbering() is called to update the DOF numberings
* structure. * and mappings on all rank due to the new mesh structure.
* *
* \param[in] tryRepartition If this parameter is true, repartitioning may be * \param[in] tryRepartition If this parameter is true, repartitioning
* done. This depends on several other parameters. If * may be done. This depends on several other
* the parameter is false, the mesh is only checked * parameters. If the parameter is false, the
* and adapted but never repartitioned. * mesh is only checked and adapted but never
* repartitioned.
*/ */
void checkMeshChange(bool tryRepartition = true); void checkMeshChange(bool tryRepartition = true);
/** \brief /** \brief
* Test, if the mesh consists of macro elements only. The mesh partitioning of * Test, if the mesh consists of macro elements only. The mesh partitioning
* the parallelization works for macro meshes only and would fail, if the mesh * of the parallelization works for macro meshes only and would fail, if the
* is already refined in some way. Therefore, this function will exit the program * mesh is already refined in some way. Therefore, this function will exit
* if it finds a non macro element in the mesh. * the program if it finds a non macro element in the mesh.
*/ */
void testForMacroMesh(); void testForMacroMesh();
/// Set for each element on the partitioning level the number of leaf elements. /// Set for each element on the partitioning level the number of
/// leaf elements.
void setInitialElementWeights(); void setInitialElementWeights();
inline virtual string getName() inline virtual string getName()
...@@ -154,7 +157,8 @@ namespace AMDiS { ...@@ -154,7 +157,8 @@ namespace AMDiS {
return periodicDof; return periodicDof;
} }
/// Returns for a global dof index its periodic mapping for a given boundary type. /// Returns for a global dof index its periodic mapping for a given
/// boundary type.
inline int getPeriodicMapping(int globalDofIndex, BoundaryType type) inline int getPeriodicMapping(int globalDofIndex, BoundaryType type)
{ {
FUNCNAME("MeshDistributor::getPeriodicMapping()"); FUNCNAME("MeshDistributor::getPeriodicMapping()");
...@@ -167,7 +171,8 @@ namespace AMDiS { ...@@ -167,7 +171,8 @@ namespace AMDiS {
} }
/// For a given global DOF index, this function returns the set of periodic /// For a given global DOF index, this function returns the set of periodic
/// associations, i.e., the boundary types the DOF is associated to, for this DOF. /// associations, i.e., the boundary types the DOF is associated to, for
/// this DOF.
inline std::set<BoundaryType>& getPerDofAssociations(int globalDofIndex) inline std::set<BoundaryType>& getPerDofAssociations(int globalDofIndex)
{ {
TEST_EXIT_DBG(periodicDofAssociations.count(globalDofIndex)) TEST_EXIT_DBG(periodicDofAssociations.count(globalDofIndex))
...@@ -190,8 +195,8 @@ namespace AMDiS { ...@@ -190,8 +195,8 @@ namespace AMDiS {
return (periodicDof[type].count(globalDofIndex) > 0); return (periodicDof[type].count(globalDofIndex) > 0);
} }
/// Return true, if the given DOF is owned by the rank. If false, the DOF is in /// Return true, if the given DOF is owned by the rank. If false, the DOF
/// rank's partition, but is owned by some other rank. /// is in rank's partition, but is owned by some other rank.
inline bool getIsRankDof(DegreeOfFreedom dof) inline bool getIsRankDof(DegreeOfFreedom dof)
{ {
if (isRankDof.count(dof)) if (isRankDof.count(dof))
...@@ -251,14 +256,14 @@ namespace AMDiS { ...@@ -251,14 +256,14 @@ namespace AMDiS {
void deserialize(istream &in); void deserialize(istream &in);
/** \brief /** \brief
* This function must be used if the values of a DOFVector must be synchronised * This function must be used if the values of a DOFVector must be
* over all ranks. That means, that each rank sends the values of the DOFs, which * synchronised over all ranks. That means, that each rank sends the
* are owned by the rank and lie on an interior bounday, to all other ranks also * values of the DOFs, which are owned by the rank and lie on an interior
* having these DOFs. * bounday, to all other ranks also having these DOFs.
* *
* This function must be used, for example, after the lineary system is solved, or * This function must be used, for example, after the lineary system is
* after the DOFVector is set by some user defined functions, e.g., initial * solved, or after the DOFVector is set by some user defined functions,
* solution functions. * e.g., initial solution functions.
*/ */
template<typename T> template<typename T>
void synchVector(DOFVector<T> &vec) void synchVector(DOFVector<T> &vec)
...@@ -290,9 +295,9 @@ namespace AMDiS { ...@@ -290,9 +295,9 @@ namespace AMDiS {
} }
/** \brief /** \brief
* Works in the same way as the function above defined for DOFVectors. Due to * Works in the same way as the function above defined for DOFVectors. Due
* performance, this function does not call \ref synchVector for each DOFVector, * to performance, this function does not call \ref synchVector for each
* but instead sends all values of all DOFVectors all at once. * DOFVector, but instead sends all values of all DOFVectors all at once.
*/ */
void synchVector(SystemVector &vec); void synchVector(SystemVector &vec);
...@@ -319,10 +324,8 @@ namespace AMDiS { ...@@ -319,10 +324,8 @@ namespace AMDiS {
protected: protected:
void addProblemStat(ProblemStatSeq *probStat); void addProblemStat(ProblemStatSeq *probStat);
/** \brief /// Determines the interior boundaries, i.e. boundaries between ranks, and
* Determines the interior boundaries, i.e. boundaries between ranks, and stores /// stores all information about them in \ref interiorBoundary.
* all information about them in \ref interiorBoundary.
*/
void createInteriorBoundaryInfo(); void createInteriorBoundaryInfo();
void updateInteriorBoundaryInfo(); void updateInteriorBoundaryInfo();
...@@ -333,36 +336,47 @@ namespace AMDiS { ...@@ -333,36 +336,47 @@ namespace AMDiS {
void createBoundaryDofs(); void createBoundaryDofs();
/// Removes all macro elements from the mesh that are not part of ranks partition. /// Removes all macro elements from the mesh that are not part of ranks
/// partition.
void removeMacroElements(); void removeMacroElements();
/// Updates the local and global DOF numbering after the mesh has been changed. /// Updates the local and global DOF numbering after the mesh has been
/// changed.
void updateLocalGlobalNumbering(); void updateLocalGlobalNumbering();
/** \brief /** \brief
* Creates to all dofs in rank's partition that are on a periodic boundary the * Creates to all dofs in rank's partition that are on a periodic boundary
* mapping from dof index to the other periodic dof indices. This information * the mapping from dof index to the other periodic dof indices. This
* is stored in \ref periodicDof. * information is stored in \ref periodicDof.
*/ */
void createPeriodicMap(); void createPeriodicMap();
/** \brief
* This function is called only once during the initialization when the
* whole macro mesh is available on all cores. It copies the pointers of all
* macro elements to \ref allMacroElements and stores all neighbour
* information based on macro element indices (and not pointer based) in
* \ref macroElementNeighbours. These information are then used to
* reconstruct macro elements during mesh redistribution.
*/
void createMacroElementInfo(); void createMacroElementInfo();
void updateMacroElementInfo(); void updateMacroElementInfo();
/** \brief /** \brief
* Checks for all given interior boundaries if the elements fit together on both * Checks for all given interior boundaries if the elements fit together on
* sides of the boundaries. If this is not the case, the mesh is adapted. Because * both sides of the boundaries. If this is not the case, the mesh is
* refinement of a certain element may forces the refinement of other elements, * adapted. Because refinement of a certain element may forces the
* it is not guaranteed that all rank's meshes fit together after this function * refinement of other elements, it is not guaranteed that all rank's meshes
* terminates. Hence, it must be called until a stable mesh refinement is reached. * fit together after this function terminates. Hence, it must be called
* until a stable mesh refinement is reached.
* *
* \param[in] allBound Defines a map from rank to interior boundaries which * \param[in] allBound Defines a map from rank to interior boundaries
* should be checked. * which should be checked.
* *
* \return If the mesh has been changed by this function, it returns true. * \return If the mesh has been changed by this function, it returns
* Otherwise, it returns false, i.e., the given interior boundaries * true. Otherwise, it returns false, i.e., the given interior
* fit together on both sides. * boundaries fit together on both sides.
*/ */
bool checkAndAdaptBoundary(RankToBoundMap &allBound); bool checkAndAdaptBoundary(RankToBoundMap &allBound);
...@@ -377,7 +391,8 @@ namespace AMDiS { ...@@ -377,7 +391,8 @@ namespace AMDiS {
/// stationary problem. /// stationary problem.
void setRankDofs(ProblemStatSeq *probStat); void setRankDofs(ProblemStatSeq *probStat);
/// Sets \ref isRankDof to all matrices and rhs vectors in all stationary problems. /// Sets \ref isRankDof to all matrices and rhs vectors in all
/// stationary problems.
void setRankDofs(); void setRankDofs();
/// Removes all periodic boundary condition information from all matrices and /// Removes all periodic boundary condition information from all matrices and
...@@ -450,7 +465,8 @@ namespace AMDiS { ...@@ -450,7 +465,8 @@ namespace AMDiS {
} }
protected: protected:
/// List of all stationary problems that are managed by this mesh distributor. /// List of all stationary problems that are managed by this mesh
/// distributor.
vector<ProblemStatSeq*> problemStat; vector<ProblemStatSeq*> problemStat;
/// If true, the mesh distributor is already initialized; /// If true, the mesh distributor is already initialized;
...@@ -462,11 +478,8 @@ namespace AMDiS { ...@@ -462,11 +478,8 @@ namespace AMDiS {
/// Overall number of processes. /// Overall number of processes.
int mpiSize; int mpiSize;
/** \brief /// MPI communicator collected all processes, which should be used for
* MPI communicator collected all processes, which should /// calculation. The Debug procces is not included in this communicator.
* be used for calculation. The Debug procces is not included
* in this communicator.
*/
MPI::Intracomm mpiComm; MPI::Intracomm mpiComm;
/// Name of the problem (as used in the init files) /// Name of the problem (as used in the init files)
...@@ -479,19 +492,21 @@ namespace AMDiS { ...@@ -479,19 +492,21 @@ namespace AMDiS {
Mesh *mesh; Mesh *mesh;
/** \brief /** \brief
* A refinement manager that should be used on the mesh. It is used to refine * A refinement manager that should be used on the mesh. It is used to
* elements at interior boundaries in order to fit together with elements on the * refine elements at interior boundaries in order to fit together with
* other side of the interior boundary. * elements on the other side of the interior boundary.
*/ */
RefinementManager *refineManager; RefinementManager *refineManager;
/// Info level. /// Info level.
int info; int info;
/// Pointer to a mesh partitioner that is used to partition the mesh to the ranks. /// Pointer to a mesh partitioner that is used to partition the mesh to
/// the ranks.
MeshPartitioner *partitioner; MeshPartitioner *partitioner;
/// Weights for the elements, i.e., the number of leaf elements within this element. /// Weights for the elements, i.e., the number of leaf elements within
/// this element.
map<int, double> elemWeights; map<int, double> elemWeights;
/** \brief /** \brief
...@@ -506,47 +521,48 @@ namespace AMDiS { ...@@ -506,47 +521,48 @@ namespace AMDiS {
/// Number of DOFs in the whole domain. /// Number of DOFs in the whole domain.
int nOverallDofs; int nOverallDofs;
// Data structure to store all sub-objects of all elements of the macro mesh. /// Data structure to store all sub-objects of all elements of the
/// macro mesh.
ElementObjects elObjects; ElementObjects elObjects;
// Maps to each macro element index a pointer to the corresponding element. /// Maps to each macro element index a pointer to the corresponding element.
map<int, Element*> macroElIndexMap; map<int, Element*> macroElIndexMap;
// Maps to each macro element index the type of this element. /// Maps to each macro element index the type of this element.
map<int, int> macroElIndexTypeMap; map<int, int> macroElIndexTypeMap;
/** \brief /** \brief
* Defines the interior boundaries of the domain that result from partitioning * Defines the interior boundaries of the domain that result from
* the whole mesh. Contains only the boundaries, which are owned by the rank, i.e., * partitioning the whole mesh. Contains only the boundaries, which are
* the object gives for every neighbour rank i the boundaries this rank owns and * owned by the rank, i.e., the object gives for every neighbour rank i
* shares with rank i. * the boundaries this rank owns and shares with rank i.
*/ */
InteriorBoundary myIntBoundary; InteriorBoundary myIntBoundary;
/** \brief /** \brief
* Defines the interior boundaries of the domain that result from partitioning * Defines the interior boundaries of the domain that result from
* the whole mesh. Contains only the boundaries, which are not owned by the rank, * partitioning the whole mesh. Contains only the boundaries, which are
* i.e., the object gives for every neighbour rank i the boundaries that are * not owned by the rank, i.e., the object gives for every neighbour rank
* owned by rank i and are shared with this rank. * i the boundaries that are owned by rank i and are shared with this rank.
*/ */
InteriorBoundary otherIntBoundary; InteriorBoundary otherIntBoundary;
/** \brief /** \brief
* Defines the periodic boundaries with other ranks. Periodic boundaries have * Defines the periodic boundaries with other ranks. Periodic boundaries
* no owner, as it is the case of interior boundaries. * have no owner, as it is the case of interior boundaries.
*/ */
InteriorBoundary periodicBoundary; InteriorBoundary periodicBoundary;
/** \brief /** \brief
* This map contains for each rank the list of DOFs the current rank must send * This map contains for each rank the list of DOFs the current rank must
* to exchange solution DOFs at the interior boundaries. * send to exchange solution DOFs at the interior boundaries.
*/ */
RankToDofContainer sendDofs; RankToDofContainer sendDofs;