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