From 28431695d227df3063f562ecc9e08d9b3cf75646 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 11 Aug 2010 14:09:38 +0000
Subject: [PATCH] Fixed some bugs with dual mesh. Some more code refactoring.
---
AMDiS/src/Assembler.cc | 92 ++++++++++++++++++---------------
AMDiS/src/Assembler.h | 1 +
AMDiS/src/DOFMatrix.cc | 16 ++++--
AMDiS/src/DOFVector.cc | 3 +-
AMDiS/src/ElInfo2d.cc | 55 +++++++++-----------
AMDiS/src/ElInfo2d.h | 2 +-
AMDiS/src/FirstOrderTerm.cc | 75 ++++++++++++++-------------
AMDiS/src/FirstOrderTerm.h | 14 ++---
AMDiS/src/Operator.cc | 2 +
AMDiS/src/Operator.h | 3 +-
AMDiS/src/ProblemVec.cc | 26 +++++++---
AMDiS/src/SecondOrderTerm.cc | 64 +++++++++++------------
AMDiS/src/SecondOrderTerm.h | 12 ++---
AMDiS/src/ZeroOrderAssembler.cc | 5 +-
14 files changed, 194 insertions(+), 176 deletions(-)
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index c516f5cd..b4b540bd 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -80,6 +80,7 @@ namespace AMDiS {
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
+ bool rowColFeSpaceEqual,
ElementMatrix& userMat,
double factor)
{
@@ -93,9 +94,9 @@ namespace AMDiS {
if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
if (smallElInfo == colElInfo)
- initElement(smallElInfo);
+ initElement(smallElInfo);
else
- initElement(smallElInfo, largeElInfo);
+ initElement(smallElInfo, largeElInfo);
}
if (el != lastMatEl || !operat->isOptimized()) {
@@ -118,65 +119,70 @@ namespace AMDiS {
ElementMatrix &m =
smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
- if (smallElInfo == colElInfo)
- tmpMat = m * mat;
- else
- tmpMat = mat * trans(m);
-
- mat = tmpMat;
+ if (!rowColFeSpaceEqual) {
+ if (smallElInfo == colElInfo)
+ tmpMat = m * mat;
+ else
+ tmpMat = mat * trans(m);
+
+ mat = tmpMat;
+ }
}
if (firstOrderAssemblerGrdPsi) {
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
- if (largeElInfo == rowElInfo) {
- ElementMatrix &m =
- smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-
- tmpMat = m * mat;
- } else {
- ElementMatrix &m =
- smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+ if (!rowColFeSpaceEqual) {
+ if (largeElInfo == rowElInfo) {
+ ElementMatrix &m =
+ smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+
+ tmpMat = m * mat;
+ } else {
+ ElementMatrix &m =
+ smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+
+ tmpMat = mat * trans(m);
+ }
- tmpMat = mat * trans(m);
+ mat = tmpMat;
}
-
- mat = tmpMat;
}
if (firstOrderAssemblerGrdPhi) {
firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
- if (largeElInfo == colElInfo) {
- ElementMatrix &m =
- smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-
- tmpMat = mat * trans(m);
- } else {
- ElementMatrix &m =
- smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+ if (!rowColFeSpaceEqual) {
+ if (largeElInfo == colElInfo) {
+ ElementMatrix &m =
+ smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+
+ tmpMat = mat * trans(m);
+ } else {
+ ElementMatrix &m =
+ smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+
+ tmpMat = m * mat;
+ }
- tmpMat = m * mat;
+ mat = tmpMat;
}
-
- mat = tmpMat;
}
if (zeroOrderAssembler) {
zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
-
- ElementMatrix &m =
- smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
-
- if (smallElInfo == colElInfo)
- tmpMat = m * mat;
- else
- tmpMat = mat;
-// else
-// tmpMat = mat * trans(m);
-
- mat = tmpMat;
+ if (!rowColFeSpaceEqual) {
+ ElementMatrix &m =
+ smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
+
+ if (smallElInfo == colElInfo)
+ tmpMat = m * mat;
+ else
+ tmpMat = mat * trans(m);
+
+ mat = tmpMat;
+ }
}
if (rememberElMat && &userMat != &elementMatrix)
@@ -342,7 +348,7 @@ namespace AMDiS {
if (mainEl != lastMatEl) {
set_to_zero(elementMatrix);
calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo,
- elementMatrix);
+ false, elementMatrix);
}
for (int i = 0; i < nBasFcts; i++) {
diff --git a/AMDiS/src/Assembler.h b/AMDiS/src/Assembler.h
index e6e13a21..97e7f904 100644
--- a/AMDiS/src/Assembler.h
+++ b/AMDiS/src/Assembler.h
@@ -69,6 +69,7 @@ namespace AMDiS {
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
+ bool rowColFeSpaceEqual,
ElementMatrix& userMat,
double factor = 1.0);
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 519aa4bb..10c05cf5 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -286,13 +286,16 @@ namespace AMDiS {
if (op) {
op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo,
- elementMatrix);
+ false, elementMatrix);
} else {
std::vector<Operator*>::iterator it = operators.begin();
std::vector<double*>::iterator factorIt = operatorFactor.begin();
for (; it != operators.end(); ++it, ++factorIt)
- (*it)->getElementMatrix(rowElInfo, colElInfo,
- smallElInfo, largeElInfo,
+ (*it)->getElementMatrix(rowElInfo,
+ colElInfo,
+ smallElInfo,
+ largeElInfo,
+ false,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
@@ -328,8 +331,11 @@ namespace AMDiS {
it != operators.end();
++it, ++factorIt) {
if ((*it)->getNeedDualTraverse()) {
- (*it)->getElementMatrix(mainElInfo, auxElInfo,
- smallElInfo, largeElInfo,
+ (*it)->getElementMatrix(mainElInfo,
+ auxElInfo,
+ smallElInfo,
+ largeElInfo,
+ rowFeSpace == colFeSpace,
elementMatrix,
*factorIt ? **factorIt : 1.0);
}
diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 3e4abce1..bdc099ae 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -776,7 +776,8 @@ namespace AMDiS {
double val = 0.0;
for (int k = 0; k < nBasFcts; k++)
val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(i));
- vecAtQPs[i] += val;
+
+ vecAtQPs[i] += localVec[j] * val;
}
}
}
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 1f8db69b..5dbd2c16 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -108,19 +108,11 @@ namespace AMDiS {
ElInfo2d::ElInfo2d(Mesh *aMesh)
: ElInfo(aMesh)
- {
- e1 = new WorldVector<double>;
- e2 = new WorldVector<double>;
- normal = new WorldVector<double>;
- }
+ {}
ElInfo2d::~ElInfo2d()
- {
- delete e1;
- delete e2;
- delete normal;
- }
+ {}
void ElInfo2d::fillMacroInfo(const MacroElement * mel)
@@ -585,7 +577,7 @@ namespace AMDiS {
}
}
- double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
+ double ElInfo2d::calcGrdLambda(DimVec<WorldVector<double> >& grd)
{
FUNCNAME("ElInfo2d::calcGrdLambda()");
@@ -595,31 +587,30 @@ namespace AMDiS {
int dim = mesh->getDim();
for (int i = 0; i < dimOfWorld; i++) {
- (*e1)[i] = coord[1][i] - coord[0][i];
- (*e2)[i] = coord[2][i] - coord[0][i];
+ e1[i] = coord[1][i] - coord[0][i];
+ e2[i] = coord[2][i] - coord[0][i];
}
if (dimOfWorld == 2) {
- double sdet = (*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0];
+ double sdet = e1[0] * e2[1] - e1[1] * e2[0];
adet = abs(sdet);
if (adet < 1.0E-25) {
MSG("abs(det) = %f\n", adet);
for (int i = 0; i <= dim; i++)
- for (int j = 0; j < dimOfWorld; j++)
- grd_lam[i][j] = 0.0;
+ grd[i].set(0.0);
} else {
double det1 = 1.0 / sdet;
- grd_lam[1][0] = (*e2)[1] * det1; // a11: (a_ij) = A^{-T}
- grd_lam[1][1] = -(*e2)[0] * det1; // a21
- grd_lam[2][0] = -(*e1)[1] * det1; // a12
- grd_lam[2][1] = (*e1)[0] * det1; // a22
- grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
- grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
+ grd[1][0] = e2[1] * det1; // a11: (a_ij) = A^{-T}
+ grd[1][1] = -e2[0] * det1; // a21
+ grd[2][0] = -e1[1] * det1; // a12
+ grd[2][1] = e1[0] * det1; // a22
+ grd[0][0] = -grd[1][0] - grd[2][0];
+ grd[0][1] = -grd[1][1] - grd[2][1];
}
} else {
- vectorProduct(*e1, *e2, *normal);
+ vectorProduct(e1, e2, normal);
adet = norm(normal);
@@ -627,21 +618,21 @@ namespace AMDiS {
MSG("abs(det) = %lf\n", adet);
for (int i = 0; i <= dim; i++)
for (int j = 0; j < dimOfWorld; j++)
- grd_lam[i][j] = 0.0;
+ grd[i][j] = 0.0;
} else {
- vectorProduct(*e2, *normal, grd_lam[1]);
- vectorProduct(*normal, *e1, grd_lam[2]);
+ vectorProduct(e2, normal, grd[1]);
+ vectorProduct(normal, e1, grd[2]);
double adet2 = 1.0 / (adet * adet);
for (int i = 0; i < dimOfWorld; i++) {
- grd_lam[1][i] *= adet2;
- grd_lam[2][i] *= adet2;
+ grd[1][i] *= adet2;
+ grd[2][i] *= adet2;
}
- grd_lam[0][0] = - grd_lam[1][0] - grd_lam[2][0];
- grd_lam[0][1] = - grd_lam[1][1] - grd_lam[2][1];
- grd_lam[0][2] = - grd_lam[1][2] - grd_lam[2][2];
+ grd[0][0] = -grd[1][0] - grd[2][0];
+ grd[0][1] = -grd[1][1] - grd[2][1];
+ grd[0][2] = -grd[1][2] - grd[2][2];
}
}
@@ -710,7 +701,7 @@ namespace AMDiS {
normal[0] = coord[i1][1] - coord[i0][1];
normal[1] = coord[i0][0] - coord[i1][0];
} else { // dow == 3
- WorldVector<double> e0, e1,e2, elementNormal;
+ WorldVector<double> e0, elementNormal;
e0 = coord[i1];
e0 -= coord[i0];
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index 09cb0f88..f50115c6 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/AMDiS/src/ElInfo2d.h
@@ -64,7 +64,7 @@ namespace AMDiS {
protected:
/// Temp vectors for function \ref calcGrdLambda.
- WorldVector<double> *e1, *e2, *normal;
+ WorldVector<double> e1, e2, normal;
static double mat_d1_left_val[3][3];
static mtl::dense2D<double> mat_d1_left;
diff --git a/AMDiS/src/FirstOrderTerm.cc b/AMDiS/src/FirstOrderTerm.cc
index c4c48bb2..c39a63f2 100644
--- a/AMDiS/src/FirstOrderTerm.cc
+++ b/AMDiS/src/FirstOrderTerm.cc
@@ -44,17 +44,17 @@ namespace AMDiS {
void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
+ lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
+ lb(grdLambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
- l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
+ l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
}
}
@@ -92,9 +92,9 @@ namespace AMDiS {
void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq]));
+ l1(grdLambda, Lb[iq], (*g)(coordsAtQPs[iq]));
}
void CoordsAtQP_FOT::eval(int nPoints,
@@ -130,9 +130,9 @@ namespace AMDiS {
void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, *b, Lb[iq], (*g)(coordsAtQPs[iq]));
+ lb(grdLambda, *b, Lb[iq], (*g)(coordsAtQPs[iq]));
}
void VecCoordsAtQP_FOT::eval(int nPoints,
@@ -178,13 +178,13 @@ namespace AMDiS {
void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (f) {
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
+ lb(grdLambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
} else {
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
+ lb(grdLambda, gradAtQPs[iq], Lb[iq], 1.0);
}
}
@@ -230,9 +230,9 @@ namespace AMDiS {
void VectorFct_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
+ lb(grdLambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
}
void VectorFct_FOT::eval(int nPoints,
@@ -264,9 +264,9 @@ namespace AMDiS {
void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
+ lb(grdLambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
}
void VecFctAtQP_FOT::eval(int nPoints,
@@ -322,9 +322,9 @@ namespace AMDiS {
void VecGrad_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const {
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
+ lb(grdLambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
}
void VecGrad_FOT::eval(int nPoints,
@@ -356,9 +356,9 @@ namespace AMDiS {
void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
+ lb(grdLambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
}
void FctGrad2_FOT::eval(int nPoints,
@@ -402,12 +402,12 @@ namespace AMDiS {
}
void Vec2AndGradAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
if (b)
- lb(Lambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
+ lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
else
- l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]));
+ l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]));
}
}
@@ -448,12 +448,12 @@ namespace AMDiS {
void FctVecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
if (b)
- lb(Lambda, *b, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
+ lb(grdLambda, *b, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
else
- l1(Lambda, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
+ l1(grdLambda, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
}
}
@@ -509,17 +509,17 @@ namespace AMDiS {
void Vec2AtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
+ lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else if (b) {
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
+ lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++)
- l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
+ l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
}
}
@@ -584,16 +584,16 @@ namespace AMDiS {
void Vec3FctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
+ lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
} else {
for (int iq = 0; iq < nPoints; iq++) {
if (b)
- lb(Lambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
+ lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
else
- l1(Lambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
+ l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
}
}
}
@@ -663,7 +663,7 @@ namespace AMDiS {
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
for (int i = 0; i < nVecs; i++)
@@ -671,7 +671,7 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- lb(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
+ lb(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
result[iq], 1.0);
}
}
@@ -765,14 +765,15 @@ namespace AMDiS {
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = 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];
- lb(Lambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg), result[iq], 1.0);
+ lb(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
+ result[iq], 1.0);
}
}
diff --git a/AMDiS/src/FirstOrderTerm.h b/AMDiS/src/FirstOrderTerm.h
index 6dfb0b19..457fefba 100644
--- a/AMDiS/src/FirstOrderTerm.h
+++ b/AMDiS/src/FirstOrderTerm.h
@@ -144,10 +144,10 @@ namespace AMDiS {
int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1(Lambda, Lb[iq], 1.0);
+ l1(grdLambda, Lb[iq], 1.0);
}
};
@@ -177,10 +177,10 @@ namespace AMDiS {
/// Implements FirstOrderTerm::getLb().
inline void getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1(Lambda, Lb[iq], (*factor));
+ l1(grdLambda, Lb[iq], (*factor));
}
private:
@@ -214,14 +214,14 @@ namespace AMDiS {
int nPoints,
VectorOfFixVecs<DimVec<double> >&Lb) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
if (bOne > -1) {
for (int iq = 0; iq < nPoints; iq++)
- lb_one(Lambda, Lb[iq], 1.0);
+ lb_one(grdLambda, Lb[iq], 1.0);
} else {
for (int iq = 0; iq < nPoints; iq++)
- lb(Lambda, *b, Lb[iq], 1.0);
+ lb(grdLambda, *b, Lb[iq], 1.0);
}
}
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index 85c6e04d..251d5fd6 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -115,6 +115,7 @@ namespace AMDiS {
void Operator::getElementMatrix(const ElInfo *rowElInfo, const ElInfo *colElInfo,
const ElInfo *smallElInfo, const ElInfo *largeElInfo,
+ bool rowColFeSpaceEqual,
ElementMatrix& userMat,
double factor)
{
@@ -125,6 +126,7 @@ namespace AMDiS {
assembler[myRank]->calculateElementMatrix(rowElInfo, colElInfo,
smallElInfo, largeElInfo,
+ rowColFeSpaceEqual,
userMat, factor);
}
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index aaf9cbf2..ebafa21a 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -101,7 +101,8 @@ namespace AMDiS {
const ElInfo *colElInfo,
const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
- ElementMatrix& userMat,
+ bool rowColFeSpaceEqual,
+ ElementMatrix& userMat,
double factor = 1.0);
/** \brief
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index f4e5eca3..cca99909 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -780,7 +780,7 @@ namespace AMDiS {
{
FUNCNAME("ProblemVec::dualMeshTraverseRequired()");
- TEST_EXIT(meshes.size() <= 2)("More than two mneshes are not supported yet!\n");
+ TEST_EXIT(meshes.size() <= 2)("More than two meshes are not yet supported!\n");
return (meshes.size() == 2);
}
@@ -818,11 +818,22 @@ namespace AMDiS {
traverseInfo.updateStatus();
+ if (writeAsmInfo) {
+ MSG("TraverseInfo:\n");
+ for (int i = 0; i < nComponents; i++) {
+ MSG(" component %d: difAuxSpace = %d\n", i, traverseInfo.difAuxSpace(i));
+
+ for (int j = 0; j < nComponents; j++) {
+ MSG(" component %d-%d: difAuxSpace = %d\n",
+ i, j, traverseInfo.difAuxSpace(i, j));
+ }
+ }
+ }
+
// Used to calculate the overall number of non zero entries.
int nnz = 0;
for (int i = 0; i < nComponents; i++) {
-
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
@@ -927,7 +938,7 @@ namespace AMDiS {
elInfo = dualElInfo.colElInfo;
- if (elInfo != NULL) {
+ if (elInfo != NULL) {
if (useGetBound)
basisFcts->getBound(elInfo, bound);
@@ -952,11 +963,12 @@ namespace AMDiS {
if (traverseInfo.difAuxSpace(i) && i == j)
rhs->getDOFVector(i)->assemble2(1.0, mainElInfo, auxElInfo,
- dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
+ dualElInfo.smallElInfo,
+ dualElInfo.largeElInfo, bound);
if (traverseInfo.difAuxSpace(i, j) && matrix)
matrix->assemble2(1.0, mainElInfo, auxElInfo,
- dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
+ dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
if (matrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->fillBoundaryConditions(mainElInfo, matrix);
@@ -986,8 +998,8 @@ namespace AMDiS {
matrix->getBoundaryManager()->fillBoundaryConditions(rowElInfo, matrix);
}
- if (i == j)
- rhs->getDOFVector(i)->assemble(1.0, rowElInfo, bound);
+ if (i == j)
+ rhs->getDOFVector(i)->assemble(1.0, rowElInfo, bound);
}
}
}
diff --git a/AMDiS/src/SecondOrderTerm.cc b/AMDiS/src/SecondOrderTerm.cc
index e6e7e5aa..d5a52a17 100644
--- a/AMDiS/src/SecondOrderTerm.cc
+++ b/AMDiS/src/SecondOrderTerm.cc
@@ -86,9 +86,9 @@ namespace AMDiS {
void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lalt(Lambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
}
void MatrixFct_SOT::eval(int nPoints,
@@ -192,9 +192,9 @@ namespace AMDiS {
void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
void VecAtQP_SOT::eval(int nPoints,
@@ -254,9 +254,9 @@ namespace AMDiS {
void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
}
void Vec2AtQP_SOT::eval(int nPoints,
@@ -302,9 +302,9 @@ namespace AMDiS {
void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
+ l1lt(grdLambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
}
void CoordsAtQP_SOT::eval(int nPoints,
@@ -376,9 +376,9 @@ namespace AMDiS {
void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lalt(Lambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
void MatrixGradient_SOT::eval(int nPoints,
@@ -445,9 +445,9 @@ namespace AMDiS {
void FctGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
}
void FctGradient_SOT::eval(int nPoints,
@@ -508,9 +508,9 @@ namespace AMDiS {
void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lalt(Lambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
void VecMatrixGradientAtQP_SOT::eval(int nPoints,
@@ -575,9 +575,9 @@ namespace AMDiS {
void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
}
void VecGradCoordsAtQP_SOT::eval(int nPoints,
@@ -635,9 +635,9 @@ namespace AMDiS {
void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
}
void VecAndCoordsAtQP_SOT::eval(int nPoints,
@@ -699,9 +699,9 @@ namespace AMDiS {
void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lalt(Lambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
+ lalt(grdLambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
void MatrixGradientAndCoords_SOT::eval(int nPoints,
@@ -777,9 +777,9 @@ namespace AMDiS {
void MatrixVec2_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lalt(Lambda, A, *(LALt[iq]), symmetric, (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]));
+ lalt(grdLambda, A, *(LALt[iq]), symmetric, (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]));
}
void MatrixVec2_SOT::eval(int nPoints,
@@ -871,7 +871,7 @@ namespace AMDiS {
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
for (int i = 0; i < nVecs; i++)
@@ -879,7 +879,7 @@ namespace AMDiS {
for (int i = 0; i < nGrads; i++)
gradsArg[i] = gradsAtQPs_[i][iq];
- lalt(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
+ lalt(grdLambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg),
*(LALt[iq]), symmetric_, 1.0);
}
}
@@ -1004,14 +1004,14 @@ namespace AMDiS {
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = 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];
- lalt(Lambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
+ lalt(grdLambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg),
*(LALt[iq]), symmetric_, 1.0);
}
}
@@ -1099,9 +1099,9 @@ namespace AMDiS {
void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
void VecAndGradAtQP_SOT::eval(int nPoints,
@@ -1149,9 +1149,9 @@ namespace AMDiS {
void VecGrad_SOT::getLALt(const ElInfo *elInfo, int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
+ l1lt(grdLambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
}
void VecGrad_SOT::eval(int nPoints,
@@ -1254,9 +1254,9 @@ namespace AMDiS {
int nPoints,
DimMat<double> **LALt) const
{
- const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) {
- lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*f)(vecAtQPs[iq]));
+ lalt_kl(grdLambda, xi, xj, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
}
diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h
index 63a87674..0f6b15c5 100644
--- a/AMDiS/src/SecondOrderTerm.h
+++ b/AMDiS/src/SecondOrderTerm.h
@@ -120,10 +120,10 @@ namespace AMDiS {
/// Implenetation of 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> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), 1.0);
+ l1lt(grdLambda, *(LALt[iq]), 1.0);
}
/// Implementation of SecondOrderTerm::eval().
@@ -188,9 +188,9 @@ namespace AMDiS {
/// 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> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- l1lt(Lambda, *(LALt[iq]), (*factor));
+ l1lt(grdLambda, *(LALt[iq]), (*factor));
}
/// Implenetation of SecondOrderTerm::eval().
@@ -360,9 +360,9 @@ namespace AMDiS {
/// 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> > &grdLambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++)
- lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor));
+ lalt_kl(grdLambda, xi, xj, *(LALt[iq]), (*factor));
}
/// Implenetation of SecondOrderTerm::eval().
diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc
index 017115d1..90fc7c0c 100644
--- a/AMDiS/src/ZeroOrderAssembler.cc
+++ b/AMDiS/src/ZeroOrderAssembler.cc
@@ -63,9 +63,6 @@ namespace AMDiS {
}
}
- newAssembler = new StandardZOA(op, assembler, quad);
-
-
// create new assembler
if (!optimized) {
newAssembler = new StandardZOA(op, assembler, quad);
@@ -101,7 +98,7 @@ namespace AMDiS {
std::vector<OperatorTerm*>::iterator termIt;
for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
(static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
-
+
if (symmetric) {
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
--
GitLab