Line data Source code
1 : /*
2 : * Copyright (c) 2013-2018: G-CSC, Goethe University Frankfurt
3 : * Author: Dmitry Logashenko
4 : * Based on the modules by Andreas Vogel
5 : *
6 : * This file is part of UG4.
7 : *
8 : * UG4 is free software: you can redistribute it and/or modify it under the
9 : * terms of the GNU Lesser General Public License version 3 (as published by the
10 : * Free Software Foundation) with the following additional attribution
11 : * requirements (according to LGPL/GPL v3 §7):
12 : *
13 : * (1) The following notice must be displayed in the Appropriate Legal Notices
14 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
15 : *
16 : * (2) The following notice must be displayed at a prominent place in the
17 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
18 : *
19 : * (3) The following bibliography is recommended for citation and must be
20 : * preserved in all covered files:
21 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
22 : * parallel geometric multigrid solver on hierarchically distributed grids.
23 : * Computing and visualization in science 16, 4 (2013), 151-164"
24 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
25 : * flexible software system for simulating pde based models on high performance
26 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
27 : *
28 : * This program is distributed in the hope that it will be useful,
29 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
30 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 : * GNU Lesser General Public License for more details.
32 : */
33 :
34 : #ifndef __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FRACT_FV1__
35 : #define __H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FRACT_FV1__
36 :
37 : // ug4 headers
38 : #include "lib_grid/algorithms/deg_layer_mngr.h"
39 : #include "lib_disc/spatial_disc/disc_util/conv_shape_interface.h"
40 :
41 : // plugin's internal headers
42 : #include "../convection_diffusion_base.h"
43 : #include "../convection_diffusion_sss.h"
44 :
45 : namespace ug{
46 : namespace ConvectionDiffusionPlugin{
47 :
48 : /// Discretization for the Convection-Diffusion Equation in fractures
49 : /**
50 : * This class implements the IElemDisc interface to provide element local
51 : * assemblings for the convection diffusion equation in the low-dimensional fractures.
52 : * using the degenerated elements. This module is complementary to ConvectionDiffusionFV1.
53 : * Cf. ConvectionDiffusionBase base class for the problem setting.
54 : *
55 : * \see ConvectionDiffusionBase for the problem settings
56 : * \see ConvectionDiffusionFV1 for the FV discretization in the full-dim. domain
57 : *
58 : * \tparam TDomain Domain
59 : */
60 : template<typename TDomain>
61 : class ConvectionDiffusionFractFV1 : public ConvectionDiffusionBase<TDomain>
62 : {
63 : /// Own type
64 : typedef ConvectionDiffusionFractFV1<TDomain> this_type;
65 :
66 : public:
67 :
68 : /// Base class type
69 : typedef ConvectionDiffusionBase<TDomain> base_type;
70 :
71 : /// domain type
72 : typedef typename base_type::domain_type domain_type;
73 :
74 : /// position type
75 : typedef typename base_type::position_type position_type;
76 :
77 : /// World ('full') dimension
78 : static const int dim = base_type::dim;
79 :
80 : /// Manifold ('low') dimension
81 : static const int low_dim = dim - 1;
82 :
83 : /// fracture manager type
84 : typedef DegeneratedLayerManager<dim> fract_manager_type;
85 :
86 : private:
87 :
88 : /// FV geometry for the fractures
89 : typedef DimFV1Geometry<low_dim, dim> TFractFVGeom;
90 :
91 : /// type of the sides of elements (als low-dimensional fracture elements)
92 : typedef typename fract_manager_type::side_type side_type;
93 :
94 : /// max. number of corners of non-degenerated sides
95 : static const size_t maxFractSideCorners = fract_manager_type::maxLayerSideCorners;
96 :
97 : /// abbreviation for local function: brine mass fraction
98 : static const size_t _C_ = 0;
99 :
100 : /// convection shapes type (for the upwind)
101 : typedef IConvectionShapes<dim> conv_shape_type;
102 :
103 : public:
104 : /// Constructor
105 : ConvectionDiffusionFractFV1(const char* functions, const char* subsets);
106 :
107 : /// sets the fracture manager
108 : /**
109 : * This function sets the fracture manager object to recognize the fractures.
110 : */
111 0 : void set_fract_manager (SmartPtr<fract_manager_type> fract_manager) {m_spFractManager = fract_manager;}
112 :
113 : /// set the upwind method along the fracture
114 : /**
115 : * This method sets the upwind method used to upwind the convection
116 : * along the fracture.
117 : *
118 : * \param shapes upwind method
119 : */
120 0 : void set_upwind (SmartPtr<conv_shape_type> shapes) {m_spConvShape = shapes;}
121 :
122 : // set Specific parameters for low-dimensional subdomains (fractures):
123 :
124 0 : void set_aperture(SmartPtr<CplUserData<number, dim> > user)
125 : {
126 0 : m_imAperture.set_data(user);
127 0 : }
128 0 : void set_aperture(number val)
129 : {
130 0 : set_aperture(make_sp(new ConstUserNumber<dim>(val)));
131 0 : }
132 : #ifdef UG_FOR_LUA
133 0 : void set_aperture(const char* fctName)
134 : {
135 0 : set_aperture(LuaUserDataFactory<number, dim>::create(fctName));
136 0 : }
137 0 : void set_aperture(LuaFunctionHandle fct)
138 : {
139 0 : set_aperture(make_sp(new LuaUserData<number,dim>(fct)));
140 0 : }
141 : #endif
142 :
143 0 : void set_ortho_velocity(SmartPtr<CplUserData<number, dim> > user)
144 : {
145 0 : m_imOrthoVelocity.set_data(user);
146 0 : }
147 0 : void set_ortho_velocity(number val)
148 : {
149 0 : set_ortho_velocity(make_sp(new ConstUserNumber<dim>(val)));
150 0 : }
151 : #ifdef UG_FOR_LUA
152 0 : void set_ortho_velocity(const char* fctName)
153 : {
154 0 : set_ortho_velocity(LuaUserDataFactory<number,dim>::create(fctName));
155 0 : }
156 0 : void set_ortho_velocity(LuaFunctionHandle fct)
157 : {
158 0 : set_ortho_velocity(make_sp(new LuaUserData<number,dim>(fct)));
159 0 : }
160 : #endif
161 :
162 0 : void set_ortho_diffusion(SmartPtr<CplUserData<number, dim> > user)
163 : {
164 0 : m_imOrthoDiffusion.set_data(user);
165 0 : }
166 0 : void set_ortho_diffusion(number val)
167 : {
168 0 : set_ortho_diffusion(make_sp(new ConstUserNumber<dim>(val)));
169 0 : }
170 : #ifdef UG_FOR_LUA
171 0 : void set_ortho_diffusion(const char* fctName)
172 : {
173 0 : set_ortho_diffusion(LuaUserDataFactory<number,dim>::create(fctName));
174 0 : }
175 0 : void set_ortho_diffusion(LuaFunctionHandle fct)
176 : {
177 0 : set_ortho_diffusion(make_sp(new LuaUserData<number,dim>(fct)));
178 0 : }
179 : #endif
180 :
181 : /// set singular sources and sinks
182 0 : void set_sss_manager(SmartPtr<CDSingularSourcesAndSinks<dim> > sss_mngr) {m_sss_mngr = sss_mngr;}
183 :
184 : /// get singular sources and sinks
185 0 : SmartPtr<CDSingularSourcesAndSinks<dim> > sss_manager() {return m_sss_mngr;}
186 :
187 : /// hanging nodes are not allowed in this discretization
188 0 : virtual bool use_hanging() const {return false;}
189 :
190 : private:
191 :
192 : /// registers the interface function
193 : void register_all_funcs ();
194 :
195 : /// auxiliary class for registering functions
196 : struct RegisterLocalDiscr
197 : {
198 : this_type * m_pThis;
199 :
200 0 : RegisterLocalDiscr(this_type * pThis) : m_pThis(pThis) {}
201 :
202 0 : template<typename TElem> void operator () (TElem &) {m_pThis->register_func<TElem> ();}
203 : };
204 :
205 : /// registers the interface function for one element type
206 : template <typename TElem> void register_func ();
207 :
208 : private:
209 :
210 : //---- General interface assembling functions: ----
211 :
212 : /// check the grid and the type of trial space for each function used
213 : virtual void prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid);
214 :
215 : /// prepares the loop over all elements
216 : /**
217 : * This method prepares the loop over all elements. It resizes the Position
218 : * array for the corner coordinates and schedules the local ip positions
219 : * at the data imports.
220 : */
221 : template <typename TElem>
222 : void prep_elem_loop(const ReferenceObjectID roid, const int si);
223 :
224 : /// finishes the loop over all elements
225 : template <typename TElem>
226 : void fsh_elem_loop ();
227 :
228 : /// prepares the element for assembling
229 : template <typename TElem>
230 : void prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const position_type vCornerCoords[]);
231 :
232 : /// assembles the local stiffness matrix using a finite volume scheme
233 : template <typename TElem>
234 : void add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
235 :
236 : /// assembles the local mass matrix using a finite volume scheme
237 : template <typename TElem>
238 : void add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
239 :
240 : /// assembles the stiffness part of the local defect
241 : template <typename TElem>
242 : void add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
243 :
244 : /// assembles the stiffness part of the local defect explicit reaction, reaction_rate and source
245 : template <typename TElem>
246 : void add_def_A_expl_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
247 :
248 : /// assembles the mass part of the local defect
249 : template <typename TElem>
250 : void add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const position_type vCornerCoords[]);
251 :
252 : /// assembles the local right hand side
253 : template <typename TElem>
254 : void add_rhs_elem(LocalVector& d, GridObject* elem, const position_type vCornerCoords[]);
255 :
256 : /// computes the linearized defect w.r.t. to the fracture velocity
257 : template <typename TElem>
258 : void lin_def_velocity(const LocalVector& u,
259 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
260 : const size_t nip);
261 :
262 : /// computes the linearized defect w.r.t. to the orthogonal velocity
263 : template <typename TElem>
264 : void lin_def_ortho_velocity(const LocalVector& u,
265 : std::vector<std::vector<number> > vvvLinDef[],
266 : const size_t nip);
267 :
268 : /// computes the linearized defect w.r.t. to the fracture diffusion
269 : template <typename TElem>
270 : void lin_def_diffusion(const LocalVector& u,
271 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
272 : const size_t nip);
273 :
274 : /// computes the linearized defect w.r.t. to the orthogonal diffusion
275 : template <typename TElem>
276 : void lin_def_ortho_diffusion(const LocalVector& u,
277 : std::vector<std::vector<number> > vvvLinDef[],
278 : const size_t nip);
279 :
280 : /// computes the linearized defect w.r.t. to the fracture flux
281 : template <typename TElem>
282 : void lin_def_flux(const LocalVector& u,
283 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
284 : const size_t nip);
285 :
286 : /// computes the linearized defect w.r.t. to the orthogonal flux
287 : template <typename TElem>
288 : void lin_def_ortho_flux(const LocalVector& u,
289 : std::vector<std::vector<number> > vvvLinDef[],
290 : const size_t nip);
291 :
292 : /// computes the linearized defect w.r.t. to the reaction
293 : template <typename TElem>
294 : void lin_def_reaction(const LocalVector& u,
295 : std::vector<std::vector<number> > vvvLinDef[],
296 : const size_t nip);
297 :
298 : /// computes the linearized defect w.r.t. to the reaction
299 : template <typename TElem>
300 : void lin_def_reaction_rate(const LocalVector& u,
301 : std::vector<std::vector<number> > vvvLinDef[],
302 : const size_t nip);
303 :
304 : /// computes the linearized defect w.r.t to the source term
305 : template <typename TElem>
306 : void lin_def_source(const LocalVector& u,
307 : std::vector<std::vector<number> > vvvLinDef[],
308 : const size_t nip);
309 :
310 : /// computes the linearized defect w.r.t to the vector source term
311 : template <typename TElem>
312 : void lin_def_vector_source(const LocalVector& u,
313 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
314 : const size_t nip);
315 :
316 : /// computes the linearized defect w.r.t to the vector source term
317 : template <typename TElem>
318 : void lin_def_ortho_vector_source(const LocalVector& u,
319 : std::vector<std::vector<number> > vvvLinDef[],
320 : const size_t nip);
321 :
322 : /// computes the linearized defect w.r.t to the mass scale term
323 : template <typename TElem>
324 : void lin_def_mass_scale(const LocalVector& u,
325 : std::vector<std::vector<number> > vvvLinDef[],
326 : const size_t nip);
327 :
328 : /// computes the linearized defect w.r.t to the mass scale term
329 : template <typename TElem>
330 : void lin_def_mass(const LocalVector& u,
331 : std::vector<std::vector<number> > vvvLinDef[],
332 : const size_t nip);
333 :
334 : //---- Assembling functions for fractures: ----
335 :
336 : template <typename TElem>
337 : inline void fract_add_jac_A_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
338 :
339 : template <typename TElem>
340 : inline void fract_bulk_add_jac_A_elem(LocalMatrix& J, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
341 :
342 : template <typename TElem>
343 : inline void fract_add_def_A_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
344 :
345 : template <typename TElem>
346 : inline void fract_bulk_add_def_A_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
347 :
348 : template <typename TElem>
349 : inline void fract_add_def_A_expl_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
350 :
351 : template <typename TElem>
352 : inline void fract_bulk_add_def_A_expl_elem(LocalVector& d, const LocalVector& u, TElem* elem, const position_type vCornerCoords[]);
353 :
354 : template <typename TElem>
355 : inline void fract_add_rhs_elem(LocalVector& b, TElem* elem, const position_type vCornerCoords[]);
356 :
357 : template <typename TElem>
358 : inline void fract_bulk_add_rhs_elem(LocalVector& b, TElem* elem, const position_type vCornerCoords[]);
359 :
360 : template <typename TElem>
361 : inline void add_sss_def_elem(LocalVector& d, const LocalVector& u, TElem* pElem, size_t co, number flux);
362 :
363 : template <typename TElem>
364 : inline void add_sss_jac_elem(LocalMatrix& J, const LocalVector& u, TElem* pElem, size_t co, number flux);
365 :
366 : private:
367 : /// returns the updated convection shapes
368 : const conv_shape_type& get_updated_conv_shapes (bool computeDeriv);
369 :
370 : protected:
371 :
372 : // Standard parameters of the convection-diffusion plugin
373 : using base_type::m_imDiffusion;
374 : using base_type::m_imVelocity;
375 : using base_type::m_imFlux;
376 : using base_type::m_imSource;
377 : using base_type::m_imSourceExpl;
378 : using base_type::m_imVectorSource;
379 : using base_type::m_imReactionRate;
380 : using base_type::m_imReactionRateExpl;
381 : using base_type::m_imReaction;
382 : using base_type::m_imReactionExpl;
383 : using base_type::m_imMassScale;
384 : using base_type::m_imMass;
385 :
386 : using base_type::m_exGrad;
387 : using base_type::m_exValue;
388 :
389 : // Specific parameters for low-dimensional subdomains (fractures):
390 : DataImport<number, dim> m_imAperture; ///< the fracture width (constant per element)
391 : DataImport<number, dim> m_imOrthoDiffusion; ///< diffusion through the fracture-bulk interface
392 : DataImport<number, dim> m_imOrthoVelocity; ///< convection velocity through the fracture-bulk interface
393 : DataImport<number, dim> m_imOrthoFlux; ///< flux through the fracture-bulk interface
394 : DataImport<number, dim> m_imOrthoVectorSource; ///< vector through the fracture-bulk interface
395 :
396 : private:
397 :
398 : // singular sources and sinks manager
399 : SmartPtr<CDSingularSourcesAndSinks<dim> > m_sss_mngr;
400 :
401 : protected:
402 :
403 : // Parameters of the discretization:
404 : SmartPtr<fract_manager_type> m_spFractManager; ///< degenerated fracture manager (may be SPNULL)
405 : SmartPtr<conv_shape_type> m_spConvShape; ///< method to compute the upwind shapes
406 :
407 : private:
408 :
409 : // Temporary data used in the assembling
410 : TFractFVGeom * m_pFractGeo; ///< FV geometry object of fracture elements
411 : size_t m_numFractCo; ///< number of corners of the fracture side
412 : side_type * m_innerFractSide; ///< inner side of the fracture element
413 : side_type * m_outerFractSide; ///< outer side of the fracture element
414 : size_t m_innerFractSideIdx; ///< index of the inner side in the ref. elem.
415 : size_t m_outerFractSideIdx; ///< index of the outer side in the ref. elem.
416 : size_t m_innerSideCo[maxFractSideCorners]; ///< inner side corner idx -> elem. corner idx
417 : size_t m_outerSideCo[maxFractSideCorners]; ///< outer side corner idx -> elem. corner idx
418 : size_t m_assCo [2 * maxFractSideCorners]; ///< correspondence of the corners of the sides
419 : };
420 :
421 : } // end ConvectionDiffusionPlugin
422 : } // end namespace ug
423 :
424 : #endif /*__H__UG__LIB_DISC__CONVECTION_DIFFUSION__CONVECTION_DIFFUSION_FRACT_FV1__*/
425 :
426 : /* End of File */
|
