Moved to new module 'dune-solvers'

[[Imported from SVN: r2373]]

dune-solvers/boxconstraint.hh

+#ifndef DUNE_BOX_CONSTRAINT_HH
+#define DUNE_BOX_CONSTRAINT_HH
+
+#include <limits>
+
+#include <dune/common/fixedarray.hh>
+#include <dune/common/fvector.hh>
+
+template <class T, int dim>
+class BoxConstraint {
+public:
+
+    /** \brief Default constructor, makes the obstacle as large as possible */
+    BoxConstraint() {
+        for (int i=0; i<dim; i++) {
+            val[2*i]   = -std::numeric_limits<T>::max();
+            val[2*i+1] =  std::numeric_limits<T>::max();
+        }
+    }
+
+    /** \brief Set obstacle values to +/- std::numeric_limits<T>::max() */
+    void clear() {
+        for (int i=0; i<dim; i++) {
+            val[2*i]   = -std::numeric_limits<T>::max();
+            val[2*i+1] =  std::numeric_limits<T>::max();
+        }
+    }
+
+    //! Subtract vector from box
+    BoxConstraint<T,dim>& operator-= (const Dune::FieldVector<T, dim>& v) 
+    {
+        for (int i=0; i<dim; i++) {
+            val[2*i]   -= v[i];
+            val[2*i+1] -= v[i];
+        }
+        return *this;
+    }
+
+    //! Access to the lower obstacles
+    T& lower(size_t i) {return val[2*i];}
+
+    //! Const access to the lower obstacles
+    const T& lower(size_t i) const {return val[2*i];}
+
+    //! Access to the upper obstacles
+    T& upper(size_t i) {return val[2*i+1];}
+
+    //! Const access to the upper obstacles
+    const T& upper(size_t i) const {return val[2*i+1];}
+
+    //! Send box constraint to an output stream
+    friend
+    std::ostream& operator<< (std::ostream& s, const BoxConstraint<T,dim>& v)
+    {
+        for (int i=0; i<dim; i++)
+            s << "Direction: " << i <<",  val[0] = " << v.val[2*i] 
+              << ", val[1] = " << v.val[2*i+1] << std::endl;
+        
+        return s;
+    }
+
+protected:
+
+    Dune::array<T, 2*dim> val;
+
+};
+
+#endif

dune-solvers/transferoperators/mandelobsrestrictor.cc

+#include <dune-solvers/transferoperators/compressedmultigridtransfer.hh>
+
+template <class DiscFuncType>
+void MandelObstacleRestrictor<DiscFuncType>::
+restrict(const std::vector<BoxConstraint<typename DiscFuncType::field_type,blocksize> >&f, 
+         std::vector<BoxConstraint<typename DiscFuncType::field_type,blocksize> > &t,
+         const Dune::BitSetVector<1>& fHasObstacle,
+         const Dune::BitSetVector<1>& tHasObstacle,
+         const MultigridTransfer<DiscFuncType>& transfer,
+         const Dune::BitSetVector<blocksize>& critical)
+{
+    if (!(dynamic_cast<const CompressedMultigridTransfer<DiscFuncType>*>(&transfer)))
+    	DUNE_THROW(Dune::NotImplemented, "Implementation assumes transfer to be of type CompressedMultigridTransfer!");
+
+    const Dune::BCRSMatrix<Dune::FieldMatrix<field_type, 1, 1> >& transferMat = 
+        dynamic_cast<const CompressedMultigridTransfer<DiscFuncType>*>(&transfer)->getMatrix();
+    
+    t.resize(transferMat.M());
+    for (size_t i=0; i<t.size(); i++)
+        t[i].clear();
+    
+    // The Mandel restriction computes the maximum $m$ of all nodal base coefficients
+    // in the support of a given coarse grid base function \f$ \phi_c \f$.  If the 
+    // coefficient of \f$ \phi_c \f$ is smaller than \f$ m \f$ , it is set to \f$ m \f$.
+    
+    // We use the multigrid transfer matrix to find all relevant fine grid dofs
+    typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> >::row_type RowType;
+    typedef typename RowType::ConstIterator ColumnIterator;
+    
+    // Loop over all coarse grid dofs
+    for (int i=0; i<transferMat.N(); i++) {
+        
+        const RowType& row = transferMat[i];
+        
+        ColumnIterator cIt    = row.begin();
+        ColumnIterator cEndIt = row.end();
+        
+        // Loop over all fine grid base functions in the support of tObsIt
+        for(; cIt!=cEndIt; ++cIt) {
+            
+            // Sort out spurious entries due to numerical dirt
+            if ((*cIt)[0][0] < 1e-3)
+                continue;
+            
+            for (int j=0; j<blocksize; j++)
+                if (!critical[i][j]) {
+                    t[cIt.index()].lower(j) = std::max(t[cIt.index()].lower(j), f[i].lower(j));
+                    t[cIt.index()].upper(j) = std::min(t[cIt.index()].upper(j), f[i].upper(j));
+                }
+
+        }
+        
+    }
+
+}

dune-solvers/transferoperators/mandelobsrestrictor.hh

+#ifndef MANDEL_OBSTACLE_RESTRICTOR_HH
+#define MANDEL_OBSTACLE_RESTRICTOR_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/ag-common/obstaclerestrictor.hh>
+#include <dune-solvers/transferoperators/multigridtransfer.hh>
+#include "boxconstraint.hh"
+
+template <class DiscFuncType>
+class MandelObstacleRestrictor : public ObstacleRestrictor<DiscFuncType>
+{
+    typedef typename DiscFuncType::field_type field_type;
+    
+    enum {blocksize = DiscFuncType::block_type::dimension};
+
+public:
+    
+    virtual void restrict(const std::vector<BoxConstraint<field_type,blocksize> >& f, 
+                          std::vector<BoxConstraint<field_type,blocksize> >& t,
+                          const Dune::BitSetVector<1>& fHasObstacle,
+                          const Dune::BitSetVector<1>& tHasObstacle,
+                          const MultigridTransfer<DiscFuncType>& transfer,
+                          const Dune::BitSetVector<blocksize>& critical);
+};
+
+// include implementation
+#include "mandelobsrestrictor.cc"
+
+#endif

dune-solvers/transferoperators/obstaclerestrictor.hh

+#ifndef DUNE_OBSTACLE_RESTRICTOR_HH
+#define DUNE_OBSTACLE_RESTRICTOR_HH
+
+#include <vector>
+
+#include <dune-solvers/transferoperators/multigridtransfer.hh>
+#include <dune/ag-common/boxconstraint.hh>
+
+/** 
+ * \todo Do we need the template parameter?
+ */
+template <class DiscFuncType>
+class ObstacleRestrictor
+{
+
+    enum {blocksize = DiscFuncType::block_type::dimension};
+
+    typedef typename DiscFuncType::field_type field_type;
+
+public:
+    
+    virtual void restrict(const std::vector<BoxConstraint<field_type,blocksize> >& f, 
+                          std::vector<BoxConstraint<field_type,blocksize> >& t,
+                          const Dune::BitSetVector<1>& fHasObstacle,
+                          const Dune::BitSetVector<1>& tHasObstacle,
+                          const MultigridTransfer<DiscFuncType>& transfer,
+                          const Dune::BitSetVector<blocksize>& critical) = 0;
+    
+    
+};
+
+#endif