Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Mesh/Gradients.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 computing gradients on unstructured meshes for tetrahedra
9 : : \details Functions computing gradients using linear finite element shape
10 : : functions on unstructured meshes for tetrahedra.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <cstddef>
15 : :
16 : : #include "Exception.hpp"
17 : : #include "Gradients.hpp"
18 : : #include "Vector.hpp"
19 : : #include "Around.hpp"
20 : :
21 : : namespace tk {
22 : :
23 : : std::array< tk::real, 3 >
24 : 882 : nodegrad( std::size_t node,
25 : : const std::array< std::vector< tk::real >, 3 >& coord,
26 : : const std::vector< std::size_t >& inpoel,
27 : : const std::pair< std::vector< std::size_t >,
28 : : std::vector< std::size_t > >& esup,
29 : : const tk::Fields& U,
30 : : ncomp_t c )
31 : : // *****************************************************************************
32 : : // Compute gradient at a mesh node
33 : : //! \param[in] node Node id at which to compute gradient
34 : : //! \param[in] coord Mesh node coordinates
35 : : //! \param[in] inpoel Mesh element connectivity
36 : : //! \param[in] esup Linked lists storing elements surrounding points, see
37 : : //! tk::genEsup()
38 : : //! \param[in] U Field vector whose component gradient to compute
39 : : //! \param[in] c Scalar component to compute gradient of
40 : : //! \return Gradient of U(c) at mesh node
41 : : // *****************************************************************************
42 : : {
43 [ - + ][ - - ]: 882 : Assert( c < U.nprop(), "Indexing out of field data" );
[ - - ][ - - ]
44 : :
45 : 882 : const auto& x = coord[0];
46 : 882 : const auto& y = coord[1];
47 : 882 : const auto& z = coord[2];
48 : :
49 : : // storage for gradient and volume at the mesh node
50 : 882 : std::array< tk::real, 3 > g{{ 0.0, 0.0, 0.0 }};
51 : 882 : tk::real vol = 0.0;
52 : :
53 : : // loop over cells surrounding mesh node
54 [ + + ]: 7378 : for (auto e : tk::Around(esup,node)) {
55 : : // access node IDs
56 : 6496 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
57 : 6496 : inpoel[e*4+2], inpoel[e*4+3] }};
58 : :
59 : : // compute element Jacobi determinant
60 : : const std::array< tk::real, 3 >
61 : 6496 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
62 : 6496 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
63 : 6496 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
64 : 6496 : const auto J = tk::triple( ba, ca, da ); // J = 6V
65 [ - + ][ - - ]: 6496 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
66 : :
67 : : // shape function derivatives, nnode*ndim [4][3]
68 : : std::array< std::array< tk::real, 3 >, 4 > grad;
69 : 6496 : grad[1] = tk::crossdiv( ca, da, J );
70 : 6496 : grad[2] = tk::crossdiv( da, ba, J );
71 : 6496 : grad[3] = tk::crossdiv( ba, ca, J );
72 [ + + ]: 25984 : for (std::size_t i=0; i<3; ++i)
73 : 19488 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
74 : :
75 : : // access field data for scalar component c at nodes of element
76 [ + - ]: 6496 : auto u = U.extract( c, 0, N );
77 : :
78 : : // compute nodal volume: every element contributes their volume / 4
79 : 6496 : vol += 5.0*J/120.0;
80 : :
81 : : // compute gradient over element weighed by cell volume / 4 and sum to node
82 [ + + ]: 25984 : for (std::size_t j=0; j<3; ++j) {
83 [ + + ]: 97440 : for (std::size_t i=0; i<4; ++i) {
84 : 77952 : g[j] += grad[i][j] * u[i] * 5.0*J/120.0;
85 : : }
86 : : }
87 : : }
88 : :
89 : : // divide components of nodal gradient by nodal volume
90 [ + + ]: 3528 : for (std::size_t j=0; j<3; ++j) g[j] /= vol;
91 : :
92 : 882 : return g;
93 : : }
94 : :
95 : : std::array< tk::real, 3 >
96 : 343 : edgegrad( const std::array< std::vector< tk::real >, 3 >& coord,
97 : : const std::vector< std::size_t >& inpoel,
98 : : const std::vector< std::size_t >& esued,
99 : : const tk::Fields& U,
100 : : ncomp_t c )
101 : : // *****************************************************************************
102 : : // Compute gradient at a mesh edge
103 : : //! \param[in] coord Mesh node coordinates
104 : : //! \param[in] inpoel Mesh element connectivity
105 : : //! \param[in] esued List of elements surrounding edge, see tk::genEsued()
106 : : //! \param[in] U Field vector whose component gradient to compute
107 : : //! \param[in] c Scalar component to compute gradient of
108 : : //! \return Gradient of U(c) at mesh edge
109 : : // *****************************************************************************
110 : : {
111 [ - + ][ - - ]: 343 : Assert( c < U.nprop(), "Indexing out of field data" );
[ - - ][ - - ]
112 : :
113 : 343 : const auto& x = coord[0];
114 : 343 : const auto& y = coord[1];
115 : 343 : const auto& z = coord[2];
116 : :
117 : : // storage for gradient and volume at the mesh edge
118 : 343 : std::array< tk::real, 3 > g{{ 0.0, 0.0, 0.0 }};
119 : 343 : tk::real vol = 0.0;
120 : :
121 : : // loop over elements surrounding edge
122 [ + + ]: 1351 : for (auto e : esued) {
123 : : // access node IDs
124 : 1008 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
125 : 1008 : inpoel[e*4+2], inpoel[e*4+3] }};
126 : :
127 : : // compute element Jacobi determinant
128 : : const std::array< tk::real, 3 >
129 : 1008 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
130 : 1008 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
131 : 1008 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
132 : 1008 : const auto J = tk::triple( ba, ca, da ); // J = 6V
133 [ - + ][ - - ]: 1008 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
134 : :
135 : : // shape function derivatives, nnode*ndim [4][3]
136 : : std::array< std::array< tk::real, 3 >, 4 > grad;
137 : 1008 : grad[1] = tk::crossdiv( ca, da, J );
138 : 1008 : grad[2] = tk::crossdiv( da, ba, J );
139 : 1008 : grad[3] = tk::crossdiv( ba, ca, J );
140 [ + + ]: 4032 : for (std::size_t i=0; i<3; ++i)
141 : 3024 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
142 : :
143 : : // access field data for scalar component c at nodes of element
144 [ + - ]: 1008 : auto u = U.extract( c, 0, N );
145 : :
146 : : // compute edge volume: every element contributes their volume / 6
147 : 1008 : vol += J/36.0;
148 : :
149 : : // compute gradient over element weighed by cell volume / 6 and sum to edge
150 [ + + ]: 4032 : for (std::size_t j=0; j<3; ++j) {
151 [ + + ]: 15120 : for (std::size_t i=0; i<4; ++i) {
152 : 12096 : g[j] += grad[i][j] * u[i] * J/36.0;
153 : : }
154 : : }
155 : : }
156 : :
157 : : // divide components of gradient by edge volume
158 [ + + ]: 1372 : for (std::size_t j=0; j<3; ++j) g[j] /= vol;
159 : :
160 : 343 : return g;
161 : : }
162 : :
163 : : } // tk::
|
