Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Basis.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Functions for computing the Dubiner basis functions in DG methods
9 : : \details This file contains functionality for computing the basis functions
10 : : and relating coordinates transformation functions used in discontinuous
11 : : Galerkin methods for variaous orders of numerical representation. The basis
12 : : functions chosen for the DG method are the Dubiner basis, which are
13 : : Legendre polynomials modified for tetrahedra, which are defined only on
14 : : the reference/master tetrahedron.
15 : : \see [1] https://doi.org/10.1007/BF01060030
16 : : \see [2] https://doi.org/10.1093/imamat/hxh111
17 : : */
18 : : // *****************************************************************************
19 : :
20 : : #include <array>
21 : : #include "QuinoaConfig.hpp"
22 : :
23 : : #include "Basis.hpp"
24 : : #include "Vector.hpp"
25 : : #include "Mass.hpp"
26 : :
27 : : std::array< tk::real, 3 >
28 : 133326316 : tk::eval_gp ( const std::size_t igp,
29 : : const std::array< std::array< tk::real, 3>, 3 >& coordfa,
30 : : const std::array< std::vector< tk::real >, 2 >& coordgp )
31 : : // *****************************************************************************
32 : : // Compute the coordinates of quadrature points for face integral in physical
33 : : // space
34 : : //! \param[in] igp Index of quadrature points
35 : : //! \param[in] coordfa Array of nodal coordinates for face element
36 : : //! \param[in] coordgp Array of coordinates for quadrature points in reference
37 : : //! space
38 : : //! \return Array of coordinates for quadrature points in physical space
39 : : // *****************************************************************************
40 : : {
41 : : // Barycentric coordinates for the triangular face
42 : 133326316 : auto shp1 = 1.0 - coordgp[0][igp] - coordgp[1][igp];
43 : : auto shp2 = coordgp[0][igp];
44 : : auto shp3 = coordgp[1][igp];
45 : :
46 : : // Transformation of the quadrature point from the 2D reference/master
47 : : // element to physical space, to obtain its physical (x,y,z) coordinates.
48 : 133326316 : return {{ coordfa[0][0]*shp1 + coordfa[1][0]*shp2 + coordfa[2][0]*shp3,
49 : 133326316 : coordfa[0][1]*shp1 + coordfa[1][1]*shp2 + coordfa[2][1]*shp3,
50 : 133326316 : coordfa[0][2]*shp1 + coordfa[1][2]*shp2 + coordfa[2][2]*shp3 }};
51 : : }
52 : :
53 : : std::array< tk::real, 3 >
54 : 65325203 : tk::eval_gp ( const std::size_t igp,
55 : : const std::array< std::array< tk::real, 3>, 4 >& coord,
56 : : const std::array< std::vector< tk::real >, 3 >& coordgp )
57 : : // *****************************************************************************
58 : : // Compute the coordinates of quadrature points for volume integral in
59 : : // physical space
60 : : //! \param[in] igp Index of quadrature points
61 : : //! \param[in] coord Array of nodal coordinates for tetrahedron element
62 : : //! \param[in] coordgp Array of coordinates for quadrature points in reference space
63 : : //! \return Array of coordinates for quadrature points in physical space
64 : : // *****************************************************************************
65 : : {
66 : : // Barycentric coordinates for the tetradedron element
67 : 65325203 : auto shp1 = 1.0 - coordgp[0][igp] - coordgp[1][igp] - coordgp[2][igp];
68 : : auto shp2 = coordgp[0][igp];
69 : : auto shp3 = coordgp[1][igp];
70 : : auto shp4 = coordgp[2][igp];
71 : :
72 : : // Transformation of the quadrature point from the reference/master
73 : : // element to physical space, to obtain its physical (x,y,z) coordinates.
74 : : return {{
75 : 65325203 : coord[0][0]*shp1 + coord[1][0]*shp2 + coord[2][0]*shp3 + coord[3][0]*shp4,
76 : 65325203 : coord[0][1]*shp1 + coord[1][1]*shp2 + coord[2][1]*shp3 + coord[3][1]*shp4,
77 : 65325203 : coord[0][2]*shp1 + coord[1][2]*shp2 + coord[2][2]*shp3 + coord[3][2]*shp4 }};
78 : : }
79 : :
80 : : std::array< std::vector<tk::real>, 3 >
81 : 1861455 : tk::eval_dBdxi( const std::size_t ndof,
82 : : const std::array< tk::real, 3 >& coordgp )
83 : : // *****************************************************************************
84 : : // Compute the derivatives of Dubiner basis wrt. reference coordinates
85 : : //! \param[in] ndof Number of degrees of freedom
86 : : //! \param[in] coordgp Coordinates in ref element where derivatives are needed
87 : : //! \return Array of the derivatives of basis functions
88 : : // *****************************************************************************
89 : : {
90 : : // Initialize array
91 : : std::array< std::vector< tk::real >, 3 > dBdxi;
92 [ + + ]: 7445820 : for (std::size_t idir=0; idir<3; ++idir) {
93 [ + - ]: 5584365 : dBdxi[idir].resize(ndof, 0.0);
94 : : }
95 : :
96 : : // high-order basis
97 [ + - ]: 1861455 : if (ndof > 1) {
98 : 1861455 : dBdxi[0][0] = 0.0;
99 : 1861455 : dBdxi[1][0] = 0.0;
100 : 1861455 : dBdxi[2][0] = 0.0;
101 : :
102 : 1861455 : dBdxi[0][1] = 2.0;
103 : 1861455 : dBdxi[1][1] = 1.0;
104 : 1861455 : dBdxi[2][1] = 1.0;
105 : :
106 : 1861455 : dBdxi[0][2] = 0.0;
107 : 1861455 : dBdxi[1][2] = 3.0;
108 : 1861455 : dBdxi[2][2] = 1.0;
109 : :
110 : 1861455 : dBdxi[0][3] = 0.0;
111 : 1861455 : dBdxi[1][3] = 0.0;
112 : 1861455 : dBdxi[2][3] = 4.0;
113 : :
114 [ + - ]: 1861455 : if (ndof > 4) {
115 : 1861455 : dBdxi[0][4] = 12.0 * coordgp[0] + 6.0 * coordgp[1]
116 : 1861455 : + 6.0 * coordgp[2] - 6.0;
117 : 1861455 : dBdxi[1][4] = 6.0 * coordgp[0] + 2.0 * coordgp[1]
118 : 1861455 : + 2.0 * coordgp[2] - 2.0;
119 : 1861455 : dBdxi[2][4] = 6.0 * coordgp[0] + 2.0 * coordgp[1]
120 : 1861455 : + 2.0 * coordgp[2] - 2.0;
121 : :
122 : 1861455 : dBdxi[0][5] = 10.0 * coordgp[1] + 2.0 * coordgp[2] - 2.0;
123 : 1861455 : dBdxi[1][5] = 10.0 * coordgp[0] + 10.0 * coordgp[1]
124 : 1861455 : + 6.0 * coordgp[2] - 6.0;
125 : 1861455 : dBdxi[2][5] = 2.0 * coordgp[0] + 6.0 * coordgp[1]
126 : 1861455 : + 2.0 * coordgp[2] - 2.0;
127 : :
128 : 1861455 : dBdxi[0][6] = 12.0 * coordgp[2] - 2.0;
129 : 1861455 : dBdxi[1][6] = 6.0 * coordgp[2] - 1.0;
130 : 1861455 : dBdxi[2][6] = 12.0 * coordgp[0] + 6.0 * coordgp[1]
131 : 1861455 : + 12.0 * coordgp[2] - 7.0;
132 : :
133 : 1861455 : dBdxi[0][7] = 0.0;
134 : 1861455 : dBdxi[1][7] = 20.0 * coordgp[1] + 8.0 * coordgp[2] - 8.0;
135 : 1861455 : dBdxi[2][7] = 8.0 * coordgp[1] + 2.0 * coordgp[2] - 2.0;
136 : :
137 : 1861455 : dBdxi[0][8] = 0.0;
138 : 1861455 : dBdxi[1][8] = 18.0 * coordgp[2] - 3.0;
139 : 1861455 : dBdxi[2][8] = 18.0 * coordgp[1] + 12.0 * coordgp[2] - 7.0;
140 : :
141 : 1861455 : dBdxi[0][9] = 0.0;
142 : 1861455 : dBdxi[1][9] = 0.0;
143 : 1861455 : dBdxi[2][9] = 30.0 * coordgp[2] - 10.0;
144 : : }
145 : : }
146 : :
147 : 1861455 : return dBdxi;
148 : : }
149 : :
150 : : std::array< std::vector<tk::real>, 3 >
151 : 8266911 : tk::eval_dBdx_p1( const std::size_t ndof,
152 : : const std::array< std::array< tk::real, 3 >, 3 >& jacInv )
153 : : // *****************************************************************************
154 : : // Compute the derivatives of basis functions for DG(P1)
155 : : //! \param[in] ndof Number of degrees of freedom
156 : : //! \param[in] jacInv Array of the inverse of Jacobian
157 : : //! \return Array of the derivatives of basis functions
158 : : // *****************************************************************************
159 : : {
160 : : // The derivatives of the basis functions dB/dx are easily calculated
161 : : // via a transformation to the reference space as,
162 : : // dB/dx = dB/dxi . dxi/dx,
163 : : // where, x = (x,y,z) are the physical coordinates, and
164 : : // xi = (xi, eta, zeta) are the reference coordinates.
165 : : // The matrix dxi/dx is the inverse of the Jacobian of transformation
166 : : // and the matrix vector product has to be calculated. This follows.
167 : :
168 : : std::array< std::vector<tk::real>, 3 > dBdx;
169 [ + - ][ + - ]: 8266911 : dBdx[0].resize( ndof, 0 );
170 [ + - ][ + - ]: 8266911 : dBdx[1].resize( ndof, 0 );
171 [ + - ]: 8266911 : dBdx[2].resize( ndof, 0 );
172 : :
173 : : auto db2dxi1 = 2.0;
174 : : auto db2dxi2 = 1.0;
175 : : auto db2dxi3 = 1.0;
176 : :
177 : : auto db3dxi1 = 0.0;
178 : : auto db3dxi2 = 3.0;
179 : : auto db3dxi3 = 1.0;
180 : :
181 : : auto db4dxi1 = 0.0;
182 : : auto db4dxi2 = 0.0;
183 : : auto db4dxi3 = 4.0;
184 : :
185 : 8266911 : dBdx[0][1] = db2dxi1 * jacInv[0][0]
186 : 8266911 : + db2dxi2 * jacInv[1][0]
187 : 8266911 : + db2dxi3 * jacInv[2][0];
188 : :
189 : 8266911 : dBdx[1][1] = db2dxi1 * jacInv[0][1]
190 : 8266911 : + db2dxi2 * jacInv[1][1]
191 : 8266911 : + db2dxi3 * jacInv[2][1];
192 : :
193 : 8266911 : dBdx[2][1] = db2dxi1 * jacInv[0][2]
194 : 8266911 : + db2dxi2 * jacInv[1][2]
195 : 8266911 : + db2dxi3 * jacInv[2][2];
196 : :
197 : 8266911 : dBdx[0][2] = db3dxi1 * jacInv[0][0]
198 : 8266911 : + db3dxi2 * jacInv[1][0]
199 : 8266911 : + db3dxi3 * jacInv[2][0];
200 : :
201 : 8266911 : dBdx[1][2] = db3dxi1 * jacInv[0][1]
202 : 8266911 : + db3dxi2 * jacInv[1][1]
203 : 8266911 : + db3dxi3 * jacInv[2][1];
204 : :
205 : 8266911 : dBdx[2][2] = db3dxi1 * jacInv[0][2]
206 : 8266911 : + db3dxi2 * jacInv[1][2]
207 : 8266911 : + db3dxi3 * jacInv[2][2];
208 : :
209 : 8266911 : dBdx[0][3] = db4dxi1 * jacInv[0][0]
210 : 8266911 : + db4dxi2 * jacInv[1][0]
211 : 8266911 : + db4dxi3 * jacInv[2][0];
212 : :
213 : 8266911 : dBdx[1][3] = db4dxi1 * jacInv[0][1]
214 : 8266911 : + db4dxi2 * jacInv[1][1]
215 : 8266911 : + db4dxi3 * jacInv[2][1];
216 : :
217 : 8266911 : dBdx[2][3] = db4dxi1 * jacInv[0][2]
218 : 8266911 : + db4dxi2 * jacInv[1][2]
219 : 8266911 : + db4dxi3 * jacInv[2][2];
220 : :
221 : 8266911 : return dBdx;
222 : : }
223 : :
224 : : void
225 : 14707176 : tk::eval_dBdx_p2( const std::size_t igp,
226 : : const std::array< std::vector< tk::real >, 3 >& coordgp,
227 : : const std::array< std::array< tk::real, 3 >, 3 >& jacInv,
228 : : std::array< std::vector<tk::real>, 3 >& dBdx )
229 : : // *****************************************************************************
230 : : // Compute the derivatives of Dubiner basis function for DG(P2)
231 : : //! \param[in] igp Index of quadrature points
232 : : //! \param[in] coordgp Gauss point coordinates for tetrahedron element
233 : : //! \param[in] jacInv Array of the inverse of Jacobian
234 : : //! \param[in,out] dBdx Array of the derivatives of basis function
235 : : // *****************************************************************************
236 : : {
237 : 14707176 : auto db5dxi1 = 12.0 * coordgp[0][igp] + 6.0 * coordgp[1][igp]
238 : 14707176 : + 6.0 * coordgp[2][igp] - 6.0;
239 : 14707176 : auto db5dxi2 = 6.0 * coordgp[0][igp] + 2.0 * coordgp[1][igp]
240 : 14707176 : + 2.0 * coordgp[2][igp] - 2.0;
241 : : auto db5dxi3 = 6.0 * coordgp[0][igp] + 2.0 * coordgp[1][igp]
242 : : + 2.0 * coordgp[2][igp] - 2.0;
243 : :
244 : 14707176 : auto db6dxi1 = 10.0 * coordgp[1][igp] + 2.0 * coordgp[2][igp] - 2.0;
245 : 14707176 : auto db6dxi2 = 10.0 * coordgp[0][igp] + 10.0 * coordgp[1][igp]
246 : 14707176 : + 6.0 * coordgp[2][igp] - 6.0;
247 : 14707176 : auto db6dxi3 = 2.0 * coordgp[0][igp] + 6.0 * coordgp[1][igp]
248 : 14707176 : + 2.0 * coordgp[2][igp] - 2.0;
249 : :
250 : 14707176 : auto db7dxi1 = 12.0 * coordgp[2][igp] - 2.0;
251 : 14707176 : auto db7dxi2 = 6.0 * coordgp[2][igp] - 1.0;
252 : : auto db7dxi3 = 12.0 * coordgp[0][igp] + 6.0 * coordgp[1][igp]
253 : 14707176 : + 12.0 * coordgp[2][igp] - 7.0;
254 : :
255 : : auto db8dxi1 = 0;
256 : 14707176 : auto db8dxi2 = 20.0 * coordgp[1][igp] + 8.0 * coordgp[2][igp] - 8.0;
257 : 14707176 : auto db8dxi3 = 8.0 * coordgp[1][igp] + 2.0 * coordgp[2][igp] - 2.0;
258 : :
259 : : auto db9dxi1 = 0;
260 : 14707176 : auto db9dxi2 = 18.0 * coordgp[2][igp] - 3.0;
261 : 14707176 : auto db9dxi3 = 18.0 * coordgp[1][igp] + 12.0 * coordgp[2][igp] - 7.0;
262 : :
263 : : auto db10dxi1 = 0;
264 : : auto db10dxi2 = 0;
265 : 14707176 : auto db10dxi3 = 30.0 * coordgp[2][igp] - 10.0;
266 : :
267 : 14707176 : dBdx[0][4] = db5dxi1 * jacInv[0][0]
268 : 14707176 : + db5dxi2 * jacInv[1][0]
269 : 14707176 : + db5dxi3 * jacInv[2][0];
270 : :
271 : 14707176 : dBdx[1][4] = db5dxi1 * jacInv[0][1]
272 : 14707176 : + db5dxi2 * jacInv[1][1]
273 : 14707176 : + db5dxi3 * jacInv[2][1];
274 : :
275 : 14707176 : dBdx[2][4] = db5dxi1 * jacInv[0][2]
276 : 14707176 : + db5dxi2 * jacInv[1][2]
277 : 14707176 : + db5dxi3 * jacInv[2][2];
278 : :
279 : 14707176 : dBdx[0][5] = db6dxi1 * jacInv[0][0]
280 : 14707176 : + db6dxi2 * jacInv[1][0]
281 : 14707176 : + db6dxi3 * jacInv[2][0];
282 : :
283 : 14707176 : dBdx[1][5] = db6dxi1 * jacInv[0][1]
284 : 14707176 : + db6dxi2 * jacInv[1][1]
285 : 14707176 : + db6dxi3 * jacInv[2][1];
286 : :
287 : 14707176 : dBdx[2][5] = db6dxi1 * jacInv[0][2]
288 : 14707176 : + db6dxi2 * jacInv[1][2]
289 : 14707176 : + db6dxi3 * jacInv[2][2];
290 : :
291 : 14707176 : dBdx[0][6] = db7dxi1 * jacInv[0][0]
292 : 14707176 : + db7dxi2 * jacInv[1][0]
293 : 14707176 : + db7dxi3 * jacInv[2][0];
294 : :
295 : 14707176 : dBdx[1][6] = db7dxi1 * jacInv[0][1]
296 : 14707176 : + db7dxi2 * jacInv[1][1]
297 : 14707176 : + db7dxi3 * jacInv[2][1];
298 : :
299 : 14707176 : dBdx[2][6] = db7dxi1 * jacInv[0][2]
300 : 14707176 : + db7dxi2 * jacInv[1][2]
301 : 14707176 : + db7dxi3 * jacInv[2][2];
302 : :
303 : 14707176 : dBdx[0][7] = db8dxi1 * jacInv[0][0]
304 : 14707176 : + db8dxi2 * jacInv[1][0]
305 : 14707176 : + db8dxi3 * jacInv[2][0];
306 : :
307 : 14707176 : dBdx[1][7] = db8dxi1 * jacInv[0][1]
308 : 14707176 : + db8dxi2 * jacInv[1][1]
309 : 14707176 : + db8dxi3 * jacInv[2][1];
310 : :
311 : 14707176 : dBdx[2][7] = db8dxi1 * jacInv[0][2]
312 : 14707176 : + db8dxi2 * jacInv[1][2]
313 : 14707176 : + db8dxi3 * jacInv[2][2];
314 : :
315 : 14707176 : dBdx[0][8] = db9dxi1 * jacInv[0][0]
316 : 14707176 : + db9dxi2 * jacInv[1][0]
317 : 14707176 : + db9dxi3 * jacInv[2][0];
318 : :
319 : 14707176 : dBdx[1][8] = db9dxi1 * jacInv[0][1]
320 : 14707176 : + db9dxi2 * jacInv[1][1]
321 : 14707176 : + db9dxi3 * jacInv[2][1];
322 : :
323 : 14707176 : dBdx[2][8] = db9dxi1 * jacInv[0][2]
324 : 14707176 : + db9dxi2 * jacInv[1][2]
325 : 14707176 : + db9dxi3 * jacInv[2][2];
326 : :
327 : 14707176 : dBdx[0][9] = db10dxi1 * jacInv[0][0]
328 : 14707176 : + db10dxi2 * jacInv[1][0]
329 : 14707176 : + db10dxi3 * jacInv[2][0];
330 : :
331 : 14707176 : dBdx[1][9] = db10dxi1 * jacInv[0][1]
332 : 14707176 : + db10dxi2 * jacInv[1][1]
333 : 14707176 : + db10dxi3 * jacInv[2][1];
334 : :
335 : 14707176 : dBdx[2][9] = db10dxi1 * jacInv[0][2]
336 : 14707176 : + db10dxi2 * jacInv[1][2]
337 : 14707176 : + db10dxi3 * jacInv[2][2];
338 : 14707176 : }
339 : :
340 : : std::vector< tk::real >
341 : 289902982 : tk::eval_basis( const std::size_t ndof,
342 : : const tk::real xi,
343 : : const tk::real eta,
344 : : const tk::real zeta )
345 : : // *****************************************************************************
346 : : // Compute the Dubiner basis functions
347 : : //! \param[in] ndof Number of degrees of freedom
348 : : //! \param[in] xi,eta,zeta Coordinates for quadrature points in reference space
349 : : //! \return Vector of basis functions
350 : : // *****************************************************************************
351 : : {
352 : : // Array of basis functions
353 : 289902982 : std::vector< tk::real > B( ndof, 1.0 );
354 : :
355 [ + + ]: 289902982 : if ( ndof > 1 ) // DG(P1)
356 : : {
357 : 229157811 : B[1] = 2.0 * xi + eta + zeta - 1.0;
358 : 229157811 : B[2] = 3.0 * eta + zeta - 1.0;
359 : 229157811 : B[3] = 4.0 * zeta - 1.0;
360 : :
361 [ + + ]: 229157811 : if( ndof > 4 ) // DG(P2)
362 : : {
363 : 65740965 : B[4] = 6.0 * xi * xi + eta * eta + zeta * zeta
364 : 65740965 : + 6.0 * xi * eta + 6.0 * xi * zeta + 2.0 * eta * zeta
365 : 65740965 : - 6.0 * xi - 2.0 * eta - 2.0 * zeta + 1.0;
366 : 65740965 : B[5] = 5.0 * eta * eta + zeta * zeta
367 : 65740965 : + 10.0 * xi * eta + 2.0 * xi * zeta + 6.0 * eta * zeta
368 : 65740965 : - 2.0 * xi - 6.0 * eta - 2.0 * zeta + 1.0;
369 : 65740965 : B[6] = 6.0 * zeta * zeta + 12.0 * xi * zeta + 6.0 * eta * zeta - 2.0 * xi
370 : 65740965 : - eta - 7.0 * zeta + 1.0;
371 : 65740965 : B[7] = 10.0 * eta * eta + zeta * zeta + 8.0 * eta * zeta
372 : 65740965 : - 8.0 * eta - 2.0 * zeta + 1.0;
373 : 65740965 : B[8] = 6.0 * zeta * zeta + 18.0 * eta * zeta - 3.0 * eta - 7.0 * zeta
374 : 65740965 : + 1.0;
375 : 65740965 : B[9] = 15.0 * zeta * zeta - 10.0 * zeta + 1.0;
376 : : }
377 : : }
378 : :
379 : 289902982 : return B;
380 : : }
381 : :
382 : : std::vector< tk::real >
383 [ + + ]: 455145443 : tk::eval_state ( ncomp_t ncomp,
384 : : const std::size_t ndof,
385 : : const std::size_t ndof_el,
386 : : const std::size_t e,
387 : : const Fields& U,
388 : : const std::vector< tk::real >& B )
389 : : // *****************************************************************************
390 : : // Compute the state variables for the tetrahedron element
391 : : //! \param[in] ncomp Number of scalar components in this PDE system
392 : : //! \param[in] ndof Maximum number of degrees of freedom
393 : : //! \param[in] ndof_el Number of degrees of freedom for the local element
394 : : //! \param[in] e Index for the tetrahedron element
395 : : //! \param[in] U Solution vector at recent time step
396 : : //! \param[in] B Vector of basis functions
397 : : //! \return Vector of state variable for tetrahedron element
398 : : // *****************************************************************************
399 : : {
400 : : // This is commented for now because that when p0/p1 adaptive with limiter
401 : : // applied, the size of basis will be 10. However, ndof_el will be 4 which
402 : : // leads to a size mismatch in limiter function.
403 : : //Assert( B.size() == ndof_el, "Size mismatch" );
404 : :
405 [ + + ]: 455145443 : if (U.empty()) return {};
406 : :
407 : : // Array of state variable for tetrahedron element
408 : 294420305 : std::vector< tk::real > state( ncomp, 0.0 );
409 : :
410 [ + + ]: 1578054840 : for (ncomp_t c=0; c<ncomp; ++c)
411 : : {
412 [ + + ]: 1283634535 : auto mark = c*ndof;
413 [ + + ]: 1283634535 : state[c] = U( e, mark );
414 : :
415 [ + + ]: 1283634535 : if(ndof_el > 1) // Second order polynomial solution
416 : : {
417 : 979207047 : state[c] += U( e, mark+1 ) * B[1]
418 : 979207047 : + U( e, mark+2 ) * B[2]
419 : 979207047 : + U( e, mark+3 ) * B[3];
420 : : }
421 : :
422 [ + + ]: 1283634535 : if(ndof_el > 4) // Third order polynomial solution
423 : : {
424 : 174496433 : state[c] += U( e, mark+4 ) * B[4]
425 : 174496433 : + U( e, mark+5 ) * B[5]
426 : 174496433 : + U( e, mark+6 ) * B[6]
427 : 174496433 : + U( e, mark+7 ) * B[7]
428 : 174496433 : + U( e, mark+8 ) * B[8]
429 : 174496433 : + U( e, mark+9 ) * B[9];
430 : : }
431 : : }
432 : :
433 : : return state;
434 : : }
435 : :
436 : : std::vector< std::vector< tk::real > >
437 : 0 : tk::DubinerToTaylor( ncomp_t ncomp,
438 : : const std::size_t e,
439 : : const std::size_t ndof,
440 : : const tk::Fields& U,
441 : : const std::vector< std::size_t >& inpoel,
442 : : const tk::UnsMesh::Coords& coord )
443 : : // *****************************************************************************
444 : : // Transform the solution with Dubiner basis to the solution with Taylor basis
445 : : //! \param[in] ncomp Number of scalar components in this PDE system
446 : : //! \param[in] e Id of element whose solution is to be limited
447 : : //! \param[in] ndof Maximum number of degrees of freedom
448 : : //! \param[in] U High-order solution vector with Dubiner basis
449 : : //! \param[in] inpoel Element connectivity
450 : : //! \param[in] coord Array of nodal coordinates
451 : : //! \return High-order solution vector with Taylor basis
452 : : // *****************************************************************************
453 : : {
454 : : std::vector< std::vector< tk::real > >
455 [ - - ][ - - ]: 0 : unk(ncomp, std::vector<tk::real>(ndof, 0.0));
456 : :
457 : : const auto& cx = coord[0];
458 : : const auto& cy = coord[1];
459 : : const auto& cz = coord[2];
460 : :
461 : : std::array< std::vector< tk::real >, 3 > center;
462 [ - - ]: 0 : center[0].resize(1, 0.25);
463 [ - - ]: 0 : center[1].resize(1, 0.25);
464 [ - - ]: 0 : center[2].resize(1, 0.25);
465 : :
466 : : // Evaluate the cell center solution
467 [ - - ]: 0 : for(ncomp_t icomp = 0; icomp < ncomp; icomp++)
468 : : {
469 : 0 : auto mark = icomp * ndof;
470 : 0 : unk[icomp][0] = U(e, mark);
471 : : }
472 : :
473 : : // Evaluate the first order derivative
474 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
475 : 0 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
476 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
477 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
478 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
479 : 0 : }};
480 : :
481 : : auto jacInv =
482 : 0 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
483 : :
484 : : // Compute the derivatives of basis function for DG(P1)
485 [ - - ]: 0 : auto dBdx = tk::eval_dBdx_p1( ndof, jacInv );
486 : :
487 [ - - ]: 0 : if(ndof > 4) {
488 [ - - ]: 0 : tk::eval_dBdx_p2(0, center, jacInv, dBdx);
489 : : }
490 : :
491 [ - - ]: 0 : for(ncomp_t icomp = 0; icomp < ncomp; icomp++)
492 : : {
493 : 0 : auto mark = icomp * ndof;
494 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
495 : : {
496 : 0 : unk[icomp][idir+1] = 0;
497 [ - - ]: 0 : for(std::size_t idof = 1; idof < ndof; idof++)
498 : 0 : unk[icomp][idir+1] += U(e, mark+idof) * dBdx[idir][idof];
499 : : }
500 : : }
501 : :
502 : : // Evaluate the second order derivative if DGP2 is applied
503 : : // The basic idea of the computation follows
504 : : // d2Udx2 = /sum u_i * (d2B_i/dx2)
505 : : // where d2B_i/dx2 = d( dB_i/dxi * dxi/dx ) / dxi * dxi/dx
506 [ - - ]: 0 : if(ndof > 4)
507 : : {
508 : : // Matrix to store the second order derivatives of basis functions in
509 : : // reference domain
510 : 0 : tk::real d2Bdxi2[6][6] =
511 : : { { 12.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
512 : : { 2.0, 10.0, 0.0, 20.0, 0.0, 0.0 },
513 : : { 2.0, 2.0, 12.0, 2.0, 12.0, 30.0 },
514 : : { 6.0, 10.0, 0.0, 0.0, 0.0, 0.0 },
515 : : { 6.0, 2.0, 12.0, 0.0, 0.0, 0.0 },
516 : : { 2.0, 6.0, 6.0, 8.0, 18.0, 0.0 } };
517 : :
518 : : // Transform matrix to convert the second order derivatives of basis
519 : : // function in reference domain to the one in physical domain
520 : : tk::real d2xdxi2[6][6];
521 : :
522 : 0 : d2xdxi2[0][0] = jacInv[0][0] * jacInv[0][0];
523 : 0 : d2xdxi2[0][1] = jacInv[1][0] * jacInv[1][0];
524 : 0 : d2xdxi2[0][2] = jacInv[2][0] * jacInv[2][0];
525 : 0 : d2xdxi2[0][3] = jacInv[0][0] * jacInv[1][0] * 2.0;
526 : 0 : d2xdxi2[0][4] = jacInv[0][0] * jacInv[2][0] * 2.0;
527 : 0 : d2xdxi2[0][5] = jacInv[1][0] * jacInv[2][0] * 2.0;
528 : :
529 : 0 : d2xdxi2[1][0] = jacInv[0][1] * jacInv[0][1];
530 : 0 : d2xdxi2[1][1] = jacInv[1][1] * jacInv[1][1];
531 : 0 : d2xdxi2[1][2] = jacInv[2][1] * jacInv[2][1];
532 : 0 : d2xdxi2[1][3] = jacInv[0][1] * jacInv[1][1] * 2.0;
533 : 0 : d2xdxi2[1][4] = jacInv[0][1] * jacInv[2][1] * 2.0;
534 : 0 : d2xdxi2[1][5] = jacInv[1][1] * jacInv[2][1] * 2.0;
535 : :
536 : 0 : d2xdxi2[2][0] = jacInv[0][2] * jacInv[0][2];
537 : 0 : d2xdxi2[2][1] = jacInv[1][2] * jacInv[1][2];
538 : 0 : d2xdxi2[2][2] = jacInv[2][2] * jacInv[2][2];
539 : 0 : d2xdxi2[2][3] = jacInv[0][2] * jacInv[1][2] * 2.0;
540 : 0 : d2xdxi2[2][4] = jacInv[0][2] * jacInv[2][2] * 2.0;
541 : 0 : d2xdxi2[2][5] = jacInv[1][2] * jacInv[2][2] * 2.0;
542 : :
543 : 0 : d2xdxi2[3][0] = jacInv[0][0] * jacInv[0][1];
544 : 0 : d2xdxi2[3][1] = jacInv[1][0] * jacInv[1][1];
545 : 0 : d2xdxi2[3][2] = jacInv[2][0] * jacInv[2][1];
546 : 0 : d2xdxi2[3][3] = jacInv[0][0] * jacInv[1][1] + jacInv[1][0] * jacInv[0][1];
547 : 0 : d2xdxi2[3][4] = jacInv[0][0] * jacInv[2][1] + jacInv[2][0] * jacInv[0][1];
548 : 0 : d2xdxi2[3][5] = jacInv[1][0] * jacInv[2][1] + jacInv[2][0] * jacInv[1][1];
549 : :
550 : 0 : d2xdxi2[4][0] = jacInv[0][0] * jacInv[0][2];
551 : 0 : d2xdxi2[4][1] = jacInv[1][0] * jacInv[1][2];
552 : 0 : d2xdxi2[4][2] = jacInv[2][0] * jacInv[2][2];
553 : 0 : d2xdxi2[4][3] = jacInv[0][0] * jacInv[1][2] + jacInv[1][0] * jacInv[0][2];
554 : 0 : d2xdxi2[4][4] = jacInv[0][0] * jacInv[2][2] + jacInv[2][0] * jacInv[0][2];
555 : 0 : d2xdxi2[4][5] = jacInv[1][0] * jacInv[2][2] + jacInv[2][0] * jacInv[1][2];
556 : :
557 : 0 : d2xdxi2[5][0] = jacInv[0][1] * jacInv[0][2];
558 : 0 : d2xdxi2[5][1] = jacInv[1][1] * jacInv[1][2];
559 : 0 : d2xdxi2[5][2] = jacInv[2][1] * jacInv[2][2];
560 : 0 : d2xdxi2[5][3] = jacInv[0][1] * jacInv[1][2] + jacInv[1][1] * jacInv[0][2];
561 : 0 : d2xdxi2[5][4] = jacInv[0][1] * jacInv[2][2] + jacInv[2][1] * jacInv[0][2];
562 : 0 : d2xdxi2[5][5] = jacInv[1][1] * jacInv[2][2] + jacInv[2][1] * jacInv[1][2];
563 : :
564 : : // Matrix to store the second order derivatives of basis functions in
565 : : // physical domain
566 : : tk::real d2Bdx2[6][6];
567 [ - - ]: 0 : for(std::size_t ibasis = 0; ibasis < 6; ibasis++) {
568 [ - - ]: 0 : for(std::size_t idir = 0; idir < 6; idir++) {
569 : 0 : d2Bdx2[idir][ibasis] = 0;
570 [ - - ]: 0 : for(std::size_t k = 0; k < 6; k++)
571 : 0 : d2Bdx2[idir][ibasis] += d2xdxi2[idir][k] * d2Bdxi2[k][ibasis];
572 : : }
573 : : }
574 : :
575 [ - - ]: 0 : for(ncomp_t icomp = 0; icomp < ncomp; icomp++)
576 : : {
577 : 0 : auto mark = icomp * ndof;
578 [ - - ]: 0 : for(std::size_t idir = 0; idir < 6; idir++)
579 : : {
580 : 0 : unk[icomp][idir+4] = 0;
581 [ - - ]: 0 : for(std::size_t ibasis = 0; ibasis < 6; ibasis++)
582 : 0 : unk[icomp][idir+4] += U(e, mark+4+ibasis) * d2Bdx2[idir][ibasis];
583 : : }
584 : : }
585 : : }
586 : 0 : return unk;
587 : : }
588 : :
589 : : void
590 : 0 : tk::TaylorToDubiner( ncomp_t ncomp,
591 : : std::size_t e,
592 : : std::size_t ndof,
593 : : const std::vector< std::size_t >& inpoel,
594 : : const tk::UnsMesh::Coords& coord,
595 : : const tk::Fields& geoElem,
596 : : std::vector< std::vector< tk::real > >& unk )
597 : : // *****************************************************************************
598 : : // Convert the solution with Taylor basis to the solution with Dubiner basis by
599 : : // projection method
600 : : //! \param[in] ncomp Number of scalar components in this PDE system
601 : : //! \param[in] e Id of element whose solution is to be limited
602 : : //! \param[in] ndof Maximum number of degrees of freedom
603 : : //! \param[in] inpoel Element connectivity
604 : : //! \param[in] coord Array of nodal coordinates
605 : : //! \param[in, out] unk High-order solution vector with Taylor basis
606 : : // *****************************************************************************
607 : : {
608 : : Assert( ncomp > 0, "Number of scalar components is incorrect" );
609 : :
610 : : // The diagonal of mass matrix
611 : 0 : std::vector< tk::real > L(ndof, 0.0);
612 : :
613 : : tk::real vol = 1.0 / 6.0;
614 : :
615 : 0 : L[0] = vol;
616 : :
617 [ - - ]: 0 : if(ndof > 1) {
618 : : Assert( (ndof == 4)||(ndof == 10),
619 : : "Mismatch in number of degrees of freedom" );
620 : 0 : L[1] = vol / 10.0;
621 : 0 : L[2] = vol * 3.0/10.0;
622 : 0 : L[3] = vol * 3.0/5.0;
623 : : }
624 : :
625 [ - - ]: 0 : if(ndof > 4) {
626 : : Assert( ndof == 10, "Mismatch in number of degrees of freedom" );
627 : 0 : L[4] = vol / 35.0;
628 : 0 : L[5] = vol / 21.0;
629 : 0 : L[6] = vol / 14.0;
630 : 0 : L[7] = vol / 7.0;
631 : 0 : L[8] = vol * 3.0/14.0;
632 : 0 : L[9] = vol * 3.0/7.0;
633 : : }
634 : :
635 : : // Coordinates of the centroid in physical domain
636 : 0 : std::array< tk::real, 3 > x_c{geoElem(e,1), geoElem(e,2), geoElem(e,3)};
637 : :
638 : : const auto& cx = coord[0];
639 : : const auto& cy = coord[1];
640 : : const auto& cz = coord[2];
641 : :
642 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
643 : 0 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
644 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
645 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
646 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
647 : 0 : }};
648 : :
649 : : // Number of quadrature points for volume integration
650 [ - - ]: 0 : auto ng = tk::NGvol(ndof);
651 : :
652 : : // arrays for quadrature points
653 : : std::array< std::vector< tk::real >, 3 > coordgp;
654 : : std::vector< tk::real > wgp;
655 : :
656 [ - - ]: 0 : coordgp[0].resize( ng );
657 [ - - ]: 0 : coordgp[1].resize( ng );
658 [ - - ]: 0 : coordgp[2].resize( ng );
659 [ - - ]: 0 : wgp.resize( ng );
660 : :
661 : : // get quadrature point weights and coordinates for triangle
662 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
663 : :
664 : : // right hand side vector
665 [ - - ][ - - ]: 0 : std::vector< tk::real > R( ncomp*ndof, 0.0 );
666 : :
667 : : // Gaussian quadrature
668 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
669 : : {
670 [ - - ]: 0 : auto wt = wgp[igp] * vol;
671 : :
672 [ - - ]: 0 : auto gp = tk::eval_gp( igp, coordel, coordgp );
673 : :
674 [ - - ]: 0 : auto B_taylor = eval_TaylorBasis( ndof, gp, x_c, coordel);
675 : :
676 : : // Compute high order solution at gauss point
677 [ - - ][ - - ]: 0 : std::vector< tk::real > state( ncomp, 0.0 );
678 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
679 : : {
680 [ - - ]: 0 : state[c] = unk[c][0];
681 : 0 : state[c] += unk[c][1] * B_taylor[1]
682 : 0 : + unk[c][2] * B_taylor[2]
683 : 0 : + unk[c][3] * B_taylor[3];
684 : :
685 [ - - ]: 0 : if(ndof > 4)
686 : 0 : state[c] += unk[c][4] * B_taylor[4] + unk[c][5] * B_taylor[5]
687 : 0 : + unk[c][6] * B_taylor[6] + unk[c][7] * B_taylor[7]
688 : 0 : + unk[c][8] * B_taylor[8] + unk[c][9] * B_taylor[9];
689 : : }
690 : :
691 [ - - ]: 0 : auto B = tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
692 : :
693 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
694 : : {
695 : 0 : auto mark = c*ndof;
696 [ - - ]: 0 : R[mark] += wt * state[c];
697 : :
698 [ - - ]: 0 : if(ndof > 1)
699 : : {
700 [ - - ]: 0 : R[mark+1] += wt * state[c] * B[1];
701 : 0 : R[mark+2] += wt * state[c] * B[2];
702 : 0 : R[mark+3] += wt * state[c] * B[3];
703 : :
704 [ - - ]: 0 : if(ndof > 4)
705 : : {
706 : 0 : R[mark+4] += wt * state[c] * B[4];
707 : 0 : R[mark+5] += wt * state[c] * B[5];
708 : 0 : R[mark+6] += wt * state[c] * B[6];
709 : 0 : R[mark+7] += wt * state[c] * B[7];
710 : 0 : R[mark+8] += wt * state[c] * B[8];
711 : 0 : R[mark+9] += wt * state[c] * B[9];
712 : : }
713 : : }
714 : : }
715 : : }
716 : :
717 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
718 : : {
719 : 0 : auto mark = c*ndof;
720 [ - - ]: 0 : for(std::size_t idof = 0; idof < ndof; idof++)
721 : 0 : unk[c][idof] = R[mark+idof] / L[idof];
722 : : }
723 : 0 : }
724 : :
725 : : std::vector< tk::real >
726 : 0 : tk::eval_TaylorBasis( const std::size_t ndof,
727 : : const std::array< tk::real, 3 >& x,
728 : : const std::array< tk::real, 3 >& x_c,
729 : : const std::array< std::array< tk::real, 3>, 4 >& coordel )
730 : : // *****************************************************************************
731 : : // Evaluate the Taylor basis at points
732 : : //! \param[in] ndof Maximum number of degrees of freedom
733 : : //! \param[in] x Nodal coordinates
734 : : //! \param[in] x_c Coordinates of the centroid
735 : : //! \param[in] coordel Array of nodal coordinates for the tetrahedron
736 : : // *****************************************************************************
737 : : {
738 : 0 : std::vector< tk::real > avg( 6, 0.0 );
739 [ - - ]: 0 : if(ndof > 4)
740 : : {
741 : : Assert( ndof == 10, "Mismatch in number of degrees of freedom" );
742 [ - - ]: 0 : auto ng = tk::NGvol(ndof);
743 : :
744 : : std::array< std::vector< tk::real >, 3 > coordgp;
745 : : std::vector< tk::real > wgp;
746 : :
747 [ - - ]: 0 : coordgp[0].resize( ng );
748 [ - - ]: 0 : coordgp[1].resize( ng );
749 [ - - ]: 0 : coordgp[2].resize( ng );
750 [ - - ]: 0 : wgp.resize( ng );
751 : :
752 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
753 : :
754 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
755 : : {
756 : : // Compute the coordinates of quadrature point at physical domain
757 [ - - ]: 0 : auto gp = tk::eval_gp( igp, coordel, coordgp );
758 : :
759 : 0 : avg[0] += wgp[igp] * (gp[0] - x_c[0]) * (gp[0] - x_c[0]) * 0.5;
760 : 0 : avg[1] += wgp[igp] * (gp[1] - x_c[1]) * (gp[1] - x_c[1]) * 0.5;
761 : 0 : avg[2] += wgp[igp] * (gp[2] - x_c[2]) * (gp[2] - x_c[2]) * 0.5;
762 : 0 : avg[3] += wgp[igp] * (gp[0] - x_c[0]) * (gp[1] - x_c[1]);
763 : 0 : avg[4] += wgp[igp] * (gp[0] - x_c[0]) * (gp[2] - x_c[2]);
764 : 0 : avg[5] += wgp[igp] * (gp[1] - x_c[1]) * (gp[2] - x_c[2]);
765 : : }
766 : : }
767 : :
768 [ - - ][ - - ]: 0 : std::vector< tk::real > B( ndof, 1.0 );
769 : :
770 [ - - ]: 0 : if(ndof > 1) {
771 : : Assert( (ndof == 4)||(ndof == 10) ,
772 : : "Mismatch in number of degrees of freedom" );
773 : 0 : B[1] = x[0] - x_c[0];
774 : 0 : B[2] = x[1] - x_c[1];
775 : 0 : B[3] = x[2] - x_c[2];
776 : : }
777 : :
778 [ - - ]: 0 : if(ndof > 4) {
779 : 0 : B[4] = B[1] * B[1] * 0.5 - avg[0];
780 : 0 : B[5] = B[2] * B[2] * 0.5 - avg[1];
781 : 0 : B[6] = B[3] * B[3] * 0.5 - avg[2];
782 : 0 : B[7] = B[1] * B[2] - avg[3];
783 : 0 : B[8] = B[1] * B[3] - avg[4];
784 : 0 : B[9] = B[2] * B[3] - avg[5];
785 : : }
786 : :
787 : 0 : return B;
788 : : }
789 : :
790 : : // -----------------------------------------------------------------------------
791 : : // Functions for reference element Taylor basis and related Xforms
792 : : // -----------------------------------------------------------------------------
793 : :
794 : : std::vector< std::vector< tk::real > >
795 : 0 : tk::DubinerToTaylorRefEl( ncomp_t ncomp,
796 : : const std::size_t e,
797 : : const std::size_t ndof,
798 : : const std::size_t ndof_el,
799 : : const std::vector< std::vector< tk::real > >& mtInv,
800 : : const tk::Fields& U )
801 : : // *****************************************************************************
802 : : // Transform the solution from Dubiner basis to Taylor basis
803 : : //! \param[in] ncomp Number of scalar components in this PDE system
804 : : //! \param[in] e Id of element whose solution is to be limited
805 : : //! \param[in] ndof Maximum number of degrees of freedom
806 : : //! \param[in] ndof_el Local number of degrees of freedom for the element
807 : : //! \param[in] mtInv Inverse of Taylor mass matrix
808 : : //! \param[in] U High-order solution vector with Dubiner basis
809 : : //! \return High-order solution vector with Taylor basis (ref element)
810 : : // *****************************************************************************
811 : : {
812 : : auto vol = 1.0/6.0;
813 : :
814 : : // 1. Get rhs for L2-projection
815 : : // Quadrature setup
816 : 0 : auto ng = tk::NGvol(ndof_el);
817 : : std::array< std::vector< real >, 3 > coordgp;
818 : : std::vector< real > wgp;
819 [ - - ]: 0 : coordgp[0].resize( ng );
820 [ - - ]: 0 : coordgp[1].resize( ng );
821 [ - - ]: 0 : coordgp[2].resize( ng );
822 [ - - ]: 0 : wgp.resize( ng );
823 [ - - ]: 0 : GaussQuadratureTet( ng, coordgp, wgp );
824 : :
825 : : // Gaussian quadrature
826 : : std::vector< std::vector< tk::real > >
827 [ - - ][ - - ]: 0 : R(ncomp, std::vector<tk::real>(ndof_el, 0.0));
[ - - ]
828 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
829 : : {
830 : : // Dubiner basis functions
831 [ - - ]: 0 : auto B = eval_basis( ndof_el, coordgp[0][igp], coordgp[1][igp],
832 [ - - ]: 0 : coordgp[2][igp] );
833 : : // Taylor basis functions
834 [ - - ]: 0 : auto Bt = eval_TaylorBasisRefEl(ndof_el, coordgp[0][igp], coordgp[1][igp],
835 [ - - ]: 0 : coordgp[2][igp]);
836 : :
837 [ - - ]: 0 : auto state = tk::eval_state(ncomp, ndof, ndof_el, e, U, B);
838 : :
839 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
840 [ - - ]: 0 : for (std::size_t id=0; id<ndof_el; ++id) {
841 : 0 : R[c][id] += wgp[igp] * vol * state[c] * Bt[id];
842 : : }
843 : : }
844 : : }
845 : :
846 : :
847 : : // 2. Get Taylor solution by premultiplying by mass matrix inverse
848 : : std::vector< std::vector< tk::real > >
849 [ - - ][ - - ]: 0 : unk(ncomp, std::vector<tk::real>(ndof_el, 0.0));
850 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
851 [ - - ]: 0 : for (std::size_t id=0; id<ndof_el; ++id) {
852 [ - - ]: 0 : for (std::size_t jd=0; jd<ndof_el; ++jd) {
853 : 0 : unk[c][id] += mtInv[id][jd] * R[c][jd];
854 : : }
855 : : }
856 : : }
857 : :
858 : 0 : return unk;
859 : : }
860 : :
861 : : void
862 : 0 : tk::TaylorToDubinerRefEl( ncomp_t ncomp,
863 : : const std::size_t ndof,
864 : : std::vector< std::vector< tk::real > >& unk )
865 : : // *****************************************************************************
866 : : // Transform the solution from Taylor to Dubiner basis
867 : : //! \param[in] ncomp Number of scalar components in this PDE system
868 : : //! \param[in] ndof Number of degrees of freedom
869 : : //! \param[in,out] unk High-order solution vector with Taylor basis that gets
870 : : //! transformed to solution with Dubiner basis
871 : : // *****************************************************************************
872 : : {
873 : : auto vol = 1.0/6.0;
874 : :
875 : 0 : auto M = massMatrixDubiner(ndof, vol);
876 : :
877 : : // 1. Get rhs for L2-projection
878 : : // Quadrature setup
879 [ - - ]: 0 : auto ng = tk::NGvol(ndof);
880 : : std::array< std::vector< real >, 3 > coordgp;
881 : : std::vector< real > wgp;
882 [ - - ]: 0 : coordgp[0].resize( ng );
883 [ - - ]: 0 : coordgp[1].resize( ng );
884 [ - - ]: 0 : coordgp[2].resize( ng );
885 [ - - ]: 0 : wgp.resize( ng );
886 [ - - ]: 0 : GaussQuadratureTet( ng, coordgp, wgp );
887 : :
888 : : // Gaussian quadrature
889 : : std::vector< std::vector< tk::real > >
890 [ - - ][ - - ]: 0 : R(ncomp, std::vector<tk::real>(ndof, 0.0));
[ - - ]
891 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
892 : : {
893 : : // Dubiner basis functions
894 [ - - ]: 0 : auto B = eval_basis( ndof, coordgp[0][igp], coordgp[1][igp],
895 [ - - ]: 0 : coordgp[2][igp] );
896 : : // Taylor basis functions
897 [ - - ]: 0 : auto Bt = eval_TaylorBasisRefEl(ndof, coordgp[0][igp], coordgp[1][igp],
898 [ - - ]: 0 : coordgp[2][igp]);
899 : :
900 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
901 : : real state(0.0);
902 [ - - ]: 0 : for (std::size_t id=0; id<ndof; ++id) {
903 : 0 : state += unk[c][id] * Bt[id];
904 : : }
905 [ - - ]: 0 : for (std::size_t id=0; id<ndof; ++id) {
906 : 0 : R[c][id] += wgp[igp] * vol * state * B[id];
907 : : }
908 : : }
909 : : }
910 : :
911 : : // 2. Get Dubiner solution by premultiplying by mass matrix inverse
912 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
913 [ - - ]: 0 : for (std::size_t id=0; id<ndof; ++id) {
914 : 0 : unk[c][id] = R[c][id] / M[id];
915 : : }
916 : : }
917 : 0 : }
918 : :
919 : : std::vector< tk::real >
920 : 2844 : tk::eval_TaylorBasisRefEl( std::size_t ndof, tk::real x, tk::real y,
921 : : tk::real z )
922 : : // *****************************************************************************
923 : : // Evaluate the Taylor basis at a point in the reference element
924 : : //! \param[in] ndof Number of degrees of freedom
925 : : //! \param[in] x Xi coordinate of point in reference element
926 : : //! \param[in] y Eta coordinate of point in reference element
927 : : //! \param[in] z Zeta coordinate of point in reference element
928 : : // *****************************************************************************
929 : : {
930 : : // Get averages required for P2 basis functions
931 : 2844 : std::vector< tk::real > avg( 6, 0.0 );
932 [ + + ]: 2844 : if(ndof > 4)
933 : : {
934 [ + - ]: 1694 : auto ng = tk::NGvol(ndof);
935 : : std::array< std::vector< tk::real >, 3 > coordgp;
936 : : std::vector< tk::real > wgp;
937 [ + - ]: 1694 : coordgp[0].resize( ng );
938 [ + - ]: 1694 : coordgp[1].resize( ng );
939 [ + - ]: 1694 : coordgp[2].resize( ng );
940 [ + - ]: 1694 : wgp.resize( ng );
941 [ + - ]: 1694 : tk::GaussQuadratureTet( ng, coordgp, wgp );
942 : :
943 [ + + ]: 20328 : for (std::size_t igp=0; igp<ng; ++igp)
944 : : {
945 : 18634 : avg[0] += wgp[igp] * (coordgp[0][igp] - 0.25) * (coordgp[0][igp] - 0.25) * 0.5;
946 : 18634 : avg[1] += wgp[igp] * (coordgp[1][igp] - 0.25) * (coordgp[1][igp] - 0.25) * 0.5;
947 : 18634 : avg[2] += wgp[igp] * (coordgp[2][igp] - 0.25) * (coordgp[2][igp] - 0.25) * 0.5;
948 : 18634 : avg[3] += wgp[igp] * (coordgp[0][igp] - 0.25) * (coordgp[1][igp] - 0.25);
949 : 18634 : avg[4] += wgp[igp] * (coordgp[0][igp] - 0.25) * (coordgp[2][igp] - 0.25);
950 : 18634 : avg[5] += wgp[igp] * (coordgp[1][igp] - 0.25) * (coordgp[2][igp] - 0.25);
951 : : }
952 : : }
953 : :
954 : : // Get Taylor basis functions
955 [ + - ][ - - ]: 2844 : std::vector< tk::real > B( ndof, 1.0 );
956 [ + - ]: 2844 : if(ndof > 1) {
957 : 2844 : B[1] = x - 0.25;
958 : 2844 : B[2] = y - 0.25;
959 : 2844 : B[3] = z - 0.25;
960 [ + + ]: 2844 : if(ndof > 4) {
961 : 1694 : B[4] = B[1] * B[1] * 0.5 - avg[0];
962 : 1694 : B[5] = B[2] * B[2] * 0.5 - avg[1];
963 : 1694 : B[6] = B[3] * B[3] * 0.5 - avg[2];
964 : 1694 : B[7] = B[1] * B[2] - avg[3];
965 : 1694 : B[8] = B[1] * B[3] - avg[4];
966 : 1694 : B[9] = B[2] * B[3] - avg[5];
967 : : }
968 : : }
969 : :
970 : 2844 : return B;
971 : : }
972 : :
973 : : std::vector< std::vector< tk::real > >
974 : 801 : tk::invMassMatTaylorRefEl( std::size_t dof )
975 : : // *****************************************************************************
976 : : // Obtain inverse mass matrix for Taylor basis in reference element
977 : : //! \param[in] dof Number of degrees of freedom
978 : : //! \return Inverse mass matrix
979 : : // *****************************************************************************
980 : : {
981 : : // Get Taylor mass matrix
982 : 801 : auto Mt = massMatrixTaylorRefEl(dof);
983 : :
984 : : // Only invert if DGP2
985 [ + + ]: 801 : if (dof > 4) {
986 : : double mtInv[10*10];
987 [ + + ]: 1694 : for (std::size_t i=0; i<Mt.size(); ++i) {
988 [ + + ]: 16940 : for (std::size_t j=0; j<Mt[i].size(); ++j) {
989 : 15400 : std::size_t idx = 10*i+j;
990 : 15400 : mtInv[idx] = Mt[i][j];
991 : : }
992 : : }
993 : : lapack_int ipiv[10];
994 : : // LU-factorization for inversion
995 [ + - ]: 154 : lapack_int info1 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, 10, 10, mtInv, 10, ipiv);
996 [ - + ][ - - ]: 154 : if (info1 != 0) Throw("Taylor mass matrix is singular");
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
997 : : // Inversion
998 [ + - ]: 154 : lapack_int info2 = LAPACKE_dgetri(LAPACK_ROW_MAJOR, 10, mtInv, 10, ipiv);
999 [ + - ][ - - ]: 154 : if (info2 != 0) Throw("Error while inverting Taylor mass matrix");
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1000 : :
1001 : : // Get 2D vector from 1D array mass matrix inverse
1002 [ + + ]: 1694 : for (std::size_t i=0; i<Mt.size(); ++i) {
1003 [ + + ]: 16940 : for (std::size_t j=0; j<Mt[i].size(); ++j) {
1004 : 15400 : std::size_t idx = 10*i+j;
1005 : 15400 : Mt[i][j] = mtInv[idx];
1006 : : }
1007 : : }
1008 : : }
1009 : :
1010 : 801 : return Mt;
1011 : : }
1012 : :
1013 : : std::vector< std::vector< tk::real > >
1014 : 801 : tk::massMatrixTaylorRefEl(std::size_t dof)
1015 : : // *****************************************************************************
1016 : : // Obtain mass matrix for Taylor basis in reference element
1017 : : //! \param[in] dof Number of degrees of freedom
1018 : : //! \return Mass matrix
1019 : : // *****************************************************************************
1020 : : {
1021 : : std::vector< std::vector< tk::real > >
1022 [ + - ][ + - ]: 801 : Mt(dof, std::vector<tk::real>(dof,0.0));
1023 : :
1024 : : // Mt(1,1)
1025 : : tk::real vol = 1.0/6.0;
1026 : 801 : Mt[0][0] = vol;
1027 : :
1028 : : // Mt(i,j) for i,j > 1
1029 [ + + ]: 801 : if (dof > 1) {
1030 : : // Quadrature information
1031 [ + - ]: 384 : auto ng = tk::NGvol(dof);
1032 : : std::array< std::vector< tk::real >, 3 > coordgp;
1033 : : std::vector< tk::real > wgp;
1034 [ + - ]: 384 : coordgp[0].resize( ng );
1035 [ + - ]: 384 : coordgp[1].resize( ng );
1036 [ + - ]: 384 : coordgp[2].resize( ng );
1037 [ + - ]: 384 : wgp.resize( ng );
1038 [ + - ]: 384 : tk::GaussQuadratureTet( ng, coordgp, wgp );
1039 : :
1040 [ + + ]: 3228 : for (std::size_t igp=0; igp<ng; ++igp)
1041 : : {
1042 [ + - ]: 2844 : auto Bt = eval_TaylorBasisRefEl(dof, coordgp[0][igp], coordgp[1][igp],
1043 [ + - ][ + - ]: 5688 : coordgp[2][igp]);
1044 [ + + ]: 21540 : for (std::size_t id=1; id<dof; ++id) {
1045 [ + + ]: 166260 : for (std::size_t jd=1; jd<dof; ++jd) {
1046 : 147564 : Mt[id][jd] += vol*wgp[igp]*Bt[id]*Bt[jd];
1047 : : }
1048 : : }
1049 : : }
1050 : : }
1051 : :
1052 : 801 : return Mt;
1053 : : }
|
