Commit 844f5cc3 authored by Andreas Nüßing's avatar Andreas Nüßing

[MEEGDriver] add statistics method

We add a new method to the MEEGDriver interface, namely `statistics`.
The exact semantics of the method depends on the implementation, but
it should write out statistics about the driver to the datatree.
For the fitted driver, this currently means the volume and surface area
of the different discrete compartments.
For the unfitted driver, the semantics is similar as it writes out the
volume and surface area of the different domains.
parent c97477a5
#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
......@@ -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