Commit 29b6099b authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Not a bug fix, a feature fix.

parent ec1e26c5
...@@ -38,15 +38,13 @@ namespace AMDiS { ...@@ -38,15 +38,13 @@ namespace AMDiS {
{ {
FUNCNAME("Assembler::calculateElementMatrix()"); FUNCNAME("Assembler::calculateElementMatrix()");
if (remember && (factor != 1.0 || operat->uhOld)) { if (remember && (factor != 1.0 || operat->uhOld))
rememberElMat = true; rememberElMat = true;
}
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) { if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
initElement(elInfo); initElement(elInfo);
}
if (el != lastMatEl || !operat->isOptimized()) { if (el != lastMatEl || !operat->isOptimized()) {
if (rememberElMat) if (rememberElMat)
...@@ -84,16 +82,14 @@ namespace AMDiS { ...@@ -84,16 +82,14 @@ namespace AMDiS {
{ {
FUNCNAME("Assembler::calculateElementMatrix()"); FUNCNAME("Assembler::calculateElementMatrix()");
if (remember && ((factor != 1.0) || (operat->uhOld))) { if (remember && ((factor != 1.0) || (operat->uhOld)))
rememberElMat = true; rememberElMat = true;
}
Element *el = smallElInfo->getElement(); Element *el = smallElInfo->getElement();
lastVecEl = lastMatEl = NULL; lastVecEl = lastMatEl = NULL;
if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) { if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
initElement(smallElInfo, largeElInfo); initElement(smallElInfo, largeElInfo);
}
if (el != lastMatEl || !operat->isOptimized()) { if (el != lastMatEl || !operat->isOptimized()) {
if (rememberElMat) if (rememberElMat)
...@@ -149,15 +145,13 @@ namespace AMDiS { ...@@ -149,15 +145,13 @@ namespace AMDiS {
{ {
FUNCNAME("Assembler::calculateElementVector()"); FUNCNAME("Assembler::calculateElementVector()");
if (remember && factor != 1.0) { if (remember && factor != 1.0)
rememberElVec = true; rememberElVec = true;
}
Element *el = elInfo->getElement(); Element *el = elInfo->getElement();
if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) { if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
initElement(elInfo); initElement(elInfo);
}
if (el != lastVecEl || !operat->isOptimized()) { if (el != lastVecEl || !operat->isOptimized()) {
if (rememberElVec) if (rememberElVec)
...@@ -174,9 +168,9 @@ namespace AMDiS { ...@@ -174,9 +168,9 @@ namespace AMDiS {
ElementVector& vec = rememberElVec ? elementVector : userVec; ElementVector& vec = rememberElVec ? elementVector : userVec;
if (operat->uhOld && remember) { if (operat->uhOld && remember) {
matVecAssemble(elInfo, vec); matVecAssemble(elInfo, vec);
if (rememberElVec) { if (rememberElVec)
userVec += factor * elementVector; userVec += factor * elementVector;
}
return; return;
} }
...@@ -184,6 +178,7 @@ namespace AMDiS { ...@@ -184,6 +178,7 @@ namespace AMDiS {
firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec); firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
if (zeroOrderAssembler) if (zeroOrderAssembler)
zeroOrderAssembler->calculateElementVector(elInfo, vec); zeroOrderAssembler->calculateElementVector(elInfo, vec);
if (rememberElVec) if (rememberElVec)
userVec += factor * elementVector; userVec += factor * elementVector;
} }
...@@ -197,15 +192,13 @@ namespace AMDiS { ...@@ -197,15 +192,13 @@ namespace AMDiS {
{ {
FUNCNAME("Assembler::calculateElementVector()"); FUNCNAME("Assembler::calculateElementVector()");
if (remember && factor != 1.0) { if (remember && factor != 1.0)
rememberElVec = true; rememberElVec = true;
}
Element *el = mainElInfo->getElement(); Element *el = mainElInfo->getElement();
if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) { if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
initElement(auxElInfo); initElement(auxElInfo);
}
if (el != lastVecEl || !operat->isOptimized()) { if (el != lastVecEl || !operat->isOptimized()) {
if (rememberElVec) if (rememberElVec)
...@@ -228,9 +221,8 @@ namespace AMDiS { ...@@ -228,9 +221,8 @@ namespace AMDiS {
matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec); matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);
} }
if (rememberElVec) { if (rememberElVec)
userVec += factor * elementVector; userVec += factor * elementVector;
}
return; return;
} }
......
...@@ -32,7 +32,7 @@ namespace AMDiS { ...@@ -32,7 +32,7 @@ namespace AMDiS {
grdTmpVec1[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0); grdTmpVec1[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
grdTmpVec2[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0); grdTmpVec2[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
} }
}; }
BasisFunction::~BasisFunction() BasisFunction::~BasisFunction()
{ {
......
...@@ -40,30 +40,27 @@ namespace AMDiS { ...@@ -40,30 +40,27 @@ namespace AMDiS {
Element *el = el_info->getElement(); Element *el = el_info->getElement();
signed char mark; signed char mark;
if (el->getChild(0)) if (el->getChild(0)) {
{ /****************************************************************************/
/****************************************************************************/ /* interior node of the tree */
/* interior node of the tree */ /****************************************************************************/
/****************************************************************************/ mark = max(el->getChild(0)->getMark(), el->getChild(1)->getMark());
mark = max(el->getChild(0)->getMark(), el->getChild(1)->getMark()); el->setMark(std::min(mark + 1, 0));
el->setMark(std::min(mark + 1, 0)); } else {
} /****************************************************************************/
else /* leaf node of the tree */
{ /****************************************************************************/
/****************************************************************************/ if (el->getMark() < 0)
/* leaf node of the tree */ el->setMark(el->getMark() - 1);
/****************************************************************************/ }
if (el->getMark() < 0) el->setMark(el->getMark() - 1);
}
return 0; return 0;
} }
void CoarseningManager::spreadCoarsenMark() void CoarseningManager::spreadCoarsenMark()
{ {
traversePtr = this; traversePtr = this;
mesh->traverse(-1, mesh->traverse(-1, Mesh::CALL_EVERY_EL_POSTORDER, spreadCoarsenMarkFunction);
Mesh::CALL_EVERY_EL_POSTORDER,
spreadCoarsenMarkFunction);
} }
...@@ -93,7 +90,6 @@ namespace AMDiS { ...@@ -93,7 +90,6 @@ namespace AMDiS {
Flag CoarseningManager::coarsenMesh(Mesh *aMesh) Flag CoarseningManager::coarsenMesh(Mesh *aMesh)
{ {
int n_elements; int n_elements;
Flag flag = Mesh::CALL_EVERY_EL_POSTORDER | Mesh::FILL_NEIGH;
ElInfo *el_info; ElInfo *el_info;
mesh = aMesh; mesh = aMesh;
...@@ -102,18 +98,23 @@ namespace AMDiS { ...@@ -102,18 +98,23 @@ namespace AMDiS {
spreadCoarsenMark(); spreadCoarsenMark();
stack = NEW TraverseStack; stack = new TraverseStack;
do { do {
doMore = false; doMore = false;
el_info = stack->traverseFirst(mesh, -1, flag); el_info = stack->traverseFirst(mesh, -1,
Mesh::CALL_EVERY_EL_POSTORDER | Mesh::FILL_NEIGH);
while (el_info) { while (el_info) {
coarsenFunction(el_info); int idx = el_info->getElement()->getIndex();
// if (idx != 2288 && idx != 2283)
coarsenFunction(el_info);
el_info = stack->traverseNext(el_info); el_info = stack->traverseNext(el_info);
} }
} while (doMore); } while (doMore);
DELETE stack; delete stack;
cleanUpAfterCoarsen(); cleanUpAfterCoarsen();
......
...@@ -87,15 +87,15 @@ namespace AMDiS { ...@@ -87,15 +87,15 @@ namespace AMDiS {
/****************************************************************************/ /****************************************************************************/
/* get new dof on el at the midpoint of the coarsening edge */ /* get new dof on el at the midpoint of the coarsening edge */
/****************************************************************************/ /****************************************************************************/
if (!el->getDOF(node+2)) { if (!el->getDOF(node + 2)) {
el->setDOF(node+2, mesh->getDOF(EDGE)); el->setDOF(node + 2, mesh->getDOF(EDGE));
if (neigh) { if (neigh) {
neigh->setDOF(node+2, const_cast<int*>( el->getDOF(node+2))); neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node+2)));
} }
} }
} }
if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) { if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) {
coarsenList->addDOFParents(n_neigh); coarsenList->addDOFParents(n_neigh);
} }
...@@ -132,7 +132,7 @@ namespace AMDiS { ...@@ -132,7 +132,7 @@ namespace AMDiS {
int CoarseningManager2d::coarsenFunction(ElInfo *el_info) int CoarseningManager2d::coarsenFunction(ElInfo *el_info)
{ {
Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( el_info->getElement())); Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>(el_info->getElement()));
DegreeOfFreedom *edge[2]; DegreeOfFreedom *edge[2];
int n_neigh, bound = 0; int n_neigh, bound = 0;
RCNeighbourList coarse_list(2); RCNeighbourList coarse_list(2);
...@@ -165,11 +165,11 @@ namespace AMDiS { ...@@ -165,11 +165,11 @@ namespace AMDiS {
/****************************************************************************/ /****************************************************************************/
if (el->getDOF(0,0) < el->getDOF(1,0)) { if (el->getDOF(0,0) < el->getDOF(1,0)) {
edge[0] = const_cast<int*>( el->getDOF(0)); edge[0] = const_cast<int*>(el->getDOF(0));
edge[1] = const_cast<int*>( el->getDOF(1)); edge[1] = const_cast<int*>(el->getDOF(1));
} else { } else {
edge[1] = const_cast<int*>( el->getDOF(0)); edge[1] = const_cast<int*>(el->getDOF(0));
edge[0] = const_cast<int*>( el->getDOF(1)); edge[0] = const_cast<int*>(el->getDOF(1));
} }
coarse_list.setElement(0, el, true); coarse_list.setElement(0, el, true);
......
...@@ -89,7 +89,7 @@ namespace AMDiS { ...@@ -89,7 +89,7 @@ namespace AMDiS {
ERROR_EXIT("TODO\n"); ERROR_EXIT("TODO\n");
} }
void DOFAdmin::freeDOFIndex(int dof) { void DOFAdmin::freeDOFIndex(int dof) {
FUNCNAME("DOFAdmin::freeDOFIndex()"); FUNCNAME("DOFAdmin::freeDOFIndex()");
TEST_EXIT_DBG(usedCount > 0)("no dofs in use\n"); TEST_EXIT_DBG(usedCount > 0)("no dofs in use\n");
...@@ -98,19 +98,17 @@ namespace AMDiS { ...@@ -98,19 +98,17 @@ namespace AMDiS {
std::list<DOFIndexedBase*>::iterator di; std::list<DOFIndexedBase*>::iterator di;
std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end(); std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end();
for (di = dofIndexedList.begin(); di != end; ++di) { for (di = dofIndexedList.begin(); di != end; ++di)
(*di)->freeDOFContent(dof); (*di)->freeDOFContent(dof);
}
std::list<DOFContainer*>::iterator dc; std::list<DOFContainer*>::iterator dc;
std::list<DOFContainer*>::iterator dcend = dofContainerList.end(); std::list<DOFContainer*>::iterator dcend = dofContainerList.end();
for (dc = dofContainerList.begin(); dc != dcend; ++dc) { for (dc = dofContainerList.begin(); dc != dcend; ++dc)
(*dc)->freeDOFIndex(dof); (*dc)->freeDOFIndex(dof);
}
dofFree[dof] = true; dofFree[dof] = true;
if (static_cast<int>(firstHole) > dof) if (static_cast<int>(firstHole) > dof)
firstHole = dof; firstHole = dof;
...@@ -118,8 +116,6 @@ namespace AMDiS { ...@@ -118,8 +116,6 @@ namespace AMDiS {
holeCount++; holeCount++;
} }
/****************************************************************************/
int DOFAdmin::getDOFIndex() int DOFAdmin::getDOFIndex()
{ {
FUNCNAME("DOFAdmin::getDOFIndex()"); FUNCNAME("DOFAdmin::getDOFIndex()");
...@@ -160,9 +156,6 @@ namespace AMDiS { ...@@ -160,9 +156,6 @@ namespace AMDiS {
return(dof); return(dof);
} }
/****************************************************************************/
void DOFAdmin::enlargeDOFLists(int minsize) void DOFAdmin::enlargeDOFLists(int minsize)
{ {
FUNCNAME("DOFAdmin::enlargeDOFLists()"); FUNCNAME("DOFAdmin::enlargeDOFLists()");
...@@ -256,9 +249,6 @@ namespace AMDiS { ...@@ -256,9 +249,6 @@ namespace AMDiS {
ERROR("container not in list\n"); ERROR("container not in list\n");
} }
/****************************************************************************/
void DOFAdmin::compress(std::vector<DegreeOfFreedom> &new_dof) void DOFAdmin::compress(std::vector<DegreeOfFreedom> &new_dof)
{ {
FUNCNAME("DOFAdmin::compress()"); FUNCNAME("DOFAdmin::compress()");
...@@ -269,16 +259,14 @@ namespace AMDiS { ...@@ -269,16 +259,14 @@ namespace AMDiS {
if (holeCount < 1) return; if (holeCount < 1) return;
// vector to mark used dofs // vector to mark used dofs
for (int i = 0; i < size; i++) { for (int i = 0; i < size; i++)
new_dof[i] = -1; new_dof[i] = -1;
}
// mark used dofs // mark used dofs
DOFIteratorBase it(this, USED_DOFS); DOFIteratorBase it(this, USED_DOFS);
for (it.reset(); !it.end(); ++it) { for (it.reset(); !it.end(); ++it)
new_dof[it.getDOFIndex()] = 1; new_dof[it.getDOFIndex()] = 1;
}
int n = 0, last = 0; int n = 0, last = 0;
for (int i = 0; i < size; i++) { /* create a MONOTONE compress */ for (int i = 0; i < size; i++) { /* create a MONOTONE compress */
if (new_dof[i] == 1) { if (new_dof[i] == 1) {
...@@ -290,13 +278,12 @@ namespace AMDiS { ...@@ -290,13 +278,12 @@ namespace AMDiS {
TEST_EXIT_DBG(n == usedCount)("count %d != usedCount %d\n", n, usedCount); TEST_EXIT_DBG(n == usedCount)("count %d != usedCount %d\n", n, usedCount);
// mark used dofs in compressed dofFree // mark used dofs in compressed dofFree
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++)
dofFree[i] = false; dofFree[i] = false;
}
// mark unused dofs in compressed dofFree // mark unused dofs in compressed dofFree
for (int i = n; i < size; i++) { for (int i = n; i < size; i++)
dofFree[i] = true; dofFree[i] = true;
}
firstHole = n; firstHole = n;
holeCount = 0; holeCount = 0;
...@@ -313,18 +300,14 @@ namespace AMDiS { ...@@ -313,18 +300,14 @@ namespace AMDiS {
std::list<DOFIndexedBase*>::iterator di; std::list<DOFIndexedBase*>::iterator di;
std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end(); std::list<DOFIndexedBase*>::iterator end = dofIndexedList.end();
for (di = dofIndexedList.begin(); di != end; ++di) { for (di = dofIndexedList.begin(); di != end; ++di)
(*di)->compressDOFIndexed(first, last, new_dof); (*di)->compressDOFIndexed(first, last, new_dof);
};
std::list<DOFContainer*>::iterator dc; std::list<DOFContainer*>::iterator dc;
std::list<DOFContainer*>::iterator endc = dofContainerList.end(); std::list<DOFContainer*>::iterator endc = dofContainerList.end();
for (dc = dofContainerList.begin(); dc != endc; dc++) { for (dc = dofContainerList.begin(); dc != endc; dc++)
(*dc)->compressDOFContainer(n, new_dof); (*dc)->compressDOFContainer(n, new_dof);
};
return;
} }
void DOFAdmin::setNumberOfDOFs(int i,int v) { void DOFAdmin::setNumberOfDOFs(int i,int v) {
......
...@@ -98,12 +98,11 @@ namespace AMDiS { ...@@ -98,12 +98,11 @@ namespace AMDiS {
typedef traits::range_generator<major, Matrix>::type cursor_type; typedef traits::range_generator<major, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type; typedef traits::range_generator<nz, cursor_type>::type icursor_type;
for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) { std::cout.precision(10);
for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
Msg::print(" (%3d,%3d,%20.17lf)", row(*icursor), col(*icursor), value(*icursor)); std::cout << row(*icursor) << " " << col(*icursor) << " " << value(*icursor) << "\n";
Msg::print("\n");
}
} }
bool DOFMatrix::symmetric() bool DOFMatrix::symmetric()
...@@ -201,7 +200,6 @@ namespace AMDiS { ...@@ -201,7 +200,6 @@ namespace AMDiS {
} }
} }
// === Add element matrix to the global matrix using the indices mapping. === // === Add element matrix to the global matrix using the indices mapping. ===
DegreeOfFreedom row, col; DegreeOfFreedom row, col;
...@@ -219,6 +217,7 @@ namespace AMDiS { ...@@ -219,6 +217,7 @@ namespace AMDiS {
for (int j = 0; j < nCol; j++) { // for all columns for (int j = 0; j < nCol; j++) { // for all columns
col = colIndices[j]; col = colIndices[j];
entry = elMat[i][j]; entry = elMat[i][j];
if (add) if (add)
ins[row][col]+= sign * entry; ins[row][col]+= sign * entry;
else else
...@@ -233,33 +232,8 @@ namespace AMDiS { ...@@ -233,33 +232,8 @@ namespace AMDiS {
} }
void DOFMatrix::freeDOFContent(int index) void DOFMatrix::freeDOFContent(int index)
{ {}
if (matrix.nnz() == 0) return;
using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits= mtl::traits;
typedef base_matrix_type Matrix;
traits::row<Matrix>::type row(matrix);
traits::col<Matrix>::type col(matrix);
typedef traits::range_generator<major, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
cursor_type cursor = begin<major>(matrix);
// Jump directly to corresponding row or column
cursor+= index;
// Requires structural symmetry !!!
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) {
int my_row= row(*icursor), my_col= col(*icursor);
// Not very efficient (but general)
matrix.lvalue(my_row, my_col) = 0.0; // Need to call crop somewhere !!!! Peter
matrix.lvalue(my_col, my_row) = 0.0;
}
}
// Should work as before
void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound) void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
{ {
FUNCNAME("DOFMatrix::assemble()"); FUNCNAME("DOFMatrix::assemble()");
...@@ -268,13 +242,9 @@ namespace AMDiS { ...@@ -268,13 +242,9 @@ namespace AMDiS {
std::vector<Operator*>::iterator it = operators.begin(); std::vector<Operator*>::iterator it = operators.begin();
std::vector<double*>::iterator factorIt = operatorFactor.begin(); std::vector<double*>::iterator factorIt = operatorFactor.begin();
for (; it != operators.end(); ++it, ++factorIt) { for (; it != operators.end(); ++it, ++factorIt)
if ((*it)->getNeedDualTraverse() == false) { if ((*it)->getNeedDualTraverse() == false)
(*it)->getElementMatrix(elInfo, (*it)->getElementMatrix(elInfo, elementMatrix, *factorIt ? **factorIt : 1.0);
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
}