Commit ed9107af authored by Praetorius, Simon's avatar Praetorius, Simon

some renaming and started to implement boundary terms

parent 633b910f
Pipeline #904 passed with stage
in 13 minutes and 48 seconds
......@@ -68,7 +68,7 @@ namespace AMDiS
void AdaptInfo::printTimeErrorLowInfo() const
{
for (size_t i = 0; i < scalContents.size(); i++)
for (std::size_t i = 0; i < scalContents.size(); i++)
{
std::cout << " Time error estimate ["<<i<<"] = "
<< getTimeEstCombined(i) << "\n"
......
......@@ -87,7 +87,7 @@ namespace AMDiS
/// Destructor.
virtual ~AdaptInfo()
{
for (size_t i = 0; i < scalContents.size(); i++)
for (std::size_t i = 0; i < scalContents.size(); i++)
delete scalContents[i];
}
......@@ -97,7 +97,7 @@ namespace AMDiS
/// Returns whether space tolerance is reached.
virtual bool spaceToleranceReached() const
{
for (size_t i = 0; i < scalContents.size(); i++)
for (std::size_t i = 0; i < scalContents.size(); i++)
{
if (!(scalContents[i]->est_sum < scalContents[i]->spaceTolerance))
return false;
......@@ -118,7 +118,7 @@ namespace AMDiS
/// Returns whether time tolerance is reached.
virtual bool timeToleranceReached() const
{
for (size_t i = 0; i < scalContents.size(); i++)
for (std::size_t i = 0; i < scalContents.size(); i++)
if (!(getTimeEstCombined(i) < scalContents[i]->timeTolerance))
return false;
......@@ -137,7 +137,7 @@ namespace AMDiS
/// Returns whether time error is under its lower bound.
virtual bool timeErrorLow() const
{
for (size_t i = 0; i < scalContents.size(); i++)
for (std::size_t i = 0; i < scalContents.size(); i++)
if (!(getTimeEstCombined(i) < scalContents[i]->timeErrLow))
return false;
......@@ -145,7 +145,7 @@ namespace AMDiS
}
/// Returns the time estimation as a combination
/// of maximal and integral time error
double getTimeEstCombined(size_t i) const
double getTimeEstCombined(std::size_t i) const
{
return
scalContents[i]->est_t_max * scalContents[i]->fac_max +
......@@ -304,7 +304,7 @@ namespace AMDiS
double getEstSum(int index) const
{
AMDIS_FUNCNAME_DBG("AdaptInfo::getEstSum()");
test_exit_dbg(static_cast<size_t>(index) < scalContents.size(),
test_exit_dbg(static_cast<std::size_t>(index) < scalContents.size(),
"Wrong index for adaptInfo!\n");
return scalContents[index]->est_sum;
......@@ -320,7 +320,7 @@ namespace AMDiS
double getEstMax(int index) const
{
AMDIS_FUNCNAME_DBG("AdaptInfo::getEstSum()");
test_exit_dbg(static_cast<size_t>(index) < scalContents.size(),
test_exit_dbg(static_cast<std::size_t>(index) < scalContents.size(),
"Wrong index for adaptInfo!\n");
return scalContents[index]->est_max;
......
......@@ -56,7 +56,7 @@ namespace AMDiS
{
vtkWriter->clear();
// copy dofvector to vertex data
for_<0, nComponents>([this, &solutions](const auto _i)
forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
{
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
......@@ -71,7 +71,7 @@ namespace AMDiS
{
vtkWriter->clear();
// copy dofvector to vertex data
for_<0, nComponents>([this, &solutions](const auto _i)
forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
{
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
......@@ -106,16 +106,16 @@ namespace AMDiS
std::vector<Dune::FieldVector<double,1> > shapeValues;
size_t nVertices = element.subEntities(dim);
for (size_t i = 0; i < nVertices; ++i) {
std::size_t nVertices = element.subEntities(dim);
for (std::size_t i = 0; i < nVertices; ++i) {
auto const& v = element.template subEntity<dim>(i);
auto pos = refElement.position(i, dim);
localBasis.evaluateFunction(pos, shapeValues);
size_t idx = indexSet.index(v);
std::size_t idx = indexSet.index(v);
data[idx] = 0.0;
for (size_t j = 0; j < shapeValues.size(); ++j) {
for (std::size_t j = 0; j < shapeValues.size(); ++j) {
const auto global_idx = localIndexSet.index(j);
data[idx] += dofvector[global_idx] * shapeValues[j];
}
......
......@@ -91,7 +91,7 @@ namespace AMDiS
};
// convert string to vector
template <class V, size_t dim>
template <class V, std::size_t dim>
struct Convert<std::array<V, dim>>
{
using T = std::array<V, dim>;
......@@ -102,7 +102,7 @@ namespace AMDiS
boost::char_separator<char> sep(",; ");
Tokenizer tokens(valStr, sep);
size_t i = 0;
std::size_t i = 0;
for (auto token : tokens)
{
test_exit(i < dim, "Vector data exceeds array dimension!");
......@@ -166,7 +166,7 @@ namespace AMDiS
{
if (from.empty())
return;
size_t start_pos = 0;
std::size_t start_pos = 0;
while ((start_pos = str.find(from, start_pos)) != std::string::npos)
{
str.replace(start_pos, from.length(), to);
......
......@@ -5,7 +5,7 @@
#include <cmath>
#include <cfloat>
#include <boost/math/special_functions/pow.hpp>
#include <boost/math/special_functions/pow.hpp>
#include <dune/amdis/common/ScalarTypes.hpp>
......@@ -30,14 +30,14 @@ namespace AMDiS
return a*a;
}
template <size_t p, class T,
template <std::size_t p, class T,
class = std::enable_if_t<Concepts::Arithmetic<T>> >
constexpr auto pow(T v)
{
return boost::math::pow<p>(v);
}
/// Implementation of the minimum of two values \f$ min(a,b)\f$ of any type
/// Implementation of the minimum of two values \f$ min(a,b)\f$ of any type
/// supporting the `>` relation.
template <class T0, class T1>
constexpr auto min(T0 a, T1 b)
......@@ -46,7 +46,7 @@ namespace AMDiS
}
/// Implementation of the maximum of two values \f$ max(a,b)\f$ of any type
/// Implementation of the maximum of two values \f$ max(a,b)\f$ of any type
/// supporting the `>` relation.
template <class T0, class T1>
constexpr auto max(T0 a, T1 b)
......
......@@ -107,10 +107,10 @@ namespace AMDiS
/// Calls \ref zot(), \ref for() or \ref sot(), depending on template
/// parameter \p Order.
template <size_t Order, class... Args>
template <std::size_t Order, class... Args>
static shared_ptr<Self> create(Args&&... args)
{
return create(index_<Order>{}, std::forward<Args>(args)...);
return create(index_<Order>, std::forward<Args>(args)...);
}
......@@ -186,19 +186,19 @@ namespace AMDiS
template <class... Args>
static shared_ptr<Self> create(index_<0>, Args&&... args)
static shared_ptr<Self> create(index_t<0>, Args&&... args)
{
return zot(std::forward<Args>(args)...);
}
template <class... Args>
static shared_ptr<Self> create(index_<1>, Args&&... args)
static shared_ptr<Self> create(index_t<1>, Args&&... args)
{
return fot(std::forward<Args>(args)...);
}
template <class... Args>
static shared_ptr<Self> create(index_<2>, Args&&... args)
static shared_ptr<Self> create(index_t<2>, Args&&... args)
{
return sot(std::forward<Args>(args)...);
}
......
......@@ -194,7 +194,7 @@ namespace AMDiS
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
Dune::FieldVector<double,dim> const& quadPos = quad[iq].position();
......@@ -209,8 +209,8 @@ namespace AMDiS
else
colLocalBasis.evaluateFunction(quadPos, colShapeValues);
for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : zeroOrder)
......@@ -243,7 +243,7 @@ namespace AMDiS
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -253,7 +253,7 @@ namespace AMDiS
std::vector<Dune::FieldVector<double,1> > rowShapeValues;
rowLocalBasis.evaluateFunction(quadPos, rowShapeValues);
for (size_t i = 0; i < size(elementvector); ++i) {
for (std::size_t i = 0; i < size(elementvector); ++i) {
const int local_i = rowView.tree().localIndex(i);
for (auto* operatorTerm : zeroOrder)
elementvector[local_i] += operatorTerm->evalZot(iq, rowShapeValues[i]) * factor;
......@@ -284,7 +284,7 @@ namespace AMDiS
for (auto* operatorTerm : firstOrderGrdPhi)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -304,11 +304,11 @@ namespace AMDiS
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dow> > colGradients(colReferenceGradients.size());
for (size_t i = 0; i < colGradients.size(); ++i)
for (std::size_t i = 0; i < colGradients.size(); ++i)
jacobian.mv(colReferenceGradients[i][0], colGradients[i]);
for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : firstOrderGrdPhi)
......@@ -342,7 +342,7 @@ namespace AMDiS
for (auto* operatorTerm : firstOrderGrdPsi)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -359,14 +359,14 @@ namespace AMDiS
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
for (size_t i = 0; i < rowGradients.size(); ++i)
for (std::size_t i = 0; i < rowGradients.size(); ++i)
jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
std::vector<Dune::FieldVector<double,1> > colShapeValues;
colLocalBasis.evaluateFunction(quadPos, colShapeValues);
for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : firstOrderGrdPsi)
......@@ -399,7 +399,7 @@ namespace AMDiS
for (auto* operatorTerm : firstOrderGrdPsi)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -416,10 +416,10 @@ namespace AMDiS
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dow> > rowGradients(rowReferenceGradients.size());
for (size_t i = 0; i < rowGradients.size(); ++i)
for (std::size_t i = 0; i < rowGradients.size(); ++i)
jacobian.mv(rowReferenceGradients[i][0], rowGradients[i]);
for (size_t i = 0; i < size(elementvector); ++i) {
for (std::size_t i = 0; i < size(elementvector); ++i) {
const int local_i = rowView.tree().localIndex(i);
for (auto* operatorTerm : firstOrderGrdPsi)
elementvector[local_i] += operatorTerm->evalFot2(iq, rowGradients[i]) * factor;
......@@ -453,7 +453,7 @@ namespace AMDiS
// TODO: currently only the implementation for equal fespaces
assert( psiDegree == phiDegree );
for (size_t iq = 0; iq < quad.size(); ++iq) {
for (std::size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -470,11 +470,11 @@ namespace AMDiS
// Compute the shape function gradients on the real element
std::vector<Dune::FieldVector<double,dow> > gradients(referenceGradients.size());
for (size_t i = 0; i < gradients.size(); ++i)
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
for (size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (size_t j = 0; j < num_cols(elementMatrix); ++j) {
for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
const int local_i = rowView.tree().localIndex(i);
const int local_j = colView.tree().localIndex(j);
for (auto* operatorTerm : secondOrder)
......@@ -513,7 +513,7 @@ namespace AMDiS
template <class MeshView>
template <class Term>
Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
size_t i,
std::size_t i,
FirstOrderType firstOrderType)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent>;
......@@ -538,7 +538,7 @@ namespace AMDiS
template <class MeshView>
template <class Term>
Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term,
size_t i, size_t j)
std::size_t i, std::size_t j)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent>;
secondOrder.push_back(new OpTerm(term, {i,j}));
......
......@@ -32,19 +32,19 @@ namespace AMDiS
public:
virtual void init(Element const& element, PointList const& points) = 0;
virtual double evalZot(size_t iq,
virtual double evalZot(std::size_t iq,
Dune::FieldVector<double,1> const& test,
Dune::FieldVector<double,1> const trial = 1.0) const = 0;
virtual double evalFot1(size_t iq,
virtual double evalFot1(std::size_t iq,
Dune::FieldVector<double,1> const& test,
Dune::FieldVector<double,dow> const& grad_trial) const = 0;
virtual double evalFot2(size_t iq,
virtual double evalFot2(std::size_t iq,
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,1> const trial = 1.0) const = 0;
virtual double evalSot(size_t iq,
virtual double evalSot(std::size_t iq,
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,dow> const& grad_trial) const = 0;
......@@ -78,32 +78,32 @@ namespace AMDiS
// cache term evaluation
values.resize(points.size());
for (size_t iq = 0; iq < points.size(); ++iq)
for (std::size_t iq = 0; iq < points.size(); ++iq)
values[iq] = term[iq];
}
virtual double evalZot(size_t iq,
virtual double evalZot(std::size_t iq,
Dune::FieldVector<double,1> const& test,
Dune::FieldVector<double,1> const trial = 1.0) const override
{
return this->evalZotImpl(_cat{}, traits, values[iq], test, trial);
}
virtual double evalFot1(size_t iq,
virtual double evalFot1(std::size_t iq,
Dune::FieldVector<double,1> const& test,
Dune::FieldVector<double,dow> const& grad_trial) const override
{
return this->evalFotImpl(_cat{}, traits, values[iq], grad_trial, test);
}
virtual double evalFot2(size_t iq,
virtual double evalFot2(std::size_t iq,
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,1> const trial = 1.0) const override
{
return this->evalFotImpl(_cat{}, traits, values[iq], grad_test, trial);
}
virtual double evalSot(size_t iq,
virtual double evalSot(std::size_t iq,
Dune::FieldVector<double,dow> const& grad_test,
Dune::FieldVector<double,dow> const& grad_trial) const override
{
......@@ -119,7 +119,7 @@ namespace AMDiS
Term term;
Traits traits;
using value_type = std::decay_t< decltype( std::declval<Term>()[std::declval<size_t>()] ) >;
using value_type = std::decay_t< decltype( std::declval<Term>()[std::size_t(0)] ) >;
using _cat = ValueCategory_t<value_type>;
std::vector<value_type> values; // NOTE: maybe caching is not necessary here, since cached already in the term
......
......@@ -59,8 +59,8 @@ namespace AMDiS
decltype(auto) getOldSolution() const { return *oldSolution; }
/// Returns the I'th component of \ref oldSolution.
template <size_t I = 0>
decltype(auto) getOldSolution(const index_<I> _i = {})
template <std::size_t I = 0>
decltype(auto) getOldSolution(const index_t<I> _i = {})
{
return (*oldSolution)[_i];
}
......
......@@ -26,6 +26,7 @@
#include <dune/amdis/ProblemStatTraits.hpp>
#include <dune/amdis/StandardProblemIteration.hpp>
#include <dune/amdis/common/OptionalDelete.hpp>
#include <dune/amdis/common/Timer.hpp>
#include <dune/amdis/common/TupleUtility.hpp>
#include <dune/amdis/common/TypeDefs.hpp>
......@@ -63,7 +64,7 @@ namespace AMDiS
static constexpr int nComponents = Traits::nComponents;
template <size_t I>
template <std::size_t I>
using FeSpace = std::tuple_element_t<I, FeSpaces>;
using WorldVector = typename Codim0::Geometry::GlobalCoordinate;
......@@ -87,10 +88,16 @@ namespace AMDiS
: name(name)
{
Parameters::get(name + "->names", componentNames);
for (size_t i = componentNames.size(); i < nComponents; ++i)
for (std::size_t i = componentNames.size(); i < nComponents; ++i)
componentNames.push_back("solution[" + std::to_string(i) + "]");
}
ProblemStatSeq(std::string name, Mesh& mesh)
: ProblemStatSeq(name)
{
this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
}
/**
* \brief Initialisation of the problem.
......@@ -181,8 +188,8 @@ namespace AMDiS
auto getSolution() const { return solution; }
/// Return the i'th \ref solution component
template <size_t I = 0>
decltype(auto) getSolution(const index_<I> _i = {})
template <std::size_t I = 0>
decltype(auto) getSolution(const index_t<I> _i = {})
{
return (*solution)[_i];
}
......@@ -202,12 +209,11 @@ namespace AMDiS
}
/// Return the \ref mesh
auto getMesh() { return mesh; }
auto const& getMesh() { return mesh; }
void setMesh(std::shared_ptr<Mesh> const& mesh_)
void setMesh(Mesh& mesh_)
{
mesh = mesh_;
meshView = std::make_shared<MeshView>(mesh->leafGridView());
mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
createFeSpaces();
createMatricesAndVectors();
......@@ -215,16 +221,25 @@ namespace AMDiS
}
/// Return the \ref meshView
auto getMeshView() { return meshView; }
/// Return the gridView of the leaf-level
auto leafGridView() { return mesh->leafGridView(); }
/// Return the gridView of levle `level`
auto levelGridView(int level) { return mesh->levelGridView(level); }
/// Return the \ref feSpaces
auto getFeSpaces() { return feSpaces; }
/// Return the I'th \ref feSpaces component
template <size_t I = 0>
FeSpace<I> const& getFeSpace(const index_<I> = {}) const
template <std::size_t I = 0>
FeSpace<I> const& getFeSpace(const index_t<I> = {}) const
{
return std::get<I>(*feSpaces);
}
template <std::size_t I = 0>
FeSpace<I>& getFeSpace(const index_t<I> = {})
{
return std::get<I>(*feSpaces);
}
......@@ -251,7 +266,6 @@ namespace AMDiS
"No mesh name specified for '", name, "->mesh'!");
mesh = MeshCreator<Mesh>::create(meshName);
meshView = std::make_shared<MeshView>(mesh->leafGridView());
msg("Create mesh:");
msg("#elements = " , mesh->size(0));
......@@ -262,7 +276,7 @@ namespace AMDiS
void createFeSpaces()
{
feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(*meshView));
feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(leafGridView()));
}
void createMatricesAndVectors()
......@@ -287,14 +301,14 @@ namespace AMDiS
void createFileWriter()
{
filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", *meshView, componentNames);
filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", leafGridView(), componentNames);
}
protected: // sub-methods to assemble DOFMatrix
template <class LhsData, class RhsData, class Elements>
void assemble(LhsData lhs, RhsData rhs, Elements const& elements);
template <class LhsData, class RhsData, class GV>
void assemble(LhsData lhs, RhsData rhs, GV const& gridView);
template <class RowView, class ColView>
bool getElementMatrix(RowView const& rowView,
......@@ -334,10 +348,7 @@ namespace AMDiS
std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
/// Name of the mesh
std::string meshName;
/// A gridView object
std::shared_ptr<MeshView> meshView;
std::string meshName = "none";
/// Pointer to the meshes for the different problem components
std::vector<Mesh*> componentMeshes;
......@@ -396,7 +407,7 @@ namespace AMDiS
DOFMatrixType& matrix;
std::list<std::shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors;
bool assemble;
bool assemble = true;
std::pair<int,int> row_col = {R, C};
};
......@@ -409,7 +420,7 @@ namespace AMDiS
DOFVectorType& vector;
std::list<std::shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors;
bool assemble;
bool assemble = true;
int row = I;
};
......@@ -438,6 +449,12 @@ namespace AMDiS
, StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
{}
template <class Grid>
ProblemStat(std::string nameStr, Grid& grid)
: ProblemStatType(nameStr, grid)
, StandardProblemIteration(dynamic_cast<ProblemStatBase&>(*this))
{}
/**
* \brief Determines the execution order of the single adaption steps.
*
......@@ -449,7 +466,7 @@ namespace AMDiS
**/
virtual Flag oneIteration(AdaptInfo& adaptInfo, Flag toDo = FULL_ITERATION) override
{
for (size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
for (std::size_t i = 0; i < ProblemStatType::getNumComponents(); ++i)
if (adaptInfo.spaceToleranceReached(i))
adaptInfo.allowRefinement(false, i);
else
......
......@@ -258,28 +258,27 @@ namespace AMDiS
{
Timer t;
For<0, nComponents>::loop([this](auto const _r) {
this->getFeSpace(_r).update(this->getMeshView());
forEach(range_<0, nComponents>, [this](auto const _r) {
this->getFeSpace(_r).update(this->leafGridView());
});
size_t nnz = 0;
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const _r)
std::size_t nnz = 0;
forEach(range_<0, nComponents>, [&,this](auto const _r)
{
msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", _r, "]");
static const int R = decltype(_r)::value;
msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", R, "]");
bool asmVector = asmVector_ && (!vectorAssembled[int(_r)] || vectorChanging[int(_r)]);
bool asmVector = asmVector_ && (!vectorAssembled[R] || vectorChanging[R]);
if (asmVector) {
rhs->compress(_r);
rhs->getVector(_r) = 0.0;
}
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector](auto const _c)
forEach(range_<0, nComponents>, [&,this](auto const _c)
{
static const int C = decltype(_c)::value;
static const int R = decltype(_r)::value;
index_<R> _r = {};
using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<R, C>;
using VectorData = typename ProblemStatSeq<Traits>::template VectorData<R>;
......@@ -311,7 +310,8 @@ namespace AMDiS
// assemble the DOFMatrix block and the corresponding rhs vector, of r==c
dofmatrix.init(assembleMatrix);
this->assemble(mat, vec, elements(*meshView));
this->assemble(mat, vec, this->leafGridView());
// TODO: assemble boundary conditions
dofmatrix.finish();
if (assembleMatrix)
......@@ -350,9 +350,9 @@ namespace AMDiS
template <class Traits>
template <class LhsData, class RhsData, class Elements>
template <class LhsData, class RhsData, class GV>
void ProblemStatSeq<Traits>::assemble(LhsData lhs, RhsData rhs,
Elements const& elements)
GV const& gridView)
{
auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
auto const& colFeSpace = lhs.matrix.getColFeSpace();
......@@ -373,7 +373,7 @@ namespace AMDiS
auto colLocalView = colFeSpace.localView();
auto colIndexSet = colFeSpace.localIndexSet();
for (auto const& element : elements) {
for (auto const& element : elements(gridView)) {
rowLocalView.bind(element);
rowIndexSet.bind(rowLocalView);
......@@ -426,6 +426,24 @@ namespace AMDiS
add = add || add_op;
}
#if 0
auto const& element = rowLocalView.element();