Line data Source code
1 : /*
2 : * Copyright (c) 2018-: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Naegel
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 : #include "convection_diffusion_stab_fe.h"
34 :
35 :
36 : #include "lib_disc/spatial_disc/disc_util/fe_geom.h"
37 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
38 : #include "lib_disc/local_finite_element/lagrange/lagrange.h"
39 : #include "lib_disc/local_finite_element/lagrange/lagrangep1.h"
40 : #include "lib_disc/quadrature/gauss/gauss_quad.h"
41 :
42 : #include "lib_grid/algorithms/volume_calculation.h" // For computing bounding box
43 :
44 : namespace ug{
45 : namespace ConvectionDiffusionPlugin{
46 :
47 : ////////////////////////////////////////////////////////////////////////////////
48 : // general
49 : ////////////////////////////////////////////////////////////////////////////////
50 :
51 : template<typename TDomain>
52 0 : ConvectionDiffusionStabFE<TDomain>::
53 : ConvectionDiffusionStabFE(const char* functions, const char* subsets)
54 : : base_type(functions,subsets),
55 0 : m_bQuadOrderUserDef(false), m_stabParamM(0.0), m_stabParamA(0.0)
56 : {
57 : this->clear_add_fct();
58 0 : }
59 :
60 :
61 : template<typename TDomain>
62 0 : ConvectionDiffusionStabFE<TDomain>::
63 : ConvectionDiffusionStabFE(const char* functions, const char* subsets, double stabM)
64 : : base_type(functions,subsets),
65 0 : m_bQuadOrderUserDef(false), m_stabParamM(stabM), m_stabParamA(0.0)
66 : {
67 : this->clear_add_fct();
68 0 : }
69 :
70 : template<typename TDomain>
71 0 : ConvectionDiffusionStabFE<TDomain>::
72 : ConvectionDiffusionStabFE(const char* functions, const char* subsets, double stabM, double stabA)
73 : : base_type(functions,subsets),
74 0 : m_bQuadOrderUserDef(false), m_stabParamM(stabM), m_stabParamA(stabA)
75 : {
76 : this->clear_add_fct();
77 0 : }
78 :
79 :
80 : template<typename TDomain>
81 0 : void ConvectionDiffusionStabFE<TDomain>::
82 : set_quad_order(size_t order)
83 : {
84 0 : m_quadOrder = order;
85 0 : m_bQuadOrderUserDef = true;
86 0 : }
87 :
88 : template<typename TDomain>
89 0 : void ConvectionDiffusionStabFE<TDomain>::
90 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
91 : {
92 : // check number of fcts
93 0 : if(vLfeID.size() != 1)
94 0 : UG_THROW("ConvectionDiffusionStab: Wrong number of functions given. "
95 : "Need exactly "<<1);
96 :
97 : // check that not ADAPTIVE
98 0 : if(vLfeID[0].order() < 1)
99 0 : UG_THROW("ConvectionDiffusionStab: Adaptive order not implemented.");
100 :
101 : // set order
102 0 : m_lfeID = vLfeID[0];
103 0 : if(!m_bQuadOrderUserDef) m_quadOrder = 2*m_lfeID.order()+1;
104 :
105 0 : register_all_funcs(m_lfeID, m_quadOrder);
106 0 : }
107 :
108 : template<typename TDomain>
109 0 : bool ConvectionDiffusionStabFE<TDomain>::
110 : use_hanging() const
111 : {
112 0 : return false;
113 : }
114 :
115 : ////////////////////////////////////////////////////////////////////////////////
116 : // Assembling functions
117 : ////////////////////////////////////////////////////////////////////////////////
118 :
119 : template<typename TDomain>
120 : template<typename TElem, typename TFEGeom>
121 0 : void ConvectionDiffusionStabFE<TDomain>::
122 : prep_elem_loop(const ReferenceObjectID roid, const int si)
123 : {
124 :
125 : // request geometry
126 0 : TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
127 :
128 : // prepare geometry for type and order
129 : try{
130 0 : geo.update_local(roid, m_lfeID, m_quadOrder);
131 0 : }UG_CATCH_THROW("ConvectionDiffusion::prep_elem_loop:"
132 : " Cannot update Finite Element Geometry.");
133 :
134 : // set local positions
135 : // static const int refDim = TElem::dim;
136 : // m_imDiffusion.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
137 :
138 0 : }
139 :
140 : template<typename TDomain>
141 : template<typename TElem, typename TFEGeom>
142 0 : void ConvectionDiffusionStabFE<TDomain>::
143 : fsh_elem_loop()
144 0 : {}
145 :
146 : template<typename TDomain>
147 : template<typename TElem, typename TFEGeom>
148 0 : void ConvectionDiffusionStabFE<TDomain>::
149 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
150 : {
151 : // request geometry
152 0 : TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
153 :
154 : try{
155 : geo.update(elem, vCornerCoords);
156 : }
157 0 : UG_CATCH_THROW("ConvectionDiffusion::prep_elem:"
158 : " Cannot update Finite Element Geometry.");
159 :
160 : // set global positions for rhs
161 : // m_imDiffusion. set_global_ips(geo.global_ips(), geo.num_ip());
162 :
163 0 : }
164 :
165 :
166 : /*!
167 : * Assemble contribution: \sigma \int_{tau} \sum h_i^2 \partial_i u^k \partial_i v
168 : */
169 : template<typename TDomain>
170 : template<typename TElem, typename TFEGeom>
171 0 : void ConvectionDiffusionStabFE<TDomain>::
172 : add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
173 : {
174 : const TFEGeom& geo = // Geometry
175 0 : GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
176 : /* const number scale = // Scaling factor
177 : ElementDiameterSq<GridObject, TDomain>(*elem, *this->domain()); */
178 :
179 0 : MathVector<dim> vBoundingBox[2]; // Element bounding box
180 0 : CalculateBoundingBox<dim>(NumVertices((TElem*)elem), vCornerCoords, vBoundingBox[0], vBoundingBox[1]);
181 :
182 : MathVector<dim> gradU; // Temporary
183 :
184 : // Loop integration points.
185 0 : for (size_t ip = 0; ip < geo.num_ip(); ++ip)
186 : {
187 : // Loop trial space
188 : VecSet(gradU, 0.0);
189 0 : for (size_t psh = 0; psh < geo.num_sh(); ++psh)
190 : VecScaleAppend(gradU, u(_C_, psh), geo.global_grad(ip, psh)); // * scale
191 :
192 : // Scale gradient
193 0 : for(size_t i = 0; i < dim; ++i)
194 : {
195 0 : double hi = vBoundingBox[1][i] - vBoundingBox[0][i];
196 0 : gradU[i] *= (hi*hi);
197 : }
198 :
199 0 : for (size_t psh = 0; psh < geo.num_sh(); ++psh)
200 : {
201 0 : d(_C_, psh) +=
202 0 : m_stabParamM * VecDot(geo.global_grad(ip, psh), gradU) * geo.weight(ip);
203 : }
204 : }
205 :
206 0 : }
207 :
208 :
209 :
210 : template<typename TDomain>
211 : template<typename TElem, typename TFEGeom>
212 0 : void ConvectionDiffusionStabFE<TDomain>::
213 : add_jac_X_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], double m_stabParam)
214 : {
215 : const TFEGeom& geo = // request geometry
216 0 : GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
217 : /* const number scale = // Scaling factor
218 : ElementDiameterSq<GridObject, TDomain>(*elem, *this->domain());*/
219 :
220 0 : MathVector<dim> vBoundingBox[2];
221 0 : CalculateBoundingBox<dim>(NumVertices((TElem*)elem), vCornerCoords, vBoundingBox[0], vBoundingBox[1]);
222 :
223 : // rescale using bounding box
224 0 : for (size_t ip = 0; ip < geo.num_ip(); ++ip){
225 0 : for (size_t psh = 0; psh < geo.num_sh(); ++psh){
226 :
227 : // Rescaled gradient (using bounding box)
228 : MathVector<dim> grad = geo.global_grad(ip, psh); // *scale
229 0 : for(size_t i = 0; i < dim; ++i)
230 : {
231 0 : double hi = vBoundingBox[1][i] - vBoundingBox[0][i];
232 0 : grad[i] *= (hi*hi);
233 : }
234 :
235 0 : for (size_t psh2 = 0; psh2 < geo.num_sh(); ++psh2)
236 : {
237 0 : J(_C_, psh, _C_, psh2) += //* scale
238 0 : m_stabParam* VecDot(grad, geo.global_grad(ip, psh2)) * geo.weight(ip);
239 : }
240 : }
241 : }
242 :
243 0 : }
244 :
245 : template<typename TDomain>
246 : template<typename TElem, typename TFEGeom>
247 0 : void ConvectionDiffusionStabFE<TDomain>::
248 : add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
249 : {
250 0 : if (m_stabParamM != 0.0) add_jac_X_elem<TElem,TFEGeom>(J, u, elem,vCornerCoords, m_stabParamM);
251 0 : }
252 :
253 : template<typename TDomain>
254 : template<typename TElem, typename TFEGeom>
255 0 : void ConvectionDiffusionStabFE<TDomain>::
256 : add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
257 : {
258 0 : if (m_stabParamA != 0.0) add_jac_X_elem<TElem,TFEGeom>(J, u, elem,vCornerCoords, m_stabParamA);
259 0 : }
260 :
261 :
262 :
263 :
264 :
265 : ////////////////////////////////////////////////////////////////////////////////
266 : // register assemble functions
267 : ////////////////////////////////////////////////////////////////////////////////
268 :
269 : #ifdef UG_DIM_1
270 : template<>
271 0 : void ConvectionDiffusionStabFE<Domain1d>::
272 : register_all_funcs(const LFEID& lfeid, const int quadOrder)
273 : {
274 : // RegularEdge
275 0 : register_func<RegularEdge, DimFEGeometry<dim> >();
276 0 : }
277 : #endif
278 :
279 : #ifdef UG_DIM_2
280 : template<>
281 0 : void ConvectionDiffusionStabFE<Domain2d>::
282 : register_all_funcs(const LFEID& lfeid, const int quadOrder)
283 : {
284 0 : register_func<RegularEdge, DimFEGeometry<dim, 1> >();
285 :
286 : const int order = lfeid.order();
287 0 : if(quadOrder != 2*order+1 || lfeid.type() != LFEID::LAGRANGE)
288 : {
289 0 : register_func<Triangle, DimFEGeometry<dim> >();
290 0 : register_func<Quadrilateral, DimFEGeometry<dim> >();
291 0 : return;
292 : }
293 :
294 : // special compiled cases
295 :
296 : // Triangle
297 0 : switch(order)
298 : {
299 0 : case 1: {typedef FEGeometry<Triangle, dim, LagrangeLSFS<ReferenceTriangle, 1>, GaussQuadrature<ReferenceTriangle, 3> > FEGeom;
300 0 : register_func<Triangle, FEGeom >(); break;}
301 0 : case 2: {typedef FEGeometry<Triangle, dim, LagrangeLSFS<ReferenceTriangle, 2>, GaussQuadrature<ReferenceTriangle, 5> > FEGeom;
302 0 : register_func<Triangle, FEGeom >(); break;}
303 0 : case 3: {typedef FEGeometry<Triangle, dim, LagrangeLSFS<ReferenceTriangle, 3>, GaussQuadrature<ReferenceTriangle, 7> > FEGeom;
304 0 : register_func<Triangle, FEGeom >(); break;}
305 0 : default: register_func<Triangle, DimFEGeometry<dim> >(); break;
306 : }
307 :
308 : // Quadrilateral
309 0 : switch(order) {
310 0 : case 1: {typedef FEGeometry<Quadrilateral, dim, LagrangeLSFS<ReferenceQuadrilateral, 1>, GaussQuadrature<ReferenceQuadrilateral, 3> > FEGeom;
311 0 : register_func<Quadrilateral, FEGeom >(); break;}
312 0 : case 2: {typedef FEGeometry<Quadrilateral, dim, LagrangeLSFS<ReferenceQuadrilateral, 2>, GaussQuadrature<ReferenceQuadrilateral, 7> > FEGeom;
313 0 : register_func<Quadrilateral, FEGeom >(); break;}
314 0 : case 3: {typedef FEGeometry<Quadrilateral, dim, LagrangeLSFS<ReferenceQuadrilateral, 3>, GaussQuadrature<ReferenceQuadrilateral, 11> > FEGeom;
315 0 : register_func<Quadrilateral, FEGeom >(); break;}
316 0 : default: register_func<Quadrilateral, DimFEGeometry<dim> >(); break;
317 : }
318 : }
319 : #endif
320 :
321 : #ifdef UG_DIM_3
322 : template<>
323 0 : void ConvectionDiffusionStabFE<Domain3d>::
324 : register_all_funcs(const LFEID& lfeid, const int quadOrder)
325 : {
326 0 : register_func<RegularEdge, DimFEGeometry<dim, 1> >();
327 0 : register_func<Triangle, DimFEGeometry<dim, 2> >();
328 0 : register_func<Quadrilateral, DimFEGeometry<dim, 2> >();
329 :
330 : const int order = lfeid.order();
331 0 : if(quadOrder != 2*order+1 || lfeid.type() != LFEID::LAGRANGE)
332 : {
333 0 : register_func<Tetrahedron, DimFEGeometry<dim> >();
334 0 : register_func<Prism, DimFEGeometry<dim> >();
335 0 : register_func<Pyramid, DimFEGeometry<dim> >();
336 0 : register_func<Hexahedron, DimFEGeometry<dim> >();
337 0 : register_func<Octahedron, DimFEGeometry<dim> >();
338 0 : return;
339 : }
340 :
341 : // special compiled cases
342 :
343 : // Tetrahedron
344 0 : switch(order)
345 : {
346 0 : case 1: {typedef FEGeometry<Tetrahedron, dim, LagrangeLSFS<ReferenceTetrahedron, 1>, GaussQuadrature<ReferenceTetrahedron, 3> > FEGeom;
347 0 : register_func<Tetrahedron, FEGeom >(); break;}
348 0 : case 2: {typedef FEGeometry<Tetrahedron, dim, LagrangeLSFS<ReferenceTetrahedron, 2>, GaussQuadrature<ReferenceTetrahedron, 5> > FEGeom;
349 0 : register_func<Tetrahedron, FEGeom >(); break;}
350 0 : case 3: {typedef FEGeometry<Tetrahedron, dim, LagrangeLSFS<ReferenceTetrahedron, 3>, GaussQuadrature<ReferenceTetrahedron, 7> > FEGeom;
351 0 : register_func<Tetrahedron, FEGeom >(); break;}
352 0 : default: register_func<Tetrahedron, DimFEGeometry<dim> >(); break;
353 : }
354 :
355 : // Prism
356 : switch(order) {
357 0 : default: register_func<Prism, DimFEGeometry<dim> >(); break;
358 : }
359 :
360 : // Pyramid
361 : switch(order)
362 : {
363 0 : default: register_func<Pyramid, DimFEGeometry<dim> >(); break;
364 : }
365 :
366 : // Octahedron
367 : switch(order)
368 : {
369 0 : default: register_func<Octahedron, DimFEGeometry<dim> >(); break;
370 : }
371 :
372 : // Hexahedron
373 0 : switch(order)
374 : {
375 0 : case 1: {typedef FEGeometry<Hexahedron, dim, LagrangeLSFS<ReferenceHexahedron, 1>, GaussQuadrature<ReferenceHexahedron, 3> > FEGeom;
376 0 : register_func<Hexahedron, FEGeom >(); break;}
377 0 : case 2: {typedef FEGeometry<Hexahedron, dim, LagrangeLSFS<ReferenceHexahedron, 2>, GaussQuadrature<ReferenceHexahedron, 7> > FEGeom;
378 0 : register_func<Hexahedron, FEGeom >(); break;}
379 0 : case 3: {typedef FEGeometry<Hexahedron, dim, LagrangeLSFS<ReferenceHexahedron, 3>, GaussQuadrature<ReferenceHexahedron, 11> > FEGeom;
380 0 : register_func<Hexahedron, FEGeom >(); break;}
381 0 : default: register_func<Hexahedron, DimFEGeometry<dim> >(); break;
382 : }
383 : }
384 : #endif
385 :
386 : template <typename TDomain>
387 : template <typename TElem, typename TFEGeom>
388 0 : void ConvectionDiffusionStabFE<TDomain>::register_func()
389 : {
390 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
391 : typedef this_type T;
392 : // static const int refDim = reference_element_traits<TElem>::dim;
393 :
394 : this->clear_add_fct(id);
395 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFEGeom>);
396 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFEGeom>);
397 : this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop<TElem, TFEGeom>);
398 : this->set_add_jac_A_elem_fct(id, &T::template add_jac_A_elem<TElem, TFEGeom>);
399 : this->set_add_jac_M_elem_fct(id, &T::template add_jac_M_elem<TElem, TFEGeom>);
400 : this->set_add_def_A_elem_fct(id, &T::template add_def_A_elem<TElem, TFEGeom>);
401 : this->set_add_def_M_elem_fct(id, &T::template add_def_M_elem<TElem, TFEGeom>);
402 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFEGeom>);
403 :
404 : // error estimator parts
405 : this->set_prep_err_est_elem_loop(id, &T::template prep_err_est_elem_loop<TElem, TFEGeom>);
406 : this->set_prep_err_est_elem(id, &T::template prep_err_est_elem<TElem, TFEGeom>);
407 : this->set_compute_err_est_A_elem(id, &T::template compute_err_est_A_elem<TElem, TFEGeom>);
408 : this->set_compute_err_est_M_elem(id, &T::template compute_err_est_M_elem<TElem, TFEGeom>);
409 : this->set_compute_err_est_rhs_elem(id, &T::template compute_err_est_rhs_elem<TElem, TFEGeom>);
410 : this->set_fsh_err_est_elem_loop(id, &T::template fsh_err_est_elem_loop<TElem, TFEGeom>);
411 :
412 : // set computation of linearized defect w.r.t velocity
413 : /* m_imDiffusion. set_fct(id, this, &T::template lin_def_diffusion<TElem, TFEGeom>);
414 : m_imVelocity. set_fct(id, this, &T::template lin_def_velocity<TElem, TFEGeom>);
415 : m_imFlux. set_fct(id, this, &T::template lin_def_flux<TElem, TFEGeom>);
416 : m_imReactionRate. set_fct(id, this, &T::template lin_def_reaction_rate<TElem, TFEGeom>);
417 : m_imReaction. set_fct(id, this, &T::template lin_def_reaction<TElem, TFEGeom>);
418 : m_imSource. set_fct(id, this, &T::template lin_def_source<TElem, TFEGeom>);
419 : m_imVectorSource. set_fct(id, this, &T::template lin_def_vector_source<TElem, TFEGeom>);
420 : m_imMassScale. set_fct(id, this, &T::template lin_def_mass_scale<TElem, TFEGeom>);
421 : m_imMass. set_fct(id, this, &T::template lin_def_mass<TElem, TFEGeom>);
422 :
423 : // exports
424 : m_exValue-> template set_fct<T,refDim>(id, this, &T::template ex_value<TElem, TFEGeom>);
425 : m_exGrad-> template set_fct<T,refDim>(id, this, &T::template ex_grad<TElem, TFEGeom>);*/
426 0 : }
427 :
428 : ////////////////////////////////////////////////////////////////////////////////
429 : // explicit template instantiations
430 : ////////////////////////////////////////////////////////////////////////////////
431 :
432 : #ifdef UG_DIM_1
433 : template class ConvectionDiffusionStabFE<Domain1d>;
434 : #endif
435 : #ifdef UG_DIM_2
436 : template class ConvectionDiffusionStabFE<Domain2d>;
437 : #endif
438 : #ifdef UG_DIM_3
439 : template class ConvectionDiffusionStabFE<Domain3d>;
440 : #endif
441 :
442 : } // end namespace ConvectionDiffusionPlugin
443 : } // namespace ug
444 :
|
