Commit 0d049e2e authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Work on residual estimator for multiple meshes in problem definition.

parent cc0849b7
......@@ -299,9 +299,8 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
uhOldLoc2[i] = 0.0;
for (int j = 0; j < nBasFcts; j++) {
for (int j = 0; j < nBasFcts; j++)
uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
}
}
......@@ -312,9 +311,8 @@ namespace AMDiS {
for (int i = 0; i < nBasFcts; i++) {
double val = 0.0;
for (int j = 0; j < nBasFcts; j++) {
for (int j = 0; j < nBasFcts; j++)
val += elementMatrix[i][j] * uhOldLoc2[j];
}
vec[i] += val;
}
......@@ -322,6 +320,7 @@ namespace AMDiS {
delete [] uhOldLoc2;
}
void Assembler::initElement(const ElInfo *smallElInfo,
const ElInfo *largeElInfo,
Quadrature *quad)
......@@ -336,6 +335,7 @@ namespace AMDiS {
zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
}
void Assembler::checkQuadratures()
{
if (secondOrderAssembler) {
......
......@@ -165,6 +165,7 @@ namespace AMDiS {
}
}
void DualTraverse::fillSubElInfo(ElInfo *elInfo1,
ElInfo *elInfo2,
ElInfo *elInfoSmall,
......@@ -173,19 +174,36 @@ namespace AMDiS {
if (!fillSubElemMat)
return;
// mtl::dense2D<double>& subCoordsMat =
// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat());
// mtl::dense2D<double>& subCoordsMat_so =
// const_cast<mtl::dense2D<double>&>(elInfoSmall->getSubElemCoordsMat_so());
if (elInfo1 == elInfoSmall) {
// stack1.getCoordsInElem(elInfo2, basisFcts, subCoordsMat);
// stack1.getCoordsInElem_so(elInfo2, basisFcts, subCoordsMat_so);
if (elInfo1 == elInfoSmall)
stack1.fillRefinementPath(*elInfoSmall, *elInfo2);
} else {
// stack2.getCoordsInElem(elInfo1, basisFcts, subCoordsMat);
// stack2.getCoordsInElem_so(elInfo1, basisFcts, subCoordsMat_so);
stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
}
else
stack2.fillRefinementPath(*elInfoSmall, *elInfo1);
}
int DualTraverse::getFace(DualElInfo *dualElInfo, int largeFace)
{
FUNCNAME("DualTraverse::getFace()");
TEST_EXIT_DBG(dualElInfo)("No dual element info object!\n");
ElInfo *largeElInfo = dualElInfo->largeElInfo;
ElInfo *smallElInfo = dualElInfo->smallElInfo;
TEST_EXIT_DBG(largeElInfo->getLevel() <= smallElInfo->getLevel())
("Should not happen!\n");
if (largeElInfo->getLevel() == smallElInfo->getLevel())
return largeFace;
int refinementPathLength = smallElInfo->getRefinementPathLength();
unsigned long refinementPath = smallElInfo->getRefinementPath();
TEST_EXIT_DBG(refinementPathLength != 0)
("Refinement path is empty. This should not happen!\n");
std::cout << "REFINEMENTPATHLENGTH = " << refinementPathLength << std::endl;
return -1;
}
}
......@@ -115,6 +115,16 @@ namespace AMDiS {
fillSubElemMat = b;
basisFcts = fcts;
}
/** \brief
* Checks if the small element has an edge/face which is part of a given edge/face
* of the large element. If this is the case, it returns the local number of the
* small edge/face, and -1 otherwise.
*
* \param[in] dualElInfo Dual element info with large and small element infos.
* \param[in] largeFace A local edge/face number on the large element.
*/
static int getFace(DualElInfo *dualElInfo, int largeFace);
protected:
/** \brief
......
......@@ -227,6 +227,19 @@ namespace AMDiS {
return parametric;
}
/// Returns \ref refinementPath
inline unsigned long getRefinementPath() const
{
return refinementPath;
}
/// Get \ref refinementPathLength
inline int getRefinementPathLength() const
{
return refinementPathLength;
}
virtual void getSubElemCoordsMat(mtl::dense2D<double>& mat,
int basisFctsDegree) const {}
......
......@@ -716,12 +716,9 @@ namespace AMDiS {
if (refinementPath & (1 << i)) {
tmpMat = mat_d1_right * mat;
mat = tmpMat;
// mat *= mat_d1_right;
} else {
tmpMat = mat_d1_left * mat;
mat = tmpMat;
// mat *= mat_d1_left;
}
}
}
......
......@@ -132,6 +132,8 @@ namespace AMDiS {
{
FUNCNAME("OperatorTerm::getGradientsAtQPs()");
ERROR_EXIT("Not yet tested!\n");
TEST_EXIT(smallElInfo->getMesh() == vec->getFESpace()->getMesh() ||
largeElInfo->getMesh() == vec->getFESpace()->getMesh())
("There is something wrong!\n");
......@@ -1026,7 +1028,12 @@ namespace AMDiS {
AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
bool symm)
: SecondOrderTerm(af->getDegree()), vec(dv), f(af), divFct(divAf), symmetric(symm)
: SecondOrderTerm(af->getDegree()),
vec(dv),
f(af),
divFct(divAf),
gradAtQPs(NULL),
symmetric(symm)
{
setSymmetric(symmetric);
......@@ -1163,13 +1170,12 @@ namespace AMDiS {
}
}
void MatrixVec2_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void MatrixVec2_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++)
result[iq] += A * grdUhAtQP[iq] * (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]);
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
result[iq] += A * grdUhAtQP[iq] * (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]);
}
General_SOT::General_SOT(std::vector<DOFVectorBase<double>*> vecs,
......@@ -1302,6 +1308,14 @@ namespace AMDiS {
gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
}
void MatrixGradient_SOT::initElement(const ElInfo* smallElInfo,
const ElInfo* largeElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
gradAtQPs = getGradientsAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
}
void FctGradient_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -1944,17 +1958,15 @@ namespace AMDiS {
}
}
void MatrixFct_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void MatrixFct_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*matrixFct)(vecAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
}
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*matrixFct)(vecAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
}
......@@ -1979,13 +1991,12 @@ namespace AMDiS {
}
}
void Matrix_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void Matrix_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP)
for (int iq = 0; iq < nPoints; iq++)
result[iq] += matrix * grdUhAtQP[iq];
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++)
result[iq] += matrix * grdUhAtQP[iq];
}
void MatrixGradient_SOT::eval(int nPoints,
......@@ -2014,19 +2025,24 @@ namespace AMDiS {
}
}
void MatrixGradient_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void MatrixGradient_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*f)(gradAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
FUNCNAME("MatrixGradient_SOT::weakEval()");
TEST_EXIT_DBG(f)("No function f!\n");
TEST_EXIT_DBG(gradAtQPs)("Operator was not initialized correctly!\n");
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*f)(gradAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
}
void VecAtQP_SOT::eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
......@@ -2047,15 +2063,13 @@ namespace AMDiS {
}
}
void VecAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2079,15 +2093,13 @@ namespace AMDiS {
}
}
void Vec2AtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void Vec2AtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2111,15 +2123,13 @@ namespace AMDiS {
}
}
void CoordsAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void CoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2143,15 +2153,13 @@ namespace AMDiS {
}
}
void FctGradient_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void FctGradient_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2175,15 +2183,13 @@ namespace AMDiS {
}
}
void VecAndGradientAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecAndGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2213,16 +2219,14 @@ namespace AMDiS {
}
}
void VecMatrixGradientAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecMatrixGradientAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
}
......@@ -2246,15 +2250,13 @@ namespace AMDiS {
}
}
void VecAndCoordsAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecAndCoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2284,16 +2286,14 @@ namespace AMDiS {
}
}
void MatrixGradientAndCoords_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void MatrixGradientAndCoords_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
int nPoints = grdUhAtQP.size();
WorldMatrix<double> A;
for (int iq = 0; iq < nPoints; iq++) {
A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
}
......@@ -2317,15 +2317,13 @@ namespace AMDiS {
}
}
void VecGradCoordsAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecGradCoordsAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2489,15 +2487,13 @@ namespace AMDiS {
}
}
void CoordsAtQP_IJ_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void CoordsAtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
result[iq][xi] += grdUhAtQP[iq][xj] * factor;
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*g)(coordsAtQPs[iq]);
result[iq][xi] += grdUhAtQP[iq][xj] * factor;
}
}
......@@ -2533,15 +2529,13 @@ namespace AMDiS {
}
}
void VecAtQP_IJ_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecAtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
result[iq][xi] += grdUhAtQP[iq][xj] * factor;
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq]);
result[iq][xi] += grdUhAtQP[iq][xj] * factor;
}
}
......@@ -2688,15 +2682,13 @@ namespace AMDiS {
}
}
void VecAndGradAtQP_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void VecAndGradAtQP_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
int nPoints = grdUhAtQP.size();
for (int iq = 0; iq < nPoints; iq++) {
double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
......@@ -2791,27 +2783,25 @@ namespace AMDiS {
}
}
void General_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void General_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
int nPoints = grdUhAtQP.size();
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
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++)
gradsArg[i] = gradsAtQPs_[i][iq];
A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
result[iq] += A * grdUhAtQP[iq];
}
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++)
gradsArg[i] = gradsAtQPs_[i][iq];
A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
result[iq] += A * grdUhAtQP[iq];
}
}
......@@ -3015,27 +3005,25 @@ namespace AMDiS {
}
}
void GeneralParametric_SOT::weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
void GeneralParametric_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
std::vector<WorldVector<double> > &result) const
{
int nPoints = grdUhAtQP.size();
int nVecs = static_cast<int>(vecs_.size());
int nGrads = static_cast<int>(grads_.size());
std::vector<double> vecsArg(nVecs);
std::vector<WorldVector<double> > gradsArg(nGrads);
if (grdUhAtQP) {