From 92f468303351fc2c2c3438d07fbd40c099c38f67 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Tue, 23 Jun 2009 15:16:27 +0000
Subject: [PATCH] Periodic boundary condition work for higher elements.
---
AMDiS/src/ElInfo2d.cc | 99 +++++++++++++++++++++++++++++++---------
AMDiS/src/PeriodicBC.cc | 23 ++--------
AMDiS/src/ProblemScal.cc | 6 +--
AMDiS/src/TimedObject.h | 39 ++++++++--------
4 files changed, 103 insertions(+), 64 deletions(-)
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 6f1accf3..fbd38c3f 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -77,31 +77,82 @@ namespace AMDiS {
Element *nb = const_cast<Element*>(neighbour_[i]);
int edgeNo = oppVertex_[i] = mel->getOppVertex(i);
- if (nb->getFirstChild() && (edgeNo != 2)) { // make nb nearest el.
+
+
+ if (nb->getFirstChild() && (edgeNo != 2)) {
+
+ // Search for the next neighbour. In many cases, the neighbour element
+ // may be refinemed in a way, such that there is no new vertex on the
+ // common edge. This situation is shown in the following picture:
+ //
+ // /|\
+ // / | \
+ // / | \
+ // /\ | \
+ // / \ | \
+ // / \| \
+ // -------------
+ //
+ // nb el
+ //
+ // Note that we know (because of the last if statement), that the
+ // neighbour element has children and the common edge is not the
+ // refinement edge, which has always the number 2, of our element.
+
+
if (edgeNo == 0) {
+ // The situation is as follows:
+ //
+ // -------
+ // \ /|\
+ // \ / | \
+ // \/ | \
+ // \ | \
+ // \ | \
+ // \| \
+ // -------
+ //
+ // nb el
+ // That means, the edge 0 of the same level neighbour is the common
+ // edge, i.e., the direct neighbour is the second child of the same
+ // level neighbour.
+
nb = neighbour_[i] = nb->getSecondChild();
} else {
+ // The situation is as shown in the picture above. So the next
+ // neighbour is the first child of the same level neighbour element.
nb = neighbour_[i] = nb->getFirstChild();
}
+ // In both cases the opp vertex number is 2, as one can see in the
+ // pictures above.
oppVertex_[i] = 2;
if (fill_opp_coords) {
if (nb->isNewCoordSet()) {
oppCoord_[i] = *(nb->getNewCoord());
} else {
- oppCoord_[i] = (macroNeighbour->coord[0] + macroNeighbour->coord[1]) * 0.5;
- }
+ // In both cases, that are shown in the pictures above, the opp
+ // vertex of the neighbour edge is the midpoint of the vertex 0
+ // and vertex 1 of the same level neighbour element.
+ oppCoord_[i] = (macroNeighbour->coord[0] +
+ macroNeighbour->coord[1]) * 0.5;
+ }
switch (i) {
case 0:
- if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
+ // The common edge is the edge 0 of this element.
+
+ switch (edgeNo) {
+ case 1:
neighbourCoord_[i][0] = macroNeighbour->coord[2];
neighbourCoord_[i][1] = macroNeighbour->coord[0];
- } else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
+ break;
+ case 0:
neighbourCoord_[i][0] = macroNeighbour->coord[1];
neighbourCoord_[i][1] = macroNeighbour->coord[2];
- } else {
+ break;
+ default:
ERROR_EXIT("Should not happen!\n");
}
@@ -109,13 +160,17 @@ namespace AMDiS {
break;
case 1:
- if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
+ // The commonedge is the edge 1 of this element.
+ switch (edgeNo) {
+ case 0:
neighbourCoord_[i][0] = macroNeighbour->coord[1];
neighbourCoord_[i][1] = macroNeighbour->coord[2];
- } else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
+ break;
+ case 1:
neighbourCoord_[i][0] = macroNeighbour->coord[2];
neighbourCoord_[i][1] = macroNeighbour->coord[0];
- } else {
+ break;
+ default:
ERROR_EXIT("Should not happen!\n");
}
@@ -123,15 +178,10 @@ namespace AMDiS {
break;
case 2:
- if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
- neighbourCoord_[i][0] = macroNeighbour->coord[2];
- neighbourCoord_[i][1] = macroNeighbour->coord[1];
- } else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
- neighbourCoord_[i][0] = macroNeighbour->coord[0];
- neighbourCoord_[i][1] = macroNeighbour->coord[2];
- } else {
- ERROR_EXIT("Should not happen!\n");
- }
+ // I've deleted here some code, be I think that this case is not
+ // possible. If an error occurs in this line, please check AMDiS
+ // revision <= 476 at the same position.
+ ERROR_EXIT("Should not happen!\n");
break;
default:
@@ -146,17 +196,24 @@ namespace AMDiS {
} else {
std::cout << " Neighbour " << j << ": not existing" << std::endl;
}
- std::cout << " OppVertex " << j << ": " << static_cast<int>(mel->getOppVertex(j)) << std::endl;
- std::cout << std::endl;
+ std::cout << " OppVertex " << j << ": "
+ << static_cast<int>(mel->getOppVertex(j))
+ << std::endl << std::endl;
}
ERROR_EXIT("should not happen!\n");
break;
}
}
} else {
+
+ // In this case, we know that the common edge is the refinement edge.
+ // This makes everything much more simpler, because we know that the
+ // next neighbour is equal to the samel level neighbour. If the same
+ // level neighbour would be refinement, also this element must to be
+ // refinement, because they share the refinement edge.
+
if (fill_opp_coords) {
oppCoord_[i] = macroNeighbour->coord[edgeNo];
-
neighbourCoord_[i] = macroNeighbour->coord;
}
}
diff --git a/AMDiS/src/PeriodicBC.cc b/AMDiS/src/PeriodicBC.cc
index f63282ef..a252749d 100644
--- a/AMDiS/src/PeriodicBC.cc
+++ b/AMDiS/src/PeriodicBC.cc
@@ -205,32 +205,19 @@ namespace AMDiS {
masterMatrix_ = NULL;
}
- matrix->print();
-
using namespace mtl;
DOFAdmin* admin = rowFESpace->getAdmin();
std::vector<int> dofMap(admin->getUsedSize());
- for (int i = 0; i < admin->getUsedSize(); i++) {
+ for (int i = 0; i < admin->getUsedSize(); i++)
dofMap[i] = (*associated)[i];
- std::cout << "map " << i << " to " << dofMap[i] << std::endl;
- }
// Compute reorder matrix (newRow and newCol yields transposed!!!)
matrix::traits::reorder<>::type R= matrix::reorder(dofMap);
- DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), B, D, E, TR;
-
- A*= 0.5;
- // Half of entries with decreased row index + half of the strict lower origing
- B= strict_upper(R) * A + strict_lower(A);
- // Add half of entries with same row and decreased column index + half of increased columns
- // Remark regarding column permutation: trans(strict_upper(trans(X))) == strict_lower(X)
- // B+= bands(trans(R), 0, 1) * (A * strict_lower(R) + lower(A)); -> make this work one day!!!
- TR= trans(R);
- D= bands(TR, 0, 1);
- E= A * strict_lower(R) + lower(A);
- B+= D * E;
- swap(A, B);
+ DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), C;
+
+ C = R * A * trans(R) + A;
+ A = 0.5 * C;
}
void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 4583f070..c71a268e 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -586,6 +586,8 @@ namespace AMDiS {
systemMatrix->removeRowsWithDBC(systemMatrix->getApplyDBCs());
+ systemMatrix->finishInsertion();
+
// TODO: ExitMatrix should be called after finishInsertion!
if (systemMatrix->getBoundaryManager())
systemMatrix->getBoundaryManager()->exitMatrix(systemMatrix);
@@ -594,8 +596,6 @@ namespace AMDiS {
if (solution->getBoundaryManager())
solution->getBoundaryManager()->exitVector(solution);
- systemMatrix->finishInsertion();
-
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n", TIME_USED(first,clock()));
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -606,8 +606,6 @@ namespace AMDiS {
#endif
createPrecon();
-
- systemMatrix->print();
}
void ProblemScal::writeResidualMesh(AdaptInfo *adaptInfo, const std::string name)
diff --git a/AMDiS/src/TimedObject.h b/AMDiS/src/TimedObject.h
index df73b4f4..7731144b 100644
--- a/AMDiS/src/TimedObject.h
+++ b/AMDiS/src/TimedObject.h
@@ -24,10 +24,6 @@
namespace AMDiS {
- // ===========================================================================
- // ===== class TimedObject ===================================================
- // ===========================================================================
-
/** \brief
* This class can be used as base class for time dependent objects where
* different objects refer to the same time. Herefore a pointer to
@@ -37,24 +33,25 @@ namespace AMDiS {
class TimedObject
{
public:
- /** \brief
- * Constructor.
- */
- TimedObject() : timePtr(NULL) {};
-
- /** \brief
- * Sets the time pointer.
- */
- inline void setTimePtr(double *timePtr_) { timePtr = timePtr_; };
-
- /** \brief
- * Returns the time pointer.
- */
- inline double *getTimePtr() { return timePtr; };
+ /// Constructor.
+ TimedObject()
+ : timePtr(NULL)
+ {}
+
+ /// Sets the time pointer.
+ inline void setTimePtr(double *timePtr_)
+ {
+ timePtr = timePtr_;
+ }
+
+ /// Returns the time pointer.
+ inline double *getTimePtr()
+ {
+ return timePtr;
+ }
+
protected:
- /** \brief
- * Pointer to the externally managed time value.
- */
+ /// Pointer to the externally managed time value.
double *timePtr;
};
--
GitLab