/* Last edited on 2005-02-04 16:46:08 by stolfi */ /* Datapoint interpolation with 2D dyadic splines. */ /* This program takes a set of points {p[0..n-1]}, corresponding data values {f[0..n-1]}, and derivatives, and computes a dyadic spline {h} such that {h(p[k]) = f[k]} for all {k}. Each sample point must be a vertex of the {d}-dimensional dyadic multigrid, specified by the index {c[k]} of its superior cell. Note that this input also specifies a rank {r[k]} for each point {p[k]}. The spline {h} is defined over the smallest finite dyadic grid {G} where every data point {p[k]} is a complete vertex in layer {r[k]} of {G}. The interpolant has the form h = SUM { u[i]*bas[i] : i = 0..m-1 } (1) where {bas[0..m-1]} consists of all the maximal (or minimal) polynomial tents of continuity order {c} and degree {g} over {G} that are centered on the vertices {p[0..n-1]} and have ranks {r[0..n-1]}. The coefficients {u[i]} are found by the ordered interpolation algorithm (see C.Cardoso's thesis). */ #include #include #include #include #include #include #include #include #include #include #include #include /* Max dimension of grid's domain: */ #define MAX_DDIM (4) /* Max dimension of function's range space: */ #define MAX_FDIM (4) /* TYPES */ int main(int argc, char **argv) { tin_options_t *o = tin_get_options(argn, argc); int n = o->nPts; int dDim = o->dDim; int fDim = o->fDim; int c = o->cont; int g = o->deg; dg_Node *rootG = tin_make_trivial(); int i; /* List of loci (vertices + depths) for interpolation: */ dg_Locus_vec L = dg_Locus_vec_new(n); /* Number of function values per tent: */ int tSize = ipow(c + 1, dDim); /* Number of function values per data point: */ int fSize = fDim*tSize; /* The function values: */ double f[fSize*n]; /* {f[i*fSize + k*fDim + r]} is derivative {k} of f-coord {r} of point {i} */ for (i = 0; i < n; i++) { /* Read a sampling point and its function values: */ tap_read(&(L.el[i]), &(f[i*fSize])); /* Make sure it is a complete vertex in {G}: */ tin_make_locus_complete(rootG, L.el[i]); } /* Sort vertices by increasing depth: */ tin_sort_loci(&L); /* Build tent basis for those vertices */ dg_Tent_vec_t bas = dg_tent_basis(dDim, DG_TK_H, c, g, L); assert(bas.nel == tDim*n); /* Coefficient vector: */ double u[fDim*bas.nel]; /* {u[i*fSize + k*fDim + r]} is coeff in f-coord {r} of tent {k} assoc to point {i} */ /* Interpolation algorithm: */ for (i = 0; i < L.nel; i++) { /* Get the locus (vertex) {Li}: */ dg_Locus Li = L.el[i]; /* Get the root-relative bbox {B} of its superior cell: */ dg_Interval B[dDim]; dg_cell_box_root_relative(dDim, Li.cell, B); /* Extract the lower corner {x} (root-relative vertex coords): */ double x[dDim]; for (r=0; r < dDim; r++) { x[i] = B[i].end[0]; } /* Get the function values {fi} specified for it: */ double *fi = &(f[i*fSize]); /* Add values and derivatives {s[...]} of all tents for all prev loci: */ double s[fSize]; for (kr = 0; kr < fSize; kr++) { s[kr] = 0; } for (j = 0; j < i; j++) { for (tix = 0; tix < tSize; tix++) { dg_Tent t = bas.el[j*tSize + tix]; assert(t.tix = tix); assert(t.cell = Li.cell); /* Coefficients for this tent: */ double *ujt = &(u[j*fSize + tix*fDim]); /* Get value and all derivs of {t} at {x}: */ double v[tDim]; dg_tent_eval_diff_root_relative(dDim, DG_TK_H, c, g, t, x, v); /* Accumulate on {s}: */ for (k = 0; k < tDim, k++) { for (r = 0; r < fDim; r++) { s[k*fDim + r] += ujt[r] * v[k*fDim + r]; } } } } /* Assume that the tents are cardinal with respect to derivative: */ /* Subtract {s} from specified values {fi}, result is coeffs {u}: */ for (kr = 0; kr < fSize; kr++) { ui[kr] = fi[kr] - s[kr]; } } } dg_Locus_vec_t tin_locus_sort(dg_CellIndex *p, int n) { /* Fills {ix[0..n-1]} with a permutation of {0..n-1} such that {i \leq j} implies {rank(p[ix[i]]) \leq rank(p[ix[j]])}. */ int ix[n]; for (i = 0; i < n; i++) { ix[i] = i; } auto int cmp_depth(int x, int y); /* Compares cell indices {c[x]} with {c[y]} in depth order. */ int cmp_depth(int x, int y) { if (c[x] < c[y]) { return -1; } else if (c[x] > c[y]) { return +1; } else { return 0; } } isrt_heapsort(ix, n, cmp_depth, +1); } /* Number of cells in grid: */ long int minCells = ipow(2, o->minDepth) - 1; long int maxCells = ipow(2, o->maxDepth) - 1; long int totCells = o->nCells; if (totCells < minCells) { totCells = minCells; } if (totCells > maxCells) { totCells = maxCells; } /* Cell list and heap */ int nCells = 0; dg_Node *node[totCells]; /* Cells are {node[0..nCells-1]}. */ /* Heap of indices into {node}, sorted by badness: */ int nLeaves = 0; int ix[totCells]; /* {node[ix[0]]} is the largest leaf cell. */ /* Start with a trivial tree: */ dg_Node *root = tin_make_trivial_tree(); node[0] = root; nCells = 1; ix[0] = 0; nLeaves = 1; /* Break leaves until we have {totCells} cells: */ while (nCells < totCells) { /* Get leaf cell with largest badness: */ int ixw = ihp_heap_pop(ix, &nLeaves, cmp_badness, +1); dg_Node *w = node[ixw]; /* Split it and add children to cell list: */ dg_node_split(w); for (ich = 0; ich < 2; ich++) { dg_Node *ch = w->ch[ich]; int ixch = nCells; nCells++; node[ixch] = ch; ihp_heap_insert(ix, &nLeaves, ixch, cmp_badness, +1); } } /* Find complete vertices: */ dg_Locus_vec_t vert = dg_find_critical_loci(DDIM, root, TRUE, o->wide); /* Build tent basis. */ dg_Tent_vec_t bas = dg_tent_basis(DDIM, DG_PK_B, 0, 1, vert) int SDIM = bas.nel; /* Write basis? */ /* tin_write_basis(bas, o->outName); */ /* Compute the matrix {H[i,j] = (lower half only). */ smt_matrix_t H = smt_matrix_new(SDIM*15); /* Initial guess. */ for (i = 0; i < SDIM; i++) { dg_Tent ti = bas.el[i]; for (j = 0; j <= i; j++) { dg_Tent tj = bas.el[j]; double Hij = dg_tent_dot(DDIM, DG_BP_K, 0, 1, ti,tj); smt_append(&H, i, j, Hih); } } dg_interval rootBox[DDIM]; /* absolute bounds of root cell */ /* Compute the right-hand-side vector(s): */ double_vec b = double_vec_new(FDIM * SDIM); for (i = 0; i < SDIM; i++); { dg_Tent ti = bas.el[i]; /* Position of the tent's reference cell (indices in grid layer): */ dg_GridPos refPos[DDIM]; dg_cell_position(DDIM, ti.cell, refPos); /* Position of the integration cell {C} for the {integrand} function: */ dg_GridPos cellPos[DDIM]; /* Absolute bounds of cell {C}: */ dg_Interval cellBox[DDIM]; /* Work vectors for {integrand}: */ double y[DDIM]; /* Point coords relative to the tent's reference cell. */ double x[DDIM]; /* Absolute coordinates. */ auto void integrand(double *z); void integrand(double *z, double *f) { /* {z} is relative to the current cell {C}. */ int i; /* Compute coordinates {y} relative to the tent's reference cell: */ for (i = 0; i < DDIM; i++) { double zi = z[i] y[i] = zi + (cellPos[i] - refPos[i]); x[i] = (1-zi)*cellBox[i][0] + zi*cellBox[i][1]; } double tx = dg_tent_eval_cell_relative(DDIM, DG_PK_B, 0, 1, ti, y); tin_eval_function(func, DDIM, x, FDIM, f); for (i = 0; i < FDIM; i++) { f[i] *= t; } } int nc = dg_tent_domain_size(DDIM, DG_PK_B, 0, 1, ti); int ic; for (k = 0; k < FDIM; k++) { b.el[i*FDIM+k] = 0; } for (ic = 0; ic < nc; ic++) { /* Get index and position of the cell number {ic} in the domain of {ti}: */ dg_CellIndex k = dg_tent_domain_cell(DDIM, DG_PK_B, 0, 1, ti, ic); dg_cell_position(ddim, k, cellPos); /* Get the cell's box in absolute coordinates: */ dg_cell_box_root_relative(DDIM, k, cellBox); dg_box_box_map(DDIM, cellBox, rootBox, cellBox); /* Integrate over that cell and accumulate: */ double res[FDIM]; gauss_integrate(DDIM, FDIM, o->gaussOrder, integrand, res); for (k = 0; k < FDIM; k++) { b[i*FDIM+k] += res[k]; } } } /* Allocate auxiliary vectors: */ double_vec u = double_vec_new(FDIM * SDIM); double_vec uj = double_vec_new(SDIM); // Column {j} of {u}. double_vec bj = double_vec_new(SDIM); // Column {j} of {b}. double_vec y = double_vec_new(SDIM); if (cholesky) { smt_matrix_t C = smt_cholesky(H, 0.0); tin_Cholesky_solve(C, FDIM, b, u, bj, uj, y, TRUE); } /* else if (conjGrad) { smt_matrix_t G = SOApprox_GetBasisMatrix(o->matName, FALSE); SOApprox_GuessSol(v); SOApprox_ConjugateGradientSolve(G, b, v, u, bj, uj, TRUE); } */ else if (o->gaussSeidel) { tin_guess_solution(v); tin_insert_upper_half(H); tin_Gauss_Seidel_solve ( H, b, o->omega, o->maxIter, o->absTol, o->relTol, v, u, bj, uj, TRUE ); } else { assert(FALSE , "unspecified/invalid solution method"); } /* Output final solution: */ double fMax[FDIM], eMax[FDIM]; /* Compute the spline: */ tin_build_spline(bas, u, root); /* The basis approximation: */ auto void appr (double *x, dg_Dim d, double *f, dg_Dim fDim); void appr (double *x, dg_Dim d, double *f, dg_Dim fDim) { dg_Node *nd = tin_find_leaf_cell(root, x); bz_Coeffs_t *bz = (bz_Coeffs_t *)nd->bz; bz_eval(bz, x, f); } tin_print_max_error_values(func, appr, fMax, eMax); /* 2D Plotting */ if((o->plot)&&(DDIM == 2)&&(FDIM == 1)) { /* Builds error function: */ Basis error_bas = SOFunctionRef_vec_new(2); double_vec errc = double_vec_new(2); errc.el[0] = 1; errc.el[1] = -1; error_bas.el[0] = (SOFunction *)f; error_bas.el[1] = (SOFunction *)g; SOFunction *e = (SOFunction *)SOLinCombFunction_Make(DDIM, FDIM,"Errorf.bas", error_bas, errc); double fPlotMin = 0, fPlotMax = 0; double gPlotMin = 0, gPlotMax = 0; double ePlotMin = 0, ePlotMax = 0; /* Plots the true solution {f}, the approximation {g}, and the error {e}: */ SOApprox_PlotSingleFunction(f, &NoMap, tree, txtcat(o->outName, "-sol"), &(o->plt), &fPlotMin, &fPlotMax); SOApprox_PlotSingleFunction(g, &NoMap, tree, txtcat(o->outName, "-apr"), &(o->plt), &gPlotMin, &gPlotMax); /* Cortes 1D. */ double initP[2] = {0.0, 0.5}, finP[2] = {1.0, 0.5}, D1min, D1max; SO1DPlot_SingleFunction (TRUE, NULL, 300, 200, g, &NoMap, DDIM, initP, finP, -0.5, 12.0, 8, 0.10, "corteX", &D1min, &D1max); initP[0] = 0.25; initP[1] = 0.0; finP[0] = 0.25; finP[1] = SQRTHALF; SO1DPlot_SingleFunction (TRUE, NULL, 300, 200, g, &NoMap, DDIM, initP, finP, -0.5, 12.0, 8, 0.10, "corteY", &D1min, &D1max); o->plt.fRange = o->plt.fRange / 10.0; o->plt.fStep = o->plt.fStep / 10.0; SOApprox_PlotSingleFunction(e, &NoMap, tree, txtcat(o->outName, "-err"), &(o->plt), &ePlotMin, &ePlotMax); fprintf(stderr, "observed SOLUTION function extrema:"); fprintf(stderr, " min = %16.12f max = %16.12f\n", fPlotMin, fPlotMax); fprintf(stderr, "\n"); fprintf(stderr, "observed APPROXIMATION function extrema:"); fprintf(stderr, " min = %16.12f max = %16.12f\n", gPlotMin, gPlotMax); fprintf(stderr, "\n"); fprintf(stderr, "observed ERROR function extrema:"); fprintf(stderr, " min = %16.12f max = %16.12f\n", ePlotMin, ePlotMax); fprintf(stderr, "\n"); fprintf(stderr, "Basis with %d tents; Grid with %d leaf cells.\n", SDIM, leaf_cells); double distance = 0.0; FILE *rd2 = open_read("../SOFindTentBasis/regular10.tree"); SOGrid_Tree *tree10 = SOGrid_Tree_read(rd2); SOFunction_Dot(e, &NoMap, e, &NoMap, (SOFunction *)NULL, tree10, &distance, TRUE); distance = sqrt(distance); fprintf(stderr, "observed DISTANCE ( |f-g| = sqrt() ): %.16g \n", distance); } return 0; } #define PPUSAGE SOParams_SetUsage tin_options_t *tin_get_options(int argn, char **argc) { tin_options_t *o = (tin_options_t *)notnull(malloc(sizeof(tin_options_t)), "no mem"); SOParams_T *pp = SOParams_NewT(stderr, argn, argc); PPUSAGE(pp, "SOComputeApprox \\\n"); PPUSAGE(pp, " -funcName NAME \\\n"); PPUSAGE(pp, " -matName NAME \\\n"); PPUSAGE(pp, " -minRank NUM \\\n"); PPUSAGE(pp, " -maxRank NUM \\\n"); PPUSAGE(pp, " -diff_limit NUM \\\n"); PPUSAGE(pp, " [ -wide ] \\\n"); PPUSAGE(pp, " [ -cholesky | \\\n"); PPUSAGE(pp, " -gaussSeidel [-omega NUM] [-tol NUM] [-maxIter NUM] | \\\n"); PPUSAGE(pp, " -conjGrad ] \\\n"); PPUSAGE(pp, " [ -gaussOrder NUM ] \\\n"); PPUSAGE(pp, " -outName NAME \\\n"); /* Plotting options instructions. */ PPUSAGE(pp, " [ -plot ] \\\n"); PPUSAGE(pp, SOPlotParams_FunctionHelp " \n"); SOParams_GetKeyword(pp, "-funcName"); o->funcName = SOParams_GetNext(pp); SOParams_KeywordPresent(pp, "-matName"); o->matName = SOParams_GetNext(pp); SOParams_GetKeyword(pp, "-outName"); o->outName = SOParams_GetNext(pp); cholesky = SOParams_KeywordPresent(pp, "-cholesky"); o->gaussSeidel = SOParams_KeywordPresent(pp, "-gaussSeidel"); conjGrad = SOParams_KeywordPresent(pp, "-conjGrad"); if (! (cholesky || o->gaussSeidel || conjGrad)) { SOParams_Error(pp, "must specify a linear solution method"); } else if ((cholesky + o->gaussSeidel + conjGrad) > 1) { SOParams_Error(pp, "please specify only one linear solution method"); } SOParams_KeywordPresent(pp, "-minRank"); o->minRank = SOParams_GetNextInt(pp, 0, INT_MAX); SOParams_KeywordPresent(pp, "-maxRank"); o->maxRank = SOParams_GetNextInt(pp, 0, INT_MAX); SOParams_KeywordPresent(pp, "-diff_limit"); o->diff_limit = SOParams_GetNextDouble(pp, 0, MAXDOUBLE); o->wide = SOParams_KeywordPresent(pp, "-wide"); if (o->gaussSeidel) { if (SOParams_KeywordPresent(pp, "-omega")) { o->omega = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); } else { o->omega = 1.0; } if (SOParams_KeywordPresent(pp, "-maxIter")) { o->maxIter = SOParams_GetNextInt(pp, 0, INT_MAX); } else { o->maxIter = 20; } if (SOParams_KeywordPresent(pp, "-absTol")) { o->absTol = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); } else { o->absTol = 0.0; } if (SOParams_KeywordPresent(pp, "-relTol")) { o->relTol = SOParams_GetNextDouble(pp, 0.0, MAXDOUBLE); } else { o->relTol = 0.0; } } SOParams_GetKeyword(pp, "-gaussOrder"); o->gaussOrder = SOParams_GetNextInt(pp, 1, INT_MAX); /* Plotting options: */ o->plot = SOParams_KeywordPresent(pp, "-plot"); o->plt = SOPlotParams_FunctionParse(pp); SOParams_Finish(pp); return o; } static PSStream *ps; static dg_Rank maxPlotDepth; void tin_plot2D( dg_Dim d, dg_Locus E, dg_Rank r, bz_Patch *b, dg_NodeStar *na ) { int m = b->m; /* Dimension of item. */ int n = b->n; /* Dimension of range (coords plus func values). */ fprintf(stderr, "+ tin_plot2D(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r); fprintf(stderr, "item shape\n"); bz_print(stderr, b, "%6.2f"); affirm(m <= 2, "bad domain dimension"); double fMax = 1.0; double fStep = 0.2; double fSync = fStep/2; int kMin = pswr_inf_isoline(fSync, fStep, -fMax); int kMax = pswr_sup_isoline(fSync, fStep, +fMax); int N; double *R, *G, *B; pswr_make_color_table( fSync, fStep, kMin, kMax, 0.000, 0.500, 1.000, 1.000, 1.000, 1.000, 1.000, 0.167, 0.000, &N, &R, &G, &B); switch(m) { case 0: { double *c = &(b->c[0]); dg_plot_2D_vertex(ps, c, n); } break; case 1: affirm(b->g[0] == 1, "bad patch degree"); { double *c0 = &(b->c[0]); double *c1 = &(b->c[n]); dg_plot_2D_edge(ps, c0, c1, n); } break; case 2: affirm((b->g[0] == 1) && (b->g[1] == 1), "bad patch degree"); { double *c00 = &(b->c[0]); double *c01 = &(b->c[n]); double *c10 = &(b->c[2*n]); double *c11 = &(b->c[3*n]); int subDp = maxPlotDepth - r + 1; if (subDp < 0) { subDp = 0; } if (subDp > 3) { subDp = 3; } dg_plot_2D_face ( ps, c00, c01, c10, c11, n, subDp, fSync, fStep, kMin, kMax, R,G,B ); } break; default: affirm(FALSE, "bad item dimension"); } fprintf(stderr, "- tin_plot2D\n"); } void tin_test_plot(char *name, bz_Degree g, bz_RDim n, dg_Rank treeDp, dg_Rank plotDp) { bool epsformat = FALSE; bz_Patch b = bz_make_test_patch(DDIM, n, g); dg_Node *grid_root = tin_make_grid(DDIM, treeDp, &b); dg_Locus E0 = dg_locus(0, DG_ROOT_CELL); dg_NodeStar na = dg_NodeRef_vec_new(1); dg_Interval bbox[MAX_RDIM]; affirm(b.m == 2, "wrong dimension"); fprintf(stderr, "grid cells = %d\n", grid_root->nodes); fprintf(stderr, "root shape =\n"); bz_print(stderr, &b, "%6.2f"); /* bz_Patch_dev_print(stderr, sh->dev, b.m, "%.2f"); */ /* fprintf(stderr, "\n"); */ maxPlotDepth = plotDp; bz_compute_bbox(&b, bbox); ps = tin_init_plot(name, epsformat, "a3", bbox); if (DEBUG_SHAPE) { tin_slow_plot(ps, &b); } else { na.el[0] = grid_root; dg_enum_locus_leaves(DDIM, E0, 0, &b, &na, plotDp, tin_plot2D); } tin_finish_plot(ps, epsformat); } PSStream *tin_init_plot ( char *name, bool epsformat, char *paperSize, dg_Interval bbox[] ) { double hpt, vpt; if (epsformat) { hpt = 5.0*72.0; vpt = hpt * 4.0/3.0; } else { pswr_get_paper_dimensions(paperSize, &hpt, &vpt); } PSStream *ps = pswr_new_stream(epsformat, name, NULL, paperSize, hpt+6, vpt+8); pswr_set_page_layout(ps, 3.0, 4.0, TRUE, 0.0, 0.0, (epsformat ? 0 : 1), 1,1); double dx = 0.1*(bbox[0].end[BHI] - bbox[0].end[BLO])/2; double dy = 0.1*(bbox[1].end[BHI] - bbox[1].end[BLO])/2; pswr_new_picture ( ps, bbox[0].end[BLO]-dx, bbox[0].end[BHI]+dx, bbox[1].end[BLO]-dy, bbox[1].end[BHI]+dy ); pswr_set_pen(ps, 0,0,0, 0.15, 0,0); if (! epsformat) { pswr_add_caption(ps, "Bezier-distorted grid", 0.0); } pswr_set_pen(ps, 1,0,0, 0.25, 0,0); return ps; } void tin_finish_plot ( PSStream *ps, bool epsformat ) { pswr_close_stream(ps); } #define x_well_1 (2.5) #define y_well_1 (4.5) #define x_well_2 (4.0) #define y_well_2 (5.0) #define y_fault (2.5) #define r_fault (0.5) #define w_fault (1.5) double tin_tol(double x, double y) { double d1 = tin_dsoft(x - x_well_1, y - y_well_1, 0.01); // double d1 = 100000.00; double d2 = tin_dsoft(x - x_well_2, y - y_well_2, 0.01); // double d2 = 100000.00; double yfx = y_fault + r_fault*sin(w_fault*x); double d3 = tin_dsoft(0.0, y - yfx, 0.01); // double d3 = 100000.00; double maxr = 0.25 * 3.0/(1.0/d1 + 1.0/d2 + 1.0/d3); return maxr; } double tin_inside(double x, double y, double rad) { double dctr = hypot(x - 0.5, y - 0.5); return dctr - rad < 0.4; } double tin_dsoft(double dx, double dy, double eps) { return sqrt(dx*dx + dy*dy + eps*eps); } dg_Node *tin_make_grid(dg_Dim d, dg_Rank treeDp, bz_Patch *b) { dg_Interval bbox[MAX_RDIM]; double ctr[MAX_RDIM]; int n = b->n; auto bool omit(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p); /* Tells {dg_make_grid} when to discard a cell. Basically, when it does not intersect a specified circle. */ auto bool split(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p); /* Tells {dg_make_grid} when to split a cell. Basically, when the variation in the coordinates and function values within the cell is larger than some tolerance, which depends on the cell's center. */ bool omit(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p) { // int i; double r2 = 0.0; // if (r >= treeDp) return FALSE; // bz_compute_bbox(b, bbox); // for (i = 0; i < d; i++) // { dg_Interval bi = bbox[i]; // double ri = (bi.hi - bi.lo)/2; // ctr[i] = (bi.lo + bi.hi)/2; // r2 += ri*ri; // } // return ! tin_inside(ctr[0], ctr[1], sqrt(r2)); fprintf(stderr, " omit k = %lld\n", k); return FALSE; } bool split(dg_CellIndex k, dg_Rank r, bz_Patch *b, dg_Node *p) { int i; double r2 = 0.0; double rad, radMax; if (r >= treeDp) return FALSE; bz_compute_bbox(b, bbox); for (i = 0; i < n; i++) { dg_Interval bi = bbox[i]; double ri = (bi.end[BHI] - bi.end[BLO])/2; /* Half-width of interval. */ ctr[i] = (bi.end[BLO] + bi.end[BHI])/2; r2 += ri*ri; } rad = sqrt(r2); radMax = tin_tol(ctr[0], ctr[1]); fprintf(stderr, " split k = %lld rad = %7.4f radMax = %7.4f\n", k, rad, radMax); return rad >= radMax; } return dg_make_grid(d, DG_ROOT_CELL, 0, b, NULL, &omit, &split); } void tin_slow_plot(PSStream *ps, bz_Patch *b) { int BSTEPS = 10; double R = 1.0, G = 0.9, B = 0.6; int f; fprintf(stderr, "+tin_slow_plot dim = %d\n", b->m); /* Plot faces */ for (f = 0; f < NF; f++) { if (dg_face_dimension(f, b->m) == 2) { bz_Patch *t = tin_do_extract_face(b, f); tin_slow_plot_face(ps, t, BSTEPS, R, G, B); free(t->c); free(t); } } /* Plot edges */ pswr_set_pen(ps, 0.0, 0.0, 0.0, 0.15, 0.0, 0.0); for (f = 0; f < NF; f++) { if (dg_face_dimension(f, b->m) == 1) { bz_Patch *t = tin_do_extract_face(b, f); tin_slow_plot_edge(ps, t, BSTEPS); free(t->c); free(t); } } /* Plot vertices */ pswr_set_pen(ps, 0.0, 0.0, 0.0, 0.15, 0.0, 0.0); for (f = 0; f < NF; f++) if (dg_face_dimension(f, b->m) == 0) { bz_Patch *t = tin_do_extract_face(b, f); tin_slow_plot_vertex(ps, t, BSTEPS); free(t->c); free(t); } fprintf(stderr, "-tin_slow_plot\n"); } void tin_slow_plot_vertex(PSStream *ps, bz_Patch *b, int steps) { fprintf(stderr, "+tin_slow_plot_vertex dim = %d\n", b->m); affirm(b->m == 0, "invalid vertex patch"); pswr_set_fill_color(ps, 0.0,0.0,0.0); pswr_dot(ps, b->c[0], b->c[1], 0.25, TRUE, FALSE); fprintf(stderr, "-tin_slow_plot_vertex\n"); } void tin_slow_plot_edge(PSStream *ps, bz_Patch *b, int steps) { fprintf(stderr, "+tin_slow_plot_edge dim = %d\n", b->m); affirm(b->m == 1, "invalid edge patch"); double p[2], q[2]; int i; fprintf(stderr, "shape =\n"); bz_print(stderr, b, "%6.2f"); for (i = 0; i <= steps; i++) { double ui = ((double)i)/((double)steps); bz_eval(b, &ui, p); if (i > 0) { pswr_segment(ps, q[0], q[1], p[0], p[1]); } q[0] = p[0]; q[1] = p[1]; } fprintf(stderr, "-tin_slow_plot_edge\n"); } void tin_slow_plot_face(PSStream *ps, bz_Patch *b, int steps, double R, double G, double B) { int i, j; fprintf(stderr, "+tin_slow_plot_face dim = %d\n", b->m); affirm(b->m == 2, "invalid face patch"); fprintf(stderr, "shape =\n"); bz_print(stderr, b, "%6.2f"); double u[2], v[2], p[MAX_RDIM], q[MAX_RDIM], r[MAX_RDIM], s[MAX_RDIM]; for (i = 1; i <= steps; i++) { for (j = 0; j <= steps; j++) { u[0] = ((double)i)/((double)steps); u[1] = ((double)j)/((double)steps); bz_eval(b, u, p); v[0] = ((double)i-1)/((double)steps); v[1] = ((double)j)/((double)steps); bz_eval(b, v, r); if (j > 0) { double xm = (p[0]+q[0]+r[0]+s[0])/4.0; double ym = (p[1]+q[1]+r[1]+s[1])/4.0; pswr_set_fill_color(ps, R,G,B); pswr_triangle(ps, p[0], p[1], q[0], q[1], xm, ym, TRUE, FALSE); pswr_triangle(ps, p[0], p[1], r[0], r[1], xm, ym, TRUE, FALSE); pswr_triangle(ps, s[0], s[1], q[0], q[1], xm, ym, TRUE, FALSE); pswr_triangle(ps, s[0], s[1], r[0], r[1], xm, ym, TRUE, FALSE); } q[0] = p[0]; q[1] = p[1]; s[0] = r[0]; s[1] = r[1]; } } fprintf(stderr, "-tin_slow_plot_face\n"); } bz_Patch *tin_do_extract_face(bz_Patch *b, dg_FaceIndex f) { dg_Dim bm = b->m, tm = 0; dg_SignedDir loc[DDIM]; dg_face_signature(f, bm, loc); bz_Patch *t = (bz_Patch *)notnull(malloc(sizeof(bz_Patch)), "no mem"); bz_Degree tg[DDIM]; int i; /* Extract degree sequence: */ for(i = 0; i < bm; i++) { if (loc[i] == SMD) { tg[tm] = b->g[i]; tm++; } } (*t) = bz_Patch_new(tm, b->n, tg); bz_get_face(b, loc, t); return t; } void tin_shatter_node ( SOGrid_Dim DDIM, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, double *maxhi, int minRank, int maxRank, double diff_limit, SOFunction *f, int *leaf_cells ); void tin_shatter_node ( SOGrid_Dim DDIM, SOGrid_Node *p, SOGrid_Rank r, double *lo, double *hi, double *maxhi, int minRank, int maxRank, double diff_limit, SOFunction *f, int *leaf_cells ) { bool split;//, border = FALSE; int i, j; double mid_val, maxdiff=0, diff, lo_child[DDIM], hi_child[DDIM]; //, abs_mid_val; double pnt_mid[DDIM], pnt_lo[DDIM], pnt_hi[DDIM], pnt_lo_val, pnt_hi_val; for(i = 0; i < DDIM; i++) pnt_mid[i] = (double) (hi[i] + lo[i]) / 2.0; f->m->eval(f, pnt_mid, &mid_val); for(j = 0; j < DDIM; j++) { for(i = 0; i < DDIM; i++) if(i != j){ pnt_lo[i] = pnt_mid[i]; pnt_hi[i] = pnt_mid[i]; } pnt_lo[j] = lo[j]; f->m->eval(f, pnt_lo, &pnt_lo_val); pnt_hi[j] = hi[j]; f->m->eval(f, pnt_hi, &pnt_hi_val); // printf(" ###### DEBUG-> pnt_lo[0]: %f, pnt_lo[1]: %f \n\n", pnt_lo[0], pnt_lo[1]); // printf(" ###### DEBUG-> pnt_hi[0]: %f, pnt_hi[1]: %f \n\n", pnt_hi[0], pnt_hi[1]); diff = fabs(mid_val - (pnt_lo_val + pnt_hi_val)/2.0); if(diff > maxdiff) maxdiff = diff; // if(lo[j] == 0.0 || hi[j] == maxhi[j]){ border = TRUE; abs_mid_val = fabs(mid_val); } } // printf(" ###### DEBUG maxdiff: %.16g diff_limit: %.16g \n\n", maxdiff, diff_limit); split = (maxdiff > diff_limit);// || (border && abs_mid_val > 2.0 * diff_limit)); if (r < minRank){split = TRUE;} if (r >= maxRank){split = FALSE;} if (split) { *leaf_cells = *leaf_cells + 1; if( p->c[LO] == NULL && p->c[HI] == NULL )SOGrid_Node_split(p); /* Compute low and high coordinates of low child. */ for (i = 0; i < DDIM; i++) { lo_child[i] = lo[i]; hi_child[i] = (i == (r%DDIM) ? (lo[i]+hi[i])/2 : hi[i]); } tin_shatter_node(DDIM, p->c[LO], r+1, lo_child, hi_child, maxhi, minRank, maxRank, diff_limit, f, leaf_cells); /* Compute low and high coordinates of high child. */ for (i = 0; i < DDIM; i++) { lo_child[i] = (i == (r%DDIM) ? (lo[i]+hi[i])/2 : lo[i]); hi_child[i] = hi[i]; } tin_shatter_node(DDIM, p->c[HI], r+1, lo_child, hi_child, maxhi, minRank, maxRank, diff_limit, f, leaf_cells); } } /* ====================================================================== */ void SetGenDotMatrices ( MatEntry_vec *evaldot, smt_matrix_t Meval, unsigned int offset, unsigned int maxindices ); void SetDotMatrices ( unsigned int dim, unsigned int offset, Basis bas, MatEntry_vec *evaldot, MatEntry_vec *Hev_ents ); void SetGenDotMatrices ( MatEntry_vec *evaldot, smt_matrix_t Meval, unsigned int offset, unsigned int maxindices ) { unsigned int i, index, indexi; for(i = 0; i < maxindices; i++) { evaldot[i] = MatEntry_vec_new(Meval.cols); evaldot[i].nel = 0; } i = 0; while(i < Meval.ents.nel) { index = Meval.ents.el[i].row; while(Meval.ents.el[i].row == index) { indexi = Meval.ents.el[i].row - offset; evaldot[indexi].el[evaldot[indexi].nel].row = Meval.ents.el[i].row; evaldot[indexi].el[evaldot[indexi].nel].col = Meval.ents.el[i].col; evaldot[indexi].el[evaldot[indexi].nel].va = Meval.ents.el[i].va; evaldot[indexi].nel++; i++; } } } void SetDotMatrices ( unsigned int dim, unsigned int offset, Basis bas, MatEntry_vec *evaldot, MatEntry_vec *Hev_ents ) { int i, j, k, hi_ind, low_ind; SOTentFunction *itf, *jtf; Hev_ents->nel = 0; for(i = 0; i < dim; i++) for(j = 0; j < dim; j++) { itf = (SOTentFunction *)bas.el[i]; jtf = (SOTentFunction *)bas.el[j]; if(itf->d->index > jtf->d->index){ hi_ind = itf->d->index; low_ind = jtf->d->index; } else { hi_ind = jtf->d->index; low_ind = itf->d->index; } k = 0; while(k < evaldot[hi_ind - offset].nel) { if(evaldot[hi_ind - offset].el[k].col == low_ind) { Hev_ents->el[Hev_ents->nel].row = i; Hev_ents->el[Hev_ents->nel].col = j; Hev_ents->el[Hev_ents->nel].va = evaldot[hi_ind - offset].el[k].va; Hev_ents->nel++; } k++; } } }