Commit e22abd25 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* This and that

parent 1500f6ae
......@@ -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++) {
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++) {
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++) {
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++) {
} 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);
}
......@@ -1602,41 +1511,27 @@ 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));
if (rowFESpace == colFESpace) {
elMat->colIndices = elMat->rowIndices;
} else {
colFESpace->getBasisFcts()->getLocalIndicesVec(element,
colAdmin