Line data Source code
1 : /*
2 : * Copyright (c) 2010-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 : #include "convection_diffusion_fe.h"
34 :
35 : #include "lib_disc/spatial_disc/disc_util/fe_geom.h"
36 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
37 : #include "lib_disc/local_finite_element/lagrange/lagrange.h"
38 : #include "lib_disc/local_finite_element/lagrange/lagrangep1.h"
39 : #include "lib_disc/quadrature/gauss/gauss_quad.h"
40 :
41 : namespace ug{
42 : namespace ConvectionDiffusionPlugin{
43 :
44 : ////////////////////////////////////////////////////////////////////////////////
45 : // general
46 : ////////////////////////////////////////////////////////////////////////////////
47 :
48 : template<typename TDomain>
49 0 : ConvectionDiffusionFE<TDomain>::
50 : ConvectionDiffusionFE(const char* functions, const char* subsets)
51 : : ConvectionDiffusionBase<TDomain>(functions,subsets),
52 0 : m_bQuadOrderUserDef(false)
53 : {
54 : this->clear_add_fct();
55 0 : }
56 :
57 : template<typename TDomain>
58 0 : void ConvectionDiffusionFE<TDomain>::set_quad_order(size_t order)
59 : {
60 0 : m_quadOrder = order;
61 0 : m_bQuadOrderUserDef = true;
62 0 : }
63 :
64 : template<typename TDomain>
65 0 : void ConvectionDiffusionFE<TDomain>::
66 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
67 : {
68 : // check number of fcts
69 0 : UG_COND_THROW(vLfeID.size() != 1, "ConvectionDiffusion:" <<
70 : "Wrong number of functions given. Need exactly "<<1);
71 :
72 : // check that not ADAPTIVE
73 0 : UG_COND_THROW(vLfeID[0].order() < 0, "ConvectionDiffusion: "<<
74 : "Adaptive order not implemented.");
75 :
76 : // set order
77 0 : m_lfeID = vLfeID[0];
78 0 : if(!m_bQuadOrderUserDef) m_quadOrder = 2*m_lfeID.order()+1;
79 :
80 0 : register_all_funcs(m_lfeID, m_quadOrder);
81 0 : }
82 :
83 : template<typename TDomain>
84 0 : bool ConvectionDiffusionFE<TDomain>::
85 : use_hanging() const
86 0 : { return false; }
87 :
88 : ////////////////////////////////////////////////////////////////////////////////
89 : // Assembling functions
90 : ////////////////////////////////////////////////////////////////////////////////
91 :
92 : template<typename TDomain>
93 : template<typename TElem, typename TFEGeom>
94 0 : void ConvectionDiffusionFE<TDomain>::
95 : prep_elem_loop(const ReferenceObjectID roid, const int si)
96 : {
97 : // Check for explicit terms.
98 0 : UG_COND_THROW(m_imSourceExpl.data_given() || m_imReactionExpl.data_given() || m_imReactionRateExpl.data_given(),
99 : "ConvectionDiffusionFE: Explicit terms not implemented.");
100 :
101 : // Request geometry.
102 0 : TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
103 :
104 : // Prepare geometry for type and order.
105 : try{
106 0 : geo.update_local(roid, m_lfeID, m_quadOrder);
107 0 : } UG_CATCH_THROW("ConvectionDiffusion::prep_elem_loop:" <<
108 : " Cannot update Finite Element Geometry.");
109 :
110 : // Set local positions.
111 : static const int refDim = TElem::dim;
112 0 : m_imDiffusion.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
113 0 : m_imVelocity.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
114 0 : m_imFlux.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
115 0 : m_imSource.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
116 0 : m_imVectorSource.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
117 0 : m_imReactionRate.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
118 0 : m_imReaction.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
119 0 : m_imMassScale.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
120 0 : m_imMass.template set_local_ips<refDim>(geo.local_ips(), geo.num_ip(), false);
121 0 : }
122 :
123 : template<typename TDomain>
124 : template<typename TElem, typename TFEGeom>
125 0 : void ConvectionDiffusionFE<TDomain>::
126 : fsh_elem_loop()
127 0 : {}
128 :
129 : template<typename TDomain>
130 : template<typename TElem, typename TFEGeom>
131 0 : void ConvectionDiffusionFE<TDomain>::
132 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
133 : {
134 : // request geometry
135 0 : TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
136 :
137 : try{
138 : geo.update(elem, vCornerCoords);
139 : }
140 0 : UG_CATCH_THROW("ConvectionDiffusion::prep_elem:"
141 : " Cannot update Finite Element Geometry.");
142 :
143 : // set global positions for rhs
144 0 : m_imDiffusion. set_global_ips(geo.global_ips(), geo.num_ip());
145 0 : m_imVelocity. set_global_ips(geo.global_ips(), geo.num_ip());
146 0 : m_imFlux. set_global_ips(geo.global_ips(), geo.num_ip());
147 0 : m_imSource. set_global_ips(geo.global_ips(), geo.num_ip());
148 0 : m_imVectorSource.set_global_ips(geo.global_ips(), geo.num_ip());
149 0 : m_imReactionRate.set_global_ips(geo.global_ips(), geo.num_ip());
150 0 : m_imReaction. set_global_ips(geo.global_ips(), geo.num_ip());
151 0 : m_imMassScale. set_global_ips(geo.global_ips(), geo.num_ip());
152 0 : m_imMass. set_global_ips(geo.global_ips(), geo.num_ip());
153 0 : }
154 :
155 : template<typename TDomain>
156 : template<typename TElem, typename TFEGeom>
157 0 : void ConvectionDiffusionFE<TDomain>::
158 : add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
159 : {
160 : // request geometry
161 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
162 :
163 : MathVector<dim> v, Dgrad;
164 :
165 : // loop integration points
166 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
167 : {
168 : // loop trial space
169 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
170 : {
171 : // Diffusion
172 0 : if(m_imDiffusion.data_given())
173 : MatVecMult(Dgrad, m_imDiffusion[ip], geo.global_grad(ip, j));
174 : else
175 : VecSet(Dgrad, 0.0);
176 :
177 : // Convection
178 0 : if(m_imVelocity.data_given())
179 0 : VecScaleAppend(Dgrad, -1*geo.shape(ip,j), m_imVelocity[ip]);
180 :
181 : // no explicit dependency on m_imFlux
182 :
183 : // loop test space
184 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
185 : {
186 : // compute integrand
187 : number integrand = VecDot(Dgrad, geo.global_grad(ip, i));
188 :
189 : // Reaction
190 0 : if(m_imReactionRate.data_given())
191 0 : integrand += m_imReactionRate[ip] * geo.shape(ip, j) * geo.shape(ip, i);
192 :
193 : // no explicit dependency on m_imReaction
194 :
195 : // multiply by weight
196 0 : integrand *= geo.weight(ip);
197 :
198 : // add to local matrix
199 0 : J(_C_, i, _C_, j) += integrand;
200 : }
201 : }
202 : }
203 0 : }
204 :
205 :
206 : template<typename TDomain>
207 : template<typename TElem, typename TFEGeom>
208 0 : void ConvectionDiffusionFE<TDomain>::
209 : add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
210 : {
211 : // request geometry
212 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
213 :
214 0 : if(!m_imMassScale.data_given()) return;
215 :
216 : // loop integration points
217 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
218 : {
219 : // loop test space
220 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
221 : {
222 : // loop trial space
223 0 : for(size_t j= 0; j < geo.num_sh(); ++j)
224 : {
225 : // add to local matrix
226 0 : J(_C_, i, _C_, j) += geo.shape(ip, i) *geo.shape(ip, j)
227 0 : * geo.weight(ip) * m_imMassScale[ip];
228 :
229 : // no explicit dependency on m_imMass
230 : }
231 : }
232 : }
233 : }
234 :
235 :
236 : template<typename TDomain>
237 : template<typename TElem, typename TFEGeom>
238 0 : void ConvectionDiffusionFE<TDomain>::
239 : add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
240 : {
241 : // request geometry
242 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
243 :
244 : number integrand, shape_u;
245 : MathMatrix<dim,dim> D;
246 : MathVector<dim> v, Dgrad_u, grad_u;
247 :
248 : // loop integration points
249 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
250 : {
251 : // get current u and grad_u
252 : VecSet(grad_u, 0.0);
253 : shape_u = 0.0;
254 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
255 : {
256 : VecScaleAppend(grad_u, u(_C_,j), geo.global_grad(ip, j));
257 0 : shape_u += u(_C_,j) * geo.shape(ip, j);
258 : }
259 :
260 : // Diffusion
261 0 : if(m_imDiffusion.data_given())
262 : MatVecMult(Dgrad_u, m_imDiffusion[ip], grad_u);
263 : else
264 : VecSet(Dgrad_u, 0.0);
265 :
266 : // Convection
267 0 : if(m_imVelocity.data_given())
268 0 : VecScaleAppend(Dgrad_u, -1*shape_u, m_imVelocity[ip]);
269 :
270 : // Convection
271 0 : if(m_imFlux.data_given())
272 : VecScaleAppend(Dgrad_u, 1.0, m_imFlux[ip]);
273 :
274 : // loop test spaces
275 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
276 : {
277 : // compute integrand
278 : integrand = VecDot(Dgrad_u, geo.global_grad(ip, i));
279 :
280 : // add Reaction Rate
281 0 : if(m_imReactionRate.data_given())
282 0 : integrand += m_imReactionRate[ip] * shape_u * geo.shape(ip, i);
283 :
284 : // add Reaction
285 0 : if(m_imReaction.data_given())
286 0 : integrand += m_imReaction[ip] * geo.shape(ip, i);
287 :
288 : // multiply by integration weight
289 0 : integrand *= geo.weight(ip);
290 :
291 : // add to local defect
292 0 : d(_C_, i) += integrand;
293 : }
294 : }
295 0 : }
296 :
297 :
298 : template<typename TDomain>
299 : template<typename TElem, typename TFEGeom>
300 0 : void ConvectionDiffusionFE<TDomain>::
301 : add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
302 : {
303 : // request geometry
304 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
305 :
306 : number shape_u = 0.0;
307 :
308 0 : if(!m_imMassScale.data_given() && !m_imMass.data_given()) return;
309 :
310 : // loop integration points
311 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
312 : {
313 : // compute value of current solution at ip
314 0 : if(m_imMassScale.data_given()){
315 : shape_u = 0.0;
316 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
317 0 : shape_u += u(_C_,j) * geo.shape(ip, j);
318 : }
319 :
320 : // loop test spaces
321 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
322 : {
323 : // compute contribution
324 : number val = 0.0;
325 :
326 : // add MassScaling
327 0 : if(m_imMassScale.data_given())
328 0 : val += shape_u * m_imMassScale[ip];
329 :
330 : // add Maxx
331 0 : if(m_imMass.data_given())
332 0 : val += m_imMass[ip];
333 :
334 : // add to local defect
335 0 : d(_C_, i) += val * geo.shape(ip, i) * geo.weight(ip);
336 : }
337 : }
338 : };
339 :
340 : template<typename TDomain>
341 : template<typename TElem, typename TFEGeom>
342 0 : void ConvectionDiffusionFE<TDomain>::
343 : add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[])
344 : {
345 : // request geometry
346 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
347 :
348 : // skip if no source present
349 0 : if(!m_imSource.data_given() && !m_imVectorSource.data_given()) return;
350 :
351 : // loop integration points
352 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
353 : {
354 : // loop test spaces
355 :
356 : // only do this if (volume) source is given
357 0 : if(m_imSource.data_given()) {
358 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
359 : {
360 : // add contribution to local defect
361 0 : d(_C_, i) += m_imSource[ip] * geo.shape(ip, i) * geo.weight(ip);
362 : }
363 : }
364 :
365 : // only do this if vector source is given
366 0 : if(m_imVectorSource.data_given()) {
367 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
368 : {
369 : // add contribution to local defect
370 0 : d(_C_, i) += geo.weight(ip) * VecDot(m_imVectorSource[ip], geo.global_grad(ip, i));
371 : }
372 : }
373 : }
374 : }
375 :
376 : // ////////////////////////////////
377 : // error estimation (begin) ///
378 :
379 : // prepares the loop over all elements of one type for the computation of the error estimator
380 : template<typename TDomain>
381 : template<typename TElem, typename TFEGeom>
382 0 : void ConvectionDiffusionFE<TDomain>::
383 : prep_err_est_elem_loop(const ReferenceObjectID roid, const int si)
384 : {
385 : // get the error estimator data object and check that it is of the right type
386 : // we check this at this point in order to be able to dispense with this check later on
387 : // (i.e. in prep_err_est_elem and compute_err_est_A_elem())
388 0 : UG_COND_THROW(this->m_spErrEstData.get() == NULL,
389 : "No ErrEstData object has been given to this ElemDisc!");
390 :
391 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
392 :
393 0 : UG_COND_THROW(!err_est_data, "Dynamic cast to SideAndElemErrEstData failed." << std::endl
394 : << "Make sure you handed the correct type of ErrEstData to this discretization.");
395 :
396 :
397 : // set local positions
398 : static const int refDim = TElem::dim;
399 :
400 : // get local IPs
401 : size_t numSideIPs, numElemIPs;
402 : const MathVector<refDim>* sideIPs;
403 : const MathVector<refDim>* elemIPs;
404 : try
405 : {
406 : numSideIPs = err_est_data->num_all_side_ips(roid);
407 0 : numElemIPs = err_est_data->num_elem_ips(roid);
408 0 : sideIPs = err_est_data->template side_local_ips<refDim>(roid);
409 0 : elemIPs = err_est_data->template elem_local_ips<refDim>(roid);
410 :
411 0 : if (!sideIPs || !elemIPs) return; // are NULL if TElem is not of the same dim as TDomain
412 : }
413 0 : UG_CATCH_THROW("Integration points for error estimator cannot be set.");
414 :
415 : // set local IPs in imports
416 0 : m_imDiffusion.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
417 0 : m_imVelocity.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
418 0 : m_imFlux.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
419 0 : m_imSource.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
420 0 : m_imVectorSource.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
421 0 : m_imReactionRate.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
422 0 : m_imReaction.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
423 0 : m_imMassScale.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
424 0 : m_imMass.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
425 :
426 : // store values of shape functions in local IPs
427 : LagrangeP1<typename reference_element_traits<TElem>::reference_element_type> trialSpace
428 0 : = Provider<LagrangeP1<typename reference_element_traits<TElem>::reference_element_type> >::get();
429 :
430 0 : m_shapeValues.resize(numElemIPs, numSideIPs, trialSpace.num_sh());
431 0 : for (size_t ip = 0; ip < numElemIPs; ip++)
432 0 : trialSpace.shapes(m_shapeValues.shapesAtElemIP(ip), elemIPs[ip]);
433 0 : for (size_t ip = 0; ip < numSideIPs; ip++)
434 0 : trialSpace.shapes(m_shapeValues.shapesAtSideIP(ip), sideIPs[ip]);
435 : }
436 :
437 : template<typename TDomain>
438 : template<typename TElem, typename TFEGeom>
439 0 : void ConvectionDiffusionFE<TDomain>::
440 : prep_err_est_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
441 : {
442 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
443 :
444 : // request geometry
445 0 : TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
446 :
447 : try{
448 : geo.update(elem, vCornerCoords);
449 : }
450 0 : UG_CATCH_THROW("ConvectionDiffusion::prep_elem:"
451 : " Cannot update Finite Element Geometry.");
452 :
453 : // roid
454 0 : ReferenceObjectID roid = elem->reference_object_id();
455 :
456 : // set global positions
457 : size_t numSideIPs, numElemIPs;
458 : const MathVector<dim>* sideIPs;
459 : const MathVector<dim>* elemIPs;
460 :
461 : try
462 : {
463 : numSideIPs = err_est_data->num_all_side_ips(roid);
464 0 : numElemIPs = err_est_data->num_elem_ips(roid);
465 :
466 0 : sideIPs = err_est_data->all_side_global_ips(elem, vCornerCoords);
467 0 : elemIPs = err_est_data->elem_global_ips(elem, vCornerCoords);
468 : }
469 0 : UG_CATCH_THROW("Global integration points for error estimator cannot be set.");
470 :
471 0 : m_imDiffusion. set_global_ips(&sideIPs[0], numSideIPs);
472 0 : m_imVelocity. set_global_ips(&sideIPs[0], numSideIPs);
473 0 : m_imFlux. set_global_ips(&sideIPs[0], numSideIPs);
474 0 : m_imSource. set_global_ips(&elemIPs[0], numElemIPs);
475 0 : m_imVectorSource. set_global_ips(&sideIPs[0], numSideIPs);
476 0 : m_imReactionRate. set_global_ips(&elemIPs[0], numElemIPs);
477 0 : m_imReaction. set_global_ips(&elemIPs[0], numElemIPs);
478 0 : m_imMassScale. set_global_ips(&elemIPs[0], numElemIPs);
479 0 : m_imMass. set_global_ips(&elemIPs[0], numElemIPs);
480 0 : }
481 :
482 : // computes the error estimator contribution (stiffness part) for one element
483 : template<typename TDomain>
484 : template<typename TElem, typename TFEGeom>
485 0 : void ConvectionDiffusionFE<TDomain>::
486 : compute_err_est_A_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
487 : {
488 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
489 :
490 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
491 0 : UG_COND_THROW(err_est_data->surface_view().get() == NULL, "Error estimator has NULL surface view.");
492 :
493 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
494 :
495 : // Request geometry.
496 0 : static const TFEGeom& geo = GeomProvider<TFEGeom>::get();
497 :
498 :
499 : // SIDE TERMS //
500 :
501 : // get the sides of the element
502 : // We have to cast elem to a pointer of type SideAndElemErrEstData::elem_type
503 : // for the SideAndElemErrEstData::operator() to work properly.
504 : // This cannot generally be achieved by casting to TElem*, since this method is also registered for
505 : // lower-dimensional types TElem, and must therefore be compilable, even if it is never EVER to be executed.
506 : // The way we achieve this here, is by calling associated_elements_sorted() which has an implementation for
507 : // all possible types. Whatever comes out of it is of course complete nonsense if (and only if)
508 : // SideAndElemErrEstData::elem_type != TElem. To be on the safe side, we throw an error if the number of
509 : // entries in the list is not as it should be.
510 :
511 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::side_type>::secure_container side_list;
512 0 : pErrEstGrid->associated_elements_sorted(side_list, (TElem*) elem);
513 :
514 0 : UG_COND_THROW (side_list.size() != (size_t) ref_elem_type::numSides,
515 : "Mismatch of numbers of sides in 'ConvectionDiffusionFE::compute_err_est_elem'");
516 :
517 : // Some auxiliary variables.
518 : MathVector<dim> fluxDensity, gradC, normal;
519 :
520 : // FIXME: The computation of the gradient has to be reworked.
521 : // In the case of P1 shape functions, it is valid. For Q1 shape functions, however,
522 : // the gradient is not constant (but bilinear) on the element - and along the sides.
523 : // We cannot use the FVGeom here. Instead, we need to calculate the gradient in each IP!
524 :
525 : // Calculate grad(u) as average (over scvf).
526 : VecSet(gradC, 0.0);
527 0 : for(size_t ii = 0; ii < geo.num_ip(); ++ii)
528 : {
529 0 : for (size_t j=0; j<m_shapeValues.num_sh(); j++)
530 0 : VecScaleAppend(gradC, u(_C_,j), geo.global_grad(ii, j));
531 : }
532 0 : VecScale(gradC, gradC, (1.0/geo.num_ip()));
533 :
534 : // Calculate flux through the sides.
535 : size_t passedIPs = 0;
536 0 : for (size_t side=0; side < (size_t) ref_elem_type::numSides; side++)
537 : {
538 : // normal on side
539 0 : SideNormal<ref_elem_type,dim>(normal, side, vCornerCoords);
540 0 : VecNormalize(normal, normal);
541 :
542 : try
543 : {
544 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
545 : {
546 0 : size_t ip = passedIPs + sip;
547 :
548 : VecSet(fluxDensity, 0.0);
549 :
550 : // diffusion //
551 0 : if (m_imDiffusion.data_given())
552 : MatVecScaleMultAppend(fluxDensity, -1.0, m_imDiffusion[ip], gradC);
553 :
554 : // convection //
555 0 : if (m_imVelocity.data_given())
556 : {
557 : number val = 0.0;
558 0 : for (size_t sh = 0; sh < m_shapeValues.num_sh(); sh++)
559 0 : val += u(_C_,sh) * m_shapeValues.shapeAtSideIP(sh,sip);
560 :
561 : VecScaleAppend(fluxDensity, val, m_imVelocity[ip]);
562 : }
563 :
564 : // general flux //
565 0 : if (m_imFlux.data_given())
566 : VecAppend(fluxDensity, m_imFlux[ip]);
567 :
568 0 : (*err_est_data)(side_list[side],sip) += scale * VecDot(fluxDensity, normal);
569 : }
570 :
571 0 : passedIPs += err_est_data->num_side_ips(side_list[side]);
572 : }
573 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
574 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
575 : }
576 :
577 : // VOLUME TERMS //
578 :
579 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::elem_type>::secure_container elem_list;
580 : pErrEstGrid->associated_elements_sorted(elem_list, (TElem*) elem);
581 0 : UG_COND_THROW (elem_list.size() != 1, "Mismatch of numbers of sides in 'ConvectionDiffusionFE::compute_err_est_elem'");
582 :
583 : try
584 : {
585 0 : for (size_t ip = 0; ip < err_est_data->num_elem_ips(elem->reference_object_id()); ip++)
586 : {
587 : number total = 0.0;
588 :
589 : // diffusion // TODO ONLY FOR (PIECEWISE) CONSTANT DIFFUSION TENSOR SO FAR!
590 : // div(D*grad(c))
591 : // nothing to do, as c is piecewise linear and div(D*grad(c)) disappears
592 : // if D is diagonal and c bilinear, this should also vanish (confirm this!)
593 :
594 : // convection // TODO ONLY FOR (PIECEWISE) CONSTANT OR DIVERGENCE-FREE
595 : // VELOCITY FIELDS SO FAR!
596 : // div(v*c) = div(v)*c + v*grad(c) -- gradC has been calculated above
597 0 : if (m_imVelocity.data_given())
598 0 : total += VecDot(m_imVelocity[ip], gradC);
599 :
600 : // general flux // TODO ONLY FOR DIVERGENCE-FREE FLUX FIELD SO FAR!
601 : // nothing to do
602 :
603 : // reaction //
604 0 : if (m_imReactionRate.data_given())
605 : {
606 : number val = 0.0;
607 0 : for (size_t sh = 0; sh < geo.num_sh(); sh++)
608 0 : val += u(_C_,sh) * m_shapeValues.shapeAtElemIP(sh,ip);
609 :
610 0 : total += m_imReactionRate[ip] * val;
611 : }
612 :
613 0 : if (m_imReaction.data_given())
614 : {
615 0 : total += m_imReaction[ip];
616 : }
617 :
618 0 : (*err_est_data)(elem_list[0],ip) += scale * total;
619 : }
620 : }
621 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
622 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
623 0 : }
624 :
625 : // computes the error estimator contribution (mass part) for one element
626 : template<typename TDomain>
627 : template<typename TElem, typename TFEGeom>
628 0 : void ConvectionDiffusionFE<TDomain>::
629 : compute_err_est_M_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
630 : {
631 : // note: mass parts only enter volume term
632 :
633 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
634 0 : UG_COND_THROW(err_est_data->surface_view().get() == NULL,
635 : "Error estimator has NULL surface view.");
636 :
637 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
638 :
639 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::elem_type>::secure_container elem_list;
640 0 : pErrEstGrid->associated_elements_sorted(elem_list, (TElem*) elem);
641 0 : UG_COND_THROW (elem_list.size() != 1,
642 : "Mismatch of numbers of sides in 'ConvectionDiffusionFE::compute_err_est_elem'");
643 :
644 : // request geometry
645 0 : static const TFEGeom& geo = GeomProvider<TFEGeom>::get();
646 :
647 : // loop integration points
648 : try
649 : {
650 0 : for (size_t ip = 0; ip < err_est_data->num_elem_ips(elem->reference_object_id()); ip++)
651 : {
652 : number total = 0.0;
653 :
654 : // mass scale //
655 0 : if (m_imMassScale.data_given())
656 : {
657 : number val = 0.0;
658 0 : for (size_t sh = 0; sh < geo.num_sh(); sh++)
659 0 : val += u(_C_,sh) * m_shapeValues.shapeAtElemIP(sh,ip);
660 :
661 0 : total += m_imMassScale[ip] * val;
662 : }
663 :
664 : // mass //
665 0 : if (m_imMass.data_given())
666 : {
667 0 : total += m_imMass[ip];
668 : }
669 :
670 0 : (*err_est_data)(elem_list[0],ip) += scale * total;
671 : }
672 : }
673 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
674 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
675 0 : }
676 :
677 : // computes the error estimator contribution (rhs part) for one element
678 : template<typename TDomain>
679 : template<typename TElem, typename TFEGeom>
680 0 : void ConvectionDiffusionFE<TDomain>::
681 : compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
682 : {
683 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
684 :
685 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
686 0 : UG_COND_THROW(err_est_data->surface_view().get() == NULL, "Error estimator has NULL surface view.");
687 :
688 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
689 :
690 : // SIDE TERMS //
691 :
692 : // get the sides of the element
693 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::side_type>::secure_container side_list;
694 0 : pErrEstGrid->associated_elements_sorted(side_list, (TElem*) elem);
695 :
696 0 : UG_COND_THROW(side_list.size() != (size_t) ref_elem_type::numSides,
697 : "Mismatch of numbers of sides in 'ConvectionDiffusionFE::compute_err_est_elem'");
698 :
699 : // loop sides
700 : size_t passedIPs = 0;
701 0 : for (size_t side = 0; side < (size_t) ref_elem_type::numSides; side++)
702 : {
703 : // normal on side
704 : MathVector<dim> normal;
705 0 : SideNormal<ref_elem_type,dim>(normal, side, vCornerCoords);
706 0 : VecNormalize(normal, normal);
707 :
708 : try
709 : {
710 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
711 : {
712 0 : size_t ip = passedIPs + sip;
713 :
714 : // vector source //
715 0 : if (m_imVectorSource.data_given())
716 0 : (*err_est_data)(side_list[side],sip) += scale * VecDot(m_imVectorSource[ip], normal);
717 : }
718 :
719 0 : passedIPs += err_est_data->num_side_ips(side_list[side]);
720 : }
721 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
722 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
723 : }
724 :
725 : // VOLUME TERMS //
726 :
727 0 : if (!m_imSource.data_given()) return;
728 :
729 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::elem_type>::secure_container elem_list;
730 : pErrEstGrid->associated_elements_sorted(elem_list, (TElem*) elem);
731 :
732 0 : UG_COND_THROW (elem_list.size() != 1,
733 : "Mismatch of numbers of sides in 'ConvectionDiffusionFE::compute_err_est_elem'");
734 :
735 : // source //
736 : try
737 : {
738 0 : for (size_t ip = 0; ip < err_est_data->num_elem_ips(elem->reference_object_id()); ip++)
739 0 : (*err_est_data)(elem_list[0],ip) += scale * m_imSource[ip];
740 : }
741 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
742 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
743 : }
744 :
745 : // postprocesses the loop over all elements of one type in the computation of the error estimator
746 : template<typename TDomain>
747 : template<typename TElem, typename TFEGeom>
748 0 : void ConvectionDiffusionFE<TDomain>::
749 : fsh_err_est_elem_loop()
750 : {
751 : // finish the element loop in the same way as the actual discretization
752 : this->template fsh_elem_loop<TElem, TFEGeom> ();
753 0 : };
754 :
755 : // error estimation (end) ///
756 : // /////////////////////////////////
757 :
758 :
759 : // computes the linearized defect w.r.t to the velocity
760 : template<typename TDomain>
761 : template <typename TElem, typename TFEGeom>
762 0 : void ConvectionDiffusionFE<TDomain>::
763 : lin_def_velocity(const LocalVector& u,
764 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
765 : const size_t nip)
766 : {
767 : // request geometry
768 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
769 :
770 : // loop integration points
771 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
772 : {
773 : // get current u and grad_u
774 : number shape_u = 0.0;
775 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
776 0 : shape_u += u(_C_,j) * geo.shape(ip, j);
777 :
778 : // loop test spaces
779 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
780 : {
781 : // add to local defect
782 0 : VecScale(vvvLinDef[ip][_C_][i], geo.global_grad(ip, i),
783 0 : (-1)* geo.weight(ip) * shape_u);
784 : }
785 : }
786 0 : }
787 :
788 : // computes the linearized defect w.r.t to the flux
789 : template<typename TDomain>
790 : template <typename TElem, typename TFEGeom>
791 0 : void ConvectionDiffusionFE<TDomain>::
792 : lin_def_flux(const LocalVector& u,
793 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
794 : const size_t nip)
795 : {
796 : // request geometry
797 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
798 :
799 : // loop integration points
800 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
801 : {
802 : // loop test spaces
803 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
804 : {
805 : // add to local defect
806 0 : VecScale(vvvLinDef[ip][_C_][i], geo.global_grad(ip, i),
807 : (-1)* geo.weight(ip));
808 : }
809 : }
810 0 : }
811 :
812 : // computes the linearized defect w.r.t to the velocity
813 : template<typename TDomain>
814 : template <typename TElem, typename TFEGeom>
815 0 : void ConvectionDiffusionFE<TDomain>::
816 : lin_def_diffusion(const LocalVector& u,
817 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
818 : const size_t nip)
819 : {
820 : // request geometry
821 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
822 :
823 : MathVector<dim> grad_u;
824 :
825 : // loop integration points
826 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
827 : {
828 : // get current u and grad_u
829 : VecSet(grad_u, 0.0);
830 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
831 : VecScaleAppend(grad_u, u(_C_,j), geo.global_grad(ip, j));
832 :
833 : // loop test spaces
834 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
835 : {
836 0 : for(size_t k = 0; k < (size_t)dim; ++k)
837 0 : for(size_t j = 0; j < (size_t)dim; ++j)
838 0 : (vvvLinDef[ip][_C_][i])(k,j) = grad_u[j] * geo.global_grad(ip, i)[k]
839 0 : * geo.weight(ip);
840 : }
841 : }
842 0 : }
843 :
844 : // computes the linearized defect w.r.t to the reaction
845 : template<typename TDomain>
846 : template <typename TElem, typename TFEGeom>
847 0 : void ConvectionDiffusionFE<TDomain>::
848 : lin_def_reaction(const LocalVector& u,
849 : std::vector<std::vector<number> > vvvLinDef[],
850 : const size_t nip)
851 : {
852 : // request geometry
853 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
854 :
855 : // loop integration points
856 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
857 : {
858 : // loop test spaces
859 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
860 : {
861 : // compute contribution
862 0 : const number val = geo.shape(ip, i) * geo.weight(ip);
863 :
864 : // add to local defect
865 0 : vvvLinDef[ip][_C_][i] = val;
866 : }
867 : }
868 0 : }
869 :
870 : // computes the linearized defect w.r.t to the reaction
871 : template<typename TDomain>
872 : template <typename TElem, typename TFEGeom>
873 0 : void ConvectionDiffusionFE<TDomain>::
874 : lin_def_reaction_rate(const LocalVector& u,
875 : std::vector<std::vector<number> > vvvLinDef[],
876 : const size_t nip)
877 : {
878 : // request geometry
879 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
880 :
881 : // loop integration points
882 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
883 : {
884 : // compute value of current solution at ip
885 : number shape_u = 0.0;
886 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
887 0 : shape_u += u(_C_,j) * geo.shape(ip, j);
888 :
889 : // loop test spaces
890 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
891 : {
892 : // compute contribution
893 0 : const number val = shape_u * geo.shape(ip, i) * geo.weight(ip);
894 :
895 : // add to local defect
896 0 : vvvLinDef[ip][_C_][i] = val;
897 : }
898 : }
899 0 : }
900 :
901 :
902 : // computes the linearized defect w.r.t to the source
903 : template<typename TDomain>
904 : template <typename TElem, typename TFEGeom>
905 0 : void ConvectionDiffusionFE<TDomain>::
906 : lin_def_source(const LocalVector& u,
907 : std::vector<std::vector<number> > vvvLinDef[],
908 : const size_t nip)
909 : {
910 : // request geometry
911 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
912 :
913 : // loop integration points
914 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
915 : {
916 : // loop test spaces
917 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
918 : {
919 : // add contribution to local defect
920 0 : vvvLinDef[ip][_C_][i] = geo.shape(ip, i) * geo.weight(ip);
921 : }
922 : }
923 0 : }
924 :
925 : // computes the linearized defect w.r.t to the "vector source"
926 : template<typename TDomain>
927 : template <typename TElem, typename TFEGeom>
928 0 : void ConvectionDiffusionFE<TDomain>::
929 : lin_def_vector_source(const LocalVector& u,
930 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
931 : const size_t nip)
932 : {
933 : // request geometry
934 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
935 :
936 : // loop integration points
937 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
938 : {
939 : // loop test spaces
940 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
941 : {
942 : // add to local defect
943 0 : VecScale(vvvLinDef[ip][_C_][i], geo.global_grad(ip, i),
944 : geo.weight(ip));
945 : }
946 : }
947 0 : }
948 :
949 : // computes the linearized defect w.r.t to the mass scale
950 : template<typename TDomain>
951 : template <typename TElem, typename TFEGeom>
952 0 : void ConvectionDiffusionFE<TDomain>::
953 : lin_def_mass_scale(const LocalVector& u,
954 : std::vector<std::vector<number> > vvvLinDef[],
955 : const size_t nip)
956 : {
957 : // request geometry
958 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
959 :
960 : // loop integration points
961 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
962 : {
963 : // compute value of current solution at ip
964 : number shape_u = 0.0;
965 0 : for(size_t j = 0; j < geo.num_sh(); ++j)
966 0 : shape_u += u(_C_,j) * geo.shape(ip, j);
967 :
968 : // loop test spaces
969 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
970 : {
971 : // compute contribution
972 0 : const number val = shape_u * geo.shape(ip, i) * geo.weight(ip);
973 :
974 : // add to local defect
975 0 : vvvLinDef[ip][_C_][i] = val;
976 : }
977 : }
978 0 : }
979 :
980 : // computes the linearized defect w.r.t to the mass scale
981 : template<typename TDomain>
982 : template <typename TElem, typename TFEGeom>
983 0 : void ConvectionDiffusionFE<TDomain>::
984 : lin_def_mass(const LocalVector& u,
985 : std::vector<std::vector<number> > vvvLinDef[],
986 : const size_t nip)
987 : {
988 : // request geometry
989 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
990 :
991 : // loop integration points
992 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
993 : {
994 : // loop test spaces
995 0 : for(size_t i = 0; i < geo.num_sh(); ++i)
996 : {
997 : // compute contribution
998 0 : const number val = geo.shape(ip, i) * geo.weight(ip);
999 :
1000 : // add to local defect
1001 0 : vvvLinDef[ip][_C_][i] = val;
1002 : }
1003 : }
1004 0 : }
1005 :
1006 : // computes the linearized defect w.r.t to the velocity
1007 : template<typename TDomain>
1008 : template <typename TElem, typename TFEGeom>
1009 0 : void ConvectionDiffusionFE<TDomain>::
1010 : ex_value(number vValue[],
1011 : const MathVector<dim> vGlobIP[],
1012 : number time, int si,
1013 : const LocalVector& u,
1014 : GridObject* elem,
1015 : const MathVector<dim> vCornerCoords[],
1016 : const MathVector<TFEGeom::dim> vLocIP[],
1017 : const size_t nip,
1018 : bool bDeriv,
1019 : std::vector<std::vector<number> > vvvDeriv[])
1020 : {
1021 : // request geometry
1022 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
1023 :
1024 : // reference element
1025 : typedef typename reference_element_traits<TElem>::reference_element_type
1026 : ref_elem_type;
1027 :
1028 : // reference dimension
1029 : static const int refDim = reference_element_traits<TElem>::dim;
1030 :
1031 : // reference object id
1032 : static const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
1033 :
1034 : // FE ip
1035 0 : if(vLocIP == geo.local_ips())
1036 : {
1037 : // Loop ips
1038 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
1039 : {
1040 : // compute concentration at ip
1041 0 : vValue[ip] = 0.0;
1042 0 : for(size_t sh = 0; sh < geo.num_sh(); ++sh)
1043 0 : vValue[ip] += u(_C_, sh) * geo.shape(ip, sh);
1044 :
1045 : // compute derivative w.r.t. to unknowns iff needed
1046 0 : if(bDeriv)
1047 0 : for(size_t sh = 0; sh < geo.num_sh(); ++sh)
1048 0 : vvvDeriv[ip][_C_][sh] = geo.shape(ip, sh);
1049 : }
1050 : }
1051 : // general case
1052 : else
1053 : {
1054 : // request for trial space
1055 : try{
1056 : const LocalShapeFunctionSet<refDim>& rTrialSpace
1057 0 : = LocalFiniteElementProvider::get<refDim>(roid, m_lfeID);
1058 :
1059 : // number of shape functions
1060 0 : const size_t numSH = rTrialSpace.num_sh();
1061 :
1062 : // storage for shape function at ip
1063 0 : std::vector<number> vShape(numSH);
1064 :
1065 : // loop ips
1066 0 : for(size_t ip = 0; ip < nip; ++ip)
1067 : {
1068 : // evaluate at shapes at ip
1069 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
1070 :
1071 : // compute concentration at ip
1072 0 : vValue[ip] = 0.0;
1073 0 : for(size_t sh = 0; sh < numSH; ++sh)
1074 0 : vValue[ip] += u(_C_, sh) * vShape[sh];
1075 :
1076 : // compute derivative w.r.t. to unknowns iff needed
1077 : // \todo: maybe store shapes directly in vvvDeriv
1078 0 : if(bDeriv)
1079 0 : for(size_t sh = 0; sh < numSH; ++sh)
1080 0 : vvvDeriv[ip][_C_][sh] = vShape[sh];
1081 : }
1082 :
1083 0 : }
1084 0 : UG_CATCH_THROW("ConvectionDiffusion::ex_value: trial space missing.");
1085 : }
1086 0 : }
1087 :
1088 : template<typename TDomain>
1089 : template <typename TElem, typename TFEGeom>
1090 0 : void ConvectionDiffusionFE<TDomain>::
1091 : ex_grad(MathVector<dim> vValue[],
1092 : const MathVector<dim> vGlobIP[],
1093 : number time, int si,
1094 : const LocalVector& u,
1095 : GridObject* elem,
1096 : const MathVector<dim> vCornerCoords[],
1097 : const MathVector<TFEGeom::dim> vLocIP[],
1098 : const size_t nip,
1099 : bool bDeriv,
1100 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[])
1101 : {
1102 : // request geometry
1103 0 : const TFEGeom& geo = GeomProvider<TFEGeom>::get(m_lfeID, m_quadOrder);
1104 :
1105 : // reference element
1106 : typedef typename reference_element_traits<TElem>::reference_element_type
1107 : ref_elem_type;
1108 :
1109 : // reference dimension
1110 : static const int refDim = reference_element_traits<TElem>::dim;
1111 :
1112 : // reference object id
1113 : static const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
1114 :
1115 : // FE
1116 0 : if(vLocIP == geo.local_ips())
1117 : {
1118 : // Loop ip
1119 0 : for(size_t ip = 0; ip < geo.num_ip(); ++ip)
1120 : {
1121 0 : VecSet(vValue[ip], 0.0);
1122 0 : for(size_t sh = 0; sh < geo.num_sh(); ++sh)
1123 : VecScaleAppend(vValue[ip], u(_C_, sh), geo.global_grad(ip, sh));
1124 :
1125 0 : if(bDeriv)
1126 0 : for(size_t sh = 0; sh < geo.num_sh(); ++sh)
1127 0 : vvvDeriv[ip][_C_][sh] = geo.global_grad(ip, sh);
1128 : }
1129 : }
1130 : // general case
1131 : else
1132 : {
1133 : // request for trial space
1134 : try{
1135 : const LocalShapeFunctionSet<refDim>& rTrialSpace
1136 0 : = LocalFiniteElementProvider::get<refDim>(roid, m_lfeID);
1137 :
1138 : // number of shape functions
1139 0 : const size_t numSH = rTrialSpace.num_sh();
1140 :
1141 : // storage for shape function at ip
1142 0 : std::vector<MathVector<refDim> > vLocGrad(numSH);
1143 : MathVector<refDim> locGrad;
1144 :
1145 : // Reference Mapping
1146 : MathMatrix<dim, refDim> JTInv;
1147 0 : ReferenceMapping<ref_elem_type, dim> mapping(vCornerCoords);
1148 :
1149 : // loop ips
1150 0 : for(size_t ip = 0; ip < nip; ++ip)
1151 : {
1152 : // evaluate at shapes at ip
1153 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
1154 :
1155 : // compute grad at ip
1156 : VecSet(locGrad, 0.0);
1157 0 : for(size_t sh = 0; sh < numSH; ++sh)
1158 : VecScaleAppend(locGrad, u(_C_, sh), vLocGrad[sh]);
1159 :
1160 : // compute global grad
1161 0 : mapping.jacobian_transposed_inverse(JTInv, vLocIP[ip]);
1162 0 : MatVecMult(vValue[ip], JTInv, locGrad);
1163 :
1164 : // compute derivative w.r.t. to unknowns iff needed
1165 0 : if(bDeriv)
1166 0 : for(size_t sh = 0; sh < numSH; ++sh)
1167 0 : MatVecMult(vvvDeriv[ip][_C_][sh], JTInv, vLocGrad[sh]);
1168 : }
1169 0 : }
1170 0 : UG_CATCH_THROW("ConvectionDiffusion::ex_grad: trial space missing.");
1171 : }
1172 0 : };
1173 :
1174 : ////////////////////////////////////////////////////////////////////////////////
1175 : // register assemble functions
1176 : ////////////////////////////////////////////////////////////////////////////////
1177 :
1178 : #ifdef UG_DIM_1
1179 : template<>
1180 0 : void ConvectionDiffusionFE<Domain1d>::
1181 : register_all_funcs(const LFEID& lfeid, const int quadOrder)
1182 : {
1183 : // RegularEdge
1184 0 : register_func<RegularEdge, DimFEGeometry<dim> >();
1185 0 : }
1186 : #endif
1187 :
1188 : #ifdef UG_DIM_2
1189 : template<>
1190 0 : void ConvectionDiffusionFE<Domain2d>::
1191 : register_all_funcs(const LFEID& lfeid, const int quadOrder)
1192 : {
1193 0 : register_func<RegularEdge, DimFEGeometry<dim, 1> >();
1194 :
1195 : const int order = lfeid.order();
1196 0 : if(quadOrder != 2*order+1 || lfeid.type() != LFEID::LAGRANGE)
1197 : {
1198 0 : register_func<Triangle, DimFEGeometry<dim> >();
1199 0 : register_func<Quadrilateral, DimFEGeometry<dim> >();
1200 0 : return;
1201 : }
1202 :
1203 : // special compiled cases
1204 :
1205 : // Triangle
1206 0 : switch(order)
1207 : {
1208 0 : case 1: {typedef FEGeometry<Triangle, dim, LagrangeLSFS<ReferenceTriangle, 1>, GaussQuadrature<ReferenceTriangle, 3> > FEGeom;
1209 0 : register_func<Triangle, FEGeom >(); break;}
1210 0 : case 2: {typedef FEGeometry<Triangle, dim, LagrangeLSFS<ReferenceTriangle, 2>, GaussQuadrature<ReferenceTriangle, 5> > FEGeom;
1211 0 : register_func<Triangle, FEGeom >(); break;}
1212 0 : case 3: {typedef FEGeometry<Triangle, dim, LagrangeLSFS<ReferenceTriangle, 3>, GaussQuadrature<ReferenceTriangle, 7> > FEGeom;
1213 0 : register_func<Triangle, FEGeom >(); break;}
1214 0 : default: register_func<Triangle, DimFEGeometry<dim> >(); break;
1215 : }
1216 :
1217 : // Quadrilateral
1218 0 : switch(order) {
1219 0 : case 1: {typedef FEGeometry<Quadrilateral, dim, LagrangeLSFS<ReferenceQuadrilateral, 1>, GaussQuadrature<ReferenceQuadrilateral, 3> > FEGeom;
1220 0 : register_func<Quadrilateral, FEGeom >(); break;}
1221 0 : case 2: {typedef FEGeometry<Quadrilateral, dim, LagrangeLSFS<ReferenceQuadrilateral, 2>, GaussQuadrature<ReferenceQuadrilateral, 7> > FEGeom;
1222 0 : register_func<Quadrilateral, FEGeom >(); break;}
1223 0 : case 3: {typedef FEGeometry<Quadrilateral, dim, LagrangeLSFS<ReferenceQuadrilateral, 3>, GaussQuadrature<ReferenceQuadrilateral, 11> > FEGeom;
1224 0 : register_func<Quadrilateral, FEGeom >(); break;}
1225 0 : default: register_func<Quadrilateral, DimFEGeometry<dim> >(); break;
1226 : }
1227 : }
1228 : #endif
1229 :
1230 : #ifdef UG_DIM_3
1231 : template<>
1232 0 : void ConvectionDiffusionFE<Domain3d>::
1233 : register_all_funcs(const LFEID& lfeid, const int quadOrder)
1234 : {
1235 0 : register_func<RegularEdge, DimFEGeometry<dim, 1> >();
1236 0 : register_func<Triangle, DimFEGeometry<dim, 2> >();
1237 0 : register_func<Quadrilateral, DimFEGeometry<dim, 2> >();
1238 :
1239 : const int order = lfeid.order();
1240 0 : if(quadOrder != 2*order+1 || lfeid.type() != LFEID::LAGRANGE)
1241 : {
1242 0 : register_func<Tetrahedron, DimFEGeometry<dim> >();
1243 0 : register_func<Prism, DimFEGeometry<dim> >();
1244 0 : register_func<Pyramid, DimFEGeometry<dim> >();
1245 0 : register_func<Hexahedron, DimFEGeometry<dim> >();
1246 0 : register_func<Octahedron, DimFEGeometry<dim> >();
1247 0 : return;
1248 : }
1249 :
1250 : // special compiled cases
1251 :
1252 : // Tetrahedron
1253 0 : switch(order)
1254 : {
1255 0 : case 1: {typedef FEGeometry<Tetrahedron, dim, LagrangeLSFS<ReferenceTetrahedron, 1>, GaussQuadrature<ReferenceTetrahedron, 3> > FEGeom;
1256 0 : register_func<Tetrahedron, FEGeom >(); break;}
1257 0 : case 2: {typedef FEGeometry<Tetrahedron, dim, LagrangeLSFS<ReferenceTetrahedron, 2>, GaussQuadrature<ReferenceTetrahedron, 5> > FEGeom;
1258 0 : register_func<Tetrahedron, FEGeom >(); break;}
1259 0 : case 3: {typedef FEGeometry<Tetrahedron, dim, LagrangeLSFS<ReferenceTetrahedron, 3>, GaussQuadrature<ReferenceTetrahedron, 7> > FEGeom;
1260 0 : register_func<Tetrahedron, FEGeom >(); break;}
1261 0 : default: register_func<Tetrahedron, DimFEGeometry<dim> >(); break;
1262 : }
1263 :
1264 : // Prism
1265 : switch(order) {
1266 0 : default: register_func<Prism, DimFEGeometry<dim> >(); break;
1267 : }
1268 :
1269 : // Pyramid
1270 : switch(order)
1271 : {
1272 0 : default: register_func<Pyramid, DimFEGeometry<dim> >(); break;
1273 : }
1274 :
1275 : // Octahedron
1276 : switch(order)
1277 : {
1278 0 : default: register_func<Octahedron, DimFEGeometry<dim> >(); break;
1279 : }
1280 :
1281 : // Hexahedron
1282 0 : switch(order)
1283 : {
1284 0 : case 1: {typedef FEGeometry<Hexahedron, dim, LagrangeLSFS<ReferenceHexahedron, 1>, GaussQuadrature<ReferenceHexahedron, 3> > FEGeom;
1285 0 : register_func<Hexahedron, FEGeom >(); break;}
1286 0 : case 2: {typedef FEGeometry<Hexahedron, dim, LagrangeLSFS<ReferenceHexahedron, 2>, GaussQuadrature<ReferenceHexahedron, 7> > FEGeom;
1287 0 : register_func<Hexahedron, FEGeom >(); break;}
1288 0 : case 3: {typedef FEGeometry<Hexahedron, dim, LagrangeLSFS<ReferenceHexahedron, 3>, GaussQuadrature<ReferenceHexahedron, 11> > FEGeom;
1289 0 : register_func<Hexahedron, FEGeom >(); break;}
1290 0 : default: register_func<Hexahedron, DimFEGeometry<dim> >(); break;
1291 : }
1292 : }
1293 : #endif
1294 :
1295 : template <typename TDomain>
1296 : template <typename TElem, typename TFEGeom>
1297 0 : void ConvectionDiffusionFE<TDomain>::register_func()
1298 : {
1299 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1300 : typedef this_type T;
1301 : static const int refDim = reference_element_traits<TElem>::dim;
1302 :
1303 : this->clear_add_fct(id);
1304 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFEGeom>);
1305 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFEGeom>);
1306 : this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop<TElem, TFEGeom>);
1307 : this->set_add_jac_A_elem_fct(id, &T::template add_jac_A_elem<TElem, TFEGeom>);
1308 : this->set_add_jac_M_elem_fct(id, &T::template add_jac_M_elem<TElem, TFEGeom>);
1309 : this->set_add_def_A_elem_fct(id, &T::template add_def_A_elem<TElem, TFEGeom>);
1310 : this->set_add_def_M_elem_fct(id, &T::template add_def_M_elem<TElem, TFEGeom>);
1311 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFEGeom>);
1312 :
1313 : // error estimator parts
1314 : this->set_prep_err_est_elem_loop(id, &T::template prep_err_est_elem_loop<TElem, TFEGeom>);
1315 : this->set_prep_err_est_elem(id, &T::template prep_err_est_elem<TElem, TFEGeom>);
1316 : this->set_compute_err_est_A_elem(id, &T::template compute_err_est_A_elem<TElem, TFEGeom>);
1317 : this->set_compute_err_est_M_elem(id, &T::template compute_err_est_M_elem<TElem, TFEGeom>);
1318 : this->set_compute_err_est_rhs_elem(id, &T::template compute_err_est_rhs_elem<TElem, TFEGeom>);
1319 : this->set_fsh_err_est_elem_loop(id, &T::template fsh_err_est_elem_loop<TElem, TFEGeom>);
1320 :
1321 : // set computation of linearized defect w.r.t velocity
1322 0 : m_imDiffusion. set_fct(id, this, &T::template lin_def_diffusion<TElem, TFEGeom>);
1323 0 : m_imVelocity. set_fct(id, this, &T::template lin_def_velocity<TElem, TFEGeom>);
1324 0 : m_imFlux. set_fct(id, this, &T::template lin_def_flux<TElem, TFEGeom>);
1325 0 : m_imReactionRate. set_fct(id, this, &T::template lin_def_reaction_rate<TElem, TFEGeom>);
1326 0 : m_imReaction. set_fct(id, this, &T::template lin_def_reaction<TElem, TFEGeom>);
1327 0 : m_imSource. set_fct(id, this, &T::template lin_def_source<TElem, TFEGeom>);
1328 0 : m_imVectorSource. set_fct(id, this, &T::template lin_def_vector_source<TElem, TFEGeom>);
1329 0 : m_imMassScale. set_fct(id, this, &T::template lin_def_mass_scale<TElem, TFEGeom>);
1330 0 : m_imMass. set_fct(id, this, &T::template lin_def_mass<TElem, TFEGeom>);
1331 :
1332 : // exports
1333 0 : m_exValue-> template set_fct<T,refDim>(id, this, &T::template ex_value<TElem, TFEGeom>);
1334 0 : m_exGrad-> template set_fct<T,refDim>(id, this, &T::template ex_grad<TElem, TFEGeom>);
1335 0 : }
1336 :
1337 : ////////////////////////////////////////////////////////////////////////////////
1338 : // explicit template instantiations
1339 : ////////////////////////////////////////////////////////////////////////////////
1340 :
1341 : #ifdef UG_DIM_1
1342 : template class ConvectionDiffusionFE<Domain1d>;
1343 : #endif
1344 : #ifdef UG_DIM_2
1345 : template class ConvectionDiffusionFE<Domain2d>;
1346 : #endif
1347 : #ifdef UG_DIM_3
1348 : template class ConvectionDiffusionFE<Domain3d>;
1349 : #endif
1350 :
1351 : } // end namespace ConvectionDiffusionPlugin
1352 : } // namespace ug
1353 :
|
