From 561814e24d08ab26a7554d4a21f71cce36638092 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 20 Nov 2008 08:35:47 +0000
Subject: [PATCH] * This and that
---
AMDiS/src/DOFVector.cc | 121 ++++++++++---------------
AMDiS/src/DOFVector.h | 16 ++++
AMDiS/src/DOFVector.hh | 150 ++++++++++++-------------------
AMDiS/src/ElInfo3d.cc | 59 ++++++++++++
AMDiS/src/ElInfo3d.h | 15 +---
AMDiS/src/ProblemScal.cc | 2 -
AMDiS/src/SolutionDataStorage.h | 21 +++++
AMDiS/src/SolutionDataStorage.hh | 39 +++++++-
8 files changed, 239 insertions(+), 184 deletions(-)
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 9bfb7d97..90ba3a08 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -149,12 +149,12 @@ namespace AMDiS {
DOFVector<WorldVector<double> > *result = grad;
- if(!result) {
- if(vec && vec->getFESpace() != feSpace) {
+ if (!result) {
+ if (vec && vec->getFESpace() != feSpace) {
DELETE vec;
vec = NULL;
}
- if(!vec) {
+ if (!vec) {
vec = NEW DOFVector<WorldVector<double> >(feSpace, "gradient");
}
result = vec;
@@ -165,36 +165,27 @@ namespace AMDiS {
DOFVector<double> volume(feSpace, "volume");
volume.set(0.0);
- Mesh *mesh = feSpace->getMesh();
-
- int dim = mesh->getDim();
-
const BasisFunction *basFcts = feSpace->getBasisFcts();
-
- DOFAdmin *admin = feSpace->getAdmin();
-
- int numPreDOFs = admin->getNumberOfPreDOFs(0);
-
+ int nPreDOFs = feSpace->getAdmin()->getNumberOfPreDOFs(0);
DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));
+ WorldVector<double> grd;
// traverse mesh
+ Mesh *mesh = feSpace->getMesh();
TraverseStack stack;
-
- Flag fillFlag =
- Mesh::CALL_LEAF_EL | Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
-
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);
while (elInfo) {
double det = elInfo->getDet();
const DegreeOfFreedom **dof = elInfo->getElement()->getDOF();
const double *localUh = getLocalVector(elInfo->getElement(), NULL);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
- WorldVector<double> grd;
basFcts->evalGrdUh(bary, grdLambda, localUh, &grd);
for (int i = 0; i < dim + 1; i++) {
- DegreeOfFreedom dofIndex = dof[i][numPreDOFs];
+ DegreeOfFreedom dofIndex = dof[i][nPreDOFs];
(*result)[dofIndex] += grd * det;
volume[dofIndex] += det;
}
@@ -207,7 +198,7 @@ namespace AMDiS {
for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt) {
if (*volIt != 0.0) {
- *grdIt *= 1.0/(*volIt);
+ *grdIt *= 1.0 / (*volIt);
}
}
@@ -232,12 +223,9 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
- Element *el = elInfo->getElement();
- int dim = elInfo->getMesh()->getDim();
int dow = Global::getGeo(WORLD);
- const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
- const BasisFunction *basFcts = feSpace->getBasisFcts();
- int numPoints = quadrature->getNumPoints();
+ int myRank = omp_get_thread_num();
+ int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldVector<double> *grd = NULL;
WorldVector<double> *result;
@@ -247,22 +235,22 @@ namespace AMDiS {
if (grd)
delete [] grd;
- grd = new WorldVector<double>[numPoints];
- for (int i = 0; i < numPoints; i++) {
+ grd = new WorldVector<double>[nPoints];
+ for (int i = 0; i < nPoints; i++) {
grd[i].set(0.0);
}
result = grd;
}
- double *localVec = new double[nBasFcts];
- getLocalVector(el, localVec);
+ double *localVec = localVectors[myRank];
+ getLocalVector(elInfo->getElement(), localVec);
- DimVec<double> grd1(dim, NO_INIT);
+ DimVec<double> &grd1 = *grdTmp[myRank];
int parts = Global::getGeo(PARTS, dim);
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (quadFast) {
- for (int i = 0; i < numPoints; i++) {
+ for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < dim + 1; j++) {
grd1[j] = 0.0;
}
@@ -281,10 +269,10 @@ namespace AMDiS {
}
}
} else {
- // DimVec<double> grdPhi(dim, DEFAULT_VALUE, 0.0);
- DimVec<double>* grdPhi = grdPhis[omp_get_thread_num()];
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
+ DimVec<double>* grdPhi = grdPhis[myRank];
- for (int i = 0; i < numPoints; i++) {
+ for (int i = 0; i < nPoints; i++) {
for (int j = 0; j < dim + 1; j++) {
grd1[j] = 0.0;
}
@@ -305,8 +293,6 @@ namespace AMDiS {
}
}
- delete [] localVec;
-
return const_cast<const WorldVector<double>*>(result);
}
@@ -330,31 +316,25 @@ namespace AMDiS {
Element *el = elInfo->getElement();
- int dim = elInfo->getMesh()->getDim();
+ int myRank = omp_get_thread_num();
int dow = Global::getGeo(WORLD);
-
- const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
- const BasisFunction *basFcts = feSpace->getBasisFcts();
-
- int numPoints = quadrature->getNumPoints();
- int i, j, k, l, iq;
+ int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static WorldMatrix<double> *vec = NULL;
-
WorldMatrix<double> *result;
if (d2AtQPs) {
result = d2AtQPs;
} else {
if(vec) delete [] vec;
- vec = new WorldMatrix<double>[numPoints];
- for(i = 0; i < numPoints; i++) {
+ vec = new WorldMatrix<double>[nPoints];
+ for (int i = 0; i < nPoints; i++) {
vec[i].set(0.0);
}
result = vec;
}
- double *localVec = new double[nBasFcts];
+ double *localVec = localVectors[myRank];
getLocalVector(el, localVec);
DimMat<double> D2Tmp(dim, DEFAULT_VALUE, 0.0);
@@ -362,56 +342,53 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (quadFast) {
- for (iq = 0; iq < numPoints; iq++) {
- for (k = 0; k < parts; k++)
- for (l = 0; l < parts; l++)
+ for (int iq = 0; iq < nPoints; iq++) {
+ for (int k = 0; k < parts; k++)
+ for (int l = 0; l < parts; l++)
D2Tmp[k][l] = 0.0;
- for (i = 0; i < nBasFcts; i++) {
- for (k = 0; k < parts; k++)
- for (l = 0; l < parts; l++)
- D2Tmp[k][l] += localVec[i]* quadFast->getSecDer(iq, i, k, l);
+ for (int i = 0; i < nBasFcts; i++) {
+ for (int k = 0; k < parts; k++)
+ for (int l = 0; l < parts; l++)
+ D2Tmp[k][l] += localVec[i] * quadFast->getSecDer(iq, i, k, l);
}
- for (i = 0; i < dow; i++)
- for (j = 0; j < dow; j++) {
+ for (int i = 0; i < dow; i++)
+ for (int j = 0; j < dow; j++) {
result[iq][i][j] = 0.0;
- for (k = 0; k < parts; k++)
- for (l = 0; l < parts; l++)
+ for (int k = 0; k < parts; k++)
+ for (int l = 0; l < parts; l++)
result[iq][i][j] += grdLambda[k][i]*grdLambda[l][j]*D2Tmp[k][l];
}
}
} else {
- // DimMat<double> D2Phi(dim, NO_INIT);
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
DimMat<double>* D2Phi = D2Phis[omp_get_thread_num()];
- for (iq = 0; iq < numPoints; iq++) {
- for (k = 0; k < parts; k++)
- for (l = 0; l < parts; l++)
+ for (int iq = 0; iq < nPoints; iq++) {
+ for (int k = 0; k < parts; k++)
+ for (int l = 0; l < parts; l++)
D2Tmp[k][l] = 0.0;
- for (i = 0; i < nBasFcts; i++) {
+ for (int i = 0; i < nBasFcts; i++) {
WARNING("not tested after index correction\n");
(*(basFcts->getD2Phi(i)))(quad->getLambda(iq), *D2Phi);
- for (k = 0; k < parts; k++)
- for (l = 0; l < parts; l++)
+ for (int k = 0; k < parts; k++)
+ for (int l = 0; l < parts; l++)
D2Tmp[k][l] += localVec[i] * (*D2Phi)[k][l];
}
- for (i = 0; i < dow; i++)
- for (j = 0; j < dow; j++) {
+ for (int i = 0; i < dow; i++)
+ for (int j = 0; j < dow; j++) {
result[iq][i][j] = 0.0;
- for (k = 0; k < parts; k++)
- for (l = 0; l < parts; l++)
+ for (int k = 0; k < parts; k++)
+ for (int l = 0; l < parts; l++)
result[iq][i][j] += grdLambda[k][i] * grdLambda[l][j] * D2Tmp[k][l];
}
}
}
- delete [] localVec;
-
-
return const_cast<const WorldMatrix<double>*>(result);
}
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index c54d0f0a..7bfa7ed1 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -231,6 +231,17 @@ namespace AMDiS {
*/
std::vector<DegreeOfFreedom*> localIndices;
+ /** \brief
+ * Are used to store temporary local values of an element.
+ * Thread safe.
+ */
+ std::vector<T*> localVectors;
+
+ /** \brief
+ * Temporarly used in \ref getGrdAtQPs. Thread safe.
+ */
+ std::vector<DimVec<double>*> grdTmp;
+
/** \brief
* Temporarly used in \ref getGrdAtQPs. Thread safe.
*/
@@ -240,6 +251,11 @@ namespace AMDiS {
* Temporarly used in \ref getD2AtQPs. Thread safe.
*/
std::vector<DimMat<double>*> D2Phis;
+
+ /** \brief
+ * Dimension of the mesh this DOFVectorBase belongs to
+ */
+ int dim;
};
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 5a920133..192f2fba 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -29,15 +29,19 @@ namespace AMDiS {
boundaryManager(NULL)
{
nBasFcts = feSpace->getBasisFcts()->getNumber();
- int dim = feSpace->getMesh()->getDim();
+ dim = feSpace->getMesh()->getDim();
localIndices.resize(omp_get_overall_max_threads());
+ localVectors.resize(omp_get_overall_max_threads());
grdPhis.resize(omp_get_overall_max_threads());
+ grdTmp.resize(omp_get_overall_max_threads());
D2Phis.resize(omp_get_overall_max_threads());
for (int i = 0; i < omp_get_overall_max_threads(); i++) {
- localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
+ localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
+ localVectors[i] = GET_MEMORY(T, this->nBasFcts);
grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
+ grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
}
}
@@ -47,7 +51,9 @@ namespace AMDiS {
{
for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts);
+ FREE_MEMORY(localVectors[i], T, this->nBasFcts);
DELETE grdPhis[i];
+ DELETE grdTmp[i];
DELETE D2Phis[i];
}
}
@@ -540,26 +546,19 @@ namespace AMDiS {
return 0;
}
- // das hat Andreas eingefuegt: Integral...
-
template<typename T>
double DOFVector<T>::Int(Quadrature* q) const
{
- FUNCNAME("DOFVector::Int");
+ FUNCNAME("DOFVector::Int()");
Mesh* mesh = feSpace->getMesh();
- int deg;
- dim = mesh->getDim();
-
- if (!q)
- {
- deg = 2*feSpace->getBasisFcts()->getDegree();
- q = Quadrature::provideQuadrature(dim, deg);
- }
+ if (!q) {
+ int deg = 2 * feSpace->getBasisFcts()->getDegree();
+ q = Quadrature::provideQuadrature(this->dim, deg);
+ }
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *q, INIT_PHI);
-
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
@@ -591,28 +590,20 @@ namespace AMDiS {
return 0;
}
- // ... und die L1-Norm ...
-
template<typename T>
double DOFVector<T>::L1Norm(Quadrature* q) const
{
- FUNCNAME("DOFVector::L1Norm");
+ FUNCNAME("DOFVector::L1Norm()");
- Mesh* mesh = feSpace->getMesh();
-
- int deg;
- dim = mesh->getDim();
-
- if (!q)
- {
- deg = 2*feSpace->getBasisFcts()->getDegree();
- q = Quadrature::provideQuadrature(dim, deg);
- }
+ if (!q) {
+ int deg = 2 * feSpace->getBasisFcts()->getDegree();
+ q = Quadrature::provideQuadrature(this->dim, deg);
+ }
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);
-
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
+ Mesh* mesh = feSpace->getMesh();
mesh->traverse(-1,
Mesh::CALL_LEAF_EL|
@@ -643,23 +634,15 @@ namespace AMDiS {
}
- // bis hierhin gehen Andreas Ergaenzungen...
-
template<typename T>
double DOFVector<T>::L2NormSquare(Quadrature* q) const
{
- FUNCNAME("DOFVector::L2NormSquare");
-
- Mesh* mesh = feSpace->getMesh();
+ FUNCNAME("DOFVector::L2NormSquare()");
- int deg;
- dim = mesh->getDim();
-
- if (!q)
- {
- deg = 2*feSpace->getBasisFcts()->getDegree();
- q = Quadrature::provideQuadrature(dim, deg);
- }
+ if (!q) {
+ int deg = 2 * feSpace->getBasisFcts()->getDegree();
+ q = Quadrature::provideQuadrature(this->dim, deg);
+ }
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q,
@@ -667,8 +650,8 @@ namespace AMDiS {
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
-
- mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, L2Norm_fct);
+ Mesh* mesh = feSpace->getMesh();
+ mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
return norm;
}
@@ -722,24 +705,19 @@ namespace AMDiS {
template<typename T>
double DOFVector<T>::H1NormSquare(Quadrature *q) const
{
- FUNCNAME("DOFVector::H1NormSquare");
- int deg;
-
- Mesh *mesh = feSpace->getMesh();
- dim = mesh->getDim();
+ FUNCNAME("DOFVector::H1NormSquare()");
- if (!q)
- {
- deg = 2*feSpace->getBasisFcts()->getDegree()-2;
- q = Quadrature::provideQuadrature(dim, deg);
- }
- quad_fast =
- FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
- *q,
- INIT_GRD_PHI);
+ if (!q) {
+ int deg = 2 * feSpace->getBasisFcts()->getDegree() - 2;
+ q = Quadrature::provideQuadrature(this->dim, deg);
+ }
+
+ quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
+ *q, INIT_GRD_PHI);
norm = 0.0;
- traverseVector = const_cast<DOFVector<T>*>(this);
+ traverseVector = const_cast<DOFVector<T>*>(this);
+ Mesh *mesh = feSpace->getMesh();
mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS |
Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
@@ -752,10 +730,9 @@ namespace AMDiS {
void DOFVector<T>::compressDOFIndexed(int first, int last,
std::vector<DegreeOfFreedom> &newDOF)
{
- int i, j;
- for(i = first; i <= last; i++) {
- if((j = newDOF[i]) >= 0) {
- vec[j] = vec[i];
+ for (int i = first; i <= last; i++) {
+ if (newDOF[i] >= 0) {
+ vec[newDOF[i]] = vec[i];
}
}
}
@@ -1252,15 +1229,9 @@ namespace AMDiS {
TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
("invalid basis functions");
- Element *el = elInfo->getElement();
-
- const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
const BasisFunction *basFcts = feSpace->getBasisFcts();
-
- int numPoints = quadrature->getNumPoints();
-
+ int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();
static T *localvec = NULL;
-
T *result;
if (vecAtQPs) {
@@ -1268,17 +1239,17 @@ namespace AMDiS {
} else {
if (localvec)
delete [] localvec;
- localvec = new T[numPoints];
- for (int i = 0; i < numPoints; i++) {
+ localvec = new T[nPoints];
+ for (int i = 0; i < nPoints; i++) {
localvec[i] = 0.0;
}
result = localvec;
}
- T *localVec = new T[nBasFcts];
- getLocalVector(el, localVec);
+ T *localVec = localVectors[omp_get_thread_num()];
+ getLocalVector(elInfo->getElement(), localVec);
- for (int i = 0; i < numPoints; i++) {
+ for (int i = 0; i < nPoints; i++) {
result[i] = 0.0;
for (int j = 0; j < nBasFcts; j++) {
result[i] += localVec[j] *
@@ -1287,9 +1258,7 @@ namespace AMDiS {
((*(basFcts->getPhi(j)))(quad->getLambda(i))));
}
}
-
- delete [] localVec;
-
+
return const_cast<const T*>(result);
}
@@ -1299,23 +1268,18 @@ namespace AMDiS {
{
FUNCNAME("DOFVector::DoubleWell()");
- Mesh* mesh = feSpace->getMesh();
-
- int deg = 0;
- dim = mesh->getDim();
-
if (!q) {
- deg = 2 * feSpace->getBasisFcts()->getDegree();
- q = Quadrature::provideQuadrature(dim, deg);
+ int deg = 2 * feSpace->getBasisFcts()->getDegree();
+ q = Quadrature::provideQuadrature(this->dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q, INIT_PHI);
-
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
-
- mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, DoubleWell_fct);
+ Mesh* mesh = feSpace->getMesh();
+ mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET,
+ DoubleWell_fct);
return norm;
}
@@ -1323,15 +1287,13 @@ namespace AMDiS {
template<typename T>
int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
{
- double det = elinfo->getDet();
- const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL,
- quad_fast, NULL);
- int numPoints = quad_fast->getNumPoints();
+ const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
+ int nPoints = quad_fast->getNumPoints();
double normT = 0.0;
- for (int iq = 0; iq < numPoints; iq++) {
- normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq])*sqr(1.0 - uh_vec[iq]);
+ for (int iq = 0; iq < nPoints; iq++) {
+ normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);
}
- norm += det * normT;
+ norm += elinfo->getDet() * normT;
return 0;
}
diff --git a/AMDiS/src/ElInfo3d.cc b/AMDiS/src/ElInfo3d.cc
index 4b8589ab..30b28c54 100644
--- a/AMDiS/src/ElInfo3d.cc
+++ b/AMDiS/src/ElInfo3d.cc
@@ -609,4 +609,63 @@ namespace AMDiS {
}
}
+ void ElInfo3d::getRefSimplexCoords(DimMat<double> *coords) const
+ {
+ (*coords)[0][0] = 1.0;
+ (*coords)[1][0] = 0.0;
+ (*coords)[2][0] = 0.0;
+ (*coords)[3][0] = 0.0;
+
+ (*coords)[0][1] = 0.0;
+ (*coords)[1][1] = 1.0;
+ (*coords)[2][1] = 0.0;
+ (*coords)[3][1] = 0.0;
+
+ (*coords)[0][2] = 0.0;
+ (*coords)[1][2] = 0.0;
+ (*coords)[2][2] = 1.0;
+ (*coords)[3][2] = 0.0;
+
+ (*coords)[0][3] = 0.0;
+ (*coords)[1][3] = 0.0;
+ (*coords)[2][3] = 0.0;
+ (*coords)[3][3] = 1.0;
+ }
+
+ void ElInfo3d::getSubElementCoords(DimMat<double> *coords,
+ int iChild) const
+ {
+ FUNCNAME("ElInfo3d::getSubElementCoords()");
+
+ ERROR_EXIT("Not yet\n");
+
+ 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;
+ double c3 = ((*coords)[3][0] + (*coords)[3][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][3] = c0;
+ (*coords)[1][3] = c1;
+ (*coords)[2][3] = c2;
+ (*coords)[3][3] = c3;
+ }
+
}
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index af403f2a..adc7931c 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/AMDiS/src/ElInfo3d.h
@@ -119,20 +119,9 @@ namespace AMDiS {
orientation = o;
}
- void getRefSimplexCoords(DimMat<double> *coords) const
- {
- FUNCNAME("ElInfo3d::getRefSimplexCoords()");
-
- ERROR_EXIT("NOT YET\n");
- }
+ void getRefSimplexCoords(DimMat<double> *coords) const;
- void getSubElementCoords(DimMat<double> *coords,
- int iChild) const
- {
- FUNCNAME("ElInfo3d::getSubElementCoords()");
-
- ERROR_EXIT("NOT YET\n");
- }
+ void getSubElementCoords(DimMat<double> *coords, int iChild) const;
protected:
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index 06fc4606..b1841b8f 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -425,8 +425,6 @@ namespace AMDiS {
// create finite element space
int degree = 1;
GET_PARAMETER(0, name_ + "->polynomial degree" ,"%d", °ree);
- std::cout << degree << std::endl;
- WAIT_REALLY;
feSpace_ = FiniteElemSpace::provideFESpace(NULL,
Lagrange::getLagrange(mesh->getDim(), degree),
mesh,
diff --git a/AMDiS/src/SolutionDataStorage.h b/AMDiS/src/SolutionDataStorage.h
index d54309c6..bbf67af0 100644
--- a/AMDiS/src/SolutionDataStorage.h
+++ b/AMDiS/src/SolutionDataStorage.h
@@ -3,6 +3,7 @@
#include <iostream>
#include <vector>
+#include <map>
#include "DOFVector.h"
#include "SystemVector.h"
#include "Parameters.h"
@@ -21,6 +22,12 @@ namespace AMDiS {
typedef std::vector<FiniteElemSpace*> type;
};
+
+ class DataContainer
+ {
+ public:
+ };
+
template<typename T>
class SolutionDataStorage
{
@@ -85,6 +92,16 @@ namespace AMDiS {
return poped;
}
+ void addContainer(std::string name);
+
+ void push(std::string name, WorldVector<double> value);
+
+ void pop(std::string name, WorldVector<double> &value);
+
+ void reset(std::string name);
+
+ void clear(std::string name);
+
protected:
/** \brief
* Deletes the fe space and all it content, i.e., also the dof admin,
@@ -172,6 +189,10 @@ namespace AMDiS {
* Counts the memory usage of all solutions stored in memory.
*/
int memoryUsage;
+
+ std::map<std::string, std::vector<WorldVector<double > > > containers;
+
+ std::map<std::string, int> containersPos;
};
}
diff --git a/AMDiS/src/SolutionDataStorage.hh b/AMDiS/src/SolutionDataStorage.hh
index 41586dfe..09fc5bbc 100644
--- a/AMDiS/src/SolutionDataStorage.hh
+++ b/AMDiS/src/SolutionDataStorage.hh
@@ -52,8 +52,6 @@ namespace AMDiS {
solutions.push_back(solution);
timestamps.push_back(timestamp);
- std::cout << "STACK SIZE = " << solutions.size() << std::endl;
-
lastPos++;
// memoryUsage += solution->calcMemoryUsage();
}
@@ -128,7 +126,7 @@ namespace AMDiS {
feSpaces.clear();
timestamps.clear();
- lastPos = -1;
+ lastPos = -1;
}
template<typename T>
@@ -140,4 +138,39 @@ namespace AMDiS {
}
}
+ template<typename T>
+ void SolutionDataStorage<T>::addContainer(std::string name)
+ {
+ std::vector<WorldVector<double> > emptyVec;
+ containers.insert(std::pair<std::string, std::vector<WorldVector<double> > >(name, emptyVec) );
+ containersPos.insert(std::pair<std::string, int>(name, -1));
+ }
+
+ template<typename T>
+ void SolutionDataStorage<T>::push(std::string name, WorldVector<double> value)
+ {
+ containers[name].push_back(value);
+ containersPos[name] = containers[name].size() - 1;
+ }
+
+ template<typename T>
+ void SolutionDataStorage<T>::pop(std::string name, WorldVector<double> &value)
+ {
+ value = containers[name][containersPos[name]];
+ containersPos[name]--;
+ }
+
+ template<typename T>
+ void SolutionDataStorage<T>::reset(std::string name)
+ {
+ containersPos[name] = containers[name].size() - 1;
+ }
+
+ template<typename T>
+ void SolutionDataStorage<T>::clear(std::string name)
+ {
+ containers[name].clear();
+ containersPos[name] = -1;
+ }
+
}
--
GitLab