← Back to branch summary

~fdm-lab/fidlab/main

« back to all changes in this revision

Viewing changes to examples/gsl2d_dir_0.cpp

  • Committer: Apostol Faliagas
  • Date: 2019-02-12 21:52:42 UTC
  • Revision ID: apostol.faliagas@gmail.com-20190212215242-mua2h96aikfvog8y
Published version 0.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples/gsl2d_dir_0.cpp
 
1
#include <iostream>
 
2
#include "pde2d.hpp"
 
3
 
 
4
using namespace fd;
 
5
using namespace fd2;
 
6
 
 
7
const char *project_info =
 
8
R"(//!  Project: Laplace equation with Dirichlet BCs, No. 0
 
9
//!
 
10
//!    y ^
 
11
//!      |
 
12
//!      |  [1] b3
 
13
//!    1 +----------+
 
14
//!      |          |
 
15
//!   [2]|          |[0]
 
16
//!   b2 |          | b1
 
17
//!      |          |
 
18
//!    0 +----------+-----> x
 
19
//!      0  [3] b4  1
 
20
//!
 
21
//! BCs:
 
22
//! b1: $u=f1=u_exact$
 
23
//! b3: $u=f3=u_exact$
 
24
//! b2: $u=f2=u_exact$
 
25
//! b4: $u=f4=u_exact$
 
26
//!
 
27
//! PDE:
 
28
//!     $u_{xx}+u_{yy}=4(x^2+y^2+1)e^{x^2+y^2}$
 
29
//!
 
30
//!     Exact solution:
 
31
//!     $u=e^{x^2+y^2}$
 
32
//!
 
33
//! Boundary discretization method: inward 4-point differencing-scheme (INWARD_DIFF4)
 
34
//! Linear solver: GSL LU
 
35
)";
 
36
 
 
37
int main()
 
38
{
 
39
        try
 
40
        {
 
41
                xBoundary b1(1., 0., 1.);
 
42
                xBoundary b2(0., 1., 0.);
 
43
                yBoundary b3(1., 1., 0.);
 
44
                yBoundary b4(0., 0., 1.);
 
45
 
 
46
                Mesh mesh(50, 50);
 
47
                mesh.set_boundaries({b3, b2, b4 ,b1});
 
48
 
 
49
                DirichletBC bc1(b1, "f1");
 
50
                DirichletBC bc2(b2, "f2");
 
51
                DirichletBC bc3(b3, "f3");
 
52
                DirichletBC bc4(b4, "f4");
 
53
 
 
54
                // Part 1: declaration of fields and functions
 
55
                Field u("u");
 
56
                Function f("f"), uex("exact_solution"), f1("f1"), f2("f2"), f3("f3"), f4("f4");
 
57
 
 
58
                // Part 2: declaration and customization of PDE
 
59
                PDE pde(&mesh);
 
60
                pde.name = "PDE2d";
 
61
 
 
62
                // Part 3: set PDE fields, BCs and functions
 
63
                pde.set_field(u);
 
64
//              pde.set_BC(bc1);
 
65
                pde.set_BC(bc2);
 
66
//              pde.set_BC(bc3, u);
 
67
//              pde.set_BC(bc4, u);
 
68
                pde.set_BCs({bc1, /*bc2, */bc3, bc4});
 
69
                pde.set_functions({f, f1, f2, f3, f4, uex});
 
70
                #if 0
 
71
                // This is the old way
 
72
                pde.get_function("f").implementation = " {return 4*(x*x+y*y+1)*exp(x*x+y*y);}";
 
73
                pde.get_function("f1").implementation = " {return exact_solution(x,y);}";
 
74
                pde.get_function("f2").implementation = " {return exact_solution(x,y);}";
 
75
                pde.get_function("f3").implementation = " {return exact_solution(x,y);}";
 
76
                pde.get_function("f4").implementation = " {return exact_solution(x,y);}";
 
77
                pde.get_function("exact_solution").implementation = " {return exp(x*x+y*y);}";
 
78
                #else
 
79
                f.set_implementation(" {return 4*(x*x+y*y+1)*exp(x*x+y*y);}");
 
80
                f1.set_implementation("exact_solution(x,y)");
 
81
                f2.set_implementation("exact_solution(x,y)");
 
82
                f3.set_implementation("exact_solution(x,y)");
 
83
                f4.set_implementation("exact_solution(x,y)");
 
84
                uex.set_implementation("exp(x*x+y*y)");
 
85
                #endif // 0
 
86
 
 
87
                // Part 4: set PDE equations
 
88
                pde("Dirichlet") = Dxx[u] + Dyy[u] - f;
 
89
 
 
90
                std::cout << "PDE";
 
91
                std::cout << (pde.linear_BCs() ? " has linear BCs" : " has nonlinear BCs") << '\n';
 
92
                std::cout << "PDE";
 
93
                std::cout << (pde.linear_eqs() ? " has linear eqs" : " has nonlinear eqs") << '\n';
 
94
                std::cout << "PDE";
 
95
                std::cout << (pde.linear() ? " is linear." : " is nonlinear.") << '\n';
 
96
 
 
97
                // Part 5: set methods
 
98
                pde.set_library("gsl");
 
99
                pde.set_matrix_method("LU");
 
100
 
 
101
                // Part 6: output
 
102
                pde.executable = "gsl2d_dir_0";
 
103
                pde.compiler = "clang++";
 
104
                pde.project_info = project_info;
 
105
                pde.set_output("code/gsl2d_dir_0");
 
106
                pde.generate_code();
 
107
                pde.generate_makefile();
 
108
        }
 
109
        catch (const std::string& msg)
 
110
        {
 
111
                std::cout << msg << std::endl;
 
112
                exit(1);
 
113
        }
 
114
        catch (const char *msg)
 
115
        {
 
116
                std::cout << "Caught const char* exception!\n";
 
117
                std::cout << msg << std::endl;
 
118
                exit(1);
 
119
        }
 
120
 
 
121
        return 0;
 
122
}