Skip to content

GitLab

Commit 44e4efd1 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* hope, it works, no more bugs at all

parent 44677772
Hide whitespace changes
Inline Side-by-side
Showing with 162 additions and 140 deletions
AMDiS/src/CoarseningManager.cc
View file @ 44e4efd1
......@@ -107,11 +107,10 @@ namespace AMDiS {
do {
doMore = false;
el_info = stack->traverseFirst(mesh, -1, flag);
while (el_info)
{
coarsenFunction(el_info);
el_info = stack->traverseNext(el_info);
}
while (el_info) {
coarsenFunction(el_info);
el_info = stack->traverseNext(el_info);
}
} while (doMore);
DELETE stack;
......
AMDiS/src/ComponentTraverseInfo.h
View file @ 44e4efd1
......@@ -70,7 +70,11 @@ namespace AMDiS {
}
const FiniteElemSpace *getAuxFESpace(int i) {
return ((i <= static_cast<int>(auxFESpaces.size())) ? auxFESpaces[i] : NULL);
return ((i < static_cast<int>(auxFESpaces.size())) ? auxFESpaces[i] : NULL);
}
void setAuxFESpace(const FiniteElemSpace* fe, int pos) {
auxFESpaces[pos] = fe;
}
int getStatus() {
......@@ -164,6 +168,22 @@ namespace AMDiS {
return vectorComponents[row].getAuxFESpace(0);
}
void fakeFESpace(const FiniteElemSpace *feOld,
const FiniteElemSpace *feNew)
{
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (matrixComponents[i][j].getAuxFESpace(0) == feOld) {
matrixComponents[i][j].setAuxFESpace(feNew, 0);
}
}
if (vectorComponents[i].getAuxFESpace(0) == feOld) {
vectorComponents[i].setAuxFESpace(feNew, 0);
}
}
}
protected:
int nComponents;
......
AMDiS/src/DOFMatrix.cc
View file @ 44e4efd1
......@@ -441,11 +441,11 @@ namespace AMDiS {
std::vector<Operator*>::iterator it = operators.begin();
std::vector<double*>::iterator factorIt = operatorFactor.begin();
for (; it != operators.end(); ++it, ++factorIt) {
if ((*it)->getNeedDualTraverse() == false) {
(*it)->getElementMatrix(elInfo,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
if ((*it)->getNeedDualTraverse() == false) {
(*it)->getElementMatrix(elInfo,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
}
addElementMatrix(factor, *elementMatrix, bound);
......
AMDiS/src/DOFVector.cc
View file @ 44e4efd1
......@@ -16,6 +16,12 @@ namespace AMDiS {
(const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseRestr(this, &list, n);
break;
case COARSE_INTERPOL:
if (feSpace == NULL || !feSpace) {
ERROR_EXIT("ERR1\n");
}
if (feSpace->getBasisFcts() == NULL || !(feSpace->getBasisFcts())) {
ERROR_EXIT("ERR2\n");
}
(const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n);
break;
default:
......
AMDiS/src/DOFVector.h
View file @ 44e4efd1
......@@ -498,7 +498,7 @@ namespace AMDiS {
*/
inline double H1Norm(Quadrature* q = NULL) const {
return sqrt(H1NormSquare());
};
}
/** \brief
* Calculates square of H1 norm of this DOFVector
......
AMDiS/src/DOFVector.hh
View file @ 44e4efd1
......@@ -744,6 +744,7 @@ namespace AMDiS {
DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
{
feSpace = rhs.feSpace;
DOFVectorBase<T>::feSpace = rhs.feSpace;
this->nBasFcts = rhs.nBasFcts;
vec = rhs.vec;
this->elementVector = new ElementVector(this->nBasFcts);
......
AMDiS/src/Operator.h
View file @ 44e4efd1
......@@ -1575,20 +1575,14 @@ namespace AMDiS {
AbstractFunction<double, double> *af,
WorldVector<double> *wv);
/** \brief
* Implementation of \ref OperatorTerm::initElement().
*/
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
/** \brief
* Implements FirstOrderTerm::getLb().
*/
/// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
/** \brief
* Implements FirstOrderTerm::eval().
*/
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
......@@ -1597,23 +1591,16 @@ namespace AMDiS {
double factor) const;
protected:
/** \brief
* DOFVector to be evaluated at quadrature points.
*/
/// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec;
/** \brief
* v at quadrature points.
*/
/// v at quadrature points.
double *vecAtQPs;
/** \brief
* Function f.
*/
/// Function f.
AbstractFunction<double, double> *f;
/**
*/
///
WorldVector<double> *b;
};
......
AMDiS/src/ProblemVec.cc
View file @ 44e4efd1
......@@ -729,7 +729,6 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
// The DOFMatrix which should be assembled (or not, if assembleMatrix
......@@ -780,18 +779,18 @@ namespace AMDiS {
} else if (traverseInfo.getStatus(i) == SingleComponentInfo::EQ_SPACES_WITH_DIF_AUX) {
// Row fe space and col fe space are both equal, but right hand side has at
// least one another aux fe space. In this case, do not assemble the rhs,
// this will be done afterwards.
// least one another aux fe space.
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
assembleOnDifMeshes2(componentSpaces[1], componentSpaces[0],
assembleOnDifMeshes2(componentSpaces[i],
traverseInfo.getAuxFESpace(i, j),
assembleFlag,
NULL,
rhs->getDOFVector(1));
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
} else {
ERROR_EXIT("Possible? If yes, not yet implemented!\n");
......@@ -803,12 +802,12 @@ namespace AMDiS {
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
assembleOnDifMeshes2(componentSpaces[i],
traverseInfo.getAuxFESpace(i, j),
assembleFlag,
assembleMatrix ? matrix : NULL,
NULL);
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
} else if (traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_NO_AUX ||
traverseInfo.getStatus(i, j) == SingleComponentInfo::DIF_SPACES_WITH_AUX) {
......@@ -888,7 +887,8 @@ namespace AMDiS {
void ProblemVec::addMatrixOperator(Operator *op,
int i, int j,
double *factor,
double *estFactor)
double *estFactor,
bool fake)
{
FUNCNAME("ProblemVec::addMatrixOperator()");
......@@ -900,34 +900,41 @@ namespace AMDiS {
setBoundaryConditionMap((*systemMatrix)[i][i]->getBoundaryManager()->
getBoundaryConditionMap());
}
(*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces());
(*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
if ((op->getAuxFESpaces())[k] != componentSpaces[i] ||
(op->getAuxFESpaces())[k] != componentSpaces[j]) {
op->setNeedDualTraverse(true);
break;
}
}
if (!fake) {
traverseInfo.getMatrix(i, j).setAuxFESpaces(op->getAuxFESpaces());
for (int k = 0; k < static_cast<int>(op->getAuxFESpaces().size()); k++) {
if ((op->getAuxFESpaces())[k]->getMesh() != componentSpaces[i]->getMesh() ||
(op->getAuxFESpaces())[k]->getMesh() != componentSpaces[j]->getMesh()) {
op->setNeedDualTraverse(true);
break;
}
}
}
}
void ProblemVec::addVectorOperator(Operator *op, int i,
double *factor,
double *estFactor)
double *estFactor,
bool fake)
{
FUNCNAME("ProblemVec::addVectorOperator()");
rhs->getDOFVector(i)->addOperator(op, factor, estFactor);
traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces());
for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
if ((op->getAuxFESpaces())[j] != componentSpaces[i]) {
op->setNeedDualTraverse(true);
break;
}
}
if (!fake) {
traverseInfo.getVector(i).setAuxFESpaces(op->getAuxFESpaces());
for (int j = 0; j < static_cast<int>(op->getAuxFESpaces().size()); j++) {
if ((op->getAuxFESpaces())[j]->getMesh() != componentSpaces[i]->getMesh()) {
op->setNeedDualTraverse(true);
break;
}
}
}
}
void ProblemVec::addDirichletBC(BoundaryType type, int system,
......
AMDiS/src/ProblemVec.h
View file @ 44e4efd1
......@@ -216,12 +216,14 @@ namespace AMDiS {
/// Adds an Operator to \ref A.
void addMatrixOperator(Operator *op, int i, int j,
double *factor = NULL,
double *estFactor = NULL);
double *estFactor = NULL,
bool fake = false);
/// Adds an Operator to \ref rhs.
void addVectorOperator(Operator *op, int i,
double *factor = NULL,
double *estFactor = NULL);
double *estFactor = NULL,
bool fake = false);
/// Adds dirichlet boundary conditions.
virtual void addDirichletBC(BoundaryType type, int system,
......
AMDiS/src/Quadrature.h
View file @ 44e4efd1
......@@ -84,19 +84,13 @@ namespace AMDiS {
{}
public:
/** \brief
* Copy constructor
*/
/// Copy constructor
Quadrature(const Quadrature&);
/** \brief
* Destructor
*/
/// Destructor
~Quadrature();
/** \brief
* Returns a Quadrature for dimension dim exact for degree degree.
*/
/// Returns a Quadrature for dimension dim exact for degree degree.
static Quadrature *provideQuadrature(int dim, int degree);
/** \brief
......@@ -116,42 +110,40 @@ namespace AMDiS {
*/
inline const std::string& getName() {
return name;
};
}
/** \brief
* Returns \ref n_points
*/
/// Returns \ref n_points
inline int getNumPoints() const {
return n_points;
};
}
/** \brief
* Returns \ref w[p]
*/
inline double getWeight(int p) const {
return w[p];
};
}
/** \brief
* Returns \ref w.
*/
inline double* getWeight() const {
return w;
};
}
/** \brief
* Returns \ref dim
*/
inline int getDim() const {
return dim;
};
}
/** \brief
* Returns \ref degree
*/
inline int getDegree() const {
return degree;
};
}
/** \brief
* Returns a pointer to a vector storing the values of a doubled valued
......@@ -184,7 +176,7 @@ namespace AMDiS {
*/
inline double getLambda(int a, int b) const {
return (lambda ? (*lambda)[a][b] : 0.0);
};
}
/** \brief
* Returns \ref lambda[a] which is a DimVec<double> containing the
......@@ -192,14 +184,14 @@ namespace AMDiS {
*/
inline const DimVec<double>& getLambda(int a) const {
return (*lambda)[a];
};
}
/** \brief
* Returns \ref lambda which is a VectorOfFixvecs<DimVec<double> >.
*/
VectorOfFixVecs<DimVec<double> > *getLambda() const {
return lambda;
};
}
public:
......@@ -224,9 +216,7 @@ namespace AMDiS {
*/
int dim;
/** \brief
* Number of quadrature points
*/
/// Number of quadrature points
int n_points;
/** \brief
......@@ -378,7 +368,7 @@ namespace AMDiS {
D2Phi(NULL),
quadrature(quad),
basisFunctions(basFcts)
{};
{}
/** \brief
* Copy constructor
......@@ -424,21 +414,21 @@ namespace AMDiS {
ERROR_EXIT("invalid flag\n");
return false;
};
}
/** \brief
* Returns \ref quadrature
*/
inline const Quadrature* getQuadrature() const {
return quadrature;
};
}
/** \brief
* Returns \ref max_points
*/
inline int getMaxQuadPoints() {
return max_points;
};
}
/** \brief
* Returns (*\ref D2Phi)[q][i][j][m]
......@@ -455,63 +445,63 @@ namespace AMDiS {
*/
inline const double getGradient(int q, int i ,int j) const {
return (grdPhi) ? (*grdPhi)[q][i][j] : 0.0;
};
}
/** \brief
* Returns (*\ref grdPhi)[q]
*/
inline VectorOfFixVecs<DimVec<double> >* getGradient(int q) const {
return (grdPhi) ? &((*grdPhi)[q]) : NULL;
};
}
/** \brief
* Returns \ref phi[q][i]
*/
inline const double getPhi(int q, int i) const {
return (phi) ? phi[q][i] : 0.0;
};
}
/** \brief
* Returns \ref phi[q]
*/
inline const double *getPhi(int q) const {
return (phi) ? phi[q] : NULL;
};
}
/** \brief
* Returns \ref quadrature ->integrateStdSimplex(f)
*/
inline double integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f) {
return quadrature->integrateStdSimplex(f);
};
}
/** \brief
* Returns \ref quadrature ->getNumPoints()
*/
inline int getNumPoints() const {
return quadrature->getNumPoints();
};
}
/** \brief
* Returns \ref quadrature ->getWeight(p)
*/
inline double getWeight(int p) const {
return quadrature->getWeight(p);
};
}
/** \brief
* Returns \ref quadrature ->getDim()
*/
inline int getDim() const {
return quadrature->getDim();
};
}
/** \brief
* Returns \ref quadrature ->getDegree()
*/
inline int getDegree() const {
return quadrature->getDegree();
};
}
/** \brief
* Returns \ref quadrature ->grdFAtQp(f, vec)
......@@ -522,7 +512,7 @@ namespace AMDiS {
WorldVector<double>* vec) const
{
return quadrature->grdFAtQp(f, vec);
};
}
/** \brief
* Returns \ref quadrature ->fAtQp(f, vec)
......@@ -531,28 +521,28 @@ namespace AMDiS {
DimVec<double> >& f, double *vec) const
{
return quadrature->fAtQp(f, vec);
};
}
/** \brief
* Returns \ref quadrature ->getLambda(a,b)
*/
inline double getLambda(int a,int b) const {
return quadrature->getLambda(a,b);
};
}
/** \brief
* Returns \ref quadrature ->getLambda(a)
*/
inline const DimVec<double>& getLambda(int a) const {
return quadrature->getLambda(a);
};
}
/** \brief
* Returns \ref basisFunctions
*/
inline BasisFunction* getBasisFunctions() const {
return basisFunctions;
};
}
protected:
/** \brief
......
AMDiS/src/ResidualEstimator.cc
View file @ 44e4efd1
......@@ -178,13 +178,14 @@ namespace AMDiS {
#if 1
void ResidualEstimator::estimateElement(ElInfo *elInfo)
{
{
FUNCNAME("ResidualEstimator::estimateElement()");
TEST_EXIT_DBG(nSystems > 0)("no system set\n");
double val = 0.0;
std::vector<Operator*>::iterator it;
std::vector<double*>::iterator itfac;
Element *el = elInfo->getElement();
double det = elInfo->getDet();
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
......@@ -201,17 +202,20 @@ namespace AMDiS {
continue;
// init assemblers
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it) {
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
++it, ++itfac) {
if (**itfac != 0.0) {
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
}
}
if (C0) {
for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd();
++it) {
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
}
}
......@@ -238,20 +242,23 @@ namespace AMDiS {
}
if (C0) {
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it) {
if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
uhQP = GET_MEMORY(double, nPoints);
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
}
if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
grdUh_qp = new WorldVector<double>[nPoints];
uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
}
if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) {
D2uhqp = new WorldMatrix<double>[nPoints];
uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
++it, ++itfac) {
if (**itfac != 0.0) {
if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
uhQP = GET_MEMORY(double, nPoints);
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
}
if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
grdUh_qp = new WorldVector<double>[nPoints];
uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
}
if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) {
D2uhqp = new WorldMatrix<double>[nPoints];