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