Plugins
Loading...
Searching...
No Matches
fract_ddf_fv1_impl.h
Go to the documentation of this file.
1/*
2 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
3 * SPDX-License-Identifier: EUPL-1.2
4 * SPDX-FileContributor: Dmitry Logashenko, Sabine Stichel
5 * SPDX-FileContributor: Goethe Universität Frankfurt
6 * SPDX-FileType: SOURCE
7 *
8 * This file is part of d3f++.
9 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
10 */
11
12
13/*
14 * fract_ddf_fv1_impl.h - class member implementation of the discretization of
15 * the FV-discretization of the density-driven flow in fractured porous media.
16 *
17 * Based on the modules by A. Vogel.
18 */
19
20/* UG4 headers: */
21#include "lib_disc/spatial_disc/disc_util/geom_provider.h"
22#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
23#include "lib_disc/spatial_disc/user_data/const_user_data.h"
24#include "lib_disc/spatial_disc/disc_util/conv_shape.h"
25#ifdef UG_FOR_LUA
26#include "bindings/lua/lua_user_data.h"
27#endif
28
29namespace ug{
30namespace d3f{
31
33// Local discretization interface
35
37template<typename TDomain>
38void FractDDF_FV1<TDomain>::prepare_setting
39(
40 const std::vector<LFEID> & vLfeID,
41 bool bNonRegular
42)
43{
44// check the grid
45 if (bNonRegular)
46 UG_THROW ("ERROR in FractDDF_FV1:"
47 " The discretization does not support hanging nodes.\n");
48
49// check number of the components
50 if (vLfeID.size () != 2)
51 UG_THROW ("FractDDF_FV1: The density driven flow requires two coponents in every vertex.");
52
53// check whether these are the LagrangeP1 elements
54 if (vLfeID[0] != LFEID(LFEID::LAGRANGE, dim, 1) || vLfeID[1] != LFEID(LFEID::LAGRANGE, dim, 1))
55 UG_THROW ("FractDDF_FV1: This discretization works with the LagrangeP1-elements only.");
56}
57
59// assembling
61
63template<typename TDomain>
64void FractDDF_FV1<TDomain>::prep_assemble_loop ()
65{
66 if (m_spVolStabData.valid())
67 // set the "volume extensions" to 0
68 m_spVolStabData->reset_flux();
69
70 if (m_sss_mngr.valid ())
71 {
72 // reset the markers of the singular sources and sinks
73 m_sss_mngr->init_all_point_sss ();
74 m_sss_mngr->init_all_line_sss ();
75 }
76}
77
79template<typename TDomain>
80template<typename TElem>
81void FractDDF_FV1<TDomain>::prepare_element_loop
82(
83 ReferenceObjectID roid,
84 int si
85)
86{
87// check whether we are in a degenerated fracture
88 if (! m_spFractManager.valid())
89 m_isFracture = false;
90 else
91 {
92 if (! m_spFractManager->is_closed ())
93 UG_THROW ("FractDDF_FV1: Fracture manager not closed");
94 m_isFracture = m_spFractManager->contains (si);
95 }
96
97// check the imports
98 if (!m_imDensityCo.data_given())
99 UG_THROW ("FractDDF_FV1: Missing Import 'density at corners'.");
100 if (!m_imDensityIP.data_given())
101 UG_THROW ("FractDDF_FV1: Missing Import 'density at IPs'.");
102 if (!m_imViscosityCo.data_given())
103 UG_THROW ("FractDDF_FV1: Missing Import 'viscosity at corners'.");
104 if (!m_imViscosityIP.data_given())
105 UG_THROW ("FractDDF_FV1: Missing Import 'viscosity at IP's'.");
106 if (!m_imConstGravity.data_given())
107 UG_THROW ("FractDDF_FV1: Missing Import 'gravity'.");
108 if (!m_imPorosity.data_given())
109 UG_THROW ("FractDDF_FV1: Missing Import 'porosity'.");
110 if(m_spVolStabData.valid() && !m_imOldDensityCo.data_given())
111 UG_THROW ("FractDDF_FV1: Missing Import 'density at corners'.");
112
113// get the gravity
114 if (m_imConstGravity.constant())
115 (*m_imConstGravity.user_data()) (m_Gravity, MathVector<dim>(), 0.0, 0);
116 else
117 UG_THROW ("FractDDF_FV1: Gravity must be constant.");
118
119// call the specific routine:
120 if (! m_isFracture)
121 this->template bulk_prepare_element_loop<TElem> (roid, si);
122 else
123 this->template fract_prepare_element_loop<TElem> (roid, si);
124
125// check if the volume stabilization is active
126 m_bVolStabDataActive = (m_spVolStabData.valid() && m_spVolStabData->flux_is_active());
127}
128
130template<typename TDomain>
131template<typename TElem>
132void FractDDF_FV1<TDomain>::bulk_prepare_element_loop
133(
134 ReferenceObjectID roid,
135 int si
136)
137{
138 typedef FV1Geometry<TElem, dim> TBulkFVGeom;
139
140// check the imports
141 if (!m_imPermeability.data_given())
142 UG_THROW ("FractDDF_FV1: Missing Import 'full-dim. permeability'.");
143 if (!m_imDiffusion.data_given())
144 UG_THROW ("FractDDF_FV1: Missing Import 'full_dim. diffusion'.");
145
146// check, that upwind has been set
147 if (m_spUpwind.invalid())
148 UG_THROW("FractDDF_FV1: Upwind has not been set.");
149
150 TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
151 static const int refDim = TElem::dim;
152
153// connect to the own export velocity (not implemented yet for the bulk medium)
154 m_imFractDarcyVelIP.set_data(SPNULL);
155
156// set local IP coordinates
157 const MathVector<refDim>* vSCVFip = bulk_geo.scvf_local_ips();
158 size_t numSCVFip = bulk_geo.num_scvf_ips();
159 m_imDensityIP.template set_local_ips<refDim> (vSCVFip, numSCVFip, false);
160 m_imViscosityIP.template set_local_ips<refDim> (vSCVFip, numSCVFip, false);
161
162 const MathVector<refDim>* vSCVip = bulk_geo.scv_local_ips();
163 size_t numSCVip = bulk_geo.num_scv_ips();
164 m_imDensityCo.template set_local_ips<refDim> (vSCVip, numSCVip, false);
165 m_imViscosityCo.template set_local_ips<refDim> (vSCVip, numSCVip, false);
166 if (m_spVolStabData.valid ())
167 m_imOldDensityCo.template set_local_ips<refDim> (vSCVip, numSCVip, 1, false);
168
169 const MathVector<refDim>* coe_local = bulk_geo.coe_local ();
170 m_imPorosity.template set_local_ips<refDim> (coe_local, 1, false);
171 m_imPermeability.template set_local_ips<refDim> (coe_local, 1, false);
172 m_imDiffusion.template set_local_ips<refDim> (coe_local, 1, false);
173
174// init upwind for element type
175 if(! m_spUpwind->template set_geometry_type<TBulkFVGeom>(bulk_geo))
176 UG_THROW("FractDDF_FV1: Cannot init upwind for bulk element type.");
177}
178
180template<typename TDomain>
181template<typename TElem>
182void FractDDF_FV1<TDomain>::fract_prepare_element_loop
183(
184 ReferenceObjectID roid,
185 int si
186)
187{
188 typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
189 static const int refDim = TElem::dim; // dimensionality of the element, not the side!
190 ref_elem_type& rRefElem = Provider<ref_elem_type>::get ();
191
192// initialize the pointer to the FV geometry for fracture elements
193 m_pFractGeo = &GeomProvider<TFractFVGeom>::get (LFEID (LFEID::LAGRANGE, low_dim, 1), 1);
194
195// connect to the own export fracture velocity
196 m_imFractDarcyVelIP.set_data(m_exFractDarcyVel);
197
198// check the imports
199 if (!m_imAperture.data_given())
200 UG_THROW ("FractDDF_FV1: Missing Import 'fracture width (aperture)'.");
201 if (!m_imFractPermeability.data_given())
202 UG_THROW ("FractDDF_FV1: Missing Import 'permeability along fracture'.");
203 if (!m_imPermeability.data_given() && !m_imOrthoPermeability.data_given())
204 UG_THROW ("FractDDF_FV1: Missing Import 'permeability' (for the fract.-bulk interface interaction).");
205 if (!m_imFractDiffusion.data_given())
206 UG_THROW ("FractDDF_FV1: Missing Import 'diffusion along fracture'.");
207 if (!m_imDiffusion.data_given() && !m_imOrthoDiffusion.data_given())
208 UG_THROW ("FractDDF_FV1: Missing Import 'diffusion' (for the fract.-bulk interface interaction).");
209
210// check, that upwind has been set
211 if (m_spUpwind.invalid())
212 UG_THROW("FractDDF_FV1: Upwind has not been set.");
213
214// set up local ip coordinates for corner import parameters
215// REMARK: Note that for the fracture elements, values of the corner import
216// parameters are indexed not by scv (as for the normal elements) but by
217// the corner indices in the reference element
218 m_imDensityCo.template set_local_ips<refDim> (rRefElem.corners(), ref_elem_type::numCorners, false);
219 m_imViscosityCo.template set_local_ips<refDim> (rRefElem.corners(), ref_elem_type::numCorners, false);
220 if(m_spVolStabData.valid())
221 m_imOldDensityCo.template set_local_ips<refDim> (rRefElem.corners(), ref_elem_type::numCorners, 1, false);
222}
223
225template<typename TDomain>
226template<typename TElem>
230
232template<typename TDomain>
233template<typename TElem>
234void FractDDF_FV1<TDomain>::prepare_element
235(
236 const LocalVector & u,
237 GridObject * elem,
238 ReferenceObjectID roid, // id of reference element used for assembling
239 const position_type vCornerCoords []
240)
241{
242 TElem * pElem = static_cast<TElem*> (elem);
243
244// call the specific routine:
245 if (! m_isFracture)
246 this->template bulk_prepare_element<TElem> (u, pElem, vCornerCoords);
247 else
248 this->template fract_prepare_element<TElem> (u, pElem, vCornerCoords);
249}
250
252template<typename TDomain>
253template<typename TElem>
254void FractDDF_FV1<TDomain>::bulk_prepare_element
255(
256 const LocalVector & u,
257 TElem * elem,
258 const position_type vCornerCoords []
259)
260{
261 typedef FV1Geometry<TElem, dim> TBulkFVGeom;
262
263// update the FV geometry for this element
264 TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
265 try
266 {
267 bulk_geo.update (elem, vCornerCoords, &(this->subset_handler()));
268 }
269 UG_CATCH_THROW("FractDDF_FV1: Cannot update the Finite Volume Geometry for a bulk element.");
270
271// set global positions for user data
272 const MathVector<dim>* vSCVFip = bulk_geo.scvf_global_ips();
273 size_t numSCVFip = bulk_geo.num_scvf_ips();
274 m_imDensityIP.set_global_ips (vSCVFip, numSCVFip);
275 m_imViscosityIP.set_global_ips (vSCVFip, numSCVFip);
276
277 const MathVector<dim>* vSCVip = bulk_geo.scv_global_ips();
278 size_t numSCVip = bulk_geo.num_scv_ips();
279 m_imDensityCo.set_global_ips (vSCVip, numSCVip);
280 m_imViscosityCo.set_global_ips (vSCVip, numSCVip);
281 if(m_spVolStabData.valid())
282 m_imOldDensityCo.set_global_ips (vSCVip, numSCVip);
283
284 const MathVector<dim>* coe_global = bulk_geo.coe_global ();
285 m_imPorosity.set_global_ips (coe_global, 1);
286 m_imPermeability.set_global_ips (coe_global, 1);
287 m_imDiffusion.set_global_ips (coe_global, 1);
288}
289
291template<typename TDomain>
292template<typename TElem>
293void FractDDF_FV1<TDomain>::fract_prepare_element
294(
295 const LocalVector & u,
296 TElem * elem,
297 const position_type vCornerCoords []
298)
299{
300 typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
301 ref_elem_type& rRefElem = Provider<ref_elem_type>::get ();
302
303// get the non-degenerated sides of the fracture element
304 try
305 {
306 m_spFractManager->get_layer_sides
307 (elem,
308 m_numFractCo, m_innerFractSide, m_innerFractSideIdx, m_innerSideCo,
309 m_outerFractSide, m_outerFractSideIdx, m_outerSideCo,
310 m_assCo);
311 }
312 UG_CATCH_THROW("FractDDF_FV1: Cannot find orientation of a fracture element.");
313
314// compute the FV geometry of the inner side
315 MathVector<dim> vSideCornerCoords [maxFractSideCorners];
316 try
317 {
318 for (size_t co = 0; co < m_numFractCo; co++)
319 vSideCornerCoords [co] = vCornerCoords [m_innerSideCo [co]];
320 m_pFractGeo->update (m_innerFractSide, vSideCornerCoords, &(this->subset_handler()));
321 }
322 UG_CATCH_THROW("FractDDF_FV1: Cannot update the Finite Volume Geometry for a fracture element.");
323 size_t numSCVFip = m_pFractGeo->num_scvf_ips ();
324
325// convert local coordinates of the side into the local coordinates of the element (for the input parameters)
326 MathVector<dim> vSideLocCornerCoords [maxFractSideCorners];
327 try
328 {
329 DimReferenceMapping<low_dim, dim>& rMapping
330 = ReferenceMappingProvider::get<low_dim, dim> (m_innerFractSide->reference_object_id ());
331 for (size_t co = 0; co < m_numFractCo; co++)
332 vSideLocCornerCoords [co] = rRefElem.corner (m_innerSideCo [co]);
333 rMapping.update (vSideLocCornerCoords);
334
335 rMapping.local_to_global (m_elem_loc_coe, *(m_pFractGeo->coe_local ()));
336 rMapping.local_to_global (m_elem_loc_scvf, m_pFractGeo->scvf_local_ips (), numSCVFip);
337 }
338 UG_CATCH_THROW("FractDDF_FV1: Cannot transform local side coordinates to local element coordinates in a fracture element.");
339
340// set local IP coordinates
341 m_imDensityIP.template set_local_ips<dim> (m_elem_loc_scvf, numSCVFip);
342 m_imViscosityIP.template set_local_ips<dim> (m_elem_loc_scvf, numSCVFip);
343
344 m_imAperture.template set_local_ips<dim> (&m_elem_loc_coe, 1);
345 m_imPorosity.template set_local_ips<dim> (&m_elem_loc_coe, 1);
346 m_imFractPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
347 if (m_imOrthoPermeability.data_given ())
348 m_imOrthoPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
349 else
350 m_imPermeability.template set_local_ips<dim> (&m_elem_loc_coe, 1);
351 m_imFractDiffusion.template set_local_ips<dim> (&m_elem_loc_coe, 1);
352 if (m_imOrthoDiffusion.data_given ())
353 m_imOrthoDiffusion.template set_local_ips<dim> (&m_elem_loc_coe, 1);
354 else
355 m_imDiffusion.template set_local_ips<dim> (&m_elem_loc_coe, 1);
356
357// set global positions for user data
358 const MathVector<dim>* vSCVFip = m_pFractGeo->scvf_global_ips();
359 m_imDensityIP.set_global_ips (vSCVFip, numSCVFip);
360 m_imViscosityIP.set_global_ips (vSCVFip, numSCVFip);
361
362 const MathVector<dim>* coe_global = m_pFractGeo->coe_global ();
363 m_imAperture.set_global_ips (coe_global, 1);
364 m_imPorosity.set_global_ips (coe_global, 1);
365 m_imFractPermeability.set_global_ips (coe_global, 1);
366 if (m_imOrthoPermeability.data_given ())
367 m_imOrthoPermeability.set_global_ips (coe_global, 1);
368 else
369 m_imPermeability.set_global_ips (coe_global, 1);
370 m_imFractDiffusion.set_global_ips (coe_global, 1);
371 if (m_imOrthoDiffusion.data_given ())
372 m_imOrthoDiffusion.set_global_ips (coe_global, 1);
373 else
374 m_imDiffusion.set_global_ips (coe_global, 1);
375
376// set the correct subdomain for the bulk-fracture permeability and diffusion interface
377//TODO: This would be a good idea to get the (orthogonal) permeability and diffusion from the interface
378// if (! m_imPermeability.data()->constant ()) // constant data are the same everywhere!
379// m_imPermeability.data()->set_subset
380// (this->subset_handler().get_subset_index (m_outerFractSide));
381//END TODO
382
383// set up global ip coordinates for corner import parameters
384// REMARK: Note that for the fracture elements, values of the corner import
385// parameters are indexed not by scv (as for the normal elements) but by
386// the corner indices in the reference element
387 m_imDensityCo.set_global_ips (vCornerCoords, ref_elem_type::numCorners);
388 m_imViscosityCo.set_global_ips (vCornerCoords, ref_elem_type::numCorners);
389 if(m_spVolStabData.valid())
390 m_imOldDensityCo.set_global_ips (vCornerCoords, ref_elem_type::numCorners);
391
392// set global and local positions for velocity import
393 m_imFractDarcyVelIP.template set_local_ips<dim> (m_elem_loc_scvf, numSCVFip);
394 m_imFractDarcyVelIP.set_global_ips(vSCVFip, numSCVFip);
395
396// compute the outer normal for the fracture side and the orthogonal gravity
397 SideNormal<ref_elem_type, dim> (m_unitOuterNormal, m_outerFractSideIdx, vCornerCoords);
398 const number outerNormalNorm = VecLength (m_unitOuterNormal);
399 if (outerNormalNorm < 1e-32)
400 UG_THROW ("FractDDF_FV1: Cannot get the normal to a fracture.")
401 m_unitOuterNormal /= outerNormalNorm;
402 m_orthGravity = VecDot (m_Gravity, m_unitOuterNormal);
403
404// init upwind for element type
405 if(! m_spUpwind->template set_geometry_type<TFractFVGeom>(*m_pFractGeo))
406 UG_THROW("FractDDF_FV1: Cannot init upwind for fracture element type.");
407}
408
410template<typename TDomain>
411template<typename TFVGeom, typename TConsGravity>
412void FractDDF_FV1<TDomain>::compute_ip_Darcy_velocity
413(
414 MathVector<dim>& Vel,
415 size_t ip,
416 const TFVGeom& geo,
417 TConsGravity& ConsGravityMethod,
418 MathVector<TFVGeom::dim> vConsGravity[],
419 number vPressure[],
420 number Viscosity
421)
422{
423 const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
424
425// Compute rho * g (as the consistent gravity force)
426 ConsGravityMethod.template compute<dim>
427 (Vel, scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(), vConsGravity);
428
429// The pressure-gradient part:
430 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
431 VecScaleAppend (Vel, - vPressure[sh], scvf.global_grad(sh));
432
433// The viscosity factor:
434 VecScale (Vel, Vel, 1 / Viscosity);
435}
436
438template<typename TDomain>
439template<typename TElem>
440void FractDDF_FV1<TDomain>::ass_dA_elem
441(
442 LocalVector & d,
443 const LocalVector & u,
444 GridObject * elem,
445 const position_type vCornerCoords []
446)
447{
448 TElem * pElem = static_cast<TElem*> (elem);
449
450// call the specific routine:
451 if (! m_isFracture)
452 {
453 this->template bulk_ass_dA_elem<TElem> (d, u, pElem, vCornerCoords);
454 }
455 else
456 {
457 this->template fract_ass_dA_elem<TElem> (d, u, pElem, vCornerCoords);
458 this->template fract_bulk_ass_dA_elem<TElem> (d, u, pElem, vCornerCoords);
459 }
460}
461
463template<typename TDomain>
464template<typename TElem>
465void FractDDF_FV1<TDomain>::bulk_ass_dA_elem
466(
467 LocalVector & d,
468 const LocalVector & u,
469 TElem * pElem,
470 const position_type vCornerCoords []
471)
472{
473 typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
474 typedef FV1Geometry<TElem, dim> TBulkFVGeom;
475 typedef StdLinConsistentGravityX<ref_elem_type::dim> TConsGravity;
476
477 MathMatrix<dim,dim> Diffusion[TBulkFVGeom::numSCVF];
478 MathVector<dim> Dgrad_c_ip;
479 MathVector<dim> Vel[TBulkFVGeom::numSCVF], DarcyVel[TBulkFVGeom::numSCVF];
480 MathVector<dim> grad_c_ip;
481
482// Get finite volume geometry
483 const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
484
485 const size_t numSh = bulk_geo.num_sh();
486 const size_t numScvf = bulk_geo.num_scvf();
487
488// Consistent gravity and its derivative at corners
489 TConsGravity ConsGravityMethod;
490 MathVector<ref_elem_type::dim> vConsGravity [ref_elem_type::numCorners];
491
492// Prepare the consistent gravity
493 try
494 {
495 ConsGravityMethod.template prepare<dim>
496 (vConsGravity, numSh, vCornerCoords, m_imDensityCo.values(), m_Gravity);
497 }
498 UG_CATCH_THROW ("FractDDF_FV1::bulk_ass_dA_elem: Cannot prepare Consistent Gravity.");
499
500// Get the corner pressure
501 number vPressure [ref_elem_type::numCorners];
502 for (size_t sh = 0; sh < numSh; sh++)
503 vPressure[sh] = u(_P_, sh);
504
505// Compute Diffusion - Dispersion
506 MatScale(Diffusion[0], m_imPorosity[0], m_imDiffusion[0]);
507 for (size_t ip = 1; ip < numScvf; ip++)
508 Diffusion[ip] = Diffusion[0];
509 //TODO add Dispersion
510
511// Compute the Darcy velocity at the ips and the convection shapes
512 for (size_t ip = 0; ip < numScvf; ip++)
513 {
514 this->template compute_ip_Darcy_velocity<TBulkFVGeom, TConsGravity>
515 (Vel[ip], ip, bulk_geo, ConsGravityMethod, vConsGravity, vPressure,
516 m_imViscosityIP[ip]);
517 MatVecMult (DarcyVel[ip], m_imPermeability[0], Vel[ip]);
518 }
519
520// compute upwind shapes for transport equation
521 if(!m_spUpwind->update(&bulk_geo, DarcyVel, Diffusion, false))
522 UG_THROW("FractDDF_FV1::bulk_ass_dA_elem: Cannot compute convection shapes.");
523
524// get a const (!!) reference to the upwind
525 const IConvectionShapes<dim>& convShape
526 = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
527
528// Loop Sub Control Volume Faces (SCVF)
529 for (size_t ip = 0; ip < numScvf; ip++)
530 {
531 // Get current SCVF
532 const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf(ip);
533
534 // The ip-concentration and the gradient of c at the ip:
535 VecSet (grad_c_ip, 0.0);
536 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
537 VecScaleAppend (grad_c_ip, u(_C_, sh), scvf.global_grad(sh));
538
540 // Transport Equation
542
543 number flux = 0;
544
545 // Compute Convective Flux
546 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
547 flux += convShape(ip, sh) * u(_C_, sh);
548
549 // Compute Diffusive Flux
550 MatVecMult(Dgrad_c_ip, Diffusion[ip], grad_c_ip);
551 const number diffFlux = VecDot(Dgrad_c_ip, scvf.normal());
552
553 // Sum total flux
554 flux -= diffFlux;
555 if(!m_BoussinesqTransport) flux *= m_imDensityIP[ip];
556
557 // Add contribution to transport equation
558 d(_C_,scvf.from()) += flux;
559 d(_C_,scvf.to()) -= flux;
560
562 // Flow Equation
564
565 // Compute flux
566 flux = VecDot(DarcyVel[ip], scvf.normal());
567 if(!m_BoussinesqFlow) flux *= m_imDensityIP[ip];
568
569 // Add contribution to flow equation
570 d(_P_,scvf.from()) += flux;
571 d(_P_,scvf.to()) -= flux;
572
573 // Volume stabilization: Add the continuity equation to the transport equation
574 if(m_bVolStabDataActive)
575 {
576 m_spVolStabData->stiff(pElem->vertex(scvf.from())) -= flux;
577 m_spVolStabData->stiff(pElem->vertex(scvf.to())) += flux;
578 //d(_C_,scvf.from()) -= flux * u(_C_,scvf.from());
579 //d(_C_,scvf.to()) += flux * u(_C_,scvf.to());
580 }
581 }
582
583// Assemble the singular sources and sinks
584 if (m_sss_mngr.valid () && (m_sss_mngr->num_points () != 0 || m_sss_mngr->num_lines () != 0))
585 {
586 typedef typename domain_type::position_accessor_type t_pos_accessor;
587 typedef typename DDSingularSourcesAndSinks<dim>::template
588 point_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_pnt_sss_iter;
589 typedef typename DDSingularSourcesAndSinks<dim>::template
590 line_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_lin_sss_iter;
591
592 t_pos_accessor& aaPos = this->domain()->position_accessor ();
593 Grid& grid = (Grid&) *this->domain()->grid ();
594
595 for(size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
596 {
597 size_t co = bulk_geo.scv(ip).node_id ();
598
599 // point sources
600 for (t_pnt_sss_iter pnt (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
601 ! pnt.is_over (); ++pnt)
602 {
603 FVPointSourceOrSink<dim, dd_point_sss_data<dim> > * pnt_sss = *pnt;
604 if (! pnt_sss->marked_for (pElem, co))
605 continue;
606 pnt_sss->compute (pnt_sss->position (), this->time (), -1); //TODO: set the subset id instead of -1
607 ass_sss_dA_elem (d, u, pElem, co, pnt_sss->intensity (), pnt_sss->concentration ());
608 }
609
610 // line sources
611 for (t_lin_sss_iter line (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
612 ! line.is_over (); ++line)
613 {
614 FVLineSourceOrSink<dim, dd_line_sss_data<dim> > * line_sss = *line;
615 number len = VecDistance (line.seg_start (), line.seg_end ());
616 line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
617 ass_sss_dA_elem (d, u, pElem, co, line_sss->intensity () * len, line_sss->concentration ());
618 }
619 }
620 }
621}
622
624template<typename TDomain>
625template<typename TElem>
626void FractDDF_FV1<TDomain>::fract_ass_dA_elem
627(
628 LocalVector & d,
629 const LocalVector & u,
630 TElem * pElem,
631 const position_type vCornerCoords []
632)
633{
634 const number half_fr_width = m_imAperture[0] / 2;
635
636 MathMatrix<dim,dim> Diffusion[TFractFVGeom::maxNumSCVF];
637 MathVector<dim> Dgrad_c_ip;
638 MathVector<dim> grad_c_ip;
639
640 const size_t numSh = m_pFractGeo->num_sh();
641 const size_t numScvf = m_pFractGeo->num_scvf();
642
643// Get the corner values
644 number vConcentration [maxFractSideCorners];
645 for (size_t sh = 0; sh < numSh; sh++)
646 {
647 size_t co = m_innerSideCo[sh];
648 vConcentration[sh] = u(_C_, co);
649 }
650
651// Compute Diffusion - Dispersion
652 MatSet (Diffusion[0], 0);
653 MatDiagSet (Diffusion[0], m_imPorosity[0] * m_imFractDiffusion[0]);
654 for (size_t ip = 1; ip < numScvf; ip++)
655 Diffusion[ip] = Diffusion[0];
656 //TODO add Dispersion
657
658// compute upwind shapes for transport equation
659 if(!m_spUpwind->update(m_pFractGeo, m_imFractDarcyVelIP.values(), Diffusion, false))
660 UG_THROW("FractDDF_FV1::fract_ass_dA_elem: Cannot compute convection shapes.");
661
662// get a const (!!) reference to the upwind
663 const IConvectionShapes<dim>& convShape
664 = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
665
666// Loop Sub Control Volume Faces (SCVF)
667 for (size_t ip = 0; ip < numScvf; ip++)
668 {
669 // Get current SCVF
670 const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf(ip);
671
672 // The ip-concentration and the gradient of c at the ip:
673 VecSet (grad_c_ip, 0.0);
674 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
675 VecScaleAppend (grad_c_ip, vConcentration[sh], scvf.global_grad(sh));
676
678 // Transport Equation
680
681 number flux = 0;
682
683 // Compute Convective Flux
684 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
685 flux += convShape(ip, sh) * vConcentration[sh];
686
687 // Compute Diffusive Flux
688 MatVecMult(Dgrad_c_ip, Diffusion[ip], grad_c_ip);
689 const number diffFlux = VecDot(Dgrad_c_ip, scvf.normal());
690
691 // Sum total flux
692 flux = (flux - diffFlux) * half_fr_width;
693 if(!m_BoussinesqTransport) flux *= m_imDensityIP[ip];
694
695 // Add contribution to transport equation
696 d(_C_, m_innerSideCo[scvf.from()]) += flux;
697 d(_C_, m_innerSideCo[scvf.to()]) -= flux;
698
700 // Flow Equation
702
703 // Compute flux
704 flux = VecDot(m_imFractDarcyVelIP[ip], scvf.normal()) * half_fr_width;
705
706 if(!m_BoussinesqFlow) flux *= m_imDensityIP[ip];
707
708 // Add contribution to flow equation
709 d(_P_, m_innerSideCo[scvf.from()]) += flux;
710 d(_P_, m_innerSideCo[scvf.to()]) -= flux;
711
712 // Volume stabilization: Add the continuity equation to the transport equation
713 if(m_bVolStabDataActive)
714 {
715 m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.from()])) -= flux;
716 m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.to()])) += flux;
717 //d(_C_,m_innerSideCo[scvf.from()]) -= flux * u(_C_,m_innerSideCo[scvf.from()]);
718 //d(_C_,m_innerSideCo[scvf.to()] ) += flux * u(_C_,m_innerSideCo[scvf.to()]);
719 }
720 }
721
722// Assemble the singular sources and sinks
723 if (m_sss_mngr.valid () && m_sss_mngr->num_lines () != 0)
724 {
725 typedef typename domain_type::position_accessor_type t_pos_accessor;
726 typedef typename DDSingularSourcesAndSinks<dim>::template
727 line_iterator<side_type,t_pos_accessor,TFractFVGeom> t_lin_sss_iter;
728
729 t_pos_accessor& aaPos = this->domain()->position_accessor ();
730 Grid& grid = (Grid&) *this->domain()->grid ();
731
732 for(size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
733 {
734 // Get the corner of the face
735 size_t side_co = m_pFractGeo->scv(ip).node_id ();
736 // Get associated node of the element (not side!)
737 size_t co = m_innerSideCo [m_pFractGeo->scv(ip).node_id ()];
738
739 // line sources (that correspond to the point sources)
740 for (t_lin_sss_iter line (m_sss_mngr.get (), m_innerFractSide, grid, aaPos, *m_pFractGeo, side_co);
741 ! line.is_over (); ++line)
742 {
743 FVLineSourceOrSink<dim, dd_line_sss_data<dim> > * line_sss = *line;
744 if (! line_sss->marked_for (m_innerFractSide, side_co))
745 continue;
746 line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
747 ass_sss_dA_elem (d, u, pElem, co, line_sss->intensity () / 2, line_sss->concentration ());
748 /* Remark: "/ 2" because the source is taken into account twise. */
749 }
750 }
751 }
752}
753
755template<typename TDomain>
756template<typename TElem>
757void FractDDF_FV1<TDomain>::fract_bulk_ass_dA_elem
758(
759 LocalVector & d,
760 const LocalVector & u,
761 TElem * pElem,
762 const position_type vCornerCoords []
763)
764{
765// the permeability of the fracture-bulk interface
766 number orthPermeability;
767 if (m_imOrthoPermeability.data_given ())
768 orthPermeability = m_imOrthoPermeability[0];
769 else
770 {
771 MathVector<dim> orthPermeabilityVec;
772 MatVecMult (orthPermeabilityVec, m_imPermeability[0], m_unitOuterNormal);
773 orthPermeability = VecDot (orthPermeabilityVec, m_unitOuterNormal);
774 }
775
776// the diffusion of the fracture-bulk interface
777 number orthDiffusion;
778 if (m_imOrthoDiffusion.data_given())
779 orthDiffusion = m_imPorosity[0] * m_imOrthoDiffusion[0];
780 else
781 {
782 MathVector<dim> orthDiffusionVec;
783 MatVecMult (orthDiffusionVec, m_imDiffusion[0], m_unitOuterNormal);
784 orthDiffusion = m_imPorosity[0] * VecDot (orthDiffusionVec, m_unitOuterNormal);
785 }
786
787// the "height" of the degenerated element
788 const number half_fr_width = m_imAperture[0] / 2;
789
790// loop over the corners of the inner side
791 for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
792 {
793 // Get current SCV
794 const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
795
796 // Get associated node of the element (not side!)
797 const int co = m_innerSideCo [scv.node_id()];
798
799 // Get the corner values
800 const number orthC_f = u(_C_, co);
801 const number orthC_m = u(_C_, m_assCo[co]);
802 const number orthP_f = u(_P_, co);
803 const number orthP_m = u(_P_, m_assCo[co]);
804 const number fractDensity = m_imDensityCo[co];
805 const number orthDensity = m_imDensityCo [m_assCo[co]];
806 const number orthViscosity = m_imViscosityCo [m_assCo[co]];
807
808 // Flux velocity between the fracture and the bulk medium
809 const number orthVelocity = (m_orthGravity * (orthDensity - fractDensity)
810 - (orthP_m - orthP_f) / half_fr_width)
811 * orthPermeability / orthViscosity;
812
814 number flux;
815
817 // Transport Equation
819
820 flux = orthDiffusion * (orthC_m - orthC_f) / half_fr_width;
821 /* We use the full upwind here: */
822 flux -= orthVelocity * ((orthVelocity >= 0)? orthC_f : orthC_m);
823 flux *= scv.volume();
824 if (! m_BoussinesqTransport) flux *= orthDensity;
825 d(_C_, m_assCo[co]) += flux;
826 d(_C_, co) -= flux;
827
829 // Flow Equation
831
832 flux = orthVelocity * scv.volume();
833 if (! m_BoussinesqFlow) flux *= orthDensity;
834 d(_P_, m_assCo[co]) -= flux;
835 d(_P_, co) += flux;
836
837 // Volume stabilization: Add the continuity equation to the transport equation
838 if(m_bVolStabDataActive)
839 {
840 m_spVolStabData->stiff(pElem->vertex(m_assCo[co])) += flux;
841 m_spVolStabData->stiff(pElem->vertex(co)) -= flux;
842 //d(_C_,m_assCo[co]) += flux * u(_C_,m_assCo[co]);
843 //d(_C_,co ) -= flux * u(_C_,co);
844 }
845 }
846}
847
849template<typename TDomain>
850template<typename TFVGeom, typename TConsGravity, size_t maxCorners>
851void FractDDF_FV1<TDomain>::compute_J_ip_Darcy_velocity
852(
853 MathVector<dim>& Vel,
854 MathVector<dim> Vel_c[],
855 MathVector<dim> Vel_p[],
856 size_t ip,
857 const TFVGeom& geo,
858 TConsGravity& ConsGravityMethod,
859 MathVector<TFVGeom::dim> vConsGravity_c[][maxCorners],
860 number vPressure[],
861 number Viscosity,
862 number Viscosity_c[]
863)
864{
865 const typename TFVGeom::SCVF& scvf = geo.scvf(ip);
866 const number InvVisco = 1 / Viscosity;
867 const size_t numSh = geo.num_sh();
868
869// Compute the derivatives of rho g - grad p, multiplied by the viscosity
870 for (size_t sh = 0; sh < numSh; sh++)
871 {
872 ConsGravityMethod.template compute<dim>
873 (Vel_c[sh], scvf.local_ip(), scvf.JTInv(), scvf.local_grad_vector(),
874 vConsGravity_c[sh]);
875 VecScale (Vel_c[sh], Vel_c[sh], InvVisco);
876
877 VecScale (Vel_p[sh], scvf.global_grad(sh), -InvVisco);
878 }
879
880// Take into account the derivative of the viscosity
881 if (Viscosity_c != NULL)
882 for (size_t sh = 0; sh < numSh; sh++)
883 VecScaleAppend (Vel_c[sh], -Viscosity_c[sh] * InvVisco, Vel);
884}
885
887template<typename TDomain>
888template<typename TElem>
889void FractDDF_FV1<TDomain>::ass_JA_elem
890(
891 LocalMatrix & J,
892 const LocalVector & u,
893 GridObject * elem,
894 const position_type vCornerCoords []
895)
896{
897 TElem * pElem = static_cast<TElem*> (elem);
898
899// call the specific routine:
900 if (! m_isFracture)
901 this->template bulk_ass_JA_elem<TElem> (J, u, pElem, vCornerCoords);
902 else
903 {
904 this->template fract_ass_JA_elem<TElem> (J, u, pElem, vCornerCoords);
905 this->template fract_bulk_ass_JA_elem<TElem> (J, u, pElem, vCornerCoords);
906 }
907}
908
910template<typename TDomain>
911template<typename TElem>
912void FractDDF_FV1<TDomain>::bulk_ass_JA_elem
913(
914 LocalMatrix & J,
915 const LocalVector & u,
916 TElem * pElem,
917 const position_type vCornerCoords []
918)
919{
920 typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
921 typedef FV1Geometry<TElem, dim> TBulkFVGeom;
922 typedef StdLinConsistentGravityX<ref_elem_type::dim> TConsGravity;
923 static const size_t numCorners = ref_elem_type::numCorners;
924
925 // Get finite volume geometry
926 const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
927
928 const size_t numSh = bulk_geo.num_sh();
929 const size_t numScvf = bulk_geo.num_scvf();
930
931 MathMatrix<dim,dim> Diffusion[TBulkFVGeom::numSCVF];
932 MathVector<dim> Dgrad;
933 MathVector<dim> Vel[TBulkFVGeom::numSCVF], DarcyVel[TBulkFVGeom::numSCVF];
934 MathVector<dim> Vel_c[numCorners], Vel_p[numCorners];
935 MathVector<dim> vDDarcyVel_c[numCorners], vDDarcyVel_p[numCorners];
936 MathVector<dim> grad_c_ip;
937 number vDFlux_c [numCorners], vDFlux_p [numCorners];
938 number vViscosity_c [numCorners], * pViscosity_c;
939
940// Consistent gravity and its derivative at corners
941 TConsGravity ConsGravityMethod;
942 MathVector<ref_elem_type::dim> vConsGravity [numCorners];
943 MathVector<ref_elem_type::dim> vConsGravity_c [numCorners][numCorners];
944
945// Prepare the consistent gravity
946 try
947 {
948 ConsGravityMethod.template prepare<dim>
949 (vConsGravity, numSh, vCornerCoords, m_imDensityCo.values(), m_Gravity);
950
951 number co_density [numCorners];
952 memset (co_density, 0, ref_elem_type::numCorners * sizeof (number));
953 for (size_t sh = 0; sh < numSh; sh++)
954 { //TODO: use scv, not merely shape id here!
955 co_density[sh] = (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv (sh, _C_, sh)); // we assume the completely local dependence of \f$\rho\f$ of \f$c\f$
956 ConsGravityMethod.template prepare<dim>
957 (vConsGravity_c[sh], numSh, vCornerCoords, co_density, m_Gravity);
958 co_density[sh] = 0.0;
959 }
960 }
961 UG_CATCH_THROW ("FractDDF_FV1::bulk_ass_JA_elem: Cannot prepare Consistent Gravity or its derivatives.");
962
963// Get the corner pressure
964 number vPressure [numCorners];
965 for (size_t sh = 0; sh < numSh; sh++)
966 vPressure[sh] = u(_P_, sh);
967
968// Compute Diffusion - Dispersion
969 MatScale(Diffusion[0], m_imPorosity[0], m_imDiffusion[0]);
970 for (size_t ip = 1; ip < numScvf; ip++)
971 Diffusion[ip] = Diffusion[0];
972 //TODO add Dispersion
973
974// Compute the Darcy velocity at the ips and the convection shapes
975 for (size_t ip = 0; ip < numScvf; ip++)
976 {
977 this->template compute_ip_Darcy_velocity<TBulkFVGeom, TConsGravity>
978 (Vel[ip], ip, bulk_geo, ConsGravityMethod, vConsGravity, vPressure,
979 m_imViscosityIP[ip]);
980 MatVecMult (DarcyVel[ip], m_imPermeability[0], Vel[ip]);
981 }
982
983// compute upwind shapes for transport equation
984 if(!m_spUpwind->update(&bulk_geo, DarcyVel, Diffusion, true))
985 UG_THROW("FractDDF_FV1::bulk_ass_dA_elem: Cannot compute convection shapes.");
986
987// get a const (!!) reference to the upwind
988 const IConvectionShapes<dim>& convShape
989 = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
990
991// Loop Sub Control Volume Faces (SCVF)
992 for (size_t ip = 0; ip < numScvf; ip++)
993 {
994 // Get current SCVF
995 const typename TBulkFVGeom::SCVF& scvf = bulk_geo.scvf(ip);
996
997 // The ip-concentration and the gradient of c at the ip:
998 VecSet (grad_c_ip, 0.0);
999 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1000 VecScaleAppend (grad_c_ip, u(_C_, sh), scvf.global_grad(sh));
1001
1002 // Compute the derivative of the viscosity
1003 if (m_imViscosityIP.constant ())
1004 pViscosity_c = NULL;
1005 else
1006 {
1007 for (size_t sh = 0; sh < numSh; sh++)
1008 vViscosity_c[sh] = (m_imViscosityIP.constant() ? 0.0 : m_imViscosityIP.deriv(ip, _C_, sh));
1009 pViscosity_c = vViscosity_c;
1010 }
1011
1012 // Compute derivatives of the Darcy velocity
1013 this->template compute_J_ip_Darcy_velocity <TBulkFVGeom, TConsGravity, numCorners>
1014 (Vel[ip], Vel_c, Vel_p, ip, bulk_geo, ConsGravityMethod, vConsGravity_c, vPressure,
1015 m_imViscosityIP[ip], pViscosity_c);
1016 for (size_t sh = 0; sh < numSh; sh++)
1017 {
1018 MatVecMult (vDDarcyVel_c[sh], m_imPermeability[0], Vel_c[sh]);
1019 MatVecMult (vDDarcyVel_p[sh], m_imPermeability[0], Vel_p[sh]);
1020 }
1021
1023 // Transport Equation
1025
1026 // Loop Shape Functions
1027 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1028 {
1029 // Compute Derivative of Convective Flux
1030 vDFlux_c[sh] = convShape(ip, sh);
1031 vDFlux_p[sh] = 0.0;
1032
1033 // Derivative w.r.t. Velocity
1034 for(size_t sh1 = 0; sh1 < scvf.num_sh(); ++sh1)
1035 {
1036 vDFlux_c[sh] += u(_C_, sh1) * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_c[sh]);
1037 vDFlux_p[sh] += u(_C_, sh1) * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_p[sh]);
1038 }
1039
1040 //TODO: add derivative of Dispersion
1041
1042 // Add Derivative of Diffusive Flux
1043 MatVecMult(Dgrad, Diffusion[ip], scvf.global_grad(sh));
1044 vDFlux_c[sh] -= VecDot(Dgrad, scvf.normal());
1045 }
1046
1047 // Handle density in case of full equation
1048 if(!m_BoussinesqTransport)
1049 {
1050 // Convective Flux
1051 number flux = 0;
1052 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1053 flux += convShape(ip, sh) * u(_C_, sh);
1054
1055 // Diffusive Flux
1056 MatVecMult(Dgrad, Diffusion[ip], grad_c_ip);
1057 flux -= VecDot(Dgrad, scvf.normal());
1058
1059 // Derivative of product
1060 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1061 {
1062 vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
1063 (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, sh)) * flux;
1064 vDFlux_p[sh] *= m_imDensityIP[ip];
1065 }
1066 }
1067
1068 // Add Flux contribution
1069 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1070 {
1071 J(_C_, scvf.from(), _C_, sh) += vDFlux_c[sh];
1072 J(_C_, scvf.to(), _C_, sh) -= vDFlux_c[sh];
1073 J(_C_, scvf.from(), _P_, sh) += vDFlux_p[sh];
1074 J(_C_, scvf.to(), _P_, sh) -= vDFlux_p[sh];
1075 }
1076
1077
1079 // Flow Equation
1081
1082 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1083 {
1084 vDFlux_c[sh] = VecDot(vDDarcyVel_c[sh], scvf.normal());
1085 vDFlux_p[sh] = VecDot(vDDarcyVel_p[sh], scvf.normal());
1086 }
1087
1088 number flux = VecDot(DarcyVel[ip], scvf.normal());
1089
1090 if(!m_BoussinesqFlow)
1091 {
1092 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1093 {
1094 vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
1095 (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, sh)) * flux;
1096 vDFlux_p[sh] *= m_imDensityIP[ip];
1097 }
1098
1099 flux *= m_imDensityIP[ip];
1100 }
1101
1102 // Add Flux contribution
1103 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1104 {
1105 J(_P_, scvf.from(), _C_, sh) += vDFlux_c[sh];
1106 J(_P_, scvf.to(), _C_, sh) -= vDFlux_c[sh];
1107 J(_P_, scvf.from(), _P_, sh) += vDFlux_p[sh];
1108 J(_P_, scvf.to(), _P_, sh) -= vDFlux_p[sh];
1109 }
1110
1111 // Volume stabilization: Add the continuity equation to the transport equation
1112 // REMARK: We do not compute derivatives of the defect of the conf. eq. here because this defect tends to 0.
1113 if(m_bVolStabDataActive)
1114 {
1115 m_spVolStabData->stiff(pElem->vertex(scvf.from())) -= flux;
1116 m_spVolStabData->stiff(pElem->vertex(scvf.to())) += flux;
1117 //J(_C_, scvf.from(), _C_, scvf.from()) -= flux;
1118 //J(_C_, scvf.to(), _C_, scvf.to()) += flux;
1119 }
1120 }
1121
1122// Assemble the singular sources and sinks
1123 if (m_sss_mngr.valid () && (m_sss_mngr->num_points () != 0 || m_sss_mngr->num_lines () != 0))
1124 {
1125 typedef typename domain_type::position_accessor_type t_pos_accessor;
1126 typedef typename DDSingularSourcesAndSinks<dim>::template
1127 point_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_pnt_sss_iter;
1128 typedef typename DDSingularSourcesAndSinks<dim>::template
1129 line_iterator<TElem,t_pos_accessor,TBulkFVGeom> t_lin_sss_iter;
1130
1131 t_pos_accessor& aaPos = this->domain()->position_accessor ();
1132 Grid& grid = (Grid&) *this->domain()->grid ();
1133
1134 for(size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
1135 {
1136 size_t co = bulk_geo.scv(ip).node_id ();
1137
1138 // point sources
1139 for (t_pnt_sss_iter pnt (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
1140 ! pnt.is_over (); ++pnt)
1141 {
1142 FVPointSourceOrSink<dim, dd_point_sss_data<dim> > * pnt_sss = *pnt;
1143 if (! pnt_sss->marked_for (pElem, co))
1144 continue;
1145 pnt_sss->compute (pnt_sss->position (), this->time (), -1); //TODO: set the subset id instead of -1
1146 ass_sss_JA_elem (J, u, pElem, co, pnt_sss->intensity (), pnt_sss->concentration ());
1147 }
1148
1149 // line sources
1150 for (t_lin_sss_iter line (m_sss_mngr.get (), pElem, grid, aaPos, bulk_geo, co);
1151 ! line.is_over (); ++line)
1152 {
1153 FVLineSourceOrSink<dim, dd_line_sss_data<dim> > * line_sss = *line;
1154 number len = VecDistance (line.seg_start (), line.seg_end ());
1155 line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
1156 ass_sss_JA_elem (J, u, pElem, co, line_sss->intensity () * len, line_sss->concentration ());
1157 }
1158 }
1159 }
1160}
1161
1163template<typename TDomain>
1164template<typename TElem>
1165void FractDDF_FV1<TDomain>::fract_ass_JA_elem
1166(
1167 LocalMatrix & J,
1168 const LocalVector & u,
1169 TElem * pElem,
1170 const position_type vCornerCoords []
1171)
1172{
1173 number half_fr_width = m_imAperture[0] / 2;
1174
1175 const size_t numSh = m_pFractGeo->num_sh();
1176 const size_t numScvf = m_pFractGeo->num_scvf();
1177
1178 MathMatrix<dim,dim> Diffusion[TFractFVGeom::maxNumSCVF];
1179 MathVector<dim> Dgrad;
1180 MathVector<dim> grad_c_ip;
1181 number vDFlux_c [maxFractSideCorners], vDFlux_p [maxFractSideCorners];
1182
1183// Get the corner values
1184 number vConcentration [maxFractSideCorners];
1185 for (size_t sh = 0; sh < numSh; sh++)
1186 {
1187 size_t co = m_innerSideCo[sh]; // co: index of the corner in the (degenerated) full-dim. element
1188 vConcentration[sh] = u(_C_, co);
1189 }
1190
1191// Compute Diffusion - Dispersion
1192 MatSet (Diffusion[0], 0);
1193 MatDiagSet (Diffusion[0], m_imPorosity[0] * m_imFractDiffusion[0]);
1194 for (size_t ip = 1; ip < numScvf; ip++)
1195 Diffusion[ip] = Diffusion[0];
1196 //TODO add Dispersion
1197
1198// compute upwind shapes for transport equation
1199 if(!m_spUpwind->update(m_pFractGeo, m_imFractDarcyVelIP.values(), Diffusion, true))
1200 UG_THROW("FractDDF_FV1::fract_ass_dA_elem: Cannot compute convection shapes.");
1201
1202// get a const (!!) reference to the upwind
1203 const IConvectionShapes<dim>& convShape
1204 = *const_cast<const IConvectionShapes<dim>*>(m_spUpwind.get());
1205
1206// Loop Sub Control Volume Faces (SCVF)
1207 for (size_t ip = 0; ip < numScvf; ip++)
1208 {
1209 // Get current SCVF
1210 const typename TFractFVGeom::SCVF& scvf = m_pFractGeo->scvf(ip);
1211
1212 // The ip-concentration and the gradient of c at the ip:
1213 VecSet (grad_c_ip, 0.0);
1214 for (size_t sh = 0; sh < scvf.num_sh(); sh++)
1215 VecScaleAppend (grad_c_ip, vConcentration[sh], scvf.global_grad(sh));
1216
1217 const MathVector<dim>* vDDarcyVel_c = m_imFractDarcyVelIP.deriv(ip, _C_);
1218 const MathVector<dim>* vDDarcyVel_p = m_imFractDarcyVelIP.deriv(ip, _P_);
1219
1221 // Transport Equation
1223
1224 // Loop Shape Functions
1225 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1226 {
1227 size_t co = m_innerSideCo[sh]; // co: index of the corner in the (degenerated) full-dim. element
1228
1229 // Compute Derivative of Convective Flux
1230 vDFlux_c[sh] = convShape(ip, sh);
1231 vDFlux_p[sh] = 0.0;
1232
1233 // Derivative w.r.t. Velocity
1234 for(size_t sh1 = 0; sh1 < scvf.num_sh(); ++sh1)
1235 {
1236 vDFlux_c[sh] += vConcentration[sh1] * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_c[co]);
1237 vDFlux_p[sh] += vConcentration[sh1] * VecDot(convShape.D_vel(ip, sh1), vDDarcyVel_p[co]);
1238 }
1239 //TODO: add derivative of Dispersion
1240
1241 // Add Derivative of Diffusive Flux
1242 MatVecMult(Dgrad, Diffusion[ip], scvf.global_grad(sh));
1243 vDFlux_c[sh] -= VecDot(Dgrad, scvf.normal());
1244 }
1245
1246 // Handle density in case of full equation
1247 if(!m_BoussinesqTransport)
1248 {
1249 // Convective Flux
1250 number flux = 0;
1251 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1252 flux += convShape(ip, sh) * vConcentration[sh];
1253
1254 // Diffusive Flux
1255 MatVecMult(Dgrad, Diffusion[ip], grad_c_ip);
1256 flux -= VecDot(Dgrad, scvf.normal());
1257
1258 // Derivative of product
1259 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1260 {
1261 vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
1262 (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, m_innerSideCo[sh])) * flux;
1263 vDFlux_p[sh] *= m_imDensityIP[ip];
1264 }
1265 }
1266
1267 // Multiply by the fracture width
1268 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1269 {
1270 vDFlux_c[sh] *= half_fr_width; vDFlux_p[sh] *= half_fr_width;
1271 }
1272
1273 // Add Flux contribution
1274 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1275 {
1276 size_t co = m_innerSideCo[sh];
1277 size_t co_from = m_innerSideCo[scvf.from()];
1278 size_t co_to = m_innerSideCo[scvf.to()];
1279
1280 J(_C_, co_from, _C_, co) += vDFlux_c[sh];
1281 J(_C_, co_to, _C_, co) -= vDFlux_c[sh];
1282 J(_C_, co_from, _P_, co) += vDFlux_p[sh];
1283 J(_C_, co_to, _P_, co) -= vDFlux_p[sh];
1284 }
1285
1286
1288 // Flow Equation
1290
1291 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1292 {
1293 size_t co = m_innerSideCo[sh];
1294 vDFlux_c[sh] = VecDot(vDDarcyVel_c[co], scvf.normal());
1295 vDFlux_p[sh] = VecDot(vDDarcyVel_p[co], scvf.normal());
1296 }
1297
1298 number flux = VecDot(m_imFractDarcyVelIP[ip], scvf.normal());
1299
1300 if(!m_BoussinesqFlow)
1301 {
1302 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1303 {
1304 vDFlux_c[sh] = m_imDensityIP[ip] * vDFlux_c[sh] +
1305 (m_imDensityIP.constant() ? 0.0 : m_imDensityIP.deriv(ip, _C_, m_innerSideCo[sh])) * flux;
1306 vDFlux_p[sh] *= m_imDensityIP[ip];
1307 }
1308
1309 flux *= m_imDensityIP[ip];
1310 }
1311
1312 // Multiply by the fracture width
1313 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1314 {
1315 vDFlux_c[sh] *= half_fr_width; vDFlux_p[sh] *= half_fr_width;
1316 }
1317
1318 // Add Flux contribution
1319 for(size_t sh = 0; sh < scvf.num_sh(); ++sh)
1320 {
1321 size_t co = m_innerSideCo[sh];
1322 size_t co_from = m_innerSideCo[scvf.from()];
1323 size_t co_to = m_innerSideCo[scvf.to()];
1324
1325 J(_P_, co_from, _C_, co) += vDFlux_c[sh];
1326 J(_P_, co_to, _C_, co) -= vDFlux_c[sh];
1327 J(_P_, co_from, _P_, co) += vDFlux_p[sh];
1328 J(_P_, co_to, _P_, co) -= vDFlux_p[sh];
1329 }
1330
1331 // Volume stabilization: Add the continuity equation to the transport equation
1332 // REMARK: We do not compute derivatives of the defect of the conf. eq. here because this defect tends to 0.
1333 if(m_bVolStabDataActive)
1334 {
1335 flux *= half_fr_width;
1336 m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.from()])) -= flux;
1337 m_spVolStabData->stiff(pElem->vertex(m_innerSideCo[scvf.to()])) += flux;
1338 //J(_C_, m_innerSideCo[scvf.from()], _C_, m_innerSideCo[scvf.from()]) -= flux;
1339 //J(_C_, m_innerSideCo[scvf.to()], _C_, m_innerSideCo[scvf.to()] ) += flux;
1340 }
1341 }
1342
1343// Assemble the singular sources and sinks
1344 if (m_sss_mngr.valid () && m_sss_mngr->num_lines () != 0)
1345 {
1346 typedef typename domain_type::position_accessor_type t_pos_accessor;
1347 typedef typename DDSingularSourcesAndSinks<dim>::template
1348 line_iterator<side_type,t_pos_accessor,TFractFVGeom> t_lin_sss_iter;
1349
1350 t_pos_accessor& aaPos = this->domain()->position_accessor ();
1351 Grid& grid = (Grid&) *this->domain()->grid ();
1352
1353 for(size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1354 {
1355 // Get the corner of the face
1356 size_t side_co = m_pFractGeo->scv(ip).node_id ();
1357 // Get associated node of the element (not side!)
1358 size_t co = m_innerSideCo [m_pFractGeo->scv(ip).node_id ()];
1359
1360 // line sources (that correspond to the point sources)
1361 for (t_lin_sss_iter line (m_sss_mngr.get (), m_innerFractSide, grid, aaPos, *m_pFractGeo, side_co);
1362 ! line.is_over (); ++line)
1363 {
1364 FVLineSourceOrSink<dim, dd_line_sss_data<dim> > * line_sss = *line;
1365 if (! line_sss->marked_for (m_innerFractSide, side_co))
1366 continue;
1367 line_sss->compute (line.seg_start (), this->time (), -1); //TODO: set the subset id instead of -1
1368 ass_sss_JA_elem (J, u, pElem, co, line_sss->intensity () / 2, line_sss->concentration ());
1369 /* Remark: "/ 2" because the source is taken into account twise. */
1370 }
1371 }
1372 }
1373}
1374
1376template<typename TDomain>
1377template<typename TElem>
1378void FractDDF_FV1<TDomain>::fract_bulk_ass_JA_elem
1379(
1380 LocalMatrix & J,
1381 const LocalVector & u,
1382 TElem * pElem,
1383 const position_type vCornerCoords []
1384)
1385{
1386// the permeability of the fracture-bulk interface
1387 number orthPermeability;
1388 if (m_imOrthoPermeability.data_given ())
1389 orthPermeability = m_imOrthoPermeability[0];
1390 else
1391 {
1392 MathVector<dim> orthPermeabilityVec;
1393 MatVecMult (orthPermeabilityVec, m_imPermeability[0], m_unitOuterNormal);
1394 orthPermeability = VecDot (orthPermeabilityVec, m_unitOuterNormal);
1395 }
1396
1397// the diffusion of the fracture-bulk interface
1398 number orthDiffusion;
1399 if (m_imOrthoDiffusion.data_given())
1400 orthDiffusion = m_imPorosity[0] * m_imOrthoDiffusion[0];
1401 else
1402 {
1403 MathVector<dim> orthDiffusionVec;
1404 MatVecMult (orthDiffusionVec, m_imDiffusion[0], m_unitOuterNormal);
1405 orthDiffusion = m_imPorosity[0] * VecDot (orthDiffusionVec, m_unitOuterNormal);
1406 }
1407
1408// the "height" of the degenerated element
1409 const number half_fr_width = m_imAperture[0] / 2;
1410
1411// loop over the corners of the inner side
1412 for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1413 {
1414 // Get current SCV
1415 const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1416 number s = scv.volume ();
1417
1418 // Get associated node of the element (not side!)
1419 const int co = m_innerSideCo [scv.node_id()];
1420
1421 // Get the corner values
1422 const number orthC_f = u(_C_, co);
1423 const number orthC_m = u(_C_, m_assCo[co]);
1424 const number orthP_f = u(_P_, co);
1425 const number orthP_m = u(_P_, m_assCo[co]);
1426
1427 const number fractDensity = m_imDensityCo[co];
1428 const number orthDensity = m_imDensityCo [m_assCo[co]];
1429 number D_fractDensity, D_orthDensity;
1430 if (! m_imDensityCo.constant())
1431 {
1432 D_fractDensity = m_imDensityCo.deriv (co, _C_, co);
1433 D_orthDensity = m_imDensityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
1434 }
1435 else D_fractDensity = D_orthDensity = 0;
1436
1437 const number orthViscosity = m_imViscosityCo [m_assCo[co]];
1438 const number D_orthViscosity = (m_imViscosityCo.constant())? 0.0
1439 : m_imViscosityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
1440
1441 // Flux velocity between the fracture and the bulk medium
1442
1443 number orthVelocity = (m_orthGravity * (orthDensity - fractDensity)
1444 - (orthP_m - orthP_f) / half_fr_width)
1445 * orthPermeability;
1446
1447 number D_orthVelocity [2], D_orthVelocity_fr [2]; // two components: w.r.t. _C_ and _P_ resp.
1448
1449 D_orthVelocity [_C_] = (m_orthGravity * D_orthDensity
1450 * orthPermeability * orthViscosity
1451 - orthVelocity * D_orthViscosity) / orthViscosity / orthViscosity;
1452
1453 D_orthVelocity_fr [_C_] = - m_orthGravity * D_fractDensity
1454 * orthPermeability / orthViscosity;
1455
1456 D_orthVelocity [_P_] = - orthPermeability / half_fr_width / orthViscosity;
1457
1458 D_orthVelocity_fr [_P_] = - D_orthVelocity [_P_];
1459
1460 orthVelocity /= orthViscosity;
1461
1463 number flux, D_flux [2], D_flux_fr [2]; // two components of the derivatives: w.r.t. _C_ and _P_ resp.
1464
1466 // Transport Equation
1468
1469 D_flux [_C_] = orthDiffusion / half_fr_width;
1470 D_flux_fr [_C_] = - D_flux [_C_];
1471 /* We use the full upwind here: */
1472 if (orthVelocity >= 0)
1473 {
1474 D_flux [_C_] -= D_orthVelocity [_C_] * orthC_f;
1475 D_flux_fr [_C_] -= orthVelocity + D_orthVelocity_fr [_C_] * orthC_f;
1476 D_flux [_P_] = - D_orthVelocity [_P_] * orthC_f;
1477 D_flux_fr [_P_] = - D_orthVelocity_fr [_P_] * orthC_f;
1478 }
1479 else
1480 {
1481 D_flux [_C_] -= orthVelocity + D_orthVelocity [_C_] * orthC_m;
1482 D_flux_fr [_C_] -= D_orthVelocity_fr [_C_] * orthC_m;
1483 D_flux [_P_] = - D_orthVelocity [_P_] * orthC_m;
1484 D_flux_fr [_P_] = - D_orthVelocity_fr [_P_] * orthC_m;
1485 }
1486
1487 if (! m_BoussinesqTransport)
1488 {
1489 flux = orthDiffusion * (orthC_m - orthC_f) / half_fr_width;
1490 /* We use the full upwind here: */
1491 flux -= orthVelocity * ((orthVelocity >= 0)? orthC_f : orthC_m);
1492
1493 D_flux [_C_] = D_flux [_C_] * orthDensity + flux * D_orthDensity;
1494 D_flux_fr [_C_] *= orthDensity;
1495 D_flux [_P_] *= orthDensity;
1496 D_flux_fr [_P_] *= orthDensity;
1497 }
1498
1499 J(_C_, m_assCo [co], _C_, m_assCo [co]) += D_flux [_C_] * s;
1500 J(_C_, m_assCo [co], _P_, m_assCo [co]) += D_flux [_P_] * s;
1501 J(_C_, m_assCo [co], _C_, co) += D_flux_fr [_C_] * s;
1502 J(_C_, m_assCo [co], _P_, co) += D_flux_fr [_P_] * s;
1503
1504 J(_C_, co, _C_, m_assCo [co]) -= D_flux [_C_] * s;
1505 J(_C_, co, _P_, m_assCo [co]) -= D_flux [_P_] * s;
1506 J(_C_, co, _C_, co) -= D_flux_fr [_C_] * s;
1507 J(_C_, co, _P_, co) -= D_flux_fr [_P_] * s;
1508
1510 // Flow Equation
1512
1513 if (! m_BoussinesqFlow)
1514 {
1515 D_flux [_C_] = (D_orthVelocity [_C_] * orthDensity
1516 + orthVelocity * D_orthDensity) * s;
1517 D_flux_fr [_C_] = D_orthVelocity_fr [_C_] * orthDensity * s;
1518 }
1519 else
1520 {
1521 D_flux [_C_] = D_orthVelocity [_C_] * s;
1522 D_flux_fr [_C_] = D_orthVelocity_fr [_C_] * s;
1523 }
1524 J(_P_, m_assCo [co], _C_, m_assCo [co]) -= D_flux [_C_];
1525 J(_P_, co, _C_, m_assCo [co]) += D_flux [_C_];
1526 J(_P_, m_assCo [co], _C_, co) -= D_flux_fr [_C_];
1527 J(_P_, co, _C_, co) += D_flux_fr [_C_];
1528
1529 if (! m_BoussinesqFlow) s *= orthDensity;
1530
1531 flux = D_orthVelocity [_P_] * s;
1532 J(_P_, m_assCo [co], _P_, m_assCo [co]) -= flux;
1533 J(_P_, co, _P_, m_assCo [co]) += flux;
1534
1535 flux = D_orthVelocity_fr [_P_] * s;
1536 J(_P_, m_assCo [co], _P_, co) -= flux;
1537 J(_P_, co, _P_, co) += flux;
1538
1539 // Volume stabilization: Add the continuity equation to the transport equation
1540 // REMARK: We do not compute derivatives of the defect of the conf. eq. here because this defect tends to 0.
1541 if(m_bVolStabDataActive)
1542 {
1543 flux = orthVelocity * s;
1544 m_spVolStabData->stiff(pElem->vertex(m_assCo[co])) += flux;
1545 m_spVolStabData->stiff(pElem->vertex(co)) -= flux;
1546 //J(_C_, m_assCo[co], _C_, m_assCo[co]) += flux;
1547 //J(_C_, co, _C_, co ) -= flux;
1548 }
1549 }
1550}
1551
1553template<typename TDomain>
1554template<typename TElem>
1555void FractDDF_FV1<TDomain>::ass_sss_dA_elem
1556(
1557 LocalVector& d,
1558 const LocalVector& u,
1559 TElem * pElem,
1560 size_t co,
1561 number intensity,
1562 number concentration
1563)
1564{
1565 if (intensity > 0)
1566 {
1567 // source
1568 number density = 1000 + 000 * concentration; //FIXME!!
1569 d(_C_, co) -= m_BoussinesqTransport? intensity * concentration : intensity * density * concentration;
1570 d(_P_, co) -= m_BoussinesqFlow? intensity : intensity * density;
1571 if(m_bVolStabDataActive)
1572 m_spVolStabData->stiff(pElem->vertex(co)) -= m_BoussinesqFlow? intensity : intensity * density;
1573 }
1574 else
1575 {
1576 // sink
1577 d(_C_, co) -= m_BoussinesqTransport? intensity * u(_C_, co) : intensity * m_imDensityCo[co] * u(_C_, co);
1578 d(_P_, co) -= m_BoussinesqFlow? intensity : intensity * m_imDensityCo[co];
1579 if(m_bVolStabDataActive)
1580 m_spVolStabData->stiff(pElem->vertex(co))
1581 -= m_BoussinesqFlow? intensity : intensity * m_imDensityCo[co];
1582 }
1583}
1584
1586template<typename TDomain>
1587template<typename TElem>
1588void FractDDF_FV1<TDomain>::ass_sss_JA_elem
1589(
1590 LocalMatrix& J,
1591 const LocalVector& u,
1592 TElem * pElem,
1593 size_t co,
1594 number intensity,
1595 number concentration
1596)
1597{
1598 if (intensity > 0)
1599 {
1600 // source: no contribution to the maxrix
1601 number density = 1000 + 000 * concentration; //FIXME!!
1602 if(m_bVolStabDataActive)
1603 m_spVolStabData->stiff(pElem->vertex(co))
1604 -= m_BoussinesqFlow? intensity : intensity * density;
1605 }
1606 else
1607 {
1608 // sink: add the intensity to the matrix
1609 if (! m_BoussinesqTransport)
1610 {
1611 if (!m_imDensityCo.constant ())
1612 J(_C_, co, _C_, co)
1613 -= intensity * (m_imDensityCo[co] + m_imDensityCo.deriv (co, _C_, co) * u (_C_, co));
1614 else
1615 J(_C_, co, _C_, co) -= intensity * m_imDensityCo[co];
1616 }
1617 else
1618 J(_C_, co, _C_, co) -= intensity;
1619
1620 if (! m_BoussinesqFlow && ! m_imDensityCo.constant ())
1621 J(_P_, co, _C_, co) -= intensity * m_imDensityCo.deriv (co, _C_, co);
1622
1623 if(m_bVolStabDataActive)
1624 m_spVolStabData->stiff(pElem->vertex(co))
1625 -= m_BoussinesqFlow? intensity : intensity * m_imDensityCo[co];
1626 }
1627}
1628
1630template<typename TDomain>
1631template<typename TElem>
1632void FractDDF_FV1<TDomain>::ass_dM_elem
1633(
1634 LocalVector & d,
1635 const LocalVector & u,
1636 GridObject * elem,
1637 const position_type vCornerCoords []
1638)
1639{
1640 TElem * pElem = static_cast<TElem*> (elem);
1641
1642// call the specific routine:
1643 if (! m_isFracture)
1644 this->template bulk_ass_dM_elem<TElem> (d, u, pElem, vCornerCoords);
1645 else
1646 this->template fract_ass_dM_elem<TElem> (d, u, pElem, vCornerCoords);
1647}
1648
1650template<typename TDomain>
1651template<typename TElem>
1652void FractDDF_FV1<TDomain>::bulk_ass_dM_elem
1653(
1654 LocalVector & d,
1655 const LocalVector & u,
1656 TElem * pElem,
1657 const position_type vCornerCoords []
1658)
1659{
1660 typedef FV1Geometry<TElem, dim> TBulkFVGeom;
1661
1662// Get finite volume geometry
1663 const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
1664
1665// The porosity is constant per element
1666 number porosity = m_imPorosity [0];
1667
1668// Loop Sub Control Volumes (SCV)
1669 for (size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
1670 {
1671 // Get current SCV
1672 const typename TBulkFVGeom::SCV& scv = bulk_geo.scv(ip);
1673
1674 // Get associated node
1675 const int co = scv.node_id();
1676
1677 // Add to local matrix
1678 if(m_BoussinesqTransport)
1679 d(_C_,co) += porosity * u(_C_,co) * scv.volume();
1680 else
1681 d(_C_,co) += porosity * m_imDensityCo[ip] * u(_C_,co) * scv.volume();
1682
1683 if(m_BoussinesqFlow)
1684 d(_P_,co) += porosity * scv.volume();
1685 else
1686 {
1687 d(_P_,co) += porosity * m_imDensityCo[ip] * scv.volume();
1688
1689 if (m_bVolStabDataActive)
1690 {
1691 // Volume stabilization: Add the continuity equation to the transport equation
1692 // REMARK: We assume the implicit Euler scheme
1693 m_spVolStabData->mass(pElem->vertex(co))
1694 -= porosity * (m_imDensityCo[ip] - m_imOldDensityCo[ip])
1695 * scv.volume();
1696 }
1697 }
1698 }
1699}
1700
1702template<typename TDomain>
1703template<typename TElem>
1704void FractDDF_FV1<TDomain>::fract_ass_dM_elem
1705(
1706 LocalVector & d,
1707 const LocalVector & u,
1708 TElem * pElem,
1709 const position_type vCornerCoords []
1710)
1711{
1712// The porosity (constant per element) and fracture width ('/ 2' because every fracture is a boundary of 2 elements)
1713 number porosity = m_imPorosity[0] * m_imAperture[0] / 2;
1714
1715// Loop Sub Control Volumes (SCV)
1716 for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1717 {
1718 // Get current SCV
1719 const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1720
1721 // Get associated node of the element (not side!)
1722 const int co = m_innerSideCo [scv.node_id()];
1723
1724 // Add to local matrix
1725 if(m_BoussinesqTransport)
1726 d(_C_,co) += porosity * u(_C_,co) * scv.volume();
1727 else
1728 d(_C_,co) += porosity * m_imDensityCo[co] * u(_C_,co) * scv.volume();
1729
1730 if(m_BoussinesqFlow)
1731 d(_P_,co) += porosity * scv.volume();
1732 else
1733 {
1734 d(_P_,co) += porosity * m_imDensityCo[co] * scv.volume();
1735
1736 if (m_bVolStabDataActive)
1737 {
1738 // Volume stabilization: Add the continuity equation to the transport equation
1739 // REMARK: We assume the implicit Euler scheme
1740 m_spVolStabData->mass(pElem->vertex(co))
1741 -= porosity * (m_imDensityCo[co] - m_imOldDensityCo[co])
1742 * scv.volume();
1743 }
1744 }
1745 }
1746}
1747
1749template<typename TDomain>
1750template<typename TElem>
1751void FractDDF_FV1<TDomain>::ass_JM_elem
1752(
1753 LocalMatrix & J,
1754 const LocalVector & u,
1755 GridObject * elem,
1756 const position_type vCornerCoords []
1757)
1758{
1759 TElem * pElem = static_cast<TElem*> (elem);
1760
1761// call the specific routine:
1762 if (! m_isFracture)
1763 this->template bulk_ass_JM_elem<TElem> (J, u, pElem, vCornerCoords);
1764 else
1765 this->template fract_ass_JM_elem<TElem> (J, u, pElem, vCornerCoords);
1766}
1767
1769template<typename TDomain>
1770template<typename TElem>
1771void FractDDF_FV1<TDomain>::bulk_ass_JM_elem
1772(
1773 LocalMatrix & J,
1774 const LocalVector & u,
1775 TElem * pElem,
1776 const position_type vCornerCoords []
1777)
1778{
1779 typedef FV1Geometry<TElem, dim> TBulkFVGeom;
1780
1781// get finite volume geometry
1782 const TBulkFVGeom& bulk_geo = GeomProvider<TBulkFVGeom>::get();
1783
1784// The porosity is constant per element
1785 number porosity = m_imPorosity [0];
1786
1787// loop Sub Control Volumes (SCV)
1788 for (size_t ip = 0; ip < bulk_geo.num_scv(); ip++)
1789 {
1790 // get current SCV
1791 const typename TBulkFVGeom::SCV& scv = bulk_geo.scv(ip);
1792
1793 // get associated node
1794 const int co = scv.node_id();
1795
1796 // Add to local matrix
1797 if(m_BoussinesqTransport)
1798 J(_C_, co, _C_, co) += porosity * scv.volume();
1799 else
1800 J(_C_, co, _C_, co) +=
1801 porosity * scv.volume() *
1802 (m_imDensityCo[ip] + (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(ip, _C_, co)) * u(_C_, co));
1803
1804 if(!m_BoussinesqFlow)
1805 {
1806 J(_P_, co, _C_, co) += porosity * (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(ip, _C_, co)) * scv.volume();
1807
1808 if (m_bVolStabDataActive)
1809 {
1810 // Volume stabilization: Add the continuity equation to the transport equation
1811 // REMARK: We assume the implicit Euler scheme
1812 m_spVolStabData->mass(pElem->vertex(co))
1813 -= porosity * scv.volume() * (m_imDensityCo[ip] - m_imOldDensityCo[ip]);
1814
1815 /*TODO: Do we need this derivative? What we do is subtracting the defect of the flow
1816 * equation multiplied by _C_ from the transport equation. We differentiate the transport
1817 * equation, but not this artificial sink. Thus, we should skip this derivative. Otherwise
1818 * the contribution of the sink will not converge to 0 as soon as we approach the solution
1819 * of the system. However, in ug3, this derivative is present.
1820 *
1821 if(!m_imDensityScv.constant())
1822 m_spVolStabData->mass(pElem->vertex(co)) -=
1823 porosity * scv.volume() * m_imDensityCo.deriv(ip, _C_, co)*u(_C_, co);
1824 */
1825 }
1826 }
1827 //else
1828 //J(_P_, co, _C_, co) += 0;
1829
1830 // Remark: Other summands are zero
1831 //J(_C_, co, _P_, co) += 0;
1832 //J(_P_, co, _P_, co) += 0;
1833 }
1834}
1835
1837template<typename TDomain>
1838template<typename TElem>
1839void FractDDF_FV1<TDomain>::fract_ass_JM_elem
1840(
1841 LocalMatrix & J,
1842 const LocalVector & u,
1843 TElem * pElem,
1844 const position_type vCornerCoords []
1845)
1846{
1847// The porosity (constant per element) and fracture width ('/ 2' because every fracture is a boundary of 2 elements)
1848 number porosity = m_imPorosity[0] * m_imAperture[0] / 2;
1849
1850// loop Sub Control Volumes (SCV)
1851 for (size_t ip = 0; ip < m_pFractGeo->num_scv(); ip++)
1852 {
1853 // get current SCV
1854 const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
1855
1856 // Get associated node of the element (not side!)
1857 const int co = m_innerSideCo [scv.node_id()];
1858
1859 // Add to local matrix
1860 if(m_BoussinesqTransport)
1861 J(_C_, co, _C_, co) += porosity * scv.volume();
1862 else
1863 J(_C_, co, _C_, co) +=
1864 porosity * scv.volume() *
1865 (m_imDensityCo[ip] + (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _C_, co)) * u(_C_, co));
1866
1867 if(!m_BoussinesqFlow)
1868 {
1869 J(_P_, co, _C_, co) += porosity * (m_imDensityCo.constant() ? 0.0 : m_imDensityCo.deriv(co, _C_, co)) * scv.volume();
1870
1871 if (m_bVolStabDataActive)
1872 {
1873 // Volume stabilization: Add the continuity equation to the transport equation
1874 // REMARK: We assume the implicit Euler scheme
1875 m_spVolStabData->mass(pElem->vertex(co))
1876 -= porosity * scv.volume() * (m_imDensityCo[co] - m_imOldDensityCo[co]);
1877
1878 /*TODO: Do we need this derivative? What we do is subtracting the defect of the flow
1879 * equation multiplied by _C_ from the transport equation. We differentiate the transport
1880 * equation, but not this artificial sink. Thus, we should skip this derivative. Otherwise
1881 * the contribution of the sink will not converge to 0 as soon as we approach the solution
1882 * of the system. However, in ug3, this derivative is present.
1883 *
1884 if(!m_imDensityScv.constant())
1885 m_spVolStabData->mass(pElem->vertex(co)) -=
1886 porosity * scv.volume() * m_imDensityCo.deriv(co, _C_, co)*u(_C_, co);
1887 */
1888 }
1889 }
1890 //else
1891 //J(_P_, co, _C_, co) += 0;
1892
1893 // Remark: Other summands are zero
1894 //J(_C_, co, _P_, co) += 0;
1895 //J(_P_, co, _P_, co) += 0;
1896 }
1897}
1898
1900template<typename TDomain>
1901template<typename TElem>
1902void FractDDF_FV1<TDomain>::ass_rhs_elem
1903(
1904 LocalVector & b,
1905 GridObject * elem,
1906 const position_type vCornerCoords []
1907)
1908{
1909}
1910
1912// implementation of the export parameters
1914
1916template<typename TDomain>
1917template<typename TElem>
1918void FractDDF_FV1<TDomain>::
1919ex_brine
1920(
1921 number vValue[],
1922 const MathVector<dim> vGlobIP[],
1923 number time, int si,
1924 const LocalVector& u,
1925 GridObject* elem,
1926 const MathVector<dim> vCornerCoords[],
1927 const MathVector<dim> vLocIP[],
1928 const size_t nip,
1929 bool bDeriv,
1930 std::vector<std::vector<number> > vvvDeriv[]
1931)
1932{
1933// reference element
1934 typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
1935
1936// number of shape functions
1937 static const size_t numSH = ref_elem_type::numCorners;
1938
1939// get trial space
1940 LagrangeP1<ref_elem_type> rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
1941
1942// storage for shape function at ip
1943 number vShape[numSH];
1944
1945// loop ips
1946 for(size_t ip = 0; ip < nip; ++ip)
1947 {
1948 // evaluate at shapes at ip
1949 rTrialSpace.shapes(vShape, vLocIP[ip]);
1950
1951 // compute concentration at ip
1952 vValue[ip] = 0.0;
1953 for(size_t sh = 0; sh < numSH; ++sh)
1954 vValue[ip] += u(_C_, sh) * vShape[sh];
1955
1956 // compute derivative w.r.t. to unknowns iff needed
1957 if(bDeriv)
1958 for(size_t sh = 0; sh < numSH; ++sh)
1959 {
1960 vvvDeriv[ip][_C_][sh] = vShape[sh];
1961 vvvDeriv[ip][_P_][sh] = 0.0;
1962 }
1963 }
1964}
1965
1967template<typename TDomain>
1968template<typename TElem>
1969void FractDDF_FV1<TDomain>::
1970ex_pressure
1971(
1972 number vValue[],
1973 const MathVector<dim> vGlobIP[],
1974 number time, int si,
1975 const LocalVector& u,
1976 GridObject* elem,
1977 const MathVector<dim> vCornerCoords[],
1978 const MathVector<dim> vLocIP[],
1979 const size_t nip,
1980 bool bDeriv,
1981 std::vector<std::vector<number> > vvvDeriv[]
1982)
1983{
1984// reference element
1985 typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
1986
1987// number of shape functions
1988 static const size_t numSH = ref_elem_type::numCorners;
1989
1990// get trial space
1991 LagrangeP1<ref_elem_type> rTrialSpace = Provider<LagrangeP1<ref_elem_type> >::get();
1992
1993// storage for shape function at ip
1994 number vShape[numSH];
1995
1996// loop ips
1997 for(size_t ip = 0; ip < nip; ++ip)
1998 {
1999 // evaluate at shapes at ip
2000 rTrialSpace.shapes(vShape, vLocIP[ip]);
2001
2002 // compute concentration at ip
2003 vValue[ip] = 0.0;
2004 for(size_t sh = 0; sh < numSH; ++sh)
2005 vValue[ip] += u(_P_, sh) * vShape[sh];
2006
2007 // compute derivative w.r.t. to unknowns iff needed
2008 if(bDeriv)
2009 for(size_t sh = 0; sh < numSH; ++sh)
2010 {
2011 vvvDeriv[ip][_C_][sh] = 0;
2012 vvvDeriv[ip][_P_][sh] = vShape[sh];
2013 }
2014 }
2015}
2016
2018template<typename TDomain>
2019template<typename TElem>
2020void FractDDF_FV1<TDomain>::
2021ex_darcy_fract
2022(
2023 MathVector<dim> vValue[],
2024 const MathVector<dim> vGlobIP[],
2025 number time, int si,
2026 const LocalVector& u,
2027 GridObject* elem,
2028 const MathVector<dim> vCornerCoords[],
2029 const MathVector<dim> vLocIP[],
2030 const size_t nip,
2031 bool bDeriv,
2032 std::vector<std::vector<MathVector<dim> > > vvvDeriv[]
2033)
2034{
2035 if (m_isFracture)
2036 {
2037 typedef StdLinConsistentGravity<low_dim> TConsGravity;
2038
2039 // this evaluation is implemented only for the standard set of the IPs
2040 if (nip != m_pFractGeo->num_scvf())
2041 UG_THROW ("FractDDF_FV1: The Darcy velocity export parameter is only implemented for the standard set of IPs.");
2042
2043 // number of shape functions
2044 const size_t numSh = m_pFractGeo->num_sh();
2045
2046 MathVector<dim> Vel;
2047 MathVector<dim> Vel_c[maxFractSideCorners], Vel_p[maxFractSideCorners];
2048 number vViscosity_c [maxFractSideCorners], * pViscosity_c;
2049
2050 // Get the corner values
2051 number vDensity [maxFractSideCorners];
2052 number vPressure [maxFractSideCorners];
2053 for (size_t sh = 0; sh < numSh; sh++)
2054 {
2055 size_t co = m_innerSideCo[sh];
2056 vDensity[sh] = m_imDensityCo[co];
2057 vPressure[sh] = u(_P_, co);
2058 }
2059
2060 // Consistent gravity and its derivative at corners
2061 TConsGravity ConsGravityMethod;
2062 MathVector<low_dim> vConsGravity [maxFractSideCorners];
2063 MathVector<low_dim> vConsGravity_c [maxFractSideCorners][maxFractSideCorners];
2064
2065 // Prepare the consistent gravity
2066 try
2067 {
2068 ConsGravityMethod.template prepare<dim>
2069 (vConsGravity, numSh, m_pFractGeo->corners(), vDensity, m_Gravity);
2070
2071 if(bDeriv)
2072 {
2073 number co_density [maxFractSideCorners];
2074 memset (co_density, 0, numSh * sizeof (number));
2075 for (size_t sh = 0; sh < numSh; sh++)
2076 { //TODO: use scv, not merely shape id here!
2077 size_t co = m_innerSideCo[sh];
2078 co_density[sh] = (m_imDensityCo.constant()) ? 0.0 : m_imDensityCo.deriv (co, _C_, co); // we assume the completely local dependence of \f$\rho\f$ of \f$c\f$
2079 ConsGravityMethod.template prepare<dim>
2080 (vConsGravity_c[sh], numSh, m_pFractGeo->corners(), co_density, m_Gravity);
2081 co_density[sh] = 0.0;
2082 }
2083 }
2084 }
2085 UG_CATCH_THROW ("FractDDF_FV1::fract_ass_dA_elem: Cannot prepare Consistent Gravity.");
2086
2087 // Loop the ips
2088 for (size_t ip = 0; ip < nip; ip++)
2089 {
2090 // Darcy velocity to be filled
2091 MathVector<dim>& DarcyVel = vValue[ip];
2092
2093 // Compute the Darcy velocity at the ips and the convection shapes
2094 this->template compute_ip_Darcy_velocity<TFractFVGeom, TConsGravity>
2095 (Vel, ip, *m_pFractGeo, ConsGravityMethod, vConsGravity, vPressure,
2096 m_imViscosityIP[ip]);
2097 VecScale (DarcyVel, Vel, m_imFractPermeability[0]);
2098
2099 // Compute the derivatives
2100 if(bDeriv)
2101 {
2102 // Derivatives to be filled
2103 MathVector<dim>* DarcyVel_c = &vvvDeriv[ip][_C_][0]; // array indexed by the shape
2104 MathVector<dim>* DarcyVel_p = &vvvDeriv[ip][_P_][0]; // array indexed by the shape
2105
2106 // Clear all the derivatives (as we fill up not all of them):
2107 for (size_t co = 0; co < vvvDeriv[ip][_C_].size(); co++)
2108 DarcyVel_c[co] = 0.0;
2109 for (size_t co = 0; co < vvvDeriv[ip][_P_].size(); co++)
2110 DarcyVel_p[co] = 0.0;
2111
2112 // Compute the derivative of the viscosity
2113 if (m_imViscosityIP.constant ())
2114 pViscosity_c = NULL;
2115 else
2116 {
2117 for (size_t sh = 0; sh < numSh; sh++)
2118 {
2119 size_t co = m_innerSideCo[sh];
2120 vViscosity_c[sh] = m_imViscosityIP.deriv(ip, _C_, co);
2121 }
2122 pViscosity_c = vViscosity_c;
2123 }
2124
2125 // Compute derivatives of the Darcy velocity
2126 this->template compute_J_ip_Darcy_velocity <TFractFVGeom, TConsGravity, maxFractSideCorners>
2127 (Vel, Vel_c, Vel_p, ip, *m_pFractGeo, ConsGravityMethod, vConsGravity_c, vPressure,
2128 m_imViscosityIP[ip], pViscosity_c);
2129 for (size_t sh = 0; sh < numSh; sh++)
2130 {
2131 size_t co = m_innerSideCo[sh];
2132 VecScale (DarcyVel_c[co], Vel_c[sh], m_imFractPermeability[0]);
2133 VecScale (DarcyVel_p[co], Vel_p[sh], m_imFractPermeability[0]);
2134 }
2135 }
2136 }
2137 }
2138 else
2139 UG_THROW ("FractDDF_FV1: The Darcy velocity export parameter is currently implemented only for fractures.");
2140}
2141
2143template<typename TDomain>
2144template<typename TElem>
2145void FractDDF_FV1<TDomain>::
2146ex_darcy_ortho_fract
2147(
2148 number vValue[],
2149 const MathVector<dim> vGlobIP[],
2150 number time, int si,
2151 const LocalVector& u,
2152 GridObject* elem,
2153 const MathVector<dim> vCornerCoords[],
2154 const MathVector<dim> vLocIP[],
2155 const size_t nip,
2156 bool bDeriv,
2157 std::vector<std::vector<number> > vvvDeriv[]
2158)
2159{
2160 if (m_isFracture)
2161 {
2162 // this evaluation is implemented only for the standard set of the IPs
2163 if (nip != m_pFractGeo->num_scv())
2164 UG_THROW ("FractDDF_FV1: The Darcy velocity export parameter is only implemented for the standard set of IPs.");
2165
2166 const number half_fr_width = m_imAperture[0] / 2;
2167
2168 // the permeability of the fracture-bulk interface
2169 number orthPermeability;
2170 if (m_imOrthoPermeability.data_given ())
2171 orthPermeability = m_imOrthoPermeability[0];
2172 else
2173 {
2174 MathVector<dim> orthPermeabilityVec;
2175 MatVecMult (orthPermeabilityVec, m_imPermeability[0], m_unitOuterNormal);
2176 orthPermeability = VecDot (orthPermeabilityVec, m_unitOuterNormal);
2177 }
2178
2179 // loop over the corners of the inner side
2180 for (size_t ip = 0; ip < nip; ip++)
2181 {
2182 // Get current SCV
2183 const typename TFractFVGeom::SCV& scv = m_pFractGeo->scv(ip);
2184 number& orthVelocity = vValue[ip];
2185
2186 // Get associated node of the element (not side!)
2187 const int co = m_innerSideCo [scv.node_id()];
2188
2189 // Get the corner values
2190 const number orthP_f = u(_P_, co);
2191 const number orthP_m = u(_P_, m_assCo[co]);
2192 const number fractDensity = m_imDensityCo[co];
2193 const number orthDensity = m_imDensityCo [m_assCo[co]];
2194 const number orthViscosity = m_imViscosityCo [m_assCo[co]];
2195
2196 // Flux velocity between the fracture and the bulk medium
2197 orthVelocity = (m_orthGravity * (orthDensity - fractDensity)
2198 - (orthP_m - orthP_f) / half_fr_width)
2199 * orthPermeability;
2200
2201 // Compute the derivatives
2202 if(bDeriv)
2203 {
2204 // Derivatives to be filled
2205 number* DarcyVel_c = &vvvDeriv[ip][_C_][0]; // array indexed by the shape
2206 number* DarcyVel_p = &vvvDeriv[ip][_P_][0]; // array indexed by the shape
2207
2208 // Clear all the derivatives (as we fill up not all of them):
2209 for (size_t i = 0; i < vvvDeriv[ip][_C_].size(); i++) DarcyVel_c[i] = 0;
2210 for (size_t i = 0; i < vvvDeriv[ip][_P_].size(); i++) DarcyVel_p[i] = 0;
2211
2212 // Compute the derivatives:
2213 number D_fractDensity, D_orthDensity;
2214 if (! m_imDensityCo.constant())
2215 {
2216 D_fractDensity = m_imDensityCo.deriv (co, _C_, co);
2217 D_orthDensity = m_imDensityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
2218 }
2219 else D_fractDensity = D_orthDensity = 0;
2220
2221 const number D_orthViscosity = (m_imViscosityCo.constant())? 0.0
2222 : m_imViscosityCo.deriv (m_assCo[co], _C_, m_assCo[co]);
2223
2224 number D_orthVelocity [2], D_orthVelocity_fr [2]; // two components: w.r.t. _C_ and _P_ resp.
2225
2226 D_orthVelocity [_C_] = (m_orthGravity * D_orthDensity
2227 * orthPermeability * orthViscosity
2228 - orthVelocity * D_orthViscosity) / orthViscosity / orthViscosity;
2229
2230 D_orthVelocity_fr [_C_] = - m_orthGravity * D_fractDensity
2231 * orthPermeability / orthViscosity;
2232
2233 D_orthVelocity [_P_] = - orthPermeability / half_fr_width / orthViscosity;
2234
2235 D_orthVelocity_fr [_P_] = - D_orthVelocity [_P_];
2236
2237 DarcyVel_c[co] = D_orthVelocity_fr[_C_];
2238 DarcyVel_p[co] = D_orthVelocity_fr[_P_];
2239
2240 DarcyVel_c[m_assCo[co]] = D_orthVelocity[_C_];
2241 DarcyVel_p[m_assCo[co]] = D_orthVelocity[_P_];
2242 }
2243
2244 // The viscosity:
2245 orthVelocity /= orthViscosity;
2246 }
2247 }
2248 else
2249 UG_THROW ("FractDDF_FV1: The orthogonal Darcy velocity export parameter is implemented only for fractures.");
2250}
2251
2253// register assembling functions
2255
2257template<typename TDomain>
2258template<typename TElem> // the element to register for
2259void FractDDF_FV1<TDomain>::register_loc_discr_func ()
2260{
2261 static const ReferenceObjectID id = geometry_traits<TElem>::REFERENCE_OBJECT_ID;
2262 static const int refDim = reference_element_traits<TElem>::dim;
2263
2264 this->clear_add_fct(id);
2265
2266 this->set_prep_elem_loop_fct(id, & this_type::template prepare_element_loop<TElem>);
2267 this->set_prep_elem_fct (id, & this_type::template prepare_element<TElem>);
2268 this->set_fsh_elem_loop_fct (id, & this_type::template finish_element_loop<TElem>);
2269 this->set_add_jac_A_elem_fct(id, & this_type::template ass_JA_elem<TElem>);
2270 this->set_add_jac_M_elem_fct(id, & this_type::template ass_JM_elem<TElem>);
2271 this->set_add_def_A_elem_fct(id, & this_type::template ass_dA_elem<TElem>);
2272 this->set_add_def_M_elem_fct(id, & this_type::template ass_dM_elem<TElem>);
2273 this->set_add_rhs_elem_fct (id, & this_type::template ass_rhs_elem<TElem>);
2274
2275 m_exBrine->template set_fct<this_type, refDim> (id, this, &this_type::template ex_brine<TElem>);
2276 m_exPressure->template set_fct<this_type, refDim> (id, this, &this_type::template ex_pressure<TElem>);
2277 m_exFractDarcyVel->template set_fct<this_type,refDim>(id, this, &this_type::template ex_darcy_fract<TElem>);
2278 m_exOrthoFractDarcyVel->template set_fct<this_type,refDim>(id, this, &this_type::template ex_darcy_ortho_fract<TElem>);
2279}
2280
2284template <typename TDomain>
2285FractDDF_FV1<TDomain>::FractDDF_FV1
2286(
2287 const char* functions,
2288 const char* subsets
2289)
2290: IElemDisc<TDomain> (functions, subsets)
2291{
2292 init();
2293};
2294
2295template <typename TDomain>
2296FractDDF_FV1<TDomain>::FractDDF_FV1
2297(
2298 const std::vector<std::string>& vFct,
2299 const std::vector<std::string>& vSubset
2300)
2301: IElemDisc<TDomain> (vFct, vSubset)
2302{
2303 init();
2304};
2305
2306template <typename TDomain>
2307void FractDDF_FV1<TDomain>::init()
2308{
2309 set_boussinesq (false);
2310
2311 m_imDensityCo.set_comp_lin_defect (false);
2312 m_imOldDensityCo.set_comp_lin_defect (false);
2313 m_imDensityIP.set_comp_lin_defect (false),
2314 m_imViscosityCo.set_comp_lin_defect (false);
2315 m_imViscosityIP.set_comp_lin_defect (false);
2316 m_imConstGravity.set_comp_lin_defect (false);
2317 m_imPorosity.set_comp_lin_defect (false);
2318 m_imPermeability .set_comp_lin_defect(false);
2319 m_imDiffusion.set_comp_lin_defect (false);
2320 m_imAperture.set_comp_lin_defect (false);
2321 m_imFractPermeability.set_comp_lin_defect (false);
2322 m_imOrthoPermeability.set_comp_lin_defect (false);
2323 m_imFractDiffusion.set_comp_lin_defect (false);
2324 m_imOrthoDiffusion.set_comp_lin_defect (false);
2325 m_spUpwind = make_sp(new ConvectionShapesFullUpwind<dim>);
2326 m_imFractDarcyVelIP.set_comp_lin_defect (false);
2327
2328 m_spVolStabData = SPNULL;
2329
2330 std::string functions;
2331 for(size_t i = 0; i < this->symb_fcts().size(); ++i)
2332 {
2333 if(i > 0) functions.append(",");
2334 functions.append(this->symb_fcts()[i]);
2335 }
2336
2337 m_exBrine = make_sp (new DataExport<number, dim> (functions.c_str()));
2338 m_exPressure = make_sp (new DataExport<number, dim> (functions.c_str()));
2339 m_exFractDarcyVel = make_sp(new DataExport<MathVector<dim>, dim>(functions.c_str()));
2340 m_exOrthoFractDarcyVel = make_sp(new DataExport<number, dim>(functions.c_str()));
2341
2342// check number of functions
2343 if (this->num_fct () != 2)
2344 UG_THROW ("Wrong number of functions: The ElemDisc 'FractDDF_FV1'"
2345 " needs exactly 2 symbolic function"
2346 " (one for the mass fraction and one for the pressure).");
2347
2348// register imports
2349
2350 this->register_import (m_imDensityCo);
2351 this->register_import (m_imOldDensityCo);
2352 this->register_import (m_imDensityIP);
2353 this->register_import (m_imViscosityCo);
2354 this->register_import (m_imViscosityIP);
2355 this->register_import (m_imConstGravity);
2356
2357 this->register_import (m_imPorosity);
2358
2359 this->register_import (m_imPermeability);
2360 this->register_import (m_imDiffusion);
2361
2362 this->register_import (m_imAperture);
2363 this->register_import (m_imFractPermeability);
2364 this->register_import (m_imOrthoPermeability);
2365 this->register_import (m_imFractDiffusion);
2366 this->register_import (m_imOrthoDiffusion);
2367
2368 this->register_import(m_imFractDarcyVelIP);
2369
2370// register the local assembling functions
2371 boost::mpl::for_each<AssembleElemList> (RegisterLocalDiscr (this));
2372
2373}
2374
2375} // namespace d3f
2376} // end namespace ug
2377
2378/* End of File */
s
parameterString s
Definition Biogas.lua:2
ug::BaseLSFS::shapes
void shapes(shape_type *vShape, const MathVector< dim > &x) const
ug::DimReferenceMapping::local_to_global
virtual void local_to_global(MathVector< worldDim > &globPos, const MathVector< dim > &locPos) const=0
ug::DimReferenceMapping::update
virtual void update(const MathVector< worldDim > *vCornerCoord)=0
ug::FVPointSourceOrSink::position
const MathVector< dim > & position()
ug::GeomProvider::get
static TGeom & get()
ug::IConvectionShapes::D_vel
const MathVector< dim > & D_vel(size_t scvf, size_t sh) const
ug::Provider::get
static TClass & get()
ug::d3f::DDSingularSourcesAndSinks
Definition d3f_sss.h:317
ug::d3f::FractDDF_FV1::prep_assemble_loop
virtual void prep_assemble_loop()
called once bevore assembling
Definition fract_ddf_fv1_impl.h:64
ug::d3f::FractDDF_FV1::ex_darcy_ortho_fract
void ex_darcy_ortho_fract(number vValue[], const MathVector< dim > vGlobIP[], number time, int si, const LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< dim > vLocIP[], const size_t nip, bool bDeriv, std::vector< std::vector< number > > vvvDeriv[])
export parameter for the orthogonal Darcy velocity
Definition fract_ddf_fv1_impl.h:2147
ug::d3f::FractDDF_FV1::fract_ass_JM_elem
void fract_ass_JM_elem(LocalMatrix &J, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the mass matrix of a time-dependent problem on a fracture element
Definition fract_ddf_fv1_impl.h:1840
ug::d3f::FractDDF_FV1::bulk_ass_dM_elem
void bulk_ass_dM_elem(LocalVector &d, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the mass part of the defect of a time-dependent problem on a bulk element
Definition fract_ddf_fv1_impl.h:1653
ug::d3f::FractDDF_FV1::ass_JA_elem
void ass_JA_elem(LocalMatrix &J, const LocalVector &u, GridObject *elem, const position_type vCornerCoords[])
computes the local stiffness matrix
Definition fract_ddf_fv1_impl.h:890
ug::d3f::FractDDF_FV1::prepare_setting
virtual void prepare_setting(const std::vector< LFEID > &vLfeID, bool bNonRegular)
check type of the grid and the trial space
Definition fract_ddf_fv1_impl.h:39
ug::d3f::FractDDF_FV1::bulk_ass_dA_elem
void bulk_ass_dA_elem(LocalVector &d, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the stiffness part of the local defect on a bulk element
Definition fract_ddf_fv1_impl.h:466
ug::d3f::FractDDF_FV1::FractDDF_FV1
FractDDF_FV1(const char *functions, const char *subsets)
class constructor
Definition fract_ddf_fv1_impl.h:2286
ug::d3f::FractDDF_FV1::register_loc_discr_func
void register_loc_discr_func()
registers the local assembler functions for a given element
Definition fract_ddf_fv1_impl.h:2259
ug::d3f::FractDDF_FV1::ass_dA_elem
void ass_dA_elem(LocalVector &d, const LocalVector &u, GridObject *elem, const position_type vCornerCoords[])
computes the stiffness part of the local defect
Definition fract_ddf_fv1_impl.h:441
ug::d3f::FractDDF_FV1::fract_bulk_ass_JA_elem
void fract_bulk_ass_JA_elem(LocalMatrix &J, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the local stiffness matrix of the fracture-bulk interaction terms on a fracture element
Definition fract_ddf_fv1_impl.h:1379
ug::d3f::FractDDF_FV1::ex_brine
void ex_brine(number vValue[], const MathVector< dim > vGlobIP[], number time, int si, const LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< dim > vLocIP[], const size_t nip, bool bDeriv, std::vector< std::vector< number > > vvvDeriv[])
export parameter for the concentration (to compute density and viscosity)
Definition fract_ddf_fv1_impl.h:1920
ug::d3f::FractDDF_FV1::bulk_prepare_element
void bulk_prepare_element(const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
prepares a given bulk element for assembling
Definition fract_ddf_fv1_impl.h:255
ug::d3f::FractDDF_FV1::bulk_ass_JA_elem
void bulk_ass_JA_elem(LocalMatrix &J, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the local stiffness matrix on a bulk element
Definition fract_ddf_fv1_impl.h:913
ug::d3f::FractDDF_FV1::fract_prepare_element_loop
void fract_prepare_element_loop(ReferenceObjectID roid, int si)
prepares the loop over the elements: the 'fracture' version
Definition fract_ddf_fv1_impl.h:183
ug::d3f::FractDDF_FV1::bulk_prepare_element_loop
void bulk_prepare_element_loop(ReferenceObjectID roid, int si)
prepares the loop over the elements: the 'bulk' version
Definition fract_ddf_fv1_impl.h:133
ug::d3f::FractDDF_FV1::bulk_ass_JM_elem
void bulk_ass_JM_elem(LocalMatrix &J, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the mass matrix of a time-dependent problem on a bulk element
Definition fract_ddf_fv1_impl.h:1772
ug::d3f::FractDDF_FV1::compute_ip_Darcy_velocity
void compute_ip_Darcy_velocity(MathVector< dim > &Vel, size_t ip, const TFVGeom &bulk_geo, TConsGravity &ConsGravityMethod, MathVector< TFVGeom::dim > vConsGravity[], number vPressure[], number Viscosity)
computes the Darcy velocity (not scaled with the permeability)
Definition fract_ddf_fv1_impl.h:413
ug::d3f::FractDDF_FV1::fract_bulk_ass_dA_elem
void fract_bulk_ass_dA_elem(LocalVector &d, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the stiffness fracture-bulk interaction terms of the local defect on a fracture element
Definition fract_ddf_fv1_impl.h:758
ug::d3f::FractDDF_FV1::finish_element_loop
void finish_element_loop()
finalizes the loop over the elements
Definition fract_ddf_fv1_impl.h:227
ug::d3f::FractDDF_FV1::fract_ass_JA_elem
void fract_ass_JA_elem(LocalMatrix &J, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the local stiffness matrix on a fracture element
Definition fract_ddf_fv1_impl.h:1166
ug::d3f::FractDDF_FV1::ass_sss_JA_elem
void ass_sss_JA_elem(LocalMatrix &J, const LocalVector &u, TElem *pElem, size_t co, number intensity, number concentration)
assembles a singular source or sink in the jacobian
Definition fract_ddf_fv1_impl.h:1589
ug::d3f::FractDDF_FV1::compute_J_ip_Darcy_velocity
void compute_J_ip_Darcy_velocity(MathVector< dim > &Vel, MathVector< dim > Vel_c[], MathVector< dim > Vel_p[], size_t ip, const TFVGeom &geo, TConsGravity &ConsGravityMethod, MathVector< TFVGeom::dim > vConsGravity_c[][maxCorners], number vPressure[], number Viscosity, number Viscosity_c[])
computes the derivatives of the Darcy velocity (not scaled with the permeability)
Definition fract_ddf_fv1_impl.h:852
ug::d3f::FractDDF_FV1::ass_rhs_elem
void ass_rhs_elem(LocalVector &d, GridObject *elem, const position_type vCornerCoords[])
computes the right-hand side due to the sources
Definition fract_ddf_fv1_impl.h:1903
ug::d3f::FractDDF_FV1::fract_ass_dA_elem
void fract_ass_dA_elem(LocalVector &d, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the stiffness part of the local defect on a fracture element
Definition fract_ddf_fv1_impl.h:627
ug::d3f::FractDDF_FV1::prepare_element_loop
void prepare_element_loop(ReferenceObjectID roid, int si)
prepares the loop over the elements: checks whether the parameters are set, ...
Definition fract_ddf_fv1_impl.h:82
ug::d3f::FractDDF_FV1::position_type
base_type::position_type position_type
position type
Definition fract_ddf_fv1.h:78
ug::d3f::FractDDF_FV1::ass_JM_elem
void ass_JM_elem(LocalMatrix &J, const LocalVector &u, GridObject *elem, const position_type vCornerCoords[])
computes the mass matrix of a time-dependent problem
Definition fract_ddf_fv1_impl.h:1752
ug::d3f::FractDDF_FV1::ex_darcy_fract
void ex_darcy_fract(MathVector< dim > vValue[], const MathVector< dim > vGlobIP[], number time, int si, const LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< dim > vLocIP[], const size_t nip, bool bDeriv, std::vector< std::vector< MathVector< dim > > > vvvDeriv[])
export parameter for the Darcy velocity in the fracture
Definition fract_ddf_fv1_impl.h:2022
ug::d3f::FractDDF_FV1::ass_dM_elem
void ass_dM_elem(LocalVector &d, const LocalVector &u, GridObject *elem, const position_type vCornerCoords[])
computes the mass part of the defect of a time-dependent problem
Definition fract_ddf_fv1_impl.h:1633
ug::d3f::FractDDF_FV1::ex_pressure
void ex_pressure(number vValue[], const MathVector< dim > vGlobIP[], number time, int si, const LocalVector &u, GridObject *elem, const MathVector< dim > vCornerCoords[], const MathVector< dim > vLocIP[], const size_t nip, bool bDeriv, std::vector< std::vector< number > > vvvDeriv[])
export parameter for the pressure
Definition fract_ddf_fv1_impl.h:1971
ug::d3f::FractDDF_FV1::init
void init()
Definition fract_ddf_fv1_impl.h:2307
ug::d3f::FractDDF_FV1::fract_prepare_element
void fract_prepare_element(const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
prepares a given fracture element for assembling
Definition fract_ddf_fv1_impl.h:294
ug::d3f::FractDDF_FV1::fract_ass_dM_elem
void fract_ass_dM_elem(LocalVector &d, const LocalVector &u, TElem *elem, const position_type vCornerCoords[])
computes the mass part of the defect of a time-dependent problem on a fracture element
Definition fract_ddf_fv1_impl.h:1705
ug::d3f::FractDDF_FV1::prepare_element
void prepare_element(const LocalVector &u, GridObject *elem, ReferenceObjectID roid, const position_type vCornerCoords[])
prepares a given element for assembling
Definition fract_ddf_fv1_impl.h:235
ug::d3f::FractDDF_FV1::ass_sss_dA_elem
void ass_sss_dA_elem(LocalVector &d, const LocalVector &u, TElem *pElem, size_t co, number intensity, number concentration)
assembles a singular source or sink in the defect
Definition fract_ddf_fv1_impl.h:1556
ug::geometry_traits
const_user_data.h
conv_shape.h
get
function get(x, y, z, t)
Viscosity
function util d3f parse Viscosity(ViscosityDesc, w)
fv1_geom.h
geom_provider.h
subset_handler
SmartPtr< TSubsetHandler > subset_handler()
grid
SmartPtr< TGrid > grid()
ug::MatSet
void MatSet(matrix_t &mInOut, typename matrix_t::value_type s)
ug::MatDiagSet
void MatDiagSet(matrix_t &mInOut, typename matrix_t::value_type s)
ug::MatScale
void MatScale(matrix_t &mOut, typename matrix_t::value_type s, const matrix_t &m)
SPNULL
const NullSmartPtr SPNULL
UG_CATCH_THROW
#define UG_CATCH_THROW(msg)
UG_THROW
#define UG_THROW(msg)
number
double number
ug::MatVecMult
void MatVecMult(vector_t_out &vOut, const matrix_t &m, const vector_t_in &v)
ug::VecLength
vector_t::value_type VecLength(const vector_t &v)
ug::VecScaleAppend
void VecScaleAppend(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1)
ug::VecDistance
vector_t::value_type VecDistance(const vector_t &v1, const vector_t &v2)
ug::VecScale
void VecScale(vector_t &vOut, const vector_t &v, typename vector_t::value_type s)
ug::VecDot
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
lua_user_data.h
ug::ReferenceObjectID
ReferenceObjectID
ug::VecSet
void VecSet(vector_t &dest, number alpha, const std::vector< size_t > vIndex)
make_sp
SmartPtr< T, FreePolicy > make_sp(T *inst)
ug::d3f::FractDDF_FV1::RegisterLocalDiscr
Definition fract_ddf_fv1.h:601
ug::reference_element_traits