Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/CGPDE.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 common to ALECG
9 : : \details Functions common to ALECG.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <array>
14 : : #include <vector>
15 : : #include <unordered_map>
16 : :
17 : : #include "Vector.hpp"
18 : : #include "DerivedData.hpp"
19 : : #include "Exception.hpp"
20 : : #include "Around.hpp"
21 : : #include "Fields.hpp"
22 : : #include "CGPDE.hpp"
23 : : #include "FunctionPrototypes.hpp"
24 : :
25 : : namespace inciter {
26 : : namespace cg {
27 : :
28 : : std::vector< tk::real >
29 [ - - ]: 0 : solinc( tk::ncomp_t ncomp,
30 : : const std::vector< EOS >& mat_blk, tk::real x, tk::real y,
31 : : tk::real z, tk::real t, tk::real dt, tk::InitializeFn solution )
32 : : // *****************************************************************************
33 : : // Evaluate the increment from t to t+dt of the analytical solution at (x,y,z)
34 : : // for all components
35 : : //! \param[in] ncomp Number of scalar components in this PDE system
36 : : //! \param[in] x X coordinate where to evaluate the solution
37 : : //! \param[in] y Y coordinate where to evaluate the solution
38 : : //! \param[in] z Z coordinate where to evaluate the solution
39 : : //! \param[in] t Time where to evaluate the solution increment starting from
40 : : //! \param[in] dt Time increment at which evaluate the solution increment to
41 : : //! \param[in] solution Function used to evaluate the solution
42 : : //! \return Increment in values of all components evaluated at (x,y,z,t+dt)
43 : : // *****************************************************************************
44 : : {
45 : 0 : auto st1 = solution( ncomp, mat_blk, x, y, z, t );
46 [ - - ]: 0 : auto st2 = solution( ncomp, mat_blk, x, y, z, t+dt );
47 : :
48 : : std::transform( begin(st1), end(st1), begin(st2), begin(st2),
49 : 0 : []( tk::real s, tk::real& d ){ return d -= s; } );
50 : :
51 : 0 : return st2;
52 : : }
53 : :
54 : : std::unordered_map< int,
55 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >
56 : 4760 : bnorm(
57 : : const std::map< int, std::vector< std::size_t > >& bface,
58 : : const std::vector< std::size_t >& triinpoel,
59 : : const std::array< std::vector< tk::real >, 3 >& coord,
60 : : const std::vector< std::size_t >& gid,
61 : : const std::unordered_map< int, std::unordered_set< std::size_t > >& bcnodes )
62 : : // *****************************************************************************
63 : : //! Compute boundary point normals
64 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
65 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (local ids)
66 : : //! \param[in] coord Mesh node coordinates
67 : : //! \param[in] gid Local->global node id map
68 : : //! \param[in] bcnodes Local node ids associated to side set ids at which BCs
69 : : //! are set that require normals
70 : : //! \return Face normals in boundary points, Inner key: local node id, value:
71 : : //! unit normal and inverse distance square between face centroids and points,
72 : : //! outer key: side set id
73 : : // *****************************************************************************
74 : : {
75 : : const auto& x = coord[0];
76 : : const auto& y = coord[1];
77 : : const auto& z = coord[2];
78 : :
79 : : // Lambda to compute the inverse distance squared between boundary face
80 : : // centroid and boundary point. Here p is the global node id and g is the
81 : : // geometry of the boundary face, see tk::geoFaceTri().
82 : 4039074 : auto invdistsq = [&]( const tk::Fields& g, std::size_t p ){
83 : 4039074 : return 1.0 / ( (g(0,4) - x[p])*(g(0,4) - x[p]) +
84 : 4039074 : (g(0,5) - y[p])*(g(0,5) - y[p]) +
85 : 4039074 : (g(0,6) - z[p])*(g(0,6) - z[p]) );
86 : 2380 : };
87 : :
88 : : // Compute boundary point normals on all side sets summing inverse distance
89 : : // weighted face normals to points. This is only a partial sum at shared
90 : : // boundary points in parallel. Inner key: global node id, value: normals and
91 : : // inverse distance square, outer key, side set id.
92 : : std::unordered_map< int,
93 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > norm;
94 [ + + ]: 7117 : for (const auto& [ setid, faceids ] : bface) { // for all side sets
95 [ + + ]: 1420857 : for (auto f : faceids) { // for all side set triangles
96 : : tk::UnsMesh::Face
97 [ + - ]: 1416120 : face{ triinpoel[f*3+0], triinpoel[f*3+1], triinpoel[f*3+2] };
98 : 1416120 : std::array< tk::real, 3 > fx{ x[face[0]], x[face[1]], x[face[2]] };
99 : 1416120 : std::array< tk::real, 3 > fy{ y[face[0]], y[face[1]], y[face[2]] };
100 : 1416120 : std::array< tk::real, 3 > fz{ z[face[0]], z[face[1]], z[face[2]] };
101 [ + - ]: 1416120 : auto g = tk::geoFaceTri( fx, fy, fz );
102 [ + + ]: 5664480 : for (auto p : face) { // for all 3 nodes of a boundary triangle face
103 [ + + ]: 13991574 : for (const auto& [s,nodes] : bcnodes) { // for all bnd nodes w/ normals
104 [ + + ]: 9743214 : if (setid == s) { // only contribute to side set we operate on
105 : : auto i = nodes.find(p);
106 [ + - ]: 4039074 : if (i != end(nodes)) { // only if user set bc on node
107 [ + - ]: 4039074 : tk::real r = invdistsq(g,p);
108 [ + - ]: 4039074 : auto& n = norm[s][gid[p]]; // associate set id and global node id
109 : 4039074 : n[0] += r*g(0,1);
110 : 4039074 : n[1] += r*g(0,2);
111 : 4039074 : n[2] += r*g(0,3);
112 : 4039074 : n[3] += r;
113 : : }
114 : : }
115 : : }
116 : : }
117 : : }
118 : : }
119 : :
120 : 2380 : return norm;
121 : : }
122 : :
123 : : } // cg::
124 : : } // inciter::
|
