Commit 3e164d33 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

New umfpack parameter; Several bugfixes for multiple meshes.

parent 3ee7ae23
...@@ -20,6 +20,7 @@ AMDIS_LIB = -L$(AMDIS_DIR)/lib -lamdis -lcompositeFEM ...@@ -20,6 +20,7 @@ AMDIS_LIB = -L$(AMDIS_DIR)/lib -lamdis -lcompositeFEM
PNG_LIB = -lpng PNG_LIB = -lpng
UMFPACK_LIB = -L$(AMDIS_DIR)/lib/AMD/Lib -L$(AMDIS_DIR)/lib/UMFPACK/Lib UMFPACK_LIB = -L$(AMDIS_DIR)/lib/AMD/Lib -L$(AMDIS_DIR)/lib/UMFPACK/Lib
#UMFPACK_LIB = -L/u/witkowski/test/AMD/Lib -L/u/witkowski/test/UMFPACK/Lib
ifeq ($(strip $(USE_SERVER)), mars) ifeq ($(strip $(USE_SERVER)), mars)
UMFPACK_LIB += -lmkl -lumfpack -lamd UMFPACK_LIB += -lmkl -lumfpack -lamd
else else
......
...@@ -83,7 +83,7 @@ namespace AMDiS { ...@@ -83,7 +83,7 @@ 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();
...@@ -103,10 +103,11 @@ namespace AMDiS { ...@@ -103,10 +103,11 @@ namespace AMDiS {
return; return;
} }
} }
ElementMatrix& mat = rememberElMat ? elementMatrix : userMat; ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
if (secondOrderAssembler) { if (secondOrderAssembler) {
ERROR_EXIT("Da muss i noch ma testen!\n");
secondOrderAssembler->calculateElementMatrix(smallElInfo, mat); secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
ElementMatrix m(nRow, nRow); ElementMatrix m(nRow, nRow);
...@@ -134,17 +135,17 @@ namespace AMDiS { ...@@ -134,17 +135,17 @@ namespace AMDiS {
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree()); smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
ElementMatrix tmpMat(nRow, nRow); ElementMatrix tmpMat(nRow, nRow);
if (smallElInfo == colElInfo) if (smallElInfo == colElInfo)
tmpMat = m * mat; tmpMat = m * mat;
else else
tmpMat = mat * trans(m); tmpMat = mat * trans(m);
mat = tmpMat; mat = tmpMat;
} }
if (rememberElMat) if (rememberElMat && &userMat != &elementMatrix)
userMat += factor * elementMatrix; userMat += factor * elementMatrix;
} }
...@@ -178,7 +179,7 @@ namespace AMDiS { ...@@ -178,7 +179,7 @@ namespace AMDiS {
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;
} }
...@@ -189,9 +190,10 @@ namespace AMDiS { ...@@ -189,9 +190,10 @@ namespace AMDiS {
zeroOrderAssembler->calculateElementVector(elInfo, vec); zeroOrderAssembler->calculateElementVector(elInfo, vec);
if (rememberElVec) if (rememberElVec)
userVec += factor * elementVector; userVec += factor * elementVector;
} }
void Assembler::calculateElementVector(const ElInfo *mainElInfo, void Assembler::calculateElementVector(const ElInfo *mainElInfo,
const ElInfo *auxElInfo, const ElInfo *auxElInfo,
const ElInfo *smallElInfo, const ElInfo *smallElInfo,
...@@ -207,8 +209,8 @@ namespace AMDiS { ...@@ -207,8 +209,8 @@ namespace AMDiS {
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(smallElInfo, largeElInfo);
if (el != lastVecEl || !operat->isOptimized()) { if (el != lastVecEl || !operat->isOptimized()) {
if (rememberElVec) if (rememberElVec)
set_to_zero(elementVector); set_to_zero(elementVector);
...@@ -223,32 +225,30 @@ namespace AMDiS { ...@@ -223,32 +225,30 @@ namespace AMDiS {
ElementVector& vec = rememberElVec ? elementVector : userVec; ElementVector& vec = rememberElVec ? elementVector : userVec;
if (operat->uhOld && remember) { if (operat->uhOld && remember) {
if (smallElInfo->getLevel() == largeElInfo->getLevel())
if (smallElInfo->getLevel() == largeElInfo->getLevel()) {
matVecAssemble(auxElInfo, vec); matVecAssemble(auxElInfo, vec);
} else { else
matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec); matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);
}
if (rememberElVec) if (rememberElVec)
userVec += factor * elementVector; userVec += factor * elementVector;
return; return;
} }
ERROR_EXIT("Not yet implemented!\n"); ERROR_EXIT("Not yet implemented!\n");
// if (firstOrderAssemblerGrdPsi) {
// firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
// }
// if (zeroOrderAssembler) {
// zeroOrderAssembler->calculateElementVector(elInfo, vec);
// }
// if (rememberElVec) {
// axpy(factor, *elementVector, *userVec);
// }
#if 0
if (firstOrderAssemblerGrdPsi)
firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
if (zeroOrderAssembler)
zeroOrderAssembler->calculateElementVector(elInfo, vec);
if (rememberElVec)
axpy(factor, *elementVector, *userVec);
#endif
} }
void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec) void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
{ {
FUNCNAME("Assembler::matVecAssemble()"); FUNCNAME("Assembler::matVecAssemble()");
...@@ -258,20 +258,23 @@ namespace AMDiS { ...@@ -258,20 +258,23 @@ namespace AMDiS {
operat->uhOld->getLocalVector(el, uhOldLoc); operat->uhOld->getLocalVector(el, uhOldLoc);
if (el != lastMatEl) if (el != lastMatEl) {
set_to_zero(elementMatrix);
calculateElementMatrix(elInfo, elementMatrix); calculateElementMatrix(elInfo, elementMatrix);
}
for (int i = 0; i < nRow; i++) { for (int i = 0; i < nRow; i++) {
double val = 0.0; double val = 0.0;
for (int j = 0; j < nRow; j++) { for (int j = 0; j < nRow; j++)
val += elementMatrix[i][j] * uhOldLoc[j]; val += elementMatrix[i][j] * uhOldLoc[j];
}
vec[i] += val; vec[i] += val;
} }
delete [] uhOldLoc; delete [] uhOldLoc;
} }
void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo, void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
const ElInfo *smallElInfo, const ElInfo *largeElInfo, const ElInfo *smallElInfo, const ElInfo *largeElInfo,
ElementVector& vec) ElementVector& vec)
...@@ -289,35 +292,22 @@ namespace AMDiS { ...@@ -289,35 +292,22 @@ namespace AMDiS {
const BasisFunction *basFcts = rowFESpace->getBasisFcts(); const BasisFunction *basFcts = rowFESpace->getBasisFcts();
int nBasFcts = basFcts->getNumber(); int nBasFcts = basFcts->getNumber();
double *uhOldLoc = new double[nBasFcts]; std::vector<double> uhOldLoc(nBasFcts);
double *uhOldLoc2 = new double[nBasFcts];
operat->uhOld->getLocalVector(auxEl, uhOldLoc);
mtl::dense2D<double> m(nRow, nRow); operat->uhOld->getLocalVector(auxEl, &(uhOldLoc[0]));
smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
for (int i = 0; i < nBasFcts; i++) {
uhOldLoc2[i] = 0.0;
for (int j = 0; j < nBasFcts; j++)
uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
}
if (mainEl != lastMatEl) { if (mainEl != lastMatEl) {
set_to_zero(elementMatrix);
calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo,
elementMatrix); elementMatrix);
} }
for (int i = 0; i < nBasFcts; i++) { for (int i = 0; i < nBasFcts; i++) {
double val = 0.0; double val = 0.0;
for (int j = 0; j < nBasFcts; j++) for (int j = 0; j < nBasFcts; j++)
val += elementMatrix[i][j] * uhOldLoc2[j]; val += elementMatrix[i][j] * uhOldLoc[j];
vec[i] += val; vec[i] += val;
} }
delete [] uhOldLoc;
delete [] uhOldLoc2;
} }
......
...@@ -22,8 +22,10 @@ namespace AMDiS { ...@@ -22,8 +22,10 @@ namespace AMDiS {
status = SingleComponentInfo::EQ_SPACES_NO_AUX; status = SingleComponentInfo::EQ_SPACES_NO_AUX;
} else { } else {
status = SingleComponentInfo::EQ_SPACES_WITH_AUX; status = SingleComponentInfo::EQ_SPACES_WITH_AUX;
for (int i = 0; i < static_cast<int>(auxFeSpaces.size()); i++) {
if (auxFeSpaces[i]->getMesh() != rowFeSpace->getMesh()) { for (std::set<const FiniteElemSpace*>::iterator it = auxFeSpaces.begin();
it != auxFeSpaces.end(); ++it) {
if ((*it)->getMesh() != rowFeSpace->getMesh()) {
status = SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX; status = SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX;
break; break;
} }
...@@ -34,9 +36,10 @@ namespace AMDiS { ...@@ -34,9 +36,10 @@ namespace AMDiS {
status = SingleComponentInfo::DIF_SPACES_NO_AUX; status = SingleComponentInfo::DIF_SPACES_NO_AUX;
} else { } else {
status = SingleComponentInfo::DIF_SPACES_WITH_AUX; status = SingleComponentInfo::DIF_SPACES_WITH_AUX;
for (int i = 0; i < static_cast<int>(auxFeSpaces.size()); i++) { for (std::set<const FiniteElemSpace*>::iterator it = auxFeSpaces.begin();
if (auxFeSpaces[i]->getMesh() != rowFeSpace->getMesh() && it != auxFeSpaces.end(); ++it) {
auxFeSpaces[i]->getMesh() != colFeSpace->getMesh()) { if ((*it)->getMesh() != rowFeSpace->getMesh() &&
(*it)->getMesh() != colFeSpace->getMesh()) {
status = SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX; status = SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX;
break; break;
} }
...@@ -72,12 +75,12 @@ namespace AMDiS { ...@@ -72,12 +75,12 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) { for (int i = 0; i < nComponents; i++) {
if (matrixComponents[row][i].getColFeSpace() != rowFeSpace) if (matrixComponents[row][i].getColFeSpace() != rowFeSpace)
return matrixComponents[row][i].getColFeSpace(); return matrixComponents[row][i].getColFeSpace();
if (matrixComponents[row][i].getAuxFeSpace(0) != rowFeSpace) if (matrixComponents[row][i].getAuxFeSpace() != rowFeSpace)
return matrixComponents[row][i].getAuxFeSpace(0); return matrixComponents[row][i].getAuxFeSpace();
} }
if (vectorComponents[row].getAuxFeSpace(0) != rowFeSpace) if (vectorComponents[row].getAuxFeSpace() != rowFeSpace)
return vectorComponents[row].getAuxFeSpace(0); return vectorComponents[row].getAuxFeSpace();
return NULL; return NULL;
} }
......
...@@ -22,7 +22,7 @@ ...@@ -22,7 +22,7 @@
#ifndef AMDIS_COMPONENTTRAVERSEINFO_H #ifndef AMDIS_COMPONENTTRAVERSEINFO_H
#define AMDIS_COMPONENTTRAVERSEINFO_H #define AMDIS_COMPONENTTRAVERSEINFO_H
#include <vector> #include <set>
#include "Global.h" #include "Global.h"
#include "FiniteElemSpace.h" #include "FiniteElemSpace.h"
...@@ -34,7 +34,6 @@ namespace AMDiS { ...@@ -34,7 +34,6 @@ namespace AMDiS {
SingleComponentInfo() SingleComponentInfo()
: rowFeSpace(NULL), : rowFeSpace(NULL),
colFeSpace(NULL), colFeSpace(NULL),
auxFeSpaces(0),
status(0) status(0)
{} {}
...@@ -44,14 +43,14 @@ namespace AMDiS { ...@@ -44,14 +43,14 @@ namespace AMDiS {
colFeSpace = col; colFeSpace = col;
} }
void setAuxFeSpaces(std::vector<const FiniteElemSpace*> feSpaces) void setAuxFeSpaces(std::set<const FiniteElemSpace*> feSpaces)
{ {
auxFeSpaces = feSpaces; auxFeSpaces = feSpaces;
} }
void addAuxFeSpace(const FiniteElemSpace *fe) void addAuxFeSpace(const FiniteElemSpace *fe)
{ {
auxFeSpaces.push_back(fe); auxFeSpaces.insert(fe);
} }
bool hasFeSpace() bool hasFeSpace()
...@@ -76,14 +75,16 @@ namespace AMDiS { ...@@ -76,14 +75,16 @@ namespace AMDiS {
return colFeSpace; return colFeSpace;
} }
const FiniteElemSpace *getAuxFeSpace(int i) const FiniteElemSpace *getAuxFeSpace()
{ {
return ((i < static_cast<int>(auxFeSpaces.size())) ? auxFeSpaces[i] : NULL); FUNCNAME("SingleComponentInfo::getAuxFeSpace()");
}
void setAuxFeSpace(const FiniteElemSpace* fe, int pos) TEST_EXIT_DBG(auxFeSpaces.size() <= 1)("More than one aux FE space!\n");
{
auxFeSpaces[pos] = fe; if (auxFeSpaces.size() == 1)
return (*(auxFeSpaces.begin()));
return NULL;
} }
int getStatus() int getStatus()
...@@ -96,7 +97,7 @@ namespace AMDiS { ...@@ -96,7 +97,7 @@ namespace AMDiS {
FiniteElemSpace *colFeSpace; FiniteElemSpace *colFeSpace;
std::vector<const FiniteElemSpace*> auxFeSpaces; std::set<const FiniteElemSpace*> auxFeSpaces;
/// Status of the component, see the possible status flags below. /// Status of the component, see the possible status flags below.
int status; int status;
...@@ -177,12 +178,12 @@ namespace AMDiS { ...@@ -177,12 +178,12 @@ namespace AMDiS {
const FiniteElemSpace* getAuxFeSpace(int row, int col) const FiniteElemSpace* getAuxFeSpace(int row, int col)
{ {
return matrixComponents[row][col].getAuxFeSpace(0); return matrixComponents[row][col].getAuxFeSpace();
} }
const FiniteElemSpace* getAuxFeSpace(int row) const FiniteElemSpace* getAuxFeSpace(int row)
{ {
return vectorComponents[row].getAuxFeSpace(0); return vectorComponents[row].getAuxFeSpace();
} }
/** \brief /** \brief
......
...@@ -808,9 +808,11 @@ namespace AMDiS { ...@@ -808,9 +808,11 @@ namespace AMDiS {
return; return;
set_to_zero(this->elementVector); set_to_zero(this->elementVector);
bool addVector = false;
if (op) { if (op) {
op->getElementVector(elInfo, this->elementVector); op->getElementVector(elInfo, this->elementVector);
addVector = true;
} else { } else {
std::vector<Operator*>::iterator it; std::vector<Operator*>::iterator it;
std::vector<double*>::iterator factorIt; std::vector<double*>::iterator factorIt;
...@@ -818,12 +820,15 @@ namespace AMDiS { ...@@ -818,12 +820,15 @@ namespace AMDiS {
for (it = this->operators.begin(), factorIt = this->operatorFactor.begin(); for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
it != this->operators.end(); it != this->operators.end();
++it, ++factorIt) ++it, ++factorIt)
if ((*it)->getNeedDualTraverse() == false) if ((*it)->getNeedDualTraverse() == false) {
(*it)->getElementVector(elInfo, this->elementVector, (*it)->getElementVector(elInfo, this->elementVector,
*factorIt ? **factorIt : 1.0); *factorIt ? **factorIt : 1.0);
addVector = true;
}
} }
addElementVector(factor, this->elementVector, bound, elInfo); if (addVector)
addElementVector(factor, this->elementVector, bound, elInfo);
} }
template<> template<>
...@@ -839,27 +844,25 @@ namespace AMDiS { ...@@ -839,27 +844,25 @@ namespace AMDiS {
return; return;
set_to_zero(this->elementVector); set_to_zero(this->elementVector);
bool addVector = false;
if (op) {
ERROR_EXIT("TODO"); TEST_EXIT(!op)("Not yet implemented!\n");
// op->getElementVector(mainElInfo, this->elementVector);
} else { std::vector<Operator*>::iterator it;
std::vector<Operator*>::iterator it; std::vector<double*>::iterator factorIt;
std::vector<double*>::iterator factorIt; for (it = this->operators.begin(), factorIt = this->operatorFactor.begin();
for (it = this->operators.begin(), factorIt = this->operatorFactor.begin(); it != this->operators.end();
it != this->operators.end(); ++it, ++factorIt)
++it, ++factorIt) { if ((*it)->getNeedDualTraverse()) {
if ((*it)->getNeedDualTraverse()) { (*it)->getElementVector(mainElInfo, auxElInfo,
(*it)->getElementVector(mainElInfo, auxElInfo, smallElInfo, largeElInfo,
smallElInfo, largeElInfo, this->elementVector,
this->elementVector, *factorIt ? **factorIt : 1.0);
*factorIt ? **factorIt : 1.0); addVector = true;
}
} }
}
if (addVector)
addElementVector(factor, this->elementVector, bound, mainElInfo); addElementVector(factor, this->elementVector, bound, mainElInfo);
} }
......
...@@ -19,28 +19,16 @@ namespace AMDiS { ...@@ -19,28 +19,16 @@ namespace AMDiS {
{0.0, 0.0, 1.0}}; {0.0, 0.0, 1.0}};
mtl::dense2D<double> ElInfo2d::mat_d1(mat_d1_val); mtl::dense2D<double> ElInfo2d::mat_d1(mat_d1_val);
/*
double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5}, double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 1.0, 0.5},
{0.0, 0.0, 0.5}, {0.0, 0.0, 0.5},
{1.0, 0.0, 0.0}}; {1.0, 0.0, 0.0}};
*/
double ElInfo2d::mat_d1_left_val[3][3] = {{0.0, 0.0, 1.0},
{1.0, 0.0, 0.0},
{0.5, 0.5, 0.0}};
mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val); mtl::dense2D<double> ElInfo2d::mat_d1_left(mat_d1_left_val);
/*
double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5}, double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 0.0, 0.5},
{1.0, 0.0, 0.5}, {1.0, 0.0, 0.5},
{0.0, 1.0, 0.0}}; {0.0, 1.0, 0.0}};
*/
double ElInfo2d::mat_d1_right_val[3][3] = {{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0},
{0.5, 0.5, 0.0}};
mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val); mtl::dense2D<double> ElInfo2d::mat_d1_right(mat_d1_right_val);
......
...@@ -24,6 +24,8 @@ namespace AMDiS { ...@@ -24,6 +24,8 @@ namespace AMDiS {
FUNCNAME("Estimator::estimate()"); FUNCNAME("Estimator::estimate()");
bool dualTraverse = false; bool dualTraverse = false;
/*
for (unsigned int i = 0; i < matrix.size(); i++) { for (unsigned int i = 0; i < matrix.size(); i++) {
TEST_EXIT(traverseInfo.getStatus(row, i) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX) TEST_EXIT(traverseInfo.getStatus(row, i) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX)
("Not yet implemented!\n"); ("Not yet implemented!\n");
...@@ -51,7 +53,11 @@ namespace AMDiS { ...@@ -51,7 +53,11 @@ namespace AMDiS {
auxFeSpace->getBasisFcts()->getDegree()) auxFeSpace->getBasisFcts()->getDegree())
("Mh, do you really want to do this? Think about it ...\n"); ("Mh, do you really want to do this? Think about it ...\n");
} }
*/