grad2DCurv

 Returns a 2D curvilinear mimetic gradient
 ----------------------------------------------------------------------------
 SPDX-License-Identifier: GPL-3.0-or-later
 © 2008-2024 San Diego State University Research Foundation (SDSURF).
 See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
 ----------------------------------------------------------------------------
 Get the determinant of the jacobian and the metrics

0001 function G = grad2DCurv(k, X, Y)
0002 % Returns a 2D curvilinear mimetic gradient
0003 % ----------------------------------------------------------------------------
0004 % SPDX-License-Identifier: GPL-3.0-or-later
0005 % © 2008-2024 San Diego State University Research Foundation (SDSURF).
0006 % See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
0007 % ----------------------------------------------------------------------------
0008     % Get the determinant of the jacobian and the metrics
0009     [J, Xe, Xn, Ye, Yn] = jacobian2D(k, X, Y);
0010     
0011     % Dimensions of nodal grid
0012     [n, m] = size(X);
0013     
0014     % Make them surfaces so they can be interpolated
0015     J = reshape(J, m, n)';
0016     Xe = reshape(Xe, m, n)';
0017     Xn = reshape(Xn, m, n)';
0018     Ye = reshape(Ye, m, n)';
0019     Yn = reshape(Yn, m, n)';
0020     
0021     % Logical grids
0022     [Xl, Yl] = meshgrid(1:m, 1:n);
0023     
0024     % Interpolate the metrics on the logical grid for positions u and v
0025     Ju = interp2(Xl, Yl, J, (Xl(1:end-1, :)+Xl(2:end, :))/2,...
0026                                             (Yl(1:end-1, :)+Yl(2:end, :))/2);
0027     Jv = interp2(Xl, Yl, J, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,...
0028                                             (Yl(:, 1:end-1)+Yl(:, 2:end))/2);
0029     Xev = interp2(Xl, Yl, Xe, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,...
0030                                             (Yl(:, 1:end-1)+Yl(:, 2:end))/2);
0031     Xnv = interp2(Xl, Yl, Xn, (Xl(:, 1:end-1)+Xl(:, 2:end))/2,...
0032                                             (Yl(:, 1:end-1)+Yl(:, 2:end))/2);
0033     Yeu = interp2(Xl, Yl, Ye, (Xl(1:end-1, :)+Xl(2:end, :))/2,...
0034                                             (Yl(1:end-1, :)+Yl(2:end, :))/2);
0035     Ynu = interp2(Xl, Yl, Yn, (Xl(1:end-1, :)+Xl(2:end, :))/2,...
0036                                             (Yl(1:end-1, :)+Yl(2:end, :))/2);
0037     
0038     % Convert metrics to diagonal matrices so they can be multiplied by the
0039     % logical operators
0040     Ju = spdiags(1./reshape(Ju', [], 1), 0, numel(Ju), numel(Ju));
0041     Jv = spdiags(1./reshape(Jv', [], 1), 0, numel(Jv), numel(Jv));
0042     Xev = spdiags(reshape(Xev', [], 1), 0, numel(Xev), numel(Xev));
0043     Xnv = spdiags(reshape(Xnv', [], 1), 0, numel(Xnv), numel(Xnv));
0044     Yeu = spdiags(reshape(Yeu', [], 1), 0, numel(Yeu), numel(Yeu));
0045     Ynu = spdiags(reshape(Ynu', [], 1), 0, numel(Ynu), numel(Ynu));
0046     
0047     % Construct 2D uniform mimetic gradient operator (d/de, d/dn)
0048     G = grad2D(k, m-1, 1, n-1, 1);
0049     Ge = G(1:m*(n-1), :);
0050     Gn = G(m*(n-1)+1:end, :);
0051     
0052     % Apply transformation
0053     Gx = Ju*(Ynu*Ge-Yeu*GI2(Gn, m-1, n-1, 'Gn'));
0054     Gy = Jv*(-Xnv*GI2(Ge, m-1, n-1, 'Ge')+Xev*Gn);
0055     
0056     % Final 2D curvilinear mimetic gradient operator (d/dx, d/dy)
0057     G = [Gx; Gy];
0058 end