|
0
|
1 //Solves |
|
|
2 // |
|
|
3 // nabla^2 u = f on Omega |
|
|
4 // u = g on dOmega |
|
|
5 // |
|
|
6 //Input data is read from interior.matrix and boundary.matrix |
|
|
7 //files. Each row in interior.matrix corresponds to an interior |
|
|
8 //point of Omega, and similarly for boundary.matrix. Although it's not |
|
|
9 //used, also need a normal vector at each point of the boundary as |
|
|
10 //given by normals.matrix |
|
|
11 |
|
|
12 #include <iostream> |
|
|
13 #include <fstream> |
|
|
14 #include <list> |
|
|
15 #include <cmath> |
|
|
16 |
|
|
17 #include <boost/shared_ptr.hpp> |
|
|
18 |
|
|
19 #include "include/linalg.hpp" |
|
|
20 #include "include/error.hpp" |
|
|
21 #include "include/rbf.hpp" |
|
|
22 #include "include/diff_op.hpp" |
|
|
23 #include "include/bvp.hpp" |
|
|
24 #include "include/interpolator.hpp" |
|
|
25 #include "include/utils.hpp" |
|
|
26 #include "include/func.hpp" |
|
|
27 |
|
|
28 #define RBF_TYPE inverse_multiquadric //Replace by template later. |
|
|
29 |
|
|
30 using namespace std; |
|
|
31 using namespace linalg; |
|
|
32 using namespace rbf; |
|
|
33 using namespace bvp; |
|
|
34 |
|
|
35 //Some arbitrary interior function for Poisson equation (change to |
|
|
36 //"return 0*p(1);", for example, for Laplace equation). |
|
|
37 class intr_func : public realfunc{ |
|
|
38 public: |
|
|
39 double at(const point& p) const{return 0*(p(1)+p(2));}; |
|
|
40 }; |
|
|
41 |
|
|
42 //Boundary function for Poisson equation |
|
|
43 class bdry_func : public realfunc{ |
|
|
44 public: |
|
|
45 //sin(theta)^3 |
|
|
46 double at(const point& p) const{return |
|
|
47 std::pow(std::sin(std::atan2(p(1),p(2))),3);}; |
|
|
48 }; |
|
|
49 |
|
|
50 void output_result(const interpolator<RBF_TYPE>& u); |
|
|
51 |
|
|
52 int main(){ |
|
|
53 using linalg::vector; |
|
|
54 |
|
|
55 try{ |
|
|
56 gsl_set_error_handler(&error_handling::errorHandler); |
|
|
57 |
|
|
58 //Define the domain from given data |
|
|
59 shared_ptr<domain> Omega(new domain("data/interior.matrix", |
|
|
60 "data/boundary.matrix", |
|
|
61 "data/normals.matrix")); |
|
|
62 |
|
|
63 //Define the interior and boundary operators |
|
|
64 shared_ptr<Laplacian> L(new Laplacian); |
|
|
65 shared_ptr<dirichlet_op> B(new dirichlet_op); |
|
|
66 |
|
|
67 //rhs of BVP |
|
|
68 intr_func f; |
|
|
69 bdry_func g; |
|
|
70 |
|
|
71 //Define the Poisson boundary value problem. |
|
|
72 shared_ptr<linear_BVP2> Poisson_BVP (new linear_BVP2(Omega, L, B, f, g)); |
|
|
73 |
|
|
74 //Define interpolator and solve the BVP. |
|
|
75 interpolator<RBF_TYPE> u(Poisson_BVP); |
|
|
76 |
|
|
77 //Write result to file |
|
|
78 output_result(u); |
|
|
79 |
|
|
80 } |
|
|
81 catch(error& exc){ |
|
|
82 utils::show_exception(exc); |
|
|
83 } |
|
|
84 |
|
|
85 return 0; |
|
|
86 } |
|
|
87 |
|
|
88 void output_result(const interpolator<RBF_TYPE>& u) |
|
|
89 { |
|
|
90 //Output result as a 41x41 meshgrid (to be plotted, say, with |
|
|
91 //Octave; show_sol script included to work with this data). |
|
|
92 matrix sol(41,41); |
|
|
93 |
|
|
94 double r, th; |
|
|
95 const double pi = atan(1)*4; |
|
|
96 point p(2); |
|
|
97 th = 0; |
|
|
98 for(size_t i = 1; i <= 41; i++){ |
|
|
99 r = 0; |
|
|
100 for(size_t j = 1; j <= 41; j++){ |
|
|
101 p(1) = r*cos(th); |
|
|
102 p(2) = r*sin(th); |
|
|
103 sol(i,j) = u(p); //Evaluate interpolator at point p |
|
|
104 r += 1.0/40; |
|
|
105 } |
|
|
106 |
|
|
107 th += 2.0*pi/40; |
|
|
108 } |
|
|
109 |
|
|
110 |
|
|
111 ofstream ofs("data/sol.matrix"); |
|
|
112 ofs << sol; |
|
|
113 } |
