diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 7686e186e1237207efc3aba1987f0a38570c6948..df62642d505711613c8fb1b586f03cdfd5a33367 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -498,7 +498,7 @@ namespace AMDiS {
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
- double psival;
+ double psival;
double *phival = GET_MEMORY(double, nCol);
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
@@ -544,7 +544,7 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
for (int j = 0; j < nCol; j++) {
- (*mat)[i][j] += quadrature->getWeight(iq)*c[iq]*psival*phival[j];
+ (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
}
}
}
@@ -1191,12 +1191,14 @@ namespace AMDiS {
Quad2::~Quad2()
{
- int nPoints = quadrature->getNumPoints();
- for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
- for (int j = 0; j < nPoints; j++) {
- DELETE tmpLALt[i][j];
+ if (!firstCall) {
+ int nPoints = quadrature->getNumPoints();
+ for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
+ for (int j = 0; j < nPoints; j++) {
+ DELETE tmpLALt[i][j];
+ }
+ DELETE [] tmpLALt[i];
}
- DELETE [] tmpLALt[i];
}
}
@@ -1219,7 +1221,7 @@ namespace AMDiS {
}
DimMat<double> **LALt = tmpLALt[myRank];
-
+
for (int i = 0; i < nPoints; i++) {
LALt[i]->set(0.0);
}
@@ -1227,8 +1229,8 @@ namespace AMDiS {
for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
}
-
- VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
+
+ VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
if (symmetric) {
for (int iq = 0; iq < nPoints; iq++) {
@@ -1239,10 +1241,11 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
(*mat)[i][i] += quadrature->getWeight(iq) *
- ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));
+ xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[i]);
for (int j = i + 1; j < nCol; j++) {
- double val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
+ double val = quadrature->getWeight(iq) *
+ xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
@@ -1257,8 +1260,8 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
for (int j = 0; j < nCol; j++) {
- (*mat)[i][j] += quadrature->getWeight(iq) *
- ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
+ (*mat)[i][j] += quadrature->getWeight(iq) *
+ xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
}
}
}
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 798d1069d302232ef6329fae44a136567bdd1941..02f1cb693ebac20856886f21b44034cc8fb95f5d 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -440,9 +440,6 @@ namespace AMDiS {
}
Operator *operat = op ? op : operators[0];
-
- // !!! Do no combine the next two lines into one !!!
- // ElementMatrix *elementMatrix = NULL;
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
if (op) {
@@ -452,17 +449,14 @@ namespace AMDiS {
::std::vector<double*>::iterator factorIt;
for (it = operators.begin(), factorIt = operatorFactor.begin();
it != operators.end();
- ++it, ++factorIt)
- {
- (*it)->getElementMatrix(elInfo,
- elementMatrix,
- *factorIt ? **factorIt : 1.0);
- }
-
+ ++it, ++factorIt) {
+ (*it)->getElementMatrix(elInfo,
+ elementMatrix,
+ *factorIt ? **factorIt : 1.0);
+ }
}
addElementMatrix(factor, *elementMatrix, bound);
- // DELETE elementMatrix;
}
Flag DOFMatrix::getAssembleFlag()
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index aacea37647cb9d61f56fd1f46ff03596edc8b0a7..3ec0cde47873b7515f71fb28f1eacd1053cb2f93 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -71,9 +71,9 @@ namespace AMDiS {
DOFAdmin *admin = feSpace->getAdmin();
// count number of nodes and dofs per node
- ::std::vector<int> numNodeDOFs;
- ::std::vector<int> numNodePreDOFs;
- ::std::vector<DimVec<double>*> bary;
+ std::vector<int> numNodeDOFs;
+ std::vector<int> numNodePreDOFs;
+ std::vector<DimVec<double>*> bary;
int numNodes = 0;
int numDOFs = 0;
@@ -102,7 +102,7 @@ namespace AMDiS {
}
// traverse mesh
- ::std::vector<bool> visited(getUsedSize(), false);
+ std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
@@ -213,10 +213,10 @@ namespace AMDiS {
}
template<>
- const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo,
- const Quadrature *quad,
+ const WorldVector<double> *DOFVectorBase<double>::getGrdAtQPs(const ElInfo *elInfo,
+ const Quadrature *quad,
const FastQuadrature *quadFast,
- WorldVector<double> *grdAtQPs) const
+ WorldVector<double> *grdAtQPs) const
{
FUNCNAME("DOFVector<double>::getGrdAtQPs()");
@@ -231,17 +231,12 @@ namespace AMDiS {
("invalid basis functions");
Element *el = elInfo->getElement();
-
int dim = elInfo->getMesh()->getDim();
int dow = Global::getGeo(WORLD);
-
const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
-
int numPoints = quadrature->getNumPoints();
-
static WorldVector<double> *grd = NULL;
-
WorldVector<double> *result;
if (grdAtQPs) {
@@ -589,9 +584,9 @@ namespace AMDiS {
}
// count number of nodes and dofs per node
- ::std::vector<int> numNodeDOFs;
- ::std::vector<int> numNodePreDOFs;
- ::std::vector<DimVec<double>*> bary;
+ std::vector<int> numNodeDOFs;
+ std::vector<int> numNodePreDOFs;
+ std::vector<DimVec<double>*> bary;
int numNodes = 0;
int numDOFs = 0;
@@ -619,7 +614,7 @@ namespace AMDiS {
}
// traverse mesh
- ::std::vector<bool> visited(getUsedSize(), false);
+ std::vector<bool> visited(getUsedSize(), false);
TraverseStack stack;
Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
@@ -661,8 +656,8 @@ namespace AMDiS {
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
- ::std::vector<DegreeOfFreedom>::iterator it;
- ::std::vector<DegreeOfFreedom>::iterator end = vec.end();
+ std::vector<DegreeOfFreedom>::iterator it;
+ std::vector<DegreeOfFreedom>::iterator end = vec.end();
DegreeOfFreedom pos = 0;
for (it = vec.begin(); it != end; ++it, ++pos) {
if (*it == dof) *it = pos;
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 5ae335c960dc61a83eaca0bce0c8827562cb8525..3a53bf06079145177cd868d5510a374c821ec62a 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -78,7 +78,7 @@ namespace AMDiS {
nBasFcts(0)
{};
- DOFVectorBase(const FiniteElemSpace *f, ::std::string n);
+ DOFVectorBase(const FiniteElemSpace *f, std::string n);
virtual ~DOFVectorBase();
@@ -121,28 +121,28 @@ namespace AMDiS {
operatorEstFactor.push_back(estFactor);
};
- inline ::std::vector<double*>::iterator getOperatorFactorBegin() {
+ inline std::vector<double*>::iterator getOperatorFactorBegin() {
return operatorFactor.begin();
};
- inline ::std::vector<double*>::iterator getOperatorFactorEnd() {
+ inline std::vector<double*>::iterator getOperatorFactorEnd() {
return operatorFactor.end();
};
- inline ::std::vector<double*>::iterator getOperatorEstFactorBegin() {
+ inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
return operatorEstFactor.begin();
};
- inline ::std::vector<double*>::iterator getOperatorEstFactorEnd() {
+ inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
return operatorEstFactor.end();
};
- inline ::std::vector<Operator*>::iterator getOperatorsBegin() {
+ inline std::vector<Operator*>::iterator getOperatorsBegin() {
return operators.begin();
};
- inline ::std::vector<Operator*>::iterator getOperatorsEnd() {
+ inline std::vector<Operator*>::iterator getOperatorsEnd() {
return operators.end();
};
@@ -156,22 +156,22 @@ namespace AMDiS {
*/
T evalUh(const DimVec<double>& lambda, DegreeOfFreedom* ind);
- inline ::std::vector<Operator*>& getOperators() {
+ inline std::vector<Operator*>& getOperators() {
return operators;
};
- inline ::std::vector<double*>& getOperatorFactor() {
+ inline std::vector<double*>& getOperatorFactor() {
return operatorFactor;
};
- inline ::std::vector<double*>& getOperatorEstFactor() {
+ inline std::vector<double*>& getOperatorEstFactor() {
return operatorEstFactor;
};
/** \brief
* Returns \ref name
*/
- inline const ::std::string& getName() const {
+ inline const std::string& getName() const {
return name;
};
@@ -192,7 +192,7 @@ namespace AMDiS {
/** \brief
*
*/
- ::std::string name;
+ std::string name;
/** \brief
*
@@ -202,17 +202,17 @@ namespace AMDiS {
/** \brief
*
*/
- ::std::vector<Operator*> operators;
+ std::vector<Operator*> operators;
/** \brief
*
*/
- ::std::vector<double*> operatorFactor;
+ std::vector<double*> operatorFactor;
/** \brief
*
*/
- ::std::vector<double*> operatorEstFactor;
+ std::vector<double*> operatorEstFactor;
/** \brief
*
@@ -228,7 +228,7 @@ namespace AMDiS {
* Are used to store temporary local dofs of an element.
* Thread safe.
*/
- ::std::vector<DegreeOfFreedom* > localIndices;
+ std::vector<DegreeOfFreedom* > localIndices;
};
@@ -314,12 +314,12 @@ namespace AMDiS {
/** \brief
* Constructs a DOFVector with name n belonging to FiniteElemSpace f
*/
- DOFVector(const FiniteElemSpace* f, ::std::string n);
+ DOFVector(const FiniteElemSpace* f, std::string n);
/** \brief
* Initialization.
*/
- void init(const FiniteElemSpace* f, ::std::string n);
+ void init(const FiniteElemSpace* f, std::string n);
/** \brief
* Copy Constructor
@@ -340,14 +340,14 @@ namespace AMDiS {
/** \brief
* Returns iterator to the begin of \ref vec
*/
- typename ::std::vector<T>::iterator begin() {
+ typename std::vector<T>::iterator begin() {
return vec.begin();
};
/** \brief
* Returns iterator to the end of \ref vec
*/
- typename ::std::vector<T>::iterator end() {
+ typename std::vector<T>::iterator end() {
return vec.end();
};
@@ -356,7 +356,7 @@ namespace AMDiS {
* DOFIndexedBase::compress()
*/
virtual void compressDOFIndexed(int first, int last,
- ::std::vector<DegreeOfFreedom> &newDof);
+ std::vector<DegreeOfFreedom> &newDof);
/** \brief
* Sets \ref refineInter to b
@@ -399,7 +399,7 @@ namespace AMDiS {
/** \brief
* Returns \ref vec
*/
- ::std::vector<T>& getVector() {
+ std::vector<T>& getVector() {
return vec;
};
@@ -580,13 +580,13 @@ namespace AMDiS {
// ===== Serializable implementation =====
- void serialize(::std::ostream &out) {
+ void serialize(std::ostream &out) {
unsigned int size = vec.size();
out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
};
- void deserialize(::std::istream &in) {
+ void deserialize(std::istream &in) {
unsigned int size;
in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
vec.resize(size);
@@ -626,7 +626,7 @@ namespace AMDiS {
/** \brief
* Name of this DOFVector
*/
- ::std::string name;
+ std::string name;
/** \brief
* FiniteElemSpace of the vector
@@ -636,7 +636,7 @@ namespace AMDiS {
/** \brief
* Data container
*/
- ::std::vector<T> vec;
+ std::vector<T> vec;
/** \brief
* Specifies whether interpolation should be performed after refinement
@@ -709,7 +709,7 @@ namespace AMDiS {
* Calls constructor of DOFVector<DegreeOfFreedom> and registers itself
* as DOFContainer at DOFAdmin
*/
- DOFVectorDOF(const FiniteElemSpace* feSpace_, ::std::string name_)
+ DOFVectorDOF(const FiniteElemSpace* feSpace_, std::string name_)
: DOFVector<DegreeOfFreedom>(feSpace_, name_)
{
feSpace->getAdmin()->addDOFContainer(this);
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index bb09516283aa3d5990e126917bd21aee995311b6..d9b3938f104e433ef1229ef2ec181f4f87549529 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -18,7 +18,7 @@
namespace AMDiS {
template<typename T>
- DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, ::std::string n)
+ DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
: feSpace(f),
name(n),
elementVector(NULL),
@@ -42,7 +42,7 @@ namespace AMDiS {
}
template<typename T>
- DOFVector<T>::DOFVector(const FiniteElemSpace* f, ::std::string n)
+ DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
: DOFVectorBase<T>(f, n),
refineInter(false),
coarsenOperation(NO_OPERATION)
@@ -51,7 +51,7 @@ namespace AMDiS {
}
template<typename T>
- void DOFVector<T>::init(const FiniteElemSpace* f, ::std::string n)
+ void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
{
name = n;
feSpace = f;
@@ -244,7 +244,7 @@ namespace AMDiS {
T m;
Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
- m = ::std::min(m, *vecIterator);
+ m = std::min(m, *vecIterator);
}
return m;
@@ -264,7 +264,7 @@ namespace AMDiS {
T m;
Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
- m = ::std::max(m, *vecIterator);
+ m = std::max(m, *vecIterator);
}
return m;
@@ -338,7 +338,7 @@ namespace AMDiS {
++rowIterator, ++vecIterator) {
sum = 0;
row = &(a[rowIterator.getDOFIndex()]);
- for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
+ for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
jcol = colIterator->col;
@@ -367,7 +367,7 @@ namespace AMDiS {
++rowIterator, ++xIterator) {
ax = alpha * (*xIterator);
row = &(a[rowIterator.getDOFIndex()]);
- for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
+ for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
jcol = colIterator->col;
@@ -726,7 +726,7 @@ namespace AMDiS {
template<typename T>
void DOFVector<T>::compressDOFIndexed(int first, int last,
- ::std::vector<DegreeOfFreedom> &newDOF)
+ std::vector<DegreeOfFreedom> &newDOF)
{
int i, j;
for(i = first; i <= last; i++) {
@@ -754,8 +754,8 @@ namespace AMDiS {
if (op) {
op->getElementVector(elInfo, this->elementVector);
} else {
- ::std::vector<Operator*>::iterator it;
- ::std::vector<double*>::iterator factorIt;
+ std::vector<Operator*>::iterator it;
+ std::vector<double*>::iterator factorIt;
for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
it != this->operators.end();
++it, ++factorIt) {
@@ -776,7 +776,7 @@ namespace AMDiS {
Flag DOFVectorBase<T>::getAssembleFlag()
{
Flag fillFlag(0);
- ::std::vector<Operator*>::iterator op;
+ std::vector<Operator*>::iterator op;
for(op = this->operators.begin(); op != this->operators.end(); ++op) {
fillFlag |= (*op)->getFillFlag();
}
@@ -962,7 +962,7 @@ namespace AMDiS {
double sum = 0;
if (!add)
*vecIterator = 0.0;
- for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
+ for(std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
@@ -1002,7 +1002,7 @@ namespace AMDiS {
!rowIterator.end();
++rowIterator, ++xIterator) {
T ax = (*xIterator);
- for(::std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
+ for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
colIterator != rowIterator->end();
colIterator++) {
int jcol = colIterator->col;
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 485cd385e3aafb8b90995262a486b287382fc1fd..f8c90468c51c211c9b9dcba6c2be124deaf2028d 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -14,6 +14,21 @@
namespace AMDiS {
+ ElInfo2d::ElInfo2d(Mesh *aMesh)
+ : ElInfo(aMesh)
+ {
+ e1 = NEW WorldVector<double>;
+ e2 = NEW WorldVector<double>;
+ normal = NEW WorldVector<double>;
+ }
+
+ ElInfo2d::~ElInfo2d()
+ {
+ DELETE e1;
+ DELETE e2;
+ DELETE normal;
+ }
+
void ElInfo2d::fillMacroInfo(const MacroElement * mel)
{
FUNCNAME("ElInfo::fillMacroInfo()");
@@ -394,21 +409,19 @@ namespace AMDiS {
double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
{
FUNCNAME("ElInfo2d::calcGrdLambda()");
-
- WorldVector<double> e1, e2;
- double adet = 0.0;
- const WorldVector<double> v0 = coord_[0];
testFlag(Mesh::FILL_COORDS);
+
+ double adet = 0.0;
int dim = mesh_->getDim();
for (int i = 0; i < dimOfWorld; i++) {
- e1[i] = coord_[1][i] - v0[i];
- e2[i] = coord_[2][i] - v0[i];
+ (*e1)[i] = coord_[1][i] - coord_[0][i];
+ (*e2)[i] = coord_[2][i] - coord_[0][i];
}
if (dimOfWorld == 2) {
- double sdet = e1[0] * e2[1] - e1[1] * e2[0];
+ double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
adet = abs(sdet);
if (adet < 1.0E-25) {
@@ -418,24 +431,18 @@ namespace AMDiS {
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[1][0] = (*e2)[1] * det1; // a11: (a_ij) = A^{-T}
+ grd_lam[1][1] = -(*e2)[0] * det1; // a21
+ grd_lam[2][0] = -(*e1)[1] * det1; // a12
+ grd_lam[2][1] = (*e1)[0] * det1; // 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;
-
- vectorProduct(e1, e2, normal);
+ } else {
+ vectorProduct(*e1, *e2, *normal);
- adet = norm(&normal);
+ adet = norm(normal);
if (adet < 1.0E-15) {
MSG("abs(det) = %lf\n", adet);
@@ -443,8 +450,8 @@ namespace AMDiS {
for (int j = 0; j < dimOfWorld; j++)
grd_lam[i][j] = 0.0;
} else {
- vectorProduct(e2, normal, grd_lam[1]);
- vectorProduct(normal, e1, grd_lam[2]);
+ vectorProduct(*e2, *normal, grd_lam[1]);
+ vectorProduct(*normal, *e1, grd_lam[2]);
double adet2 = 1.0 / (adet * adet);
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index 217196af3fdfda4aab7a066d9fc63fc2b7ac84c3..67b7291d4c8bdc443a1ac2e7c14b1304c5cda583 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/AMDiS/src/ElInfo2d.h
@@ -43,7 +43,9 @@ namespace AMDiS {
/** \brief
* Constructor. Calls ElInfo's protected Constructor.
*/
- ElInfo2d(Mesh* aMesh) : ElInfo(aMesh) {};
+ ElInfo2d(Mesh* aMesh);
+
+ ~ElInfo2d();
/** \brief
* 2-dimensional realisation of ElInfo's fillElInfo method.
@@ -74,6 +76,12 @@ namespace AMDiS {
* 2-dimensional realisation of ElInfo's getElementNormal method.
*/
double getElementNormal( WorldVector<double> &normal) const;
+
+ protected:
+ /** \brief
+ * Temp vectors for function \ref calcGrdLambda.
+ */
+ WorldVector<double> *e1, *e2, *normal;
};
}
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index be835f0cf6e3138a1d40aab3bcea42f1f32b4d5a..ca9a691ecbc4b02882125bf527279abd8f223311 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -449,15 +449,15 @@ namespace AMDiS {
/** \brief
* Calls the corresponding constructor of VectorOfFixVecs
*/
- DimMat(int dim, InitType initType=NO_INIT)
- : Matrix<T>(dim+1, dim+1)
+ DimMat(int dim, InitType initType = NO_INIT)
+ : Matrix<T>(dim + 1, dim + 1)
{};
/** \brief
* Calls the corresponding constructor of VectorOfFixVecs
*/
DimMat(int dim, InitType initType, const T& ini)
- : Matrix<T>(dim+1, dim+1)
+ : Matrix<T>(dim + 1, dim + 1)
{
TEST_EXIT_DBG(initType == DEFAULT_VALUE)
("wrong initType or wrong initializer\n");
@@ -468,7 +468,7 @@ namespace AMDiS {
* Calls the corresponding constructor of VectorOfFixVecs
*/
DimMat(int dim, InitType initType, T* ini)
- : Matrix<T>(dim+1, dim+1)
+ : Matrix<T>(dim + 1, dim + 1)
{
TEST_EXIT_DBG(initType == VALUE_LIST)("wrong initType or wrong initializer\n");
setValues(ini);
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 09e9c0fa59f4e995873d144a0501ef51ae9febde..bea2474a917cbe2c0dcd3fc9b7c7417b2ab35dde 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -757,10 +757,8 @@ namespace AMDiS {
}
// center
- if (dimOfPosition == 3) {
- if (degree == 4) {
- return &sortedCenterDeg4;
- }
+ if ((dimOfPosition == 3) && (degree == 4)) {
+ return &sortedCenterDeg4;
}
ERROR_EXIT("should not be reached\n");
diff --git a/AMDiS/src/MacroReader.cc b/AMDiS/src/MacroReader.cc
index c6733b0efc50e664721caa95607062b2947711b8..1ae037d23f572f03be4a844cda92b56c39e82fb6 100644
--- a/AMDiS/src/MacroReader.cc
+++ b/AMDiS/src/MacroReader.cc
@@ -349,35 +349,31 @@ namespace AMDiS {
TEST_EXIT(pmesh)("no mesh\n");
- int i;
- int dim = pmesh->getDim();
-
+ int dim = pmesh->getDim();
mesh = pmesh;
mesh->setNumberOfElements(ne);
mesh->setNumberOfLeaves(ne);
mesh->setNumberOfVertices(nv);
-
- for(i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
MacroElement *newMacro = NEW MacroElement(mesh->getDim());
mel.push_back(newMacro);
mesh->addMacroElement(mel[i]);
}
dof = GET_MEMORY(DegreeOfFreedom*, nv);
-
- coords = NEW WorldVector<double>[nv];
-
+ coords = NEW WorldVector<double>[nv];
mel_vertex = GET_MEMORY(int*, ne);
- for (i = 0; i < ne; i++) {
+
+ for (int i = 0; i < ne; i++) {
mel_vertex[i]=GET_MEMORY(int, mesh->getGeo(VERTEX));
}
- for (i = 0; i < nv; i++)
+ for (int i = 0; i < nv; i++)
dof[i] = mesh->getDOF(VERTEX);
- for (i = 0; i < ne; i++) {
+ for (int i = 0; i < ne; i++) {
mel[i]->element = mesh->createNewElement();
(mel)[i]->index = i;
@@ -387,8 +383,6 @@ namespace AMDiS {
}
neigh_set = false;
bound_set = false;
-
- return;
}
void MacroInfo::clear(int ne, int nv)
@@ -448,10 +442,9 @@ namespace AMDiS {
static int get_key_no(const char *key)
{
- int i;
-
- for (i = 0; i < N_KEYS; i++)
- if (!strcmp(keys[i], key)) return(i);
+ for (int i = 0; i < N_KEYS; i++)
+ if (!strcmp(keys[i], key))
+ return(i);
return(-1);
}
@@ -460,8 +453,8 @@ namespace AMDiS {
static const char *read_key(const char *line)
{
- static char key[100];
- char *k = key;
+ static char key[100];
+ char *k = key;
while (isspace(*line))
line++;
@@ -835,7 +828,6 @@ namespace AMDiS {
mel[i]->boundary[k] = DIRICHLET;
}
}
- return;
}
@@ -938,7 +930,7 @@ namespace AMDiS {
FUNCNAME("MacroReader::computeNeighbours()");
int dim = mesh->getDim();
- FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
+ FixVec<DegreeOfFreedom*, DIMEN> dof(dim, NO_INIT);
for (int i = 0; i < mesh->getNumberOfLeaves(); i++) {
for (int k = 0; k < mesh->getGeo(NEIGH); k++) {
@@ -963,7 +955,7 @@ namespace AMDiS {
for (int l = 0; l < dim; l++)
dof[l] = const_cast<DegreeOfFreedom*>(mesh->getMacroElement(i)->
getElement()->
- getDOF((k+l+1)%(dim+1)));
+ getDOF((k + l + 1) % (dim + 1)));
}
int j = 0;
@@ -978,14 +970,29 @@ namespace AMDiS {
}
}
- TEST_EXIT(j < mesh->getNumberOfLeaves())
- ("could not find neighbour %d of element %d\n", k, i);
+ if (j >= mesh->getNumberOfLeaves()) {
+ std::cout << "----------- ERROR ------------" << std::endl;
+ std::cout << "Cannot find neighbour " << k << " of element " << i << std::endl;
+ std::cout << " dim = " << dim << std::endl;
+ std::cout << " coords of element = ";
+ for (int l = 0; l <= dim; l++) {
+ std::cout << mesh->getMacroElement(i)->getCoord(l);
+ if (l < dim) {
+ std::cout << " / ";
+ }
+ }
+ std::cout << std::endl;
+ std::cout << " dofs = ";
+ for (int l = 0; l < dim; l++) {
+ std::cout << *(dof[l]) << " ";
+ }
+ std::cout << std::endl;
+
+ ERROR_EXIT("\n");
+ }
}
}
}
-
-
- return;
}
@@ -1153,8 +1160,6 @@ namespace AMDiS {
} else {
mesh->setMaxEdgeNeigh(dim-1);
}
-
- return;
}
/*
@@ -1181,8 +1186,6 @@ namespace AMDiS {
ERROR_EXIT("mesh->feSpace\n");
}
}
-
- return;
}
/****************************************************************************/
@@ -1474,8 +1477,8 @@ namespace AMDiS {
("neighbour index %d < element index %d\n", nei->getIndex(),
mel_index);
- if (nei->getIndex() < mel_index) return false;
-
+ if (nei->getIndex() < mel_index)
+ return false;
edge_no = Tetrahedron::edgeOfDOFs[j][k];
@@ -1547,8 +1550,7 @@ namespace AMDiS {
void MacroReader::fillMelBoundary(Mesh *mesh, MacroElement *mel,
FixVec<BoundaryType ,NEIGH> ind)
{
- int i;
- for(i=0; i < mesh->getGeo(NEIGH); i++) {
+ for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
mel->boundary[i] = ind[i];
}
}
@@ -1558,17 +1560,14 @@ namespace AMDiS {
::std::deque<MacroElement*>& macro_elements,
FixVec<int,NEIGH> ind)
{
- int k;
int dim = mel->element->getMesh()->getDim();
- for (k = 0; k < Global::getGeo(NEIGH, dim); k++)
- {
- if (ind[k] >= 0)
- mel->neighbour[k] = macro_elements[ind[k]];
- else
- mel->neighbour[k] = NULL;
- }
- return;
+ for (int k = 0; k < Global::getGeo(NEIGH, dim); k++) {
+ if (ind[k] >= 0)
+ mel->neighbour[k] = macro_elements[ind[k]];
+ else
+ mel->neighbour[k] = NULL;
+ }
}
@@ -1732,82 +1731,79 @@ namespace AMDiS {
int vertices = mesh->getGeo(VERTEX);
- if (!mel || test[mel->getIndex()]==1)
- {
- return;
- }
- else
- {
- test[mel->getIndex()]=1;
+ if (!mel || test[mel->getIndex()]==1) {
+ return;
+ } else {
+ test[mel->getIndex()]=1;
- laengstekante(mel->coord,l,v);
-
- if (v[1] == mesh->getGeo(VERTEX)) /*nur eine laengste Kante*/
- {
- umbvk(mesh,mel,v[0],ele);
-
- for (i=0; i < vertices; i++) {
- recumb(mesh,
- mel->neighbour[i],
- mel,
- test,
- l[0],
- 0,
- ele,
- umbvk);
- }
- return;
+ laengstekante(mel->coord,l,v);
+
+ if (v[1] == mesh->getGeo(VERTEX)) /*nur eine laengste Kante*/
+ {
+ umbvk(mesh,mel,v[0],ele);
+
+ for (i=0; i < vertices; i++) {
+ recumb(mesh,
+ mel->neighbour[i],
+ mel,
+ test,
+ l[0],
+ 0,
+ ele,
+ umbvk);
}
- else
- {
- if (ka == 1)
- {
- if (fabs((l[0]-lg)/lg) < 0.0000001)
- {
- for (i=0; i < vertices; i++)
- {
- if (mel->neighbour[i] == macroalt)
- {
- k=i;
- }
- }
-
- umbvk(mesh,mel,k,ele);
-
- for (i=0; i < vertices; i++) {
- recumb(mesh,
- mel->neighbour[i],
- mel,
- test,
- l[0],
- 0,
- ele,
- umbvk);
+ return;
+ }
+ else
+ {
+ if (ka == 1)
+ {
+ if (fabs((l[0]-lg)/lg) < 0.0000001)
+ {
+ for (i=0; i < vertices; i++)
+ {
+ if (mel->neighbour[i] == macroalt)
+ {
+ k=i;
+ }
}
- return;
+
+ umbvk(mesh,mel,k,ele);
+
+ for (i=0; i < vertices; i++) {
+ recumb(mesh,
+ mel->neighbour[i],
+ mel,
+ test,
+ l[0],
+ 0,
+ ele,
+ umbvk);
}
- }
-
- n=const_cast<MacroElement*>(mel->getNeighbour(v[0]));
- /*Nachbar an fiktiv verlaengerter Kante*/
- umbvk(mesh,mel,v[0],ele);
-
- recumb(mesh, n, mel,test,l[0],1,ele,umbvk);
- for (i=0; i < vertices; i++) {
- recumb(mesh,
- mel->neighbour[i],
- mel,
- test,
- l[0],
- 0,
- ele,
- umbvk);
+ return;
+ }
}
- return;
+
+ n=const_cast<MacroElement*>(mel->getNeighbour(v[0]));
+ /*Nachbar an fiktiv verlaengerter Kante*/
+ umbvk(mesh,mel,v[0],ele);
+
+ recumb(mesh, n, mel,test,l[0],1,ele,umbvk);
+ for (i=0; i < vertices; i++) {
+ recumb(mesh,
+ mel->neighbour[i],
+ mel,
+ test,
+ l[0],
+ 0,
+ ele,
+ umbvk);
}
- }
+ return;
+ }
+ }
}
-
+
// berechnet aus Koord. der Eckpunkte eines Elementes
// die Laenge der laengsten Kante
@@ -1939,8 +1935,6 @@ namespace AMDiS {
mesh->traverse(-1, Mesh::CALL_EVERY_EL_INORDER, basicNodeFct);
WAIT_REALLY;
}
-
- return;
}
int MacroReader::basicCheckFct(ElInfo* elinfo)
diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h
index e2f61122b7f0bd002e3b43f3a31031e98dca572f..db9bea122a9b7031d0ba6592d830ab44dbde6843 100644
--- a/AMDiS/src/MatrixVector.h
+++ b/AMDiS/src/MatrixVector.h
@@ -232,26 +232,26 @@ namespace AMDiS {
};
void print() const {
- ::std::cout << this->size << " vector: " << ::std::endl;
+ std::cout << this->size << " vector: " << std::endl;
for (int i = 0; i < size; i++) {
- ::std::cout << this->valArray[i] << " ";
+ std::cout << this->valArray[i] << " ";
}
- ::std::cout << ::std::endl;
+ std::cout << std::endl;
};
// ===== Serializable implementation =====
- void serialize(::std::ostream &out) {
+ void serialize(std::ostream &out) {
out.write(reinterpret_cast<const char*>(&size), sizeof(int));
out.write(reinterpret_cast<const char*>(valArray), size * sizeof(T));
};
- void deserialize(::std::istream &in) {
+ void deserialize(std::istream &in) {
in.read(reinterpret_cast<char*>(&size), sizeof(int));
in.read(reinterpret_cast<char*>(valArray), size * sizeof(T));
};
- ::std::string getTypeName() const {
+ std::string getTypeName() const {
return "Vector";
};
@@ -305,13 +305,11 @@ namespace AMDiS {
/** \brief
* Assignement operator.
*/
- inline const Matrix<T>& operator=(const T& scal)
- {
+ inline const Matrix<T>& operator=(const T& scal) {
return static_cast<const Matrix<T>&>(Vector<T>::operator=(scal));
};
- inline bool operator==(const Matrix<T>& rhs) const
- {
+ inline bool operator==(const Matrix<T>& rhs) const {
if (rows != rhs.getNumRows()) return false;
if (cols != rhs.getNumCols()) return false;
return Vector<T>::operator == (rhs);
@@ -320,48 +318,42 @@ namespace AMDiS {
/** \brief
* Comparison operator.
*/
- inline bool operator!=(const Matrix<T>& rhs) const
- {
- return !(*this==rhs);
+ inline bool operator!=(const Matrix<T>& rhs) const {
+ return !(*this == rhs);
}
/** \brief
* Acces to i-th matrix row.
*/
- inline T *operator[](int i)
- {
+ inline T *operator[](int i) {
return this->valArray + cols * i;
};
/** \brief
* Acces to i-th matrix row for constant matrices.
*/
- inline const T *operator[](int i) const
- {
+ inline const T *operator[](int i) const {
return this->valArray + cols * i;
};
/** \brief
* Returns \ref rows.
*/
- inline int getNumRows() const
- {
+ inline int getNumRows() const {
return rows;
};
/** \brief
* Return \ref cols.
*/
- inline int getNumCols() const
- {
+ inline int getNumCols() const {
return cols;
};
/** \brief
* Returns \ref rows.
*/
- inline int getSize() const
- {
+ inline int getSize() const {
return rows;
};
@@ -373,12 +365,12 @@ namespace AMDiS {
};
void print() const {
- ::std::cout << this->rows << " x " << this->cols << " matrix: " << ::std::endl;
+ std::cout << this->rows << " x " << this->cols << " matrix: " << std::endl;
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
- ::std::cout << this->valArray[i * cols + j] << " ";
+ std::cout << this->valArray[i * cols + j] << " ";
}
- ::std::cout << ::std::endl;
+ std::cout << std::endl;
}
};
@@ -407,22 +399,44 @@ namespace AMDiS {
for (resultIt = result.begin(), mIt = m.begin();
resultIt != result.end();
- ++resultIt)
- {
- *resultIt = 0;
- for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) {
- *resultIt += *mIt * *vIt;
- }
+ ++resultIt) {
+ *resultIt = 0;
+ for (vIt = v.begin(); vIt != v.end(); ++vIt, ++mIt) {
+ *resultIt += *mIt * *vIt;
}
+ }
return result;
};
/** \brief
- * Matrix vector multiplication.
+ * Computes x(Ay)^T, with x and y vectors, and A a matrix.
*/
template<typename T>
- inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m)
+ inline T xAy(const Vector<T>& x, const Matrix<T>& a, const Vector<T>& y)
{
+ TEST_EXIT_DBG(a.getNumCols() == x.getSize())("A and x not compatible\n");
+ TEST_EXIT_DBG(a.getNumCols() == y.getSize())("A and y not compatible\n");
+
+ T result = 0;
+ T tmp = 0;
+ T *aIt, *xIt, *yIt;
+
+ for (aIt = a.begin(), xIt = x.begin(); aIt != a.end(); ++xIt) {
+ tmp = 0;
+ for (yIt = y.begin(); yIt != y.end(); ++yIt, ++aIt) {
+ tmp += *aIt * *yIt;
+ }
+ result += *xIt * tmp;
+ }
+
+ return result;
+ }
+
+ /** \brief
+ * Matrix vector multiplication.
+ */
+ template<typename T>
+ inline const Vector<T>& operator*=(const Vector<T>& v, const Matrix<T>& m) {
return mv(m, v, v);
};
@@ -430,8 +444,7 @@ namespace AMDiS {
* Matrix vector multiplication.
*/
template<typename T>
- inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v)
- {
+ inline Vector<T> operator*(const Matrix<T>& m, const Vector<T>& v) {
Vector<T> result(v.getSize());
return mv(m, v, result);
};
@@ -441,17 +454,16 @@ namespace AMDiS {
* Scalar product.
*/
template<typename T>
- inline double operator*(const Vector<T>& v1, const Vector<T>& v2)
+ inline double operator*(const Vector<T>& v1, const Vector<T>& v2)
{
- double result = 0;
+ double result = 0.0;
T *v1It, *v2It;
for (v1It = v1.begin(), v2It = v2.begin();
v1It != v1.end();
- ++v1It, ++v2It)
- {
- result += *v1It * *v2It;
- }
+ ++v1It, ++v2It) {
+ result += *v1It * *v2It;
+ }
return result;
};
@@ -484,13 +496,13 @@ namespace AMDiS {
Vector<T>& result)
{
TEST_EXIT_DBG(v.getSize() == result.getSize())("invalid size\n");
+
T *vIt, *resultIt;
for (vIt = v.begin(), resultIt = result.begin();
vIt != v.end();
- ++vIt, ++resultIt)
- {
- *resultIt = scal * *vIt;
- }
+ ++vIt, ++resultIt) {
+ *resultIt = scal * *vIt;
+ }
return result;
};
@@ -525,44 +537,38 @@ namespace AMDiS {
T *xIt, *yIt;
for (xIt = x.begin(), yIt = y.begin();
xIt != x.end();
- ++xIt, ++yIt)
- {
- *yIt += a * *xIt;
- }
+ ++xIt, ++yIt) {
+ *yIt += a * *xIt;
+ }
return y;
};
template<typename T>
- inline const Vector<T>& operator*=(Vector<T>& v, const T& scal)
- {
+ inline const Vector<T>& operator*=(Vector<T>& v, const T& scal) {
return mult(scal, v, v);
};
template<typename T>
- inline Vector<T> operator*(const Vector<T>& v, const T& scal)
- {
+ inline Vector<T> operator*(const Vector<T>& v, const T& scal) {
Vector<T> result = v;
result *= scal;
return result;
};
template<typename T>
- inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal)
- {
+ inline const Vector<T>& operator+(const Vector<T>& v1, const T& scal) {
Vector<T> result(v1.getSize());
return add(v1, scal, result);
};
template<typename T>
- inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2)
- {
+ inline const Vector<T>& operator+=(Vector<T>& v1, const Vector<T>& v2) {
return add(v1, v2, v1);
};
template<typename T>
- inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2)
- {
+ inline Vector<T> operator+(const Vector<T>& v1, const Vector<T>& v2) {
Vector<T> result = v1;
result += v2;
return result;
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index c99b5a0863a1233b685a6d0bdb2589775f02f7fa..de40e63af3ef8f2caa1c0ff63daff96410f1a4cd 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -348,7 +348,8 @@ namespace AMDiS {
* Constructor.
*/
FactorLaplace_SOT(double f)
- : SecondOrderTerm(0) {
+ : SecondOrderTerm(0)
+ {
factor = new double;
*factor = f;
@@ -359,7 +360,8 @@ namespace AMDiS {
* Constructor.
*/
FactorLaplace_SOT(double *fptr)
- : SecondOrderTerm(0), factor(fptr)
+ : SecondOrderTerm(0),
+ factor(fptr)
{
setSymmetric(true);
};
@@ -369,10 +371,8 @@ namespace AMDiS {
*/
inline void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- //double det = elInfo->getDet();
- int iq;
- for(iq = 0; iq < numPoints; iq++)
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ for (int iq = 0; iq < numPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*factor));
};
@@ -381,7 +381,7 @@ namespace AMDiS {
* Implenetation of SecondOrderTerm::eval().
*/
void eval(int numPoints,
- const double *,
+ const double *,
const WorldVector<double> *,
const WorldMatrix<double> *D2UhAtQP,
double *result,
@@ -519,7 +519,8 @@ namespace AMDiS {
* Constructor
*/
Matrix_SOT(WorldMatrix<double> mat)
- : SecondOrderTerm(0), matrix(mat)
+ : SecondOrderTerm(0),
+ matrix(mat)
{
symmetric = matrix.isSymmetric();
setSymmetric(symmetric);
@@ -582,7 +583,9 @@ namespace AMDiS {
* Constructor.
*/
FactorIJ_SOT(int x_i, int x_j, double f)
- : SecondOrderTerm(0), xi(x_i), xj(x_j)
+ : SecondOrderTerm(0),
+ xi(x_i),
+ xj(x_j)
{
factor = new double;
*factor = f;
@@ -594,7 +597,10 @@ namespace AMDiS {
* Constructor.
*/
FactorIJ_SOT(int x_i, int x_j, double *fptr)
- : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr)
+ : SecondOrderTerm(0),
+ xi(x_i),
+ xj(x_j),
+ factor(fptr)
{
setSymmetric(xi == xj);
};