diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index b9486fb210dee9fc42cffb24be522c1d5d1f3328..27561fd4e80af1a68e198b67dee353aef5257ae1 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -1100,32 +1100,41 @@ namespace AMDiS {
Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
- {}
+ {
+ q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(),
+ owner->getColFESpace()->getBasisFcts(),
+ quadrature);
+ tmpLALt.resize(omp_get_max_threads());
+ for (int i = 0; i < omp_get_max_threads(); i++) {
+ tmpLALt[i] = NEW DimMat<double>*;
+ *(tmpLALt[i]) = NEW DimMat<double>(dim, NO_INIT);
+ }
+ }
+
+ Pre2::~Pre2()
+ {
+ for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
+ DELETE *(tmpLALt[i]);
+ DELETE tmpLALt[i];
+ }
+ }
void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
- DimMat<double> **LALt = NEW DimMat<double>*;
- *LALt=NEW DimMat<double>(dim, NO_INIT);
const int **nEntries;
const int *k, *l;
const double *values;
- double val;
- if (firstCall) {
- q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(),
- owner->getColFESpace()->getBasisFcts(),
- quadrature);
- firstCall = false;
- }
-
- LALt[0]->set(0.0);
int myRank = omp_get_thread_num();
+ DimMat<double> **LALt = tmpLALt[myRank];
+ DimMat<double> &tmpMat = *LALt[0];
+ tmpMat.set(0.0);
for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
}
- (*LALt[0]) *= elInfo->getDet();
+ tmpMat *= elInfo->getDet();
nEntries = q11->getNumberEntries();
@@ -1134,9 +1143,9 @@ namespace AMDiS {
k = q11->getKVec(i, i);
l = q11->getLVec(i, i);
values = q11->getValVec(i, i);
- val = 0.0;
+ double val = 0.0;
for (int m = 0; m < nEntries[i][i]; m++) {
- val += values[m] * (*LALt[0])[k[m]][l[m]];
+ val += values[m] * tmpMat[k[m]][l[m]];
}
(*mat)[i][i] += val;
@@ -1146,7 +1155,7 @@ namespace AMDiS {
values = q11->getValVec(i, j);
val = 0.0;
for (int m = 0; m < nEntries[i][j]; m++) {
- val += values[m] * (*LALt[0])[k[m]][l[m]];
+ val += values[m] * tmpMat[k[m]][l[m]];
}
(*mat)[i][j] += val;
(*mat)[j][i] += val;
@@ -1158,17 +1167,14 @@ namespace AMDiS {
k = q11->getKVec(i, j);
l = q11->getLVec(i, j);
values = q11->getValVec(i, j);
- val = 0.0;
+ double val = 0.0;
for (int m = 0; m < nEntries[i][j]; m++) {
- val += values[m] * (*LALt[0])[k[m]][l[m]];
+ val += values[m] * tmpMat[k[m]][l[m]];
}
(*mat)[i][j] += val;
}
}
}
-
- DELETE *LALt;
- DELETE LALt;
}
Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad)
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index fa01e0e5008a55ce2a5f540b77cea2edb0e90a85..df347ea93d915840d4e146941aac179c31081c76 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -859,6 +859,11 @@ namespace AMDiS {
*/
Pre2(Operator *op, Assembler *assembler, Quadrature *quad);
+ /** \brief
+ * Destructor.
+ */
+ ~Pre2();
+
/** \brief
* Implements SubAssembler::calculateElementMatrix().
*/
@@ -867,7 +872,7 @@ namespace AMDiS {
/** \brief
* Implements SubAssembler::calculateElementVector().
*/
- void calculateElementVector(const ElInfo *, ElementVector */*vec*/) {
+ void calculateElementVector(const ElInfo *, ElementVector *) {
ERROR_EXIT("should not be called\n");
};
@@ -878,6 +883,11 @@ namespace AMDiS {
*/
const Q11PsiPhi *q11;
+ /** \brief
+ * Thread safe temporary vector of DimMats for calculation in calculateElementMatrix().
+ */
+ ::std::vector< DimMat<double>** > tmpLALt;
+
friend class SecondOrderAssembler;
};
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 4198906a9b44fdd1a32e18ce65818aa35b5afdfd..91cca45d7ff797849873d1ce21314bfd20bcd62f 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -298,33 +298,33 @@ namespace AMDiS {
if (j!=static_cast<int>(matrix[a].size())) matrix[a][j].col=c;
}
- double *DOFMatrix::addSparseDOFEntry(double sign, int irow,
- int jcol, double entry,
- bool add)
+ void DOFMatrix::addSparseDOFEntry(double sign, int irow,
+ int jcol, double entry,
+ bool add)
{
FUNCNAME("DOFMatrix::addSparseDOFEntry()");
MatrixRow *row = &(matrix[irow]);
if (add && !entry)
- return NULL;
-
- double *result = NULL;
+ return;
- int i, freeCol = -1, rowSize = static_cast<int>( row->size());
+ int freeCol = -1;
+ int rowSize = static_cast<int>( row->size());
TEST_EXIT_DBG(jcol >= 0 &&
jcol < colFESpace->getAdmin()->getUsedSize())
("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1);
// first entry is diagonal entry
- if (rowFESpace==colFESpace)
+ if (rowFESpace == colFESpace)
if (rowSize == 0) {
MatEntry newEntry = {irow, 0.0};
row->push_back(newEntry);
rowSize = 1;
}
+ int i;
// search jcol
for (i = 0; i < rowSize; i++) {
// jcol found ?
@@ -350,22 +350,17 @@ namespace AMDiS {
if (!add)
(*row)[i].entry = 0.0;
(*row)[i].entry += sign * entry;
- result = &((*row)[i].entry);
} else {
- if(freeCol == -1) {
+ if (freeCol == -1) {
MatEntry newEntry = {jcol, sign * entry};
row->push_back(newEntry);
- result = &((*row)[row->size() - 1].entry);
} else {
(*row)[freeCol].col = jcol;
if (!add)
(*row)[freeCol].entry = 0.0;
(*row)[freeCol].entry += sign * entry;
- result = &((*row)[freeCol].entry);
}
- }
-
- return result;
+ }
}
void DOFMatrix::addMatEntry(int row, MatEntry entry)
@@ -527,7 +522,7 @@ namespace AMDiS {
entry += a[rowIndex][j].entry * b[logIndex][physIndex].entry;
}
}
- if(entry != 0.0) {
+ if (entry != 0.0) {
addSparseDOFEntry(1.0, rowIndex, i, entry);
}
}
@@ -552,7 +547,7 @@ namespace AMDiS {
double entry = aRowIt->entry * bRowIt->entry;
- if(entry != 0.0) {
+ if (entry != 0.0) {
addSparseDOFEntry(1.0, aCol, bCol, entry);
}
}
@@ -592,14 +587,14 @@ namespace AMDiS {
// add x contributions to this row
for(i=0; i < static_cast<int>((*xIterator).size()); i++) {
colIndex = (*xIterator)[i].col;
- if(colIndex >= 0) {
+ if (colIndex >= 0) {
addSparseDOFEntry(a, rowIndex, colIndex, (*xIterator)[i].entry);
}
}
// add y contributions to this row
for(i=0; i < static_cast<int>((*yIterator).size()); i++) {
colIndex = (*yIterator)[i].col;
- if(colIndex >= 0) {
+ if (colIndex >= 0) {
addSparseDOFEntry(1.0, rowIndex, colIndex, (*yIterator)[i].entry);
}
}
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index 27f0b99396aa7e5a5cf24e2fc341d39bbf26d8ab..9a105bb1e6a8fd3789049dc1ead01b4a75829a2f 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -515,9 +515,9 @@ namespace AMDiS {
* Creates an entry with logical indices irow, icol if there is no entry
* yet. Than sign * entry is added to the value at this logical indices
*/
- double *addSparseDOFEntry(double sign,
- int irow, int jcol, double entry,
- bool add = true);
+ void addSparseDOFEntry(double sign,
+ int irow, int jcol, double entry,
+ bool add = true);
inline double *hasSparseDOFEntry(int irow, int jcol) {
::std::vector<MatEntry>::iterator it;
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 9f3f51e76a71999fd51d3c54408745aba8bde871..6095bf018a62dba06496976e758a53b50eb3dd64 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -178,7 +178,7 @@ namespace AMDiS {
int dim = mesh_->getDim();
int posIndex = DIM_OF_INDEX(pos, dim);
int offset = indexOffset[dim - 1][posIndex];
-
+
return boundary_[offset + i];
}
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 4b6a21e08bd0ac243c25c31ced387a8e36bfda03..5ac0f2759b955eb4bf881aeb7efc9644319ccca5 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -190,7 +190,7 @@ namespace AMDiS {
/** \brief
* Get ElInfo's \ref boundary_[i]
*/
- virtual BoundaryType getBoundary(int i) const {
+ inline BoundaryType getBoundary(int i) const {
return boundary_[i];
};
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 8865b711fe214bd6bfdc8687bc0c2dc0062e45c5..88b49dd11d4c55c7d6c851a87abf0266132ebf5c 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -314,8 +314,6 @@ namespace AMDiS {
double det = 0.0;
- int myRank = omp_get_thread_num();
-
WorldVector<double> *e0 = &tmpWorldVecs[0];
WorldVector<double> *e1 = &tmpWorldVecs[1];
WorldVector<double> *e2 = &tmpWorldVecs[2];
@@ -362,13 +360,10 @@ namespace AMDiS {
const int (*cvg)[4] = NULL; /* cvg = child_vertex[el_type] */
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;
- FixVec<Element*, NEIGH> *neigh_local;
- Flag fill_opp_coords;
Mesh *mesh = elinfo_old->getMesh();
TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n");
@@ -378,7 +373,7 @@ namespace AMDiS {
fillFlag_ = fillFlag__local;
parent_ = el_old;
level_ = elinfo_old->level_ + 1;
- int el_type_local = ( dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>( elinfo_old)))->getType();
+ int el_type_local = (dynamic_cast<const ElInfo3d*>(elinfo_old))->getType();
el_type = (el_type_local + 1) % 3;
TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
@@ -409,8 +404,10 @@ namespace AMDiS {
if (fillFlag__local.isSet(Mesh::FILL_NEIGH) ||
fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
- neigh_local = &neighbour_;
- neigh_old = &elinfo_old->neighbour_;
+ FixVec<Element*, NEIGH> *neigh_local = &neighbour_;
+ const FixVec<Element*, NEIGH> *neigh_old = &elinfo_old->neighbour_;
+
+ Flag fill_opp_coords;
fill_opp_coords.setFlags(fillFlag__local & Mesh::FILL_OPP_COORDS);
/*----- nb[0] is other child --------------------------------------------*/
diff --git a/AMDiS/src/Element.h b/AMDiS/src/Element.h
index e0f371300d5013cdb8861a96ac4a6ebe7ddad27f..8c3c5026983706e21c7f38a77654e4a01702d4cb 100644
--- a/AMDiS/src/Element.h
+++ b/AMDiS/src/Element.h
@@ -110,22 +110,21 @@ namespace AMDiS {
/** \brief
* Returns \ref child[0]
*/
- virtual Element* getFirstChild() const {
+ inline Element* getFirstChild() const {
return child[0];
};
/** \brief
* Returns \ref child[1]
*/
- virtual Element* getSecondChild() const {
+ inline Element* getSecondChild() const {
return child[1];
};
/** \brief
* Returns \ref child[i], i=0,1
*/
- virtual Element* getChild(int i) const {
- FUNCNAME("Element::getChild()");
+ inline Element* getChild(int i) const {
TEST_EXIT_DBG(i==0 || i==1)("i must be 0 or 1\n");
return child[i];
};
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 473a0751df9d6d2bde7b4af3fa8966ab628b5a97..09e9c0fa59f4e995873d144a0501ef51ae9febde 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -790,26 +790,23 @@ namespace AMDiS {
el_info->testFlag(Mesh::FILL_BOUND);
- DimVec<int> parts(dim, NO_INIT);
- for (int i = 0; i < dim + 1; i++)
- parts[i] = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
-
- int index = 0;
-
// boundaries
+ int index = 0;
int offset = 0;
BoundaryType boundaryType;
for (int i = dim - 1; i > 0; i--)
- offset += parts[i];
+ offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
for (int i = 0; i < dim; i++) {
- for (int j = offset; j < offset + parts[i]; j++) {
+ int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
+ for (int j = offset; j < jto; j++) {
boundaryType = el_info->getBoundary(j);
- for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
+ int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
+ for (int k = 0; k < kto; k++) {
result[index++] = boundaryType;
}
}
- offset -= parts[i+1];
+ offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
}
// interior nodes in the center
diff --git a/AMDiS/src/OEMSolver.hh b/AMDiS/src/OEMSolver.hh
index b10d5d817f005dd16b89054282de78557bb46041..0bb86e0214291aceabf69c6411a9c83211aad563 100644
--- a/AMDiS/src/OEMSolver.hh
+++ b/AMDiS/src/OEMSolver.hh
@@ -96,7 +96,7 @@ namespace AMDiS {
}
if (res <= tolerance) {
- INFO(info,6)("finished successfully with %d iterations\n", iter);
+ INFO(info,6)("finished successfully with %d iterations\n");
return(1);
}
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index e662bfd8f8eecdbd8df47a4ab04abe04d4a5e400..ab73a1f8a3915ee10d48c933f47413d1d77b4a1d 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -122,22 +122,23 @@ namespace AMDiS {
}
void OperatorTerm::l1lt(const DimVec<WorldVector<double> >& Lambda,
- DimMat<double>& LALt,
- double factor)
+ DimMat<double>& LALt,
+ double factor)
{
- int i, j, k;
- static const int dimOfWorld = Global::getGeo(WORLD);
- int dim = LALt.getNumRows() - 1;
+ static const int dimOfWorld = Global::getGeo(WORLD);
+ int dim = LALt.getNumRows() - 1;
double val = 0.0;
- for (i = 0; i <= dim; i++) {
- for (val = k = 0; k < dimOfWorld; k++)
+ for (int i = 0; i <= dim; i++) {
+ val = 0.0;
+ for (int k = 0; k < dimOfWorld; k++)
val += Lambda[i][k] * Lambda[i][k];
val *= factor;
LALt[i][i] += val;
- for (j = i+1; j <= dim; j++) {
- for (val = k = 0; k < dimOfWorld; k++)
+ for (int j = i + 1; j <= dim; j++) {
+ val = 0.0;
+ for (int k = 0; k < dimOfWorld; k++)
val += Lambda[i][k] * Lambda[j][k];
val *= factor;
LALt[i][j] += val;
diff --git a/AMDiS/src/PeriodicBC.cc b/AMDiS/src/PeriodicBC.cc
index 35325323c30dce5c3961ad7fd2b8582f5628f10b..ed60cceab00b6fff8b74efc4bbab84e2ea5402b2 100644
--- a/AMDiS/src/PeriodicBC.cc
+++ b/AMDiS/src/PeriodicBC.cc
@@ -209,9 +209,8 @@ namespace AMDiS {
TEST_EXIT_DBG(matrix)("no matrix\n");
if (matrix == masterMatrix_) {
- const BasisFunction *basFcts = rowFESpace->getBasisFcts();
- int num = basFcts->getNumber();
- FREE_MEMORY(neighIndices_, DegreeOfFreedom, num);
+ FREE_MEMORY(neighIndices_, DegreeOfFreedom,
+ rowFESpace->getBasisFcts()->getNumber());
masterMatrix_ = NULL;
}
@@ -221,13 +220,15 @@ 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++) {
col = (*rowIt)[colIndex].col;
- if(col == DOFMatrix::UNUSED_ENTRY) continue;
- if(col == DOFMatrix::NO_MORE_ENTRIES) break;
+ if (col == DOFMatrix::UNUSED_ENTRY)
+ continue;
+ if (col == DOFMatrix::NO_MORE_ENTRIES)
+ break;
newCol = (*associated_)[col];
newEntryPtr = matrix->hasSparseDOFEntry(newRow, newCol);
diff --git a/AMDiS/src/Preconditioner.h b/AMDiS/src/Preconditioner.h
index 21dc117e7d3967a09637c5d844fc59e1618c1a60..6855691f60ae6f6ae7f86c6234d94f603a4f7e3e 100644
--- a/AMDiS/src/Preconditioner.h
+++ b/AMDiS/src/Preconditioner.h
@@ -241,7 +241,7 @@ namespace AMDiS {
int size = scalPrecons.getSize();
#ifdef _OPENMP
-#pragma omp parallel for
+#pragma omp parallel for num_threads(size)
#endif
for (i = 0; i < size; i++) {
scalPrecons[i]->init();
@@ -255,7 +255,7 @@ namespace AMDiS {
int i;
int size = scalPrecons.getSize();
#ifdef _OPENMP
-#pragma omp parallel for
+#pragma omp parallel for num_threads(size)
#endif
for (i = 0; i < size; i++) {
scalPrecons[i]->precon(x->getDOFVector(i));
@@ -269,7 +269,7 @@ namespace AMDiS {
int i;
int size = scalPrecons.getSize();
#ifdef _OPENMP
-#pragma omp parallel for
+#pragma omp parallel for num_threads(size)
#endif
for (i = 0; i < size; i++) {
scalPrecons[i]->exit();
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 7434c9ad3bcb45a8444eaf39c910cce453122009..094e7053c6c212298a329eabdbf57a3cffaca87d 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -701,6 +701,8 @@ namespace AMDiS {
#pragma omp parallel for
#endif
for (i = 0; i < numComponents_; i++) {
+ const BasisFunction *basisFcts = componentSpaces_[i]->getBasisFcts();
+
for (int j = 0; j < numComponents_; j++) {
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
@@ -734,7 +736,7 @@ namespace AMDiS {
} else {
BoundaryType *bound = NULL;
if (useGetBound_) {
- bound = GET_MEMORY(BoundaryType, componentSpaces_[i]->getBasisFcts()->getNumber());
+ bound = GET_MEMORY(BoundaryType, basisFcts->getNumber());
}
TraverseStack stack;
@@ -742,7 +744,7 @@ namespace AMDiS {
while (elInfo) {
if (useGetBound_) {
- componentSpaces_[i]->getBasisFcts()->getBound(elInfo, bound);
+ basisFcts->getBound(elInfo, bound);
}
if (assembleMatrix) {
@@ -767,7 +769,7 @@ namespace AMDiS {
matrix->getBoundaryManager()->exitMatrix(matrix);
if (useGetBound_) {
- FREE_MEMORY(bound, BoundaryType, componentSpaces_[i]->getBasisFcts()->getNumber());
+ FREE_MEMORY(bound, BoundaryType, basisFcts->getNumber());
}
}
diff --git a/AMDiS/src/RefinementManager.cc b/AMDiS/src/RefinementManager.cc
index 6f3b9ca2aaea352390f1b5a0f2a1e13fa66e265a..1e2abc776cd1ca86e0b3f9277bd65c0e9d0add37 100644
--- a/AMDiS/src/RefinementManager.cc
+++ b/AMDiS/src/RefinementManager.cc
@@ -15,9 +15,7 @@ namespace AMDiS {
RefinementManager* RefinementManager::traversePtr = NULL;
bool RefinementManager::doMoreRecursiveRefine = false;
int RefinementManager::callRefineInterpol = 0;
-
RCNeighbourList* RefinementManager3d::refList = NULL;
-
int RefinementManager::globalMark = 0;
int RefinementManager::globalRefineFct(ElInfo* elinfo)
@@ -32,7 +30,9 @@ namespace AMDiS {
{
FUNCNAME("RefinementManager::globalRefine()");
- if(mark <= 0) return static_cast<Flag>(0);
+ if (mark <= 0)
+ return static_cast<Flag>(0);
+
globalMark = mark;
aMesh->traverse(-1,
Mesh::CALL_LEAF_EL |
@@ -50,7 +50,6 @@ namespace AMDiS {
FUNCNAME("RefinementManager::refineMesh()");
mesh = aMesh;
-
int n_elements = mesh->getNumberOfLeaves();
ElInfo *el_info;
diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc
index a4223e9cc0e8040521e9159f48120d173a082298..7ac4d19c3564996299c67bf6e66bdd15a4fdfec8 100644
--- a/AMDiS/src/RefinementManager2d.cc
+++ b/AMDiS/src/RefinementManager2d.cc
@@ -17,8 +17,9 @@ namespace AMDiS {
ElInfo* RefinementManager2d::refineFunction(ElInfo* el_info)
{
- FUNCNAME("RefinementManager::refineFunction");
- int n_neigh;
+ FUNCNAME("RefinementManager::refineFunction()");
+
+ int n_neigh;
bool bound = false;
DegreeOfFreedom *edge[2];
@@ -76,7 +77,7 @@ namespace AMDiS {
::std::map<int, VertexVector*>::iterator it;
::std::map<int, VertexVector*>::iterator end = mesh->getPeriodicAssociations().end();
- while(edge[0] != NULL) {
+ while (edge[0] != NULL) {
periodicList = refineList->periodicSplit(edge,
next_edge,
&n_neigh,
@@ -134,8 +135,6 @@ namespace AMDiS {
/* and now refine the patch */
/****************************************************************************/
- //newDOF = refinePatch(edge, refineList, n_neigh, bound);
-
DELETE refineList;
return(el_info);
@@ -206,63 +205,10 @@ namespace AMDiS {
/****************************************************************************/
/* first refine the element */
/****************************************************************************/
-
- // check for periodic boundary
- // DegreeOfFreedom *periodicDOF[3] = {NULL, NULL, NULL};
- // RCNeighbourList *periodicRefineList = NULL;
- // int n_neigh_periodic;
-
- // if (neigh &&
- // (neigh->getDOF(0) != el->getDOF(1)) &&
- // (neigh->getDOF(0) != el->getDOF(0)))
- // {
- // periodicDOF[0] = mesh->getDOF(VERTEX);
- // mesh->incrementNumberOfVertices(1);
- // mesh->incrementNumberOfEdges(1);
- // if (mesh->getNumberOfDOFs(EDGE)) {
- // periodicDOF[1] = mesh->getDOF(EDGE);
- // periodicDOF[2] = mesh->getDOF(EDGE);
- // }
-
- // DegreeOfFreedom *edge[2] = {
- // const_cast<DegreeOfFreedom*>(el->getDOF(0)),
- // const_cast<DegreeOfFreedom*>(el->getDOF(1))
- // };
- // DegreeOfFreedom *periodicEdge[2] = {
- // const_cast<DegreeOfFreedom*>(neigh->getDOF(0)),
- // const_cast<DegreeOfFreedom*>(neigh->getDOF(1))
- // };
-
- // periodicRefineList = refineList->periodicSplit(edge,
- // &n_neigh,
- // periodicEdge,
- // &n_neigh_periodic);
-
- // // dynamic_cast<LeafDataPeriodic*>(neigh->getElementData(PERIODIC))->
- // // setNewDOF(dof[0][0]);
- // // dynamic_cast<LeafDataPeriodic*>(el->getElementData(PERIODIC))->
- // // setNewDOF(periodicDOF[0][0]);
-
- // // MSG("el %d new dof %d periodic dof %d\n",
- // // el->getIndex(),
- // // dof[0][0],
- // // periodicDOF[0][0]);
- // }
-
bisectTriangle(el, dof);
- // MSG("el %d -> %d %d\n",
- // el->getIndex(),
- // el->getChild(0)->getIndex(),
- // el->getChild(1)->getIndex());
-
if (neigh) {
- // if(periodicDOF[0]) {
- // dof[0] = periodicDOF[0];
- // dof[1] = periodicDOF[1];
- // dof[2] = periodicDOF[2];
- // } else {
DegreeOfFreedom *tmp = dof[1];
/****************************************************************************/
/* there is a neighbour; refine it also, but first exchange dof[1] and */
@@ -270,19 +216,11 @@ namespace AMDiS {
/****************************************************************************/
dof[1] = dof[2];
dof[2] = tmp;
- // }
bisectTriangle(neigh, dof);
-
- // MSG("neigh %d -> %d %d\n",
- // neigh->getIndex(),
- // neigh->getChild(0)->getIndex(),
- // neigh->getChild(1)->getIndex());
+ } else {
+ newCoords = true;
}
- else
- {
- newCoords = true;
- }
/****************************************************************************/
/* if there are functions to interpolate data to the finer grid, do so */
@@ -296,39 +234,30 @@ namespace AMDiS {
::std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
for(it = admin->beginDOFIndexed(); it != end; it++) {
(*it)->refineInterpol(*refineList, n_neigh);
- // if(periodicRefineList) {
- // (*it)->refineInterpol(*periodicRefineList, n_neigh_periodic);
- // }
}
}
- if (!mesh->queryCoarseDOFs())
- {
+ if (!mesh->queryCoarseDOFs()) {
+ /****************************************************************************/
+ /* if there should be no dof information on interior leaf elements remove */
+ /* dofs from edges and the centers of parents */
+ /****************************************************************************/
+ if (mesh->getNumberOfDOFs(EDGE)) {
+
+ int node = mesh->getNode(EDGE);
+
/****************************************************************************/
- /* if there should be no dof information on interior leaf elements remove */
- /* dofs from edges and the centers of parents */
+ /* the only DOF that can be freed is that in the refinement edge; all other*/
+ /* DOFs are handed on the children */
/****************************************************************************/
- if (mesh->getNumberOfDOFs(EDGE)) {
-
- int node = mesh->getNode(EDGE);
-
- /****************************************************************************/
- /* the only DOF that can be freed is that in the refinement edge; all other*/
- /* DOFs are handed on the children */
- /****************************************************************************/
-
- mesh->freeDOF(const_cast<int*>( el->getDOF(node+2)), EDGE);
- }
- if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) {
- refineList->removeDOFParents(n_neigh);
- // if(periodicRefineList) {
- // periodicRefineList->removeDOFParents(n_neigh_periodic);
- // }
- }
+
+ mesh->freeDOF(const_cast<int*>( el->getDOF(node+2)), EDGE);
}
-
- // if(periodicRefineList) DELETE periodicRefineList;
+ if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) {
+ refineList->removeDOFParents(n_neigh);
+ }
+ }
return dof[0][0];
}
@@ -358,7 +287,8 @@ namespace AMDiS {
el->setFirstChild(child[0]);
el->setSecondChild(child[1]);
- if (newMark > 0) doMoreRecursiveRefine = true;
+ if (newMark > 0)
+ doMoreRecursiveRefine = true;
/****************************************************************************/
/* vertex 2 is the newest vertex */
@@ -408,7 +338,7 @@ namespace AMDiS {
}
if (mesh->getNumberOfDOFs(CENTER)) {
- int node = mesh->getNode(CENTER);
+ int node = mesh->getNode(CENTER);
/****************************************************************************/
/* there are dofs at the barycenter of the triangles */
@@ -424,9 +354,7 @@ namespace AMDiS {
int *n_neigh)
{
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) {
@@ -436,7 +364,7 @@ namespace AMDiS {
/****************************************************************************/
opp_vertex = (*el_info)->getOppVertex(2);
- neigh_info = stack->traverseNeighbour2d(*el_info, 2);
+ ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(),
1));
neigh_info = refineFunction(neigh_info);
diff --git a/AMDiS/src/RefinementManager3d.cc b/AMDiS/src/RefinementManager3d.cc
index 57c0aadd277fd49fd29b60a57a842b547bec2e21..1a8c94e1ce85fee24f9d0a5119f36c16ce93aa80 100644
--- a/AMDiS/src/RefinementManager3d.cc
+++ b/AMDiS/src/RefinementManager3d.cc
@@ -22,8 +22,8 @@ namespace AMDiS {
DegreeOfFreedom *edge[2])
{
Tetrahedron *el = dynamic_cast<Tetrahedron*>(const_cast<Element*>( ref_list->getElement(index))), *child[2];
- int i, node;
- int el_type = ref_list->getType(index);
+ int i, node;
+ int el_type = ref_list->getType(index);
child[0] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
child[1] = dynamic_cast<Tetrahedron*>(mesh->createNewElement(el));
diff --git a/AMDiS/src/SystemVector.h b/AMDiS/src/SystemVector.h
index b6d4b95b65cf545d170b2053a29807f6cadde5ac..75f90eaaa650f844a140f0b4144ede1887c26dfe 100644
--- a/AMDiS/src/SystemVector.h
+++ b/AMDiS/src/SystemVector.h
@@ -168,8 +168,8 @@ namespace AMDiS {
* Returns sum of used vector sizes.
*/
inline int getUsedSize() const {
- int i, totalSize = 0, size = vectors.getSize();
- for(i = 0; i < size; i++) {
+ int totalSize = 0, size = vectors.getSize();
+ for (int i = 0; i < size; i++) {
totalSize += vectors[i]->getUsedSize();
}
return totalSize;
@@ -250,8 +250,8 @@ namespace AMDiS {
*/
inline SystemVector& operator=(const SystemVector& rhs) {
TEST_EXIT_DBG(rhs.vectors.getSize() == vectors.getSize())("invalied sizes\n");
- int i, size = vectors.getSize();
- for(i = 0; i < size; i++) {
+ int size = vectors.getSize();
+ for (int i = 0; i < size; i++) {
(*(vectors[i])) = (*(rhs.getDOFVector(i)));
}
return *this;
@@ -476,7 +476,7 @@ namespace AMDiS {
TEST_EXIT_DBG(size == matrix.getNumCols())("incompatible sizes\n");
#ifdef _OPENMP
-#pragma omp parallel for schedule(static, 1)
+#pragma omp parallel for schedule(static, 1) num_threads(size)
#endif
for (i = 0; i < size; i++) {
if (!add)
@@ -506,7 +506,7 @@ namespace AMDiS {
int i;
#ifdef _OPENMP
-#pragma omp parallel for schedule(static, 1)
+#pragma omp parallel for schedule(static, 1) num_threads(size)
#endif
for (i = 0; i < size; i++) {
axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));