forked from FESTUNG/FESTUNG
-
Notifications
You must be signed in to change notification settings - Fork 0
/
assembleMatEdgePhiIntPhiIntNu.m
169 lines (163 loc) · 7.25 KB
/
assembleMatEdgePhiIntPhiIntNu.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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
% Assembles two matrices containing integrals over edges of products of two
% basis functions from the interior of each element with a component of the
% edge normal.
%===============================================================================
%> @file assembleMatEdgePhiIntPhiIntNu.m
%>
%> @brief Assembles two matrices containing integrals over edges of products of
%> two basis functions from the interior of each element, multiplied by
%> the corresponding component of the unit normal.
%===============================================================================
%>
%> @brief Assembles two matrices @f$\mathsf{{Q}}^m_\mathrm{N}, m\in\{1,2\}@f$
%> containing integrals over edges of products of two basis functions
%> from the interior of each element with a component of the
%> edge normal.
%>
%> The matrix @f$\mathsf{{Q}}^m_\mathrm{N} \in \mathbb{R}^{KN\times KN}@f$
%> is block diagonal and defined as
%> @f[
%> [\mathsf{{Q}}^m_\mathrm{N}]_{(k-1)N+i,(k-1)N+j} =
%> \sum_{E_{kn} \in \partial T_k \cap \mathcal{E}_N}
%> \int_{E_{kn}} \nu_{kn}^m \varphi_{ki} \varphi_{kj} \mathrm{d}s \,.
%> @f]
%> with @f$\nu_{kn}@f$ the @f$m@f$-th component (@f$m\in\{1,2\}@f$) of the edge
%> normal.
%> All other entries are zero.
%> To allow for vectorization, the assembly is reformulated as
%> @f[
%> \mathsf{{Q}}^m_\mathrm{N} = \sum_{n=1}^3
%> \begin{bmatrix}
%> \delta_{E_{1n}\in\mathcal{E}_\mathrm{N}} & & \\
%> & ~\ddots~ & \\
%> & & \delta_{E_{Kn}\in\mathcal{E}_\mathrm{N}}
%> \end{bmatrix} \circ \begin{bmatrix}
%> \nu^m_{1n} | E_{1n} | & & \\
%> & ~\ddots~ & \\
%> & & \nu^m_{Kn} | E_{Kn} |
%> \end{bmatrix} \otimes [\hat{\mathsf{{S}}}^\mathrm{diag}]_{:,:,n}\;,
%> @f]
%> where @f$\delta_{E_{kn}\in\mathcal{E}_\mathrm{N}}@f$ denotes the Kronecker
%> delta, @f$\circ@f$ denotes the Hadamard product, and @f$\otimes@f$ denotes
%> the Kronecker product.
%>
%> The entries of matrix
%> @f$\hat{\mathsf{{S}}}^\mathrm{diag}\in\mathbb{R}^{N\times N\times3}@f$
%> are given by
%> @f[
%> [\hat{\mathsf{{S}}}^\mathrm{diag}]_{i,j,n} =
%> \int_0^1 \hat{\varphi}_i \circ \hat{\mathbf{\gamma}}_n(s)
%> \hat{\varphi}_j\circ \hat{\mathbf{\gamma}}_n(s) \mathrm{d}s \,,
%> @f]
%> where the mapping @f$\hat{\mathbf{\gamma}}_n@f$ is defined in
%> <code>gammaMap()</code>.
%>
%> It is essentially the same as the diagonal part of
%> <code>assembleMatEdgePhiPhiNu()</code>.
%>
%> @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 markE0Tbdr <code>logical</code> arrays that mark each triangles
%> (boundary) edges on which the matrix blocks should be
%> assembled @f$[K \times 3]@f$
%> @param refEdgePhiIntPhiInt Local matrix
%> @f$\hat{\mathsf{{S}}}^\text{diag}@f$ as provided
%> by <code>integrateRefEdgePhiIntPhiInt()</code>.
%> @f$[N \times N \times 3]@f$
%> @param areaNuE0Tbdr (optional) argument to provide precomputed values
%> for the products of <code>markE0Tbdr</code>,
%> <code>g.areaE0T</code>, and <code>g.nuE0T</code>
%> @f$[3 \times 2 \text{ cell}]@f$
%> @retval ret The assembled matrices @f$[2 \times 1 \text{ cell}]@f$
%>
%> This file is part of FESTUNG
%>
%> @copyright 2014-2015 Florian Frank, Balthasar Reuter, Vadym Aizinger
%> Modified by Hennes Hajduk, 2016-04-06
%>
%> @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 <https://www.gnu.org/licenses/>.
%> @endparblock
%
function ret = assembleMatEdgePhiIntPhiIntNu(g, markE0Tbdr, refEdgePhiIntPhiInt, areaNuE0Tbdr)
% Extract dimensions
K = g.numT; N = size(refEdgePhiIntPhiInt, 1);
% Check function arguments that are directly used
validateattributes(markE0Tbdr, {'logical'}, {'size', [K 3]}, mfilename, 'markE0Tbdr');
validateattributes(refEdgePhiIntPhiInt, {'numeric'}, {'size', [N N 3]}, mfilename, 'refEdgePhiIntPhiInt');
if nargin > 3
ret = assembleMatEdgePhiIntPhiIntNu_withAreaNuE0Tbdr(g, refEdgePhiIntPhiInt, areaNuE0Tbdr);
elseif isfield(g, 'areaNuE0T')
ret = assembleMatEdgePhiIntPhiIntNu_noAreaNuE0Tbdr_withAreaNuE0T(g, markE0Tbdr, refEdgePhiIntPhiInt);
else
ret = assembleMatEdgePhiIntPhiIntNu_noAreaNuE0Tbdr_noAreaNuE0T(g, markE0Tbdr, refEdgePhiIntPhiInt);
end
end % function
%
%===============================================================================
%> @brief Helper function for the case that assembleMatEdgePhiIntPhiIntNu()
%> was called with a precomputed field areaNuE0Tbdr.
%
function ret = assembleMatEdgePhiIntPhiIntNu_withAreaNuE0Tbdr(g, refEdgePhiIntPhiInt, areaNuE0Tbdr)
% Extract dimensions
K = g.numT; N = size(refEdgePhiIntPhiInt, 1);
% Assemble matrices
ret = cell(2,1);
ret{1} = sparse(K*N, K*N);
ret{2} = sparse(K*N, K*N);
for n = 1 : 3
ret{1} = ret{1} + kron(spdiags(areaNuE0Tbdr{n,1}, 0,K,K), refEdgePhiIntPhiInt(:,:,n));
ret{2} = ret{2} + kron(spdiags(areaNuE0Tbdr{n,2}, 0,K,K), refEdgePhiIntPhiInt(:,:,n));
end % for
end % function
%
%===============================================================================
%> @brief Helper function for the case that assembleMatEdgePhiIntPhiIntNu()
%> was called without a precomputed field areaNuE0Tbdr and parameter g provides
%> a precomputed field areaNuE0T.
%
function ret = assembleMatEdgePhiIntPhiIntNu_noAreaNuE0Tbdr_withAreaNuE0T(g, markE0Tbdr, refEdgePhiIntPhiInt)
% Extract dimensions
K = g.numT; N = size(refEdgePhiIntPhiInt, 1);
% Assemble matrices
ret = cell(2,1);
ret{1} = sparse(K*N, K*N);
ret{2} = sparse(K*N, K*N);
for n = 1 : 3
ret{1} = ret{1} + kron(spdiags(markE0Tbdr(:,n).*g.areaNuE0T{n,1}, 0,K,K), refEdgePhiIntPhiInt(:,:,n));
ret{2} = ret{2} + kron(spdiags(markE0Tbdr(:,n).*g.areaNuE0T{n,2}, 0,K,K), refEdgePhiIntPhiInt(:,:,n));
end % for
end % function
%
%===============================================================================
%> @brief Helper function for the case that assembleMatEdgePhiIntPhiIntNu()
%> was called without a precomputed field areaNuE0Tbdr and parameter g provides
%> no precomputed field areaNuE0T.
%
function ret = assembleMatEdgePhiIntPhiIntNu_noAreaNuE0Tbdr_noAreaNuE0T(g, markE0Tbdr, refEdgePhiIntPhiInt)
% Extract dimensions
K = g.numT; N = size(refEdgePhiIntPhiInt, 1);
% Assemble matrices
ret = cell(2,1);
ret{1} = sparse(K*N, K*N);
ret{2} = sparse(K*N, K*N);
for n = 1 : 3
QNkn = markE0Tbdr(:,n) .* g.areaE0T(:,n);
ret{1} = ret{1} + kron(spdiags(QNkn .* g.nuE0T(:,n,1), 0,K,K), refEdgePhiIntPhiInt(:,:,n));
ret{2} = ret{2} + kron(spdiags(QNkn .* g.nuE0T(:,n,2), 0,K,K), refEdgePhiIntPhiInt(:,:,n));
end % for
end % function