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

* This and that

parent ecbb2779
......@@ -336,7 +336,7 @@ namespace AMDiS {
if (!optimized) {
newAssembler = NEW Stand0(op, assembler, quad);
} else {
if(pwConst) {
if (pwConst) {
newAssembler = NEW Pre0(op, assembler, quad);
} else {
newAssembler = NEW Quad0(op, assembler, quad);
......@@ -493,12 +493,9 @@ namespace AMDiS {
void Stand0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
double val;
const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
double psival;
double *phival = GET_MEMORY(double, nCol);
int nPoints = quadrature->getNumPoints();
double *c = GET_MEMORY(double, nPoints);
......@@ -523,10 +520,10 @@ namespace AMDiS {
}
for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
(*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];
double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
......@@ -542,7 +539,7 @@ namespace AMDiS {
}
for (int i = 0; i < nRow; i++) {
psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
for (int j = 0; j < nCol; j++) {
(*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
}
......@@ -1141,7 +1138,6 @@ namespace AMDiS {
}
tmpMat *= elInfo->getDet();
nEntries = q11->getNumberEntries();
if (symmetric) {
......@@ -1274,7 +1270,6 @@ namespace AMDiS {
void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
{
double val;
DimVec<double> grdPsi(dim, NO_INIT);
VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);
......@@ -1309,7 +1304,7 @@ namespace AMDiS {
(grdPsi * ((*LALt[iq]) * grdPhi[i]));
for (int j = i + 1; j < nCol; j++) {
val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j]));
double val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j]));
(*mat)[i][j] += val;
(*mat)[j][i] += val;
}
......@@ -1517,13 +1512,10 @@ namespace AMDiS {
// create sub assemblers
secondOrderAssembler =
SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
firstOrderAssemblerGrdPsi =
FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
firstOrderAssemblerGrdPhi =
FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
zeroOrderAssembler =
ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
}
......
......@@ -53,13 +53,9 @@ namespace AMDiS {
class Q0Psi;
class Q1Psi;
class Q2Psi;
// class Operator;
// class OperatorTerm;
template<typename T> class DOFVectorBase;
// enum FirstOrderType;
// ============================================================================
// ===== class SubAssembler ===================================================
// ============================================================================
......@@ -1121,17 +1117,17 @@ namespace AMDiS {
/** \brief
* SubAssembler for the first order terms (grdPsi)
*/
FirstOrderAssembler *firstOrderAssemblerGrdPsi;
FirstOrderAssembler *firstOrderAssemblerGrdPsi;
/** \brief
* SubAssembler for the first order terms (grdPhi)
*/
FirstOrderAssembler *firstOrderAssemblerGrdPhi;
FirstOrderAssembler *firstOrderAssemblerGrdPhi;
/** \brief
* SubAssembler for the zero order terms
*/
ZeroOrderAssembler *zeroOrderAssembler;
ZeroOrderAssembler *zeroOrderAssembler;
bool remember;
......@@ -1196,13 +1192,13 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
StandardAssembler(Operator *op,
StandardAssembler(Operator *op,
Quadrature *quad2,
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0,
const FiniteElemSpace *rowFESpace_,
const FiniteElemSpace *colFESpace_=NULL);
const FiniteElemSpace *rowFESpace,
const FiniteElemSpace *colFESpace = NULL);
};
// ============================================================================
......@@ -1223,13 +1219,13 @@ namespace AMDiS {
/** \brief
* Constructor.
*/
OptimizedAssembler(Operator *op,
OptimizedAssembler(Operator *op,
Quadrature *quad2,
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0,
const FiniteElemSpace *rowFESpace_,
const FiniteElemSpace *colFESpace_=NULL);
const FiniteElemSpace *rowFESpace,
const FiniteElemSpace *colFESpace = NULL);
};
}
......
......@@ -59,8 +59,7 @@ namespace AMDiS {
MacroElement *macroNeighbour = mel->getNeighbour(i);
if (macroNeighbour) {
neighbour_[i] = macroNeighbour->getElement();
neighbour_[i] = macroNeighbour->getElement();
Element *nb = const_cast<Element*>(neighbour_[i]);
int edgeNo = oppVertex_[i] = mel->getOppVertex(i);
......@@ -110,6 +109,19 @@ namespace AMDiS {
break;
default:
std::cout << "------------- Error --------------" << std::endl;
std::cout << " Element index = " << element_->getIndex() << "\n\n";
for (int j = 0; j < neighbours; j++) {
if (mel->getNeighbour(j)) {
std::cout << " Neighbour " << j << ": "
<< mel->getNeighbour(j)->getElement()->getIndex()
<< std::endl;
} else {
std::cout << " Neighbour " << j << ": not existing" << std::endl;
}
std::cout << " OppVertex " << j << ": " << static_cast<int>(mel->getOppVertex(j)) << std::endl;
std::cout << std::endl;
}
ERROR_EXIT("should not happen!\n");
break;
}
......
......@@ -261,7 +261,7 @@ namespace AMDiS {
void VecAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
Quadrature* quad)
{
vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
}
......@@ -450,8 +450,7 @@ namespace AMDiS {
void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
lalt(Lambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
}
}
......@@ -502,8 +501,7 @@ namespace AMDiS {
void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int numPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
lalt(Lambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
}
......@@ -511,17 +509,15 @@ namespace AMDiS {
void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
}
}
void VecAtQP_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
if(b)
for (int iq = 0; iq < numPoints; iq++) {
if (b)
lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
else
l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
......@@ -530,8 +526,7 @@ namespace AMDiS {
void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq]));
}
}
......@@ -540,22 +535,19 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >& Lb) const
{
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++)
{
lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq]));
}
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq]));
}
}
void VectorGradient_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
if(f) {
for(iq = 0; iq < numPoints; iq++) {
if (f) {
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
}
} else {
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
}
}
......@@ -563,16 +555,14 @@ namespace AMDiS {
void VectorFct_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
}
}
void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
}
}
......@@ -596,8 +586,7 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
}
}
......@@ -615,13 +604,9 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f1)(vecAtQPs1[iq]) *
(*f2)(vecAtQPs2[iq]) *
uhAtQP[iq];
for (int iq = 0; iq < numPoints; iq++) {
result[iq] += fac * (*f1)(vecAtQPs1[iq]) *
(*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
}
}
......@@ -638,11 +623,8 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs1[iq], vecAtQPs2[iq]) *
for (int iq = 0; iq < numPoints; iq++) {
result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) *
uhAtQP[iq];
}
}
......@@ -660,11 +642,8 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs[iq], coordsAtQPs[iq]) *
for (int iq = 0; iq < numPoints; iq++) {
result[iq] += fac * (*f)(vecAtQPs[iq], coordsAtQPs[iq]) *
uhAtQP[iq];
}
}
......@@ -683,11 +662,8 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(gradAtQPs[iq], coordsAtQPs[iq]) *
for (int iq = 0; iq < numPoints; iq++) {
result[iq] += fac * (*f)(gradAtQPs[iq], coordsAtQPs[iq]) *
uhAtQP[iq];
}
}
......@@ -705,8 +681,7 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]) *
......@@ -714,8 +689,6 @@ namespace AMDiS {
}
}
// bis hierhin
void VecAndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
......@@ -729,8 +702,7 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs[iq], gradAtQPs[iq]) *
......@@ -752,8 +724,7 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs[iq], gradAtQPs[iq]) *
......@@ -774,8 +745,7 @@ namespace AMDiS {
double *result,
double fac) const
{
int iq;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) *
......@@ -826,23 +796,22 @@ namespace AMDiS {
double *result,
double factor) const
{
int i, j, dow = Global::getGeo(WORLD);
int iq;
int dow = Global::getGeo(WORLD);
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
double resultQP = 0.0;
WorldMatrix<double> A = (*matrixFct)(vecAtQPs[iq]);
if(D2UhAtQP) {
for(i=0; i < dow; i++) {
for(j=0; j < dow; j++) {
if (D2UhAtQP) {
for (int i = 0; i < dow; i++) {
for (int j = 0; j < dow; j++) {
resultQP += A[i][j] * D2UhAtQP[iq][j][i];
}
}
}
if(grdUhAtQP) {
if (grdUhAtQP) {
resultQP += (*divFct)(A) * grdUhAtQP[iq];
}
......@@ -854,10 +823,9 @@ namespace AMDiS {
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
int iq;
if(grdUhAtQP) {
if (grdUhAtQP) {
WorldMatrix<double> A;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
A = (*matrixFct)(vecAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
......@@ -872,14 +840,13 @@ namespace AMDiS {
double *result,
double factor) const
{
int i, j, dow = Global::getGeo(WORLD);
int iq;
int dow = Global::getGeo(WORLD);
if(D2UhAtQP) {
for(iq = 0; iq < numPoints; iq++) {
if (D2UhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double resultQP = 0.0;
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++) {
resultQP += matrix[i][j] * D2UhAtQP[iq][j][i];
}
}
......@@ -892,9 +859,8 @@ namespace AMDiS {
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if(grdUhAtQP) {
int iq;
for(iq = 0; iq < numPoints; iq++) {
if (grdUhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
result[iq] += matrix * grdUhAtQP[iq];
}
}
......@@ -909,23 +875,22 @@ namespace AMDiS {
double *result,
double factor) const
{
int i, j, dow = Global::getGeo(WORLD);
int iq;
int dow = Global::getGeo(WORLD);
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
double resultQP = 0.0;
WorldMatrix<double> A = (*f)(gradAtQPs[iq]);
if(D2UhAtQP) {
for(i=0; i < dow; i++) {
for(j=0; j < dow; j++) {
if (D2UhAtQP) {
for (int i = 0; i < dow; i++) {
for (int j = 0; j < dow; j++) {
resultQP += A[i][j] * D2UhAtQP[iq][j][i];
}
}
}
if(grdUhAtQP) {
if (grdUhAtQP) {
resultQP += (*divFct)(A) * grdUhAtQP[iq];
}
......@@ -937,10 +902,9 @@ namespace AMDiS {
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if(grdUhAtQP) {
int iq;
if (grdUhAtQP) {
WorldMatrix<double> A;
for(iq = 0; iq < numPoints; iq++) {
for (int iq = 0; iq < numPoints; iq++) {
A = (*f)(gradAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
......@@ -954,14 +918,13 @@ namespace AMDiS {
double *result,
double fac) const
{
int i, dow = Global::getGeo(WORLD);
int iq;
int dow = Global::getGeo(WORLD);
if(D2UhAtQP) {
for(iq = 0; iq < numPoints; iq++) {
if (D2UhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
double resultQP = 0.0;
for(i=0; i < dow; i++) {
for (int i = 0; i < dow; i++) {
resultQP += D2UhAtQP[iq][i][i];
}
result[iq] += fac * factor * resultQP;
......@@ -973,9 +936,8 @@ namespace AMDiS {
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if(grdUhAtQP) {
int iq;
for(iq = 0; iq < numPoints; iq++) {
if (grdUhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
......@@ -989,14 +951,13 @@ namespace AMDiS {
double *result,
double f) const
{
int i, dow = Global::getGeo(WORLD);
int iq;
int dow = Global::getGeo(WORLD);
if(D2UhAtQP) {
for(iq = 0; iq < numPoints; iq++) {
if (D2UhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
double resultQP = 0.0;
for(i=0; i < dow; i++) {
for (int i = 0; i < dow; i++) {
resultQP += D2UhAtQP[iq][i][i];
}
result[iq] += factor * f * resultQP;
......@@ -1008,9 +969,8 @@ namespace AMDiS {
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if(grdUhAtQP) {
int iq;
for(iq = 0; iq < numPoints; iq++) {
if (grdUhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
......@@ -1024,15 +984,13 @@ namespace AMDiS {
double *result,
double fac) const
{
int i, dow = Global::getGeo(WORLD);
int iq;
int dow = Global::getGeo(WORLD);
if(D2UhAtQP) {
for(iq = 0; iq < numPoints; iq++) {
if (D2UhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*f)(gradAtQPs[iq]);
double resultQP = 0.0;
for(i=0; i < dow; i++) {
for (int i = 0; i < dow; i++) {
resultQP += D2UhAtQP[iq][i][i];
}
result[iq] += fac * factor * resultQP;
......@@ -1044,9 +1002,8 @@ namespace AMDiS {
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if(grdUhAtQP) {
int iq;
for(iq = 0; iq < numPoints; iq++) {
if (grdUhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*f)(gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
......@@ -1060,15 +1017,13 @@ namespace AMDiS {
double *result,
double fac) const
{
int i, dow = Global::getGeo(WORLD);