Commit 306f294f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Removed libpng dependent code. Added HAVE_PNG compiler flag, but not yet done in configure.

parent bbd936aa
#include <algorithm> #include <algorithm>
#include <png.h>
#include <boost/numeric/mtl/mtl.hpp> #include <boost/numeric/mtl/mtl.hpp>
#include "DOFMatrix.h" #include "DOFMatrix.h"
#include "QPsiPhi.h" #include "QPsiPhi.h"
......
...@@ -260,12 +260,14 @@ namespace AMDiS { ...@@ -260,12 +260,14 @@ namespace AMDiS {
#endif #endif
} }
#ifdef HAVE_PNG
if (writePngFormat) { if (writePngFormat) {
PngWriter pngWriter(dataCollectors[0]); PngWriter pngWriter(dataCollectors[0]);
pngWriter.writeFile(fn + ".png", pngType); pngWriter.writeFile(fn + ".png", pngType);
MSG("PNG image file written to %s\n", (fn + ".png").c_str()); MSG("PNG image file written to %s\n", (fn + ".png").c_str());
} }
#endif
for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++) for (int i = 0; i < static_cast<int>(dataCollectors.size()); i++)
delete dataCollectors[i]; delete dataCollectors[i];
......
...@@ -836,7 +836,8 @@ namespace AMDiS { ...@@ -836,7 +836,8 @@ namespace AMDiS {
vec2AtQPs = subAssembler->getVectorAtQPs(vec2, elInfo, quad); 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(); const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
...@@ -1076,6 +1077,75 @@ namespace AMDiS { ...@@ -1076,6 +1077,75 @@ namespace AMDiS {
auxFESpaces.push_back(dv->getFESpace()); 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, General_SOT::General_SOT(std::vector<DOFVectorBase<double>*> vecs,
std::vector<DOFVectorBase<double>*> grads, std::vector<DOFVectorBase<double>*> grads,
TertiaryAbstractFunction<WorldMatrix<double>, TertiaryAbstractFunction<WorldMatrix<double>,
......
...@@ -406,20 +406,14 @@ namespace AMDiS { ...@@ -406,20 +406,14 @@ namespace AMDiS {
AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div, AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
bool sym = false); bool sym = false);
/** \brief /// Implementation of \ref OperatorTerm::initElement().
* Implementation of \ref OperatorTerm::initElement().
*/
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler, void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL); Quadrature *quad = NULL);
/** \brief /// Implements SecondOrderTerm::getLALt().
* Implements SecondOrderTerm::getLALt().
*/
void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const; void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
/** \brief /// Implenetation of SecondOrderTerm::eval().
* Implenetation of SecondOrderTerm::eval().
*/
void eval(int nPoints, void eval(int nPoints,
const double *uhAtQP, const double *uhAtQP,
const WorldVector<double> *grdUhAtQP, const WorldVector<double> *grdUhAtQP,
...@@ -427,37 +421,25 @@ namespace AMDiS { ...@@ -427,37 +421,25 @@ namespace AMDiS {
double *result, double *result,
double factor) const; double factor) const;
/** \brief /// Implenetation of SecondOrderTerm::weakEval().
* Implenetation of SecondOrderTerm::weakEval().
*/
void weakEval(int nPoints, void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP, const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const; WorldVector<double> *result) const;
protected: protected:
/** \brief /// DOFVector to be evaluated at quadrature points.
* DOFVector to be evaluated at quadrature points.
*/
DOFVectorBase<double>* vec; 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; double* vecAtQPs;
/** \brief /// Function for A.
* Function for A.
*/
AbstractFunction<WorldMatrix<double>, double>* matrixFct; AbstractFunction<WorldMatrix<double>, double>* matrixFct;
/** \brief ///
*
*/
AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct; 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; bool symmetric;
}; };
...@@ -1195,6 +1177,55 @@ namespace AMDiS { ...@@ -1195,6 +1177,55 @@ namespace AMDiS {
bool symmetric; 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 class General_SOT : public SecondOrderTerm
{ {
public: public:
...@@ -1209,21 +1240,13 @@ namespace AMDiS { ...@@ -1209,21 +1240,13 @@ namespace AMDiS {
WorldMatrix<double> > *divFct, WorldMatrix<double> > *divFct,
bool symmetric); bool symmetric);
/** \brief /// Implementation of \ref OperatorTerm::initElement().
* Implementation of \ref OperatorTerm::initElement(). void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL);
*/
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; void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
/** \brief /// Implenetation of SecondOrderTerm::eval().
* Implenetation of SecondOrderTerm::eval().
*/
void eval(int nPoints, void eval(int nPoints,
const double *uhAtQP, const double *uhAtQP,
const WorldVector<double> *grdUhAtQP, const WorldVector<double> *grdUhAtQP,
...@@ -1231,9 +1254,7 @@ namespace AMDiS { ...@@ -1231,9 +1254,7 @@ namespace AMDiS {
double *result, double *result,
double factor) const; double factor) const;
/** \brief /// Implenetation of SecondOrderTerm::weakEval().
* Implenetation of SecondOrderTerm::weakEval().
*/
void weakEval(int nPoints, void weakEval(int nPoints,
const WorldVector<double> *grdUhAtQP, const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const; WorldVector<double> *result) const;
......
#ifdef HAVE_PNG
#include <float.h> #include <float.h>
#include <png.h> #include <png.h>
...@@ -122,3 +124,5 @@ namespace AMDiS { ...@@ -122,3 +124,5 @@ namespace AMDiS {
return 0; return 0;
} }
} }
#endif
...@@ -22,6 +22,8 @@ ...@@ -22,6 +22,8 @@
#ifndef AMDIS_PNGWRITER_H #ifndef AMDIS_PNGWRITER_H
#define AMDIS_PNGWRITER_H #define AMDIS_PNGWRITER_H
#ifdef HAVE_PNG
#include "BasisFunction.h" #include "BasisFunction.h"
#include "DataCollector.h" #include "DataCollector.h"
#include "FileWriter.h" #include "FileWriter.h"
...@@ -53,3 +55,5 @@ namespace AMDiS { ...@@ -53,3 +55,5 @@ namespace AMDiS {
} }
#endif #endif
#endif
...@@ -764,7 +764,12 @@ namespace AMDiS { ...@@ -764,7 +764,12 @@ namespace AMDiS {
matrix->getBoundaryManager()->exitMatrix(matrix); matrix->getBoundaryManager()->exitMatrix(matrix);
if (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 // And now assemble boundary conditions on the vectors
......
...@@ -184,9 +184,8 @@ namespace AMDiS { ...@@ -184,9 +184,8 @@ namespace AMDiS {
if (firstCall) { if (firstCall) {
tmpLALt[myRank] = new DimMat<double>*[nPoints]; tmpLALt[myRank] = new DimMat<double>*[nPoints];
for (int j = 0; j < nPoints; j++) { for (int j = 0; j < nPoints; j++)
tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT); tmpLALt[myRank][j] = new DimMat<double>(dim, NO_INIT);
}
const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts(); const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI); psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
...@@ -197,21 +196,17 @@ namespace AMDiS { ...@@ -197,21 +196,17 @@ namespace AMDiS {
DimMat<double> **LALt = tmpLALt[myRank]; DimMat<double> **LALt = tmpLALt[myRank];
for (int i = 0; i < nPoints; i++) { for (int i = 0; i < nPoints; i++)
LALt[i]->set(0.0); 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); (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; VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
if (symmetric) { if (symmetric) {
// === Symmetric assembling. ===
TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n"); TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
...@@ -232,7 +227,8 @@ namespace AMDiS { ...@@ -232,7 +227,8 @@ namespace AMDiS {
} }
} }
} }
} else { /* non symmetric assembling */ } else {
// === Non symmetric assembling. ===
for (int iq = 0; iq < nPoints; iq++) { for (int iq = 0; iq < nPoints; iq++) {
(*LALt[iq]) *= elInfo->getDet(); (*LALt[iq]) *= elInfo->getDet();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment