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_FE__
34 : #define __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FE__
35 :
36 : // library intern headers
37 : #include "../convection_diffusion_base.h"
38 : #include "lib_disc/spatial_disc/disc_util/conv_shape_interface.h"
39 :
40 : namespace ug{
41 : namespace ConvectionDiffusionPlugin{
42 :
43 : // \ingroup lib_disc_elem_disc
44 : /// \addtogroup convection_diffusion
45 : /// \{
46 :
47 : /// Discretization for the Convection-Diffusion Equation
48 : /**
49 : * This is a FE discretization of the convection-diffusion equation.
50 : * Cf. ConvectionDiffusionBase base class for the problem setting.
51 : *
52 : * \see ConvectionDiffusionBase for the problem settings
53 : *
54 : * \tparam TDomain Domain
55 : */
56 : template< typename TDomain>
57 : class ConvectionDiffusionFE : public ConvectionDiffusionBase<TDomain>
58 : {
59 : private:
60 : /// Base class type
61 : typedef ConvectionDiffusionBase<TDomain> base_type;
62 :
63 : /// Own type
64 : typedef ConvectionDiffusionFE<TDomain> this_type;
65 :
66 : /// error estimator type
67 : typedef SideAndElemErrEstData<TDomain> err_est_type;
68 :
69 : public:
70 : /// World dimension
71 : static const int dim = base_type::dim;
72 :
73 : public:
74 : /// Constructor
75 : ConvectionDiffusionFE(const char* functions, const char* subsets);
76 :
77 : /// sets the quad order
78 : void set_quad_order(size_t order);
79 :
80 : private:
81 : /// prepares the loop over all elements
82 : /**
83 : * This method prepares the loop over all elements. It resizes the Position
84 : * array for the corner coordinates and schedules the local ip positions
85 : * at the data imports.
86 : */
87 : template <typename TElem, typename TFEGeom>
88 : void prep_elem_loop(const ReferenceObjectID roid, const int si);
89 :
90 : /// prepares the element for assembling
91 : /**
92 : * This methods prepares an element for the assembling. The Positions of
93 : * the Element Corners are read and the Finite Volume Geometry is updated.
94 : * The global ip positions are scheduled at the data imports.
95 : */
96 : template <typename TElem, typename TFEGeom>
97 : void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[]);
98 :
99 : /// finishes the loop over all elements
100 : template <typename TElem, typename TFEGeom>
101 : void fsh_elem_loop();
102 :
103 : /// assembles the local stiffness matrix using a finite volume scheme
104 : template <typename TElem, typename TFEGeom>
105 : void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
106 :
107 : /// assembles the local mass matrix using a finite volume scheme
108 : template <typename TElem, typename TFEGeom>
109 : void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
110 :
111 : /// assembles the stiffness part of the local defect
112 : template <typename TElem, typename TFEGeom>
113 : void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
114 :
115 : /// assembles the mass part of the local defect
116 : template <typename TElem, typename TFEGeom>
117 : void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
118 :
119 : /// assembles the local right hand side
120 : template <typename TElem, typename TFEGeom>
121 : void add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[]);
122 :
123 : /// prepares the loop over all elements of one type for the computation of the error estimator
124 : template <typename TElem, typename TFEGeom>
125 : void prep_err_est_elem_loop(const ReferenceObjectID roid, const int si);
126 :
127 : /// prepares the element for assembling the error estimator
128 : template <typename TElem, typename TFEGeom>
129 : void prep_err_est_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[]);
130 :
131 : /// computes the error estimator contribution for one element
132 : template <typename TElem, typename TFEGeom>
133 : void compute_err_est_A_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
134 :
135 : /// computes the error estimator contribution for one element
136 : template <typename TElem, typename TFEGeom>
137 : void compute_err_est_M_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
138 :
139 : /// computes the error estimator contribution for one element
140 : template <typename TElem, typename TFEGeom>
141 : void compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale);
142 :
143 : /// postprocesses the loop over all elements of one type in the computation of the error estimator
144 : template <typename TElem, typename TFEGeom>
145 : void fsh_err_est_elem_loop();
146 :
147 : protected:
148 : /// computes the linearized defect w.r.t to the velocity
149 : template <typename TElem, typename TFEGeom>
150 : void lin_def_velocity(const LocalVector& u,
151 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
152 : const size_t nip);
153 :
154 : /// computes the linearized defect w.r.t to the velocity
155 : template <typename TElem, typename TFEGeom>
156 : void lin_def_diffusion(const LocalVector& u,
157 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
158 : const size_t nip);
159 :
160 : /// computes the linearized defect w.r.t to the flux
161 : template <typename TElem, typename TFEGeom>
162 : void lin_def_flux(const LocalVector& u,
163 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
164 : const size_t nip);
165 :
166 : /// computes the linearized defect w.r.t to the reaction
167 : template <typename TElem, typename TFEGeom>
168 : void lin_def_reaction(const LocalVector& u,
169 : std::vector<std::vector<number> > vvvLinDef[],
170 : const size_t nip);
171 :
172 : /// computes the linearized defect w.r.t to the reaction
173 : template <typename TElem, typename TFEGeom>
174 : void lin_def_reaction_rate(const LocalVector& u,
175 : std::vector<std::vector<number> > vvvLinDef[],
176 : const size_t nip);
177 :
178 : /// computes the linearized defect w.r.t to the source term
179 : template <typename TElem, typename TFEGeom>
180 : void lin_def_source(const LocalVector& u,
181 : std::vector<std::vector<number> > vvvLinDef[],
182 : const size_t nip);
183 :
184 : /// computes the linearized defect w.r.t to the vector source term
185 : template <typename TElem, typename TFEGeom>
186 : void lin_def_vector_source(const LocalVector& u,
187 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
188 : const size_t nip);
189 :
190 : /// computes the linearized defect w.r.t to the mass scale term
191 : template <typename TElem, typename TFEGeom>
192 : void lin_def_mass_scale(const LocalVector& u,
193 : std::vector<std::vector<number> > vvvLinDef[],
194 : const size_t nip);
195 :
196 : /// computes the linearized defect w.r.t to the mass scale term
197 : template <typename TElem, typename TFEGeom>
198 : void lin_def_mass(const LocalVector& u,
199 : std::vector<std::vector<number> > vvvLinDef[],
200 : const size_t nip);
201 :
202 : private:
203 : /// abbreviation for the local solution
204 : static const size_t _C_ = 0;
205 :
206 : using base_type::m_imDiffusion;
207 : using base_type::m_imVelocity;
208 : using base_type::m_imFlux;
209 : using base_type::m_imSource;
210 : using base_type::m_imSourceExpl;
211 : using base_type::m_imVectorSource;
212 : using base_type::m_imReactionRate;
213 : using base_type::m_imReactionRateExpl;
214 : using base_type::m_imReaction;
215 : using base_type::m_imReactionExpl;
216 : using base_type::m_imMassScale;
217 : using base_type::m_imMass;
218 :
219 : using base_type::m_exGrad;
220 : using base_type::m_exValue;
221 :
222 : protected:
223 : /// computes the concentration
224 : template <typename TElem, typename TFEGeom>
225 : void ex_value(number vValue[],
226 : const MathVector<dim> vGlobIP[],
227 : number time, int si,
228 : const LocalVector& u,
229 : GridObject* elem,
230 : const MathVector<dim> vCornerCoords[],
231 : const MathVector<TFEGeom::dim> vLocIP[],
232 : const size_t nip,
233 : bool bDeriv,
234 : std::vector<std::vector<number> > vvvDeriv[]);
235 :
236 : /// computes the gradient of the concentration
237 : template <typename TElem, typename TFEGeom>
238 : void ex_grad(MathVector<dim> vValue[],
239 : const MathVector<dim> vGlobIP[],
240 : number time, int si,
241 : const LocalVector& u,
242 : GridObject* elem,
243 : const MathVector<dim> vCornerCoords[],
244 : const MathVector<TFEGeom::dim> vLocIP[],
245 : const size_t nip,
246 : bool bDeriv,
247 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[]);
248 :
249 : public:
250 : /// type of trial space for each function used
251 : virtual void prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid);
252 :
253 : /// returns if hanging nodes are needed
254 : virtual bool use_hanging() const;
255 :
256 : protected:
257 : /// current integration order
258 : bool m_bQuadOrderUserDef;
259 : int m_quadOrder;
260 :
261 : /// current shape function set
262 : LFEID m_lfeID;
263 :
264 : /// register utils
265 : /// \{
266 : void register_all_funcs(const LFEID& lfeid, const int quadOrder);
267 : template <typename TElem, typename TFEGeom> void register_func();
268 : /// \}
269 :
270 : private:
271 : /// struct holding values of shape functions in IPs
272 0 : struct ShapeValues
273 : {
274 : public:
275 0 : void resize(std::size_t nEip, std::size_t nSip, std::size_t _nSh)
276 : {
277 0 : nSh = _nSh;
278 0 : elemVals.resize(nEip);
279 0 : sideVals.resize(nSip);
280 0 : for (std::size_t i = 0; i < nEip; i++) elemVals[i].resize(nSh);
281 0 : for (std::size_t i = 0; i < nSip; i++) sideVals[i].resize(nSh);
282 0 : }
283 0 : number& shapeAtElemIP(std::size_t sh, std::size_t ip) {return elemVals[ip][sh];}
284 0 : number& shapeAtSideIP(std::size_t sh, std::size_t ip) {return sideVals[ip][sh];}
285 0 : number* shapesAtElemIP(std::size_t ip) {return &elemVals[ip][0];}
286 0 : number* shapesAtSideIP(std::size_t ip) {return &sideVals[ip][0];}
287 0 : std::size_t num_sh() {return nSh;}
288 : private:
289 : std::size_t nSh;
290 : std::vector<std::vector<number> > elemVals;
291 : std::vector<std::vector<number> > sideVals;
292 : } m_shapeValues;
293 :
294 : };
295 :
296 : // end group convection_diffusion
297 : /// \}
298 :
299 : } // end ConvectionDiffusionPlugin
300 : } // end namespace ug
301 :
302 :
303 : #endif /*__H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FE__*/
|
