From e22abd25aab0a1ad2c3d41989b61196fc89dfc54 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 7 May 2008 10:24:03 +0000
Subject: [PATCH] * This and that
---
AMDiS/src/Assembler.cc | 369 ++++++++++++++------------------------
AMDiS/src/Assembler.h | 2 +-
AMDiS/src/MemoryManager.h | 9 +-
AMDiS/src/Operator.cc | 122 ++-----------
AMDiS/src/Quadrature.h | 67 +++----
5 files changed, 174 insertions(+), 395 deletions(-)
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 373fa100..a3fab99f 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -60,7 +60,7 @@ namespace AMDiS {
// check if all terms are symmetric
symmetric = true;
- for (int i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
if (!(terms[i])->isSymmetric()) {
symmetric = false;
break;
@@ -181,29 +181,20 @@ namespace AMDiS {
// calculate new values
const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();
- //double* uhLoc = GET_MEMORY(double, basFcts->getNumber());
- //vec->getLocalVector(elInfo->getElement(), uhLoc);
-
if (opt && !quad && sameFESpaces) {
if (psiFast->getBasisFunctions() == basFcts) {
- //psiFast->uhAtQp(uhLoc, values);
vec->getVecAtQPs(elInfo, NULL, psiFast, values);
} else if(phiFast->getBasisFunctions() == basFcts) {
- //phiFast->uhAtQp(uhLoc, values);
vec->getVecAtQPs(elInfo, NULL, phiFast, values);
} else {
vec->getVecAtQPs(elInfo, localQuad, NULL, values);
}
} else {
- //localQuad->uhAtQp(basFcts, uhLoc, values);
vec->getVecAtQPs(elInfo, localQuad, NULL, values);
}
-
- //FREE_MEMORY(uhLoc, double, basFcts->getNumber());
-
+
valuesAtQPs[vec]->valid = true;
- // return values
return values;
}
@@ -518,9 +509,9 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
- (*mat)[i][i] += quadrature->getWeight(iq)*c[iq]*psival*phival[i];
+ (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
for (int j = i + 1; j < nCol; j++) {
- val = quadrature->getWeight(iq)*c[iq]*psival*phival[j];
+ val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
@@ -544,8 +535,8 @@ namespace AMDiS {
}
}
- FREE_MEMORY(c, double, numPoints);
FREE_MEMORY(phival, double, nCol);
+ FREE_MEMORY(c, double, numPoints);
}
void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
@@ -858,9 +849,10 @@ namespace AMDiS {
for (int j = 0; j < nCol; j++) {
k = q10->getKVec(i, j);
values = q10->getValVec(i, j);
- int m;
- for (val = m = 0; m < nEntries[i][j]; m++)
- val += values[m]*Lb[0][k[m]];
+ val = 0.0;
+ for (int m = 0; m < nEntries[i][j]; m++) {
+ val += values[m] * Lb[0][k[m]];
+ }
(*mat)[i][j] += val;
}
}
@@ -1045,9 +1037,10 @@ namespace AMDiS {
for (int j = 0; j < nCol; j++) {
l = q01->getLVec(i, j);
values = q01->getValVec(i, j);
- int m;
- for (val = m = 0; m < nEntries[i][j]; m++)
- val += values[m]*Lb[0][l[m]];
+ val = 0.0;
+ for (int m = 0; m < nEntries[i][j]; m++) {
+ val += values[m] * Lb[0][l[m]];
+ }
(*mat)[i][j] += val;
}
}
@@ -1059,10 +1052,9 @@ namespace AMDiS {
const int *k;
const double *values;
- int i, m;
double val;
- if(firstCall) {
+ if (firstCall) {
q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -1074,30 +1066,26 @@ namespace AMDiS {
const int *nEntries = q1->getNumberEntries();
Lb[0].set(0.0);
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<FirstOrderTerm*>(terms[i]))->getLb(elInfo, 1, Lb);
}
Lb[0] *= elInfo->getDet();
- for (i = 0; i < nRow; i++) {
- k = q1->getKVec(i);
+ for (int i = 0; i < nRow; i++) {
+ k = q1->getKVec(i);
values = q1->getValVec(i);
- for (val = m = 0; m < nEntries[i]; m++)
- val += values[m]*Lb[0][k[m]];
+ val = 0.0;
+ for (int m = 0; m < nEntries[i]; m++) {
+ val += values[m] * Lb[0][k[m]];
+ }
(*vec)[i] += val;
}
-
- //DELETE [] Lb;
}
Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
- {
- // q11 = Q11PsiPhi::provideQ11PsiPhi(assembler->getRowFESpace()->getBasisFcts(),
- // assembler->getColFESpace()->getBasisFcts(),
- // quadrature);
- }
+ {}
void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
@@ -1106,10 +1094,9 @@ namespace AMDiS {
const int **nEntries;
const int *k, *l;
const double *values;
- int i, j, m;
double val;
- if(firstCall) {
+ if (firstCall) {
q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(),
owner->getColFESpace()->getBasisFcts(),
quadrature);
@@ -1118,7 +1105,7 @@ namespace AMDiS {
LALt[0]->set(0.0);
- for(i=0; i < static_cast<int>( terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>( terms.size()); i++) {
(static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, 1, LALt);
}
@@ -1127,32 +1114,38 @@ namespace AMDiS {
nEntries = q11->getNumberEntries();
if (symmetric) {
- for (i = 0; i < nRow; i++) {
- k = q11->getKVec(i, i);
- l = q11->getLVec(i, i);
+ for (int i = 0; i < nRow; i++) {
+ k = q11->getKVec(i, i);
+ l = q11->getLVec(i, i);
values = q11->getValVec(i, i);
- for (val = m = 0; m < nEntries[i][i]; m++)
- val += values[m]*(*LALt[0])[k[m]][l[m]];
+ val = 0.0;
+ for (int m = 0; m < nEntries[i][i]; m++) {
+ val += values[m] * (*LALt[0])[k[m]][l[m]];
+ }
(*mat)[i][i] += val;
- for (j = i+1; j < nCol; j++) {
- k = q11->getKVec(i, j);
- l = q11->getLVec(i, j);
+
+ for (int j = i + 1; j < nCol; j++) {
+ k = q11->getKVec(i, j);
+ l = q11->getLVec(i, j);
values = q11->getValVec(i, j);
- for (val = m = 0; m < nEntries[i][j]; m++)
- val += values[m]*(*LALt[0])[k[m]][l[m]];
+ val = 0.0;
+ for (int m = 0; m < nEntries[i][j]; m++) {
+ val += values[m] * (*LALt[0])[k[m]][l[m]];
+ }
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
- }
- else { /* A not symmetric or psi != phi */
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++) {
- k = q11->getKVec(i, j);
- l = q11->getLVec(i, j);
+ } else { /* A not symmetric or psi != phi */
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
+ k = q11->getKVec(i, j);
+ l = q11->getLVec(i, j);
values = q11->getValVec(i, j);
- for (val = m = 0; m < nEntries[i][j]; m++)
- val += values[m]*(*LALt[0])[k[m]][l[m]];
+ val = 0.0;
+ for (int m = 0; m < nEntries[i][j]; m++) {
+ val += values[m] * (*LALt[0])[k[m]][l[m]];
+ }
(*mat)[i][j] += val;
}
}
@@ -1162,67 +1155,16 @@ namespace AMDiS {
DELETE LALt;
}
- // void Pre2::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // FUNCNAME("Pre2::calculateElementVector");
- // ERROR_EXIT("should not be called\n");
- // }
-
- // void Pre2::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // DimMat<double> LALt(dim, NO_INIT);
- // const int *nEntries;
- // const int *k, *l;
- // const double *values;
- // int i, m;
- // double val;
-
- // LALt.set(0.0);
-
- // for(i=0; i < static_cast<int>( terms->size()); i++) {
- // (static_cast<SecondOrderTerm*>((*terms)[i])->eval(elInfo, 0, &LALt);
- // }
-
- // nEntries = q2->getNumberEntries();
-
- // for (i = 0; i < nRow; i++) {
- // k = q2->getKVec(i);
- // l = q2->getLVec(i);
- // values = q2->getValVec(i);
- // for (val = m = 0; m < nEntries[i]; m++)
- // val += values[m]*LALt[k[m]][l[m]];
- // vec[i] += val;
- // }
- // }
-
Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, true)
- {
- // if(!psiFast) {
- // psiFast = FastQuadrature::provideFastQuadrature(
- // assembler->getRowFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_GRD_PHI);
- // } else {
- // psiFast->init(INIT_GRD_PHI);
- // }
- // if(!phiFast) {
- // phiFast = FastQuadrature::provideFastQuadrature(
- // assembler->getColFESpace()->getBasisFcts(),
- // *quadrature,
- // INIT_PHI);
- // } else {
- // phiFast->init(INIT_GRD_PHI);
- // }
- }
+ {}
void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
double val;
VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;
- int iq, i, j;
- if(firstCall) {
+ if (firstCall) {
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
basFcts = owner->getColFESpace()->getBasisFcts();
@@ -1233,26 +1175,28 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
- for(i=0;i<nPoints;i++) LALt[i]=NEW DimMat<double>(dim, NO_INIT);
- for (iq = 0; iq < nPoints; iq++) {
+ for (int i = 0; i < nPoints;i++) {
+ LALt[i] = NEW DimMat<double>(dim, NO_INIT);
+ }
+ for (int iq = 0; iq < nPoints; iq++) {
LALt[iq]->set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, nPoints, LALt);
}
if (symmetric) {
- for (iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
grdPhi = phiFast->getGradient(iq);
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
(*mat)[i][i] += quadrature->getWeight(iq) *
((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));
- for (j = i+1; j < nCol; j++) {
+ for (int j = i + 1; j < nCol; j++) {
val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
(*mat)[i][j] += val;
(*mat)[j][i] += val;
@@ -1261,14 +1205,14 @@ namespace AMDiS {
}
}
else { /* non symmetric assembling */
- for (iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
grdPsi = psiFast->getGradient(iq);
grdPhi = phiFast->getGradient(iq);
- for (i = 0; i < nRow; i++) {
- for (j = 0; j < nCol; j++) {
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq) *
((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
}
@@ -1276,16 +1220,12 @@ namespace AMDiS {
}
}
- for(i=0;i<nPoints;i++) DELETE LALt[i];
+ for (int i = 0;i < nPoints; i++)
+ DELETE LALt[i];
+
DELETE [] LALt;
}
- // void Quad2::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // FUNCNAME("Quad2::calculateElementVector");
- // ERROR_EXIT("should not be called\n");
- // }
-
Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad)
: SecondOrderAssembler(op, assembler, quad, false)
{}
@@ -1295,7 +1235,6 @@ namespace AMDiS {
double val;
DimVec<double> grdPsi(dim, NO_INIT);
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
- int iq, i, j;
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
@@ -1303,46 +1242,46 @@ namespace AMDiS {
int nPoints = quadrature->getNumPoints();
DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
- for (iq = 0; iq < nPoints; iq++) {
- LALt[iq]=NEW DimMat<double>(dim,NO_INIT);
+ for (int iq = 0; iq < nPoints; iq++) {
+ LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
LALt[iq]->set(0.0);
}
- for(i=0; i < static_cast<int>(terms.size()); i++) {
+
+ for (int i = 0; i < static_cast<int>(terms.size()); i++) {
(static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, nPoints, LALt);
}
if (symmetric) {
- for (iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
- for(i=0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++) {
grdPhi[i] = (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq));
}
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
grdPsi = (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq));
(*mat)[i][i] += quadrature->getWeight(iq) *
(grdPsi * ((*LALt[iq]) * grdPhi[i]));
- for (j = i+1; j < nCol; j++) {
+ for (int j = i + 1; j < nCol; j++) {
val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j]));
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
}
}
- }
- else { /* non symmetric assembling */
- for (iq = 0; iq < nPoints; iq++) {
+ } else { /* non symmetric assembling */
+ for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();
- for(i=0; i < nCol; i++) {
+ for (int i = 0; i < nCol; i++) {
grdPhi[i] = (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq));
}
- for (i = 0; i < nRow; i++) {
+ for (int i = 0; i < nRow; i++) {
grdPsi = (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq));
- for (j = 0; j < nCol; j++) {
+ for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq) *
(grdPsi * ((*LALt[iq]) * grdPhi[j]));
}
@@ -1350,16 +1289,12 @@ namespace AMDiS {
}
}
- for(iq=0;iq<nPoints;iq++) DELETE LALt[iq];
+ for (int iq = 0; iq < nPoints; iq++)
+ DELETE LALt[iq];
+
DELETE [] LALt;
}
- // void Stand2::calculateElementVector(const ElInfo *elInfo, double *vec)
- // {
- // FUNCNAME("Stand2::calculateElementVector");
- // ERROR_EXIT("should not be called\n");
- // }
-
Assembler::Assembler(Operator *op,
const FiniteElemSpace *rowFESpace_,
const FiniteElemSpace *colFESpace_)
@@ -1377,9 +1312,7 @@ namespace AMDiS {
lastVecEl(NULL),
lastTraverseId(-1)
- {
- //if(op->uhOld) rememberElMat = true;
- }
+ {}
void Assembler::calculateElementMatrix(const ElInfo *elInfo,
ElementMatrix *userMat,
@@ -1439,7 +1372,7 @@ namespace AMDiS {
{
FUNCNAME("Assembler::calculateElementVector()");
- if(remember && factor != 1.0) {
+ if (remember && factor != 1.0) {
rememberElVec = true;
}
@@ -1452,100 +1385,76 @@ namespace AMDiS {
checkQuadratures();
- if((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
+ if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
initElement(elInfo);
}
- if(el != lastVecEl || !operat->isOptimized()) {
- if(rememberElVec) {
+ if (el != lastVecEl || !operat->isOptimized()) {
+ if (rememberElVec) {
elementVector->set(0.0);
}
lastVecEl = el;
} else {
- if(rememberElVec) {
+ if (rememberElVec) {
axpy(factor, *elementVector, *userVec);
- //*userVec += *elementVector * factor;
- //operat->addElementVector(elementVector, userVec, factor);
return;
}
}
ElementVector *vec = rememberElVec ? elementVector : userVec;
- if(operat->uhOld && remember) {
+ if (operat->uhOld && remember) {
matVecAssemble(elInfo, vec);
- if(rememberElVec) {
+ if (rememberElVec) {
axpy(factor, *elementVector, *userVec);
- //*userVec += *elementVector * factor;
- //operat->addElementVector(elementVector, userVec, factor);
}
return;
}
- //if(secondOrderAssembler)
- // secondOrderAssembler->calculateElementVector(elInfo, vec);
- if(firstOrderAssemblerGrdPsi)
+ if (firstOrderAssemblerGrdPsi)
firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
- //if(firstOrderAssemblerGrdPhi)
- // firstOrderAssemblerGrdPhi->calculateElementVector(elInfo, vec);
- if(zeroOrderAssembler)
+ if (zeroOrderAssembler)
zeroOrderAssembler->calculateElementVector(elInfo, vec);
- if(rememberElVec) {
+ if (rememberElVec) {
axpy(factor, *elementVector, *userVec);
- //*userVec += *elementVector * factor;
- //operat->addElementVector(elementVector, userVec, factor);
}
-
- // MSG("\n");
- // for(int i=0; i < 3; i++) {
- // MSG("%e\n", (*vec)[i]);
- // }
-
}
void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
{
FUNCNAME("Assembler::matVecAssemble()");
- int i, j;
-
Element *el = elInfo->getElement();
const BasisFunction *basFcts = rowFESpace->getBasisFcts();
const double *uhOldLoc = operat->uhOld->getLocalVector(el, NULL);
int n = basFcts->getNumber();
- if(el != lastMatEl) {
+ if (el != lastMatEl) {
calculateElementMatrix(elInfo, NULL);
}
double val;
- for (i = 0; i < n; i++) {
- val = 0;
- for (j = 0; j < n; j++) {
+ for (int i = 0; i < n; i++) {
+ val = 0.0;
+ for (int j = 0; j < n; j++) {
val += (*elementMatrix)[i][j]*uhOldLoc[j];
}
(*vec)[i] += val;
}
-
- // MSG("\n");
- // MSG("uh loc:\n");
- // for (i = 0; i < n; i++) {
- // MSG("%e\n", uhOldLoc[i]);
- // }
}
void Assembler::initElement(const ElInfo *elInfo, Quadrature *quad)
{
checkQuadratures();
- if(secondOrderAssembler)
+ if (secondOrderAssembler)
secondOrderAssembler->initElement(elInfo, quad);
- if(firstOrderAssemblerGrdPsi)
+ if (firstOrderAssemblerGrdPsi)
firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
- if(firstOrderAssemblerGrdPhi)
+ if (firstOrderAssemblerGrdPhi)
firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
- if(zeroOrderAssembler)
+ if (zeroOrderAssembler)
zeroOrderAssembler->initElement(elInfo, quad);
}
@@ -1561,13 +1470,13 @@ namespace AMDiS {
bool opt = (rowFESpace_ == colFESpace_);
// create sub assemblers
- secondOrderAssembler =
+ secondOrderAssembler =
SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
- firstOrderAssemblerGrdPsi =
+ firstOrderAssemblerGrdPsi =
FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
- firstOrderAssemblerGrdPhi =
+ firstOrderAssemblerGrdPhi =
FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
- zeroOrderAssembler =
+ zeroOrderAssembler =
ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
}
@@ -1583,13 +1492,13 @@ namespace AMDiS {
remember = false;
// create sub assemblers
- secondOrderAssembler =
+ secondOrderAssembler =
SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
- firstOrderAssemblerGrdPsi =
+ firstOrderAssemblerGrdPsi =
FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
- firstOrderAssemblerGrdPhi =
+ firstOrderAssemblerGrdPhi =
FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
- zeroOrderAssembler =
+ zeroOrderAssembler =
ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
}
@@ -1601,42 +1510,28 @@ namespace AMDiS {
}
elMat->set(0.0);
-
- DOFAdmin *rowAdmin = rowFESpace->getAdmin();
- DOFAdmin *colAdmin = colFESpace->getAdmin();
-
-
+
Element *element = elInfo->getElement();
- /*
- elMat->rowIndices =
- rowFESpace->getBasisFcts()->getLocalIndices(element,
- rowAdmin,
- NULL);
- */
-
- /*
- elMat->colIndices =
- colFESpace->getBasisFcts()->getLocalIndices(element,
- colAdmin,
- NULL);
- */
-
-
rowFESpace->getBasisFcts()->getLocalIndicesVec(element,
- rowAdmin,
+ rowFESpace->getAdmin(),
&(elMat->rowIndices));
- colFESpace->getBasisFcts()->getLocalIndicesVec(element,
- colAdmin,
- &(elMat->colIndices));
+ if (rowFESpace == colFESpace) {
+ elMat->colIndices = elMat->rowIndices;
+ } else {
+ colFESpace->getBasisFcts()->getLocalIndicesVec(element,
+ colFESpace->getAdmin(),
+ &(elMat->colIndices));
+ }
+
return elMat;
}
ElementVector *Assembler::initElementVector(ElementVector *elVec,
const ElInfo *elInfo)
{
- if(!elVec) {
+ if (!elVec) {
elVec = NEW ElementVector(nRow);
}
@@ -1655,40 +1550,38 @@ namespace AMDiS {
}
void Assembler::checkQuadratures()
- {
- int dim;
-
- if(secondOrderAssembler) {
+ {
+ if (secondOrderAssembler) {
// create quadrature
- if(!secondOrderAssembler->getQuadrature()) {
- dim = rowFESpace->getMesh()->getDim();
+ if (!secondOrderAssembler->getQuadrature()) {
+ int dim = rowFESpace->getMesh()->getDim();
int degree = operat->getQuadratureDegree(2);
Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
secondOrderAssembler->setQuadrature(quadrature);
}
}
- if(firstOrderAssemblerGrdPsi) {
+ if (firstOrderAssemblerGrdPsi) {
// create quadrature
- if(!firstOrderAssemblerGrdPsi->getQuadrature()) {
- dim = rowFESpace->getMesh()->getDim();
+ if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
+ int dim = rowFESpace->getMesh()->getDim();
int degree = operat->getQuadratureDegree(1, GRD_PSI);
Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
}
}
- if(firstOrderAssemblerGrdPhi) {
+ if (firstOrderAssemblerGrdPhi) {
// create quadrature
- if(!firstOrderAssemblerGrdPhi->getQuadrature()) {
- dim = rowFESpace->getMesh()->getDim();
+ if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
+ int dim = rowFESpace->getMesh()->getDim();
int degree = operat->getQuadratureDegree(1, GRD_PHI);
Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
}
}
- if(zeroOrderAssembler) {
+ if (zeroOrderAssembler) {
// create quadrature
- if(!zeroOrderAssembler->getQuadrature()) {
- dim = rowFESpace->getMesh()->getDim();
+ if (!zeroOrderAssembler->getQuadrature()) {
+ int dim = rowFESpace->getMesh()->getDim();
int degree = operat->getQuadratureDegree(0);
Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
zeroOrderAssembler->setQuadrature(quadrature);
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index ca874eb4..c3cf09d4 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -70,7 +70,7 @@ namespace AMDiS {
* \brief
* Base class for SecondOrderAssembler, FirstOrderAssembler,
* ZeroOrderAssembler. The task of a SubAssembler is to assemble a list of
- * terms of a spezial order and add their contributions to a DOFMatrix or a
+ * terms of a special order and add their contributions to a DOFMatrix or a
* DOFVector. An Assembler can consist of up to four SubAssemblers: one
* SecondOrderAssembler for second order terms, one ZeroOrderAssembler for
* terms of order zero, and two FirstOrderAssemblers. One for terms with
diff --git a/AMDiS/src/MemoryManager.h b/AMDiS/src/MemoryManager.h
index d6daea39..c2045de6 100644
--- a/AMDiS/src/MemoryManager.h
+++ b/AMDiS/src/MemoryManager.h
@@ -251,13 +251,16 @@ namespace AMDiS {
*/
static T* getMemory(size_t s, const char* filename, int line) {
FUNCNAME("MemoryManager<T>::getMemory()");
- if(!singleton) init();
+ if (!singleton)
+ init();
byteCount += s;
singleton->instanceCount += s/sizeof(T);
singleton->typedByteCount += s;
- if(printInfo) {
- if(filename && line)
+ if (printInfo) {
+ if (filename && line) {
MSG("FILE: %s, LINE: %d - ", filename, line);
+ }
+
Msg::print("%s::getMemory : %d instances, total : %d instances\n",
singleton->getName().c_str(), s/sizeof(T),
static_cast<int>(singleton->instanceCount));
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 2d1eaed2..2ac46ec0 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -239,8 +239,8 @@ namespace AMDiS {
Quadrature *quad1GrdPhi,
Quadrature *quad0)
{
- if(!assembler)
- if(optimized) {
+ if (!assembler)
+ if (optimized) {
assembler =
NEW OptimizedAssembler(this,
quad2, quad1GrdPsi, quad1GrdPhi, quad0,
@@ -597,13 +597,12 @@ namespace AMDiS {
}
void VecAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- if(f) {
- for(iq = 0; iq < numPoints; iq++) {
+ if (f) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq]);
}
} else {
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += vecAtQPs[iq];
}
}
@@ -625,8 +624,7 @@ namespace AMDiS {
// Andreas ergaenzt
void MultVecAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
}
}
@@ -649,8 +647,7 @@ namespace AMDiS {
}
void Vec2AtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
}
}
@@ -672,8 +669,7 @@ namespace AMDiS {
}
void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
}
}
@@ -696,8 +692,7 @@ namespace AMDiS {
void FctGradientCoords_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
}
}
@@ -719,8 +714,7 @@ namespace AMDiS {
}
void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
}
}
@@ -744,8 +738,7 @@ namespace AMDiS {
// bis hierhin
void VecAndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
}
}
@@ -768,8 +761,7 @@ namespace AMDiS {
void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const
{
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
}
}
@@ -791,8 +783,7 @@ namespace AMDiS {
}
void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
}
}
@@ -814,8 +805,7 @@ namespace AMDiS {
}
void FctGradient_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
- int iq;
- for(iq = 0; iq < numPoints; iq++) {
+ for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(gradAtQPs[iq]);
}
}
@@ -1387,90 +1377,6 @@ namespace AMDiS {
return assembler;
}
-
- // Matrix<double> **OperatorVec::createElementMatrix()
- // {
- // int i;
- // Matrix<double> **newElementMatrix = GET_MEMORY(Matrix<double>*, nRow);
- // for(i=0; i < nRow; i++) {
- // newElementMatrix[i] = NEW Matrix<double>[nCol](vecSize, vecSize);
- // }
-
- // return newElementMatrix;
- // }
-
- // void OperatorVec::freeElementMatrix(Matrix<double> **elMat) {
- // int i;
- // for(i = 0; i < nRow; i++) {
- // DELETE [] elMat[i];
- // }
- // FREE_MEMORY(elMat, Matrix<double>*, nRow);
- // }
-
- // Vector<double> *OperatorVec::createElementVector()
- // {
- // return NEW Vector<double>[nRow](vecSize);
- // }
-
- // void OperatorVec::freeElementVector(Vector<double> *elVec)
- // {
- // DELETE [] elVec;
- // }
-
-
-
- // void OperatorVec::getElementMatrix(ElInfo *elInfo,
- // Matrix<double> **userMat,
- // double factor = 1.0)
- // {
- // int i, j, k, l;
- // Operator *op = NULL;
- // double **elMat = NULL;
-
- // for(i = 0; i < operators.getNumRows(); i++) {
- // for(j = 0; j < operators.getNumCols(); j++) {
- // op = operators[i][j];
- // if(op) {
- // if(!elMat) elMat = op->createElementMatrix();
- // op->resetElementMatrix(elMat);
- // op->getElementMatrix(elInfo, elMat, factor);
- // for(k = 0; k < nRow; k++) {
- // for(l = 0; l < nCol; l++) {
- // userMat[i][j][k][l] += elMat[k][l];
- // }
- // }
- // }
- // }
- // }
-
- // if(elMat) op->freeElementMatrix(elMat);
- // }
-
- // void OperatorVec::getElementVector(const ElInfo *elInfo,
- // Vector<double> *userVec,
- // double factor = 1.0)
- // {
- // int i, j, k;
- // Operator *op = NULL;
- // double *elVec = NULL;
-
- // for(i = 0; i < operators.getNumRows(); i++) {
- // for(j = 0; j < operators.getNumCols(); j++) {
- // op = operators[i][j];
- // if(op) {
- // if(!elVec) elVec = op->createElementVector();
- // op->resetElementVector(elVec);
- // op->getElementVector(elInfo, elVec, factor);
- // for(k = 0; k < nRow; k++) {
- // userVec[i][k] += elVec[k];
- // }
- // }
- // }
- // }
-
- // if(elVec) op->freeElementVector(elVec);
- // }
-
void CoordsAtQP_IJ_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index 33720912..32c17f62 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -114,32 +114,44 @@ namespace AMDiS {
/** \brief
* Returns \ref name
*/
- inline const ::std::string& getName() { return name; };
+ inline const ::std::string& getName() {
+ return name;
+ };
/** \brief
* Returns \ref n_points
*/
- inline int getNumPoints() const { return n_points;};
+ inline int getNumPoints() const {
+ return n_points;
+ };
/** \brief
* Returns \ref w[p]
*/
- inline double getWeight(int p) const {return w[p];};
+ inline double getWeight(int p) const {
+ return w[p];
+ };
/** \brief
* Returns \ref w.
*/
- inline double* getWeight() const { return w; };
+ inline double* getWeight() const {
+ return w;
+ };
/** \brief
* Returns \ref dim
*/
- inline int getDim() const { return dim; };
+ inline int getDim() const {
+ return dim;
+ };
/** \brief
* Returns \ref degree
*/
- inline int getDegree() const { return degree; };
+ inline int getDegree() const {
+ return degree;
+ };
/** \brief
* Returns a pointer to a vector storing the values of a doubled valued
@@ -170,7 +182,7 @@ namespace AMDiS {
* Returns \ref lambda[a][b] which is the b-th coordinate entry of the a-th
* quadrature point
*/
- inline double getLambda(int a,int b) const {
+ inline double getLambda(int a, int b) const {
return (lambda ? (*lambda)[a][b] : 0.0);
};
@@ -185,45 +197,10 @@ namespace AMDiS {
/** \brief
* Returns \ref lambda which is a VectorOfFixvecs<DimVec<double> >.
*/
- VectorOfFixVecs<DimVec<double> > *getLambda() const { return lambda; };
-
- /** \brief
- * The function returns a pointer ptr to a vector of length
- * \ref n_points storing the values of \f$ u_h \f$ at all
- * quadrature points, i.e.
- * \f[ ptr[l] = u_h(lambda[l]) \f] where
- * l = 0, . . . , n_points - 1. vec is an optional memory pointer
- */
- // const double *uhAtQp(const BasisFunction *basFcts,
- // const double *uh_loc, double *vec) const;
-
-
- /** \brief
- * The function returns a pointer ptr to a vector of length
- * \ref n_points WorldVectors storing \f$ \nabla u_h \f$ at all
- * quadrature points, i.e.
- * \f[ ptr[l][i] = u_{h,xi}(lambda[l]) \f]
- * where l = 0, ... , n_points - 1, and i = 0, ... ,
- * DIM_OF_WORLD - 1; vec is an optional memory pointer
- */
- // const WorldVector<double> *grdUhAtQp(const BasisFunction *basFcts,
- // const DimVec<WorldVector<double> >& Lambda,
- // const double *uh_loc,
- // WorldVector<double> *vec) const;
+ VectorOfFixVecs<DimVec<double> > *getLambda() const {
+ return lambda;
+ };
- /** \brief
- * The function returns a pointer ptr to a vector of length
- * \ref n_points of WorldMatrices storing D2uh at all quadrature
- * points, i.e.
- * \f[ ptr[l][i][j] = u_{h,x_ix_j}(lambda[l]) \f]
- * where l = 0, ... , n_points - 1, and i, j = 0, ... ,
- * DIM_OF_WORLD - 1; vec is an optional memory pointer
- */
- // const WorldMatrix<double> *D2UhAtQp(const BasisFunction *basFcts,
- // const DimVec<WorldVector<double> >& Lambda,
- // const double *uh_loc,
- // WorldMatrix<double> *vec) const;
-
public:
/** \brief
--
GitLab