Skip to content
Snippets Groups Projects
Commit 052c87c5 authored by Jö Fahlke's avatar Jö Fahlke
Browse files

ConstraintsAssemblerHelper: inline

parent f5a01aa3
No related branches found
No related tags found
No related merge requests found
......@@ -3,151 +3,114 @@
namespace PPS {
//! construct constraints
/**
* \code
* #include <dune/pdelab/constraints/common/constraints.hh>
* \endcode
* \tparam P Type implementing a constraints parameter tree
* \tparam GFS Type implementing the model GridFunctionSpace
* \tparam CG Type implementing the model
* GridFunctionSpace::ConstraintsContainer::Type
* \tparam isFunction bool to identify old-style parameters, which were implemented the Dune::PDELab::FunctionInterface
*/
template<typename P, typename GFS, typename CG, bool isFunction>
struct ConstraintsAssemblerHelper
template<typename P, typename GFS, typename CG>
void constraints(const P& p, const GFS& gfs, CG& cg,
const bool verbose = false)
{
//! construct constraints from given boundary condition function
/**
* \code
* #include <dune/pdelab/constraints/common/constraints.hh>
* \endcode
* \tparam P Type implementing a constraints parameter tree
* \tparam GFS Type implementing the model GridFunctionSpace
* \tparam CG Type implementing the model
* GridFunctionSpace::ConstraintsContainer::Type
*
* \param p The parameter object
* \param gfs The gridfunctionspace
* \param cg The constraints container
* \param verbose Print information about the constaints at the end
*/
static void
assemble(const P& p, const GFS& gfs, CG& cg, const bool verbose)
{
// get some types
using ES = typename GFS::Traits::EntitySet;
using Element = typename ES::Traits::Element;
using Intersection = typename ES::Traits::Intersection;
ES es = gfs.entitySet();
// clear global constraints
cg.clear();
using ES = typename GFS::Traits::EntitySet;
using Element = typename ES::Traits::Element;
using Intersection = typename ES::Traits::Intersection;
// make local function space
using LFS = Dune::PDELab::LocalFunctionSpace<GFS>;
LFS lfs_e(gfs);
Dune::PDELab::LFSIndexCache<LFS> lfs_cache_e(lfs_e);
LFS lfs_f(gfs);
Dune::PDELab::LFSIndexCache<LFS> lfs_cache_f(lfs_f);
ES es = gfs.entitySet();
// get index set
auto& is = es.indexSet();
// make local function space
using LFS = Dune::PDELab::LocalFunctionSpace<GFS>;
LFS lfs_e(gfs);
Dune::PDELab::LFSIndexCache<LFS> lfs_cache_e(lfs_e);
LFS lfs_f(gfs);
Dune::PDELab::LFSIndexCache<LFS> lfs_cache_f(lfs_f);
// loop once over the grid
for (const auto& element : elements(es))
{
// get index set
auto& is = es.indexSet();
auto id = is.uniqueIndex(element);
// loop once over the grid
for (const auto& element : elements(es))
{
// bind local function space to element
lfs_e.bind(element);
auto id = is.uniqueIndex(element);
using CL = typename CG::LocalTransformation;
// bind local function space to element
lfs_e.bind(element);
CL cl_self;
using CL = typename CG::LocalTransformation;
using ElementWrapper = Dune::PDELab::ElementGeometry<Element>;
using IntersectionWrapper = Dune::PDELab::IntersectionGeometry<Intersection>;
CL cl_self;
Dune::TypeTree::applyToTreePair(p,lfs_e,Dune::PDELab::VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
using ElementWrapper = Dune::PDELab::ElementGeometry<Element>;
using IntersectionWrapper = Dune::PDELab::IntersectionGeometry<Intersection>;
// iterate over intersections and call metaprogram
unsigned int intersection_index = 0;
for (const auto& intersection : intersections(es,element))
{
Dune::TypeTree::applyToTreePair(p,lfs_e,Dune::PDELab::VolumeConstraints<ElementWrapper,CL>(ElementWrapper(element),cl_self));
auto intersection_data = classifyIntersection(es,intersection);
auto intersection_type = std::get<0>(intersection_data);
auto& outside_element = std::get<1>(intersection_data);
// iterate over intersections and call metaprogram
unsigned int intersection_index = 0;
for (const auto& intersection : intersections(es,element))
{
switch (intersection_type) {
auto intersection_data = classifyIntersection(es,intersection);
auto intersection_type = std::get<0>(intersection_data);
auto& outside_element = std::get<1>(intersection_data);
case Dune::PDELab::IntersectionType::skeleton:
case Dune::PDELab::IntersectionType::periodic:
{
auto idn = is.uniqueIndex(outside_element);
switch (intersection_type) {
if(id > idn){
// bind local function space to element in neighbor
lfs_f.bind(outside_element);
case Dune::PDELab::IntersectionType::skeleton:
case Dune::PDELab::IntersectionType::periodic:
{
auto idn = is.uniqueIndex(outside_element);
CL cl_neighbor;
if(id > idn){
// bind local function space to element in neighbor
lfs_f.bind(outside_element);
Dune::TypeTree::applyToTreePair(lfs_e,lfs_f,Dune::PDELab::SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
CL cl_neighbor;
if (!cl_neighbor.empty())
{
lfs_cache_f.update();
cg.import_local_transformation(cl_neighbor,lfs_cache_f);
}
Dune::TypeTree::applyToTreePair(lfs_e,lfs_f,Dune::PDELab::SkeletonConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self,cl_neighbor));
if (!cl_neighbor.empty())
{
lfs_cache_f.update();
cg.import_local_transformation(cl_neighbor,lfs_cache_f);
}
break;
}
case Dune::PDELab::IntersectionType::boundary:
Dune::TypeTree::applyToTreePair(p,lfs_e,Dune::PDELab::BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
}
break;
}
case Dune::PDELab::IntersectionType::processor:
Dune::TypeTree::applyToTree(lfs_e,Dune::PDELab::ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
break;
case Dune::PDELab::IntersectionType::boundary:
Dune::TypeTree::applyToTreePair(p,lfs_e,Dune::PDELab::BoundaryConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
break;
}
++intersection_index;
}
case Dune::PDELab::IntersectionType::processor:
Dune::TypeTree::applyToTree(lfs_e,Dune::PDELab::ProcessorConstraints<IntersectionWrapper,CL>(IntersectionWrapper(intersection,intersection_index),cl_self));
break;
if (!cl_self.empty())
{
lfs_cache_e.update();
cg.import_local_transformation(cl_self,lfs_cache_e);
}
++intersection_index;
}
if (!cl_self.empty())
{
lfs_cache_e.update();
cg.import_local_transformation(cl_self,lfs_cache_e);
}
// print result
if(verbose){
std::cout << "constraints:" << std::endl;
}
std::cout << cg.size() << " constrained degrees of freedom" << std::endl;
// print result
if(verbose){
std::cout << "constraints:" << std::endl;
for (const auto& col : cg)
{
std::cout << col.first << ": "; // col index
for (const auto& row : col.second)
std::cout << "(" << row.first << "," << row.second << ") "; // row index , value
std::cout << std::endl;
}
std::cout << cg.size() << " constrained degrees of freedom" << std::endl;
for (const auto& col : cg)
{
std::cout << col.first << ": "; // col index
for (const auto& row : col.second)
std::cout << "(" << row.first << "," << row.second << ") "; // row index , value
std::cout << std::endl;
}
}
}; // end ConstraintsAssemblerHelper
template<typename P, typename GFS, typename CG>
void constraints(const P& p, const GFS& gfs, CG& cg,
const bool verbose = false)
{
// clear global constraints
cg.clear();
PPS::ConstraintsAssemblerHelper<P, GFS, CG, false>::assemble(p,gfs,cg,verbose);
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment