pointconstraintsmaker.h

Go to the documentation of this file.

#ifndef POINTCONSTRAINTSMAKER_H_
#define POINTCONSTRAINTSMAKER_H_

#include "userdefineddofconstraints.h"

namespace DOpE
{
        template<template<int, int> class DH, int dopedim, int dealdim = dopedim>
        class PointConstraints : public UserDefinedDoFConstraints<DH,dopedim,dealdim>
        {
        public:
          PointConstraints(const std::vector<dealii::Point<dealdim> > &c_points,
                           const std::vector<std::vector<bool> > &c_comps)
        : UserDefinedDoFConstraints<DH,dopedim,dealdim>(), c_points_(c_points), c_comps_(c_comps)
        {
          if(c_points_.size() != c_comps_.size())
            throw DOpEException("Number of Entries not matching!","PointConstraints::PointConstraints");
        }
          
          virtual void MakeStateDoFConstraints(
            const DOpEWrapper::DoFHandler<dealdim, DH> & dof_handler,
            dealii::ConstraintMatrix& constraint_matrix) const;
        
          virtual void
            MakeControlDoFConstraints(
              const DOpEWrapper::DoFHandler<dopedim, DH> & /*dof_handler*/,
              dealii::ConstraintMatrix& /*dof_constraints*/) const {}

        private:
          
          const std::vector<Point<dealdim> > & c_points_;
          const std::vector<std::vector<bool> > & c_comps_;
        };

        template<template<int, int> class DH, int dopedim, int dealdim>
          void PointConstraints<DH, dopedim, dealdim>::MakeStateDoFConstraints(
            const DOpEWrapper::DoFHandler<dealdim, DH > & dof_handler,
            dealii::ConstraintMatrix& constraint_matrix) const
        {
          std::vector<dealii::Point<dealdim> > support_points(dof_handler.n_dofs());
          STHInternals::MapDoFsToSupportPoints(this->GetMapping(),dof_handler, support_points);
          for(unsigned int i = 0; i < c_points_.size(); i++)
          {       
            std::vector<bool> selected_dofs(dof_handler.n_dofs());
#if DEAL_II_MAJOR_VERSION >= 8
            //Newer dealii Versions have changed the interface
            dealii::ComponentMask components(c_comps_[i]);
            DoFTools::extract_dofs(dof_handler,components,selected_dofs);
#else //Less than 8.0
#if DEAL_II_MAJOR_VERSION >= 7
#if DEAL_II_MINOR_VERSION >= 3
            //Newer dealii Versions have changed the interface
            dealii::ComponentMask components(c_comps_[i]);
            DoFTools::extract_dofs(dof_handler,components,selected_dofs);
#else //Less than 7.3 
            DoFTools::extract_dofs(dof_handler,c_comps_[i],selected_dofs);
#endif
#else //Less than 7.0
            DoFTools::extract_dofs(dof_handler,c_comps_[i],selected_dofs);
#endif
#endif

            bool found = false;
            for(unsigned int p = 0; p < support_points.size(); p++)
            {
              if(c_points_[i].distance(support_points[p]) < 1.e-12)
              {
                found = true;
                if(selected_dofs[p] == true)
                {
                  constraint_matrix.add_line(p);
                }
                //Note that you cannot stop the loop here, since more than 
                //one dof may be located at a given point!
              }
            }
            if(!found)
            {
              throw DOpEException("Points not found!","PointConstraints::MakeStateDoFConstraints");
            }
          }
        }

}

#endif /* POINTCONSTRAINTSMAKER_H_ */