Commit f0eafd9a authored by Andreas Nüßing's avatar Andreas Nüßing

Merge branch 'feature/meeg-add-statistics' into 'master'

MEEGDriver add statistics

See merge request duneuro/duneuro!49
parents ef4d2970 844f5cc3
......@@ -257,7 +257,7 @@ namespace duneuro
return deformedPositions;
}
template <class HostGrid, class L = char>
template <class HostGrid, class L = std::size_t>
std::unique_ptr<Dune::GeometryGrid<HostGrid,
DeformationFunction<typename HostGrid::LeafGridView>>>
create_geometry_adapted_grid(std::unique_ptr<HostGrid> hostGrid,
......@@ -297,7 +297,7 @@ namespace duneuro
std::unique_ptr<GeometryGrid> geometryGrid;
std::unique_ptr<GridType> grid;
std::unique_ptr<std::vector<char>> labels;
std::unique_ptr<std::vector<std::size_t>> labels;
};
template <int dim>
......@@ -325,7 +325,7 @@ namespace duneuro
}
grid->createEnd();
auto sgGv = grid->leafGridView();
auto subLabels = Dune::Std::make_unique<std::vector<char>>(grid->size(0));
auto subLabels = Dune::Std::make_unique<std::vector<std::size_t>>(grid->size(0));
for (const auto& subElement : Dune::elements(sgGv)) {
(*subLabels)[sgGv.indexSet().index(subElement)] =
data.labels[ggGv.indexSet().index(grid->template getHostEntity<0>(subElement))] - 1;
......@@ -337,12 +337,11 @@ namespace duneuro
std::shared_ptr<VolumeConductor<typename GeometryAdaptedGrid<dim>::GridType>>
make_geometry_adapted_volume_conductor(
std::unique_ptr<typename GeometryAdaptedGrid<dim>::GridType> grid,
std::unique_ptr<std::vector<char>> labels, const std::vector<double>& conductivities,
std::unique_ptr<std::vector<std::size_t>> labels, const std::vector<double>& conductivities,
const Dune::ParameterTree& config)
{
using VC = VolumeConductor<typename GeometryAdaptedGrid<dim>::GridType>;
using Tensor = typename VC::TensorType;
using MappingType = typename VC::MappingType;
std::vector<Tensor> tensors;
for (auto value : conductivities) {
Tensor t;
......@@ -354,10 +353,7 @@ namespace duneuro
tensors.push_back(t);
}
typename VC::GridView gv = grid->leafGridView();
return std::make_shared<VC>(
std::move(grid),
Dune::Std::make_unique<MappingType>(
IndirectEntityMapping<typename VC::GridView, Tensor, char>(gv, tensors, *labels)));
return std::make_shared<VC>(std::move(grid), *labels, tensors);
}
#if HAVE_NIFTI
......@@ -365,7 +361,7 @@ namespace duneuro
struct GeometryAdaptedGridReader {
static GeometryAdaptedGrid<dim> read(const Dune::ParameterTree& config)
{
std::vector<char> labels;
std::vector<std::size_t> labels;
std::array<unsigned int, dim> cells;
NiftiImageReader::read<dim>(config.get<std::string>("filename"), std::back_inserter(labels),
cells);
......@@ -389,7 +385,7 @@ namespace duneuro
}
grid->createEnd();
auto sgGv = grid->leafGridView();
auto subLabels = Dune::Std::make_unique<std::vector<char>>(grid->size(0));
auto subLabels = Dune::Std::make_unique<std::vector<std::size_t>>(grid->size(0));
for (const auto& subElement : Dune::elements(sgGv)) {
(*subLabels)[sgGv.indexSet().index(subElement)] =
labels[ggGv.indexSet().index(grid->template getHostEntity<0>(subElement))] - 1;
......
......@@ -13,66 +13,6 @@
namespace duneuro
{
template <class GV, class T, class I = std::size_t>
class IndirectEntityMapping
{
public:
typedef T TensorType;
typedef typename GV::template Codim<0>::Entity EntityType;
IndirectEntityMapping(const GV& gridView, const std::vector<T>& tensors,
const std::vector<I>& indexToTensor)
: gridView_(gridView), mapper_(gridView), tensors_(tensors), indexToTensor_(indexToTensor)
{
assert(tensors.size() > 0);
assert(indexToTensor_.size() == static_cast<std::size_t>(gridView_.size(0)));
}
const T& operator()(const EntityType& e) const
{
auto index = mapper_.index(e);
assert(index < indexToTensor_.size());
auto tensorIndex = indexToTensor_[index];
assert(tensorIndex < tensors_.size());
return tensors_[tensorIndex];
}
private:
// copy the grid view so that the mapper is always valid
GV gridView_;
Dune::SingleCodimSingleGeomTypeMapper<GV, 0> mapper_;
std::vector<T> tensors_;
// maps an entity index to the tensor index
std::vector<I> indexToTensor_;
};
template <class GV, class T>
class DirectEntityMapping
{
public:
typedef T TensorsType;
typedef typename GV::template Codim<0>::Entity EntityType;
DirectEntityMapping(const GV& gridView, const std::vector<T>& tensors)
: gridView_(gridView), mapper_(gridView_), tensors_(tensors)
{
assert(gridView_.size(0) == tensors_.size());
}
const T& operator()(const EntityType& e) const
{
auto index = mapper_.index(e);
assert(index < tensors_.size());
return tensors_[index];
}
private:
// copy the grid view so that the mapper is always valid
GV gridView_;
Dune::SingleCodimSingleGeomTypeMapper<GV, 0> mapper_;
std::vector<T> tensors_;
};
template <class G>
class VolumeConductor
{
......@@ -82,12 +22,31 @@ namespace duneuro
typedef typename G::ctype ctype;
typedef typename G::template Codim<0>::Entity EntityType;
typedef Dune::FieldMatrix<ctype, dim, dim> TensorType;
typedef std::function<const TensorType&(const EntityType&)> MappingType;
typedef typename G::LeafGridView GridView;
VolumeConductor(std::unique_ptr<G> grid, std::unique_ptr<MappingType> mapper)
: grid_(std::move(grid)), mapper_(std::move(mapper)), gridView_(grid_->leafGridView())
VolumeConductor(std::unique_ptr<G> grid, std::vector<std::size_t> labels,
std::vector<TensorType> tensors)
: grid_(std::move(grid))
, labels_(labels)
, tensors_(tensors)
, gridView_(grid_->leafGridView())
, elementMapper_(gridView_)
{
// check if we are given one label for each element
if (labels.size() != elementMapper_.size()) {
DUNE_THROW(Dune::Exception, "number of labels ("
<< labels.size()
<< ") has to match the number of elements ("
<< elementMapper_.size() << ")");
}
// check if all labels are valid
for (const auto& l : labels) {
if (l >= tensors_.size()) {
DUNE_THROW(Dune::Exception,
"label " << l << " is not a valid index into the tensors vector of size "
<< tensors_.size());
}
}
}
const G& grid() const
......@@ -107,23 +66,25 @@ namespace duneuro
const TensorType& tensor(const EntityType& entity) const
{
return (*mapper_)(entity);
return tensors_[label(entity)];
}
G* releaseGrid()
std::size_t label(const EntityType& entity) const
{
return grid_.release();
return labels_[elementMapper_.index(entity)];
}
MappingType* releaseMapping()
G* releaseGrid()
{
return mapper_.release();
return grid_.release();
}
private:
std::unique_ptr<G> grid_;
std::unique_ptr<MappingType> mapper_;
const std::vector<std::size_t> labels_;
const std::vector<TensorType> tensors_;
GridView gridView_;
Dune::SingleCodimSingleGeomTypeMapper<GridView, 0> elementMapper_;
};
}
......
#ifndef DUNEURO_VOLUME_CONDUCTOR_STATISTICS_HH
#define DUNEURO_VOLUME_CONDUCTOR_STATISTICS_HH
#include <map>
#include <dune/geometry/quadraturerules.hh>
#include <duneuro/common/volume_conductor.hh>
namespace duneuro
{
namespace VolumeConductorStatisticsDetail
{
template <class G>
typename G::ctype computeVolume(const G& geometry, std::size_t order)
{
using Real = typename G::ctype;
const int dim = G::mydimension;
const auto& rule = Dune::QuadratureRules<Real, dim>::rule(geometry.type(), order);
Real result = 0.0;
for (const auto& qp : rule) {
result += qp.weight() * geometry.integrationElement(qp.position());
}
return result;
}
}
struct VolumeConductorStatistics {
std::map<std::size_t, double> domainToVolume;
std::map<std::pair<std::size_t, std::size_t>, double> interfaceToVolume;
};
template <class G>
VolumeConductorStatistics
computeVolumeConductorStatistics(const VolumeConductor<G>& volumeConductor)
{
VolumeConductorStatistics result;
Dune::SingleCodimSingleGeomTypeMapper<typename VolumeConductor<G>::GridView, 0> elementMapper(
volumeConductor.gridView());
for (const auto& element : Dune::elements(volumeConductor.gridView())) {
std::size_t insideLabel = volumeConductor.label(element);
result.domainToVolume[insideLabel] +=
VolumeConductorStatisticsDetail::computeVolume(element.geometry(), 2);
auto insideIndex = elementMapper.index(element);
for (const auto& intersection : Dune::intersections(volumeConductor.gridView(), element)) {
std::size_t outsideLabel;
if (intersection.neighbor()) {
const auto& outside = intersection.outside();
if (insideIndex > elementMapper.index(outside)) {
continue;
}
outsideLabel = volumeConductor.label(outside);
} else {
outsideLabel = std::numeric_limits<std::size_t>::max();
}
std::size_t minLabel = std::min(insideLabel, outsideLabel);
std::size_t maxLabel = std::max(insideLabel, outsideLabel);
result.interfaceToVolume[std::make_pair(minLabel, maxLabel)] +=
VolumeConductorStatisticsDetail::computeVolume(intersection.geometry(), 2);
}
}
return result;
}
}
#endif // DUNEURO_VOLUME_CONDUCTOR_STATISTICS_HH
......@@ -24,7 +24,6 @@ namespace duneuro
{
public:
typedef typename VolumeConductor<G>::TensorType TensorType;
typedef typename VolumeConductor<G>::MappingType MappingType;
typedef typename G::LeafGridView GV;
typedef typename G::ctype ctype;
enum { dim = G::dimension };
......@@ -93,10 +92,8 @@ namespace duneuro
timer.stop();
dataTree.set("time_reordering_labels", timer.lastElapsed());
dataTree.set("time", timer.elapsed());
return std::make_shared<VolumeConductor<G>>(
std::move(grid),
std::unique_ptr<MappingType>(new MappingType(
IndirectEntityMapping<GV, TensorType>(gv, data.tensors, reordered_labels))));
return std::make_shared<VolumeConductor<G>>(std::move(grid), reordered_labels,
data.tensors);
} else if (data.labels.size() > 0) {
std::vector<std::size_t> reordered_labels(gv.size(0));
if (std::size_t(mapper.size()) != reordered_labels.size()) {
......@@ -134,10 +131,7 @@ namespace duneuro
tensors.push_back(t);
}
dataTree.set("time", timer.elapsed());
return std::make_shared<VolumeConductor<G>>(
std::move(grid),
std::unique_ptr<MappingType>(new MappingType(
IndirectEntityMapping<GV, TensorType>(gv, tensors, reordered_labels))));
return std::make_shared<VolumeConductor<G>>(std::move(grid), reordered_labels, tensors);
} else {
DUNE_THROW(Dune::Exception, "you have to provide labels or tensors");
}
......@@ -188,10 +182,7 @@ namespace duneuro
timer.stop();
dataTree.set("time_reordering_indices", timer.lastElapsed());
dataTree.set("time", timer.elapsed());
return std::make_shared<VolumeConductor<G>>(
std::move(grid),
std::unique_ptr<MappingType>(new MappingType(
IndirectEntityMapping<GV, TensorType>(gv, tensors, indexToTensor))));
return std::make_shared<VolumeConductor<G>>(std::move(grid), indexToTensor, tensors);
} else if (extension == "dgf") {
timer.start();
typedef Dune::SingleCodimSingleGeomTypeMapper<GV, 0> Mapper;
......@@ -224,10 +215,8 @@ namespace duneuro
timer.stop();
dataTree.set("time_reordering_tensors", timer.lastElapsed());
dataTree.set("time", timer.elapsed());
return std::make_shared<VolumeConductor<G>>(
std::unique_ptr<G>(gptr.release()),
std::unique_ptr<MappingType>(new MappingType(
IndirectEntityMapping<GV, TensorType>(gv, tensors, indexToTensor))));
return std::make_shared<VolumeConductor<G>>(std::unique_ptr<G>(gptr.release()),
indexToTensor, tensors);
} else if (extension == "geo") {
// assuming cauchy grid format
Dune::GridFactory<G> factory;
......@@ -238,13 +227,18 @@ namespace duneuro
Mapper mapper(gv);
std::vector<TensorType> tensors;
CauchyTensorReader<G>::read(tensorFilename, std::back_inserter(tensors));
if (tensors.size() != mapper.size()) {
DUNE_THROW(Dune::Exception,
"grid file contains a different number of elements than the tensors file ("
<< mapper.size() << " != " << tensors.size() << ")");
}
std::vector<TensorType> reorderedTensors(tensors.size());
std::vector<std::size_t> labels(tensors.size());
for (const auto& element : Dune::elements(gv)) {
labels.push_back(labels.size());
reorderedTensors[mapper.index(element)] = tensors[factory.insertionIndex(element)];
}
return std::make_shared<VolumeConductor<G>>(
std::move(grid), std::unique_ptr<MappingType>(new MappingType(
DirectEntityMapping<GV, TensorType>(gv, reorderedTensors))));
return std::make_shared<VolumeConductor<G>>(std::move(grid), labels, reorderedTensors);
} else {
DUNE_THROW(Dune::IOError, "cannot infer file type from extension \"" << extension << "\"");
}
......
......@@ -21,6 +21,7 @@
#include <duneuro/common/matrix_utilities.hh>
#include <duneuro/common/stl.hh>
#include <duneuro/common/volume_conductor.hh>
#include <duneuro/common/volume_conductor_statistics.hh>
#include <duneuro/common/volume_conductor_storage.hh>
#include <duneuro/eeg/cg_source_model_factory.hh>
#include <duneuro/eeg/conforming_eeg_forward_solver.hh>
......@@ -405,6 +406,21 @@ namespace duneuro
return projectedGlobalElectrodes_;
}
virtual void statistics(DataTree dataTree) const override
{
auto volumeConductorStatistics =
computeVolumeConductorStatistics(*(volumeConductorStorage_.get()));
auto sub = dataTree.sub("volume_conductor");
for (const auto& dtv : volumeConductorStatistics.domainToVolume) {
sub.set("volume_label_" + std::to_string(dtv.first), dtv.second);
}
for (const auto& itv : volumeConductorStatistics.interfaceToVolume) {
sub.set("surface_labels_" + std::to_string(itv.first.first) + "_"
+ std::to_string(itv.first.second),
itv.second);
}
}
private:
Dune::ParameterTree config_;
typename Traits::VCStorage volumeConductorStorage_;
......
......@@ -141,6 +141,15 @@ namespace duneuro
virtual std::vector<CoordinateType> getProjectedElectrodes() const = 0;
/**
* \brief obtain different statistics of the driver
*
* The results will be stored in the dataTree object. The entries depend on the implementation
* and might contain information such as the volume of the different compartments or the number
* of entities of different codimensions
*/
virtual void statistics(DataTree dataTree) const = 0;
virtual ~MEEGDriverInterface()
{
}
......
......@@ -27,6 +27,7 @@
#include <duneuro/io/vtk_functors.hh>
#include <duneuro/meeg/meeg_driver_interface.hh>
#include <duneuro/meeg/unfitted_meeg_driver_data.hh>
#include <duneuro/udg/subtriangulation_statistics.hh>
namespace duneuro
{
......@@ -347,6 +348,20 @@ namespace duneuro
return projectedGlobalElectrodes_;
}
virtual void statistics(DataTree dataTree) const override
{
auto subtriangulationStatistics = computeSubtriangulationStatistics(*subTriangulation_);
auto sub = dataTree.sub("subtriangulation");
for (const auto& dtv : subtriangulationStatistics.domainToVolume) {
sub.set("volume_label_" + std::to_string(dtv.first), dtv.second);
}
for (const auto& itv : subtriangulationStatistics.interfaceToVolume) {
sub.set("surface_labels_" + std::to_string(itv.first.first) + "_"
+ std::to_string(itv.first.second),
itv.second);
}
}
private:
void checkElectrodes() const
{
......
#ifndef DUNEURO_SUBTRIANGULATION_STATISTICS_HH
#define DUNEURO_SUBTRIANGULATION_STATISTICS_HH
#include <map>
#include <dune/geometry/quadraturerules.hh>
#include <dune/udg/pdelab/subtriangulation.hh>
#include <duneuro/common/volume_conductor.hh>
#include <duneuro/meeg/unfitted_meeg_driver.hh>
#include <duneuro/meeg/unfitted_meeg_driver_data.hh>
namespace duneuro
{
namespace SubtriangulationStatisticsDetail
{
template <class G>
typename G::ctype computeVolume(const G& geometry, std::size_t order)
{
using Real = typename G::ctype;
const int dim = G::mydimension;
const auto& rule = Dune::QuadratureRules<Real, dim>::rule(geometry.type(), order);
Real result = 0.0;
for (const auto& qp : rule) {
result += qp.weight() * geometry.integrationElement(qp.position());
}
return result;
}
}
struct SubtriangulationStatistics {
std::map<std::size_t, double> domainToVolume;
std::map<std::pair<std::size_t, std::size_t>, double> interfaceToVolume;
};
template <int dim>
struct SubTriangulationTraits;
template <class ST>
SubtriangulationStatistics computeSubtriangulationStatistics(const ST& subTriangulation)
{
using STTraits = SubTriangulationTraits<ST::dim>;
SubtriangulationStatistics result;
Dune::PDELab::UnfittedSubTriangulation<typename STTraits::GridView> ust(
subTriangulation.gridView(), subTriangulation);
for (const auto& element : Dune::elements(subTriangulation.gridView())) {
ust.create(element);
for (const auto& part : ust) {
result.domainToVolume[part.domainIndex()] +=
SubtriangulationStatisticsDetail::computeVolume(part.geometry(), 2);
}
for (auto it = ust.ibegin(); it != ust.iend(); ++it) {
if (it->insideDomainIndex() > it->outsideDomainIndex() && it->outsideDomainIndex() >= 0)
continue;
std::size_t outsideLabel;
if (it->outsideDomainIndex() < 0)
outsideLabel = std::numeric_limits<std::size_t>::max();
else
outsideLabel = it->outsideDomainIndex();
result.interfaceToVolume[std::make_pair(it->insideDomainIndex(), outsideLabel)] +=
SubtriangulationStatisticsDetail::computeVolume(it->geometry(), 2);
}
}
return result;
}
}
#endif // DUNEURO_SUBTRIANGULATION_STATISTICS_HH
......@@ -40,28 +40,6 @@ namespace duneuro
}
return result;
}
std::map<std::size_t, double> domainVolumes() const
{
Dune::PDELab::UnfittedSubTriangulation<typename STTraits::GridView> ust(fundamentalGridView_,
*subTriangulation_);
std::map<std::size_t, double> volume;
for (const auto& element : Dune::elements(fundamentalGridView_)) {
ust.create(element);
for (const auto& part : ust) {
const auto& geometry = part.geometry();
const auto& quadRule =
Dune::QuadratureRules<typename STTraits::GridView::ctype, dim>::rule(
geometry.type(), geometry.type().isCube() ? 2 : 0);
for (const auto& quadPoint : quadRule) {
volume[part.domainIndex()] +=
quadPoint.weight() * geometry.integrationElement(quadPoint.position());
}
}
}
return volume;
}
private:
UnfittedMEEGDriverData<dim> data_;
std::unique_ptr<typename STTraits::Grid> grid_;
......
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