diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 711aec0a60874c3ea0ffb96b98e906b08a5dfe67..44413fe2f81345cda4c6500785b22bc787ad4939 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -84,6 +84,8 @@ namespace AMDiS {
void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *userMat,
double factor)
{
@@ -122,14 +124,24 @@ namespace AMDiS {
ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;
- if (secondOrderAssembler)
- secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, mat);
- if (firstOrderAssemblerGrdPsi)
- firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo, mat);
- if (firstOrderAssemblerGrdPhi)
- firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo, mat);
+ if (secondOrderAssembler) {
+ secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo, mat);
+ }
+
+ if (firstOrderAssemblerGrdPsi) {
+ firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo, mat);
+ }
+
+ if (firstOrderAssemblerGrdPhi) {
+ firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo, mat);
+ }
+
if (zeroOrderAssembler) {
- zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, mat);
+ zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo, mat);
}
if (rememberElMat && userMat) {
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index 6b9170d0e48c5fa3d6e0f23ba57888a4bf6c168d..cd0edc0a5563270fb9be484ae3e62511d073c329 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -68,7 +68,8 @@ namespace AMDiS {
const FiniteElemSpace *rowFESpace,
const FiniteElemSpace *colFESpace = NULL);
- virtual ~Assembler() {};
+ virtual ~Assembler()
+ {}
ElementMatrix *initElementMatrix(ElementMatrix *elMat,
const ElInfo *rowElInfo,
@@ -87,6 +88,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *userMat,
double factor = 1.0);
@@ -102,49 +105,49 @@ namespace AMDiS {
*/
inline const FiniteElemSpace* getRowFESpace() {
return rowFESpace;
- };
+ }
/** \brief
* Returns \ref colFESpace.
*/
inline const FiniteElemSpace* getColFESpace() {
return colFESpace;
- };
+ }
/** \brief
* Returns \ref nRow.
*/
inline int getNRow() {
return nRow;
- };
+ }
/** \brief
* Returns \ref nCol.
*/
inline int getNCol() {
return nCol;
- };
+ }
/** \brief
* Sets \ref rememberElMat.
*/
inline void rememberElementMatrix(bool rem) {
rememberElMat = rem;
- };
+ }
/** \brief
* Sets \ref rememberElVec.
*/
inline void rememberElementVector(bool rem) {
rememberElVec = rem;
- };
+ }
/** \brief
* Returns \ref zeroOrderAssembler.
*/
inline ZeroOrderAssembler* getZeroOrderAssembler() {
return zeroOrderAssembler;
- };
+ }
/** \brief
* Returns \ref firstOrderAssemblerGrdPsi or \ref firstOrderAssemblerGrdPhi
@@ -155,21 +158,21 @@ namespace AMDiS {
return (type == GRD_PSI) ?
firstOrderAssemblerGrdPsi :
firstOrderAssemblerGrdPhi;
- };
+ }
/** \brief
* Returns \ref secondOrderAssembler.
*/
inline SecondOrderAssembler* getSecondOrderAssembler() {
return secondOrderAssembler;
- };
+ }
/** \brief
* Returns \ref operat;
*/
inline Operator* getOperator() {
return operat;
- };
+ }
/** \brief
* Initialisation for the given ElInfo. The call is deligated to
@@ -223,7 +226,7 @@ namespace AMDiS {
lastVecEl = lastMatEl = NULL;
lastTraverseId = ElInfo::traverseId[omp_get_thread_num()];
}
- };
+ }
/** \brief
* Checks whether quadratures for sub assemblers are already set.
diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index e11450a3ae524399930e1753aa80ddbf726de14a..80aaa3f0c5f4062bfbad0f942dc44ceb9c5feda6 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -276,21 +276,21 @@ namespace AMDiS {
*/
inline BasFctType *getPhi(int i) const {
return (*phi)[i];
- };
+ }
/** \brief
* Returns the gradient of the i-th local basis function
*/
inline GrdBasFctType *getGrdPhi(int i) const {
return (*grdPhi)[i];
- };
+ }
/** \brief
* Returns the second derivative of the i-th local basis function
*/
inline D2BasFctType *getD2Phi(int i) const {
return (*d2Phi)[i];
- };
+ }
/** \brief
* Approximates the L2 scalar products of a given function with all basis
@@ -314,45 +314,53 @@ namespace AMDiS {
*/
virtual void l2ScpFctBas(Quadrature*,
AbstractFunction<double, WorldVector<double> >* /*f*/,
- DOFVector<double>* /*fh*/) {};
+ DOFVector<double>* /*fh*/)
+ {}
/** \brief
* WorldVector<double> valued l2ScpFctBas function
*/
virtual void l2ScpFctBas(Quadrature* ,
AbstractFunction<WorldVector<double>, WorldVector<double> >* /*f*/,
- DOFVector<WorldVector<double> >* /*fh*/) {};
+ DOFVector<WorldVector<double> >* /*fh*/)
+ {}
/** \brief
* Interpolates a DOFIndexed<double> after refinement
*/
- virtual void refineInter(DOFIndexed<double> *, RCNeighbourList*, int){};
+ virtual void refineInter(DOFIndexed<double> *, RCNeighbourList*, int)
+ {}
/** \brief
* Interpolates a DOFIndexed<double> after coarsening
*/
- virtual void coarseInter(DOFIndexed<double> *, RCNeighbourList*, int){};
+ virtual void coarseInter(DOFIndexed<double> *, RCNeighbourList*, int)
+ {}
/** \brief
* Restricts a DOFIndexed<double> after coarsening
*/
- virtual void coarseRestr(DOFIndexed<double> *, RCNeighbourList*, int){};
+ virtual void coarseRestr(DOFIndexed<double> *, RCNeighbourList*, int)
+ {}
/** \brief
* Interpolates a DOFVector<WorldVector<double> > after refinement
*/
- virtual void refineInter(DOFVector<WorldVector<double> >*, RCNeighbourList*, int){};
+ virtual void refineInter(DOFVector<WorldVector<double> >*, RCNeighbourList*, int)
+ {}
/** \brief
* Interpolates a DOFVector<WorldVector<double> > after coarsening
*/
- virtual void coarseInter(DOFVector<WorldVector<double> >*, RCNeighbourList*, int){};
+ virtual void coarseInter(DOFVector<WorldVector<double> >*, RCNeighbourList*, int)
+ {}
/** \brief
* Restricts a DOFVector<WorldVector<double> > after coarsening
*/
- virtual void coarseRestr(DOFVector<WorldVector<double> >*, RCNeighbourList*, int){};
+ virtual void coarseRestr(DOFVector<WorldVector<double> >*, RCNeighbourList*, int)
+ {}
/** \brief
* Returns local dof indices of the element for the given fe space.
@@ -362,7 +370,7 @@ namespace AMDiS {
DegreeOfFreedom*) const
{
return NULL;
- };
+ }
/** \brief
* Returns local dof indices of the element for the given fe space.
@@ -370,7 +378,7 @@ namespace AMDiS {
virtual void getLocalIndicesVec(const Element*,
const DOFAdmin*,
Vector<DegreeOfFreedom>*) const
- {};
+ {}
/** \brief
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index ec34e98bff01e3817f4141c0812846b71fdc0733..ae0d0ec85f16ad5b7545993d34097a123d5b5371 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -250,14 +250,6 @@ namespace AMDiS {
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
if (condition && condition->isDirichlet()) {
- /* MatrixRow *matrixRow = &(matrix[row]);
- if (coupleMatrix) {
- matrixRow->resize(0);
- } else {
- matrixRow->resize(1);
- ((*matrixRow)[0]).col = row;
- ((*matrixRow)[0]).entry = 1.0;
- } */
applyDBCs.insert(static_cast<int>(row));
} else {
for (int j = 0; j < nCol; j++) { // for all columns
@@ -312,7 +304,7 @@ namespace AMDiS {
("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1);
// first entry is diagonal entry
- if ((rowSize == 0) && (rowFESpace == colFESpace)) {
+ if ((rowSize == 0) && (colFESpace == rowFESpace)) {
MatEntry newEntry = {irow, 0.0};
row->push_back(newEntry);
rowSize = 1;
@@ -458,11 +450,13 @@ namespace AMDiS {
*factorIt ? **factorIt : 1.0);
}
}
-
+
addElementMatrix(factor, *elementMatrix, bound);
}
- void DOFMatrix::assemble(double factor, ElInfo *rowElInfo, ElInfo *colElInfo,
+ void DOFMatrix::assemble(double factor,
+ ElInfo *rowElInfo, ElInfo *colElInfo,
+ ElInfo *smallElInfo, ElInfo *largeElInfo,
const BoundaryType *bound, Operator *op)
{
FUNCNAME("DOFMatrix::assemble()");
@@ -476,19 +470,30 @@ namespace AMDiS {
initElementMatrix(elementMatrix, rowElInfo, colElInfo);
if (op) {
- op->getElementMatrix(rowElInfo, colElInfo, elementMatrix);
+ op->getElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo,
+ elementMatrix);
} else {
std::vector<Operator*>::iterator it;
std::vector<double*>::iterator factorIt;
for (it = operators.begin(), factorIt = operatorFactor.begin();
it != operators.end();
++it, ++factorIt) {
- (*it)->getElementMatrix(rowElInfo, colElInfo,
- elementMatrix,
- *factorIt ? **factorIt : 1.0);
+ // If both levels are equal, we may use the standard assembler for
+ // one element.
+ if (rowElInfo->getLevel() == colElInfo->getLevel()) {
+ (*it)->getElementMatrix(rowElInfo,
+ elementMatrix,
+ *factorIt ? **factorIt : 1.0);
+ } else {
+ (*it)->getElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo,
+ elementMatrix,
+ *factorIt ? **factorIt : 1.0);
+ }
}
}
-
+
addElementMatrix(factor, *elementMatrix, bound);
}
@@ -671,6 +676,23 @@ namespace AMDiS {
}
}
+ int DOFMatrix::getNumCols() const {
+ int max = 0;
+ int rowSize = static_cast<int>(matrix.size());
+ for (int i = 0; i < rowSize; i++) {
+ int colSize = static_cast<int>(matrix[i].size());
+ for (int j = 0; j < colSize; j++) {
+ if (matrix[i][j].col > max) {
+ max = matrix[i][j].col;
+ }
+ }
+ }
+
+ // Add one to the maximum value, because the indeces start with 0.
+ return max + 1;
+ }
+
+
void DOFMatrix::createPictureFile(const char* filename, int dim)
{
png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
diff --git a/AMDiS/src/DOFMatrix.h b/AMDiS/src/DOFMatrix.h
index e95ab45373299fc6f2a94848e01917928fb24532..0285814aba27f92d34f35f77849af442594cca79 100644
--- a/AMDiS/src/DOFMatrix.h
+++ b/AMDiS/src/DOFMatrix.h
@@ -88,11 +88,11 @@ namespace AMDiS {
public:
MatEntryValueLess(const double& value)
: value_(value)
- {};
+ {}
bool operator()(const MatEntry& m) {
return (fabs(m.entry) < value_);
- };
+ }
};
@@ -106,11 +106,11 @@ namespace AMDiS {
public:
MatEntryColLess(const int& col)
: col_(col)
- {};
+ {}
bool operator()(const MatEntry& m) {
return (fabs(m.col) < col_);
- };
+ }
};
@@ -121,7 +121,7 @@ namespace AMDiS {
struct CmpMatEntryCol : public std::binary_function<MatEntry, MatEntry, bool> {
bool operator()(const MatEntry& m1, const MatEntry m2) const {
return m1.col < m2.col;
- };
+ }
};
@@ -132,7 +132,7 @@ namespace AMDiS {
struct CmpMatEntryValue : public std::binary_function<MatEntry, MatEntry, bool> {
bool operator()(const MatEntry& m1, const MatEntry m2) const {
return m1.entry < m2.entry;
- };
+ }
};
@@ -143,7 +143,7 @@ namespace AMDiS {
struct CmpMatEntryAbsValueLess : public std::binary_function<MatEntry, MatEntry, bool> {
bool operator()(const MatEntry& m1, const MatEntry m2) const {
return fabs(m1.entry) < fabs(m2.entry);
- };
+ }
};
/** \brief
@@ -153,7 +153,7 @@ namespace AMDiS {
struct CmpMatEntryAbsValueGreater : public std::binary_function<MatEntry, MatEntry, bool> {
bool operator()(const MatEntry& m1, const MatEntry m2) const {
return fabs(m1.entry) > fabs(m2.entry);
- };
+ }
};
// ============================================================================
@@ -183,7 +183,7 @@ namespace AMDiS {
Iterator(DOFIndexed< std::vector<MatEntry> > *c,
DOFIteratorType type)
: DOFIterator< std::vector<MatEntry> >(c, type)
- {};
+ {}
};
/** \brief
@@ -237,14 +237,14 @@ namespace AMDiS {
*/
std::vector< std::vector<MatEntry> >::iterator begin() {
return matrix.begin();
- };
+ }
/** \brief
* Returns an iterator to the end of matrix rows (\ref matrix.end())
*/
std::vector< std::vector<MatEntry> >::iterator end() {
return matrix.end();
- };
+ }
/** \brief
* Used by DOFAdmin to compress the DOFMatrix
@@ -262,28 +262,28 @@ namespace AMDiS {
*/
inline bool entryUsed(DegreeOfFreedom a,int b) const {
return (((matrix[a])[b]).col >= 0);
- };
+ }
/** \brief
* Returns true if matrix[a][b] has entry \ref NO_MORE_ENTRIES
*/
inline bool noMoreEntries(DegreeOfFreedom a,int b) const {
return (NO_MORE_ENTRIES == ((matrix[a])[b]).col);
- };
+ }
/** \brief
* Returns \ref coupleMatrix.
*/
inline bool isCoupleMatrix() {
return coupleMatrix;
- };
+ }
/** \brief
* Returns \ref coupleMatrix.
*/
inline void setCoupleMatrix(bool c) {
coupleMatrix = c;
- };
+ }
/** \brief
* Matrix-matrix multiplication.
@@ -310,31 +310,31 @@ namespace AMDiS {
operators.push_back(op);
operatorFactor.push_back(factor);
operatorEstFactor.push_back(estFactor);
- };
+ }
inline std::vector<double*>::iterator getOperatorFactorBegin() {
return operatorFactor.begin();
- };
+ }
inline std::vector<double*>::iterator getOperatorFactorEnd() {
return operatorFactor.end();
- };
+ }
inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
return operatorEstFactor.begin();
- };
+ }
inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
return operatorEstFactor.end();
- };
+ }
inline std::vector<Operator*>::iterator getOperatorsBegin() {
return operators.begin();
- };
+ }
inline std::vector<Operator*>::iterator getOperatorsEnd() {
return operators.end();
- };
+ }
Flag getAssembleFlag();
@@ -365,7 +365,9 @@ namespace AMDiS {
const BoundaryType *bound, Operator *op = NULL);
- void assemble(double factor, ElInfo *rowElInfo, ElInfo *colElInfo,
+ void assemble(double factor,
+ ElInfo *rowElInfo, ElInfo *colElInfo,
+ ElInfo *smallElInfo, ElInfo *largeElInfo,
const BoundaryType *bound, Operator *op = NULL);
@@ -382,25 +384,25 @@ namespace AMDiS {
*/
std::vector< std::vector<MatEntry> >& getMatrix() {
return matrix;
- };
+ }
void setMatrix(std::vector< std::vector<MatEntry> > m) {
matrix = m;
- };
+ }
/** \brief
* Returns \ref matrix[n]
*/
const std::vector<MatEntry>& getRow(int n) const {
return matrix[n];
- };
+ }
/** \brief
* Returns \ref matrix[n]
*/
std::vector<MatEntry>& getRow(int n) {
return matrix[n];
- };
+ }
/** \brief
* Returns whether restriction should be performed after coarsening
@@ -408,35 +410,35 @@ namespace AMDiS {
*/
virtual bool coarseRestrict() {
return false;
- };
+ }
/** \brief
* Returns const \ref rowFESpace
*/
const FiniteElemSpace* getRowFESpace() const {
return rowFESpace;
- };
+ }
/** \brief
* Returns const \ref colFESpace
*/
const FiniteElemSpace* getColFESpace() const {
return colFESpace;
- };
+ }
/** \brief
* Returns const \ref rowFESpace
*/
const FiniteElemSpace* getFESpace() const {
return rowFESpace;
- };
+ }
/** \brief
* Returns number of rows (\ref matrix.size())
*/
inline int getSize() const {
return matrix.size();
- };
+ }
/** \brief
* Returns the number of used rows (equal to number of used DOFs in
@@ -444,34 +446,20 @@ namespace AMDiS {
*/
inline int getUsedSize() const {
return rowFESpace->getAdmin()->getUsedSize();
- };
+ }
/** \brief
* Returns number of cols. For that, the function iteratos over all
* rows and searchs for the entry with the highest col number.
*/
- int getNumCols() const {
- int max = 0;
- int rowSize = static_cast<int>(matrix.size());
- for (int i = 0; i < rowSize; i++) {
- int colSize = static_cast<int>(matrix[i].size());
- for (int j = 0; j < colSize; j++) {
- if (matrix[i][j].col > max) {
- max = matrix[i][j].col;
- }
- }
- }
-
- // Add one to the maximum value, because the indeces start with 0.
- return max + 1;
- };
+ int getNumCols() const;
/** \brief
* Returns \ref name
*/
inline const std::string& getName() const {
return name;
- };
+ }
/** \brief
* Resizes \ref matrix to n rows
@@ -479,7 +467,7 @@ namespace AMDiS {
inline void resize(int n) {
TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
matrix.resize(n);
- };
+ }
/** \brief
* Returns \ref matrix[i] which is the i-th row
@@ -488,7 +476,7 @@ namespace AMDiS {
TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
("Illegal matrix index %d.\n",i);
return matrix[i];
- };
+ }
/** \brief
* Returns \ref matrix[i] which is the i-th row
@@ -497,14 +485,14 @@ namespace AMDiS {
TEST_EXIT_DBG((i >= 0) && (i < (static_cast<int>(matrix.size()))))
("Illegal vector index %d.\n", i);
return matrix[i];
- };
+ }
/** \brief
* Access to \ref matrix[a][b].entry
*/
inline double physAcc(DegreeOfFreedom a, DegreeOfFreedom b) const {
return matrix[a][b].entry;
- };
+ }
/** \brief
* Access to \ref matrix[a][b].col
@@ -513,7 +501,7 @@ namespace AMDiS {
DegreeOfFreedom b) const
{
return matrix[a][b].col;
- };
+ }
/** \brief
* Returns physical column index of logical index b in row a
@@ -526,7 +514,7 @@ namespace AMDiS {
return j;
return -1;
- };
+ }
/** \brief
* Returns value at logical indices a,b
@@ -546,7 +534,7 @@ namespace AMDiS {
DegreeOfFreedom c)
{
matrix[a][b].col = c;
- };
+ }
/** \brief
* Creates an entry with logical indices irow, icol if there is no entry
@@ -566,7 +554,7 @@ namespace AMDiS {
return &(it->entry);
};
return NULL;
- };
+ }
void addMatEntry(int row, MatEntry entry);
@@ -600,19 +588,19 @@ namespace AMDiS {
inline std::vector<Operator*> getOperators() {
return operators;
- };
+ }
inline std::vector<double*> getOperatorFactor() {
return operatorFactor;
- };
+ }
inline std::vector<double*> getOperatorEstFactor() {
return operatorEstFactor;
- };
+ }
inline BoundaryManager* getBoundaryManager() const {
return boundaryManager;
- };
+ }
std::set<int>* getApplyDBCs() {
return &applyDBCs;
@@ -620,7 +608,7 @@ namespace AMDiS {
inline void setBoundaryManager(BoundaryManager *bm) {
boundaryManager = bm;
- };
+ }
void createPictureFile(const char* filename, int dim);
@@ -637,7 +625,7 @@ namespace AMDiS {
out.write(reinterpret_cast<const char*>(&vecSize), sizeof(unsigned int));
out.write(reinterpret_cast<const char*>(&(matrix[i][0])), vecSize * sizeof(MatEntry));
}
- };
+ }
void deserialize(std::istream &in) {
unsigned int matrixSize, vecSize;
@@ -649,7 +637,7 @@ namespace AMDiS {
matrix[i].resize(vecSize);
in.read(reinterpret_cast<char*>(&(matrix[i][0])), vecSize * sizeof(MatEntry));
}
- };
+ }
int memsize();
diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index 214df2cbe9f66a4067238f123ab108054ec384fe..b890b8d0c19bc2ab6438e0c5f12bbb2e36612f94 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -132,9 +132,9 @@ namespace AMDiS {
*elInfo2;
// update rest
- double newRest = rest - 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
+ rest -= 1.0 / (1 << ((*elInfoSmall)->getLevel() - (*elInfoLarge)->getLevel()));
- if (newRest < 1e-32) {
+ if (rest < 1e-32) {
// large element finished -> increment both elements
rest = 1.0;
inc1 = true;
@@ -145,19 +145,6 @@ namespace AMDiS {
// increment only small element
inc1 = (*elInfo1 == *elInfoSmall) ? true : false;
inc2 = (*elInfo2 == *elInfoSmall) ? true : false;
-
- // If rest is 1, than this is the first element pair with this large
- // element. Therefore, fill the displacement stack of the small element
- // traverse stack.
- if (rest == 1.0) {
- if (*elInfo1 == *elInfoSmall) {
- stack1.startDisplacementCalculation((*elInfo2)->getLevel());
- } else {
- stack2.startDisplacementCalculation((*elInfo1)->getLevel());
- }
- }
-
- rest = newRest;
}
}
@@ -169,30 +156,17 @@ namespace AMDiS {
if (!fillSubElemMat)
return;
- VectorOfFixVecs<DimVec<double> > *subCoords =
- NEW VectorOfFixVecs<DimVec<double> >(elInfo1->getMesh()->getDim(),
- elInfo1->getMesh()->getDim() + 1,
- NO_INIT);
-
DimMat<double> *subCoordsMat = elInfoSmall->getSubElemCoordsMat();
if (!subCoordsMat) {
subCoordsMat = NEW DimMat<double>(elInfo1->getMesh()->getDim());
}
if (elInfo1 == elInfoSmall) {
- stack1.getCoordsInElem(elInfo2, subCoords);
+ stack1.getCoordsInElem(elInfo2, subCoordsMat);
} else {
- stack2.getCoordsInElem(elInfo1, subCoords);
+ stack2.getCoordsInElem(elInfo1, subCoordsMat);
}
- for (int i = 0; i < elInfo1->getMesh()->getDim() + 1; i++) {
- for (int j = 0; j < elInfo1->getMesh()->getDim() + 1; j++) {
- (*subCoordsMat)[j][i] = (*subCoords)[i][j];
- }
- }
-
elInfoSmall->setSubElemCoordsMat(subCoordsMat);
-
- DELETE subCoords;
}
}
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 3a9dc0ed3b416da979e172be5c7eb1dd626c9565..3951b8296f15bec2385e599302ec5c64d7161f66 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -477,9 +477,9 @@ namespace AMDiS {
return(0.0);
}
- virtual void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const = 0;
+ virtual void getRefSimplexCoords(DimMat<double> *coords) const = 0;
- virtual void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
+ virtual void getSubElementCoords(DimMat<double> *coords,
int iChild) const = 0;
protected:
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 9b6fbeb6635751cd87b78f7918c1106eccbfc87c..9ae0c0a03074771091a5111bdf3cf955f25ec7b7 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -272,24 +272,24 @@ namespace AMDiS {
return;
}
- void ElInfo1d::getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const
+ void ElInfo1d::getRefSimplexCoords(DimMat<double> *coords) const
{
(*coords)[0][0] = 1.0;
- (*coords)[0][1] = 0.0;
-
(*coords)[1][0] = 0.0;
+
+ (*coords)[0][1] = 0.0;
(*coords)[1][1] = 1.0;
}
- void ElInfo1d::getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
+ void ElInfo1d::getSubElementCoords(DimMat<double> *coords,
int iChild) const
{
if (iChild == 0) {
- (*coords)[1][0] = ((*coords)[0][0] + (*coords)[1][0]) * 0.5;
- (*coords)[1][1] = ((*coords)[0][1] + (*coords)[1][1]) * 0.5;
+ (*coords)[0][1] = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
+ (*coords)[1][1] = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
} else {
- (*coords)[0][0] = ((*coords)[0][0] + (*coords)[1][0]) * 0.5;
- (*coords)[0][1] = ((*coords)[0][1] + (*coords)[1][1]) * 0.5;
+ (*coords)[0][0] = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
+ (*coords)[1][0] = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
}
}
diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h
index 728118ff7474724a9acb8d3e467b1a70322609ce..619bf3e698a2b9ce09566dfb5c0e360c8faddd9d 100644
--- a/AMDiS/src/ElInfo1d.h
+++ b/AMDiS/src/ElInfo1d.h
@@ -81,9 +81,9 @@ namespace AMDiS {
return (i + 1) % 2;
}
- void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const;
+ void getRefSimplexCoords(DimMat<double> *coords) const;
- void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
+ void getSubElementCoords(DimMat<double> *coords,
int iChild) const;
};
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index f3f33d18cd344159077ba2fb0cf932241008a803..af4bebfaf4e0cfeacab270186377a82b134dde17 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -591,4 +591,50 @@ namespace AMDiS {
return(det);
}
+
+ void ElInfo2d::getRefSimplexCoords(DimMat<double> *coords) const
+ {
+ (*coords)[0][0] = 1.0;
+ (*coords)[1][0] = 0.0;
+ (*coords)[2][0] = 0.0;
+
+ (*coords)[0][1] = 0.0;
+ (*coords)[1][1] = 1.0;
+ (*coords)[2][1] = 0.0;
+
+ (*coords)[0][2] = 0.0;
+ (*coords)[1][2] = 0.0;
+ (*coords)[2][2] = 1.0;
+ }
+
+ void ElInfo2d::getSubElementCoords(DimMat<double> *coords,
+ int iChild) const
+ {
+ double c0 = ((*coords)[0][0] + (*coords)[0][1]) * 0.5;
+ double c1 = ((*coords)[1][0] + (*coords)[1][1]) * 0.5;
+ double c2 = ((*coords)[2][0] + (*coords)[2][1]) * 0.5;
+
+ if (iChild == 0) {
+ (*coords)[0][1] = (*coords)[0][0];
+ (*coords)[1][1] = (*coords)[1][0];
+ (*coords)[2][1] = (*coords)[2][0];
+
+ (*coords)[0][0] = (*coords)[0][2];
+ (*coords)[1][0] = (*coords)[1][2];
+ (*coords)[2][0] = (*coords)[2][2];
+ } else {
+ (*coords)[0][1] = (*coords)[0][2];
+ (*coords)[1][1] = (*coords)[1][2];
+ (*coords)[2][1] = (*coords)[2][2];
+
+ (*coords)[0][0] = (*coords)[0][1];
+ (*coords)[1][0] = (*coords)[1][1];
+ (*coords)[2][0] = (*coords)[2][1];
+ }
+
+ (*coords)[0][2] = c0;
+ (*coords)[1][2] = c1;
+ (*coords)[2][2] = c2;
+ }
+
}
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index 48ca680b7ce1deb168e79807f7d649995efc4c9b..056111ff740e5de193e377e536368b02b68f5a58 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/AMDiS/src/ElInfo2d.h
@@ -75,22 +75,12 @@ namespace AMDiS {
/** \brief
* 2-dimensional realisation of ElInfo's getElementNormal method.
*/
- double getElementNormal( WorldVector<double> &normal) const;
+ double getElementNormal(WorldVector<double> &normal) const;
- void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const
- {
- FUNCNAME("ElInfo2d::getRefSimplexCoords()");
+ void getRefSimplexCoords(DimMat<double> *coords) const;
- ERROR_EXIT("NOT YET\n");
- }
-
- void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
- int iChild) const
- {
- FUNCNAME("ElInfo2d::getSubElementCoords()");
-
- ERROR_EXIT("NOT YET\n");
- }
+ void getSubElementCoords(DimMat<double> *coords,
+ int iChild) const;
protected:
/** \brief
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index 1be858be72fd5f04fdac98f8751f00ff1a9271ad..af403f2a5a27fa935bcc5c406aa57437f8d7c021 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/AMDiS/src/ElInfo3d.h
@@ -49,7 +49,7 @@ namespace AMDiS {
: ElInfo(aMesh)
{
tmpWorldVecs.resize(4);
- };
+ }
/** \brief
* Assignment operator
@@ -96,37 +96,37 @@ namespace AMDiS {
*/
inline unsigned char getType() const {
return el_type;
- };
+ }
/** \brief
* get ElInfo's \ref orientation
*/
inline signed char getOrientation() const {
return orientation;
- };
+ }
/** \brief
* set ElInfo's \ref el_type to t
*/
inline void setType(unsigned char t) {
el_type = t;
- };
+ }
/** \brief
* set ElInfo's \ref orientation to o
*/
inline void setOrientation(signed char o) {
orientation = o;
- };
+ }
- void getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const
+ void getRefSimplexCoords(DimMat<double> *coords) const
{
FUNCNAME("ElInfo3d::getRefSimplexCoords()");
ERROR_EXIT("NOT YET\n");
}
- void getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
+ void getSubElementCoords(DimMat<double> *coords,
int iChild) const
{
FUNCNAME("ElInfo3d::getSubElementCoords()");
diff --git a/AMDiS/src/FirstOrderAssembler.h b/AMDiS/src/FirstOrderAssembler.h
index af3b7f2963b09a805c6a18d0cb880385a0a9e919..ee10b3c8faa2ad9d83014fe237c6391b8bbd795a 100644
--- a/AMDiS/src/FirstOrderAssembler.h
+++ b/AMDiS/src/FirstOrderAssembler.h
@@ -98,6 +98,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -136,6 +138,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -173,6 +177,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -208,6 +214,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -245,6 +253,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -293,6 +303,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index a86aa0da0c239184f0459f7fe67d66d2e4563e4b..6739111080c33f262e9a7e8eed6d4a7e83481571 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -16,12 +16,12 @@
namespace AMDiS {
- std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- int Lagrange::nBasFctsDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- std::vector<BasFctType*> Lagrange::phiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- std::vector<GrdBasFctType*> Lagrange::grdPhiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- std::vector<D2BasFctType*> Lagrange::D2PhiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
+ std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ int Lagrange::nBasFctsDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ std::vector<BasFctType*> Lagrange::phiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ std::vector<GrdBasFctType*> Lagrange::grdPhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ std::vector<D2BasFctType*> Lagrange::D2PhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
std::list<Lagrange*> Lagrange::allBasFcts;
@@ -1615,50 +1615,47 @@ namespace AMDiS {
cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);
- if (typ == 0)
- {
- switch(lr_set)
- {
- case 1:
- cdi = el->getChild(1)->getDOF(node0+1, n0);
- TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n");
- (*drv)[cdi] =
- (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
- + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
- + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
- break;
- case 2:
- cdi = el->getChild(1)->getDOF(node0+2, n0);
- TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n");
- (*drv)[cdi] =
- (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
- + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
- + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
- break;
- }
- }
- else
- {
- switch(lr_set)
- {
- case 1:
- cdi = el->getChild(1)->getDOF(node0+2, n0);
- TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n");
- (*drv)[cdi] =
- (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
- + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
- + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
- break;
- case 2:
- cdi = el->getChild(1)->getDOF(node0+1, n0);
- TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n");
- (*drv)[cdi] =
- (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
- + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
- + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
- break;
- }
- }
+ if (typ == 0) {
+ switch(lr_set)
+ {
+ case 1:
+ cdi = el->getChild(1)->getDOF(node0+1, n0);
+ TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n");
+ (*drv)[cdi] =
+ (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
+ + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
+ break;
+ case 2:
+ cdi = el->getChild(1)->getDOF(node0+2, n0);
+ TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n");
+ (*drv)[cdi] =
+ (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
+ + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
+ break;
+ }
+ } else {
+ switch(lr_set)
+ {
+ case 1:
+ cdi = el->getChild(1)->getDOF(node0+2, n0);
+ TEST_EXIT_DBG(cdi == cd[18])("cdi != cd[18]\n");
+ (*drv)[cdi] =
+ (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
+ + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
+ break;
+ case 2:
+ cdi = el->getChild(1)->getDOF(node0+1, n0);
+ TEST_EXIT_DBG(cdi == cd[17])("cdi != cd[17]\n");
+ (*drv)[cdi] =
+ (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]])
+ + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
+ + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
+ break;
+ }
+ }
}
return;
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index 36a0165ee3a7742eac535240c09a8eda344d86f4..a977e2fac6f394aa75aa2549631cdb56ae94d7f3 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -67,7 +67,7 @@ namespace AMDiS {
* degree. Multiple instantiation of identical basis functions is avoided
* by rembering once created basis functions in \ref allBasFcts.
*/
- static Lagrange* getLagrange(int dim_, int degree_);
+ static Lagrange* getLagrange(int dim, int degree);
/** \brief
* Implements BasisFunction::interpol
@@ -94,7 +94,8 @@ namespace AMDiS {
/** \brief
* Implements BasisFunction::getDOFIndices
*/
- const DegreeOfFreedom* getDOFIndices(const Element*, const DOFAdmin&,
+ const DegreeOfFreedom* getDOFIndices(const Element*,
+ const DOFAdmin&,
DegreeOfFreedom*) const;
/** \brief
@@ -124,7 +125,7 @@ namespace AMDiS {
{
if (refineInter_fct)
(*refineInter_fct)(drv, list, n, this);
- };
+ }
/** \brief
* Implements BasisFunction::coarseRestrict
@@ -135,7 +136,7 @@ namespace AMDiS {
{
if (coarseRestr_fct)
(*coarseRestr_fct)(drv, list, n, this);
- };
+ }
/** \brief
* Implements BasisFunction::coarseInter
@@ -185,7 +186,7 @@ namespace AMDiS {
* Recursive calculation of coordinates. Used by \ref setBary
*/
void createCoords(int* coordInd, int numCoords, int dimIndex, int rest,
- DimVec<double>* vec=NULL);
+ DimVec<double>* vec = NULL);
/** \brief
* Used by \ref setBary
@@ -205,7 +206,8 @@ namespace AMDiS {
/** \brief
* Used by \ref getDOFIndices and \ref getVec
*/
- int* orderOfPositionIndices(const Element* el, GeoIndex position,
+ int* orderOfPositionIndices(const Element* el,
+ GeoIndex position,
int positionIndex) const;
/** \brief
@@ -223,12 +225,12 @@ namespace AMDiS {
/** \name static dim-degree-arrays
* \{
*/
- static std::vector<DimVec<double>* > baryDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- static DimVec<int>* ndofDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- static int nBasFctsDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- static std::vector<BasFctType*> phiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- static std::vector<GrdBasFctType*> grdPhiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
- static std::vector<D2BasFctType*> D2PhiDimDegree[MAX_DIM+1][MAX_DEGREE+1];
+ static std::vector<DimVec<double>* > baryDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ static DimVec<int>* ndofDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ static int nBasFctsDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ static std::vector<BasFctType*> phiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ static std::vector<GrdBasFctType*> grdPhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
+ static std::vector<D2BasFctType*> D2PhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
/** \} */
/** \brief
@@ -378,83 +380,83 @@ namespace AMDiS {
*/
inline double operator()(const DimVec<double>& lambda) const {
return func(lambda, vertices);
- };
+ }
/** \name basis functions for different degrees
* \{
*/
- // ====== Lagrange0 ================================================
+ // ====== Lagrange, degree = 0 =====================================
// center
inline static double phi0c(const DimVec<double>&, int*) {
return 1.0;
- };
+ }
- // ====== Lagrange1 ================================================
+ // ====== Lagrange, degree = 1 =====================================
// vertex
inline static double phi1v(const DimVec<double>& lambda, int* vertices) {
return lambda[vertices[0]];
- };
+ }
- // ====== Lagrange2 ================================================
+ // ====== Lagrange, degree = 2 =====================================
// vertex
inline static double phi2v(const DimVec<double>& lambda, int* vertices) {
return lambda[vertices[0]] * (2.0 * lambda[vertices[0]] - 1.0);
- };
+ }
// edge
inline static double phi2e(const DimVec<double>& lambda, int* vertices) {
return (4.0 * lambda[vertices[0]] * lambda[vertices[1]]);
- };
+ }
- // ====== Lagrange3 ================================================
+ // ====== Lagrange, degree = 3 =====================================
// vertex
inline static double phi3v(const DimVec<double>& lambda, int* vertices) {
return (4.5 * (lambda[vertices[0]] - 1.0) * lambda[vertices[0]] + 1.0) *
lambda[vertices[0]];
- };
+ }
// edge
inline static double phi3e(const DimVec<double>& lambda, int* vertices) {
return (13.5 * lambda[vertices[0]] - 4.5) *
lambda[vertices[0]] * lambda[vertices[1]];
- };
+ }
// face
inline static double phi3f(const DimVec<double>& lambda, int* vertices) {
return 27.0 * lambda[vertices[0]] * lambda[vertices[1]] *
lambda[vertices[2]];
- };
+ }
- // ====== Lagrange4 ================================================
+ // ====== Lagrange, degree = 4 ======================================
// vertex
inline static double phi4v(const DimVec<double>& lambda, int* vertices) {
return (((32.0 * lambda[vertices[0]] - 48.0) * lambda[vertices[0]] + 22.0)
* lambda[vertices[0]] - 3.0) * lambda[vertices[0]] / 3.0;
- };
+ }
// edge
inline static double phi4e0(const DimVec<double>& lambda, int* vertices) {
return ((128.0 * lambda[vertices[0]] - 96.0) * lambda[vertices[0]] + 16.0)
* lambda[vertices[0]] * lambda[vertices[1]] / 3.0;
- };
+ }
inline static double phi4e1(const DimVec<double>& lambda, int* vertices) {
return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] *
(4.0 * lambda[vertices[1]] - 1.0) * lambda[vertices[1]] * 4.0;
- };
+ }
// face
inline static double phi4f(const DimVec<double>& lambda, int* vertices) {
return (4.0 * lambda[vertices[0]] - 1.0) * lambda[vertices[0]] *
lambda[vertices[1]] * lambda[vertices[2]] * 32.0;
- };
+ }
// center
inline static double phi4c(const DimVec<double>& lambda, int* vertices) {
return 256.0 * lambda[vertices[0]] * lambda[vertices[1]] *
lambda[vertices[2]] * lambda[vertices[3]];
- };
+ }
};
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 676451a4bc21c5b03a7936f29f59692a51ebb8d2..f4ab8bcb4b3990e7d24f67bbe636666dc25b7359 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -220,6 +220,8 @@ namespace AMDiS {
void Operator::getElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *userMat,
double factor)
{
@@ -230,6 +232,7 @@ namespace AMDiS {
}
assembler[myRank]->calculateElementMatrix(rowElInfo, colElInfo,
+ smallElInfo, largeElInfo,
userMat, factor);
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 73540d4cbadee09efea8a719d2ab533443b03ea2..11210f56a5b3ee0651ac28cf7c7fccdc88cb6a1b 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -70,12 +70,13 @@ namespace AMDiS {
OperatorTerm(int degree_)
: properties(0),
degree(degree_)
- {};
+ {}
/** \brief
* Destructor.
*/
- virtual ~OperatorTerm() {};
+ virtual ~OperatorTerm()
+ {}
/** \brief
* Virtual method. It's called by SubAssembler::initElement() for
@@ -85,7 +86,7 @@ namespace AMDiS {
virtual void initElement(const ElInfo*,
SubAssembler*,
Quadrature *quad = NULL)
- {};
+ {}
/** \brief
* Specifies whether the matrix of the term is symmetric
@@ -97,7 +98,7 @@ namespace AMDiS {
*/
inline bool isPWConst() {
return (degree == 0);
- };
+ }
/** \brief
* Returns true, if the term has a symmetric matrix,
@@ -110,7 +111,7 @@ namespace AMDiS {
*/
inline int getDegree() {
return degree;
- };
+ }
/** \brief
* Evaluation of the OperatorTerm at all quadrature points.
@@ -178,7 +179,7 @@ namespace AMDiS {
val *= factor;
Lb[i] += val;
}
- };
+ }
protected:
/** \brief
@@ -230,12 +231,15 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
- SecondOrderTerm(int deg) : OperatorTerm(deg) {};
+ SecondOrderTerm(int deg)
+ : OperatorTerm(deg)
+ {}
/** \brief
* Destructor.
*/
- virtual ~SecondOrderTerm() {};
+ virtual ~SecondOrderTerm()
+ {}
/** \brief
* Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
@@ -271,7 +275,7 @@ namespace AMDiS {
: SecondOrderTerm(0)
{
setSymmetric(true);
- };
+ }
/** \brief
* Implenetation of SecondOrderTerm::getLALt().
@@ -283,7 +287,7 @@ namespace AMDiS {
for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), 1.0);
}
- };
+ }
/** \brief
* Implementation of SecondOrderTerm::eval().
@@ -306,7 +310,7 @@ namespace AMDiS {
result[iq] += factor * resultQP;
}
}
- };
+ }
/** \brief
* Implenetation of SecondOrderTerm::weakEval().
@@ -320,8 +324,8 @@ namespace AMDiS {
result[iq] += grdUhAtQP[iq];
}
}
- };
- };
+ }
+ };
// ============================================================================
@@ -345,7 +349,7 @@ namespace AMDiS {
*factor = f;
setSymmetric(true);
- };
+ }
/** \brief
* Constructor.
@@ -355,7 +359,7 @@ namespace AMDiS {
factor(fptr)
{
setSymmetric(true);
- };
+ }
/** \brief
* Implements SecondOrderTerm::getLALt().
@@ -365,7 +369,7 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < numPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*factor));
- };
+ }
/** \brief
@@ -389,7 +393,7 @@ namespace AMDiS {
result[iq] += resultQP * f * (*factor);
}
}
- };
+ }
/** \brief
* Implenetation of SecondOrderTerm::weakEval().
@@ -403,7 +407,7 @@ namespace AMDiS {
axpy(*factor, grdUhAtQP[iq], result[iq]);
}
}
- };
+ }
private:
/** \brief
@@ -439,7 +443,7 @@ namespace AMDiS {
symmetric(sym)
{
setSymmetric(symmetric);
- };
+ }
/** \brief
* Implementation of \ref OperatorTerm::initElement().
@@ -513,7 +517,7 @@ namespace AMDiS {
{
symmetric = matrix.isSymmetric();
setSymmetric(symmetric);
- };
+ }
/** \brief
* Implements SecondOrderTerm::getLALt().
@@ -524,7 +528,7 @@ namespace AMDiS {
int iq;
for(iq = 0; iq < numPoints; iq++)
lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
- };
+ }
/** \brief
* Implenetation of SecondOrderTerm::eval().
@@ -580,7 +584,7 @@ namespace AMDiS {
*factor = f;
setSymmetric(xi == xj);
- };
+ }
/** \brief
* Constructor.
@@ -592,7 +596,7 @@ namespace AMDiS {
factor(fptr)
{
setSymmetric(xi == xj);
- };
+ }
/** \brief
* Implements SecondOrderTerm::getLALt().
@@ -603,7 +607,7 @@ namespace AMDiS {
int iq;
for(iq = 0; iq < numPoints; iq++)
lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor));
- };
+ }
/** \brief
* Implenetation of SecondOrderTerm::eval().
@@ -621,7 +625,7 @@ namespace AMDiS {
result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac;
}
}
- };
+ }
/** \brief
* Implenetation of SecondOrderTerm::weakEval().
@@ -636,7 +640,7 @@ namespace AMDiS {
result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
}
}
- };
+ }
private:
/** \brief
@@ -3781,6 +3785,8 @@ namespace AMDiS {
virtual void getElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *userMat,
double factor = 1.0);
@@ -3799,14 +3805,14 @@ namespace AMDiS {
*/
inline const FiniteElemSpace *getRowFESpace() {
return rowFESpace;
- };
+ }
/** \brief
* Returns \ref colFESpace
*/
inline const FiniteElemSpace *getColFESpace() {
return colFESpace;
- };
+ }
/** \brief
* Sets \ref uhOld.
@@ -3818,7 +3824,7 @@ namespace AMDiS {
*/
inline const DOFVectorBase<double> *getUhOld() {
return uhOld;
- };
+ }
/** \brief
* Returns \ref assembler
@@ -3835,28 +3841,28 @@ namespace AMDiS {
*/
inline const bool isMatrixOperator() {
return type.isSet(MATRIX_OPERATOR);
- };
+ }
/** \brief
* Returns whether this is a vector operator
*/
inline const bool isVectorOperator() {
return type.isSet(VECTOR_OPERATOR);
- };
+ }
/** \brief
* Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
*/
inline void setFillFlag(Flag f) {
fillFlag = f;
- };
+ }
/** \brief
* Returns \ref fillFlag
*/
inline Flag getFillFlag() {
return fillFlag;
- };
+ }
/** \brief
* Initialization of the needed SubAssemblers using the given quadratures.
diff --git a/AMDiS/src/Preconditioner.h b/AMDiS/src/Preconditioner.h
index 651fea5d70225408fa0fb78bfe5738e705d9a2f4..32fe39ee5308edf4fe2b4cece8db29a9a01690b8 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 num_threads(size)
+#pragma omp parallel for num_threads(min(size, omp_get_max_threads()))
#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 num_threads(size)
+#pragma omp parallel for num_threads(min(size, omp_get_max_threads()))
#endif
for (i = 0; i < size; i++) {
scalPrecons[i]->precon(x->getDOFVector(i));
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 65ddfab53658b68e6bda7793cc8ad1dd47ff53fa..bdfb390bcb61762e962d07c3ef8769ea18b9346f 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -1033,6 +1033,7 @@ namespace AMDiS {
ElInfo *rowElInfo, *colElInfo;
ElInfo *largeElInfo, *smallElInfo;
+ dualTraverse.setFillSubElemMat(true);
bool cont = dualTraverse.traverseFirst(rowMesh, colMesh, -1, -1,
assembleFlag, assembleFlag,
&rowElInfo, &colElInfo,
@@ -1040,33 +1041,18 @@ namespace AMDiS {
while (cont) {
Element *rowElem = rowElInfo->getElement();
Element *colElem = colElInfo->getElement();
-
- if (rowElInfo->getLevel() != colElInfo->getLevel()) {
- if (smallElInfo == rowElInfo) {
- std::cout << "Row = small\n";
- ERROR_EXIT("NOCH EIN PAAR GEDANKEN MACHEN!\n");
- }
-
- if (useGetBound_) {
- basisFcts->getBound(rowElInfo, bound);
- }
-
- if (matrix) {
- matrix->assemble(1.0, rowElInfo, colElInfo, bound);
- }
- } else {
- if (useGetBound_) {
- basisFcts->getBound(rowElInfo, bound);
- }
- if (matrix) {
- matrix->assemble(1.0, rowElInfo, bound);
-
- if (matrix->getBoundaryManager()) {
- matrix->getBoundaryManager()->
- fillBoundaryConditions(rowElInfo, matrix);
- }
- }
+ if (useGetBound_) {
+ basisFcts->getBound(rowElInfo, bound);
+ }
+
+ if (matrix) {
+ matrix->assemble(1.0, rowElInfo, colElInfo, smallElInfo, largeElInfo, bound);
+
+ if (matrix->getBoundaryManager()) {
+ matrix->getBoundaryManager()->
+ fillBoundaryConditions(rowElInfo, matrix);
+ }
}
if (vector) {
diff --git a/AMDiS/src/SecondOrderAssembler.h b/AMDiS/src/SecondOrderAssembler.h
index 732ec3cecf66e7f73dcf1379ad7ca008d9dfa11a..e31e2d161eaedb20fda5372325d9353ece931fda 100644
--- a/AMDiS/src/SecondOrderAssembler.h
+++ b/AMDiS/src/SecondOrderAssembler.h
@@ -83,6 +83,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -122,6 +124,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -170,6 +174,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
diff --git a/AMDiS/src/SubAssembler.h b/AMDiS/src/SubAssembler.h
index a74420dbb7467aaa8b4e7b8ed0b0c32fadb2ccf9..0870472ddfc03b8b87e314f884793915b6478e17 100644
--- a/AMDiS/src/SubAssembler.h
+++ b/AMDiS/src/SubAssembler.h
@@ -68,6 +68,8 @@ namespace AMDiS {
virtual void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat) = 0;
/** \brief
diff --git a/AMDiS/src/SystemVector.h b/AMDiS/src/SystemVector.h
index a9bc7811c2d8fc9ff083ee22c2c34d0c9b503eb9..9a1042ae39704ecc8342ac360f65b43c88143c96 100644
--- a/AMDiS/src/SystemVector.h
+++ b/AMDiS/src/SystemVector.h
@@ -511,7 +511,7 @@ namespace AMDiS {
TEST_EXIT_DBG(size == matrix.getNumCols())("incompatible sizes\n");
#ifdef _OPENMP
-#pragma omp parallel for schedule(static, 1) num_threads(size)
+#pragma omp parallel for schedule(static, 1) num_threads(min(size, omp_get_max_threads()))
#endif
for (i = 0; i < size; i++) {
if (!add)
@@ -541,7 +541,7 @@ namespace AMDiS {
int i;
#ifdef _OPENMP
-#pragma omp parallel for schedule(static, 1) num_threads(size)
+#pragma omp parallel for schedule(static, 1) num_threads(min(size, omp_get_max_threads()))
#endif
for (i = 0; i < size; i++) {
axpy(a, *(x.getDOFVector(i)), *(y.getDOFVector(i)));
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 20981eb35444cd852abd3ceaec5d022983c841b3..6f8383ae0d222048055d5fc344cf39147588d52b 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -1092,17 +1092,15 @@ namespace AMDiS {
}
void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo,
- VectorOfFixVecs<DimVec<double> > *coords)
+ DimMat<double> *coords)
{
- int upperLevel = static_cast<int>(upperElInfo->getLevel());
- int lowerLevel = static_cast<int>(elinfo_stack[stack_used]->getLevel());
- int stackPos = stack_used;
+ int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo->getLevel();
+
upperElInfo->getRefSimplexCoords(coords);
- while (lowerLevel > upperLevel) {
- upperElInfo->getSubElementCoords(coords, elinfo_stack[stackPos]->getIChild());
- stackPos--;
- lowerLevel--;
+ for (int i = 1; i <= levelDif; i++) {
+ upperElInfo->getSubElementCoords(coords,
+ elinfo_stack[stack_used - levelDif + i]->getIChild());
}
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index f85388ee7a0909608b1b7ed08c3b6203192f8389..ac9b2a714dd674c8684c269cfbfd90ef48e41309 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -150,7 +150,7 @@ namespace AMDiS {
void update();
void getCoordsInElem(const ElInfo *upperElInfo,
- VectorOfFixVecs<DimVec<double> > *coords);
+ DimMat<double> *coords);
/** \brief
* Starts the calculation of the displacement relative to given level.
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 5361c772ded2ca742ed0845b5c45a511adb325f9..e3579a361e2890fd9002d2f92b016543e41c18eb 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -147,10 +147,33 @@ namespace AMDiS {
void StandardZOA::calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
- const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts(); // grosses Element
- const BasisFunction *phi = owner->getColFESpace()->getBasisFcts(); // kleines Element
+ FUNCNAME("StandardZOA::calculateElementMatrix()");
+
+ TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");
+
+ const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
+ const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
+ DimMat<double> *m = smallElInfo->getSubElemCoordsMat();
+ // m->print();
+ // WAIT_REALLY;
+
+ // The basis function on the larger element is defined as a linear
+ // combination of basis functions of the smaller element. At the moment,
+ // this is supported only for linear lagrange basis functions.
+ //
+ // If b0 and b1 are the basis functions of the larger element restricted to
+ // the small element, and s0 and s1 are the two basis function of the small
+ // element, b0 and b1 are expressed as follows:
+ //
+ // b0 = c00 * s0 + c01 * s1;
+ // b1 = c10 * s0 + c11 * s1;
+ //
+ // The constants are defined by the subElement Matrix.
+
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
@@ -166,26 +189,19 @@ namespace AMDiS {
}
for (int iq = 0; iq < nPoints; iq++) {
- c[iq] *= colElInfo->getDet();
-
- for (int i = 0; i < 2; i++) { // iteration ueber col-phi basis funktionen
- for (int j = 0; j < 2; j++) { // iteration ueber row-psi basis funktionen
- if (j == 0) {
- double val = quadrature->getWeight(iq) * c[iq];
- val *= (*(phi->getPhi(i)))(quadrature->getLambda(iq)) * (*(phi->getPhi(0)))(quadrature->getLambda(iq));
-
- (*mat)[j][i] += val;
+ c[iq] *= smallElInfo->getDet();
- val = quadrature->getWeight(iq) * c[iq];
- val *= (*(phi->getPhi(i)))(quadrature->getLambda(iq)) * 0.5 * (*(phi->getPhi(1)))(quadrature->getLambda(iq));
+ for (int i = 0; i < nCol; i++) {
+ for (int j = 0; j < nRow; j++) {
+ double val = quadrature->getWeight(iq) * c[iq] * (*(phi->getPhi(i)))(quadrature->getLambda(iq));
- (*mat)[j][i] += val;
- } else {
- double val = quadrature->getWeight(iq) * c[iq];
- val *= (*(phi->getPhi(i)))(quadrature->getLambda(iq)) * 0.5 * (*(phi->getPhi(1)))(quadrature->getLambda(iq));
-
- (*mat)[j][i] += val;
+ double tmpval = 0.0;
+ for (int k = 0; k < nCol; k++) {
+ tmpval += (*m)[j][k] * (*(phi->getPhi(k)))(quadrature->getLambda(iq));
}
+ val *= tmpval;
+
+ (*mat)[j][i] += val;
}
}
}
diff --git a/AMDiS/src/ZeroOrderAssembler.h b/AMDiS/src/ZeroOrderAssembler.h
index 4f9c65efa322de6bb4cba7f40a406d02e6e61da3..aa33fcce4344a792ddcb0085dc996ed74cda35a4 100644
--- a/AMDiS/src/ZeroOrderAssembler.h
+++ b/AMDiS/src/ZeroOrderAssembler.h
@@ -83,6 +83,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat);
/** \brief
@@ -117,6 +119,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");
@@ -155,6 +159,8 @@ namespace AMDiS {
void calculateElementMatrix(const ElInfo *rowElInfo,
const ElInfo *colElInfo,
+ const ElInfo *smallElInfo,
+ const ElInfo *largeElInfo,
ElementMatrix *mat)
{
ERROR_EXIT("CEM not yet\n");