Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__
34 : #define __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__
35 :
36 : // ug4 headers
37 : #include "lib_disc/spatial_disc/disc_util/conv_shape_interface.h"
38 :
39 : // plugin's internal headers
40 : #include "../convection_diffusion_base.h"
41 : #include "../convection_diffusion_sss.h"
42 :
43 : namespace ug{
44 : namespace ConvectionDiffusionPlugin{
45 :
46 : // \ingroup lib_disc_elem_disc
47 : /// \addtogroup convection_diffusion
48 : /// \{
49 :
50 : /// FV Discretization for the Convection-Diffusion Equation
51 : /**
52 : * This is the simple (upwinded) FV discretization of the convection-diffusion equation.
53 : * Cf. ConvectionDiffusionBase base class for the problem setting.
54 : *
55 : * \see ConvectionDiffusionBase for the problem settings
56 : *
57 : * \tparam TDomain Domain
58 : */
59 : template< typename TDomain>
60 : class ConvectionDiffusionFV1 : public ConvectionDiffusionBase<TDomain>
61 : {
62 : private:
63 : /// Base class type
64 : typedef ConvectionDiffusionBase<TDomain> base_type;
65 :
66 : /// Own type
67 : typedef ConvectionDiffusionFV1<TDomain> this_type;
68 :
69 : /// error estimator type
70 : typedef SideAndElemErrEstData<TDomain> err_est_type;
71 :
72 : public:
73 : /// World dimension
74 : static const int dim = base_type::dim;
75 :
76 : public:
77 : /// Constructor
78 : ConvectionDiffusionFV1(const char* functions, const char* subsets);
79 :
80 : /// set the upwind method
81 : /**
82 : * This method sets the upwind method used to upwind the convection.
83 : *
84 : * \param shapes upwind method
85 : */
86 : void set_upwind(SmartPtr<IConvectionShapes<dim> > shapes);
87 :
88 : /// set singular sources and sinks
89 0 : void set_sss_manager(SmartPtr<CDSingularSourcesAndSinks<dim> > sss_mngr) {m_sss_mngr = sss_mngr;}
90 :
91 : /// get singular sources and sinks
92 0 : SmartPtr<CDSingularSourcesAndSinks<dim> > sss_manager() {return m_sss_mngr;}
93 :
94 : /// activates/deactivates the condensed scvf ip's for the FV scheme
95 0 : void set_condensed_FV(bool condensed) {m_bCondensedFV = condensed;}
96 :
97 : /// returns the 'condensed scvf ip' flag
98 0 : bool condensed_FV() {return m_bCondensedFV;}
99 :
100 : private:
101 : /// prepares assembling
102 : virtual void prep_assemble_loop();
103 :
104 : /// prepares the loop over all elements
105 : /**
106 : * This method prepares the loop over all elements. It resizes the Position
107 : * array for the corner coordinates and schedules the local ip positions
108 : * at the data imports.
109 : */
110 : template <typename TElem, typename TFVGeom>
111 : void prep_elem_loop(const ReferenceObjectID roid, const int si);
112 :
113 : /// prepares the element for assembling
114 : /**
115 : * This methods prepares an element for the assembling. The Positions of
116 : * the Element Corners are read and the Finite Volume Geometry is updated.
117 : * The global ip positions are scheduled at the data imports.
118 : */
119 : template <typename TElem, typename TFVGeom>
120 : void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[]);
121 :
122 : /// finishes the loop over all elements
123 : template <typename TElem, typename TFVGeom>
124 : void fsh_elem_loop();
125 :
126 : /// assembles the local stiffness matrix using a finite volume scheme
127 : template <typename TElem, typename TFVGeom>
128 : void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
129 :
130 : /// assembles the local mass matrix using a finite volume scheme
131 : template <typename TElem, typename TFVGeom>
132 : void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
133 :
134 : /// assembles the stiffness part of the local defect
135 : template <typename TElem, typename TFVGeom>
136 : void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
137 :
138 : /// assembles the stiffness part of the local defect explicit reaction, reaction_rate and source
139 : template <typename TElem, typename TFVGeom>
140 : void add_def_A_expl_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
141 :
142 : /// assembles the mass part of the local defect
143 : template <typename TElem, typename TFVGeom>
144 : void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
145 :
146 : /// assembles the local right hand side
147 : template <typename TElem, typename TFVGeom>
148 : void add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[]);
149 :
150 : /// prepares the loop over all elements of one type for the computation of the error estimator
151 : template <typename TElem, typename TFVGeom>
152 : void prep_err_est_elem_loop(const ReferenceObjectID roid, const int si);
153 :
154 : /// prepares the element for assembling the error estimator
155 : template <typename TElem, typename TFVGeom>
156 : void prep_err_est_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
157 :
158 : /// computes the error estimator contribution for one element
159 : template <typename TElem, typename TFVGeom>
160 : void compute_err_est_A_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
161 :
162 : /// computes the error estimator contribution for one element
163 : template <typename TElem, typename TFVGeom>
164 : void compute_err_est_M_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
165 :
166 : /// computes the error estimator contribution for one element
167 : template <typename TElem, typename TFVGeom>
168 : void compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
169 :
170 : /// postprocesses the loop over all elements of one type in the computation of the error estimator
171 : template <typename TElem, typename TFVGeom>
172 : void fsh_err_est_elem_loop();
173 :
174 : private:
175 :
176 : /// adds contributions of a singular source or sink to the matrix
177 : template<typename TElem, typename TFVGeom>
178 : void add_sss_jac_elem
179 : (
180 : LocalMatrix& J, ///< the matrix to update
181 : const LocalVector& u, ///< current solution
182 : GridObject* elem, ///< the element
183 : const TFVGeom& geo, ///< the FV geometry for that element
184 : size_t i, ///< index of the SCV
185 : number flux ///< flux through source/sink (premultiplied by the length for lines)
186 : );
187 :
188 : /// adds contributions of a singular source or sink to the defect
189 : template<typename TElem, typename TFVGeom>
190 : void add_sss_def_elem
191 : (
192 : LocalVector& d, ///< the defect to update
193 : const LocalVector& u, ///< current solution
194 : GridObject* elem, ///< the element
195 : const TFVGeom& geo, ///< the FV geometry for that element
196 : size_t i, ///< index of the SCV
197 : number flux ///< flux through source/sink (premultiplied by the length for lines)
198 : );
199 :
200 : protected:
201 : /// computes the linearized defect w.r.t to the velocity
202 : template <typename TElem, typename TFVGeom>
203 : void lin_def_velocity(const LocalVector& u,
204 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
205 : const size_t nip);
206 :
207 : /// computes the linearized defect w.r.t to the diffusion
208 : template <typename TElem, typename TFVGeom>
209 : void lin_def_diffusion(const LocalVector& u,
210 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
211 : const size_t nip);
212 :
213 : /// computes the linearized defect w.r.t to the flux
214 : template <typename TElem, typename TFVGeom>
215 : void lin_def_flux(const LocalVector& u,
216 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
217 : const size_t nip);
218 :
219 : /// computes the linearized defect w.r.t to the reaction
220 : template <typename TElem, typename TFVGeom>
221 : void lin_def_reaction(const LocalVector& u,
222 : std::vector<std::vector<number> > vvvLinDef[],
223 : const size_t nip);
224 :
225 : /// computes the linearized defect w.r.t to the reaction
226 : template <typename TElem, typename TFVGeom>
227 : void lin_def_reaction_rate(const LocalVector& u,
228 : std::vector<std::vector<number> > vvvLinDef[],
229 : const size_t nip);
230 :
231 : /// computes the linearized defect w.r.t to the source term
232 : template <typename TElem, typename TFVGeom>
233 : void lin_def_source(const LocalVector& u,
234 : std::vector<std::vector<number> > vvvLinDef[],
235 : const size_t nip);
236 :
237 : /// computes the linearized defect w.r.t to the vector source term
238 : template <typename TElem, typename TFVGeom>
239 : void lin_def_vector_source(const LocalVector& u,
240 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
241 : const size_t nip);
242 :
243 : /// computes the linearized defect w.r.t to the mass scale term
244 : template <typename TElem, typename TFVGeom>
245 : void lin_def_mass_scale(const LocalVector& u,
246 : std::vector<std::vector<number> > vvvLinDef[],
247 : const size_t nip);
248 :
249 : /// computes the linearized defect w.r.t to the mass scale term
250 : template <typename TElem, typename TFVGeom>
251 : void lin_def_mass(const LocalVector& u,
252 : std::vector<std::vector<number> > vvvLinDef[],
253 : const size_t nip);
254 :
255 : private:
256 : /// abbreviation for the local solution
257 : static const size_t _C_ = 0;
258 :
259 : /// singular sources and sinks manager
260 : SmartPtr<CDSingularSourcesAndSinks<dim> > m_sss_mngr;
261 :
262 : using base_type::m_imDiffusion;
263 : using base_type::m_imVelocity;
264 : using base_type::m_imFlux;
265 : using base_type::m_imSource;
266 : using base_type::m_imSourceExpl;
267 : using base_type::m_imVectorSource;
268 : using base_type::m_imReactionRate;
269 : using base_type::m_imReactionRateExpl;
270 : using base_type::m_imReaction;
271 : using base_type::m_imReactionExpl;
272 : using base_type::m_imMassScale;
273 : using base_type::m_imMass;
274 :
275 : using base_type::m_exGrad;
276 : using base_type::m_exValue;
277 :
278 : protected:
279 : /// method to compute the upwind shapes
280 : SmartPtr<IConvectionShapes<dim> > m_spConvShape;
281 :
282 : /// returns the updated convection shapes
283 : typedef IConvectionShapes<dim> conv_shape_type;
284 : const IConvectionShapes<dim>& get_updated_conv_shapes(const FVGeometryBase& geo, bool compute_deriv);
285 :
286 : /// computes the concentration
287 : template <typename TElem, typename TFVGeom>
288 : void ex_value(number vValue[],
289 : const MathVector<dim> vGlobIP[],
290 : number time, int si,
291 : const LocalVector& u,
292 : GridObject* elem,
293 : const MathVector<dim> vCornerCoords[],
294 : const MathVector<TFVGeom::dim> vLocIP[],
295 : const size_t nip,
296 : bool bDeriv,
297 : std::vector<std::vector<number> > vvvDeriv[]);
298 :
299 : /// computes the gradient of the concentration
300 : template <typename TElem, typename TFVGeom>
301 : void ex_grad(MathVector<dim> vValue[],
302 : const MathVector<dim> vGlobIP[],
303 : number time, int si,
304 : const LocalVector& u,
305 : GridObject* elem,
306 : const MathVector<dim> vCornerCoords[],
307 : const MathVector<TFVGeom::dim> vLocIP[],
308 : const size_t nip,
309 : bool bDeriv,
310 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[]);
311 :
312 : public:
313 : /// type of trial space for each function used
314 : virtual void prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid);
315 :
316 : /// returns if hanging nodes are needed
317 : virtual bool use_hanging() const;
318 :
319 : protected:
320 : /// current regular grid flag
321 : bool m_bNonRegularGrid;
322 :
323 : /// if to use the 'condensed' FV scvf ip's
324 : bool m_bCondensedFV;
325 :
326 : /// register utils
327 : /// \{
328 : void register_all_funcs(bool bHang);
329 : template <typename TElem> void register_func_for_(bool bHang);
330 : template <typename TElem, typename TFVGeom> void register_func();
331 : struct RegisterLocalDiscr
332 : {
333 : this_type * m_pThis; bool m_bHang;
334 0 : RegisterLocalDiscr(this_type * pThis, bool bHang) : m_pThis(pThis), m_bHang(bHang) {}
335 : template< typename TElem > void operator() (TElem &)
336 0 : {m_pThis->register_func_for_<TElem> (m_bHang);}
337 : };
338 : /// \}
339 :
340 : private:
341 : /// struct holding values of shape functions in IPs
342 0 : struct ShapeValues
343 : {
344 : public:
345 0 : void resize(std::size_t nEip, std::size_t nSip, std::size_t _nSh)
346 : {
347 0 : nSh = _nSh;
348 0 : elemVals.resize(nEip);
349 0 : sideVals.resize(nSip);
350 0 : for (std::size_t i = 0; i < nEip; i++) elemVals[i].resize(nSh);
351 0 : for (std::size_t i = 0; i < nSip; i++) sideVals[i].resize(nSh);
352 0 : }
353 0 : number& shapeAtElemIP(std::size_t sh, std::size_t ip) {return elemVals[ip][sh];}
354 0 : number& shapeAtSideIP(std::size_t sh, std::size_t ip) {return sideVals[ip][sh];}
355 0 : number* shapesAtElemIP(std::size_t ip) {return &elemVals[ip][0];}
356 0 : number* shapesAtSideIP(std::size_t ip) {return &sideVals[ip][0];}
357 0 : std::size_t num_sh() {return nSh;}
358 : private:
359 : std::size_t nSh;
360 : std::vector<std::vector<number> > elemVals;
361 : std::vector<std::vector<number> > sideVals;
362 : } m_shapeValues;
363 : };
364 :
365 : // end group convection_diffusion
366 : /// \}
367 :
368 : } // end ConvectionDiffusionPlugin
369 : } // end namespace ug
370 :
371 :
372 : #endif /*__H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FV1__*/
|
