DUNE-ACFEM (2.5.1)

functionoperations.hh
1 #ifndef __FUNCTIONEXPRESSION_HH__
2 #define __FUNCTIONEXPRESSION_HH__
3 
4 #include "../expressions/expressionoperations.hh"
5 
6 namespace Dune
7 {
8 
9  namespace ACFem
10  {
11 
20  template<class Operation>
21  struct FunctionExpressionOperation
22  {};
23 
24  template<class Operation>
25  struct UnaryFunctionExpressionOperation
26  : public FunctionExpressionOperation<Operation>
27  {};
28 
29  template<class Operation>
30  struct BinaryFunctionExpressionOperation
31  : public FunctionExpressionOperation<Operation>
32  {};
33 
34  template<>
35  struct BinaryFunctionExpressionOperation<PlusOperation>
36  : public FunctionExpressionOperation<PlusOperation>
37  {
38  template<class LeftFunctionType, class RightFunctionType,
39  class PointType, class RangeType>
40  static void evaluate (const LeftFunctionType& f, const RightFunctionType& g,
41  const PointType &x, RangeType &ret)
42  {
43  RangeType tmp;
44 
45  f.evaluate(x, ret);
46  g.evaluate(x, tmp);
47 
48  //std::cout << ret << " + " << tmp << std::endl;
49 
50  ret += tmp;
51  }
52 
53  template<class LeftFunctionType,
54  class RightFunctionType,
55  class PointType, class JacobianRangeType,
56  class LeftPWConst = std::false_type,
57  class RightPWConst = std::false_type>
58  static void jacobian (const LeftFunctionType& f, const RightFunctionType& g,
59  const PointType &x, JacobianRangeType &ret,
60  const LeftPWConst& fConst = LeftPWConst(),
61  const RightPWConst& gConst = RightPWConst())
62  {
63  if (fConst) {
64  g.jacobian(x, ret);
65  return;
66  }
67 
68  if (gConst) {
69  f.jacobian(x, ret);
70  return;
71  }
72 
73  JacobianRangeType tmp;
74 
75  f.jacobian(x, ret);
76  g.jacobian(x, tmp);
77 
78  ret += tmp;
79  }
80 
81  template<class LeftFunctionType, class RightFunctionType,
82  class PointType, class HessianRangeType,
83  class LeftPWConst = std::false_type,
84  class RightPWConst = std::false_type>
85  static void hessian (const LeftFunctionType& f, const RightFunctionType& g,
86  const PointType &x, HessianRangeType &ret,
87  const LeftPWConst& fConst = LeftPWConst(),
88  const RightPWConst& gConst = RightPWConst())
89  {
90  if (fConst) {
91  g.hessian(x, ret);
92  return;
93  }
94 
95  if (gConst) {
96  f.hessian(x, ret);
97  return;
98  }
99 
100  HessianRangeType tmp;
101 
102  f.hessian(x, ret);
103  g.hessian(x, tmp);
104 
105  //ret += tmp; Dune bug: does not work
106  for (int i = 0; i < HessianRangeType::dimension; ++i) {
107  ret[i] += tmp[i];
108  }
109  }
110 
111  };
112 
113  template<>
114  struct BinaryFunctionExpressionOperation<MinusOperation>
115  : public FunctionExpressionOperation<MinusOperation>
116  {
117  template<class LeftFunctionType, class RightFunctionType, class PointType, class RangeType>
118  static void evaluate(const LeftFunctionType& f, const RightFunctionType& g, const PointType &x, RangeType &ret)
119  {
120  RangeType tmp;
121 
122  f.evaluate(x, ret);
123  g.evaluate(x, tmp);
124 
125  //std::cout << ret << " - " << tmp << std::endl;
126 
127  ret -= tmp;
128  }
129 
130  template<class LeftFunctionType,
131  class RightFunctionType,
132  class PointType, class JacobianRangeType,
133  class LeftPWConst = std::false_type,
134  class RightPWConst = std::false_type>
135  static void jacobian(const LeftFunctionType& f, const RightFunctionType& g,
136  const PointType &x, JacobianRangeType &ret,
137  const LeftPWConst& fConst = LeftPWConst(),
138  const RightPWConst& gConst = RightPWConst())
139  {
140  if (fConst) {
141  g.jacobian(x, ret);
142  return;
143  }
144 
145  if (gConst) {
146  f.jacobian(x, ret);
147  ret *= -1.0;
148  return;
149  }
150 
151  JacobianRangeType tmp;
152 
153  f.jacobian(x, ret);
154  g.jacobian(x, tmp);
155 
156  ret -= tmp;
157  }
158 
159  template<class LeftFunctionType, class RightFunctionType,
160  class PointType, class HessianRangeType,
161  class LeftPWConst = std::false_type,
162  class RightPWConst = std::false_type>
163  static void hessian(const LeftFunctionType& f, const RightFunctionType& g,
164  const PointType &x, HessianRangeType &ret,
165  const LeftPWConst& fConst = LeftPWConst(),
166  const RightPWConst& gConst = RightPWConst())
167  {
168  if (fConst) {
169  g.hessian(x, ret);
170  return;
171  }
172 
173  if (gConst) {
174  f.hessian(x, ret);
175  return;
176  }
177 
178  HessianRangeType tmp;
179 
180  f.hessian(x, ret);
181  g.hessian(x, tmp);
182 
183  //ret -= tmp; Dune bug: does not work
184  for (int i = 0; i < HessianRangeType::dimension; ++i) {
185  ret[i] -= tmp[i];
186  }
187  }
188 
189  };
190 
192  template<>
193  struct BinaryFunctionExpressionOperation<MultiplyOperation>
194  : public FunctionExpressionOperation<MultiplyOperation>
195  {
196  // This is the scalar product case
197  template<class Result, class Matrix, class Factor, int numRows, int dimFactor>
198  struct MulEqHelper
199  {
200  //
201  static void evaluate(Result& res, const Matrix& val, const Factor& factor)
202  {
203  static_assert(Result::rows == 1,
204  "Wrong result type for jacobian of multiplication");
205 
206  // This is the case where res is a FieldMatrix with one row, so:
207  val.mtv(factor, res[0]);
208  }
209  };
210 
211  // factor is a scalar
212  template<class Result, class Matrix, class Factor, int numRows>
213  struct MulEqHelper<Result, Matrix, Factor, numRows, 1>
214  {
215  static void evaluate(Result& res, const Matrix& val, const Factor& factor)
216  {
217  res = val;
218  res *= factor;
219  }
220  };
221 
222  // Matrix has one rows, factor is a vector, Kronecker product
223  template<class Result, class Matrix, class Factor, int dimFactor>
224  struct MulEqHelper<Result, Matrix, Factor, 1, dimFactor>
225  {
226  static void evaluate(Result& res, const Matrix& val, const Factor& factor)
227  {
228  for (int i = 0; i < dimFactor; ++i) {
229  for (int j = 0; j < Matrix::cols; ++j) {
230  res[i][j] = val[0][j] * factor[i];
231  }
232  }
233  }
234  };
235 
236  // This is just the 1-1 diagonal case which HAS to be
237  // specialized explicitly in order to avoid ambuigities
238  template<class Result, class Matrix, class Factor>
239  struct MulEqHelper<Result, Matrix, Factor, 1, 1>
240  {
241  static void evaluate(Result& res, const Matrix& val, const Factor& factor)
242  {
243  res = val;
244  res *= factor;
245  }
246  };
247 
248  template<class Result, class Matrix, class Factor>
249  static void callMulEq(Result& res, const Matrix& val, const Factor& factor)
250  {
251  //std::cerr << __PRETTY_FUNCTION__ << std::endl;
252  MulEqHelper<Result, Matrix, Factor, Matrix::rows, Factor::dimension>::evaluate(res, val, factor);
253  }
254 
255  template<class Left, class Right, class Point, class Result,
256  unsigned dimF, unsigned dimG, unsigned dimResult>
257  struct MulEvalHelper;
258 
259  // Scalar product case
260  template<class Left, class Right, class Point, class Result, unsigned dimArg>
261  struct MulEvalHelper<Left, Right, Point, Result, dimArg, dimArg, 1>
262  {
263  static void evaluate(const Left &f, const Right& g, const Point& x, Result& res)
264  {
265  typename Left::RangeType fval;
266  typename Right::RangeType gval;
267 
268  f.evaluate(x, fval);
269  g.evaluate(x, gval);
270 
271  res = fval * gval; // scalar product
272  }
273  };
274 
275  // Scalar * Vector
276  template<class Left, class Right, class Point, class Result, unsigned dimResult>
277  struct MulEvalHelper<Left, Right, Point, Result, 1, dimResult, dimResult>
278  {
279  static void evaluate(const Left &f, const Right& g, const Point& x, Result& res)
280  {
281  typename Left::RangeType fval;
282 
283  f.evaluate(x, fval);
284  g.evaluate(x, res);
285 
286  res *= fval; // S-mult
287  }
288  };
289 
290  // Vector * Scalar
291  template<class Left, class Right, class Point, class Result, unsigned dimResult>
292  struct MulEvalHelper<Left, Right, Point, Result, dimResult, 1, dimResult>
293  {
294  static void evaluate(const Left &f, const Right& g, const Point& x, Result& res)
295  {
296  typename Right::RangeType gval;
297 
298  f.evaluate(x, res);
299  g.evaluate(x, gval);
300 
301  res *= gval; // S-mult
302  }
303  };
304 
305  // Scalar * Scalar
306  template<class Left, class Right, class Point, class Result>
307  struct MulEvalHelper<Left, Right, Point, Result, 1, 1, 1>
308  {
309  static void evaluate(const Left &f, const Right& g, const Point& x, Result& res)
310  {
311  typename Right::RangeType gval;
312 
313  f.evaluate(x, res);
314  g.evaluate(x, gval);
315 
316  res *= gval; // S-mult
317  }
318  };
319 
320  template<class Left, class Right, class Point, class Result>
321  static void callMulEval(const Left& f, const Right& g, const Point& x, Result& res)
322  {
323  enum {
324  leftDimRange = Left::FunctionSpaceType::dimRange,
325  rightDimRange = Right::FunctionSpaceType::dimRange,
326  };
327  typedef
328  MulEvalHelper<Left, Right, Point, Result, leftDimRange, rightDimRange, Result::dimension>
329  HelperType;
330 
331  HelperType::evaluate(f, g, x, res);
332  }
333 
334  template<class ResultType, class LeftType, class RightType, unsigned LeftDim, unsigned RightDim>
335  struct HessianHelper
336  {
337  static void evaluate(ResultType& result, const LeftType& hessian, const RightType& value)
338  {
339  // do nothing
340  }
341  };
342 
343  template<class TensorType, class ScalarType, unsigned OuterDim>
344  struct HessianHelper<TensorType, TensorType, ScalarType, OuterDim, 1>
345  {
346  static void evaluate(TensorType& result, const TensorType& hessian, const ScalarType& value)
347  {
348  result = hessian;
349  for (int k = 0; k < TensorType::dimension; ++k) {
350  result[k] *= value;
351  }
352  }
353  };
354 
355  template<class Left, class Right, class Result>
356  static void callHessianHelper(Result& result, const Left& left, const Right& right)
357  {
358  enum {
359  leftDim = Left::dimension,
360  rightDim = Right::dimension
361  };
362  HessianHelper<Result, Left, Right, leftDim, rightDim>::evaluate(result, left, right);
363  }
364 
365  template<class LeftFunctionType, class RightFunctionType, class PointType, class RangeType>
366  static void evaluate(const LeftFunctionType& f, const RightFunctionType& g, const PointType &x, RangeType &res)
367  {
368  callMulEval(f, g, x, res);
369  }
370 
371  template<class LeftFunctionType, class RightFunctionType,
372  class PointType, class JacobianRangeType,
373  class LeftPWConst = std::false_type,
374  class RightPWConst = std::false_type>
375  static void jacobian(const LeftFunctionType& f, const RightFunctionType& g,
376  const PointType &x, JacobianRangeType &ret,
377  const LeftPWConst& fConst = LeftPWConst(),
378  const RightPWConst& gConst = RightPWConst())
379  {
380  typename LeftFunctionType::RangeType fval;
381  typename RightFunctionType::RangeType gval;
382 
383  typename LeftFunctionType::JacobianRangeType fjac;
384  typename RightFunctionType::JacobianRangeType gjac;
385 
386  if (fConst) {
387  f.evaluate(x, fval);
388  g.jacobian(x, gjac);
389  callMulEq(ret, gjac, fval);
390  return;
391  }
392 
393  if (gConst) {
394  g.evaluate(x, gval);
395  f.jacobian(x, fjac);
396  callMulEq(ret, fjac, gval);
397  return;
398  }
399 
400  f.evaluate(x, fval);
401  g.evaluate(x, gval);
402 
403  f.jacobian(x, fjac);
404  g.jacobian(x, gjac);
405 
406  callMulEq(ret, fjac, gval);
407  JacobianRangeType tmp;
408  callMulEq(tmp, gjac, fval);
409  ret += tmp;
410  }
411 
412  template<class LeftFunctionType, class RightFunctionType,
413  class PointType, class HessianRangeType,
414  class LeftPWConst = std::false_type,
415  class RightPWConst = std::false_type>
416  static void hessian(const LeftFunctionType& f, const RightFunctionType& g,
417  const PointType &x, HessianRangeType &ret,
418  const LeftPWConst& fConst = LeftPWConst(),
419  const RightPWConst& gConst = RightPWConst())
420  {
421  const int dimDomain = LeftFunctionType::DomainType::dimension;
422 
423  typename LeftFunctionType::RangeType fval;
424  typename RightFunctionType::RangeType gval;
425 
426  f.evaluate(x, fval);
427  g.evaluate(x, gval);
428 
429  typename LeftFunctionType::JacobianRangeType fjac;
430  typename RightFunctionType::JacobianRangeType gjac;
431 
432  f.jacobian(x, fjac);
433  g.jacobian(x, gjac);
434 
435  typename LeftFunctionType::HessianRangeType fhes;
436  typename RightFunctionType::HessianRangeType ghes;
437 
438  f.hessian(x, fhes);
439  g.hessian(x, ghes);
440 
441  if (LeftFunctionType::RangeType::dimension == 1) {
442  callHessianHelper(ret, ghes, fval);
443  if (fConst) {
444  return;
445  }
446  for (int k = 0; k < RightFunctionType::RangeType::dimension; ++k) {
447  for (int i = 0; i < dimDomain; ++i) {
448  for (int j = 0; j < dimDomain; ++j) {
449  ret[k][i][j] += fjac[0][i]*gjac[k][j] + fjac[0][j] * gjac[k][i];
450  ret[k][i][j] += fhes[0][i][j] * gval[k];
451  }
452  }
453  }
454  } else if (RightFunctionType::RangeType::dimension == 1) {
455  callHessianHelper(ret, fhes, gval);
456  if (gConst) {
457  return;
458  }
459  for (int k = 0; k < LeftFunctionType::RangeType::dimension; ++k) {
460  for (int i = 0; i < dimDomain; ++i) {
461  for (int j = 0; j < dimDomain; ++j) {
462  ret[k][i][j] += gjac[0][i]*fjac[k][j] + gjac[0][j]*fjac[k][i];
463  ret[k][i][j] += ghes[0][i][j] * fval[k];
464  }
465  }
466  }
467  } else {
468  // Both vectors, equal size assumed. This is -- consequently
469  // -- the scalar-product case.
470  const unsigned dimDomain = LeftFunctionType::DomainType::dimension;
471  const unsigned dimRange = LeftFunctionType::RangeType::dimension;
472 
473  ret[0] = 0.;
474  for (unsigned alpha = 0; alpha < dimRange; ++alpha) {
475  if (!gConst) {
476  ghes[alpha] *= fval[alpha];
477  ret[0] += ghes[alpha];
478  }
479  if (!fConst) {
480  fhes[alpha] *= gval[alpha];
481  ret[0] += fhes[alpha];
482  }
483 
484  if (fConst || gConst) {
485  continue;
486  }
487 
488  for (unsigned j = 0; j < dimDomain; ++j) {
489  ret[0][j][j] += 2.0 * fjac[alpha][j] * gjac[alpha][j];
490  for (unsigned k = j+1; k < dimDomain; ++k) {
491  typename LeftFunctionType::RangeFieldType tmp;
492  tmp = fjac[alpha][j]*gjac[alpha][k] + fjac[alpha][k]*gjac[alpha][j];
493  ret[0][j][k] += tmp;
494  ret[0][k][j] += tmp;
495  }
496  }
497  }
498  }
499  }
500  };
501 
502  // Division. f must be __SCALAR__, we do not support matrix inversion.
503  template<>
504  struct UnaryFunctionExpressionOperation<InvertOperation>
505  : public FunctionExpressionOperation<InvertOperation>
506  {
507  template<class FunctionType, class PointType, class RangeType>
508  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
509  {
510  static_assert(FunctionType::RangeType::dimension == 1,
511  "Operand must be scalar for inversion");
512 
513  f.evaluate(x, ret);
514  ret = 1.0 / ret;
515  }
516 
517  template<class FunctionType, class PointType, class JacobianRangeType>
518  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
519  {
520  static_assert(FunctionType::RangeType::dimension == 1,
521  "Operand must be scalar for inversion");
522 
523  typename FunctionType::RangeType fval;
524 
525  f.evaluate(x, fval);
526  f.jacobian(x, ret);
527 
528  ret /= - fval * fval;
529  }
530 
531  template<class FunctionType, class PointType, class HessianRangeType>
532  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
533  {
534  static_assert(FunctionType::RangeType::dimension == 1, "Operand must be scalar for division");
535 
536  typename FunctionType::RangeType fval;
537  f.evaluate(x, fval);
538  typename FunctionType::RangeFieldType factor;
539  factor = fval * fval;
540 
541  typename FunctionType::JacobianRangeType fjac;
542  f.jacobian(x, fjac);
543 
544  f.hessian(x, ret);
545 
546  ret[0] /= - factor;
547  factor *= fval;
548 
549  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
550  ret[0][i][i] += 2.0 * fjac[0][i] * fjac[0][i] / factor;
551  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
552  typename FunctionType::RangeFieldType tmp;
553 
554  tmp = 2.0 * fjac[0][i] * fjac[0][j] / factor;
555  ret[0][i][j] += tmp;
556  ret[0][j][i] += tmp;
557  }
558  }
559  }
560 
561  };
562 
564  template<>
565  struct BinaryFunctionExpressionOperation<SMultiplyOperation>
566  : public FunctionExpressionOperation<SMultiplyOperation>
567  {
568  template<class RangeFieldType, class FunctionType, class PointType, class RangeType>
569  static void evaluate (const RangeFieldType& s, const FunctionType& f, const PointType &x, RangeType &ret)
570  {
571  f.evaluate(x, ret);
572 
573  ret *= s;
574  }
575 
576  template<class RangeFieldType, class FunctionType, class PointType, class JacobianRangeType>
577  static void jacobian (const RangeFieldType& s, const FunctionType& f, const PointType &x, JacobianRangeType &ret)
578  {
579  f.jacobian(x, ret);
580 
581  ret *= s;
582  }
583 
584  template<class RangeFieldType, class FunctionType, class PointType, class HessianRangeType>
585  static void hessian (const RangeFieldType& s, const FunctionType& f, const PointType &x, HessianRangeType &ret)
586  {
587  f.hessian(x, ret);
588 
589  ret *= s;
590  }
591  };
592 
593  // Identity operation, quite simple
594  template<>
595  struct UnaryFunctionExpressionOperation<IdentityOperation>
596  : public FunctionExpressionOperation<IdentityOperation>
597  {
598  template<class FunctionType, class PointType, class RangeType>
599  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
600  {
601  f.evaluate(x, ret);
602  }
603 
604  template<class FunctionType, class PointType, class JacobianRangeType>
605  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
606  {
607  f.jacobian(x, ret);
608  }
609 
610  template<class FunctionType, class PointType, class HessianRangeType>
611  static void hessian (const FunctionType& f, const PointType &x, HessianRangeType &ret)
612  {
613  f.hessian(x, ret);
614  }
615 
616  };
617 
618  // Unary minus, quite simple
619  template<>
620  struct UnaryFunctionExpressionOperation<MinusOperation>
621  : public FunctionExpressionOperation<MinusOperation>
622  {
623  template<class FunctionType, class PointType, class RangeType>
624  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
625  {
626  f.evaluate(x, ret);
627  ret *= -1.0;
628  }
629 
630  template<class FunctionType, class PointType, class JacobianRangeType>
631  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
632  {
633  f.jacobian(x, ret);
634  ret *= -1.0;
635  }
636 
637  template<class FunctionType, class PointType, class HessianRangeType>
638  static void hessian (const FunctionType& f, const PointType &x, HessianRangeType &ret)
639  {
640  f.hessian(x, ret);
641  for (int k = 0; k < HessianRangeType::dimension; ++k) {
642  ret[k] *= -1.0;
643  }
644  }
645 
646  };
647 
648  // Apply other functions, componentwise. One could in principle
649  // also implement composition, but that would be a binary expression.
650 
651  template<>
652  struct UnaryFunctionExpressionOperation<ExpOperation>
653  : public FunctionExpressionOperation<ExpOperation>
654  {
655  typedef UnaryFunctionExpressionOperation<ExpOperation> ThisType;
656 
657  template<class FunctionType, class PointType, class RangeType>
658  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
659  {
660  f.evaluate(x, ret);
661  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
662  ret[i] = exp(ret[i]);
663  }
664  }
665 
666  template<class FunctionType, class PointType, class JacobianRangeType>
667  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
668  {
669  typename FunctionType::RangeType val;
670  ThisType::evaluate(f, x, val);
671 
672  f.jacobian(x, ret);
673  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
674  ret[i] *= val[i];
675  }
676  }
677 
678  template<class FunctionType, class PointType, class HessianRangeType>
679  static void hessian (const FunctionType& f, const PointType &x, HessianRangeType &ret)
680  {
681  typename FunctionType::RangeType val;
682  ThisType::evaluate(f, x, val);
683 
684  typename FunctionType::JacobianRangeType fjac;
685  f.jacobian(x, fjac);
686  f.hessian(x, ret);
687 
688  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
689  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
690  ret[k][i][i] += fjac[k][i] * fjac[k][i];
691  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
692  typename FunctionType::RangeFieldType tmp;
693  tmp = fjac[k][i] * fjac[k][j];
694  ret[k][i][j] += tmp;
695  ret[k][j][i] += tmp;
696  }
697  }
698  ret[k] *= val;
699  }
700  }
701  };
702 
703  template<>
704  struct UnaryFunctionExpressionOperation<LogOperation>
705  : public FunctionExpressionOperation<LogOperation>
706  {
707  typedef UnaryFunctionExpressionOperation<LogOperation> ThisType;
708 
709  template<class FunctionType, class PointType, class RangeType>
710  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
711  {
712  f.evaluate(x, ret);
713  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
714  ret[i] = log(ret[i]);
715  }
716  }
717 
718  template<class FunctionType, class PointType, class JacobianRangeType>
719  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
720  {
721  typename FunctionType::RangeType fval;
722  f.evaluate(x, fval);
723 
724  f.jacobian(x, ret);
725  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
726  ret[i] /= fval[i];
727  }
728  }
729 
730  template<class FunctionType, class PointType, class HessianRangeType>
731  static void hessian (const FunctionType& f, const PointType &x, HessianRangeType &ret)
732  {
733  typename FunctionType::RangeType fval;
734  f.evaluate(x, fval);
735 
736  typename FunctionType::JacobianRangeType fjac;
737  f.jacobian(x, fjac);
738  f.hessian(x, ret);
739 
740  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
741  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
742  ret[k][i][i] -= fjac[k][i] * fjac[k][i] / fval[k];
743  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
744  typename FunctionType::RangeFieldType tmp;
745  tmp = fjac[k][i] * fjac[k][j] / fval[k];
746  ret[k][i][j] -= tmp;
747  ret[k][j][i] -= tmp;
748  }
749  }
750  ret[k] /= fval[k];
751  }
752  }
753  };
754 
755  template<>
756  struct UnaryFunctionExpressionOperation<SqrtOperation>
757  : public FunctionExpressionOperation<SqrtOperation>
758  {
759  typedef UnaryFunctionExpressionOperation<SqrtOperation> ThisType;
760 
761  template<class FunctionType, class PointType, class RangeType>
762  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
763  {
764  f.evaluate(x, ret);
765  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
766  ret[i] = sqrt(ret[i]);
767  }
768  }
769 
770  template<class FunctionType, class PointType, class JacobianRangeType>
771  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
772  {
773  typename FunctionType::RangeType val;
774  ThisType::evaluate(f, x, val);
775 
776  f.jacobian(x, ret);
777  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
778  ret[i] /= 2.0*val[i];
779  }
780  }
781 
782  template<class FunctionType, class PointType, class HessianRangeType>
783  static void hessian (const FunctionType& f, const PointType &x, HessianRangeType &ret)
784  {
785  typename FunctionType::RangeType fval;
786  f.evaluate(x, fval);
787 
788  typename FunctionType::JacobianRangeType fjac;
789  f.jacobian(x, fjac);
790  f.hessian(x, ret);
791 
792  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
793  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
794  ret[k][i][i] -= fjac[k][i] * fjac[k][i] / 2.0 / fval[k];
795  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
796  typename FunctionType::RangeFieldType tmp;
797  tmp = fjac[k][i] * fjac[k][j] / 2.0 / fval[k];
798  ret[k][i][j] -= tmp;
799  ret[k][j][i] -= tmp;
800  }
801  }
802  ret[k] /= 2.0 / sqrt(fval[k]);
803  }
804  }
805  };
806 
807  template<>
808  struct UnaryFunctionExpressionOperation<SquareOperation>
809  : public FunctionExpressionOperation<SquareOperation>
810  {
811  typedef UnaryFunctionExpressionOperation<SquareOperation> ThisType;
812 
813  template<class FunctionType, class PointType, class RangeType>
814  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
815  {
816  typename FunctionType::RangeType tmp;
817  f.evaluate(x, tmp);
818 
819  ret = tmp * tmp;
820  }
821 
822  template<class FunctionType, class PointType, class JacobianRangeType>
823  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
824  {
825  typename FunctionType::RangeType fval;
826  f.evaluate(x, fval);
827  fval *= 2.;
828 
829  typename FunctionType::JacobianRangeType fjac;
830  f.jacobian(x, fjac);
831 
832  fjac.mtv(fval, ret[0]);
833  }
834 
835  template<class FunctionType, class PointType, class HessianRangeType>
836  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
837  {
838  typename FunctionType::RangeType fval;
839  f.evaluate(x, fval);
840 
841  typename FunctionType::JacobianRangeType fjac;
842  f.jacobian(x, fjac);
843 
844  typename FunctionType::HessianRangeType fhes;
845  f.hessian(x, fhes);
846 
847  ret = 0;
848  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
849  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
850  ret[0][i][i] += 2.0 * (fjac[k][i] * fjac[k][i] + fval[k] * fhes[k][i][i]);
851  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
852  typename FunctionType::RangeFieldType tmp;
853  tmp = 2.0 * (fjac[k][i]*fjac[k][j] + fval[k] * fhes[k][i][j]);
854  ret[0][i][j] += tmp;
855  ret[0][j][i] += tmp;
856  }
857  }
858  }
859  }
860  };
861 
862  template<>
863  struct UnaryFunctionExpressionOperation<SinOperation>
864  : public FunctionExpressionOperation<SinOperation>
865  {
866  typedef UnaryFunctionExpressionOperation<SinOperation> ThisType;
867 
868  template<class FunctionType, class PointType, class RangeType>
869  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
870  {
871  f.evaluate(x, ret);
872  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
873  ret[i] = std::sin(ret[i]);
874  }
875  }
876 
877  template<class FunctionType, class PointType, class JacobianRangeType>
878  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
879  {
880  typename FunctionType::RangeType fval;
881  f.evaluate(x, fval);
882 
883  f.jacobian(x, ret);
884  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
885  ret[k] *= std::cos(fval[k]);
886  }
887  }
888 
889  template<class FunctionType, class PointType, class HessianRangeType>
890  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
891  {
892  typename FunctionType::RangeType fval;
893  f.evaluate(x, fval);
894 
895  typename FunctionType::JacobianRangeType fjac;
896  f.jacobian(x, fjac);
897  f.hessian(x, ret);
898 
899  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
900  typename FunctionType::RangeFieldType cosval, sinval;
901  cosval = std::cos(fval[k]);
902  sinval = std::sin(fval[k]);
903  ret[k] *= cosval;
904  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
905  ret[k][i][i] -= fjac[k][i] * fjac[k][i] * sinval;
906  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
907  typename FunctionType::RangeFieldType tmp;
908  tmp = fjac[k][i] * fjac[k][j] * sinval;
909  ret[k][i][j] -= tmp;
910  ret[k][j][i] -= tmp;
911  }
912  }
913  }
914  }
915  };
916 
917  template<>
918  struct UnaryFunctionExpressionOperation<CosOperation>
919  : public FunctionExpressionOperation<CosOperation>
920  {
921  typedef UnaryFunctionExpressionOperation<CosOperation> ThisType;
922 
923  template<class FunctionType, class PointType, class RangeType>
924  static void evaluate (const FunctionType& f, const PointType &x, RangeType &ret)
925  {
926  f.evaluate(x, ret);
927  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
928  ret[i] = std::cos(ret[i]);
929  }
930  }
931 
932  template<class FunctionType, class PointType, class JacobianRangeType>
933  static void jacobian (const FunctionType& f, const PointType &x, JacobianRangeType &ret)
934  {
935  typename FunctionType::RangeType fval;
936  f.evaluate(x, fval);
937 
938  f.jacobian(x, ret);
939  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
940  ret[k] *= -std::sin(fval[k]);
941  }
942  }
943 
944  template<class FunctionType, class PointType, class HessianRangeType>
945  static void hessian (const FunctionType& f, const PointType &x, HessianRangeType &ret)
946  {
947  typename FunctionType::RangeType fval;
948  f.evaluate(x, fval);
949 
950  typename FunctionType::JacobianRangeType fjac;
951  f.jacobian(x, fjac);
952  f.hessian(x, ret);
953 
954  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
955  typename FunctionType::RangeFieldType cosval, sinval;
956  cosval = std::cos(fval[k]);
957  sinval = std::sin(fval[k]);
958  ret[k] *= -sinval;
959  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
960  ret[k][i][i] += fjac[k][i] * fjac[k][i] * cosval;
961  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
962  typename FunctionType::RangeFieldType tmp;
963  tmp = fjac[k][i] * fjac[k][j] * cosval;
964  ret[k][i][j] -= tmp;
965  ret[k][j][i] -= tmp;
966  }
967  }
968  }
969  }
970  };
971 
972  template<>
973  struct UnaryFunctionExpressionOperation<TanOperation>
974  : public FunctionExpressionOperation<TanOperation>
975  {
976  typedef UnaryFunctionExpressionOperation<TanOperation> ThisType;
977 
978  template<class FunctionType, class PointType, class RangeType>
979  static void evaluate(const FunctionType& f, const PointType &x, RangeType &ret)
980  {
981  f.evaluate(x, ret);
982  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
983  ret[i] = tan(ret[i]);
984  }
985  }
986 
987  template<class FunctionType, class PointType, class JacobianRangeType>
988  static void jacobian(const FunctionType& f, const PointType &x, JacobianRangeType &ret)
989  {
990  typename FunctionType::RangeType fval;
991  f.evaluate(x, fval);
992 
993  f.jacobian(x, ret);
994  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
995  // cos should be declard "pure", so evaluating twice should
996  // be ok
997  ret[k] /= (std::cos(fval[k]) * std::cos(fval[k]));
998  }
999  }
1000 
1001  template<class FunctionType, class PointType, class HessianRangeType>
1002  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
1003  {
1004  typename FunctionType::RangeType fval;
1005  f.evaluate(x, fval);
1006 
1007  typename FunctionType::JacobianRangeType fjac;
1008  f.jacobian(x, fjac);
1009  f.hessian(x, ret);
1010 
1011  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1012  typename FunctionType::RangeFieldType D1, D2;
1013  D1 = 1.0 / (std::cos(fval[k]) * std::cos(fval[k]));
1014  D2 = 2.0 * D1 * tan(fval[k]);
1015  ret[k] *= D1;
1016  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
1017  ret[k][i][i] += fjac[k][i] * fjac[k][i] * D1;
1018  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
1019  typename FunctionType::RangeFieldType tmp;
1020  tmp = fjac[k][i] * fjac[k][j] * D2;
1021  ret[k][i][j] += tmp;
1022  ret[k][j][i] += tmp;
1023  }
1024  }
1025  }
1026  }
1027  };
1028 
1029  template<>
1030  struct UnaryFunctionExpressionOperation<AtanOperation>
1031  : public FunctionExpressionOperation<AtanOperation>
1032  {
1033  typedef UnaryFunctionExpressionOperation<AtanOperation> ThisType;
1034 
1035  template<class FunctionType, class PointType, class RangeType>
1036  static void evaluate(const FunctionType& f, const PointType &x, RangeType &ret)
1037  {
1038  f.evaluate(x, ret);
1039  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
1040  ret[i] = atan(ret[i]);
1041  }
1042  }
1043 
1044  template<class FunctionType, class PointType, class JacobianRangeType>
1045  static void jacobian(const FunctionType& f, const PointType &x, JacobianRangeType &ret)
1046  {
1047  typename FunctionType::RangeType fval;
1048  f.evaluate(x, fval);
1049 
1050  f.jacobian(x, ret);
1051  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1052  // cos should be declard "pure", so evaluating twice should
1053  // be ok
1054  ret[k] /= fval[k] * fval[k] + 1.;
1055  }
1056  }
1057 
1058  template<class FunctionType, class PointType, class HessianRangeType>
1059  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
1060  {
1061  typename FunctionType::RangeType fval;
1062  f.evaluate(x, fval);
1063 
1064  typename FunctionType::JacobianRangeType fjac;
1065  f.jacobian(x, fjac);
1066  f.hessian(x, ret);
1067 
1068  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1069  typename FunctionType::RangeFieldType D1, D2;
1070  D1 = 1.0 / (fval[k] * fval[k] + 1.);
1071  D2 = -2. * fval[k] * D1 * D1;
1072  ret[k] *= D1;
1073  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
1074  ret[k][i][i] += fjac[k][i] * fjac[k][i] * D1;
1075  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
1076  typename FunctionType::RangeFieldType tmp;
1077  tmp = fjac[k][i] * fjac[k][j] * D2;
1078  ret[k][i][j] += tmp;
1079  ret[k][j][i] += tmp;
1080  }
1081  }
1082  }
1083  }
1084  };
1085 
1086  template<>
1087  struct UnaryFunctionExpressionOperation<AsinOperation>
1088  : public FunctionExpressionOperation<AsinOperation>
1089  {
1090  typedef UnaryFunctionExpressionOperation<AsinOperation> ThisType;
1091 
1092  template<class FunctionType, class PointType, class RangeType>
1093  static void evaluate(const FunctionType& f, const PointType &x, RangeType &ret)
1094  {
1095  f.evaluate(x, ret);
1096  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
1097  ret[i] = std::asin(ret[i]);
1098  }
1099  }
1100 
1101  template<class FunctionType, class PointType, class JacobianRangeType>
1102  static void jacobian(const FunctionType& f, const PointType &x, JacobianRangeType &ret)
1103  {
1104  typename FunctionType::RangeType fval;
1105  f.evaluate(x, fval);
1106 
1107  f.jacobian(x, ret);
1108  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1109  // cos should be declard "pure", so evaluating twice should
1110  // be ok
1111  ret[k] /= std::sqrt(1. - fval[k] * fval[k]);
1112  }
1113  }
1114 
1115  template<class FunctionType, class PointType, class HessianRangeType>
1116  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
1117  {
1118  typename FunctionType::RangeType fval;
1119  f.evaluate(x, fval);
1120 
1121  typename FunctionType::JacobianRangeType fjac;
1122  f.jacobian(x, fjac);
1123  f.hessian(x, ret);
1124 
1125  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1126  typename FunctionType::RangeFieldType D1, D2;
1127  D1 = 1.0 / std::sqrt(1. - fval[k] * fval[k]);
1128  D2 = fval[k] * D1 * D1 * D1;
1129  ret[k] *= D1;
1130  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
1131  ret[k][i][i] += fjac[k][i] * fjac[k][i] * D1;
1132  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
1133  typename FunctionType::RangeFieldType tmp;
1134  tmp = fjac[k][i] * fjac[k][j] * D2;
1135  ret[k][i][j] += tmp;
1136  ret[k][j][i] += tmp;
1137  }
1138  }
1139  }
1140  }
1141  };
1142 
1143  template<>
1144  struct UnaryFunctionExpressionOperation<AcosOperation>
1145  : public FunctionExpressionOperation<AcosOperation>
1146  {
1147  typedef UnaryFunctionExpressionOperation<AcosOperation> ThisType;
1148 
1149  template<class FunctionType, class PointType, class RangeType>
1150  static void evaluate(const FunctionType& f, const PointType &x, RangeType &ret)
1151  {
1152  f.evaluate(x, ret);
1153  for (int i = 0; i < FunctionType::RangeType::dimension; ++i) {
1154  ret[i] = std::acos(ret[i]);
1155  }
1156  }
1157 
1158  template<class FunctionType, class PointType, class JacobianRangeType>
1159  static void jacobian(const FunctionType& f, const PointType &x, JacobianRangeType &ret)
1160  {
1161  typename FunctionType::RangeType fval;
1162  f.evaluate(x, fval);
1163 
1164  f.jacobian(x, ret);
1165  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1166  // cos should be declard "pure", so evaluating twice should
1167  // be ok
1168  ret[k] /= -std::sqrt(1. - fval[k] * fval[k]);
1169  }
1170  }
1171 
1172  template<class FunctionType, class PointType, class HessianRangeType>
1173  static void hessian(const FunctionType& f, const PointType &x, HessianRangeType &ret)
1174  {
1175  typename FunctionType::RangeType fval;
1176  f.evaluate(x, fval);
1177 
1178  typename FunctionType::JacobianRangeType fjac;
1179  f.jacobian(x, fjac);
1180  f.hessian(x, ret);
1181 
1182  for (int k = 0; k < FunctionType::RangeType::dimension; ++k) {
1183  typename FunctionType::RangeFieldType D1, D2;
1184  D1 = -1.0 / std::sqrt(1. - fval[k] * fval[k]);
1185  D2 = fval[k] * D1 * D1 * D1;
1186  ret[k] *= D1;
1187  for (int i = 0; i < FunctionType::DomainType::dimension; ++i) {
1188  ret[k][i][i] += fjac[k][i] * fjac[k][i] * D1;
1189  for (int j = i+1; j < FunctionType::DomainType::dimension; ++j) {
1190  typename FunctionType::RangeFieldType tmp;
1191  tmp = fjac[k][i] * fjac[k][j] * D2;
1192  ret[k][i][j] += tmp;
1193  ret[k][j][i] += tmp;
1194  }
1195  }
1196  }
1197  }
1198  };
1199 
1201 
1203 
1204  } // namespace ACFem
1205 
1206 } // namespace Dune
1207 
1208 #endif
Dune::ACFem::log
static UnaryGridFunctionExpression< LogOperation, Function > log(const Fem::Function< typename Function::FunctionSpaceType, Function > &f_)
Component-wise logarithm, log(f)(x) = [log(f(x)_0),...,log(f(x)_N].
Definition: gridfunctionexpression.hh:2844
Dune::ACFem::IdentityOperation
Identity, i.e. just wrap the object.
Definition: expressionoperations.hh:284
Dune::ACFem::InvertOperation
Inversion of an object.
Definition: expressionoperations.hh:275
Dune::ACFem::MultiplyOperation
Multiplication of two objects.
Definition: expressionoperations.hh:248
Dune::ACFem::SMultiplyOperation
Multiplication by scalars from the left.
Definition: expressionoperations.hh:257
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 16, 22:29, 2024)