2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
6
/// Rheolef is free software; you can redistribute it and/or modify
7
/// it under the terms of the GNU General Public License as published by
8
/// the Free Software Foundation; either version 2 of the License, or
9
/// (at your option) any later version.
11
/// Rheolef is distributed in the hope that it will be useful,
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14
/// GNU General Public License for more details.
16
/// You should have received a copy of the GNU General Public License
17
/// along with Rheolef; if not, write to the Free Software
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20
/// =========================================================================
23
// \int_\Gamma weight(x) dx = mes(\Gamma) with weight
25
// where Gamma is a boundary of the [0,1]^d cube
28
// form_mass_bdr_tst -app P1 -weight 0 -I../data carre-v2 left
30
#include "rheolef/rheolef.h"
31
using namespace rheolef;
34
typedef Float (*function_type)(const point&);
36
Float weight_1 (const point& x) { return 1; }
37
Float weight_x (const point& x) { return x[0]; }
38
Float weight_y (const point& x) { return x[1]; }
39
Float weight_z (const point& x) { return x[2]; }
40
Float weight_x2 (const point& x) { return x[0]*x[0]; }
41
Float weight_xy (const point& x) { return x[0]*x[1]; }
42
Float weight_y2 (const point& x) { return x[1]*x[1]; }
43
Float weight_xz (const point& x) { return x[0]*x[2]; }
44
Float weight_yz (const point& x) { return x[1]*x[2]; }
45
Float weight_z2 (const point& x) { return x[2]*x[2]; }
47
const char * domain_name [6] = {
48
"left", // y=0, xz varies
49
"back", // x=0, yz varies (3d)
50
"bottom", // z=0, xy varies
52
"right", // y=1, xz varies
53
"front", // x=1, yz varies (3d)
54
"top", // z=1, xy varies
57
domain_index (const string& name)
59
for (size_t i = 0; i < 6; i++)
60
if (domain_name[i] == name) return i;
61
error_macro ("unexpected domain name `" << name << "'");
66
function_type function;
74
list_type weight_list_2d[]={
76
// left back bottom right front top
77
// x=0 - y=0 x=1 - y=1
80
{"1", weight_1, {1, -1, 1, 1, -1, 1 } },
81
{"x", weight_x, {0, -1, 0.5, 1, -1, 0.5 } },
82
{"x2",weight_x2, {0, -1, T(1)/3, 1, -1, T(1)/3} },
83
{"y", weight_y, {0.5, -1, 0, 0.5, -1, 1 } },
84
{"xy",weight_xy, {0, -1, 0, 0.5, -1, 0.5 } },
85
{"y2",weight_y2, {T(1)/3, -1, 0, T(1)/3, -1, 1 } }
87
const size_t max_weight_2d = sizeof(weight_list_2d)/sizeof(list_type);
91
list_type weight_list_3d[]={
93
// left back bottom right front top
94
// y=0 x=0 z=0 y=1 x=1 z=1
95
// dx.dz dy.dz dx.dy dx.dz dy.dz dx.dy
97
{"1", weight_1, {1, 1, 1, 1, 1, 1 } },
98
{"x", weight_x, {0.5, 0, 0.5, 0.5, 1, 0.5 } },
99
{"x2",weight_x2, {T(1)/3, 0, T(1)/3, T(1)/3, 1, T(1)/3} },
100
{"y", weight_y, {0, 0.5, 0.5, 1, 0.5, 0.5 } },
101
{"xy",weight_xy, {0, 0, 0.25, 0.5, 0.5, 0.25 } },
102
{"y2",weight_y2, {0, T(1)/3, T(1)/3, 1, T(1)/3, T(1)/3} },
103
{"z", weight_z, {0.5, 0.5, 0, 0.5, 0.5, 1 } },
104
{"xz",weight_xz, {0.25, 0, 0, 0.25, 0.5, 0.5 } },
105
{"yz",weight_yz, {0, 0.25, 0, 0.5, 0.25, 0.5 } },
106
{"z2",weight_z2, {T(1)/3, T(1)/3, 0, T(1)/3, T(1)/3, 1 } }
108
const size_t max_weight_3d = sizeof(weight_list_3d)/sizeof(list_type);
110
// axisymmetric case: TODO recompute the good values because of swap columns
112
list_type weight_list_rz []={
114
// left back bottom right front top
115
// y=0 x=0 z=0 y=1 x=1 z=1
116
// x.dx 0 - x.dx dy -
118
{"1", weight_1, {0.5, 0, -1, 0.5, 1, -1 } },
119
{"x", weight_x, {T(1)/3, 0, -1, T(1)/3, 1, -1 } },
120
{"x2",weight_x2, {0.25, 0, -1, 0.25, 1, -1 } },
121
{"y", weight_y, {0, 0, -1, 0.5, 0.5, -1 } },
122
{"xy",weight_xy, {0, 0, -1, T(1)/3, 0.5, -1 } },
123
{"y2",weight_y2, {0, 0, -1, 0.5, T(1)/3, -1 } }
125
const size_t max_weight_rz = sizeof(weight_list_rz)/sizeof(list_type);
128
weight_index (size_t dim, const string& name, const string& coord_sys)
130
if (coord_sys == "cartesian" && dim == 3) {
132
for (size_t i = 0; i < max_weight_3d; i++)
133
if (weight_list_3d[i].name == name) return i;
134
error_macro ("unexpected 3d weight `" << name << "'");
135
return max_weight_3d;
137
} else if (coord_sys == "cartesian" && dim == 2) {
139
for (size_t i = 0; i < max_weight_2d; i++)
140
if (weight_list_2d[i].name == name) return i;
141
error_macro ("unexpected 2d weight `" << name << "'");
142
return max_weight_2d;
144
} else { // axisymmetric rz
146
for (size_t i = 0; i < max_weight_rz; i++)
147
if (weight_list_rz[i].name == name) return i;
148
error_macro ("unexpected weight `" << name << "' in axisymmetric system");
149
return max_weight_rz;
153
get_function (size_t dim, string weight, string coord_sys)
155
size_t weight_idx = weight_index (dim, weight, coord_sys);
156
if (coord_sys == "cartesian" && dim == 3) {
157
return weight_list_3d [weight_idx].function;
158
} else if (coord_sys == "cartesian" && dim == 2) {
159
return weight_list_2d [weight_idx].function;
161
return weight_list_rz [weight_idx].function;
165
get_expected_value (size_t dim, string weight, string domain_name, string coord_sys)
167
size_t weight_idx = weight_index (dim, weight, coord_sys);
168
size_t domain_idx = domain_index (domain_name);
169
if (coord_sys == "cartesian" && dim == 3) {
170
return weight_list_3d [weight_idx].expect [domain_idx];
171
} else if (coord_sys == "cartesian" && dim == 2) {
172
return weight_list_2d [weight_idx].expect [domain_idx];
174
return weight_list_rz [weight_idx].expect [domain_idx];
179
is_on_axis (const space& W)
181
const point& xmin = W.get_geo().xmin();
182
const point& xmax = W.get_geo().xmax();
183
return (xmin[0] == Float(0) && xmax[0] == Float(0));
189
dcerr << "form_mass_bdr_tst: usage: form_mass_bdr_tst"
194
<< " [-weight string]"
197
<< "example:" << endl
198
<< " form_mass_bdr_tst -app P1 -I../data carre top -weight x" << endl;
201
int main(int argc, char**argv)
203
environment distributed (argc, argv);
204
typedef geo::size_type size_type;
205
Float epsilon = 1e-15;
207
// load geometry and options
210
string approx = "P1";
212
bool mesh_done = false;
213
geo::size_type n_dom = 0;
215
string coord_sys = "cartesian";
218
if (argc <= 1) usage() ;
221
for (int i=1 ; i < argc ; i++) {
223
if (argv [i][0] == '-' && argv [i][1] == 'I')
224
append_dir_to_rheo_path (argv[i]+2) ;
225
else if (strcmp(argv[i], "-app") == 0) approx = argv[++i];
226
else if (strcmp(argv[i], "-weight") == 0) weight = argv[++i];
227
else if (strcmp(argv[i], "-rz") == 0) coord_sys = "rz";
228
else if (strcmp(argv[i], "-") == 0) {
230
// input geo on standard input
231
if (mesh_done) usage() ;
232
dcerr << " ! mass: read geo on stdin" << endl ;
236
} else if (!mesh_done) {
239
omega = geo(argv[i]);
248
warning_macro ("main [1]");
250
warning_macro("no mesh specified.");
254
warning_macro("no weight specified");
257
if (dom_name == "") {
258
warning_macro("no boundary domain specified");
261
size_t dim = omega.dimension();
262
warning_macro ("main [2]");
264
omega.set_coordinate_system(coord_sys);
267
// build boundary domain
269
warning_macro ("main [3]");
270
domain gamma(omega[dom_name]) ;
274
warning_macro ("main [4]");
275
space V(omega, approx);
276
space W(gamma, approx);
280
warning_macro ("main [5]");
283
string sys_coord = omega.coordinate_system();
285
string sys_coord = "cartesian";
286
function_type weight_fct = get_function(dim,weight, sys_coord);
287
field weight_h = interpolate(V,weight_fct);
291
warning_macro ("main [6]");
292
form m_h (V, V, "mass", gamma) ;
293
// ---------------------------------------------------
294
// first check : compute boundary area with i-th weight
295
// ---------------------------------------------------
296
warning_macro ("main [7]");
297
Float mes_gamma_1 = dot(one, m_h*weight_h);
298
cout << "mes1(gamma) = " << double(mes_gamma_1) << endl;
299
// ---------------------------------------------------
300
// 2nd check : use trace
301
// ---------------------------------------------------
303
bool axi_and_on_axis = (coord_sys == "rz" && is_on_axis(W));
305
bool axi_and_on_axis = false;
306
if (!axi_and_on_axis) {
310
warning_macro ("main [8]");
311
field weight_bdr_bug = interpolate(W,weight_fct);
315
warning_macro ("main [9]");
316
form m_gamma_h (W, W, "mass") ;
317
warning_macro ("main [9a]");
318
form m_trace_h (V, W, "mass") ;
320
// deduce field on W by solving
321
// m_gamma_h*weight_bdr = m_trace_h*weight_h
323
warning_macro ("main [10]");
324
field one_bdr (W, 1.);
326
solver sa (m_gamma_h.uu);
328
weight_bdr.u = sa.solve(m_trace_h.uu*weight_h.u
329
+ m_trace_h.ub*weight_h.b - m_gamma_h.ub*weight_bdr.b);
331
// compute boundary area with i-th weight
333
warning_macro ("main [11]");
334
Float mes_gamma_2 = dot(one_bdr, m_gamma_h*weight_bdr);
335
cout << "mes2(gamma) = " << double(mes_gamma_2) << endl;
336
// ---------------------------------------------------
337
// 3nd check : compare two results
338
// ---------------------------------------------------
339
warning_macro ("main [12]");
340
if (fabs(mes_gamma_1-mes_gamma_2) > sqrt(epsilon)) {
341
dcerr << "boundary mass problem detected." << endl;
345
// ---------------------------------------------------
346
// 4th check : compare with exact data
347
// ---------------------------------------------------
348
warning_macro ("main [13]");
349
Float expect = get_expected_value (dim, weight, dom_name, coord_sys);
350
cout << "exact = " << expect << endl;
351
if (fabs(mes_gamma_1-expect) <= sqrt(epsilon)) {
352
dcerr << "ok" << endl;
354
dcerr << "error = " << mes_gamma_1-expect << endl;
357
warning_macro ("main [14]");
361
<< "weight = " << weight_h.u << ";" << endl
362
<< "one_bdr = " << one_bdr.u << ";" << endl
363
<< "m_bdr = " << m_gamma_h.uu << ";" << endl
364
<< "m_trace = " << m_trace_h.uu << ";" << endl
365
<< "weight_bdr = " << weight_bdr.u << ";" << endl
added
removed
nfem/ptst/form_mass_bdr_tst.cc
