forked from FESTUNG/FESTUNG
-
Notifications
You must be signed in to change notification settings - Fork 0
/
applySlopeLimiterTaylorStrict.m
144 lines (136 loc) · 6.88 KB
/
applySlopeLimiterTaylorStrict.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
% Applies the stricter hierarchical vertex based slope limiter to a discrete
% function given in Taylor basis.
%===============================================================================
%> @file applySlopeLimiterTaylorStrict.m
%>
%> @brief Applies the stricter hierarchical vertex based slope limiter to a
%> discrete function given in Taylor basis.
%===============================================================================
%>
%> @brief Applies the stricter hierarchical vertex based slope limiter to a
%> discrete function given in Taylor basis.
%>
%> This is a modification of the hierarchical vertex based limiter (cf.
%> <code>applySlopeLimiterTaylorHierarchicalVertex()</code>).
%> Here, the full reconstructions instead of linear reconstruction is used.
%>
%> Let @f$\mathcal{A}_q:=\{\mathbf{a}\in\mathbb{N}_0^2\,\big|\,|\mathbf{a}|=q\}@f$
%> be the set of all two-dimensional multi-indices of order @f$q@f$
%> (see <code>phiTaylorPhy()</code> for properties of multi-indices).
%> We determine the correction factor @f$\alpha_{ke}^{(q)}@f$ for each order
%> @f$q \le p@f$ (@f$p@f$ being the polynomial degree of the solution)
%> by computing correction factors using full reconstructions of
%> derivatives,
%> @f[
%> \forall\,\mathbf{a}\in\mathcal{A}_{q-1}\,,\quad
%> c_{k,\mathbf{a},i} \;:=\; \sum_{0 \le |\mathbf{b}| < p - q}
%> C_{k,I(\mathbf{a}+\mathbf{b})} \, \phi_{k,I(\mathbf{b})}(\mathbf{x}_{ki})\,,
%> @f]
%> where the indices of the corresponding degrees of freedom are given by
%> @f$I(\mathbf{a})@f$ (cf. <code>mult2ind()</code>).
%> The correction factor is defined as
%> @f[
%> \alpha_{ke}^{(q)} = \min_{\mathbf{a}\in\mathcal{A}_{q-1}} \, \alpha_{k\mathbf{a}}^{(q)} \,,
%> \qquad\text{with}\quad
%> \alpha_{k\mathbf{a}}^{(q)} \;:=\; \min_i \left\{
%> \begin{array}{llr}
%> (c_{k,\mathbf{a},i}^\mathrm{max} - c_{k,\mathbf{a},c})\big/(c_{k,\mathbf{a},i} - c_{k,\mathbf{a},c})&\text{if}\;&c_{k,\mathbf{a},i} > c_{k,\mathbf{a},i}^\mathrm{max}\\
%> 1 &\text{if}\;&c_{k,\mathbf{a},i}^\mathrm{min}\le c_{k,\mathbf{a},i}\le c_{k,\mathbf{a},i}^\mathrm{max}\\
%> (c_{k,\mathbf{a},i}^\mathrm{min} - c_{k,\mathbf{a},c})\big/(c_{k,\mathbf{a},i} - c_{k,\mathbf{a},c})&\text{if}\;&c_{k,\mathbf{a},i} < c_{k,\mathbf{a},i}^\mathrm{min}
%> \end{array}\right\}\;,
%> @f]
%> where @f$c_{k,\mathbf{a},i}^\mathrm{min}@f$, @f$c_{k,\mathbf{a},i}^\mathrm{max}@f$
%> are minimum- and maximum values of above reconstruction as determined by
%> <code>computeMinMaxV0TElementPatch()</code>, @f$c_{k,\mathbf{a},i}@f$
%> are the function values in the vertices
%> @f$\mathbf{x}_{ki}, i\in\{1,2,3\}@f$, and @f$c_{k,\mathbf{a},c}@f$
%> is the function value in the centroid @f$\mathbf{x}_{kc}@f$ of element
%> @f$T_k@f$.
%>
%> We begin with the highest-order degrees of freedom, compute each
%> correction coefficient and apply it to all higher order degrees of
%> freedom.
%> The limited solution becomes then
%> @f[
%> c_h(\mathbf{x}) = C_{k,1} \,\phi_{k1}
%> + \alpha_{ke}^{(1)} C_{k,2} \,\phi_{k2} (\mathbf{x})
%> + \alpha_{ke}^{(1)} C_{k,3} \,\phi_{k3} (\mathbf{x})
%> + \sum_{i=4}^{N} (\alpha_{ke}^{(1)} \cdots
%> \alpha_{ke}^{(|\mathbf{a}_i|)} ) C_{ki} \, \phi_{ki}(\mathbf{x}) \,.
%> @f]
%>
%> For details to the stricter hierarchical vertex based limiter, see
%> @ref RAWFK2016.
%>
%> @param g The lists describing the geometric and topological
%> properties of a triangulation (see
%> <code>generateGridData()</code>)
%> @f$[1 \times 1 \text{ struct}]@f$
%> @param dataTaylor The representation matrix of the unlimited function
%> @f$c_h@f$. @f$[K \times N]@f$
%> @param markV0TbdrD <code>logical</code> arrays that mark each triangles
%> (Dirichlet boundary) vertices on which additional
%> function values should be considered during the
%> slope limiting routine. @f$[K \times 3]@f$
%> @param dataV0T The function values for (Dirichlet boundary) vertices
%> specified by <code>markV0TbdrD</code>. @f$[K \times 3]@f$
%> @param basesOnQuad A struct containing precomputed values of (Taylor) basis
%> functions on quadrature points. Must provide at
%> least phiTaylorV0T.
%> @retval dataTaylorLim The representation matrix of the limited function
%> @f$\mathsf{\Phi}^\mathrm{Taylor}c_h@f$. @f$[K \times N]@f$
%>
%> This file is part of FESTUNG
%>
%> @copyright 2014-2016 Florian Frank, Balthasar Reuter, Vadym Aizinger
%>
%> @par License
%> @parblock
%> This program is free software: you can redistribute it and/or modify
%> it under the terms of the GNU General Public License as published by
%> the Free Software Foundation, either version 3 of the License, or
%> (at your option) any later version.
%>
%> This program is distributed in the hope that it will be useful,
%> but WITHOUT ANY WARRANTY; without even the implied warranty of
%> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%> GNU General Public License for more details.
%>
%> You should have received a copy of the GNU General Public License
%> along with this program. If not, see <http:https://www.gnu.org/licenses/>.
%> @endparblock
%
function dataTaylorLim = applySlopeLimiterTaylorStrict(g, dataTaylor, markV0TbdrD, dataV0T, basesOnQuad)
% Number of elements and polynomial degree
[K, N] = size(dataTaylor);
p = (sqrt(8*N+1)-3)/2;
% Check function arguments that are directly used
validateattributes(dataTaylor, {'numeric'}, {'size', [g.numT NaN]}, mfilename, 'dataTaylor');
assert(size(dataTaylor, 2) >= 3, 'Number of local degrees of freedom in dataTaylor does not correspond to p>=1')
validateattributes(basesOnQuad, {'struct'}, {}, mfilename, 'basesOnQuad')
validateattributes(basesOnQuad.phiTaylorV0T, {'numeric'}, {'size', [K 3 N]}, mfilename, 'phiTaylorV0T');
% Initialize limited coefficients
dataTaylorLim = dataTaylor;
for ord = p : -1 : 1
alphaOrd = ones(K, 1);
indDOF = ord * (ord + 1) / 2 + 1 : N;
% Number of necessary reconstructions is equal to polynomial degree
for i = 1 : ord
% Determine multi-indices for reconstruction, shift them by multi-index
% of the centroid, and compute corresponding linear indices
pReconstruction = p - ord + 1;
mult = bsxfun(@plus, [ord - i, i - 1], multiindex(pReconstruction));
ind = mult2ind(mult);
% Compute limiter parameters for each element
valV0T = computeFuncDiscAtPoints(dataTaylorLim(:, ind), basesOnQuad.phiTaylorV0T(:, :, 1 : size(ind, 1)));
if ord > 1
alphaTmp = computeVertexBasedLimiter(g, dataTaylorLim(:, ind(1)), valV0T, markV0TbdrD, valV0T);
else
alphaTmp = computeVertexBasedLimiter(g, dataTaylorLim(:, ind(1)), valV0T, markV0TbdrD, dataV0T);
end
alphaOrd = min(alphaOrd, alphaTmp);
end
% Apply limiter
dataTaylorLim(:, indDOF) = bsxfun(@times, alphaOrd, dataTaylorLim(:, indDOF));
end
end % function