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 : :
22 : : #include "Basis.hpp"
23 : :
24 : : std::array< tk::real, 3 >
25 : 176152853 : tk::eval_gp ( const std::size_t igp,
26 : : const std::array< std::array< tk::real, 3>, 3 >& coordfa,
27 : : const std::array< std::vector< tk::real >, 2 >& coordgp )
28 : : // *****************************************************************************
29 : : // Compute the coordinates of quadrature points for face integral in physical
30 : : // space
31 : : //! \param[in] igp Index of quadrature points
32 : : //! \param[in] coordfa Array of nodal coordinates for face element
33 : : //! \param[in] coordgp Array of coordinates for quadrature points in reference
34 : : //! space
35 : : //! \return Array of coordinates for quadrature points in physical space
36 : : // *****************************************************************************
37 : : {
38 : : // Barycentric coordinates for the triangular face
39 : 176152853 : auto shp1 = 1.0 - coordgp[0][igp] - coordgp[1][igp];
40 : : auto shp2 = coordgp[0][igp];
41 : : auto shp3 = coordgp[1][igp];
42 : :
43 : : // Transformation of the quadrature point from the 2D reference/master
44 : : // element to physical space, to obtain its physical (x,y,z) coordinates.
45 : 176152853 : return {{ coordfa[0][0]*shp1 + coordfa[1][0]*shp2 + coordfa[2][0]*shp3,
46 : 176152853 : coordfa[0][1]*shp1 + coordfa[1][1]*shp2 + coordfa[2][1]*shp3,
47 : 176152853 : coordfa[0][2]*shp1 + coordfa[1][2]*shp2 + coordfa[2][2]*shp3 }};
48 : : }
49 : :
50 : : std::array< tk::real, 3 >
51 : 78647977 : tk::eval_gp ( const std::size_t igp,
52 : : const std::array< std::array< tk::real, 3>, 4 >& coord,
53 : : const std::array< std::vector< tk::real >, 3 >& coordgp )
54 : : // *****************************************************************************
55 : : // Compute the coordinates of quadrature points for volume integral in
56 : : // physical space
57 : : //! \param[in] igp Index of quadrature points
58 : : //! \param[in] coord Array of nodal coordinates for tetrahedron element
59 : : //! \param[in] coordgp Array of coordinates for quadrature points in reference space
60 : : //! \return Array of coordinates for quadrature points in physical space
61 : : // *****************************************************************************
62 : : {
63 : : // Barycentric coordinates for the tetradedron element
64 : 78647977 : auto shp1 = 1.0 - coordgp[0][igp] - coordgp[1][igp] - coordgp[2][igp];
65 : : auto shp2 = coordgp[0][igp];
66 : : auto shp3 = coordgp[1][igp];
67 : : auto shp4 = coordgp[2][igp];
68 : :
69 : : // Transformation of the quadrature point from the reference/master
70 : : // element to physical space, to obtain its physical (x,y,z) coordinates.
71 : : return {{
72 : 78647977 : coord[0][0]*shp1 + coord[1][0]*shp2 + coord[2][0]*shp3 + coord[3][0]*shp4,
73 : 78647977 : coord[0][1]*shp1 + coord[1][1]*shp2 + coord[2][1]*shp3 + coord[3][1]*shp4,
74 : 78647977 : coord[0][2]*shp1 + coord[1][2]*shp2 + coord[2][2]*shp3 + coord[3][2]*shp4 }};
75 : : }
76 : :
77 : : std::array< std::vector<tk::real>, 3 >
78 : 12109623 : tk::eval_dBdx_p1( const std::size_t ndof,
79 : : const std::array< std::array< tk::real, 3 >, 3 >& jacInv )
80 : : // *****************************************************************************
81 : : // Compute the derivatives of basis functions for DG(P1)
82 : : //! \param[in] ndof Number of degrees of freedom
83 : : //! \param[in] jacInv Array of the inverse of Jacobian
84 : : //! \return Array of the derivatives of basis functions
85 : : // *****************************************************************************
86 : : {
87 : : // The derivatives of the basis functions dB/dx are easily calculated
88 : : // via a transformation to the reference space as,
89 : : // dB/dx = dB/dxi . dxi/dx,
90 : : // where, x = (x,y,z) are the physical coordinates, and
91 : : // xi = (xi, eta, zeta) are the reference coordinates.
92 : : // The matrix dxi/dx is the inverse of the Jacobian of transformation
93 : : // and the matrix vector product has to be calculated. This follows.
94 : :
95 : : std::array< std::vector<tk::real>, 3 > dBdx;
96 [ + - ][ + - ]: 12109623 : dBdx[0].resize( ndof, 0 );
97 [ + - ][ + - ]: 12109623 : dBdx[1].resize( ndof, 0 );
98 [ + - ]: 12109623 : dBdx[2].resize( ndof, 0 );
99 : :
100 : : auto db2dxi1 = 2.0;
101 : : auto db2dxi2 = 1.0;
102 : : auto db2dxi3 = 1.0;
103 : :
104 : : auto db3dxi1 = 0.0;
105 : : auto db3dxi2 = 3.0;
106 : : auto db3dxi3 = 1.0;
107 : :
108 : : auto db4dxi1 = 0.0;
109 : : auto db4dxi2 = 0.0;
110 : : auto db4dxi3 = 4.0;
111 : :
112 : 12109623 : dBdx[0][1] = db2dxi1 * jacInv[0][0]
113 : 12109623 : + db2dxi2 * jacInv[1][0]
114 : 12109623 : + db2dxi3 * jacInv[2][0];
115 : :
116 : 12109623 : dBdx[1][1] = db2dxi1 * jacInv[0][1]
117 : 12109623 : + db2dxi2 * jacInv[1][1]
118 : 12109623 : + db2dxi3 * jacInv[2][1];
119 : :
120 : 12109623 : dBdx[2][1] = db2dxi1 * jacInv[0][2]
121 : 12109623 : + db2dxi2 * jacInv[1][2]
122 : 12109623 : + db2dxi3 * jacInv[2][2];
123 : :
124 : 12109623 : dBdx[0][2] = db3dxi1 * jacInv[0][0]
125 : 12109623 : + db3dxi2 * jacInv[1][0]
126 : 12109623 : + db3dxi3 * jacInv[2][0];
127 : :
128 : 12109623 : dBdx[1][2] = db3dxi1 * jacInv[0][1]
129 : 12109623 : + db3dxi2 * jacInv[1][1]
130 : 12109623 : + db3dxi3 * jacInv[2][1];
131 : :
132 : 12109623 : dBdx[2][2] = db3dxi1 * jacInv[0][2]
133 : 12109623 : + db3dxi2 * jacInv[1][2]
134 : 12109623 : + db3dxi3 * jacInv[2][2];
135 : :
136 : 12109623 : dBdx[0][3] = db4dxi1 * jacInv[0][0]
137 : 12109623 : + db4dxi2 * jacInv[1][0]
138 : 12109623 : + db4dxi3 * jacInv[2][0];
139 : :
140 : 12109623 : dBdx[1][3] = db4dxi1 * jacInv[0][1]
141 : 12109623 : + db4dxi2 * jacInv[1][1]
142 : 12109623 : + db4dxi3 * jacInv[2][1];
143 : :
144 : 12109623 : dBdx[2][3] = db4dxi1 * jacInv[0][2]
145 : 12109623 : + db4dxi2 * jacInv[1][2]
146 : 12109623 : + db4dxi3 * jacInv[2][2];
147 : :
148 : 12109623 : return dBdx;
149 : : }
150 : :
151 : : void
152 : 20852475 : tk::eval_dBdx_p2( const std::size_t igp,
153 : : const std::array< std::vector< tk::real >, 3 >& coordgp,
154 : : const std::array< std::array< tk::real, 3 >, 3 >& jacInv,
155 : : std::array< std::vector<tk::real>, 3 >& dBdx )
156 : : // *****************************************************************************
157 : : // Compute the derivatives of Dubiner basis function for DG(P2)
158 : : //! \param[in] igp Index of quadrature points
159 : : //! \param[in] coordgp Gauss point coordinates for tetrahedron element
160 : : //! \param[in] jacInv Array of the inverse of Jacobian
161 : : //! \param[in,out] dBdx Array of the derivatives of basis function
162 : : // *****************************************************************************
163 : : {
164 : 20852475 : auto db5dxi1 = 12.0 * coordgp[0][igp] + 6.0 * coordgp[1][igp]
165 : 20852475 : + 6.0 * coordgp[2][igp] - 6.0;
166 : 20852475 : auto db5dxi2 = 6.0 * coordgp[0][igp] + 2.0 * coordgp[1][igp]
167 : 20852475 : + 2.0 * coordgp[2][igp] - 2.0;
168 : : auto db5dxi3 = 6.0 * coordgp[0][igp] + 2.0 * coordgp[1][igp]
169 : : + 2.0 * coordgp[2][igp] - 2.0;
170 : :
171 : 20852475 : auto db6dxi1 = 10.0 * coordgp[1][igp] + 2.0 * coordgp[2][igp] - 2.0;
172 : 20852475 : auto db6dxi2 = 10.0 * coordgp[0][igp] + 10.0 * coordgp[1][igp]
173 : 20852475 : + 6.0 * coordgp[2][igp] - 6.0;
174 : 20852475 : auto db6dxi3 = 2.0 * coordgp[0][igp] + 6.0 * coordgp[1][igp]
175 : 20852475 : + 2.0 * coordgp[2][igp] - 2.0;
176 : :
177 : 20852475 : auto db7dxi1 = 12.0 * coordgp[2][igp] - 2.0;
178 : 20852475 : auto db7dxi2 = 6.0 * coordgp[2][igp] - 1.0;
179 : : auto db7dxi3 = 12.0 * coordgp[0][igp] + 6.0 * coordgp[1][igp]
180 : 20852475 : + 12.0 * coordgp[2][igp] - 7.0;
181 : :
182 : : auto db8dxi1 = 0;
183 : 20852475 : auto db8dxi2 = 20.0 * coordgp[1][igp] + 8.0 * coordgp[2][igp] - 8.0;
184 : 20852475 : auto db8dxi3 = 8.0 * coordgp[1][igp] + 2.0 * coordgp[2][igp] - 2.0;
185 : :
186 : : auto db9dxi1 = 0;
187 : 20852475 : auto db9dxi2 = 18.0 * coordgp[2][igp] - 3.0;
188 : 20852475 : auto db9dxi3 = 18.0 * coordgp[1][igp] + 12.0 * coordgp[2][igp] - 7.0;
189 : :
190 : : auto db10dxi1 = 0;
191 : : auto db10dxi2 = 0;
192 : 20852475 : auto db10dxi3 = 30.0 * coordgp[2][igp] - 10.0;
193 : :
194 : 20852475 : dBdx[0][4] = db5dxi1 * jacInv[0][0]
195 : 20852475 : + db5dxi2 * jacInv[1][0]
196 : 20852475 : + db5dxi3 * jacInv[2][0];
197 : :
198 : 20852475 : dBdx[1][4] = db5dxi1 * jacInv[0][1]
199 : 20852475 : + db5dxi2 * jacInv[1][1]
200 : 20852475 : + db5dxi3 * jacInv[2][1];
201 : :
202 : 20852475 : dBdx[2][4] = db5dxi1 * jacInv[0][2]
203 : 20852475 : + db5dxi2 * jacInv[1][2]
204 : 20852475 : + db5dxi3 * jacInv[2][2];
205 : :
206 : 20852475 : dBdx[0][5] = db6dxi1 * jacInv[0][0]
207 : 20852475 : + db6dxi2 * jacInv[1][0]
208 : 20852475 : + db6dxi3 * jacInv[2][0];
209 : :
210 : 20852475 : dBdx[1][5] = db6dxi1 * jacInv[0][1]
211 : 20852475 : + db6dxi2 * jacInv[1][1]
212 : 20852475 : + db6dxi3 * jacInv[2][1];
213 : :
214 : 20852475 : dBdx[2][5] = db6dxi1 * jacInv[0][2]
215 : 20852475 : + db6dxi2 * jacInv[1][2]
216 : 20852475 : + db6dxi3 * jacInv[2][2];
217 : :
218 : 20852475 : dBdx[0][6] = db7dxi1 * jacInv[0][0]
219 : 20852475 : + db7dxi2 * jacInv[1][0]
220 : 20852475 : + db7dxi3 * jacInv[2][0];
221 : :
222 : 20852475 : dBdx[1][6] = db7dxi1 * jacInv[0][1]
223 : 20852475 : + db7dxi2 * jacInv[1][1]
224 : 20852475 : + db7dxi3 * jacInv[2][1];
225 : :
226 : 20852475 : dBdx[2][6] = db7dxi1 * jacInv[0][2]
227 : 20852475 : + db7dxi2 * jacInv[1][2]
228 : 20852475 : + db7dxi3 * jacInv[2][2];
229 : :
230 : 20852475 : dBdx[0][7] = db8dxi1 * jacInv[0][0]
231 : 20852475 : + db8dxi2 * jacInv[1][0]
232 : 20852475 : + db8dxi3 * jacInv[2][0];
233 : :
234 : 20852475 : dBdx[1][7] = db8dxi1 * jacInv[0][1]
235 : 20852475 : + db8dxi2 * jacInv[1][1]
236 : 20852475 : + db8dxi3 * jacInv[2][1];
237 : :
238 : 20852475 : dBdx[2][7] = db8dxi1 * jacInv[0][2]
239 : 20852475 : + db8dxi2 * jacInv[1][2]
240 : 20852475 : + db8dxi3 * jacInv[2][2];
241 : :
242 : 20852475 : dBdx[0][8] = db9dxi1 * jacInv[0][0]
243 : 20852475 : + db9dxi2 * jacInv[1][0]
244 : 20852475 : + db9dxi3 * jacInv[2][0];
245 : :
246 : 20852475 : dBdx[1][8] = db9dxi1 * jacInv[0][1]
247 : 20852475 : + db9dxi2 * jacInv[1][1]
248 : 20852475 : + db9dxi3 * jacInv[2][1];
249 : :
250 : 20852475 : dBdx[2][8] = db9dxi1 * jacInv[0][2]
251 : 20852475 : + db9dxi2 * jacInv[1][2]
252 : 20852475 : + db9dxi3 * jacInv[2][2];
253 : :
254 : 20852475 : dBdx[0][9] = db10dxi1 * jacInv[0][0]
255 : 20852475 : + db10dxi2 * jacInv[1][0]
256 : 20852475 : + db10dxi3 * jacInv[2][0];
257 : :
258 : 20852475 : dBdx[1][9] = db10dxi1 * jacInv[0][1]
259 : 20852475 : + db10dxi2 * jacInv[1][1]
260 : 20852475 : + db10dxi3 * jacInv[2][1];
261 : :
262 : 20852475 : dBdx[2][9] = db10dxi1 * jacInv[0][2]
263 : 20852475 : + db10dxi2 * jacInv[1][2]
264 : 20852475 : + db10dxi3 * jacInv[2][2];
265 : 20852475 : }
266 : :
267 : : std::vector< tk::real >
268 : 382846871 : tk::eval_basis( const std::size_t ndof,
269 : : const tk::real xi,
270 : : const tk::real eta,
271 : : const tk::real zeta )
272 : : // *****************************************************************************
273 : : // Compute the Dubiner basis functions
274 : : //! \param[in] ndof Number of degrees of freedom
275 : : //! \param[in] xi,eta,zeta Coordinates for quadrature points in reference space
276 : : //! \return Vector of basis functions
277 : : // *****************************************************************************
278 : : {
279 : : // Array of basis functions
280 : 382846871 : std::vector< tk::real > B( ndof, 1.0 );
281 : :
282 [ + + ]: 382846871 : if ( ndof > 1 ) // DG(P1)
283 : : {
284 : 248971324 : B[1] = 2.0 * xi + eta + zeta - 1.0;
285 : 248971324 : B[2] = 3.0 * eta + zeta - 1.0;
286 : 248971324 : B[3] = 4.0 * zeta - 1.0;
287 : :
288 [ + + ]: 248971324 : if( ndof > 4 ) // DG(P2)
289 : : {
290 : 89726295 : B[4] = 6.0 * xi * xi + eta * eta + zeta * zeta
291 : 89726295 : + 6.0 * xi * eta + 6.0 * xi * zeta + 2.0 * eta * zeta
292 : 89726295 : - 6.0 * xi - 2.0 * eta - 2.0 * zeta + 1.0;
293 : 89726295 : B[5] = 5.0 * eta * eta + zeta * zeta
294 : 89726295 : + 10.0 * xi * eta + 2.0 * xi * zeta + 6.0 * eta * zeta
295 : 89726295 : - 2.0 * xi - 6.0 * eta - 2.0 * zeta + 1.0;
296 : 89726295 : B[6] = 6.0 * zeta * zeta + 12.0 * xi * zeta + 6.0 * eta * zeta - 2.0 * xi
297 : 89726295 : - eta - 7.0 * zeta + 1.0;
298 : 89726295 : B[7] = 10.0 * eta * eta + zeta * zeta + 8.0 * eta * zeta
299 : 89726295 : - 8.0 * eta - 2.0 * zeta + 1.0;
300 : 89726295 : B[8] = 6.0 * zeta * zeta + 18.0 * eta * zeta - 3.0 * eta - 7.0 * zeta
301 : 89726295 : + 1.0;
302 : 89726295 : B[9] = 15.0 * zeta * zeta - 10.0 * zeta + 1.0;
303 : : }
304 : : }
305 : :
306 : 382846871 : return B;
307 : : }
308 : :
309 : : std::vector< tk::real >
310 [ + + ]: 611291486 : tk::eval_state ( ncomp_t ncomp,
311 : : ncomp_t offset,
312 : : const std::size_t ndof,
313 : : const std::size_t ndof_el,
314 : : const std::size_t e,
315 : : const Fields& U,
316 : : const std::vector< tk::real >& B,
317 : : const std::array< std::size_t, 2 >& VarRange )
318 : : // *****************************************************************************
319 : : // Compute the state variables for the tetrahedron element
320 : : //! \param[in] ncomp Number of scalar components in this PDE system
321 : : //! \param[in] offset Offset this PDE system operates from
322 : : //! \param[in] ndof Maximum number of degrees of freedom
323 : : //! \param[in] ndof_el Number of degrees of freedom for the local element
324 : : //! \param[in] e Index for the tetrahedron element
325 : : //! \param[in] U Solution vector at recent time step
326 : : //! \param[in] B Vector of basis functions
327 : : //! \param[in] VarRange Range of the variables to be evaluated
328 : : //! \return Vector of state variable for tetrahedron element
329 : : // *****************************************************************************
330 : : {
331 : : // This is commented for now because that when p0/p1 adaptive with limiter
332 : : // applied, the size of basis will be 10. However, ndof_el will be 4 which
333 : : // leads to a size mismatch in limiter function.
334 : : //Assert( B.size() == ndof_el, "Size mismatch" );
335 : :
336 [ + + ]: 611291486 : if (U.empty()) return {};
337 : :
338 : : // Array of state variable for tetrahedron element
339 : 360719054 : std::vector< tk::real > state( ncomp, 0.0 );
340 : :
341 [ + + ]: 1465842178 : for (ncomp_t c=VarRange[0]; c<=VarRange[1]; ++c)
342 : : {
343 [ + + ]: 1105123124 : auto mark = c*ndof;
344 [ + + ]: 1105123124 : state[c] = U( e, mark, offset );
345 : :
346 [ + + ]: 1105123124 : if(ndof_el > 1) //DG(P1)
347 : : {
348 : 801647181 : state[c] += U( e, mark+1, offset ) * B[1]
349 : 801647181 : + U( e, mark+2, offset ) * B[2]
350 : 801647181 : + U( e, mark+3, offset ) * B[3];
351 : : }
352 : :
353 [ + + ]: 1105123124 : if(ndof_el > 4) //DG(P2)
354 : : {
355 : 254490903 : state[c] += U( e, mark+4, offset ) * B[4]
356 : 254490903 : + U( e, mark+5, offset ) * B[5]
357 : 254490903 : + U( e, mark+6, offset ) * B[6]
358 : 254490903 : + U( e, mark+7, offset ) * B[7]
359 : 254490903 : + U( e, mark+8, offset ) * B[8]
360 : 254490903 : + U( e, mark+9, offset ) * B[9];
361 : : }
362 : : }
363 : :
364 : : return state;
365 : : }
366 : :
367 : : std::vector< std::vector< tk::real > >
368 : 3604 : tk::DubinerToTaylor( ncomp_t ncomp,
369 : : ncomp_t offset,
370 : : const std::size_t e,
371 : : const std::size_t ndof,
372 : : const tk::Fields& U,
373 : : const std::vector< std::size_t >& inpoel,
374 : : const tk::UnsMesh::Coords& coord )
375 : : // *****************************************************************************
376 : : // Transform the solution with Dubiner basis to the solution with Taylor basis
377 : : //! \param[in] ncomp Number of scalar components in this PDE system
378 : : //! \param[in] offset Index for equation systems
379 : : //! \param[in] e Id of element whose solution is to be limited
380 : : //! \param[in] ndof Maximum number of degrees of freedom
381 : : //! \param[in] U High-order solution vector with Dubiner basis
382 : : //! \param[in] inpoel Element connectivity
383 : : //! \param[in] coord Array of nodal coordinates
384 : : //! \return High-order solution vector with Taylor basis
385 : : // *****************************************************************************
386 : : {
387 : : std::vector< std::vector< tk::real > >
388 [ + - ][ + - ]: 7208 : unk(ncomp, std::vector<tk::real>(ndof, 0.0));
389 : :
390 : : const auto& cx = coord[0];
391 : : const auto& cy = coord[1];
392 : : const auto& cz = coord[2];
393 : :
394 : : std::array< std::vector< tk::real >, 3 > center;
395 [ + - ]: 3604 : center[0].resize(1, 0.25);
396 [ + - ]: 3604 : center[1].resize(1, 0.25);
397 [ + - ]: 3604 : center[2].resize(1, 0.25);
398 : :
399 : : // Evaluate the cell center solution
400 [ + + ]: 21624 : for(ncomp_t icomp = 0; icomp < ncomp; icomp++)
401 : : {
402 : 18020 : auto mark = icomp * ndof;
403 : 18020 : unk[icomp][0] = U(e, mark, offset);
404 : : }
405 : :
406 : : // Evaluate the first order derivative
407 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
408 : 3604 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
409 : 3604 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
410 : 3604 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
411 : 3604 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
412 : 3604 : }};
413 : :
414 : : auto jacInv =
415 : 3604 : tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
416 : :
417 : : // Compute the derivatives of basis function for DG(P1)
418 [ + - ]: 3604 : auto dBdx = tk::eval_dBdx_p1( ndof, jacInv );
419 : :
420 [ + - ]: 3604 : if(ndof > 4) {
421 [ - + ]: 3604 : tk::eval_dBdx_p2(0, center, jacInv, dBdx);
422 : : }
423 : :
424 [ + + ]: 21624 : for(ncomp_t icomp = 0; icomp < ncomp; icomp++)
425 : : {
426 : 18020 : auto mark = icomp * ndof;
427 [ + + ]: 72080 : for(std::size_t idir = 0; idir < 3; idir++)
428 : : {
429 : 54060 : unk[icomp][idir+1] = 0;
430 [ + + ]: 540600 : for(std::size_t idof = 1; idof < ndof; idof++)
431 : 486540 : unk[icomp][idir+1] += U(e, mark+idof, offset) * dBdx[idir][idof];
432 : : }
433 : : }
434 : :
435 : : // Evaluate the second order derivative if DGP2 is applied
436 : : // The basic idea of the computation follows
437 : : // d2Udx2 = /sum u_i * (d2B_i/dx2)
438 : : // where d2B_i/dx2 = d( dB_i/dxi * dxi/dx ) / dxi * dxi/dx
439 [ + - ]: 3604 : if(ndof > 4)
440 : : {
441 : : // Matrix to store the second order derivatives of basis functions in
442 : : // reference domain
443 : 3604 : tk::real d2Bdxi2[6][6] =
444 : : { { 12.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
445 : : { 2.0, 10.0, 0.0, 20.0, 0.0, 0.0 },
446 : : { 2.0, 2.0, 12.0, 2.0, 12.0, 30.0 },
447 : : { 6.0, 10.0, 0.0, 0.0, 0.0, 0.0 },
448 : : { 6.0, 2.0, 12.0, 0.0, 0.0, 0.0 },
449 : : { 2.0, 6.0, 6.0, 8.0, 18.0, 0.0 } };
450 : :
451 : : // Transform matrix to convert the second order derivatives of basis
452 : : // function in reference domain to the one in physical domain
453 : : tk::real d2xdxi2[6][6];
454 : :
455 : 3604 : d2xdxi2[0][0] = jacInv[0][0] * jacInv[0][0];
456 : 3604 : d2xdxi2[0][1] = jacInv[1][0] * jacInv[1][0];
457 : 3604 : d2xdxi2[0][2] = jacInv[2][0] * jacInv[2][0];
458 : 3604 : d2xdxi2[0][3] = jacInv[0][0] * jacInv[1][0] * 2.0;
459 : 3604 : d2xdxi2[0][4] = jacInv[0][0] * jacInv[2][0] * 2.0;
460 : 3604 : d2xdxi2[0][5] = jacInv[1][0] * jacInv[2][0] * 2.0;
461 : :
462 : 3604 : d2xdxi2[1][0] = jacInv[0][1] * jacInv[0][1];
463 : 3604 : d2xdxi2[1][1] = jacInv[1][1] * jacInv[1][1];
464 : 3604 : d2xdxi2[1][2] = jacInv[2][1] * jacInv[2][1];
465 : 3604 : d2xdxi2[1][3] = jacInv[0][1] * jacInv[1][1] * 2.0;
466 : 3604 : d2xdxi2[1][4] = jacInv[0][1] * jacInv[2][1] * 2.0;
467 : 3604 : d2xdxi2[1][5] = jacInv[1][1] * jacInv[2][1] * 2.0;
468 : :
469 : 3604 : d2xdxi2[2][0] = jacInv[0][2] * jacInv[0][2];
470 : 3604 : d2xdxi2[2][1] = jacInv[1][2] * jacInv[1][2];
471 : 3604 : d2xdxi2[2][2] = jacInv[2][2] * jacInv[2][2];
472 : 3604 : d2xdxi2[2][3] = jacInv[0][2] * jacInv[1][2] * 2.0;
473 : 3604 : d2xdxi2[2][4] = jacInv[0][2] * jacInv[2][2] * 2.0;
474 : 3604 : d2xdxi2[2][5] = jacInv[1][2] * jacInv[2][2] * 2.0;
475 : :
476 : 3604 : d2xdxi2[3][0] = jacInv[0][0] * jacInv[0][1];
477 : 3604 : d2xdxi2[3][1] = jacInv[1][0] * jacInv[1][1];
478 : 3604 : d2xdxi2[3][2] = jacInv[2][0] * jacInv[2][1];
479 : 3604 : d2xdxi2[3][3] = jacInv[0][0] * jacInv[1][1] + jacInv[1][0] * jacInv[0][1];
480 : 3604 : d2xdxi2[3][4] = jacInv[0][0] * jacInv[2][1] + jacInv[2][0] * jacInv[0][1];
481 : 3604 : d2xdxi2[3][5] = jacInv[1][0] * jacInv[2][1] + jacInv[2][0] * jacInv[1][1];
482 : :
483 : 3604 : d2xdxi2[4][0] = jacInv[0][0] * jacInv[0][2];
484 : 3604 : d2xdxi2[4][1] = jacInv[1][0] * jacInv[1][2];
485 : 3604 : d2xdxi2[4][2] = jacInv[2][0] * jacInv[2][2];
486 : 3604 : d2xdxi2[4][3] = jacInv[0][0] * jacInv[1][2] + jacInv[1][0] * jacInv[0][2];
487 : 3604 : d2xdxi2[4][4] = jacInv[0][0] * jacInv[2][2] + jacInv[2][0] * jacInv[0][2];
488 : 3604 : d2xdxi2[4][5] = jacInv[1][0] * jacInv[2][2] + jacInv[2][0] * jacInv[1][2];
489 : :
490 : 3604 : d2xdxi2[5][0] = jacInv[0][1] * jacInv[0][2];
491 : 3604 : d2xdxi2[5][1] = jacInv[1][1] * jacInv[1][2];
492 : 3604 : d2xdxi2[5][2] = jacInv[2][1] * jacInv[2][2];
493 : 3604 : d2xdxi2[5][3] = jacInv[0][1] * jacInv[1][2] + jacInv[1][1] * jacInv[0][2];
494 : 3604 : d2xdxi2[5][4] = jacInv[0][1] * jacInv[2][2] + jacInv[2][1] * jacInv[0][2];
495 : 3604 : d2xdxi2[5][5] = jacInv[1][1] * jacInv[2][2] + jacInv[2][1] * jacInv[1][2];
496 : :
497 : : // Matrix to store the second order derivatives of basis functions in
498 : : // physical domain
499 : : tk::real d2Bdx2[6][6];
500 [ + + ]: 25228 : for(std::size_t ibasis = 0; ibasis < 6; ibasis++) {
501 [ + + ]: 151368 : for(std::size_t idir = 0; idir < 6; idir++) {
502 : 129744 : d2Bdx2[idir][ibasis] = 0;
503 [ + + ]: 908208 : for(std::size_t k = 0; k < 6; k++)
504 : 778464 : d2Bdx2[idir][ibasis] += d2xdxi2[idir][k] * d2Bdxi2[k][ibasis];
505 : : }
506 : : }
507 : :
508 [ + + ]: 21624 : for(ncomp_t icomp = 0; icomp < ncomp; icomp++)
509 : : {
510 : 18020 : auto mark = icomp * ndof;
511 [ + + ]: 126140 : for(std::size_t idir = 0; idir < 6; idir++)
512 : : {
513 : 108120 : unk[icomp][idir+4] = 0;
514 [ + + ]: 756840 : for(std::size_t ibasis = 0; ibasis < 6; ibasis++)
515 : 648720 : unk[icomp][idir+4] += U(e, mark+4+ibasis, offset) * d2Bdx2[idir][ibasis];
516 : : }
517 : : }
518 : : }
519 : 3604 : return unk;
520 : : }
521 : :
522 : : void
523 : 3604 : tk::TaylorToDubiner( ncomp_t ncomp,
524 : : std::size_t e,
525 : : std::size_t ndof,
526 : : const std::vector< std::size_t >& inpoel,
527 : : const tk::UnsMesh::Coords& coord,
528 : : const tk::Fields& geoElem,
529 : : std::vector< std::vector< tk::real > >& unk )
530 : : // *****************************************************************************
531 : : // Convert the solution with Taylor basis to the solution with Dubiner basis by
532 : : // projection method
533 : : //! \param[in] ncomp Number of scalar components in this PDE system
534 : : //! \param[in] e Id of element whose solution is to be limited
535 : : //! \param[in] ndof Maximum number of degrees of freedom
536 : : //! \param[in] inpoel Element connectivity
537 : : //! \param[in] coord Array of nodal coordinates
538 : : //! \param[in, out] unk High-order solution vector with Taylor basis
539 : : // *****************************************************************************
540 : : {
541 : : Assert( ncomp > 0, "Number of scalar components is incorrect" );
542 : :
543 : : // The diagonal of mass matrix
544 : 3604 : std::vector< tk::real > L(ndof, 0.0);
545 : :
546 : : tk::real vol = 1.0 / 6.0;
547 : :
548 : 3604 : L[0] = vol;
549 : :
550 [ + - ]: 3604 : if(ndof > 1) {
551 : : Assert( (ndof == 4)||(ndof == 10),
552 : : "Mismatch in number of degrees of freedom" );
553 : 3604 : L[1] = vol / 10.0;
554 : 3604 : L[2] = vol * 3.0/10.0;
555 : 3604 : L[3] = vol * 3.0/5.0;
556 : : }
557 : :
558 [ + - ]: 3604 : if(ndof > 4) {
559 : : Assert( ndof == 10, "Mismatch in number of degrees of freedom" );
560 : 3604 : L[4] = vol / 35.0;
561 : 3604 : L[5] = vol / 21.0;
562 : 3604 : L[6] = vol / 14.0;
563 : 3604 : L[7] = vol / 7.0;
564 : 3604 : L[8] = vol * 3.0/14.0;
565 : 3604 : L[9] = vol * 3.0/7.0;
566 : : }
567 : :
568 : : // Coordinates of the centroid in physical domain
569 : 3604 : std::array< tk::real, 3 > x_c{geoElem(e,1,0), geoElem(e,2,0), geoElem(e,3,0)};
570 : :
571 : : const auto& cx = coord[0];
572 : : const auto& cy = coord[1];
573 : : const auto& cz = coord[2];
574 : :
575 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
576 : 3604 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
577 : 3604 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
578 : 3604 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
579 : 3604 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
580 : 3604 : }};
581 : :
582 : : // Number of quadrature points for volume integration
583 [ + - ]: 3604 : auto ng = tk::NGvol(ndof);
584 : :
585 : : // arrays for quadrature points
586 : : std::array< std::vector< tk::real >, 3 > coordgp;
587 : : std::vector< tk::real > wgp;
588 : :
589 [ + - ]: 3604 : coordgp[0].resize( ng );
590 [ + - ]: 3604 : coordgp[1].resize( ng );
591 [ + - ]: 3604 : coordgp[2].resize( ng );
592 [ + - ]: 3604 : wgp.resize( ng );
593 : :
594 : : // get quadrature point weights and coordinates for triangle
595 [ + - ]: 3604 : tk::GaussQuadratureTet( ng, coordgp, wgp );
596 : :
597 : : // right hand side vector
598 [ + - ][ - - ]: 3604 : std::vector< tk::real > R( ncomp*ndof, 0.0 );
599 : :
600 : : // Gaussian quadrature
601 [ + + ]: 43248 : for (std::size_t igp=0; igp<ng; ++igp)
602 : : {
603 [ + - ]: 39644 : auto wt = wgp[igp] * vol;
604 : :
605 [ + - ]: 39644 : auto gp = tk::eval_gp( igp, coordel, coordgp );
606 : :
607 [ + - ]: 39644 : auto B_taylor = eval_TaylorBasis( ndof, gp, x_c, coordel);
608 : :
609 : : // Compute high order solution at gauss point
610 [ + - ][ - - ]: 39644 : std::vector< tk::real > state( ncomp, 0.0 );
611 [ + + ]: 237864 : for (ncomp_t c=0; c<ncomp; ++c)
612 : : {
613 [ + - ]: 198220 : state[c] = unk[c][0];
614 : 198220 : state[c] += unk[c][1] * B_taylor[1]
615 : 198220 : + unk[c][2] * B_taylor[2]
616 : 198220 : + unk[c][3] * B_taylor[3];
617 : :
618 [ + - ]: 198220 : if(ndof > 4)
619 : 198220 : state[c] += unk[c][4] * B_taylor[4] + unk[c][5] * B_taylor[5]
620 : 198220 : + unk[c][6] * B_taylor[6] + unk[c][7] * B_taylor[7]
621 : 198220 : + unk[c][8] * B_taylor[8] + unk[c][9] * B_taylor[9];
622 : : }
623 : :
624 [ + - ]: 39644 : auto B = tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
625 : :
626 [ + + ]: 237864 : for (ncomp_t c=0; c<ncomp; ++c)
627 : : {
628 : 198220 : auto mark = c*ndof;
629 [ + - ]: 198220 : R[mark] += wt * state[c];
630 : :
631 [ + - ]: 198220 : if(ndof > 1)
632 : : {
633 [ + - ]: 198220 : R[mark+1] += wt * state[c] * B[1];
634 : 198220 : R[mark+2] += wt * state[c] * B[2];
635 : 198220 : R[mark+3] += wt * state[c] * B[3];
636 : :
637 [ + - ]: 198220 : if(ndof > 4)
638 : : {
639 : 198220 : R[mark+4] += wt * state[c] * B[4];
640 : 198220 : R[mark+5] += wt * state[c] * B[5];
641 : 198220 : R[mark+6] += wt * state[c] * B[6];
642 : 198220 : R[mark+7] += wt * state[c] * B[7];
643 : 198220 : R[mark+8] += wt * state[c] * B[8];
644 : 198220 : R[mark+9] += wt * state[c] * B[9];
645 : : }
646 : : }
647 : : }
648 : : }
649 : :
650 [ + + ]: 21624 : for (ncomp_t c=0; c<ncomp; ++c)
651 : : {
652 : 18020 : auto mark = c*ndof;
653 [ + + ]: 198220 : for(std::size_t idof = 0; idof < ndof; idof++)
654 : 180200 : unk[c][idof] = R[mark+idof] / L[idof];
655 : : }
656 : 3604 : }
657 : :
658 : : std::vector< tk::real >
659 : 54060 : tk::eval_TaylorBasis( const std::size_t ndof,
660 : : const std::array< tk::real, 3 >& x,
661 : : const std::array< tk::real, 3 >& x_c,
662 : : const std::array< std::array< tk::real, 3>, 4 >& coordel )
663 : : // *****************************************************************************
664 : : // Evaluate the Taylor basis at points
665 : : //! \param[in] ndof Maximum number of degrees of freedom
666 : : //! \param[in] x Nodal coordinates
667 : : //! \param[in] x_c Coordinates of the centroid
668 : : //! \param[in] coordel Array of nodal coordinates for the tetrahedron
669 : : // *****************************************************************************
670 : : {
671 : 54060 : std::vector< tk::real > avg( 6, 0.0 );
672 [ + - ]: 54060 : if(ndof > 4)
673 : : {
674 : : Assert( ndof == 10, "Mismatch in number of degrees of freedom" );
675 [ + - ]: 54060 : auto ng = tk::NGvol(ndof);
676 : :
677 : : std::array< std::vector< tk::real >, 3 > coordgp;
678 : : std::vector< tk::real > wgp;
679 : :
680 [ + - ]: 54060 : coordgp[0].resize( ng );
681 [ + - ]: 54060 : coordgp[1].resize( ng );
682 [ + - ]: 54060 : coordgp[2].resize( ng );
683 [ + - ]: 54060 : wgp.resize( ng );
684 : :
685 [ + - ]: 54060 : tk::GaussQuadratureTet( ng, coordgp, wgp );
686 : :
687 [ + + ]: 648720 : for (std::size_t igp=0; igp<ng; ++igp)
688 : : {
689 : : // Compute the coordinates of quadrature point at physical domain
690 [ + - ]: 594660 : auto gp = tk::eval_gp( igp, coordel, coordgp );
691 : :
692 : 594660 : avg[0] += wgp[igp] * (gp[0] - x_c[0]) * (gp[0] - x_c[0]) * 0.5;
693 : 594660 : avg[1] += wgp[igp] * (gp[1] - x_c[1]) * (gp[1] - x_c[1]) * 0.5;
694 : 594660 : avg[2] += wgp[igp] * (gp[2] - x_c[2]) * (gp[2] - x_c[2]) * 0.5;
695 : 594660 : avg[3] += wgp[igp] * (gp[0] - x_c[0]) * (gp[1] - x_c[1]);
696 : 594660 : avg[4] += wgp[igp] * (gp[0] - x_c[0]) * (gp[2] - x_c[2]);
697 : 594660 : avg[5] += wgp[igp] * (gp[1] - x_c[1]) * (gp[2] - x_c[2]);
698 : : }
699 : : }
700 : :
701 [ + - ][ - - ]: 54060 : std::vector< tk::real > B( ndof, 1.0 );
702 : :
703 [ + - ]: 54060 : if(ndof > 1) {
704 : : Assert( (ndof == 4)||(ndof == 10) ,
705 : : "Mismatch in number of degrees of freedom" );
706 : 54060 : B[1] = x[0] - x_c[0];
707 : 54060 : B[2] = x[1] - x_c[1];
708 : 54060 : B[3] = x[2] - x_c[2];
709 : : }
710 : :
711 [ + - ]: 54060 : if(ndof > 4) {
712 : 54060 : B[4] = B[1] * B[1] * 0.5 - avg[0];
713 : 54060 : B[5] = B[2] * B[2] * 0.5 - avg[1];
714 : 54060 : B[6] = B[3] * B[3] * 0.5 - avg[2];
715 : 54060 : B[7] = B[1] * B[2] - avg[3];
716 : 54060 : B[8] = B[1] * B[3] - avg[4];
717 : 54060 : B[9] = B[2] * B[3] - avg[5];
718 : : }
719 : :
720 : 54060 : return B;
721 : : }
|
