Line data Source code
1 : /*
2 : * Copyright (c) 2013-2018: G-CSC, Goethe University Frankfurt
3 : * Author: Dmitry Logashenko
4 : * Based on the modules by Andreas Vogel
5 : *
6 : * This file is part of UG4.
7 : *
8 : * UG4 is free software: you can redistribute it and/or modify it under the
9 : * terms of the GNU Lesser General Public License version 3 (as published by the
10 : * Free Software Foundation) with the following additional attribution
11 : * requirements (according to LGPL/GPL v3 §7):
12 : *
13 : * (1) The following notice must be displayed in the Appropriate Legal Notices
14 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
15 : *
16 : * (2) The following notice must be displayed at a prominent place in the
17 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
18 : *
19 : * (3) The following bibliography is recommended for citation and must be
20 : * preserved in all covered files:
21 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
22 : * parallel geometric multigrid solver on hierarchically distributed grids.
23 : * Computing and visualization in science 16, 4 (2013), 151-164"
24 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
25 : * flexible software system for simulating pde based models on high performance
26 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
27 : *
28 : * This program is distributed in the hope that it will be useful,
29 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
30 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 : * GNU Lesser General Public License for more details.
32 : */
33 :
34 : #include "convection_diffusion_fractfv1.h"
35 :
36 : #include "lib_disc/spatial_disc/disc_util/geom_provider.h"
37 : #include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
38 : #include "lib_disc/spatial_disc/disc_util/hfv1_geom.h"
39 : #include "lib_disc/spatial_disc/disc_util/conv_shape.h"
40 :
41 : namespace ug{
42 : namespace ConvectionDiffusionPlugin{
43 :
44 : ////////////////////////////////////////////////////////////////////////////////
45 : // general
46 : ////////////////////////////////////////////////////////////////////////////////
47 :
48 : template<typename TDomain>
49 0 : ConvectionDiffusionFractFV1<TDomain>::
50 : ConvectionDiffusionFractFV1 (const char* functions, const char* subsets)
51 : : ConvectionDiffusionBase<TDomain> (functions, subsets),
52 0 : m_spConvShape (new ConvectionShapesNoUpwind<dim>)
53 : {
54 : // initialize the fracture-specific input-export parameters that are not in the base
55 : m_imAperture.set_comp_lin_defect (false);
56 :
57 : // register the fracture-specific input-export parameters that are not in the base
58 0 : this->register_import (m_imAperture);
59 0 : this->register_import (m_imOrthoDiffusion);
60 0 : this->register_import (m_imOrthoVelocity);
61 0 : this->register_import (m_imOrthoFlux);
62 0 : this->register_import (m_imOrthoVectorSource);
63 :
64 : // register the element assembling functions
65 : register_all_funcs ();
66 0 : }
67 :
68 : ////////////////////////////////////////////////////////////////////////////////
69 : // Local discretization interface
70 : ////////////////////////////////////////////////////////////////////////////////
71 :
72 : /// checks the grid and the shape functions
73 : template<typename TDomain>
74 0 : void ConvectionDiffusionFractFV1<TDomain>::prepare_setting
75 : (
76 : const std::vector<LFEID> & vLfeID,
77 : bool bNonRegular
78 : )
79 : {
80 : // check the grid
81 0 : if (bNonRegular)
82 0 : UG_THROW ("ERROR in ConvectionDiffusionFractFV1::prepare_setting:"
83 : " This discretization does not support hanging nodes.\n");
84 :
85 : // check number of the components
86 0 : if (vLfeID.size () != 1)
87 0 : UG_THROW ("ConvectionDiffusionFractFV1::prepare_setting: Wrong number of functions given. Need exactly 1.");
88 :
89 : // check whether these are the LagrangeP1 elements
90 : if (vLfeID[0] != LFEID(LFEID::LAGRANGE, dim, 1))
91 0 : UG_THROW ("ConvectionDiffusionFractFV1::prepare_setting: ConvectionDiffusion FV Scheme only implemented for 1st order.");
92 0 : }
93 :
94 : ////////////////////////////////////////////////////////////////////////////////
95 : // Assembling functions
96 : ////////////////////////////////////////////////////////////////////////////////
97 :
98 : /// computes and returns the upwind shapes for the given velocity
99 : template<typename TDomain>
100 : const typename ConvectionDiffusionFractFV1<TDomain>::conv_shape_type&
101 0 : ConvectionDiffusionFractFV1<TDomain>::get_updated_conv_shapes
102 : (
103 : bool computeDeriv ///< whether to compute the derivatives of the shapes
104 : )
105 : {
106 0 : if(m_imVelocity.data_given())
107 : {
108 : // get diffusion at ips
109 : const MathMatrix<dim, dim>* vDiffusion = NULL;
110 0 : if (m_imDiffusion.data_given ()) vDiffusion = m_imDiffusion.values ();
111 :
112 : // update convection shapes
113 0 : if (!m_spConvShape->update (m_pFractGeo, m_imVelocity.values (), vDiffusion, computeDeriv))
114 : {
115 : UG_LOG("ERROR in 'ConvectionDiffusionFractFV1::get_updated_conv_shapes': "
116 : "Cannot compute convection shapes.\n");
117 : }
118 : }
119 :
120 : // return a const (!!) reference to the upwind
121 0 : return *const_cast<const IConvectionShapes<dim>*>(m_spConvShape.get());
122 : }
123 :
124 : /// prepares the loop over the elements: checks whether the parameters are set, ...
125 : template<typename TDomain>
126 : template<typename TElem>
127 0 : void ConvectionDiffusionFractFV1<TDomain>::prep_elem_loop
128 : (
129 : ReferenceObjectID roid, ///< only elements with this roid are looped over
130 : int si ///< and only in this subdomain
131 : )
132 : {
133 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
134 :
135 : // check the imports
136 0 : if (!m_imAperture.data_given())
137 0 : UG_THROW ("ConvectionDiffusionFractFV1::prep_elem_loop: Missing Import 'fracture width (aperture)'.");
138 :
139 : // check whether we are in a degenerated fracture
140 0 : if (! m_spFractManager.valid())
141 0 : UG_THROW ("ConvectionDiffusionFractFV1::prep_elem_loop: No fracture manager specified.");
142 0 : if (! m_spFractManager->is_closed ())
143 0 : UG_THROW ("ConvectionDiffusionFractFV1::prep_elem_loop: Fracture manager not closed.");
144 0 : if (! m_spFractManager->contains (si))
145 0 : UG_THROW ("ConvectionDiffusionFractFV1::prep_elem_loop:"
146 : " The fract. discretization is used for subset " << si << " which is not a degenerated fracture.");
147 :
148 : // check, that upwind has been set
149 0 : if(m_spConvShape.invalid())
150 0 : UG_THROW("ConvectionDiffusionFractFV1::prep_elem_loop: Upwind has not been set.");
151 :
152 : // initialize the pointer to the FV geometry for fracture elements
153 0 : m_pFractGeo = &GeomProvider<TFractFVGeom>::get (LFEID (LFEID::LAGRANGE, low_dim, 1), 1);
154 :
155 : // set up local ip coordinates for corner import parameters
156 : // REMARK: Note that for the fracture elements, values of the corner import
157 : // parameters are indexed not by scv (as for the normal elements) but by
158 : // the corner indices in the reference element
159 0 : }
160 :
161 : template<typename TDomain>
162 : template<typename TElem>
163 0 : void ConvectionDiffusionFractFV1<TDomain>::fsh_elem_loop ()
164 0 : {}
165 :
166 : template<typename TDomain>
167 : template<typename TElem>
168 0 : void ConvectionDiffusionFractFV1<TDomain>::prep_elem
169 : (
170 : const LocalVector & u, ///< local solution at the dofs associated with elem
171 : GridObject * elem, ///< element to prepare
172 : ReferenceObjectID roid, // id of reference element used for assembling
173 : const position_type vCornerCoords [] ///< coordinates of the corners of the element
174 : )
175 : {
176 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
177 :
178 : static const int refDim = TElem::dim; // dimensionality of the element, not the side!
179 : TElem * pElem = static_cast<TElem*> (elem);
180 0 : ref_elem_type& rRefElem = Provider<ref_elem_type>::get ();
181 :
182 : // get the non-degenerated sides of the fracture element
183 : try
184 : {
185 : m_spFractManager->get_layer_sides
186 0 : (pElem,
187 0 : m_numFractCo, m_innerFractSide, m_innerFractSideIdx, m_innerSideCo,
188 0 : m_outerFractSide, m_outerFractSideIdx, m_outerSideCo,
189 0 : m_assCo);
190 : }
191 0 : UG_CATCH_THROW("ConvectionDiffusionFractFV1::prep_elem: Cannot find orientation of a fracture element.");
192 :
193 : // compute the FV geometry of the inner side
194 0 : position_type vSideCornerCoords [maxFractSideCorners]; // global side corner coords
195 : try
196 : {
197 0 : for (size_t co = 0; co < m_numFractCo; co++)
198 0 : vSideCornerCoords [co] = vCornerCoords [m_innerSideCo [co]];
199 0 : m_pFractGeo->update (m_innerFractSide, vSideCornerCoords, &(this->subset_handler()));
200 : }
201 0 : UG_CATCH_THROW("ConvectionDiffusionFractFV1::prep_elem: Cannot update the Finite Volume Geometry for a fracture element.");
202 0 : size_t numSCVFip = m_pFractGeo->num_scvf_ips ();
203 : size_t numSCVip = m_pFractGeo->num_scv_ips ();
204 :
205 : // convert local coordinates of the side into the local coordinates of the element (for the input parameters)
206 0 : position_type vSideLocCornerCoords [maxFractSideCorners]; // local side corner coords
207 0 : position_type vSideLocSCVFipCoords [TFractFVGeom::maxNumSCVF]; // local scvf ip's in a fracture element
208 0 : position_type vSideLocSCVipCoords [TFractFVGeom::maxNumSCV]; // local scv ip's in a fracture element
209 : position_type elemLocCOE; // local coordinates of the mass center of a fracture element
210 : try
211 : {
212 : DimReferenceMapping<low_dim, dim>& rMapping
213 0 : = ReferenceMappingProvider::get<low_dim, dim> (m_innerFractSide->reference_object_id ());
214 0 : for (size_t co = 0; co < m_numFractCo; co++)
215 0 : vSideLocCornerCoords [co] = rRefElem.corner (m_innerSideCo [co]);
216 0 : rMapping.update (vSideLocCornerCoords);
217 :
218 0 : rMapping.local_to_global (elemLocCOE, *(m_pFractGeo->coe_local ()));
219 0 : rMapping.local_to_global (vSideLocSCVFipCoords, m_pFractGeo->scvf_local_ips (), numSCVFip);
220 0 : rMapping.local_to_global (vSideLocSCVipCoords, m_pFractGeo->scv_local_ips (), numSCVip);
221 : }
222 0 : UG_CATCH_THROW("ConvectionDiffusionFractFV1::prep_elem: Cannot transform local side coordinates to local element coordinates in a fracture element.");
223 :
224 :
225 : // set local positions
226 :
227 0 : m_imDiffusion.template set_local_ips<refDim> (vSideLocSCVFipCoords, numSCVFip);
228 0 : m_imVelocity.template set_local_ips<refDim> (vSideLocSCVFipCoords, numSCVFip);
229 0 : m_imFlux.template set_local_ips<refDim> (vSideLocSCVFipCoords, numSCVFip);
230 0 : m_imVectorSource.template set_local_ips<refDim> (vSideLocSCVFipCoords, numSCVFip);
231 :
232 0 : m_imSource.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
233 0 : m_imReactionRate.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
234 0 : m_imReaction.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
235 0 : m_imReactionRateExpl.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
236 0 : m_imReactionExpl.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
237 0 : m_imSourceExpl.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
238 0 : m_imMassScale.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
239 0 : m_imMass.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
240 :
241 0 : m_imOrthoDiffusion.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
242 0 : m_imOrthoVelocity.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
243 0 : m_imOrthoFlux.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
244 0 : m_imOrthoVectorSource.template set_local_ips<refDim> (vSideLocSCVipCoords, numSCVip);
245 :
246 0 : m_imAperture.template set_local_ips<dim> (&elemLocCOE, 1);
247 :
248 : // set global positions
249 :
250 0 : const position_type* vSCVFip = m_pFractGeo->scvf_global_ips ();
251 0 : m_imDiffusion. set_global_ips (vSCVFip, numSCVFip);
252 0 : m_imVelocity. set_global_ips (vSCVFip, numSCVFip);
253 0 : m_imFlux. set_global_ips (vSCVFip, numSCVFip);
254 0 : m_imVectorSource. set_global_ips (vSCVFip, numSCVFip);
255 :
256 0 : const position_type* vSCVip = m_pFractGeo->scv_global_ips ();
257 0 : m_imSource. set_global_ips (vSCVip, numSCVip);
258 0 : m_imReactionRate. set_global_ips (vSCVip, numSCVip);
259 0 : m_imReactionRateExpl. set_global_ips (vSCVip, numSCVip);
260 0 : m_imReactionExpl. set_global_ips (vSCVip, numSCVip);
261 0 : m_imSourceExpl. set_global_ips (vSCVip, numSCVip);
262 0 : m_imReaction. set_global_ips (vSCVip, numSCVip);
263 0 : m_imMassScale. set_global_ips (vSCVip, numSCVip);
264 0 : m_imMass. set_global_ips (vSCVip, numSCVip);
265 :
266 0 : const position_type* coe_global = m_pFractGeo->coe_global ();
267 0 : m_imAperture.set_global_ips (coe_global, 1);
268 :
269 : // init upwind for element type
270 0 : if(m_spConvShape.valid())
271 0 : if(!m_spConvShape->template set_geometry_type<TFractFVGeom> (*m_pFractGeo))
272 0 : UG_THROW("ConvectionDiffusionFractFV1::prep_elem:"
273 : " Cannot init upwind for element type.");
274 0 : }
275 :
276 : /// computes the local stiffness matrix
277 : template<typename TDomain>
278 : template<typename TElem>
279 0 : void ConvectionDiffusionFractFV1<TDomain>::add_jac_A_elem
280 : (
281 : LocalMatrix & J, ///< the local matrix to update
282 : const LocalVector & u, ///< current approximation of the solution
283 : GridObject * elem, ///< grid element
284 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
285 : )
286 : {
287 : TElem * pElem = static_cast<TElem*> (elem);
288 :
289 0 : this->template fract_add_jac_A_elem<TElem> (J, u, pElem, vCornerCoords);
290 0 : this->template fract_bulk_add_jac_A_elem<TElem> (J, u, pElem, vCornerCoords);
291 0 : }
292 :
293 : /// assembles a singular source or sink in the jacobian
294 : template<typename TDomain>
295 : template<typename TElem>
296 : void ConvectionDiffusionFractFV1<TDomain>::add_sss_jac_elem
297 : (
298 : LocalMatrix& J, ///< the matrix to update
299 : const LocalVector& u, ///< current solution
300 : TElem* pElem, ///< the element
301 : size_t co, ///< corner for the source/sink
302 : number flux ///< flux through the source/sink (premultiplied by the length for lines)
303 : )
304 : {
305 0 : if (flux < 0.0)
306 : // sink
307 0 : J(_C_, co, _C_, co) -= flux;
308 : }
309 :
310 : /// computes the local stiffness matrix on a fracture element
311 : template<typename TDomain>
312 : template<typename TElem>
313 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_add_jac_A_elem
314 : (
315 : LocalMatrix & J, ///< the local matrix to update
316 : const LocalVector & u, ///< current approximation of the solution
317 : TElem * pElem, ///< grid element
318 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
319 : )
320 : {
321 0 : number half_fr_width = m_imAperture[0] / 2;
322 :
323 0 : const size_t numSh = m_pFractGeo->num_sh ();
324 : const size_t numScvf = m_pFractGeo->num_scvf ();
325 :
326 : // Diffusion and Velocity Term
327 0 : if (m_imDiffusion.data_given () || m_imVelocity.data_given ())
328 : {
329 : // get conv shapes
330 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes (false);
331 :
332 : // loop Sub Control Volume Faces (SCVF)
333 0 : for (size_t ip = 0; ip < numScvf; ++ip)
334 : {
335 : // get current SCVF
336 0 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
337 :
338 : ////////////////////////////////////////////////////
339 : // Diffusive Term
340 : ////////////////////////////////////////////////////
341 0 : if (m_imDiffusion.data_given ())
342 : {
343 : // Diff. Tensor times Gradient
344 : MathVector<dim> Dgrad;
345 :
346 : // loop shape functions
347 0 : for (size_t sh = 0; sh < numSh; ++sh)
348 : {
349 : // Compute Diffusion Tensor times Gradient
350 : MatVecMult (Dgrad, m_imDiffusion[ip], scvf.global_grad (sh));
351 :
352 : // Compute flux at IP
353 0 : const number D_diff_flux = VecDot (Dgrad, scvf.normal ()) * half_fr_width;
354 :
355 0 : J(_C_, m_innerSideCo[scvf.from()], _C_, m_innerSideCo[sh]) -= D_diff_flux;
356 0 : J(_C_, m_innerSideCo[scvf.to() ], _C_, m_innerSideCo[sh]) += D_diff_flux;
357 : }
358 : }
359 :
360 : ////////////////////////////////////////////////////
361 : // Convective Term
362 : ////////////////////////////////////////////////////
363 0 : if (m_imVelocity.data_given ())
364 : {
365 : // Add Flux contribution
366 0 : for (size_t sh = 0; sh < numSh; ++sh)
367 : {
368 0 : const number D_conv_flux = convShape (ip, sh) * half_fr_width;
369 :
370 : // Add flux term to local matrix
371 0 : J(_C_, m_innerSideCo[scvf.from()], _C_, m_innerSideCo[sh]) += D_conv_flux;
372 0 : J(_C_, m_innerSideCo[scvf.to() ], _C_, m_innerSideCo[sh]) -= D_conv_flux;
373 : }
374 : }
375 :
376 : // no explicit dependency on flux import
377 : }
378 : }
379 :
380 : // Reaction rate
381 0 : if(m_imReactionRate.data_given())
382 : {
383 : // loop Sub Control Volume (SCV)
384 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
385 : {
386 : // get current SCV
387 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
388 :
389 : // get associated node
390 0 : const int co = m_innerSideCo [scv.node_id ()];
391 :
392 : // Add to local matrix
393 0 : J(_C_, co, _C_, co) += m_imReactionRate[ip] * scv.volume () * half_fr_width;
394 : }
395 : }
396 :
397 : // Reaction term does not explicitly depend on the associated unknown function
398 :
399 : // Assemble the singular sources and sinks
400 0 : if (m_sss_mngr.valid () && m_sss_mngr->num_lines () != 0)
401 : {
402 : typedef typename domain_type::position_accessor_type t_pos_accessor;
403 : typedef typename CDSingularSourcesAndSinks<dim>::template
404 : line_iterator<side_type,t_pos_accessor,TFractFVGeom> t_lin_sss_iter;
405 :
406 0 : t_pos_accessor& aaPos = this->domain()->position_accessor ();
407 0 : Grid& grid = (Grid&) *this->domain()->grid ();
408 :
409 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
410 : {
411 : // Get the corner of the face
412 : size_t side_co = m_pFractGeo->scv(ip).node_id ();
413 : // Get associated node of the element (not side!)
414 0 : size_t co = m_innerSideCo [m_pFractGeo->scv(ip).node_id ()];
415 :
416 : // line sources (that correspond to the point sources)
417 0 : for (t_lin_sss_iter line (m_sss_mngr.get (), m_innerFractSide, grid, aaPos, *m_pFractGeo, side_co);
418 0 : ! line.is_over (); ++line)
419 : {
420 : FVLineSourceOrSink<dim, cd_line_sss_data<dim> > * line_sss = *line;
421 0 : if (! line_sss->marked_for (m_innerFractSide, side_co))
422 0 : continue;
423 : line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
424 0 : add_sss_jac_elem (J, u, pElem, co, line_sss->flux () / 2);
425 : /* Remark: "/ 2" because the source is taken into account twice. */
426 : }
427 : }
428 : }
429 0 : }
430 :
431 : /// computes the local stiffness matrix of the fracture-bulk interaction terms on a fracture element
432 : template<typename TDomain>
433 : template<typename TElem>
434 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_bulk_add_jac_A_elem
435 : (
436 : LocalMatrix & J, ///< the local matrix to update
437 : const LocalVector & u, ///< current approximation of the solution
438 : TElem * pElem, ///< grid element
439 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
440 : )
441 : {
442 0 : number half_fr_width = m_imAperture[0] / 2;
443 :
444 : // loop over the corners of the inner side
445 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
446 : {
447 : // Get current SCV
448 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
449 : number s = scv.volume ();
450 :
451 : // Get associated node of the element (not side!)
452 0 : const int co = m_innerSideCo [scv.node_id()];
453 :
454 : // Assemble Diffusion and Convection
455 : number D_flux = 0, D_flux_fr = 0;
456 :
457 0 : if (m_imOrthoDiffusion.data_given ())
458 : {
459 0 : D_flux = m_imOrthoDiffusion[ip] / half_fr_width;
460 0 : D_flux_fr = - D_flux;
461 : }
462 :
463 0 : if (m_imOrthoVelocity.data_given ())
464 : {
465 : /* We use the full upwind here: */
466 0 : if (m_imOrthoVelocity[ip] >= 0)
467 0 : D_flux_fr -= m_imOrthoVelocity[ip];
468 : else
469 0 : D_flux -= m_imOrthoVelocity[ip];
470 : }
471 :
472 0 : J(_C_, m_assCo [co], _C_, m_assCo [co]) += D_flux * s;
473 0 : J(_C_, m_assCo [co], _C_, co ) += D_flux_fr * s;
474 :
475 0 : J(_C_, co, _C_, m_assCo [co]) -= D_flux * s;
476 0 : J(_C_, co, _C_, co ) -= D_flux_fr * s;
477 : }
478 0 : }
479 :
480 :
481 : /// computes the mass matrix of a time-dependent problem
482 : template<typename TDomain>
483 : template<typename TElem>
484 0 : void ConvectionDiffusionFractFV1<TDomain>::add_jac_M_elem
485 : (
486 : LocalMatrix & J, ///< the local matrix to update
487 : const LocalVector & u, ///< current approximation of the solution
488 : GridObject * elem, ///< grid element
489 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
490 : )
491 : {
492 0 : if (!m_imMassScale.data_given ()) return;
493 :
494 0 : number half_fr_width = m_imAperture[0] / 2;
495 :
496 : // loop Sub Control Volumes (SCV)
497 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
498 : {
499 : // get current SCV
500 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
501 :
502 : // get associated node
503 0 : const int co = m_innerSideCo [scv.node_id ()];
504 :
505 : // Add to local matrix
506 0 : J(_C_, co, _C_, co) += scv.volume() * m_imMassScale[ip] * half_fr_width;
507 : }
508 :
509 : // m_imMass part does not explicitly depend on associated unknown function
510 : }
511 :
512 : /// computes the stiffness part of the local defect
513 : template<typename TDomain>
514 : template<typename TElem>
515 0 : void ConvectionDiffusionFractFV1<TDomain>::add_def_A_elem
516 : (
517 : LocalVector & d, ///< the local defect to update
518 : const LocalVector & u, ///< current approximation of the solution
519 : GridObject * elem, ///< grid element
520 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
521 : )
522 : {
523 : TElem * pElem = static_cast<TElem*> (elem);
524 :
525 0 : this->template fract_add_def_A_elem<TElem> (d, u, pElem, vCornerCoords);
526 0 : this->template fract_bulk_add_def_A_elem<TElem> (d, u, pElem, vCornerCoords);
527 0 : }
528 :
529 : /// assembles a singular source or sink in the defect
530 : template<typename TDomain>
531 : template<typename TElem>
532 : void ConvectionDiffusionFractFV1<TDomain>::add_sss_def_elem
533 : (
534 : LocalVector& d, ///< the defect to update
535 : const LocalVector& u, ///< current solution
536 : TElem * pElem, ///< the element
537 : size_t co, ///< corner for the source/sink
538 : number flux ///< flux through the source/sink (premultiplied by the length for lines)
539 : )
540 : {
541 0 : if (flux > 0.0)
542 : // source
543 0 : d(_C_, co) -= flux;
544 : else
545 : // sink
546 0 : d(_C_, co) -= flux * u(_C_, co);
547 : }
548 :
549 : /// computes the stiffness part of the local defect on a fracture element
550 : template<typename TDomain>
551 : template<typename TElem>
552 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_add_def_A_elem
553 : (
554 : LocalVector & d, ///< the local defect to update
555 : const LocalVector & u, ///< current approximation of the solution
556 : TElem * pElem, ///< grid element
557 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
558 : )
559 : {
560 0 : number half_fr_width = m_imAperture[0] / 2;
561 :
562 : // Diff. Tensor times Gradient
563 : MathVector<dim> Dgrad;
564 :
565 0 : const size_t numSh = m_pFractGeo->num_sh ();
566 : const size_t numScvf = m_pFractGeo->num_scvf ();
567 :
568 : // Diffusion and Velocity Term
569 0 : if (m_imDiffusion.data_given () || m_imVelocity.data_given () || m_imFlux.data_given ())
570 : {
571 : // get conv shapes
572 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes (false);
573 :
574 : // loop Sub Control Volume Faces (SCVF)
575 0 : for(size_t ip = 0; ip < numScvf; ++ip)
576 : {
577 : // get current SCVF
578 0 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
579 :
580 : ////////////////////////////////////////////////////
581 : // Diffusive Term
582 : ////////////////////////////////////////////////////
583 0 : if(m_imDiffusion.data_given ())
584 : {
585 : // to compute D \nabla c
586 : MathVector<dim> Dgrad_c, grad_c;
587 :
588 : // compute gradient and shape at ip
589 : VecSet (grad_c, 0.0);
590 0 : for(size_t sh = 0; sh < numSh; ++sh)
591 0 : VecScaleAppend (grad_c,
592 : u (_C_, m_innerSideCo[sh]), scvf.global_grad (sh));
593 :
594 : // scale by diffusion tensor
595 : MatVecMult (Dgrad_c, m_imDiffusion[ip], grad_c);
596 :
597 : // Compute flux
598 0 : const number diff_flux = VecDot (Dgrad_c, scvf.normal ()) * half_fr_width;
599 :
600 : // Add to local defect
601 0 : d(_C_, m_innerSideCo [scvf.from()]) -= diff_flux;
602 0 : d(_C_, m_innerSideCo [scvf.to() ]) += diff_flux;
603 : }
604 :
605 : ////////////////////////////////////////////////////
606 : // Convective Term
607 : ////////////////////////////////////////////////////
608 0 : if(m_imVelocity.data_given())
609 : {
610 : // sum up convective flux using convection shapes
611 : number conv_flux = 0.0;
612 0 : for(size_t sh = 0; sh < numSh; ++sh)
613 0 : conv_flux += u (_C_, m_innerSideCo[sh]) * convShape (ip, sh);
614 0 : conv_flux *= half_fr_width;
615 :
616 : // add to local defect
617 0 : d(_C_, m_innerSideCo [scvf.from()]) += conv_flux;
618 0 : d(_C_, m_innerSideCo [scvf.to() ]) -= conv_flux;
619 : }
620 :
621 : /////////////////////////////////////////////////////
622 : // Flux Term
623 : /////////////////////////////////////////////////////
624 0 : if(m_imFlux.data_given())
625 : {
626 : // sum up flux
627 0 : const number flux = VecDot (m_imFlux[ip], scvf.normal ()) * half_fr_width;
628 :
629 : // add to local defect
630 0 : d(_C_, m_innerSideCo [scvf.from()]) += flux;
631 0 : d(_C_, m_innerSideCo [scvf.to() ]) -= flux;
632 : }
633 : }
634 : }
635 :
636 : // Reaction rate
637 0 : if (m_imReactionRate.data_given ())
638 : {
639 : // loop Sub Control Volumes (SCV)
640 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
641 : {
642 : // get current SCV
643 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
644 :
645 : // get associated node
646 0 : const int co = m_innerSideCo [scv.node_id ()];
647 :
648 : // Add to local defect
649 0 : d(_C_, co) += u(_C_, co) * m_imReactionRate[ip] * scv.volume () * half_fr_width;
650 : }
651 : }
652 :
653 : // Reaction term
654 0 : if (m_imReaction.data_given ())
655 : {
656 : // loop Sub Control Volumes (SCV)
657 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
658 : {
659 : // get current SCV
660 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
661 :
662 : // get associated node
663 0 : const int co = m_innerSideCo [scv.node_id ()];
664 :
665 : // Add to local defect
666 0 : d(_C_, co) += m_imReaction[ip] * scv.volume () * half_fr_width;
667 : }
668 : }
669 :
670 : // Assemble the singular sources and sinks
671 0 : if (m_sss_mngr.valid () && m_sss_mngr->num_lines () != 0)
672 : {
673 : typedef typename domain_type::position_accessor_type t_pos_accessor;
674 : typedef typename CDSingularSourcesAndSinks<dim>::template
675 : line_iterator<side_type,t_pos_accessor,TFractFVGeom> t_lin_sss_iter;
676 :
677 0 : t_pos_accessor& aaPos = this->domain()->position_accessor ();
678 0 : Grid& grid = (Grid&) *this->domain()->grid ();
679 :
680 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
681 : {
682 : // Get the corner of the face
683 : size_t side_co = m_pFractGeo->scv(ip).node_id ();
684 : // Get associated node of the element (not side!)
685 0 : size_t co = m_innerSideCo [m_pFractGeo->scv(ip).node_id ()];
686 :
687 : // line sources (that correspond to the point sources)
688 0 : for (t_lin_sss_iter line (m_sss_mngr.get (), m_innerFractSide, grid, aaPos, *m_pFractGeo, side_co);
689 0 : ! line.is_over (); ++line)
690 : {
691 : FVLineSourceOrSink<dim, cd_line_sss_data<dim> > * line_sss = *line;
692 0 : if (! line_sss->marked_for (m_innerFractSide, side_co))
693 0 : continue;
694 : line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
695 0 : add_sss_def_elem (d, u, pElem, co, line_sss->flux () / 2);
696 : /* Remark: "/ 2" because the source is taken into account twice. */
697 : }
698 : }
699 : }
700 0 : }
701 :
702 : /// computes the stiffness fracture-bulk interaction terms of the local defect on a fracture element
703 : template<typename TDomain>
704 : template<typename TElem>
705 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_bulk_add_def_A_elem
706 : (
707 : LocalVector & d, ///< the local defect to update
708 : const LocalVector & u, ///< current approximation of the solution
709 : TElem * pElem, ///< grid element
710 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
711 : )
712 : {
713 : number orthC_f, orthC_m;
714 :
715 0 : number half_fr_width = m_imAperture[0] / 2;
716 :
717 : // loop over the corners of the inner side
718 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
719 : {
720 : // Get current SCV
721 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
722 :
723 : // Get associated node of the element (not side!)
724 0 : const int co = m_innerSideCo [scv.node_id()];
725 :
726 : // Get the corner values
727 0 : orthC_f = u(_C_, co);
728 0 : orthC_m = u(_C_, m_assCo[co]);
729 :
730 : // Assemble Diffusion, Convection and the Flux
731 : number flux = 0;
732 :
733 0 : if (m_imOrthoDiffusion.data_given ())
734 0 : flux = m_imOrthoDiffusion[ip] * (orthC_m - orthC_f) / half_fr_width;
735 :
736 0 : if (m_imOrthoVelocity.data_given ())
737 : {
738 0 : number orthVelocity = m_imOrthoVelocity[ip];
739 : /* We use the full upwind here: */
740 0 : flux -= orthVelocity * ((orthVelocity >= 0)? orthC_f : orthC_m);
741 : }
742 :
743 0 : if (m_imOrthoFlux.data_given ())
744 0 : flux -= m_imOrthoFlux[ip];
745 :
746 0 : flux *= scv.volume ();
747 0 : d(_C_, m_assCo[co]) += flux;
748 0 : d(_C_, co) -= flux;
749 : }
750 0 : }
751 :
752 : /// computes the stiffness part of the local defect
753 : template<typename TDomain>
754 : template<typename TElem>
755 0 : void ConvectionDiffusionFractFV1<TDomain>::add_def_A_expl_elem
756 : (
757 : LocalVector & d, ///< the local defect to update
758 : const LocalVector & u, ///< current approximation of the solution
759 : GridObject * elem, ///< grid element
760 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
761 : )
762 : {
763 : TElem * pElem = static_cast<TElem*> (elem);
764 :
765 0 : this->template fract_add_def_A_expl_elem<TElem> (d, u, pElem, vCornerCoords);
766 : this->template fract_bulk_add_def_A_expl_elem<TElem> (d, u, pElem, vCornerCoords);
767 0 : }
768 :
769 : /// computes the stiffness part of the local defect on a fracture element
770 : template<typename TDomain>
771 : template<typename TElem>
772 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_add_def_A_expl_elem
773 : (
774 : LocalVector & d, ///< the local defect to update
775 : const LocalVector & u, ///< current approximation of the solution
776 : TElem * pElem, ///< grid element
777 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
778 : )
779 : {
780 0 : number half_fr_width = m_imAperture[0] / 2;
781 :
782 : // Reaction rate
783 0 : if (m_imReactionRateExpl.data_given ())
784 : {
785 : // loop Sub Control Volumes (SCV)
786 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
787 : {
788 : // get current SCV
789 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
790 :
791 : // get associated node
792 0 : const int co = m_innerSideCo [scv.node_id ()];
793 :
794 : // Add to local defect
795 0 : d(_C_, co) += u(_C_, co) * m_imReactionRateExpl[ip] * scv.volume () * half_fr_width;
796 : }
797 : }
798 :
799 : // reaction
800 0 : if (m_imReactionExpl.data_given ())
801 : {
802 : // loop Sub Control Volumes (SCV)
803 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
804 : {
805 : // get current SCV
806 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
807 :
808 : // get associated node
809 0 : const int co = m_innerSideCo [scv.node_id ()];
810 :
811 : // Add to local defect
812 0 : d(_C_, co) += m_imReactionExpl[ip] * scv.volume () * half_fr_width;
813 : }
814 : }
815 :
816 : // source
817 0 : if (m_imSourceExpl.data_given ())
818 : {
819 : // loop Sub Control Volumes (SCV)
820 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
821 : {
822 : // get current SCV
823 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
824 :
825 : // get associated node
826 0 : const int co = m_innerSideCo [scv.node_id ()];
827 :
828 : // Add to local rhs
829 0 : d(_C_, co) -= m_imSourceExpl[ip] * scv.volume () * half_fr_width;
830 : }
831 : }
832 0 : }
833 :
834 : /// computes the stiffness fracture-bulk interaction terms of the local defect on a fracture element
835 : template<typename TDomain>
836 : template<typename TElem>
837 : void ConvectionDiffusionFractFV1<TDomain>::fract_bulk_add_def_A_expl_elem
838 : (
839 : LocalVector & d, ///< the local defect to update
840 : const LocalVector & u, ///< current approximation of the solution
841 : TElem * pElem, ///< grid element
842 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
843 : )
844 : {
845 : }
846 :
847 : /// computes the mass part of the defect of a time-dependent problem
848 : template<typename TDomain>
849 : template<typename TElem>
850 0 : void ConvectionDiffusionFractFV1<TDomain>::add_def_M_elem
851 : (
852 : LocalVector & d, ///< the local defect to update
853 : const LocalVector & u, ///< current approximation of the solution
854 : GridObject * elem, ///< grid element
855 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
856 : )
857 : {
858 0 : if (!m_imMassScale.data_given () && !m_imMass.data_given ()) return;
859 :
860 0 : number half_fr_width = m_imAperture[0] / 2;
861 :
862 : // loop Sub Control Volumes (SCV)
863 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
864 : {
865 : // get current SCV
866 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
867 :
868 : // get associated node
869 0 : const int co = m_innerSideCo [scv.node_id ()];
870 :
871 : // mass value
872 : number val = 0;
873 :
874 : // multiply by scaling
875 0 : if (m_imMassScale.data_given ())
876 0 : val += m_imMassScale[ip] * u(_C_, co);
877 :
878 : // add mass
879 0 : if (m_imMass.data_given ())
880 0 : val += m_imMass[ip];
881 :
882 : // Add to local defect
883 0 : d(_C_, co) += val * scv.volume () * half_fr_width;
884 : }
885 : }
886 :
887 : /// computes the right-hand side due to the sources
888 : template<typename TDomain>
889 : template<typename TElem>
890 0 : void ConvectionDiffusionFractFV1<TDomain>::add_rhs_elem
891 : (
892 : LocalVector & d, ///< the right-hand side to update
893 : GridObject * elem, ///< grid element
894 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
895 : )
896 : {
897 : TElem * pElem = static_cast<TElem*> (elem);
898 :
899 0 : this->template fract_add_rhs_elem<TElem> (d, pElem, vCornerCoords);
900 0 : this->template fract_bulk_add_rhs_elem<TElem> (d, pElem, vCornerCoords);
901 0 : }
902 :
903 : /// computes the right-hand side due to the sources on a fracture element
904 : template<typename TDomain>
905 : template<typename TElem>
906 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_add_rhs_elem
907 : (
908 : LocalVector & d, ///< the local defect to update
909 : TElem * pElem, ///< grid element
910 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
911 : )
912 : {
913 0 : number half_fr_width = m_imAperture[0] / 2;
914 :
915 0 : if (m_imSource.data_given ())
916 : // loop Sub Control Volumes (SCV)
917 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
918 : {
919 : // get current SCV
920 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
921 :
922 : // get associated node
923 0 : const int co = m_innerSideCo [scv.node_id ()];
924 :
925 : // Add to local rhs
926 0 : d(_C_, co) += m_imSource[ip] * scv.volume () * half_fr_width;
927 : }
928 :
929 0 : if (m_imVectorSource.data_given ())
930 : // loop Sub Control Volume Fraces (SCVF)
931 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scvf (); ++ip)
932 : {
933 : // get current SCVF
934 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
935 :
936 : // compute the flux
937 0 : number flux = VecDot (m_imVectorSource[ip], scvf.normal ()) * half_fr_width;
938 :
939 : // Add to local rhs
940 0 : d(_C_, m_innerSideCo [scvf.from()]) -= flux;
941 0 : d(_C_, m_innerSideCo [scvf.to() ]) += flux;
942 : }
943 0 : }
944 :
945 : /// computes the fracture-bulk interaction terms of the local rhs on a fracture element
946 : template<typename TDomain>
947 : template<typename TElem>
948 0 : void ConvectionDiffusionFractFV1<TDomain>::fract_bulk_add_rhs_elem
949 : (
950 : LocalVector & d, ///< the local defect to update
951 : TElem * pElem, ///< grid element
952 : const position_type vCornerCoords [] ///< corner coordinates of the grid element
953 : )
954 : {
955 0 : if (m_imOrthoVectorSource.data_given ())
956 : // loop Sub Control Volumes (SCV)
957 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
958 : {
959 : // Get current SCV
960 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
961 :
962 : // Get associated node of the element (not side!)
963 0 : const int co = m_innerSideCo [scv.node_id()];
964 :
965 0 : number flux = m_imOrthoVectorSource[ip] * scv.volume ();
966 0 : d(_C_, m_assCo[co]) += flux;
967 0 : d(_C_, co) -= flux;
968 : }
969 0 : }
970 :
971 : /// computes the linearized defect w.r.t to the fracture velocity
972 : template<typename TDomain>
973 : template<typename TElem>
974 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_velocity
975 : (
976 : const LocalVector& u,
977 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
978 : const size_t nip
979 : )
980 : {
981 : // reset the values for the linearized defect
982 0 : for(size_t ip = 0; ip < nip; ++ip)
983 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
984 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
985 : vvvLinDef[ip][c][sh] = 0.0;
986 :
987 : // get conv shapes
988 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes (true);
989 :
990 : // loop Sub Control Volume Faces (SCVF)
991 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scvf (); ++ip)
992 : {
993 : // get current SCVF
994 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
995 :
996 : // sum up contributions of convection shapes
997 : MathVector<dim> linDefect;
998 : VecSet(linDefect, 0.0);
999 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1000 0 : VecScaleAppend (linDefect, u (_C_, m_innerSideCo[sh]), convShape.D_vel (ip, sh));
1001 :
1002 : // add parts for both sides of scvf
1003 0 : vvvLinDef[ip][_C_][m_innerSideCo [scvf.from()]] += linDefect;
1004 0 : vvvLinDef[ip][_C_][m_innerSideCo [scvf.to()] ] -= linDefect;
1005 : }
1006 0 : }
1007 :
1008 : /// computes the linearized defect w.r.t to the orthogonal velocity
1009 : template<typename TDomain>
1010 : template<typename TElem>
1011 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_ortho_velocity
1012 : (
1013 : const LocalVector& u,
1014 : std::vector<std::vector<number> > vvvLinDef[],
1015 : const size_t nip
1016 : )
1017 : {
1018 : // reset the values for the linearized defect
1019 0 : for(size_t ip = 0; ip < nip; ++ip)
1020 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1021 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1022 0 : vvvLinDef[ip][c][sh] = 0.0;
1023 :
1024 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1025 : {
1026 : // Get current SCV
1027 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1028 :
1029 : // Get associated node of the element (not side!)
1030 0 : const int co = m_innerSideCo [scv.node_id()];
1031 :
1032 : // Get the corner values
1033 0 : number orthC_f = u(_C_, co);
1034 0 : number orthC_m = u(_C_, m_assCo[co]);
1035 :
1036 : // Compute the derivative of the flux
1037 0 : number D_flux_vel = - ((m_imOrthoVelocity[ip] >= 0)? orthC_f : orthC_m) * scv.volume ();
1038 :
1039 : // add parts for both sides of the fracture
1040 0 : vvvLinDef[ip][_C_][m_assCo[co]] += D_flux_vel;
1041 0 : vvvLinDef[ip][_C_][co ] -= D_flux_vel;
1042 : }
1043 0 : }
1044 :
1045 : /// computes the linearized defect w.r.t to the fracture diffusion
1046 : template<typename TDomain>
1047 : template<typename TElem>
1048 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_diffusion
1049 : (
1050 : const LocalVector& u,
1051 : std::vector<std::vector<MathMatrix<dim,dim> > > vvvLinDef[],
1052 : const size_t nip
1053 : )
1054 : {
1055 : // reset the values for the linearized defect
1056 0 : for(size_t ip = 0; ip < nip; ++ip)
1057 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1058 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1059 : vvvLinDef[ip][c][sh] = 0.0;
1060 :
1061 : // get conv shapes
1062 0 : const IConvectionShapes<dim>& convShape = get_updated_conv_shapes (true);
1063 :
1064 : // loop Sub Control Volume Faces (SCVF)
1065 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scvf (); ++ip)
1066 : {
1067 : // get current SCVF
1068 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
1069 :
1070 : // compute gradient at ip
1071 : MathVector<dim> grad_c;
1072 : VecSet (grad_c, 0.0);
1073 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1074 0 : VecScaleAppend (grad_c, u (_C_, m_innerSideCo[sh]), scvf.global_grad (sh));
1075 :
1076 : // compute the lin defect at this ip
1077 : MathMatrix<dim,dim> linDefect;
1078 :
1079 : // part coming from $-\nabla u * \vec{n}
1080 0 : for(size_t k=0; k < (size_t) dim; ++k)
1081 0 : for(size_t j = 0; j < (size_t)dim; ++j)
1082 0 : linDefect(j,k) = (scvf.normal())[j] * grad_c[k];
1083 :
1084 : // add contribution from convection shapes
1085 0 : if(convShape.non_zero_deriv_diffusion ())
1086 0 : for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1087 0 : MatAdd (linDefect,
1088 : convShape.D_diffusion (ip, sh), u (_C_, m_innerSideCo[sh]));
1089 :
1090 : // add contributions
1091 0 : vvvLinDef[ip][_C_][m_innerSideCo[scvf.from()]] -= linDefect;
1092 0 : vvvLinDef[ip][_C_][m_innerSideCo[scvf.to()] ] += linDefect;
1093 : }
1094 0 : }
1095 :
1096 : /// computes the linearized defect w.r.t to the orthogonal diffusion
1097 : template<typename TDomain>
1098 : template<typename TElem>
1099 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_ortho_diffusion
1100 : (
1101 : const LocalVector& u,
1102 : std::vector<std::vector<number> > vvvLinDef[],
1103 : const size_t nip
1104 : )
1105 : {
1106 : // reset the values for the linearized defect
1107 0 : for(size_t ip = 0; ip < nip; ++ip)
1108 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1109 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1110 0 : vvvLinDef[ip][c][sh] = 0.0;
1111 :
1112 : number orthC_f, orthC_m;
1113 :
1114 0 : number half_fr_width = m_imAperture[0] / 2;
1115 :
1116 : // loop over the corners of the inner side
1117 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1118 : {
1119 : // Get current SCV
1120 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1121 :
1122 : // Get associated node of the element (not side!)
1123 0 : const int co = m_innerSideCo [scv.node_id()];
1124 :
1125 : // Get the corner values
1126 0 : orthC_f = u(_C_, co);
1127 0 : orthC_m = u(_C_, m_assCo[co]);
1128 :
1129 0 : number linDefect = (orthC_m - orthC_f) * scv.volume () / half_fr_width;
1130 :
1131 : // add contributions
1132 0 : vvvLinDef[ip][_C_][m_assCo[co]] += linDefect;
1133 0 : vvvLinDef[ip][_C_][co ] -= linDefect;
1134 : }
1135 0 : }
1136 :
1137 : /// computes the linearized defect w.r.t to the fracture flux
1138 : template<typename TDomain>
1139 : template<typename TElem>
1140 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_flux
1141 : (
1142 : const LocalVector& u,
1143 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
1144 : const size_t nip
1145 : )
1146 : {
1147 : // reset the values for the linearized defect
1148 0 : for(size_t ip = 0; ip < nip; ++ip)
1149 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1150 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1151 : vvvLinDef[ip][c][sh] = 0.0;
1152 :
1153 : // loop Sub Control Volume Faces (SCVF)
1154 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scvf(); ++ip)
1155 : {
1156 : // get current SCVF
1157 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
1158 :
1159 : // add parts for both sides of scvf
1160 0 : vvvLinDef[ip][_C_][m_innerSideCo [scvf.from()]] += scvf.normal();
1161 0 : vvvLinDef[ip][_C_][m_innerSideCo [scvf.to() ]] -= scvf.normal();
1162 : }
1163 0 : }
1164 :
1165 : /// computes the linearized defect w.r.t to the orthogonal flux
1166 : template<typename TDomain>
1167 : template<typename TElem>
1168 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_ortho_flux
1169 : (
1170 : const LocalVector& u,
1171 : std::vector<std::vector<number> > vvvLinDef[],
1172 : const size_t nip
1173 : )
1174 : {
1175 : // reset the values for the linearized defect
1176 0 : for(size_t ip = 0; ip < nip; ++ip)
1177 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1178 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1179 0 : vvvLinDef[ip][c][sh] = 0.0;
1180 :
1181 : // loop over the corners of the inner side
1182 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1183 : {
1184 : // Get current SCV
1185 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1186 :
1187 : // Get associated node of the element (not side!)
1188 0 : const int co = m_innerSideCo [scv.node_id()];
1189 :
1190 : // add parts for both sides of the fracture
1191 0 : vvvLinDef[ip][_C_][m_assCo[co]] += scv.volume ();
1192 0 : vvvLinDef[ip][_C_][co ] -= scv.volume ();
1193 : }
1194 0 : }
1195 :
1196 : /// computes the linearized defect w.r.t to the reaction source
1197 : template<typename TDomain>
1198 : template<typename TElem>
1199 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_reaction
1200 : (
1201 : const LocalVector& u,
1202 : std::vector<std::vector<number> > vvvLinDef[],
1203 : const size_t nip
1204 : )
1205 : {
1206 : // reset the values for the linearized defect
1207 0 : for(size_t ip = 0; ip < nip; ++ip)
1208 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1209 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1210 0 : vvvLinDef[ip][c][sh] = 0.0;
1211 :
1212 0 : number half_fr_width = m_imAperture[0] / 2;
1213 :
1214 : // loop over the corners of the inner side
1215 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1216 : {
1217 : // Get current SCV
1218 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1219 :
1220 : // Get associated node of the element (not side!)
1221 0 : const int co = m_innerSideCo [scv.node_id()];
1222 :
1223 : // add parts for both sides of the fracture
1224 0 : vvvLinDef[ip][_C_][co] = scv.volume () * half_fr_width;
1225 : }
1226 0 : }
1227 :
1228 : /// computes the linearized defect w.r.t to the reaction rate
1229 : template<typename TDomain>
1230 : template<typename TElem>
1231 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_reaction_rate
1232 : (
1233 : const LocalVector& u,
1234 : std::vector<std::vector<number> > vvvLinDef[],
1235 : const size_t nip
1236 : )
1237 : {
1238 : // reset the values for the linearized defect
1239 0 : for(size_t ip = 0; ip < nip; ++ip)
1240 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1241 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1242 0 : vvvLinDef[ip][c][sh] = 0.0;
1243 :
1244 0 : number half_fr_width = m_imAperture[0] / 2;
1245 :
1246 : // loop over the corners of the inner side
1247 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1248 : {
1249 : // Get current SCV
1250 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1251 :
1252 : // Get associated node of the element (not side!)
1253 0 : const int co = m_innerSideCo [scv.node_id()];
1254 :
1255 : // add parts for both sides of the fracture
1256 0 : vvvLinDef[ip][_C_][co] = u(_C_, co) * scv.volume () * half_fr_width;
1257 : }
1258 0 : }
1259 :
1260 : /// computes the linearized defect w.r.t to the source
1261 : template<typename TDomain>
1262 : template<typename TElem>
1263 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_source
1264 : (
1265 : const LocalVector& u,
1266 : std::vector<std::vector<number> > vvvLinDef[],
1267 : const size_t nip
1268 : )
1269 : {
1270 : // reset the values for the linearized defect
1271 0 : for(size_t ip = 0; ip < nip; ++ip)
1272 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1273 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1274 0 : vvvLinDef[ip][c][sh] = 0.0;
1275 :
1276 0 : number half_fr_width = m_imAperture[0] / 2;
1277 :
1278 : // loop over the corners of the inner side
1279 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1280 : {
1281 : // Get current SCV
1282 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1283 :
1284 : // Get associated node of the element (not side!)
1285 0 : const int co = m_innerSideCo [scv.node_id()];
1286 :
1287 : // add parts for both sides of the fracture
1288 0 : vvvLinDef[ip][_C_][co] = scv.volume () * half_fr_width;
1289 : }
1290 0 : }
1291 :
1292 : /// computes the linearized defect w.r.t to the fracture flux
1293 : template<typename TDomain>
1294 : template<typename TElem>
1295 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_vector_source
1296 : (
1297 : const LocalVector& u,
1298 : std::vector<std::vector<MathVector<dim> > > vvvLinDef[],
1299 : const size_t nip
1300 : )
1301 : {
1302 : // reset the values for the linearized defect
1303 0 : for(size_t ip = 0; ip < nip; ++ip)
1304 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1305 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1306 : vvvLinDef[ip][c][sh] = 0.0;
1307 :
1308 : // loop Sub Control Volume Faces (SCVF)
1309 0 : for(size_t ip = 0; ip < m_pFractGeo->num_scvf(); ++ip)
1310 : {
1311 : // get current SCVF
1312 : const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf (ip);
1313 :
1314 : // add parts for both sides of scvf
1315 0 : vvvLinDef[ip][_C_][m_innerSideCo [scvf.from()]] -= scvf.normal();
1316 0 : vvvLinDef[ip][_C_][m_innerSideCo [scvf.to() ]] += scvf.normal();
1317 : }
1318 0 : }
1319 :
1320 : /// computes the linearized defect w.r.t to the orthogonal flux
1321 : template<typename TDomain>
1322 : template<typename TElem>
1323 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_ortho_vector_source
1324 : (
1325 : const LocalVector& u,
1326 : std::vector<std::vector<number> > vvvLinDef[],
1327 : const size_t nip
1328 : )
1329 : {
1330 : // reset the values for the linearized defect
1331 0 : for(size_t ip = 0; ip < nip; ++ip)
1332 0 : for(size_t c = 0; c < vvvLinDef[ip].size(); ++c)
1333 0 : for(size_t sh = 0; sh < vvvLinDef[ip][c].size(); ++sh)
1334 0 : vvvLinDef[ip][c][sh] = 0.0;
1335 :
1336 : // loop over the corners of the inner side
1337 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1338 : {
1339 : // Get current SCV
1340 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1341 :
1342 : // Get associated node of the element (not side!)
1343 0 : const int co = m_innerSideCo [scv.node_id()];
1344 :
1345 : // add parts for both sides of the fracture
1346 0 : vvvLinDef[ip][_C_][m_assCo[co]] -= scv.volume ();
1347 0 : vvvLinDef[ip][_C_][co ] += scv.volume ();
1348 : }
1349 0 : }
1350 :
1351 : // computes the linearized defect w.r.t to the mass scale
1352 : template<typename TDomain>
1353 : template <typename TElem>
1354 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_mass_scale
1355 : (
1356 : const LocalVector& u,
1357 : std::vector<std::vector<number> > vvvLinDef[],
1358 : const size_t nip
1359 : )
1360 : {
1361 0 : number half_fr_width = m_imAperture[0] / 2;
1362 :
1363 : // loop Sub Control Volumes (SCV)
1364 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
1365 : {
1366 : // get current SCV
1367 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
1368 :
1369 : // get associated node
1370 0 : const int co = m_innerSideCo [scv.node_id ()];
1371 :
1372 : // set lin defect
1373 0 : vvvLinDef[co][_C_][co] = u(_C_, co) * scv.volume () * half_fr_width;
1374 : }
1375 0 : }
1376 :
1377 : // computes the linearized defect w.r.t to the mass scale
1378 : template<typename TDomain>
1379 : template <typename TElem>
1380 0 : void ConvectionDiffusionFractFV1<TDomain>::lin_def_mass
1381 : (
1382 : const LocalVector& u,
1383 : std::vector<std::vector<number> > vvvLinDef[],
1384 : const size_t nip
1385 : )
1386 : {
1387 0 : number half_fr_width = m_imAperture[0] / 2;
1388 :
1389 : // loop Sub Control Volumes (SCV)
1390 0 : for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ++ip)
1391 : {
1392 : // get current SCV
1393 : const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv (ip);
1394 :
1395 : // get associated node
1396 0 : const int co = m_innerSideCo [scv.node_id ()];
1397 :
1398 : // set lin defect
1399 0 : vvvLinDef[co][_C_][co] = scv.volume () * half_fr_width;
1400 : }
1401 0 : }
1402 :
1403 : /// registers the local assembler functions for a given element
1404 : template<typename TDomain>
1405 : template<typename TElem>
1406 0 : void ConvectionDiffusionFractFV1<TDomain>::register_func ()
1407 : {
1408 : ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
1409 : typedef this_type T;
1410 :
1411 : this->clear_add_fct(id);
1412 :
1413 : this->set_prep_elem_loop_fct (id, &T::template prep_elem_loop<TElem>);
1414 : this->set_prep_elem_fct ( id, &T::template prep_elem<TElem>);
1415 : this->set_fsh_elem_loop_fct ( id, &T::template fsh_elem_loop<TElem>);
1416 : this->set_add_jac_A_elem_fct (id, &T::template add_jac_A_elem<TElem>);
1417 : this->set_add_jac_M_elem_fct (id, &T::template add_jac_M_elem<TElem>);
1418 : this->set_add_def_A_elem_fct (id, &T::template add_def_A_elem<TElem>);
1419 : this->set_add_def_A_expl_elem_fct(id, &T::template add_def_A_expl_elem<TElem>);
1420 : this->set_add_def_M_elem_fct (id, &T::template add_def_M_elem<TElem>);
1421 : this->set_add_rhs_elem_fct ( id, &T::template add_rhs_elem<TElem>);
1422 :
1423 : // set computation of linearized defect w.r.t velocity, diffusion etc.
1424 0 : m_imDiffusion.set_fct (id, this, &T::template lin_def_diffusion<TElem>);
1425 0 : m_imOrthoDiffusion.set_fct (id, this, &T::template lin_def_ortho_diffusion<TElem>);
1426 0 : m_imVelocity.set_fct (id, this, &T::template lin_def_velocity<TElem>);
1427 0 : m_imOrthoVelocity.set_fct (id, this, &T::template lin_def_ortho_velocity<TElem>);
1428 0 : m_imFlux.set_fct (id, this, &T::template lin_def_flux<TElem>);
1429 0 : m_imOrthoFlux.set_fct (id, this, &T::template lin_def_ortho_flux<TElem>);
1430 0 : m_imReactionRate.set_fct (id, this, &T::template lin_def_reaction_rate<TElem>);
1431 0 : m_imReaction.set_fct (id, this, &T::template lin_def_reaction<TElem>);
1432 0 : m_imSource.set_fct (id, this, &T::template lin_def_source<TElem>);
1433 0 : m_imVectorSource.set_fct (id, this, &T::template lin_def_vector_source<TElem>);
1434 0 : m_imOrthoVectorSource.set_fct (id, this, &T::template lin_def_ortho_vector_source<TElem>);
1435 0 : m_imMassScale.set_fct (id, this, &T::template lin_def_mass_scale<TElem>);
1436 0 : m_imMass.set_fct (id, this, &T::template lin_def_mass<TElem>);
1437 0 : }
1438 :
1439 : /// registers the interface functions for all types of the elements
1440 : template<typename TDomain>
1441 0 : void ConvectionDiffusionFractFV1<TDomain>::register_all_funcs ()
1442 : {
1443 : typedef typename domain_traits<dim>::DimElemList AssembleElemList;
1444 :
1445 0 : boost::mpl::for_each<AssembleElemList> (RegisterLocalDiscr (this));
1446 0 : }
1447 :
1448 : ////////////////////////////////////////////////////////////////////////////////
1449 : // explicit template instantiations
1450 : ////////////////////////////////////////////////////////////////////////////////
1451 :
1452 : #ifdef UG_DIM_2
1453 : template class ConvectionDiffusionFractFV1<Domain2d>;
1454 : #endif
1455 : #ifdef UG_DIM_3
1456 : template class ConvectionDiffusionFractFV1<Domain3d>;
1457 : #endif
1458 :
1459 : } // end namespace ConvectionDiffusionPlugin
1460 : } // namespace ug
1461 :
1462 : /* End of File */
|
