Line data Source code
1 : /*
2 : * Copyright (c) 2026: CEMSE, KAUST
3 : * Author: Dmitry Logashenko
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__CONSISTENT_GRAVITY_LINKER__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__CONSISTENT_GRAVITY_LINKER__
35 :
36 : #include <vector>
37 :
38 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
39 : #include "lib_disc/spatial_disc/disc_util/consistent_gravity.h"
40 : #include "linker.h"
41 : #ifdef UG_FOR_LUA
42 : #include "bindings/lua/lua_user_data.h"
43 : #endif
44 :
45 : namespace ug{
46 :
47 :
48 : ////////////////////////////////////////////////////////////////////////////////
49 : // Consistent Gravity linker
50 : ////////////////////////////////////////////////////////////////////////////////
51 :
52 : /// Linker for the consistent gravity (according to P. Frolkovic)
53 : /**
54 : * This linker computes the s.c. consistent gravity \f$[ \rho \mathbf{g} ]_{consistent} \f$, according to Frolkovic,
55 : * where
56 : * <ul>
57 : * <li> \f$ \rho \f$ the density
58 : * <li> \f$ \mathbf{g} \f$ constant gravity
59 : * </ul>
60 : * is the input parameter. ( \f$ \mathbf{g} \f$ should be constant.)
61 : *
62 : * References:
63 : * <ul>
64 : * <li> P. Frolkovic, P. Knabner, Consistent Velocity Approximations in Finite
65 : * Element or Volume Discretizations of Density Driven Flow, In: Computational
66 : * Methods in WaterResources XI, Vol. 1 (A.A. Aldama et al., eds.),
67 : * Computational Mechanics Publication, Southhampten, 1996, p. 93-100
68 : * </li>
69 : * </ul>
70 : */
71 : template <int dim>
72 : class ConsistentGravityLinker
73 : : public StdDataLinker<ConsistentGravityLinker<dim>, MathVector<dim>, dim>
74 : {
75 : /// Base class type
76 : typedef StdDataLinker<ConsistentGravityLinker<dim>, MathVector<dim>, dim> base_type;
77 :
78 : /// The world dimension
79 : static const int world_dim = dim;
80 :
81 : public:
82 0 : ConsistentGravityLinker() :
83 : m_spDensity(NULL), m_spDDensity(NULL),
84 0 : m_gravity(0), m_noDerivatives(false)
85 : {
86 : this->set_num_input(1);
87 0 : }
88 :
89 :
90 0 : inline void evaluate (MathVector<dim>& value,
91 : const MathVector<dim>& globIP,
92 : number time, int si) const
93 : {
94 0 : UG_THROW("ConsistentGravityLinker: Element required for evaluation.")
95 : }
96 :
97 : template <int refDim>
98 0 : inline void evaluate(MathVector<dim> vValue[],
99 : const MathVector<dim> vGlobIP[],
100 : number time, int si,
101 : GridObject* elem,
102 : const MathVector<dim> vCornerCoords[],
103 : const MathVector<refDim> vLocalIP[],
104 : const size_t nip,
105 : LocalVector* u,
106 : const MathMatrix<refDim, dim>* vJT = NULL) const
107 : {
108 0 : const ReferenceObjectID roid = elem->reference_object_id();
109 : const DimReferenceElement<refDim>& refElem = ReferenceElementProvider::get<refDim>(roid);
110 : const size_t nco = refElem.num(ROID_VERTEX);
111 :
112 : const LocalShapeFunctionSet<refDim>& trialSpace =
113 0 : LocalFiniteElementProvider::get<refDim>(roid, LFEID(LFEID::LAGRANGE, refDim, 1));
114 :
115 : // consistent gravity
116 : StdLinConsistentGravity<refDim> ConsGravityMethod;
117 0 : std::vector<MathVector<refDim> > vConsGravity(nco);
118 0 : std::vector<MathVector<refDim> > vLocalGrad(nco);
119 :
120 : // coefficients at the integration points
121 : MathVector<dim> gravity; // we assume that the gravity is constant in the elem
122 0 : std::vector<number> vDensity(nco);
123 :
124 0 : (*m_spDensity)(&(vDensity[0]), vCornerCoords, time, si,
125 : elem, vCornerCoords, refElem.corners(), nco, u, NULL);
126 :
127 : try
128 : {
129 : ConsGravityMethod.template prepare<dim>
130 0 : (&(vConsGravity[0]), nco, vCornerCoords, &(vDensity[0]), m_gravity);
131 : }
132 0 : UG_CATCH_THROW ("ConsistentDarcyVelLinker: Cannot prepare the consistent gravity.");
133 :
134 : // get jacobians of the transformation mapping if not passed
135 : std::vector<MathMatrix<refDim, dim> > vJT_;
136 0 : if(vJT == NULL)
137 : {
138 : DimReferenceMapping<refDim, dim>& mapping
139 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
140 0 : vJT_.resize(nip);
141 0 : mapping.jacobian_transposed(&(vJT_[0]), vLocalIP, nip);
142 : vJT = &(vJT_[0]);
143 : }
144 :
145 0 : for(size_t ip = 0; ip < nip; ++ip)
146 : {
147 : // get the local gradient (assuming the Lagrange-1 basis functions)
148 0 : trialSpace.grads(&(vLocalGrad[0]), vLocalIP[ip]);
149 :
150 : // get the inverse Jacobian
151 : MathMatrix<dim,refDim> JTInv;
152 0 : RightInverse(JTInv, vJT[ip]);
153 :
154 : // compute [rho*g]_consistent
155 : ConsGravityMethod.template compute<dim>
156 0 : (vValue[ip], vLocalIP[ip], JTInv, &(vLocalGrad[0]), &(vConsGravity[0]));
157 : }
158 0 : }
159 :
160 : template <int refDim>
161 0 : void prepare_dim_elem(GridObject* elem,
162 : const ReferenceObjectID roid,
163 : const MathVector<dim> vCornerCoords[])
164 : {
165 : const DimReferenceElement<refDim>& refElem = ReferenceElementProvider::get<refDim>(roid);
166 :
167 0 : m_densityCornerS = m_spDensity->template register_local_ip_series<refDim>
168 0 : (refElem.corners(), refElem.num(0), this->time_point(), false);
169 0 : m_spDensity->set_global_ips(m_densityCornerS, vCornerCoords, refElem.num(0));
170 0 : }
171 :
172 0 : virtual void prepare_element(GridObject* elem, const MathVector<dim> vCornerCoords[])
173 : {
174 0 : const ReferenceObjectID roid = elem->reference_object_id();
175 : const ReferenceElement& refElem = ReferenceElementProvider::get(roid);
176 : const int ref_dim = refElem.dimension();
177 0 : switch(ref_dim)
178 : {
179 0 : case 1: this->template prepare_dim_elem<1> (elem, roid, vCornerCoords); break;
180 0 : case 2: this->template prepare_dim_elem<2> (elem, roid, vCornerCoords); break;
181 0 : case 3: this->template prepare_dim_elem<3> (elem, roid, vCornerCoords); break;
182 0 : default: UG_THROW("ConsistentDarcyVelLinker: Ref. dimension " << ref_dim << " not supported.");
183 : }
184 0 : }
185 :
186 : template <int refDim>
187 0 : void eval_and_deriv(MathVector<dim> vValue[],
188 : const MathVector<dim> vGlobIP[],
189 : number time, int si,
190 : GridObject* elem,
191 : const MathVector<dim> vCornerCoords[],
192 : const MathVector<refDim> vLocalIP[],
193 : const size_t nip,
194 : LocalVector* u,
195 : bool bDeriv,
196 : int s,
197 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
198 : const MathMatrix<refDim, dim>* vJT = NULL) const
199 : {
200 0 : const ReferenceObjectID roid = elem->reference_object_id();
201 : const DimReferenceElement<refDim>& refElem = ReferenceElementProvider::get<refDim>(roid);
202 : const size_t nco = refElem.num(ROID_VERTEX);
203 :
204 : const LocalShapeFunctionSet<refDim>& trialSpace =
205 0 : LocalFiniteElementProvider::get<refDim>(roid, LFEID(LFEID::LAGRANGE, refDim, 1));
206 :
207 : // consistent gravity
208 : StdLinConsistentGravity<refDim> ConsGravityMethod;
209 0 : std::vector<MathVector<refDim> > vConsGravity(nco);
210 :
211 : // inverse jacobians of the transformation at the integration points
212 0 : std::vector<MathMatrix<dim,refDim> > vJTInv(nip);
213 :
214 : // local gradients at the integration points
215 0 : std::vector<MathVector<refDim> > vLocalGrad(nco);
216 :
217 : // get the data of the ip series
218 0 : const number* vDensity = m_spDensity->values(m_densityCornerS);
219 :
220 : // prepare the consistent gravity term
221 : try
222 : {
223 : ConsGravityMethod.template prepare<dim>
224 0 : (&(vConsGravity[0]), nco, vCornerCoords, vDensity, m_gravity);
225 : }
226 0 : UG_CATCH_THROW ("ConsistentGravityLinker: Cannot prepare the consistent gravity.");
227 :
228 : // get jacobians of the transformation mapping if not passed
229 : std::vector<MathMatrix<refDim, dim> > vJT_;
230 0 : if(vJT == NULL)
231 : {
232 : DimReferenceMapping<refDim, dim>& mapping
233 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
234 0 : vJT_.resize(nip);
235 0 : mapping.jacobian_transposed(&(vJT_[0]), vLocalIP, nip);
236 : vJT = &(vJT_[0]);
237 : }
238 :
239 0 : for(size_t ip = 0; ip < nip; ++ip)
240 : {
241 : // get the local gradient (assuming the Lagrange-1 basis functions)
242 0 : trialSpace.grads(&(vLocalGrad[0]), vLocalIP[ip]);
243 :
244 : // get the inverse Jacobian
245 0 : RightInverse(vJTInv[ip], vJT[ip]);
246 :
247 : // compute [rho*g]_consistent
248 : ConsGravityMethod.template compute<dim>
249 0 : (vValue[ip], vLocalIP[ip], vJTInv[ip], &(vLocalGrad[0]), &(vConsGravity[0]));
250 : }
251 :
252 0 : if(!bDeriv)
253 : return;
254 :
255 : // Compute the derivatives at all ips
256 :
257 0 : this->set_zero(vvvDeriv, nip);
258 :
259 0 : if(m_noDerivatives || this->zero_derivative() || (!m_spDDensity.valid()) || m_spDDensity->zero_derivative())
260 0 : return;
261 :
262 : // prepare derivatives of the primary function at the corners
263 0 : std::vector<std::vector<MathVector<refDim> > > vvDConsGravity(nco);
264 : try
265 : {
266 0 : std::vector<number> DCoVal(nco);
267 0 : DCoVal.assign(nco, 0);
268 0 : for (size_t co = 0; co < nco; co++)
269 : {
270 0 : DCoVal[co] = 1;
271 0 : vvDConsGravity[co].resize(nco);
272 : ConsGravityMethod.template prepare<dim>
273 0 : (&(vvDConsGravity[co][0]), nco, vCornerCoords, &(DCoVal[0]), m_gravity);
274 0 : DCoVal[co] = 0;
275 : }
276 0 : }
277 0 : UG_CATCH_THROW ("ConsistentGravityLinker: Cannot prepare Consistent Gravity for its derivatives.");
278 :
279 : // compute the derivatives of the density
280 0 : if(m_noDerivatives || m_spDDensity->zero_derivative())
281 : return;
282 :
283 0 : for(size_t fct = 0; fct < m_spDDensity->num_fct(); ++fct) // fct = concentration, pressure, ...
284 : {
285 : // get common fct id for this function
286 : const size_t commonFct = this->input_common_fct(_RHO_, fct);
287 0 : if(this->num_sh(commonFct) != nco)
288 0 : UG_THROW ("ConsistentGravityLinker: Number of shapes mismatch.");
289 :
290 0 : for(size_t ip = 0; ip < nip; ++ip) // we derive the cons. grav. at this ip
291 : {
292 0 : trialSpace.grads(&(vLocalGrad[0]), vLocalIP[ip]);
293 :
294 0 : for(size_t co = 0; co < nco; ++co) // w.r.t to the DoF at this corner (shape idx.)
295 : {
296 0 : const number* vDDensity = m_spDDensity->deriv(m_densityCornerS, co, fct);
297 0 : MathVector<dim>& deriv = vvvDeriv[ip][commonFct][co];
298 :
299 0 : ConsGravityMethod.template compute<dim>(deriv,
300 : vLocalIP[ip], vJTInv[ip], &(vLocalGrad[0]), &(vvDConsGravity[co][0]));
301 0 : VecScale(deriv, deriv, vDDensity[co]);
302 :
303 : /* Note that here we assume that the density has the local dependence
304 : * on the arguments (concentration, pressure, ...): The density at
305 : * corner co depends only on the DoFs at that corner. However, this
306 : * excludes any differential operators, interpolations etc.
307 : */
308 : }
309 : }
310 : }
311 0 : }
312 :
313 : public:
314 :
315 : /// set density import
316 0 : void set_density(SmartPtr<CplUserData<number, dim> > data)
317 : {
318 0 : m_spDensity = data;
319 0 : m_spDDensity = data.template cast_dynamic<DependentUserData<number, dim> >();
320 0 : base_type::set_input(_RHO_, data, data);
321 0 : }
322 :
323 0 : void set_density(number val)
324 : {
325 0 : set_density(make_sp(new ConstUserNumber<dim>(val)));
326 0 : }
327 :
328 : /// set gravity vector
329 : void set_gravity(MathVector<dim> g)
330 : {
331 : m_gravity = g;
332 : }
333 :
334 : /// set gravity from an array
335 0 : void set_gravity(const std::vector<number>& vGravity)
336 : {
337 0 : if(vGravity.size() != dim)
338 0 : UG_THROW("ConsistentGravityLinker: Illegal dimension of the specified gravity vector.");
339 0 : for(size_t i = 0; i < dim; i++) m_gravity[i] = vGravity[i];
340 0 : }
341 :
342 : /// set gravity in the z direction
343 0 : void set_gravity(number g)
344 : {
345 : m_gravity = 0;
346 0 : m_gravity[dim-1] = g;
347 0 : }
348 :
349 : protected:
350 : /// import density
351 : static const size_t _RHO_ = 0;
352 : SmartPtr<CplUserData<number, dim> > m_spDensity;
353 : SmartPtr<DependentUserData<number, dim> > m_spDDensity;
354 : size_t m_densityCornerS;
355 :
356 : /// constant gravity
357 : MathVector<dim> m_gravity;
358 :
359 :
360 : public:
361 0 : void set_no_derivatives(bool v) { m_noDerivatives = v; }
362 :
363 : protected:
364 : // disable the derivatives
365 : bool m_noDerivatives;
366 : };
367 :
368 : } // end namespace ug
369 :
370 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__CONSISTENT_GRAVITY_LINKER__ */
|
