Skip to content

Commit a34bf0e

Browse files
gthelergtheler
committed
HAVE_GSL in config.h
1 parent 1307e20 commit a34bf0e

3 files changed

Lines changed: 66 additions & 65 deletions

File tree

‎config.h.ac‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#undef HAVE_ASPRINTF
22
#undef HAVE_CLOCK_GETTIME
33
#undef HAVE_GETRUSAGE
4+
#undef HAVE_GSL
45
#undef HAVE_PETSC
56
#undef HAVE_SLEPC
67
#undef HAVE_SUNDIALS

‎src/math/fit.c‎

Lines changed: 32 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -25,9 +25,9 @@
2525

2626
int feenox_fit_f(const gsl_vector *via, void *arg, gsl_vector *f);
2727
int feenox_fit_df(const gsl_vector *via, void *arg, gsl_matrix *J);
28-
int feenox_fit_in_range(fit_t *fit);
29-
void feenox_fit_update_x(fit_t *fit, size_t j);
30-
void feenox_fit_update_vias(fit_t *fit, const gsl_vector *via);
28+
int feenox_fit_in_range(fit_t *this);
29+
void feenox_fit_update_x(fit_t *this, size_t j);
30+
void feenox_fit_update_vias(fit_t *this, const gsl_vector *via);
3131
void feenox_fit_print_state(const size_t iter, void *arg, const gsl_multifit_nlinear_workspace *w);
3232

3333

@@ -139,32 +139,32 @@ int feenox_instruction_fit(void *arg) {
139139

140140
}
141141

