diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc index 03e16e3e7d7a88170b943c5e950eb8750ad3e43b..b9470a61c556d227cdc33c3d3eb0abbb42fbfea2 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 5c80ac11052cd7fdc35bc35d6fd6e0712db017c4..2ddefba00fbd770e3ce86d4b8f33346768244f70 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 795a3b6bc3d8955bd858249b394512e57b894b28..01360e1105af9e7ad1c60cdc273d7b72d93d507c 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 f07620c2b986b560d2b70081f7dea2fbd096520d..e876ab19884d5a49e661e06832bfb25250a0cc77 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 d8680438061f82be376d61ba92c1d7839621ce70..0217373e318caf2ef15dfe592d062bc03dedffd3 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 cd11943820649c64fd17a34b1ae35fa2fbef6cf5..3d99e78b83c5bb3fa23ad685c1f6320fc6328aff 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;