/* See dgtent.h */ /* Last edited on 2004-08-24 14:49:28 by stolfi */ #include #include #include #include #include #include #include #include #include /* INTERNAL PROTOTYPES */ /* IMPLEMENTATIONS */ dg_TentIndex dg_tent_index ( dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_PulseIndex pix[] ) { int i; dg_TentIndex tix = 0; switch(kind) { case DG_TK_H: { affirm(g == 2*c + 1, "bad H tent degree"); int M = c+1, P = 1; for (i = 0; i < d; i++) { affirm(pix[i] > 0, "negative H pulse index"); affirm(pix[i] < M, "bad H pulse index"); tix += pix[i]*P; P *= M; } } break; case DG_TK_B: { affirm(g > 0, "bad B tent degree"); int M = g - c, P = 1; for (i = 0; i < d; i++) { affirm(pix[i] > 0, "negative B pulse index"); affirm(pix[i] < M, "bad B pulse index"); tix += pix[i]*P; P *= M; } } break; case DG_TK_D: { affirm(g == 2*c + 1, "bad D tent degree"); int M = c + 1, P = 1; for (i = 0; i < d; i++) { affirm(pix[i] > 0, "negative D pulse index"); affirm(pix[i] < M, "bad D pulse index"); tix += pix[i]*P; P *= M; M -= pix[i]; } } break; default: affirm(FALSE, "bad tent kind"); } return tix; } dg_TextIndex dg_tent_index_max ( dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g ) { int i; dg_TentIndex tix = 0; switch(kind) { case DG_TK_H: { affirm(g == 2*c + 1, "bad H tent degree"); tix = ipow(c + 1, d) - 1; } break; case DG_TK_B: { affirm(g > 0, "bad B tent degree"); tix = ipow(g - c, d) - 1; } break; case DG_TK_D: { affirm(g == 2*c + 1, "bad D tent degree"); tix = c*ipow(c + 1, d - 1); } break; default: affirm(FALSE, "bad tent kind"); } return tix; } void dg_tent_get_pulse_indices ( dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_TentIndex tix, dg_PulseIndex pix[] ) { int i; switch(kind) { case DG_TK_H: { affirm(g == 2*c + 1, "bad H tent degree"); int M = c + 1; for (i = 0; i < d; i++) { pix[i] = tix % M; tix /= M; } affirm(tix == 0, "bad H tent index"); } break; case DG_TK_B: { affirm(g > 0, "bad B tent degree"); int M = g - c; for (i = 0; i < d; i++) { pix[i] = tix % M; tix /= M; } affirm(tix == 0, "bad B tent index"); } break; case DG_TK_D: { affirm(g == 2*c + 1, "bad D tent degree"); int M = c + 1; for (i = 0; i < d; i++) { pix[i] = tix % M; tix /= M; M -= pix[i]; } affirm(tix == 0, "bad D tent index"); } break; default: affirm(FALSE, "bad tent kind"); } } dg_Locus dg_tent_center ( dg_Dim d, dg_Cont c, dg_Degree g, dg_CellIndex k, dg_PulseIndex pix[] ) { dg_Axes NrmE = 0; int i; for (i = 0; i < d; i++) { if (pix[i] <= c) { NrmE = dg_axis_include(i, NrmE); } } return dg_locus(NrmE, k); } double dg_tent_eval_cell_relative ( dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_TentIndex tix, double *z, dg_GridSize sz[] ) { affirm(g >= 0, "negative degree"); affirm(g <= DG_MAX_B_PULSE_DEGREE, "degree too high"); affirm(g >= 2*c + 1, "overconstrained pulses not implemented yet"); dg_PulseIndex pix[DG_MAX_GRID_DIM]; dg_tent_indices(d, kind, c, g, tix, pix); double v = 1.0; int i; for (i = 0; i < d; i++) { v *= dg_mother_pulse_sum(kind, c, g, pix[i], sz[i], z[i]); } return v; } double dg_tent_eval_root_relative ( dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_Tent t, double *x ) { /* Get the the cell's rank and position: */ dg_Rank rank = dg_cell_rank(t.cell); static dg_GridPos pos[DG_MAX_GRID_DIM]; dg_cell_position(d, t.cell, pos); /* Get the the size of level {rank}: */ static dg_GridSize sz[DG_MAX_GRID_DIM]; dg_grid_size(d, rank, sz); int i; /* Compute the coords {z} of {x} rel to the tent's ref cell {C0}: */ static double z[DG_MAX_GRID_DIM]; for (i = 0; i < d; i++) { z[i] = x[i]*(double)(sz[i]) - pos[i]; } /* Now evaluate the tent function: */ return dg_tent_eval_cell_relative(d, kind, c, g, t.tix, z, sz); } void dg_tent_to_bezier ( dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_TentIndex tix, double **bp ) { /* Data for component mother pulses, unfolded, per axis: */ dg_PulseKind pkind[DG_MAX_GRID_DIM]; /* Pulse kind. */ bz_Degrees gg; /* Degrees of pulse ({g} on all axes). */ dg_PulseIndex pix[DG_MAX_GRID_DIM]; /* Pulse index. */ int wd[DG_MAX_GRID_DIM]; /* Size of tent support. */ double **bt[DG_MAX_GRID_DIM]; /* Bézier coeffs per interval. */ int k, i; dg_tent_indices(d, kind, c, g, tix, pix); /* Number of Bézier coeffs per patch: */ int NB = 1; for (i = 0; i < d; i++) { gg[i] = g; NB *= (g+1); switch (kind) { case DG_TK_H: case DG_TK_D: pkind[i] = DG_PK_H; break; case DG_TK_B: pkind[i] = DG_PK_B; break; default: affirm(FALSE, "bad tent kind"); } wd[i] = 2; } /* Bézier coeffs for a single patch: */ double *bc[NB]; /* ?? unneded ? */ /* Get mother pulses along each axis: */ for (i = 0; i < d; i ++) { pkind[i] = dg_tent_get_pulse_kind(d, kind, t.tix, ???); pix[i] = dg_tent wd[i] = dg_mother_pulse_supp_count(pkind[i], c, g, /* Make sure that the mother pulse in direction {i} spans only 2 intervals. */ for (s = 0; s < 1; s++) { bt[i] = (double *)notnull(malloc()); /* Get Bézier coeffs {bt[i][s,pix]} of interval {s} of unfolded mother spline */ /* */ } } for (k = 0; k < NC; k++) { /* Decompose {k} into displacements {s[0..d-1]} */ int bix; for (bix = 0; bix < NB; bix++) { for (i = 0; i < d; i++) { /* Extract index {pix[i]} from {bix} */ /* Fetch Bezier coeff {pix[i]} of */ dg_tent_get_bezier_patch(d, kind, c, g, t, k, bc) } } /* PRINTOUT: */ void dg_tent_print(FILE *wr, dg_Dim d, dg_TentKind kind, dg_Cont c, dg_Degree g, dg_Tent t) { dg_PulseIndex pix[DG_MAX_GRID_DIM]; dg_tent_indices(d, c, g, t.tix, pix); dg_tent_kind_print(wr, kind); fprintf(wr, ":"); dg_cell_print(wr, d, t.cell); fprintf(wr, "["); int i; for (i = 0; i < d; i++) { fprintf(wr, "%d", (int)pix[i]); } fprintf(wr, "]"); } void dg_tent_kind_print(FILE *wr, dg_TentKind kind) { switch(kind) { case DG_TK_H: fprintf(wr, "H"); break; case DG_TK_B: fprintf(wr, "B"); break; case DG_TK_D: fprintf(wr, "D"); break; default: affirm(FALSE, "bad tent kind"); } } /* Manipulation of tent vectors: */ dg_Tent_vec_t dg_Tent_vec_new(nat nel) { /* This is not a macro only because gcc does not allow cast of struct: */ vec_t v = vec_new(nel, sizeof(dg_Tent)); dg_Tent_vec_t r; r.nel = v.nel; r.el = (dg_Tent *)v.el; return r; }