B ]t\@sdZddlZddlZddlmZddlmZddlm Z dZ Gddde Z Gd d d e Z Gd d d e ZGd ddZGdddZGdddeZGdddeZGdddeZGddde Zd-ddZddZdd Zd!d"Zd#d$Zd%d&Zd'd(Zd)d*Zd+d,ZdS).z( Interpolation inside triangular grids. N) Triangulation) TriFinder) TriAnalyzer)TriInterpolatorLinearTriInterpolatorCubicTriInterpolatorc@s4eZdZdZd ddZdZdZd dd Zd d ZdS)ra Abstract base class for classes used to perform interpolation on triangular grids. Derived classes implement the following methods: - ``__call__(x, y)`` , where x, y are array_like point coordinates of the same shape, and that returns a masked array of the same shape containing the interpolated z-values. - ``gradient(x, y)`` , where x, y are array_like point coordinates of the same shape, and that returns a list of 2 masked arrays of the same shape containing the 2 derivatives of the interpolator (derivatives of interpolated z values with respect to x and y). NcCs~t|tstd||_t||_|jj|jjjkr>td|dk rXt|t sXtd|pd|j |_ d|_ d|_ d|_dS)NzExpected a Triangulation objectz=z array must have same length as triangulation x and y arrayszExpected a TriFinder objectg?) isinstancer ValueError_triangulationnpasarray_zshapexrZ get_trifinder _trifinder_unit_x_unit_y _tri_renum)self triangulationz trifinderrZ tri_analyzerZcompressed_trianglesZ compressed_xZ compressed_yZ tri_renumZ node_renumZ node_maskrrrrs"    zCubicTriInterpolator.__init__cCs|j||ddddS)N)r)r-r.r)r2)rrr,rrrr;s zCubicTriInterpolator.__call__cCs|j||dddS)N)r r!)r-r.)r2)rrr,rrrr<s zCubicTriInterpolator.gradientc Cs|j|}||||}|j|}tj|j|dd}|dkrN|j|||S|dkr||} |j || ||} |dkr| ddddfS| ddddfSn t d|dS)Nr)axisr)r r!r rzInvalid return_key: ) rA_get_alpha_vecrCr expand_dimsrErGget_function_values _get_jacobianget_function_derivativesr ) rr0r-rr,tris_ptsalphaeccdofJr rrrr*s     z,CubicTriInterpolator._interpolate_single_keycCs`|dkr&|dkrtdt||d}n2|dkr8t|}n |dkrJt|}ntd||S)ap Computes and returns nodal dofs according to kind Parameters ---------- kind: {'min_E', 'geom', 'user'} Choice of the _DOF_estimator subclass to perform the gradient estimation. dz: tuple of array_likes (dzdx, dzdy), optional Used only if *kind*=user; in this case passed to the :class:`_DOF_estimator_user`. Returns ------- dof : array_like, shape (npts,2) Estimation of the gradient at triangulation nodes (stored as degree of freedoms of reduced-HCT triangle elements). userNzQFor a CubicTriInterpolator with *kind*='user', a valid *dz* argument is expected.)r>Zgeomr=zTCubicTriInterpolator *kind* proposed: {0}; should be one of: 'user', 'geom', 'min_E')r _DOF_estimator_user_DOF_estimator_geom_DOF_estimator_min_Er$compute_dof_from_df)rrHr>ZTErrrrDs  z!CubicTriInterpolator._compute_dofcCs.|jd}|dddddf|dddddf}|dddddf|dddddf}tj||gdd}t|}tj||gdd|dddddf}t||} t| } t|tt||} t| | } td| ddddf| ddddfg| ddddfg| ddddfgg} | S)aT Fast (vectorized) function to compute barycentric coordinates alpha. Parameters ---------- x, y : array-like of dim 1 (shape (nx,)) Coordinates of the points whose points barycentric coordinates are requested tris_pts : array like of dim 3 (shape: (nx,3,2)) Coordinates of the containing triangles apexes. Returns ------- alpha : array of dim 2 (shape (nx,3)) Barycentric coordinates of the points inside the containing triangles. rNrrr)rI)ndimr stack_transpose_vectorized_prod_vectorized_pseudo_inv22sym_vectorizedrK_to_matrix_vectorized)rr,rOrYabZabTZabZOMZmetricZ metric_invZCovarksirPrrrrJs ,,(  Rz#CubicTriInterpolator._get_alpha_veccCst|dddddf|dddddf}t|dddddf|dddddf}t|dddf|dddfg|dddf|dddfgg}|S)a Fast (vectorized) function to compute triangle jacobian matrix. Parameters ---------- tris_pts : array like of dim 3 (shape: (nx,3,2)) Coordinates of the containing triangles apexes. Returns ------- J : array of dim 3 (shape (nx,2,2)) Barycentric coordinates of the points inside the containing triangles. J[itri,:,:] is the jacobian matrix at apex 0 of the triangle itri, so that the following (matrix) relationship holds: [dz/dksi] = [J] x [dz/dx] with x: global coordinates ksi: element parametric coordinates in triangle first apex local basis. Nrrr)r arrayr^)rOr_r`rSrrrrMs 22 $z"CubicTriInterpolator._get_jacobiancCs"tj|dddddf|dddddfdd}tj|dddddf|dddddfdd}tj|dddddf|dddddfdd}tt||ddddf}tt||ddddf}tt||ddddf}t|||g|||g|||ggS)a Computes triangle eccentricities Parameters ---------- tris_pts : array like of dim 3 (shape: (nx,3,2)) Coordinates of the triangles apexes. Returns ------- ecc : array like of dim 2 (shape: (nx,3)) The so-called eccentricity parameters [1] needed for HCT triangular element. Nrr)rIr)r rKr\r[r^)rOr_r`cZdot_aZdot_bZdot_crrrrB1s666 z0CubicTriInterpolator._compute_tri_eccentricities)r=NN)N)r4r5r6r7rr;rr8r<r9r*rD staticmethodrJrMrBrrrrr%s^ ( # ( rc@seZdZdZeddddddddddg ddddddddddg ddddddddddg dddddd d ddd g ddddd dddddg d ddd ddddddg dddddddd d d g d dddd dddddg dddd ddddddg g Zeddddddddddg ddddddddddg ddddddddddg d ddd d dddddg d dddddddddg ddddddddddg dddddddddd g ddddddddddg dddddddddd g g Zed ddd ddddddg ddddddddddg ddddddddddg ddddddddddg ddddddddddg ddddddddddg ddddddddddg d dddddddddg ddddddddddg g Zeddddddddddg d dddddddddg ddddddddddg d dddd dddddg ddddddddddg ddddddddddg ddddddddddg ddddddddddg ddddddddddg g Z eddgddgddgd d gd d gddggZ edddgdddgdddgdddgdddgddd gdddgdddgddd gg Z dZ ejdddgdddgdddgdddgdddgdddgdddgdddgdddgg ej dZejdgej ddZedddgdddgdddggZed dgd dggZedd gdd ggZddZddZddZd d!Zd"d#Zd*d%d&Zd'd(Zd)S)+rFaN Implementation of reduced HCT triangular element with explicit shape functions. Computes z, dz, d2z and the element stiffness matrix for bending energy: E(f) = integral( (d2z/dx2 + d2z/dy2)**2 dA) *** Reference for the shape functions: *** [1] Basis functions for general Hsieh-Clough-Tocher _triangles, complete or reduced. Michel Bernadou, Kamal Hassan International Journal for Numerical Methods in Engineering. 17(5):784 - 789. 2.01 *** Element description: *** 9 dofs: z and dz given at 3 apex C1 (conform) gg@gпg?g?g?gg@g?ggg?gg?gg gqq?gqq?gqq?g98?)rg"@g@c Cstj|dddddf}t|| dd}t|| dd}|ddddf}|ddddf}|ddddf} ||} ||} | | } t| |g| |g| | g| | g| |g| |g| | g| |g| |g||| gg } t|j| }|t|ddddft|j| 7}|t|ddddft|j| 7}|t|ddddft|j | 7}t|d|dd}t||ddddfS)a Parameters ---------- alpha : is a (N x 3 x 1) array (array of column-matrices) of barycentric coordinates, ecc : is a (N x 3 x 1) array (array of column-matrices) of triangle eccentricities, dofs : is a (N x 1 x 9) arrays (arrays of row-matrices) of computed degrees of freedom. Returns ------- Returns the N-array of interpolated function values. r)rINrr) r argmin_roll_vectorizedr^r\M_scalar_vectorizedM0M1M2)rrPrQdofssubtriraErr,rx_sqy_sqz_sqVprodsrrrrLs*0* z'_ReducedHCT_Element.get_function_valuescCstj|dddddf}t|| dd}t|| dd}|ddddf}|ddddf} |ddddf} ||} | | } | | } td| d| gd| dgdd| gd || d || | gd || | d || gd || | | gd | | | g| d | | g| d || | g|| | | || | | gg }t|t|j|ddd }t|j|}|t|ddddft|j |7}|t|ddddft|j |7}|t|ddddft|j |7}t|d |dd}t||}t |}t|t |}|S) a Parameters ---------- *alpha* is a (N x 3 x 1) array (array of column-matrices of barycentric coordinates) *J* is a (N x 2 x 2) array of jacobian matrices (jacobian matrix at triangle first apex) *ecc* is a (N x 3 x 1) array (array of column-matrices of triangle eccentricities) *dofs* is a (N x 1 x 9) arrays (arrays of row-matrices) of computed degrees of freedom. Returns ------- Returns the values of interpolated function derivatives [dz/dx, dz/dy] in global coordinates at locations alpha, as a column-matrices of shape (N x 2 x 1). r)rINrrgg@ggg@) block_sizerIrf)r rgrhr^r\_extract_submatrices rotate_dVrirjrkrlrm_safe_inv22_vectorizedr[)rrPrSrQrnrorarprr,rrqrrrsZdVruZdsdksidfdksiJ_invZdfdxrrrrNsD  $  z,_ReducedHCT_Element.get_function_derivativesc Cs2|||}t||}||}t||}t|S)a Parameters ---------- *alpha* is a (N x 3 x 1) array (array of column-matrices) of barycentric coordinates *J* is a (N x 2 x 2) array of jacobian matrices (jacobian matrix at triangle first apex) *ecc* is a (N x 3 x 1) array (array of column-matrices) of triangle eccentricities *dofs* is a (N x 1 x 9) arrays (arrays of row-matrices) of computed degrees of freedom. Returns ------- Returns the values of interpolated function 2nd-derivatives [d2z/dx2, d2z/dy2, d2z/dxdy] in global coordinates at locations alpha, as a column-matrices of shape (N x 3 x 1). )get_d2Sidksij2r\get_Hrot_from_Jr[) rrPrSrQrnd2sdksi2Zd2fdksi2H_rotZd2fdx2rrrget_function_hessianss     z)_ReducedHCT_Element.get_function_hessiansc Cstj|dddddf}t|| dd}t|| dd}|ddddf}|ddddf}|ddddf}td|d|d|gd|ddgdd|dgd|d|d |d|d|gd|d |d|d|d|gd|d |dd |gd|dd|gdd|d|gdd|d |d |gd |d ||||gg } t| t|j|d dd } t|j| } | t|ddddft|j | 7} | t|ddddft|j | 7} | t|ddddft|j | 7} t| d |dd} | S) a Parameters ---------- *alpha* is a (N x 3 x 1) array (array of column-matrices) of barycentric coordinates *ecc* is a (N x 3 x 1) array (array of column-matrices) of triangle eccentricities Returns ------- Returns the arrays d2sdksi2 (N x 3 x 1) Hessian of shape functions expressed in covariante coordinates in first apex basis. r)rINrrg@gg@g@grf)rwrI) r rgrhr^r\rx rotate_d2Vrirjrkrlrm) rrPrQrorarprr,rZd2Vrurrrrr}+s8  $$ z"_ReducedHCT_Element.get_d2Sidksij2cCst|d}t|j|}t|j|}tj|ddgtjd}d|ddddf<d|ddddf<d|ddddf<||ddddddf<||ddddddf<||ddd dd df<|j|d d \}}tj|ddgtjd} |j} |j } xzt |j D]l} t | | ddf| |d} t| d } | | }|| |}t||}| |tt||jt|7} qWttt|| |} t|| S) a Parameters ---------- *J* is a (N x 2 x 2) array of jacobian matrices (jacobian matrix at triangle first apex) *ecc* is a (N x 3 x 1) array (array of column-matrices) of triangle eccentricities Returns ------- Returns the element K matrices for bending energy expressed in GLOBAL nodal coordinates. K_ij = integral [ (d2zi/dx2 + d2zi/dy2) * (d2zj/dx2 + d2zj/dy2) dA] tri_J is needed to rotate dofs from local basis to global basis rre)rrNrfT) return_arear)r r&r\J0_to_J1J0_to_J2zerosr#r~gauss_w gauss_ptsrangen_gaussZtiler+rKr}rpr[rj)rrSrQnJ1J2ZDOF_rotrareaKweightsZptsZigaussrPZweightZ d2Skdksi2Zd2Skdx2rrrget_bending_matricesWs2         z(_ReducedHCT_Element.get_bending_matricesFc Cst|}|ddddf}|ddddf}|ddddf}|ddddf}t||||||g||||||gd||d||||||gg}|s|Sd|ddddf|ddddf|ddddf|ddddf} || fSdS)as Parameters ---------- *J* is a (N x 2 x 2) array of jacobian matrices (jacobian matrix at triangle first apex) Returns ------- Returns H_rot used to rotate Hessian from local basis of first apex, to global coordinates. if *return_area* is True, returns also the triangle area (0.5*det(J)) Nrrrg?)rzr^) rrSrr|ZJi00ZJi11ZJi10ZJi01rrrrrr~s*Lz#_ReducedHCT_Element.get_Hrot_from_Jc Cst|d}tj|tjd}tj|tjd }ddddddg}dd d g} t||d d dfd|d d dfdd||d d dfd|d d dfdd||d d dfd|d d dfddg g} tj|d dgtjd} tt| t| } t| | } |||}t |t |||}t | t |||}t | t |||}|t ||| }tj |dd }t||d d d d df }| t |dg|d d dd d f}tj t |t |d}||||fS)aG Builds K and F for the following elliptic formulation: minimization of curvature energy with value of function at node imposed and derivatives 'free'. Builds the global Kff matrix in cco format. Builds the full Ff vec Ff = - Kfc x Uc Parameters ---------- *J* is a (N x 2 x 2) array of jacobian matrices (jacobian matrix at triangle first apex) *ecc* is a (N x 3 x 1) array (array of column-matrices) of triangle eccentricities *triangles* is a (N x 3) array of nodes indexes. *Uc* is (N x 3) array of imposed displacements at nodes Returns ------- (Kff_rows, Kff_cols, Kff_vals) Kff matrix in coo format - Duplicate (row, col) entries must be summed. Ff: force vector - dim npts * 3 r)rrrrrrfrNre)rI)r) r r&arangeint32onesr^r\r[rr%Zix_rKbincount)rrSrQ trianglesUcZntriZ vec_rangeZ c_indicesZf_dofZc_dofZ f_dof_indicesZexpand_indicesZ f_row_indicesZ f_col_indicesZK_elemKff_valsKff_rowsKff_colsZKfc_elemZUc_elemZFf_elemZ Ff_indicesFfrrrget_Kff_and_Ffs.  **2   &z"_ReducedHCT_Element.get_Kff_and_FfN)F)r4r5r6r7r rbrirkrlrmryrrr#rrrrprrrLrNrr}rr~rrrrrrFRs  "%;,5 rFc@s4eZdZdZddZddZddZedd Zd S) _DOF_estimatoraj Abstract base class for classes used to perform estimation of a function first derivatives, and deduce the dofs for a CubicTriInterpolator using a reduced HCT element formulation. Derived classes implement compute_df(self,**kwargs), returning np.vstack([dfx,dfy]).T where : dfx, dfy are the estimation of the 2 gradient coordinates. cKs^t|tstd|j|_|j|_|j|_|j|_|j|j |_|_ |j f||_ | dS)Nz&Expected a CubicTriInterpolator object) rrr r@rAr rr?rr compute_dzr>rX)rZ interpolatorkwargsrrrrs z_DOF_estimator.__init__cKstdS)N)r3)rrrrrr sz_DOF_estimator.compute_dzcCs6t|j}|j|j}|j|j}||||}|S)za Computes reduced-HCT elements degrees of freedom, knowing the gradient. )rrMrArr?r> get_dof_vec)rrStri_ztri_dzZtri_dofrrrrX s    z"_DOF_estimator.compute_dof_from_dfc Csr|jd}tj|dgtjd}ttj|}ttj|}t|tj|dddddfdd}t|tj|dddddfdd}t|tj|dddddfdd} t |ddddf|ddddf| ddddfg|ddddf|ddddf| ddddfgg} ||ddddd f<| dddf|dddd d f<| dddf|ddddd f<|S) a Computes the dof vector of a triangle, knowing the value of f, df and of the local Jacobian at each node. *tri_z*: array of shape (3,) of f nodal values *tri_dz*: array of shape (3,2) of df/dx, df/dy nodal values *J*: Jacobian matrix in local basis of apex 0 Returns dof array of shape (9,) so that for each apex iapex: dof[iapex*3+0] = f(Ai) dof[iapex*3+1] = df(Ai).(AiAi+) dof[iapex*3+2] = df(Ai).(AiAi-)] rre)rNr)rIrrrfr) rr rr#r\rFrrrKr^) rrrSZnptrRrrZcol0Zcol1Zcol2r{rrrrs   &&&28""z_DOF_estimator.get_dof_vecN) r4r5r6r7rrrXrdrrrrrrs   rc@seZdZdZddZdS)rUz5 dz is imposed by user / Accounts for scaling if any cCs,|\}}||j}||j}t||gjS)N)rrr vstackT)rr>r r!rrrr:s  z_DOF_estimator_user.compute_dzN)r4r5r6r7rrrrrrU8srUc@s(eZdZdZddZddZddZdS) rVz? Fast 'geometric' approximation, recommended for large arrays. c Cs|}|}tjt|jt|d}t|}t|}xhtdD]\}|dd|f|dddf|dd|f<|dd|f|dddf|dd|f<qJWtjt|jt|d}tjt|jt|d}||} ||} t| | gj S)a self.df is computed as weighted average of _triangles sharing a common node. On each triangle itri f is first assumed linear (= ~f), which allows to compute d~f[itri] Then the following approximation of df nodal values is then proposed: f[ipt] = SUM ( w[itri] x d~f[itri] , for itri sharing apex ipt) The weighted coeff. w[itri] are proportional to the angle of the triangle itri at apex ipt )rrfNrr) compute_geom_weightscompute_geom_gradsr rr%r? empty_likerrr) rZ el_geom_wZ el_geom_gradZ w_node_sumZdfx_el_wZdfy_el_wZiapexZ dfx_node_sumZ dfy_node_sumZ dfx_estimZ dfy_estimrrrrCs   ,0z_DOF_estimator_geom.compute_dzc CsPtt|jddg}|j}x*tdD]}|dd|dddf}|dd|ddddf}|dd|ddddf}t|dddf|dddf|dddf|dddf}t|dddf|dddf|dddf|dddf}tt||tj d} dt| d|dd|f<q*W|S)z Builds the (nelems x 3) weights coeffs of _triangles angles, renormalized so that np.sum(weights, axis=1) == np.ones(nelems) rrfNrg?g?) r rr&r?rArZarctan2absmodZpi) rrrOZiptZp0Zp1Zp2Zalpha1Zalpha2ZanglerrrrdsDD"z(_DOF_estimator_geom.compute_geom_weightsc Cs~|j}|j|j}|dddddf|dddddf}|dddddf|dddddf}t||g}t|}|dddf|dddf}|dddf|dddf}t||gj} t| } | dddf|ddddf| dddf|ddddf| dddf<| dddf|ddddf| dddf|ddddf| dddf<| S)z Compute the (global) gradient component of f assumed linear (~f). returns array df of shape (nelems,2) df[ielem].dM[ielem] = dz[ielem] i.e. df = dz x dM = dM.T^-1 x dz Nrrr) rArr?r Zdstackrzrrr) rrOZtris_fZdM1ZdM2ZdMZdM_invZdZ1ZdZ2ZdZZdfrrrrzs ,,   PPz&_DOF_estimator_geom.compute_geom_gradsN)r4r5r6r7rrrrrrrrVAs!rVc@s eZdZdZddZddZdS)rWz The 'smoothest' approximation, df is computed through global minimization of the bending energy: E(f) = integral[(d2z/dx2 + d2z/dy2 + 2 d2z/dxdy)**2 dA] cCs|j|_t||dS)N)rCrVr)rZ Interpolatorrrrrsz_DOF_estimator_min_E.__init__cCst|}t|}t}t|j}|j}|j }|j |j }| ||||\}} } } d} | j d} t | || | | fd}|t|| || d\}}tj||| }||krtd|}tj|jj ddgtjd}|ddd|dddf<|d dd|ddd f<|S) z{ Elliptic solver for bending energy minimization. Uses a dedicated 'toy' sparse Jacobi PCG solver. g|=r)r)Ar`x0tolzIn TriCubicInterpolator initialization, PCG sparse solver did not converge after 1000 iterations. `geom` approximation is used instead of `min_E`r)rNr)rVrr r%rFrrMrArCr?rrr_Sparse_Matrix_coo compress_csc_cglinalgnormdotwarningswarnr(r@r#)rZdz_initZUf0Zreference_elementrSZeccsrrrrrrrZn_dofZKff_cooZUferrZerr0r>rrrrs.       z_DOF_estimator_min_E.compute_dzN)r4r5r6r7rrrrrrrWsrWc@sHeZdZddZddZddZddZd d Zd d Ze d dZ dS)rcCsF|\|_|_tj|tjd|_tj|tjd|_tj|tjd|_dS)a4 Creates a sparse matrix in coo format *vals*: arrays of values of non-null entries of the matrix *rows*: int arrays of rows of non-null entries of the matrix *cols*: int arrays of cols of non-null entries of the matrix *shape*: 2-tuple (n,m) of matrix shape )rN) rmr r r#valsrrowscols)rrrrrrrrrs z_Sparse_Matrix_coo.__init__cCs2|j|jfksttj|j|j||j|jdS)z Dot product of self by a vector *V* in sparse-dense to dense format *V* dense vector of shape (self.m,) )rZ minlength)rrAssertionErrorr rrrr)rrtrrrrsz_Sparse_Matrix_coo.dotcCsRtj|j|j|jddd\}}}|j||_|j||_tj||jd|_dS)zV Compress rows, cols, vals / summing duplicates. Sort for csc format. T)r1return_inverse)rN)r uniquerrrrr)r_rindicesrrrrs   z_Sparse_Matrix_coo.compress_csccCsRtj|j|j|jddd\}}}|j||_|j||_tj||jd|_dS)zV Compress rows, cols, vals / summing duplicates. Sort for csr format. T)r1r)rN)r rrrrrr)rrrrrrr compress_csrs   z_Sparse_Matrix_coo.compress_csrcCs\tj|j|jgtjd}|jj}x6t|D]*}||j||j |f|j|7<q*W|S)zb Returns a dense matrix representing self. Mainly for debugging purposes. )r) r rrrr#rr&rrr)rr/Znvalsirrrto_denses *z_Sparse_Matrix_coo.to_densecCs |S)N)r__str__)rrrrr sz_Sparse_Matrix_coo.__str__cCs>|j|jk}tjt|j|jtjd}|j|||j|<|S)zF Returns the (dense) vector of the diagonal elements. )r)rrr rminrr#r)rZin_diagdiagrrrr s z_Sparse_Matrix_coo.diagN) r4r5r6rrrrrrpropertyrrrrrrs    r绽|=cCs>|j}|j|kst|j|ks"ttj|}|j}t|dk|d}|dkrZt |}n|}|| |} | |} t |} d} t | | } d}xt t | ||kr||kr| | | } | | }| t | |}| ||} | |} | }t | | } ||| }| |} |d7}qWtj| ||}||fS)ah Use Preconditioned Conjugate Gradient iteration to solve A x = b A simple Jacobi (diagonal) preconditionner is used. Parameters ---------- A: _Sparse_Matrix_coo *A* must have been compressed before by compress_csc or compress_csr method. b: array Right hand side of the linear system. Returns ------- x: array. The converged solution. err: float The absolute error np.linalg.norm(A.dot(x) - b) Other parameters ---------------- x0: array. Starting guess for the solution. tol: float. Tolerance to achieve. The algorithm terminates when the relative residual is below tol. maxiter: integer. Maximum number of iterations. Iteration will stop after maxiter steps even if the specified tolerance has not been achieved. gư>Ngrr) r&rrrr rrrwhererrZsqrtr)rr`rrmaxiterrZb_normZkvecrrwpZbetaZrhokrrPZrhooldrrrrrs8!    $      rcCs^|jdkst|jdddks$tt|}|ddddf|ddddf}||ddddf|ddddf}t|dt|k}t|rd|}n t|jd}d||||<|ddddf||ddddf<|ddddf ||ddddf<|ddddf ||ddddf<|ddddf||ddddf<|S) z Inversion of arrays of (2,2) matrices, returns 0 for rank-deficient matrices. *M* : array of (2,2) matrices to inverse, shape (n,2,2) rfN)rrrrg:0yE>g?)rYrrr rrallr)riM_invprod1deltarank2Z delta_invrrrrzs $(  $&&$rzc Csj|jdkst|jdddks$tt|}|ddddf|ddddf}||ddddf|ddddf}t|dt|k}t|r6|ddddf||ddddf<|ddddf ||ddddf<|ddddf ||ddddf<|ddddf||ddddf<n0||}||ddf|||ddf<||ddf |||ddf<||ddf |||ddf<||ddf|||ddf<|}||ddf||ddf}t|dk}d||d |}||ddf|||ddf<||ddf|||ddf<||ddf|||ddf<||ddf|||ddf<|S) a Inversion of arrays of (2,2) SYMMETRIC matrices; returns the (Moore-Penrose) pseudo-inverse for rank-deficient matrices. In case M is of rank 1, we have M = trace(M) x P where P is the orthogonal projection on Im(M), and we return trace(M)^-1 x P == M / trace(M)**2 In case M is of rank 0, we return the null matrix. *M* : array of (2,2) matrices to inverse, shape (n,2,2) rfrN)rrrrg:0yE>g?r)rYrrr rrr) rirrrrZrank01ZtrZtr_zerosZ sq_tr_invrrrr]s2  $( $&&(r]cCs|j}|j}t|dkstt|dks,t|d|dks@tt|}t|d|d|df}tt||dtjf|dtjddfdS)zm Matrix product between arrays of matrices, or a matrix and an array of matrices (*M1* and *M2*) rrrr.N)rlenrrr sum transposenewaxis)rlrmZsh1Zsh2Zndim1Zt1_indexrrrr\sr\cCs|ddtjtjf|S)z6 Scalar product between scalars and matrices. N)r r)ZscalarrirrrrjsrjcCst|dddgS)z4 Transposition of an array of matrices *M*. rrr)r r)rirrrr[sr[c Cs&|dks t|j}|dkst|j}|dks0t|j}|dd\}}|d|jdks\ttj|dtjd}t|} |dkrxt|D]<} x6t|D]*} ||| | || f| dd| | f<qWqWnT|dkr"xHt|D]<} x6t|D]*} ||| | | |f| dd| | f<qWqW| S)z Rolls an array of matrices along an axis according to an array of indices *roll_indices* *axis* can be either 0 (rolls rows) or 1 (rolls columns). )rrrfrrNr)r)rrYrr rrrr) riZ roll_indicesrIrYZ ndim_rollshrrcZ vec_indicesZM_rolliricrrrrhs&    0 .rhc Cst|ttfsttdd|Ds(ttdd|D}t||ddksVtt|}|d}t|dd}|j}|j d||g}tj ||d}xBt |D]6}x0t |D]$} t||| |dd|| f<qWqW|S)z Builds an array of matrices from individuals np.arrays of identical shapes. *M*: ncols-list of nrows-lists of shape sh. Returns M_res np.array of shape (sh, nrow, ncols) so that: M_res[...,i,j] = M[i][j] css|]}t|ttfVqdS)N)rtuplelist).0itemrrr $sz(_to_matrix_vectorized..cSsg|] }t|qSr)r)rrrrr %sz)_to_matrix_vectorized..r)rN) rrrrrr r rrrr(r) riZc_vecrrcZM00dtrZM_retZirowZicolrrrr^s (r^c Cs|jdkst|dkst|j\}}|dkr>|jd||g}n|dkrV|jd||g}|j}tj||d}|dkrxt|D].} |||| ddf|dd| ddf<q|WnD|dkrx:t|D].} |dd||| f|dddd| f<qW|S)z Extracts selected blocks of a matrices *M* depending on parameters *block_indices* and *block_size*. Returns the array of extracted matrices *Mres* so that: M_res[...,ir,:] = M[(block_indices*block_size+ir), :] r)rrr)rN)rYrrrr r(r) riZ block_indicesrwrIrrcrrZM_resrrrrrrx3s   0.rx)Nrr)r7rZnumpyr Zmatplotlib.trirZmatplotlib.tri.trifinderrZmatplotlib.tri.tritoolsr__all__objectrrrrFrrUrVrWrrrzr]r\rjr[rhr^rxrrrrs<   ^8/%D T:H y/