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 : 1104211 : solinc( tk::ncomp_t system, tk::ncomp_t ncomp, tk::real x, tk::real y,
30 : : tk::real z, tk::real t, tk::real dt, tk::InitializeFn solution )
31 : : // *****************************************************************************
32 : : // Evaluate the increment from t to t+dt of the analytical solution at (x,y,z)
33 : : // for all components
34 : : //! \param[in] system Equation system index, i.e., which equation system we
35 : : //! operate on among the systems of PDEs configured by the user
36 : : //! \param[in] ncomp Number of scalar components in this PDE system
37 : : //! \param[in] x X coordinate where to evaluate the solution
38 : : //! \param[in] y Y coordinate where to evaluate the solution
39 : : //! \param[in] z Z coordinate where to evaluate the solution
40 : : //! \param[in] t Time where to evaluate the solution increment starting from
41 : : //! \param[in] dt Time increment at which evaluate the solution increment to
42 : : //! \param[in] solution Function used to evaluate the solution
43 : : //! \return Increment in values of all components evaluated at (x,y,z,t+dt)
44 : : // *****************************************************************************
45 : : {
46 [ + - ]: 2208422 : auto st1 = solution( system, ncomp, x, y, z, t );
47 [ + - ]: 1104211 : auto st2 = solution( system, ncomp, x, y, z, t+dt );
48 : :
49 : : std::transform( begin(st1), end(st1), begin(st2), begin(st2),
50 : 5795926 : []( tk::real s, tk::real& d ){ return d -= s; } );
51 : :
52 : 2208422 : return st2;
53 : : }
54 : :
55 : : std::unordered_map< int,
56 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >
57 : 2268 : bnorm(
58 : : const std::map< int, std::vector< std::size_t > >& bface,
59 : : const std::vector< std::size_t >& triinpoel,
60 : : const std::array< std::vector< tk::real >, 3 >& coord,
61 : : const std::vector< std::size_t >& gid,
62 : : const std::unordered_map< int, std::unordered_set< std::size_t > >& bcnodes )
63 : : // *****************************************************************************
64 : : //! Compute boundary point normals
65 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
66 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (local ids)
67 : : //! \param[in] coord Mesh node coordinates
68 : : //! \param[in] gid Local->global node id map
69 : : //! \param[in] bcnodes Local node ids associated to side set ids at which BCs
70 : : //! are set that require normals
71 : : //! \return Face normals in boundary points, Inner key: local node id, value:
72 : : //! unit normal and inverse distance square between face centroids and points,
73 : : //! outer key: side set id
74 : : // *****************************************************************************
75 : : {
76 : 2268 : const auto& x = coord[0];
77 : 2268 : const auto& y = coord[1];
78 : 2268 : const auto& z = coord[2];
79 : :
80 : : // Lambda to compute the inverse distance squared between boundary face
81 : : // centroid and boundary point. Here p is the global node id and g is the
82 : : // geometry of the boundary face, see tk::geoFaceTri().
83 : 3633768 : auto invdistsq = [&]( const tk::Fields& g, std::size_t p ){
84 : 3633768 : return 1.0 / ( (g(0,4,0) - x[p])*(g(0,4,0) - x[p]) +
85 : 3633768 : (g(0,5,0) - y[p])*(g(0,5,0) - y[p]) +
86 : 3633768 : (g(0,6,0) - z[p])*(g(0,6,0) - z[p]) );
87 : 2268 : };
88 : :
89 : : // Compute boundary point normals on all side sets summing inverse distance
90 : : // weighted face normals to points. This is only a partial sum at shared
91 : : // boundary points in parallel. Inner key: global node id, value: normals and
92 : : // inverse distance square, outer key, side set id.
93 : : std::unordered_map< int,
94 : 2268 : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > norm;
95 [ + + ]: 7373 : for (const auto& [ setid, faceids ] : bface) { // for all side sets
96 [ + + ]: 1337983 : for (auto f : faceids) { // for all side set triangles
97 : : tk::UnsMesh::Face
98 : 1332878 : face{ triinpoel[f*3+0], triinpoel[f*3+1], triinpoel[f*3+2] };
99 : 1332878 : std::array< tk::real, 3 > fx{ x[face[0]], x[face[1]], x[face[2]] };
100 : 1332878 : std::array< tk::real, 3 > fy{ y[face[0]], y[face[1]], y[face[2]] };
101 : 1332878 : std::array< tk::real, 3 > fz{ z[face[0]], z[face[1]], z[face[2]] };
102 [ + - ]: 2665756 : auto g = tk::geoFaceTri( fx, fy, fz );
103 [ + + ]: 5331512 : for (auto p : face) { // for all 3 nodes of a boundary triangle face
104 [ + + ]: 13010514 : for (const auto& [s,nodes] : bcnodes) { // for all bnd nodes w/ normals
105 [ + + ]: 9011880 : if (setid == s) { // only contribute to side set we operate on
106 [ + - ]: 3633768 : auto i = nodes.find(p);
107 [ + - ]: 3633768 : if (i != end(nodes)) { // only if user set bc on node
108 [ + - ]: 3633768 : tk::real r = invdistsq(g,p);
109 [ + - ][ + - ]: 3633768 : auto& n = norm[s][gid[p]]; // associate set id and global node id
110 [ + - ]: 3633768 : n[0] += r*g(0,1,0);
111 [ + - ]: 3633768 : n[1] += r*g(0,2,0);
112 [ + - ]: 3633768 : n[2] += r*g(0,3,0);
113 : 3633768 : n[3] += r;
114 : : }
115 : : }
116 : : }
117 : : }
118 : : }
119 : : }
120 : :
121 : 4536 : return norm;
122 : : }
123 : :
124 : : } // cg::
125 : : } // inciter::
|
