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 : #include "convection_diffusion_fv.h"
34 :
35 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
36 : #include "lib_disc/spatial_disc/disc_util/fvho_geom.h"
37 :
38 : namespace ug{
39 : namespace ConvectionDiffusionPlugin{
40 :
41 : ////////////////////////////////////////////////////////////////////////////////
42 : // general
43 : ////////////////////////////////////////////////////////////////////////////////
44 :
45 : template<typename TDomain>
46 0 : ConvectionDiffusionFV<TDomain>::
47 : ConvectionDiffusionFV(const char* functions, const char* subsets)
48 : : ConvectionDiffusionBase<TDomain>(functions,subsets),
49 0 : m_bQuadOrderUserDef(false)
50 : {
51 : this->clear_add_fct();
52 0 : }
53 :
54 : template<typename TDomain>
55 0 : void ConvectionDiffusionFV<TDomain>::
56 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
57 : {
58 : // check grid
59 0 : if(bNonRegularGrid)
60 0 : UG_THROW("ConvectionDiffusion: Only regular grid implemented.");
61 :
62 : // check number
63 0 : if(vLfeID.size() != 1)
64 0 : UG_THROW("ConvectionDiffusion: Wrong number of functions given. "
65 : "Need exactly "<<1);
66 :
67 : // check that Lagrange
68 0 : if(vLfeID[0].type() != LFEID::LAGRANGE)
69 0 : UG_THROW("ConvectionDiffusion: Only Lagrange supported.");
70 :
71 : // check that not ADAPTIVE
72 0 : if(vLfeID[0].order() < 1)
73 0 : UG_THROW("ConvectionDiffusion: Adaptive order not implemented.");
74 :
75 : // set order
76 0 : m_lfeID = vLfeID[0];
77 0 : if(!m_bQuadOrderUserDef) {
78 0 : m_quadOrder = m_lfeID.order()+1;
79 : }
80 :
81 0 : register_all_funcs(m_lfeID, m_quadOrder);
82 0 : }
83 :
84 : template<typename TDomain>
85 0 : bool ConvectionDiffusionFV<TDomain>::
86 : use_hanging() const
87 : {
88 0 : return true;
89 : }
90 :
91 : template<typename TDomain>
92 0 : void ConvectionDiffusionFV<TDomain>::
93 : set_quad_order(size_t order)
94 : {
95 0 : m_quadOrder = order; m_bQuadOrderUserDef = true;
96 0 : }
97 :
98 : ////////////////////////////////////////////////////////////////////////////////
99 : // Assembling functions
100 : ////////////////////////////////////////////////////////////////////////////////
101 :
102 : template<typename TDomain>
103 : template<typename TElem, typename TFVGeom>
104 0 : void ConvectionDiffusionFV<TDomain>::
105 : prep_elem_loop(const ReferenceObjectID roid, const int si)
106 : {
107 0 : if( m_imSourceExpl.data_given() ||
108 0 : m_imReactionExpl.data_given() ||
109 : m_imReactionRateExpl.data_given())
110 0 : UG_THROW("ConvectionDiffusionFV: Explicit terms not implemented.");
111 :
112 : // request geometry
113 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
114 :
115 : try{
116 0 : geo.update_local(roid, m_lfeID, m_quadOrder);
117 : }
118 0 : UG_CATCH_THROW("ConvectionDiffusion::prep_elem_loop:"
119 : " Cannot update Finite Volume Geometry.");
120 :
121 : // set local positions
122 : if(!TFVGeom::usesHangingNodes)
123 : {
124 : static const int refDim = TElem::dim;
125 : const MathVector<refDim>* vSCVFip = geo.scvf_local_ips();
126 : const size_t numSCVFip = geo.num_scvf_ips();
127 : const MathVector<refDim>* vSCVip = geo.scv_local_ips();
128 : const size_t numSCVip = geo.num_scv_ips();
129 0 : m_imDiffusion.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
130 0 : m_imVelocity.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
131 0 : m_imFlux.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
132 0 : m_imSource.template set_local_ips<refDim>(vSCVip,numSCVip, false);
133 0 : m_imReactionRate.template set_local_ips<refDim>(vSCVip,numSCVip, false);
134 0 : m_imReaction.template set_local_ips<refDim>(vSCVip,numSCVip, false);
135 0 : m_imMassScale.template set_local_ips<refDim>(vSCVip,numSCVip, false);
136 0 : m_imMass.template set_local_ips<refDim>(vSCVip,numSCVip, false);
137 : }
138 0 : }
139 :
140 : template<typename TDomain>
141 : template<typename TElem, typename TFVGeom>
142 0 : void ConvectionDiffusionFV<TDomain>::
143 : fsh_elem_loop()
144 0 : {}
145 :
146 : template<typename TDomain>
147 : template<typename TElem, typename TFVGeom>
148 0 : void ConvectionDiffusionFV<TDomain>::
149 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
150 : {
151 : // request geometry
152 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
153 :
154 : try{
155 0 : geo.update(elem, vCornerCoords, &(this->subset_handler()));
156 : }
157 0 : UG_CATCH_THROW("ConvectionDiffusion::prep_elem:"
158 : " Cannot update Finite Volume Geometry.");
159 :
160 : // set local positions
161 : const size_t numSCVFip = geo.num_scvf_ips();
162 : const size_t numSCVip = geo.num_scv_ips();
163 : if(TFVGeom::usesHangingNodes)
164 : {
165 : static const int refDim = TElem::dim;
166 : const MathVector<refDim>* vSCVFip = geo.scvf_local_ips();
167 : const MathVector<refDim>* vSCVip = geo.scv_local_ips();
168 : m_imDiffusion.template set_local_ips<refDim>(vSCVFip,numSCVFip);
169 : m_imVelocity.template set_local_ips<refDim>(vSCVFip,numSCVFip);
170 : m_imFlux.template set_local_ips<refDim>(vSCVFip,numSCVFip);
171 : m_imSource.template set_local_ips<refDim>(vSCVip,numSCVip);
172 : m_imReactionRate.template set_local_ips<refDim>(vSCVip,numSCVip);
173 : m_imReaction.template set_local_ips<refDim>(vSCVip,numSCVip);
174 : m_imMassScale.template set_local_ips<refDim>(vSCVip,numSCVip);
175 : m_imMass.template set_local_ips<refDim>(vSCVip,numSCVip);
176 : }
177 :
178 : // set global positions
179 : const MathVector<dim>* vSCVFip = geo.scvf_global_ips();
180 : const MathVector<dim>* vSCVip = geo.scv_global_ips();
181 0 : m_imDiffusion. set_global_ips(vSCVFip, numSCVFip);
182 0 : m_imVelocity. set_global_ips(vSCVFip, numSCVFip);
183 0 : m_imFlux. set_global_ips(vSCVFip, numSCVFip);
184 0 : m_imSource. set_global_ips(vSCVip, numSCVip);
185 0 : m_imReactionRate. set_global_ips(vSCVip, numSCVip);
186 0 : m_imReaction. set_global_ips(vSCVip, numSCVip);
187 0 : m_imMassScale. set_global_ips(vSCVip, numSCVip);
188 0 : m_imMass. set_global_ips(vSCVip, numSCVip);
189 0 : }
190 :
191 : template<typename TDomain>
192 : template<typename TElem, typename TFVGeom>
193 0 : void ConvectionDiffusionFV<TDomain>::
194 : add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
195 : {
196 : // request geometry
197 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
198 :
199 : // Diff. Tensor times Gradient
200 : MathVector<dim> Dgrad;
201 :
202 : // Diffusion and Velocity Term
203 0 : if(m_imDiffusion.data_given() || m_imVelocity.data_given())
204 : {
205 : // loop Sub Control Volume Faces (SCVF)
206 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
207 : {
208 : // get current SCVF
209 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
210 :
211 : // loop integration points
212 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
213 : {
214 : ////////////////////////////////////////////////////
215 : // Diffusive Term
216 : ////////////////////////////////////////////////////
217 0 : if(m_imDiffusion.data_given())
218 : {
219 : // loop shape functions
220 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
221 : {
222 : // Compute Diffusion Tensor times Gradient
223 : MatVecMult(Dgrad, m_imDiffusion[ip], scvf.global_grad(i, sh));
224 :
225 : // Compute flux at IP
226 : const number D_diff_flux = VecDot(Dgrad, scvf.normal());
227 :
228 : // Add flux term to local matrix
229 0 : J(_C_, scvf.from(), _C_, sh) -= D_diff_flux * scvf.weight(i);
230 0 : J(_C_, scvf.to() , _C_, sh) += D_diff_flux * scvf.weight(i);
231 : }
232 : }
233 :
234 : ////////////////////////////////////////////////////
235 : // Convective Term
236 : ////////////////////////////////////////////////////
237 0 : if(m_imVelocity.data_given())
238 : {
239 : // Add Flux contribution
240 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
241 : {
242 0 : const number D_conv_flux = VecDot(m_imVelocity[ip], scvf.normal())
243 : * scvf.shape(i, sh);
244 :
245 : // Add fkux term to local matrix
246 0 : J(_C_, scvf.from(), _C_, sh) += D_conv_flux * scvf.weight(i);
247 0 : J(_C_, scvf.to(), _C_, sh) -= D_conv_flux * scvf.weight(i);
248 : }
249 : }
250 :
251 : // no explicit dependency on flux import
252 : } // end loop ip
253 : } // end loop scvf
254 : } // end data given
255 :
256 : ////////////////////////////////////////////////////
257 : // Reaction Term
258 : ////////////////////////////////////////////////////
259 :
260 : // if no data for reaction, return
261 0 : if(!m_imReactionRate.data_given()) return;
262 :
263 : // loop Sub Control Volume (SCV)
264 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
265 : {
266 : // get current SCV
267 : const typename TFVGeom::SCV& scv = geo.scv(s);
268 :
269 : // get associated node
270 0 : const int co = scv.node_id();
271 :
272 : // loop shapes
273 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
274 : {
275 : // reset integral
276 : number integral = 0;
277 :
278 : // loop integration points
279 0 : for(size_t i = 0; i < scv.num_ip(); ++i)
280 : {
281 0 : integral += m_imReactionRate[ip+i] * scv.shape(i, sh) * scv.weight(i);
282 : }
283 :
284 : // Add to local matrix
285 0 : J(_C_, co, _C_, sh) += integral;
286 : }
287 :
288 : // increase ip offset
289 0 : ip += scv.num_ip();
290 : }
291 :
292 : // no explicit dependency in m_imReaction
293 : }
294 :
295 :
296 : template<typename TDomain>
297 : template<typename TElem, typename TFVGeom>
298 0 : void ConvectionDiffusionFV<TDomain>::
299 : add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
300 : {
301 : // request geometry
302 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
303 :
304 0 : if(!m_imMassScale.data_given()) return;
305 :
306 : // loop Sub Control Volumes (SCV)
307 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
308 : {
309 : // get current SCV
310 : const typename TFVGeom::SCV& scv = geo.scv(s);
311 :
312 : // get associated node
313 0 : const int co = scv.node_id();
314 :
315 : // loop shapes
316 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
317 : {
318 : // reset integral
319 : number integral = 0;
320 :
321 : // loop integration points
322 0 : for(size_t i = 0; i < scv.num_ip(); ++i)
323 0 : integral += scv.shape(i, sh) * scv.weight(i) * m_imMassScale[ip+i];
324 :
325 : // no explicit dependency in m_imMass
326 :
327 : // Add to local matrix
328 0 : J(_C_, co, _C_, sh) += integral;
329 : }
330 :
331 : // increase ip offset
332 0 : ip += scv.num_ip();
333 : }
334 : }
335 :
336 :
337 : template<typename TDomain>
338 : template<typename TElem, typename TFVGeom>
339 0 : void ConvectionDiffusionFV<TDomain>::
340 : add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
341 : {
342 : // request geometry
343 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
344 :
345 0 : if(m_imDiffusion.data_given() || m_imVelocity.data_given() || m_imFlux.data_given())
346 : {
347 : // loop Sub Control Volume Faces (SCVF)
348 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
349 : {
350 : // get current SCVF
351 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
352 :
353 : // the flux of the scvf
354 : number flux = 0;
355 :
356 : // loop integration points
357 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
358 : {
359 : number fluxIP = 0;
360 :
361 : /////////////////////////////////////////////////////
362 : // Diffusive Term
363 : /////////////////////////////////////////////////////
364 0 : if(m_imDiffusion.data_given())
365 : {
366 : // to compute D \nabla c
367 : MathVector<dim> Dgrad_c, grad_c;
368 :
369 : // compute gradient and shape at ip
370 : VecSet(grad_c, 0.0);
371 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
372 : VecScaleAppend(grad_c, u(_C_,sh), scvf.global_grad(i, sh));
373 :
374 : // scale by diffusion tensor
375 : MatVecMult(Dgrad_c, m_imDiffusion[ip], grad_c);
376 :
377 : // Compute flux
378 0 : fluxIP = -VecDot(Dgrad_c, scvf.normal());
379 : }
380 :
381 : /////////////////////////////////////////////////////
382 : // Convective Term
383 : /////////////////////////////////////////////////////
384 0 : if(m_imVelocity.data_given())
385 : {
386 : // sum up solution
387 : number solIP = 0;
388 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
389 0 : solIP += u(_C_, sh) * scvf.shape(i, sh);
390 :
391 : // add convective flux
392 0 : fluxIP += solIP * VecDot(m_imVelocity[ip], scvf.normal());
393 : }
394 :
395 : /////////////////////////////////////////////////////
396 : // Flux Term
397 : /////////////////////////////////////////////////////
398 0 : if(m_imFlux.data_given())
399 : {
400 : // add flux
401 0 : fluxIP += VecDot(m_imFlux[ip], scvf.normal());
402 : }
403 :
404 : // sum flux
405 0 : flux += fluxIP * scvf.weight(i);
406 : } // end loop ip
407 :
408 : // no multiplication with volume is needed, since already contained
409 : // in the normal
410 :
411 : // add to local defect
412 0 : d(_C_, scvf.from()) += flux;
413 0 : d(_C_, scvf.to() ) -= flux;
414 :
415 : } // end loop scvf
416 : } // end data switch
417 :
418 : // reaction rate
419 0 : if(m_imReactionRate.data_given())
420 : {
421 : // loop Sub Control Volumes (SCV)
422 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
423 : {
424 : // get current SCV
425 : const typename TFVGeom::SCV& scv = geo.scv(s);
426 :
427 : // reset integral
428 : number integral = 0;
429 :
430 : // loop integration points
431 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
432 : {
433 : // compute solution at ip
434 : number solIP = 0;
435 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
436 0 : solIP += u(_C_, sh) * scv.shape(i, sh);
437 :
438 : // add to integral-sum
439 0 : integral += m_imReactionRate[ip] * solIP * scv.weight(i);
440 : }
441 :
442 : // get associated node
443 0 : const int co = scv.node_id();
444 :
445 : // Add to local defect
446 0 : d(_C_, co) += integral;
447 : }
448 : }
449 :
450 : // reaction rate
451 0 : if(m_imReaction.data_given())
452 : {
453 : // loop Sub Control Volumes (SCV)
454 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
455 : {
456 : // get current SCV
457 : const typename TFVGeom::SCV& scv = geo.scv(s);
458 :
459 : // reset integral
460 : number integral = 0;
461 :
462 : // loop integration points
463 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
464 : {
465 : // add to integral-sum
466 0 : integral += m_imReaction[ip] * scv.weight(i);
467 : }
468 :
469 : // get associated node
470 0 : const int co = scv.node_id();
471 :
472 : // Add to local defect
473 0 : d(_C_, co) += integral;
474 : }
475 : }
476 0 : }
477 :
478 :
479 : template<typename TDomain>
480 : template<typename TElem, typename TFVGeom>
481 0 : void ConvectionDiffusionFV<TDomain>::
482 : add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
483 : {
484 : // request geometry
485 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
486 :
487 0 : if(!m_imMassScale.data_given() && !m_imMass.data_given()) return;
488 :
489 : // loop Sub Control Volumes (SCV)
490 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
491 : {
492 : // get current SCV
493 : const typename TFVGeom::SCV& scv = geo.scv(s);
494 :
495 : // reset integral
496 : number integral = 0;
497 :
498 : // loop integration points
499 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
500 : {
501 : // compute value
502 : number val = 0.0;
503 :
504 : // multiply by scaling
505 0 : if(m_imMassScale.data_given()){
506 :
507 : number solIP = 0;
508 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
509 0 : solIP += u(_C_, sh) * scv.shape(i, sh);
510 :
511 0 : val += m_imMassScale[ip] * solIP;
512 : }
513 :
514 : // add mass
515 0 : if(m_imMass.data_given())
516 0 : val += m_imMass[ip];
517 :
518 : // add to integral-sum
519 0 : integral += val * scv.weight(i);
520 : }
521 :
522 : // get associated node
523 0 : const int co = scv.node_id();
524 :
525 : // Add to local defect
526 0 : d(_C_, co) += integral;
527 : }
528 : }
529 :
530 :
531 : template<typename TDomain>
532 : template<typename TElem, typename TFVGeom>
533 0 : void ConvectionDiffusionFV<TDomain>::
534 : add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[])
535 : {
536 : // if zero data given, return
537 0 : if(!m_imSource.data_given()) return;
538 :
539 : // request geometry
540 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
541 :
542 : // loop Sub Control Volumes (SCV)
543 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
544 : {
545 : // get current SCV
546 : const typename TFVGeom::SCV& scv = geo.scv(s);
547 :
548 : // reset integral
549 : number integral = 0;
550 :
551 : // loop integration points
552 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
553 : {
554 : // add to integral-sum
555 0 : integral += m_imSource[ip] * scv.weight(i);
556 : }
557 :
558 : // get associated node
559 0 : const int co = scv.node_id();
560 :
561 : // Add to local defect
562 0 : d(_C_, co) += integral;
563 : }
564 : }
565 :
566 :
567 : // computes the linearized defect w.r.t to the velocity
568 : template<typename TDomain>
569 : template <typename TElem, typename TFVGeom>
570 0 : void ConvectionDiffusionFV<TDomain>::
571 : lin_def_velocity(const LocalVector& u,
572 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
573 : const size_t nip)
574 : {
575 : // request geometry
576 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
577 :
578 : // reset the values for the linearized defect
579 0 : for(size_t ip = 0; ip < nip; ++ip)
580 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
581 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
582 : vvvLinDef[ip][c][sh] = 0.0;
583 :
584 : // loop Sub Control Volume Faces (SCVF)
585 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
586 : {
587 : // get current SCVF
588 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
589 :
590 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
591 : {
592 : // compute shape at ip
593 : number solIP = 0.0;
594 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
595 0 : solIP += u(_C_,sh) * scvf.shape(i, sh);
596 :
597 : // add parts for both sides of scvf
598 0 : VecScale(vvvLinDef[ip][_C_][scvf.from()], scvf.normal(), solIP * scvf.weight(i));
599 0 : VecScale(vvvLinDef[ip][_C_][scvf.to() ], scvf.normal(), -solIP * scvf.weight(i));
600 : }
601 : }
602 0 : }
603 :
604 : // computes the linearized defect w.r.t to the flux
605 : template<typename TDomain>
606 : template <typename TElem, typename TFVGeom>
607 0 : void ConvectionDiffusionFV<TDomain>::
608 : lin_def_flux(const LocalVector& u,
609 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
610 : const size_t nip)
611 : {
612 : // request geometry
613 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
614 :
615 : // reset the values for the linearized defect
616 0 : for(size_t ip = 0; ip < nip; ++ip)
617 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
618 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
619 : vvvLinDef[ip][c][sh] = 0.0;
620 :
621 : // loop Sub Control Volume Faces (SCVF)
622 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
623 : {
624 : // get current SCVF
625 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
626 :
627 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
628 : {
629 : // add parts for both sides of scvf
630 0 : VecScale(vvvLinDef[ip][_C_][scvf.from()], scvf.normal(), scvf.weight(i));
631 0 : VecScale(vvvLinDef[ip][_C_][scvf.to() ], scvf.normal(), -scvf.weight(i));
632 : }
633 : }
634 0 : }
635 :
636 : // computes the linearized defect w.r.t to the velocity
637 : template<typename TDomain>
638 : template <typename TElem, typename TFVGeom>
639 0 : void ConvectionDiffusionFV<TDomain>::
640 : lin_def_diffusion(const LocalVector& u,
641 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
642 : const size_t nip)
643 : {
644 : // request geometry
645 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
646 :
647 : // reset the values for the linearized defect
648 0 : for(size_t ip = 0; ip < nip; ++ip)
649 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
650 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
651 : vvvLinDef[ip][c][sh] = 0.0;
652 :
653 : // loop Sub Control Volume Faces (SCVF)
654 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
655 : {
656 : // get current SCVF
657 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
658 :
659 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
660 : {
661 : // compute gradient at ip
662 : MathVector<dim> gradIP; VecSet(gradIP, 0.0);
663 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
664 : VecScaleAppend(gradIP, u(_C_,sh), scvf.global_grad(i, sh));
665 :
666 : // part coming from -\nabla u * \vec{n}
667 0 : for(size_t k=0; k < (size_t)dim; ++k)
668 0 : for(size_t j = 0; j < (size_t)dim; ++j)
669 : {
670 0 : const number val = (scvf.normal())[k] * gradIP[j];
671 :
672 0 : vvvLinDef[ip][_C_][scvf.from()](k,j) = -val * scvf.weight(i);
673 0 : vvvLinDef[ip][_C_][scvf.to() ](k,j) = val * scvf.weight(i);
674 : }
675 : }
676 : }
677 0 : }
678 :
679 : // computes the linearized defect w.r.t to the reaction rate
680 : template<typename TDomain>
681 : template <typename TElem, typename TFVGeom>
682 0 : void ConvectionDiffusionFV<TDomain>::
683 : lin_def_reaction_rate(const LocalVector& u,
684 : std::vector<std::vector<number> > vvvLinDef[],
685 : const size_t nip)
686 : {
687 : // request geometry
688 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
689 :
690 : // loop Sub Control Volumes (SCV)
691 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
692 : {
693 : // get current SCV
694 : const typename TFVGeom::SCV& scv = geo.scv(s);
695 :
696 : // get associated node
697 0 : const int co = scv.node_id();
698 :
699 : // loop integration points
700 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
701 : {
702 : // compute solution at ip
703 : number solIP = 0;
704 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
705 0 : solIP += u(_C_, sh) * scv.shape(i, sh);
706 :
707 : // set lin defect
708 0 : vvvLinDef[ip][_C_][co] = solIP * scv.weight(i);
709 : }
710 : }
711 0 : }
712 :
713 : // computes the linearized defect w.r.t to the reaction
714 : template<typename TDomain>
715 : template <typename TElem, typename TFVGeom>
716 0 : void ConvectionDiffusionFV<TDomain>::
717 : lin_def_reaction(const LocalVector& u,
718 : std::vector<std::vector<number> > vvvLinDef[],
719 : const size_t nip)
720 : {
721 : // request geometry
722 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
723 :
724 : // loop Sub Control Volumes (SCV)
725 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
726 : {
727 : // get current SCV
728 : const typename TFVGeom::SCV& scv = geo.scv(s);
729 :
730 : // get associated node
731 0 : const int co = scv.node_id();
732 :
733 : // loop integration points
734 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
735 : {
736 : // set lin defect
737 0 : vvvLinDef[ip][_C_][co] = scv.weight(i);
738 : }
739 : }
740 0 : }
741 :
742 : // computes the linearized defect w.r.t to the source
743 : template<typename TDomain>
744 : template <typename TElem, typename TFVGeom>
745 0 : void ConvectionDiffusionFV<TDomain>::
746 : lin_def_source(const LocalVector& u,
747 : std::vector<std::vector<number> > vvvLinDef[],
748 : const size_t nip)
749 : {
750 : // request geometry
751 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
752 :
753 : // loop Sub Control Volumes (SCV)
754 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
755 : {
756 : // get current SCV
757 : const typename TFVGeom::SCV& scv = geo.scv(s);
758 :
759 : // get associated node
760 0 : const int co = scv.node_id();
761 :
762 : // loop integration points
763 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
764 : {
765 : // set lin defect
766 0 : vvvLinDef[ip][_C_][co] = scv.weight(i);
767 : }
768 : }
769 0 : }
770 :
771 : // computes the linearized defect w.r.t to the mass scale
772 : template<typename TDomain>
773 : template <typename TElem, typename TFVGeom>
774 0 : void ConvectionDiffusionFV<TDomain>::
775 : lin_def_mass_scale(const LocalVector& u,
776 : std::vector<std::vector<number> > vvvLinDef[],
777 : const size_t nip)
778 : {
779 : // request geometry
780 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
781 :
782 : // loop Sub Control Volumes (SCV)
783 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
784 : {
785 : // get current SCV
786 : const typename TFVGeom::SCV& scv = geo.scv(s);
787 :
788 : // get associated node
789 0 : const int co = scv.node_id();
790 :
791 : // loop integration points
792 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
793 : {
794 : // compute solution at ip
795 : number solIP = 0;
796 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
797 0 : solIP += u(_C_, sh) * scv.shape(i, sh);
798 :
799 : // set lin defect
800 0 : vvvLinDef[ip][_C_][co] = solIP * scv.weight(i);
801 : }
802 : }
803 0 : }
804 :
805 : // computes the linearized defect w.r.t to the mass scale
806 : template<typename TDomain>
807 : template <typename TElem, typename TFVGeom>
808 0 : void ConvectionDiffusionFV<TDomain>::
809 : lin_def_mass(const LocalVector& u,
810 : std::vector<std::vector<number> > vvvLinDef[],
811 : const size_t nip)
812 : {
813 : // request geometry
814 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
815 :
816 : // loop Sub Control Volumes (SCV)
817 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
818 : {
819 : // get current SCV
820 : const typename TFVGeom::SCV& scv = geo.scv(s);
821 :
822 : // get associated node
823 0 : const int co = scv.node_id();
824 :
825 : // loop integration points
826 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
827 : {
828 : // set lin defect
829 0 : vvvLinDef[ip][_C_][co] = scv.weight(i);
830 : }
831 : }
832 0 : }
833 :
834 : // computes the linearized defect w.r.t to the velocity
835 : template<typename TDomain>
836 : template <typename TElem, typename TFVGeom>
837 0 : void ConvectionDiffusionFV<TDomain>::
838 : ex_value(number vValue[],
839 : const MathVector<dim> vGlobIP[],
840 : number time, int si,
841 : const LocalVector& u,
842 : GridObject* elem,
843 : const MathVector<dim> vCornerCoords[],
844 : const MathVector<TFVGeom::dim> vLocIP[],
845 : const size_t nip,
846 : bool bDeriv,
847 : std::vector<std::vector<number> > vvvDeriv[])
848 : {
849 : // request geometry
850 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
851 :
852 : // reference element
853 : typedef typename reference_element_traits<TElem>::reference_element_type
854 : ref_elem_type;
855 :
856 : // reference object id
857 : static const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
858 :
859 : // FV SCVF ip
860 0 : if(vLocIP == geo.scvf_local_ips())
861 : {
862 : // loop Sub Control Volume Faces (SCVF)
863 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
864 : {
865 : // get current SCVF
866 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
867 :
868 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
869 : {
870 : // compute concentration at ip
871 0 : vValue[ip] = 0.0;
872 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
873 0 : vValue[ip] += u(_C_, sh) * scvf.shape(i, sh);
874 :
875 : // compute derivative w.r.t. to unknowns iff needed
876 0 : if(bDeriv)
877 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
878 0 : vvvDeriv[ip][_C_][sh] = scvf.shape(i, sh);
879 : }
880 : }
881 : }
882 : // FV SCV ip
883 0 : else if(vLocIP == geo.scv_local_ips())
884 : {
885 : // loop Sub Control Volumes (SCV)
886 0 : for(size_t s = 0, ip = 0; s < geo.num_scv(); ++s)
887 : {
888 : // get current SCV
889 : const typename TFVGeom::SCV& scv = geo.scv(s);
890 :
891 : // loop integration points
892 0 : for(size_t i = 0; i < scv.num_ip(); ++i, ++ip)
893 : {
894 : // compute solution at ip
895 0 : vValue[ip] = 0.0;
896 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
897 0 : vValue[ip] += u(_C_, sh) * scv.shape(i, sh);
898 :
899 : // compute derivative w.r.t. to unknowns iff needed
900 0 : if(bDeriv)
901 0 : for(size_t sh = 0; sh < scv.num_sh(); ++sh)
902 0 : vvvDeriv[ip][_C_][sh] = scv.shape(i, sh);
903 : }
904 : }
905 : }
906 : // general case
907 : else
908 : {
909 : // request for trial space
910 : try{
911 : const LocalShapeFunctionSet<TFVGeom::dim>& rTrialSpace
912 0 : = LocalFiniteElementProvider::get<TFVGeom::dim>(roid, m_lfeID);
913 :
914 : // storage for shape function at ip
915 0 : std::vector<number> vShape(rTrialSpace.num_sh());
916 :
917 : // loop ips
918 0 : for(size_t ip = 0; ip < nip; ++ip)
919 : {
920 : // evaluate at shapes at ip
921 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
922 :
923 : // compute concentration at ip
924 0 : vValue[ip] = 0.0;
925 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
926 0 : vValue[ip] += u(_C_, sh) * vShape[sh];
927 :
928 : // compute derivative w.r.t. to unknowns iff needed
929 : // \todo: maybe store shapes directly in vvvDeriv
930 0 : if(bDeriv)
931 0 : for(size_t sh = 0; sh < vShape.size(); ++sh){
932 : UG_ASSERT(_C_ < vvvDeriv[ip].size(), _C_<<", "<<vvvDeriv[ip].size());
933 : UG_ASSERT(sh < vvvDeriv[ip][_C_].size(), sh<<", "<<vvvDeriv[ip][_C_].size());
934 0 : vvvDeriv[ip][_C_][sh] = vShape[sh];
935 : }
936 : }
937 0 : }
938 0 : UG_CATCH_THROW("ConvectionDiffusion::ex_value: trial space missing.");
939 : }
940 0 : }
941 :
942 : // computes the linearized defect w.r.t to the velocity
943 : template<typename TDomain>
944 : template <typename TElem, typename TFVGeom>
945 0 : void ConvectionDiffusionFV<TDomain>::
946 : ex_grad(MathVector<dim> vValue[],
947 : const MathVector<dim> vGlobIP[],
948 : number time, int si,
949 : const LocalVector& u,
950 : GridObject* elem,
951 : const MathVector<dim> vCornerCoords[],
952 : const MathVector<TFVGeom::dim> vLocIP[],
953 : const size_t nip,
954 : bool bDeriv,
955 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[])
956 : {
957 : // request geometry
958 0 : const TFVGeom& geo = GeomProvider<TFVGeom>::get(m_lfeID, m_quadOrder);
959 :
960 : // reference element
961 : typedef typename reference_element_traits<TElem>::reference_element_type
962 : ref_elem_type;
963 :
964 : // reference dimension
965 : static const int refDim = TElem::dim;
966 :
967 : // reference object id
968 : static const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
969 :
970 : // FV SCVF ip
971 0 : if(vLocIP == geo.scvf_local_ips())
972 : {
973 : // loop Sub Control Volume Faces (SCVF)
974 0 : for(size_t s = 0, ip = 0; s < geo.num_scvf(); ++s)
975 : {
976 : // get current SCVF
977 : const typename TFVGeom::SCVF& scvf = geo.scvf(s);
978 :
979 0 : for(size_t i = 0; i < scvf.num_ip(); ++i, ++ip)
980 : {
981 0 : VecSet(vValue[ip], 0.0);
982 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
983 : VecScaleAppend(vValue[ip], u(_C_, sh), scvf.global_grad(i, sh));
984 :
985 0 : if(bDeriv)
986 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
987 0 : vvvDeriv[ip][_C_][sh] = scvf.global_grad(i, sh);
988 : }
989 : }
990 : }
991 : // general case
992 : else
993 : {
994 : // request for trial space
995 : try{
996 : const LocalShapeFunctionSet<TFVGeom::dim>& rTrialSpace
997 0 : = LocalFiniteElementProvider::get<TFVGeom::dim>(roid, m_lfeID);
998 :
999 : // storage for shape function at ip
1000 0 : const size_t numSH = rTrialSpace.num_sh();
1001 0 : std::vector<MathVector<refDim> > vLocGrad(numSH);
1002 : MathVector<refDim> locGrad;
1003 :
1004 : // Reference Mapping
1005 : MathMatrix<dim, refDim> JTInv;
1006 0 : ReferenceMapping<ref_elem_type, dim> mapping(vCornerCoords);
1007 :
1008 : // loop ips
1009 0 : for(size_t ip = 0; ip < nip; ++ip)
1010 : {
1011 : // evaluate at shapes at ip
1012 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
1013 :
1014 : // compute grad at ip
1015 : VecSet(locGrad, 0.0);
1016 0 : for(size_t sh = 0; sh < numSH; ++sh)
1017 : VecScaleAppend(locGrad, u(_C_, sh), vLocGrad[sh]);
1018 :
1019 : // compute global grad
1020 0 : mapping.jacobian_transposed_inverse(JTInv, vLocIP[ip]);
1021 0 : MatVecMult(vValue[ip], JTInv, locGrad);
1022 :
1023 : // compute derivative w.r.t. to unknowns iff needed
1024 0 : if(bDeriv)
1025 0 : for(size_t sh = 0; sh < numSH; ++sh)
1026 0 : MatVecMult(vvvDeriv[ip][_C_][sh], JTInv, vLocGrad[sh]);
1027 : }
1028 :
1029 0 : }
1030 0 : UG_CATCH_THROW("ConvectionDiffusion::ex_grad: trial space missing");
1031 : }
1032 0 : };
1033 :
1034 : ////////////////////////////////////////////////////////////////////////////////
1035 : // register assemble functions
1036 : ////////////////////////////////////////////////////////////////////////////////
1037 :
1038 : #ifdef UG_DIM_1
1039 : template<>
1040 0 : void ConvectionDiffusionFV<Domain1d>::
1041 : register_all_funcs(const LFEID& lfeID, const int quadOrder)
1042 : {
1043 : // const int order = lfeID.order();
1044 : typedef DimFVGeometry<dim> FVGeom;
1045 0 : register_func<RegularEdge, FVGeom >();
1046 0 : }
1047 : #endif
1048 :
1049 : #ifdef UG_DIM_2
1050 : template<>
1051 0 : void ConvectionDiffusionFV<Domain2d>::
1052 : register_all_funcs(const LFEID& lfeID, const int quadOrder)
1053 : {
1054 : const int order = lfeID.order();
1055 0 : if(quadOrder == order+1 && lfeID.type() == LFEID::LAGRANGE)
1056 : {
1057 : // RegularEdge
1058 0 : switch(order)
1059 : {
1060 0 : case 1: {typedef FVGeometry<1, RegularEdge, dim> FVGeom;
1061 0 : register_func<RegularEdge, FVGeom >(); break;}
1062 0 : case 2: {typedef FVGeometry<2, RegularEdge, dim> FVGeom;
1063 0 : register_func<RegularEdge, FVGeom >(); break;}
1064 0 : case 3: {typedef FVGeometry<3, RegularEdge, dim> FVGeom;
1065 0 : register_func<RegularEdge, FVGeom >(); break;}
1066 0 : default: {typedef DimFVGeometry<dim, 1> FVGeom;
1067 0 : register_func<RegularEdge, FVGeom >(); break;}
1068 : }
1069 :
1070 : // Triangle
1071 0 : switch(order)
1072 : {
1073 0 : case 1: {typedef FVGeometry<1, Triangle, dim> FVGeom;
1074 0 : register_func<Triangle, FVGeom >(); break;}
1075 0 : case 2: {typedef FVGeometry<2, Triangle, dim> FVGeom;
1076 0 : register_func<Triangle, FVGeom >(); break;}
1077 0 : case 3: {typedef FVGeometry<3, Triangle, dim> FVGeom;
1078 0 : register_func<Triangle, FVGeom >(); break;}
1079 0 : default: {typedef DimFVGeometry<dim> FVGeom;
1080 0 : register_func<Triangle, FVGeom >(); break;}
1081 : }
1082 :
1083 : // Quadrilateral
1084 0 : switch(order) {
1085 0 : case 1: {typedef FVGeometry<1, Quadrilateral, dim> FVGeom;
1086 0 : register_func<Quadrilateral, FVGeom >(); break;}
1087 0 : case 2: {typedef FVGeometry<2, Quadrilateral, dim> FVGeom;
1088 0 : register_func<Quadrilateral, FVGeom >(); break;}
1089 0 : case 3: {typedef FVGeometry<3, Quadrilateral, dim> FVGeom;
1090 0 : register_func<Quadrilateral, FVGeom >(); break;}
1091 0 : default: {typedef DimFVGeometry<dim> FVGeom;
1092 0 : register_func<Quadrilateral, FVGeom >(); break;}
1093 : }
1094 : }
1095 : else
1096 : {
1097 0 : register_func<RegularEdge, DimFVGeometry<dim, 1> >();
1098 :
1099 : typedef DimFVGeometry<dim> FVGeom;
1100 0 : register_func<Triangle, FVGeom >();
1101 0 : register_func<Quadrilateral, FVGeom >();
1102 : }
1103 0 : }
1104 : #endif
1105 :
1106 : #ifdef UG_DIM_3
1107 : template<>
1108 0 : void ConvectionDiffusionFV<Domain3d>::
1109 : register_all_funcs(const LFEID& lfeID, const int quadOrder)
1110 : {
1111 : const int order = lfeID.order();
1112 0 : if(quadOrder == order+1 && lfeID.type() == LFEID::LAGRANGE)
1113 : {
1114 : // RegularEdge
1115 0 : switch(order)
1116 : {
1117 0 : case 1: {typedef FVGeometry<1, RegularEdge, dim> FVGeom;
1118 0 : register_func<RegularEdge, FVGeom >(); break;}
1119 0 : case 2: {typedef FVGeometry<2, RegularEdge, dim> FVGeom;
1120 0 : register_func<RegularEdge, FVGeom >(); break;}
1121 0 : case 3: {typedef FVGeometry<3, RegularEdge, dim> FVGeom;
1122 0 : register_func<RegularEdge, FVGeom >(); break;}
1123 0 : default: {typedef DimFVGeometry<dim, 1> FVGeom;
1124 0 : register_func<RegularEdge, FVGeom >(); break;}
1125 : }
1126 :
1127 : // Tetrahedron
1128 0 : switch(order)
1129 : {
1130 0 : case 1: {typedef FVGeometry<1, Tetrahedron, dim> FVGeom;
1131 0 : register_func<Tetrahedron, FVGeom >(); break;}
1132 0 : case 2: {typedef FVGeometry<2, Tetrahedron, dim> FVGeom;
1133 0 : register_func<Tetrahedron, FVGeom >(); break;}
1134 0 : case 3: {typedef FVGeometry<3, Tetrahedron, dim> FVGeom;
1135 0 : register_func<Tetrahedron, FVGeom >(); break;}
1136 0 : default: {typedef DimFVGeometry<dim> FVGeom;
1137 0 : register_func<Tetrahedron, FVGeom >(); break;}
1138 : }
1139 :
1140 : // Prism
1141 0 : switch(order) {
1142 0 : case 1: {typedef FVGeometry<1, Prism, dim> FVGeom;
1143 0 : register_func<Prism, FVGeom >(); break;}
1144 0 : default: {typedef DimFVGeometry<dim> FVGeom;
1145 0 : register_func<Prism, FVGeom >(); break;}
1146 : }
1147 :
1148 : // Hexahedron
1149 0 : switch(order)
1150 : {
1151 0 : case 1: {typedef FVGeometry<1, Hexahedron, dim> FVGeom;
1152 0 : register_func<Hexahedron, FVGeom >(); break;}
1153 0 : case 2: {typedef FVGeometry<2, Hexahedron, dim> FVGeom;
1154 0 : register_func<Hexahedron, FVGeom >(); break;}
1155 0 : case 3: {typedef FVGeometry<3, Hexahedron, dim> FVGeom;
1156 0 : register_func<Hexahedron, FVGeom >(); break;}
1157 0 : default: {typedef DimFVGeometry<dim> FVGeom;
1158 0 : register_func<Hexahedron, FVGeom >(); break;}
1159 : }
1160 : }
1161 : else
1162 : {
1163 : typedef DimFVGeometry<dim> FVGeom;
1164 0 : register_func<Tetrahedron, FVGeom >();
1165 0 : register_func<Prism, FVGeom >();
1166 0 : register_func<Hexahedron, FVGeom >();
1167 : }
1168 :
1169 0 : }
1170 : #endif
1171 :
1172 : template<typename TDomain>
1173 : template<typename TElem, typename TFVGeom>
1174 0 : void ConvectionDiffusionFV<TDomain>::
1175 : register_func()
1176 : {
1177 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1178 : typedef this_type T;
1179 : static const int refDim = reference_element_traits<TElem>::dim;
1180 :
1181 : this->clear_add_fct(id);
1182 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFVGeom>);
1183 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFVGeom>);
1184 : this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop<TElem, TFVGeom>);
1185 : this->set_add_jac_A_elem_fct(id, &T::template add_jac_A_elem<TElem, TFVGeom>);
1186 : this->set_add_jac_M_elem_fct(id, &T::template add_jac_M_elem<TElem, TFVGeom>);
1187 : this->set_add_def_A_elem_fct(id, &T::template add_def_A_elem<TElem, TFVGeom>);
1188 : this->set_add_def_M_elem_fct(id, &T::template add_def_M_elem<TElem, TFVGeom>);
1189 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFVGeom>);
1190 :
1191 : // set computation of linearized defect w.r.t velocity
1192 0 : m_imDiffusion. set_fct(id, this, &T::template lin_def_diffusion<TElem, TFVGeom>);
1193 0 : m_imVelocity. set_fct(id, this, &T::template lin_def_velocity<TElem, TFVGeom>);
1194 0 : m_imFlux. set_fct(id, this, &T::template lin_def_flux<TElem, TFVGeom>);
1195 0 : m_imReactionRate.set_fct(id, this, &T::template lin_def_reaction_rate<TElem, TFVGeom>);
1196 0 : m_imReaction. set_fct(id, this, &T::template lin_def_reaction<TElem, TFVGeom>);
1197 0 : m_imSource. set_fct(id, this, &T::template lin_def_source<TElem, TFVGeom>);
1198 0 : m_imMassScale.set_fct(id, this, &T::template lin_def_mass_scale<TElem, TFVGeom>);
1199 0 : m_imMass. set_fct(id, this, &T::template lin_def_mass<TElem, TFVGeom>);
1200 :
1201 : // exports
1202 0 : m_exValue-> template set_fct<T,refDim>(id, this, &T::template ex_value<TElem, TFVGeom>);
1203 0 : m_exGrad->template set_fct<T,refDim>(id, this, &T::template ex_grad<TElem, TFVGeom>);
1204 0 : }
1205 :
1206 : ////////////////////////////////////////////////////////////////////////////////
1207 : // explicit template instantiations
1208 : ////////////////////////////////////////////////////////////////////////////////
1209 :
1210 : #ifdef UG_DIM_1
1211 : template class ConvectionDiffusionFV<Domain1d>;
1212 : #endif
1213 : #ifdef UG_DIM_2
1214 : template class ConvectionDiffusionFV<Domain2d>;
1215 : #endif
1216 : #ifdef UG_DIM_3
1217 : template class ConvectionDiffusionFV<Domain3d>;
1218 : #endif
1219 :
1220 : } // end namespace ConvectionDiffusionPlugin
1221 : } // namespace ug
1222 :
|
