DUNE-ACFEM (2.5.1)

functionoperations.hh
1#ifndef __FUNCTIONEXPRESSION_HH__
2#define __FUNCTIONEXPRESSION_HH__
3
4#include "../expressions/expressionoperations.hh"
5
6namespace 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.111.3 (Nov 12, 23:30, 2024)