FiniteElementProblem.hh 7.09 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
//=============================================================================
//
//  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;
48
49

  virtual double max_feasible_step( const double* _x, const double* _v) { return DBL_MAX;}
David Bommes's avatar
David Bommes committed
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
};


// 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)];

185
      triplets_.clear();
David Bommes's avatar
David Bommes committed
186
187
188
189
190
      element_.eval_hessian(x_, instances_.c(i), triplets_);

      for(unsigned int j=0; j<triplets_.size(); ++j)
      {
        // add re-indexed Triplet
Max Lyon's avatar
Max Lyon committed
191
192
        _triplets.push_back(Triplet( (int)instances_.index(i,triplets_[j].row()),
                                     (int)instances_.index(i,triplets_[j].col()),
David Bommes's avatar
David Bommes committed
193
194
195
196
197
                                     triplets_[j].value()                   ));
      }
    }
  }

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
  virtual double max_feasible_step( const double* _x, const double* _v)
  {
    double t = DBL_MAX;

    for(unsigned int i=0; i<instances_.size(); ++i)
    {
      // get local x vector and v vector
      for(unsigned int j=0; j<NV; ++j)
      {
        x_[j] = _x[instances_.index(i,j)];
        g_[j] = _v[instances_.index(i,j)];
      }
      t = std::min(t, element_.max_feasible_step(x_, g_, instances_.c(i)));
    }

    return t;
  }


David Bommes's avatar
David Bommes committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
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
239
  FiniteElementProblem(const unsigned int _n);
David Bommes's avatar
David Bommes committed
240

241
  void add_set(FiniteElementSetBase* _fe_set);
David Bommes's avatar
David Bommes committed
242

243
  void clear_sets();
David Bommes's avatar
David Bommes committed
244

245
  std::vector<double>& x();
David Bommes's avatar
David Bommes committed
246
247

  // problem definition
248
  virtual int    n_unknowns   (                                );
David Bommes's avatar
David Bommes committed
249

250
  virtual void   initial_x    (       double* _x               );
David Bommes's avatar
David Bommes committed
251

252
  virtual double eval_f       ( const double* _x );
David Bommes's avatar
David Bommes committed
253

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

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

258
259
260
  // return largest value t such that _x+t*_v is feasible
  virtual double max_feasible_step( const double* _x, const double* _v);

David Bommes's avatar
David Bommes committed
261

262
  virtual void   store_result ( const double* _x );
David Bommes's avatar
David Bommes committed
263
264

  // advanced properties (ToDo better handling)
265
266
  virtual bool   constant_gradient() const;
  virtual bool   constant_hessian()  const;
David Bommes's avatar
David Bommes committed
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289

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
//=============================================================================