G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsFuncData.h
Go to the documentation of this file.
1 
15 #pragma once
16 
18 #include<gsCore/gsBoundary.h>
19 
20 namespace gismo
21 {
22 
23 template <typename T> class gsFunctionSet;
24 
42 namespace util {
43 
44 // Adaptor to compute Hessian
45 template <typename Derived>
46 gsMatrix<typename Derived::Scalar> secDerToHessian(const gsEigen::DenseBase<Derived> & secDers,
47  const index_t dim)
48 {
49  index_t sz = dim*(dim+1)/2;
50  auto ders = secDers.reshaped(sz, secDers.size() / sz );
51  gsMatrix<typename Derived::Scalar> hessian(dim*dim, ders.cols() );
52 
53  switch ( dim )
54  {
55  case 1:
56  hessian = secDers.transpose();
57  break;
58  case 2:
59  hessian.row(0)=ders.row(0);//0,0
60  hessian.row(1)=//1,0
61  hessian.row(2)=ders.row(2);//0,1
62  hessian.row(3)=ders.row(1);//1,1
63  break;
64  case 3:
65  hessian.row(0)=ders.row(0);//0,0
66  hessian.row(3)=//0,1
67  hessian.row(1)=ders.row(3);//1,0
68  hessian.row(6)=//0,2
69  hessian.row(2)=ders.row(4);//2,0
70  hessian.row(4)=ders.row(1);//1,1
71  hessian.row(7)=//1,2
72  hessian.row(5)=ders.row(5);//2,1
73  hessian.row(8)=ders.row(2);//2,2
74  break;
75  default:
76  sz = 0;
77  for (index_t k=0; k!=dim; ++k ) // for all rows
78  {
79  hessian.row((dim+1)*k) = ders.row(k);
80  for (index_t l=k+1; l<dim; ++l ) // for all cols
81  hessian.row(dim*k+l) =
82  hessian.row(dim*l+k) = ders.row(dim + sz++);
83  }
84  break;
85  }
86  return hessian;
87 }
88 
89 template <typename Derived>
90 void hessianToSecDer (const gsEigen::DenseBase<Derived> & hessian,
91  const index_t dim,
93 {
94  GISMO_ASSERT(hessian.cols() == dim, "single Hessian implemented");
95  secDers.resize(dim*(dim+1)/2, hessian.cols() / dim );
96  switch ( dim )
97  {
98  case 1:
99  secDers=hessian.transpose();
100  break;
101  case 2:
102  secDers.at(0)=hessian(0,0);
103  secDers.at(1)=hessian(1,1);
104  secDers.at(2)=hessian(1,0);//==hessian(0,1));
105  break;
106  case 3:
107  secDers.at(0)=hessian(0,0);
108  secDers.at(1)=hessian(1,1);
109  secDers.at(2)=hessian(2,2);
110  secDers.at(3)=hessian(0,1);//==hessian(1,0));
111  secDers.at(4)=hessian(0,2);//==hessian(2,0));
112  secDers.at(5)=hessian(1,2);//==hessian(2,1));
113  break;
114  default:
115  GISMO_ERROR("NO_IMPLEMENTATION");
116  break;
117  }
118 }
119 
120 }//namespace util
121 
122 template <typename T>
123 class gsFuncData
124 {
125  friend class gsFunctionSet<T>;
126 
127 public:
128  typedef typename gsMatrix<T>::constColumn constColumn;
129  // Types for returning quick access to data in matrix format
130  typedef gsAsConstMatrix<T, -1, -1> matrixView;
131  typedef gsEigen::Transpose<typename matrixView::Base> matrixTransposeView;
132 
133 public:
134  mutable unsigned flags;
135  index_t patchId; // move to mapdata
136 
137  gsMatrix<index_t> actives;
138 
142  std::vector<gsMatrix<T> > values;
143 
144  gsMatrix<T> curls;
145  gsMatrix<T> divs;
146  gsMatrix<T> laplacians;
147 
151  std::pair<short_t, short_t> dim;
152 
153 public:
159  explicit gsFuncData(unsigned flags = 0, int patch = 0)
160  : flags(flags), patchId(patch)
161  { }
162 
163 public:
164 
170  void addFlags (unsigned newFlags)
171  { flags = flags|newFlags; }
172 
173 
174  int maxDeriv() const
175  {
176  if (flags & (NEED_LAPLACIAN|NEED_DERIV2|NEED_HESSIAN) )
177  return 2;
178  else if (flags & (NEED_DERIV|NEED_GRAD_TRANSFORM|NEED_CURL|NEED_DIV) )
179  return 1;
180  else if (flags & (NEED_VALUE) )
181  return 0;
182  return -1;
183  }
184 
189  unsigned bytesUsed() const
190  { /*to do*/
191  return 1;
192  }
193 
195  void clear()
196  {
197  flags = 0;
198  patchId = -1;
199  actives.clear();
200  values.clear();
201  curls.clear();
202  divs.clear();
203  laplacians.clear();
204  //dim;
205  }
206 
208  void swap(gsFuncData & other)
209  {
210  std::swap(flags , other.flags );
211  std::swap(patchId, other.patchId);
212  std::swap(dim, other.dim);
213  actives .swap(other.actives );
214  values .swap(other.values );
215  curls .swap(other.curls );
216  divs .swap(other.divs );
217  laplacians.swap(other.laplacians);
218  }
219 
220 public:
221 
222  inline const gsMatrix<index_t> & allActives() const
223  {
224  GISMO_ASSERT(flags & NEED_ACTIVE,
225  "actives are not computed unless the NEED_ACTIVE flag is set.");
226  GISMO_ASSERT(0!=actives.size(), "actives were not computed.");
227  return actives;
228  }
229 
230  inline const gsMatrix<T> & allValues() const
231  {
232  GISMO_ASSERT(flags & NEED_ACTIVE,
233  "values are not computed unless the NEED_ACTIVE flag is set.");
234  GISMO_ASSERT(0!=values.size(), "values were not computed.");
235  return values.front();
236  }
237 
238  inline gsMatrix<index_t>::constColumn active(index_t point = 0) const
239  {
240  GISMO_ASSERT(flags & NEED_ACTIVE,
241  "actives are not computed unless the NEED_ACTIVE flag is set.");
242  return actives.col(point);
243  }
244 
245  inline matrixView eval (index_t point) const
246  {
247  GISMO_ASSERT(flags & NEED_VALUE,
248  "values are not computed unless the NEED_VALUE flag is set.");
249  return values[0].reshapeCol(point, dim.second, values[0].rows()/dim.second);
250  }
251 
252  inline matrixView deriv (index_t point) const
253  {
254  GISMO_ASSERT(flags & NEED_DERIV,
255  "derivs are not computed unless the NEED_DERIV flag is set.");
256  return values[1].reshapeCol(point, derivSize(), values[1].rows()/derivSize());
257  }
258 
259  inline matrixView deriv2 (index_t point) const
260  {
261  GISMO_ASSERT(flags & NEED_DERIV2,
262  "deriv2s are not computed unless the NEED_DERIV2 flag is set.");
263  return values[2].reshapeCol(point, deriv2Size(), values[2].rows()/deriv2Size());
264  }
265 
266  inline matrixView curl (index_t point) const
267  {
268  GISMO_ASSERT(flags & NEED_CURL,
269  "curls are not computed unless the NEED_CURL flag is set.");
270  return curls.reshapeCol(point, dim.second, curls.rows()/dim.second );
271  }
272 
273  inline matrixView div (index_t point) const
274  {
275  GISMO_ASSERT(flags & NEED_DIV,
276  "divs are not computed unless the NEED_DIV flag is set.");
277  return divs.reshapeCol(point, divSize(), divs.rows()/divSize() );
278  }
279 
280  inline matrixView laplacian (index_t point) const
281  {
283  "laplacians are not computed unless the NEED_LAPLACIAN flag is set.");
284  return laplacians.reshapeCol(point, dim.second, laplacians.rows()/dim.second );
285  }
286 
287 
288  inline matrixTransposeView jacobian(index_t point, index_t func = 0) const
289  {
290  GISMO_ASSERT(flags & (NEED_DERIV),
291  "jacobian access needs the computation of derivs: set the NEED_JACOBIAN flag.");
292  return gsAsConstMatrix<T, Dynamic, Dynamic>(&values[1].coeffRef(func*derivSize(),point),dim.first,dim.second).transpose();
293  }
294 
295  inline gsMatrix<T> hessian(index_t point, index_t func = 0) const
296  {
297  GISMO_ASSERT(flags & NEED_HESSIAN,
298  "hessian access needs the computation of 2nd derivs: set the NEED_HESSIAN flag.");
299  gsMatrix<T> res(dim.first,dim.first);
300  const index_t dsz = dim.first*(dim.first+1) / 2;
301  res = util::secDerToHessian(values[2].block(func*dsz,point,dsz,1), dim.first);
302  res.resize(dim.first,dim.first);
303  return res;
304  }
305 
306 //protected:
307 
309  int derivSize () const {return dim.first*dim.second;}
311  int deriv2Size() const {return dim.second*dim.first*(dim.first+1) / 2; }
313  int divSize () const {return dim.second/dim.first;}
314 };
315 
316 
323 template <typename T>
324 class gsMapData : public gsFuncData<T>
325 {
326 public:
327  typedef gsFuncData<T> Base;
328  typedef typename Base::constColumn constColumn;
329  typedef typename Base::matrixView matrixView;
330  typedef typename Base::matrixTransposeView matrixTransposeView;
331 
332 public:
337  explicit gsMapData(unsigned flags = 0)
338  : Base(flags), side(boundary::none)
339  { }
340 
341 public:
342  using Base::flags;
343  using Base::values;
344  using Base::dim;
345 
346  boxSide side;
347 
349 
350  gsMatrix<T> measures;
353  gsMatrix<T> normals;
354  gsMatrix<T> outNormals; // only for the boundary
355 
356 public:
357  inline constColumn point(const index_t point) const { return points.col(point);}
358 
359  inline T measure(const index_t point) const
360  {
361  GISMO_ASSERT(flags & NEED_MEASURE,
362  "measures are not computed unless the NEED_MEASURE flag is set.");
363  return measures(0,point);
364  }
365 
366  inline matrixView fundForm(const index_t point) const
367  {
369  "fundForms are not computed unless the NEED_2ND_FFORM flag is set.");
370  return fundForms.reshapeCol(point, dim.first, dim.first);
371  }
372 
373  inline constColumn normal(const index_t point) const
374  {
375  GISMO_ASSERT(flags & NEED_NORMAL,
376  "normals are not computed unless the NEED_NORMAL flag is set.");
377  return normals.col(point);
378  }
379 
380  inline constColumn outNormal(const index_t point) const
381  {
383  "normals are not computed unless the NEED_NORMAL flag is set.");
384  return outNormals.col(point);
385  }
386 
387  inline matrixTransposeView jacobians() const
388  {
389  GISMO_ASSERT(flags & NEED_DERIV,
390  "jacobian access needs the computation of derivs: set the NEED_DERIV flag." << this->maxDeriv() );
391  return gsAsConstMatrix<T, Dynamic, Dynamic>(&values[1].coeffRef(0,0), dim.first,dim.second*values[1].cols()).transpose();
392  }
393 };
394 
395 } // namespace gismo
396 
397 namespace std
398 {
399 
400 // Specializations of std::swap
401 template <typename T>
402 void swap(gismo::gsFuncData<T> & f1, gismo::gsFuncData<T> & f2)
403 {
404  f1.swap(f2);
405 }
406 
407 } // namespace std
Div operator.
Definition: gsForwardDeclarations.h:78
gsMatrix< T > fundForms
Second fundumental forms.
Definition: gsFuncData.h:351
gsMapData(unsigned flags=0)
Main constructor.
Definition: gsFuncData.h:337
Gradient of the object.
Definition: gsForwardDeclarations.h:73
Provides structs and classes related to interfaces and boundaries.
Struct that defines the boundary sides and corners and types of a geometric object.
Definition: gsBoundary.h:55
Outward normal on the boundary.
Definition: gsForwardDeclarations.h:86
the gsMapData is a cache of pre-computed function (map) values.
Definition: gsFuncData.h:324
Second derivatives.
Definition: gsForwardDeclarations.h:80
#define index_t
Definition: gsConfig.h:32
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition: gsMatrix.h:38
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gsMatrix< T > jacInvTr
Inverse of the Jacobian matrix (transposed)
Definition: gsFuncData.h:352
Laplacian.
Definition: gsForwardDeclarations.h:83
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition: gsMatrix.h:211
Normal vector of the object.
Definition: gsForwardDeclarations.h:85
Hessian matrix.
Definition: gsForwardDeclarations.h:82
Active function ids.
Definition: gsForwardDeclarations.h:84
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition: gsFuncData.h:23
gsAsMatrix< T, Dynamic, Dynamic > reshapeCol(index_t c, index_t n, index_t m)
Returns column c of the matrix resized to n x m matrix This function assumes that the matrix is size ...
Definition: gsMatrix.h:231
Curl operator.
Definition: gsForwardDeclarations.h:79
Struct which represents a certain side of a box.
Definition: gsBoundary.h:84
Creates a mapped object or data pointer to a const matrix without copying data.
Definition: gsLinearAlgebra.h:127
Second fundamental form.
Definition: gsForwardDeclarations.h:87
The density of the measure pull back.
Definition: gsForwardDeclarations.h:76
EIGEN_STRONG_INLINE curl_expr< T > curl(const gsFeVariable< T > &u)
The curl of a finite element variable.
Definition: gsExpressions.h:4500
Value of the object.
Definition: gsForwardDeclarations.h:72
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
This is the main header file that collects wrappers of Eigen for linear algebra.
gsMatrix< T > points
input (parametric) points
Definition: gsFuncData.h:348
Gradient transformation matrix.
Definition: gsForwardDeclarations.h:77