From 869d35d6991d5c7350b2283bb4620e189a9098e5 Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 12 Apr 2010 12:36:45 +0000
Subject: [PATCH] Bugfix in usage of FOT operators in multi-mesh method.
---
AMDiS/src/Assembler.cc | 6 ++++--
AMDiS/src/ElInfo1d.cc | 2 +-
AMDiS/src/ElInfo2d.cc | 22 ++++++++++++++++++++--
AMDiS/src/FirstOrderAssembler.cc | 16 ++++++++--------
AMDiS/src/FirstOrderTerm.cc | 11 ++++++++++-
AMDiS/src/FirstOrderTerm.h | 6 ++++++
6 files changed, 49 insertions(+), 14 deletions(-)
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 03e16e3e..b9470a61 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -120,9 +120,11 @@ namespace AMDiS {
}
if (firstOrderAssemblerGrdPsi) {
+ // std::cout << "PSI!" << std::endl;
+
firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
- if (largeElInfo == colElInfo) {
+ if (largeElInfo == rowElInfo) {
ElementMatrix &m =
smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
@@ -140,7 +142,7 @@ namespace AMDiS {
if (firstOrderAssemblerGrdPhi) {
firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
- if (largeElInfo == rowElInfo) {
+ if (largeElInfo == colElInfo) {
ElementMatrix &m =
smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
diff --git a/AMDiS/src/ElInfo1d.cc b/AMDiS/src/ElInfo1d.cc
index 5c80ac11..2ddefba0 100644
--- a/AMDiS/src/ElInfo1d.cc
+++ b/AMDiS/src/ElInfo1d.cc
@@ -353,7 +353,7 @@ namespace AMDiS {
using namespace mtl;
- if (subElemMatrices[degree].count(refinementPath) == 0) {
+ if (subElemGradMatrices[degree].count(refinementPath) == 0) {
dense2D<double> mat(mat_d1);
for (int i = 0; i < refinementPathLength; i++)
diff --git a/AMDiS/src/ElInfo2d.cc b/AMDiS/src/ElInfo2d.cc
index 795a3b6b..01360e11 100644
--- a/AMDiS/src/ElInfo2d.cc
+++ b/AMDiS/src/ElInfo2d.cc
@@ -731,11 +731,29 @@ namespace AMDiS {
using namespace mtl;
- if (subElemMatrices[degree].count(refinementPath) == 0) {
+ if (subElemGradMatrices[degree].count(refinementPath) == 0) {
dense2D<double> mat(mat_d1);
+ dense2D<double> tmpMat(num_rows(mat), num_rows(mat));
+
+ double test_left[3][3] = {{0.0, 0.0, 0.5},
+ {-0.5, -0.5, 0.0},
+ {1.0, 0.0, 0.0}};
+ double test_right[3][3] = {{0.0, 0.0, 0.5},
+ {0.5, -0.5, 0.0},
+ {0.0, 1.0, 0.0}};
+
+ mtl::dense2D<double> mat_left(test_left);
+ mtl::dense2D<double> mat_right(test_right);
+
for (int i = 0; i < refinementPathLength; i++)
- mat *= 0.5;
+ if (refinementPath & (1 << i)) {
+ tmpMat = mat_right * mat;
+ mat = tmpMat;
+ } else {
+ tmpMat = mat_left * mat;
+ mat = tmpMat;
+ }
subElemGradMatrices[1][refinementPath] = mat;
}
diff --git a/AMDiS/src/FirstOrderAssembler.cc b/AMDiS/src/FirstOrderAssembler.cc
index f07620c2..e876ab19 100644
--- a/AMDiS/src/FirstOrderAssembler.cc
+++ b/AMDiS/src/FirstOrderAssembler.cc
@@ -85,15 +85,15 @@ namespace AMDiS {
dynamic_cast<FirstOrderAssembler*>(new Stand01(op, assembler, quad));
} else {
if (pwConst) {
- newAssembler =
- (type == GRD_PSI) ?
- dynamic_cast<FirstOrderAssembler*>(new Pre10(op, assembler, quad)) :
- dynamic_cast<FirstOrderAssembler*>(new Pre01(op, assembler, quad));
+ newAssembler =
+ (type == GRD_PSI) ?
+ dynamic_cast<FirstOrderAssembler*>(new Pre10(op, assembler, quad)) :
+ dynamic_cast<FirstOrderAssembler*>(new Pre01(op, assembler, quad));
} else {
- newAssembler =
- (type == GRD_PSI) ?
- dynamic_cast<FirstOrderAssembler*>(new Quad10(op, assembler, quad)) :
- dynamic_cast<FirstOrderAssembler*>(new Quad01(op, assembler, quad));
+ newAssembler =
+ (type == GRD_PSI) ?
+ dynamic_cast<FirstOrderAssembler*>(new Quad10(op, assembler, quad)) :
+ dynamic_cast<FirstOrderAssembler*>(new Quad01(op, assembler, quad));
}
}
diff --git a/AMDiS/src/FirstOrderTerm.cc b/AMDiS/src/FirstOrderTerm.cc
index d8680438..0217373e 100644
--- a/AMDiS/src/FirstOrderTerm.cc
+++ b/AMDiS/src/FirstOrderTerm.cc
@@ -33,6 +33,14 @@ namespace AMDiS {
vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
}
+ void VecAtQP_FOT::initElement(const ElInfo* smallElInfo,
+ const ElInfo* largeElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vecAtQPs = getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
+ }
+
void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const
{
@@ -64,7 +72,8 @@ namespace AMDiS {
double factor = (*f)(vecAtQPs[iq]);
double resultQP = 0.0;
for (int i = 0; i < dow; i++)
- resultQP += grdUhAtQP[iq][i];
+ resultQP += grdUhAtQP[iq][i];
+
result[iq] += fac * factor * resultQP;
}
}
diff --git a/AMDiS/src/FirstOrderTerm.h b/AMDiS/src/FirstOrderTerm.h
index cd119438..3d99e78b 100644
--- a/AMDiS/src/FirstOrderTerm.h
+++ b/AMDiS/src/FirstOrderTerm.h
@@ -266,6 +266,12 @@ namespace AMDiS {
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
+ /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
+ void initElement(const ElInfo* smallElInfo,
+ const ElInfo* largeElInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
/// Implements FirstOrderTerm::getLb().
void getLb(const ElInfo *elInfo, int nPoints,
VectorOfFixVecs<DimVec<double> >& Lb) const;
--
GitLab