Skip to content
Snippets Groups Projects
Commit 8c1c1d09 authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Make compile and clean-up

parent 68978f68
No related branches found
No related tags found
No related merge requests found
...@@ -2,7 +2,6 @@ ...@@ -2,7 +2,6 @@
#define DUNE_GEOMETRIC_ERROR_ESTIMATOR_HH #define DUNE_GEOMETRIC_ERROR_ESTIMATOR_HH
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
/** \brief Error estimator which marks an element for refinement if it fulfills a geometric condition /** \brief Error estimator which marks an element for refinement if it fulfills a geometric condition
*/ */
...@@ -13,64 +12,49 @@ class GeometricEstimator { ...@@ -13,64 +12,49 @@ class GeometricEstimator {
enum {dimworld = GridType::dimensionworld}; enum {dimworld = GridType::dimensionworld};
public: public:
static void estimate(GridType& grid,
bool (*refinementCondition)(const Dune::FieldVector<double, dim>& pos)) {
void estimate(GridType& grid, for (const auto& e : elements(grid.leafGridView())) {
bool (*refinementCondition)(const Dune::FieldVector<double, GridType::dimension>& pos) ) {
typedef typename GridType::template Codim<0>::LeafIterator LeafIterator;
LeafIterator lIt = grid.template leafbegin<0>();
LeafIterator lEndIt = grid.template leafend<0>();
for (; lIt!=lEndIt; ++lIt) {
// Compute element's center of mass // Compute element's center of mass
Dune::FieldVector<double, GridType::dimension> center(0); const auto& geometry = e.geometry();
int numCorners = lIt->template count<GridType::dimension>(); auto center = geometry.corner(0);
for (int i=0; i<numCorners; i++) for (int i=1; i<geometry.corners(); i++)
center += lIt->geometry().corner(i); center += geometry.corner(i);
center /= (double)numCorners; center /= (double)geometry.corners();
// Mark if the condition is fulfilled // Mark if the condition is fulfilled
if ( (*refinementCondition)(center) ) if ((*refinementCondition)(center) )
grid.mark(1, *lIt); grid.mark(1, e);
} }
} }
void estimate(GridType& grid, template <class BlockVectorType>
static void estimate(GridType& grid,
bool (*refinementCondition)(const Dune::FieldVector<double, dimworld>& pos), bool (*refinementCondition)(const Dune::FieldVector<double, dimworld>& pos),
const Dune::BlockVector<Dune::FieldVector<double, dimworld> >& deformation) { const BlockVectorType& deformation) {
const typename GridType::Traits::LeafIndexSet& indexSet = grid.leafIndexSet(); const auto& indexSet = grid.leafIndexSet();
typedef typename GridType::template Codim<0>::LeafIterator LeafIterator;
LeafIterator lIt = grid.template leafbegin<0>(); for (const auto& e : elements(grid.leafGridView())) {
LeafIterator lEndIt = grid.template leafend<0>();
for (; lIt!=lEndIt; ++lIt) {
// Compute element's center of mass // Compute element's center of mass
Dune::FieldVector<double, dimworld> center(0); Dune::FieldVector<double, dimworld> center(0);
int numCorners = lIt->template count<GridType::dimension>(); const auto& geometry = e.geometry();
for (int i=0; i<numCorners; i++) { for (int i=0; i<geometry.corners(); i++)
center += lIt->geometry()[i] + deformation[indexSet.template subIndex<dim>(*lIt,i)]; center += geometry.corner(i) + deformation[indexSet.template subIndex<dim>(e,i)];
}
center /= numCorners; center /= (double) geometry.corners();
// Mark if the condition is fulfilled // Mark if the condition is fulfilled
if ( (*refinementCondition)(center) > 0) if ( (*refinementCondition)(center) > 0)
grid.mark(1, lIt); grid.mark(1, e);
} }
} }
}; };
#endif #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment