38 size_t order = gamma.order();
41 for (
size_type iv = 0, nv = gamma.n_vertex(); iv < nv; iv++) {
45 std::vector<size_type> dis_inod;
46 for (
size_type iedg = 0, nedg = gamma.size(1); iedg < nedg; iedg++) {
47 const geo_element& E = gamma.get_geo_element (1, iedg);
48 gamma.dis_inod (E, dis_inod);
49 const point& x0 = node[dis_inod[0]];
50 const point& x1 = node[dis_inod[1]];
53 point xi = x0 + ((1.0*i)/order)*(x1 - x0);
58 for (
size_type ie = 0, ne = gamma.size(2); ie < ne; ie++) {
59 const geo_element& K = gamma.get_geo_element (2, ie);
60 gamma.dis_inod (K, dis_inod);
63 const point& x0 = node[dis_inod[0]];
64 const point& x1 = node[dis_inod[1]];
65 const point& x2 = node[dis_inod[2]];
69 point hat_x ((1.0*i)/order, (1.0*j)/order);
70 point xi = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
77 const point& x0 = node[dis_inod[0]];
78 const point& x1 = node[dis_inod[1]];
79 const point& x2 = node[dis_inod[2]];
80 const point& x3 = node[dis_inod[3]];
84 point hat_x (-1+2.*i/order, -1+2.*j/order);
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;
132 std::vector<size_type> dis_inod;
133 for (
size_type iedg = 0, nedg = omega.size(1); iedg < nedg; iedg++) {
134 const geo_element& E = omega.get_geo_element(1, iedg);
135 omega.dis_inod (E, dis_inod);
136 const point& x0 = node[dis_inod[0]];
137 const point& x1 = node[dis_inod[1]];
140 point xi = x0 + ((1.0*i)/order)*(x1 - x0);
141 if ((
d == 2 && bdry_sid_mark [iedg]) || (
d == 3 && bdry_edg_mark [iedg])) {
144 node[dis_inod[loc_inod]] = xi;
149 for (
size_type ifac = 0, nfac = omega.size(2); ifac < nfac; ifac++) {
150 const geo_element& S = omega.get_geo_element(2, ifac);
151 omega.dis_inod (S, dis_inod);
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");
169 const point& x0 = node[dis_inod[0]];
170 const point& x1 = node[dis_inod[1]];
171 const point& x2 = node[dis_inod[2]];
174 for (
size_type j = 1; i+j < order; j++) {
176 point hat_x ((1.0*i)/order, (1.0*j)/order);
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]) {
186 node[dis_inod[loc_inod]] = xi;
196 default:
error_macro (
"element variant `"<<S.
name()<<
"' not yet supported");
201 for (
size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
203 omega.dis_inod (K, dis_inod);
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);
233 const point& x0 = node[dis_inod[0]];
234 const point& x1 = node[dis_inod[1]];
235 const point& x2 = node[dis_inod[2]];
238 for (
size_type j = 1; i+j < order; j++) {
240 point hat_x ((1.0*i)/order, (1.0*j)/order);
245 x = (1 - hat_x[0] - hat_x[1])*x0 + hat_x[0]*x1 + hat_x[1]*x2;
247 node[dis_inod[loc_inod]] = x;
253 if (is_on_bdry || has_bdry_edge) {
254 const point& x0 = node[dis_inod[0]];
255 const point& x1 = node[dis_inod[1]];
256 const point& x2 = node[dis_inod[2]];
257 const point& x3 = node[dis_inod[3]];
261 }
else if (has_bdry_edge) {
265 for (
size_type j = 1; i+j < order; j++) {
266 for (
size_type k = 1; i+j+k < order; k++) {
268 point hat_x ((1.0*i)/order, (1.0*j)/order, (1.0*k)/order);
269 node[dis_inod[loc_inod]] = F (hat_x);
274 const point& x0 = node[dis_inod[0]];
275 const point& x1 = node[dis_inod[1]];
276 const point& x2 = node[dis_inod[2]];
277 const point& x3 = node[dis_inod[3]];
279 for (
size_type j = 1; i+j < order; j++) {
280 for (
size_type k = 1; i+j+k < order; k++) {
282 point hat_x ((1.0*i)/order, (1.0*j)/order, (1.0*k)/order);
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;
284 node[dis_inod[loc_inod]] = x;
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;
299 const point& x0 = node[dis_inod[shift[0]]];
300 const point& x1 = node[dis_inod[shift[1]]];
301 const point& x2 = node[dis_inod[shift[2]]];
302 const point& x3 = node[dis_inod[shift[3]]];
309 coord [2] = order - i;
310 coord [3] = order - j;
313 point hat_x (-1+2.0*i1/order, -1+2.0*j1/order);
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);
323 node[dis_inod[loc_inod]] = x;
328 default:
error_macro (
"element variant `"<<K.
name()<<
"' not yet supported");
332 omega_out.set_nodes (node);