FiniteElementProblem.hh 6.45 KB
Newer Older
David Bommes's avatar
David Bommes committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
//=============================================================================
//
//  CLASS FiniteElementProblem
//
//=============================================================================


#ifndef COMISO_FINITEELEMENTPROBLEM_HH
#define COMISO_FINITEELEMENTPROBLEM_HH


//== INCLUDES =================================================================

#include <CoMISo/Config/CoMISoDefines.hh>
#include "NProblemInterface.hh"


//== FORWARDDECLARATIONS ======================================================

//== NAMESPACES ===============================================================

namespace COMISO { 

//== CLASS DEFINITION =========================================================

	      
/** \class FiniteElementProblem

    Brief Description.
  
    A more elaborate description follows.
*/


class FiniteElementSetBase
{
public:

  typedef Eigen::Triplet<double> Triplet;

  virtual ~FiniteElementSetBase() {}

  virtual double eval_f( const double* _x ) = 0;

  virtual void accumulate_gradient( const double* _x , double* _g) = 0;

  virtual void accumulate_hessian ( const double* _x , std::vector<Triplet>& _triplets) = 0;
};


// Concepts

//class FiniteElementType
//{
//  // define dimensions of variables and constants
//  const static int NV = 3;
//  const static int NC = 1;
//
//  typedef Eigen::Matrix<size_t,NV,1> VecI;
//  typedef Eigen::Matrix<double,NV,1> VecV;
//  typedef Eigen::Matrix<double,NC,1> VecC;
//  typedef Eigen::Triplet<double> Triplet;
//
//  inline double eval_f       (const VecV& _x, const VecC& _c) const;
//  inline void   eval_gradient(const VecV& _x, const VecC& _c, VecV& _g) const;
//  inline void   eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets) const;
//};


//class FiniteElementInstancesType
//{
//
//  typedef Eigen::Matrix<double,NC,1> VecC;
//
//  size_t size() const;
//
//  size_t index( const size_t _instance, const size_t _nr) const;
//
//  const VecC& c( const size_t _instance) const;
//
//};

template<class FiniteElementType>
class FiniteElementInstancesVector
{
public:
  const static int NV = FiniteElementType::NV;
  const static int NC = FiniteElementType::NC;

  typedef typename FiniteElementType::VecI VecI;
  typedef typename FiniteElementType::VecC VecC;

  void add_element(const VecI& _indices, const VecC& _constants)
  {
    indices_.push_back(_indices);
    constants_.push_back( _constants);
  }

  void clear_elements()
  {
    indices_.clear();
    constants_.clear();
  }

  size_t size() const
  {
    return indices_.size();
  }

  size_t index( const size_t _instance, const size_t _nr) const
  {
    return indices_[_instance][_nr];
  }

  const VecC& c( const size_t _instance) const
  {
    return constants_[_instance];
  }

private:
  std::vector<VecI> indices_;
  std::vector<VecC> constants_;
};


template<class FiniteElementType, class FiniteElementInstancesType = FiniteElementInstancesVector<FiniteElementType> >
class FiniteElementSet : public FiniteElementSetBase
{
public:

  // export dimensions of element
  const static int NV = FiniteElementType::NV;
  const static int NC = FiniteElementType::NC;

  typedef typename FiniteElementType::VecI VecI;
  typedef typename FiniteElementType::VecV VecV;
  typedef typename FiniteElementType::VecC VecC;

  // access element for setting constants per element etc.
  FiniteElementType& element() { return element_;}

  // access instances for adding etc.
  FiniteElementInstancesType& instances() { return instances_;}

  virtual double eval_f( const double* _x )
  {
    double f(0.0);
    for(unsigned int i=0; i<instances_.size(); ++i)
    {
      // get local x vector
      for(unsigned int j=0; j<NV; ++j)
        x_[j] = _x[instances_.index(i,j)];

      f += element_.eval_f(x_, instances_.c(i));
    }
    return f;
  }

  virtual void accumulate_gradient( const double* _x , double* _g)
  {
    for(unsigned int i=0; i<instances_.size(); ++i)
    {
      // get local x vector
      for(unsigned int j=0; j<NV; ++j)
        x_[j] = _x[instances_.index(i,j)];

      element_.eval_gradient(x_, instances_.c(i), g_);

      // accumulate into global gradient
      for(unsigned int j=0; j<NV; ++j)
        _g[instances_.index(i,j)] += g_[j];
    }
  }

  virtual void accumulate_hessian ( const double* _x , std::vector<Triplet>& _triplets)
  {
    for(unsigned int i=0; i<instances_.size(); ++i)
    {
      // get local x vector
      for(unsigned int j=0; j<NV; ++j)
        x_[j] = _x[instances_.index(i,j)];

      element_.eval_hessian(x_, instances_.c(i), triplets_);

      for(unsigned int j=0; j<triplets_.size(); ++j)
      {
        // add re-indexed Triplet
        _triplets.push_back(Triplet( instances_.index(i,triplets_[j].row()),
                                     instances_.index(i,triplets_[j].col()),
                                     triplets_[j].value()                   ));
      }
    }
  }

private:

  FiniteElementType element_;

  FiniteElementInstancesType instances_;

  VecV x_;
  VecV g_;

  std::vector<Triplet> triplets_;
};


class COMISODLLEXPORT FiniteElementProblem : public NProblemInterface
{
public:

//  typedef Eigen::SparseMatrix<double,Eigen::ColMajor> SMatrixNP;

  typedef FiniteElementSetBase::Triplet Triplet;

  /// Default constructor
217
  FiniteElementProblem(const unsigned int _n);
David Bommes's avatar
David Bommes committed
218 219

  /// Destructor
220
  virtual ~FiniteElementProblem();
David Bommes's avatar
David Bommes committed
221

222
  void add_set(FiniteElementSetBase* _fe_set);
David Bommes's avatar
David Bommes committed
223

224
  void clear_sets();
David Bommes's avatar
David Bommes committed
225

226
  std::vector<double>& x();
David Bommes's avatar
David Bommes committed
227 228

  // problem definition
229
  virtual int    n_unknowns   (                                );
David Bommes's avatar
David Bommes committed
230

231
  virtual void   initial_x    (       double* _x               );
David Bommes's avatar
David Bommes committed
232

233
  virtual double eval_f       ( const double* _x );
David Bommes's avatar
David Bommes committed
234

235
  virtual void   eval_gradient( const double* _x, double*    _g);
David Bommes's avatar
David Bommes committed
236

237
  virtual void   eval_hessian ( const double* _x, SMatrixNP& _H);
David Bommes's avatar
David Bommes committed
238 239


240
  virtual void   store_result ( const double* _x );
David Bommes's avatar
David Bommes committed
241 242

  // advanced properties (ToDo better handling)
243 244
  virtual bool   constant_gradient() const;
  virtual bool   constant_hessian()  const;
David Bommes's avatar
David Bommes committed
245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267

private:

  // number of unknowns
  unsigned int n_;

  // current/initial solution
  std::vector<double> x_;

  // finite element sets
  std::vector<FiniteElementSetBase*> fe_sets_;

  // vector of triplets (avoid re-allocation)
  std::vector<Triplet> triplets_;
};


//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_FINITEELEMENTPROBLEM_HH defined
//=============================================================================