From 4dea0e15042a06a00e82ca90023a664c772b9602 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 13 May 2009 08:08:01 +0000
Subject: [PATCH] Better solution to create element matrices for assembling on
different meshes.
---
AMDiS/bin/Makefile.am | 1 -
AMDiS/src/Assembler.cc | 50 ++++++++++++--------
AMDiS/src/DOFIterator.h | 31 ++++++------
AMDiS/src/DOFVector.cc | 30 ++++++------
AMDiS/src/DOFVector.h | 6 ++-
AMDiS/src/DualTraverse.cc | 15 ++++--
AMDiS/src/ElInfo.cc | 11 ++---
AMDiS/src/ElInfo.h | 46 +++++++++---------
AMDiS/src/ElInfo1d.cc | 65 ++++++++++++++------------
AMDiS/src/ElInfo1d.h | 13 ++++--
AMDiS/src/ElInfo2d.cc | 32 +------------
AMDiS/src/ElInfo2d.h | 7 ---
AMDiS/src/ElInfo3d.cc | 2 +
AMDiS/src/ElInfo3d.h | 7 ---
AMDiS/src/FirstOrderAssembler.cc | 26 ++++-------
AMDiS/src/Operator.cc | 78 ++++++++++++-------------------
AMDiS/src/Operator.h | 56 +++++++---------------
AMDiS/src/ProblemVec.cc | 34 +++++---------
AMDiS/src/SecondOrderAssembler.cc | 38 +++++++--------
AMDiS/src/Traverse.cc | 19 ++++----
AMDiS/src/Traverse.h | 3 +-
AMDiS/src/ZeroOrderAssembler.cc | 5 +-
22 files changed, 245 insertions(+), 330 deletions(-)
diff --git a/AMDiS/bin/Makefile.am b/AMDiS/bin/Makefile.am
index 04ad127f..22e25027 100644
--- a/AMDiS/bin/Makefile.am
+++ b/AMDiS/bin/Makefile.am
@@ -16,7 +16,6 @@ if USE_PARALLEL_AMDIS
$(PARALLEL_DIR)/MeshStructure_ED.h \
$(PARALLEL_DIR)/ParallelError.h $(PARALLEL_DIR)/ParallelError.hh \
$(PARALLEL_DIR)/ParallelProblem.h $(PARALLEL_DIR)/ParallelProblem.cc \
- $(PARALLEL_DIR)/ParallelDomainProblem.h $(PARALLEL_DIR)/ParallelDomainProblem.cc \
$(PARALLEL_DIR)/ParMetisPartitioner.h $(PARALLEL_DIR)/ParMetisPartitioner.cc \
$(PARALLEL_DIR)/PartitionElementData.h
PARALLEL_INCLUDES = -I/work/home7/witkowsk/local/include \
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 50a4f276..6f816dab 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -109,22 +109,35 @@ namespace AMDiS {
ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
- TEST_EXIT(zeroOrderAssembler)("Not yet implemented for not zero order assembler!\n");
-
- if (secondOrderAssembler)
+ if (secondOrderAssembler) {
secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
- if (firstOrderAssemblerGrdPsi)
+
+ ElementMatrix m(nRow, nRow);
+ smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
+ ElementMatrix tmpMat(nRow, nRow);
+ tmpMat = m * mat;
+ mat = tmpMat;
+ }
+
+ if (firstOrderAssemblerGrdPsi) {
+ ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
- if (firstOrderAssemblerGrdPhi)
+ }
+
+ if (firstOrderAssemblerGrdPhi) {
+ ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
- if (zeroOrderAssembler)
- zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
+ }
- const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
- ElementMatrix tmpMat(nRow, nRow);
- tmpMat = m * mat;
- mat = tmpMat;
+ if (zeroOrderAssembler) {
+ zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
+ ElementMatrix m(nRow, nRow);
+ smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
+ ElementMatrix tmpMat(nRow, nRow);
+ tmpMat = m * mat;
+ mat = tmpMat;
+ }
if (rememberElMat)
userMat += factor * elementMatrix;
@@ -240,23 +253,20 @@ namespace AMDiS {
FUNCNAME("Assembler::matVecAssemble()");
Element *el = elInfo->getElement();
- const BasisFunction *basFcts = rowFESpace->getBasisFcts();
- int nBasFcts = basFcts->getNumber();
- double *uhOldLoc = new double[nBasFcts];
+ double *uhOldLoc = new double[nRow];
operat->uhOld->getLocalVector(el, uhOldLoc);
if (el != lastMatEl)
calculateElementMatrix(elInfo, elementMatrix);
- for (int i = 0; i < nBasFcts; i++) {
+ for (int i = 0; i < nRow; i++) {
double val = 0.0;
- for (int j = 0; j < nBasFcts; j++) {
+ for (int j = 0; j < nRow; j++) {
val += elementMatrix[i][j] * uhOldLoc[j];
}
vec[i] += val;
- }
-
+ }
delete [] uhOldLoc;
}
@@ -283,7 +293,8 @@ namespace AMDiS {
operat->uhOld->getLocalVector(auxEl, uhOldLoc);
- const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
+ mtl::dense2D<double> m(nRow, nRow);
+ smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
for (int i = 0; i < nBasFcts; i++) {
uhOldLoc2[i] = 0.0;
@@ -306,7 +317,6 @@ namespace AMDiS {
vec[i] += val;
}
-
delete [] uhOldLoc;
delete [] uhOldLoc2;
}
diff --git a/AMDiS/src/DOFIterator.h b/AMDiS/src/DOFIterator.h
index 6be65ae1..a61d3e01 100644
--- a/AMDiS/src/DOFIterator.h
+++ b/AMDiS/src/DOFIterator.h
@@ -72,14 +72,13 @@ namespace AMDiS {
virtual void reset() {
position = 0;
dofFreeIterator = dofFree->begin();
- if(dofFreeIterator == dofFree->end()) {
+ if (dofFreeIterator == dofFree->end())
return;
- }
+
goToBeginOfIteratedObject();
- if(type != ALL_DOFS) {
- if(*dofFreeIterator == (type == USED_DOFS))
+ if (type != ALL_DOFS)
+ if (*dofFreeIterator == (type == USED_DOFS))
operator++();
- }
}
/** \brief
@@ -196,14 +195,20 @@ namespace AMDiS {
virtual void decObjectIterator() {}
protected:
- DOFAdmin *dofAdmin; /**< DOFAdmin which contains
- * the dofFree vector
- */
- int position; /**< current position index */
- std::vector<bool> *dofFree; /**< stores which DOFs are used
- */
- std::vector<bool>::iterator dofFreeIterator; /**< iterator for dofFree */
- const DOFIteratorType type; /**< type of this iterator */
+ /// DOFAdmin which contains the dofFree vector.
+ DOFAdmin *dofAdmin;
+
+ /// Current position index.
+ int position;
+
+ /// Stores which DOFs are used.
+ std::vector<bool> *dofFree;
+
+ /// Iterator for dofFree.
+ std::vector<bool>::iterator dofFreeIterator;
+
+ /// Type of this iterator.
+ const DOFIteratorType type;
};
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 38be5c24..f90020f5 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -1,5 +1,4 @@
#include <boost/numeric/mtl/mtl.hpp>
-
#include "DOFVector.h"
#include "Traverse.h"
#include "DualTraverse.h"
@@ -314,10 +313,8 @@ namespace AMDiS {
TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
- if (quad && quadFast) {
- TEST_EXIT_DBG(quad == quadFast->getQuadrature())
- ("quad != quadFast->quadrature\n");
- }
+ if (quad && quadFast)
+ TEST_EXIT_DBG(quad == quadFast->getQuadrature())("quad != quadFast->quadrature\n");
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
@@ -335,29 +332,29 @@ namespace AMDiS {
delete [] grd;
grd = new WorldVector<double>[nPoints];
- for (int i = 0; i < nPoints; i++) {
+ for (int i = 0; i < nPoints; i++)
grd[i].set(0.0);
- }
+
result = grd;
}
double *localVec = localVectors[myRank];
getLocalVector(largeElInfo->getElement(), localVec);
- const mtl::dense2D<double> m = smallElInfo->getSubElemCoordsMat();
+
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
+ mtl::dense2D<double> m(nBasFcts, nBasFcts);
+ smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();
-
- const BasisFunction *basFcts = feSpace->getBasisFcts();
DimVec<double>* grdPhi = grdPhis[myRank];
DimVec<double> tmp(dim, DEFAULT_VALUE, 0.0);
for (int i = 0; i < nPoints; i++) {
- for (int j = 0; j < dim + 1; j++) {
+ for (int j = 0; j < dim + 1; j++)
grd1[j] = 0.0;
- }
for (int j = 0; j < nBasFcts; j++) {
@@ -368,16 +365,14 @@ namespace AMDiS {
*grdPhi += tmp;
}
- for (int k = 0; k < parts; k++) {
+ for (int k = 0; k < parts; k++)
grd1[k] += (*grdPhi)[k] * localVec[j];
- }
}
for (int l = 0; l < dow; l++) {
result[i][l] = 0.0;
- for (int k = 0; k < parts; k++) {
+ for (int k = 0; k < parts; k++)
result[i][l] += grdLambda[k][l] * grd1[k];
- }
}
}
@@ -804,7 +799,8 @@ namespace AMDiS {
double *localVec = localVectors[omp_get_thread_num()];
getLocalVector(largeElInfo->getElement(), localVec);
- const mtl::dense2D<double> m = smallElInfo->getSubElemCoordsMat();
+ mtl::dense2D<double> m(nBasFcts, nBasFcts);
+ smallElInfo->getSubElemCoordsMat(m, basFcts->getDegree());
for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index ac16ebd9..123c8c49 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -290,11 +290,11 @@ namespace AMDiS {
{}
DOFVector<T> *create() {
- return NEW DOFVector<T>(feSpace, "");
+ return new DOFVector<T>(feSpace, "");
}
void free(DOFVector<T> *vec) {
- DELETE vec;
+ delete vec;
}
private:
@@ -493,12 +493,14 @@ namespace AMDiS {
void interpol(DOFVector<T> *v, double factor);
+ /// Writes the data of the DOFVector to an output stream.
void serialize(std::ostream &out) {
unsigned int size = vec.size();
out.write(reinterpret_cast<const char*>(&size), sizeof(unsigned int));
out.write(reinterpret_cast<const char*>(&(vec[0])), size * sizeof(T));
}
+ /// Reads data of an DOFVector from an input stream.
void deserialize(std::istream &in) {
unsigned int size;
in.read(reinterpret_cast<char*>(&size), sizeof(unsigned int));
diff --git a/AMDiS/src/DualTraverse.cc b/AMDiS/src/DualTraverse.cc
index aed16a83..95811279 100644
--- a/AMDiS/src/DualTraverse.cc
+++ b/AMDiS/src/DualTraverse.cc
@@ -170,12 +170,19 @@ namespace AMDiS {
if (!fillSubElemMat)
return;
- mtl::dense2D<double>& subCoordsMat =
- const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
+// mtl::dense2D<double>& subCoordsMat =
+// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
+// mtl::dense2D<double>& subCoordsMat_so =
+// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());
+
if (elInfo1 == elInfoSmall) {
- stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
+// stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
+// stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
+ stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
} else {
- stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
+// stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
+// stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
+ stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
}
}
}
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 7431e21c..bfc2430c 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -29,7 +29,8 @@ namespace AMDiS {
neighbourCoord_(mesh_->getDim(), NO_INIT),
oppVertex_(mesh_->getDim(), NO_INIT),
grdLambda(mesh_->getDim(), NO_INIT),
- subElemCoordsMat(2, 2)
+ refinementPath(0),
+ refinementPathLength(0)
{
projection_.set(NULL);
@@ -46,14 +47,12 @@ namespace AMDiS {
WorldVector<double>& w) const
{
- int dim = l.getSize() - 1;
-
double c = l[0];
for (int j = 0; j < dimOfWorld; j++)
w[j] = c * coord_[0][j];
- int vertices = Global::getGeo(VERTEX, dim);
+ int vertices = Global::getGeo(VERTEX, l.getSize() - 1);
for (int i = 1; i < vertices; i++) {
c = l[i];
@@ -62,10 +61,6 @@ namespace AMDiS {
}
}
- /****************************************************************************/
- /* compute volume of an element */
- /****************************************************************************/
-
double ElInfo::calcDet() const
{
testFlag(Mesh::FILL_COORDS);
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 4a9fe862..45b3bc92 100644
--- a/AMDiS/src/ElInfo.h
+++ b/AMDiS/src/ElInfo.h
@@ -205,13 +205,11 @@ namespace AMDiS {
return parametric_;
}
- /** \brief
- * Returns the element transformation matrix \ref subElemCoordsMat .
- * This value is set only during dual traverse.
- */
- inline const mtl::dense2D<double>& getSubElemCoordsMat() const {
- return subElemCoordsMat;
- }
+ virtual void getSubElemCoordsMat(mtl::dense2D<double>& mat,
+ int basisFctsDegree) const {}
+
+ virtual void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
+ int basisFctsDegree) const {}
/** \} */
@@ -279,14 +277,18 @@ namespace AMDiS {
parametric_ = param;
}
- /** \brief
- * Sets the element transformation matrix \ref subElemCoordsMat .
- * This value is used only during dual traverse.
- */
-// inline void setSubElemCoordsMat(DimMat<double> *mat) {
-// subElemCoordsMat = mat;
-// }
-
+ /// Set \ref refinementPath
+ inline void setRefinementPath(unsigned long rPath)
+ {
+ refinementPath = rPath;
+ }
+
+ /// Set \ref refinementPathLength
+ inline void setRefinementPathLength(int length)
+ {
+ refinementPathLength = length;
+ }
+
/** \} */
@@ -372,14 +374,6 @@ namespace AMDiS {
return(0.0);
}
- virtual void getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const = 0;
-
- virtual void getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const = 0;
-
-
protected:
/// Pointer to the current mesh
Mesh *mesh_;
@@ -471,6 +465,10 @@ namespace AMDiS {
/// Stores the world dimension.
int dimOfWorld;
+ unsigned long refinementPath;
+
+ int refinementPathLength;
+
/** \brief
* This is a transformation matrix used during dual traverse. It is set, if
* the current element is the smaller element of an element pair in the traverse.
@@ -480,6 +478,8 @@ namespace AMDiS {
*/
mtl::dense2D<double> subElemCoordsMat;
+ mtl::dense2D<double> subElemCoordsMat_so;
+
public:
/** \brief
* child_vertex[el_type][child][i] = father's local vertex index of new
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 060a038f..0e631c07 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -40,6 +40,10 @@ namespace AMDiS {
{1.0, 0.75, 0.0},
{0.0, 0.375, 1.0}};
mtl::dense2D<double> ElInfo1d::mat_d2_right(mat_d2_right_val);
+
+ double ElInfo1d::test1_val[2][2] = {{1.0, 0.0},
+ {0.0, 1.0}};
+ mtl::dense2D<double> ElInfo1d::test2_val(test1_val);
void ElInfo1d::fillMacroInfo(const MacroElement * mel)
@@ -122,11 +126,10 @@ namespace AMDiS {
if (adet2 < 1.0E-15) {
MSG("det*det = %lf\n", adet2);
- for (int i = 0; i <= 1; i++)
- grd_lam[i] = 0.0;
+ grd_lam[0] = grd_lam[1] = 0.0;
} else {
grd_lam[1] = e * (1.0 / adet2);
- grd_lam[0] = grd_lam[1] * (-1.0);
+ grd_lam[0] = grd_lam[1] * -1.0;
}
return sqrt(adet2);
@@ -296,48 +299,52 @@ namespace AMDiS {
}
}
- void ElInfo1d::getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const
+ void ElInfo1d::getSubElemCoordsMat(mtl::dense2D<double>& mat, int basisFctsDegree) const
{
- FUNCNAME("ElInfo1d::getRefSimplexCoords()");
-
- switch (basisFcts->getDegree()) {
+ switch (basisFctsDegree) {
case 1:
- coords = mat_d1;
+ mat = mat_d1;
+
+ for (int i = 0; i < refinementPathLength; i++) {
+ if (refinementPath & (1 << i))
+ mat *= mat_d1_left;
+ else
+ mat *= mat_d1_right;
+ }
+
break;
+
case 2:
- coords = mat_d2;
+ mat = mat_d2;
+
+ for (int i = 0; i < refinementPathLength; i++) {
+ if (refinementPath & (1 << i))
+ mat *= mat_d2_left;
+ else
+ mat *= mat_d2_right;
+ }
+
break;
+
default:
- ERROR_EXIT("Not implemented for basis function with degree %d!\n",\
- basisFcts->getDegree());
+ ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
}
}
- void ElInfo1d::getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const
+ void ElInfo1d::getSubElemCoordsMat_so(mtl::dense2D<double>& mat, int basisFctsDegree) const
{
- FUNCNAME("ElInfo1d::getSubElementCoords()");
-
- switch (basisFcts->getDegree()) {
+ switch (basisFctsDegree) {
case 1:
- if (iChild == 0)
- coords *= mat_d1_left;
- else
- coords *= mat_d1_right;
+ mat = test2_val;
- break;
- case 2:
- if (iChild == 0)
- coords *= mat_d2_left;
- else
- coords *= mat_d2_right;
+ for (int i = 0; i < refinementPathLength; i++)
+ mat *= 0.5;
break;
default:
- ERROR_EXIT("Not yet implemented!\n");
+ ERROR_EXIT("Not supported for basis function degree: %d\n", basisFctsDegree);
}
}
+
}
diff --git a/AMDiS/src/ElInfo1d.h b/AMDiS/src/ElInfo1d.h
index 11294a3a..e91edfba 100644
--- a/AMDiS/src/ElInfo1d.h
+++ b/AMDiS/src/ElInfo1d.h
@@ -61,12 +61,11 @@ namespace AMDiS {
return (i + 1) % 2;
}
- void getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const;
+ void getSubElemCoordsMat(mtl::dense2D<double>& mat,
+ int basisFctsDegree) const;
- void getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const;
+ void getSubElemCoordsMat_so(mtl::dense2D<double>& mat,
+ int basisFctsDegree) const;
protected:
static double mat_d1_val[2][2];
@@ -93,6 +92,10 @@ namespace AMDiS {
static mtl::dense2D<double> mat_d2_right;
+ static double test1_val[2][2];
+
+ static mtl::dense2D<double> test2_val;
+
};
}
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 5524fa4b..8bcea08e 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -58,9 +58,8 @@ namespace AMDiS {
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
int vertices = mesh_->getGeo(VERTEX);
- for (int i = 0; i < vertices; i++) {
+ for (int i = 0; i < vertices; i++)
coord_[i] = mel->coord[i];
- }
}
int neighbours = mesh_->getGeo(NEIGH);
@@ -619,33 +618,4 @@ namespace AMDiS {
return det;
}
- void ElInfo2d::getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const
- {
- TEST_EXIT(basisFcts->getDegree() == 1)("Wrong basis function degree!\n");
-
- coords = mat_d1;
- }
-
- void ElInfo2d::getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const
- {
- FUNCNAME("ElInfo1d::getSubElementCoords()");
-
- switch (basisFcts->getDegree()) {
- case 1:
- if (iChild == 0)
- coords *= mat_d1_left;
- else
- coords *= mat_d1_right;
-
- break;
- case 2:
- break;
- default:
- ERROR_EXIT("Not yet implemented!\n");
- }
- }
-
}
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index 29a1ca26..797597ef 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/AMDiS/src/ElInfo2d.h
@@ -58,13 +58,6 @@ namespace AMDiS {
/// 2-dimensional realisation of ElInfo's getElementNormal method.
double getElementNormal(WorldVector<double> &normal) const;
- void getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const;
-
- void getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const;
-
protected:
/// Temp vectors for function \ref calcGrdLambda.
WorldVector<double> *e1, *e2, *normal;
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index d8ae8dab..55567f35 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -609,6 +609,7 @@ namespace AMDiS {
}
}
+#if 0
void ElInfo3d::getRefSimplexCoords(const BasisFunction *basisFcts,
mtl::dense2D<double>& coords) const
{
@@ -670,4 +671,5 @@ namespace AMDiS {
// (*coords)[3][3] = c3;
}
+#endif
}
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index 9cd217f2..b4412a08 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/AMDiS/src/ElInfo3d.h
@@ -87,13 +87,6 @@ namespace AMDiS {
orientation = o;
}
- void getRefSimplexCoords(const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords) const;
-
- void getSubElementCoords(const BasisFunction *basisFcts,
- int iChild,
- mtl::dense2D<double>& coords) const;
-
protected:
/// \ref el 's type. Is Filled automatically by the traversal routines.
unsigned char el_type;
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index 7df4bcb2..dcb19c3c 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -38,19 +38,17 @@ namespace AMDiS {
bool optimized)
{
std::vector<SubAssembler*> *subAssemblers =
- optimized ?
- (type == GRD_PSI ?
- &optimizedSubAssemblersGrdPsi :
- &optimizedSubAssemblersGrdPhi) :
- (type == GRD_PSI ?
- &standardSubAssemblersGrdPsi :
- &standardSubAssemblersGrdPhi);
-
+ optimized ?
+ (type == GRD_PSI ?
+ &optimizedSubAssemblersGrdPsi :
+ &optimizedSubAssemblersGrdPhi) :
+ (type == GRD_PSI ?
+ &standardSubAssemblersGrdPsi :
+ &standardSubAssemblersGrdPhi);
+
int myRank = omp_get_thread_num();
std::vector<OperatorTerm*> opTerms
- = (type == GRD_PSI) ?
- op->firstOrderGrdPsi[myRank] :
- op->firstOrderGrdPhi[myRank];
+ = (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank];
// check if a assembler is needed at all
if (opTerms.size() == 0)
@@ -62,16 +60,12 @@ namespace AMDiS {
// check if a new assembler is needed
for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
-
std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
sort(assTerms.begin(), assTerms.end());
- if ((opTerms == assTerms) &&
- ((*subAssemblers)[i]->getQuadrature() == quad)) {
-
+ if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
- }
}
// check if all terms are pw_const
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 427bf073..79ca035a 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -227,11 +227,9 @@ namespace AMDiS {
{
const int dimOfWorld = Global::getGeo(WORLD);
int dim = LALt.getNumRows() - 1;
-
- double val = 0.0;
for (int i = 0; i <= dim; i++) {
- val = 0.0;
+ double val = 0.0;
for (int k = 0; k < dimOfWorld; k++)
val += Lambda[i][k] * Lambda[i][k];
val *= factor;
@@ -307,9 +305,8 @@ namespace AMDiS {
{
int myRank = omp_get_thread_num();
- if (!assembler[myRank]) {
+ if (!assembler[myRank])
initAssembler(myRank, NULL, NULL, NULL, NULL);
- }
assembler[myRank]->calculateElementMatrix(elInfo, userMat, factor);
}
@@ -321,9 +318,8 @@ namespace AMDiS {
{
int myRank = omp_get_thread_num();
- if (!assembler[myRank]) {
+ if (!assembler[myRank])
initAssembler(myRank, NULL, NULL, NULL, NULL);
- }
assembler[myRank]->calculateElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo,
@@ -336,9 +332,8 @@ namespace AMDiS {
{
int myRank = omp_get_thread_num();
- if (!assembler[myRank]) {
+ if (!assembler[myRank])
initAssembler(myRank, NULL, NULL, NULL, NULL);
- }
assembler[myRank]->calculateElementVector(elInfo, userVec, factor);
}
@@ -350,9 +345,8 @@ namespace AMDiS {
{
int myRank = omp_get_thread_num();
- if (!assembler[myRank]) {
+ if (!assembler[myRank])
initAssembler(myRank, NULL, NULL, NULL, NULL);
- }
assembler[myRank]->calculateElementVector(mainElInfo, auxElInfo,
smallElInfo, largeElInfo,
@@ -370,15 +364,13 @@ namespace AMDiS {
#endif
{
if (optimized) {
- assembler[rank] =
- NEW OptimizedAssembler(this,
- quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFESpace, colFESpace);
+ assembler[rank] = new OptimizedAssembler(this,
+ quad2, quad1GrdPsi, quad1GrdPhi, quad0,
+ rowFESpace, colFESpace);
} else {
- assembler[rank] =
- NEW StandardAssembler(this,
- quad2, quad1GrdPsi, quad1GrdPhi, quad0,
- rowFESpace, colFESpace);
+ assembler[rank] = new StandardAssembler(this,
+ quad2, quad1GrdPsi, quad1GrdPhi, quad0,
+ rowFESpace, colFESpace);
}
}
}
@@ -1194,83 +1186,72 @@ namespace AMDiS {
void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
- }
}
void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
- }
}
void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
- }
}
void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
- }
}
void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
- }
}
void FctGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
- }
}
void VecAndGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
- }
}
void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
- }
}
void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt)
const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
- }
}
void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
lalt(Lambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
- }
}
void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
- }
}
void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
@@ -2924,9 +2905,9 @@ namespace AMDiS {
Quadrature *quad)
{
int size = static_cast<int>(vecs.size());
- for (int i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++)
gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
- }
+
vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
}
@@ -2937,9 +2918,9 @@ namespace AMDiS {
std::vector<WorldVector<double>*> arg(size);
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++)
arg[i] = &(gradsAtQPs[i][iq]);
- }
+
C[iq] += (*f)(vecAtQPs[iq], arg);
}
}
@@ -2955,9 +2936,9 @@ namespace AMDiS {
std::vector<WorldVector<double>*> arg(size);
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++)
arg[i] = &(gradsAtQPs[i][iq]);
- }
+
result[iq] += fac * (*f)(vecAtQPs[iq], arg) * uhAtQP[iq];
}
}
@@ -2982,9 +2963,8 @@ namespace AMDiS {
void VecGrad_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
- }
}
void VecGrad_FOT::eval(int nPoints,
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 530a2744..6378e98b 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -157,11 +157,8 @@ namespace AMDiS {
DimMat<double>& LALt,
double factor);
- /** \brief
- * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the
- * identity.
- */
- static void l1lt(const DimVec<WorldVector<double> >& Lambda,
+ /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
+ static void l1lt(const DimVec<WorldVector<double> >& Lambda,
DimMat<double>& LALt,
double factor);
@@ -171,10 +168,7 @@ namespace AMDiS {
DimVec<double>& Lb,
double factor);
- /** \brief
- * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in
- * each component.
- */
+ /// Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in each component.
static void l1(const DimVec<WorldVector<double> >& Lambda,
DimVec<double>& Lb,
double factor)
@@ -185,9 +179,9 @@ namespace AMDiS {
for (int i = 0; i <= dim; i++) {
double val = 0.0;
- for (int j = 0; j < dimOfWorld; j++) {
+ for (int j = 0; j < dimOfWorld; j++)
val += Lambda[i][j];
- }
+
val *= factor;
Lb[i] += val;
}
@@ -219,10 +213,6 @@ namespace AMDiS {
friend class Operator;
};
- // ============================================================================
- // ===== class SecondOrderTerm ===============================================
- // ============================================================================
-
/**
* \ingroup Assembler
*
@@ -237,22 +227,16 @@ namespace AMDiS {
: OperatorTerm(deg)
{}
- /** \brief
- * Destructor.
- */
+ /// Destructor.
virtual ~SecondOrderTerm()
{}
- /** \brief
- * Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
- */
+ /// Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
virtual void getLALt(const ElInfo *elInfo,
int nPoints,
DimMat<double> **result) const = 0;
- /** \brief
- * Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
- */
+ /// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
virtual void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const = 0;
@@ -275,24 +259,19 @@ namespace AMDiS {
setSymmetric(true);
}
- /** \brief
- * Implenetation of SecondOrderTerm::getLALt().
- */
+ /// Implenetation of SecondOrderTerm::getLALt().
inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
l1lt(Lambda, *(LALt[iq]), 1.0);
- }
}
- /** \brief
- * Implementation of SecondOrderTerm::eval().
- */
+ /// Implementation of SecondOrderTerm::eval().
inline void eval(int nPoints,
- const double * , // uhAtQP
- const WorldVector<double> * , // grdUhAtQP
+ const double *,
+ const WorldVector<double> *,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const
@@ -310,17 +289,14 @@ namespace AMDiS {
}
}
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if (grdUhAtQP) {
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
result[iq] += grdUhAtQP[iq];
- }
}
}
};
@@ -1467,7 +1443,7 @@ namespace AMDiS {
FactorSimple_FOT(double f)
: FirstOrderTerm(0)
{
- factor = NEW double;
+ factor = new double;
*factor = f;
}
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 43c563ea..e59fe7c9 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -687,6 +687,7 @@ 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
@@ -695,18 +696,12 @@ namespace AMDiS {
// If the matrix was assembled before and it is marked to be assembled
// only once, it will not be assembled.
- if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
+ if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j])
assembleMatrix = false;
- }
// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
- if (!matrix) {
+ if (!matrix || !assembleMatrix)
assembleMatrix = false;
- }
-
- if (!asmMatrix) {
- assembleMatrix = false;
- }
// If the matrix should not be assembled, the rhs vector has to be considered.
// This will be only done, if i == j. So, if both is not true, we can jump
@@ -1055,23 +1050,20 @@ namespace AMDiS {
#pragma omp barrier
#endif
while (elInfo) {
- if (useGetBound) {
+ if (useGetBound)
basisFcts->getBound(elInfo, bound);
- }
if (matrix) {
tmpMatrix->assemble(1.0, elInfo, bound);
// Take the matrix boundary manager from the public matrix,
// but assemble the boundary conditions on the thread private matrix.
- if (matrix->getBoundaryManager()) {
- matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpMatrix);
- }
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpMatrix);
}
- if (vector) {
+ if (vector)
tmpVector->assemble(1.0, elInfo, bound, NULL);
- }
elInfo = stack.traverseNext(elInfo);
}
@@ -1098,9 +1090,8 @@ namespace AMDiS {
#pragma omp master
#endif
{
- if (matrix) {
+ if (matrix)
matrix->startInsertion();
- }
}
if (matrix) {
@@ -1160,15 +1151,12 @@ namespace AMDiS {
if (matrix) {
matrix->assemble(1.0, rowElInfo, colElInfo, smallElInfo, largeElInfo, bound);
- if (matrix->getBoundaryManager()) {
- matrix->getBoundaryManager()->
- fillBoundaryConditions(rowElInfo, matrix);
- }
+ if (matrix->getBoundaryManager())
+ matrix->getBoundaryManager()->fillBoundaryConditions(rowElInfo, matrix);
}
- if (vector) {
+ if (vector)
vector->assemble(1.0, rowElInfo, bound);
- }
cont = dualTraverse.traverseNext(&rowElInfo, &colElInfo,
&smallElInfo, &largeElInfo);
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index 603b9f19..1a9e6ed8 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -36,11 +36,11 @@ namespace AMDiS {
SecondOrderAssembler *newAssembler;
std::vector<SubAssembler*> *subAssemblers =
- optimized ?
- &optimizedSubAssemblers :
- &standardSubAssemblers;
+ optimized ?
+ &optimizedSubAssemblers :
+ &standardSubAssemblers;
- std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
+ std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
sort(opTerms.begin(), opTerms.end());
@@ -50,11 +50,8 @@ namespace AMDiS {
sort(assTerms.begin(), assTerms.end());
- if ((opTerms == assTerms) &&
- ((*subAssemblers)[i]->getQuadrature() == quad)) {
-
+ if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
- }
}
// check if all terms are pw_const
@@ -68,12 +65,12 @@ namespace AMDiS {
// create new assembler
if (!optimized) {
- newAssembler = NEW Stand2(op, assembler, quad);
+ newAssembler = new Stand2(op, assembler, quad);
} else {
if (pwConst) {
- newAssembler = NEW Pre2(op, assembler, quad);
+ newAssembler = new Pre2(op, assembler, quad);
} else {
- newAssembler = NEW Quad2(op, assembler, quad);
+ newAssembler = new Quad2(op, assembler, quad);
}
}
@@ -89,8 +86,8 @@ namespace AMDiS {
quadrature);
tmpLALt.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) {
- tmpLALt[i] = NEW DimMat<double>*;
- *(tmpLALt[i]) = NEW DimMat<double>(dim, NO_INIT);
+ tmpLALt[i] = new DimMat<double>*;
+ *(tmpLALt[i]) = new DimMat<double>(dim, NO_INIT);
}
}
@@ -186,9 +183,9 @@ namespace AMDiS {
int myRank = omp_get_thread_num();
if (firstCall) {
- tmpLALt[myRank] = NEW DimMat<double>*[nPoints];
+ tmpLALt[myRank] = new DimMat<double>*[nPoints];
for (int j = 0; j < nPoints; j++) {
- tmpLALt[myRank][j] = NEW DimMat<double>(dim, NO_INIT);
+ tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
}
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
@@ -262,9 +259,9 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
- DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
+ DimMat<double> **LALt = new DimMat<double>*[nPoints];
for (int iq = 0; iq < nPoints; iq++) {
- LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
+ LALt[iq] = new DimMat<double>(dim, NO_INIT);
LALt[iq]->set(0.0);
}
@@ -280,12 +277,12 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
- for (int i = 0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++)
(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
- }
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
+
mat[i][i] += quadrature->getWeight(iq) *
(grdPsi * ((*LALt[iq]) * grdPhi[i]));
@@ -301,9 +298,8 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
- for (int i = 0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++)
(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
- }
for (int i = 0; i < nRow; i++) {
(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index b972ec70..6dee2f61 100644
--- a/AMDiS/src/Traverse.cc
+++ b/AMDiS/src/Traverse.cc
@@ -1068,18 +1068,17 @@ namespace AMDiS {
dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
}
- void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo,
- const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords)
+ void TraverseStack::fillRefinementPath(ElInfo& elInfo, const ElInfo& upperElInfo)
{
- int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo->getLevel();
+ int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo.getLevel();
+ unsigned long rPath = 0;
- upperElInfo->getRefSimplexCoords(basisFcts, coords);
-
- for (int i = 1; i <= levelDif; i++)
- upperElInfo->getSubElementCoords(basisFcts,
- elinfo_stack[stack_used - levelDif + i]->getIChild(),
- coords);
+ for (int i = 1; i <= levelDif; i++)
+ if (elinfo_stack[stack_used - levelDif + i]->getIChild())
+ rPath = rPath | (1 << (i - 1));
+
+ elInfo.setRefinementPath(rPath);
+ elInfo.setRefinementPathLength(levelDif);
}
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 56c86ead..c95695b1 100644
--- a/AMDiS/src/Traverse.h
+++ b/AMDiS/src/Traverse.h
@@ -114,8 +114,7 @@ namespace AMDiS {
/// Only for 3d: Calls update of all ElInfo3d onjects in \ref elinfo_stack
void update();
- void getCoordsInElem(const ElInfo *upperElInfo, const BasisFunction *basisFcts,
- mtl::dense2D<double>& coords);
+ void fillRefinementPath(ElInfo& elInfo, const ElInfo& upperElInfo);
/// Is used for parallel mesh traverse.
inline void setMyThreadId(int myThreadId) {
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index b93e7b73..e3d26969 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -26,8 +26,10 @@ namespace AMDiS {
Quadrature *quad,
bool optimized)
{
+ int myRank = omp_get_thread_num();
+
// check if a assembler is needed at all
- if (op->zeroOrder.size() == 0)
+ if (op->zeroOrder[myRank].size() == 0)
return NULL;
ZeroOrderAssembler *newAssembler;
@@ -35,7 +37,6 @@ namespace AMDiS {
std::vector<SubAssembler*> *subAssemblers =
optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
- int myRank = omp_get_thread_num();
std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
sort(opTerms.begin(), opTerms.end());
--
GitLab