PolyFEM
Loading...
Searching...
No Matches
FrankeProblem.cpp
Go to the documentation of this file.
1#include "FrankeProblem.hpp"
2
3#include <iostream>
4
5namespace polyfem
6{
7 namespace problem
8 {
9 namespace
10 {
11 template <typename T>
12 T franke_fun(T x, T y)
13 {
14 auto cx2 = (9 * x - 2) * (9 * x - 2);
15 auto cy2 = (9 * y - 2) * (9 * y - 2);
16
17 auto cx1 = (9 * x + 1) * (9 * x + 1);
18 auto cx7 = (9 * x - 7) * (9 * x - 7);
19
20 auto cy3 = (9 * y - 3) * (9 * y - 3);
21 auto cx4 = (9 * x - 4) * (9 * x - 4);
22
23 auto cy7 = (9 * y - 7) * (9 * y - 7);
24
25 return (3. / 4.) * exp(-(1. / 4.) * cx2 - (1. / 4.) * cy2) + (3. / 4.) * exp(-(1. / 49.) * cx1 - (9. / 10.) * y - 1. / 10.) + (1. / 2.) * exp(-(1. / 4.) * cx7 - (1. / 4.) * cy3) - (1. / 5.) * exp(-cx4 - cy7);
26 }
27
28 template <typename T>
29 T franke_fun(T x, T y, T z)
30 {
31 auto cx2 = (9 * x - 2) * (9 * x - 2);
32 auto cy2 = (9 * y - 2) * (9 * y - 2);
33 auto cz2 = (9 * z - 2) * (9 * z - 2);
34
35 auto cx1 = (9 * x + 1) * (9 * x + 1);
36 auto cx7 = (9 * x - 7) * (9 * x - 7);
37
38 auto cy3 = (9 * y - 3) * (9 * y - 3);
39 auto cx4 = (9 * x - 4) * (9 * x - 4);
40 auto cy7 = (9 * y - 7) * (9 * y - 7);
41
42 auto cz5 = (9 * z - 5) * (9 * z - 5);
43
44 return 3. / 4. * exp(-1. / 4. * cx2 - 1. / 4. * cy2 - 1. / 4. * cz2) + 3. / 4. * exp(-1. / 49. * cx1 - 9. / 10. * y - 1. / 10. - 9. / 10. * z - 1. / 10.) + 1. / 2. * exp(-1. / 4. * cx7 - 1. / 4. * cy3 - 1. / 4. * cz5) - 1. / 5. * exp(-cx4 - cy7 - cz5);
45 }
46
47 template <typename T>
48 T franke_fun_old(T x, T y, T z)
49 {
50 auto cx2 = (9 * x - 2) * (9 * x - 2);
51 auto cy2 = (9 * y - 2) * (9 * y - 2);
52 auto cz2 = (9 * z - 2) * (9 * z - 2);
53
54 auto cx1 = (9 * x + 1) * (9 * x + 1);
55 auto cx7 = (9 * x - 7) * (9 * x - 7);
56
57 auto cy3 = (9 * y - 3) * (9 * y - 3);
58 auto cx4 = (9 * x - 4) * (9 * x - 4);
59 auto cy7 = (9 * y - 7) * (9 * y - 7);
60
61 auto cz5 = (9 * y - 5) * (9 * y - 5);
62
63 return 3. / 4. * exp(-1. / 4. * cx2 - 1. / 4. * cy2 - 1. / 4. * cz2) + 3. / 4. * exp(-1. / 49. * cx1 - 9. / 10. * y - 1. / 10. - 9. / 10. * z - 1. / 10.) + 1. / 2. * exp(-1. / 4. * cx7 - 1. / 4. * cy3 - 1. / 4. * cz5) - 1. / 5. * exp(-cx4 - cy7 - cz5);
64 }
65 } // namespace
66
67 FrankeProblem::FrankeProblem(const std::string &name)
68 : ProblemWithSolution(name)
69 {
70 }
71
72 VectorNd FrankeProblem::eval_fun(const VectorNd &pt, const double t) const
73 {
74 VectorNd res(1);
75
76 if (pt.size() == 2)
77 res(0) = franke_fun(pt(0), pt(1)) * t;
78 else if (pt.size() == 3)
79 res(0) = franke_fun(pt(0), pt(1), pt(2)) * t;
80 else
81 assert(false);
82
83 return res;
84 }
85
86 AutodiffGradPt FrankeProblem::eval_fun(const AutodiffGradPt &pt, const double t) const
87 {
88 AutodiffGradPt res(1);
89
90 if (pt.size() == 2)
91 res(0) = franke_fun(pt(0), pt(1)) * t;
92 else if (pt.size() == 3)
93 res(0) = franke_fun(pt(0), pt(1), pt(2)) * t;
94 else
95 assert(false);
96
97 return res;
98 }
99
100 AutodiffHessianPt FrankeProblem::eval_fun(const AutodiffHessianPt &pt, const double t) const
101 {
102 AutodiffHessianPt res(1);
103
104 if (pt.size() == 2)
105 res(0) = franke_fun(pt(0), pt(1)) * t;
106 else if (pt.size() == 3)
107 res(0) = franke_fun(pt(0), pt(1), pt(2)) * t;
108 else
109 assert(false);
110
111 return res;
112 }
113
115
116 FrankeProblemOld::FrankeProblemOld(const std::string &name)
117 : ProblemWithSolution(name)
118 {
119 }
120
121 VectorNd FrankeProblemOld::eval_fun(const VectorNd &pt, const double t) const
122 {
123 VectorNd res(1);
124
125 if (pt.size() == 2)
126 res(0) = franke_fun(pt(0), pt(1)) * t;
127 else if (pt.size() == 3)
128 res(0) = franke_fun_old(pt(0), pt(1), pt(2)) * t;
129 else
130 assert(false);
131
132 return res;
133 }
134
135 AutodiffGradPt FrankeProblemOld::eval_fun(const AutodiffGradPt &pt, const double t) const
136 {
137 AutodiffGradPt res(1);
138
139 if (pt.size() == 2)
140 res(0) = franke_fun(pt(0), pt(1)) * t;
141 else if (pt.size() == 3)
142 res(0) = franke_fun_old(pt(0), pt(1), pt(2)) * t;
143 else
144 assert(false);
145
146 return res;
147 }
148
149 AutodiffHessianPt FrankeProblemOld::eval_fun(const AutodiffHessianPt &pt, const double t) const
150 {
151 AutodiffHessianPt res(1);
152
153 if (pt.size() == 2)
154 res(0) = franke_fun(pt(0), pt(1)) * t;
155 else if (pt.size() == 3)
156 res(0) = franke_fun_old(pt(0), pt(1), pt(2)) * t;
157 else
158 assert(false);
159
160 return res;
161 }
162 } // namespace problem
163} // namespace polyfem
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
x
int x
Definition SplineBasis3d.cpp:55
polyfem::problem::FrankeProblem::eval_fun
VectorNd eval_fun(const VectorNd &pt, const double t) const override
Definition FrankeProblem.cpp:72
polyfem::problem::FrankeProblem::FrankeProblem
FrankeProblem(const std::string &name)
Definition FrankeProblem.cpp:67
polyfem::problem::FrankeProblemOld::eval_fun
VectorNd eval_fun(const VectorNd &pt, const double t) const override
Definition FrankeProblem.cpp:121
polyfem::problem::FrankeProblemOld::FrankeProblemOld
FrankeProblemOld(const std::string &name)
Definition FrankeProblem.cpp:116
generate_rotation_mtx.T
T
Definition generate_rotation_mtx.py:43
polyfem::AutodiffHessianPt
Eigen::Matrix< AutodiffScalarHessian, Eigen::Dynamic, 1, 0, 3, 1 > AutodiffHessianPt
Definition AutodiffTypes.hpp:14
polyfem::VectorNd
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > VectorNd
Definition Types.hpp:11
polyfem::AutodiffGradPt
Eigen::Matrix< AutodiffScalarGrad, Eigen::Dynamic, 1, 0, 3, 1 > AutodiffGradPt
Definition AutodiffTypes.hpp:13