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

[VolumeConductor] use always indirect mapping

In order to prepare for providing statistics about the volume conductor,
we remove the direct mapping to be able to access element labels.
Most of the meshes (but some legacy files) provide labeled elements and
we thus have a direct access to these labels.
For the legacy formats, this leads to a slight increase in time consumption
when accessing the tensor of an element, but compared to other tasks, this
should be negligible.
Besides being able to access the labels, this also gets rid of the type
erasure in the volume conductor, leading to way more readable code.
parent ef4d2970
......@@ -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_;
};
}
......
......@@ -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 << "\"");
}
......
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