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_fv1.h"
34 :
35 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
36 : #include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
37 : #include "lib_disc/spatial_disc/disc_util/hfv1_geom.h"
38 : #include "lib_disc/spatial_disc/disc_util/conv_shape.h"
39 :
40 : namespace ug{
41 : namespace ConvectionDiffusionPlugin{
42 :
43 : DebugID DID_CONV_DIFF_FV1("CONV_DIFF_FV1");
44 :
45 : ////////////////////////////////////////////////////////////////////////////////
46 : // general
47 : ////////////////////////////////////////////////////////////////////////////////
48 :
49 : template<typename TDomain>
50 0 : ConvectionDiffusionFV1<TDomain>::
51 : ConvectionDiffusionFV1(const char* functions, const char* subsets)
52 : : ConvectionDiffusionBase<TDomain>(functions,subsets),
53 0 : m_spConvShape(new ConvectionShapesNoUpwind<dim>),
54 0 : m_bNonRegularGrid(false), m_bCondensedFV(false)
55 : {
56 : register_all_funcs(m_bNonRegularGrid);
57 0 : }
58 :
59 : template<typename TDomain>
60 0 : void ConvectionDiffusionFV1<TDomain>::
61 : prepare_setting(const std::vector<LFEID>& vLfeID, bool bNonRegularGrid)
62 : {
63 : // check number
64 0 : if(vLfeID.size() != 1)
65 0 : UG_THROW("ConvectionDiffusion: Wrong number of functions given. "
66 : "Need exactly "<<1);
67 :
68 0 : if(vLfeID[0].order() != 1 || vLfeID[0].type() != LFEID::LAGRANGE)
69 0 : UG_THROW("ConvectionDiffusion FV Scheme only implemented for 1st order.");
70 :
71 : // remember
72 0 : m_bNonRegularGrid = bNonRegularGrid;
73 :
74 : // update assemble functions
75 : register_all_funcs(m_bNonRegularGrid);
76 0 : }
77 :
78 : template<typename TDomain>
79 0 : bool ConvectionDiffusionFV1<TDomain>::
80 : use_hanging() const
81 : {
82 0 : return true;
83 : }
84 :
85 : ////////////////////////////////////////////////////////////////////////////////
86 : // Assembling functions
87 : ////////////////////////////////////////////////////////////////////////////////
88 :
89 : template<typename TDomain>
90 0 : void ConvectionDiffusionFV1<TDomain>::
91 : prep_assemble_loop()
92 : {
93 0 : if (m_sss_mngr.valid ())
94 : // reset the markers of the point sources and sinks (as there are no marks for the lines)
95 : m_sss_mngr->init_all_point_sss ();
96 0 : }
97 :
98 : template<typename TDomain>
99 : template<typename TElem, typename TFVGeom>
100 0 : void ConvectionDiffusionFV1<TDomain>::
101 : prep_elem_loop(const ReferenceObjectID roid, const int si)
102 : {
103 :
104 : // check, that upwind has been set
105 : //TODO: It is really used only if we have the convection
106 0 : if(m_spConvShape.invalid())
107 0 : UG_THROW("ConvectionDiffusionFV1::prep_elem_loop:"
108 : " Upwind has not been set.");
109 :
110 : // set local positions
111 : if(!TFVGeom::usesHangingNodes)
112 : {
113 : static const int refDim = TElem::dim;
114 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get();
115 : const MathVector<refDim>* vSCVFip = geo.scvf_local_ips();
116 : const size_t numSCVFip = geo.num_scvf_ips();
117 : const MathVector<refDim>* vSCVip = geo.scv_local_ips();
118 : const size_t numSCVip = geo.num_scv_ips();
119 0 : m_imDiffusion.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
120 0 : m_imVelocity.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
121 0 : m_imFlux.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
122 0 : m_imSource.template set_local_ips<refDim>(vSCVip,numSCVip, false);
123 0 : m_imVectorSource.template set_local_ips<refDim>(vSCVFip,numSCVFip, false);
124 0 : m_imReactionRate.template set_local_ips<refDim>(vSCVip,numSCVip, false);
125 0 : m_imReaction.template set_local_ips<refDim>(vSCVip,numSCVip, false);
126 0 : m_imReactionRateExpl.template set_local_ips<refDim>(vSCVip,numSCVip, false);
127 0 : m_imReactionExpl.template set_local_ips<refDim>(vSCVip,numSCVip, false);
128 0 : m_imSourceExpl.template set_local_ips<refDim>(vSCVip,numSCVip, false);
129 0 : m_imMassScale.template set_local_ips<refDim>(vSCVip,numSCVip, false);
130 0 : m_imMass.template set_local_ips<refDim>(vSCVip,numSCVip, false);
131 :
132 : // init upwind for element type
133 0 : if(!m_spConvShape->template set_geometry_type<TFVGeom>(geo))
134 0 : UG_THROW("ConvectionDiffusionFV1::prep_elem_loop:"
135 : " Cannot init upwind for element type.");
136 : }
137 0 : }
138 :
139 : template<typename TDomain>
140 : template<typename TElem, typename TFVGeom>
141 0 : void ConvectionDiffusionFV1<TDomain>::
142 : fsh_elem_loop()
143 0 : {}
144 :
145 : template<typename TDomain>
146 : template<typename TElem, typename TFVGeom>
147 0 : void ConvectionDiffusionFV1<TDomain>::
148 : prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, const MathVector<dim> vCornerCoords[])
149 : {
150 : // Update Geometry for this element
151 0 : static TFVGeom& geo = GeomProvider<TFVGeom>::get();
152 :
153 : try{
154 : UG_DLOG(DID_CONV_DIFF_FV1, 2, ">>OCT_DISC_DEBUG: " << "convection_diffusion_fv1.cpp: " << "prep_elem(): update(): "<< roid << std::endl);
155 0 : geo.update(elem, vCornerCoords, &(this->subset_handler()));
156 0 : }UG_CATCH_THROW("ConvectionDiffusionFV1::prep_elem:"
157 : " Cannot update Finite Volume Geometry.");
158 :
159 : // set local positions
160 : if(TFVGeom::usesHangingNodes)
161 : {
162 : const int refDim = TElem::dim;
163 0 : const MathVector<refDim>* vSCVFip = geo.scvf_local_ips();
164 : const size_t numSCVFip = geo.num_scvf_ips();
165 : const MathVector<refDim>* vSCVip = geo.scv_local_ips();
166 : const size_t numSCVip = geo.num_scv_ips();
167 0 : m_imDiffusion.template set_local_ips<refDim>(vSCVFip,numSCVFip);
168 0 : m_imVelocity.template set_local_ips<refDim>(vSCVFip,numSCVFip);
169 0 : m_imFlux.template set_local_ips<refDim>(vSCVFip,numSCVFip);
170 0 : m_imSource.template set_local_ips<refDim>(vSCVip,numSCVip);
171 0 : m_imVectorSource.template set_local_ips<refDim>(vSCVFip,numSCVFip);
172 0 : m_imReactionRate.template set_local_ips<refDim>(vSCVip,numSCVip);
173 0 : m_imReaction.template set_local_ips<refDim>(vSCVip,numSCVip);
174 0 : m_imReactionRateExpl.template set_local_ips<refDim>(vSCVip,numSCVip);
175 0 : m_imReactionExpl.template set_local_ips<refDim>(vSCVip,numSCVip);
176 0 : m_imSourceExpl.template set_local_ips<refDim>(vSCVip,numSCVip);
177 0 : m_imMassScale.template set_local_ips<refDim>(vSCVip,numSCVip);
178 0 : m_imMass.template set_local_ips<refDim>(vSCVip,numSCVip);
179 :
180 0 : if(m_spConvShape.valid())
181 0 : if(!m_spConvShape->template set_geometry_type<TFVGeom>(geo))
182 0 : UG_THROW("ConvectionDiffusionFV1::prep_elem:"
183 : " Cannot init upwind for element type.");
184 : }
185 :
186 : // set global positions
187 0 : const MathVector<dim>* vSCVFip = geo.scvf_global_ips();
188 : const size_t numSCVFip = geo.num_scvf_ips();
189 : const MathVector<dim>* vSCVip = geo.scv_global_ips();
190 : const size_t numSCVip = geo.num_scv_ips();
191 0 : m_imDiffusion. set_global_ips(vSCVFip, numSCVFip);
192 0 : m_imVelocity. set_global_ips(vSCVFip, numSCVFip);
193 0 : m_imFlux. set_global_ips(vSCVFip, numSCVFip);
194 0 : m_imSource. set_global_ips(vSCVip, numSCVip);
195 0 : m_imVectorSource. set_global_ips(vSCVFip, numSCVFip);
196 0 : m_imReactionRate. set_global_ips(vSCVip, numSCVip);
197 0 : m_imReactionRateExpl. set_global_ips(vSCVip, numSCVip);
198 0 : m_imReactionExpl. set_global_ips(vSCVip, numSCVip);
199 0 : m_imSourceExpl. set_global_ips(vSCVip, numSCVip);
200 0 : m_imReaction. set_global_ips(vSCVip, numSCVip);
201 0 : m_imMassScale. set_global_ips(vSCVip, numSCVip);
202 0 : m_imMass. set_global_ips(vSCVip, numSCVip);
203 0 : }
204 :
205 : template <class TVector>
206 : static TVector CalculateCenter(GridObject* o, const TVector* coords)
207 : {
208 : TVector v;
209 : VecSet(v, 0);
210 :
211 : size_t numCoords = 0;
212 : switch(o->base_object_id()){
213 : case VERTEX: numCoords = 1; break;
214 : case EDGE: numCoords = static_cast<Edge*>(o)->num_vertices(); break;
215 : case FACE: numCoords = static_cast<Face*>(o)->num_vertices(); break;
216 : case VOLUME: numCoords = static_cast<Volume*>(o)->num_vertices(); break;
217 : default: UG_THROW("Unknown element type."); break;
218 : }
219 :
220 : for(size_t i = 0; i < numCoords; ++i)
221 : VecAdd(v, v, coords[i]);
222 :
223 : if(numCoords > 0)
224 : VecScale(v, v, 1. / (number)numCoords);
225 :
226 : return v;
227 : }
228 :
229 : template<typename TDomain>
230 : template<typename TElem, typename TFVGeom>
231 : void ConvectionDiffusionFV1<TDomain>::
232 : add_sss_jac_elem
233 : (
234 : LocalMatrix& J, ///< the matrix to update
235 : const LocalVector& u, ///< current solution
236 : GridObject* elem, ///< the element
237 : const TFVGeom& geo, ///< the FV geometry for that element
238 : size_t co, ///< index of the SCV
239 : number flux ///< flux through source/sink (premultiplied by the length for lines)
240 : )
241 : {
242 :
243 0 : if (flux < 0.0)
244 : // sink
245 0 : J(_C_, co, _C_, co) -= flux;
246 : }
247 :
248 : template<typename TDomain>
249 : template<typename TElem, typename TFVGeom>
250 0 : void ConvectionDiffusionFV1<TDomain>::
251 : add_jac_A_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
252 : {
253 : // get finite volume geometry
254 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
255 :
256 : // Diff. Tensor times Gradient
257 : MathVector<dim> Dgrad;
258 :
259 : // get conv shapes
260 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo, false);
261 :
262 : // Diffusion and Velocity Term
263 0 : if(m_imDiffusion.data_given() || m_imVelocity.data_given())
264 : {
265 : // loop Sub Control Volume Faces (SCVF)
266 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
267 : {
268 : // get current SCVF
269 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
270 :
271 : ////////////////////////////////////////////////////
272 : // Diffusive Term
273 : ////////////////////////////////////////////////////
274 0 : if(m_imDiffusion.data_given())
275 : {
276 : #ifdef UG_ENABLE_DEBUG_LOGS
277 : // DID_CONV_DIFF_FV1
278 : number D_diff_flux_sum = 0.0;
279 : #endif
280 :
281 : // loop shape functions
282 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
283 : {
284 : // Compute Diffusion Tensor times Gradient
285 : MatVecMult(Dgrad, m_imDiffusion[ip], scvf.global_grad(sh));
286 :
287 : // Compute flux at IP
288 : const number D_diff_flux = VecDot(Dgrad, scvf.normal());
289 : UG_DLOG(DID_CONV_DIFF_FV1, 2, ">>OCT_DISC_DEBUG: " << "convection_diffusion_fv1.cpp: " << "add_jac_A_elem(): " << "sh # " << sh << " ; normalSize scvf # " << ip << ": " << VecLength(scvf.normal()) << "; \t from "<< scvf.from() << "; to " << scvf.to() << "; D_diff_flux: " << D_diff_flux << "; scvf.global_grad(sh): " << scvf.global_grad(sh) << std::endl);
290 :
291 : // Add flux term to local matrix // HIER MATRIXINDIZES!!!
292 : UG_ASSERT((scvf.from() < J.num_row_dof(_C_)) && (scvf.to() < J.num_col_dof(_C_)),
293 : "Bad local dof-index on element with object-id " << elem->base_object_id()
294 : << " with center: " << CalculateCenter(elem, vCornerCoords));
295 :
296 0 : J(_C_, scvf.from(), _C_, sh) -= D_diff_flux;
297 0 : J(_C_, scvf.to() , _C_, sh) += D_diff_flux;
298 :
299 : #ifdef UG_ENABLE_DEBUG_LOGS
300 : // DID_CONV_DIFF_FV1
301 : D_diff_flux_sum += D_diff_flux;
302 : #endif
303 : }
304 :
305 : UG_DLOG(DID_CONV_DIFF_FV1, 2, "D_diff_flux_sum = " << D_diff_flux_sum << std::endl << std::endl);
306 : }
307 :
308 : ////////////////////////////////////////////////////
309 : // Convective Term
310 : ////////////////////////////////////////////////////
311 0 : if(m_imVelocity.data_given())
312 : {
313 : // Add Flux contribution
314 0 : for(size_t sh = 0; sh < convShape.num_sh(); ++sh)
315 : {
316 : const number D_conv_flux = convShape(ip, sh);
317 :
318 : // Add flux term to local matrix
319 0 : J(_C_, scvf.from(), _C_, sh) += D_conv_flux;
320 0 : J(_C_, scvf.to(), _C_, sh) -= D_conv_flux;
321 : }
322 : }
323 :
324 : // no explicit dependency on flux import
325 : }
326 : }
327 :
328 : //UG_LOG("Local Matrix is: \n"<<J<<"\n");
329 :
330 : ////////////////////////////////////////////////////
331 : // Reaction Term (using lumping)
332 : ////////////////////////////////////////////////////
333 :
334 0 : if(m_imReactionRate.data_given())
335 : {
336 : // loop Sub Control Volume (SCV)
337 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
338 : {
339 : // get current SCV
340 : const typename TFVGeom::SCV& scv = geo.scv(ip);
341 :
342 : // get associated node
343 0 : const int co = scv.node_id();
344 :
345 : // Add to local matrix
346 0 : J(_C_, co, _C_, co) += m_imReactionRate[ip] * scv.volume();
347 : }
348 : }
349 :
350 : // reaction term does not explicitly depend on the associated unknown function
351 :
352 : ////////////////////////////////
353 : // Singular sources and sinks
354 : ////////////////////////////////
355 :
356 0 : if (m_sss_mngr.valid () && (m_sss_mngr->num_points () != 0 || m_sss_mngr->num_lines () != 0))
357 : {
358 : typedef typename TDomain::position_accessor_type t_pos_accessor;
359 : typedef typename CDSingularSourcesAndSinks<dim>::template
360 : point_iterator<TElem,t_pos_accessor,TFVGeom> t_pnt_sss_iter;
361 : typedef typename CDSingularSourcesAndSinks<dim>::template
362 : line_iterator<TElem,t_pos_accessor,TFVGeom> t_lin_sss_iter;
363 :
364 0 : t_pos_accessor& aaPos = this->domain()->position_accessor();
365 0 : Grid& grid = (Grid&) *this->domain()->grid();
366 :
367 0 : for(size_t co = 0; co < geo.num_sh(); co++)
368 : {
369 :
370 : // point sources
371 0 : for (t_pnt_sss_iter pnt (m_sss_mngr.get(), (TElem *) elem, grid, aaPos, geo, co);
372 0 : ! pnt.is_over (); ++pnt)
373 : {
374 : FVPointSourceOrSink<dim, cd_point_sss_data<dim> > * pnt_sss = *pnt;
375 0 : if (! pnt_sss->marked_for (elem, co))
376 0 : continue;
377 : pnt_sss->compute (pnt_sss->position (), this->time (), -1); //TODO: set the subset id instead of -1
378 : this->template add_sss_jac_elem<TElem, TFVGeom> (J, u, elem, geo, co,
379 : pnt_sss->flux ());
380 : }
381 :
382 : // line sources
383 0 : for (t_lin_sss_iter line (m_sss_mngr.get(), (TElem *) elem, grid, aaPos, geo, co);
384 0 : ! line.is_over (); ++line)
385 : {
386 : FVLineSourceOrSink<dim, cd_line_sss_data<dim> > * line_sss = *line;
387 0 : number len = VecDistance (line.seg_start (), line.seg_end ());
388 : line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
389 0 : this->template add_sss_jac_elem<TElem, TFVGeom> (J, u, elem, geo, co,
390 : line_sss->flux () * len);
391 : }
392 : }
393 : }
394 0 : }
395 :
396 :
397 : template<typename TDomain>
398 : template<typename TElem, typename TFVGeom>
399 0 : void ConvectionDiffusionFV1<TDomain>::
400 : add_jac_M_elem(LocalMatrix& J, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
401 : {
402 : // get finite volume geometry
403 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
404 :
405 0 : if(!m_imMassScale.data_given()) return;
406 :
407 : // loop Sub Control Volumes (SCV)
408 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
409 : {
410 : // get current SCV
411 : const typename TFVGeom::SCV& scv = geo.scv(ip);
412 :
413 : // get associated node
414 0 : const int co = scv.node_id();
415 :
416 : // Add to local matrix
417 0 : J(_C_, co, _C_, co) += scv.volume() * m_imMassScale[ip];
418 : }
419 :
420 : // m_imMass part does not explicitly depend on associated unknown function
421 : }
422 :
423 :
424 : template<typename TDomain>
425 : template<typename TElem, typename TFVGeom>
426 : void ConvectionDiffusionFV1<TDomain>::
427 : add_sss_def_elem
428 : (
429 : LocalVector& d, ///< the defect to update
430 : const LocalVector& u, ///< current solution
431 : GridObject* elem, ///< the element
432 : const TFVGeom& geo, ///< the FV geometry for that element
433 : size_t co, ///< index of the SCV
434 : number flux ///< flux through source/sink (premultiplied by the length for lines)
435 : )
436 : {
437 0 : if (flux > 0.0)
438 : // source
439 0 : d(_C_, co) -= flux;
440 : else
441 : // sink
442 0 : d(_C_, co) -= flux * u(_C_, co);
443 : }
444 :
445 : template<typename TDomain>
446 : template<typename TElem, typename TFVGeom>
447 0 : void ConvectionDiffusionFV1<TDomain>::
448 : add_def_A_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
449 : {
450 : // get finite volume geometry
451 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
452 :
453 : // get conv shapes
454 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo, false);
455 :
456 0 : if(m_imDiffusion.data_given() || m_imVelocity.data_given() || m_imFlux.data_given())
457 : {
458 : // loop Sub Control Volume Faces (SCVF)
459 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
460 : {
461 : // get current SCVF
462 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
463 :
464 : /////////////////////////////////////////////////////
465 : // Diffusive Term
466 : /////////////////////////////////////////////////////
467 0 : if(m_imDiffusion.data_given())
468 : {
469 : // to compute D \nabla c
470 : MathVector<dim> Dgrad_c, grad_c;
471 :
472 : // compute gradient and shape at ip
473 : VecSet(grad_c, 0.0);
474 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
475 : VecScaleAppend(grad_c, u(_C_,sh), scvf.global_grad(sh));
476 :
477 : // scale by diffusion tensor
478 : MatVecMult(Dgrad_c, m_imDiffusion[ip], grad_c);
479 :
480 : // Compute flux
481 : const number diff_flux = VecDot(Dgrad_c, scvf.normal());
482 :
483 : // Add to local defect
484 0 : d(_C_, scvf.from()) -= diff_flux;
485 0 : d(_C_, scvf.to() ) += diff_flux;
486 : }
487 :
488 : /////////////////////////////////////////////////////
489 : // Convective Term
490 : /////////////////////////////////////////////////////
491 0 : if(m_imVelocity.data_given())
492 : {
493 : // sum up convective flux using convection shapes
494 : number conv_flux = 0.0;
495 0 : for(size_t sh = 0; sh < convShape.num_sh(); ++sh)
496 0 : conv_flux += u(_C_, sh) * convShape(ip, sh);
497 :
498 : // add to local defect
499 0 : d(_C_, scvf.from()) += conv_flux;
500 0 : d(_C_, scvf.to() ) -= conv_flux;
501 : }
502 :
503 : /////////////////////////////////////////////////////
504 : // Flux Term
505 : /////////////////////////////////////////////////////
506 0 : if(m_imFlux.data_given())
507 : {
508 : // sum up flux
509 : const number flux = VecDot(m_imFlux[ip], scvf.normal());
510 :
511 : // add to local defect
512 0 : d(_C_, scvf.from()) += flux;
513 0 : d(_C_, scvf.to() ) -= flux;
514 : }
515 :
516 : }
517 : }
518 :
519 : // reaction rate
520 0 : if(m_imReactionRate.data_given())
521 : {
522 : // loop Sub Control Volumes (SCV)
523 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
524 : {
525 : // get current SCV
526 : const typename TFVGeom::SCV& scv = geo.scv(ip);
527 :
528 : // get associated node
529 0 : const int co = scv.node_id();
530 :
531 : // Add to local defect
532 0 : d(_C_, co) += u(_C_, co) * m_imReactionRate[ip] * scv.volume();
533 : }
534 : }
535 :
536 : // reaction term
537 0 : if(m_imReaction.data_given())
538 : {
539 : // loop Sub Control Volumes (SCV)
540 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
541 : {
542 : // get current SCV
543 : const typename TFVGeom::SCV& scv = geo.scv(ip);
544 :
545 : // get associated node
546 0 : const int co = scv.node_id();
547 :
548 : // Add to local defect
549 0 : d(_C_, co) += m_imReaction[ip] * scv.volume();
550 : }
551 : }
552 :
553 : ////////////////////////////////
554 : // Singular sources and sinks
555 : ////////////////////////////////
556 :
557 0 : if (m_sss_mngr.valid () && (m_sss_mngr->num_points () != 0 || m_sss_mngr->num_lines () != 0))
558 : {
559 : typedef typename TDomain::position_accessor_type t_pos_accessor;
560 : typedef typename CDSingularSourcesAndSinks<dim>::template
561 : point_iterator<TElem,t_pos_accessor,TFVGeom> t_pnt_sss_iter;
562 : typedef typename CDSingularSourcesAndSinks<dim>::template
563 : line_iterator<TElem,t_pos_accessor,TFVGeom> t_lin_sss_iter;
564 :
565 0 : t_pos_accessor& aaPos = this->domain()->position_accessor();
566 0 : Grid& grid = (Grid&) *this->domain()->grid();
567 :
568 : //for(size_t i = 0; i < geo.num_scv(); i++)
569 0 : for(size_t co = 0; co < geo.num_sh(); co++)
570 : {
571 :
572 : // point sources
573 0 : for (t_pnt_sss_iter pnt (m_sss_mngr.get(), (TElem *) elem, grid, aaPos, geo, co);
574 0 : ! pnt.is_over (); ++pnt)
575 : {
576 : FVPointSourceOrSink<dim, cd_point_sss_data<dim> > * pnt_sss = *pnt;
577 0 : if (! pnt_sss->marked_for (elem, co))
578 0 : continue;
579 : pnt_sss->compute (pnt_sss->position (), this->time (), -1); //TODO: set the subset id instead of -1
580 : this->template add_sss_def_elem<TElem, TFVGeom> (d, u, elem, geo, co,
581 : pnt_sss->flux ());
582 : }
583 :
584 : // line sources
585 0 : for (t_lin_sss_iter line (m_sss_mngr.get(), (TElem *) elem, grid, aaPos, geo, co);
586 0 : ! line.is_over (); ++line)
587 : {
588 : FVLineSourceOrSink<dim, cd_line_sss_data<dim> > * line_sss = *line;
589 0 : number len = VecDistance (line.seg_start (), line.seg_end ());
590 : line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
591 0 : this->template add_sss_def_elem<TElem, TFVGeom> (d, u, elem, geo, co,
592 : line_sss->flux () * len);
593 : }
594 : }
595 : }
596 0 : }
597 :
598 : template<typename TDomain>
599 : template<typename TElem, typename TFVGeom>
600 0 : void ConvectionDiffusionFV1<TDomain>::
601 : add_def_A_expl_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
602 : {
603 : // get finite volume geometry
604 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
605 :
606 : // reaction rate
607 0 : if(m_imReactionRateExpl.data_given())
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 : const typename TFVGeom::SCV& scv = geo.scv(ip);
614 :
615 : // get associated node
616 0 : const int co = scv.node_id();
617 :
618 : // Add to local defect
619 0 : d(_C_, co) += u(_C_, co) * m_imReactionRateExpl[ip] * scv.volume();
620 : }
621 : }
622 :
623 : // reaction
624 0 : if(m_imReactionExpl.data_given())
625 : {
626 : // loop Sub Control Volumes (SCV)
627 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
628 : {
629 : // get current SCV
630 : const typename TFVGeom::SCV& scv = geo.scv(ip);
631 :
632 : // get associated node
633 0 : const int co = scv.node_id();
634 :
635 : // Add to local defect
636 0 : d(_C_, co) += m_imReactionExpl[ip] * scv.volume();
637 : }
638 : }
639 :
640 0 : if(m_imSourceExpl.data_given())
641 : {
642 : // loop Sub Control Volumes (SCV)
643 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
644 : {
645 : // get current SCV
646 : const typename TFVGeom::SCV& scv = geo.scv(ip);
647 :
648 : // get associated node
649 0 : const int co = scv.node_id();
650 :
651 : // Add to local rhs
652 0 : d(_C_, co) -= m_imSourceExpl[ip] * scv.volume();
653 : }
654 : }
655 0 : }
656 :
657 : template<typename TDomain>
658 : template<typename TElem, typename TFVGeom>
659 0 : void ConvectionDiffusionFV1<TDomain>::
660 : add_def_M_elem(LocalVector& d, const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
661 : {
662 : // get finite volume geometry
663 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
664 :
665 0 : if(!m_imMassScale.data_given() && !m_imMass.data_given()) return;
666 :
667 : // loop Sub Control Volumes (SCV)
668 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
669 : {
670 : // get current SCV
671 : const typename TFVGeom::SCV& scv = geo.scv(ip);
672 :
673 : // get associated node
674 0 : const int co = scv.node_id();
675 :
676 : // mass value
677 : number val = 0.0;
678 :
679 : // multiply by scaling
680 0 : if(m_imMassScale.data_given())
681 0 : val += m_imMassScale[ip] * u(_C_, co);
682 :
683 : // add mass
684 0 : if(m_imMass.data_given())
685 0 : val += m_imMass[ip];
686 :
687 : // Add to local defect
688 0 : d(_C_, co) += val * scv.volume();
689 : }
690 : }
691 :
692 :
693 : template<typename TDomain>
694 : template<typename TElem, typename TFVGeom>
695 0 : void ConvectionDiffusionFV1<TDomain>::
696 : add_rhs_elem(LocalVector& d, GridObject* elem, const MathVector<dim> vCornerCoords[])
697 : {
698 : // get finite volume geometry
699 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
700 :
701 : // loop Sub Control Volumes (SCV)
702 0 : if ( m_imSource.data_given() ) {
703 0 : for ( size_t ip = 0; ip < geo.num_scv(); ++ip ) {
704 : // get current SCV
705 : const typename TFVGeom::SCV& scv = geo.scv( ip );
706 :
707 : // get associated node
708 0 : const int co = scv.node_id();
709 :
710 : // Add to local rhs
711 0 : d(_C_, co) += m_imSource[ip] * scv.volume();
712 : //UG_LOG("d(_C_, co) = " << d(_C_, co) << "; \t ip " << ip << "; \t co " << co << "; \t scv_vol " << scv.volume() << "; \t m_imSource[ip] " << m_imSource[ip] << std::endl);
713 : }
714 : }
715 :
716 : // loop Sub Control Volumes (SCVF)
717 0 : if ( m_imVectorSource.data_given() ) {
718 0 : for ( size_t ip = 0; ip < geo.num_scvf(); ++ip ) {
719 : // get current SCVF
720 : const typename TFVGeom::SCVF& scvf = geo.scvf( ip );
721 :
722 : // Add to local rhs
723 : number flux = VecDot(m_imVectorSource[ip], scvf.normal());
724 0 : d(_C_, scvf.from()) -= flux;
725 0 : d(_C_, scvf.to() ) += flux;
726 : }
727 : }
728 0 : }
729 :
730 :
731 : // ////////////////////////////////
732 : // error estimation (begin) ///
733 :
734 : // prepares the loop over all elements of one type for the computation of the error estimator
735 : template<typename TDomain>
736 : template<typename TElem, typename TFVGeom>
737 0 : void ConvectionDiffusionFV1<TDomain>::
738 : prep_err_est_elem_loop(const ReferenceObjectID roid, const int si)
739 : {
740 : // get the error estimator data object and check that it is of the right type
741 : // we check this at this point in order to be able to dispense with this check later on
742 : // (i.e. in prep_err_est_elem and compute_err_est_A_elem())
743 0 : if (this->m_spErrEstData.get() == NULL)
744 : {
745 0 : UG_THROW("No ErrEstData object has been given to this ElemDisc!");
746 : }
747 :
748 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
749 :
750 0 : if (!err_est_data)
751 : {
752 0 : UG_THROW("Dynamic cast to SideAndElemErrEstData failed."
753 : << std::endl << "Make sure you handed the correct type of ErrEstData to this discretization.");
754 : }
755 :
756 :
757 : // check that upwind has been set
758 0 : if (m_spConvShape.invalid())
759 0 : UG_THROW("ConvectionDiffusionFV1::prep_err_est_elem_loop: "
760 : "Upwind has not been set.");
761 :
762 : // set local positions
763 : if (!TFVGeom::usesHangingNodes)
764 : {
765 : static const int refDim = TElem::dim;
766 :
767 : // get local IPs
768 : size_t numSideIPs, numElemIPs;
769 : const MathVector<refDim>* sideIPs;
770 : const MathVector<refDim>* elemIPs;
771 : try
772 : {
773 : numSideIPs = err_est_data->num_all_side_ips(roid);
774 0 : numElemIPs = err_est_data->num_elem_ips(roid);
775 0 : sideIPs = err_est_data->template side_local_ips<refDim>(roid);
776 0 : elemIPs = err_est_data->template elem_local_ips<refDim>(roid);
777 :
778 0 : if (!sideIPs || !elemIPs) return; // are NULL if TElem is not of the same dim as TDomain
779 : }
780 0 : UG_CATCH_THROW("Integration points for error estimator cannot be set.");
781 :
782 : // set local IPs in imports
783 0 : m_imDiffusion.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
784 0 : m_imVelocity.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
785 0 : m_imFlux.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
786 0 : m_imSource.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
787 0 : m_imVectorSource.template set_local_ips<refDim>(sideIPs, numSideIPs, false);
788 0 : m_imReactionRate.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
789 0 : m_imReaction.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
790 0 : m_imMassScale.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
791 0 : m_imMass.template set_local_ips<refDim>(elemIPs, numElemIPs, false);
792 :
793 : // init upwind for element type
794 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get();
795 0 : if (!m_spConvShape->template set_geometry_type<TFVGeom>(geo))
796 0 : UG_THROW("ConvectionDiffusionFV1::prep_err_est_elem_loop: "
797 : "Cannot init upwind for element type.");
798 :
799 : // store values of shape functions in local IPs
800 : LagrangeP1<typename reference_element_traits<TElem>::reference_element_type> trialSpace
801 0 : = Provider<LagrangeP1<typename reference_element_traits<TElem>::reference_element_type> >::get();
802 :
803 0 : m_shapeValues.resize(numElemIPs, numSideIPs, trialSpace.num_sh());
804 0 : for (size_t ip = 0; ip < numElemIPs; ip++)
805 0 : trialSpace.shapes(m_shapeValues.shapesAtElemIP(ip), elemIPs[ip]);
806 0 : for (size_t ip = 0; ip < numSideIPs; ip++)
807 0 : trialSpace.shapes(m_shapeValues.shapesAtSideIP(ip), sideIPs[ip]);
808 : }
809 : }
810 :
811 : template<typename TDomain>
812 : template<typename TElem, typename TFVGeom>
813 0 : void ConvectionDiffusionFV1<TDomain>::
814 : prep_err_est_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[])
815 : {
816 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
817 :
818 : // update geometry for this element
819 0 : static TFVGeom& geo = GeomProvider<TFVGeom>::get();
820 : try
821 : {
822 0 : geo.update(elem, vCornerCoords, &(this->subset_handler()));
823 : }
824 0 : UG_CATCH_THROW("ConvectionDiffusionFV1::prep_err_est_elem: Cannot update Finite Volume Geometry.");
825 :
826 : // roid
827 0 : ReferenceObjectID roid = elem->reference_object_id();
828 :
829 : // set local positions
830 : if (TFVGeom::usesHangingNodes)
831 : {
832 : static const int refDim = TElem::dim;
833 :
834 : size_t numSideIPs, numElemIPs;
835 : const MathVector<refDim>* sideIPs;
836 : const MathVector<refDim>* elemIPs;
837 : try
838 : {
839 : numSideIPs = err_est_data->num_all_side_ips(roid);
840 0 : numElemIPs = err_est_data->num_elem_ips(roid);
841 0 : sideIPs = err_est_data->template side_local_ips<refDim>(roid);
842 0 : elemIPs = err_est_data->template elem_local_ips<refDim>(roid);
843 :
844 0 : if (!sideIPs || !elemIPs) return; // are NULL if TElem is not of the same dim as TDomain
845 : }
846 0 : UG_CATCH_THROW("Integration points for error estimator cannot be set.");
847 :
848 0 : m_imDiffusion.template set_local_ips<refDim>(sideIPs, numSideIPs);
849 0 : m_imVelocity.template set_local_ips<refDim>(sideIPs, numSideIPs);
850 0 : m_imFlux.template set_local_ips<refDim>(sideIPs, numSideIPs);
851 0 : m_imSource.template set_local_ips<refDim>(elemIPs, numElemIPs);
852 0 : m_imVectorSource.template set_local_ips<refDim>(sideIPs, numSideIPs);
853 0 : m_imReactionRate.template set_local_ips<refDim>(elemIPs, numElemIPs);
854 0 : m_imReaction.template set_local_ips<refDim>(elemIPs, numElemIPs);
855 0 : m_imMassScale.template set_local_ips<refDim>(elemIPs, numElemIPs);
856 0 : m_imMass.template set_local_ips<refDim>(elemIPs, numElemIPs);
857 :
858 : // init upwind for element type
859 0 : TFVGeom& geo = GeomProvider<TFVGeom>::get();
860 0 : if (!m_spConvShape->template set_geometry_type<TFVGeom>(geo))
861 0 : UG_THROW("ConvectionDiffusionFV1::prep_err_est_elem_loop: "
862 : "Cannot init upwind for element type.");
863 :
864 : // store values of shape functions in local IPs
865 : LagrangeP1<typename reference_element_traits<TElem>::reference_element_type> trialSpace
866 0 : = Provider<LagrangeP1<typename reference_element_traits<TElem>::reference_element_type> >::get();
867 :
868 0 : m_shapeValues.resize(numElemIPs, numSideIPs, trialSpace.num_sh());
869 0 : for (size_t ip = 0; ip < numElemIPs; ip++)
870 0 : trialSpace.shapes(m_shapeValues.shapesAtElemIP(ip), elemIPs[ip]);
871 0 : for (size_t ip = 0; ip < numSideIPs; ip++)
872 0 : trialSpace.shapes(m_shapeValues.shapesAtSideIP(ip), sideIPs[ip]);
873 : }
874 :
875 : // set global positions
876 : size_t numSideIPs, numElemIPs;
877 : const MathVector<dim>* sideIPs;
878 : const MathVector<dim>* elemIPs;
879 :
880 : try
881 : {
882 : numSideIPs = err_est_data->num_all_side_ips(roid);
883 0 : numElemIPs = err_est_data->num_elem_ips(roid);
884 :
885 0 : sideIPs = err_est_data->all_side_global_ips(elem, vCornerCoords);
886 0 : elemIPs = err_est_data->elem_global_ips(elem, vCornerCoords);
887 : }
888 0 : UG_CATCH_THROW("Global integration points for error estimator cannot be set.");
889 :
890 0 : m_imDiffusion. set_global_ips(&sideIPs[0], numSideIPs);
891 0 : m_imVelocity. set_global_ips(&sideIPs[0], numSideIPs);
892 0 : m_imFlux. set_global_ips(&sideIPs[0], numSideIPs);
893 0 : m_imSource. set_global_ips(&elemIPs[0], numElemIPs);
894 0 : m_imVectorSource. set_global_ips(&sideIPs[0], numSideIPs);
895 0 : m_imReactionRate. set_global_ips(&elemIPs[0], numElemIPs);
896 0 : m_imReaction. set_global_ips(&elemIPs[0], numElemIPs);
897 0 : m_imMassScale. set_global_ips(&elemIPs[0], numElemIPs);
898 0 : m_imMass. set_global_ips(&elemIPs[0], numElemIPs);
899 : }
900 :
901 : // computes the error estimator contribution (stiffness part) for one element
902 : template<typename TDomain>
903 : template<typename TElem, typename TFVGeom>
904 0 : void ConvectionDiffusionFV1<TDomain>::
905 : compute_err_est_A_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
906 : {
907 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
908 :
909 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
910 :
911 0 : if (err_est_data->surface_view().get() == NULL) {UG_THROW("Error estimator has NULL surface view.");}
912 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
913 :
914 : // request geometry
915 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
916 :
917 :
918 : // SIDE TERMS //
919 :
920 : // get the sides of the element
921 : // We have to cast elem to a pointer of type SideAndElemErrEstData::elem_type
922 : // for the SideAndElemErrEstData::operator() to work properly.
923 : // This cannot generally be achieved by casting to TElem*, since this method is also registered for
924 : // lower-dimensional types TElem, and must therefore be compilable, even if it is never EVER to be executed.
925 : // The way we achieve this here, is by calling associated_elements_sorted() which has an implementation for
926 : // all possible types. Whatever comes out of it is of course complete nonsense if (and only if)
927 : // SideAndElemErrEstData::elem_type != TElem. To be on the safe side, we throw an error if the number of
928 : // entries in the list is not as it should be.
929 :
930 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::side_type>::secure_container side_list;
931 0 : pErrEstGrid->associated_elements_sorted(side_list, (TElem*) elem);
932 0 : if (side_list.size() != (size_t) ref_elem_type::numSides)
933 0 : UG_THROW ("Mismatch of numbers of sides in 'ConvectionDiffusionFV1::compute_err_est_elem'");
934 :
935 : // some help variables
936 : MathVector<dim> fluxDensity, gradC, normal;
937 :
938 : // FIXME: The computation of the gradient has to be reworked.
939 : // In the case of P1 shape functions, it is valid. For Q1 shape functions, however,
940 : // the gradient is not constant (but bilinear) on the element - and along the sides.
941 : // We cannot use the FVGeom here. Instead, we need to calculate the gradient in each IP!
942 :
943 : // calculate grad u (take grad from first scvf ip (grad u is constant on the entire element))
944 : /* if (geo.num_scvf() < 1) {UG_THROW("Element has no SCVFs!");}
945 : const typename TFVGeom::SCVF& scvf = geo.scvf(0);
946 :
947 : VecSet(gradC, 0.0);
948 : for (size_t j=0; j<m_shapeValues.num_sh(); j++)
949 : VecScaleAppend(gradC, u(_C_,j), scvf.global_grad(j));*/
950 :
951 : // calculate grad u as average (over scvf)
952 : VecSet(gradC, 0.0);
953 0 : for(size_t ii = 0; ii < geo.num_scvf(); ++ii)
954 : {
955 : const typename TFVGeom::SCVF& scvf = geo.scvf(ii);
956 0 : for (size_t j=0; j<m_shapeValues.num_sh(); j++)
957 : VecScaleAppend(gradC, u(_C_,j), scvf.global_grad(j));
958 : }
959 0 : VecScale(gradC, gradC, (1.0/geo.num_scvf()));
960 :
961 : /*VecSet(gradC, 0.0);
962 : for(size_t ii = 0; ii < geo.num_scv(); ++ii)
963 : {
964 : const typename TFVGeom::SCV& scv = geo.scv(ii);
965 : for (size_t j=0; j<m_shapeValues.num_sh(); j++)
966 : VecScaleAppend(gradC, u(_C_,j), scv.global_grad(j));
967 : }
968 : VecScale(gradC, gradC, (1.0/geo.num_scvf()));
969 : */
970 :
971 : // calculate flux through the sides
972 : size_t passedIPs = 0;
973 0 : for (size_t side=0; side < (size_t) ref_elem_type::numSides; side++)
974 : {
975 : // normal on side
976 0 : SideNormal<ref_elem_type,dim>(normal, side, vCornerCoords);
977 0 : VecNormalize(normal, normal);
978 :
979 : try
980 : {
981 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
982 : {
983 0 : size_t ip = passedIPs + sip;
984 :
985 : VecSet(fluxDensity, 0.0);
986 :
987 : // diffusion //
988 0 : if (m_imDiffusion.data_given())
989 : MatVecScaleMultAppend(fluxDensity, -1.0, m_imDiffusion[ip], gradC);
990 :
991 : // convection //
992 0 : if (m_imVelocity.data_given())
993 : {
994 : number val = 0.0;
995 0 : for (size_t sh = 0; sh < m_shapeValues.num_sh(); sh++)
996 0 : val += u(_C_,sh) * m_shapeValues.shapeAtSideIP(sh,sip);
997 :
998 : VecScaleAppend(fluxDensity, val, m_imVelocity[ip]);
999 : }
1000 :
1001 : // general flux //
1002 0 : if (m_imFlux.data_given())
1003 : VecAppend(fluxDensity, m_imFlux[ip]);
1004 :
1005 0 : (*err_est_data)(side_list[side],sip) += scale * VecDot(fluxDensity, normal);
1006 : }
1007 :
1008 0 : passedIPs += err_est_data->num_side_ips(side_list[side]);
1009 : }
1010 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
1011 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
1012 : }
1013 :
1014 : // VOLUME TERMS //
1015 :
1016 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::elem_type>::secure_container elem_list;
1017 : pErrEstGrid->associated_elements_sorted(elem_list, (TElem*) elem);
1018 0 : if (elem_list.size() != 1)
1019 0 : UG_THROW ("Mismatch of numbers of sides in 'ConvectionDiffusionFV1::compute_err_est_elem'");
1020 :
1021 : try
1022 : {
1023 0 : for (size_t ip = 0; ip < err_est_data->num_elem_ips(elem->reference_object_id()); ip++)
1024 : {
1025 : number total = 0.0;
1026 :
1027 : // diffusion // TODO ONLY FOR (PIECEWISE) CONSTANT DIFFUSION TENSOR SO FAR!
1028 : // div(D*grad(c))
1029 : // nothing to do, as c is piecewise linear and div(D*grad(c)) disappears
1030 : // if D is diagonal and c bilinear, this should also vanish (confirm this!)
1031 :
1032 : // convection // TODO ONLY FOR (PIECEWISE) CONSTANT OR DIVERGENCE-FREE
1033 : // VELOCITY FIELDS SO FAR!
1034 : // div(v*c) = div(v)*c + v*grad(c) -- gradC has been calculated above
1035 0 : if (m_imVelocity.data_given())
1036 0 : total += VecDot(m_imVelocity[ip], gradC);
1037 :
1038 : // general flux // TODO ONLY FOR DIVERGENCE-FREE FLUX FIELD SO FAR!
1039 : // nothing to do
1040 :
1041 : // reaction //
1042 0 : if (m_imReactionRate.data_given())
1043 : {
1044 : number val = 0.0;
1045 0 : for (size_t sh = 0; sh < geo.num_sh(); sh++)
1046 0 : val += u(_C_,sh) * m_shapeValues.shapeAtElemIP(sh,ip);
1047 :
1048 0 : total += m_imReactionRate[ip] * val;
1049 : }
1050 :
1051 0 : if (m_imReaction.data_given())
1052 : {
1053 0 : total += m_imReaction[ip];
1054 : }
1055 :
1056 0 : (*err_est_data)(elem_list[0],ip) += scale * total;
1057 : }
1058 : }
1059 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
1060 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
1061 0 : }
1062 :
1063 : // computes the error estimator contribution (mass part) for one element
1064 : template<typename TDomain>
1065 : template<typename TElem, typename TFVGeom>
1066 0 : void ConvectionDiffusionFV1<TDomain>::
1067 : compute_err_est_M_elem(const LocalVector& u, GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
1068 : {
1069 : // note: mass parts only enter volume term
1070 :
1071 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
1072 :
1073 0 : if (err_est_data->surface_view().get() == NULL) {UG_THROW("Error estimator has NULL surface view.");}
1074 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
1075 :
1076 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::elem_type>::secure_container elem_list;
1077 0 : pErrEstGrid->associated_elements_sorted(elem_list, (TElem*) elem);
1078 0 : if (elem_list.size() != 1)
1079 0 : UG_THROW ("Mismatch of numbers of sides in 'ConvectionDiffusionFV1::compute_err_est_elem'");
1080 :
1081 : // request geometry
1082 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1083 :
1084 : // loop integration points
1085 : try
1086 : {
1087 0 : for (size_t ip = 0; ip < err_est_data->num_elem_ips(elem->reference_object_id()); ip++)
1088 : {
1089 : number total = 0.0;
1090 :
1091 : // mass scale //
1092 0 : if (m_imMassScale.data_given())
1093 : {
1094 : number val = 0.0;
1095 0 : for (size_t sh = 0; sh < geo.num_sh(); sh++)
1096 0 : val += u(_C_,sh) * m_shapeValues.shapeAtElemIP(sh,ip);
1097 :
1098 0 : total += m_imMassScale[ip] * val;
1099 : }
1100 :
1101 : // mass //
1102 0 : if (m_imMass.data_given())
1103 : {
1104 0 : total += m_imMass[ip];
1105 : }
1106 :
1107 0 : (*err_est_data)(elem_list[0],ip) += scale * total;
1108 : }
1109 : }
1110 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
1111 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
1112 0 : }
1113 :
1114 : // computes the error estimator contribution (rhs part) for one element
1115 : template<typename TDomain>
1116 : template<typename TElem, typename TFVGeom>
1117 0 : void ConvectionDiffusionFV1<TDomain>::
1118 : compute_err_est_rhs_elem(GridObject* elem, const MathVector<dim> vCornerCoords[], const number& scale)
1119 : {
1120 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
1121 :
1122 0 : err_est_type* err_est_data = dynamic_cast<err_est_type*>(this->m_spErrEstData.get());
1123 :
1124 0 : if (err_est_data->surface_view().get() == NULL) {UG_THROW("Error estimator has NULL surface view.");}
1125 0 : MultiGrid* pErrEstGrid = (MultiGrid*) (err_est_data->surface_view()->subset_handler()->multi_grid());
1126 :
1127 : // SIDE TERMS //
1128 :
1129 : // get the sides of the element
1130 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::side_type>::secure_container side_list;
1131 0 : pErrEstGrid->associated_elements_sorted(side_list, (TElem*) elem);
1132 0 : if (side_list.size() != (size_t) ref_elem_type::numSides)
1133 0 : UG_THROW ("Mismatch of numbers of sides in 'ConvectionDiffusionFV1::compute_err_est_elem'");
1134 :
1135 : // loop sides
1136 : size_t passedIPs = 0;
1137 0 : for (size_t side = 0; side < (size_t) ref_elem_type::numSides; side++)
1138 : {
1139 : // normal on side
1140 : MathVector<dim> normal;
1141 0 : SideNormal<ref_elem_type,dim>(normal, side, vCornerCoords);
1142 0 : VecNormalize(normal, normal);
1143 :
1144 : try
1145 : {
1146 0 : for (size_t sip = 0; sip < err_est_data->num_side_ips(side_list[side]); sip++)
1147 : {
1148 0 : size_t ip = passedIPs + sip;
1149 :
1150 : // vector source //
1151 0 : if (m_imVectorSource.data_given())
1152 0 : (*err_est_data)(side_list[side],sip) += scale * VecDot(m_imVectorSource[ip], normal);
1153 : }
1154 :
1155 0 : passedIPs += err_est_data->num_side_ips(side_list[side]);
1156 : }
1157 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
1158 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
1159 : }
1160 :
1161 : // VOLUME TERMS //
1162 :
1163 0 : if (!m_imSource.data_given()) return;
1164 :
1165 : typename MultiGrid::traits<typename SideAndElemErrEstData<TDomain>::elem_type>::secure_container elem_list;
1166 : pErrEstGrid->associated_elements_sorted(elem_list, (TElem*) elem);
1167 0 : if (elem_list.size() != 1)
1168 0 : UG_THROW ("Mismatch of numbers of sides in 'ConvectionDiffusionFV1::compute_err_est_elem'");
1169 :
1170 : // source //
1171 : try
1172 : {
1173 0 : for (size_t ip = 0; ip < err_est_data->num_elem_ips(elem->reference_object_id()); ip++)
1174 0 : (*err_est_data)(elem_list[0],ip) += scale * m_imSource[ip];
1175 : }
1176 0 : UG_CATCH_THROW("Values for the error estimator could not be assembled at every IP." << std::endl
1177 : << "Maybe wrong type of ErrEstData object? This implementation needs: SideAndElemErrEstData.");
1178 : }
1179 :
1180 : // postprocesses the loop over all elements of one type in the computation of the error estimator
1181 : template<typename TDomain>
1182 : template<typename TElem, typename TFVGeom>
1183 0 : void ConvectionDiffusionFV1<TDomain>::
1184 : fsh_err_est_elem_loop()
1185 : {
1186 : // finish the element loop in the same way as the actual discretization
1187 : this->template fsh_elem_loop<TElem, TFVGeom> ();
1188 0 : };
1189 :
1190 : // error estimation (end) ///
1191 : // /////////////////////////////////
1192 :
1193 : // computes the linearized defect w.r.t to the velocity
1194 : template<typename TDomain>
1195 : template <typename TElem, typename TFVGeom>
1196 0 : void ConvectionDiffusionFV1<TDomain>::
1197 : lin_def_velocity(const LocalVector& u,
1198 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
1199 : const size_t nip)
1200 : {
1201 : // get finite volume geometry
1202 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1203 :
1204 : // get conv shapes
1205 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo, true);
1206 :
1207 : // reset the values for the linearized defect
1208 0 : for(size_t ip = 0; ip < nip; ++ip)
1209 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1210 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1211 : vvvLinDef[ip][c][sh] = 0.0;
1212 :
1213 : // loop Sub Control Volume Faces (SCVF)
1214 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
1215 : {
1216 : // get current SCVF
1217 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
1218 :
1219 : // sum up contributions of convection shapes
1220 : MathVector<dim> linDefect;
1221 : VecSet(linDefect, 0.0);
1222 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1223 : VecScaleAppend(linDefect, u(_C_,sh), convShape.D_vel(ip, sh));
1224 :
1225 : // add parts for both sides of scvf
1226 0 : vvvLinDef[ip][_C_][scvf.from()] += linDefect;
1227 : vvvLinDef[ip][_C_][scvf.to()] -= linDefect;
1228 : }
1229 0 : }
1230 :
1231 : // computes the linearized defect w.r.t to the velocity
1232 : template<typename TDomain>
1233 : template <typename TElem, typename TFVGeom>
1234 0 : void ConvectionDiffusionFV1<TDomain>::
1235 : lin_def_diffusion(const LocalVector& u,
1236 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
1237 : const size_t nip)
1238 : {
1239 : // get finite volume geometry
1240 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1241 :
1242 : // get conv shapes
1243 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes(geo, true);
1244 :
1245 : // reset the values for the linearized defect
1246 0 : for(size_t ip = 0; ip < nip; ++ip)
1247 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1248 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1249 : vvvLinDef[ip][c][sh] = 0.0;
1250 :
1251 : // loop Sub Control Volume Faces (SCVF)
1252 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
1253 : {
1254 : // get current SCVF
1255 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
1256 :
1257 : // compute gradient at ip
1258 : MathVector<dim> grad_u; VecSet(grad_u, 0.0);
1259 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1260 : VecScaleAppend(grad_u, u(_C_,sh), scvf.global_grad(sh));
1261 :
1262 : // compute the lin defect at this ip
1263 : MathMatrix<dim,dim> linDefect;
1264 :
1265 : // part coming from -\nabla u * \vec{n}
1266 0 : for(size_t k=0; k < (size_t)dim; ++k)
1267 0 : for(size_t j = 0; j < (size_t)dim; ++j)
1268 0 : linDefect(j,k) = (scvf.normal())[j] * grad_u[k];
1269 :
1270 : // add contribution from convection shapes
1271 0 : if(convShape.non_zero_deriv_diffusion())
1272 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1273 : MatAdd(linDefect, convShape.D_diffusion(ip, sh), u(_C_, sh));
1274 :
1275 : // add contributions
1276 0 : vvvLinDef[ip][_C_][scvf.from()] -= linDefect;
1277 : vvvLinDef[ip][_C_][scvf.to() ] += linDefect;
1278 : }
1279 0 : }
1280 :
1281 : template<typename TDomain>
1282 : template <typename TElem, typename TFVGeom>
1283 0 : void ConvectionDiffusionFV1<TDomain>::
1284 : lin_def_flux(const LocalVector& u,
1285 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
1286 : const size_t nip)
1287 : {
1288 : // get finite volume geometry
1289 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1290 :
1291 : // reset the values for the linearized defect
1292 0 : for(size_t ip = 0; ip < nip; ++ip)
1293 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1294 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1295 : vvvLinDef[ip][c][sh] = 0.0;
1296 :
1297 : // loop Sub Control Volume Faces (SCVF)
1298 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
1299 : {
1300 : // get current SCVF
1301 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
1302 :
1303 : // add parts for both sides of scvf
1304 0 : vvvLinDef[ip][_C_][scvf.from()] += scvf.normal();
1305 : vvvLinDef[ip][_C_][scvf.to()] -= scvf.normal();
1306 : }
1307 0 : }
1308 :
1309 : // computes the linearized defect w.r.t to the reaction rate
1310 : template<typename TDomain>
1311 : template <typename TElem, typename TFVGeom>
1312 0 : void ConvectionDiffusionFV1<TDomain>::
1313 : lin_def_reaction_rate(const LocalVector& u,
1314 : std::vector<std::vector<number> > vvvLinDef[],
1315 : const size_t nip)
1316 : {
1317 : // get finite volume geometry
1318 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1319 :
1320 : // loop Sub Control Volumes (SCV)
1321 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
1322 : {
1323 : // get current SCV
1324 : const typename TFVGeom::SCV& scv = geo.scv(ip);
1325 :
1326 : // get associated node
1327 0 : const int co = scv.node_id();
1328 :
1329 : // set lin defect
1330 0 : vvvLinDef[ip][_C_][co] = u(_C_, co) * scv.volume();
1331 : }
1332 0 : }
1333 :
1334 : // computes the linearized defect w.r.t to the reaction
1335 : template<typename TDomain>
1336 : template <typename TElem, typename TFVGeom>
1337 0 : void ConvectionDiffusionFV1<TDomain>::
1338 : lin_def_reaction(const LocalVector& u,
1339 : std::vector<std::vector<number> > vvvLinDef[],
1340 : const size_t nip)
1341 : {
1342 : // get finite volume geometry
1343 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1344 :
1345 : // loop Sub Control Volumes (SCV)
1346 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
1347 : {
1348 : // get current SCV
1349 : const typename TFVGeom::SCV& scv = geo.scv(ip);
1350 :
1351 : // get associated node
1352 0 : const int co = scv.node_id();
1353 :
1354 : // set lin defect
1355 0 : vvvLinDef[ip][_C_][co] = scv.volume();
1356 : }
1357 0 : }
1358 :
1359 : // computes the linearized defect w.r.t to the source
1360 : template<typename TDomain>
1361 : template <typename TElem, typename TFVGeom>
1362 0 : void ConvectionDiffusionFV1<TDomain>::
1363 : lin_def_source(const LocalVector& u,
1364 : std::vector<std::vector<number> > vvvLinDef[],
1365 : const size_t nip)
1366 : {
1367 : // get finite volume geometry
1368 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1369 :
1370 : // loop Sub Control Volumes (SCV)
1371 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
1372 : {
1373 : // get current SCV
1374 : const typename TFVGeom::SCV& scv = geo.scv(ip);
1375 :
1376 : // get associated node
1377 0 : const int co = scv.node_id();
1378 :
1379 : // set lin defect
1380 0 : vvvLinDef[ip][_C_][co] = scv.volume();
1381 : }
1382 0 : }
1383 :
1384 : // computes the linearized defect w.r.t to the vector source
1385 : // (in analogy to velocity)
1386 : template<typename TDomain>
1387 : template <typename TElem, typename TFVGeom>
1388 0 : void ConvectionDiffusionFV1<TDomain>::
1389 : lin_def_vector_source(const LocalVector& u,
1390 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
1391 : const size_t nip)
1392 : {
1393 : // get finite volume geometry
1394 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1395 :
1396 : // reset the values for the linearized defect
1397 0 : for(size_t ip = 0; ip < nip; ++ip)
1398 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1399 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1400 : vvvLinDef[ip][c][sh] = 0.0;
1401 :
1402 : // loop Sub Control Volumes Faces (SCVF)
1403 0 : for ( size_t ip = 0; ip < geo.num_scvf(); ++ip ) {
1404 : // get current SCVF
1405 : const typename TFVGeom::SCVF& scvf = geo.scvf( ip );
1406 :
1407 : // add parts for both sides of scvf
1408 0 : vvvLinDef[ip][_C_][scvf.from()] -= scvf.normal();
1409 : vvvLinDef[ip][_C_][scvf.to()] += scvf.normal();
1410 : }
1411 0 : }
1412 :
1413 : // computes the linearized defect w.r.t to the mass scale
1414 : template<typename TDomain>
1415 : template <typename TElem, typename TFVGeom>
1416 0 : void ConvectionDiffusionFV1<TDomain>::
1417 : lin_def_mass_scale(const LocalVector& u,
1418 : std::vector<std::vector<number> > vvvLinDef[],
1419 : const size_t nip)
1420 : {
1421 : // get finite volume geometry
1422 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1423 :
1424 : // loop Sub Control Volumes (SCV)
1425 0 : for(size_t co = 0; co < geo.num_scv(); ++co)
1426 : {
1427 : // get current SCV
1428 : const typename TFVGeom::SCV& scv = geo.scv(co);
1429 :
1430 : // Check associated node
1431 : UG_ASSERT(co == scv.node_id(), "Only one shape per SCV");
1432 :
1433 : // set lin defect
1434 0 : vvvLinDef[co][_C_][co] = u(_C_, co) * scv.volume();
1435 : }
1436 0 : }
1437 :
1438 : // computes the linearized defect w.r.t to the mass scale
1439 : template<typename TDomain>
1440 : template <typename TElem, typename TFVGeom>
1441 0 : void ConvectionDiffusionFV1<TDomain>::
1442 : lin_def_mass(const LocalVector& u,
1443 : std::vector<std::vector<number> > vvvLinDef[],
1444 : const size_t nip)
1445 : {
1446 : // get finite volume geometry
1447 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1448 :
1449 : // loop Sub Control Volumes (SCV)
1450 0 : for(size_t co = 0; co < geo.num_scv(); ++co)
1451 : {
1452 : // get current SCV
1453 : const typename TFVGeom::SCV& scv = geo.scv(co);
1454 :
1455 : // Check associated node
1456 : UG_ASSERT(co == scv.node_id(), "Only one shape per SCV");
1457 :
1458 : // set lin defect
1459 0 : vvvLinDef[co][_C_][co] = scv.volume();
1460 : }
1461 0 : }
1462 :
1463 : // computes the linearized defect w.r.t to the velocity
1464 : template<typename TDomain>
1465 : template <typename TElem, typename TFVGeom>
1466 0 : void ConvectionDiffusionFV1<TDomain>::
1467 : ex_value(number vValue[],
1468 : const MathVector<dim> vGlobIP[],
1469 : number time, int si,
1470 : const LocalVector& u,
1471 : GridObject* elem,
1472 : const MathVector<dim> vCornerCoords[],
1473 : const MathVector<TFVGeom::dim> vLocIP[],
1474 : const size_t nip,
1475 : bool bDeriv,
1476 : std::vector<std::vector<number> > vvvDeriv[])
1477 : {
1478 : // get finite volume geometry
1479 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1480 :
1481 : // reference element
1482 : typedef typename reference_element_traits<TElem>::reference_element_type
1483 : ref_elem_type;
1484 :
1485 : // number of shape functions
1486 : static const size_t numSH = ref_elem_type::numCorners;
1487 :
1488 :
1489 : // FV1 SCVF ip
1490 0 : if(vLocIP == geo.scvf_local_ips())
1491 : {
1492 : // Loop Sub Control Volume Faces (SCVF)
1493 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
1494 : {
1495 : // Get current SCVF
1496 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
1497 :
1498 : // compute concentration at ip
1499 0 : vValue[ip] = 0.0;
1500 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1501 0 : vValue[ip] += u(_C_, sh) * scvf.shape(sh);
1502 :
1503 : // compute derivative w.r.t. to unknowns iff needed
1504 0 : if(bDeriv)
1505 : {
1506 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1507 0 : vvvDeriv[ip][_C_][sh] = scvf.shape(sh);
1508 :
1509 : // do not forget that number of DoFs (== vvvDeriv[ip][_C_])
1510 : // might be > scvf.num_sh() in case of hanging nodes!
1511 0 : size_t ndof = vvvDeriv[ip][_C_].size();
1512 0 : for (size_t sh = scvf.num_sh(); sh < ndof; ++sh)
1513 0 : vvvDeriv[ip][_C_][sh] = 0.0;
1514 : }
1515 : }
1516 : }
1517 : // FV1 SCV ip
1518 0 : else if(vLocIP == geo.scv_local_ips())
1519 : {
1520 : // Loop Sub Control Volumes (SCV)
1521 0 : for(size_t ip = 0; ip < geo.num_scv(); ++ip)
1522 : {
1523 : // Get current SCV
1524 : const typename TFVGeom::SCV& scv = geo.scv(ip);
1525 :
1526 : // get corner of SCV
1527 : const size_t co = scv.node_id();
1528 :
1529 : // solution at ip
1530 0 : vValue[ip] = u(_C_, co);
1531 :
1532 : // set derivatives if needed
1533 0 : if(bDeriv)
1534 : {
1535 0 : size_t ndof = vvvDeriv[ip][_C_].size();
1536 0 : for(size_t sh = 0; sh < ndof; ++sh)
1537 0 : vvvDeriv[ip][_C_][sh] = (sh==co) ? 1.0 : 0.0;
1538 : }
1539 : }
1540 : }
1541 : // general case
1542 : else
1543 : {
1544 : // get trial space
1545 0 : LagrangeP1<ref_elem_type>& rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
1546 :
1547 : // storage for shape function at ip
1548 : number vShape[numSH];
1549 :
1550 : // loop ips
1551 0 : for(size_t ip = 0; ip < nip; ++ip)
1552 : {
1553 : // evaluate at shapes at ip
1554 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
1555 :
1556 : // compute concentration at ip
1557 0 : vValue[ip] = 0.0;
1558 0 : for(size_t sh = 0; sh < numSH; ++sh)
1559 0 : vValue[ip] += u(_C_, sh) * vShape[sh];
1560 :
1561 : // compute derivative w.r.t. to unknowns iff needed
1562 : // \todo: maybe store shapes directly in vvvDeriv
1563 0 : if(bDeriv)
1564 : {
1565 0 : for(size_t sh = 0; sh < numSH; ++sh)
1566 0 : vvvDeriv[ip][_C_][sh] = vShape[sh];
1567 :
1568 : // beware of hanging nodes!
1569 0 : size_t ndof = vvvDeriv[ip][_C_].size();
1570 0 : for (size_t sh = numSH; sh < ndof; ++sh)
1571 0 : vvvDeriv[ip][_C_][sh] = 0.0;
1572 : }
1573 : }
1574 : }
1575 0 : }
1576 :
1577 : // computes the linearized defect w.r.t to the velocity
1578 : template<typename TDomain>
1579 : template <typename TElem, typename TFVGeom>
1580 0 : void ConvectionDiffusionFV1<TDomain>::
1581 : ex_grad(MathVector<dim> vValue[],
1582 : const MathVector<dim> vGlobIP[],
1583 : number time, int si,
1584 : const LocalVector& u,
1585 : GridObject* elem,
1586 : const MathVector<dim> vCornerCoords[],
1587 : const MathVector<TFVGeom::dim> vLocIP[],
1588 : const size_t nip,
1589 : bool bDeriv,
1590 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[])
1591 : {
1592 : // Get finite volume geometry
1593 0 : static const TFVGeom& geo = GeomProvider<TFVGeom>::get();
1594 :
1595 : // reference element
1596 : typedef typename reference_element_traits<TElem>::reference_element_type
1597 : ref_elem_type;
1598 :
1599 : // reference dimension
1600 : static const int refDim = ref_elem_type::dim;
1601 :
1602 : // number of shape functions
1603 : static const size_t numSH = ref_elem_type::numCorners;
1604 :
1605 :
1606 : // FV1 SCVF ip
1607 0 : if(vLocIP == geo.scvf_local_ips())
1608 : {
1609 : // Loop Sub Control Volume Faces (SCVF)
1610 0 : for(size_t ip = 0; ip < geo.num_scvf(); ++ip)
1611 : {
1612 : // Get current SCVF
1613 : const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
1614 :
1615 0 : VecSet(vValue[ip], 0.0);
1616 :
1617 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1618 : VecScaleAppend(vValue[ip], u(_C_, sh), scvf.global_grad(sh));
1619 :
1620 0 : if(bDeriv)
1621 : {
1622 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1623 0 : vvvDeriv[ip][_C_][sh] = scvf.global_grad(sh);
1624 :
1625 : // beware of hanging nodes!
1626 0 : size_t ndof = vvvDeriv[ip][_C_].size();
1627 0 : for (size_t sh = scvf.num_sh(); sh < ndof; ++sh)
1628 : vvvDeriv[ip][_C_][sh] = 0.0;
1629 : }
1630 : }
1631 : }
1632 : // general case
1633 : else
1634 : {
1635 : // get trial space
1636 0 : LagrangeP1<ref_elem_type>& rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
1637 :
1638 : // storage for shape function at ip
1639 0 : MathVector<refDim> vLocGrad[numSH];
1640 : MathVector<refDim> locGrad;
1641 :
1642 : // Reference Mapping
1643 : MathMatrix<dim, refDim> JTInv;
1644 0 : ReferenceMapping<ref_elem_type, dim> mapping(vCornerCoords);
1645 :
1646 : // loop ips
1647 0 : for(size_t ip = 0; ip < nip; ++ip)
1648 : {
1649 : // evaluate at shapes at ip
1650 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
1651 :
1652 : // compute grad at ip
1653 : VecSet(locGrad, 0.0);
1654 0 : for(size_t sh = 0; sh < numSH; ++sh)
1655 0 : VecScaleAppend(locGrad, u(_C_, sh), vLocGrad[sh]);
1656 :
1657 : // compute global grad
1658 0 : mapping.jacobian_transposed_inverse(JTInv, vLocIP[ip]);
1659 0 : MatVecMult(vValue[ip], JTInv, locGrad);
1660 :
1661 : // compute derivative w.r.t. to unknowns iff needed
1662 0 : if(bDeriv)
1663 : {
1664 0 : for(size_t sh = 0; sh < numSH; ++sh)
1665 0 : MatVecMult(vvvDeriv[ip][_C_][sh], JTInv, vLocGrad[sh]);
1666 :
1667 : // beware of hanging nodes!
1668 0 : size_t ndof = vvvDeriv[ip][_C_].size();
1669 0 : for (size_t sh = numSH; sh < ndof; ++sh)
1670 : vvvDeriv[ip][_C_][sh] = 0.0;
1671 : }
1672 : }
1673 : }
1674 0 : };
1675 :
1676 : ////////////////////////////////////////////////////////////////////////////////
1677 : // upwind
1678 : ////////////////////////////////////////////////////////////////////////////////
1679 :
1680 : template<typename TDomain>
1681 0 : void ConvectionDiffusionFV1<TDomain>::
1682 0 : set_upwind(SmartPtr<IConvectionShapes<dim> > shapes) {m_spConvShape = shapes;}
1683 :
1684 : // computes the linearized defect w.r.t to the velocity
1685 : template<typename TDomain>
1686 : const typename ConvectionDiffusionFV1<TDomain>::conv_shape_type&
1687 0 : ConvectionDiffusionFV1<TDomain>::
1688 : get_updated_conv_shapes(const FVGeometryBase& geo, bool compute_deriv)
1689 : {
1690 : // compute upwind shapes for transport equation
1691 : // \todo: we should move this computation into the preparation part of the
1692 : // disc, to only compute the shapes once, reusing them several times.
1693 0 : if(m_imVelocity.data_given())
1694 : {
1695 : // get diffusion at ips
1696 : const MathMatrix<dim, dim>* vDiffusion = NULL;
1697 0 : if(m_imDiffusion.data_given()) vDiffusion = m_imDiffusion.values();
1698 :
1699 : // update convection shapes
1700 0 : if(!m_spConvShape->update(&geo, m_imVelocity.values(), vDiffusion, compute_deriv))
1701 : {
1702 : UG_LOG("ERROR in 'ConvectionDiffusionFV1::get_updated_conv_shapes': "
1703 : "Cannot compute convection shapes.\n");
1704 : }
1705 : }
1706 :
1707 : // return a const (!!) reference to the upwind
1708 0 : return *const_cast<const IConvectionShapes<dim>*>(m_spConvShape.get());
1709 : }
1710 :
1711 :
1712 : ////////////////////////////////////////////////////////////////////////////////
1713 : // register assemble functions
1714 : ////////////////////////////////////////////////////////////////////////////////
1715 :
1716 : /**
1717 : * Registers the assembling functions for all the element types
1718 : */
1719 : template<typename TDomain>
1720 0 : void ConvectionDiffusionFV1<TDomain>::
1721 : register_all_funcs(bool bHang)
1722 : {
1723 : // list of element types for assembling
1724 : typedef typename domain_traits<dim>::AllElemList AssembleElemList;
1725 :
1726 : // register the local assembling functions
1727 0 : boost::mpl::for_each<AssembleElemList> (RegisterLocalDiscr (this, bHang));
1728 0 : }
1729 :
1730 : /**
1731 : * Registers the assembling functions for an element type.
1732 : */
1733 : template<typename TDomain>
1734 : template<typename TElem>
1735 : void
1736 0 : ConvectionDiffusionFV1<TDomain>::
1737 : register_func_for_(bool bHang)
1738 : {
1739 : // switch assemble functions
1740 0 : if(!bHang)
1741 : {
1742 0 : if(!m_bCondensedFV)
1743 0 : register_func<TElem, FV1Geometry<TElem, dim> >();
1744 : else
1745 0 : register_func<TElem, FV1CondensedGeometry<TElem, dim> >();
1746 : }
1747 : else
1748 : {
1749 0 : if(m_bCondensedFV)
1750 0 : UG_THROW("ConvectionDiffusionFV1: Condensed FV not supported for hanging nodes.");
1751 0 : register_func<TElem, HFV1Geometry<TElem, dim> >();
1752 : }
1753 0 : }
1754 :
1755 : template<typename TDomain>
1756 : template<typename TElem, typename TFVGeom>
1757 0 : void ConvectionDiffusionFV1<TDomain>::
1758 : register_func()
1759 : {
1760 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1761 : typedef this_type T;
1762 : static const int refDim = reference_element_traits<TElem>::dim;
1763 :
1764 : this->clear_add_fct(id);
1765 : this->set_prep_elem_loop_fct(id, &T::template prep_elem_loop<TElem, TFVGeom>);
1766 : this->set_prep_elem_fct( id, &T::template prep_elem<TElem, TFVGeom>);
1767 : this->set_fsh_elem_loop_fct( id, &T::template fsh_elem_loop<TElem, TFVGeom>);
1768 : this->set_add_jac_A_elem_fct(id, &T::template add_jac_A_elem<TElem, TFVGeom>);
1769 : this->set_add_jac_M_elem_fct(id, &T::template add_jac_M_elem<TElem, TFVGeom>);
1770 : this->set_add_def_A_elem_fct(id, &T::template add_def_A_elem<TElem, TFVGeom>);
1771 : this->set_add_def_A_expl_elem_fct(id, &T::template add_def_A_expl_elem<TElem, TFVGeom>);
1772 : this->set_add_def_M_elem_fct(id, &T::template add_def_M_elem<TElem, TFVGeom>);
1773 : this->set_add_rhs_elem_fct( id, &T::template add_rhs_elem<TElem, TFVGeom>);
1774 :
1775 : // error estimator parts
1776 : this->set_prep_err_est_elem_loop(id, &T::template prep_err_est_elem_loop<TElem, TFVGeom>);
1777 : this->set_prep_err_est_elem(id, &T::template prep_err_est_elem<TElem, TFVGeom>);
1778 : this->set_compute_err_est_A_elem(id, &T::template compute_err_est_A_elem<TElem, TFVGeom>);
1779 : this->set_compute_err_est_M_elem(id, &T::template compute_err_est_M_elem<TElem, TFVGeom>);
1780 : this->set_compute_err_est_rhs_elem(id, &T::template compute_err_est_rhs_elem<TElem, TFVGeom>);
1781 : this->set_fsh_err_est_elem_loop(id, &T::template fsh_err_est_elem_loop<TElem, TFVGeom>);
1782 :
1783 : // set computation of linearized defect w.r.t velocity
1784 0 : m_imDiffusion.set_fct(id, this, &T::template lin_def_diffusion<TElem, TFVGeom>);
1785 0 : m_imVelocity. set_fct(id, this, &T::template lin_def_velocity<TElem, TFVGeom>);
1786 0 : m_imFlux.set_fct(id, this, &T::template lin_def_flux<TElem, TFVGeom>);
1787 0 : m_imReactionRate. set_fct(id, this, &T::template lin_def_reaction_rate<TElem, TFVGeom>);
1788 0 : m_imReaction. set_fct(id, this, &T::template lin_def_reaction<TElem, TFVGeom>);
1789 0 : m_imSource. set_fct(id, this, &T::template lin_def_source<TElem, TFVGeom>);
1790 0 : m_imVectorSource.set_fct(id, this, &T::template lin_def_vector_source<TElem, TFVGeom>);
1791 0 : m_imMassScale.set_fct(id, this, &T::template lin_def_mass_scale<TElem, TFVGeom>);
1792 0 : m_imMass. set_fct(id, this, &T::template lin_def_mass<TElem, TFVGeom>);
1793 :
1794 : // exports
1795 0 : m_exValue-> template set_fct<T,refDim>(id, this, &T::template ex_value<TElem, TFVGeom>);
1796 0 : m_exGrad->template set_fct<T,refDim>(id, this, &T::template ex_grad<TElem, TFVGeom>);
1797 0 : }
1798 :
1799 : ////////////////////////////////////////////////////////////////////////////////
1800 : // explicit template instantiations
1801 : ////////////////////////////////////////////////////////////////////////////////
1802 :
1803 : #ifdef UG_DIM_1
1804 : template class ConvectionDiffusionFV1<Domain1d>;
1805 : #endif
1806 : #ifdef UG_DIM_2
1807 : template class ConvectionDiffusionFV1<Domain2d>;
1808 : #endif
1809 : #ifdef UG_DIM_3
1810 : template class ConvectionDiffusionFV1<Domain3d>;
1811 : #endif
1812 :
1813 : } // end namespace ConvectionDiffusionPlugin
1814 : } // namespace ug
1815 :
|
