diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index 796639ee104df4f7fbda01e4f84f8abe87ec579a..66f1fe5ef28a3c70bc5936e3835f908fa4afe5d6 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -134,8 +134,16 @@ $(SOURCE_DIR)/TFQMR.h \ $(SOURCE_DIR)/TFQMR.hh \
$(SOURCE_DIR)/VecSymSolver.h \ $(SOURCE_DIR)/VecSymSolver.hh \
$(SOURCE_DIR)/UmfPackSolver.h \ $(SOURCE_DIR)/UmfPackSolver.hh \
$(SOURCE_DIR)/Lagrange.h $(SOURCE_DIR)/Line.h \
-$(SOURCE_DIR)/MacroElement.h $(SOURCE_DIR)/MacroWriter.h $(SOURCE_DIR)/Markings.h $(SOURCE_DIR)/Markings.hh $(SOURCE_DIR)/MemoryManager.h $(SOURCE_DIR)/Mesh.h $(SOURCE_DIR)/ODirSolver.h $(SOURCE_DIR)/ODirSolver.hh $(SOURCE_DIR)/OEMSolver.h \
-$(SOURCE_DIR)/OEMSolver.hh $(SOURCE_DIR)/OResSolver.h $(SOURCE_DIR)/OResSolver.hh $(SOURCE_DIR)/Parameters.h $(SOURCE_DIR)/Parametric.h $(SOURCE_DIR)/Preconditioner.h \
+$(SOURCE_DIR)/MacroElement.h $(SOURCE_DIR)/MacroWriter.h \
+$(SOURCE_DIR)/Markings.h $(SOURCE_DIR)/Markings.hh \
+$(SOURCE_DIR)/MemoryManager.h \
+$(SOURCE_DIR)/Mesh.h \
+$(SOURCE_DIR)/ODirSolver.h $(SOURCE_DIR)/ODirSolver.hh \
+$(SOURCE_DIR)/OEMSolver.h $(SOURCE_DIR)/OEMSolver.hh \
+$(SOURCE_DIR)/OResSolver.h $(SOURCE_DIR)/OResSolver.hh \
+$(SOURCE_DIR)/Parameters.h \
+$(SOURCE_DIR)/Parametric.h \
+$(SOURCE_DIR)/Preconditioner.h \
$(SOURCE_DIR)/Quadrature.h \
$(SOURCE_DIR)/RCNeighbourList.h $(SOURCE_DIRe)/RefinementManager.h $(SOURCE_DIR)/RefinementManager1d.h $(SOURCE_DIR)/RefinementManager2d.h $(SOURCE_DIR)/RefinementManager3d.h $(SOURCE_DIR)/TecPlotWriter.h $(SOURCE_DIR)/TecPlotWriter.hh $(SOURCE_DIR)/Tetrahedron.h \
$(SOURCE_DIR)/Traverse.h $(SOURCE_DIR)/Triangle.h $(SOURCE_DIR)/NonLinSolver.h $(SOURCE_DIR)/NonLinSolver.hh $(SOURCE_DIR)/ProblemInstat.h $(SOURCE_DIR)/ProblemInstat.cc $(SOURCE_DIR)/ProblemTimeInterface.h $(SOURCE_DIR)/ProblemNonLin.h $(SOURCE_DIR)/ProblemNonLin.cc \
diff --git a/AMDiS/bin/Makefile.in b/AMDiS/bin/Makefile.in
index 939c3ef65b2ca0ee719e0ca8c28b1abfa1be827e..5a344dab4e41f227649b97a7a777955a184b6e82 100644
--- a/AMDiS/bin/Makefile.in
+++ b/AMDiS/bin/Makefile.in
@@ -548,8 +548,16 @@ $(SOURCE_DIR)/TFQMR.h \ $(SOURCE_DIR)/TFQMR.hh \
$(SOURCE_DIR)/VecSymSolver.h \ $(SOURCE_DIR)/VecSymSolver.hh \
$(SOURCE_DIR)/UmfPackSolver.h \ $(SOURCE_DIR)/UmfPackSolver.hh \
$(SOURCE_DIR)/Lagrange.h $(SOURCE_DIR)/Line.h \
-$(SOURCE_DIR)/MacroElement.h $(SOURCE_DIR)/MacroWriter.h $(SOURCE_DIR)/Markings.h $(SOURCE_DIR)/Markings.hh $(SOURCE_DIR)/MemoryManager.h $(SOURCE_DIR)/Mesh.h $(SOURCE_DIR)/ODirSolver.h $(SOURCE_DIR)/ODirSolver.hh $(SOURCE_DIR)/OEMSolver.h \
-$(SOURCE_DIR)/OEMSolver.hh $(SOURCE_DIR)/OResSolver.h $(SOURCE_DIR)/OResSolver.hh $(SOURCE_DIR)/Parameters.h $(SOURCE_DIR)/Parametric.h $(SOURCE_DIR)/Preconditioner.h \
+$(SOURCE_DIR)/MacroElement.h $(SOURCE_DIR)/MacroWriter.h \
+$(SOURCE_DIR)/Markings.h $(SOURCE_DIR)/Markings.hh \
+$(SOURCE_DIR)/MemoryManager.h \
+$(SOURCE_DIR)/Mesh.h \
+$(SOURCE_DIR)/ODirSolver.h $(SOURCE_DIR)/ODirSolver.hh \
+$(SOURCE_DIR)/OEMSolver.h $(SOURCE_DIR)/OEMSolver.hh \
+$(SOURCE_DIR)/OResSolver.h $(SOURCE_DIR)/OResSolver.hh \
+$(SOURCE_DIR)/Parameters.h \
+$(SOURCE_DIR)/Parametric.h \
+$(SOURCE_DIR)/Preconditioner.h \
$(SOURCE_DIR)/Quadrature.h \
$(SOURCE_DIR)/RCNeighbourList.h $(SOURCE_DIRe)/RefinementManager.h $(SOURCE_DIR)/RefinementManager1d.h $(SOURCE_DIR)/RefinementManager2d.h $(SOURCE_DIR)/RefinementManager3d.h $(SOURCE_DIR)/TecPlotWriter.h $(SOURCE_DIR)/TecPlotWriter.hh $(SOURCE_DIR)/Tetrahedron.h \
$(SOURCE_DIR)/Traverse.h $(SOURCE_DIR)/Triangle.h $(SOURCE_DIR)/NonLinSolver.h $(SOURCE_DIR)/NonLinSolver.hh $(SOURCE_DIR)/ProblemInstat.h $(SOURCE_DIR)/ProblemInstat.cc $(SOURCE_DIR)/ProblemTimeInterface.h $(SOURCE_DIR)/ProblemNonLin.h $(SOURCE_DIR)/ProblemNonLin.cc \
diff --git a/AMDiS/compositeFEM/src/SubPolytope.cc b/AMDiS/compositeFEM/src/SubPolytope.cc
index dca19c6bf8f72f81528661b8a8386237fb849987..5f9c8de7ce22393ad645b0804effc32de3f07606 100644
--- a/AMDiS/compositeFEM/src/SubPolytope.cc
+++ b/AMDiS/compositeFEM/src/SubPolytope.cc
@@ -3,8 +3,7 @@
#include "SubElInfo.h"
#include "SubPolytope.h"
-bool
-SubPolytope::checkIntPoints()
+bool SubPolytope::checkIntPoints()
{
////////////////////////////////////////////////////////////////////////////
//
@@ -15,13 +14,14 @@ SubPolytope::checkIntPoints()
// false - no
////////////////////////////////////////////////////////////////////////////
- int i,j;
int zeroCounter;
- for (i = 0; i < numIntPoints; i++) {
+ for (int i = 0; i < numIntPoints; i++) {
zeroCounter = 0;
- for (j = 0; j < dim+1; j++) {
- if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) { zeroCounter++; }
+ for (int j = 0; j < dim+1; j++) {
+ if (fabs((*intPoints)[i][j]) <= 1.e-15 ) {
+ zeroCounter++;
+ }
}
/**
@@ -29,21 +29,27 @@ SubPolytope::checkIntPoints()
*
* "Inner" points on edges aren't equal to 0.0 in any component.
*/
- if ( dim == 1 && zeroCounter != 0 ) { return false; }
+ if (dim == 1 && zeroCounter != 0 ) {
+ return false;
+ }
/**
* Dimension 2
*
* "Inner" points on edges are equal to 0.0 in exactly 1 component.
*/
- if ( dim == 2 && zeroCounter != 1 ) { return false; }
+ if (dim == 2 && zeroCounter != 1) {
+ return false;
+ }
/**
* Dimension 3
*
* "Inner" points on edges are equal to 0.0 in exactly 2 components.
*/
- if ( dim == 3 && zeroCounter != 2 ) { return false; }
+ if (dim == 3 && zeroCounter != 2) {
+ return false;
+ }
}
return true;
@@ -57,7 +63,7 @@ SubPolytope::SubPolytope(const ElInfo *elInfo_,
intPoints(intPoints_),
numIntPoints(numIntPoints_)
{
- FUNCNAME("SubPolytope::SubPolytope");
+ FUNCNAME("SubPolytope::SubPolytope()");
dim = (*intPoints_)[0].getSize() - 1;
@@ -87,8 +93,7 @@ SubPolytope::SubPolytope(const ElInfo *elInfo_,
case 3:
if (numIntPoints == 3) {
createSubElementPolytopeIsSubElement2D3D();
- }
- else {
+ } else {
createSubElementsForPolytope3D(indexElVertInPol_);
}
break;
@@ -98,8 +103,8 @@ SubPolytope::SubPolytope(const ElInfo *elInfo_,
}
}
-void
-SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
+
+void SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
{
//**************************************************************************
// The intersection of the one-dimensional element (interval) divides
@@ -107,7 +112,7 @@ SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
// subelement to take.
//**************************************************************************
- FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement1D");
+ FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement1D()");
TEST_EXIT(dim == 1 && numIntPoints == 1)("invalid call of this routine\n");
@@ -158,8 +163,7 @@ SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
}
-void
-SubPolytope::createSubElementPolytopeIsSubElement2D3D()
+void SubPolytope::createSubElementPolytopeIsSubElement2D3D()
{
//**************************************************************************
// The intersection with element produced dim intersection points.
@@ -176,7 +180,7 @@ SubPolytope::createSubElementPolytopeIsSubElement2D3D()
// components.
//**************************************************************************
- FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement2D3D");
+ FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement2D3D()");
TEST_EXIT((dim == 2 && numIntPoints == 2) ||
(dim == 3 && numIntPoints == 3))
@@ -190,10 +194,11 @@ SubPolytope::createSubElementPolytopeIsSubElement2D3D()
* Get the vertex which - with the intersection points intPoints - forms
* a subelement of element.
*/
- int i,j;
- for (i = 0; i < numIntPoints; i++) {
- for (j = 0; j < dim+1; j++) {
- if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) { vertex[j] = 0.0; };
+ for (int i = 0; i < numIntPoints; i++) {
+ for (int j = 0; j < dim+1; j++) {
+ if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) {
+ vertex[j] = 0.0;
+ };
}
}
@@ -206,7 +211,7 @@ SubPolytope::createSubElementPolytopeIsSubElement2D3D()
* not to be an intersection point.
*/
(*subElVertices)[0] = vertex;
- for (i = 0; i < numIntPoints; i++) {
+ for (int i = 0; i < numIntPoints; i++) {
(*subElVertices)[i+1] = (*intPoints)[i];
}
@@ -222,22 +227,21 @@ SubPolytope::createSubElementPolytopeIsSubElement2D3D()
DELETE subElVertices;
}
-int
-SubPolytope::getIndexSecondFaceIntPoint0(int indexFirstFace, int dim)
-{
- int i;
- for (i = 0; i < dim+1; i++) {
+int SubPolytope::getIndexSecondFaceIntPoint0(int indexFirstFace, int dim)
+{
+ for (int i = 0; i < dim+1; i++) {
if ( fabs((*intPoints)[0][i]) <= 1.e-15 && i != indexFirstFace ) {
return i;
}
}
+
ERROR_EXIT("couldn't determine the second face for IntPoint0\n");
return -1;
}
-void
-SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
+
+void SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
{
//**************************************************************************
// The intersection with element produced four intersection points. Thus the
@@ -292,14 +296,13 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
int indexEdge[2]; // For a vertex lying on an edge indexEdge
// stores the indices of the two barycentric
// coordinates which are not equal to zero.
- int indexA;
- int indexB;
- int indexS_0;
- int indexS_1;
- int indexS_2;
- int indexS_3;
- int indexSecondFaceIntPoint0;
- int i,j,k;
+ int indexA = 0;
+ int indexB = 0;
+ int indexS_0 = 0;
+ int indexS_1 = 0;
+ int indexS_2 = 0;
+ int indexS_3 = 0;
+ int indexSecondFaceIntPoint0 = 0;
/**
* Get the second vertex of element lying in subpolytope.
@@ -312,10 +315,10 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
*/
// Get the edges including the intersection points.
- for (i = 0; i < numIntPoints; i++) {
- k = 0;
- for (j = 0; j < dim+1; j++) {
- if ( fabs((*intPoints)[i][j]) > 1.e-15 ) {
+ for (int i = 0; i < numIntPoints; i++) {
+ int k = 0;
+ for (int j = 0; j < dim+1; j++) {
+ if (fabs((*intPoints)[i][j]) > 1.e-15 ) {
indexEdge[k] = j;
k++;
}
@@ -327,9 +330,9 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
// Get the vertex of element adjacent with indexElVertInPol1 whose
// common edge with indexElVertInPol1 doesn't contain an
// intersection point, and store it in indexElVertInPol2.
- for (i = 0; i < dim+1; i++) {
- if ( intPointOnEdge[indexElVertInPol1][i] == false &&
- i != indexElVertInPol1 ) {
+ for (int i = 0; i < dim+1; i++) {
+ if (intPointOnEdge[indexElVertInPol1][i] == false &&
+ i != indexElVertInPol1 ) {
indexElVertInPol2 = i;
break;
}
@@ -343,23 +346,21 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
* face opposite to A.
*/
- if ( fabs((*intPoints)[0][indexElVertInPol1]) <= 1.e-15 ) {
+ if (fabs((*intPoints)[0][indexElVertInPol1]) <= 1.e-15) {
// (*intPoints)[0] lies in the face opposite to vertex
// /ref indexElVertInPol1.
indexA = indexElVertInPol2;
indexB = indexElVertInPol1;
- }
- else if ( fabs((*intPoints)[0][indexElVertInPol2]) <= 1.e-15 ) {
+ } else if (fabs((*intPoints)[0][indexElVertInPol2]) <= 1.e-15) {
// (*intPoints)[0] lies in the face opposite to vertex
// /ref indexElVertInPol2.
indexA = indexElVertInPol1;
indexB = indexElVertInPol2;
- }
- else {
+ } else {
ERROR_EXIT("couldn't determine A and B\n");
}
@@ -370,7 +371,7 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
// (*intPoints)[0] is a neighbour of A (A has been constructed this way).
indexS_0 = 0;
- if ( fabs((*intPoints)[1][indexB]) <= 1.e-15 ) {
+ if (fabs((*intPoints)[1][indexB]) <= 1.e-15) {
// (*intPoints)[1] lies in the face opposite to B, thus is a neighbour
// of A.
@@ -379,14 +380,13 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
- if ( fabs((*intPoints)[2][indexSecondFaceIntPoint0]) <= 1.e-15 ) {
+ if (fabs((*intPoints)[2][indexSecondFaceIntPoint0]) <= 1.e-15) {
// (*intPoints)[2] is neighbour of (*intPoints)[0]
indexS_2 = 3;
indexS_3 = 2;
- }
- else {
+ } else {
// (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
// plane
@@ -394,28 +394,26 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
indexS_2 = 2;
indexS_3 = 3;
}
- }
- else if ( fabs((*intPoints)[1][indexA]) <= 1.e-15 ) {
+ } else if (fabs((*intPoints)[1][indexA]) <= 1.e-15) {
// (*intPoints)[1] lies in the face opposite to A
indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
- if ( fabs((*intPoints)[1][indexSecondFaceIntPoint0]) <= 1.e-15 ) {
+ if (fabs((*intPoints)[1][indexSecondFaceIntPoint0]) <= 1.e-15) {
// (*intPoints)[1] is neighbour of (*intPoints)[0], but isn't
// neighbour of A
indexS_3 = 1;
- if ( fabs((*intPoints)[2][indexB]) <= 1.e-15 ) {
+ if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
// (*intPoints)[2] is neighbour of (*intPoints)[0] and neighbour of A
indexS_1 = 2;
indexS_2 = 3;
- }
- else {
+ } else {
// (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
// plane
@@ -423,22 +421,20 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
indexS_1 = 3;
indexS_2 = 2;
}
- }
- else {
+ } else {
// (*intPoints)[1] isn't neighbour of (*intPoints)[0], thus lies opposite
// to (*intPoints)[0] in the intersection plane
indexS_2 = 1;
- if ( fabs((*intPoints)[2][indexB]) <= 1.e-15 ) {
+ if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
// (*intPoints)[2] is neighbour of A
indexS_1 = 2;
indexS_3 = 3;
- }
- else {
+ } else {
// (*intPoints)[2] isn't neighbour of A
@@ -446,8 +442,7 @@ SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
indexS_3 = 2;
}
}
- }
- else {
+ } else {
ERROR_EXIT("IntPoint1 isn't either part of the face opposite to A nor part of the face opposite to B\n");
}
diff --git a/AMDiS/src/AbstractFunction.h b/AMDiS/src/AbstractFunction.h
index decfd6e8a046f79c2d133a45b668dc10f908950b..bb5925fd5f0bb660baa56896bcdd5c768f9a80e0 100644
--- a/AMDiS/src/AbstractFunction.h
+++ b/AMDiS/src/AbstractFunction.h
@@ -93,7 +93,9 @@ namespace AMDiS {
/** \brief
* Returns \ref degree_.
*/
- inline int getDegree() const { return degree_; };
+ inline int getDegree() const {
+ return degree_;
+ };
/** \brief
* Deligates the evaluation to overriden method f.
@@ -132,7 +134,9 @@ namespace AMDiS {
/** \brief
* Returns \ref degree_.
*/
- inline int getDegree() const { return degree_; };
+ inline int getDegree() const {
+ return degree_;
+ };
/** \brief
* function evaluation.
@@ -174,7 +178,9 @@ namespace AMDiS {
/** \brief
* Returns \ref degree_.
*/
- inline int getDegree() const { return degree_; };
+ inline int getDegree() const {
+ return degree_;
+ };
/** \brief
* function evaluation.
diff --git a/AMDiS/src/AdaptBase.h b/AMDiS/src/AdaptBase.h
index f754bbb7576fd6c2bcbef6238479bc80353e247d..83660464d99cda495369e6ebaa0a9b0ab2ecd3c1 100644
--- a/AMDiS/src/AdaptBase.h
+++ b/AMDiS/src/AdaptBase.h
@@ -69,7 +69,9 @@ namespace AMDiS {
/** \brief
* Returns \ref name_
*/
- inline const ::std::string& getName() const { return name_; };
+ inline const ::std::string& getName() const {
+ return name_;
+ };
/** \brief
* Returns \ref problemIteration_
@@ -85,7 +87,9 @@ namespace AMDiS {
/** \brief
* Returns \ref adaptInfo_
*/
- inline AdaptInfo *getAdaptInfo() { return adaptInfo_; };
+ inline AdaptInfo *getAdaptInfo() {
+ return adaptInfo_;
+ };
/** \brief
* Returns \ref problemTime_
@@ -101,7 +105,9 @@ namespace AMDiS {
/** \brief
* Returns \ref initialAdaptInfo_
*/
- inline AdaptInfo *getInitialAdaptInfo() { return initialAdaptInfo_; };
+ inline AdaptInfo *getInitialAdaptInfo() {
+ return initialAdaptInfo_;
+ };
protected:
/** \brief
diff --git a/AMDiS/src/AdaptInfo.h b/AMDiS/src/AdaptInfo.h
index 9aa4063ff2f725c563b01206b332f5cafd03e0b2..feb90de9dbd6d283a1e363c7e5521a7e70649997 100644
--- a/AMDiS/src/AdaptInfo.h
+++ b/AMDiS/src/AdaptInfo.h
@@ -197,8 +197,7 @@ namespace AMDiS {
* Destructor.
*/
virtual ~AdaptInfo() {
- int i;
- for(i = 0; i < scalContents.getSize(); i++) {
+ for (int i = 0; i < scalContents.getSize(); i++) {
delete scalContents[i];
}
};
@@ -219,9 +218,9 @@ namespace AMDiS {
* Returns whether space tolerance is reached.
*/
virtual bool spaceToleranceReached() {
- int i, size = scalContents.getSize();
- for(i = 0; i < size; i++) {
- if(!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance)) {
+ int size = scalContents.getSize();
+ for (int i = 0; i < size; i++) {
+ if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance)) {
return false;
}
}
@@ -232,7 +231,7 @@ namespace AMDiS {
* Returns whether space tolerance of component i is reached.
*/
virtual bool spaceToleranceReached(int i) {
- if(!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance)) {
+ if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance)) {
return false;
} else {
return true;
@@ -243,10 +242,9 @@ namespace AMDiS {
* Returns whether time tolerance is reached.
*/
virtual bool timeToleranceReached() {
- //if(timestep <= minTimestep) return true;
- int i, size = scalContents.getSize();
- for(i = 0; i < size; i++) {
- if(!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance)) {
+ int size = scalContents.getSize();
+ for (int i = 0; i < size; i++) {
+ if (!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance)) {
return false;
}
}
@@ -257,7 +255,7 @@ namespace AMDiS {
* Returns whether time tolerance of component i is reached.
*/
virtual bool timeToleranceReached(int i) {
- if(!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance)) {
+ if (!(scalContents[i]->est_t_sum < scalContents[i]->timeTolerance)) {
return false;
} else {
return true;
@@ -485,8 +483,11 @@ namespace AMDiS {
*/
inline double setTime(double t) {
time = t;
- if (time > endTime) time = endTime;
- if (time < startTime) time = startTime;
+ if (time > endTime)
+ time = endTime;
+ if (time < startTime)
+ time = startTime;
+
return time;
};
@@ -509,9 +510,13 @@ namespace AMDiS {
*/
inline double setTimestep(double t) {
timestep = t;
- if (timestep > maxTimestep) timestep = maxTimestep;
- if (timestep < minTimestep) timestep = minTimestep;
- if (time + timestep > endTime) timestep = endTime - time;
+ if (timestep > maxTimestep)
+ timestep = maxTimestep;
+ if (timestep < minTimestep)
+ timestep = minTimestep;
+ if (time + timestep > endTime)
+ timestep = endTime - time;
+
return timestep;
};
@@ -763,8 +768,14 @@ namespace AMDiS {
*/
int maxSolverIterations;
+ /** \brief
+ *
+ */
double solverTolerance;
+ /** \brief
+ *
+ */
double solverResidual;
/** \brief
diff --git a/AMDiS/src/AdaptInstationary.cc b/AMDiS/src/AdaptInstationary.cc
index 3d83ab1e358373d375b282b65d1df12c0f4835bf..d07de2be8a064d535b9b328b94b0247a802e0721 100644
--- a/AMDiS/src/AdaptInstationary.cc
+++ b/AMDiS/src/AdaptInstationary.cc
@@ -112,7 +112,7 @@ namespace AMDiS {
adaptInfo_->incTimestepIteration();
- if(!fixedTimestep_ &&
+ if (!fixedTimestep_ &&
!adaptInfo_->timeToleranceReached() &&
!adaptInfo_->getTimestep() <= adaptInfo_->getMinTimestep())
{
@@ -126,8 +126,8 @@ namespace AMDiS {
do {
problemIteration_->beginIteration(adaptInfo_);
- if(problemIteration_->oneIteration(adaptInfo_, FULL_ITERATION)) {
- if(!fixedTimestep_ &&
+ if (problemIteration_->oneIteration(adaptInfo_, FULL_ITERATION)) {
+ if (!fixedTimestep_ &&
!adaptInfo_->timeToleranceReached() &&
!adaptInfo_->getTimestep() <= adaptInfo_->getMinTimestep())
{
@@ -148,7 +148,7 @@ namespace AMDiS {
!adaptInfo_->getTimestep() <= adaptInfo_->getMinTimestep() &&
adaptInfo_->getTimestepIteration() <= adaptInfo_->getMaxTimestepIteration());
- if(!fixedTimestep_ && adaptInfo_->timeErrorLow()) {
+ if (!fixedTimestep_ && adaptInfo_->timeErrorLow()) {
adaptInfo_->setTimestep(adaptInfo_->getTimestep() *time_delta_2);
}
}
diff --git a/AMDiS/src/AdaptInstationary.h b/AMDiS/src/AdaptInstationary.h
index 1e08cfb0548156f3d9819161ffee7bbe23625da0..3be1793e96c41c4d81252b261215c6cf2f429419 100644
--- a/AMDiS/src/AdaptInstationary.h
+++ b/AMDiS/src/AdaptInstationary.h
@@ -149,12 +149,12 @@ namespace AMDiS {
/** \brief
* Parameter \f$ \delta_1 \f$ used in time step reduction
*/
- double time_delta_1;
+ double time_delta_1;
/** \brief
* Parameter \f$ \delta_2 \f$ used in time step enlargement
*/
- double time_delta_2;
+ double time_delta_2;
/** \brief
* If this parameter is 1 and the instationary problem is stable, hence the number
diff --git a/AMDiS/src/AdaptStationary.cc b/AMDiS/src/AdaptStationary.cc
index 458ee4a1dc13b4ae8bd62651b8bb15a0256396eb..62a297d153258dd0d2e69217db19c6314a117669 100644
--- a/AMDiS/src/AdaptStationary.cc
+++ b/AMDiS/src/AdaptStationary.cc
@@ -28,19 +28,19 @@ namespace AMDiS {
}
// adaption loop
- while(!adaptInfo_->spaceToleranceReached() &&
- (adaptInfo_->getSpaceIteration() < adaptInfo_->getMaxSpaceIteration() ||
- adaptInfo_->getMaxSpaceIteration() < 0) )
- {
- problemIteration_->beginIteration(adaptInfo_);
- Flag adapted = problemIteration_->oneIteration(adaptInfo_, FULL_ITERATION);
- problemIteration_->endIteration(adaptInfo_);
-
- if (!adapted)
- break;
-
- adaptInfo_->incSpaceIteration();
- }
+ while (!adaptInfo_->spaceToleranceReached() &&
+ (adaptInfo_->getSpaceIteration() < adaptInfo_->getMaxSpaceIteration() ||
+ adaptInfo_->getMaxSpaceIteration() < 0) ) {
+
+ problemIteration_->beginIteration(adaptInfo_);
+ Flag adapted = problemIteration_->oneIteration(adaptInfo_, FULL_ITERATION);
+ problemIteration_->endIteration(adaptInfo_);
+
+ if (!adapted)
+ break;
+
+ adaptInfo_->incSpaceIteration();
+ }
return 0;
}
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 9405f7745fb89415c438f6e84c288ea84dd8b267..74dbfa1f84618683674b7ac4d63b0d6564359c95 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -37,15 +37,13 @@ namespace AMDiS {
opt(optimized),
firstCall(true)
{
- int i;
-
const BasisFunction *psi = assembler->rowFESpace->getBasisFcts();
const BasisFunction *phi = assembler->colFESpace->getBasisFcts();
nRow = psi->getNumber();
nCol = phi->getNumber();
- switch(order) {
+ switch (order) {
case 0:
terms = op->zeroOrder;
break;
@@ -62,33 +60,26 @@ namespace AMDiS {
// check if all terms are symmetric
symmetric = true;
- for(i=0; i < static_cast<int>(terms.size()); i++) {
- if(!(terms[i])->isSymmetric()) {
+ for (int i=0; i < static_cast<int>(terms.size()); i++) {
+ if (!(terms[i])->isSymmetric()) {
symmetric = false;
break;
}
}
dim = assembler->rowFESpace->getMesh()->getDim();
- // // create quadrature
- // if(!quadrature) {
- // dim = assembler->rowFESpace->getMesh()->getDim();
- // int degree = op->getQuadratureDegree(order, type);
- // quadrature = Quadrature::provideQuadrature(dim, degree);
- // }
}
FastQuadrature *SubAssembler::updateFastQuadrature(FastQuadrature *quadFast,
const BasisFunction *psi,
Flag updateFlag)
{
- if(!quadFast) {
- quadFast =
- FastQuadrature::provideFastQuadrature(psi,
- *quadrature,
- updateFlag);
+ if (!quadFast) {
+ quadFast = FastQuadrature::provideFastQuadrature(psi,
+ *quadrature,
+ updateFlag);
} else {
- if(!quadFast->initialized(updateFlag))
+ if (!quadFast->initialized(updateFlag))
quadFast->init(updateFlag);
}
@@ -103,19 +94,19 @@ namespace AMDiS {
// set values at QPs invalid
::std::map<const DOFVectorBase<double>*, ValuesAtQPs*>::iterator it1;
- for(it1 = valuesAtQPs.begin(); it1 != valuesAtQPs.end(); ++it1) {
+ for (it1 = valuesAtQPs.begin(); it1 != valuesAtQPs.end(); ++it1) {
((*it1).second)->valid = false;
}
// set gradients at QPs invalid
::std::map<const DOFVectorBase<double>*, GradientsAtQPs*>::iterator it2;
- for(it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2) {
+ for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2) {
((*it2).second)->valid = false;
}
// calls initElement of each term
::std::vector<OperatorTerm*>::iterator it;
- for(it = terms.begin(); it != terms.end(); ++it) {
+ for (it = terms.begin(); it != terms.end(); ++it) {
(*it)->initElement(elInfo, this, quad);
}
}
@@ -128,19 +119,20 @@ namespace AMDiS {
const int numPoints = localQuad->getNumPoints();
// already calculated for this element ?
- if(coordsValid) {
+ if (coordsValid) {
return coordsAtQPs;
}
// not yet calcualted !
- if(coordsAtQPs) DELETE [] coordsAtQPs;
+ if (coordsAtQPs)
+ DELETE [] coordsAtQPs;
// allocate memory
coordsAtQPs = NEW WorldVector<double>[numPoints];
// set new values
WorldVector<double>* k = NULL;
int l;
- for(k=&(coordsAtQPs[0]), l=0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
+ for (k = &(coordsAtQPs[0]), l = 0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
elInfo->coordToWorld(localQuad->getLambda(l), k);
}
@@ -160,12 +152,12 @@ namespace AMDiS {
TEST_EXIT(vec)("no dof vector!\n");
- if(valuesAtQPs[vec] && valuesAtQPs[vec]->valid)
+ if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid)
return valuesAtQPs[vec]->values.getValArray();
Quadrature *localQuad = quad ? quad : quadrature;
- if(!valuesAtQPs[vec]) {
+ if (!valuesAtQPs[vec]) {
valuesAtQPs[vec] = new ValuesAtQPs;
}
valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());
@@ -176,10 +168,10 @@ namespace AMDiS {
(vec->getFESpace() == owner->rowFESpace) ||
(vec->getFESpace() == owner->colFESpace);
- if(opt && !quad && sameFESpaces) {
+ if (opt && !quad && sameFESpaces) {
const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
const BasisFunction *phi = owner->colFESpace->getBasisFcts();
- if(vec->getFESpace()->getBasisFcts() == psi) {
+ if (vec->getFESpace()->getBasisFcts() == psi) {
psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
} else if(vec->getFESpace()->getBasisFcts() == phi) {
phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);
@@ -192,8 +184,8 @@ namespace AMDiS {
//double* uhLoc = GET_MEMORY(double, basFcts->getNumber());
//vec->getLocalVector(elInfo->getElement(), uhLoc);
- if(opt && !quad && sameFESpaces) {
- if(psiFast->getBasisFunctions() == basFcts) {
+ if (opt && !quad && sameFESpaces) {
+ if (psiFast->getBasisFunctions() == basFcts) {
//psiFast->uhAtQp(uhLoc, values);
vec->getVecAtQPs(elInfo, NULL, psiFast, values);
} else if(phiFast->getBasisFunctions() == basFcts) {
@@ -225,7 +217,7 @@ namespace AMDiS {
TEST_EXIT(vec)("no dof vector!\n");
- if(gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid)
+ if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid)
return gradientsAtQPs[vec]->values.getValArray();
Quadrature *localQuad = quad ? quad : quadrature;
@@ -244,8 +236,8 @@ namespace AMDiS {
(vec->getFESpace() == owner->rowFESpace) ||
(vec->getFESpace() == owner->colFESpace);
- if(opt && !quad && sameFESpaces) {
- if(vec->getFESpace()->getBasisFcts() == psi) {
+ if (opt && !quad && sameFESpaces) {
+ if (vec->getFESpace()->getBasisFcts() == psi) {
psiFast = updateFastQuadrature(psiFast, psi, INIT_GRD_PHI);
} else if(vec->getFESpace()->getBasisFcts() == phi) {
phiFast = updateFastQuadrature(phiFast, phi, INIT_GRD_PHI);
@@ -253,33 +245,20 @@ namespace AMDiS {
}
// calculate new values
- //const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();
- //double* uhLoc = GET_MEMORY(double, basFcts->getNumber());
- //vec->getLocalVector(elInfo->getElement(), uhLoc);
-
- if(opt && !quad && sameFESpaces) {
+ if (opt && !quad && sameFESpaces) {
if(psiFast->getBasisFunctions() == basFcts) {
- //psiFast->grdUhAtQp(Lambda, uhLoc, values);
vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
} else {
- //phiFast->grdUhAtQp(Lambda, uhLoc, values);
vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
}
} else {
- //localQuad->grdUhAtQp(basFcts,
- // Lambda,
- // uhLoc,
- // values);
vec->getGrdAtQPs(elInfo, localQuad, NULL, values);
}
-
- //FREE_MEMORY(uhLoc, double, basFcts->getNumber());
-
+
gradientsAtQPs[vec]->valid = true;
- // return values
return values;
}
@@ -311,10 +290,8 @@ namespace AMDiS {
Quadrature *quad,
bool optimized)
{
- int i;
-
// check if a assembler is needed at all
- if(op->zeroOrder.size() == 0) {
+ if (op->zeroOrder.size() == 0) {
return NULL;
}
@@ -330,31 +307,31 @@ namespace AMDiS {
sort(opTerms.begin(), opTerms.end());
// check if a new assembler is needed
- if(quad) {
- for(i=0; i < static_cast<int>( subAssemblers->size()); i++) {
+ if (quad) {
+ for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
::std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
sort(assTerms.begin(), assTerms.end());
- if((opTerms == assTerms) &&
- ((*subAssemblers)[i]->getQuadrature() == quad))
- {
- return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
- }
+ if ((opTerms == assTerms) &&
+ ((*subAssemblers)[i]->getQuadrature() == quad)) {
+
+ return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
+ }
}
}
// check if all terms are pw_const
bool pwConst = true;
- for(i=0; i < static_cast<int>( op->zeroOrder.size()); i++) {
- if(!op->zeroOrder[i]->isPWConst()) {
+ for (int i = 0; i < static_cast<int>( op->zeroOrder.size()); i++) {
+ if (!op->zeroOrder[i]->isPWConst()) {
pwConst = false;
break;
}
}
// create new assembler
- if(!optimized) {
+ if (!optimized) {
newAssembler = NEW Stand0(op, assembler, quad);
} else {
if(pwConst) {
@@ -375,8 +352,6 @@ namespace AMDiS {
FirstOrderType type,
bool optimized)
{
- int i;
-
::std::vector<SubAssembler*> *subAssemblers =
optimized ?
(type == GRD_PSI ?
@@ -390,7 +365,7 @@ namespace AMDiS {
= (type == GRD_PSI) ? op->firstOrderGrdPsi : op->firstOrderGrdPhi;
// check if a assembler is needed at all
- if(opTerms.size() == 0) {
+ if (opTerms.size() == 0) {
return NULL;
}
@@ -399,38 +374,36 @@ namespace AMDiS {
FirstOrderAssembler *newAssembler;
// check if a new assembler is needed
- for(i=0; i < static_cast<int>( subAssemblers->size()); i++) {
+ for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
::std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
- sort(assTerms.begin(), assTerms.end());
+ sort(assTerms.begin(), assTerms.end());
-
+ if ((opTerms == assTerms) &&
+ ((*subAssemblers)[i]->getQuadrature() == quad)) {
- if((opTerms == assTerms) &&
- ((*subAssemblers)[i]->getQuadrature() == quad))
- {
- return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
- }
+ return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
+ }
}
// check if all terms are pw_const
bool pwConst = true;
- for(i=0; i < static_cast<int>( opTerms.size()); i++) {
- if(!(opTerms[i])->isPWConst()) {
+ for (int i = 0; i < static_cast<int>( opTerms.size()); i++) {
+ if (!(opTerms[i])->isPWConst()) {
pwConst = false;
break;
}
}
// create new assembler
- if(!optimized) {
+ if (!optimized) {
newAssembler =
(type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>(NEW Stand10(op, assembler, quad)) :
dynamic_cast<FirstOrderAssembler*>(NEW Stand01(op, assembler, quad));
} else {
- if(pwConst) {
+ if (pwConst) {
newAssembler =
(type == GRD_PSI) ?
dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
@@ -453,8 +426,6 @@ namespace AMDiS {
Quadrature *quad,
bool optimized)
{
- int i;
-
// check if a assembler is needed at all
if(op->secondOrder.size() == 0) {
return NULL;
@@ -472,32 +443,32 @@ namespace AMDiS {
sort(opTerms.begin(), opTerms.end());
// check if a new assembler is needed
- for(i=0; i < static_cast<int>( subAssemblers->size()); i++) {
+ for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
::std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
sort(assTerms.begin(), assTerms.end());
- if((opTerms == assTerms) &&
- ((*subAssemblers)[i]->getQuadrature() == quad))
- {
- return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
- }
+ if ((opTerms == assTerms) &&
+ ((*subAssemblers)[i]->getQuadrature() == quad)) {
+
+ return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
+ }
}
// check if all terms are pw_const
bool pwConst = true;
- for(i=0; i < static_cast<int>( op->secondOrder.size()); i++) {
- if(!op->secondOrder[i]->isPWConst()) {
+ for (int i = 0; i < static_cast<int>( op->secondOrder.size()); i++) {
+ if (!op->secondOrder[i]->isPWConst()) {
pwConst = false;
break;
}
}
// create new assembler
- if(!optimized) {
+ if (!optimized) {
newAssembler = NEW Stand2(op, assembler, quad);
} else {
- if(pwConst) {
+ if (pwConst) {
newAssembler = NEW Pre2(op, assembler, quad);
} else {
newAssembler = NEW Quad2(op, assembler, quad);
@@ -515,7 +486,6 @@ namespace AMDiS {
void Stand0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- int iq, i, j;
double val;
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
@@ -527,28 +497,28 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] = 0.0;
}
::std::vector<OperatorTerm*>::iterator termIt;
- for(termIt = terms.begin(); termIt != terms.end(); ++termIt) {
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
}
if (symmetric) {
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
- for(i=0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++) {
phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
}
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
(*mat)[i][i] += quadrature->getWeight(iq)*c[iq]*psival*phival[i];
- for (j = i+1; j < nCol; j++) {
+ for (int j = i + 1; j < nCol; j++) {
val = quadrature->getWeight(iq)*c[iq]*psival*phival[j];
(*mat)[i][j] += val;
(*mat)[j][i] += val;
@@ -556,17 +526,17 @@ namespace AMDiS {
}
}
} else { /* non symmetric assembling */
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] *= elInfo->getDet();
// calculate phi at QPs only once!
- for(i=0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++) {
phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
}
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
- for (j = 0; j < nCol; j++) {
+ for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq)*c[iq]*psival*phival[j];
}
}
@@ -578,28 +548,25 @@ namespace AMDiS {
void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
- double psi;
- int iq, i;
-
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] = 0.0;
}
::std::vector<OperatorTerm*>::iterator termIt;
- for(termIt = terms.begin(); termIt != terms.end(); ++termIt) {
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] *= elInfo->getDet();
- for (i = 0; i < nRow; i++) {
- psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
+ for (int i = 0; i < nRow; i++) {
+ double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
(quadrature->getLambda(iq));
- (*vec)[i] += quadrature->getWeight(iq)*c[iq]*psi;
+ (*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
}
}
FREE_MEMORY(c, double, numPoints);
@@ -608,29 +575,13 @@ namespace AMDiS {
Quad0::Quad0(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
{
- // if(!psiFast) {
- // psiFast =
- // FastQuadrature::provideFastQuadrature(assembler->
- // getRowFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_PHI);
- // }
- // if(!phiFast) {
- // phiFast =
- // FastQuadrature::provideFastQuadrature(assembler->
- // getColFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_PHI);
- // }
}
void Quad0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
const double *psi, *phi;
- int iq, i, j;
- double val;
- if(firstCall) {
+ if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -641,39 +592,39 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] = 0.0;
}
::std::vector<OperatorTerm*>::iterator termIt;
- for(termIt = terms.begin(); termIt != terms.end(); ++termIt) {
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
}
if (symmetric) {
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] *= elInfo->getDet();
psi = psiFast->getPhi(iq);
phi = phiFast->getPhi(iq);
- for (i = 0; i < nRow; i++) {
- (*mat)[i][i] += quadrature->getWeight(iq)*c[iq]*psi[i]*phi[i];
- for (j = i+1; j < nCol; j++) {
- val = quadrature->getWeight(iq)*c[iq]*psi[i]*phi[j];
+ for (int i = 0; i < nRow; i++) {
+ (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
+ for (int j = i + 1; j < nCol; j++) {
+ double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
}
} else { /* non symmetric assembling */
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] *= elInfo->getDet();
psi = psiFast->getPhi(iq);
phi = phiFast->getPhi(iq);
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++) {
- (*mat)[i][j] += quadrature->getWeight(iq)*c[iq]*psi[i]*phi[j];
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
+ (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
}
}
}
@@ -683,10 +634,7 @@ namespace AMDiS {
void Quad0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
- const double *psi;
- int iq, i;
-
- if(firstCall) {
+ if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -697,21 +645,21 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, numPoints);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] = 0.0;
}
::std::vector<OperatorTerm*>::iterator termIt;
- for(termIt = terms.begin(); termIt != terms.end(); ++termIt) {
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt) {
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
c[iq] *= elInfo->getDet();
- psi = psiFast->getPhi(iq);
- for (i = 0; i < nRow; i++) {
- (*vec)[i] += quadrature->getWeight(iq)*c[iq]*psi[i];
+ const double *psi = psiFast->getPhi(iq);
+ for (int i = 0; i < nRow; i++) {
+ (*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
}
}
FREE_MEMORY(c, double, numPoints);
@@ -720,19 +668,13 @@ namespace AMDiS {
Pre0::Pre0(Operator *op, Assembler *assembler, Quadrature *quad)
: ZeroOrderAssembler(op, assembler, quad, true)
{
- // q00 = Q00PsiPhi::provideQ00PsiPhi(assembler->getRowFESpace()->getBasisFcts(),
- // assembler->getColFESpace()->getBasisFcts(),
- // quadrature);
- // q0 = Q0Psi::provideQ0Psi(assembler->getRowFESpace()->getBasisFcts(),
- // quadrature);
}
void Pre0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- int i, j;
- double val, *c = GET_MEMORY(double, 1);
+ double *c = GET_MEMORY(double, 1);
- if(firstCall) {
+ if (firstCall) {
q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -742,24 +684,24 @@ namespace AMDiS {
}
c[0] = 0.0;
- for(i=0; i < static_cast<int>( terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>( terms.size()); i++) {
(static_cast<ZeroOrderTerm*>((terms[i])))->getC(elInfo, 1, c);
}
c[0] *= elInfo->getDet();
if (symmetric) {
- for (i = 0; i < nRow; i++) {
- (*mat)[i][i] += c[0]*q00->getValue(i,i);
- for (j = i + 1; j < nCol; j++) {
- val = c[0]*q00->getValue(i,j);
+ for (int i = 0; i < nRow; i++) {
+ (*mat)[i][i] += c[0] * q00->getValue(i,i);
+ for (int j = i + 1; j < nCol; j++) {
+ double val = c[0] * q00->getValue(i, j);
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
} else {
- for (i = 0; i < nRow; i++)
- for (j = 0; j < nCol; j++)
+ for (int i = 0; i < nRow; i++)
+ for (int j = 0; j < nCol; j++)
(*mat)[i][j] += c[0]*q00->getValue(i,j);
}
@@ -768,10 +710,9 @@ namespace AMDiS {
void Pre0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
- int i;
double *c = GET_MEMORY(double, 1);;
- if(firstCall) {
+ if (firstCall) {
q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -783,13 +724,13 @@ namespace AMDiS {
::std::vector<OperatorTerm*>::iterator termIt;
c[0] = 0.0;
- for(termIt = terms.begin(); termIt != terms.end(); ++termIt) {
+ for (termIt = terms.begin(); termIt != terms.end(); ++termIt) {
(static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
}
c[0] *= elInfo->getDet();
- for (i = 0; i < nRow; i++)
+ for (int i = 0; i < nRow; i++)
(*vec)[i] += c[0] * q0->getValue(i);
FREE_MEMORY(c, double, 1);
@@ -804,7 +745,6 @@ namespace AMDiS {
{
DimVec<double> grdPsi(dim, NO_INIT);
double *phival = GET_MEMORY(double, nCol);
- int iq, i, j;
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
@@ -812,72 +752,34 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
- // DimVec<double> *Lb = NEW DimVec<double>[numPoints](dim, NO_INIT);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq].set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, numPoints, Lb);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- for(i=0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++) {
phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
}
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
grdPsi = (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq));
- for (j = 0; j < nCol; j++) {
- (*mat)[i][j] +=
- quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
+ for (int j = 0; j < nCol; j++) {
+ (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
}
}
}
- // DELETE [] Lb;
FREE_MEMORY(phival, double, nCol);
}
- // void Stand10::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // DimVec<double> Lb(dim, NO_INIT);
- // DimVec<double> grdPsi(dim, NO_INIT);
- // int iq, i;
-
- // const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
-
- // for (iq = 0; iq < quadrature->getNumPoints(); iq++) {
- // Lb.set(0.0);
- // for(i=0; i < static_cast<int>( terms->size()); i++) {
- // (static_cast<FirstOrderTerm*>(((*terms)[i]))->eval(elInfo, 0, &Lb);
- // }
-
- // Lb *= elInfo->getDet();
-
- // for (i = 0; i < nRow; i++) {
- // grdPsi = (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq));
- // vec[i] += quadrature->getWeight(iq) * (Lb * grdPsi);
- // }
- // }
- // }
-
Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
{
- // if(!psiFast) {
- // psiFast = FastQuadrature::provideFastQuadrature(
- // assembler->getRowFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_GRD_PHI);
- // }
- // if(!phiFast) {
- // phiFast = FastQuadrature::provideFastQuadrature(
- // assembler->getColFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_PHI);
- // }
}
@@ -885,9 +787,8 @@ namespace AMDiS {
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
const double *phi;
- int iq, i, j;
- if(firstCall) {
+ if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -898,71 +799,42 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
- // DimVec<double> *Lb = NEW DimVec<double>[numPoints](dim, NO_INIT);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq].set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, numPoints, Lb);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq] *= elInfo->getDet();
phi = phiFast->getPhi(iq);
grdPsi = psiFast->getGradient(iq);
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++)
- (*mat)[i][j] +=
- quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++)
+ (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
}
}
- // DELETE [] Lb;
}
- // void Quad10::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // DimVec<double> Lb(dim, NO_INIT);
- // VectorOfFixVecs<DimVec<double> > *grdPsi;
- // int iq, i;
-
- // for (iq = 0; iq < quadrature->getNumPoints(); iq++) {
- // Lb.set(0.0);
- // for(i=0; i < static_cast<int>( terms->size()); i++) {
- // (static_cast<FirstOrderTerm*>(((*terms)[i]))->eval(elInfo, 0, &Lb);
- // }
-
- // Lb *= elInfo->getDet();
-
- // grdPsi = psiFast->getGradient(iq);
-
- // for (i = 0; i < nRow; i++) {
- // vec[i] += quadrature->getWeight(iq) * (Lb * (*grdPsi)[i]);
- // }
- // }
- // }
Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
{
- // q10 = Q10PsiPhi::provideQ10PsiPhi(assembler->getRowFESpace()->getBasisFcts(),
- // assembler->getColFESpace()->getBasisFcts(),
- // quadrature);
- // q1 = Q1Psi::provideQ1Psi(assembler->getRowFESpace()->getBasisFcts(),
- // quadrature);
}
+
void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
- //DimVec<double> *Lb = NEW DimVec<double>[1](dim, NO_INIT);
const int *k;
const double *values;
- int i, j, m;
double val;
- if(firstCall) {
+ if (firstCall) {
q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -974,59 +846,34 @@ namespace AMDiS {
const int **nEntries = q10->getNumberEntries();
Lb[0].set(0.0);
- for(i=0; i < static_cast<int>( terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>( terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, 1, Lb);
}
Lb[0] *= elInfo->getDet();
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++) {
- k = q10->getKVec(i, j);
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
+ k = q10->getKVec(i, j);
values = q10->getValVec(i, j);
+ int m;
for (val = m = 0; m < nEntries[i][j]; m++)
val += values[m]*Lb[0][k[m]];
(*mat)[i][j] += val;
}
}
-
- // DELETE [] Lb;
}
- // void Pre10::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // DimVec<double> Lb(dim, NO_INIT);
- // const int *nEntries = q1->getNumberEntries();
- // const int *k;
- // const double *values;
- // int i, m;
- // double val;
-
- // Lb.set(0.0);
- // for(i=0; i < static_cast<int>( terms->size()); i++) {
- // (static_cast<FirstOrderTerm*>(((*terms)[i]))->eval(elInfo, 0, &Lb);
- // }
-
- // Lb *= elInfo->getDet();
-
- // for (i = 0; i < nRow; i++) {
- // k = q1->getKVec(i);
- // values = q1->getValVec(i);
- // for (val = m = 0; m < nEntries[i]; m++)
- // val += values[m]*Lb[k[m]];
- // vec[i] += val;
- // }
- // }
Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
{}
+
void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
double psival;
- int iq, i, j;
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
@@ -1036,54 +883,49 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
// DimVec<double> *Lb = NEW DimVec<double>[numPoints](dim, NO_INIT);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq].set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, numPoints, Lb);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- for(i=0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++) {
grdPhi[i] = (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq));
}
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
- for (j = 0; j < nCol; j++)
- (*mat)[i][j] +=
- quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
+ for (int j = 0; j < nCol; j++)
+ (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
}
- }
-
- // DELETE [] Lb;
+ }
}
void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
DimVec<double> grdPsi(dim, NO_INIT);
- int iq, i;
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
- // DimVec<double> *Lb = NEW DimVec<double>[numPoints](dim, NO_INIT);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq].set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, numPoints, Lb);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
grdPsi = (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq));
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
}
@@ -1093,27 +935,13 @@ namespace AMDiS {
Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad)
: FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
{
- // if(!psiFast) {
- // psiFast = FastQuadrature::provideFastQuadrature(
- // assembler->getRowFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_PHI);
- // }
- // if(!phiFast) {
- // phiFast = FastQuadrature::provideFastQuadrature(
- // assembler->getColFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_GRD_PHI);
- // }
}
void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > *grdPhi;
- const double *psi;
- int iq, i, j;
- if(firstCall) {
+ if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -1124,37 +952,32 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
- // DimVec<double> *Lb = NEW DimVec<double>[numPoints](dim, NO_INIT);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq].set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, numPoints, Lb);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq] *= elInfo->getDet();
- psi = psiFast->getPhi(iq);
+ const double *psi = psiFast->getPhi(iq);
grdPhi = phiFast->getGradient(iq);
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++)
- (*mat)[i][j] +=
- quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++)
+ (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
}
}
-
- // DELETE [] Lb;
}
void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
VectorOfFixVecs<DimVec<double> > *grdPsi;
- int iq, i;
- if(firstCall) {
+ if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -1165,26 +988,24 @@ namespace AMDiS {
int numPoints = quadrature->getNumPoints();
VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
- // DimVec<double> *Lb = NEW DimVec<double>[numPoints](dim, NO_INIT);
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq].set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, numPoints, Lb);
}
- for (iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
Lb[iq] *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
}
}
- // DELETE [] Lb;
}
Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad)
@@ -1195,14 +1016,12 @@ namespace AMDiS {
void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
- //DimVec<double> *Lb = NEW DimVec<double>[1](dim, NO_INIT);
const int *l;
const double *values;
- int i, j, m;
double val;
- if(firstCall) {
+ if (firstCall) {
q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -1214,28 +1033,27 @@ namespace AMDiS {
const int **nEntries = q01->getNumberEntries();
Lb[0].set(0.0);
- for(i=0; i < static_cast<int>( terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>( terms.size()); i++) {
(static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, 1, Lb);
}
Lb[0] *= elInfo->getDet();
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++) {
- l = q01->getLVec(i, j);
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
+ l = q01->getLVec(i, j);
values = q01->getValVec(i, j);
+ int m;
for (val = m = 0; m < nEntries[i][j]; m++)
val += values[m]*Lb[0][l[m]];
(*mat)[i][j] += val;
}
}
- // DELETE [] Lb;
}
void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
{
VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);
- //DimVec<double> *Lb = NEW DimVec<double>[1](dim, NO_INIT);
const int *k;
const double *values;
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index e5d8131ebcbd787ba1c2de442ea3039b42fe32ae..ca874eb4d68bcd89c314105cfb22aed39a060cd8 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -117,7 +117,9 @@ namespace AMDiS {
/** \brief
* Returns \ref terms
*/
- inline ::std::vector<OperatorTerm*> *getTerms() { return &terms; };
+ inline ::std::vector<OperatorTerm*> *getTerms() {
+ return &terms;
+ };
/** \brief
* Returns \ref quadrature.
@@ -168,12 +170,16 @@ namespace AMDiS {
/** \brief
* Returns \ref psiFast.
*/
- const FastQuadrature *getPsiFast() const { return psiFast; };
+ const FastQuadrature *getPsiFast() const {
+ return psiFast;
+ };
/** \brief
* Returns \ref phiFast.
*/
- const FastQuadrature *getPhiFast() const { return phiFast; };
+ const FastQuadrature *getPhiFast() const {
+ return phiFast;
+ };
protected:
/** \brief
@@ -270,8 +276,14 @@ namespace AMDiS {
*/
::std::vector<OperatorTerm*> terms;
+ /** \brief
+ *
+ */
bool opt;
+ /** \brief
+ *
+ */
bool firstCall;
friend class Assembler;
@@ -915,32 +927,44 @@ namespace AMDiS {
/** \brief
* Returns \ref rowFESpace.
*/
- inline const FiniteElemSpace* getRowFESpace() { return rowFESpace; };
+ inline const FiniteElemSpace* getRowFESpace() {
+ return rowFESpace;
+ };
/** \brief
* Returns \ref colFESpace.
*/
- inline const FiniteElemSpace* getColFESpace() { return colFESpace; };
+ inline const FiniteElemSpace* getColFESpace() {
+ return colFESpace;
+ };
/** \brief
* Returns \ref nRow.
*/
- inline int getNRow() { return nRow; };
+ inline int getNRow() {
+ return nRow;
+ };
/** \brief
* Returns \ref nCol.
*/
- inline int getNCol() { return nCol; };
+ inline int getNCol() {
+ return nCol;
+ };
/** \brief
* Sets \ref rememberElMat.
*/
- inline void rememberElementMatrix(bool rem) { rememberElMat = rem; };
+ inline void rememberElementMatrix(bool rem) {
+ rememberElMat = rem;
+ };
/** \brief
* Sets \ref rememberElVec.
*/
- inline void rememberElementVector(bool rem) { rememberElVec = rem; };
+ inline void rememberElementVector(bool rem) {
+ rememberElVec = rem;
+ };
/** \brief
* Returns \ref zeroOrderAssembler.
@@ -972,7 +996,9 @@ namespace AMDiS {
/** \brief
* Returns \ref operat;
*/
- inline Operator* getOperator() { return operat; };
+ inline Operator* getOperator() {
+ return operat;
+ };
/** \brief
* Initialisation for the given ElInfo. The call is deligated to
@@ -1022,7 +1048,7 @@ namespace AMDiS {
* Checks whether this is a new travese.
*/
inline void checkForNewTraverse() {
- if(lastTraverseId != ElInfo::traverseId) {
+ if (lastTraverseId != ElInfo::traverseId) {
lastVecEl = lastMatEl = NULL;
lastTraverseId = ElInfo::traverseId;
}
diff --git a/AMDiS/src/Boundary.h b/AMDiS/src/Boundary.h
index 022d3153601da6543e6179c2ac114a81abe752f9..a0610a8b3cff75630dae696c1c5132a830d080d3 100644
--- a/AMDiS/src/Boundary.h
+++ b/AMDiS/src/Boundary.h
@@ -34,10 +34,10 @@ namespace AMDiS {
* Contains boundary constants
*/
typedef enum {
- INTERIOR=0, /**< interior => no boundary (b = 0) */
- DIRICHLET=1, /**< dirichlet boundary (b > 0) */
- NEUMANN=-1, /**< neumann boundary (-100 < b < 0) */
- ROBIN=-100 /**< robin boundary ( b <= -100) */
+ INTERIOR = 0, /**< interior => no boundary (b = 0) */
+ DIRICHLET = 1, /**< dirichlet boundary (b > 0) */
+ NEUMANN = -1, /**< neumann boundary (-100 < b < 0) */
+ ROBIN = -100 /**< robin boundary (b <= -100) */
} BoundaryConstants;
/** \brief
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 1ad36621bdd501e8d638c3888e9babcdee881c9b..bf92f4e0eafabe6679479265f268d6020bc4f030 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -27,12 +27,17 @@ namespace AMDiS {
projection_(mesh_->getDim(), NO_INIT),
oppCoord_(mesh_->getDim(), NO_INIT),
neighbour_(mesh_->getDim(), NO_INIT),
+ neighbourCoord_(mesh_->getDim(), NO_INIT),
oppVertex_(mesh_->getDim(), NO_INIT),
grdLambda_(mesh_->getDim(), NO_INIT)
//parametric_(false)
{
projection_.set(NULL);
+
+ for (int i = 0; i < neighbourCoord_.getSize(); i++) {
+ neighbourCoord_[i].init(mesh_->getDim());
+ }
}
ElInfo::~ElInfo()
@@ -57,29 +62,23 @@ namespace AMDiS {
WorldVector<double> *w)
{
- static WorldVector<double> world;
- double c;
- WorldVector<double> *ret;
- WorldVector<double> v;
- int i, j;
-
int dim = l.getSize() - 1;
int dimOfWorld = Global::getGeo(WORLD);
+ static WorldVector<double> world;
+ WorldVector<double> *ret = w ? w : &world;
+ WorldVector<double> v = (*coords)[0];
+ double c = l[0];
- ret = w ? w : &world;
-
- v = (*coords)[0];
- c = l[0];
- for (j = 0; j < dimOfWorld; j++)
+ for (int j = 0; j < dimOfWorld; j++)
(*ret)[j] = c*v[j];
int vertices = Global::getGeo(VERTEX, dim);
- for (i = 1; i < vertices; i++) {
+ for (int i = 1; i < vertices; i++) {
v = (*coords)[i];
c = l[i];
- for (j = 0; j < dimOfWorld; j++)
+ for (int j = 0; j < dimOfWorld; j++)
(*ret)[j] += c*v[j];
}
return ret;
@@ -97,9 +96,8 @@ namespace AMDiS {
double ElInfo::calcDet(const FixVec<WorldVector<double>, VERTEX> &coords)
{
- FUNCNAME("ElInfo::calcDet");
+ FUNCNAME("ElInfo::calcDet()");
- int i;
int dim = coords.getSize() - 1;
int dow = Global::getGeo(WORLD);
@@ -107,17 +105,15 @@ namespace AMDiS {
double det = 0.0;
- // --> changes Christina
if (dim == 0)
return 1.0;
- // --> end: changes Christina
// dim = dim of world
WorldVector<double> e1, e2, e3, v0;
v0 = coords[0];
- for (i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++) {
e1[i] = coords[1][i] - v0[i];
if(dim > 1)
e2[i] = coords[2][i] - v0[i];
@@ -140,18 +136,18 @@ namespace AMDiS {
{
WorldVector<double> n;
- if(dim > 1) {
+ if (dim > 1) {
n[0] = e1[1]*e2[2] - e1[2]*e2[1];
n[1] = e1[2]*e2[0] - e1[0]*e2[2];
n[2] = e1[0]*e2[1] - e1[1]*e2[0];
}
- if(dim == 1) {
+ if (dim == 1) {
det = norm(&e1);
} else if (dim == 2) {
det = norm(&n);
} else if (dim == 3) {
- det = n[0]*e3[0]+n[1]*e3[1]+n[2]*e3[2];
+ det = n[0] * e3[0] + n[1] * e3[1] + n[2] * e3[2];
} else
ERROR_EXIT("not yet for problem dimension = %d",dim);
break;
@@ -171,10 +167,10 @@ namespace AMDiS {
void ElInfo::fillDetGrdLambda()
{
- if(fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
+ if (fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
det_ = calcGrdLambda(grdLambda_);
} else {
- if(fillFlag_.isSet(Mesh::FILL_DET)) {
+ if (fillFlag_.isSet(Mesh::FILL_DET)) {
det_ = calcDet();
}
}
@@ -190,15 +186,9 @@ namespace AMDiS {
int dim = mesh_->getDim();
int posIndex = DIM_OF_INDEX(pos, dim);
- int offset = indexOffset[dim-1][posIndex];
+ int offset = indexOffset[dim - 1][posIndex];
return boundary_[offset + i];
}
- // BoundaryType ElInfo::getBoundOfBoundary(int i) const {
- // WARNING("USE: getBoundary() instead of getBoundOfBoundary()\n");
- // TEST_EXIT(boundary_)("no boundary\n");
- // return (*boundary_)[i];
- // }
-
}
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 179df71aaa0d7635d0c4b293cd2f15f81560438d..3e02ca564f57b07e3bbb1c7c5ce77506e409a747 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -94,6 +94,7 @@ namespace AMDiS {
boundary_ = rhs.boundary_;
oppCoord_ = rhs.oppCoord_;
neighbour_ = rhs.neighbour_;
+ neighbourCoord_ = rhs.neighbourCoord_;
oppVertex_ = rhs.oppVertex_;
return *this;
};
@@ -167,7 +168,7 @@ namespace AMDiS {
* Get ElInfo's \ref coord_. This is a FixVec<WorldVector<double> > filled with the
* coordinates of the all vertice of element \ref el.
*/
- inline FixVec<WorldVector<double>,VERTEX>& getCoords() {
+ inline FixVec<WorldVector<double>, VERTEX>& getCoords() {
return coord_;
};
@@ -175,7 +176,7 @@ namespace AMDiS {
* Get ElInfo's \ref coord_. This is a FixVec<WorldVector<double> > filled with the
* coordinates of the all vertice of element \ref el.
*/
- inline const FixVec<WorldVector<double>,VERTEX>& getCoords() const {
+ inline const FixVec<WorldVector<double>, VERTEX>& getCoords() const {
return coord_;
};
@@ -205,6 +206,13 @@ namespace AMDiS {
return neighbour_[i];
};
+ /** \brief
+ * Get ElInfo's \ref neighbourCoord_[i]
+ */
+ inline FixVec<WorldVector<double>, VERTEX> getNeighbourCoord(int i) const {
+ return neighbourCoord_[i];
+ }
+
/** \brief
* Get ElInfo's \ref oppVertex_[i]
*/
@@ -398,7 +406,7 @@ namespace AMDiS {
* parent data parentInfo.
* pure virtual => must be overriden in sub-class.
*/
- virtual void fillElInfo(int ichild, const ElInfo *parentInfo)=0;
+ virtual void fillElInfo(int ichild, const ElInfo *parentInfo) = 0;
/** \brief
* calculates the Jacobian of the barycentric coordinates on \element and stores
@@ -425,7 +433,8 @@ namespace AMDiS {
*/
virtual double getElementNormal(WorldVector<double> &elementNormal) const
{
- FUNCNAME("ElInfo::getElementNormal");
+ FUNCNAME("ElInfo::getElementNormal()");
+
ERROR("virtual function not implemented in this sub-class ");
return(0.0);
@@ -484,13 +493,13 @@ namespace AMDiS {
* of the children is (level of the parent + 1). level_ is filled always by
* the traversal routines.
*/
- unsigned char level_;
+ unsigned char level_;
/** \brief
* \ref coord_[i] is a WorldVector<double> storing the world coordinates of the
* i-th vertex of element \ref element.
*/
- FixVec<WorldVector<double>,VERTEX> coord_;
+ FixVec<WorldVector<double>, VERTEX> coord_;
/** \brief
* boundary_[i] is the BoundaryType of the i-th edge/face
@@ -499,7 +508,7 @@ namespace AMDiS {
* object of the i-th edge, for i=0,..,N_EDGES - 1. It is
* a pointer to NULL for an interior edge/face.
*/
- FixVec<BoundaryType ,BOUNDARY> boundary_;
+ FixVec<BoundaryType, BOUNDARY> boundary_;
/** \brief
* Vector storing pointers to projections for each face, edge, vertex.
@@ -510,20 +519,26 @@ namespace AMDiS {
* oppCoord_[i] coordinates of the i-th neighbour vertex opposite the
* common edge/face.
*/
- FixVec<WorldVector<double>,NEIGH> oppCoord_;
+ FixVec<WorldVector<double>, NEIGH> oppCoord_;
/** \brief
* neighbour_[i] pointer to the element at the edge/face with local index i.
* It is a pointer to NULL for boundary edges/faces.
*/
- FixVec<Element*,NEIGH> neighbour_;
+ FixVec<Element*, NEIGH> neighbour_;
+
+ /** \brief
+ * neighbourCoord_[i][j] are the coordinate of the j-th vertex of the i-th
+ * neighbour element with the common edge/face.
+ */
+ FixVec<FixVec<WorldVector<double>, VERTEX>, NEIGH> neighbourCoord_;
/** \brief
* oppVertex_[i] is undefined if neighbour_[i] is a pointer to NULL.
* Otherwise it is the local index of the neighbour's vertex opposite the
* common edge/face.
*/
- FixVec<unsigned char,NEIGH> oppVertex_;
+ FixVec<unsigned char, NEIGH> oppVertex_;
/** \brief
* Elements determinant.
@@ -536,7 +551,7 @@ namespace AMDiS {
DimVec<WorldVector<double> > grdLambda_;
/** \brief
- * True, if this elInfo stores parametrized information. Flase, otherwise.
+ * True, if this elInfo stores parametrized information. False, otherwise.
*/
bool parametric_;
@@ -549,7 +564,7 @@ namespace AMDiS {
static const int childVertex[3][2][4];
/** \brief
- *child_edge[el_type][child][i] = father's local edge index of new edge i.
+ * child_edge[el_type][child][i] = father's local edge index of new edge i.
* new edge 2 is half of old edge 0, new edges 4,5 are really new edges, and
* value is different: child_edge[][][4,5] = index of same edge in other
* child.
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 6a249bbf635c9b425e6e46c395c1de7716c75296..ac72bca6670dca5345d2c74a69c15749e08c1a45 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -16,76 +16,67 @@ namespace AMDiS {
void ElInfo1d::fillMacroInfo(const MacroElement * mel)
{
- FUNCNAME("ElInfo1d::fillMacroInfo");
+ FUNCNAME("ElInfo1d::fillMacroInfo()");
Element *nb;
MacroElement *mnb;
- int i;
- macroElement_ = const_cast<MacroElement*>( mel);
- element_ = const_cast<Element*>( mel->getElement());
- parent_ = NULL;
- level_ = 0;
+ macroElement_ = const_cast<MacroElement*>( mel);
+ element_ = const_cast<Element*>( mel->getElement());
+ parent_ = NULL;
+ level_ = 0;
int vertices = mesh_->getGeo(VERTEX);
if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
- fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA))
- {
- for (i = 0; i < vertices; i++)
- {
- coord_[i] = mel->coord[i];
- }
+ fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
+
+ for (int i = 0; i < vertices; i++) {
+ coord_[i] = mel->coord[i];
+ }
}
// if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
// det = elGrdLambda(*Lambda);
// }
- if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- WorldVector<double> oppC;
-
- int neighbours = mesh_->getGeo(NEIGH);
- for (i = 0; i < neighbours; i++)
- {
- nb = NULL;
- if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i))))
- {
- if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- oppC = mnb->coord[i];
- }
-
- nb = const_cast<Element*>( mnb->getElement());
-
- while (!(nb->isLeaf())) // make nb nearest element
- {
- if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- if (nb->getNewCoord(-1)) {
- oppC = *(nb->getNewCoord());
- } else {
- oppC = (mel->coord[i] + oppC) * 0.5;
- }
- }
- nb = const_cast<Element*>( nb->getChild(1-i));
- }
-
- if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- oppCoord_[i] = oppC;
- }
+ if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ WorldVector<double> oppC;
+
+ int neighbours = mesh_->getGeo(NEIGH);
+ for (int i = 0; i < neighbours; i++) {
+ nb = NULL;
+ if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ oppC = mnb->coord[i];
+ }
+
+ nb = const_cast<Element*>( mnb->getElement());
+
+ while (!(nb->isLeaf())) { // make nb nearest element
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ if (nb->getNewCoord(-1)) {
+ oppC = *(nb->getNewCoord());
+ } else {
+ oppC = (mel->coord[i] + oppC) * 0.5;
}
- neighbour_[i] = nb;
- oppVertex_[i] = nb ? i : -1;
+ }
+ nb = const_cast<Element*>( nb->getChild(1-i));
+ }
+
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ oppCoord_[i] = oppC;
}
+ }
+ neighbour_[i] = nb;
+ oppVertex_[i] = nb ? i : -1;
}
+ }
if (fillFlag_.isSet(Mesh::FILL_BOUND) ) {
- for (i = 0; i < vertices; i++)
+ for (int i = 0; i < vertices; i++)
boundary_[i] = mel->getBoundary(i);
- for(i = 0; i < element_->getGeo(PROJECTION); i++) {
+ for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
projection_[i] = mel->getProjection(i);
}
}
@@ -98,24 +89,21 @@ namespace AMDiS {
/****************************************************************************/
double ElInfo1d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
{
- FUNCNAME("ElInfo1d::calcGrdLambda\n");
-
- int i;
- double adet2;
+ FUNCNAME("ElInfo1d::calcGrdLambda()");
testFlag(Mesh::FILL_COORDS);
WorldVector<double> e;
-
- e = coord_[1]; e -= coord_[0];
- adet2 = e*e;
+ e = coord_[1];
+ e -= coord_[0];
+ double adet2 = e * e;
if (adet2 < 1.0E-15) {
MSG("det*det = %lf\n", adet2);
- for (i = 0; i <= 1; i++)
+ for (int i = 0; i <= 1; i++)
grd_lam[i] = 0.0;
} else {
- grd_lam[1] = e * (1.0/adet2);
+ grd_lam[1] = e * (1.0 / adet2);
grd_lam[0] = grd_lam[1] * (-1.0);
}
@@ -125,12 +113,11 @@ namespace AMDiS {
const int ElInfo1d::worldToCoord(const WorldVector<double>& x,
DimVec<double>* lambda) const
{
- FUNCNAME("ElInfo1d::worldToCoord");
+ FUNCNAME("ElInfo1d::worldToCoord()");
+
double lmin;
double a = coord_[0][0];
double length = (coord_[1][0] - a);
- int i, j, k;
-
int dim = mesh_->getDim();
static DimVec<double> vec(dim, NO_INIT);
@@ -147,10 +134,10 @@ namespace AMDiS {
(*lambda)[1] = (x[0]-a)/length;
(*lambda)[0] = 1.0 - (*lambda)[1];
- k = -1;
+ int k = -1;
lmin = 0.0;
- j = 0;
- for (i = 0; i <= dim; i++) {
+ int j = 0;
+ for (int i = 0; i <= dim; i++) {
if ((*lambda)[i] < -1.E-5) {
if ((*lambda)[i] < lmin) {
k = i;
@@ -170,26 +157,11 @@ namespace AMDiS {
/****************************************************************************/
double ElInfo1d::getNormal(int side, WorldVector<double> &normal) const
{
- double det;
normal = coord_[side] - coord_[(side + 1) % 2];
- det = norm(&normal);
+ double det = norm(&normal);
TEST_EXIT(det > 1.e-30)("det = 0 on side %d\n", side);
- normal *= 1.0/det;
-
- // int dow = Global::getGeo(WORLD);
- // double det = 0.0;
-
- // if(dow == 1){
- // normal[0] = side ? 1.0 : -1.0;
- // det = 1.0;
- // }
- // else{
- // normal = side ? coord_[1] - coord_[0] : coord_[0] - coord_[1];
- // det = norm(&normal);
- // normal *= 1.0 / det;
- // det = 1.0;
- // }
-
+ normal *= 1.0 / det;
+
return(det);
}
@@ -201,8 +173,9 @@ namespace AMDiS {
/****************************************************************************/
double ElInfo1d::getElementNormal( WorldVector<double> &elementNormal) const
{
- FUNCNAME("ElInfo::getElementNormal");
- double det = 0.0;
+ FUNCNAME("ElInfo::getElementNormal()");
+
+ double det = 0.0;
int dow = Global::getGeo(WORLD);
TEST_EXIT(dow = 2)(" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!");
@@ -222,8 +195,8 @@ namespace AMDiS {
void ElInfo1d::fillElInfo(int ichild, const ElInfo *elinfo_old)
{
- FUNCNAME("ElInfo1d::fillElInfo");
- int i; //, j;
+ FUNCNAME("ElInfo1d::fillElInfo()");
+
Element *nb;
Element *elem = elinfo_old->element_;
@@ -231,96 +204,81 @@ namespace AMDiS {
TEST_EXIT((element_ = const_cast<Element*>( elem->getChild(ichild))))
("missing child %d?\n", ichild);
- macroElement_ = elinfo_old->macroElement_;
- fillFlag_ = elinfo_old->fillFlag_;
- parent_ = elem;
- level_ = elinfo_old->level_ + 1;
+ macroElement_ = elinfo_old->macroElement_;
+ fillFlag_ = elinfo_old->fillFlag_;
+ parent_ = elem;
+ level_ = elinfo_old->level_ + 1;
- //int dow = Global::getGeo(WORLD);
int neighbours = mesh_->getGeo(NEIGH);
if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
- fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA))
- {
- const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elinfo_old->coord_);
+ fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- coord_[ichild] = (*old_coord)[ichild];
- if (elem->getNewCoord(-1)) {
- coord_[1-ichild] = *(elem->getNewCoord());
- } else {
- coord_[1-ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
- }
+ const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elinfo_old->coord_);
+
+ coord_[ichild] = (*old_coord)[ichild];
+ if (elem->getNewCoord(-1)) {
+ coord_[1 - ichild] = *(elem->getNewCoord());
+ } else {
+ coord_[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5;
}
+ }
// if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
// det = calcGrdLambda(*Lambda);
// }
- if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- WorldVector<double> oppC;
-
- TEST_EXIT(fillFlag_.isSet(Mesh::FILL_COORDS))
- ("FILL_OPP_COORDS only with FILL_COORDS\n");
-
- for (i = 0; i < neighbours; i++)
- {
- if (i != ichild)
- {
- nb = const_cast<Element*>( elem->getChild(1-ichild));
- // if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- if ( fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- oppC = elinfo_old->coord_[i];
- }
- }
- else
- {
- nb = const_cast<Element*>( elinfo_old->getNeighbour(i));
-
- if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- oppC = elinfo_old->oppCoord_[i];
- }
- }
+ if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ WorldVector<double> oppC;
+
+ TEST_EXIT(fillFlag_.isSet(Mesh::FILL_COORDS))
+ ("FILL_OPP_COORDS only with FILL_COORDS\n");
+
+ for (int i = 0; i < neighbours; i++) {
+ if (i != ichild) {
+ nb = const_cast<Element*>( elem->getChild(1-ichild));
+ if ( fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ oppC = elinfo_old->coord_[i];
+ }
+ } else {
+ nb = const_cast<Element*>( elinfo_old->getNeighbour(i));
- if (nb)
- {
- while (nb->getChild(0)) // make nb nearest element
- {
- if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- if (nb->getNewCoord(-1)) {
- oppC = *(nb->getNewCoord());
- } else {
- oppC = (coord_[i] + oppC) * 0.5;
- }
- }
- nb = const_cast<Element*>( nb->getChild(1-i));
- }
-
- if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- oppCoord_[i] = oppC;
- }
- }
- neighbour_[i] = nb;
- oppVertex_[i] = nb ? i : -1;
+ if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ oppC = elinfo_old->oppCoord_[i];
}
- }
+ }
- if (fillFlag_.isSet(Mesh::FILL_BOUND))
- {
- boundary_[ichild] = elinfo_old->getBoundary(ichild);
- boundary_[1-ichild] = INTERIOR;
+ if (nb) {
+ while (nb->getChild(0)) { // make nb nearest element
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ if (nb->getNewCoord(-1)) {
+ oppC = *(nb->getNewCoord());
+ } else {
+ oppC = (coord_[i] + oppC) * 0.5;
+ }
+ }
+ nb = const_cast<Element*>( nb->getChild(1-i));
+ }
- if(elinfo_old->getProjection(0) &&
- elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION)
- {
- projection_[0] = elinfo_old->getProjection(0);
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ oppCoord_[i] = oppC;
}
+ }
+ neighbour_[i] = nb;
+ oppVertex_[i] = nb ? i : -1;
}
+ }
+ if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
+ boundary_[ichild] = elinfo_old->getBoundary(ichild);
+ boundary_[1-ichild] = INTERIOR;
+
+ if (elinfo_old->getProjection(0) &&
+ elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
+ projection_[0] = elinfo_old->getProjection(0);
+ }
+ }
+
return;
}
diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h
index bc46f1694e3919e2301d11f2959a243bfe01c24b..f2eff8106fc29e816e4b878e066f3f7916b9b327 100644
--- a/AMDiS/src/ElInfo1d.h
+++ b/AMDiS/src/ElInfo1d.h
@@ -75,7 +75,9 @@ namespace AMDiS {
*/
double getElementNormal( WorldVector<double> &normal) const;
- int getSideOfNeighbour(int i) { return (i+1)%2; };
+ int getSideOfNeighbour(int i) {
+ return (i + 1) % 2;
+ };
};
}
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 0370a1f3331e46f5832a9dabfb41357ea047523b..c0979903d0bca9023fef15eb6e09e4551e38d356 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -16,28 +16,19 @@ namespace AMDiS {
void ElInfo2d::fillMacroInfo(const MacroElement * mel)
{
- FUNCNAME("ElInfo::fillMacroInfo");
-
- int i, k;
- Element *nb;
- MacroElement *mnb;
-#if 1 // #if !NEIGH_IN_EL
- bool fill_opp_coords;
-#endif
-
-
+ FUNCNAME("ElInfo::fillMacroInfo()");
+
macroElement_ = const_cast<MacroElement*>(mel);
- element_ = const_cast<Element*>(mel->getElement());
- parent_ = NULL;
- level_ = 0;
-
- int vertices = mesh_->getGeo(VERTEX);
- //int edges = mesh_->getGeo(EDGE);
+ element_ = const_cast<Element*>(mel->getElement());
+ parent_ = NULL;
+ level_ = 0;
if (fillFlag_.isSet(Mesh::FILL_COORDS) ||
fillFlag_.isSet(Mesh::FILL_DET) ||
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- for (i=0; i<vertices; i++) {
+
+ int vertices = mesh_->getGeo(VERTEX);
+ for (int i = 0; i < vertices; i++) {
coord_[i] = mel->coord[i];
}
}
@@ -48,51 +39,89 @@ namespace AMDiS {
int neighbours = mesh_->getGeo(NEIGH);
- if(fillFlag_.isSet(Mesh::FILL_OPP_COORDS) ||
- fillFlag_.isSet(Mesh::FILL_NEIGH))
- {
- fill_opp_coords = (fillFlag_.isSet(Mesh::FILL_OPP_COORDS));
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) ||
+ fillFlag_.isSet(Mesh::FILL_NEIGH)) {
+
+ bool fill_opp_coords = (fillFlag_.isSet(Mesh::FILL_OPP_COORDS));
- for (i = 0; i < neighbours; i++) {
- if ((mnb = mel->getNeighbour(i))) {
- neighbour_[i] = mel->getNeighbour(i)->getElement();
- nb = const_cast<Element*>(neighbour_[i]);
- k = oppVertex_[i] = mel->getOppVertex(i);
- if (nb->getFirstChild() && (k != 2)) { // make nb nearest el.
- nb = neighbour_[i] =
- (k==0) ?
- nb->getSecondChild() :
- nb->getFirstChild();
- k = oppVertex_[i] = 2;
- if (fill_opp_coords) {
- if (nb->getNewCoord(-1)) {
- oppCoord_[i] = *(nb->getNewCoord());
+ for (int i = 0; i < neighbours; i++) {
+ MacroElement *macroNeighbour = mel->getNeighbour(i);
+
+ if (macroNeighbour) {
+ neighbour_[i] = macroNeighbour->getElement();
+
+ Element *nb = const_cast<Element*>(neighbour_[i]);
+
+ int edgeNo = oppVertex_[i] = mel->getOppVertex(i);
+ if (nb->getFirstChild() && (edgeNo != 2)) { // make nb nearest el.
+ if (edgeNo == 0) {
+ nb = neighbour_[i] = nb->getSecondChild();
+ } else {
+ nb = neighbour_[i] = nb->getFirstChild();
+ }
+
+ oppVertex_[i] = 2;
+
+ if (fill_opp_coords) {
+ if (nb->getNewCoord(-1)) {
+ oppCoord_[i] = *(nb->getNewCoord());
+ } else {
+ oppCoord_[i] = (macroNeighbour->coord[0] + macroNeighbour->coord[1]) * 0.5;
+ }
+
+ switch (i) {
+ case 0:
+ if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
+ neighbourCoord_[i][0] = macroNeighbour->coord[2];
+ neighbourCoord_[i][1] = macroNeighbour->coord[0];
+ } else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(1))) {
+ neighbourCoord_[i][0] = macroNeighbour->coord[1];
+ neighbourCoord_[i][1] = macroNeighbour->coord[2];
} else {
- oppCoord_[i] =
- (mnb->coord[0] + mnb->coord[1]) * 0.5;
+ ERROR_EXIT("should not happen!\n");
}
- }
- } else {
- if (fill_opp_coords) {
- oppCoord_[i] = mnb->coord[k];
+
+ neighbourCoord_[i][2] = oppCoord_[i];
+ break;
+
+ case 1:
+ if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(2))) {
+ neighbourCoord_[i][0] = macroNeighbour->coord[1];
+ neighbourCoord_[i][1] = macroNeighbour->coord[2];
+ } else if (*(macroNeighbour->getElement()->getDOF(2)) == *(element_->getDOF(0))) {
+ neighbourCoord_[i][0] = macroNeighbour->coord[2];
+ neighbourCoord_[i][1] = macroNeighbour->coord[0];
+ } else {
+ ERROR_EXIT("should not happen!\n");
+ }
+
+ neighbourCoord_[i][2] = oppCoord_[i];
+ break;
+
+ default:
+ ERROR_EXIT("should not happen!\n");
+ break;
}
}
} else {
- neighbour_[i] = NULL;
+ if (fill_opp_coords) {
+ oppCoord_[i] = macroNeighbour->coord[edgeNo];
+
+ neighbourCoord_[i] = macroNeighbour->coord;
+ }
}
- }
+ } else {
+ neighbour_[i] = NULL;
+ }
}
-
- if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
+ }
- // for (i = 0; i < vertices; i++)
- // (*bound)[i] = mel->getBound(i);
-
- for (i = 0; i < element_->getGeo(BOUNDARY) ; i++) {
+ if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
+ for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
boundary_[i] = mel->getBoundary(i);
}
- for(i = 0; i < element_->getGeo(PROJECTION); i++) {
+ for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
projection_[i] = mel->getProjection(i);
}
}
@@ -105,191 +134,257 @@ namespace AMDiS {
void ElInfo2d::fillElInfo(int ichild, const ElInfo *elinfo_old)
{
- FUNCNAME("ElInfo::fillElInfo");
- Element *elem = elinfo_old->element_;
- int j;
- Element *nb;
+ FUNCNAME("ElInfo::fillElInfo()");
- Flag fill_flag = elinfo_old->fillFlag_;
- bool fill_opp_coords;
+ Element *elem = elinfo_old->element_;
+ Element *nb;
+
+ Flag fill_flag = elinfo_old->fillFlag_;
TEST_EXIT(elem->getFirstChild())("no children?\n");
- TEST_EXIT(element_ = const_cast<Element*>((ichild==0) ?
+ TEST_EXIT(element_ = const_cast<Element*>((ichild == 0) ?
elem->getFirstChild() :
elem->getSecondChild()))
("missing child %d?\n", ichild);
macroElement_ = elinfo_old->macroElement_;
fillFlag_ = fill_flag;
- parent_ = elem;
- level_ = elinfo_old->level_ + 1;
+ parent_ = elem;
+ level_ = elinfo_old->level_ + 1;
int dow = Global::getGeo(WORLD);
if (fillFlag_.isSet(Mesh::FILL_COORDS) ||
fillFlag_.isSet(Mesh::FILL_DET) ||
- fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA))
- {
- if (elem->getNewCoord(-1)) {
- coord_[2] = *(elem->getNewCoord());
- } else {
- for (j = 0; j < dow; j++)
- coord_[2][j] =
- 0.5 * (elinfo_old->coord_[0][j] + elinfo_old->coord_[1][j]);
- // projection !?!
+ fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
+
+ if (elem->getNewCoord(-1)) {
+ coord_[2] = *(elem->getNewCoord());
+ } else {
+ coord_[2].setMidpoint(elinfo_old->coord_[0], elinfo_old->coord_[1]);
+ }
+
+ if (ichild == 0) {
+ coord_[0] = elinfo_old->coord_[2];
+ coord_[1] = elinfo_old->coord_[0];
+ } else {
+ coord_[0] = elinfo_old->coord_[1];
+ coord_[1] = elinfo_old->coord_[2];
+ }
+ }
+
+ bool fill_opp_coords = (fill_flag.isSet(Mesh::FILL_OPP_COORDS));
+
+ if (fill_flag.isSet(Mesh::FILL_NEIGH) || fill_opp_coords) {
+ if (ichild == 0) {
+ // Calculation of the neighbour 2, its oppCoords and the
+ // cooresponding oppVertex.
+
+ neighbour_[2] = elinfo_old->neighbour_[1];
+ oppVertex_[2] = elinfo_old->oppVertex_[1];
+
+ if (neighbour_[2] && fill_opp_coords) {
+ oppCoord_[2] = elinfo_old->oppCoord_[1];
+ neighbourCoord_[2] = elinfo_old->neighbourCoord_[1];
}
+
+
+ // Calculation of the neighbour 1, its oppCoords and the
+ // cooresponding oppVertex.
+
+ if (elem->getFirstChild() &&
+ elem->getSecondChild()->getFirstChild() &&
+ elem->getSecondChild()->getFirstChild()) {
+
+ neighbour_[1] = elem->getSecondChild()->getSecondChild();
+ oppVertex_[1] = 2;
+
+ if (fill_opp_coords) {
+ if (elem->getSecondChild()->getNewCoord(-1)) {
+ oppCoord_[1] = *(elem->getSecondChild()->getNewCoord());
+ } else {
+ oppCoord_[1].setMidpoint(elinfo_old->coord_[1],
+ elinfo_old->coord_[2]);
+ }
- if (ichild==0) {
- for (j = 0; j < dow; j++) {
- coord_[0][j] = elinfo_old->coord_[2][j];
- coord_[1][j] = elinfo_old->coord_[0][j];
+ neighbourCoord_[1][0] = coord_[0];
+ neighbourCoord_[1][1] = coord_[2];
+ neighbourCoord_[1][2] = oppCoord_[1];
}
} else {
- for (j = 0; j < dow; j++) {
- coord_[0][j] = elinfo_old->coord_[1][j];
- coord_[1][j] = elinfo_old->coord_[2][j];
+ neighbour_[1] = elem->getSecondChild();
+ oppVertex_[1] = 0;
+
+ if (fill_opp_coords) {
+ oppCoord_[1] = elinfo_old->coord_[1];
+
+ neighbourCoord_[1][0] = elinfo_old->coord_[1];
+ neighbourCoord_[1][1] = elinfo_old->coord_[2];
+ neighbourCoord_[1][2] = coord_[2];
}
}
- }
- // if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- // det = calcGrdLambda(*Lambda);
- // }
- if (fill_flag.isSet(Mesh::FILL_NEIGH)
- || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- fill_opp_coords = (fill_flag.isSet(Mesh::FILL_OPP_COORDS));
- if (ichild==0) {
- if ((neighbour_[2] = elinfo_old->neighbour_[1])) {
+ // Calculation of the neighbour 0, its oppCoords and the
+ // cooresponding oppVertex.
+
+ nb = elinfo_old->neighbour_[2];
+ if (nb) {
+ TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n");
+ TEST_EXIT(nb->getFirstChild())("missing first child?\n");
+ TEST_EXIT(nb->getSecondChild())("missing second child?\n");
+
+ nb = nb->getSecondChild();
+
+ if (nb->getFirstChild()) {
+ oppVertex_[0] = 2;
+
if (fill_opp_coords) {
- for (j=0; j<dow; j++)
- oppCoord_[2][j] = elinfo_old->oppCoord_[1][j];
- }
- }
-
- oppVertex_[2] = elinfo_old->oppVertex_[1];
-
- if (elem->getFirstChild() &&
- elem->getSecondChild()->getFirstChild() &&
- elem->getSecondChild()->getFirstChild()) // ????
- {
- TEST_EXIT((neighbour_[1] = elem->getSecondChild()->getSecondChild()))
- ("???");
- oppVertex_[1] = 2;
- if (fill_opp_coords) {
- if (elem->getSecondChild()->getNewCoord(-1))
- oppCoord_[1] = *(elem->getSecondChild()->getNewCoord());
- else
- for (j=0; j<dow; j++)
- oppCoord_[1][j] =
- 0.5 * (elinfo_old->coord_[1][j] +
- elinfo_old->coord_[2][j]);
- // projection !?!
+ if (nb->getNewCoord(-1)) {
+ oppCoord_[0] = *(nb->getNewCoord());
+ } else {
+ oppCoord_[0].setMidpoint(elinfo_old->neighbourCoord_[2][1],
+ elinfo_old->neighbourCoord_[2][2]);
}
- } else {
- TEST_EXIT((neighbour_[1] = elem->getSecondChild()))("???");
- oppVertex_[1] = 0;
+
+ neighbourCoord_[0][0].setMidpoint(elinfo_old->neighbourCoord_[2][0],
+ elinfo_old->neighbourCoord_[2][1]);
+ neighbourCoord_[0][1] = elinfo_old->neighbourCoord_[2][1];
+ neighbourCoord_[0][2] = oppCoord_[0];
+ }
+
+ nb = nb->getFirstChild();
+ } else {
+ oppVertex_[0] = 1;
+
if (fill_opp_coords) {
- for (j=0; j<dow; j++)
- oppCoord_[1][j] = elinfo_old->coord_[1][j];
+ oppCoord_[0] = elinfo_old->oppCoord_[2];
+
+ neighbourCoord_[0][0] = elinfo_old->neighbourCoord_[2][0];
+ neighbourCoord_[0][1] = elinfo_old->neighbourCoord_[2][2];
+ neighbourCoord_[0][2].setMidpoint(elinfo_old->neighbourCoord_[2][0],
+ elinfo_old->neighbourCoord_[2][1]);
}
}
- if ((nb = elinfo_old->neighbour_[2])) {
- TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n");
- TEST_EXIT(nb->getFirstChild())("missing children?\n");
- TEST_EXIT((nb = nb->getSecondChild()))("missing getSecondChild()?\n");
-
- if (nb->getFirstChild()) {
- oppVertex_[0] = 2;
- if (fill_opp_coords) {
- if (nb->getNewCoord(-1)) {
- oppCoord_[0] = *(nb->getNewCoord());
- } else {
- for (j=0; j<dow; j++)
- oppCoord_[0][j] =
- 0.5*(elinfo_old->oppCoord_[2][j] +
- elinfo_old->coord_[0][j]);
- }
- }
- nb = nb->getFirstChild();
+ }
+
+ neighbour_[0] = nb;
+ } else { /* ichild == 1 */
+ // Calculation of the neighbour 2, its oppCoords and the
+ // cooresponding oppVertex.
+
+ neighbour_[2] = elinfo_old->neighbour_[0];
+ oppVertex_[2] = elinfo_old->oppVertex_[0];
+
+ if (neighbour_[2] && fill_opp_coords) {
+ oppCoord_[2] = elinfo_old->oppCoord_[0];
+ neighbourCoord_[2] = elinfo_old->neighbourCoord_[0];
+ }
+
+
+ // Calculation of the neighbour 0, its oppCoords and the
+ // cooresponding oppVertex.
+
+ if (elem->getFirstChild()->getFirstChild()) {
+ neighbour_[0] = elem->getFirstChild()->getFirstChild();
+ oppVertex_[0] = 2;
+
+ if (fill_opp_coords) {
+ if (elem->getFirstChild()->getNewCoord(-1)) {
+ oppCoord_[0] = *(elem->getFirstChild()->getNewCoord());
} else {
- oppVertex_[0] = 1;
- if (fill_opp_coords) {
- for (j=0; j<dow; j++)
- oppCoord_[0][j] = elinfo_old->oppCoord_[2][j];
- }
+ oppCoord_[0].setMidpoint(elinfo_old->coord_[0],
+ elinfo_old->coord_[2]);
}
+
+ neighbourCoord_[0][0] = coord_[2];
+ neighbourCoord_[0][1] = coord_[1];
+ neighbourCoord_[0][2] = oppCoord_[0];
}
-
- neighbour_[0] = nb;
- } else { /* ichild==1 */
- if ((neighbour_[2] = elinfo_old->neighbour_[0])) {
- if (fill_opp_coords) {
- for (j=0; j<dow; j++)
- oppCoord_[2][j] = elinfo_old->oppCoord_[0][j];
- }
+ } else {
+ neighbour_[0] = elem->getFirstChild();
+ oppVertex_[0] = 1;
+
+ if (fill_opp_coords) {
+ oppCoord_[0] = elinfo_old->coord_[0];
+
+ neighbourCoord_[0][0] = elinfo_old->coord_[2];
+ neighbourCoord_[0][1] = elinfo_old->coord_[0];
+ neighbourCoord_[0][2] = coord_[2];
}
-
- oppVertex_[2] = elinfo_old->oppVertex_[0];
- if (elem->getFirstChild()->getFirstChild()) {
- neighbour_[0] = elem->getFirstChild()->getFirstChild();
- oppVertex_[0] = 2;
+ }
+
+ // Calculation of the neighbour 1, its oppCoords and the
+ // cooresponding oppVertex.
+
+ nb = elinfo_old->neighbour_[2];
+ if (nb) {
+ TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n");
+ TEST((nb = nb->getFirstChild()))("missing child?\n");
+
+ if (nb->getFirstChild()) {
+ oppVertex_[1] = 2;
+
if (fill_opp_coords) {
- if (elem->getFirstChild()->getNewCoord(-1)) {
- oppCoord_[0] = *(elem->getFirstChild()->getNewCoord());
+ if (nb->getNewCoord(-1)) {
+ oppCoord_[1] = *(nb->getNewCoord());
} else {
- for (j=0; j<dow; j++)
- oppCoord_[0][j] =
- 0.5 * (elinfo_old->coord_[0][j] +
- elinfo_old->coord_[2][j]);
- // projection !?!
+ WorldVector<double> d1;
+ for (int j=0; j<dow; j++)
+ d1[j] =
+ 0.5*(elinfo_old->oppCoord_[2][j] +
+ elinfo_old->coord_[1][j]);
+
+ oppCoord_[1].setMidpoint(elinfo_old->neighbourCoord_[2][0],
+ elinfo_old->neighbourCoord_[2][2]);
+
+ if (!(d1 == oppCoord_[1])) {
+ if ((elinfo_old->coord_[0][0] == 0.0) ||
+ (elinfo_old->coord_[1][0] == 0.0) ||
+ (elinfo_old->coord_[2][0] == 0.0) ||
+ (elinfo_old->coord_[0][0] == 80.0) ||
+ (elinfo_old->coord_[1][0] == 80.0) ||
+ (elinfo_old->coord_[2][0] == 80.0)) {
+
+ } else {
+ ERROR_EXIT("ja 4\n");
+ }
+ }
}
+
+ neighbourCoord_[1][0] = elinfo_old->neighbourCoord_[2][0];
+ neighbourCoord_[1][1].setMidpoint(elinfo_old->neighbourCoord_[2][0],
+ elinfo_old->neighbourCoord_[2][1]);
+ neighbourCoord_[1][2] = oppCoord_[1];
}
+ nb = nb->getSecondChild();
+
} else {
- neighbour_[0] = elem->getFirstChild();
- oppVertex_[0] = 1;
+ oppVertex_[1] = 0;
+
if (fill_opp_coords) {
- for (j=0; j<dow; j++)
- oppCoord_[0][j] = elinfo_old->coord_[0][j];
- }
- }
- if ((nb = elinfo_old->neighbour_[2])) {
- TEST(elinfo_old->oppVertex_[2] == 2)("invalid neighbour\n");
- TEST((nb = nb->getFirstChild()))("missing child?\n");
- if (nb->getFirstChild()) {
- oppVertex_[1] = 2;
- if (fill_opp_coords) {
- if (nb->getNewCoord(-1))
- oppCoord_[1] = *(nb->getNewCoord());
- else
- for (j=0; j<dow; j++)
- oppCoord_[1][j] =
- 0.5*(elinfo_old->oppCoord_[2][j] +
- elinfo_old->coord_[1][j]);
- // projection !?!
- }
- nb = nb->getSecondChild();
- } else {
- oppVertex_[1] = 0;
- if (fill_opp_coords) {
- for (j=0; j<dow; j++)
- oppCoord_[1][j] = elinfo_old->oppCoord_[2][j];
- }
+ oppCoord_[1] = elinfo_old->oppCoord_[2];
+
+ neighbourCoord_[1][0] = elinfo_old->neighbourCoord_[2][2];
+ neighbourCoord_[1][1] = elinfo_old->neighbourCoord_[2][0];
+ neighbourCoord_[1][2].setMidpoint(elinfo_old->neighbourCoord_[2][0],
+ elinfo_old->neighbourCoord_[2][1]);
}
}
- neighbour_[1] = nb;
}
- }
-
+ neighbour_[1] = nb;
+ } // if (ichild == 0) {} else
+ } // if (fill_flag.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
+
+
if (fill_flag.isSet(Mesh::FILL_BOUND)) {
- if (elinfo_old->getBoundary(2))
+ if (elinfo_old->getBoundary(2)) {
boundary_[5] = elinfo_old->getBoundary(2);
- else
+ } else {
boundary_[5] = INTERIOR;
-
+ }
- if (ichild==0) {
+ if (ichild == 0) {
boundary_[3] = elinfo_old->getBoundary(5);
boundary_[4] = elinfo_old->getBoundary(3);
boundary_[0] = elinfo_old->getBoundary(2);
@@ -303,12 +398,12 @@ namespace AMDiS {
boundary_[2] = elinfo_old->getBoundary(0);
}
- if(elinfo_old->getProjection(0) &&
- elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION)
- {
- projection_[0] = elinfo_old->getProjection(0);
- } else { // boundary projection
- if (ichild==0) {
+ if (elinfo_old->getProjection(0) &&
+ elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
+
+ projection_[0] = elinfo_old->getProjection(0);
+ } else { // boundary projection
+ if (ichild == 0) {
projection_[0] = elinfo_old->getProjection(2);
projection_[1] = NULL;
projection_[2] = elinfo_old->getProjection(1);
@@ -327,46 +422,42 @@ namespace AMDiS {
// TEST_EXIT(Global::getGeo(WORLD) == 2)
// ("dim != dim_of_world ! use parametric elements!\n");
- int i, j;
WorldVector<double> e1, e2;
- double sdet, det1, adet = 0.0;
- const WorldVector<double> v0= coord_[0];
- double a11, a12, a21, a22;
+ double adet = 0.0;
+ const WorldVector<double> v0 = coord_[0];
testFlag(Mesh::FILL_COORDS);
int dow = Global::getGeo(WORLD);
int dim = mesh_->getDim();
- for (i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++) {
e1[i] = coord_[1][i] - v0[i];
e2[i] = coord_[2][i] - v0[i];
}
- if(Global::getGeo(WORLD) == 2) {
- sdet = e1[0] * e2[1] - e1[1] * e2[0];
+ if (dow == 2) {
+ double sdet = e1[0] * e2[1] - e1[1] * e2[0];
adet = abs(sdet);
if (adet < 1.0E-25) {
MSG("abs(det) = %f\n", adet);
- for (i = 0; i <= dim; i++)
- for (j = 0; j < dow; j++)
+ for (int i = 0; i <= dim; i++)
+ for (int j = 0; j < dow; j++)
grd_lam[i][j] = 0.0;
+ } else {
+ double det1 = 1.0 / sdet;
+ double a11 = e2[1] * det1; /* (a_ij) = A^{-T} */
+ double a21 = -e2[0] * det1;
+ double a12 = -e1[1] * det1;
+ double a22 = e1[0] * det1;
+
+ grd_lam[1][0] = a11;
+ grd_lam[1][1] = a21;
+ grd_lam[2][0] = a12;
+ grd_lam[2][1] = a22;
+ grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
+ grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
}
- else
- {
- det1 = 1.0 / sdet;
- a11 = e2[1] * det1; /* (a_ij) = A^{-T} */
- a21 = -e2[0] * det1;
- a12 = -e1[1] * det1;
- a22 = e1[0] * det1;
-
- grd_lam[1][0] = a11;
- grd_lam[1][1] = a21;
- grd_lam[2][0] = a12;
- grd_lam[2][1] = a22;
- grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
- grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
- }
} else {
WorldVector<double> normal;
@@ -376,16 +467,16 @@ namespace AMDiS {
if (adet < 1.0E-15) {
MSG("abs(det) = %lf\n", adet);
- for (i = 0; i <= dim; i++)
- for (j = 0; j < dow; j++)
+ for (int i = 0; i <= dim; i++)
+ for (int j = 0; j < dow; j++)
grd_lam[i][j] = 0.0;
} else {
vectorProduct(e2, normal, grd_lam[1]);
vectorProduct(normal, e1, grd_lam[2]);
- double adet2 = 1.0/(adet*adet);
+ double adet2 = 1.0 / (adet * adet);
- for (i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++) {
grd_lam[1][i] *= adet2;
grd_lam[2][i] *= adet2;
}
@@ -402,33 +493,32 @@ namespace AMDiS {
const int ElInfo2d::worldToCoord(const WorldVector<double>& xy,
DimVec<double>* lambda) const
{
- FUNCNAME("ElInfo::worldToCoord");
- DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
- WorldVector<double> x;
- double x0, det, det0, det1;
- int i, j;
-
- static DimVec<double> vec(mesh_->getDim(), NO_INIT);
+ FUNCNAME("ElInfo::worldToCoord()");
TEST_EXIT(lambda)("lambda must not be NULL\n");
+ DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
+ WorldVector<double> x;
+ static DimVec<double> vec(mesh_->getDim(), NO_INIT);
+
int dim = mesh_->getDim();
int dimOfWorld = Global::getGeo(WORLD);
- for (j=0; j<dimOfWorld; j++) {
- x0 = coord_[dim][j];
+ for (int j = 0; j < dimOfWorld; j++) {
+ double x0 = coord_[dim][j];
x[j] = xy[j] - x0;
- for (i=0; i < dim; i++)
+ for (int i = 0; i < dim; i++)
edge[i][j] = coord_[i][j] - x0;
}
- det = edge[0][0] * edge[1][1] - edge[0][1] * edge[1][0];
- det0 = x[0] * edge[1][1] - x[1] * edge[1][0];
- det1 = edge[0][0] * x[1] - edge[0][1] * x[0];
+ double det = edge[0][0] * edge[1][1] - edge[0][1] * edge[1][0];
+ double det0 = x[0] * edge[1][1] - x[1] * edge[1][0];
+ double det1 = edge[0][0] * x[1] - edge[0][1] * x[0];
if (abs(det) < DBL_TOL) {
ERROR("det = %le; abort\n", det);
- for (i=0; i<=dim; i++) (*lambda)[i] = 1.0/dim;
+ for (int i = 0; i <= dim; i++)
+ (*lambda)[i] = 1.0/dim;
return 0;
}
@@ -436,22 +526,14 @@ namespace AMDiS {
(*lambda)[1] = det1 / det;
(*lambda)[2] = 1.0 - (*lambda)[0] - (*lambda)[1];
- // if(dimOfWorld == 3) {
- // if(abs(x[2] - edge[0][2] * lambda[0] - edge[1][2] * lambda[1]) > 1.e-5) {
- // WARNING("system not solvable\n");
- // };
- // }
-
int k = -1;
double lmin = 0.0;
- j = 0;
- for (i = 0; i <= dim; i++) {
+ for (int i = 0; i <= dim; i++) {
if ((*lambda)[i] < -1.E-5) {
if ((*lambda)[i] < lmin) {
k = i;
lmin = (*lambda)[i];
}
- j++;
}
}
@@ -461,30 +543,29 @@ namespace AMDiS {
double ElInfo2d::getNormal(int side, WorldVector<double> &normal) const
{
- FUNCNAME("ElInfo::getNormal");
- double det = 0.0;
- int i0, i1;
+ FUNCNAME("ElInfo::getNormal()");
- int dow = Global::getGeo(WORLD);
-
- i0 = (side+1) % 3;
- i1 = (side+2) % 3;
+ int i0 = (side + 1) % 3;
+ int i1 = (side + 2) % 3;
- if (dow == 2){
+ if (Global::getGeo(WORLD) == 2){
normal[0] = coord_[i1][1] - coord_[i0][1];
normal[1] = coord_[i0][0] - coord_[i1][0];
} else { // dow == 3
WorldVector<double> e0, e1,e2, elementNormal;
- e0 = coord_[i1]; e0 -= coord_[i0];
- e1 = coord_[i1]; e1 -= coord_[side];
- e2 = coord_[i0]; e2 -= coord_[side];
+ e0 = coord_[i1];
+ e0 -= coord_[i0];
+ e1 = coord_[i1];
+ e1 -= coord_[side];
+ e2 = coord_[i0];
+ e2 -= coord_[side];
vectorProduct(e1, e2, elementNormal);
vectorProduct(elementNormal, e0, normal);
}
- det = norm(&normal);
+ double det = norm(&normal);
TEST_EXIT(det > 1.e-30)("det = 0 on face %d\n", side);
@@ -501,19 +582,16 @@ namespace AMDiS {
/****************************************************************************/
double ElInfo2d::getElementNormal( WorldVector<double> &elementNormal) const
{
- FUNCNAME("ElInfo::getElementNormal");
- double det = 0.0;
- int dow = Global::getGeo(WORLD);
- WorldVector<double> e0, e1;
+ FUNCNAME("ElInfo::getElementNormal()");
- TEST_EXIT(dow = 3)(" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!");
+ TEST_EXIT(Global::getGeo(WORLD) == 3)(" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!");
- e0 = coord_[1] - coord_[0];
- e1 = coord_[2] - coord_[0];
+ WorldVector<double> e0 = coord_[1] - coord_[0];
+ WorldVector<double> e1 = coord_[2] - coord_[0];
vectorProduct(e0, e1, elementNormal);
- det = norm(&elementNormal);
+ double det = norm(&elementNormal);
TEST_EXIT(det > 1.e-30)("det = 0");
@@ -521,13 +599,4 @@ namespace AMDiS {
return(det);
}
-
-
-
-
-
-
-
-
-
}
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 5796a62e8484b451a6f69b45f4114c611e0f72f1..138e7f04bde0f530d378ca7f94b9f56b849ff530 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -19,21 +19,20 @@ namespace AMDiS {
FUNCNAME("ElInfo3d::fillMacroInfo");
Element *nb;
MacroElement *mnb;
- int i, k;
Flag fill_opp_coords;
macroElement_ = const_cast<MacroElement*>( mel);
- element_ = const_cast<Element*>( mel->getElement());
- parent_ = NULL;
- level_ = 0;
- el_type = const_cast<MacroElement*>(mel)->getElType();
+ element_ = const_cast<Element*>( mel->getElement());
+ parent_ = NULL;
+ level_ = 0;
+ el_type = const_cast<MacroElement*>(mel)->getElType();
int vertices = mesh_->getGeo(VERTEX);
if (fillFlag_.isSet(Mesh::FILL_COORDS) ||
fillFlag_.isSet(Mesh::FILL_DET) ||
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
- for (i=0; i<vertices; i++) {
+ for (int i = 0; i < vertices; i++) {
coord_[i] = mel->coord[i];
}
}
@@ -44,41 +43,39 @@ namespace AMDiS {
int neighbours = mesh_->getGeo(NEIGH);
- if(fillFlag_.isSet(Mesh::FILL_OPP_COORDS) ||
- fillFlag_.isSet(Mesh::FILL_NEIGH)) {
+ if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS) ||
+ fillFlag_.isSet(Mesh::FILL_NEIGH)) {
fill_opp_coords.setFlags(fillFlag_ & Mesh::FILL_OPP_COORDS);
- for (i=0; i < neighbours; i++) {
+ for (int i = 0; i < neighbours; i++) {
if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
neighbour_[i] = const_cast<Element*>( mel->getNeighbour(i)->getElement());
nb = const_cast<Element*>( neighbour_[i]);
- k = oppVertex_[i] = mel->getOppVertex(i);
+ int k;
+ k = oppVertex_[i] = mel->getOppVertex(i);
if (nb->getChild(0) && (k < 2)) { /*make nb nearest element.*/
- if (k==1) {
+ if (k == 1) {
neighbour_[i] = const_cast<Element*>( nb->getChild(0));
nb = const_cast<Element*>( neighbour_[i]);
} else {
neighbour_[i] = const_cast<Element*>( nb->getChild(1));
nb = const_cast<Element*>( neighbour_[i]);
}
- k = oppVertex_[i] = 3;
- if (fill_opp_coords.isAnySet())
- {
- /* always edge between vertices 0 and 1 is bisected! */
- if (mnb->getElement()->isNewCoordSet())
- oppCoord_[i] = *(mnb->getElement()->getNewCoord());
- else
- oppCoord_[i] = (mnb->coord[0] + mnb->coord[1]) * 0.5;
- }
- }
- else {
+ k = oppVertex_[i] = 3;
+ if (fill_opp_coords.isAnySet()) {
+ /* always edge between vertices 0 and 1 is bisected! */
+ if (mnb->getElement()->isNewCoordSet())
+ oppCoord_[i] = *(mnb->getElement()->getNewCoord());
+ else
+ oppCoord_[i] = (mnb->coord[0] + mnb->coord[1]) * 0.5;
+ }
+ } else {
if (fill_opp_coords.isAnySet()) {
oppCoord_[i] = mnb->coord[k];
}
}
- }
- else {
+ } else {
neighbour_[i] = NULL;
}
}
@@ -87,23 +84,22 @@ namespace AMDiS {
//int faces = mesh->getGeo(FACE);
//int edges = mesh->getGeo(EDGE);
- if (fillFlag_.isSet(Mesh::FILL_BOUND))
- {
- for (i = 0; i < element_->getGeo(BOUNDARY); i++) {
- boundary_[i] = mel->getBoundary(i);
- }
-
- for (i = 0; i < element_->getGeo(PROJECTION); i++) {
- projection_[i] = mel->getProjection(i);
- }
+ if (fillFlag_.isSet(Mesh::FILL_BOUND)) {
+ for (int i = 0; i < element_->getGeo(BOUNDARY); i++) {
+ boundary_[i] = mel->getBoundary(i);
+ }
+
+ for (int i = 0; i < element_->getGeo(PROJECTION); i++) {
+ projection_[i] = mel->getProjection(i);
}
+ }
if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) {
WorldVector<WorldVector<double> > a;
double s;
- for (i=0; i<3; i++) {
- a[i] = mel->coord[i+1];
+ for (int i = 0; i < 3; i++) {
+ a[i] = mel->coord[i + 1];
a[i] -= mel->coord[0];
}
@@ -120,11 +116,11 @@ namespace AMDiS {
double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam) const
{
- FUNCNAME("ElInfo3d::calcGrdLambda");
+ FUNCNAME("ElInfo3d::calcGrdLambda()");
+
TEST_EXIT(Global::getGeo(WORLD) == 3)
("dim != dim_of_world ! use parametric elements!\n");
- int i, j;
WorldVector<double> e1, e2, e3;
WorldVector<double> v0;
double det, adet;
@@ -133,52 +129,49 @@ namespace AMDiS {
testFlag(Mesh::FILL_COORDS);
v0 = coord_[0];
- for (i = 0; i < 3; i++) {
+ for (int i = 0; i < 3; i++) {
e1[i] = coord_[1][i] - v0[i];
e2[i] = coord_[2][i] - v0[i];
e3[i] = coord_[3][i] - v0[i];
}
- det = e1[0] * (e2[1]*e3[2] - e2[2]*e3[1])
- - e1[1] * (e2[0]*e3[2] - e2[2]*e3[0])
- + e1[2] * (e2[0]*e3[1] - e2[1]*e3[0]);
+ det = e1[0] * (e2[1] * e3[2] - e2[2] * e3[1])
+ - e1[1] * (e2[0] * e3[2] - e2[2] * e3[0])
+ + e1[2] * (e2[0] * e3[1] - e2[1] * e3[0]);
adet = abs(det);
- if (adet < 1.0E-25)
- {
- MSG("abs(det) = %f\n",adet);
- for (i = 0; i < 4; i++)
- for (j = 0; j < 3; j++)
- grd_lam[i][j] = 0.0;
- }
- else
- {
- det = 1.0 / det;
- a11 = (e2[1]*e3[2] - e2[2]*e3[1]) * det; /* (a_ij) = A^{-T} */
- a12 = (e2[2]*e3[0] - e2[0]*e3[2]) * det;
- a13 = (e2[0]*e3[1] - e2[1]*e3[0]) * det;
- a21 = (e1[2]*e3[1] - e1[1]*e3[2]) * det;
- a22 = (e1[0]*e3[2] - e1[2]*e3[0]) * det;
- a23 = (e1[1]*e3[0] - e1[0]*e3[1]) * det;
- a31 = (e1[1]*e2[2] - e1[2]*e2[1]) * det;
- a32 = (e1[2]*e2[0] - e1[0]*e2[2]) * det;
- a33 = (e1[0]*e2[1] - e1[1]*e2[0]) * det;
-
- grd_lam[1][0] = a11;
- grd_lam[1][1] = a12;
- grd_lam[1][2] = a13;
- grd_lam[2][0] = a21;
- grd_lam[2][1] = a22;
- grd_lam[2][2] = a23;
- grd_lam[3][0] = a31;
- grd_lam[3][1] = a32;
- grd_lam[3][2] = a33;
-
- grd_lam[0][0] = -grd_lam[1][0] -grd_lam[2][0] -grd_lam[3][0];
- grd_lam[0][1] = -grd_lam[1][1] -grd_lam[2][1] -grd_lam[3][1];
- grd_lam[0][2] = -grd_lam[1][2] -grd_lam[2][2] -grd_lam[3][2];
- }
+ if (adet < 1.0E-25) {
+ MSG("abs(det) = %f\n",adet);
+ for (int i = 0; i < 4; i++)
+ for (int j = 0; j < 3; j++)
+ grd_lam[i][j] = 0.0;
+ } else {
+ det = 1.0 / det;
+ a11 = (e2[1] * e3[2] - e2[2] * e3[1]) * det; /* (a_ij) = A^{-T} */
+ a12 = (e2[2] * e3[0] - e2[0] * e3[2]) * det;
+ a13 = (e2[0] * e3[1] - e2[1] * e3[0]) * det;
+ a21 = (e1[2] * e3[1] - e1[1] * e3[2]) * det;
+ a22 = (e1[0] * e3[2] - e1[2] * e3[0]) * det;
+ a23 = (e1[1] * e3[0] - e1[0] * e3[1]) * det;
+ a31 = (e1[1] * e2[2] - e1[2] * e2[1]) * det;
+ a32 = (e1[2] * e2[0] - e1[0] * e2[2]) * det;
+ a33 = (e1[0] * e2[1] - e1[1] * e2[0]) * det;
+
+ grd_lam[1][0] = a11;
+ grd_lam[1][1] = a12;
+ grd_lam[1][2] = a13;
+ grd_lam[2][0] = a21;
+ grd_lam[2][1] = a22;
+ grd_lam[2][2] = a23;
+ grd_lam[3][0] = a31;
+ grd_lam[3][1] = a32;
+ grd_lam[3][2] = a33;
+
+ grd_lam[0][0] = -grd_lam[1][0] -grd_lam[2][0] -grd_lam[3][0];
+ grd_lam[0][1] = -grd_lam[1][1] -grd_lam[2][1] -grd_lam[3][1];
+ grd_lam[0][2] = -grd_lam[1][2] -grd_lam[2][2] -grd_lam[3][2];
+ }
return adet;
}
@@ -186,11 +179,11 @@ namespace AMDiS {
const int ElInfo3d::worldToCoord(const WorldVector<double>& xy,
DimVec<double>* lambda) const
{
- FUNCNAME("ElInfo::worldToCoord");
+ FUNCNAME("ElInfo::worldToCoord()");
+
DimVec<WorldVector<double> > edge(mesh_->getDim(), NO_INIT);
WorldVector<double> x;
double x0, det, det0, det1, det2;
- int i, j;
static DimVec<double> vec(mesh_->getDim(), NO_INIT);
@@ -207,10 +200,11 @@ namespace AMDiS {
/* ( q1z q2z q3z) (lambda3) (qz) */
/* mit qi=pi-p3, q=xy-p3 */
- for (j=0; j<dimOfWorld; j++) {
+ for (int j = 0; j < dimOfWorld; j++) {
x0 = coord_[dim][j];
x[j] = xy[j] - x0;
- for (i=0; i<dim; i++)
+
+ for (int i = 0; i < dim; i++)
edge[i][j] = coord_[i][j] - x0;
}
@@ -241,7 +235,11 @@ namespace AMDiS {
if (abs(det) < DBL_TOL) {
ERROR("det = %le; abort\n", det);
- for (i=0; i<=dim; i++) (*lambda)[i] = 1.0/dim;
+
+ for (int i = 0; i <= dim; i++) {
+ (*lambda)[i] = 1.0 / dim;
+ }
+
return 0;
}
@@ -252,14 +250,13 @@ namespace AMDiS {
int k = -1;
double lmin = 0.0;
- j = 0;
- for (i = 0; i <= dim; i++) {
+
+ for (int i = 0; i <= dim; i++) {
if ((*lambda)[i] < -1.E-5) {
if ((*lambda)[i] < lmin) {
k = i;
lmin = (*lambda)[i];
}
- j++;
}
}
@@ -273,48 +270,47 @@ namespace AMDiS {
void ElInfo3d::update()
{
- FUNCNAME("ElInfo::update");
+ FUNCNAME("ElInfo::update()");
int neighbours = mesh_->getGeo(NEIGH);
int vertices = mesh_->getGeo(VERTEX);
int dow = Global::getGeo(WORLD);
- if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- Tetrahedron *nb;
- int ineigh, ov, j, k;
- Flag fill_opp_coords = fillFlag_ & Mesh::FILL_OPP_COORDS;
-
- for (ineigh=0; ineigh < neighbours; ineigh++) {
- if ((nb = dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighbour_[ineigh])))) {
- ov = oppVertex_[ineigh];
- if (ov < 2 && nb->getFirstChild()) {
- if (fill_opp_coords != Flag(0)) {
- k = -1;
- for (j=0; j<vertices; j++)
- if (element_->getDOF(j) == nb->getDOF(1-ov)) k = j;
-
- if(k == -1) {
- for (j=0; j<vertices; j++)
- if (mesh_->associated(element_->getDOF(j, 0),
- nb->getDOF(1-ov, 0)))
- k = j;
+ if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+ Tetrahedron *nb;
+ Flag fill_opp_coords = fillFlag_ & Mesh::FILL_OPP_COORDS;
+
+ for (int ineigh = 0; ineigh < neighbours; ineigh++) {
+ if ((nb = dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighbour_[ineigh])))) {
+ int ov = oppVertex_[ineigh];
+ if (ov < 2 && nb->getFirstChild()) {
+ if (fill_opp_coords != Flag(0)) {
+ int k = -1;
+ for (int j = 0; j < vertices; j++)
+ if (element_->getDOF(j) == nb->getDOF(1 - ov))
+ k = j;
+
+ if (k == -1) {
+ for (int j = 0; j < vertices; j++) {
+ if (mesh_->associated(element_->getDOF(j, 0), nb->getDOF(1 - ov, 0))) {
+ k = j;
+ }
}
- TEST_EXIT(k >= 0)("neighbour dof not found\n");
-
- if (nb->isNewCoordSet())
- oppCoord_[ineigh] = *(nb->getNewCoord());
- else
- for (j = 0; j < dow; j++)
- oppCoord_[ineigh][j] =
- (oppCoord_[ineigh][j] + coord_[k][j]) / 2;
}
- neighbour_[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
- oppVertex_[ineigh] = 3;
+ TEST_EXIT(k >= 0)("neighbour dof not found\n");
+
+ if (nb->isNewCoordSet())
+ oppCoord_[ineigh] = *(nb->getNewCoord());
+ else
+ for (int j = 0; j < dow; j++)
+ oppCoord_[ineigh][j] = (oppCoord_[ineigh][j] + coord_[k][j]) / 2;
}
+ neighbour_[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
+ oppVertex_[ineigh] = 3;
}
}
}
+ }
}
@@ -322,17 +318,16 @@ namespace AMDiS {
{
FUNCNAME("ElInfo3d::getNormal");
double det = 0.0;
- int i, i0, i1, i2;
WorldVector<double> e0, e1, e2;
int dow = Global::getGeo(WORLD);
- if(dow == 3) {
- i0 = (face+1)%4;
- i1 = (face+2)%4;
- i2 = (face+3)%4;
+ if (dow == 3) {
+ int i0 = (face + 1) % 4;
+ int i1 = (face + 2) % 4;
+ int i2 = (face + 3) % 4;
- for (i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++) {
e0[i] = coord_[i1][i] - coord_[i0][i];
e1[i] = coord_[i2][i] - coord_[i0][i];
e2[i] = coord_[face][i] - coord_[i0][i];
@@ -341,7 +336,7 @@ namespace AMDiS {
vectorProduct(e0, e1, normal);
if ((e2*normal) < 0.0)
- for (i = 0; i < dow; i++)
+ for (int i = 0; i < dow; i++)
normal[i] = -normal[i];
det = norm(&normal);
@@ -362,32 +357,31 @@ namespace AMDiS {
void ElInfo3d::fillElInfo(int ichild, const ElInfo *elinfo_old)
{
- FUNCNAME("ElInfo3d::fillElInfo");
- int i,j,k;
- int el_type_local = 0; /* el_type in {0,1,2} */
- int ochild = 0; /* index of other child = 1-ichild */
- int *cv = NULL; /* cv = child_vertex[el_type][ichild] */
+ FUNCNAME("ElInfo3d::fillElInfo()");
+
+ int el_type_local = 0; /* el_type in {0,1,2} */
+ int ochild = 0; /* index of other child = 1-ichild */
+ int *cv = NULL; /* cv = child_vertex[el_type][ichild] */
const int (*cvg)[4] = NULL; /* cvg = child_vertex[el_type] */
- int *ce; /* ce = child_edge[el_type][ichild] */
- int iedge;
+ int *ce; /* ce = child_edge[el_type][ichild] */
Element *nb, *nbk;
const FixVec<Element*, NEIGH> *neigh_old;
Element *el_old = elinfo_old->element_;
Flag fillFlag__local = elinfo_old->fillFlag_;
DegreeOfFreedom *dof;
- int ov = -1;
+ int ov = -1;
FixVec<Element*, NEIGH> *neigh_local;
- Flag fill_opp_coords;
+ Flag fill_opp_coords;
Mesh *mesh = elinfo_old->getMesh();
TEST_EXIT(el_old->getChild(0))("missing child?\n"); /* Kuni 22.08.96 */
- element_ = const_cast<Element*>( el_old->getChild(ichild));
- macroElement_ = elinfo_old->macroElement_;
+ element_ = const_cast<Element*>( el_old->getChild(ichild));
+ macroElement_ = elinfo_old->macroElement_;
fillFlag_ = fillFlag__local;
- parent_ = el_old;
- level_ = elinfo_old->level_ + 1;
- el_type = (( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->el_type + 1) % 3;
+ parent_ = el_old;
+ level_ = elinfo_old->level_ + 1;
+ el_type = (( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->el_type + 1) % 3;
TEST_EXIT(element_)("missing child %d?\n", ichild);
@@ -395,262 +389,235 @@ namespace AMDiS {
el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
cvg = Tetrahedron::childVertex[el_type_local];
cv = const_cast<int*>( cvg[ichild]);
- ochild = 1-ichild;
+ ochild = 1 - ichild;
}
int dow = Global::getGeo(WORLD);
- if(fillFlag__local.isSet(Mesh::FILL_COORDS) ||
- fillFlag_.isSet(Mesh::FILL_DET) ||
- fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA))
- {
- for (i=0; i<3; i++) {
- for (j = 0; j < dow; j++) {
- coord_[i][j] = elinfo_old->coord_[cv[i]][j];
- }
+ if (fillFlag__local.isSet(Mesh::FILL_COORDS) ||
+ fillFlag_.isSet(Mesh::FILL_DET) ||
+ fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < dow; j++) {
+ coord_[i][j] = elinfo_old->coord_[cv[i]][j];
}
- if (el_old->getNewCoord())
- coord_[3] = *(el_old->getNewCoord());
- else
- for (j = 0; j < dow; j++)
- coord_[3][j] =
- (elinfo_old->coord_[0][j] + elinfo_old->coord_[1][j]) / 2;
}
+ if (el_old->getNewCoord()) {
+ coord_[3] = *(el_old->getNewCoord());
+ } else {
+ for (int j = 0; j < dow; j++) {
+ coord_[3][j] = (elinfo_old->coord_[0][j] + elinfo_old->coord_[1][j]) / 2;
+ }
+ }
+ }
// if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
// det = calcGrdLambda(*Lambda);
// }
- if(fillFlag__local.isSet(Mesh::FILL_NEIGH) ||
- fillFlag_.isSet(Mesh::FILL_OPP_COORDS))
- {
- neigh_local = &neighbour_;
- neigh_old = &elinfo_old->neighbour_;
- //oppVertex__local = *oppVertex_;
- fill_opp_coords.setFlags(fillFlag__local & Mesh::FILL_OPP_COORDS);
-
- /*----- nb[0] is other child --------------------------------------------*/
-
- /* if (nb = el_old->child[ochild]) old version */
- if(el_old->getChild(0) &&
- (nb = const_cast<Element*>( el_old->getChild(ochild)))) /*Kuni 22.08.96*/
- {
- if (nb->getChild(0))
- { /* go down one level for direct neighbour */
- if (fill_opp_coords.isAnySet())
- {
- if (nb->getNewCoord())
- {
- oppCoord_[0]= *(nb->getNewCoord());
- }
- else
- {
- k = cvg[ochild][1];
- for (j = 0; j < dow; j++)
- oppCoord_[0][j] =
- (elinfo_old->coord_[ochild][j] +
- elinfo_old->coord_[k][j]) / 2;
- }
- }
- (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
- oppVertex_[0] = 3;
- }
- else {
- if (fill_opp_coords.isAnySet()) {
- for (j = 0; j < dow; j++) {
- oppCoord_[0][j] = elinfo_old->coord_[ochild][j];
- }
+ if (fillFlag__local.isSet(Mesh::FILL_NEIGH) ||
+ fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
+
+ neigh_local = &neighbour_;
+ neigh_old = &elinfo_old->neighbour_;
+ fill_opp_coords.setFlags(fillFlag__local & Mesh::FILL_OPP_COORDS);
+
+ /*----- nb[0] is other child --------------------------------------------*/
+
+ /* if (nb = el_old->child[ochild]) old version */
+ if (el_old->getChild(0) &&
+ (nb = const_cast<Element*>( el_old->getChild(ochild)))) {
+
+ if (nb->getChild(0)) { /* go down one level for direct neighbour */
+ if (fill_opp_coords.isAnySet()) {
+ if (nb->getNewCoord()) {
+ oppCoord_[0]= *(nb->getNewCoord());
+ } else {
+ int k = cvg[ochild][1];
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[0][j] = (elinfo_old->coord_[ochild][j] + elinfo_old->coord_[k][j]) / 2;
}
- (*neigh_local)[0] = nb;
- oppVertex_[0] = 0;
}
}
- else {
- ERROR_EXIT("no other child");
- (*neigh_local)[0] = NULL;
+ (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
+ oppVertex_[0] = 3;
+ } else {
+ if (fill_opp_coords.isAnySet()) {
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[0][j] = elinfo_old->coord_[ochild][j];
+ }
+ }
+ (*neigh_local)[0] = nb;
+ oppVertex_[0] = 0;
}
+ } else {
+ ERROR_EXIT("no other child");
+ (*neigh_local)[0] = NULL;
+ }
- /*----- nb[1],nb[2] are childs of old neighbours nb_old[cv[i]] ----------*/
-
- for (i=1; i<3; i++)
- {
- if ((nb = const_cast<Element*>( (*neigh_old)[cv[i]])))
- {
- TEST_EXIT(nb->getChild(0))("nonconforming triangulation\n");
-
- for (k=0; k<2; k++) /* look at both childs of old neighbour */
- {
- nbk = const_cast<Element*>( nb->getChild(k));
- if (nbk->getDOF(0) == el_old->getDOF(ichild)) {
- /* opp. vertex */
- dof = const_cast<int*>(nb->getDOF(elinfo_old->oppVertex_[cv[i]]));
-
- if (dof == nbk->getDOF(1))
- {
- ov = 1;
- if (nbk->getChild(0))
- {
- if (fill_opp_coords.isAnySet())
- {
- if (nbk->getNewCoord())
- oppCoord_[i] = *(nbk->getNewCoord());
- else
- for (j = 0; j < dow; j++)
- oppCoord_[i][j] =
- (elinfo_old->oppCoord_[cv[i]][j]
- + elinfo_old->coord_[ichild][j]) / 2;
- }
- (*neigh_local)[i] = nbk->getChild(0);
- oppVertex_[i] = 3;
- break;
- }
- }
- else
- {
- //TEST_EXIT(dof == nbk->getDOF(2))("opp_vertex not found\n");
- if(dof != nbk->getDOF(2)) {
- ov = -1; break;
- }
- ov = 2;
- }
-
- if (fill_opp_coords.isAnySet())
- {
- for (j = 0; j < dow; j++)
- {
- oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
- }
- }
- (*neigh_local)[i] = nbk;
- oppVertex_[i] = ov;
- break;
- }
+ /*----- nb[1],nb[2] are childs of old neighbours nb_old[cv[i]] ----------*/
+
+ for (int i = 1; i < 3; i++) {
+ if ((nb = const_cast<Element*>( (*neigh_old)[cv[i]]))) {
+ TEST_EXIT(nb->getChild(0))("nonconforming triangulation\n");
- } /* end for k */
- //TEST_EXIT(k<2)("child not found with vertex\n");
-
- // periodic ?
- if(k == 2 || ov == -1) {
- for (k=0; k<2; k++) /* look at both childs of old neighbour */
- {
- nbk = const_cast<Element*>( nb->getChild(k));
- if (nbk->getDOF(0) == el_old->getDOF(ichild) ||
- mesh->associated(nbk->getDOF(0, 0), el_old->getDOF(ichild, 0))) {
- /* opp. vertex */
- dof = const_cast<int*>(nb->getDOF(elinfo_old->oppVertex_[cv[i]]));
-
- if (dof == nbk->getDOF(1) ||
- mesh->associated(dof[0], nbk->getDOF(1, 0)))
- {
- ov = 1;
- if (nbk->getChild(0))
- {
- if (fill_opp_coords.isAnySet())
- {
- if (nbk->getNewCoord())
- oppCoord_[i] = *(nbk->getNewCoord());
- else
- for (j = 0; j < dow; j++)
- oppCoord_[i][j] =
- (elinfo_old->oppCoord_[cv[i]][j]
- + elinfo_old->coord_[ichild][j]) / 2;
- }
- (*neigh_local)[i] = nbk->getChild(0);
- oppVertex_[i] = 3;
- break;
- }
- }
- else
- {
- TEST_EXIT(dof == nbk->getDOF(2) ||
- mesh->associated(dof[0], nbk->getDOF(2, 0)))
- ("opp_vertex not found\n");
- ov = 2;
- }
-
- if (fill_opp_coords.isAnySet())
- {
- for (j = 0; j < dow; j++)
- {
- oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
- }
- }
- (*neigh_local)[i] = nbk;
- oppVertex_[i] = ov;
- break;
+ int k;
+ for (k = 0; k < 2; k++) { /* look at both childs of old neighbour */
+ nbk = const_cast<Element*>( nb->getChild(k));
+ if (nbk->getDOF(0) == el_old->getDOF(ichild)) {
+ /* opp. vertex */
+ dof = const_cast<int*>(nb->getDOF(elinfo_old->oppVertex_[cv[i]]));
+
+ if (dof == nbk->getDOF(1)) {
+ ov = 1;
+ if (nbk->getChild(0)) {
+ if (fill_opp_coords.isAnySet()) {
+ if (nbk->getNewCoord()) {
+ oppCoord_[i] = *(nbk->getNewCoord());
+ } else {
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[i][j] = (elinfo_old->oppCoord_[cv[i]][j] + elinfo_old->coord_[ichild][j]) / 2;
}
-
- } /* end for k */
- TEST_EXIT(k<2)("child not found with vertex\n");
+ }
+ }
+ (*neigh_local)[i] = nbk->getChild(0);
+ oppVertex_[i] = 3;
+ break;
}
+ } else {
+ //TEST_EXIT(dof == nbk->getDOF(2))("opp_vertex not found\n");
+ if (dof != nbk->getDOF(2)) {
+ ov = -1;
+ break;
+ }
+ ov = 2;
}
- else
- {
- (*neigh_local)[i] = NULL;
- }
- } /* end for i */
-
- /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
-
- if (((*neigh_local)[3] = (*neigh_old)[ochild]))
- {
- oppVertex_[3] = elinfo_old->oppVertex_[ochild];
- if (fill_opp_coords.isAnySet()) {
- for (j = 0; j < dow; j++) {
- oppCoord_[3][j] = elinfo_old->oppCoord_[ochild][j];
+ if (fill_opp_coords.isAnySet()) {
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
+ }
}
+ (*neigh_local)[i] = nbk;
+ oppVertex_[i] = ov;
+ break;
}
- }
- }
+
+ } /* end for k */
+
+ //TEST_EXIT(k<2)("child not found with vertex\n");
+
+ // periodic ?
+ if (k == 2 || ov == -1) {
+ for (k = 0; k < 2; k++) { /* look at both childs of old neighbour */
+ nbk = const_cast<Element*>( nb->getChild(k));
+ if (nbk->getDOF(0) == el_old->getDOF(ichild) ||
+ mesh->associated(nbk->getDOF(0, 0), el_old->getDOF(ichild, 0))) {
+ /* opp. vertex */
+ dof = const_cast<int*>(nb->getDOF(elinfo_old->oppVertex_[cv[i]]));
+
+ if (dof == nbk->getDOF(1) ||
+ mesh->associated(dof[0], nbk->getDOF(1, 0))) {
+ ov = 1;
+ if (nbk->getChild(0)) {
+ if (fill_opp_coords.isAnySet()) {
+ if (nbk->getNewCoord()) {
+ oppCoord_[i] = *(nbk->getNewCoord());
+ } else {
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[i][j] = (elinfo_old->oppCoord_[cv[i]][j] + elinfo_old->coord_[ichild][j]) / 2;
+ }
+ }
+ }
+ (*neigh_local)[i] = nbk->getChild(0);
+ oppVertex_[i] = 3;
+ break;
+ }
+ } else {
+ TEST_EXIT(dof == nbk->getDOF(2) ||
+ mesh->associated(dof[0], nbk->getDOF(2, 0)))
+ ("opp_vertex not found\n");
+ ov = 2;
+ }
- if (fillFlag__local.isSet(Mesh::FILL_BOUND))
- {
- for (i = 0; i < 3; i++) {
- boundary_[10 + i] = elinfo_old->getBoundary(10 + cv[i]);
+ if (fill_opp_coords.isAnySet()) {
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[i][j] = elinfo_old->oppCoord_[cv[i]][j];
+ }
+ }
+ (*neigh_local)[i] = nbk;
+ oppVertex_[i] = ov;
+ break;
+ }
+
+ } /* end for k */
+ TEST_EXIT(k < 2)("child not found with vertex\n");
+ }
+ } else {
+ (*neigh_local)[i] = NULL;
}
+ } /* end for i */
+
+
+ /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
+
+ if (((*neigh_local)[3] = (*neigh_old)[ochild])) {
+ oppVertex_[3] = elinfo_old->oppVertex_[ochild];
+
+ if (fill_opp_coords.isAnySet()) {
+ for (int j = 0; j < dow; j++) {
+ oppCoord_[3][j] = elinfo_old->oppCoord_[ochild][j];
+ }
+ }
+ }
+ }
- boundary_[13] = elinfo_old->getBoundary(4);
-
- boundary_[0] = INTERIOR;
- boundary_[1] = elinfo_old->getBoundary(cv[1]);
- boundary_[2] = elinfo_old->getBoundary(cv[2]);
- boundary_[3] = elinfo_old->getBoundary(ochild);
+ if (fillFlag__local.isSet(Mesh::FILL_BOUND)) {
+ for (int i = 0; i < 3; i++) {
+ boundary_[10 + i] = elinfo_old->getBoundary(10 + cv[i]);
+ }
+
+ boundary_[13] = elinfo_old->getBoundary(4);
+
+ boundary_[0] = INTERIOR;
+ boundary_[1] = elinfo_old->getBoundary(cv[1]);
+ boundary_[2] = elinfo_old->getBoundary(cv[2]);
+ boundary_[3] = elinfo_old->getBoundary(ochild);
+
+ ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
+ for (int iedge = 0; iedge < 4; iedge++) {
+ boundary_[mesh_->getGeo(FACE) + iedge] =
+ elinfo_old->getBoundary(mesh_->getGeo(FACE)+ce[iedge]);
+ }
+ for (int iedge = 4; iedge < 6; iedge++) {
+ int i = 5 - cv[iedge - 3]; /* old vertex opposite new edge */
+ boundary_[mesh_->getGeo(FACE) + iedge] = elinfo_old->getBoundary(i);
+ }
- ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
- for (iedge=0; iedge<4; iedge++) {
- boundary_[mesh_->getGeo(FACE)+iedge] =
- elinfo_old->getBoundary(mesh_->getGeo(FACE)+ce[iedge]);
+ if (elinfo_old->getProjection(0) &&
+ elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) {
+
+ projection_[0] = elinfo_old->getProjection(0);
+ } else { // boundary projection
+ projection_[0] = NULL;
+ projection_[1] = elinfo_old->getProjection(cv[1]);
+ projection_[2] = elinfo_old->getProjection(cv[2]);
+ projection_[3] = elinfo_old->getProjection(ochild);
+
+ for (int iedge = 0; iedge < 4; iedge++) {
+ projection_[mesh_->getGeo(FACE)+iedge] =
+ elinfo_old->getProjection(mesh_->getGeo(FACE)+ce[iedge]);
}
- for (iedge=4; iedge<6; iedge++) {
- i = 5 - cv[iedge-3]; /* old vertex opposite new edge */
- boundary_[mesh_->getGeo(FACE)+iedge] = elinfo_old->getBoundary(i);
- }
-
- if(elinfo_old->getProjection(0) &&
- elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION)
- {
- projection_[0] = elinfo_old->getProjection(0);
- } else { // boundary projection
- projection_[0] = NULL;
- projection_[1] = elinfo_old->getProjection(cv[1]);
- projection_[2] = elinfo_old->getProjection(cv[2]);
- projection_[3] = elinfo_old->getProjection(ochild);
-
- for (iedge=0; iedge<4; iedge++) {
- projection_[mesh_->getGeo(FACE)+iedge] =
- elinfo_old->getProjection(mesh_->getGeo(FACE)+ce[iedge]);
- }
- for (iedge=4; iedge<6; iedge++) {
- i = 5 - cv[iedge-3]; /* old vertex opposite new edge */
- projection_[mesh_->getGeo(FACE)+iedge] = elinfo_old->getProjection(i);
- }
+ for (int iedge = 4; iedge < 6; iedge++) {
+ int i = 5 - cv[iedge - 3]; /* old vertex opposite new edge */
+ projection_[mesh_->getGeo(FACE) + iedge] = elinfo_old->getProjection(i);
}
}
+ }
-
+
if (fillFlag_.isSet(Mesh::FILL_ORIENTATION)) {
orientation =
( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elinfo_old)))->orientation
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index 348737acae0a5ab5669b3ece74c87eb24ab3fd35..6546270ca2084383bf4f8e8e031182028d785aaf 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/AMDiS/src/ElInfo3d.h
@@ -96,22 +96,30 @@ namespace AMDiS {
/** \brief
* get ElInfo's \ref el_type
*/
- inline unsigned char getType() const { return el_type; };
+ inline unsigned char getType() const {
+ return el_type;
+ };
/** \brief
* get ElInfo's \ref orientation
*/
- inline signed char getOrientation() const { return orientation; };
+ inline signed char getOrientation() const {
+ return orientation;
+ };
/** \brief
* set ElInfo's \ref el_type to t
*/
- inline void setType(unsigned char t) { el_type = t; };
+ inline void setType(unsigned char t) {
+ el_type = t;
+ };
/** \brief
* set ElInfo's \ref orientation to o
*/
- inline void setOrientation(signed char o) { orientation = o; };
+ inline void setOrientation(signed char o) {
+ orientation = o;
+ };
protected:
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index fa7985e6fa78aee84f0e9158d99ada41e93a1a9b..2df1b7ff8c605424859313d9c0f437dfe7292cbd 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -125,7 +125,7 @@ namespace AMDiS {
* Returns \ref child[i], i=0,1
*/
virtual Element* getChild(int i) const {
- FUNCNAME("Element::getChild");
+ FUNCNAME("Element::getChild()");
TEST_EXIT(i==0 || i==1)("i must be 0 or 1\n");
return child[i];
};
@@ -135,7 +135,7 @@ namespace AMDiS {
* false otherwise.
*/
inline const bool isLeaf() const {
- return (child[0]==NULL);
+ return (child[0] == NULL);
};
/** \brief
@@ -172,15 +172,15 @@ namespace AMDiS {
*/
inline double getEstimation(int row) const
{
- if(isLeaf()) {
+ if (isLeaf()) {
TEST_EXIT(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(ESTIMATABLE);
TEST_EXIT(ld)("leaf data not estimatable!\n");
- return dynamic_cast<LeafDataEstimatableInterface*>(ld)->
- getErrorEstimate(row);
- }
- else return 0.0;
+ return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
+ }
+
+ return 0.0;
};
/** \brief
@@ -188,15 +188,15 @@ namespace AMDiS {
* element and if it has leaf data and if this leaf data are coarsenable.
*/
inline double getCoarseningEstimation(int row) {
- if(isLeaf()) {
+ if (isLeaf()) {
TEST_EXIT(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(COARSENABLE);
TEST_EXIT(ld)("element data not coarsenable!\n");
- return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
- getCoarseningErrorEstimate(row);
+ return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->getCoarseningErrorEstimate(row);
}
- else return 0.0;
+
+ return 0.0;
};
/** \brief
@@ -204,7 +204,6 @@ namespace AMDiS {
*/
int getRegion() const;
-
/** \brief
* Returns local vertex number of the j-th vertex of the i-th edge
*/
@@ -218,8 +217,14 @@ namespace AMDiS {
int positionIndex,
int vertexIndex) const = 0;
+ /** \brief
+ *
+ */
virtual int getPositionOfVertex(int side, int vertex) const = 0;
+ /** \brief
+ *
+ */
virtual int getEdgeOfFace(int face, int edge) const = 0;
/** \brief
@@ -264,33 +269,46 @@ namespace AMDiS {
/** \brief
* Sets \ref child[0]
*/
- virtual void setFirstChild(Element *aChild) { child[0] = aChild; };
+ virtual void setFirstChild(Element *aChild) {
+ child[0] = aChild;
+ };
/** \brief
* Sets \ref child[1]
*/
- virtual void setSecondChild(Element *aChild) { child[1] = aChild; };
+ virtual void setSecondChild(Element *aChild) {
+ child[1] = aChild;
+ };
/** \brief
* Sets \ref elementData of Element
*/
- void setElementData(ElementData* ed) { elementData = ed; };
+ void setElementData(ElementData* ed) {
+ elementData = ed;
+ };
/** \brief
* Sets \ref newCoord of Element. Needed by refinement, if Element has a
* boundary edge on a curved boundary.
*/
- inline void setNewCoord(WorldVector<double>* coord) { newCoord = coord; };
+ inline void setNewCoord(WorldVector<double>* coord) {
+ newCoord = coord;
+ };
/** \brief
* Sets \ref mesh.
*/
- inline void setMesh(Mesh *m) { mesh = m; };
+ inline void setMesh(Mesh *m) {
+ mesh = m;
+ };
/** \brief
* Sets the pointer to the DOFs of the i-th node of Element
*/
- DegreeOfFreedom* setDOF(int i,DegreeOfFreedom* p) {dof[i]=p;return dof[i];};
+ DegreeOfFreedom* setDOF(int i, DegreeOfFreedom* p) {
+ dof[i] = p;
+ return dof[i];
+ };
/** \brief
* Checks whether Element is a leaf element and whether it has leaf data.
@@ -298,15 +316,14 @@ namespace AMDiS {
*/
inline void setEstimation(double est, int row)
{
- if(isLeaf()) {
+ if (isLeaf()) {
TEST_EXIT(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(ESTIMATABLE);
TEST_EXIT(ld)("leaf data not estimatable\n");
dynamic_cast<LeafDataEstimatableInterface*>(ld)->
setErrorEstimate(row, est);
- }
- else {
+ } else {
ERROR_EXIT("setEstimation only for leaf elements!\n");
}
};
@@ -317,15 +334,14 @@ namespace AMDiS {
*/
inline void setCoarseningEstimation(double est, int row)
{
- if(isLeaf()) {
+ if (isLeaf()) {
TEST_EXIT(elementData)("leaf element without leaf data\n");
ElementData *ld = elementData->getElementData(COARSENABLE);
TEST_EXIT(ld)("leaf data not coarsenable\n");
dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
setCoarseningErrorEstimate(row, est);
- }
- else {
+ } else {
ERROR_EXIT("setEstimation only for leaf elements!\n");
}
};
@@ -333,17 +349,24 @@ namespace AMDiS {
/** \brief
* Sets Elements \ref mark = mark + 1;
*/
- inline void incrementMark() {mark++;}
+ inline void incrementMark() {
+ mark++;
+ }
/** \brief
* Sets Elements \ref mark = mark - 1;
*/
- inline void decrementMark() {if (0<mark) mark--;};
+ inline void decrementMark() {
+ if (0 < mark)
+ mark--;
+ };
/** \brief
* Sets Element's \ref mark
*/
- inline void setMark(signed char m) {mark=m;};
+ inline void setMark(signed char m) {
+ mark = m;
+ };
/** \} */
@@ -358,7 +381,7 @@ namespace AMDiS {
* Used by Estimator for the jumps => same quadrature nodes from both sides!
*/
virtual const FixVec<int,WORLD>&
- sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
+ sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
/** \brief
* Returns a copy of itself. Needed by Mesh to create Elements by a
@@ -420,11 +443,10 @@ namespace AMDiS {
/** \brief
* Refines Element's leaf data
*/
- inline void refineElementData(Element* child1, Element* child2, int elType=0) {
- if(elementData) {
- bool remove =
- elementData->refineElementData(this, child1, child2, elType);
- if(remove) {
+ inline void refineElementData(Element* child1, Element* child2, int elType = 0) {
+ if (elementData) {
+ bool remove = elementData->refineElementData(this, child1, child2, elType);
+ if (remove) {
ElementData *tmp = elementData->getDecorated();
DELETE elementData;
elementData = tmp;
@@ -438,13 +460,13 @@ namespace AMDiS {
inline void coarsenElementData(Element* child1, Element* child2, int elType=0) {
ElementData *childData;
childData = child1->getElementData();
- if(childData) {
+ if (childData) {
childData->coarsenElementData(this, child1, child2, elType);
DELETE childData;
child1->setElementData(NULL);
}
childData = child2->getElementData();
- if(childData) {
+ if (childData) {
childData->coarsenElementData(this, child2, child1, elType);
DELETE childData;
child2->setElementData(NULL);
@@ -458,8 +480,11 @@ namespace AMDiS {
return elementData;
};
+ /** \brief
+ *
+ */
inline ElementData* getElementData(int typeID) const {
- if(elementData) {
+ if (elementData) {
return elementData->getElementData(typeID);
}
return NULL;
@@ -470,8 +495,8 @@ namespace AMDiS {
*/
bool deleteElementData(int typeID) {
FUNCNAME("Element::deleteElementData()");
- if(elementData) {
- if(elementData->isOfType(typeID)) {
+ if (elementData) {
+ if (elementData->isOfType(typeID)) {
ElementData *tmp = elementData;
elementData = elementData->getDecorated();
DELETE tmp;
@@ -490,12 +515,14 @@ namespace AMDiS {
* elementTyp is the type of this element (comes from ElInfo)
*/
bool isRefinedAtSide(int side, Element *el1, Element *el2,
- unsigned char elementTyp=255);
+ unsigned char elementTyp = 255);
/** \brief
* Returns whether Element's \ref newCoord is set
*/
- inline bool isNewCoordSet() const { return (newCoord != NULL);};
+ inline bool isNewCoordSet() const {
+ return (newCoord != NULL);
+ };
/** \brief
* Frees memory for \ref newCoord
@@ -518,7 +545,9 @@ namespace AMDiS {
/** \brief
* Sets Element's \ref index. Used by friend class Mesh.
*/
- inline void setIndex(int i) {index=i ; };
+ inline void setIndex(int i) {
+ index = i;
+ };
/** \brief
* Used by friend class Mesh while dofCompress
@@ -535,7 +564,7 @@ namespace AMDiS {
* Pointers to the two children of interior elements of the tree. Pointers
* to NULL for leaf elements.
*/
- Element *child[2];
+ Element *child[2];
/** \brief
* Vector of pointers to DOFs. These pointers must be available for elements
@@ -554,14 +583,14 @@ namespace AMDiS {
* and may be larger than the number of elements in the binary tree (the list
* of indices may have holes after coarsening).
*/
- int index;
+ int index;
/** \brief
* Marker for refinement and coarsening. if mark is positive for a leaf
* element, this element is refined mark times. if mark is negative for
* a leaf element, this element is coarsened -mark times.
*/
- signed char mark;
+ signed char mark;
/** \brief
* If the element has a boundary edge on a curved boundary, this is a pointer
@@ -576,21 +605,13 @@ namespace AMDiS {
/** \brief
* Pointer to the Mesh this element belongs to
*/
- Mesh* mesh;
+ Mesh* mesh;
/** \brief
* Pointer to Element's leaf data
*/
- ElementData* elementData;
-
- // struct ElementConnection
- // {
- // Element *connectedElement;
- // int connectedToPart; // number of edge/face of this element
- // ElementConnection* virtualConnection;
- // };
+ ElementData* elementData;
- // ::std::list<ElementConnection*> connections;
friend class Mesh;
diff --git a/AMDiS/src/ElementFileWriter.cc b/AMDiS/src/ElementFileWriter.cc
index 9d959bb87917cde500c23dc4f6e03362c40e0f17..b186b48cfb5af1c1f416da64c7d5b49d3c999832 100644
--- a/AMDiS/src/ElementFileWriter.cc
+++ b/AMDiS/src/ElementFileWriter.cc
@@ -373,36 +373,7 @@ namespace AMDiS {
while (elInfo) {
// Get element value.
val = vec[elInfo->getElement()->getIndex()];
-
- /* if ((elInfo->getCoord(0)[0] == 80.0) &&
- (elInfo->getCoord(0)[1] == 1.25) &&
- (elInfo->getCoord(1)[0] == 79.9609375) &&
- (elInfo->getCoord(1)[1] == 1.2109375) &&
- (elInfo->getCoord(2)[0] == 80.0) &&
- (elInfo->getCoord(2)[1] == 1.2109375))
- {
- ::std::cout << "------" << ::std::endl;
- ::std::cout << "val = " << val << " I = " << elInfo->getElement()->getIndex() << ::std::endl;
- ::std::cout << "PARENT-I = " << elInfo->getParent()->getIndex() << ::std::endl;
- ::std::cout << "vc = " << vc << ::std::endl;
- ::std::cout.precision(10);
- for (int i = 0; i < 3; i++) {
- ::std::cout << elInfo->getCoord(i)[0] << "/" << elInfo->getCoord(i)[1] << ::std::endl;
- }
- }
- */
-
- if (val > 0.3) {
- ::std::cout.precision(10);
- ::std::cout << "------" << ::std::endl;
- ::std::cout << "val = " << val << " I = " << elInfo->getElement()->getIndex() << ::std::endl;
- ::std::cout << "vc = " << vc << ::std::endl;
- for (int i = 0; i < 3; i++) {
- ::std::cout << elInfo->getCoord(i)[0] << "/" << elInfo->getCoord(i)[1] << ::std::endl;
- }
- }
-
-
+
// Write value for each vertex of each element.
for (int i = 0; i <= dim; i++) {
fout << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
diff --git a/AMDiS/src/Estimator.cc b/AMDiS/src/Estimator.cc
index 54bdb40bc2a1be90deb3fed956a0ccefe28bfb68..b1ce20824d5cca5d0d23ed66577413b059aedd7a 100755
--- a/AMDiS/src/Estimator.cc
+++ b/AMDiS/src/Estimator.cc
@@ -235,6 +235,10 @@ namespace AMDiS {
el = elInfo->getElement();
+ // if (el->getIndex() == 66065) {
+ // ::std::cout << "UNSER ELEMENT!" << ::std::endl;
+ // }
+
double det = elInfo->getDet();
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
@@ -505,6 +509,8 @@ namespace AMDiS {
double d = 0.5*(det + detNeigh);
+ // ::std::cout << "d: " << d << ::std::endl;
+
if (norm == NO_NORM || norm == L2_NORM)
val *= C1*h2_from_det(d, dim) * d;
else
@@ -534,6 +540,8 @@ namespace AMDiS {
est_el += val;
}
}
+
+ // ::std::cout << "Est: " << est_el << ::std::endl;
el->setEstimation(est_el, row);
diff --git a/AMDiS/src/Estimator.h b/AMDiS/src/Estimator.h
index 8f1cc13f9382e21391cbf5e3522ad10dbab434df..47e5ae286e2513bd1a5d357324699edf01f72dbe 100755
--- a/AMDiS/src/Estimator.h
+++ b/AMDiS/src/Estimator.h
@@ -295,12 +295,16 @@ namespace AMDiS {
/** \brief
* Sets \ref name
*/
- void setName(::std::string name_) { name = name_; };
+ void setName(::std::string name_) {
+ name = name_;
+ };
/** \brief
* Sets \ref row
*/
- void setRow(int r) { row = r; };
+ void setRow(int r) {
+ row = r;
+ };
protected:
/** \brief
@@ -361,11 +365,6 @@ namespace AMDiS {
virtual void exit();
- /* /\** \brief */
- /* * Implements \ref Estimator::estimate(). */
- /* *\/ */
- /* virtual double estimate(double timestep = 0.0); */
-
protected:
/** \brief
* Constant in front of element residual
@@ -387,19 +386,31 @@ namespace AMDiS {
*/
double C3;
+
int numSystems;
+
int numPoints;
+
int system;
+
int dim;
+
int degree;
+
Quadrature *quad;
+
FastQuadrature **quadFast;
+
const BasisFunction **basFcts;
double **uhEl;
+
double **uhOldEl;
+
double *uhQP;
+
double *uhOldQP;
+
double *riq;
};
diff --git a/AMDiS/src/FixVec.cc b/AMDiS/src/FixVec.cc
index b7740449f5aa9e9d97e5c9f2cde092700b49ab77..ed37d4187170e13d3b12de89accdc11d3769ae31 100644
--- a/AMDiS/src/FixVec.cc
+++ b/AMDiS/src/FixVec.cc
@@ -9,8 +9,8 @@ namespace AMDiS {
va_start(arg, size);
VectorOfFixVecs<DimVec<double> > *result =
NEW VectorOfFixVecs<DimVec<double> >(dim, size, NO_INIT);
- for(int i=0; i < size; i++) {
- for(int j=0; j < dim+1; j++) {
+ for (int i = 0; i < size; i++) {
+ for (int j = 0; j < dim + 1; j++) {
(*result)[i][j] = va_arg(arg, double);
}
}
@@ -24,7 +24,7 @@ namespace AMDiS {
va_start(arg, size);
double *result = GET_MEMORY(double, size);
- for(int i=0; i < size; i++) {
+ for (int i = 0; i < size; i++) {
result[i] = va_arg(arg, double);
}
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index 101236db39cda2ac387f0b7260d9d28b150ef7ea..2277103d14a4190a8dca20cd099b1f2d4a54ae10 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -103,7 +103,7 @@ namespace AMDiS {
FixVec(int dim, InitType initType, const T& ini)
: Vector<T>(calcSize(dim))
{
- TEST_EXIT(initType==DEFAULT_VALUE)("wrong initType or wrong initializer\n");
+ TEST_EXIT(initType == DEFAULT_VALUE)("wrong initType or wrong initializer\n");
this->set(ini);
}
@@ -112,8 +112,8 @@ namespace AMDiS {
*/
inline void init(int dim)
{
- TEST_EXIT(this->getSize() == 0)("already initialized\n");
- resize(calcSize(dim));
+ //TEST_EXIT(this->getSize() == 0)("already initialized\n");
+ this->resize(calcSize(dim));
};
/** \brief
@@ -174,12 +174,13 @@ namespace AMDiS {
* FixVec's constructors. size_ is the number of contained FixVecs. initType
* must be NO_INIT.
*/
- VectorOfFixVecs(int dim, int size, InitType initType) : size_(size)
+ VectorOfFixVecs(int dim, int size, InitType initType)
+ : size_(size)
{
- TEST_EXIT(initType==NO_INIT)("wrong initType or wrong initializer\n");
- //vec = GET_MEMORY(FixVecType*, size);
+ TEST_EXIT(initType == NO_INIT)("wrong initType or wrong initializer\n");
+
vec = GET_MEMORY(FixVecType*, size_);
- for(FixVecType** i=&vec[0]; i < &vec[size_]; i++) {
+ for (FixVecType** i=&vec[0]; i < &vec[size_]; i++) {
*i = NEW FixVecType(dim, NO_INIT);
}
}
@@ -192,10 +193,10 @@ namespace AMDiS {
VectorOfFixVecs(int dim, int s, InitType initType, const FixVecType* ini)
: size_(s)
{
- TEST_EXIT(initType==VALUE_LIST)("wrong initType or wrong initializer\n");
+ TEST_EXIT(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
+
vec = GET_MEMORY(FixVecType*, size_);
- //vec = new FixVecType*[size];
- for(FixVecType** i=&vec[0]; i < &vec[size_]; i++) {
+ for (FixVecType** i=&vec[0]; i < &vec[size_]; i++) {
*i = NEW FixVecType(ini[i]);
}
}
@@ -224,9 +225,8 @@ namespace AMDiS {
{
size_ = rhs.getSize();
vec = GET_MEMORY(FixVecType*, size_);
- int i;
- for(i=0; i < size_; i++) {
+ for (int i = 0; i < size_; i++) {
vec[i] = NEW FixVecType(*(rhs.vec[i]));
}
};
@@ -280,12 +280,16 @@ namespace AMDiS {
/** \brief
* Returns the \ref size of this VectorOfFixVecs
*/
- inline int getSize() const { return size_; };
+ inline int getSize() const {
+ return size_;
+ };
/** \brief
* Returns the size of the contained FixVecs
*/
- inline int getSizeOfFixVec() const { return vec[0]->getSize(); }
+ inline int getSizeOfFixVec() const {
+ return vec[0]->getSize();
+ }
protected:
/** \brief
@@ -520,6 +524,19 @@ namespace AMDiS {
return (*this);
};
+ /** \brief
+ * Sets the arrays value to the geometric midpoint of the points
+ * p1 and p2;
+ */
+ inline void setMidpoint(const WorldVector<T> &p1, const WorldVector<T> &p2)
+ {
+ int dow = Global::getGeo(WORLD);
+
+ for (int i = 0; i < dow; i++) {
+ this->valArray[i] = 0.5 * (p1[i] + p2[i]);
+ }
+ }
+
/** \brief
* Multplication of a matrix with a vector.
*/
diff --git a/AMDiS/src/FixVec.hh b/AMDiS/src/FixVec.hh
index a35d2d0ad82aa70d415f4f178685691476ff6ec7..b07eeb0a53e08a452f0e5a92962acfc462d7d0c4 100644
--- a/AMDiS/src/FixVec.hh
+++ b/AMDiS/src/FixVec.hh
@@ -3,14 +3,11 @@ namespace AMDiS {
template<typename T, GeoIndex d>
double absteukl(const FixVec<T,d>& a,const FixVec<T,d>& b)
{
- int i;
- double erg;
+ double erg = 0.0;
- erg=0.0;
- for (i=0; i < a.getSize() ; i++)
- {
- erg=erg+((a[i]-b[i])*(a[i]-b[i]));
- }
+ for (int i = 0; i < a.getSize() ; i++) {
+ erg = erg + ((a[i] - b[i]) * (a[i] - b[i]));
+ }
return sqrt(erg);
}
diff --git a/AMDiS/src/MacroElement.h b/AMDiS/src/MacroElement.h
index 5fd6e7ac6df32ed33bc15a38c2d7c72d37517b5e..03691f5f8691319acc1f04fb71c098b34742286d 100644
--- a/AMDiS/src/MacroElement.h
+++ b/AMDiS/src/MacroElement.h
@@ -74,72 +74,63 @@ namespace AMDiS {
/** \brief
* Returns \ref index.
*/
- inline int getIndex() const
- {
+ inline int getIndex() const {
return index;
};
/** \brief
* Returns ref projection[i].
*/
- inline Projection *getProjection(int i) const
- {
+ inline Projection *getProjection(int i) const {
return projection[i];
};
/** \brief
* Returns \ref el
*/
- inline Element* getElement() const
- {
+ inline Element* getElement() const {
return element;
};
/** \brief
* Returns the i-th neighbour of this MacroElement \ref neighbour[i]
*/
- inline MacroElement* getNeighbour(int i) const
- {
+ inline MacroElement* getNeighbour(int i) const {
return neighbour[i];
};
/** \brief
* Returns the i-th opp-vertex of this MacroElement \ref oppVertex[i]
*/
- inline char getOppVertex(int i) const
- {
+ inline char getOppVertex(int i) const {
return oppVertex[i];
};
/** \brief
* Returns \ref coord[i]
*/
- inline WorldVector<double>& getCoord(int i)
- {
+ inline WorldVector<double>& getCoord(int i) {
return coord[i];
};
/** \brief
* Returns \ref coord
*/
- inline FixVec<WorldVector<double>, VERTEX>& getCoord()
- {
+ inline FixVec<WorldVector<double>, VERTEX>& getCoord() {
return coord;
};
/** \brief
* Returns \ref boundary[i]
*/
- inline BoundaryType getBoundary(int i) const
- {
+ inline BoundaryType getBoundary(int i) const {
return boundary[i];
};
/** \brief
* Returns \ref elType
*/
- inline unsigned char getElType() const
- {
+ inline unsigned char getElType() const {
return elType;
};
@@ -154,63 +145,59 @@ namespace AMDiS {
/** \brief
* Sets \ref index
*/
- inline void setIndex(int n)
- {
+ inline void setIndex(int n) {
index = n ;
};
/** \brief
* Sets \ref element if not yet set.
*/
- inline void setElement(Element* element_)
- {
- if (!element) element=element_;
- else if ((element!=element_)) ERROR("Trying to change element in MacroElement\n");
+ inline void setElement(Element* element_) {
+ if (!element) {
+ element = element_;
+ } else {
+ if (element != element_)
+ ERROR("Trying to change element in MacroElement\n");
+ }
};
/** \brief
* Sets \ref elType
*/
- inline void setElType(unsigned char typ)
- {
- elType=typ;
+ inline void setElType(unsigned char typ) {
+ elType = typ;
};
/** \brief
* Sets \ref projection[i] = p.
*/
- inline void setProjection(int i, Projection *p)
- {
+ inline void setProjection(int i, Projection *p) {
projection[i] = p;
};
/** \brief
* Sets the i-th Neighbour to n
*/
- inline void setNeighbour(int i, MacroElement *n)
- {
+ inline void setNeighbour(int i, MacroElement *n) {
neighbour[i] = n;
};
/** \brief
* Sets the i-th opp vertex to c
*/
- inline void setOppVertex(int i, char c)
- {
- oppVertex[i]=c;
+ inline void setOppVertex(int i, char c) {
+ oppVertex[i] = c;
};
/** \brief
* Sets \ref boundary[i] to b
*/
- inline void setBoundary(int i, BoundaryType b)
- {
+ inline void setBoundary(int i, BoundaryType b) {
boundary[i] = b;
};
- inline void setCoord(int i, const WorldVector<double> &c)
- {
+ inline void setCoord(int i, const WorldVector<double> &c) {
coord[i] = c;
};
@@ -267,6 +254,9 @@ namespace AMDiS {
*/
unsigned char elType;
+ /** \brief
+ *
+ */
::std::vector<int> *deserializedNeighbourIndices_;
// friend classes
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index ba9b28f7136515747135d8c0cfe0bfd20c27bc78..65b2c182f089b5da302e96d25e2223b76386aadb 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -28,24 +28,15 @@ namespace AMDiS {
{
FUNCNAME("Mesh::readMacro()");
- ::std::deque<MacroElement*>::iterator mel;
- MacroElement *neigh;
- WorldVector<double> *coords;
- int **melVertex;
- DegreeOfFreedom **dof;
- int i, j, k;
- WorldVector<double> x_min, x_max;
- MacroInfo *macroInfo;
-
TEST_EXIT(filename)("no file specified; filename NULL pointer\n");
-
- macroInfo = NEW MacroInfo();
+
+ MacroInfo *macroInfo = NEW MacroInfo();
macroInfo->readAMDiSMacro(filename, mesh);
- mel = macroInfo->mel.begin();
- melVertex = macroInfo->mel_vertex;
- coords = macroInfo->coords;
- dof = macroInfo->dof;
+ ::std::deque<MacroElement*>::iterator mel = macroInfo->mel.begin();
+ int **melVertex = macroInfo->mel_vertex;
+ WorldVector<double> *coords = macroInfo->coords;
+ DegreeOfFreedom **dof = macroInfo->dof;
// === read periodic data =================================
if (periodicFile&&(strcmp(periodicFile, "") != 0)) {
@@ -69,7 +60,7 @@ namespace AMDiS {
PeriodicMap periodicMap;
- for (i = 0; i < n; i++) {
+ for (int i = 0; i < n; i++) {
::std::map<int, int> vertexMapEl1;
::std::map<int, int> vertexMapEl2;
@@ -78,15 +69,15 @@ namespace AMDiS {
TEST_EXIT(fscanf(file, "%d", &boundaryType) == 1)("boundaryType?\n");
TEST_EXIT(fscanf(file, "%d", &el1) == 1)("el1?\n");
- for(j = 0; j < dim; j++) {
+ for (int j = 0; j < dim; j++) {
TEST_EXIT(fscanf(file, "%d", &verticesEl1[j]) == 1)("vertEl1[%d]\n", j);
}
TEST_EXIT(fscanf(file, "%d", &el2) == 1)("el2?\n");
- for(j = 0; j < dim; j++) {
+ for (int j = 0; j < dim; j++) {
TEST_EXIT(fscanf(file, "%d", &verticesEl2[j]) == 1)("vertEl2[%d]\n", j);
}
- for(j = 0; j < dim; j++) {
- if(mode == 0) {
+ for (int j = 0; j < dim; j++) {
+ if (mode == 0) {
periodicMap.setEntry(melVertex[el1][verticesEl1[j]],
melVertex[el2][verticesEl2[j]]);
}
@@ -96,29 +87,29 @@ namespace AMDiS {
// calculate sides of periodic vertices
int sideEl1 = 0, sideEl2 = 0;
- if(dim == 1) {
+ if (dim == 1) {
sideEl1 = verticesEl1[0];
sideEl2 = verticesEl2[0];
} else {
- for(j = 0; j < dim + 1; j++) {
+ for (int j = 0; j < dim + 1; j++) {
sideEl1 += j;
sideEl2 += j;
}
- for(j = 0; j < dim; j++) {
+ for (int j = 0; j < dim; j++) {
sideEl1 -= verticesEl1[j];
sideEl2 -= verticesEl2[j];
}
}
// create periodic info
- DimVec<WorldVector<double> > periodicCoordsEl1(dim-1, NO_INIT);
- DimVec<WorldVector<double> > periodicCoordsEl2(dim-1, NO_INIT);
+ DimVec<WorldVector<double> > periodicCoordsEl1(dim - 1, NO_INIT);
+ DimVec<WorldVector<double> > periodicCoordsEl2(dim - 1, NO_INIT);
- Element *element1 = const_cast<Element*>((*(mel+el1))->getElement());
- Element *element2 = const_cast<Element*>((*(mel+el2))->getElement());
+ Element *element1 = const_cast<Element*>((*(mel + el1))->getElement());
+ Element *element2 = const_cast<Element*>((*(mel + el2))->getElement());
// for all vertices of this side
- for (j = 0; j < dim; j++) {
+ for (int j = 0; j < dim; j++) {
periodicCoordsEl1[element1->getPositionOfVertex(sideEl1, verticesEl1[j])] =
coords[melVertex[el2][vertexMapEl1[verticesEl1[j]]]];
periodicCoordsEl2[element2->getPositionOfVertex(sideEl2, verticesEl2[j])] =
@@ -164,6 +155,7 @@ namespace AMDiS {
if (mode != 0) {
VertexVector *associated = mesh->periodicAssociations[boundaryType];
+
if (!associated) {
associated = NEW VertexVector(mesh->getVertexAdmin(), "vertex vector");
mesh->periodicAssociations[boundaryType] = associated;
@@ -173,7 +165,7 @@ namespace AMDiS {
}
}
- for(j = 0; j < dim; j++) {
+ for (int j = 0; j < dim; j++) {
(*associated)[melVertex[el1][verticesEl1[j]]] =
melVertex[el2][vertexMapEl1[verticesEl1[j]]];
(*associated)[melVertex[el2][verticesEl2[j]]] =
@@ -186,7 +178,7 @@ namespace AMDiS {
FREE_MEMORY(verticesEl2, int, dim);
// change periodic vertex dofs
- for (i = 0; i < mesh->getNumberOfVertices(); i++) {
+ for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
if (periodicMap.getEntry(i) != -1) {
mesh->freeDOF(dof[i], VERTEX);
dof[i] = dof[periodicMap.getEntry(i)];
@@ -194,16 +186,18 @@ namespace AMDiS {
::std::map<BoundaryType, VertexVector*>::iterator assoc;
::std::map<BoundaryType, VertexVector*>::iterator assocEnd =
mesh->periodicAssociations.end();
- for(assoc = mesh->periodicAssociations.begin();
- assoc != assocEnd;
- ++assoc)
- {
- DegreeOfFreedom a = (*(assoc->second))[i];
- if (a != i) {
- (*(assoc->second))[i] = i;
- (*(assoc->second))[a] = periodicMap.getEntry(i);
- }
+
+ for (assoc = mesh->periodicAssociations.begin();
+ assoc != assocEnd;
+ ++assoc) {
+
+ DegreeOfFreedom a = (*(assoc->second))[i];
+ if (a != i) {
+ (*(assoc->second))[i] = i;
+ (*(assoc->second))[a] = periodicMap.getEntry(i);
}
+ }
+
}
}
@@ -212,66 +206,61 @@ namespace AMDiS {
mesh->periodicAssociations.end();
for (assoc = mesh->periodicAssociations.begin();
assoc != assocEnd;
- ++assoc)
- {
- for( i = 0; i < mesh->getNumberOfVertices(); i++) {
- if (i != (*(assoc->second))[i])
- MSG("association %d: vertex %d -> vertex %d\n",
- assoc->first,
- i, (*(assoc->second))[i]);
- }
+ ++assoc) {
+
+ for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
+ if (i != (*(assoc->second))[i])
+ MSG("association %d: vertex %d -> vertex %d\n",
+ assoc->first,
+ i, (*(assoc->second))[i]);
}
- for (i = 0; i < mesh->getNumberOfVertices(); i++) {
+ }
+
+ for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
if (periodicMap.getEntry(i) != -1) {
MSG("identification : vertex %d is now vertex %d\n", i, periodicMap.getEntry(i));
}
}
}
+
// =========================================================
- for (i = 0; i < mesh->getNumberOfMacros(); i++)
- {
- for (k = 0; k < mesh->getGeo(VERTEX); k++)
- {
- (*(mel+i))->setCoord(k, coords[melVertex[i][k]]);
- const_cast<Element*>((*(mel+i))->getElement())->
- setDOF(k,dof[melVertex[i][k]]);
- }
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
+ for (int k = 0; k < mesh->getGeo(VERTEX); k++) {
+ (*(mel + i))->setCoord(k, coords[melVertex[i][k]]);
+ const_cast<Element*>((*(mel+i))->getElement())->
+ setDOF(k,dof[melVertex[i][k]]);
}
+ }
- if (!macroInfo->neigh_set)
- {
- TEST_EXIT(!periodicFile)
- ("periodic boundary condition => element neighbours must be set\n");
+ if (!macroInfo->neigh_set) {
+ TEST_EXIT(!periodicFile)
+ ("periodic boundary condition => element neighbours must be set\n");
computeNeighbours(mesh);
- }
- else
- {
+ } else {
/****************************************************************************/
/* fill MEL oppVertex values when reading neighbour information form file */
/****************************************************************************/
- for (i = 0; i < mesh->getNumberOfMacros(); i++)
- {
- for (k = 0; k < mesh->getGeo(NEIGH); k++)
- {
- if ((neigh = const_cast<MacroElement*>(mel[i]->getNeighbour(k))))
- {
- for (j = 0; j < mesh->getGeo(NEIGH); j++)
- if (neigh->getNeighbour(j) == *(mel+i)) break;
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
+ for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
+ MacroElement *neigh = const_cast<MacroElement*>(mel[i]->getNeighbour(k));
+
+ if (neigh) {
+ int j;
+ for (j = 0; j < mesh->getGeo(NEIGH); j++)
+ if (neigh->getNeighbour(j) == *(mel + i))
+ break;
- TEST_EXIT(j < mesh->getGeo(NEIGH))
- ("el %d no neighbour of neighbour %d\n",
- mel[i]->getIndex(), neigh->getIndex());
- mel[i]->setOppVertex(k,j);
- }
- else
- {
- mel[i]->setOppVertex(k,-1);
- }
- }
+ TEST_EXIT(j < mesh->getGeo(NEIGH))("el %d no neighbour of neighbour %d\n",
+ mel[i]->getIndex(), neigh->getIndex());
+ mel[i]->setOppVertex(k, j);
+ } else {
+ mel[i]->setOppVertex(k, -1);
}
+ }
}
+ }
if (!macroInfo->bound_set) {
macroInfo->dirichletBoundary();
@@ -285,17 +274,17 @@ namespace AMDiS {
int numFaces = mesh->getGeo(FACE);
int dim = mesh->getDim();
mel = mesh->firstMacroElement();
- for(i = 0; i < mesh->getNumberOfLeaves(); i++) {
+ for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
MacroElement *macroEl = *(mel+i);
Projection *projector = macroEl->getProjection(0);
- if(projector && projector->getType() == VOLUME_PROJECTION) {
- for(j = 0; j <= dim; j++) {
+ if (projector && projector->getType() == VOLUME_PROJECTION) {
+ for (int j = 0; j <= dim; j++) {
projector->project(macroEl->getCoord(j));
}
} else {
- for(j = 0; j < mesh->getGeo(EDGE); j++) {
+ for (int j = 0; j < mesh->getGeo(EDGE); j++) {
projector = macroEl->getProjection(numFaces + j);
- if(projector) {
+ if (projector) {
int vertex0 = Global::getReferenceElement(dim)->getVertexOfEdge(j, 0);
int vertex1 = Global::getReferenceElement(dim)->getVertexOfEdge(j, 1);
projector->project(macroEl->getCoord(vertex0));
@@ -309,40 +298,42 @@ namespace AMDiS {
macroInfo->fillBoundaryInfo(mesh);
if (mesh->getNumberOfDOFs(CENTER)) {
- for (i = 0; i < mesh->getNumberOfMacros(); i++)
- {
- const_cast<Element*>(mel[i]->getElement())->
- setDOF(mesh->getNode(CENTER),mesh->getDOF(CENTER));
- }
+ for (int i = 0; i < mesh->getNumberOfMacros(); i++) {
+ const_cast<Element*>(mel[i]->getElement())->
+ setDOF(mesh->getNode(CENTER),mesh->getDOF(CENTER));
+ }
}
/****************************************************************************/
/* domain size */
/****************************************************************************/
- for (j = 0; j < Global::getGeo(WORLD); j++)
- {
- x_min[j] = 1.E30;
- x_max[j] = -1.E30;
- }
+ WorldVector<double> x_min, x_max;
- for (i = 0; i < mesh->getNumberOfVertices(); i++)
- for (j = 0; j < Global::getGeo(WORLD); j++)
- {
- x_min[j] = ::std::min(x_min[j], coords[i][j]);
- x_max[j] = ::std::max(x_max[j], coords[i][j]);
- }
+ for (int j = 0; j < Global::getGeo(WORLD); j++) {
+ x_min[j] = 1.E30;
+ x_max[j] = -1.E30;
+ }
+
+ for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
+ for (int j = 0; j < Global::getGeo(WORLD); j++) {
+ x_min[j] = ::std::min(x_min[j], coords[i][j]);
+ x_max[j] = ::std::max(x_max[j], coords[i][j]);
+ }
+ }
- for (j = 0; j < Global::getGeo(WORLD); j++)
+ for (int j = 0; j < Global::getGeo(WORLD); j++)
mesh->setDiameter(j, x_max[j] - x_min[j]);
if (check) {
checkMesh(mesh);
if (mesh->getDim() > 1) {
- char filenew[128];
- strncpy(filenew, filename, 128); filenew[127] = 0;
- strncat(filenew, ".new", 128); filenew[127] = 0;
+ char filenew[128];
+ strncpy(filenew, filename, 128);
+ filenew[127] = 0;
+ strncat(filenew, ".new", 128);
+ filenew[127] = 0;
macroTest(mesh, filenew);
}
}
@@ -356,7 +347,8 @@ namespace AMDiS {
void MacroInfo::fill(Mesh *pmesh, int ne, int nv)
{
- FUNCNAME("MacroInfo::fill");
+ FUNCNAME("MacroInfo::fill()");
+
int i;
int dim=pmesh->getDim();
diff --git a/AMDiS/src/Marker.h b/AMDiS/src/Marker.h
index f9e81158648341a40e5c19671e78700b2ae6abf0..0bb7cb883f3e1c1346852c2392731e7b9cf305e8 100644
--- a/AMDiS/src/Marker.h
+++ b/AMDiS/src/Marker.h
@@ -80,8 +80,20 @@ namespace AMDiS {
if (!maximumMarking || (newMark > oldMark)) {
el->setMark(newMark);
- if (oldMark > 0) elMarkRefine--; else if(oldMark < 0) elMarkCoarsen--;
- if (newMark > 0) elMarkRefine++; else if(newMark < 0) elMarkCoarsen++;
+
+ if (oldMark > 0) {
+ elMarkRefine--;
+ } else {
+ if (oldMark < 0)
+ elMarkCoarsen--;
+ }
+
+ if (newMark > 0) {
+ elMarkRefine++;
+ } else {
+ if (newMark < 0)
+ elMarkCoarsen++;
+ }
}
};
diff --git a/AMDiS/src/Mesh.h b/AMDiS/src/Mesh.h
index d0e7b9304dc379702165bddc74ae38993332990a..74fca0653fdbc318aad99709c232a26161032abe 100644
--- a/AMDiS/src/Mesh.h
+++ b/AMDiS/src/Mesh.h
@@ -273,13 +273,17 @@ namespace AMDiS {
* Returns the size of \ref admin which is the number of the DOFAdmins
* belonging to this mesh
*/
- const int getNumberOfDOFAdmin() const {return admin.size();};
+ const int getNumberOfDOFAdmin() const {
+ return admin.size();
+ };
/** \brief
* Returns the size of \ref macroElements which is the number of
* of macro elements of this mesh
*/
- const int getNumberOfMacros() const {return macroElements.size();};
+ const int getNumberOfMacros() const {
+ return macroElements.size();
+ };
/** \brief
* Returns a DOFAdmin which at least manages vertex DOFs
@@ -320,16 +324,6 @@ namespace AMDiS {
return macroElements.end();
};
- // /** \brief
- // * Returns the begin of \ref meshes
- // */
- // static ::std::list<Mesh*>::iterator begin() { return meshes.begin(); };
-
- // /** \brief
- // * Returns the end of \ref meshes
- // */
- // static ::std::list<Mesh*>::iterator end() { return meshes.end(); };
-
/** \} */
// ==========================================================================
@@ -340,42 +334,58 @@ namespace AMDiS {
/** \brief
* Sets \ref name of the mesh
*/
- inline void setName(const ::std::string& aName) { name=aName;};
+ inline void setName(const ::std::string& aName) {
+ name = aName;
+ };
/** \brief
* Sets \ref nVertices of the mesh
*/
- inline void setNumberOfVertices(int n) { nVertices=n; };
+ inline void setNumberOfVertices(int n) {
+ nVertices = n;
+ };
/** \brief
* Sets \ref nFaces of the mesh
*/
- inline void setNumberOfFaces(int n) { nFaces=n; };
+ inline void setNumberOfFaces(int n) {
+ nFaces = n;
+ };
/** \brief
* Increments \ref nVertices by inc
*/
- inline void incrementNumberOfVertices(int inc) { nVertices += inc; };
+ inline void incrementNumberOfVertices(int inc) {
+ nVertices += inc;
+ };
/** \brief
* Sets \ref nEdges of the mesh
*/
- inline void setNumberOfEdges(int n) { nEdges=n; };
+ inline void setNumberOfEdges(int n) {
+ nEdges = n;
+ };
/** \brief
* Increments \ref nEdges by inc
*/
- inline void incrementNumberOfEdges(int inc) { nEdges += inc; };
+ inline void incrementNumberOfEdges(int inc) {
+ nEdges += inc;
+ };
/** \brief
* Increments \ref nFaces by inc
*/
- inline void incrementNumberOfFaces(int inc) { nFaces += inc; };
+ inline void incrementNumberOfFaces(int inc) {
+ nFaces += inc;
+ };
/** \brief
* Sets \ref nLeaves of the mesh
*/
- inline void setNumberOfLeaves(int n) { nLeaves=n; };
+ inline void setNumberOfLeaves(int n) {
+ nLeaves = n;
+ };
/** \brief
* Increments \ref nLeaves by inc
@@ -387,12 +397,16 @@ namespace AMDiS {
/** \brief
* Sets \ref nElements of the mesh
*/
- inline void setNumberOfElements(int n) { nElements=n; };
+ inline void setNumberOfElements(int n) {
+ nElements = n;
+ };
/** \brief
* Increments \ref nElements by inc
*/
- inline void incrementNumberOfElements(int inc) { nElements += inc; };
+ inline void incrementNumberOfElements(int inc) {
+ nElements += inc;
+ };
/** \brief
* Sets *\ref diam to w
@@ -407,17 +421,23 @@ namespace AMDiS {
/** \brief
* Sets \ref preserveCoarseDOFs = true
*/
- inline void retainCoarseDOFs() {preserveCoarseDOFs=true;};
+ inline void retainCoarseDOFs() {
+ preserveCoarseDOFs = true;
+ };
/** \brief
* Sets \ref preserveCoarseDOFs = b
*/
- inline void setPreserveCoarseDOFs(bool b) {preserveCoarseDOFs=b;};
+ inline void setPreserveCoarseDOFs(bool b) {
+ preserveCoarseDOFs = b;
+ };
/** \brief
* Sets \ref preserveCoarseDOFs = false
*/
- inline void noCoarseDOFs() {preserveCoarseDOFs=false;};
+ inline void noCoarseDOFs() {
+ preserveCoarseDOFs = false;
+ };
/** \brief
* Sets \ref elementPrototype of the mesh
@@ -437,7 +457,9 @@ namespace AMDiS {
parametric = param;
};
- inline void setMaxEdgeNeigh(int m) { maxEdgeNeigh = m; };
+ inline void setMaxEdgeNeigh(int m) {
+ maxEdgeNeigh = m;
+ };
/** \} */
// ==========================================================================
@@ -490,7 +512,9 @@ namespace AMDiS {
/** \brief
* Clears \ref macroElements
*/
- inline void clearMacroElements() { macroElements.clear();};
+ inline void clearMacroElements() {
+ macroElements.clear();
+ };
/** \brief
* Adds a macro element to the mesh
@@ -559,7 +583,7 @@ namespace AMDiS {
* Returns FILL_ANY_?D
*/
inline static const Flag& getFillAnyFlag(int dim) {
- switch(dim) {
+ switch (dim) {
case 1:
return FILL_ANY_1D;
break;
@@ -628,22 +652,22 @@ namespace AMDiS {
/** \brief
*
*/
- static const Flag FILL_COORDS ;
+ static const Flag FILL_COORDS;
/** \brief
*
*/
- static const Flag FILL_BOUND ;
+ static const Flag FILL_BOUND;
/** \brief
*
*/
- static const Flag FILL_NEIGH ;
+ static const Flag FILL_NEIGH;
/** \brief
*
*/
- static const Flag FILL_OPP_COORDS ;
+ static const Flag FILL_OPP_COORDS;
/** \brief
*
@@ -653,24 +677,31 @@ namespace AMDiS {
/** \brief
*
*/
- static const Flag FILL_ADD_ALL ;
+ static const Flag FILL_ADD_ALL;
/** \brief
*
*/
- static const Flag FILL_ANY_1D ;
+ static const Flag FILL_ANY_1D;
/** \brief
*
*/
- static const Flag FILL_ANY_2D ;
+ static const Flag FILL_ANY_2D;
/** \brief
*
*/
- static const Flag FILL_ANY_3D ;
+ static const Flag FILL_ANY_3D;
+ /** \brief
+ *
+ */
static const Flag FILL_DET;
+
+ /** \brief
+ *
+ */
static const Flag FILL_GRD_LAMBDA;
//**************************************************************************
@@ -874,6 +905,9 @@ namespace AMDiS {
*/
static ::std::vector<DegreeOfFreedom> dof_used;
+ /** \brief
+ *
+ */
static ::std::map<DegreeOfFreedom, DegreeOfFreedom*> serializedDOFs;
/** \brief
@@ -912,10 +946,20 @@ namespace AMDiS {
MacroInfo *macroFileInfo_;
protected:
- // for findElement-Fcts
+ /** \brief
+ * for findElement-Fcts
+ */
DimVec<double> final_lambda;
+
+ /** \brief
+ *
+ */
const WorldVector<double> *g_xy0, *g_xy;
- double *g_sp;
+
+ /** \brief
+ *
+ */
+ double *g_sp;
friend class MacroInfo;
friend class MacroReader;
diff --git a/AMDiS/src/PeriodicBC.cc b/AMDiS/src/PeriodicBC.cc
index 96a14871895cdc59d3901a5c4491437b5282b133..373d66f7e94801f0584d4c8fb1c74f26120a56bd 100644
--- a/AMDiS/src/PeriodicBC.cc
+++ b/AMDiS/src/PeriodicBC.cc
@@ -16,8 +16,8 @@ namespace AMDiS {
{
::std::vector<PeriodicDOFMapping*>::iterator it;
::std::vector<PeriodicDOFMapping*>::iterator end = mappings_.end();
- for(it = mappings_.begin(); it != end; ++it) {
- if((*it)->basFcts_ == basFcts) {
+ for (it = mappings_.begin(); it != end; ++it) {
+ if ((*it)->basFcts_ == basFcts) {
return *it;
}
}
@@ -32,9 +32,9 @@ namespace AMDiS {
{
FUNCNAME("PeriodicDOFMapping::PeriodicDOFMapping()");
TEST_EXIT(basFcts_->getDim() > 1)("dim == 1\n");
- int i, num = basFcts_->getNumber();
+ int num = basFcts_->getNumber();
DimVec<double> *lambda;
- for(i = 0; i < num; i++) {
+ for (int i = 0; i < num; i++) {
lambda = basFcts_->getCoords(i);
indexOfCoords_[*lambda] = i;
}
@@ -43,8 +43,8 @@ namespace AMDiS {
PeriodicDOFMapping::~PeriodicDOFMapping()
{
::std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> >::iterator it;
- for(it = dofPermutation_.begin(); it != dofPermutation_.end(); ++it) {
- if(it->second) {
+ for (it = dofPermutation_.begin(); it != dofPermutation_.end(); ++it) {
+ if (it->second) {
FREE_MEMORY(it->second, DegreeOfFreedom, basFcts_->getNumber());
}
}
@@ -55,15 +55,12 @@ namespace AMDiS {
{
FUNCNAME("PeriodicDOFMapping::getDOFPermutation()");
- if(dofPermutation_[vertexPermutation] == NULL) {
- //MSG("new dof permutation needed\n");
-
- int i, j;
+ if (dofPermutation_[vertexPermutation] == NULL) {
int dim = basFcts_->getDim();
int num = basFcts_->getNumber();
int sum = 0;
- for(i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < dim + 1; i++) {
sum += i - vertexPermutation[i];
TEST_EXIT(vertexPermutation[i] < dim + 1)
("invalid vertexPermuation\n");
@@ -76,28 +73,15 @@ namespace AMDiS {
DegreeOfFreedom *mapping = GET_MEMORY(DegreeOfFreedom, num);
- for(i = 0; i < num; i++) {
+ for (int i = 0; i < num; i++) {
lambda = basFcts_->getCoords(i);
- for(j = 0; j < dim + 1; j++) {
+ for (int j = 0; j < dim + 1; j++) {
newLambda[vertexPermutation[j]] = (*lambda)[j];
}
mapping[i] = indexOfCoords_[newLambda];
}
dofPermutation_[vertexPermutation] = mapping;
-
- // MSG("vertex permutation: ");
- // for(i = 0; i < dim+1; i++) {
- // MSG("%d ", vertexPermutation[i]);
- // }
- // MSG("\n");
-
- // MSG("dof permutation: ");
- // for(i = 0; i < num; i++) {
- // MSG("%d ", dofPermutation_[vertexPermutation][i]);
- // }
- // MSG("\n");
- // WAIT;
}
return dofPermutation_[vertexPermutation];
@@ -108,7 +92,7 @@ namespace AMDiS {
: BoundaryCondition(type, rowFESpace_, NULL),
masterMatrix_(NULL)
{
- if(rowFESpace->getMesh()->getDim() > 1) {
+ if (rowFESpace->getMesh()->getDim() > 1) {
periodicDOFMapping_ =
PeriodicDOFMapping::providePeriodicDOFMapping(rowFESpace_->getBasisFcts());
} else {
@@ -123,9 +107,7 @@ namespace AMDiS {
{
FUNCNAME("PeriodicBC::initMatrix()");
- // MSG("begin initMatrix...\n");
-
- if(!masterMatrix_) {
+ if (!masterMatrix_) {
masterMatrix_ = matrix;
Mesh *mesh = matrix->getRowFESpace()->getMesh();
@@ -140,8 +122,6 @@ namespace AMDiS {
neighIndices_ = GET_MEMORY(DegreeOfFreedom, num);
}
-
- //MSG("end initMatrix...\n");
}
void PeriodicBC::fillBoundaryCondition(DOFMatrix *matrix,
@@ -150,11 +130,10 @@ namespace AMDiS {
const BoundaryType *localBound,
int nBasFcts)
{
- if(matrix == masterMatrix_) {
+ if (matrix == masterMatrix_) {
int dim = rowFESpace->getMesh()->getDim();
- if(dim > 1) {
- int i, j;
+ if (dim > 1) {
DOFAdmin *admin = rowFESpace->getAdmin();
FixVec<int, WORLD> elFace(dim, NO_INIT);
@@ -172,11 +151,11 @@ namespace AMDiS {
GeoIndex sideGeoIndex = INDEX_OF_DIM(dim-1, dim);
- for(side = 0; side < dim + 1; side++) {
+ for (side = 0; side < dim + 1; side++) {
- if(elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
+ if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
- for(vertex = 0; vertex < dim; vertex++) {
+ for (vertex = 0; vertex < dim; vertex++) {
index = element->getVertexOfPosition(sideGeoIndex,
side,
vertex);
@@ -188,16 +167,18 @@ namespace AMDiS {
basFcts->getLocalIndices(neigh, admin, neighIndices_);
int oppVertex = 0;
- for(i = 0; i < dim + 1; i++) {
+ for (int i = 0; i < dim + 1; i++) {
// get vertex permutation
- if(i == side) {
+ if (i == side) {
vertexPermutation[i] = 0;
} else {
DegreeOfFreedom periodicDOF =
periodicDOFs[element->getPositionOfVertex(side, i)];
- for(j = 0; j < dim + 1; j++) {
- if(neigh->getDOF(j, 0) == periodicDOF) break;
+ int j;
+ for (j = 0; j < dim + 1; j++) {
+ if (neigh->getDOF(j, 0) == periodicDOF)
+ break;
}
vertexPermutation[i] = j;
}
@@ -210,8 +191,8 @@ namespace AMDiS {
periodicDOFMapping_->getDOFPermutation(vertexPermutation);
// set associated dofs
- for(i = 0; i < num; i++) {
- if((*(basFcts->getCoords(i)))[side] == 0) {
+ for (int i = 0; i < num; i++) {
+ if ((*(basFcts->getCoords(i)))[side] == 0) {
(*associated_)[dofIndices[i]] = neighIndices_[dofPermutation[i]];
}
}
@@ -227,7 +208,7 @@ namespace AMDiS {
TEST_EXIT(matrix)("no matrix\n");
- if(matrix == masterMatrix_) {
+ if (matrix == masterMatrix_) {
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int num = basFcts->getNumber();
FREE_MEMORY(neighIndices_, DegreeOfFreedom, num);
@@ -317,10 +298,10 @@ namespace AMDiS {
int row, col, newRow, newCol;
double entry, *newEntryPtr;
- for(rowIt = matrix->begin(), row = 0; rowIt != rowEnd; ++rowIt, ++row) {
+ for( rowIt = matrix->begin(), row = 0; rowIt != rowEnd; ++rowIt, ++row) {
rowSize = static_cast<int>(rowIt->size());
newRow = (*associated_)[row];
- for(colIndex = 0; colIndex < rowSize; colIndex++) {
+ for (colIndex = 0; colIndex < rowSize; colIndex++) {
col = (*rowIt)[colIndex].col;
if(col == DOFMatrix::UNUSED_ENTRY) continue;
if(col == DOFMatrix::NO_MORE_ENTRIES) break;
@@ -328,11 +309,11 @@ namespace AMDiS {
newEntryPtr = matrix->hasSparseDOFEntry(newRow, newCol);
- if((row < newRow) ||
- ((row == newRow) && (col < newCol)) ||
- !newEntryPtr)
+ if ((row < newRow) ||
+ ((row == newRow) && (col < newCol)) ||
+ !newEntryPtr)
{
- if(!newEntryPtr) {
+ if (!newEntryPtr) {
entry = 0.5 * (*rowIt)[colIndex].entry;
(*rowIt)[colIndex].entry = entry;
matrix->addSparseDOFEntry(1.0, newRow, newCol, entry, false);
@@ -343,8 +324,6 @@ namespace AMDiS {
}
}
}
-
- //MSG("end exitMatrix...\n");
}
void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
@@ -361,12 +340,12 @@ namespace AMDiS {
Mesh *mesh = vector->getFESpace()->getMesh();
VertexVector *associated_ = mesh->getPeriodicAssociations()[boundaryType];
- for(vecIt.reset(); !vecIt.end(); ++vecIt) {
+ for (vecIt.reset(); !vecIt.end(); ++vecIt) {
dof = vecIt.getDOFIndex();
newDOF = (*associated_)[dof];
- if(dof < newDOF) {
+ if (dof < newDOF) {
// ---------- different assemblage style --------
// (*vector)[newDOF] += (*vector)[dof];
// (*vector)[dof] = 0.0;
@@ -382,8 +361,6 @@ namespace AMDiS {
(*vector)[newDOF] = entry;
}
}
-
- //MSG("end exitVector...\n");
}
}
diff --git a/AMDiS/src/PeriodicBC.h b/AMDiS/src/PeriodicBC.h
index 3783582bf304b8bc2b666697f771337783c18502..46a6f57ec597adc4f607abcbda9a1cc8c0ca8b82 100644
--- a/AMDiS/src/PeriodicBC.h
+++ b/AMDiS/src/PeriodicBC.h
@@ -41,10 +41,12 @@ namespace AMDiS {
{
public:
bool operator()(const DimVec<T> &v1, const DimVec<T> &v2) {
- int i, size = v1.getSize();
- for(i = 0; i < size; i++) {
- if(v1[i] < v2[i]) return true;
- if(v1[i] > v2[i]) return false;
+ int size = v1.getSize();
+ for (int i = 0; i < size; i++) {
+ if (v1[i] < v2[i])
+ return true;
+ if (v1[i] > v2[i])
+ return false;
}
return false;
};
@@ -106,24 +108,13 @@ namespace AMDiS {
void exitVector(DOFVectorBase<double>* vector);
- // DOFVectorDOF *getAssociatedDOFs() { return associated_; };
-
- // void updateAssociatedDOFs();
-
- // void clearPeriodicElementData();
-
- // void refinePeriodicData(LeafDataPeriodic *pd, Element *parent, int newDOF,
- // int elType);
-
- // inline LeafDataPeriodic *getPeriodicElementData(int index) {
- // return periodicElementData_[index];
- // };
-
protected:
VertexVector *associated_;
+
PeriodicDOFMapping *periodicDOFMapping_;
+
DegreeOfFreedom *neighIndices_;
- // ::std::map<int, LeafDataPeriodic*> periodicElementData_;
+
DOFMatrix *masterMatrix_;
};
diff --git a/AMDiS/src/PeriodicMap.h b/AMDiS/src/PeriodicMap.h
index dade895cc93fa0e8572ea78cb9e81e14ebc4c02e..f77e00474c44efebcd68a6d149a45a5f6d48d4f8 100644
--- a/AMDiS/src/PeriodicMap.h
+++ b/AMDiS/src/PeriodicMap.h
@@ -33,12 +33,12 @@ namespace AMDiS {
public:
void setEntry(DegreeOfFreedom key, DegreeOfFreedom entry) {
// no trivial entries!
- if(key == entry)
+ if (key == entry)
return;
// if a key equal to entry exists ...
- if(getEntry(entry) >= 0) {
- if(getEntry(entry) == key) {
+ if (getEntry(entry) >= 0) {
+ if (getEntry(entry) == key) {
return;
}
// ... let it be equal entries
@@ -48,14 +48,14 @@ namespace AMDiS {
// replace entries equal to key
::std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator it;
- for(it = periodicMap.begin(); it != periodicMap.end(); ++it) {
- if(it->second == key) {
+ for (it = periodicMap.begin(); it != periodicMap.end(); ++it) {
+ if (it->second == key) {
it->second = entry;
}
}
// if key exists already
- if(getEntry(key) >= 0) {
+ if (getEntry(key) >= 0) {
// insert new entry with old entry as key
setEntry(getEntry(key), entry);
}
@@ -67,7 +67,7 @@ namespace AMDiS {
DegreeOfFreedom getEntry(DegreeOfFreedom key) {
::std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator it;
it = periodicMap.find(key);
- if(it != periodicMap.end()) {
+ if (it != periodicMap.end()) {
return it->second;
} else {
return -1;
diff --git a/AMDiS/src/RefinementManager.h b/AMDiS/src/RefinementManager.h
index d4ccc4c619df6e56726c2212629ed70602f2e279..23a9e843db2aebc5a3fb68156fa7757719c81ba9 100644
--- a/AMDiS/src/RefinementManager.h
+++ b/AMDiS/src/RefinementManager.h
@@ -70,12 +70,6 @@ namespace AMDiS {
*/
virtual ~RefinementManager() {};
- // /** \brief
- // * Refines a single Element by increment its mark and then refine
- // * the mesh.
- // */
- // virtual void refineElement(Element *element);
-
/** \brief
* Generates new coordinates on curved boundaries. Can be overriden by
* sub classes if used.
@@ -102,25 +96,35 @@ namespace AMDiS {
/** \brief
* set \ref newCoords
*/
- inline void newCoord(bool nc) { newCoords = nc; };
+ inline void newCoord(bool nc) {
+ newCoords = nc;
+ };
- inline bool newCoord() { return newCoords; };
+ inline bool newCoord() {
+ return newCoords;
+ };
/** \brief
* Implements the refinement of el_info->el. Can be overriden by sub
* classes if used.
*/
virtual ElInfo* refineFunction(ElInfo* ) {
- FUNCNAME("RefinementManager::refineFunction\n");
+ FUNCNAME("RefinementManager::refineFunction()");
ERROR_EXIT("called for base class!\n");
return NULL;
};
- inline void setMesh(Mesh *m) { mesh = m; };
+ inline void setMesh(Mesh *m) {
+ mesh = m;
+ };
- inline void setStack(TraverseStack *s) { stack = s; };
+ inline void setStack(TraverseStack *s) {
+ stack = s;
+ };
- inline TraverseStack *getStack() { return stack; };
+ inline TraverseStack *getStack() {
+ return stack;
+ };
protected:
/** \brief
diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc
index e1c9aa81351255eff1d3a3c9d35506cd17c16e12..b42e95325e1df018c6d74f470c01eb8710c72aa0 100644
--- a/AMDiS/src/RefinementManager2d.cc
+++ b/AMDiS/src/RefinementManager2d.cc
@@ -31,39 +31,34 @@ namespace AMDiS {
refineList->setElement(0, el_info->getElement());
n_neigh = 1;
- if(el_info->getProjection(0) &&
- el_info->getProjection(0)->getType() == VOLUME_PROJECTION)
- {
- newCoords = true;
- }
+ if (el_info->getProjection(0) &&
+ el_info->getProjection(0)->getType() == VOLUME_PROJECTION) {
+ newCoords = true;
+ }
/****************************************************************************/
/* give the refinement edge the right orientation */
/****************************************************************************/
- if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0))
- {
- edge[0] = const_cast<int*>( el_info->getElement()->getDOF(0));
- edge[1] = const_cast<int*>( el_info->getElement()->getDOF(1));
- }
- else
- {
- edge[1] = const_cast<int*>( el_info->getElement()->getDOF(0));
- edge[0] = const_cast<int*>( el_info->getElement()->getDOF(1));
- }
+ if (el_info->getElement()->getDOF(0,0) < el_info->getElement()->getDOF(1, 0)) {
+ edge[0] = const_cast<int*>( el_info->getElement()->getDOF(0));
+ edge[1] = const_cast<int*>( el_info->getElement()->getDOF(1));
+ } else {
+ edge[1] = const_cast<int*>( el_info->getElement()->getDOF(0));
+ edge[0] = const_cast<int*>( el_info->getElement()->getDOF(1));
+ }
/****************************************************************************/
/* get the refinement patch */
/****************************************************************************/
- if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh))
- {
- /****************************************************************************/
- /* domain's boundary was reached while looping around the refinement edge */
- /****************************************************************************/
- getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
- bound = true;
- }
-
+ if (getRefinePatch(&el_info, edge, 0, refineList, &n_neigh)) {
+ /****************************************************************************/
+ /* domain's boundary was reached while looping around the refinement edge */
+ /****************************************************************************/
+ getRefinePatch(&el_info, edge, 1, refineList, &n_neigh);
+ bound = true;
+ }
+
// ==========================================================================
// === check for periodic boundary ==========================================
// ==========================================================================
@@ -92,17 +87,17 @@ namespace AMDiS {
newDOF =
refinePatch(edge, periodicList, n_neigh_periodic, bound);
- if(firstNewDOF == -1) {
+ if (firstNewDOF == -1) {
firstNewDOF = newDOF;
}
-
- if(lastNewDOF != -1) {
- for(it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
- if(it->second) {
- if(((*(it->second))[edge[0][0]] == last_edge[0][0] &&
- (*(it->second))[edge[1][0]] == last_edge[1][0]) ||
- ((*(it->second))[edge[0][0]] == last_edge[1][0] &&
- (*(it->second))[edge[1][0]] == last_edge[0][0]))
+
+ if (lastNewDOF != -1) {
+ for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
+ if (it->second) {
+ if (((*(it->second))[edge[0][0]] == last_edge[0][0] &&
+ (*(it->second))[edge[1][0]] == last_edge[1][0]) ||
+ ((*(it->second))[edge[0][0]] == last_edge[1][0] &&
+ (*(it->second))[edge[1][0]] == last_edge[0][0]))
{
(*(it->second))[lastNewDOF] = newDOF;
(*(it->second))[newDOF] = lastNewDOF;
@@ -119,13 +114,13 @@ namespace AMDiS {
edge[1] = next_edge[1];
}
- if(lastNewDOF != firstNewDOF) {
- for(it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
- if(it->second) {
- if(((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
- (*(it->second))[first_edge[1][0]] == last_edge[1][0]) ||
- ((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
- (*(it->second))[first_edge[1][0]] == last_edge[0][0]))
+ if (lastNewDOF != firstNewDOF) {
+ for (it = mesh->getPeriodicAssociations().begin(); it != end; ++it) {
+ if (it->second) {
+ if (((*(it->second))[first_edge[0][0]] == last_edge[0][0] &&
+ (*(it->second))[first_edge[1][0]] == last_edge[1][0]) ||
+ ((*(it->second))[first_edge[0][0]] == last_edge[1][0] &&
+ (*(it->second))[first_edge[1][0]] == last_edge[0][0]))
{
(*(it->second))[lastNewDOF] = firstNewDOF;
(*(it->second))[firstNewDOF] = lastNewDOF;
@@ -150,27 +145,25 @@ namespace AMDiS {
int RefinementManager2d::newCoordsFct(ElInfo *el_info)
{
Element *el = el_info->getElement();
- int j;
int dow = Global::getGeo(WORLD);
Projection *projector = el_info->getProjection(0);
- if(!projector || projector->getType() != VOLUME_PROJECTION) {
+ if (!projector || projector->getType() != VOLUME_PROJECTION) {
projector = el_info->getProjection(2);
}
- if (el->getFirstChild() && projector && (!el->isNewCoordSet()))
- {
- WorldVector<double> *new_coord = NEW WorldVector<double>;
-
- for (j = 0; j < dow; j++)
- (*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j])*0.5;
-
- projector->project(*new_coord);
-
- el->setNewCoord(new_coord);
- }
+ if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
+ WorldVector<double> *new_coord = NEW WorldVector<double>;
+
+ for (int j = 0; j < dow; j++)
+ (*new_coord)[j] = (el_info->getCoord(0)[j] + el_info->getCoord(1)[j])*0.5;
+
+ projector->project(*new_coord);
+ el->setNewCoord(new_coord);
+ }
+
return 0;
}
@@ -202,14 +195,13 @@ namespace AMDiS {
mesh->incrementNumberOfVertices(1);
mesh->incrementNumberOfEdges(1);
- if (mesh->getNumberOfDOFs(EDGE))
- {
- /****************************************************************************/
- /* there are two additional dofs in the refinement edge */
- /****************************************************************************/
- dof[1] = mesh->getDOF(EDGE);
- dof[2] = mesh->getDOF(EDGE);
- }
+ if (mesh->getNumberOfDOFs(EDGE)) {
+ /****************************************************************************/
+ /* there are two additional dofs in the refinement edge */
+ /****************************************************************************/
+ dof[1] = mesh->getDOF(EDGE);
+ dof[2] = mesh->getDOF(EDGE);
+ }
/****************************************************************************/
/* first refine the element */
@@ -345,11 +337,10 @@ namespace AMDiS {
void RefinementManager2d::bisectTriangle(Triangle *el,
DegreeOfFreedom* newDOFs[3])
{
- FUNCNAME("RefinementManager2d::bisectTriangle");
+ FUNCNAME("RefinementManager2d::bisectTriangle()");
TEST_EXIT(mesh)("no mesh!\n");
Triangle *child[2];
- int i_child;
child[0] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
child[1] = dynamic_cast<Triangle*>(mesh->createNewElement(el));
@@ -364,7 +355,6 @@ namespace AMDiS {
/****************************************************************************/
el->refineElementData(child[0], child[1]); // call of subclass-method
- // el->freeElementData();
el->setFirstChild(child[0]);
el->setSecondChild(child[1]);
@@ -381,11 +371,10 @@ namespace AMDiS {
/* the other vertices are handed on from the parent */
/****************************************************************************/
- for (i_child = 0; i_child < 2; i_child++)
- {
- child[i_child]->setDOF(i_child, const_cast<int*>( el->getDOF(2)));
- child[i_child]->setDOF(1-i_child, const_cast<int*>( el->getDOF(i_child)));
- }
+ for (int i_child = 0; i_child < 2; i_child++) {
+ child[i_child]->setDOF(i_child, const_cast<int*>( el->getDOF(2)));
+ child[i_child]->setDOF(1-i_child, const_cast<int*>( el->getDOF(i_child)));
+ }
/****************************************************************************/
/* there is one more leaf element, two hierachical elements and one more */
@@ -396,39 +385,37 @@ namespace AMDiS {
mesh->incrementNumberOfLeaves(1);
mesh->incrementNumberOfElements(2);
- if (mesh->getNumberOfDOFs(EDGE))
- {
- /****************************************************************************/
- /* there are dof's in the midpoint of the edges */
- /****************************************************************************/
- DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);
-
- child[0]->setDOF(4, newEdgeDOFs);
- child[1]->setDOF(3, newEdgeDOFs);
-
- /****************************************************************************/
- /* dofs handed on by the parent */
- /****************************************************************************/
- child[0]->setDOF(5, const_cast<int*>( el->getDOF(4)));
- child[1]->setDOF(5, const_cast<int*>( el->getDOF(3)));
-
- /****************************************************************************/
- /* dofs in the refinement edge */
- /****************************************************************************/
- child[0]->setDOF(3, newDOFs[1]);
- child[1]->setDOF(4, newDOFs[2]);
- }
-
- if (mesh->getNumberOfDOFs(CENTER))
- {
- int node = mesh->getNode(CENTER);
-
- /****************************************************************************/
- /* there are dofs at the barycenter of the triangles */
- /****************************************************************************/
- child[0]->setDOF(node, mesh->getDOF(CENTER));
- child[1]->setDOF(node, mesh->getDOF(CENTER));
- }
+ if (mesh->getNumberOfDOFs(EDGE)) {
+ /****************************************************************************/
+ /* there are dof's in the midpoint of the edges */
+ /****************************************************************************/
+ DegreeOfFreedom* newEdgeDOFs = mesh->getDOF(EDGE);
+
+ child[0]->setDOF(4, newEdgeDOFs);
+ child[1]->setDOF(3, newEdgeDOFs);
+
+ /****************************************************************************/
+ /* dofs handed on by the parent */
+ /****************************************************************************/
+ child[0]->setDOF(5, const_cast<int*>( el->getDOF(4)));
+ child[1]->setDOF(5, const_cast<int*>( el->getDOF(3)));
+
+ /****************************************************************************/
+ /* dofs in the refinement edge */
+ /****************************************************************************/
+ child[0]->setDOF(3, newDOFs[1]);
+ child[1]->setDOF(4, newDOFs[2]);
+ }
+
+ if (mesh->getNumberOfDOFs(CENTER)) {
+ int node = mesh->getNode(CENTER);
+
+ /****************************************************************************/
+ /* there are dofs at the barycenter of the triangles */
+ /****************************************************************************/
+ child[0]->setDOF(node, mesh->getDOF(CENTER));
+ child[1]->setDOF(node, mesh->getDOF(CENTER));
+ }
}
int RefinementManager2d::getRefinePatch(ElInfo **el_info,
@@ -436,31 +423,30 @@ namespace AMDiS {
int dir, RCNeighbourList* refineList,
int *n_neigh)
{
- FUNCNAME("RefinementManager2d::getRefinePatch");
+ FUNCNAME("RefinementManager2d::getRefinePatch()");
Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( (*el_info)->getElement()));
ElInfo *neigh_info;
- int opp_vertex = 0;
-
- if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2)
- {
- /****************************************************************************/
- /* neighbour is not compatible devisible; refine neighbour first, store the*/
- /* opp_vertex to traverse back to el */
- /****************************************************************************/
- opp_vertex = (*el_info)->getOppVertex(2);
+ int opp_vertex = 0;
- neigh_info = stack->traverseNeighbour2d(*el_info, 2);
- neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(),
- 1));
- neigh_info = refineFunction(neigh_info);
-
- /****************************************************************************/
- /* now go back to the original element and refine the patch */
- /****************************************************************************/
- *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
- TEST_EXIT(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
- }
+ if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
+ /****************************************************************************/
+ /* neighbour is not compatible devisible; refine neighbour first, store the*/
+ /* opp_vertex to traverse back to el */
+ /****************************************************************************/
+ opp_vertex = (*el_info)->getOppVertex(2);
+
+ neigh_info = stack->traverseNeighbour2d(*el_info, 2);
+ neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(),
+ 1));
+ neigh_info = refineFunction(neigh_info);
+
+ /****************************************************************************/
+ /* now go back to the original element and refine the patch */
+ /****************************************************************************/
+ *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
+ TEST_EXIT(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
+ }
if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
TEST_EXIT((*el_info)->getOppVertex(2) == 2)
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 83a39215173ddc0234f75e6670dde28449bc8334..2bf2eadbdb00d82ebb1efeae04e6850d28d3bb10 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -15,8 +15,7 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level,
Flag fill_flag)
{
- FUNCNAME("TraverseStack::traverseFirst");
- int i;
+ FUNCNAME("TraverseStack::traverseFirst()");
traverse_mesh = mesh;
traverse_level = level;
@@ -27,19 +26,19 @@ namespace AMDiS {
Flag FILL_ANY = Mesh::getFillAnyFlag(mesh->getDim());
- for (i=0; i<stack_size; i++)
+ for (int i = 0; i < stack_size; i++)
elinfo_stack[i]->setFillFlag(fill_flag & FILL_ANY);
elinfo_stack[0]->setMesh(mesh);
elinfo_stack[1]->setMesh(mesh);
if (fill_flag.isSet(Mesh::CALL_LEAF_EL_LEVEL)) {
- TEST_EXIT(level >= 0)("invalid level: %d\n",level);
+ TEST_EXIT(level >= 0)("invalid level: %d\n", level);
}
traverse_mel = NULL;
- stack_used = 0;
- el_count = 0;
+ stack_used = 0;
+ el_count = 0;
return(traverseNext(NULL));
}
@@ -47,8 +46,8 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseNext(ElInfo* elinfo_old)
{
- FUNCNAME("TraverseStack::traverseNext");
- ElInfo *elinfo=NULL;
+ FUNCNAME("TraverseStack::traverseNext()");
+ ElInfo *elinfo = NULL;
ElInfo::traverseId = id;
@@ -89,8 +88,8 @@ namespace AMDiS {
/* MSG("total element count:%d\n",el_count); */
}
- if(elinfo) {
- if(parametric) {
+ if (elinfo) {
+ if (parametric) {
elinfo = parametric->addParametricInfo(elinfo);
}
elinfo->fillDetGrdLambda();
@@ -137,13 +136,14 @@ namespace AMDiS {
elinfo_new->fillElInfo(1, elinfo);
sum+=recursive(elinfo_new);
DELETE elinfo_new;
- }
- else {
- if (el_fct!=NULL) {
- if(parametric) elinfo = parametric->addParametricInfo(elinfo);
+ } else {
+ if (el_fct != NULL) {
+ if (parametric)
+ elinfo = parametric->addParametricInfo(elinfo);
elinfo->fillDetGrdLambda();
sum+=el_fct(elinfo);
- if(parametric) elinfo = parametric->removeParametricInfo(elinfo);
+ if (parametric)
+ elinfo = parametric->removeParametricInfo(elinfo);
}
}
return (flag.isSet(Mesh::FILL_ADD_ALL))?sum:0;;
@@ -162,11 +162,13 @@ namespace AMDiS {
DELETE elinfo_new;
}
else {
- if ((elinfo->getLevel() == level)&&(el_fct!=NULL)) {
- if(parametric) elinfo = parametric->addParametricInfo(elinfo);
+ if ((elinfo->getLevel() == level) && (el_fct!=NULL)) {
+ if (parametric)
+ elinfo = parametric->addParametricInfo(elinfo);
elinfo->fillDetGrdLambda();
sum+=el_fct(elinfo);
- if(parametric) elinfo = parametric->removeParametricInfo(elinfo);
+ if (parametric)
+ elinfo = parametric->removeParametricInfo(elinfo);
}
}
return (flag.isSet(Mesh::FILL_ADD_ALL))?sum:0;
@@ -322,86 +324,86 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseLeafElement()
{
- FUNCNAME("TraverseStack::traverseLeafElement");
+ FUNCNAME("TraverseStack::traverseLeafElement()");
Element *el;
- int i;
- if (stack_used == 0) /* first call */
- {
- currentMacro = traverse_mesh->firstMacroElement();
- if(currentMacro == traverse_mesh->endOfMacroElements()) return NULL;
- traverse_mel = *currentMacro;
+ if (stack_used == 0) { /* first call */
+ currentMacro = traverse_mesh->firstMacroElement();
+ if (currentMacro == traverse_mesh->endOfMacroElements())
+ return NULL;
+ traverse_mel = *currentMacro;
+
+ stack_used = 1;
+ elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
+ info_stack[stack_used] = 0;
+ el = elinfo_stack[stack_used]->getElement();
+ if ((el == NULL) || (el->getFirstChild() == NULL)) {
+ return(elinfo_stack[stack_used]);
+ }
+ } else {
+ el = elinfo_stack[stack_used]->getElement();
+
+ /* go up in tree until we can go down again */
+ while((stack_used > 0) &&
+ ((info_stack[stack_used] >= 2) || (el->getFirstChild() == NULL)))
+ {
+ stack_used--;
+ el = elinfo_stack[stack_used]->getElement();
+ }
+
+ /* goto next macro element */
+ if (stack_used < 1) {
+ currentMacro++;
+ if (currentMacro == traverse_mesh->endOfMacroElements())
+ return NULL;
+ traverse_mel = *currentMacro;
+
stack_used = 1;
elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
info_stack[stack_used] = 0;
-
+
el = elinfo_stack[stack_used]->getElement();
if ((el == NULL) || (el->getFirstChild() == NULL)) {
return(elinfo_stack[stack_used]);
}
}
- else
- {
- el = elinfo_stack[stack_used]->getElement();
-
- /* go up in tree until we can go down again */
- while((stack_used > 0) &&
- ((info_stack[stack_used] >= 2) || (el->getFirstChild()==NULL)))
- {
- stack_used--;
- el = elinfo_stack[stack_used]->getElement();
- }
-
- /* goto next macro element */
- if (stack_used < 1) {
- currentMacro++;
- if(currentMacro == traverse_mesh->endOfMacroElements()) return NULL;
- traverse_mel = *currentMacro;
-
- stack_used = 1;
- elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
- info_stack[stack_used] = 0;
-
- el = elinfo_stack[stack_used]->getElement();
- if ((el == NULL) || (el->getFirstChild() == NULL)) {
- return(elinfo_stack[stack_used]);
- }
- }
- }
-
+ }
+
/* go down tree until leaf */
- while (el->getFirstChild())
- {
- if (stack_used >= stack_size-1) enlargeTraverseStack();
- i = info_stack[stack_used];
- el = const_cast<Element*>( ((i==0) ? el->getFirstChild() : el->getSecondChild()));
- info_stack[stack_used]++;
- elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
- stack_used++;
+ while (el->getFirstChild()) {
+ if (stack_used >= stack_size-1)
+ enlargeTraverseStack();
+
+ int i = info_stack[stack_used];
+ el = const_cast<Element*>( ((i == 0) ? el->getFirstChild() : el->getSecondChild()));
+ info_stack[stack_used]++;
+ elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
+ stack_used++;
+
+
+ TEST_EXIT(stack_used < stack_size)
+ ("stack_size=%d too small, level=(%d,%d)\n",
+ stack_size, elinfo_stack[stack_used]->getLevel());
+
+ info_stack[stack_used] = 0;
+ }
-
- TEST_EXIT(stack_used < stack_size)
- ("stack_size=%d too small, level=(%d,%d)\n",
- stack_size, elinfo_stack[stack_used]->getLevel());
-
- info_stack[stack_used] = 0;
- }
-
return elinfo_stack[stack_used];
}
ElInfo* TraverseStack::traverseLeafElementLevel()
{
- FUNCNAME("TraverseStack::traverseLeafElementLevel");
+ FUNCNAME("TraverseStack::traverseLeafElementLevel()");
ERROR_EXIT("not yet implemented\n");
+
return NULL;
}
ElInfo* TraverseStack::traverseElementLevel()
{
- FUNCNAME("TraverseStack::traverseElementLevel");
+ FUNCNAME("TraverseStack::traverseElementLevel()");
ElInfo *elInfo;
do {
@@ -550,23 +552,31 @@ namespace AMDiS {
ElInfo* TraverseStack::traverseEveryElementPreorder()
{
- FUNCNAME("TraverseStack::traverseEveryElementPreorder");
+ FUNCNAME("TraverseStack::traverseEveryElementPreorder()");
Element *el;
- int i;
- if (stack_used == 0) /* first call */
- {
- currentMacro = traverse_mesh->firstMacroElement();
- traverse_mel = *currentMacro;
- if (traverse_mel == NULL) return NULL;
+ if (stack_used == 0) { /* first call */
+ currentMacro = traverse_mesh->firstMacroElement();
+ traverse_mel = *currentMacro;
+ if (traverse_mel == NULL)
+ return NULL;
+
+ stack_used = 1;
+ ::std::cout << "<1>" << ::std::endl;
+ elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
+ ::std::cout << "<2>" << ::std::endl;
+ info_stack[stack_used] = 0;
- stack_used = 1;
- elinfo_stack[stack_used]->fillMacroInfo(traverse_mel);
- info_stack[stack_used] = 0;
-
- return(elinfo_stack[stack_used]);
+ for (int i = 0; i < 3; i++) {
+ if (elinfo_stack[stack_used]->getNeighbour(i)) {
+ ::std::cout << "MACRO TEST SET[" << i << "] = " << (elinfo_stack[stack_used]->getNeighbour(i))->getIndex() << ::std::endl;
+ }
}
+ ::std::cout << "<3>" << ::std::endl;
+
+ return(elinfo_stack[stack_used]);
+ }
el = elinfo_stack[stack_used]->getElement();
@@ -582,7 +592,8 @@ namespace AMDiS {
/* goto next macro element */
if (stack_used < 1) {
currentMacro++;
- if (currentMacro == traverse_mesh->endOfMacroElements()) return(NULL);
+ if (currentMacro == traverse_mesh->endOfMacroElements())
+ return(NULL);
traverse_mel = *currentMacro;
stack_used = 1;
@@ -595,8 +606,10 @@ namespace AMDiS {
/* go down tree */
- if (stack_used >= stack_size-1) enlargeTraverseStack();
- i = info_stack[stack_used];
+ if (stack_used >= stack_size-1)
+ enlargeTraverseStack();
+
+ int i = info_stack[stack_used];
info_stack[stack_used]++;
elinfo_stack[stack_used + 1]->fillElInfo(i, elinfo_stack[stack_used]);
@@ -991,8 +1004,7 @@ namespace AMDiS {
if ((nb < 2) && (save_stack_used > 1)) {
stack2_used = 2; /* go down one level in OLD hierarchy */
- }
- else {
+ } else {
stack2_used = 1;
}
elinfo2 = save_elinfo_stack[stack2_used];
@@ -1000,15 +1012,15 @@ namespace AMDiS {
i = traverse_mel->getOppVertex(nb);
traverse_mel = traverse_mel->getNeighbour(nb);
- if (traverse_mel == NULL) return NULL;
+ if (traverse_mel == NULL)
+ return NULL;
nb = i;
stack_used = 1;
elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>( traverse_mel));
info_stack[stack_used] = 0;
- }
- else { /* goto other child */
+ } else { /* goto other child */
stack2_used = stack_used + 1;
if (save_stack_used > stack2_used) {
@@ -1087,21 +1099,22 @@ namespace AMDiS {
TEST_EXIT(elinfo->getNeighbour(opp_vertex) == old_elinfo->getElement())
("didn't succeed !?!?!?");
}
- if (elinfo->getElement()->getFirstChild())
- {
- MSG(" looking for neighbour %d of element %d at %8X\n",
- neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
- MSG(" originally: neighbour %d of element %d at %8X\n",
- sav_neighbour, sav_index, sav_el);
- MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
- elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
- elinfo->getNeighbour(opp_vertex)->getIndex());
- MSG("got no leaf element\n");
- WAIT_REALLY;
- }
+ if (elinfo->getElement()->getFirstChild()) {
+ MSG(" looking for neighbour %d of element %d at %8X\n",
+ neighbour, old_elinfo->getElement()->getIndex(), old_elinfo->getElement());
+ MSG(" originally: neighbour %d of element %d at %8X\n",
+ sav_neighbour, sav_index, sav_el);
+ MSG(" got element %d at %8X with opp_vertex %d neigh %d\n",
+ elinfo->getElement()->getIndex(), elinfo->getElement(), opp_vertex,
+ elinfo->getNeighbour(opp_vertex)->getIndex());
+ MSG("got no leaf element\n");
+ WAIT_REALLY;
+ }
+
+ if (elinfo) {
+ if (parametric)
+ elinfo = parametric->addParametricInfo(elinfo);
- if(elinfo) {
- if(parametric) elinfo = parametric->addParametricInfo(elinfo);
elinfo->fillDetGrdLambda();
}
@@ -1110,12 +1123,11 @@ namespace AMDiS {
void TraverseStack::update()
{
- FUNCNAME("TraverseStack::update");
- int i;
+ FUNCNAME("TraverseStack::update()");
TEST_EXIT(traverse_mesh->getDim() == 3)("update only in 3d\n");
- for (i=stack_used; i > 0; i--) {
+ for (int i = stack_used; i > 0; i--) {
dynamic_cast<ElInfo3d*>( elinfo_stack[i])->update();
}
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index ed60c08480d867e5dbcdacbe135214600f92dbc7..170caba3af8054b629d930b27d5990110e159afa 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -238,7 +238,7 @@ namespace AMDiS {
/** \brief
*
*/
- ::std::vector<ElInfo*> save_elinfo_stack;
+ ::std::vector<ElInfo*> save_elinfo_stack;
/** \brief
*
diff --git a/AMDiS/src/Triangle.h b/AMDiS/src/Triangle.h
index deee73684ba8d4825a9d9b6ca60c790395753ad8..e9a2a392a792d5933c1f46aac3a5a316068bbe52 100644
--- a/AMDiS/src/Triangle.h
+++ b/AMDiS/src/Triangle.h
@@ -48,7 +48,9 @@ namespace AMDiS {
/** \brief
* implements Element::clone
*/
- inline Element *clone() { return NEW Triangle(mesh); };
+ inline Element *clone() {
+ return NEW Triangle(mesh);
+ };
/** \brief
@@ -95,15 +97,6 @@ namespace AMDiS {
}
};
-
- /** \brief
- * implements Element::sideOfChild
- */
- //int sideOfChild(int child,
- // int side,
- // bool *isVirtual,
- // unsigned char elementType=0);
-
/** \brief
* implements Element::hasSide
*/
@@ -119,19 +112,25 @@ namespace AMDiS {
* implements Element::isLine. Returns false because this element is a
* Triangle
*/
- inline bool isLine() const { return false; };
+ inline bool isLine() const {
+ return false;
+ };
/** \brief
* implements Element::isTriangle. Returns true because this element is a
* Triangle
*/
- inline bool isTriangle() const { return true; };
+ inline bool isTriangle() const {
+ return true;
+ };
/** \brief
* implements Element::isTetrahedron. Returns false because this element is a
* Triangle
*/
- inline bool isTetrahedron() const { return false; };
+ inline bool isTetrahedron() const {
+ return false;
+ };
/** \brief
* implements Element::getSideOfChild()
@@ -166,7 +165,9 @@ namespace AMDiS {
// ===== Serializable implementation =====
- ::std::string getTypeName() const { return "Triangle"; };
+ ::std::string getTypeName() const {
+ return "Triangle";
+ };
protected:
/** \brief