Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Christian Wehner
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_fvcr.h"
34 :
35 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
36 : #include "lib_disc/spatial_disc/disc_util/fvcr_geom.h"
37 : #include "lib_disc/spatial_disc/disc_util/hfvcr_geom.h"
38 : #include "lib_disc/spatial_disc/disc_util/conv_shape.h"
39 :
40 : namespace ug{
41 : namespace ConvectionDiffusionPlugin{
42 :
43 : ////////////////////////////////////////////////////////////////////////////////
44 : // general
45 : ////////////////////////////////////////////////////////////////////////////////
46 :
47 : template<typename TDomain>
48 0 : ConvectionDiffusionFVCR<TDomain>::
49 : ConvectionDiffusionFVCR(const char* functions, const char* subsets)
50 : : ConvectionDiffusionBase<TDomain>(functions,subsets),
51 0 : m_spConvShape(new ConvectionShapesNoUpwind<dim>),
52 0 : m_bNonRegularGrid(false)
53 : {
54 0 : register_all_funcs(m_bNonRegularGrid);
55 0 : }
56 :
57 : template<typename TDomain>
58 0 : void ConvectionDiffusionFVCR<TDomain>::
59 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
60 : {
61 : // check number
62 0 : if(vLfeID.size() != 1)
63 0 : UG_THROW("ConvectionDiffusion: Wrong number of functions given. "
64 : "Need exactly "<<1);
65 :
66 0 : if(vLfeID[0].order() != 1 || vLfeID[0].type() != LFEID::CROUZEIX_RAVIART)
67 0 : UG_THROW("ConvectionDiffusion FVCR Scheme only implemented for 1st order.");
68 :
69 : // remember
70 0 : m_bNonRegularGrid = bNonRegularGrid;
71 :
72 : // update assemble functions
73 0 : register_all_funcs(m_bNonRegularGrid);
74 0 : }
75 :
76 : template<typename TDomain>
77 0 : bool ConvectionDiffusionFVCR<TDomain>::
78 : use_hanging() const
79 : {
80 0 : return true;
81 : }
82 :
83 :
84 : ////////////////////////////////////////////////////////////////////////////////
85 : // Assembling functions
86 : ////////////////////////////////////////////////////////////////////////////////
87 :
88 :
89 : template<typename TDomain>
90 : template<typename TElem, typename TFVGeom>
91 0 : void ConvectionDiffusionFVCR<TDomain>::
92 : prep_elem_loop(const ReferenceObjectID roid, const int si)
93 : {
94 0 : if(m_imFlux.data_given())
95 0 : UG_THROW("ConvectionDiffusionFVCR: Flux import not implemented.");
96 :
97 0 : if( m_imSourceExpl.data_given() ||
98 0 : m_imReactionExpl.data_given() ||
99 : m_imReactionRateExpl.data_given())
100 0 : UG_THROW("ConvectionDiffusionFVCR: Explicit terms not implemented.");
101 :
102 : // set local positions
103 : if(!TFVGeom::usesHangingNodes)
104 : {
105 : static const int refDim = TElem::dim;
106 0 : static TFVGeom& geo = GeomProvider<TFVGeom>::get();
107 :
108 0 : geo.update_local_data();
109 :
110 0 : m_imDiffusion.template set_local_ips<refDim>(geo.scvf_local_ips(),
111 : geo.num_scvf_ips(), false);
112 0 : m_imVelocity.template set_local_ips<refDim>(geo.scvf_local_ips(),
113 : geo.num_scvf_ips(), false);
114 0 : m_imSource.template set_local_ips<refDim>(geo.scv_local_ips(),
115 : geo.num_scv_ips(), false);
116 0 : m_imReactionRate.template set_local_ips<refDim>(geo.scv_local_ips(),
117 : geo.num_scv_ips(), false);
118 0 : m_imReaction.template set_local_ips<refDim>(geo.scv_local_ips(),
119 : geo.num_scv_ips(), false);
120 0 : m_imMassScale.template set_local_ips<refDim>(geo.scv_local_ips(),
121 : geo.num_scv_ips(), false);
122 0 : m_imMass.template set_local_ips<refDim>(geo.scv_local_ips(),
123 : geo.num_scv_ips(), false);
124 : }
125 :
126 : // check, that upwind has been set
127 : // if(m_spConvShape.invalid())
128 : // UG_THROW("ConvectionDiffusion::prep_elem_loop:"
129 : // " Upwind has not been set.");
130 :
131 : // init upwind for element type
132 : // if(!m_spConvShape->template set_geometry_type<TFVGeom>())
133 : // UG_THROW("ConvectionDiffusion::prep_elem_loop:"
134 : // " Cannot init upwind for element type.");
135 0 : }
136 :
137 : template<typename TDomain>
138 : template<typename TElem, typename TFVGeom>
139 0 : void ConvectionDiffusionFVCR<TDomain>::
140 : fsh_elem_loop()
141 0 : {}
142 :
143 : template<typename TDomain>
144 : template<typename TElem, typename TFVGeom>
145 0 : void ConvectionDiffusionFVCR<TDomain>::
146 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
147 : {
148 : // Update Geometry for this element
149 0 : static TFVGeom& geo = GeomProvider<TFVGeom>::get();
150 :
151 : try{
152 0 : geo.update(elem, vCornerCoords, &(this->subset_handler()));
153 : }
154 0 : UG_CATCH_THROW("ConvectionDiffusion::prep_elem:"
155 : " Cannot update Finite Volume Geometry.");
156 :
157 : // set local positions
158 : if(TFVGeom::usesHangingNodes)
159 : {
160 : static const int refDim = TElem::dim;
161 0 : m_imDiffusion.template set_local_ips<refDim>(geo.scvf_local_ips(),
162 : geo.num_scvf_ips());
163 0 : m_imVelocity.template set_local_ips<refDim>(geo.scvf_local_ips(),
164 : geo.num_scvf_ips());
165 0 : m_imSource.template set_local_ips<refDim>(geo.scv_local_ips(),
166 : geo.num_scv_ips());
167 0 : m_imReactionRate.template set_local_ips<refDim>(geo.scv_local_ips(),
168 : geo.num_scv_ips());
169 0 : m_imReaction.template set_local_ips<refDim>(geo.scv_local_ips(),
170 : geo.num_scv_ips());
171 0 : m_imMassScale.template set_local_ips<refDim>(geo.scv_local_ips(),
172 : geo.num_scv_ips());
173 0 : m_imMass.template set_local_ips<refDim>(geo.scv_local_ips(),
174 : geo.num_scv_ips());
175 : /* if(m_spConvShape.valid())
176 : if(!m_spConvShape->template set_geometry_type<TFVGeom>(geo))
177 : UG_THROW("ConvectionDiffusion::prep_elem_loop:"
178 : " Cannot init upwind for element type.");*/
179 : }
180 :
181 : // set global positions
182 0 : const MathVector<dim>* vSCVFip = geo.scvf_global_ips();
183 : const size_t numSCVFip = geo.num_scvf_ips();
184 : const MathVector<dim>* vSCVip = geo.scv_global_ips();
185 : const size_t numSCVip = geo.num_scv_ips();
186 0 : m_imDiffusion. set_global_ips(vSCVFip, numSCVFip);
187 0 : m_imVelocity. set_global_ips(vSCVFip, numSCVFip);
188 : // m_imFlux. set_global_ips(vSCVFip, numSCVFip);
189 0 : m_imSource. set_global_ips(vSCVip, numSCVip);
190 : // m_imVectorSource. set_global_ips(vSCVFip, numSCVFip);
191 0 : m_imReactionRate. set_global_ips(vSCVip, numSCVip);
192 : // m_imReactionRateExpl. set_global_ips(vSCVip, numSCVip);
193 : // m_imReactionExpl. set_global_ips(vSCVip, numSCVip);
194 : // m_imSourceExpl. set_global_ips(vSCVip, numSCVip);
195 0 : m_imReaction. set_global_ips(vSCVip, numSCVip);
196 : // m_imMassScale. set_global_ips(vSCVip, numSCVip);
197 0 : m_imMass. set_global_ips(vSCVip, numSCVip);
198 0 : }
199 :
200 : template<typename TDomain>
201 : template<typename TElem, typename TFVGeom>
202 0 : void ConvectionDiffusionFVCR<TDomain>::
203 : add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
204 : {
205 : // get finite volume geometry
206 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
207 :
208 : // Diff. Tensor times Gradient
209 : MathVector<dim> Dgrad;
210 :
211 : // get conv shapes
212 : // const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo);
213 :
214 : // Diffusion and Velocity Term
215 0 : if(m_imDiffusion.data_given() || m_imVelocity.data_given())
216 : {
217 : // loop Sub Control Volume Faces (SCVF)
218 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
219 : {
220 : // get current SCVF
221 0 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
222 :
223 : ////////////////////////////////////////////////////
224 : // Diffusive Term
225 : ////////////////////////////////////////////////////
226 0 : if(m_imDiffusion.data_given())
227 : {
228 : // loop shape functions
229 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
230 : {
231 : // Compute Diffusion Tensor times Gradient
232 : MatVecMult(Dgrad, m_imDiffusion[ip], scvf.global_grad(sh));
233 :
234 : // Compute flux at IP
235 : const number D_diff_flux = VecDot(Dgrad, scvf.normal());
236 :
237 : // Add flux term to local matrix
238 0 : J(_C_, scvf.from(), _C_, sh) -= D_diff_flux;
239 0 : J(_C_, scvf.to() , _C_, sh) += D_diff_flux;
240 : }
241 : }
242 :
243 : ////////////////////////////////////////////////////
244 : // Convective Term
245 : ////////////////////////////////////////////////////
246 : /* if(m_imVelocity.data_given())
247 : {
248 : // Add Flux contribution
249 : for(size_t sh = 0; sh < convShape.num_sh(); ++sh)
250 : {
251 : const number D_conv_flux = convShape(ip, sh);
252 :
253 : // Add flux term to local matrix
254 : J(_C_, scvf.from(), _C_, sh) += D_conv_flux;
255 : J(_C_, scvf.to(), _C_, sh) -= D_conv_flux;
256 : }
257 : }*/
258 : }
259 : }
260 :
261 : // handle constrained dofs
262 : if(TFVGeom::usesHangingNodes){
263 0 : for (size_t i=0;i<geo.num_constrained_dofs();i++){
264 : const typename TFVGeom::CONSTRAINED_DOF& cd = geo.constrained_dof(i);
265 : const size_t index = cd.index();
266 0 : J(_C_,index,_C_,index) = 1;
267 0 : for (size_t j=0;j<cd.num_constraining_dofs();j++)
268 0 : J(_C_, index, _C_, cd.constraining_dofs_index(j)) = -cd.constraining_dofs_weight(j);
269 : // insert interpolation equation directly for all dofs
270 0 : for (size_t j=0;j<geo.num_scv();j++){
271 : const size_t nodeID = geo.scv(j).node_id();
272 0 : number alpha=J(_C_,nodeID,_C_,index);
273 0 : J(_C_,nodeID,_C_,index)=0;
274 0 : for (size_t k=0;k<cd.num_constraining_dofs();k++)
275 0 : J(_C_,nodeID,_C_,cd.constraining_dofs_index(k)) += alpha*cd.constraining_dofs_weight(k);
276 : }
277 : }
278 : }
279 :
280 : /* UG_LOG("Local Matrix is: \n");
281 : size_t n = geo.num_scv()+geo.num_constrained_dofs();
282 : UG_LOG("locA=[");
283 : for (size_t i=0;i<n;i++){
284 : for (size_t j=0;j<n;j++)
285 : UG_LOG(J(0,i,0,j) << " ");
286 : UG_LOG(";\n");
287 : }
288 : UG_LOG("]\n");*/
289 : ////////////////////////////////////////////////////
290 : // Reaction Term (using lumping)
291 : ////////////////////////////////////////////////////
292 :
293 : // if no data for reaction rate given, return
294 0 : if(!m_imReactionRate.data_given()) return;
295 :
296 : // loop Sub Control Volume (SCV)
297 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
298 : {
299 : // get current SCV
300 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
301 :
302 : // get associated node
303 0 : const int co = scv.node_id();
304 :
305 : // Add to local matrix
306 0 : J(_C_, co, _C_, co) += m_imReactionRate[ip] * scv.volume();
307 : }
308 :
309 : // reaction term does not explicitly depend on the associated unknown function
310 : }
311 :
312 :
313 : template<typename TDomain>
314 : template<typename TElem, typename TFVGeom>
315 0 : void ConvectionDiffusionFVCR<TDomain>::
316 : add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
317 : {
318 : // get finite volume geometry
319 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
320 :
321 0 : if(!m_imMassScale.data_given()) return;
322 :
323 : // loop Sub Control Volumes (SCV)
324 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
325 : {
326 : // get current SCV
327 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
328 :
329 : // get associated node
330 0 : const int co = scv.node_id();
331 :
332 : // Add to local matrix
333 0 : J(_C_, co, _C_, co) += scv.volume() * m_imMassScale[ip];
334 : }
335 :
336 : // m_imMass part does not explicitly depend on associated unknown function
337 : }
338 :
339 :
340 : template<typename TDomain>
341 : template<typename TElem, typename TFVGeom>
342 0 : void ConvectionDiffusionFVCR<TDomain>::
343 : add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
344 : {
345 : // get finite volume geometry
346 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
347 :
348 : // get conv shapes
349 : // const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo);
350 :
351 0 : if(m_imDiffusion.data_given() || m_imVelocity.data_given())
352 : {
353 : // loop Sub Control Volume Faces (SCVF)
354 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
355 : {
356 : // get current SCVF
357 0 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
358 :
359 : /////////////////////////////////////////////////////
360 : // Diffusive Term
361 : /////////////////////////////////////////////////////
362 0 : if(m_imDiffusion.data_given())
363 : {
364 : // to compute D \nabla c
365 : MathVector<dim> Dgrad_c, grad_c;
366 :
367 : // compute gradient and shape at ip
368 : VecSet(grad_c, 0.0);
369 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
370 : VecScaleAppend(grad_c, u(_C_,sh), scvf.global_grad(sh));
371 :
372 : // scale by diffusion tensor
373 : MatVecMult(Dgrad_c, m_imDiffusion[ip], grad_c);
374 :
375 : // Compute flux
376 : const number diff_flux = VecDot(Dgrad_c, scvf.normal());
377 :
378 : // Add to local defect
379 0 : d(_C_, scvf.from()) -= diff_flux;
380 0 : d(_C_, scvf.to() ) += diff_flux;
381 : }
382 :
383 : /////////////////////////////////////////////////////
384 : // Convective Term
385 : /////////////////////////////////////////////////////
386 : /* if(m_imVelocity.data_given())
387 : {
388 : // sum up convective flux using convection shapes
389 : number conv_flux = 0.0;
390 : for(size_t sh = 0; sh < convShape.num_sh(); ++sh)
391 : conv_flux += u(_C_, sh) * convShape(ip, sh);
392 :
393 : // add to local defect
394 : d(_C_, scvf.from()) += conv_flux;
395 : d(_C_, scvf.to() ) -= conv_flux;
396 : }*/
397 : }
398 : }
399 :
400 : // reaction rate
401 0 : if(m_imReactionRate.data_given())
402 : {
403 : // loop Sub Control Volumes (SCV)
404 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
405 : {
406 : // get current SCV
407 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
408 :
409 : // get associated node
410 0 : const int co = scv.node_id();
411 :
412 : // Add to local defect
413 0 : d(_C_, co) += u(_C_, co) * m_imReactionRate[ip] * scv.volume();
414 : }
415 : }
416 :
417 : // reaction rate
418 0 : if(m_imReaction.data_given())
419 : {
420 : // loop Sub Control Volumes (SCV)
421 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
422 : {
423 : // get current SCV
424 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
425 :
426 : // get associated node
427 0 : const int co = scv.node_id();
428 :
429 : // Add to local defect
430 0 : d(_C_, co) += m_imReaction[ip] * scv.volume();
431 : }
432 : }
433 :
434 : // handle constrained dofs, compute defect of interpolation equation
435 : if(TFVGeom::usesHangingNodes){
436 0 : for (size_t i=0;i<geo.num_constrained_dofs();i++){
437 : const typename TFVGeom::CONSTRAINED_DOF& cd = geo.constrained_dof(i);
438 : const size_t index = cd.index();
439 : number defect = u(_C_,index);
440 0 : for (size_t j=0;j<cd.num_constraining_dofs();j++)
441 0 : defect -= cd.constraining_dofs_weight(j) * u(_C_,cd.constraining_dofs_index(j));
442 0 : d(_C_,index) = defect;
443 : }
444 : }
445 0 : }
446 :
447 :
448 : template<typename TDomain>
449 : template<typename TElem, typename TFVGeom>
450 0 : void ConvectionDiffusionFVCR<TDomain>::
451 : add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
452 : {
453 : // get finite volume geometry
454 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
455 :
456 0 : if(!m_imMassScale.data_given() && !m_imMass.data_given()) return;
457 :
458 : // loop Sub Control Volumes (SCV)
459 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
460 : {
461 : // get current SCV
462 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
463 :
464 : // get associated node
465 0 : const int co = scv.node_id();
466 :
467 : // mass value
468 : number val = 0.0;
469 :
470 : // multiply by scaling
471 0 : if(m_imMassScale.data_given())
472 0 : val += m_imMassScale[ip] * u(_C_, co);
473 :
474 : // add mass
475 0 : if(m_imMass.data_given())
476 0 : val += m_imMass[ip];
477 :
478 : // Add to local defect
479 0 : d(_C_, co) += val * scv.volume();
480 : }
481 : }
482 :
483 :
484 : template<typename TDomain>
485 : template<typename TElem, typename TFVGeom>
486 0 : void ConvectionDiffusionFVCR<TDomain>::
487 : add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[])
488 : {
489 : // if zero data given, return
490 0 : if(!m_imSource.data_given()) return;
491 :
492 : // get finite volume geometry
493 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
494 :
495 : // loop Sub Control Volumes (SCV)
496 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
497 : {
498 : // get current SCV
499 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
500 :
501 : // get associated node
502 0 : const int co = scv.node_id();
503 :
504 : // Add to local rhs
505 0 : d(_C_, co) += m_imSource[ip] * scv.volume();
506 : }
507 : }
508 :
509 :
510 : // computes the linearized defect w.r.t to the velocity
511 : template<typename TDomain>
512 : template <typename TElem, typename TFVGeom>
513 0 : void ConvectionDiffusionFVCR<TDomain>::
514 : lin_def_velocity(const LocalVector& u,
515 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
516 : const size_t nip)
517 : {
518 : // get finite volume geometry
519 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
520 :
521 : // get conv shapes
522 : // const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo);
523 :
524 : // reset the values for the linearized defect
525 0 : for(size_t ip = 0; ip < nip; ++ip)
526 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
527 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
528 : vvvLinDef[ip][c][sh] = 0.0;
529 :
530 : // loop Sub Control Volume Faces (SCVF)
531 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
532 : {
533 : // get current SCVF
534 0 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
535 :
536 : // sum up contributions of convection shapes
537 : MathVector<dim> linDefect;
538 : VecSet(linDefect, 0.0);
539 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
540 : // VecScaleAppend(linDefect, u(_C_,sh), convShape.D_vel(ip, sh));
541 :
542 : // add parts for both sides of scvf
543 0 : vvvLinDef[ip][_C_][scvf.from()] += linDefect;
544 0 : vvvLinDef[ip][_C_][scvf.to()] -= linDefect;
545 : }
546 0 : }
547 :
548 : // computes the linearized defect w.r.t to the velocity
549 : template<typename TDomain>
550 : template <typename TElem, typename TFVGeom>
551 0 : void ConvectionDiffusionFVCR<TDomain>::
552 : lin_def_diffusion(const LocalVector& u,
553 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
554 : const size_t nip)
555 : {
556 : // get finite volume geometry
557 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
558 :
559 : // get conv shapes
560 : // const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo);
561 :
562 : // reset the values for the linearized defect
563 0 : for(size_t ip = 0; ip < nip; ++ip)
564 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
565 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
566 : vvvLinDef[ip][c][sh] = 0.0;
567 :
568 : // loop Sub Control Volume Faces (SCVF)
569 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
570 : {
571 : // get current SCVF
572 0 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
573 :
574 : // compute gradient at ip
575 : MathVector<dim> grad_u; VecSet(grad_u, 0.0);
576 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
577 : VecScaleAppend(grad_u, u(_C_,sh), scvf.global_grad(sh));
578 :
579 : // compute the lin defect at this ip
580 : MathMatrix<dim,dim> linDefect;
581 :
582 : // part coming from -\nabla u * \vec{n}
583 0 : for(size_t k=0; k < (size_t)dim; ++k)
584 0 : for(size_t j = 0; j < (size_t)dim; ++j)
585 0 : linDefect(j,k) = (scvf.normal())[j] * grad_u[k];
586 :
587 : // add contribution from convection shapes
588 : // if(convShape.non_zero_deriv_diffusion())
589 : // for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
590 : // MatAdd(linDefect, convShape.D_diffusion(ip, sh), u(_C_, sh));
591 :
592 : // add contributions
593 0 : vvvLinDef[ip][_C_][scvf.from()] -= linDefect;
594 : vvvLinDef[ip][_C_][scvf.to() ] += linDefect;
595 : }
596 0 : }
597 :
598 : // computes the linearized defect w.r.t to the reaction rate
599 : template<typename TDomain>
600 : template <typename TElem, typename TFVGeom>
601 0 : void ConvectionDiffusionFVCR<TDomain>::
602 : lin_def_reaction_rate(const LocalVector& u,
603 : std::vector<std::vector<number> > vvvLinDef[],
604 : const size_t nip)
605 : {
606 : // get finite volume geometry
607 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
608 :
609 : // loop Sub Control Volumes (SCV)
610 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
611 : {
612 : // get current SCV
613 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
614 :
615 : // get associated node
616 0 : const int co = scv.node_id();
617 :
618 : // set lin defect
619 0 : vvvLinDef[ip][_C_][co] = u(_C_, co) * scv.volume();
620 : }
621 0 : }
622 :
623 : // computes the linearized defect w.r.t to the reaction
624 : template<typename TDomain>
625 : template <typename TElem, typename TFVGeom>
626 0 : void ConvectionDiffusionFVCR<TDomain>::
627 : lin_def_reaction(const LocalVector& u,
628 : std::vector<std::vector<number> > vvvLinDef[],
629 : const size_t nip)
630 : {
631 : // get finite volume geometry
632 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
633 :
634 : // loop Sub Control Volumes (SCV)
635 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
636 : {
637 : // get current SCV
638 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
639 :
640 : // get associated node
641 0 : const int co = scv.node_id();
642 :
643 : // set lin defect
644 0 : vvvLinDef[ip][_C_][co] = scv.volume();
645 : }
646 0 : }
647 :
648 : // computes the linearized defect w.r.t to the source
649 : template<typename TDomain>
650 : template <typename TElem, typename TFVGeom>
651 0 : void ConvectionDiffusionFVCR<TDomain>::
652 : lin_def_source(const LocalVector& u,
653 : std::vector<std::vector<number> > vvvLinDef[],
654 : const size_t nip)
655 : {
656 : // get finite volume geometry
657 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
658 :
659 : // loop Sub Control Volumes (SCV)
660 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
661 : {
662 : // get current SCV
663 0 : const typename TFVGeom::SCV& scv = geo.scv(ip);
664 :
665 : // get associated node
666 0 : const int co = scv.node_id();
667 :
668 : // set lin defect
669 0 : vvvLinDef[ip][_C_][co] = scv.volume();
670 : }
671 0 : }
672 :
673 : // computes the linearized defect w.r.t to the mass scale
674 : template<typename TDomain>
675 : template <typename TElem, typename TFVGeom>
676 0 : void ConvectionDiffusionFVCR<TDomain>::
677 : lin_def_mass_scale(const LocalVector& u,
678 : std::vector<std::vector<number> > vvvLinDef[],
679 : const size_t nip)
680 : {
681 : // get finite volume geometry
682 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
683 :
684 : // loop Sub Control Volumes (SCV)
685 0 : for(size_t co = 0; co < geo.num_scv(); ++co)
686 : {
687 : // get current SCV
688 0 : const typename TFVGeom::SCV& scv = geo.scv(co);
689 :
690 : // Check associated node
691 : UG_ASSERT(co == scv.node_id(), "Only one shape per SCV");
692 :
693 : // set lin defect
694 0 : vvvLinDef[co][_C_][co] = u(_C_, co) * scv.volume();
695 : }
696 0 : }
697 :
698 : // computes the linearized defect w.r.t to the mass scale
699 : template<typename TDomain>
700 : template <typename TElem, typename TFVGeom>
701 0 : void ConvectionDiffusionFVCR<TDomain>::
702 : lin_def_mass(const LocalVector& u,
703 : std::vector<std::vector<number> > vvvLinDef[],
704 : const size_t nip)
705 : {
706 : // get finite volume geometry
707 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
708 :
709 : // loop Sub Control Volumes (SCV)
710 0 : for(size_t co = 0; co < geo.num_scv(); ++co)
711 : {
712 : // get current SCV
713 0 : const typename TFVGeom::SCV& scv = geo.scv(co);
714 :
715 : // Check associated node
716 : UG_ASSERT(co == scv.node_id(), "Only one shape per SCV");
717 :
718 : // set lin defect
719 0 : vvvLinDef[co][_C_][co] = scv.volume();
720 : }
721 0 : }
722 :
723 : // computes the linearized defect w.r.t to the velocity
724 : template<typename TDomain>
725 : template <typename TElem, typename TFVGeom>
726 0 : void ConvectionDiffusionFVCR<TDomain>::
727 : ex_value(number vValue[],
728 : const MathVector<dim> vGlobIP[],
729 : number time, int si,
730 : const LocalVector& u,
731 : GridObject* elem,
732 : const MathVector<dim> vCornerCoords[],
733 : const MathVector<TFVGeom::dim> vLocIP[],
734 : const size_t nip,
735 : bool bDeriv,
736 : std::vector<std::vector<number> > vvvDeriv[])
737 : {
738 : // get finite volume geometry
739 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
740 :
741 : // reference element
742 : typedef typename reference_element_traits<TElem>::reference_element_type
743 : ref_elem_type;
744 :
745 : // number of shape functions
746 : static const size_t numSH = ref_elem_type::numCorners;
747 :
748 : // CRFV SCVF ip
749 0 : if(vLocIP == geo.scvf_local_ips())
750 : {
751 : // Loop Sub Control Volume Faces (SCVF)
752 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
753 : {
754 : // Get current SCVF
755 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
756 :
757 : // compute concentration at ip
758 0 : vValue[ip] = 0.0;
759 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
760 0 : vValue[ip] += u(_C_, sh) * scvf.shape(sh);
761 :
762 : // compute derivative w.r.t. to unknowns iff needed
763 0 : if(bDeriv)
764 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
765 0 : vvvDeriv[ip][_C_][sh] = scvf.shape(sh);
766 : }
767 : }
768 : // CRFV SCV ip
769 0 : else if(vLocIP == geo.scv_local_ips())
770 : {
771 : // solution at ip
772 0 : for(size_t sh = 0; sh < numSH; ++sh)
773 0 : vValue[sh] = u(_C_, sh);
774 :
775 : // set derivatives if needed
776 0 : if(bDeriv)
777 0 : for(size_t sh = 0; sh < numSH; ++sh)
778 0 : for(size_t sh2 = 0; sh2 < numSH; ++sh2)
779 0 : vvvDeriv[sh][_C_][sh2] = (sh==sh2) ? 1.0 : 0.0;
780 : }
781 : // general case
782 : else
783 : {
784 : // get trial space
785 0 : LagrangeP1<ref_elem_type> rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
786 :
787 : // storage for shape function at ip
788 : number vShape[numSH];
789 :
790 : // loop ips
791 0 : for(size_t ip = 0; ip < nip; ++ip)
792 : {
793 : // evaluate at shapes at ip
794 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
795 :
796 : // compute concentration at ip
797 0 : vValue[ip] = 0.0;
798 0 : for(size_t sh = 0; sh < numSH; ++sh)
799 0 : vValue[ip] += u(_C_, sh) * vShape[sh];
800 :
801 : // compute derivative w.r.t. to unknowns iff needed
802 : // \todo: maybe store shapes directly in vvvDeriv
803 0 : if(bDeriv)
804 0 : for(size_t sh = 0; sh < numSH; ++sh)
805 0 : vvvDeriv[ip][_C_][sh] = vShape[sh];
806 : }
807 : }
808 0 : }
809 :
810 : // computes the linearized defect w.r.t to the velocity
811 : template<typename TDomain>
812 : template <typename TElem, typename TFVGeom>
813 0 : void ConvectionDiffusionFVCR<TDomain>::
814 : ex_grad(MathVector<dim> vValue[],
815 : const MathVector<dim> vGlobIP[],
816 : number time, int si,
817 : const LocalVector& u,
818 : GridObject* elem,
819 : const MathVector<dim> vCornerCoords[],
820 : const MathVector<TFVGeom::dim> vLocIP[],
821 : const size_t nip,
822 : bool bDeriv,
823 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[])
824 : {
825 : // Get finite volume geometry
826 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
827 :
828 : // reference element
829 : typedef typename reference_element_traits<TElem>::reference_element_type
830 : ref_elem_type;
831 :
832 : // reference dimension
833 : static const int refDim = ref_elem_type::dim;
834 :
835 : // number of shape functions
836 : static const size_t numSH = ref_elem_type::numCorners;
837 :
838 : // CRFV SCVF ip
839 0 : if(vLocIP == geo.scvf_local_ips())
840 : {
841 : // Loop Sub Control Volume Faces (SCVF)
842 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
843 : {
844 : // Get current SCVF
845 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
846 :
847 0 : VecSet(vValue[ip], 0.0);
848 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
849 : VecScaleAppend(vValue[ip], u(_C_, sh), scvf.global_grad(sh));
850 :
851 0 : if(bDeriv)
852 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
853 0 : vvvDeriv[ip][_C_][sh] = scvf.global_grad(sh);
854 : }
855 : }
856 : // general case
857 : else
858 : {
859 : // get trial space
860 0 : LagrangeP1<ref_elem_type>& rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
861 :
862 : // storage for shape function at ip
863 0 : MathVector<refDim> vLocGrad[numSH];
864 : MathVector<refDim> locGrad;
865 :
866 : // Reference Mapping
867 : MathMatrix<dim, refDim> JTInv;
868 0 : ReferenceMapping<ref_elem_type, dim> mapping(vCornerCoords);
869 :
870 : // loop ips
871 0 : for(size_t ip = 0; ip < nip; ++ip)
872 : {
873 : // evaluate at shapes at ip
874 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
875 :
876 : // compute grad at ip
877 : VecSet(locGrad, 0.0);
878 0 : for(size_t sh = 0; sh < numSH; ++sh)
879 : VecScaleAppend(locGrad, u(_C_, sh), vLocGrad[sh]);
880 :
881 : // compute global grad
882 0 : mapping.jacobian_transposed_inverse(JTInv, vLocIP[ip]);
883 0 : MatVecMult(vValue[ip], JTInv, locGrad);
884 :
885 : // compute derivative w.r.t. to unknowns iff needed
886 0 : if(bDeriv)
887 0 : for(size_t sh = 0; sh < numSH; ++sh)
888 0 : MatVecMult(vvvDeriv[ip][_C_][sh], JTInv, vLocGrad[sh]);
889 : }
890 : }
891 0 : };
892 :
893 : ////////////////////////////////////////////////////////////////////////////////
894 : // upwind
895 : ////////////////////////////////////////////////////////////////////////////////
896 :
897 : template<typename TDomain>
898 0 : void ConvectionDiffusionFVCR<TDomain>::
899 0 : set_upwind(SmartPtr<IConvectionShapes<dim> > shapes) {m_spConvShape = shapes;}
900 :
901 : // computes the linearized defect w.r.t to the velocity
902 : template<typename TDomain>
903 : const typename ConvectionDiffusionFVCR<TDomain>::conv_shape_type&
904 0 : ConvectionDiffusionFVCR<TDomain>::
905 : get_updated_conv_shapes(const FVGeometryBase& geo)
906 : {
907 : // compute upwind shapes for transport equation
908 : // \todo: we should move this computation into the preparation part of the
909 : // disc, to only compute the shapes once, reusing them several times.
910 0 : if(m_imVelocity.data_given())
911 : {
912 : // get diffusion at ips
913 : const MathMatrix<dim, dim>* vDiffusion = NULL;
914 0 : if(m_imDiffusion.data_given()) vDiffusion = m_imDiffusion.values();
915 :
916 : // update convection shapes
917 0 : if(!m_spConvShape->update(&geo, m_imVelocity.values(), vDiffusion, true))
918 : {
919 : UG_LOG("ERROR in 'ConvectionDiffusionFV1::add_jac_A_elem': "
920 : "Cannot compute convection shapes.\n");
921 : }
922 : }
923 :
924 : // return a const (!!) reference to the upwind
925 0 : return *const_cast<const IConvectionShapes<dim>*>(m_spConvShape.get());
926 : }
927 :
928 :
929 : ////////////////////////////////////////////////////////////////////////////////
930 : // register assemble functions
931 : ////////////////////////////////////////////////////////////////////////////////
932 :
933 : #ifdef UG_DIM_1
934 : template<>
935 0 : void ConvectionDiffusionFVCR<Domain1d>::
936 : register_all_funcs(bool bHang)
937 : {
938 0 : UG_THROW("Crouxeiz-Raviart only senseful in dimension >= 2");
939 : }
940 : #endif
941 :
942 : #ifdef UG_DIM_2
943 : template<>
944 0 : void ConvectionDiffusionFVCR<Domain2d>::
945 : register_all_funcs(bool bHang)
946 : {
947 : // switch assemble functions
948 0 : if(!bHang)
949 : {
950 0 : register_func<Triangle, CRFVGeometry<Triangle, dim> >();
951 0 : register_func<Quadrilateral, CRFVGeometry<Quadrilateral, dim> >();
952 : }
953 : else
954 : {
955 0 : register_func<Triangle, HCRFVGeometry<Triangle, dim> >();
956 0 : register_func<Quadrilateral, HCRFVGeometry<Quadrilateral, dim> >();
957 : }
958 0 : }
959 : #endif
960 :
961 : #ifdef UG_DIM_3
962 : template<>
963 0 : void ConvectionDiffusionFVCR<Domain3d>::
964 : register_all_funcs(bool bHang)
965 : {
966 : // switch assemble functions
967 0 : if(!bHang)
968 : {
969 0 : register_func<Tetrahedron, CRFVGeometry<Tetrahedron, dim> >();
970 0 : register_func<Prism, CRFVGeometry<Prism, dim> >();
971 0 : register_func<Pyramid, CRFVGeometry<Pyramid, dim> >();
972 0 : register_func<Hexahedron, CRFVGeometry<Hexahedron, dim> >();
973 : }
974 : else
975 : {
976 0 : register_func<Tetrahedron, HCRFVGeometry<Tetrahedron, dim> >();
977 0 : register_func<Prism, HCRFVGeometry<Prism, dim> >();
978 0 : register_func<Pyramid, HCRFVGeometry<Pyramid, dim> >();
979 0 : register_func<Hexahedron, HCRFVGeometry<Hexahedron, dim> >();
980 : }
981 0 : }
982 : #endif
983 :
984 : template<typename TDomain>
985 : template<typename TElem, typename TFVGeom>
986 0 : void ConvectionDiffusionFVCR<TDomain>::
987 : register_func()
988 : {
989 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
990 : typedef this_type T;
991 : static const int refDim = reference_element_traits<TElem>::dim;
992 :
993 : this->clear_add_fct(id);
994 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFVGeom>);
995 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFVGeom>);
996 : this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop<TElem, TFVGeom>);
997 : this->set_add_jac_A_elem_fct(id, &T::template add_jac_A_elem<TElem, TFVGeom>);
998 : this->set_add_jac_M_elem_fct(id, &T::template add_jac_M_elem<TElem, TFVGeom>);
999 : this->set_add_def_A_elem_fct(id, &T::template add_def_A_elem<TElem, TFVGeom>);
1000 : this->set_add_def_M_elem_fct(id, &T::template add_def_M_elem<TElem, TFVGeom>);
1001 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFVGeom>);
1002 :
1003 : // set computation of linearized defect w.r.t velocity
1004 0 : m_imVelocity. set_fct(id, this, &T::template lin_def_velocity<TElem, TFVGeom>);
1005 0 : m_imDiffusion.set_fct(id, this, &T::template lin_def_diffusion<TElem, TFVGeom>);
1006 0 : m_imReactionRate. set_fct(id, this, &T::template lin_def_reaction_rate<TElem, TFVGeom>);
1007 0 : m_imReaction. set_fct(id, this, &T::template lin_def_reaction<TElem, TFVGeom>);
1008 0 : m_imSource. set_fct(id, this, &T::template lin_def_source<TElem, TFVGeom>);
1009 0 : m_imMassScale.set_fct(id, this, &T::template lin_def_mass_scale<TElem, TFVGeom>);
1010 0 : m_imMass. set_fct(id, this, &T::template lin_def_mass<TElem, TFVGeom>);
1011 :
1012 : // exports
1013 0 : m_exValue-> template set_fct<T,refDim>(id, this, &T::template ex_value<TElem, TFVGeom>);
1014 0 : m_exGrad->template set_fct<T,refDim>(id, this, &T::template ex_grad<TElem, TFVGeom>);
1015 0 : }
1016 :
1017 : ////////////////////////////////////////////////////////////////////////////////
1018 : // explicit template instantiations
1019 : ////////////////////////////////////////////////////////////////////////////////
1020 :
1021 : #ifdef UG_DIM_1
1022 : template class ConvectionDiffusionFVCR<Domain1d>;
1023 : #endif
1024 : #ifdef UG_DIM_2
1025 : template class ConvectionDiffusionFVCR<Domain2d>;
1026 : #endif
1027 : #ifdef UG_DIM_3
1028 : template class ConvectionDiffusionFVCR<Domain3d>;
1029 : #endif
1030 :
1031 : } // end namespace ConvectionDiffusionPlugin
1032 : } // namespace ug
1033 :
|
