Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/SolidTerms.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 integrals for the RHS terms
9 : : of the inverse deformation equations in the multi-material
10 : : solid dynamics solver
11 : : \details This file contains routines to integrate right hand side terms
12 : : to the inverse deformation equations for the multi-material
13 : : solid dynamics solver.
14 : : */
15 : : // *****************************************************************************
16 : :
17 : : #include <array>
18 : : #include <vector>
19 : :
20 : : #include "SolidTerms.hpp"
21 : : #include "Vector.hpp"
22 : : #include "Quadrature.hpp"
23 : : #include "Reconstruction.hpp"
24 : : #include "MultiMatTerms.hpp"
25 : : #include "MultiMat/MultiMatIndexing.hpp"
26 : : #include "MultiMat/MiscMultiMatFns.hpp"
27 : : #include "Inciter/InputDeck/InputDeck.hpp"
28 : :
29 : : namespace inciter {
30 : : extern ctr::InputDeck g_inputdeck;
31 : : }
32 : :
33 : : namespace tk {
34 : :
35 : : void
36 : 0 : solidTermsVolInt(
37 : : std::size_t nmat,
38 : : const std::vector< inciter::EOS >& mat_blk,
39 : : const std::size_t ndof,
40 : : const std::size_t rdof,
41 : : const std::size_t nelem,
42 : : const std::vector< std::size_t >& inpoel,
43 : : const UnsMesh::Coords& coord,
44 : : const Fields& geoElem,
45 : : const Fields& U,
46 : : const Fields& P,
47 : : const std::vector< std::size_t >& ndofel,
48 : : const tk::real dt,
49 : : Fields& R,
50 : : int intsharp )
51 : : // *****************************************************************************
52 : : // Compute all RHS volume terms in the inverse deformation equations to
53 : : // satisfy the condition curl(g) = 0.
54 : : //! \param[in] nmat Number of materials in this PDE system
55 : : //! \param[in] mat_blk EOS material block
56 : : //! \param[in] ndof Maximum number of degrees of freedom
57 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
58 : : //! \param[in] nelem Maximum number of elements
59 : : //! \param[in] inpoel Element-node connectivity
60 : : //! \param[in] coord Array of nodal coordinates
61 : : //! \param[in] geoElem Element geometry array
62 : : //! \param[in] U Solution vector at recent time step
63 : : //! \param[in] P Vector of primitives at recent time step
64 : : //! \param[in] ndofel Vector of local number of degrees of freedom
65 : : //! \param[in] dt Delta time
66 : : //! \param[in,out] R Right-hand side vector computed
67 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
68 : : //! default 0, so that it is unused for single-material and transport.
69 : : // *****************************************************************************
70 : : {
71 : :
72 : : using inciter::velocityDofIdx;
73 : : using inciter::volfracDofIdx;
74 : : using inciter::deformDofIdx;
75 : :
76 : : const auto& solidx =
77 : : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
78 : :
79 : 0 : auto ncomp = R.nprop()/ndof;
80 : 0 : auto nprim = P.nprop()/rdof;
81 : : const auto& cx = coord[0];
82 : : const auto& cy = coord[1];
83 : : const auto& cz = coord[2];
84 : :
85 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
86 : : {
87 : :
88 : : // Relaxation coefficient
89 : 0 : tk::real eta = 1.0/(6.0*dt);
90 : :
91 : 0 : auto ng = tk::NGvol(ndofel[e]);
92 : :
93 : : // arrays for quadrature points
94 : : std::array< std::vector< real >, 3 > coordgp;
95 : : std::vector< real > wgp;
96 : :
97 [ - - ]: 0 : coordgp[0].resize( ng );
98 [ - - ]: 0 : coordgp[1].resize( ng );
99 [ - - ]: 0 : coordgp[2].resize( ng );
100 [ - - ]: 0 : wgp.resize( ng );
101 : :
102 [ - - ]: 0 : GaussQuadratureTet( ng, coordgp, wgp );
103 : :
104 : : // Extract the element coordinates
105 : : std::array< std::array< real, 3>, 4 > coordel {{
106 : 0 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
107 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
108 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
109 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}}};
110 : :
111 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
112 : : {
113 : : // Compute the coordinates of quadrature point at physical domain
114 [ - - ]: 0 : auto gp = eval_gp( igp, coordel, coordgp );
115 : :
116 : : // Compute the basis function
117 : : auto B =
118 [ - - ]: 0 : eval_basis(rdof,coordgp[0][igp],coordgp[1][igp],coordgp[2][igp]);
119 : :
120 : : // Get state
121 : : std::vector< real > state;
122 [ - - ]: 0 : state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
123 : : nmat, e, rdof, inpoel, coord, geoElem, gp, B, U, P);
124 : :
125 : : // Loop through materials
126 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
127 : : {
128 [ - - ]: 0 : if (solidx[k] > 0)
129 : : {
130 : 0 : tk::real alpha = state[inciter::volfracIdx(nmat, k)];
131 : : std::array< std::array< tk::real, 3 >, 3 > g;
132 : : // Compute the source terms
133 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
134 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
135 : 0 : g[i][j] = state[inciter::deformIdx(nmat,solidx[k],i,j)];
136 : :
137 : : // Compute rhs factor
138 : 0 : tk::real rho = state[inciter::densityIdx(nmat, k)]/alpha;
139 [ - - ]: 0 : tk::real rho0 = mat_blk[k].compute< inciter::EOS::rho0 >();
140 : 0 : tk::real rfact = eta*(rho/(rho0*tk::determinant(g))-1.0);
141 : :
142 : : // Compute the source terms
143 [ - - ][ - - ]: 0 : std::vector< real > s(9*ndof, 0.0);
144 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
145 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
146 [ - - ]: 0 : for (std::size_t idof=0; idof<ndof; ++idof)
147 : 0 : s[(i*3+j)*ndof+idof] = B[idof]*rfact*g[i][j];
148 : :
149 [ - - ]: 0 : auto wt = wgp[igp] * geoElem(e, 0);
150 : :
151 [ - - ]: 0 : update_rhs( nmat, ndof, wt, e, solidx[k], s, R );
152 : :
153 : : }
154 : : }
155 : : }
156 : : }
157 : :
158 : 0 : }
159 : :
160 : : void
161 : 0 : update_rhs( std::size_t nmat,
162 : : const std::size_t ndof,
163 : : const tk::real wt,
164 : : const std::size_t e,
165 : : const std::size_t solidx,
166 : : const std::vector< real >& s,
167 : : Fields& R )
168 : : // *****************************************************************************
169 : : // Update the rhs by adding the source term integrals
170 : : //! \param[in] nmat Number of materials in this PDE system
171 : : //! \param[in] rdof Maximum number of degrees of freedom
172 : : //! \param[in] wt Weight of gauss quadrature point
173 : : //! \param[in] e Element index
174 : : //! \param[in] solidx Material index indicator
175 : : //! \param[in] s Source term vector
176 : : //! \param[in,out] R Right-hand side vector computed
177 : : // *****************************************************************************
178 : : {
179 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
180 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
181 [ - - ]: 0 : for (std::size_t idof=0; idof<ndof; ++idof)
182 : : {
183 : : std::size_t dofId = inciter::deformDofIdx(nmat,solidx,i,j,ndof,idof);
184 : 0 : std::size_t srcId = (i*3+j)*ndof+idof;
185 : 0 : R(e, dofId) += wt * s[srcId];
186 : : }
187 : 0 : }
188 : :
189 : : } // tk::
|
