FiniteElementProblem.hh 7.12 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
185
186
187
188
189
};


// 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
Max Lyon's avatar
Max Lyon committed
190
191
        _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
192
193
194
195
196
                                     triplets_[j].value()                   ));
      }
    }
  }

197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
  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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
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
238
  FiniteElementProblem(const unsigned int _n);
David Bommes's avatar
David Bommes committed
239
240

  /// Destructor
241
  virtual ~FiniteElementProblem();
David Bommes's avatar
David Bommes committed
242

243
  void add_set(FiniteElementSetBase* _fe_set);
David Bommes's avatar
David Bommes committed
244

245
  void clear_sets();
David Bommes's avatar
David Bommes committed
246

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

  // problem definition
250
  virtual int    n_unknowns   (                                );
David Bommes's avatar
David Bommes committed
251

252
  virtual void   initial_x    (       double* _x               );
David Bommes's avatar
David Bommes committed
253

254
  virtual double eval_f       ( const double* _x );
David Bommes's avatar
David Bommes committed
255

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

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

260
261
262
  // 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
263

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

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

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