Mercurial > hg > kwantix
annotate src/main-Laplace.cpp @ 40:eaa99e09607d
Remove indentation due to namespace scope
| author | Jordi Gutiérrez Hermoso <jordigh@gmail.com> |
|---|---|
| date | Fri, 12 Mar 2010 09:21:07 -0600 |
| parents | 7f31f9e2d196 |
| children |
| rev | line source |
|---|---|
| 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; | |
|
34
7f31f9e2d196
Rename kwantxi to kwantix
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
30
diff
changeset
|
31 using namespace kwantix; |
| 0 | 32 |
| 33 //Some arbitrary interior function for Poisson equation (change to | |
| 34 //"return 0*p(1);", for example, for Laplace equation). | |
| 35 class intr_func : public realfunc{ | |
| 36 public: | |
| 37 double at(const point& p) const{return 0*(p(1)+p(2));}; | |
| 38 }; | |
| 39 | |
| 40 //Boundary function for Poisson equation | |
| 41 class bdry_func : public realfunc{ | |
| 42 public: | |
| 43 //sin(theta)^3 | |
| 44 double at(const point& p) const{return | |
| 45 std::pow(std::sin(std::atan2(p(1),p(2))),3);}; | |
| 46 }; | |
| 47 | |
| 48 void output_result(const interpolator<RBF_TYPE>& u); | |
| 49 | |
| 50 int main(){ | |
|
34
7f31f9e2d196
Rename kwantxi to kwantix
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
30
diff
changeset
|
51 using kwantix::vector; |
| 0 | 52 |
| 53 try{ | |
|
34
7f31f9e2d196
Rename kwantxi to kwantix
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
30
diff
changeset
|
54 gsl_set_error_handler(&kwantix::errorHandler); |
| 0 | 55 |
| 56 //Define the domain from given data | |
| 57 shared_ptr<domain> Omega(new domain("data/interior.matrix", | |
| 58 "data/boundary.matrix", | |
|
40
eaa99e09607d
Remove indentation due to namespace scope
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
34
diff
changeset
|
59 "data/normals.matrix")); |
| 0 | 60 |
| 61 //Define the interior and boundary operators | |
| 62 shared_ptr<Laplacian> L(new Laplacian); | |
| 63 shared_ptr<dirichlet_op> B(new dirichlet_op); | |
| 64 | |
| 65 //rhs of BVP | |
| 66 intr_func f; | |
| 67 bdry_func g; | |
| 68 | |
| 69 //Define the Poisson boundary value problem. | |
| 70 shared_ptr<linear_BVP2> Poisson_BVP (new linear_BVP2(Omega, L, B, f, g)); | |
| 71 | |
| 72 //Define interpolator and solve the BVP. | |
| 73 interpolator<RBF_TYPE> u(Poisson_BVP); | |
| 74 | |
| 75 //Write result to file | |
| 76 output_result(u); | |
| 77 | |
| 78 } | |
| 79 catch(error& exc){ | |
|
34
7f31f9e2d196
Rename kwantxi to kwantix
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents:
30
diff
changeset
|
80 kwantix::show_exception(exc); |
| 0 | 81 } |
| 82 | |
| 83 return 0; | |
| 84 } | |
| 85 | |
| 86 void output_result(const interpolator<RBF_TYPE>& u) | |
| 87 { | |
| 88 //Output result as a 41x41 meshgrid (to be plotted, say, with | |
| 89 //Octave; show_sol script included to work with this data). | |
| 90 matrix sol(41,41); | |
| 91 | |
| 92 double r, th; | |
| 93 const double pi = atan(1)*4; | |
| 94 point p(2); | |
| 95 th = 0; | |
| 96 for(size_t i = 1; i <= 41; i++){ | |
| 97 r = 0; | |
| 98 for(size_t j = 1; j <= 41; j++){ | |
| 99 p(1) = r*cos(th); | |
| 100 p(2) = r*sin(th); | |
| 101 sol(i,j) = u(p); //Evaluate interpolator at point p | |
| 102 r += 1.0/40; | |
| 103 } | |
| 104 | |
| 105 th += 2.0*pi/40; | |
| 106 } | |
| 107 | |
| 108 | |
| 109 ofstream ofs("data/sol.matrix"); | |
| 110 ofs << sol; | |
| 111 } |