Description of grad3DCurv

0001 function G = grad3DCurv(k, X, Y, Z)
0002 % Returns a 3D 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 
0009     % Get the determinant of the jacobian and the metrics
0010     [J, Xe, Xn, Xc, Ye, Yn, Yc, Ze, Zn, Zc] = jacobian3D(k, X, Y, Z);
0011     
0012     % Dimensions of nodal grid
0013     [n, m, o] = size(X);
0014     
0015     % Make them volumes so they can be interpolated
0016     J = permute(reshape(J, m, n, o), [2, 1, 3]);
0017     A = permute(reshape(Yn.*Zc-Zn.*Yc, m, n, o), [2, 1, 3]);
0018     B = permute(reshape(Zn.*Xc-Xn.*Zc, m, n, o), [2, 1, 3]);
0019     C = permute(reshape(Xn.*Yc-Yn.*Xc, m, n, o), [2, 1, 3]);
0020     D = permute(reshape(Ze.*Yc-Ye.*Zc, m, n, o), [2, 1, 3]);
0021     E = permute(reshape(Xe.*Zc-Ze.*Xc, m, n, o), [2, 1, 3]);
0022     F = permute(reshape(Ye.*Xc-Xe.*Yc, m, n, o), [2, 1, 3]);
0023     G = permute(reshape(Ye.*Zn-Ze.*Yn, m, n, o), [2, 1, 3]);
0024     H = permute(reshape(Ze.*Xn-Xe.*Zn, m, n, o), [2, 1, 3]);
0025     I = permute(reshape(Xe.*Yn-Ye.*Xn, m, n, o), [2, 1, 3]);
0026     
0027     % Logical grids
0028     [Xl, Yl, Zl] = meshgrid(1:m, 1:n, 1:o);
0029     
0030     % Interpolate the metrics on the logical grid for positions u, v and w
0031     Xl_ = (Xl(1:end-1, :, :)+Xl(2:end, :, :))/2;
0032     Xl_ = (Xl_(:, :, 1:end-1)+Xl_(:, :, 2:end))/2;
0033     Yl_ = (Yl(1:end-1, :, :)+Yl(2:end, :, :))/2;
0034     Yl_ = (Yl_(:, :, 1:end-1)+Yl_(:, :, 2:end))/2;
0035     Zl_ = (Zl(1:end-1, :, :)+Zl(2:end, :, :))/2;
0036     Zl_ = (Zl_(:, :, 1:end-1)+Zl_(:, :, 2:end))/2;
0037     Ju = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_);
0038     A = interp3(Xl, Yl, Zl, A, Xl_, Yl_, Zl_);
0039     D = interp3(Xl, Yl, Zl, D, Xl_, Yl_, Zl_);
0040     G = interp3(Xl, Yl, Zl, G, Xl_, Yl_, Zl_);
0041     
0042     Xl_ = (Xl(:, 1:end-1, :)+Xl(:, 2:end, :))/2;
0043     Xl_ = (Xl_(:, :, 1:end-1)+Xl_(:, :, 2:end))/2;
0044     Yl_ = (Yl(:, 1:end-1, :)+Yl(:, 2:end, :))/2;
0045     Yl_ = (Yl_(:, :, 1:end-1)+Yl_(:, :, 2:end))/2;
0046     Zl_ = (Zl(:, 1:end-1, :)+Zl(:, 2:end, :))/2;
0047     Zl_ = (Zl_(:, :, 1:end-1)+Zl_(:, :, 2:end))/2;
0048     Jv = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_);
0049     B = interp3(Xl, Yl, Zl, B, Xl_, Yl_, Zl_);
0050     E = interp3(Xl, Yl, Zl, E, Xl_, Yl_, Zl_);
0051     H = interp3(Xl, Yl, Zl, H, Xl_, Yl_, Zl_);
0052     
0053     Xl_ = (Xl(1:end-1, :, :)+Xl(2:end, :, :))/2;
0054     Xl_ = (Xl_(:, 1:end-1, :)+Xl_(:, 2:end, :))/2;
0055     Yl_ = (Yl(1:end-1, :, :)+Yl(2:end, :, :))/2;
0056     Yl_ = (Yl_(:, 1:end-1, :)+Yl_(:, 2:end, :))/2;
0057     Zl_ = (Zl(1:end-1, :, :)+Zl(2:end, :, :))/2;
0058     Zl_ = (Zl_(:, 1:end-1, :)+Zl_(:, 2:end, :))/2;
0059     Jw = interp3(Xl, Yl, Zl, J, Xl_, Yl_, Zl_);
0060     C = interp3(Xl, Yl, Zl, C, Xl_, Yl_, Zl_);
0061     F = interp3(Xl, Yl, Zl, F, Xl_, Yl_, Zl_);
0062     I = interp3(Xl, Yl, Zl, I, Xl_, Yl_, Zl_);
0063     
0064     % Convert metrics to diagonal matrices so they can be multiplied by the
0065     % logical operators
0066     Ju = spdiags(1./reshape(permute(Ju, [2, 1, 3]), [], 1), 0, numel(Ju), numel(Ju));
0067     Jv = spdiags(1./reshape(permute(Jv, [2, 1, 3]), [], 1), 0, numel(Jv), numel(Jv));
0068     Jw = spdiags(1./reshape(permute(Jw, [2, 1, 3]), [], 1), 0, numel(Jw), numel(Jw));
0069     A = spdiags(reshape(permute(A, [2, 1, 3]), [], 1), 0, numel(A), numel(A));
0070     B = spdiags(reshape(permute(B, [2, 1, 3]), [], 1), 0, numel(B), numel(B));
0071     C = spdiags(reshape(permute(C, [2, 1, 3]), [], 1), 0, numel(C), numel(C));
0072     D = spdiags(reshape(permute(D, [2, 1, 3]), [], 1), 0, numel(D), numel(D));
0073     E = spdiags(reshape(permute(E, [2, 1, 3]), [], 1), 0, numel(E), numel(E));
0074     F = spdiags(reshape(permute(F, [2, 1, 3]), [], 1), 0, numel(F), numel(F));
0075     G = spdiags(reshape(permute(G, [2, 1, 3]), [], 1), 0, numel(G), numel(G));
0076     H = spdiags(reshape(permute(H, [2, 1, 3]), [], 1), 0, numel(H), numel(H));
0077     I = spdiags(reshape(permute(I, [2, 1, 3]), [], 1), 0, numel(I), numel(I));
0078     
0079     % Construct 3D uniform mimetic gradient operator (d/de, d/dn, d/dc)
0080     Grad = grad3D(k, m-1, 1, n-1, 1, o-1, 1);
0081     Ge = Grad(1:m*(n-1)*(o-1), :);
0082     Gn = Grad(m*(n-1)*(o-1)+1:m*(n-1)*(o-1)+(m-1)*n*(o-1), :);
0083     Gc = Grad(m*(n-1)*(o-1)+(m-1)*n*(o-1)+1:end, :);
0084     
0085     % Apply transformation
0086     Gx = Ju*(A*Ge+D*GI13(Gn, m-1, n-1, o-1, 'Gn')+G*GI13(Gc, m-1, n-1, o-1, 'Gc'));
0087     Gy = Jv*(B*GI13(Ge, m-1, n-1, o-1, 'Ge')+E*Gn+H*GI13(Gc, m-1, n-1, o-1, 'Gcy'));
0088     Gz = Jw*(C*GI13(Ge, m-1, n-1, o-1, 'Gee')+F*GI13(Gn, m-1, n-1, o-1, 'Gnn')+I*Gc);
0089     
0090     % Final 3D curvilinear mimetic gradient operator (d/dx, d/dy, d/dz)
0091     G = [Gx; Gy; Gz];
0092 end
grad3DCurv

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
