From 62841d3cb18dc8c3b5849f83278b75f89645a0a0 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Wed, 2 May 2012 13:28:06 +0000
Subject: [PATCH] ElementFileWriter updated and some corrections in DOFVector
---
AMDiS/src/Assembler.cc | 3 +
AMDiS/src/CouplingIterationInterface.cc | 3 +-
AMDiS/src/CouplingProblemStat.h | 1 -
AMDiS/src/DOFVector.h | 7 -
AMDiS/src/DOFVector.hh | 8 +-
AMDiS/src/Debug.cc | 2 +-
AMDiS/src/ElInfo.cc | 6 +
AMDiS/src/Functors.h | 9 ++
AMDiS/src/Initfile.h | 6 +
AMDiS/src/MacroElement.h | 1 +
AMDiS/src/Mesh.cc | 6 +-
AMDiS/src/SecondOrderTerm.cc | 83 ++++++++++
AMDiS/src/SecondOrderTerm.h | 63 +++++++-
AMDiS/src/TransformDOF.h | 35 ++++
AMDiS/src/ZeroOrderTerm.cc | 56 +++++--
AMDiS/src/ZeroOrderTerm.h | 4 +-
AMDiS/src/io/ElementFileWriter.cc | 205 +++++++++++++++++-------
AMDiS/src/io/ElementFileWriter.h | 41 ++++-
18 files changed, 438 insertions(+), 101 deletions(-)
diff --git a/AMDiS/src/Assembler.cc b/AMDiS/src/Assembler.cc
index 6f2143ad..54458041 100644
--- a/AMDiS/src/Assembler.cc
+++ b/AMDiS/src/Assembler.cc
@@ -324,6 +324,9 @@ namespace AMDiS {
calculateElementMatrix(elInfo, elementMatrix);
}
+ // eventuell einfach
+ // vec = elementMatrix*uhOldLoc;
+ // schreiben
for (int i = 0; i < nRow; i++) {
double val = 0.0;
for (int j = 0; j < nCol; j++)
diff --git a/AMDiS/src/CouplingIterationInterface.cc b/AMDiS/src/CouplingIterationInterface.cc
index 717b74ef..d241d26f 100644
--- a/AMDiS/src/CouplingIterationInterface.cc
+++ b/AMDiS/src/CouplingIterationInterface.cc
@@ -43,9 +43,10 @@ namespace AMDiS {
void CouplingIterationInterface::beginIteration(AdaptInfo *adaptInfo)
{ FUNCNAME("CouplingIterationInterface::beginIteration()");
MSG("\n");
+ int nTimesteps = (adaptInfo->getNumberOfTimesteps() ? adaptInfo->getNumberOfTimesteps() : (adaptInfo->getEndTime()-adaptInfo->getStartTime())/adaptInfo->getTimestep());
MSG("begin of iteration number: %d/%d\n",
adaptInfo->getTimestepNumber() + 1,
- adaptInfo->getNumberOfTimesteps());
+ nTimesteps);
MSG("==================================================\n");
};
diff --git a/AMDiS/src/CouplingProblemStat.h b/AMDiS/src/CouplingProblemStat.h
index 6351b3a5..581907f8 100644
--- a/AMDiS/src/CouplingProblemStat.h
+++ b/AMDiS/src/CouplingProblemStat.h
@@ -166,7 +166,6 @@ namespace AMDiS {
if (initMesh && meshes[i])
refinementManager->globalRefine(meshes[i], globalRefinements);
}
-
}
void createRefCoarseManager()
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 8625f12c..4db016c5 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -934,13 +934,6 @@ namespace AMDiS {
TOut integrate_VecAndCoords(const DOFVector<T> &vec,
BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct);
- template<typename TOut, typename T>
- TOut integrate_VecAndCoords(const DOFVector<T> &vec,
- BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
- {
- return integrate(vec, fct);
- }
-
/// Computes the integral: \f$ \int f(\{v_i(\vec{x})\}_i)\,\text{d}\vec{x}\f$
double integrateGeneral(const vector<DOFVector<double>*> &vecs,
AbstractFunction<double, vector<double> > *fct);
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index b4bb8152..ffb06041 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -475,7 +475,7 @@ namespace AMDiS {
const DOFAdmin* admin = this->getFeSpace()->getAdmin();
int nBasFcts = basFct->getNumber();
std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
- mtl::dense_vector<double> fctInterpolValues(nBasFcts);
+ mtl::dense_vector<T> fctInterpolValues(nBasFcts);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
@@ -513,14 +513,14 @@ namespace AMDiS {
else
flag |= Mesh::CALL_EL_LEVEL;
- double result = 0.0;
+ T result; nullify(result);
int nPoints = quadFast->getNumPoints();
mtl::dense_vector<T> uh_vec(nPoints);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, meshLevel, flag);
while (elInfo) {
double det = elInfo->getDet();
- double normT = 0.0;
+ T normT; nullify(normT);
this->getVecAtQPs(elInfo, NULL, quadFast, uh_vec);
for (int iq = 0; iq < nPoints; iq++)
normT += quadFast->getWeight(iq) * (uh_vec[iq]);
@@ -769,7 +769,7 @@ namespace AMDiS {
template<typename TOut, typename T>
TOut integrate_VecAndCoords(const DOFVector<T> &vec,
- BinaryAbstractFunction<T, T, WorldVector<double> > *fct)
+ BinaryAbstractFunction<TOut, T, WorldVector<double> > *fct)
{
FUNCNAME("integrate_VecAndCoords()");
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index 73bfcaf7..752d692f 100644
--- a/AMDiS/src/Debug.cc
+++ b/AMDiS/src/Debug.cc
@@ -106,7 +106,7 @@ namespace AMDiS {
elInfo = stack.traverseNext(elInfo);
}
- ElementFileWriter::writeFile(vec, mesh, filename, "", level);
+ ElementFileWriter::writeFile(vec, mesh, filename, ".vtu", level);
}
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index c9cffb33..6ea21053 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/AMDiS/src/ElInfo.cc
@@ -195,6 +195,11 @@ namespace AMDiS {
void ElInfo::fillElInfo(const MacroElement *mel, int refinementPathLength, unsigned long refinementPath)
{
+ if (refinementPathLength == 0) {
+ fillMacroInfo(mel);
+ return;
+ }
+
std::vector<ElInfo*> elInfo;
elInfo.push_back(mesh->createNewElInfo());
elInfo.push_back(this);
@@ -208,6 +213,7 @@ namespace AMDiS {
}
if (i%2 == 0)
*this = *elInfo[0];
+
delete elInfo[0];
}
diff --git a/AMDiS/src/Functors.h b/AMDiS/src/Functors.h
index ecf20445..fae87d69 100644
--- a/AMDiS/src/Functors.h
+++ b/AMDiS/src/Functors.h
@@ -181,6 +181,15 @@ struct Norm2Sqr_comp3 : public TertiaryAbstractFunction<T,T,T,T>
T operator()(const T &v1, const T &v2, const T &v3) const { return sqr(v1)+sqr(v2)+sqr(v3); }
};
+template<typename T>
+struct Vec1WorldVec : public AbstractFunction<WorldVector<T>,T>
+{
+ WorldVector<T> operator()(const T &v0) const {
+ WorldVector<T> result;
+ result[0]=v0;
+ return result;
+ }
+};
template<typename T>
struct Vec2WorldVec : public BinaryAbstractFunction<WorldVector<T>,T,T>
{
diff --git a/AMDiS/src/Initfile.h b/AMDiS/src/Initfile.h
index 23c9e361..33c0c21e 100644
--- a/AMDiS/src/Initfile.h
+++ b/AMDiS/src/Initfile.h
@@ -438,6 +438,11 @@ namespace AMDiS {
initIntern();
if (debugInfo == -1)
debugInfo = singlett->getMsgInfo();
+ else {
+ int swap(debugInfo);
+ debugInfo = singlett->getMsgInfo();
+ singlett->msgInfo=swap;
+ }
std::string valStr;
int error_code = singlett->checkedGet(tag, valStr);
@@ -456,6 +461,7 @@ namespace AMDiS {
std::cout << "Parameter '" << tag << "'"
<< " initialized with: " << value << std::endl;
}
+ singlett->msgInfo = debugInfo;
}
diff --git a/AMDiS/src/MacroElement.h b/AMDiS/src/MacroElement.h
index aed36d52..ae2f988c 100644
--- a/AMDiS/src/MacroElement.h
+++ b/AMDiS/src/MacroElement.h
@@ -215,6 +215,7 @@ namespace AMDiS {
friend class MacroInfo;
friend class MacroReader;
+ friend class ElInfo;
friend class ElInfo1d;
friend class ElInfo2d;
friend class ElInfo3d;
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index a7f78628..8664e00e 100644
--- a/AMDiS/src/Mesh.cc
+++ b/AMDiS/src/Mesh.cc
@@ -1015,7 +1015,7 @@ namespace AMDiS {
out << name << "\n";
- SerUtil::serialize(out, dim);
+// SerUtil::serialize(out, dim);
SerUtil::serialize(out, nVertices);
SerUtil::serialize(out, nEdges);
SerUtil::serialize(out, nLeaves);
@@ -1095,8 +1095,8 @@ namespace AMDiS {
in.get();
int oldVal = dim;
- SerUtil::deserialize(in, dim);
- TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
+// SerUtil::deserialize(in, dim);
+// TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
SerUtil::deserialize(in, nVertices);
SerUtil::deserialize(in, nEdges);
diff --git a/AMDiS/src/SecondOrderTerm.cc b/AMDiS/src/SecondOrderTerm.cc
index 1af5bd9c..80edf261 100644
--- a/AMDiS/src/SecondOrderTerm.cc
+++ b/AMDiS/src/SecondOrderTerm.cc
@@ -1365,4 +1365,87 @@ namespace AMDiS {
}
+ // ========== Vec2AtQP_IJ_SOT ==========
+
+ Vec2AtQP_IJ_SOT::Vec2AtQP_IJ_SOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af,
+ int x_i, int x_j)
+ : SecondOrderTerm(af ? af->getDegree() : dv1->getFeSpace()->getBasisFcts()->getDegree()+dv2->getFeSpace()->getBasisFcts()->getDegree()),
+ vec1(dv1),
+ vec2(dv2),
+ f(af),
+ xi(x_i), xj(x_j)
+ {
+ setSymmetric(xi == xj);
+
+ TEST_EXIT(dv1)("No vector 1!\n");
+ TEST_EXIT(dv2)("No vector 2!\n");
+
+ auxFeSpaces.insert(dv1->getFeSpace());
+ auxFeSpaces.insert(dv2->getFeSpace());
+ }
+
+ void Vec2AtQP_IJ_SOT::initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad)
+ {
+ getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs1);
+ getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
+ }
+
+ void Vec2AtQP_IJ_SOT::getLALt(const ElInfo *elInfo,
+ vector<mtl::dense2D<double> > &LALt) const
+ {
+ const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
+ const int nPoints = static_cast<int>(LALt.size());
+
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++)
+ lalt_kl(grdLambda, xi, xj, LALt[iq], (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ lalt_kl(grdLambda, xi, xj, LALt[iq], vecAtQPs1[iq] * vecAtQPs2[iq]);
+ }
+ }
+
+ void Vec2AtQP_IJ_SOT::eval(int nPoints,
+ const mtl::dense_vector<double>& uhAtQP,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double fac)
+ {
+ if (num_rows(D2UhAtQP) > 0) {
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ result[iq] += D2UhAtQP[iq][xi][xj] * factor * fac;
+ }
+ } else {
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = vecAtQPs1[iq] * vecAtQPs2[iq];
+ result[iq] += D2UhAtQP[iq][xi][xj] * factor * fac;
+ }
+ }
+ }
+ }
+
+ void Vec2AtQP_IJ_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result)
+ {
+ int nPoints = grdUhAtQP.size();
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ result[iq][xi] += grdUhAtQP[iq][xj] * factor;
+ }
+ } else {
+ for (int iq = 0; iq < nPoints; iq++) {
+ double factor = vecAtQPs1[iq] * vecAtQPs2[iq];
+ result[iq][xi] += grdUhAtQP[iq][xj] * factor;
+ }
+ }
+ }
+
}
diff --git a/AMDiS/src/SecondOrderTerm.h b/AMDiS/src/SecondOrderTerm.h
index a77acc27..b2f67cea 100644
--- a/AMDiS/src/SecondOrderTerm.h
+++ b/AMDiS/src/SecondOrderTerm.h
@@ -361,7 +361,7 @@ namespace AMDiS {
{
public:
/// Constructor.
- VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);
+ VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af = NULL);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
@@ -411,8 +411,9 @@ namespace AMDiS {
{
public:
/// Constructor.
- Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
- BinaryAbstractFunction<double, double, double> *af);
+ Vec2AtQP_SOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af = NULL);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
@@ -1145,7 +1146,8 @@ namespace AMDiS {
{
public:
/// Constructor.
- VecAtQP_IJ_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af,
+ VecAtQP_IJ_SOT(DOFVectorBase<double> *dv,
+ AbstractFunction<double, double> *af,
int x_i, int x_j);
/// Implementation of \ref OperatorTerm::initElement().
@@ -1182,6 +1184,59 @@ namespace AMDiS {
int xi, xj;
};
+ /**
+ * \ingroup Assembler
+ *
+ * \brief
+ * SecondOrderTerm where A is a WorldMatrix having a in all positions
+ * except possibly the position IJ, multiplied with a function
+ * evaluated at the quadrature points of two given DOFVectors:
+ * \f$ \nabla \cdot f(v(\vec{x}),w(\vec{x})) A \nabla u(\vec{x}) \f$
+ */
+ class Vec2AtQP_IJ_SOT : public SecondOrderTerm
+ {
+ public:
+ /// Constructor.
+ Vec2AtQP_IJ_SOT(DOFVectorBase<double> *dv1,
+ DOFVectorBase<double> *dv2,
+ BinaryAbstractFunction<double, double, double> *af,
+ int x_i, int x_j);
+
+ /// Implementation of \ref OperatorTerm::initElement().
+ void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+
+ /// Implements SecondOrderTerm::getLALt().
+ void getLALt(const ElInfo *elInfo, vector<mtl::dense2D<double> > &LALt) const;
+
+ /// Implements SecondOrderTerm::eval().
+ void eval(int nPoints,
+ const mtl::dense_vector<double>& uhAtQP,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double factor);
+
+ /// Implements SecondOrderTerm::weakEval().
+ void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
+ std::vector<WorldVector<double> > &result);
+
+ protected:
+ /// DOFVector to be evaluated at quadrature points.
+ DOFVectorBase<double>* vec1;
+ DOFVectorBase<double>* vec2;
+
+ /// Pointer to an array containing the DOFVector evaluated at quadrature points.
+ mtl::dense_vector<double> vecAtQPs1;
+ mtl::dense_vector<double> vecAtQPs2;
+
+ /// Function evaluated at \ref vecAtQPs.
+ BinaryAbstractFunction<double, double, double> *f;
+
+ private:
+ /// Directions for the derivatives.
+ int xi, xj;
+ };
}
#endif
diff --git a/AMDiS/src/TransformDOF.h b/AMDiS/src/TransformDOF.h
index 43432486..e180186f 100644
--- a/AMDiS/src/TransformDOF.h
+++ b/AMDiS/src/TransformDOF.h
@@ -681,6 +681,7 @@ template<typename T1, typename T2, typename T3, typename T4>
inline void transformDOF(T1 val, DOFVector<T2> &vec2, DOFVector<T3> &vec3, DOFVector<T4> &result, TertiaryAbstractFunction<T4, T1, T2, T3> &tertiary_op)
{ transformDOF(val, &vec2, &vec3, &result, &tertiary_op); }
+
// ===========================================================================================
// return binary_op(vec, interpol(fct))
@@ -716,6 +717,40 @@ template<typename T> inline void transformDOFInterpolation(
DOFVector<T> &vec, AbstractFunction<T, WorldVector<double> > &fct, DOFVector<T> &result, BinaryAbstractFunction<T, T, T> &binary_op)
{ transformDOFInterpolation(&vec, &fct, &result, &binary_op); }
+
+// ====================================================================================
+
+template<typename T>
+T accumulateDOF_simple(DOFVector<T> *vec,
+ T value0,
+ BinaryAbstractFunction<T, T, T> *binary_op)
+{
+ DOFIterator<T> vecIter(vec, USED_DOFS);
+ T value = value0;
+ for(vecIter.reset(); !vecIter.end(); ++vecIter)
+ {
+ value = (*binary_op)(value, *vecIter);
+ }
+ return value;
+}
+
+template<typename TOut, typename T1, typename T2>
+TOut accumulateDOF_simple(DOFVector<T1> *vec1,
+ DOFVector<T2> *vec2,
+ TOut value0,
+ TertiaryAbstractFunction<TOut, TOut, T1, T2> *tertiary_op)
+{
+ TEST_EXIT(vec1->getFeSpace() == vec2->getFeSpace())("FeSpaces must be equal!\n");
+ DOFIterator<T1> vec1Iter(vec1, USED_DOFS);
+ DOFIterator<T2> vec2Iter(vec2, USED_DOFS);
+ TOut value = value0;
+ for(vec1Iter.reset(),vec2Iter.reset(); !vec1Iter.end(); ++vec1Iter,++vec2Iter)
+ {
+ value = (*tertiary_op)(value, *vec1Iter, *vec2Iter);
+ }
+ return value;
+}
+
}
#endif // AMDIS_TRANSFORM_DOF_H
diff --git a/AMDiS/src/ZeroOrderTerm.cc b/AMDiS/src/ZeroOrderTerm.cc
index c324829e..57d0276d 100644
--- a/AMDiS/src/ZeroOrderTerm.cc
+++ b/AMDiS/src/ZeroOrderTerm.cc
@@ -126,7 +126,10 @@ namespace AMDiS {
Vec2AtQP_ZOT::Vec2AtQP_ZOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *af)
- : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
+ : ZeroOrderTerm((af ? af->getDegree() : dv1->getFeSpace()->getBasisFcts()->getDegree()+dv2->getFeSpace()->getBasisFcts()->getDegree())),
+ vec1(dv1),
+ vec2(dv2),
+ f(af)
{
TEST_EXIT(dv1)("No first vector!\n");
TEST_EXIT(dv2)("No second vector!\n");
@@ -156,9 +159,14 @@ namespace AMDiS {
}
void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
- {
- for (int iq = 0; iq < nPoints; iq++)
- C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ {
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += vecAtQPs1[iq] * vecAtQPs2[iq];
+ }
}
void Vec2AtQP_ZOT::eval(int nPoints,
@@ -168,8 +176,13 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double fac)
{
- for (int iq = 0; iq < nPoints; iq++)
- result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq];
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq];
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] += fac * vecAtQPs1[iq] * vecAtQPs2[iq] * uhAtQP[iq];
+ }
}
@@ -179,7 +192,13 @@ namespace AMDiS {
DOFVectorBase<double> *dv2,
DOFVectorBase<double> *dv3,
TertiaryAbstractFunction<double, double, double, double> *af)
- : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), vec3(dv3), f(af)
+ : ZeroOrderTerm((af ? af->getDegree() : dv1->getFeSpace()->getBasisFcts()->getDegree()
+ + dv2->getFeSpace()->getBasisFcts()->getDegree()
+ + dv3->getFeSpace()->getBasisFcts()->getDegree())),
+ vec1(dv1),
+ vec2(dv2),
+ vec3(dv3),
+ f(af)
{
TEST_EXIT(dv1)("No first vector!\n");
TEST_EXIT(dv2)("No second vector!\n");
@@ -200,9 +219,14 @@ namespace AMDiS {
}
void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
- {
- for (int iq = 0; iq < nPoints; iq++)
- C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
+ {
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ C[iq] += vecAtQPs1[iq] * vecAtQPs2[iq] * vecAtQPs3[iq];
+ }
}
void Vec3AtQP_ZOT::eval(int nPoints,
@@ -212,9 +236,15 @@ namespace AMDiS {
mtl::dense_vector<double>& result,
double fac)
{
- for (int iq = 0; iq < nPoints; iq++)
- result[iq] +=
- fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * uhAtQP[iq];
+ if (f) {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] +=
+ fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * uhAtQP[iq];
+ } else {
+ for (int iq = 0; iq < nPoints; iq++)
+ result[iq] +=
+ fac * vecAtQPs1[iq] * vecAtQPs2[iq] * vecAtQPs3[iq] * uhAtQP[iq];
+ }
}
diff --git a/AMDiS/src/ZeroOrderTerm.h b/AMDiS/src/ZeroOrderTerm.h
index 21231e8e..986cfd53 100644
--- a/AMDiS/src/ZeroOrderTerm.h
+++ b/AMDiS/src/ZeroOrderTerm.h
@@ -190,7 +190,7 @@ namespace AMDiS {
/// Constructor.
Vec2AtQP_ZOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
- BinaryAbstractFunction<double, double, double> *f);
+ BinaryAbstractFunction<double, double, double> *f = NULL);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
@@ -242,7 +242,7 @@ namespace AMDiS {
Vec3AtQP_ZOT(DOFVectorBase<double> *dv1,
DOFVectorBase<double> *dv2,
DOFVectorBase<double> *dv3,
- TertiaryAbstractFunction<double, double, double, double> *f);
+ TertiaryAbstractFunction<double, double, double, double> *f = NULL);
/// Implementation of \ref OperatorTerm::initElement().
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
diff --git a/AMDiS/src/io/ElementFileWriter.cc b/AMDiS/src/io/ElementFileWriter.cc
index 6878a05f..e61d8800 100644
--- a/AMDiS/src/io/ElementFileWriter.cc
+++ b/AMDiS/src/io/ElementFileWriter.cc
@@ -33,21 +33,27 @@ namespace AMDiS {
: name(name_),
amdisMeshDatExt(".elem.mesh"),
vtkExt(".vtu"),
+ pvdExt(".pvd"),
writeAMDiSFormat(0),
writeVtkFormat(0),
+ writeVtkVectorFormat(0),
+ writeAs3dVector(true),
+ writeParaViewAnimation(0),
appendIndex(0),
indexLength(5),
indexDecimals(3),
tsModulo(1),
timestepNumber(-1),
mesh(mesh_),
- vec(mapvec)
+ vec(&mapvec),
+ vecs(NULL)
{
if (name != "") {
Parameters::get(name + "->output->filename", filename);
Parameters::get(name + "->output->AMDiS format", writeAMDiSFormat);
Parameters::get(name + "->output->AMDiS mesh-dat ext", amdisMeshDatExt);
Parameters::get(name + "->output->ParaView format", writeVtkFormat);
+ Parameters::get(name + "->output->ParaView animation", writeParaViewAnimation);
Parameters::get(name + "->output->append index", appendIndex);
Parameters::get(name + "->output->index length", indexLength);
Parameters::get(name + "->output->index decimals", indexDecimals);
@@ -56,6 +62,43 @@ namespace AMDiS {
}
+ ElementFileWriter::ElementFileWriter(string name_,
+ Mesh *mesh_,
+ map<int, vector<double> > &mapvec)
+ : name(name_),
+ amdisMeshDatExt(".elem.mesh"),
+ vtkExt(".vtu"),
+ pvdExt(".pvd"),
+ writeAMDiSFormat(0),
+ writeVtkFormat(0),
+ writeVtkVectorFormat(0),
+ writeAs3dVector(true),
+ writeParaViewAnimation(0),
+ appendIndex(0),
+ indexLength(5),
+ indexDecimals(3),
+ tsModulo(1),
+ timestepNumber(-1),
+ mesh(mesh_),
+ vec(NULL),
+ vecs(&mapvec)
+ {
+ if (name != "") {
+ Parameters::get(name + "->output->filename", filename);
+ Parameters::get(name + "->output->AMDiS format", writeAMDiSFormat);
+ Parameters::get(name + "->output->AMDiS mesh-dat ext", amdisMeshDatExt);
+ Parameters::get(name + "->output->ParaView format", writeVtkFormat);
+ Parameters::get(name + "->output->ParaView vector format", writeVtkVectorFormat);
+ Parameters::get(name + "->output->write vector as 3d vector", writeAs3dVector);
+ Parameters::get(name + "->output->ParaView animation", writeParaViewAnimation);
+ Parameters::get(name + "->output->append index", appendIndex);
+ Parameters::get(name + "->output->index length", indexLength);
+ Parameters::get(name + "->output->index decimals", indexDecimals);
+ Parameters::get(name + "->output->write every i-th timestep", tsModulo);
+ }
+ }
+
+
void ElementFileWriter::writeFiles(AdaptInfo *adaptInfo, bool force,
int level, Flag traverseFlag,
bool (*writeElem)(ElInfo*))
@@ -92,12 +135,16 @@ namespace AMDiS {
MSG("MeshDat file written to %s\n", (fn + amdisMeshDatExt).c_str());
}
- if (writeVtkFormat) {
+ if (writeVtkFormat || writeVtkVectorFormat) {
TEST_EXIT(mesh)("no mesh\n");
- writeVtkValues(fn, vtkExt);
+ writeVtkValues(fn, vtkExt,-1,writeVtkVectorFormat==1);
MSG("VTK file written to %s\n", (fn + vtkExt).c_str());
}
+
+ if (writeParaViewAnimation) {
+ VtkWriter::updateAnimationFile(adaptInfo, (fn + vtkExt), ¶ViewAnimationFrames, (filename + pvdExt));
+ }
}
@@ -105,10 +152,23 @@ namespace AMDiS {
Mesh *mesh,
string filename,
string postfix,
- int level)
+ int level,
+ bool writeAsVector)
{
ElementFileWriter efw("", mesh, vec);
- efw.writeVtkValues(filename, postfix, level);
+ efw.writeVtkValues(filename, postfix, level, writeAsVector);
+ }
+
+
+ void ElementFileWriter::writeFile(map<int, vector<double> > &vecs,
+ Mesh *mesh,
+ string filename,
+ string postfix,
+ int level,
+ bool writeAsVector)
+ {
+ ElementFileWriter efw("", mesh, vecs);
+ efw.writeVtkValues(filename, postfix, level, writeAsVector);
}
@@ -198,7 +258,7 @@ namespace AMDiS {
while (elInfo) {
// Get element value.
- val = vec[elInfo->getElement()->getIndex()];
+ val = (*vec)[elInfo->getElement()->getIndex()];
// Write value for each vertex of each element.
for (int i = 0; i <= dim; ++i) {
@@ -218,10 +278,13 @@ namespace AMDiS {
void ElementFileWriter::writeVtkValues(string fname,
string postfix,
- int level)
+ int level,
+ bool writeAsVector)
{
FUNCNAME("ElementFileWriter::writeVtkValues()");
+ TEST_EXIT((vec!=NULL || vecs!=NULL) && (vec==NULL || vecs==NULL))("Ether vec or vecs must be given, not both and not nothing!");
+
#if HAVE_PARALLEL_DOMAIN_AMDIS
string filename =
fname +
@@ -242,49 +305,47 @@ namespace AMDiS {
file.push(boost::iostreams::file_descriptor_sink(filename, ios::trunc));
int dim = mesh->getDim();
- int dimOfWorld = Global::getGeo(WORLD);
int vertices = mesh->getGeo(VERTEX);
- int nElements = mesh->getNumberOfLeaves();
+ int nVertices = 0;
+ int nElements = 0;
Flag traverseFlag = level == -1 ? Mesh::CALL_LEAF_EL : Mesh::CALL_EL_LEVEL;
- if (level != -1) {
- nElements = 0;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh, level, traverseFlag);
+ while (elInfo) {
+ nElements++;
+ elInfo = stack.traverseNext(elInfo);
+ }
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, level, Mesh::CALL_EL_LEVEL);
- while (elInfo) {
- nElements++;
- elInfo = stack.traverseNext(elInfo);
+ // === Write vertex coordinates (every vertex for every element). ===
+ elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
+
+ std::map<DegreeOfFreedom, std::pair<int, WorldVector<double> > > dof2Idx;
+ std::vector<DegreeOfFreedom> idx2Dof;
+ while (elInfo) {
+ // Write coordinates of all element vertices.
+ for (int i = 0; i <= dim; i++) {
+ DegreeOfFreedom idx = elInfo->getElement()->getDof(i, 0);
+ if (dof2Idx.find(idx) == dof2Idx.end()) {
+ dof2Idx[idx] = std::make_pair(nVertices, elInfo->getCoord(i));
+ idx2Dof.push_back(idx);
+ nVertices++;
+ }
}
+ elInfo = stack.traverseNext(elInfo);
}
file << "<?xml version=\"1.0\"?>\n";
file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
file << " <UnstructuredGrid>\n";
- file << " <Piece NumberOfPoints=\"" << (dim + 1) * nElements
+ file << " <Piece NumberOfPoints=\"" << nVertices
<< "\" NumberOfCells=\"" << nElements << "\">\n";
file << " <Points>\n";
file << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n";
-
- // === Write vertex coordinates (every vertex for every element). ===
- TraverseStack stack;
- ElInfo *elInfo =
- stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
-
- while (elInfo) {
- // Write coordinates of all element vertices.
- for (int i = 0; i <= dim; i++) {
- for (int j = 0; j < dimOfWorld; j++)
- file << elInfo->getCoord(i)[j] << " ";
-
- for (int j = dimOfWorld; j < 3; j++)
- file << "0.0";
-
- file << "\n";
- }
-
- elInfo = stack.traverseNext(elInfo);
+ for (int i = 0; i < nVertices; i++) {
+ DegreeOfFreedom dof = idx2Dof[i];
+ writeCoord(file, dof2Idx[dof].second);
}
file << " </DataArray>\n";
@@ -316,40 +377,60 @@ namespace AMDiS {
file << " </DataArray>\n";
file << " <DataArray type=\"Int32\" Name=\"connectivity\">\n";
- int vertCntr = 0;
- for (int i = 0; i < nElements; ++i) {
- for (int j = 0; j <= dim; ++j) {
- file << vertCntr << " ";
- ++vertCntr;
+ elInfo = stack.traverseFirst(mesh, level, traverseFlag);
+ while (elInfo) {
+ // Write coordinates of all element vertices.
+ for (int i = 0; i <= dim; i++) {
+ DegreeOfFreedom dof = elInfo->getElement()->getDof(i, 0);
+ file << (dof2Idx[dof].first) << " ";
}
file << "\n";
+ elInfo = stack.traverseNext(elInfo);
}
file << " </DataArray>\n";
-
file << " </Cells>\n";
- file << " <PointData>\n";
- file << " <DataArray type=\"Float32\" Name=\"value\" format=\"ascii\">\n";
- file.setf(ios::scientific,ios::floatfield);
+ int dataLength = (vecs != NULL ? (*(vecs->begin())).second.size() : 1);
+ int nComponents = (!writeAsVector || (vecs == NULL && vec != NULL) ? 1 : dataLength);
+ int nDataArrays = (!writeAsVector && (vec == NULL && vecs != NULL) ? dataLength : 1);
+ file << " <CellData>\n";
+ for (int i = 0; i < nDataArrays; i++) {
+ file << " <DataArray type=\"Float32\" Name=\"value"<<i<<"\" format=\"ascii\" NumberOfComponents=\""<<(writeAsVector ? std::max(3,nComponents) : nComponents)<<"\">\n";
- elInfo = stack.traverseFirst(mesh, level, traverseFlag | Mesh::FILL_COORDS);
- int vc = 0;
- while (elInfo) {
- // Get element value.
- double val = vec[elInfo->getElement()->getIndex()];
-
- // Write value for each vertex of each element.
- for (int i = 0; i <= dim; i++)
- file << (fabs(val) < 1e-40 ? 0.0 : val) << "\n";
-
- vc++;
- elInfo = stack.traverseNext(elInfo);
- }
+ file.setf(ios::scientific,ios::floatfield);
- file << " </DataArray>\n";
-
+ elInfo = stack.traverseFirst(mesh, level, traverseFlag);
+ int vc = 0;
+ while (elInfo) {
+ // Get element value.
+ int idx = elInfo->getElement()->getIndex();
+
+ for (int j = 0; j < nComponents; j++) {
+ double val = (vec != NULL ? (*vec)[idx] : ((*vecs)[idx].size()==dataLength ? (*vecs)[idx][i*nComponents+j] : 0.0));
+
+ // Write value for each vertex of each element.
+ if (fabs(val) < 1.e-40)
+ file << "0.0";
+ else
+ file << val;
- file << " </PointData>\n";
+ if (j < nComponents-1)
+ file << " ";
+ }
+ if (writeAs3dVector && writeAsVector && vecs != NULL) {
+ for (int j = nComponents; j < 3; j++)
+ file << " 0.0";
+ }
+ file << "\n";
+
+ vc++;
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ file << " </DataArray>\n";
+ }
+
+ file << " </CellData>\n";
file << " </Piece>\n";
file << " </UnstructuredGrid>\n";
file << "</VTKFile>\n";
@@ -361,3 +442,5 @@ namespace AMDiS {
#endif
}
}
+
+
diff --git a/AMDiS/src/io/ElementFileWriter.h b/AMDiS/src/io/ElementFileWriter.h
index e8915b1a..e678ab8d 100644
--- a/AMDiS/src/io/ElementFileWriter.h
+++ b/AMDiS/src/io/ElementFileWriter.h
@@ -48,6 +48,10 @@ namespace AMDiS {
ElementFileWriter(string name,
Mesh *mesh,
map<int, double> &vec);
+
+ ElementFileWriter(string name,
+ Mesh *mesh,
+ map<int, vector<double> > &vecs);
/// Implementation of FileWriterInterface::writeFiles().
void writeFiles(AdaptInfo *adaptInfo, bool force,
@@ -60,7 +64,15 @@ namespace AMDiS {
Mesh *mesh,
string filename,
string postfix = ".vtu",
- int level = -1);
+ int level = -1,
+ bool writeAsVector = false);
+
+ static void writeFile(map<int, vector<double> > &vecs,
+ Mesh *mesh,
+ string filename,
+ string postfix = ".vtu",
+ int level = -1,
+ bool writeAsVector = true);
protected:
/// Writes element data in AMDiS format (1 file !).
@@ -68,8 +80,20 @@ namespace AMDiS {
/// Writes element data in VTK format.
void writeVtkValues(string filename, string postfix,
- int level = -1);
-
+ int level = -1, bool writeAsVector = false);
+
+ /// Writes a world coordinate to a given file.
+ template<typename T>
+ void writeCoord(T &file, WorldVector<double> coord)
+ {
+ for (int i = 0; i < Global::getGeo(WORLD); i++)
+ file << " " << std::scientific << coord[i];
+ for (int i = Global::getGeo(WORLD); i < 3; i++)
+ file << " " << "0.0";
+
+ file << "\n";
+ }
+
protected:
/// Name.
string name;
@@ -82,6 +106,7 @@ namespace AMDiS {
/// VTK file extension.
string vtkExt;
+ string pvdExt;
/// 0: Don't write AMDiS files.
/// 1: Write AMDiS files.
@@ -90,7 +115,12 @@ namespace AMDiS {
/// 0: Don't write VTK files.
/// 1: Write VTK files.
int writeVtkFormat;
+ int writeVtkVectorFormat;
+ bool writeAs3dVector;
+
+ int writeParaViewAnimation;
+
/// 0: Don't append time index to filename prefix.
/// 1: Append time index to filename prefix.
int appendIndex;
@@ -111,7 +141,10 @@ namespace AMDiS {
Mesh *mesh;
/// Vector that stores the solution.
- map<int, double> vec;
+ map<int, double> *vec;
+ map<int, vector<double> > *vecs;
+
+ vector< pair<double, string> > paraViewAnimationFrames;
};
}
--
GitLab