Commit 3cfa5da5 authored by Sophia Müller's avatar Sophia Müller

continue CEM

parent 4a5f5bb5
......@@ -104,6 +104,10 @@ namespace duneuro
problem_.setElectrodePatches(electrodePatches);
}
void setImpedances(const std::vector<double>& impedances){
problem_.setImpedances(impedances);
}
private:
std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
......
......@@ -46,16 +46,13 @@ namespace duneuro
Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type
bctype(const typename Traits::IntersectionType& intersectionType,
const typename Traits::IntersectionDomainType& intersectionDomainType) const
{/*
{
for(unsigned int i = 0; i < electrodePatches_.size(); ++i){
if (electrodePatches_[i]->contains(intersectionType))
{std::cout<<"elementPatch" <<i<<" enthält Intersection"<< std::endl;
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Outflow;
break;}
else
{std::cout<<"elementPatch" <<i<<" enthält Intersection nicht"<< std::endl;
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann;}
}*/
}
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann;
}
......@@ -84,14 +81,43 @@ namespace duneuro
return 0.0;
}
typename Traits::RangeFieldType k(const typename Traits::ElementType&,
const typename Traits::DomainType&) const
{
return 0.0;
}
typename Traits::RangeFieldType o(const typename Traits::ElementType&,
const typename Traits::DomainType&) const
{
return 0.0;
}
//impedances
typename Traits::RangeFieldType z(const typename Traits::IntersectionType& intersectionType,
const typename Traits::IntersectionDomainType& intersectionDomainType)
{
for(unsigned int i = 0; i < electrodePatches_.size(); ++i){
if (electrodePatches_[i]->contains(intersectionType))
{
return impedances_[i];
break;}
}
}
void setElectrodePatches(
std::vector<std::shared_ptr<typename duneuro::SurfacePatch<GV>>>& electrodePatches) {
electrodePatches_ = electrodePatches;
}
void setImpedances(const std::vector<double>& impedances){
impedances_ =impedances;
}
private:
std::shared_ptr<const VC> volumeConductor_;
std::vector<std::shared_ptr<typename duneuro::SurfacePatch<GV>>> electrodePatches_;
std::vector<double> impedances_;
};
}
......
......@@ -11,6 +11,7 @@
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <duneuro/common/element_neighborhood_map.hh>
#include <duneuro/eeg/electrode_projection_interface.hh>
namespace duneuro
{
......@@ -71,6 +72,20 @@ namespace duneuro
}
}
//for creating SurfacePatch
template<typename ElementSearch>
ElementPatch(std::shared_ptr<ElementNeighborhoodMap<GV>> elementNeighborhoodMap,
const ElementSearch& elementSearch, const Coordinate& position, Element element,
ElementPatchInitialization initialization,
std::function<bool(Element)> elementFilter = [](const Element&) { return true; })
: elementNeighborhoodMap_(elementNeighborhoodMap)
, elementFilter_(elementFilter)
, elementMapper_(elementNeighborhoodMap_->gridView())
, vertexMapper_(elementNeighborhoodMap_->gridView())
{
initializeClosestVertex(elementSearch, position, element);
}
void extend(ElementPatchExtension extension)
{
std::vector<Element> candidates;
......@@ -216,6 +231,37 @@ namespace duneuro
}
}
}
template <typename ElementSearch>
void initializeClosestVertex(const ElementSearch& elementSearch, const Coordinate& position, Element& element)
{
const auto& geo = element.geometry();
// find closest corner
unsigned int minCorner = 0;
double minDistance = std::numeric_limits<double>::max();
for (unsigned int i = 0; i < geo.corners(); ++i) {
Coordinate tmp = position;
tmp -= geo.corner(i);
double tn = tmp.two_norm();
if (tn < minDistance) {
minDistance = tn;
minCorner = i;
}
}
// retrieve elements belonging to that corner
std::vector<Element> candidates;
elementNeighborhoodMap_->getNeighborsOfVertex(
vertexMapper_.subIndex(element, minCorner, GV::dimension),
std::back_inserter(candidates));
// filter them and push them to the list
for (const auto& e : candidates) {
if (elementFilter_(e)) {
elements_.push_back(e);
elementIndices_.insert(elementMapper_.index(e));
}
}
}
};
template <class VC, class ES>
......
......@@ -11,6 +11,7 @@
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <duneuro/common/element_patch.hh>
#include <duneuro/eeg/electrode_projection_interface.hh>
//////////////////////////////////////////////////////////////////////////
//Class for creating a SurfacePatch //
//a SurfacePatch is needed for the Complete Electode Model (CEM) //
......@@ -52,6 +53,17 @@ namespace duneuro
intersections_=elementpatch_.extractGlobalBoundaryIntersections();
}
template<typename ES>
SurfacePatch(std::shared_ptr<ElementNeighborhoodMap<GV>> elementNeighborhoodMap,
const ES& elementSearch, const Coordinate& position, const Element& element,
ElementPatchInitialization initialization)
: elementpatch_(elementNeighborhoodMap,elementSearch, position, element, initialization)
, position_(position)
{
intersections_=elementpatch_.extractGlobalBoundaryIntersections();
}
const std::vector<Intersection>& intersections() const
{
return intersections_;
......@@ -91,12 +103,12 @@ namespace duneuro
}
bool contains(const Element& element) const
{if(std::find(elements_.begin(),elements_.end(),element)==elements_.end())
{if(std::find(elements_.begin(),elements_.end(),element)!=elements_.end())
{return true;}else{return false;}
};
bool contains(const Intersection& intersection) const
{if(std::find(intersections_.begin(),intersections_.end(),intersection)==intersections_.end())
{if(std::find(intersections_.begin(),intersections_.end(),intersection)!=intersections_.end())
{return true;}else{return false;}
};
......@@ -135,12 +147,12 @@ namespace duneuro
template <class VC, class ES>
std::shared_ptr<SurfacePatch<typename VC::GridView>> make_electrode_patch(std::shared_ptr<ElementNeighborhoodMap<typename VC::GridView>> elementNeighborhoodMap,
const ES& elementSearch, const Dune::FieldVector<typename VC::ctype, VC::dim>& position, const double radius)
const ES& elementSearch, const Dune::FieldVector<typename VC::ctype, VC::dim>& position, const typename VC::GridView::template Codim<0>::Entity& element,const double radius)
{
const char* initialization = "closest_vertex";
auto patchInitialization = duneuro::elementPatchInitializationFromString(initialization);
auto surfacePatch = Dune::Std::make_unique<SurfacePatch<typename VC::GridView>>(
elementNeighborhoodMap, elementSearch, position,
elementNeighborhoodMap, elementSearch, position, element,
patchInitialization);
const char* extension = "vertex";
auto patchExtension = duneuro::elementPatchExtensionFromString(extension);
......
......@@ -19,10 +19,13 @@ namespace duneuro
template <class GV>
struct ElectrodeProjectionInterface {
using GlobalCoordinate = Dune::FieldVector<typename GV::ctype, GV::dimension>;
// using Element = typename GV::template Codim<0>::Entity;
virtual void setElectrodes(const std::vector<GlobalCoordinate>& electrodes) = 0;
virtual const ProjectedElectrode<GV>& getProjection(std::size_t i) const = 0;
virtual std::size_t size() const = 0;
//virtual const Element& element(std::size_t i) const =0;
virtual ~ElectrodeProjectionInterface()
{
......
......@@ -53,6 +53,7 @@ namespace duneuro
using SourceModelFactoryType = CGSourceModelFactory;
static void setElectrodes(SolverType& solver,
const std::vector<Dune::FieldVector<typename VC::ctype, VC::dim>>& GlobelElectrodes,
const ElectrodeProjectionInterface<typename VC::GridView>& electrodeProjection,
const std::vector<double>& impedances,
const std::shared_ptr<KDTreeElementSearch<typename VC::GridView>> elementSearch,
const typename VC::GridView& gridView,
......@@ -67,23 +68,30 @@ namespace duneuro
using SourceModelFactoryType = CGSourceModelFactory;
static void setElectrodes(SolverType& solver,
const std::vector<Dune::FieldVector<typename VC::ctype, VC::dim>>& GlobalElectrodes,
const ElectrodeProjectionInterface<typename VC::GridView>& electrodeProjection,
const std::vector<double>& impedances,
const std::shared_ptr<KDTreeElementSearch<typename VC::GridView>> elementSearch,
const typename VC::GridView& gridView,
const Dune::ParameterTree& config) {
if(impedances.size() != electrodeProjection.size()){
DUNE_THROW(Dune::Exception,
"number of impedances (" << impedances.size() << ") does not match number of electrodes ("
<< electrodeProjection.size() << ")");
}
std::vector<std::shared_ptr<typename duneuro::SurfacePatch<typename VC::GridView>>> electrodePatches;
auto radius = std::stod(config.get<std::string>("radius"));
std::cout<<radius<<std::endl;
//for (unsigned int i = 0; i < GlobalElectrodes.size(); ++i) {
//Dune::FieldVector<typename VC::ctype, VC::dim> electrode = GlobalElectrodes[i];
Dune::FieldVector<typename VC::ctype, VC::dim> electrode = GlobalElectrodes[0];
for (unsigned int i = 0; i < GlobalElectrodes.size(); ++i) {
const auto& proj = electrodeProjection.getProjection(i);
auto ePatch = duneuro::make_electrode_patch<VC, KDTreeElementSearch<typename VC::GridView>>(std::make_shared<duneuro::ElementNeighborhoodMap<typename VC::GridView>>(
duneuro::ElementNeighborhoodMap<typename VC::GridView>(gridView)),
*elementSearch, electrode, radius);
*elementSearch, proj.element.geometry().global(proj.localPosition), proj.element, radius);
electrodePatches.push_back(ePatch);
//}
}
solver.setElectrodePatches(electrodePatches);
}
solver.setImpedances(impedances);
}
};
template <class VC, ElementType et, int degree>
......@@ -93,6 +101,7 @@ namespace duneuro
using SourceModelFactoryType = DGSourceModelFactory;
static void setElectrodes(SolverType& solver,
const std::vector<Dune::FieldVector<typename VC::ctype, VC::dim>>& GlobelElectrodes,
const ElectrodeProjectionInterface<typename VC::GridView>& electrodeProjection,
const std::vector<double>& impedances,
const std::shared_ptr<KDTreeElementSearch<typename VC::GridView>> elementSearch,
const typename VC::GridView& gridView,
......@@ -229,8 +238,9 @@ namespace duneuro
const auto& proj = electrodeProjection_->getProjection(i);
projectedGlobalElectrodes_.push_back(proj.element.geometry().global(proj.localPosition));
}
//creating electrode Patches
SelectFittedSolver<solverType,typename Traits::VC,elementType, degree>::setElectrodes( *solver_, projectedGlobalElectrodes_,impedances, elementSearch_, volumeConductorStorage_.get()->gridView(), config);
SelectFittedSolver<solverType,typename Traits::VC,elementType, degree>::setElectrodes( *solver_, projectedGlobalElectrodes_,*electrodeProjection_,impedances, elementSearch_, volumeConductorStorage_.get()->gridView(), config);
}
virtual void
......
......@@ -48,20 +48,21 @@ int main(int argc, char** argv)
//extension of surfacepatch
surfacepatch.extend(patchExtension, radius);
//try if the right surfaces are in the extension
/* for(const auto& is:surfacepatch.intersections())
for(const auto& is:surfacepatch.intersections())
{ auto fgeo = is.geometry();
std::cout<< "SurfaceCenter: " << fgeo.center()<<std::endl;
auto boundaryelement = is.inside();
/* auto boundaryelement = is.inside();
auto beGeo = boundaryelement.geometry();
std::cout<< "BoundaryElement: " <<std::endl;
for (unsigned j=0; j<beGeo.corners();j++){
std::cout<< "Vertices: " << beGeo.corner(j) << std::endl;}
std::cout<< "Vertices: " << beGeo.corner(j) << std::endl;}*/
}
for(const auto& e:surfacepatch.elements())
{ auto egeo = e.geometry();
std::cout<< "BoundaryElement: " << egeo.center()<<std::endl;}*/
std::cout<< "BoundaryElement: " << egeo.center()<<std::endl;}
......@@ -77,20 +78,17 @@ for(const auto& e:surfacepatch.elements())
//test of mathod make_surface_patch
Dune::ParameterTree config;
config["extensions"] = "intersection";
config["initialization"] = "closestVertex";
//auto surfacepatch2 = duneuro::make_surface_patch(/*enm,*/ elementSearch/*, position, config, radius*/);
/*
duneuro::SurfacePatchVTKFunction<SP> vtksurfacepatch2function(surfacepatch2, "surfacepatch2");
auto celldata = std::make_shared<duneuro::SurfacePatchVTKFunction<SP>>(vtksurfacepatch2function);
Dune::VTKWriter<GV> vtkwritersurface2(gridView);
vtkwritersurface2.addCellData(celldata);
vtkwritersurface2.write("surface2_grid");*/
}
//try for-loop for finding out if intersection is in surfacepatch
// iterate over all cells in the LeafGridView
for (const auto& e : elements(gridView))
{for (const auto& is: intersections(gridView, e)){
// std::cout << is.geometry().center() << std::endl;
if(surfacepatch.contains(is)){
std::cout<<"in"<<std::endl;}
}
}
}
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