142-
void feenox_fit_update_x(fit_t *fit, size_t j) {
142+
void feenox_fit_update_x(fit_t *this, size_t j) {
143143
unsigned int i = 0;
144-
for (i = 0; i < fit->data->n_arguments; i++) {
145-
fit->x[i] = feenox_vector_get(fit->data->vector_argument[i], j);
146-
feenox_var_value(fit->function->var_argument[i]) = fit->x[i];
144+
for (i = 0; i < this->data->n_arguments; i++) {
145+
this->x[i] = feenox_vector_get(this->data->vector_argument[i], j);
146+
feenox_var_value(this->function->var_argument[i]) = this->x[i];
147147
}
148148
return;
149149
}
150150

151151

152-
void feenox_fit_update_vias(fit_t *fit, const gsl_vector *via) {
152+
void feenox_fit_update_vias(fit_t *this, const gsl_vector *via) {
153153
unsigned int i = 0;
154-
for (i = 0; i < fit->n_via; i++) {
155-
feenox_var_value(fit->via[i]) = gsl_vector_get(via, i);
154+
for (i = 0; i < this->n_via; i++) {
155+
feenox_var_value(this->via[i]) = gsl_vector_get(via, i);
156156
}
157157
return;
158158
}
159159

160-
int feenox_fit_in_range(fit_t *fit) {
160+
int feenox_fit_in_range(fit_t *this) {
161161

162-
if (fit->range_min == NULL || fit->range_max == NULL) {
162+
if (this->range_min == NULL || this->range_max == NULL) {
163163
return 1;
164164
}
165165
unsigned int i = 0;
166-
for (i = 0; i < fit->data->n_arguments; i++) {
167-
if (fit->x[i] < fit->range_min[i] || fit->x[i] > fit->range_max[i]) {
166+
for (i = 0; i < this->data->n_arguments; i++) {
167+
if (this->x[i] < this->range_min[i] || this->x[i] > this->range_max[i]) {
168168
return 0;
169169
}
170170
}
@@ -174,14 +174,14 @@ int feenox_fit_in_range(fit_t *fit) {
174174

175175

176176
int feenox_fit_f(const gsl_vector *via, void *arg, gsl_vector *f) {
177-
fit_t *fit = (fit_t *)arg;
177+
fit_t *this = (fit_t *)arg;
178178

179-
feenox_fit_update_vias(fit, via);
179+
feenox_fit_update_vias(this, via);
180180

181181
size_t j = 0;
182-
for (j = 0; j < fit->n_data; j++) {
183-
feenox_fit_update_x(fit, j);
184-
gsl_vector_set(f, j, feenox_fit_in_range(fit) ? (feenox_function_eval(fit->function, fit->x) - feenox_vector_get(fit->data->vector_value, j)) : 0);
182+
for (j = 0; j < this->n_data; j++) {
183+
feenox_fit_update_x(this, j);
184+
gsl_vector_set(f, j, feenox_fit_in_range(this) ? (feenox_function_eval(this->function, this->x) - feenox_vector_get(this->data->vector_value, j)) : 0);
185185
}
186186
// feenox_debug_print_gsl_vector(f, stdout);
187187

@@ -190,19 +190,19 @@ int feenox_fit_f(const gsl_vector *via, void *arg, gsl_vector *f) {
190190

191191

192192
int feenox_fit_df(const gsl_vector *via, void *arg, gsl_matrix *J) {
193-
fit_t *fit = (fit_t *)arg;
193+
fit_t *this = (fit_t *)arg;
194194

195-
feenox_fit_update_vias(fit, via);
195+
feenox_fit_update_vias(this, via);
196196

197197
unsigned int i = 0;
198198
size_t j = 0;
199199
int in_range = 0;
200-
for (j = 0; j < fit->n_data; j++) {
201-
feenox_fit_update_x(fit, j);
202-
in_range = feenox_fit_in_range(fit);
200+
for (j = 0; j < this->n_data; j++) {
201+
feenox_fit_update_x(this, j);
202+
in_range = feenox_fit_in_range(this);
203203

204-
for (i = 0; i < fit->n_via; i++) {
205-
gsl_matrix_set (J, j, i, (in_range) ? ((fit->gradient != NULL) ? feenox_expression_eval(&fit->gradient[i]) : feenox_expression_derivative_wrt_variable(&fit->function->algebraic_expression, fit->via[i], feenox_var_value(fit->via[i]))) : 0);
204+
for (i = 0; i < this->n_via; i++) {
205+
gsl_matrix_set (J, j, i, (in_range) ? ((this->gradient != NULL) ? feenox_expression_eval(&this->gradient[i]) : feenox_expression_derivative_wrt_variable(&this->function->algebraic_expression, this->via[i], feenox_var_value(this->via[i]))) : 0);
206206
}
207207
}
208208
// feenox_debug_print_gsl_matrix(J, stdout);
@@ -211,15 +211,15 @@ int feenox_fit_df(const gsl_vector *via, void *arg, gsl_matrix *J) {
211211
}
212212

213213
void feenox_fit_print_state(const size_t iter, void *arg, const gsl_multifit_nlinear_workspace *w) {
214-
fit_t *fit = (fit_t *)arg;
214+
fit_t *this = (fit_t *)arg;
215215

216-
if (fit->verbose) {
216+
if (this->verbose) {
217217
gsl_vector *via = gsl_multifit_nlinear_position(w);
218218
gsl_vector *f = gsl_multifit_nlinear_residual(w);
219219

220220
printf ("# iter %2zu\t", iter);
221221
unsigned int i = 0;
222-
for (i = 0; i < fit->n_via; i++) {
222+
for (i = 0; i < this->n_via; i++) {
223223
printf("%+e ", gsl_vector_get(via, i));
224224
}
225225
printf("\tres = %.4g\n", gsl_blas_dnrm2(f));
@@ -243,15 +243,15 @@ int feenox_fit_df(const gsl_vector *via, void *arg, gsl_matrix *J) {
243243
return -1; // GSL_FAILURE
244244
}
245245

246-
int feenox_fit_in_range(fit_t *fit) {
246+
int feenox_fit_in_range(fit_t *this) {
247247
return 0;
248248
}
249249

250-
void feenox_fit_update_x(fit_t *fit, size_t j) {
250+
void feenox_fit_update_x(fit_t *this, size_t j) {
251251
return;
252252
}
253253

254-
void feenox_fit_update_vias(fit_t *fit, const gsl_vector *via) {
254+
void feenox_fit_update_vias(fit_t *this, const gsl_vector *via) {
255255
return;
256256
}
257257

‎src/mesh/interpolate.c‎

Lines changed: 33 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -27,59 +27,59 @@ struct mesh_interp_params {
2727
};
2828

2929

30-
double feenox_mesh_interpolate_function_node(function_t *function, const double *x) {
30+
double feenox_mesh_interpolate_function_node(struct function_t *this, const double *x) {
3131

32-
if (function->vector_value == NULL) {
32+
if (this->vector_value == NULL) {
3333
return 0;
3434
}
3535

3636
// find the nearest node
37-
node_t *nearest_node = feenox_mesh_find_nearest_node(function->mesh, x);
37+
node_t *nearest_node = feenox_mesh_find_nearest_node(this->mesh, x);
3838

3939
// check is it is close enough to a node
40-
switch (function->mesh->dim) {
40+
switch (this->mesh->dim) {
4141
case 1:
42-
if (gsl_pow_2(fabs(x[0]-nearest_node->x[0])) < gsl_pow_2(function->multidim_threshold)) {
43-
return feenox_vector_get(function->vector_value, nearest_node->index_mesh);
42+
if (gsl_pow_2(fabs(x[0]-nearest_node->x[0])) < gsl_pow_2(this->multidim_threshold)) {
43+
return feenox_vector_get(this->vector_value, nearest_node->index_mesh);
4444
}
4545
break;
4646
case 2:
47-
if (feenox_mesh_subtract_squared_module2d(x, nearest_node->x) < gsl_pow_2(function->multidim_threshold)) {
48-
return feenox_vector_get(function->vector_value, nearest_node->index_mesh);
47+
if (feenox_mesh_subtract_squared_module2d(x, nearest_node->x) < gsl_pow_2(this->multidim_threshold)) {
48+
return feenox_vector_get(this->vector_value, nearest_node->index_mesh);
4949
}
5050
break;
5151
case 3:
52-
if (feenox_mesh_subtract_squared_module(x, nearest_node->x) < gsl_pow_2(function->multidim_threshold)) {
53-
return feenox_vector_get(function->vector_value, nearest_node->index_mesh);
52+
if (feenox_mesh_subtract_squared_module(x, nearest_node->x) < gsl_pow_2(this->multidim_threshold)) {
53+
return feenox_vector_get(this->vector_value, nearest_node->index_mesh);
5454
}
5555
break;
5656
}
5757

5858
element_t *element = NULL;
5959
double xi[3] = {0, 0, 0}; // vector with the local coordinates within the element
60-
if ((element = feenox_mesh_find_element(function->mesh, nearest_node, x)) != NULL) {
60+
if ((element = feenox_mesh_find_element(this->mesh, nearest_node, x)) != NULL) {
6161
if (feenox_mesh_interp_solve_for_r(element, x, xi) != FEENOX_OK) {
6262
return 0;
6363
}
6464
} else {
6565
// should we return the nearest node value?
66-
return feenox_vector_get(function->vector_value, nearest_node->index_mesh);
66+
return feenox_vector_get(this->vector_value, nearest_node->index_mesh);
6767
}
6868

6969
// compute the interpolation
7070
double y = 0;
71-
if (function->spatial_derivative_of == NULL) {
71+
if (this->spatial_derivative_of == NULL) {
7272

7373
for (int j = 0; j < element->type->nodes; j++) {
74-
y += element->type->h(j, xi) * feenox_vector_get(function->vector_value, element->node[j]->index_mesh);
74+
y += element->type->h(j, xi) * feenox_vector_get(this->vector_value, element->node[j]->index_mesh);
7575
}
7676

7777
} else {
7878

7979
gsl_matrix *B = feenox_fem_compute_B(element, xi);
8080
for (int j = 0; j < element->type->nodes; j++) {
81-
y += gsl_matrix_get(B, function->spatial_derivative_with_respect_to, j)
82-
* feenox_vector_get(function->spatial_derivative_of->vector_value, element->node[j]->index_mesh);
81+
y += gsl_matrix_get(B, this->spatial_derivative_with_respect_to, j)
82+
* feenox_vector_get(this->spatial_derivative_of->vector_value, element->node[j]->index_mesh);
8383
}
8484
gsl_matrix_free(B);
8585

@@ -90,7 +90,7 @@ double feenox_mesh_interpolate_function_node(function_t *function, const double
9090
}
9191

9292

93-
int feenox_mesh_interp_solve_for_r(element_t *element, const double *x, double *r) {
93+
int feenox_mesh_interp_solve_for_r(element_t *this, const double *x, double *r) {
9494

9595
int gsl_status;
9696
int m;
@@ -102,29 +102,29 @@ int feenox_mesh_interp_solve_for_r(element_t *element, const double *x, double *
102102
gsl_multiroot_function_fdf fun = {&feenox_mesh_interp_residual,
103103
&feenox_mesh_interp_jacob,
104104
&feenox_mesh_interp_residual_jacob,
105-
element->type->dim, &p};
105+
this->type->dim, &p};
106106

107-
if (element->type->id == ELEMENT_TYPE_TETRAHEDRON4 || element->type->id == ELEMENT_TYPE_TETRAHEDRON10) {
107+
if (this->type->id == ELEMENT_TYPE_TETRAHEDRON4 || this->type->id == ELEMENT_TYPE_TETRAHEDRON10) {
108108

109109
// the tetrahedron is an easy one (it's linear)
110110
// usually using this for tet10 is good enough
111111
// TODO: check error
112-
feenox_mesh_compute_r_tetrahedron(element, x, r);
112+
feenox_mesh_compute_r_tetrahedron(this, x, r);
113113

114114
} else {
115115

116-
p.element = element;
116+
p.element = this;
117117
p.x = x;
118118

119-
test = gsl_vector_calloc(element->type->dim); // guess initial cero
119+
test = gsl_vector_calloc(this->type->dim); // guess initial cero
120120

121121
// T = gsl_multiroot_fsolver_hybrids;
122122
// T = gsl_multiroot_fsolver_hybrid;
123123
// T = gsl_multiroot_fsolver_dnewton;
124124
// T = gsl_multiroot_fsolver_broyden;
125125
T = gsl_multiroot_fdfsolver_gnewton;
126126

127-
s = gsl_multiroot_fdfsolver_alloc (T, element->type->dim);
127+
s = gsl_multiroot_fdfsolver_alloc (T, this->type->dim);
128128
gsl_multiroot_fdfsolver_set (s, &fun, test);
129129

130130
do {
@@ -135,7 +135,7 @@ int feenox_mesh_interp_solve_for_r(element_t *element, const double *x, double *
135135
gsl_status = gsl_multiroot_test_residual(s->f, feenox_var_value(feenox.mesh.vars.eps));
136136
} while (gsl_status == -2 /* GSL_CONTINUE */ && iter < 10);
137137

138-
for (m = 0; m < element->type->dim; m++) {
138+
for (m = 0; m < this->type->dim; m++) {
139139
r[m] = gsl_vector_get(gsl_multiroot_fdfsolver_root(s), m);
140140
}
141141

@@ -210,9 +210,9 @@ int feenox_mesh_interp_residual_jacob(const gsl_vector *test, void *params, gsl_
210210

211211

212212

213-
int feenox_mesh_compute_r_tetrahedron(element_t *element, const double *x, double *r) {
213+
int feenox_mesh_compute_r_tetrahedron(element_t *this, const double *x, double *r) {
214214

215-
if (element->type->id == ELEMENT_TYPE_TETRAHEDRON4 || element->type->id == ELEMENT_TYPE_TETRAHEDRON10) {
215+
if (this->type->id == ELEMENT_TYPE_TETRAHEDRON4 || this->type->id == ELEMENT_TYPE_TETRAHEDRON10) {
216216

217217
// reference: eq (9.11) from AFEM.Ch09
218218

@@ -222,9 +222,9 @@ int feenox_mesh_compute_r_tetrahedron(element_t *element, const double *x, doubl
222222
double dz[5][5];
223223
for (int j = 0; j < 4; j++) {
224224
for (int j_prime = 0; j_prime < 4; j_prime++) {
225-
dx[j+1][j_prime+1] = element->node[j]->x[0] - element->node[j_prime]->x[0];
226-
dy[j+1][j_prime+1] = element->node[j]->x[1] - element->node[j_prime]->x[1];
227-
dz[j+1][j_prime+1] = element->node[j]->x[2] - element->node[j_prime]->x[2];
225+
dx[j+1][j_prime+1] = this->node[j]->x[0] - this->node[j_prime]->x[0];
226+
dy[j+1][j_prime+1] = this->node[j]->x[1] - this->node[j_prime]->x[1];
227+
dz[j+1][j_prime+1] = this->node[j]->x[2] - this->node[j_prime]->x[2];
228228
}
229229
}
230230

@@ -233,9 +233,9 @@ int feenox_mesh_compute_r_tetrahedron(element_t *element, const double *x, doubl
233233

234234
// but the nodes and coordinates start from zero
235235
// sixV01 = this->node[1]->x[0] * (this->node[2]->x[1]*this->node[3]->x[2] - this->node[3]->x[1]*this->node[2]->x[2]) + this->node[2]->x[0] * (this->node[3]->x[1]*this->node[1]->x[2] - this->node[1]->x[1]*this->node[3]->x[2]) + this->node[3]->x[0] * (this->node[1]->x[1]*this->node[2]->x[2] - this->node[2]->x[1]*this->node[1]->x[2]);
236-
double sixV02 = element->node[0]->x[0] * (element->node[3]->x[1]*element->node[2]->x[2] - element->node[2]->x[1]*element->node[3]->x[2]) + element->node[2]->x[0] * (element->node[0]->x[1]*element->node[3]->x[2] - element->node[3]->x[1]*element->node[0]->x[2]) + element->node[3]->x[0] * (element->node[2]->x[1]*element->node[0]->x[2] - element->node[0]->x[1]*element->node[2]->x[2]);
237-
double sixV03 = element->node[0]->x[0] * (element->node[1]->x[1]*element->node[3]->x[2] - element->node[3]->x[1]*element->node[1]->x[2]) + element->node[1]->x[0] * (element->node[3]->x[1]*element->node[0]->x[2] - element->node[0]->x[1]*element->node[3]->x[2]) + element->node[3]->x[0] * (element->node[0]->x[1]*element->node[1]->x[2] - element->node[1]->x[1]*element->node[0]->x[2]);
238-
double sixV04 = element->node[0]->x[0] * (element->node[2]->x[1]*element->node[1]->x[2] - element->node[1]->x[1]*element->node[2]->x[2]) + element->node[1]->x[0] * (element->node[0]->x[1]*element->node[2]->x[2] - element->node[2]->x[1]*element->node[0]->x[2]) + element->node[2]->x[0] * (element->node[1]->x[1]*element->node[0]->x[2] - element->node[0]->x[1]*element->node[1]->x[2]);
236+
double sixV02 = this->node[0]->x[0] * (this->node[3]->x[1]*this->node[2]->x[2] - this->node[2]->x[1]*this->node[3]->x[2]) + this->node[2]->x[0] * (this->node[0]->x[1]*this->node[3]->x[2] - this->node[3]->x[1]*this->node[0]->x[2]) + this->node[3]->x[0] * (this->node[2]->x[1]*this->node[0]->x[2] - this->node[0]->x[1]*this->node[2]->x[2]);
237+
double sixV03 = this->node[0]->x[0] * (this->node[1]->x[1]*this->node[3]->x[2] - this->node[3]->x[1]*this->node[1]->x[2]) + this->node[1]->x[0] * (this->node[3]->x[1]*this->node[0]->x[2] - this->node[0]->x[1]*this->node[3]->x[2]) + this->node[3]->x[0] * (this->node[0]->x[1]*this->node[1]->x[2] - this->node[1]->x[1]*this->node[0]->x[2]);
238+
double sixV04 = this->node[0]->x[0] * (this->node[2]->x[1]*this->node[1]->x[2] - this->node[1]->x[1]*this->node[2]->x[2]) + this->node[1]->x[0] * (this->node[0]->x[1]*this->node[2]->x[2] - this->node[2]->x[1]*this->node[0]->x[2]) + this->node[2]->x[0] * (this->node[1]->x[1]*this->node[0]->x[2] - this->node[0]->x[1]*this->node[1]->x[2]);
239239

240240
// back again from one
241241
// xi0 = 1.0/sixV * (sixV01 * 1 + (dy[4][2]*dz[3][2] - dy[3][2]*dz[4][2])*gsl_vector_get(x, 0) + (dx[3][2]*dz[4][2] - dx[4][2]*dz[3][2])*gsl_vector_get(x, 1) + (dx[4][2]*dy[3][2] - dx[3][2]*dy[4][2])*gsl_vector_get(x, 2));
@@ -252,7 +252,7 @@ int feenox_mesh_compute_r_tetrahedron(element_t *element, const double *x, doubl
252252
*/
253253

254254
} else {
255-
feenox_push_error_message("not for element type %d", element->type->id) ;
255+
feenox_push_error_message("not for element type %d", this->type->id) ;
256256
return FEENOX_ERROR;
257257
}
258258

0 commit comments

Comments
 (0)