30#include "rheolef/geo_element_curved_ball.h"
46 for (
size_type iedg = 0, nedg =
gamma.size(1); iedg < nedg; iedg++) {
52 size_type loc_inod = reference_element_e::ilat2loc_inod (
order, ilat(i));
62 case reference_element::t: {
68 size_type loc_inod = reference_element_t::ilat2loc_inod (
order, ilat(i, j));
70 point xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
76 case reference_element::q: {
83 size_type loc_inod = reference_element_q::ilat2loc_inod (
order, ilat(i, j));
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);
139 size_type loc_inod = reference_element_e::ilat2loc_inod (
order, ilat(i));
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");
167 case reference_element::t: {
175 size_type loc_inod = reference_element_t::ilat2loc_inod (
order, ilat(i, j));
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]) {
192 case reference_element::q: {
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);
231 case reference_element::t: {
239 size_type loc_inod = reference_element_t::ilat2loc_inod (
order, ilat(i, j));
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) {
267 size_type loc_inod = reference_element_T::ilat2loc_inod (
order, ilat(i, j, k));
269 node[
dis_inod[loc_inod]] = F (hat_x);
281 size_type loc_inod = reference_element_T::ilat2loc_inod (
order, ilat(i, j, k));
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;
291 case reference_element::q: {
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;
306 size_type loc_inod = reference_element_q::ilat2loc_inod (
order, ilat(i, j));
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"
349int 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
#define error_macro(message)
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_filled(const geo_basic< Float, sequential > &omega)
geo_basic< Float, sequential > fix_s(const geo_basic< Float, sequential > &gamma)
int main(int argc, char **argv)
geo_basic< Float, sequential > fix(const geo_basic< Float, sequential > &omega)
point project_on_unit_sphere(const point &x)
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