From bbd936aa6e4a48ce0fb362e58fd24efa5c5966bc Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 14 Oct 2009 08:04:37 +0000
Subject: [PATCH] Fixed a bug for FOT operators when using 2d meshes in 3d.
---
AMDiS/src/FixVec.h | 18 +-
AMDiS/src/Operator.cc | 405 ++++++++++++++++++++----------------------
AMDiS/src/Operator.h | 87 +++------
3 files changed, 226 insertions(+), 284 deletions(-)
diff --git a/AMDiS/src/FixVec.h b/AMDiS/src/FixVec.h
index c079aa8e..16e10025 100644
--- a/AMDiS/src/FixVec.h
+++ b/AMDiS/src/FixVec.h
@@ -277,7 +277,8 @@ namespace AMDiS {
* VectorOfFixVecs constructors.
*/
MatrixOfFixVecs(int dim, int r, int c, InitType initType)
- : rows(r), columns(c)
+ : rows(r),
+ columns(c)
{
TEST_EXIT_DBG(initType == NO_INIT)("wrong initType or wrong initializer\n");
vec = new VectorOfFixVecs<FixVecType>*[rows];
@@ -309,12 +310,14 @@ namespace AMDiS {
}
/// Returns \ref rows
- inline int getNumberOfRows() const {
+ inline int getNumberOfRows() const
+ {
return rows;
}
/// Returns \ref columns
- inline int getNumberOfColumns() const {
+ inline int getNumberOfColumns() const
+ {
return columns;
}
@@ -337,7 +340,8 @@ namespace AMDiS {
* parts of an element.
*/
template<typename T>
- class DimVec : public FixVec<T,PARTS> {
+ class DimVec : public FixVec<T,PARTS>
+ {
public:
DimVec() {}
@@ -549,8 +553,7 @@ namespace AMDiS {
return result;
}
- inline bool operator<(const WorldVector<double>& v1,
- const WorldVector<double>& v2)
+ inline bool operator<(const WorldVector<double>& v1, const WorldVector<double>& v2)
{
int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++) {
@@ -561,8 +564,7 @@ namespace AMDiS {
return false;
}
- inline bool operator==(const WorldVector<double>& v1,
- const WorldVector<double>& v2)
+ inline bool operator==(const WorldVector<double>& v1, const WorldVector<double>& v2)
{
int dow = Global::getGeo(WORLD);
for (int i = 0; i < dow; i++)
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 92551535..fb8d01e5 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -53,11 +53,10 @@ namespace AMDiS {
void OperatorTerm::setSymmetric(bool symm)
{
- if (symm) {
+ if (symm)
properties.setFlag(SYMMETRIC);
- } else {
- properties.unsetFlag(SYMMETRIC);
- }
+ else
+ properties.unsetFlag(SYMMETRIC);
}
bool OperatorTerm::isSymmetric()
@@ -145,23 +144,20 @@ namespace AMDiS {
// Both elements have the same size, so we can use the simple procedure
// to determine the gradients.
- if (vec->getFESpace()->getMesh() == smallElInfo->getMesh()) {
+ if (vec->getFESpace()->getMesh() == smallElInfo->getMesh())
return subAssembler->getGradientsAtQPs(vec, smallElInfo, quad);
- } else {
- return subAssembler->getGradientsAtQPs(vec, largeElInfo, quad);
- }
+ else
+ return subAssembler->getGradientsAtQPs(vec, largeElInfo, quad);
} else {
// The two elements are different. If the vector is defined on the mesh of the
// small element, we can still use the simple procedure to determine the gradients.
- if (vec->getFESpace()->getMesh() == largeElInfo->getMesh()) {
+ if (vec->getFESpace()->getMesh() == largeElInfo->getMesh())
return subAssembler->getGradientsAtQPs(vec, smallElInfo, largeElInfo, quad);
- } else {
+ else
return subAssembler->getGradientsAtQPs(vec, smallElInfo, quad);
- }
-
}
}
@@ -807,7 +803,7 @@ namespace AMDiS {
if (b)
lb(Lambda, *b, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
else
- l1(Lambda, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
+ l1(Lambda, Lb[iq], (*f)(coordsAtQPs[iq]) * (vecAtQPs[iq]));
}
}
@@ -1423,143 +1419,164 @@ namespace AMDiS {
void CoordsAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
- Quadrature *quad) {
+ Quadrature *quad)
+ {
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
- void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
+ 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++)
lalt(Lambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
}
- void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
+ 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++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
- void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
+ 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++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
}
- void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
+ 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++)
l1lt(Lambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
}
- void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
+ 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 {
+ 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++)
l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
}
void VecAndGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const {
+ DimMat<double> **LALt) const
+ {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
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 {
+ DimMat<double> **LALt) const
+ {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
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();
+ 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++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
}
void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int nPoints,
- DimMat<double> **LALt) const {
+ DimMat<double> **LALt) const
+ {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
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 {
+ DimMat<double> **LALt) const
+ {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
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 {
+ void VecAtQP_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++) {
if (b)
lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
else
- l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
+ l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
}
}
- void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void CoordsAtQP_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++)
l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq]));
- }
}
void VecCoordsAtQP_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, b, Lb[iq], (*g)(coordsAtQPs[iq]));
- }
}
- void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >& Lb) const
+ {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
if (f) {
- for (int iq = 0; iq < nPoints; iq++) {
- lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
} else {
- for (int iq = 0; iq < nPoints; iq++) {
- lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
}
}
- void VectorFct_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void VectorFct_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++) {
- lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
}
- void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
+ void VecFctAtQP_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++) {
- lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
}
- void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
+ void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
+ {
if (f) {
- for (int iq = 0; iq < nPoints; iq++) {
- C[iq] += (*f)(vecAtQPs[iq]);
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs[iq]);
} else {
- for (int iq = 0; iq < nPoints; iq++) {
- C[iq] += vecAtQPs[iq];
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += vecAtQPs[iq];
}
}
@@ -1571,21 +1588,19 @@ namespace AMDiS {
double fac) const
{
if (f) {
- for (int iq = 0; iq < nPoints; iq++) {
- result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
} else {
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
result[iq] += fac * vecAtQPs[iq] * uhAtQP[iq];
- }
}
}
- void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
- for (int iq = 0; iq < nPoints; iq++) {
- C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
- }
+ void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
}
void MultVecAtQP_ZOT::eval(int nPoints,
@@ -1595,16 +1610,14 @@ namespace AMDiS {
double *result,
double fac) const
{
- for (int iq = 0; iq < nPoints; iq++) {
- result[iq] += fac * (*f1)(vecAtQPs1[iq]) *
- (*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac * (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
}
- void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
- for (int iq = 0; iq < nPoints; iq++) {
- C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
- }
+ void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
}
void Vec2AtQP_ZOT::eval(int nPoints,
@@ -1614,16 +1627,14 @@ namespace AMDiS {
double *result,
double fac) const
{
- for (int iq = 0; iq < nPoints; iq++) {
- result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) *
- uhAtQP[iq];
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq];
}
- void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const {
- for (int iq = 0; iq < nPoints; iq++) {
- C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
- }
+ void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
}
void Vec3AtQP_ZOT::eval(int nPoints,
@@ -1633,20 +1644,17 @@ namespace AMDiS {
double *result,
double fac) const
{
- for (int iq = 0; iq < nPoints; iq++) {
+ for (int iq = 0; iq < nPoints; iq++)
result[iq] +=
- fac *
- (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) *
- uhAtQP[iq];
- }
+ fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * uhAtQP[iq];
}
void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
- std::vector<double> &C) const {
- for (int iq = 0; iq < nPoints; iq++) {
- C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
- }
+ std::vector<double> &C) const
+ {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
}
void VecAndCoordsAtQP_ZOT::eval(int nPoints,
@@ -2302,9 +2310,9 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
double resultQP = 0.0;
- for (int i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i];
- }
+
result[iq] += f * factor * resultQP;
}
}
@@ -2324,9 +2332,8 @@ namespace AMDiS {
result[iq] += b * grdUhAtQP[iq] * factor;
}
} else {
- for (int iq = 0; iq < nPoints; iq++) {
- result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor;
- }
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor;
}
}
}
@@ -2359,9 +2366,9 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
const WorldVector<double> &b = (*g)(coordsAtQPs[iq]);
- for (int i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++)
resultQP += b[i] * grdUhAtQP[iq][i];
- }
+
result[iq] += f * resultQP;
}
}
@@ -2369,9 +2376,8 @@ namespace AMDiS {
Assembler *Operator::getAssembler(int rank)
{
- if (!assembler[rank]) {
- initAssembler(rank, NULL, NULL, NULL, NULL);
- }
+ if (!assembler[rank])
+ initAssembler(rank, NULL, NULL, NULL, NULL);
return assembler[rank];
}
@@ -2393,9 +2399,8 @@ namespace AMDiS {
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_kl(Lambda, xi, xj, *(LALt[iq]), (*g)(coordsAtQPs[iq]));
- }
}
void CoordsAtQP_IJ_SOT::eval(int nPoints,
@@ -2556,9 +2561,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)(arg) * uhAtQP[iq];
}
}
@@ -2624,7 +2629,9 @@ namespace AMDiS {
}
}
- void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
+ void VecAndGradAtQP_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++)
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
@@ -2666,12 +2673,11 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
+
lalt(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
*(LALt[iq]), symmetric_, 1.0);
}
@@ -2694,26 +2700,21 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
- for (int i = 0; i < nVecs; i++) {
+ for (int i = 0; i < nVecs; i++)
vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- }
WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
if (D2UhAtQP) {
- for (int i = 0; i < dow; i++) {
- for (int j = 0; j < dow; j++) {
+ for (int i = 0; i < dow; i++)
+ for (int j = 0; j < dow; j++)
resultQP += A[i][j] * D2UhAtQP[iq][j][i];
- }
- }
}
- if (grdUhAtQP) {
- resultQP += (*divFct_)(A) * grdUhAtQP[iq];
- }
+ if (grdUhAtQP)
+ resultQP += (*divFct_)(A) * grdUhAtQP[iq];
result[iq] += resultQP * factor;
}
@@ -2732,12 +2733,11 @@ namespace AMDiS {
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- }
+
A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
result[iq] += A * grdUhAtQP[iq];
}
@@ -2753,12 +2753,10 @@ namespace AMDiS {
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
- for (int i = 0; i < nVecs; i++) {
- vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
- }
- for (int i = 0; i < nGrads; i++) {
- gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
+ for (int i = 0; i < nGrads; i++)
+ gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
void General_FOT::getLb(const ElInfo *elInfo,
@@ -2774,12 +2772,11 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
+
lb(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
result[iq], 1.0);
}
@@ -2803,18 +2800,16 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
- for (int i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i] * b[i];
- }
+
result[iq] += factor * resultQP;
}
}
@@ -2834,8 +2829,7 @@ namespace AMDiS {
gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
- void General_ZOT::getC(const ElInfo *, int nPoints,
- std::vector<double> &C) const
+ void General_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
@@ -2847,7 +2841,6 @@ namespace AMDiS {
for (int i = 0; i < nVecs; i++)
vecsArg[i] = vecsAtQPs_[i][iq];
for (int i = 0; i < nGrads; i++)
-
gradsArg[i] = gradsAtQPs_[i][iq];
C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
@@ -2868,19 +2861,14 @@ namespace AMDiS {
std::vector<WorldVector<double> > gradsArg(nGrads);
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
result[iq] += fac * (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg) * uhAtQP[iq];
}
}
-
-
- //************************************ parametric general Operatorterms *****************
void GeneralParametric_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
@@ -2891,16 +2879,13 @@ namespace AMDiS {
elInfo->getElementNormal(elementNormal_);
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
- for (int i = 0; i < nVecs; i++) {
- vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
- }
- for (int i = 0; i < nGrads; i++) {
- gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
+ for (int i = 0; i < nGrads; i++)
+ gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
- void GeneralParametric_SOT::getLALt(const ElInfo *elInfo,
- int nPoints,
+ void GeneralParametric_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
int nVecs = static_cast<int>(vecs_.size());
@@ -2912,12 +2897,10 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
lalt(Lambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
*(LALt[iq]), symmetric_, 1.0);
}
@@ -2940,26 +2923,22 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
- WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
+ WorldMatrix<double> A =
+ (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
if (D2UhAtQP) {
- for (int i = 0; i < dow; i++) {
- for (int j = 0; j < dow; j++) {
+ for (int i = 0; i < dow; i++)
+ for (int j = 0; j < dow; j++)
resultQP += A[i][j] * D2UhAtQP[iq][j][i];
- }
- }
}
- if (grdUhAtQP) {
- resultQP += (*divFct_)(A) * grdUhAtQP[iq];
- }
+ if (grdUhAtQP)
+ resultQP += (*divFct_)(A) * grdUhAtQP[iq];
result[iq] += resultQP * factor;
}
@@ -2978,12 +2957,11 @@ namespace AMDiS {
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- }
+
A = (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
result[iq] += A * grdUhAtQP[iq];
}
@@ -3000,12 +2978,10 @@ namespace AMDiS {
elInfo->getElementNormal(elementNormal_);
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
- for (int i = 0; i < nVecs; i++) {
- vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
- }
- for (int i = 0; i < nGrads; i++) {
- gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
- }
+ for (int i = 0; i < nVecs; i++)
+ vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
+ for (int i = 0; i < nGrads; i++)
+ gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
}
void GeneralParametric_FOT::getLb(const ElInfo *elInfo,
@@ -3020,12 +2996,11 @@ namespace AMDiS {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
- for (int i = 0; i < nVecs; i++) {
- vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nVecs; i++)
+ vecsArg[i] = vecsAtQPs_[i][iq];
+ for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- }
+
lb(Lambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg), result[iq], 1.0);
}
}
@@ -3048,25 +3023,23 @@ namespace AMDiS {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
- for (int i = 0; i < nVecs; i++) {
+ for (int i = 0; i < nVecs; i++)
vecsArg[i] = vecsAtQPs_[i][iq];
- }
- for (int i = 0; i < nGrads; i++) {
- gradsArg[i] = gradsAtQPs_[i][iq];
- }
+ for (int i = 0; i < nGrads; i++)
+ gradsArg[i] = gradsAtQPs_[i][iq];
const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg);
- for (int i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i] * b[i];
- }
+
result[iq] += factor * resultQP;
}
}
}
- void GeneralParametric_ZOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
- Quadrature *quad)
+ void GeneralParametric_ZOT::initElement(const ElInfo* elInfo, SubAssembler*
+ subAssembler, Quadrature *quad)
{
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
@@ -3074,12 +3047,10 @@ namespace AMDiS {
elInfo->getElementNormal(elementNormal_);
coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);
- for (int i = 0; i < nVecs; i++) {
+ for (int i = 0; i < nVecs; i++)
vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
- }
- for (int i = 0; i < nGrads; i++) {
+ for (int i = 0; i < nGrads; i++)
gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
- }
}
void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints,
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index cfeacab5..6065b87a 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -172,8 +172,9 @@ namespace AMDiS {
double factor)
{
const int dimOfWorld = Global::getGeo(WORLD);
-
- for (int i = 0; i <= dimOfWorld; i++) {
+ const int dim = Lambda.size() - 1;
+
+ for (int i = 0; i <= dim; i++) {
double val = 0.0;
for (int j = 0; j < dimOfWorld; j++)
@@ -192,8 +193,9 @@ namespace AMDiS {
double factor)
{
const int dimOfWorld = Global::getGeo(WORLD);
+ const int dim = Lambda.size() - 1;
- for (int i = 0; i <= dimOfWorld; i++) {
+ for (int i = 0; i <= dim; i++) {
double val = 0.0;
for (int j = 0; j < dimOfWorld; j++)
@@ -476,9 +478,7 @@ namespace AMDiS {
setSymmetric(symmetric);
}
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
{
const DimVec<WorldVector<double> >& Lambda = elInfo->getGrdLambda();
@@ -486,9 +486,7 @@ namespace AMDiS {
lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
}
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -496,22 +494,16 @@ namespace AMDiS {
double *result,
double factor) const;
- /** \brief
- * Implenetation of SecondOrderTerm::weakEval().
- */
+ /// Implenetation of SecondOrderTerm::weakEval().
void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
protected:
- /** \brief
- * Matrix stroring A.
- */
+ /// Matrix stroring A.
WorldMatrix<double> matrix;
- /** \brief
- * True, if \ref matrix is symmetric.
- */
+ /// True, if \ref matrix is symmetric.
bool symmetric;
};
@@ -543,12 +535,10 @@ namespace AMDiS {
setSymmetric(xi == xj);
}
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor));
}
@@ -1357,21 +1347,15 @@ namespace AMDiS {
: OperatorTerm(deg)
{}
- /** \brief
- * Destructor.
- */
+ /// Destructor.
virtual ~FirstOrderTerm() {}
- /** \brief
- * Evaluation of \f$ \Lambda b \f$.
- */
+ /// Evaluation of \f$ \Lambda b \f$.
virtual void getLb(const ElInfo *elInfo,
int nPoints,
VectorOfFixVecs<DimVec<double> >& result) const = 0;
- /** \brief
- * Implenetation of FirstOrderTerm::eval().
- */
+ /// Implenetation of FirstOrderTerm::eval().
void eval(int nPoints,
const double *,
const WorldVector<double> *grdUhAtQP,
@@ -1384,9 +1368,9 @@ namespace AMDiS {
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double resultQP = 0.0;
- for (int i = 0; i < dow; i++) {
+ for (int i = 0; i < dow; i++)
resultQP += grdUhAtQP[iq][i];
- }
+
result[iq] += resultQP * factor;
}
}
@@ -1409,9 +1393,7 @@ namespace AMDiS {
: FirstOrderTerm(0)
{}
- /** \brief
- * Implements FirstOrderTerm::getLb().
- */
+ /// Implements FirstOrderTerm::getLb().
inline void getLb(const ElInfo *elInfo,
int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
@@ -1491,9 +1473,7 @@ namespace AMDiS {
lb(Lambda, b, Lb[iq], 1.0);
}
- /** \brief
- * Implements FirstOrderTerm::eval().
- */
+ /// Implements FirstOrderTerm::eval().
void eval(int nPoints,
const double *,
const WorldVector<double> *grdUhAtQP,
@@ -1502,14 +1482,12 @@ namespace AMDiS {
double factor) const
{
if (grdUhAtQP)
- for(int iq = 0; iq < nPoints; iq++)
+ for (int iq = 0; iq < nPoints; iq++)
result[iq] += b * grdUhAtQP[iq] * factor;
}
protected:
- /** \brief
- * Vector which is multiplied with \f$ \nabla u(\vec{x}) \f$
- */
+ /// Vector which is multiplied with \f$ \nabla u(\vec{x}) \f$
WorldVector<double> b;
};
@@ -1532,7 +1510,8 @@ namespace AMDiS {
Quadrature *quad = NULL);
/// Implements FirstOrderTerm::getLb().
- void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const;
+ void getLb(const ElInfo *elInfo, int nPoints,
+ VectorOfFixVecs<DimVec<double> >& Lb) const;
/// Implements FirstOrderTerm::eval().
void eval(int nPoints,
@@ -2224,9 +2203,7 @@ namespace AMDiS {
/// Implements ZeroOrderTerm::getC().
void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
- /** \brief
- * Implements ZeroOrderTerm::eval().
- */
+ /// Implements ZeroOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -2235,19 +2212,13 @@ namespace AMDiS {
double fac) const;
protected:
- /** \brief
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec;
- /** \brief
- * Vector v at quadrature points.
- */
+ /// Vector v at quadrature points.
double *vecAtQPs;
- /** \brief
- * Function for c.
- */
+ /// Function for c.
AbstractFunction<double, double> *f;
};
@@ -2266,9 +2237,7 @@ namespace AMDiS {
AbstractFunction<double, double> *f1,
AbstractFunction<double, double> *f2);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
--
GitLab