Skip to content
Snippets Groups Projects
Commit cc564534 authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Fix up CurvedFunction

parent 8f71e42d
No related branches found
No related tags found
No related merge requests found
#ifndef CURVED_FUNCTION_HH #ifndef CURVED_FUNCTION_HH
#define CURVED_FUNCTION_HH #define CURVED_FUNCTION_HH
#include <cmath>
#include <dune/fufem/interval.hh> #include <dune/fufem/interval.hh>
namespace Dune { namespace Dune {
...@@ -35,16 +36,14 @@ template <class NonlinearityType> class CurvedFunction { ...@@ -35,16 +36,14 @@ template <class NonlinearityType> class CurvedFunction {
} }
void domain(Interval<double> &domain) const { void domain(Interval<double> &domain) const {
// TODO domain[0] = -M_PI;
domain[0] = 0; domain[1] = M_PI;
domain[1] = 1;
} }
void cartesian(double m, VectorType &y) const { void cartesian(double m, VectorType &y) const {
y = 0; y = 0;
y.axpy(1 - m, x); y.axpy(std::cos(m), x);
y.axpy(m, dir); y.axpy(std::sin(m), dir);
y /= y.two_norm();
} }
private: private:
...@@ -54,17 +53,18 @@ template <class NonlinearityType> class CurvedFunction { ...@@ -54,17 +53,18 @@ template <class NonlinearityType> class CurvedFunction {
VectorType const &x; VectorType const &x;
VectorType const &dir; VectorType const &dir;
/* /* If x and d were normalised, we would have
Tangential vector within the plane, with positive direction.
< (1-m)x + m*y,-m*|y|^2*x+(1-m)*|x|^2*y> cartesian(a) = cos(a) * x + sin(a) * d and
= -m(1-m)|y|^2<x,x> + m(1-m)|x|^2<y,y> (because <x,y> = 0) tangent(a) = -sin(a) * x + cos(a) * d.
= 0
Since we x and d are not normalised and the return of
cartesian() is fixed, we scale the tangent.
*/ */
void tangentialDirection(double m, VectorType &y) const { void tangentialDirection(double m, VectorType &y) const {
y = 0; y = 0;
y.axpy(-m * dir.two_norm2(), x); y.axpy(-std::sin(m) * dir.two_norm2(), x);
y.axpy((1 - m) * x.two_norm2(), dir); y.axpy(std::cos(m) * x.two_norm2(), dir);
} }
}; };
} }
......
...@@ -35,18 +35,18 @@ int main() { ...@@ -35,18 +35,18 @@ int main() {
start[0] = 0; start[0] = 0;
start[1] = 1; start[1] = 1;
double const ret1 = functionTester(J, start, 3); double const ret1 = functionTester(J, start, 5);
// Something random // Something random
start[0] = 279; start[0] = 279;
start[1] = -96; start[1] = -96;
double const ret2 = functionTester(J, start, 3); double const ret2 = functionTester(J, start, 5);
assert(std::abs(ret1 - ret2) < 1e-5); assert(std::abs(ret1 - ret2) < 1e-7);
start[0] = 0; start[0] = 0;
start[1] = 0; start[1] = 0;
double const ret3 = functionTester(J, start, 1); double const ret3 = functionTester(J, start, 1);
assert(std::abs(ret1 - ret3) < 1e-5); assert(std::abs(ret1 - ret3) < 1e-7);
} }
...@@ -41,12 +41,12 @@ int main() { ...@@ -41,12 +41,12 @@ int main() {
start[0] = 279; start[0] = 279;
start[1] = -96; start[1] = -96;
double const ret2 = functionTester(J, start, 3); double const ret2 = functionTester(J, start, 5);
assert(std::abs(ret1 - ret2) < 1e-5); assert(std::abs(ret1 - ret2) < 1e-7);
start[0] = 0; start[0] = 0;
start[1] = 0; start[1] = 0;
double const ret3 = functionTester(J, start, 1); double const ret3 = functionTester(J, start, 2);
assert(std::abs(ret1 - ret3) < 1e-5); assert(std::abs(ret1 - ret3) < 1e-7);
} }
...@@ -41,12 +41,12 @@ int main() { ...@@ -41,12 +41,12 @@ int main() {
start[0] = 279; start[0] = 279;
start[1] = -96; start[1] = -96;
double const ret2 = functionTester(J, start, 4); double const ret2 = functionTester(J, start, 5);
assert(std::abs(ret1 - ret2) < 1e-8); assert(std::abs(ret1 - ret2) < 1e-7);
start[0] = 0; start[0] = 0;
start[1] = 0; start[1] = 0;
double const ret3 = functionTester(J, start, 1); double const ret3 = functionTester(J, start, 2);
assert(std::abs(ret2 - ret3) < 1e-5); assert(std::abs(ret1 - ret3) < 1e-7);
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment