Commit 781a6094 authored by Praetorius, Simon's avatar Praetorius, Simon

some optimizations in dune-functions discretization

parent 301640ea
...@@ -14,6 +14,9 @@ ...@@ -14,6 +14,9 @@
namespace AMDiS namespace AMDiS
{ {
const Flag EQUAL_BASES = 0x01L;
const Flag EQUAL_LOCALBASES = 0x02L;
template <class GlobalBasis> template <class GlobalBasis>
class Assembler class Assembler
{ {
...@@ -87,11 +90,28 @@ namespace AMDiS ...@@ -87,11 +90,28 @@ namespace AMDiS
return globalBasis_.gridView(); return globalBasis_.gridView();
} }
public:
template <class RowNode, class ColNode>
Flag optimizationFlags(RowNode const& rowNode, ColNode const& colNode) const
{
Flag flag;
if (rowNode.treeIndex() == colNode.treeIndex())
flag |= EQUAL_BASES;
// NOTE: find a way to compare localBases directly
// if (rowNode.finiteElement().localBasis().order() == colNode.finiteElement().localBasis().order())
// flag |= EQUAL_LOCALBASES;
return flag;
}
private: private:
GlobalBasis& globalBasis_; GlobalBasis& globalBasis_;
MatrixOperators<GlobalBasis>& matrixOperators_; MatrixOperators<GlobalBasis>& matrixOperators_;
VectorOperators<GlobalBasis>& rhsOperators_; VectorOperators<GlobalBasis>& rhsOperators_;
Constraints<GlobalBasis>& constraints_; Constraints<GlobalBasis>& constraints_;
//TODO: add caching of localBases
}; };
} // end namespace AMDiS } // end namespace AMDiS
......
...@@ -24,14 +24,14 @@ void Assembler<GlobalBasis>::assemble( ...@@ -24,14 +24,14 @@ void Assembler<GlobalBasis>::assemble(
std::size_t localSize = localView.maxSize(); std::size_t localSize = localView.maxSize();
mtl::dense2D<double> elementMatrix(localSize, localSize); mtl::dense2D<double> elementMatrix(localSize, localSize);
set_to_zero(elementMatrix);
mtl::dense_vector<double> elementVector(localSize); mtl::dense_vector<double> elementVector(localSize);
set_to_zero(elementVector);
// 3. traverse grid and assemble operators on the elements // 3. traverse grid and assemble operators on the elements
for (auto const& element : elements(globalBasis_.gridView())) for (auto const& element : elements(globalBasis_.gridView()))
{ {
set_to_zero(elementMatrix);
set_to_zero(elementVector);
localView.bind(element); localView.bind(element);
// traverse type-tree of global-basis // traverse type-tree of global-basis
...@@ -53,7 +53,8 @@ void Assembler<GlobalBasis>::assemble( ...@@ -53,7 +53,8 @@ void Assembler<GlobalBasis>::assemble(
auto colLocalView = colBasis.localView(); auto colLocalView = colBasis.localView();
colLocalView.bind(element); // NOTE: Is this necessary colLocalView.bind(element); // NOTE: Is this necessary
assembleElementOperators(elementMatrix, matrix, matOp, rowLocalView, colLocalView); assembleElementOperators(elementMatrix, matrix, matOp, rowLocalView, colLocalView,
optimizationFlags(rowNode, colNode));
} }
}); });
}); });
...@@ -105,7 +106,7 @@ void Assembler<GlobalBasis>::assembleElementOperators( ...@@ -105,7 +106,7 @@ void Assembler<GlobalBasis>::assembleElementOperators(
auto assemble_operators = [&](auto const& context, auto& operator_list) { auto assemble_operators = [&](auto const& context, auto& operator_list) {
for (auto scaled : operator_list) { for (auto scaled : operator_list) {
bool add_op = scaled.op->assemble(context, localViews..., elementContainer, scaled.factor); bool add_op = scaled.op->assemble(context, elementContainer, localViews..., scaled.factor);
add = add || add_op; add = add || add_op;
} }
}; };
......
#pragma once #pragma once
#include <string> #include <string>
#include <vector> //#include <vector>
#include <array> //#include <array>
#include <memory> #include <memory>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh> #include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/io/file/vtk/vtksequencewriter.hh> #include <dune/grid/io/file/vtk/vtksequencewriter.hh>
#include <dune/geometry/referenceelements.hh> //#include <dune/geometry/referenceelements.hh>
#include <dune/typetree/childextraction.hh>
#include <dune/amdis/Initfile.hpp> #include <dune/amdis/Initfile.hpp>
#include <dune/amdis/common/Loops.hpp> #include <dune/amdis/io/FileWriterInterface.hpp>
#include <dune/amdis/utility/Filesystem.hpp> #include <dune/amdis/utility/Filesystem.hpp>
namespace AMDiS namespace AMDiS
{ {
template <class GlobalBasis, class TreePath>
template <class Traits>
class FileWriter class FileWriter
: public FileWriterInterface
{ {
private: // typedefs and static constants private: // typedefs and static constants
using MeshView = typename Traits::GridView; using GridView = typename GlobalBasis::GridView;
/// Dimension of the mesh /// Dimension of the mesh
static constexpr int dim = MeshView::dimension; static constexpr int dim = GridView::dimension;
/// Dimension of the world /// Dimension of the world
static constexpr int dow = MeshView::dimensionworld; static constexpr int dow = GridView::dimensionworld;
/// Number of problem components using Vector = DOFVector<GlobalBasis>;
static constexpr int nComponents = 1;
public:
public:
/// Constructor. /// Constructor.
FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names) FileWriter(std::string baseName,
: meshView(meshView) std::shared_ptr<GlobalBasis> const& basis,
, names(names) TreePath const& tp,
std::shared_ptr<Vector> const& vector)
: FileWriterInterface(baseName)
, basis_(basis)
, treePath_(tp)
, vector_(vector)
{ {
filename = "solution"; //int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0);
Parameters::get(base + "->filename", filename);
dir = "output";
Parameters::get(base + "->output directory", dir);
if (!filesystem::exists(dir)) //if (subSampling > 0)
error_exit("Output directory '",dir,"' does not exist!"); // vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(basis_->gridView(), subSampling);
//else
vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(basis_->gridView());
vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView); vtkSeqWriter_ = std::make_shared<Dune::VTKSequenceWriter<GridView>>(vtkWriter_, filename_, dir_, "");
vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
mode_ = Parameters::get<int>(baseName + "->ParaView mode").value_or(0);
} }
/// Implements \ref FileWriterInterface::writeFiles
virtual void writeFiles(AdaptInfo& adaptInfo, bool force) override
{
using Tree = typename GlobalBasis::LocalView::Tree;
using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
writeVertexData(typename Node::NodeTag{}, index_<Node::CHILDREN>, [&,this]()
{
if (mode_ == 0)
vtkSeqWriter_->write(adaptInfo.getTime(), Dune::VTK::ascii);
else if (mode_ == 1)
vtkSeqWriter_->write(adaptInfo.getTime(), Dune::VTK::appendedraw);
});
}
protected:
template <class W>
void writeVertexData(Dune::TypeTree::LeafNodeTag, index_t<0>, W write)
{
using namespace Dune::Functions::BasisBuilder;
auto fct = makeDiscreteFunction(basis_,treePath_,vector_);
auto p1basis = makeBasis(basis_->gridView(), lagrange<1>());
auto data = makeDOFVector(p1basis, name_);
interpolate(p1basis, data, *fct);
auto dataFct = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(p1basis,data);
vtkWriter_->addVertexData(dataFct, Dune::VTK::FieldInfo(name_, Dune::VTK::FieldInfo::Type::scalar, 1));
write();
vtkWriter_->clear();
}
template <std::size_t C, class W>
void writeVertexData(Dune::TypeTree::PowerNodeTag, index_t<C>, W write)
{
using namespace Dune::Functions::BasisBuilder;
assert( C == dow );
auto fct = makeDiscreteFunction(basis_,treePath_,vector_);
auto p1basis = makeBasis(basis_->gridView(), power<C>(lagrange<1>(), flatLexicographic()));
auto data = makeDOFVector(p1basis, name_);
interpolate(p1basis, data, *fct);
using Range = Dune::FieldVector<double,C>;
auto dataFct = Dune::Functions::makeDiscreteGlobalBasisFunction<Range>(p1basis,data);
vtkWriter_->addVertexData(dataFct, Dune::VTK::FieldInfo(name_, Dune::VTK::FieldInfo::Type::vector, C));
write();
vtkWriter_->clear();
}
template <class NodeTag, std::size_t C, class W>
void writeVertexData(NodeTag, index_t<C>, W) {}
#if 0
/// default write method for time-depended data /// default write method for time-depended data
template <class SystemVectorType> template <class SystemVectorType>
void write(double time, SystemVectorType const& solutions) void write(double time, SystemVectorType const& solutions)
...@@ -82,9 +144,6 @@ namespace AMDiS ...@@ -82,9 +144,6 @@ namespace AMDiS
vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/); vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/);
} }
protected:
template <class DOFVector, class Vector> template <class DOFVector, class Vector>
void dofVector2vertexVector(DOFVector const& dofvector, Vector& data) void dofVector2vertexVector(DOFVector const& dofvector, Vector& data)
{ {
...@@ -98,7 +157,7 @@ namespace AMDiS ...@@ -98,7 +157,7 @@ namespace AMDiS
auto localView = feSpace.localView(); auto localView = feSpace.localView();
auto localIndexSet = feSpace.localIndexSet(); auto localIndexSet = feSpace.localIndexSet();
#if 0
// copy data to P1-vector // copy data to P1-vector
for (auto const& element : elements(meshView)) { for (auto const& element : elements(meshView)) {
localView.bind(element); localView.bind(element);
...@@ -124,19 +183,21 @@ namespace AMDiS ...@@ -124,19 +183,21 @@ namespace AMDiS
} }
} }
} }
#endif
} }
#endif
private: private:
MeshView meshView; std::shared_ptr<GlobalBasis> basis_;
shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter; TreePath treePath_;
shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter; std::shared_ptr<Vector> vector_;
std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_;
std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_;
std::vector<std::string> names; // std::vector<double> data_vector;
std::array<std::vector<double>, nComponents> data_vectors;
std::string filename; // represents VTK::OutputType: 0...ascii, 1...appendedraw
std::string dir; int mode_;
}; };
} // end namespace AMDiS } // end namespace AMDiS
...@@ -23,15 +23,14 @@ namespace AMDiS ...@@ -23,15 +23,14 @@ namespace AMDiS
using ElementMatrix = Impl::ElementMatrix; using ElementMatrix = Impl::ElementMatrix;
using ElementVector = Impl::ElementVector; using ElementVector = Impl::ElementVector;
using OrderTag = std::conditional_t<(type==GRD_PHI), tag::fot_grd_phi, tag::fot_grd_psi>;
public: public:
template <class Operator> template <class Operator>
FirstOrderAssembler(Operator& op, LocalContext const& element) FirstOrderAssembler(Operator& op, LocalContext const& element)
: Super(op, element, 1, type) : Super(op, element, 1, type)
{ {
if (type == GRD_PHI) op.init(OrderTag{}, element, Super::getQuadrature());
op.initFot1(element, Super::getQuadrature());
else
op.initFot2(element, Super::getQuadrature());
} }
...@@ -100,12 +99,13 @@ namespace AMDiS ...@@ -100,12 +99,13 @@ namespace AMDiS
for (std::size_t j = 0; j < colLocalFE.size(); ++j) { for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
const int local_i = rowView.tree().localIndex(i); const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j); const int local_j = colView.tree().localIndex(j);
op.evalFot1(elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j]); op.eval(tag::fot_grd_phi{}, elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j]);
} }
} }
} }
} }
template <class Operator, class RowView, class ColView> template <class Operator, class RowView, class ColView>
void calculateElementMatrix(Operator& op, void calculateElementMatrix(Operator& op,
RowView const& rowView, ColView const& colView, RowView const& rowView, ColView const& colView,
...@@ -148,7 +148,7 @@ namespace AMDiS ...@@ -148,7 +148,7 @@ namespace AMDiS
for (std::size_t j = 0; j < colLocalFE.size(); ++j) { for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
const int local_i = rowView.tree().localIndex(i); const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j); const int local_j = colView.tree().localIndex(j);
op.evalFot2(elementMatrix[local_i][local_j], iq, factor, rowGradients[i], colShapeValues[j]); op.eval(tag::fot_grd_psi{}, elementMatrix[local_i][local_j], iq, factor, rowGradients[i], colShapeValues[j]);
} }
} }
} }
...@@ -201,13 +201,14 @@ namespace AMDiS ...@@ -201,13 +201,14 @@ namespace AMDiS
for (std::size_t j = 0; j < colLocalFE.size(); ++j) { for (std::size_t j = 0; j < colLocalFE.size(); ++j) {
for (std::size_t k = 0; k < degree(colView.tree()); ++k) { for (std::size_t k = 0; k < degree(colView.tree()); ++k) {
const int local_j = colView.tree().child(k).localIndex(j); const int local_j = colView.tree().child(k).localIndex(j);
op.evalZot(elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j][k]); op.eval(tag::fot_grd_phi{}, elementMatrix[local_i][local_j], iq, factor, rowShapeValues[i], colGradients[j][k]);
} }
} }
} }
} }
} }
template <class Operator, class RowView, class ColView> template <class Operator, class RowView, class ColView>
void calculateElementMatrix(Operator& op, void calculateElementMatrix(Operator& op,
RowView const& rowView, ColView const& colView, RowView const& rowView, ColView const& colView,
...@@ -254,7 +255,7 @@ namespace AMDiS ...@@ -254,7 +255,7 @@ namespace AMDiS
for (std::size_t i = 0; i < rowLocalFE.size(); ++i) { for (std::size_t i = 0; i < rowLocalFE.size(); ++i) {
for (std::size_t k = 0; k < degree(rowView.tree()); ++k) { for (std::size_t k = 0; k < degree(rowView.tree()); ++k) {
const int local_i = rowView.tree().child(k).localIndex(i); const int local_i = rowView.tree().child(k).localIndex(i);
op.evalZot(elementMatrix[local_i][local_j], iq, factor, rowGradients[i][k], colShapeValues[j]); op.eval(tag::fot_grd_psi{}, elementMatrix[local_i][local_j], iq, factor, rowGradients[i][k], colShapeValues[j]);
} }
} }
} }
...@@ -269,7 +270,10 @@ namespace AMDiS ...@@ -269,7 +270,10 @@ namespace AMDiS
ElementMatrix& elementMatrix, ElementMatrix& elementMatrix,
double fac, double fac,
RowNodeTag, ColNodeTag, FOT) RowNodeTag, ColNodeTag, FOT)
{} {
warning("FirstOrderAssembler::calculateElementMatrix not implemented for RowNode x ColNode");
}
template <class Operator, class RowView> template <class Operator, class RowView>
void calculateElementVector(Operator& op, void calculateElementVector(Operator& op,
...@@ -307,18 +311,21 @@ namespace AMDiS ...@@ -307,18 +311,21 @@ namespace AMDiS
for (std::size_t i = 0; i < localFE.size(); ++i) { for (std::size_t i = 0; i < localFE.size(); ++i) {
const int local_i = rowView.tree().localIndex(i); const int local_i = rowView.tree().localIndex(i);
op.evalFot2(elementVector[local_i], iq, factor, rowGradients[i]); op.eval(tag::fot_grd_psi{}, elementVector[local_i], iq, factor, rowGradients[i]);
} }
} }
} }
template <class Operator, class RowView, class NodeTag, class FOT> template <class Operator, class RowView, class NodeTag, class FOT>
void calculateElementVector(Operator& op, void calculateElementVector(Operator& op,
RowView const& rowView, RowView const& rowView,
ElementVector& elementVector, ElementVector& elementVector,
double fac, double fac,
NodeTag, FOT) NodeTag, FOT)
{} {
warning("FirstOrderAssembler::calculateElementVector not implemented for RowNode");
}
}; };
} // end namespace AMDiS } // end namespace AMDiS
...@@ -65,15 +65,15 @@ namespace AMDiS ...@@ -65,15 +65,15 @@ namespace AMDiS
void Initfile::getInternalParameters() void Initfile::getInternalParameters()
{ {
int val = 0; int val = 0;
get("level of information", val, 0); get("level of information", val);
msgInfo = val; msgInfo = val;
val = 1; val = 1;
get("parameter information", val, 0); get("parameter information", val);
paramInfo = val; paramInfo = val;
val = 0; val = 0;
get("break on missing tag", val, 0); get("break on missing tag", val);
breakOnMissingTag = val; breakOnMissingTag = val;
if (msgInfo == 0) if (msgInfo == 0)
...@@ -88,9 +88,9 @@ namespace AMDiS ...@@ -88,9 +88,9 @@ namespace AMDiS
} }
// explicit template instatiation // explicit template instatiation
template void Initfile::get(std::string, int&, int); template void Initfile::get(std::string, int&);
template void Initfile::get(std::string, double&, int); template void Initfile::get(std::string, double&);
template void Initfile::get(std::string, std::string&, int); template void Initfile::get(std::string, std::string&);
} // end namespace AMDiS } // end namespace AMDiS
...@@ -190,62 +190,30 @@ namespace AMDiS ...@@ -190,62 +190,30 @@ namespace AMDiS
static void init(std::string in); static void init(std::string in);
/** Static get routine for getting parameter-values from init-file
* initialized in init()-method.
* Cast the value to the desired type using std::stringstream.
* @param tag: The tag to look for
* @param value: The result.
* @param debugInfo: msgInfo for current parameter. (0..no printing,
* 1..print missing parameter info, 2..print parameter value) [optional]
*/
template <class T> template <class T>
static void get(std::string tag, T& value, int debugInfo = -1) static boost::optional<T> get(std::string tag)
{ {
if (debugInfo == -1)
debugInfo = singlett().getMsgInfo();
else
{
int swap(debugInfo);
debugInfo = singlett().getMsgInfo();
singlett().msgInfo=swap;
}
using path = boost::property_tree::ptree::path_type; using path = boost::property_tree::ptree::path_type;
replaceAll(tag, "->", ">"); replaceAll(tag, "->", ">");
auto tagPath = path(tag, '>'); auto tagPath = path(tag, '>');
// TODO: use boost::optional instead return singlett().pt.get_optional<T>(tagPath);
// TODO: use convert method from above
std::string valueStr = "-";
valueStr = singlett().pt.get(tagPath, valueStr);
if (valueStr != "-")
detail::Convert<T>::eval(valueStr, value);
// if (debugInfo == 2)
// {
// std::cout << "Parameter '" << tag << "'"
// << " initialized with: " << value << std::endl;
// }
singlett().msgInfo = debugInfo;
}
template <class T, class S>
static T get(std::string tag, S const& default_value)
{
T value = default_value;
Self::get(tag, value);
return value;
} }
/** \brief Static get routine for getting parameter-values from init-file
* initialized in init()-method.
* Cast the value to the desired type.
*
* \param tag: The tag to look for
* \param value: The result.
*/
template <class T> template <class T>
static T get(std::string tag) static void get(std::string tag, T& value)
{ {
T value; nullify(value); auto v = get<std::string>(tag);
Self::get(tag, value); if (v)
return value; detail::Convert<T>::eval(v.value(), value);
} }
...@@ -336,9 +304,9 @@ namespace AMDiS ...@@ -336,9 +304,9 @@ namespace AMDiS
using Parameters = Initfile;