First partially-implemented functional release
[2dsolver.git] / src / problem.c
1 #include <2dsolve.h>
2
3 struct boundary_condition *
4 boundary_condition_new (int type, float value)
5 {
6 struct boundary_condition *new;
7
8 new = xmalloc (sizeof (struct boundary_condition));
9
10 new->type = type;
11 new->value = value;
12
13 return new;
14 }
15
16 struct problem2d *
17 problem_get_default (struct geometry_info *geom)
18 {
19 struct problem2d *new;
20
21 new = xmalloc (sizeof (struct problem2d));
22
23 memset (new, 0, sizeof (struct problem2d));
24
25 new->geometry = geom;
26
27 jit_compile (new->a[A_11] = jit_expr_new ("1"));
28 jit_compile (new->a[A_22] = jit_expr_new ("1"));
29 jit_compile (new->a[A_21] = jit_expr_new ("0"));
30 jit_compile (new->a[A_12] = jit_expr_new ("0"));
31 jit_compile (new->a0 = jit_expr_new ("1"));
32 jit_compile (new->f = jit_expr_new ("-4.9"));
33
34 new->boundary_count = geom->edge_count;
35 new->boundary_list = xmalloc (geom->edge_count * sizeof (struct boundary_condition *));
36
37 memset (new->boundary_list, 0, geom->edge_count * sizeof (struct boundary_condition *));
38
39 return new;
40 }
41
42 void
43 problem_force_dirichlet (struct problem2d *problem, float value)
44 {
45 int i;
46
47 for (i = 0; i < problem->boundary_count; i++)
48 if (problem->boundary_list[i] == NULL)
49 problem->boundary_list[i] = boundary_condition_new (BOUNDARY_CONDITION_DIRICHLET, value);
50 }
51
52 struct problem2d *
53 problem_get_sample_from_matlab (const char *path)
54 {
55 struct problem2d *problem;
56 struct geometry_info *curr_info;
57
58 if ((curr_info = geometry_load_from_matlab (path, "p", "e", "t")) == NULL)
59 {
60 fprintf (stderr, "cannot open %s: %s\n", path, strerror (errno));
61 return NULL;
62 }
63
64 return problem_get_default (curr_info);
65 }
66
67 #define REALLY_HUGE 1e30
68
69 int
70 problem_is_incomplete (struct problem2d *problem)
71 {
72 int i;
73
74 if (problem->A == NULL)
75 return PROBLEM_INCOMPLETE_NO_MATRIX;
76
77 if (problem->b == NULL)
78 return PROBLEM_INCOMPLETE_NO_VECTOR;
79
80 for (i = 0; i < problem->boundary_count; i++)
81 if (problem->boundary_list[i] == NULL)
82 return PROBLEM_INCOMPLETE_NO_BOUNDS;
83
84 return 0;
85 }
86
87 static void
88 lock_the_right_way (struct problem2d *problem, int var, float value)
89 {
90 int i;
91
92 for (i = 0; i < problem->A->size2; i++)
93 gsl_matrix_set (problem->A, var, i, i == var);
94
95 gsl_vector_set (problem->b, var, value);
96 }
97
98 int
99 problem_lock (struct problem2d *problem)
100 {
101 int code;
102 int i;
103
104 code = problem_is_incomplete (problem);
105
106 if (code == PROBLEM_INCOMPLETE_NO_MATRIX ||
107 code == PROBLEM_INCOMPLETE_NO_VECTOR)
108 {
109 fprintf (stderr, "Cannot lock boundary conditions: stiffness matrix/second member uninitialized\n");
110 return -1;
111 }
112 else if (code == PROBLEM_INCOMPLETE_NO_BOUNDS)
113 {
114 fprintf (stderr, "Cannot lock boundary conditions: some boundaries are undefined\n");
115 return -1;
116 }
117
118 for (i = 0; i < problem->boundary_count; i++)
119 {
120 if (problem->boundary_list[i]->type != BOUNDARY_CONDITION_DIRICHLET)
121 {
122 fprintf (stderr, "Cannot understand Neumann conditions (yet)\n");
123 return -1;
124 }
125
126
127 lock_the_right_way (problem, problem->geometry->edge_list[i]->p1->index, problem->boundary_list[i]->value);
128 lock_the_right_way (problem, problem->geometry->edge_list[i]->p2->index, problem->boundary_list[i]->value);
129
130 /*
131 gsl_matrix_set (problem->A, problem->geometry->edge_list[i]->p1->index, problem->geometry->edge_list[i]->p1->index, REALLY_HUGE);
132 gsl_vector_set (problem->b, problem->geometry->edge_list[i]->p1->index, REALLY_HUGE * problem->boundary_list[i]->value);
133
134 gsl_matrix_set (problem->A, problem->geometry->edge_list[i]->p2->index, problem->geometry->edge_list[i]->p2->index, REALLY_HUGE);
135 gsl_vector_set (problem->b, problem->geometry->edge_list[i]->p2->index, REALLY_HUGE * problem->boundary_list[i]->value);
136 */
137
138 }
139 }
140