...
 
Commits (8)
# Where to search and which files to use
INPUT += @top_srcdir@/duneuro
EXTRA_PACKAGES += amsmath
#---------------------------------------------------------------------------
# Generate help files for Qt Creator
#---------------------------------------------------------------------------
GENERATE_QHP = YES
QCH_FILE = "../qthelp/qthelp.qch"
QHP_NAMESPACE = "duneuro"
QHP_VIRTUAL_FOLDER = "dune"
QHG_LOCATION = "/usr/bin/qhelpgenerator"
#ifndef DUNEURO_CG_CEM_SOLVER_HH
#define DUNEURO_CG_CEM_SOLVER_HH
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/localoperator/convectiondiffusionfem.hh>
#include <duneuro/common/assembler.hh>
#include <duneuro/common/convection_diffusion_FEM_CEM.hh>
#include <duneuro/common/convection_diffusion_cg_default_parameter.hh>
#include <duneuro/common/convection_diffusion_cg_cem_default_parameter.hh>
#include <duneuro/common/flags.hh>
#include <duneuro/common/kdtree.hh>
#include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/make_dof_vector.hh>
#include <duneuro/common/random.hh>
#include <duneuro/io/data_tree.hh>
#include <duneuro/common/element_neighborhood_map.hh>
namespace duneuro
{
template <class VC, ElementType elementType, unsigned int degree, class DF = double,
class RF = double, class JF = double>
struct CGCEMSolverTraits {
static const int dimension = VC::dim;
using VolumeConductor = VC;
using GridView = typename VC::GridView;
using CoordinateFieldType = typename VC::ctype;
using ElementSearch = KDTreeElementSearch<GridView>;
using ElementNeighborhoodMap2 = ElementNeighborhoodMap<GridView>;
using Problem = ConvectionDiffusionCGCEMDefaultParameter<VC>;
using DirichletExtension = Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem>;
using BoundaryCondition = Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem>;
using FunctionSpace =
typename CGFunctionSpaceTraits<VC, BoundaryCondition, degree, elementType>::Type;
using DomainDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, DF>;
using RangeDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, RF>;
using LocalOperator =
ConvectionDiffusionFEMCEM<Problem, typename FunctionSpace::FEM>;
using Assembler = GalerkinGlobalAssembler<FunctionSpace, LocalOperator, DF, RF, JF>;
using LinearSolver =
LinearProblemSolver<typename Assembler::GO, DomainDOFVector, RangeDOFVector>;
};
template <class VC, ElementType elementType, unsigned int degree, class DF = double,
class RF = double, class JF = double>
class CGCEMSolver
{
public:
using Traits = CGCEMSolverTraits<VC, elementType, degree, DF, RF, JF>;
CGCEMSolver(std::shared_ptr<const VC> volumeConductor,
std::shared_ptr<const typename Traits::ElementSearch> search,
const Dune::ParameterTree& config, DataTree dataTree = DataTree())
: volumeConductor_(volumeConductor)
, search_(search)
, problem_(volumeConductor)
, dirichletExtension_(volumeConductor->gridView(), problem_)
, boundaryCondition_(volumeConductor->gridView(), problem_)
, functionSpace_(const_cast<typename VC::GridType&>(volumeConductor->grid()),
boundaryCondition_)
// const_cast due to misused non-const reference in pdelab boilerplate
, localOperator_(problem_, config.get<unsigned int>("intorderadd", 0))
, assembler_(functionSpace_, localOperator_, elementType == ElementType::hexahedron ?
(1 << VC::dim) + 1 :
Dune::StaticPower<3, VC::dim>::power)
, linearSolver_(assembler_.getGO(), config)
{
dataTree.set("degree", degree);
dataTree.set("element_type", to_string(elementType));
}
/*
///neuer Konstruktor um zu testen, ob man die Elektroden mitnehmen kann
CGCEMSolver(std::shared_ptr<const VC> volumeConductor,
std::shared_ptr<const typename Traits::ElementSearch> search,
const Dune::ParameterTree& config, DataTree dataTree = DataTree(), int test=1)
: volumeConductor_(volumeConductor)
, search_(search)
, problem_(volumeConductor)
, dirichletExtension_(volumeConductor->gridView(), problem_)
, boundaryCondition_(volumeConductor->gridView(), problem_)
, functionSpace_(const_cast<typename VC::GridType&>(volumeConductor->grid()),
boundaryCondition_)
// const_cast due to misused non-const reference in pdelab boilerplate
, localOperator_(problem_, config.get<unsigned int>("intorderadd", 0))
, assembler_(functionSpace_, localOperator_, elementType == ElementType::hexahedron ?
(1 << VC::dim) + 1 :
Dune::StaticPower<3, VC::dim>::power)
, linearSolver_(assembler_.getGO(), config)
,test_(test)
{
dataTree.set("degree", degree);
dataTree.set("element_type", to_string(elementType));
/////////////////////////////////////////////
///implementing additional part for CEM ///
/////////////////////////////////////////////
/// reading electrodes
}*/
template <typename SolverBackend>
void solve(SolverBackend& solverBackend, const typename Traits::RangeDOFVector& rightHandSide,
typename Traits::DomainDOFVector& solution, const Dune::ParameterTree& config,
DataTree dataTree = DataTree())
{
Dune::Timer timer;
randomize_uniform(Dune::PDELab::Backend::native(solution), DF(-1.0), DF(1.0));
linearSolver_.apply(solverBackend, solution, rightHandSide, config, dataTree);
dataTree.set("time", timer.elapsed());
}
const typename Traits::FunctionSpace& functionSpace() const
{
return functionSpace_;
}
std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor() const
{
return volumeConductor_;
}
std::shared_ptr<const typename Traits::ElementSearch> elementSearch() const
{
return search_;
}
private:
std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
typename Traits::Problem problem_;
typename Traits::DirichletExtension dirichletExtension_;
typename Traits::BoundaryCondition boundaryCondition_;
typename Traits::FunctionSpace functionSpace_;
typename Traits::LocalOperator localOperator_;
typename Traits::Assembler assembler_;
typename Traits::LinearSolver linearSolver_;
template <class V>
friend struct MakeDOFVectorHelper;
};
}
#endif // DUNEURO_CG_CEM_SOLVER_HH
#ifndef DUNEURO_CONVECTION_DIFFUSION_FEM_CEM_HH
#define DUNEURO_CONVECTION_DIFFUSION_FEM_CEM_HH
#include <dune/pdelab/localoperator/convectiondiffusionfem.hh>
//#include <duneuro/common/convection_diffusion_cg_cem_default_parameter.hh>
#include<vector>
#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/common/referenceelements.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/idefault.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/finiteelement/localbasiscache.hh>
#include<dune/pdelab/localoperator/convectiondiffusionparameter.hh>
#include <duneuro/common/surface_patch.hh>
namespace duneuro{
template<typename T, typename FiniteElementMap>
class ConvectionDiffusionFEMCEM :
public Dune::PDELab::ConvectionDiffusionFEM<T,FiniteElementMap> {
public:
ConvectionDiffusionFEMCEM (T& param_, int intorderadd_=0)
: Dune::PDELab::ConvectionDiffusionFEM<T,FiniteElementMap>(param_,intorderadd_), param(param_), intorderadd(intorderadd_)
{ std::cout<<"das ist mit FEM_CEM!!!!"<<std::endl;
}
// jacobian of volume term
template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
M& mat) const
{
// Define types
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
using size_type = typename LFSU::Traits::SizeType;
// dimensions
const int dim = EG::Entity::dimension;
// Reference to cell
const auto& cell = eg.entity();
// Get geometry
auto geo = eg.geometry();
// evaluate diffusion tensor at cell center, assume it is constant over elements
auto ref_el = referenceElement(geo);
auto localcenter = ref_el.position(0,0);
auto tensor = param.A(cell,localcenter);
// Initialize vectors outside for loop
std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
// Transformation matrix
typename EG::Geometry::JacobianInverseTransposed jac;
// loop over quadrature points
auto intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
for (const auto& ip : quadratureRule(geo,intorder))
{
// evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
auto& js = cache.evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
// transform gradient to real element
jac = geo.jacobianInverseTransposed(ip.position());
for (size_type i=0; i<lfsu.size(); i++)
{
jac.mv(js[i][0],gradphi[i]);
tensor.mv(gradphi[i],Agradphi[i]);
}
// evaluate basis functions
auto& phi = cache.evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
// evaluate velocity field, sink term and source term
auto b = param.b(cell,ip.position());
auto c = param.c(cell,ip.position());
// integrate (A grad phi_j)*grad phi_i - phi_j b*grad phi_i + c*phi_j*phi_i
RF factor = ip.weight() * geo.integrationElement(ip.position());
for (size_type j=0; j<lfsu.size(); j++)
for (size_type i=0; i<lfsu.size(); i++)
mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i]-phi[j]*(b*gradphi[i])+c*phi[j]*phi[i] )*factor);
}
//calculate second term for CEM:
//latex: \Sum_{l=1}^{L} ( 1/(Z_l)*\integral_{e_l}\phi_i\phi_j dS )
}
private:
T& param;
int intorderadd;
using LocalBasisType = typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType;
Dune::PDELab::LocalBasisCache<LocalBasisType> cache;
//typename Dune::PDELab::ConvectionDiffusionFEM< T, FiniteElementMap> ConvectionDiffusionFEM;
};
}
#endif // DUNEURO_CONVECTION_DIFFUSION_FEM_CEM_HH
#ifndef CONVECTION_DIFFUSION_CG_CEM_DEFAULT_PARAMETER_HH
#define CONVECTION_DIFFUSION_CG_CEM_DEFAULT_PARAMETER_HH
#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
namespace duneuro
{
template <typename VC>
class ConvectionDiffusionCGCEMDefaultParameter
{
public:
using GV = typename VC::GridType::LeafGridView;
using Traits = Dune::PDELab::ConvectionDiffusionParameterTraits<GV, typename GV::ctype>;
explicit ConvectionDiffusionCGCEMDefaultParameter(std::shared_ptr<const VC> volumeConductor)
: volumeConductor_(volumeConductor)
{
assert(volumeConductor_);
std::cout<<"das ist cg_cem_default_parameter"<<std::endl;
}
template <class EG>
typename Traits::PermTensorType A(const EG& e, const typename Traits::DomainType& x) const
{
return A(e.entity(), x);
}
typename Traits::PermTensorType A(const typename Traits::ElementType& e,
const typename Traits::DomainType&) const
{
return volumeConductor_->tensor(e);
}
template <class EG>
typename Traits::RangeType b(const EG&, const typename Traits::DomainType&) const
{
return typename Traits::RangeType(0.0);
}
template <class EG>
typename Traits::RangeFieldType c(const EG&, const typename Traits::DomainType&) const
{
return typename Traits::RangeFieldType(0.0);
}
Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type
bctype(const typename Traits::IntersectionType&,
const typename Traits::IntersectionDomainType&) const
{
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann;
}
typename Traits::RangeFieldType o(const typename Traits::IntersectionType&,
const typename Traits::IntersectionDomainType&)
{
return 0.0;
}
typename Traits::RangeFieldType f(const typename Traits::ElementType&,
const typename Traits::DomainType&) const
{
return 0.0;
}
typename Traits::RangeFieldType g(const typename Traits::ElementType&,
const typename Traits::DomainType&) const
{
return 0.0;
}
typename Traits::RangeFieldType j(const typename Traits::IntersectionType&,
const typename Traits::IntersectionDomainType&) const
{
return 0.0;
}
private:
std::shared_ptr<const VC> volumeConductor_;
};
}
#endif // CONVECTION_DIFFUSION_CG_CEM_DEFAULT_PARAMETER_HH
......@@ -147,6 +147,25 @@ namespace duneuro
return out;
}
template <class I>
void extractGlobalBoundaryIntersections(I out) const
{
for (const auto& element : elements_) {
for (const auto& is : Dune::intersections(elementNeighborhoodMap_->gridView(), element)) {
if (is.boundary()) {
*out++ = is;
}
}
}
}
std::vector<Intersection> extractGlobalBoundaryIntersections() const
{
std::vector<Intersection> out;
extractGlobalBoundaryIntersections(std::back_inserter(out));
return out;
}
private:
std::shared_ptr<ElementNeighborhoodMap<GV>> elementNeighborhoodMap_;
std::function<bool(Element)> elementFilter_;
......@@ -297,10 +316,10 @@ namespace duneuro
};
template <class EP>
std::unique_ptr<ElementPatchVTKFunction<EP>>
std::shared_ptr<ElementPatchVTKFunction<EP>>
make_element_patch_vtk_function(const EP& patch, const std::string& name = "patch")
{
return Dune::Std::make_unique<ElementPatchVTKFunction<EP>>(patch, name);
return std::make_shared<ElementPatchVTKFunction<EP>>(patch, name);
}
}
......
......@@ -30,7 +30,7 @@ namespace duneuro
static const Dune::GeometryType::BasicType value = Dune::GeometryType::BasicType::simplex;
};
enum class FittedSolverType { cg, dg };
enum class FittedSolverType { cg, dg, cg_cem };
enum class UnfittedSolverType { cutfem, udg};
}
......
#ifndef DUNEURO_SURFACE_PATCH_HH
#define DUNEURO_SURFACE_PATCH_HH
#include <memory>
#include <vector>
#include <dune/common/parametertree.hh>
#include <dune/common/std/memory.hh>
#include <dune/grid/common/scsgmapper.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <duneuro/common/element_patch.hh>
//////////////////////////////////////////////////////////////////////////
//Class for creating a SurfacePatch //
//a SurfacePatch is needed for the Complete Electode Model (CEM) //
//creates a boundary surface patch for electrode, all boundarysurfaces //
//are included, which have a smaller/equal distance to the //
//electrode-projection //
//////////////////////////////////////////////////////////////////////////
// required parameters: //
// GridView cp. Element_Patch //
// position (of electrode point) Vector for electrode //
// radius of electrode int, for creating the Surface //
// ElementNeighborhoodMap cp. Element_Patch //
// elementSearch cp. Element_Patch //
// initialization two options: closest_vertex or //
// single_element ( Element_Patch) //
//////////////////////////////////////////////////////////////////////////
namespace duneuro
{
template <typename GV>
class SurfacePatch
{
public:
using GridView = GV;
using Coordinate = Dune::FieldVector<typename GV::ctype, GV::dimension>;
using Element = typename GV::template Codim<0>::Entity;
using Intersection = typename GV::Intersection;
using IntersectionMapper = Dune::SingleCodimSingleGeomTypeMapper<GV,1>;
using EP = typename duneuro::ElementPatch<GV>;
template<typename ES>
SurfacePatch(std::shared_ptr<ElementNeighborhoodMap<GV>> elementNeighborhoodMap,
const ES& elementSearch, const Coordinate& position,
ElementPatchInitialization initialization)
: elementpatch_(elementNeighborhoodMap,elementSearch, position, initialization)
, position_(position)
{
intersections_=elementpatch_.extractGlobalBoundaryIntersections();
}
const std::vector<Intersection>& intersections() const
{
return intersections_;
}
const std::vector<Element>& elements() const
{
return elements_;
}
//Extend and form surfacepatch, that the center of every Intersection is in the radius
//for ElectrodeSurface
void extend(ElementPatchExtension extension, const double& radius){
std::vector<Intersection> surface;
for (const auto& is:intersections_){
if(rfilter(is, radius)){
surface.push_back(is);
}
}
std::vector<Intersection> old;
extend(extension);
while(old.size() < surface.size())
{old = surface;
extend(extension);
for (const auto& is:intersections_){
if(rfilter(is, radius) && ifilter(is,surface) ){
surface.push_back(is);
}
}
}
intersections_=surface;
elements_.clear();
for(const auto& is:intersections_)
{ elements_.push_back(is.inside());}
//std::cout<< "größe Surface-Vektor: " << surface.size() << std::endl;
}
bool contains(const Element& element) const
{if(std::find(elements_.begin(),elements_.end(),element)==elements_.end())
{return true;}else{return false;}
};
private:
std::vector<Intersection> intersections_;
std::set<std::size_t> intersectionIndices_;
EP elementpatch_;
std::vector<Element> elements_;
const Coordinate position_;
//PatchExtension
void extend(ElementPatchExtension extension){
elementpatch_.extend(extension);
elements_=elementpatch_.elements();
intersections_ = elementpatch_.extractGlobalBoundaryIntersections();
}
//filter to extract just the intersections, whitch center is in radius distance to beginning position
bool rfilter(Intersection is, const double& radius)
{ auto isgeo = is.geometry();
auto dist0 = isgeo.center()-position_;
double dist= dist0.two_norm();
if(dist<=radius){return true;}else{return false;}
};
//filter to extract all intersections, that are already in vector
bool ifilter(Intersection is, std::vector<Intersection> vector)
{if(std::find(vector.begin(),vector.end(),is)==vector.end())
{return true;}else{return false;}
};
};
template <class VC, class ES>
std::unique_ptr<SurfacePatch<typename VC::GridView>> make_surface_patch(
std::shared_ptr<ElementNeighborhoodMap<typename VC::GridView>> elementNeighborhoodMap,
const ES& elementSearch, const Dune::FieldVector<typename VC::ctype, VC::dim>& position,
const Dune::ParameterTree& config)
{ const double radius = std::stod(config.get<std::string>("radius"));
const char* initialization = "closest_vertex";
auto patchInitialization = duneuro::elementPatchInitializationFromString(initialization);
auto surfacePatch = Dune::Std::make_unique<SurfacePatch<typename VC::GridView>>(
elementNeighborhoodMap, elementSearch, position,
patchInitialization);
const char* extension = "vertex";
auto patchExtension = duneuro::elementPatchExtensionFromString(extension);
surfacePatch->extend(patchExtension,radius);
return surfacePatch;
}
//write to VTKFunktion
template <class SP>
class SurfacePatchVTKFunction : public Dune::VTKFunction<typename SP::GridView>
{
using GV = typename SP::GridView;
typedef typename GV::ctype DF;
enum { n = GV::dimension };
typedef typename GV::template Codim<0>::Entity Entity;
public:
SurfacePatchVTKFunction(const SP& sp, std::string name) : sp_(sp), name_(name)
{
}
virtual int ncomps() const
{
return 1;
}
virtual double evaluate(int comp, const Entity& e, const Dune::FieldVector<DF, n>& xi) const
{
return sp_.contains(e) ? 0 : 1;
}
virtual std::string name() const
{
return name_;
}
private:
const SP& sp_;
std::string name_;
};
template <class SP>
std::shared_ptr<ElementPatchVTKFunction<SP>>
make_surface_patch_vtk_function(const SP& patch, const std::string& name = "patch")
{
return std::make_shared<SurfacePatchVTKFunction<SP>>(patch, name);
}
}
#endif // DUNEURO SURFACE_PATCH_HH
......@@ -8,10 +8,12 @@
#include <dune/common/std/memory.hh>
#include <duneuro/common/cg_solver.hh>
#include <duneuro/common/cg_cem_solver.hh>
#include <duneuro/common/cg_solver_backend.hh>
#include <duneuro/common/default_grids.hh>
#include <duneuro/common/dg_solver.hh>
#include <duneuro/common/dg_solver_backend.hh>
#include <duneuro/common/element_neighborhood_map.hh>
#include <duneuro/common/flags.hh>
#if HAVE_DUNE_SUBGRID
#include <duneuro/common/geometry_adaption.hh>
......@@ -37,6 +39,7 @@
#include <duneuro/meg/fitted_meg_transfer_matrix_solver.hh>
#include <duneuro/meg/meg_solver_factory.hh>
#include <duneuro/meg/meg_solver_interface.hh>
#include <duneuro/common/surface_patch.hh>
namespace duneuro
{
......@@ -50,6 +53,13 @@ namespace duneuro
using SourceModelFactoryType = CGSourceModelFactory;
};
template <class VC, ElementType et, int degree>
struct SelectFittedSolver<FittedSolverType::cg_cem, VC, et, degree> {
using SolverType = CGCEMSolver<VC, et, degree>;
using SolverBackendType = CGSolverBackend<SolverType, et>;
using SourceModelFactoryType = CGSourceModelFactory;
};
template <class VC, ElementType et, int degree>
struct SelectFittedSolver<FittedSolverType::dg, VC, et, degree> {
using SolverType = DGSolver<VC, et, degree>;
......@@ -70,6 +80,7 @@ namespace duneuro
using TransferMatrixRHSFactory = FittedTransferMatrixRHSFactory;
using DomainDOFVector = typename Solver::Traits::DomainDOFVector;
using ElementSearch = KDTreeElementSearch<typename VC::GridView>;
using GV = typename VC::GridView;
};
template <int dim, ElementType elementType, FittedSolverType solverType, int degree,
......@@ -170,6 +181,7 @@ namespace duneuro
return Dune::Std::make_unique<Function>(make_domain_dof_vector(*solver_, 0.0));
}
//if statement for creating an electrodePatch for cg_cem
virtual void
setElectrodes(const std::vector<typename MEEGDriverInterface<dim>::CoordinateType>& electrodes,
const Dune::ParameterTree& config) override
......@@ -183,6 +195,28 @@ namespace duneuro
const auto& proj = electrodeProjection_->getProjection(i);
projectedGlobalElectrodes_.push_back(proj.element.geometry().global(proj.localPosition));
}
auto solverType1 = config_.get<std::string>("solver_type");
//creating electrode Patches
if (solverType1 == "cg_cem")
{
electrodePatches_.clear();
for (unsigned int i = 0; i < projectedGlobalElectrodes_.size(); ++i) {
const auto& ePatch =duneuro::make_surface_patch<typename Traits::VC, typename Traits::ElementSearch>
(std::make_shared<duneuro::ElementNeighborhoodMap<typename Traits::GV>>(
duneuro::ElementNeighborhoodMap<typename Traits::GV>(volumeConductorStorage_.get()->gridView())),
*elementSearch_,
projectedGlobalElectrodes_[i],
config);
electrodePatches_.push_back(*ePatch);
}
auto surfacePatch = duneuro::make_surface_patch<typename Traits::VC, typename Traits::ElementSearch>
(std::make_shared<duneuro::ElementNeighborhoodMap<typename Traits::GV>>(
duneuro::ElementNeighborhoodMap<typename Traits::GV>(volumeConductorStorage_.get()->gridView())),
*elementSearch_,
projectedGlobalElectrodes_[0],
config);
}
}
virtual void setCoilsAndProjections(
......@@ -441,6 +475,7 @@ namespace duneuro
std::vector<typename duneuro::ElectrodeProjectionInterface<
typename Traits::VC::GridView>::GlobalCoordinate>
projectedGlobalElectrodes_;
std::vector<typename duneuro::SurfacePatch<typename Traits::GV>> electrodePatches_;
};
}
......
......@@ -11,6 +11,8 @@ template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1>;
template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1>;
template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg_cem, 1>;
#if HAVE_DUNE_SUBGRID
template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, true>;
......
......@@ -22,6 +22,8 @@ extern template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::tetrahe
duneuro::FittedSolverType::dg, 1, false>;
extern template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, false>;
extern template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg_cem, 1, false>;
#if HAVE_DUNE_SUBGRID
extern template class duneuro::FittedMEEGDriver<3, duneuro::ElementType::hexahedron,
......@@ -210,6 +212,14 @@ namespace duneuro
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else if (solverType == "cg_cem") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedMEEGDriver<3, ElementType::tetrahedron,
FittedSolverType::cg_cem, 1>>(
data.fittedData, config, dataTree);
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else if (solverType == "dg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedMEEGDriver<3, ElementType::tetrahedron,
......@@ -233,6 +243,7 @@ namespace duneuro
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else {
DUNE_THROW(Dune::Exception, "unknown solver type \"" << solverType << "\"");
}
......
dune_add_test(SOURCES test_electrode_projection.cc LINK_LIBRARIES duneuro)
dune_add_test(SOURCES test_numerical_flux.cc)
dune_add_test(SOURCES test_physical_flux.cc)
dune_add_test(SOURCES test_boundary_intersection_patch.cc)
#include <config.h>
#include <iostream>
#include <memory>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/grid/utility/hierarchicsearch.hh>
#include <dune/istl/scalarproducts.hh>
#include <duneuro/common/element_patch.hh>
#include <duneuro/common/surface_patch.hh>
//Test for creating a Surface Patch
int main(int argc, char** argv)
{
std::cout<<"done"<<std::endl;
const int dim =3;
using Grid = Dune::YaspGrid<dim>;
auto grid = Dune::StructuredGridFactory<Grid>::createCubeGrid({0, 0, 0}, {1, 2, 2}, {4, 100, 100});
using GV = Grid::LeafGridView;
auto gridView = grid->leafGridView();
using Coordinate = Dune::FieldVector<typename GV::ctype, GV::dimension>;
//set Data for creating Patch
using LeafIndexSet = GV::IndexSet;
const LeafIndexSet& set = gridView.indexSet();
using ElementSearch = Dune::HierarchicSearch<Grid,LeafIndexSet>;
auto elementSearch = ElementSearch(*grid, set);
auto position = Coordinate(1);
auto enh0 = duneuro::ElementNeighborhoodMap<GV>(gridView);
const char* initialization = "closest_vertex";
auto patchInitialization = duneuro::elementPatchInitializationFromString(initialization);
const char* extension = "vertex";
auto patchExtension = duneuro::elementPatchExtensionFromString(extension);
auto enm = std::make_shared<duneuro::ElementNeighborhoodMap<GV>>(enh0);
//creating SurfacePatch
const double radius =0.3;
using SP = duneuro::SurfacePatch<GV>;
auto surfacepatch = SP( enm, elementSearch, position, patchInitialization);
//extension of surfacepatch
surfacepatch.extend(patchExtension, radius);
//try if the right surfaces are in the extension
/* for(const auto& is:surfacepatch.intersections())
{ auto fgeo = is.geometry();
std::cout<< "SurfaceCenter: " << fgeo.center()<<std::endl;
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;}
}
for(const auto& e:surfacepatch.elements())
{ auto egeo = e.geometry();
std::cout<< "BoundaryElement: " << egeo.center()<<std::endl;}*/
//create vtkFunction for using in VTKWriter
//creating SurfacePatch vtk
duneuro::SurfacePatchVTKFunction<SP> vtksurfacepatchfunction(surfacepatch, "surfacepatch");
auto celldata = std::make_shared<duneuro::SurfacePatchVTKFunction<SP>>(vtksurfacepatchfunction);
Dune::VTKWriter<GV> vtkwritersurface(gridView);
vtkwritersurface.addCellData(celldata);
vtkwritersurface.write("surface_grid");
//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");*/
}