diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index f99efee188c8350adca5a54354ce42a38be6a574..7ea536218d5107c65e2268623a8950d89c9a4765 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -1,5 +1,4 @@
#include <algorithm>
-#include <png.h>
#include <boost/numeric/mtl/mtl.hpp>
#include "DOFMatrix.h"
#include "QPsiPhi.h"
diff --git a/AMDiS/src/FileWriter.cc b/AMDiS/src/FileWriter.cc
index 26af7361e01088986fc385e1148db32fad72ce60..d6bf9a7224a7f0265c9352309e83f71c8aaa905d 100644
--- a/AMDiS/src/FileWriter.cc
+++ b/AMDiS/src/FileWriter.cc
@@ -260,12 +260,14 @@ namespace AMDiS {
#endif
}
+#ifdef HAVE_PNG
if (writePngFormat) {
PngWriter pngWriter(dataCollectors[0]);
pngWriter.writeFile(fn + ".png", pngType);
MSG("PNG image file written to %s\n", (fn + ".png").c_str());
}
+#endif
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++)
delete dataCollectors[i];
diff --git a/AMDiS/src/Operator.cc b/AMDiS/src/Operator.cc
index fb8d01e51bb48ed3b6614284da1ee0be02e86add..af4974cc7ded415376b314f3392ea8059d67bb9f 100644
--- a/AMDiS/src/Operator.cc
+++ b/AMDiS/src/Operator.cc
@@ -836,7 +836,8 @@ namespace AMDiS {
vec2AtQPs = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
}
- void Vec2AtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
+ void Vec2AtQP_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++) {
@@ -1076,6 +1077,75 @@ namespace AMDiS {
auxFESpaces.push_back(dv->getFESpace());
}
+
+ MatrixVec2_SOT::MatrixVec2_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *f,
+ WorldMatrix<double> Af,
+ bool sym)
+ : SecondOrderTerm(f->getDegree()),
+ vec1(dv1),
+ vec2(dv2),
+ fct(f),
+ A(Af),
+ symmetric(sym)
+ {
+ FUNCNAME("MatrixVec2_SOT::MatrixVec2_SOT()");
+
+ setSymmetric(symmetric);
+
+ TEST_EXIT_DBG(dv1)("No vector!\n");
+ TEST_EXIT_DBG(dv2)("No vector!\n");
+
+ auxFESpaces.push_back(dv1->getFESpace());
+ auxFESpaces.push_back(dv2->getFESpace());
+ }
+
+ void MatrixVec2_SOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ vec1AtQPs = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
+ vec2AtQPs = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
+ }
+
+ void MatrixVec2_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, A, *(LALt[iq]), symmetric, (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]));
+ }
+
+ void MatrixVec2_SOT::eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const
+ {
+ int dow = Global::getGeo(WORLD);
+
+ for (int iq = 0; iq < nPoints; iq++) {
+ double resultQP = 0.0;
+
+ if (D2UhAtQP)
+ for (int i = 0; i < dow; i++)
+ for (int j = 0; j < dow; j++)
+ resultQP += A[i][j] * D2UhAtQP[iq][j][i];
+
+ result[iq] += resultQP * factor * (*fct)(vec1AtQPs[iq], vec2AtQPs[iq]);
+ }
+ }
+
+ void MatrixVec2_SOT::weakEval(int nPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const
+ {
+ if (grdUhAtQP)
+ 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,
std::vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<WorldMatrix<double>,
diff --git a/AMDiS/src/Operator.h b/AMDiS/src/Operator.h
index 6065b87a3c65f00bd05616602a027a565b1af8bf..b456bd3175dd63803fe6909297c22b02cf3f2faf 100644
--- a/AMDiS/src/Operator.h
+++ b/AMDiS/src/Operator.h
@@ -406,20 +406,14 @@ namespace AMDiS {
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
bool sym = false);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
+ /// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -427,37 +421,25 @@ 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
- * DOFVector to be evaluated at quadrature points.
- */
+ /// DOFVector to be evaluated at quadrature points.
DOFVectorBase<double>* vec;
- /** \brief
- * Pointer to the values of the DOFVector at quadrature points.
- */
+ /// Pointer to the values of the DOFVector at quadrature points.
double* vecAtQPs;
- /** \brief
- * Function for A.
- */
+ /// Function for A.
AbstractFunction<WorldMatrix<double>, double>* matrixFct;
- /** \brief
- *
- */
+ ///
AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct;
- /** \brief
- * True, if \ref matrixFct produces always symmetric matrices.
- */
+ /// True, if \ref matrixFct produces always symmetric matrices.
bool symmetric;
};
@@ -1195,6 +1177,55 @@ namespace AMDiS {
bool symmetric;
};
+ class MatrixVec2_SOT : public SecondOrderTerm
+ {
+ public:
+ /// Constructor.
+ MatrixVec2_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *f,
+ WorldMatrix<double> Af,
+ bool sym = false);
+
+ /// Implementation of \ref OperatorTerm::initElement().
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /// Implements SecondOrderTerm::getLALt().
+ void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
+
+ /// Implenetation of SecondOrderTerm::eval().
+ void eval(int nPoints,
+ const double *uhAtQP,
+ const WorldVector<double> *grdUhAtQP,
+ const WorldMatrix<double> *D2UhAtQP,
+ double *result,
+ double factor) const;
+
+ /// Implenetation of SecondOrderTerm::weakEval().
+ void weakEval(int nPoints,
+ const WorldVector<double> *grdUhAtQP,
+ WorldVector<double> *result) const;
+
+ protected:
+ /// DOFVector to be evaluated at quadrature points.
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ /// Pointer to the values of the DOFVector at quadrature points.
+ double* vec1AtQPs;
+ double* vec2AtQPs;
+
+ /// Function for A.
+ BinaryAbstractFunction<double, double, double> * fct;
+
+ ///
+ WorldMatrix<double> A;
+
+ /// True, if \ref matrixFct produces always symmetric matrices.
+ bool symmetric;
+ };
+
+
class General_SOT : public SecondOrderTerm
{
public:
@@ -1209,21 +1240,13 @@ namespace AMDiS {
WorldMatrix<double> > *divFct,
bool symmetric);
- /** \brief
- * Implementation of \ref OperatorTerm::initElement().
- */
- void initElement(const ElInfo*,
- SubAssembler* ,
- Quadrature *quad= NULL);
+ /// Implementation of \ref OperatorTerm::initElement().
+ void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL);
- /** \brief
- * Implements SecondOrderTerm::getLALt().
- */
+ /// Implements SecondOrderTerm::getLALt().
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
- /** \brief
- * Implenetation of SecondOrderTerm::eval().
- */
+ /// Implenetation of SecondOrderTerm::eval().
void eval(int nPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
@@ -1231,9 +1254,7 @@ 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;
diff --git a/AMDiS/src/PngWriter.cc b/AMDiS/src/PngWriter.cc
index 9cd553e59d803358c11c3931fa425102a9a093d6..63659beda242b74901b8f4c88a42e590b15c5839 100644
--- a/AMDiS/src/PngWriter.cc
+++ b/AMDiS/src/PngWriter.cc
@@ -1,3 +1,5 @@
+#ifdef HAVE_PNG
+
#include <float.h>
#include <png.h>
@@ -122,3 +124,5 @@ namespace AMDiS {
return 0;
}
}
+
+#endif
diff --git a/AMDiS/src/PngWriter.h b/AMDiS/src/PngWriter.h
index c4934e5301eebed205bc921e8fde0ca5dba7d97c..81c10446b0b40e7f969073c7dcc5e1675c8108f6 100644
--- a/AMDiS/src/PngWriter.h
+++ b/AMDiS/src/PngWriter.h
@@ -22,6 +22,8 @@
#ifndef AMDIS_PNGWRITER_H
#define AMDIS_PNGWRITER_H
+#ifdef HAVE_PNG
+
#include "BasisFunction.h"
#include "DataCollector.h"
#include "FileWriter.h"
@@ -53,3 +55,5 @@ namespace AMDiS {
}
#endif
+
+#endif
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 8109c91f13ab91fdda7ed24a7557a23ed8e3f068..b1847b6877e106e0507a2b7084c545e6c09e6315 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -764,7 +764,12 @@ namespace AMDiS {
matrix->getBoundaryManager()->exitMatrix(matrix);
if (matrix)
- nnz += matrix->getBaseMatrix().nnz();
+ nnz += matrix->getBaseMatrix().nnz();
+
+ if (matrix) {
+ std::cout << "i: " << i << " " << j << std::endl;
+ print(matrix->getBaseMatrix());
+ }
}
// And now assemble boundary conditions on the vectors
diff --git a/AMDiS/src/SecondOrderAssembler.cc b/AMDiS/src/SecondOrderAssembler.cc
index acd1e4d71ff6c889d6b275136ee62eaa7682c117..61e1d6481dd8681c77714672827823a4bf3982ef 100644
--- a/AMDiS/src/SecondOrderAssembler.cc
+++ b/AMDiS/src/SecondOrderAssembler.cc
@@ -184,9 +184,8 @@ namespace AMDiS {
if (firstCall) {
tmpLALt[myRank] = new DimMat<double>*[nPoints];
- for (int j = 0; j < nPoints; j++) {
- tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
- }
+ for (int j = 0; j < nPoints; j++)
+ tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
@@ -197,21 +196,17 @@ namespace AMDiS {
DimMat<double> **LALt = tmpLALt[myRank];
- for (int i = 0; i < nPoints; i++) {
+ for (int i = 0; i < nPoints; i++)
LALt[i]->set(0.0);
- }
-
- for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
+
+ for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
(static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
- }
-// for (int i = 0; i < nPoints; i++) {
-// LALt[i]->print();
-// }
-
VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
if (symmetric) {
+ // === Symmetric assembling. ===
+
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) {
@@ -232,7 +227,8 @@ namespace AMDiS {
}
}
}
- } else { /* non symmetric assembling */
+ } else {
+ // === Non symmetric assembling. ===
for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet();