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