2D surface approximation with additional constraint for all components (x,y,z) combined.
(
As: COO,
bs: Array,
data: Array,
u: Array,
v: Array,
uorder: int,
vorder: int,
uknots: int | Array = 50,
vknots: int | Array = 50,
uperiodic: bool = False,
vperiodic: bool = False,
penalty: int = 3,
lam: Real = 0,
)
| 1877 | |
| 1878 | @multidispatch |
| 1879 | def constrainedApproximate2D( |
| 1880 | As: COO, |
| 1881 | bs: Array, |
| 1882 | data: Array, |
| 1883 | u: Array, |
| 1884 | v: Array, |
| 1885 | uorder: int, |
| 1886 | vorder: int, |
| 1887 | uknots: int | Array = 50, |
| 1888 | vknots: int | Array = 50, |
| 1889 | uperiodic: bool = False, |
| 1890 | vperiodic: bool = False, |
| 1891 | penalty: int = 3, |
| 1892 | lam: Real = 0, |
| 1893 | ) -> Surface: |
| 1894 | """ |
| 1895 | 2D surface approximation with additional constraint for all components (x,y,z) combined. |
| 1896 | """ |
| 1897 | |
| 1898 | # process the knots |
| 1899 | uknots_ = uknots if isinstance(uknots, Array) else np.linspace(0, 1, uknots) |
| 1900 | vknots_ = vknots if isinstance(vknots, Array) else np.linspace(0, 1, vknots) |
| 1901 | |
| 1902 | # create the design matrix |
| 1903 | C = designMatrix2D( |
| 1904 | u, v, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic |
| 1905 | ).csc() |
| 1906 | |
| 1907 | # handle penalties if requested |
| 1908 | if lam: |
| 1909 | # construct the penalty grid |
| 1910 | up, vp = uniformGrid(uknots_, vknots_, uorder, vorder, uperiodic, vperiodic) |
| 1911 | |
| 1912 | # construct the derivative matrices |
| 1913 | penalties = penaltyMatrix2D( |
| 1914 | up, vp, uorder, vorder, penalty, uknots_, vknots_, uperiodic, vperiodic, |
| 1915 | ) |
| 1916 | |
| 1917 | # augment the design matrix |
| 1918 | tmp = [comb(penalty, i) * penalties[i].csc() for i in range(penalty + 1)] |
| 1919 | Lu = uknots_[-1] - uknots_[0] # v length of the parametric domain |
| 1920 | Lv = vknots_[-1] - vknots_[0] # u length of the parametric domain |
| 1921 | P = Lu * Lv / len(up) * sp.vstack(tmp) |
| 1922 | |
| 1923 | CtC = C.T @ C + lam * P.T @ P |
| 1924 | else: |
| 1925 | CtC = C.T @ C |
| 1926 | |
| 1927 | # assemble one large block diagonal matrix |
| 1928 | CtC_xyz = sp.block_diag((CtC, CtC, CtC)) |
| 1929 | C_xyz = sp.block_diag((C, C, C)) |
| 1930 | |
| 1931 | # solve the augmented normal equations |
| 1932 | LHS = sp.bmat(((CtC_xyz, As.coo().T,), (As.coo(), None))) |
| 1933 | RHS = np.concatenate((C_xyz.T @ data.flatten(order="F"), bs)) |
| 1934 | |
| 1935 | LHS_DM = DM(LHS) |
| 1936 | ldl = Linsol("mumps", "mumps", LHS_DM.sparsity()) |