Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/DGPDE.hpp
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 Partial differential equation base for discontinuous Galerkin PDEs
9 : : \details This file defines a generic partial differential equation (PDE)
10 : : class for PDEs that use discontinuous Galerkin spatial discretization.
11 : : The class uses runtime polymorphism without client-side inheritance:
12 : : inheritance is confined to the internals of the class, invisible to
13 : : client-code. The class exclusively deals with ownership enabling client-side
14 : : value semantics. Credit goes to Sean Parent at Adobe:
15 : : https://github.com/sean-parent/sean-parent.github.com/wiki/
16 : : Papers-and-Presentations.
17 : : */
18 : : // *****************************************************************************
19 : : #ifndef DGPDE_h
20 : : #define DGPDE_h
21 : :
22 : : #include <array>
23 : : #include <string>
24 : : #include <vector>
25 : : #include <memory>
26 : : #include <unordered_set>
27 : : #include <functional>
28 : :
29 : : #include "Types.hpp"
30 : : #include "Fields.hpp"
31 : : #include "FaceData.hpp"
32 : : #include "UnsMesh.hpp"
33 : : #include "Inciter/InputDeck/InputDeck.hpp"
34 : : #include "FunctionPrototypes.hpp"
35 : : #include "History.hpp"
36 : :
37 : : namespace inciter {
38 : :
39 : : extern ctr::InputDeck g_inputdeck;
40 : :
41 : : using ncomp_t = tk::ncomp_t;
42 : : using BCStateFn =
43 : : std::vector< std::tuple< std::vector< std::size_t >, tk::StateFn, tk::StateFn > >;
44 : :
45 : : //! Extract BC configuration ignoring if BC not specified
46 : : //! \note A more preferable way of catching errors such as this function
47 : : //! hides is during parsing, so that we don't even get here if BCs are
48 : : //! not correctly specified. For now we simply ignore if BCs are not
49 : : //! specified by allowing empty BC vectors from the user input.
50 : : struct ConfigBC {
51 : : BCStateFn& state; //!< BC state config: sidesets + statefn
52 : : const std::vector< tk::StateFn >& fn; //!< BC state functions
53 : : const std::vector< tk::StateFn >& gfn; //!< BC gradient functions
54 : : std::size_t c; //!< Counts BC types configured
55 : : //! Constructor
56 : 462 : ConfigBC( BCStateFn& s,
57 : : const std::vector< tk::StateFn >& f,
58 : 462 : const std::vector< tk::StateFn >& gf ) :
59 : 462 : state(s), fn(f), gfn(gf), c(0) {}
60 : : //! Function to call for each BC type
61 : 3234 : template< typename U > void operator()( brigand::type_<U> ) {
62 : 6468 : std::vector< std::size_t > cfg, v;
63 : : // collect sidesets across all meshes
64 [ + + ]: 6468 : for (const auto& ibc : g_inputdeck.get< tag::bc >()) {
65 [ + - ]: 3234 : v.insert(v.end(), ibc.get< U >().begin(), ibc.get< U >().end());
66 : : }
67 [ + + ][ + - ]: 3234 : if (v.size() > 0) cfg = v;
68 [ - + ][ - - ]: 3234 : Assert( fn.size() > c, "StateFn missing for BC type" );
[ - - ][ - - ]
69 [ + - ][ + - ]: 3234 : state.push_back( { cfg, fn[c], gfn[c] } );
70 : 3234 : ++c;
71 : 3234 : }
72 : : };
73 : :
74 : : //! \brief Extract information on inlet BCs, which have a different structure
75 : : //! than other BCs
76 : : void ConfigInletBC( BCStateFn&,
77 : : const tk::StateFn&,
78 : : const tk::StateFn& );
79 : :
80 : : //! \brief Extract information on back pressure BCs, which have a different
81 : : //! structure than other BCs
82 : : void ConfigBackPressureBC( BCStateFn&,
83 : : const tk::StateFn&,
84 : : const tk::StateFn& );
85 : :
86 : : //! State function for invalid/un-configured boundary conditions
87 : : [[noreturn]] tk::StateFn::result_type
88 : : invalidBC( ncomp_t, const std::vector< EOS >&,
89 : : const std::vector< tk::real >&, tk::real, tk::real, tk::real,
90 : : tk::real, const std::array< tk::real, 3> & );
91 : :
92 : : //! \brief Partial differential equation base for discontinuous Galerkin PDEs
93 : : //! \details This class uses runtime polymorphism without client-side
94 : : //! inheritance: inheritance is confined to the internals of the this class,
95 : : //! invisible to client-code. The class exclusively deals with ownership
96 : : //! enabling client-side value semantics. Credit goes to Sean Parent at Adobe:
97 : : //! https://github.com/sean-parent/sean-parent.github.com/wiki/
98 : : //! Papers-and-Presentations. For example client code that models a DGPDE,
99 : : //! see inciter::CompFlow.
100 : : class DGPDE {
101 : :
102 : : private:
103 : : using ncomp_t = tk::ncomp_t;
104 : :
105 : : public:
106 : : //! Default constructor taking no arguments for Charm++
107 : 0 : explicit DGPDE() = default;
108 : :
109 : : //! \brief Constructor taking an object modeling Concept.
110 : : //! \details The object of class T comes pre-constructed.
111 : : //! \param[in] x Instantiated object of type T given by the template
112 : : //! argument.
113 : : template< typename T > explicit DGPDE( T x ) :
114 : : self( std::make_unique< Model<T> >( std::move(x) ) ) {}
115 : :
116 : : //! \brief Constructor taking a function pointer to a constructor of an
117 : : //! object modeling Concept.
118 : : //! \details Passing std::function allows late execution of the constructor,
119 : : //! i.e., as late as inside this class' constructor, and thus usage from
120 : : //! a factory. Note that there are at least two different ways of using
121 : : //! this constructor:
122 : : //! - Bind T's constructor arguments and place it in std::function<T()>
123 : : //! and passing no arguments as args.... This case then instantiates the
124 : : //! model via its constructor and stores it in here.
125 : : //! - Bind a single placeholder argument to T's constructor and pass it in
126 : : //! as host's args..., which then forwards it to model's constructor. This
127 : : //! allows late binding, i.e., binding the argument only here.
128 : : //! \see See also the wrapper tk::recordModel() which does the former and
129 : : //! tk::recordModelLate() which does the latter, both defined in
130 : : //! src/Base/Factory.h.
131 : : //! \param[in] x Function pointer to a constructor of an object modeling
132 : : //! Concept.
133 : : //! \param[in] args Zero or more constructor arguments
134 : : template< typename T, typename...Args >
135 : 367 : explicit DGPDE( std::function<T(Args...)> x, Args&&... args ) :
136 : : self( std::make_unique< Model<T> >(
137 [ + - ]: 367 : std::move( x( std::forward<Args>(args)... ) ) ) ) {}
138 : :
139 : : //! Public interface to find number of primitive quantities for the diff eq
140 : 671 : std::size_t nprim() const
141 : 671 : { return self->nprim(); }
142 : :
143 : : //! Public interface to find number of materials for the diff eq
144 : : std::size_t nmat() const
145 : : { return self->nmat(); }
146 : :
147 : : //! Public interface to find Dofs for each equation in pde system
148 : 671 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
149 : 671 : { return self->numEquationDofs(numEqDof); }
150 : :
151 : : //! Public interface to find how 'stiff equations', which are the inverse
152 : : //! deformation equations because of plasticity
153 : 2684 : std::size_t nstiffeq() const
154 : 2684 : { return self->nstiffeq(); }
155 : :
156 : : //! Public interface to find how 'nonstiff equations', which are the inverse
157 : : //! deformation equations because of plasticity
158 : 1342 : std::size_t nnonstiffeq() const
159 : 1342 : { return self->nnonstiffeq(); }
160 : :
161 : : //! Public function to locate the stiff equations
162 : 0 : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
163 : 0 : { return self->setStiffEqIdx( stiffEqIdx ); }
164 : :
165 : : //! Public function to locate the nonstiff equations
166 : 0 : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
167 : 0 : { return self->setNonStiffEqIdx( nonStiffEqIdx ); }
168 : :
169 : : //! Public interface to determine elements that lie inside the IC box
170 : 671 : void IcBoxElems( const tk::Fields& geoElem,
171 : : std::size_t nielem,
172 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
173 : 671 : { self->IcBoxElems( geoElem, nielem, inbox ); }
174 : :
175 : : //! Public interface to setting the initial conditions for the diff eq
176 : 803 : void initialize(
177 : : const tk::Fields& L,
178 : : const std::vector< std::size_t >& inpoel,
179 : : const tk::UnsMesh::Coords& coord,
180 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
181 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
182 : : elemblkid,
183 : : tk::Fields& unk,
184 : : tk::real t,
185 : : const std::size_t nielem ) const
186 : 803 : { self->initialize( L, inpoel, coord, inbox, elemblkid, unk, t, nielem ); }
187 : :
188 : : //! Public interface for computing density constraint
189 : 1325 : void computeDensityConstr( std::size_t nelem,
190 : : tk::Fields& unk,
191 : : std::vector< tk::real >& densityConstr) const
192 : 1325 : { self->computeDensityConstr( nelem, unk, densityConstr); }
193 : :
194 : : //! Public interface to computing the left-hand side matrix for the diff eq
195 : 803 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const
196 : 803 : { self->lhs( geoElem, l ); }
197 : :
198 : : //! Public interface to updating the interface cells for the diff eq
199 : 38265 : void updateInterfaceCells( tk::Fields& unk,
200 : : std::size_t nielem,
201 : : std::vector< std::size_t >& ndofel,
202 : : std::vector< std::size_t >& interface ) const
203 : 38265 : { self->updateInterfaceCells( unk, nielem, ndofel, interface ); }
204 : :
205 : : //! Public interface to updating the primitives for the diff eq
206 : 38936 : void updatePrimitives( const tk::Fields& unk,
207 : : const tk::Fields& L,
208 : : const tk::Fields& geoElem,
209 : : tk::Fields& prim,
210 : : std::size_t nielem,
211 : : const std::vector< std::size_t >& ndofel ) const
212 : 38936 : { self->updatePrimitives( unk, L, geoElem, prim, nielem, ndofel ); }
213 : :
214 : : //! Public interface to cleaning up trace materials for the diff eq
215 : 38265 : void cleanTraceMaterial( tk::real t,
216 : : const tk::Fields& geoElem,
217 : : tk::Fields& unk,
218 : : tk::Fields& prim,
219 : : std::size_t nielem ) const
220 : 38265 : { self->cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
221 : :
222 : : //! Public interface to reconstructing the second-order solution
223 : 24840 : void reconstruct( tk::real t,
224 : : const tk::Fields& geoFace,
225 : : const tk::Fields& geoElem,
226 : : const inciter::FaceData& fd,
227 : : const std::map< std::size_t, std::vector< std::size_t > >&
228 : : esup,
229 : : const std::vector< std::size_t >& inpoel,
230 : : const tk::UnsMesh::Coords& coord,
231 : : tk::Fields& U,
232 : : tk::Fields& P,
233 : : const bool pref,
234 : : const std::vector< std::size_t >& ndofel ) const
235 : : {
236 : 24840 : self->reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P,
237 : 24840 : pref, ndofel );
238 : 24840 : }
239 : :
240 : : //! Public interface to limiting the second-order solution
241 : 24840 : void limit( tk::real t,
242 : : const bool pref,
243 : : const tk::Fields& geoFace,
244 : : const tk::Fields& geoElem,
245 : : const inciter::FaceData& fd,
246 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
247 : : const std::vector< std::size_t >& inpoel,
248 : : const tk::UnsMesh::Coords& coord,
249 : : const std::vector< std::size_t >& ndofel,
250 : : const std::vector< std::size_t >& gid,
251 : : const std::unordered_map< std::size_t, std::size_t >& bid,
252 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
253 : : const std::vector< std::vector<tk::real> >& pNodalExtrm,
254 : : const std::vector< std::vector<tk::real> >& mtInv,
255 : : tk::Fields& U,
256 : : tk::Fields& P,
257 : : std::vector< std::size_t >& shockmarker ) const
258 : : {
259 : 24840 : self->limit( t, pref, geoFace, geoElem, fd, esup, inpoel, coord, ndofel,
260 : 24840 : gid, bid, uNodalExtrm, pNodalExtrm, mtInv, U, P, shockmarker );
261 : 24840 : }
262 : :
263 : : //! Public interface to update the conservative variable solution
264 : 20790 : void CPL( const tk::Fields& prim,
265 : : const tk::Fields& geoElem,
266 : : const std::vector< std::size_t >& inpoel,
267 : : const tk::UnsMesh::Coords& coord,
268 : : tk::Fields& unk,
269 : : std::size_t nielem ) const
270 : : {
271 : 20790 : self->CPL( prim, geoElem, inpoel, coord, unk, nielem );
272 : 20790 : }
273 : :
274 : : //! Public interface to reset the high order solution for p-adaptive scheme
275 : 1770 : void resetAdapSol( const inciter::FaceData& fd,
276 : : tk::Fields& unk,
277 : : tk::Fields& prim,
278 : : const std::vector< std::size_t >& ndofel ) const
279 : : {
280 : 1770 : self->resetAdapSol( fd, unk, prim, ndofel );
281 : 1770 : }
282 : :
283 : : //! Public interface to getting the cell-averaged deformation gradients
284 : : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
285 : : const tk::Fields& U,
286 : : std::size_t nielem ) const
287 : : {
288 : : return self->cellAvgDeformGrad( U, nielem );
289 : : }
290 : :
291 : : //! Public interface to computing the P1 right-hand side vector
292 : 38265 : void rhs( tk::real t,
293 : : const bool pref,
294 : : const tk::Fields& geoFace,
295 : : const tk::Fields& geoElem,
296 : : const inciter::FaceData& fd,
297 : : const std::vector< std::size_t >& inpoel,
298 : : const std::vector< std::unordered_set< std::size_t > >& boxelems,
299 : : const tk::UnsMesh::Coords& coord,
300 : : const tk::Fields& U,
301 : : const tk::Fields& P,
302 : : const std::vector< std::size_t >& ndofel,
303 : : const tk::real dt,
304 : : tk::Fields& R ) const
305 : : {
306 : 38265 : self->rhs( t, pref, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
307 : 38265 : ndofel, dt, R );
308 : 38265 : }
309 : :
310 : : //! Evaluate the adaptive indicator and mark the ndof for each element
311 : 1636 : void eval_ndof( std::size_t nunk,
312 : : const tk::UnsMesh::Coords& coord,
313 : : const std::vector< std::size_t >& inpoel,
314 : : const inciter::FaceData& fd,
315 : : const tk::Fields& unk,
316 : : const tk::Fields& prim,
317 : : inciter::ctr::PrefIndicatorType indicator,
318 : : std::size_t ndof,
319 : : std::size_t ndofmax,
320 : : tk::real tolref,
321 : : std::vector< std::size_t >& ndofel ) const
322 : : {
323 : 1636 : self->eval_ndof( nunk, coord, inpoel, fd, unk, prim, indicator, ndof,
324 : 1636 : ndofmax, tolref, ndofel );
325 : 1636 : }
326 : :
327 : : //! Public interface for computing the minimum time step size
328 : 2830 : tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
329 : : const std::vector< std::size_t >& inpoel,
330 : : const inciter::FaceData& fd,
331 : : const tk::Fields& geoFace,
332 : : const tk::Fields& geoElem,
333 : : const std::vector< std::size_t >& ndofel,
334 : : const tk::Fields& U,
335 : : const tk::Fields& P,
336 : : const std::size_t nielem ) const
337 : 5660 : { return self->dt( coord, inpoel, fd, geoFace, geoElem, ndofel, U,
338 : 2830 : P, nielem ); }
339 : :
340 : : //! Public interface for elastic energy balance
341 : 0 : void balance_plastic_energy( std::size_t e,
342 : : std::vector< tk::real > x_star,
343 : : std::vector< tk::real > x,
344 : : tk::Fields& U ) const
345 [ - - ][ - - ]: 0 : { return self->balance_plastic_energy(e, x_star, x, U); }
346 : :
347 : : //! Public interface for computing stiff terms for an element
348 : 0 : void stiff_rhs( std::size_t e,
349 : : const tk::Fields& geoElem,
350 : : const std::vector< std::size_t >& inpoel,
351 : : const tk::UnsMesh::Coords& coord,
352 : : const tk::Fields& U,
353 : : const tk::Fields& P,
354 : : const std::vector< std::size_t >& ndofel,
355 : : tk::Fields& R ) const
356 : 0 : { return self->stiff_rhs( e, geoElem, inpoel, coord, U, P, ndofel, R); }
357 : :
358 : : //! Public interface to returning maps of output var functions
359 : 2650 : std::map< std::string, tk::GetVarFn > OutVarFn() const
360 : 2650 : { return self->OutVarFn(); }
361 : :
362 : : //! Public interface to returning analytic field output labels
363 : 684 : std::vector< std::string > analyticFieldNames() const
364 : 684 : { return self->analyticFieldNames(); }
365 : :
366 : : //! Public interface to returning time history field output labels
367 : 0 : std::vector< std::string > histNames() const { return self->histNames(); }
368 : :
369 : : //! Public interface to returning variable names
370 : 76 : std::vector< std::string > names() const { return self->names(); }
371 : :
372 : : //! Public interface to returning surface field output
373 : : std::vector< std::vector< tk::real > >
374 : : surfOutput( const std::map< int, std::vector< std::size_t > >& bnd,
375 : : tk::Fields& U ) const
376 : : { return self->surfOutput( bnd, U ); }
377 : :
378 : : //! Public interface to return point history output
379 : : std::vector< std::vector< tk::real > >
380 : 0 : histOutput( const std::vector< HistData >& h,
381 : : const std::vector< std::size_t >& inpoel,
382 : : const tk::UnsMesh::Coords& coord,
383 : : const tk::Fields& U,
384 : : const tk::Fields& P ) const
385 : 0 : { return self->histOutput( h, inpoel, coord, U, P ); }
386 : :
387 : : //! Public interface to returning analytic solution
388 : : tk::InitializeFn::result_type
389 : 669522 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
390 : 669522 : { return self->analyticSolution( xi, yi, zi, t ); }
391 : :
392 : : //! Public interface to returning the analytic solution for conserved vars
393 : : tk::InitializeFn::result_type
394 : 4508294 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
395 : 4508294 : { return self->solution( xi, yi, zi, t ); }
396 : :
397 : : //! Public interface to returning the specific total energy
398 : : tk::real
399 : 2553118 : sp_totalenergy( std::size_t e, const tk::Fields& unk ) const
400 : 2553118 : { return self->sp_totalenergy( e, unk ); }
401 : :
402 : : //! Copy assignment
403 : : DGPDE& operator=( const DGPDE& x )
404 : : { DGPDE tmp(x); *this = std::move(tmp); return *this; }
405 : : //! Copy constructor
406 : : DGPDE( const DGPDE& x ) : self( x.self->copy() ) {}
407 : : //! Move assignment
408 : : DGPDE& operator=( DGPDE&& ) noexcept = default;
409 : : //! Move constructor
410 : 367 : DGPDE( DGPDE&& ) noexcept = default;
411 : :
412 : : private:
413 : : //! \brief Concept is a pure virtual base class specifying the requirements
414 : : //! of polymorphic objects deriving from it
415 : : struct Concept {
416 : 367 : Concept() = default;
417 : 0 : Concept( const Concept& ) = default;
418 : 367 : virtual ~Concept() = default;
419 : : virtual Concept* copy() const = 0;
420 : : virtual std::size_t nprim() const = 0;
421 : : virtual std::size_t nmat() const = 0;
422 : : virtual void numEquationDofs(std::vector< std::size_t >&) const = 0;
423 : : virtual std::size_t nstiffeq() const = 0;
424 : : virtual std::size_t nnonstiffeq() const = 0;
425 : : virtual void setStiffEqIdx( std::vector< std::size_t >& ) const = 0;
426 : : virtual void setNonStiffEqIdx( std::vector< std::size_t >& ) const = 0;
427 : : virtual void IcBoxElems( const tk::Fields&,
428 : : std::size_t,
429 : : std::vector< std::unordered_set< std::size_t > >& ) const = 0;
430 : : virtual void initialize(
431 : : const tk::Fields&,
432 : : const std::vector< std::size_t >&,
433 : : const tk::UnsMesh::Coords&,
434 : : const std::vector< std::unordered_set< std::size_t > >&,
435 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
436 : : tk::Fields&,
437 : : tk::real,
438 : : const std::size_t nielem ) const = 0;
439 : : virtual void computeDensityConstr( std::size_t nelem,
440 : : tk::Fields& unk,
441 : : std::vector< tk::real >& densityConstr)
442 : : const = 0;
443 : : virtual void lhs( const tk::Fields&, tk::Fields& ) const = 0;
444 : : virtual void updateInterfaceCells( tk::Fields&,
445 : : std::size_t,
446 : : std::vector< std::size_t >&,
447 : : std::vector< std::size_t >& ) const = 0;
448 : : virtual void updatePrimitives( const tk::Fields&,
449 : : const tk::Fields&,
450 : : const tk::Fields&,
451 : : tk::Fields&,
452 : : std::size_t,
453 : : const std::vector< std::size_t >& ) const = 0;
454 : : virtual void cleanTraceMaterial( tk::real,
455 : : const tk::Fields&,
456 : : tk::Fields&,
457 : : tk::Fields&,
458 : : std::size_t ) const = 0;
459 : : virtual void reconstruct( tk::real,
460 : : const tk::Fields&,
461 : : const tk::Fields&,
462 : : const inciter::FaceData&,
463 : : const std::map< std::size_t,
464 : : std::vector< std::size_t > >&,
465 : : const std::vector< std::size_t >&,
466 : : const tk::UnsMesh::Coords&,
467 : : tk::Fields&,
468 : : tk::Fields&,
469 : : const bool,
470 : : const std::vector< std::size_t >& ) const = 0;
471 : : virtual void limit( tk::real,
472 : : const bool,
473 : : const tk::Fields&,
474 : : const tk::Fields&,
475 : : const inciter::FaceData&,
476 : : const std::map< std::size_t,
477 : : std::vector< std::size_t > >&,
478 : : const std::vector< std::size_t >&,
479 : : const tk::UnsMesh::Coords&,
480 : : const std::vector< std::size_t >&,
481 : : const std::vector< std::size_t >&,
482 : : const std::unordered_map< std::size_t, std::size_t >&,
483 : : const std::vector< std::vector<tk::real> >&,
484 : : const std::vector< std::vector<tk::real> >&,
485 : : const std::vector< std::vector<tk::real> >&,
486 : : tk::Fields&,
487 : : tk::Fields&,
488 : : std::vector< std::size_t >& ) const = 0;
489 : : virtual void CPL( const tk::Fields&,
490 : : const tk::Fields&,
491 : : const std::vector< std::size_t >&,
492 : : const tk::UnsMesh::Coords&,
493 : : tk::Fields&,
494 : : std::size_t ) const = 0;
495 : : virtual std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
496 : : const tk::Fields&,
497 : : std::size_t ) const = 0;
498 : : virtual void rhs( tk::real,
499 : : const bool,
500 : : const tk::Fields&,
501 : : const tk::Fields&,
502 : : const inciter::FaceData&,
503 : : const std::vector< std::size_t >&,
504 : : const std::vector< std::unordered_set< std::size_t > >&,
505 : : const tk::UnsMesh::Coords&,
506 : : const tk::Fields&,
507 : : const tk::Fields&,
508 : : const std::vector< std::size_t >&,
509 : : const tk::real,
510 : : tk::Fields& ) const = 0;
511 : : virtual void resetAdapSol( const inciter::FaceData&,
512 : : tk::Fields&,
513 : : tk::Fields&,
514 : : const std::vector< std::size_t >& )
515 : : const = 0;
516 : : virtual void eval_ndof( std::size_t,
517 : : const tk::UnsMesh::Coords&,
518 : : const std::vector< std::size_t >&,
519 : : const inciter::FaceData&,
520 : : const tk::Fields&,
521 : : const tk::Fields&,
522 : : inciter::ctr::PrefIndicatorType,
523 : : std::size_t,
524 : : std::size_t,
525 : : tk::real,
526 : : std::vector< std::size_t >& ) const = 0;
527 : : virtual tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
528 : : const std::vector< std::size_t >&,
529 : : const inciter::FaceData&,
530 : : const tk::Fields&,
531 : : const tk::Fields&,
532 : : const std::vector< std::size_t >&,
533 : : const tk::Fields&,
534 : : const tk::Fields&,
535 : : const std::size_t ) const = 0;
536 : : virtual void balance_plastic_energy( std::size_t,
537 : : std::vector< tk::real >,
538 : : std::vector< tk::real >,
539 : : tk::Fields& ) const = 0;
540 : : virtual void stiff_rhs( std::size_t,
541 : : const tk::Fields&,
542 : : const std::vector< std::size_t >&,
543 : : const tk::UnsMesh::Coords&,
544 : : const tk::Fields&,
545 : : const tk::Fields&,
546 : : const std::vector< std::size_t >&,
547 : : tk::Fields& ) const = 0;
548 : : virtual std::map< std::string, tk::GetVarFn > OutVarFn() const = 0;
549 : : virtual std::vector< std::string > analyticFieldNames() const = 0;
550 : : virtual std::vector< std::string > histNames() const = 0;
551 : : virtual std::vector< std::string > names() const = 0;
552 : : virtual std::vector< std::vector< tk::real > > surfOutput(
553 : : const std::map< int, std::vector< std::size_t > >&,
554 : : tk::Fields& ) const = 0;
555 : : virtual std::vector< std::vector< tk::real > > histOutput(
556 : : const std::vector< HistData >&,
557 : : const std::vector< std::size_t >&,
558 : : const tk::UnsMesh::Coords&,
559 : : const tk::Fields&,
560 : : const tk::Fields& ) const = 0;
561 : : virtual tk::InitializeFn::result_type analyticSolution(
562 : : tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
563 : : virtual tk::InitializeFn::result_type solution(
564 : : tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
565 : : virtual tk::real sp_totalenergy(
566 : : std::size_t, const tk::Fields& ) const = 0;
567 : : };
568 : :
569 : : //! \brief Model models the Concept above by deriving from it and overriding
570 : : //! the virtual functions required by Concept
571 : : template< typename T >
572 : : struct Model : Concept {
573 : 367 : explicit Model( T x ) : data( std::move(x) ) {}
574 [ - - ]: 0 : Concept* copy() const override { return new Model( *this ); }
575 : 671 : std::size_t nprim() const override
576 : 671 : { return data.nprim(); }
577 : 0 : std::size_t nmat() const override
578 : 0 : { return data.nmat(); }
579 : 671 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const override
580 : 671 : { data.numEquationDofs(numEqDof); }
581 : 2684 : std::size_t nstiffeq() const override
582 : 2684 : { return data.nstiffeq(); }
583 : 1342 : std::size_t nnonstiffeq() const override
584 : 1342 : { return data.nnonstiffeq(); }
585 : 0 : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const override
586 : 0 : { data.setStiffEqIdx(stiffEqIdx); }
587 : 0 : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const override
588 : 0 : { data.setNonStiffEqIdx(nonStiffEqIdx); }
589 : 671 : void IcBoxElems( const tk::Fields& geoElem,
590 : : std::size_t nielem,
591 : : std::vector< std::unordered_set< std::size_t > >& inbox )
592 : 671 : const override { data.IcBoxElems( geoElem, nielem, inbox ); }
593 : 803 : void initialize(
594 : : const tk::Fields& L,
595 : : const std::vector< std::size_t >& inpoel,
596 : : const tk::UnsMesh::Coords& coord,
597 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
598 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
599 : : elemblkid,
600 : : tk::Fields& unk,
601 : : tk::real t,
602 : : const std::size_t nielem )
603 : 803 : const override { data.initialize( L, inpoel, coord, inbox, elemblkid, unk,
604 : 803 : t, nielem ); }
605 : 1325 : void computeDensityConstr( std::size_t nelem,
606 : : tk::Fields& unk,
607 : : std::vector< tk::real >& densityConstr)
608 : : const override
609 : 1325 : { data.computeDensityConstr( nelem, unk, densityConstr ); }
610 : 803 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const override
611 : 803 : { data.lhs( geoElem, l ); }
612 : 38265 : void updateInterfaceCells( tk::Fields& unk,
613 : : std::size_t nielem,
614 : : std::vector< std::size_t >& ndofel,
615 : : std::vector< std::size_t >& interface )
616 : 38265 : const override { data.updateInterfaceCells( unk, nielem, ndofel, interface ); }
617 : 38936 : void updatePrimitives( const tk::Fields& unk,
618 : : const tk::Fields& L,
619 : : const tk::Fields& geoElem,
620 : : tk::Fields& prim,
621 : : std::size_t nielem,
622 : : const std::vector< std::size_t >& ndofel )
623 : : const override {
624 : 38936 : data.updatePrimitives( unk, L, geoElem, prim, nielem, ndofel );
625 : 38936 : }
626 : 38265 : void cleanTraceMaterial( tk::real t,
627 : : const tk::Fields& geoElem,
628 : : tk::Fields& unk,
629 : : tk::Fields& prim,
630 : : std::size_t nielem )
631 : 38265 : const override { data.cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
632 : 24840 : void reconstruct( tk::real t,
633 : : const tk::Fields& geoFace,
634 : : const tk::Fields& geoElem,
635 : : const inciter::FaceData& fd,
636 : : const std::map< std::size_t,
637 : : std::vector< std::size_t > >& esup,
638 : : const std::vector< std::size_t >& inpoel,
639 : : const tk::UnsMesh::Coords& coord,
640 : : tk::Fields& U,
641 : : tk::Fields& P,
642 : : const bool pref,
643 : : const std::vector< std::size_t >& ndofel )const override
644 : : {
645 : 24840 : data.reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P,
646 : : pref, ndofel );
647 : 24840 : }
648 : 24840 : void limit( tk::real t,
649 : : const bool pref,
650 : : const tk::Fields& geoFace,
651 : : const tk::Fields& geoElem,
652 : : const inciter::FaceData& fd,
653 : : const std::map< std::size_t, std::vector< std::size_t > >&
654 : : esup,
655 : : const std::vector< std::size_t >& inpoel,
656 : : const tk::UnsMesh::Coords& coord,
657 : : const std::vector< std::size_t >& ndofel,
658 : : const std::vector< std::size_t >& gid,
659 : : const std::unordered_map< std::size_t, std::size_t >& bid,
660 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
661 : : const std::vector< std::vector<tk::real> >& pNodalExtrm,
662 : : const std::vector< std::vector<tk::real> >& mtInv,
663 : : tk::Fields& U,
664 : : tk::Fields& P,
665 : : std::vector< std::size_t >& shockmarker ) const override
666 : : {
667 : 24840 : data.limit( t, pref, geoFace, geoElem, fd, esup, inpoel, coord, ndofel, gid,
668 : : bid, uNodalExtrm, pNodalExtrm, mtInv, U, P, shockmarker );
669 : 24840 : }
670 : 20790 : void CPL( const tk::Fields& prim,
671 : : const tk::Fields& geoElem,
672 : : const std::vector< std::size_t >& inpoel,
673 : : const tk::UnsMesh::Coords& coord,
674 : : tk::Fields& unk,
675 : : std::size_t nielem ) const override
676 : : {
677 : 20790 : data.CPL( prim, geoElem, inpoel, coord, unk, nielem );
678 : 20790 : }
679 : 1770 : void resetAdapSol( const inciter::FaceData& fd,
680 : : tk::Fields& unk,
681 : : tk::Fields& prim,
682 : : const std::vector< std::size_t >& ndofel ) const override
683 : : {
684 : 1770 : data.resetAdapSol( fd, unk, prim, ndofel );
685 : 1770 : }
686 : 0 : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
687 : : const tk::Fields& U,
688 : : std::size_t nielem ) const override
689 : : {
690 : 0 : return data.cellAvgDeformGrad( U, nielem );
691 : : }
692 : 38265 : void rhs(
693 : : tk::real t,
694 : : const bool pref,
695 : : const tk::Fields& geoFace,
696 : : const tk::Fields& geoElem,
697 : : const inciter::FaceData& fd,
698 : : const std::vector< std::size_t >& inpoel,
699 : : const std::vector< std::unordered_set< std::size_t > >& boxelems,
700 : : const tk::UnsMesh::Coords& coord,
701 : : const tk::Fields& U,
702 : : const tk::Fields& P,
703 : : const std::vector< std::size_t >& ndofel,
704 : : const tk::real dt,
705 : : tk::Fields& R ) const override
706 : : {
707 : 38265 : data.rhs( t, pref, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
708 : : ndofel, dt, R );
709 : 38265 : }
710 : 1636 : void eval_ndof( std::size_t nunk,
711 : : const tk::UnsMesh::Coords& coord,
712 : : const std::vector< std::size_t >& inpoel,
713 : : const inciter::FaceData& fd,
714 : : const tk::Fields& unk,
715 : : const tk::Fields& prim,
716 : : inciter::ctr::PrefIndicatorType indicator,
717 : : std::size_t ndof,
718 : : std::size_t ndofmax,
719 : : tk::real tolref,
720 : : std::vector< std::size_t >& ndofel ) const override
721 : 1636 : { data.eval_ndof( nunk, coord, inpoel, fd, unk, prim, indicator, ndof,
722 : 1636 : ndofmax, tolref, ndofel ); }
723 : 2830 : tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
724 : : const std::vector< std::size_t >& inpoel,
725 : : const inciter::FaceData& fd,
726 : : const tk::Fields& geoFace,
727 : : const tk::Fields& geoElem,
728 : : const std::vector< std::size_t >& ndofel,
729 : : const tk::Fields& U,
730 : : const tk::Fields& P,
731 : : const std::size_t nielem ) const override
732 : : { return data.dt( coord, inpoel, fd, geoFace, geoElem, ndofel,
733 : 2830 : U, P, nielem ); }
734 : 0 : void balance_plastic_energy( std::size_t e,
735 : : std::vector< tk::real > x_star,
736 : : std::vector< tk::real > x,
737 : : tk::Fields& U ) const override
738 [ - - ][ - - ]: 0 : { return data.balance_plastic_energy( e, x_star, x, U); }
739 : 0 : void stiff_rhs( std::size_t e,
740 : : const tk::Fields& geoElem,
741 : : const std::vector< std::size_t >& inpoel,
742 : : const tk::UnsMesh::Coords& coord,
743 : : const tk::Fields& U,
744 : : const tk::Fields& P,
745 : : const std::vector< std::size_t >& ndofel,
746 : : tk::Fields& R ) const override
747 : 0 : { return data.stiff_rhs( e, geoElem, inpoel, coord, U, P, ndofel, R ); }
748 : 2650 : std::map< std::string, tk::GetVarFn > OutVarFn() const override
749 : 2650 : { return data.OutVarFn(); }
750 : 684 : std::vector< std::string > analyticFieldNames() const override
751 : 684 : { return data.analyticFieldNames(); }
752 : 0 : std::vector< std::string > histNames() const override
753 : 0 : { return data.histNames(); }
754 : 76 : std::vector< std::string > names() const override
755 : 76 : { return data.names(); }
756 : 0 : std::vector< std::vector< tk::real > > surfOutput(
757 : : const std::map< int, std::vector< std::size_t > >& bnd,
758 : : tk::Fields& U ) const override
759 : 0 : { return data.surfOutput( bnd, U ); }
760 : 0 : std::vector< std::vector< tk::real > > histOutput(
761 : : const std::vector< HistData >& h,
762 : : const std::vector< std::size_t >& inpoel,
763 : : const tk::UnsMesh::Coords& coord,
764 : : const tk::Fields& U,
765 : : const tk::Fields& P ) const override
766 : 0 : { return data.histOutput( h, inpoel, coord, U, P ); }
767 : : tk::InitializeFn::result_type
768 : 669522 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
769 : 669522 : const override { return data.analyticSolution( xi, yi, zi, t ); }
770 : : tk::InitializeFn::result_type
771 : 4508294 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
772 : 4508294 : const override { return data.solution( xi, yi, zi, t ); }
773 : 2553118 : tk::real sp_totalenergy( std::size_t e, const tk::Fields& unk )
774 : 2553118 : const override { return data.sp_totalenergy( e, unk ); }
775 : : T data;
776 : : };
777 : :
778 : : std::unique_ptr< Concept > self; //!< Base pointer used polymorphically
779 : : };
780 : :
781 : : } // inciter::
782 : :
783 : : #endif // DGPDE_h
|
