30 #include "rheolef/geo_element_curved_ball.h"
46 for (
size_type iedg = 0, nedg =
gamma.size(1); iedg < nedg; iedg++) {
70 point xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
85 point xi = 0.25*((1 - hat_x[0])*(1 - hat_x[1])*x0
86 + (1 + hat_x[0])*(1 - hat_x[1])*x1
87 + (1 + hat_x[0])*(1 + hat_x[1])*x2
88 + (1 - hat_x[0])*(1 + hat_x[1])*x3);
94 default:
error_macro (
"element variant `"<<K.
name()<<
"' not yet supported");
98 gamma_out.set_nodes (node);
105 size_t order = omega.order();
109 const size_type sid_dim = bdry.map_dimension();
111 std::vector<bool> bdry_ver_mark (omega.n_node(),
false);
112 std::vector<bool> bdry_edg_mark (omega.size(1),
false);
113 std::vector<bool> bdry_sid_mark (omega.size(sid_dim),
false);
114 for (
size_type ioisid = 0, noisid = bdry.size(); ioisid < noisid; ioisid++) {
115 const geo_element& S = bdry.get_geo_element(sid_dim, ioisid);
116 bdry_sid_mark [S.
dis_ie()] =
true;
117 for (
size_type loc_iv = 0, loc_nv = S.
size(); loc_iv < loc_nv; loc_iv++) {
118 bdry_ver_mark [S[loc_iv]] =
true;
120 if (
d != 3)
continue;
121 for (
size_type loc_iedg = 0, loc_nedg = S.
n_subgeo(1); loc_iedg < loc_nedg; loc_iedg++) {
122 bdry_edg_mark [S.
edge(loc_iedg)] =
true;
127 for (
size_type iv = 0, nv = omega.n_vertex(); iv < nv; iv++) {
128 if (! bdry_ver_mark [iv])
continue;
133 for (
size_type iedg = 0, nedg = omega.size(1); iedg < nedg; iedg++) {
134 const geo_element& E = omega.get_geo_element(1, iedg);
141 if ((
d == 2 && bdry_sid_mark [iedg]) || (
d == 3 && bdry_edg_mark [iedg])) {
149 for (
size_type ifac = 0, nfac = omega.size(2); ifac < nfac; ifac++) {
150 const geo_element& S = omega.get_geo_element(2, ifac);
152 bool side_has_boundary_edge =
false;
154 if (! bdry_sid_mark [ifac]) {
156 for (
size_type loc_iedg = 0, loc_nedg = S.
n_subgeo (1); loc_iedg < loc_nedg; loc_iedg++) {
158 if (bdry_edg_mark [iedg]) {
160 loc_bdry_iedg = loc_iedg;
161 side_has_boundary_edge =
true;
164 check_macro (n_bdry_edg <= 1,
"unsupported side with "<<n_bdry_edg<<
" boundary edges");
178 if (side_has_boundary_edge) {
181 xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
182 if (bdry_sid_mark [ifac]) {
196 default:
error_macro (
"element variant `"<<S.
name()<<
"' not yet supported");
201 for (
size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
207 for (
size_type loc_isid = 0, loc_nsid = K.
n_subgeo (sid_dim); loc_isid < loc_nsid; loc_isid++) {
209 if (bdry_sid_mark [isid]) {
211 loc_bdry_isid = loc_isid;
214 check_macro (n_bdry_sid <= 1,
"unsupported element with "<<n_bdry_sid<<
" boundary sides");
215 bool is_on_bdry = (n_bdry_sid == 1);
219 if (
d == 3 && !is_on_bdry) {
220 for (
size_type loc_iedg = 0, loc_nedg = K.
n_subgeo (1); loc_iedg < loc_nedg; loc_iedg++) {
222 if (bdry_edg_mark [iedg]) {
224 loc_bdry_iedg = loc_iedg;
227 check_macro (n_bdry_edg <= 1,
"unsupported element with "<<n_bdry_edg<<
" boundary edges");
229 bool has_bdry_edge = (n_bdry_edg == 1);
245 x = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
253 if (is_on_bdry || has_bdry_edge) {
261 }
else if (has_bdry_edge) {
269 node[
dis_inod[loc_inod]] = F (hat_x);
283 point x = (1 - hat_x[0] - hat_x[1] - hat_x[2])*x0 + hat_x[0]*x1 + hat_x[1]*x2 + hat_x[2]*x3;
294 shift[0] = loc_bdry_isid;
295 shift[1] = (loc_bdry_isid+1)%4;
296 shift[2] = (loc_bdry_isid+2)%4;
297 shift[3] = (loc_bdry_isid+3)%4;
309 coord [2] =
order - i;
310 coord [3] =
order - j;
318 x = 0.25*( (1 - hat_x[0])*(1 - hat_x[1])*x0
319 + (1 + hat_x[0])*(1 - hat_x[1])*x1
320 + (1 + hat_x[0])*(1 + hat_x[1])*x2
321 + (1 - hat_x[0])*(1 + hat_x[1])*x3);
328 default:
error_macro (
"element variant `"<<K.
name()<<
"' not yet supported");
332 omega_out.set_nodes (node);
336 if (omega.map_dimension() < omega.dimension()) {
337 return fix_s (omega);
343 std::cerr <<
"mkgeo_ball_gmsh_fix: usage:" << std::endl
344 <<
"mkgeo_ball_gmsh_fix "
345 <<
"[-order int] < input.geo > output.geo"
349 int main(
int argc,
char**argv) {
352 argv[0] <<
": command may be used as mono-process only");
356 size_t order = std::numeric_limits<size_t>::max();
357 for (
int i = 1; i < argc; i++) {
359 if (strcmp (argv[i],
"-order") == 0) {
360 if (i == argc-1) { std::cerr << argv[0] <<
" -order: option argument missing" << std::endl;
usage(); }
361 order = atoi(argv[++i]);
369 if (
order != std::numeric_limits<size_t>::max()) {
370 omega.reset_order (
order);