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);
field::size_type size_type
see the Float page for the full documentation
see the point page for the full documentation
void set_boundary_face(size_t loc_ifac_curved)
void set_boundary_edge(size_t loc_iedg_curved)
see the disarray page for the full documentation
see the environment page for the full documentation
generic mesh with rerefence counting
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
variant_type variant() const
size_type n_subgeo(size_type subgeo_dim) const
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static const variant_type q
static const variant_type T
static const variant_type t
double Float
see the Float page for the full documentation
idiststream din
see the diststream page for the full documentation
odiststream dout(cout)
see the diststream page for the full documentation
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
geo_basic< Float, sequential > fix_s(const geo_basic< Float, sequential > &gamma)
int main(int argc, char **argv)
point project_on_unit_sphere(const point &x)
geo_basic< Float, sequential > fix_filled(const geo_basic< Float, sequential > &omega)
geo_basic< Float, sequential > fix(const geo_basic< Float, sequential > &omega)
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
This file is part of Rheolef.
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
rheolef - reference manual