Line data Source code
1 : /*
2 : * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Arne Nägel
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__FUNCTION_SPACE__CONST_GRID_FUNCTION_USER_DATA__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACE__CONST_GRID_FUNCTION_USER_DATA__
35 :
36 :
37 : #include <map>
38 : #include <string>
39 :
40 : // ug4 headers
41 :
42 : #include "common/common.h"
43 :
44 : #include "lib_grid/tools/subset_group.h"
45 :
46 : #include "lib_disc/common/function_group.h"
47 : #include "lib_disc/common/groups_util.h"
48 : #include "lib_disc/quadrature/quadrature.h"
49 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
50 : #include "lib_disc/spatial_disc/user_data/std_user_data.h"
51 : #include "lib_disc/reference_element/reference_mapping_provider.h"
52 :
53 :
54 : namespace ug{
55 :
56 : /**
57 : * Base class for the CplUserData-interface for grid functions.
58 : *
59 : * This base class provides the CplUserData interface for the GridFunctions at integration
60 : * points. It can be used in particular to specify constant (or explicit) import parameters
61 : * for discretizations. Note that this interface does NOT compute derivatives (so, assuming
62 : * that the data are constant and do not depend on the current guess of the solution). For
63 : * the current solution in discretizations, so that the derivatives w.r.t. the DoFs are
64 : * taken into account, use the GridFunction...Data classes.
65 : */
66 : template <typename TImpl, typename TData, typename TGridFunction>
67 : class StdExplicitGridFunctionData
68 : : public StdUserData<StdExplicitGridFunctionData<TImpl,TData, TGridFunction>, TData, TGridFunction::dim>
69 : {
70 : public:
71 : /// world dimension of grid function
72 : static const int dim = TGridFunction::dim;
73 :
74 : protected:
75 : SmartPtr<TGridFunction> m_spGridFct; ///< grid function
76 : size_t m_fct; ///< component of function
77 : LFEID m_lfeID;
78 :
79 : protected:
80 : /// access to implementation
81 : TImpl& getImpl() {return static_cast<TImpl&>(*this);}
82 :
83 : /// const access to implementation
84 : const TImpl& getImpl() const {return static_cast<const TImpl&>(*this);}
85 :
86 : public:
87 : /// common constructor
88 0 : StdExplicitGridFunctionData
89 : (
90 : SmartPtr<TGridFunction> spGridFct ///< grid function
91 : )
92 0 : : m_spGridFct(spGridFct)
93 0 : {}
94 :
95 : /// returns if data is constant
96 0 : virtual bool constant() const {return false;}
97 :
98 : /// returns if grid function is needed for evaluation
99 : /** (true, since local coordinates may not be sufficient)*/
100 0 : virtual bool requires_grid_fct() const {return true;}
101 :
102 : /// returns if provided data is continuous over geometric object boundaries
103 0 : virtual bool continuous() const {return getImpl().continuous(); }
104 :
105 : template <int refDim>
106 0 : void eval(LocalVector* u, GridObject* elem,
107 : const MathVector<dim> vCornerCoords[], bool bDeriv = false)
108 : {
109 : const int si = this->subset();
110 :
111 0 : for(size_t s = 0; s < this->num_series(); ++s)
112 : {
113 0 : getImpl().template evaluate<refDim>(this->values(s), this->ips(s), this->time(s), si,
114 : elem, vCornerCoords,
115 : this->template local_ips<refDim>(s), this->num_ip(s),
116 : u);
117 : }
118 0 : }
119 :
120 : template <int refDim>
121 0 : void eval(LocalVectorTimeSeries* u, GridObject* elem,
122 : const MathVector<dim> vCornerCoords[], bool bDeriv = false)
123 : {
124 : const int si = this->subset();
125 :
126 0 : for(size_t s = 0; s < this->num_series(); ++s)
127 : {
128 0 : getImpl().template evaluate<refDim>(this->values(s), this->ips(s), this->time(s), si,
129 : elem, vCornerCoords,
130 : this->template local_ips<refDim>(s), this->num_ip(s),
131 : &(u->solution(this->time_point(s))));
132 : }
133 0 : }
134 :
135 0 : virtual void compute(LocalVector* u, GridObject* elem,
136 : const MathVector<dim> vCornerCoords[], bool bDeriv = false)
137 : {
138 : UG_ASSERT(elem->base_object_id() == this->dim_local_ips(),
139 : "local ip dimension and reference element dimension mismatch.");
140 :
141 0 : switch(this->dim_local_ips())
142 : {
143 0 : case 1: eval<1>(u,elem,vCornerCoords,bDeriv); break;
144 0 : case 2: eval<2>(u,elem,vCornerCoords,bDeriv); break;
145 0 : case 3: eval<3>(u,elem,vCornerCoords,bDeriv); break;
146 0 : default: UG_THROW("StdExplicitGridFunctionData: Dimension not supported.");
147 : }
148 0 : }
149 :
150 0 : virtual void compute(LocalVectorTimeSeries* u, GridObject* elem,
151 : const MathVector<dim> vCornerCoords[], bool bDeriv = false){
152 :
153 : UG_ASSERT(elem->base_object_id() == this->dim_local_ips(),
154 : "local ip dimension and reference element dimension mismatch.");
155 :
156 0 : switch(this->dim_local_ips())
157 : {
158 0 : case 1: eval<1>(u,elem,vCornerCoords,bDeriv); break;
159 0 : case 2: eval<2>(u,elem,vCornerCoords,bDeriv); break;
160 0 : case 3: eval<3>(u,elem,vCornerCoords,bDeriv); break;
161 0 : default: UG_THROW("StdExplicitGridFunctionData: Dimension not supported.");
162 : }
163 0 : }
164 :
165 0 : virtual void operator() (TData& value,
166 : const MathVector<dim>& globIP,
167 : number time, int si) const
168 : {
169 0 : UG_THROW("StdExplicitGridFunctionData: Solution, element and local ips required "
170 : "for evaluation, but not passed. Cannot evaluate.");
171 : }
172 :
173 0 : virtual void operator() (TData vValue[],
174 : const MathVector<dim> vGlobIP[],
175 : number time, int si, const size_t nip) const
176 : {
177 0 : UG_THROW("StdExplicitGridFunctionData: Solution, element and local ips required "
178 : "for evaluation, but not passed. Cannot evaluate.");
179 : }
180 :
181 : template <int refDim>
182 : inline void evaluate(TData vValue[],
183 : const MathVector<dim> vGlobIP[],
184 : number time, int si,
185 : GridObject* elem,
186 : const MathVector<dim> vCornerCoords[],
187 : const MathVector<refDim> vLocIP[],
188 : const size_t nip,
189 : LocalVector* u,
190 : const MathMatrix<refDim, dim>* vJT = NULL) const
191 : {
192 : // forward
193 0 : getImpl().template evaluate<refDim> (vValue, vGlobIP, time, si,
194 : elem, vCornerCoords, vLocIP, nip, u, vJT);
195 : }
196 : };
197 :
198 : /**
199 : * Class for the CplUserData-interface for the interpolation of a scalar component of a grid function.
200 : *
201 : * This class provides the CplUserData interface for the interpolation of a scalar component
202 : * of a GridFunction at integration points. It can be used in particular to specify constant
203 : * (or explicit) import parameters for discretizations. Note that this interface does NOT
204 : * compute derivatives (so, assuming that the data are constant and do not depend on the
205 : * current guess of the solution). For the interpolation of the current solution in
206 : * discretizations, so that the derivatives w.r.t. the DoFs are taken into account, use
207 : * the GridFunctionNumberData class.
208 : */
209 : template<typename TGridFunction>
210 : class ExplicitGridFunctionValue
211 : : public StdExplicitGridFunctionData<ExplicitGridFunctionValue<TGridFunction>, number, TGridFunction>
212 : {
213 : typedef StdExplicitGridFunctionData<ExplicitGridFunctionValue<TGridFunction>, number, TGridFunction> base_type;
214 :
215 : using base_type::m_spGridFct;
216 :
217 : size_t m_fct; ///< component of the function
218 : LFEID m_lfeID; ///< type of the shape functions
219 :
220 : public:
221 : // world dimension of grid function
222 : static const int dim = TGridFunction::dim;
223 :
224 : // constructor
225 0 : ExplicitGridFunctionValue
226 : (
227 : SmartPtr<TGridFunction> gf, ///< grid function
228 : const char* cmp ///< its component
229 : )
230 0 : : base_type(gf)
231 : {
232 : // get function id of name
233 0 : m_fct = m_spGridFct->fct_id_by_name(cmp);
234 :
235 : // check that function exists
236 0 : if(m_fct >= m_spGridFct->num_fct())
237 0 : UG_THROW("ExplicitGridFunctionValue: Function space does not contain"
238 : " a function with name " << cmp << ".");
239 :
240 : // local finite element id
241 0 : m_lfeID = m_spGridFct->local_finite_element_id(m_fct);
242 0 : }
243 :
244 : template <int refDim>
245 0 : inline void evaluate(number vValue[],
246 : const MathVector<dim> vGlobIP[],
247 : number time, int si,
248 : GridObject* elem,
249 : const MathVector<dim> vCornerCoords[],
250 : const MathVector<refDim> vLocIP[],
251 : const size_t nip,
252 : LocalVector* u,
253 : const MathMatrix<refDim, dim>* vJT = NULL) const
254 : {
255 : // reference object id
256 0 : const ReferenceObjectID roid = elem->reference_object_id();
257 :
258 : // get trial space
259 : try
260 : {
261 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
262 0 : LocalFiniteElementProvider::get<refDim>(roid, m_lfeID);
263 :
264 : // memory for shapes
265 : std::vector<number> vShape;
266 :
267 : // get multiindices of element
268 : std::vector<DoFIndex> ind;
269 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
270 :
271 : // loop ips
272 0 : for(size_t ip = 0; ip < nip; ++ip)
273 : {
274 : // evaluate at shapes at ip
275 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
276 :
277 : // compute solution at integration point
278 0 : vValue[ip] = 0.0;
279 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
280 : {
281 0 : const number valSH = DoFRef(*m_spGridFct, ind[sh]);
282 0 : vValue[ip] += valSH * vShape[sh];
283 : }
284 : }
285 :
286 0 : }
287 0 : UG_CATCH_THROW("ExplicitGridFunctionValue: Shape Function Set missing for"
288 : " Reference Object: "<<roid<<", Trial Space: "
289 : <<m_lfeID<<", refDim="<<refDim);
290 :
291 0 : }
292 :
293 0 : bool continuous() const {return LocalFiniteElementProvider::continuous(m_lfeID);}
294 :
295 : };
296 :
297 : /**
298 : * Class for the CplUserData-interface for the interpolation of a vector part of a grid function.
299 : *
300 : * This class provides the CplUserData interface for the interpolation of dim components
301 : * of a GridFunction as one vector at integration points. It can be used in particular to
302 : * specify constant (or explicit) import parameters for discretizations. Note that this
303 : * interface does NOT compute derivatives (so, assuming that the data are constant and do
304 : * not depend on the current guess of the solution). For the interpolation of the current
305 : * solution in discretizations, so that the derivatives w.r.t. the DoFs are taken into
306 : * account, use the GridFunctionVectorData class.
307 : */
308 : template<typename TGridFunction>
309 : class ExplicitGridFunctionVector
310 : : public StdExplicitGridFunctionData<ExplicitGridFunctionVector<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction>
311 : {
312 : public:
313 : // world dimension of grid function
314 : static const int dim = TGridFunction::dim;
315 :
316 : private:
317 : typedef StdExplicitGridFunctionData<ExplicitGridFunctionVector<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction> base_type;
318 :
319 : using base_type::m_spGridFct;
320 :
321 : size_t m_vfct[dim]; ///< components of the function
322 : LFEID m_vlfeID[dim]; ///< types of the shape functions
323 :
324 : public:
325 :
326 : // constructor
327 0 : ExplicitGridFunctionVector
328 : (
329 : SmartPtr<TGridFunction> gf, ///< grid function
330 : const char* cmps ///< its component
331 : )
332 0 : : base_type(gf)
333 : {
334 : try
335 : {
336 : // get strings
337 0 : std::string fctString (cmps);
338 :
339 : // tokenize strings and select functions
340 : std::vector<std::string> tokens;
341 0 : TokenizeString(fctString, tokens, ',');
342 :
343 0 : for(size_t i = 0; i < tokens.size(); ++i)
344 0 : RemoveWhitespaceFromString(tokens[i]);
345 :
346 0 : if((int)tokens.size() != dim)
347 0 : UG_THROW("ExplicitGridFunctionVector: Needed " << dim << " components "
348 : "in symbolic function names, but given: " << cmps << '.');
349 :
350 : // get function id of name
351 0 : for(int i = 0; i < dim; ++i)
352 : {
353 0 : m_vfct[i] = m_spGridFct->fct_id_by_name(tokens[i].c_str());
354 0 : m_vlfeID[i] = m_spGridFct->local_finite_element_id(m_vfct[i]);
355 : }
356 :
357 0 : }
358 0 : UG_CATCH_THROW("ExplicitGridFunctionVector: Cannot find some symbolic function name in '" << cmps << "'.");
359 0 : }
360 :
361 : template <int refDim>
362 0 : inline void evaluate(MathVector<dim> vValue[],
363 : const MathVector<dim> vGlobIP[],
364 : number time, int si,
365 : GridObject* elem,
366 : const MathVector<dim> vCornerCoords[],
367 : const MathVector<refDim> vLocIP[],
368 : const size_t nip,
369 : LocalVector* u,
370 : const MathMatrix<refDim, dim>* vJT = NULL) const
371 : {
372 : // reference object id
373 0 : const ReferenceObjectID roid = elem->reference_object_id();
374 :
375 : // memory for shapes
376 : std::vector<number> vShape;
377 :
378 : // multiindices of element
379 : std::vector<DoFIndex> ind;
380 :
381 0 : for (int d = 0; d < dim; ++d)
382 : {
383 : try
384 : {
385 : // get trial space
386 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
387 0 : LocalFiniteElementProvider::get<refDim>(roid, m_vlfeID[d]);
388 :
389 0 : m_spGridFct->dof_indices(elem, m_vfct[d], ind);
390 :
391 : // loop ips
392 0 : for(size_t ip = 0; ip < nip; ++ip)
393 : {
394 : // evaluate at shapes at ip
395 0 : rTrialSpace.shapes(vShape, vLocIP[ip]);
396 :
397 : // compute solution at integration point
398 0 : vValue[ip][d] = 0.0;
399 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
400 : {
401 0 : const number valSH = DoFRef(*m_spGridFct, ind[sh]);
402 0 : vValue[ip][d] += valSH * vShape[sh];
403 : }
404 : }
405 : }
406 0 : UG_CATCH_THROW("ExplicitGridFunctionVector: Shape Function Set missing for"
407 : " Reference Object: " << roid << ", Trial Space: "
408 : << m_vlfeID << ", refDim=" <<refDim);
409 : }
410 0 : }
411 :
412 0 : bool continuous() const
413 : {
414 0 : for(int i = 0; i < dim; ++i)
415 0 : if(!LocalFiniteElementProvider::continuous(m_vlfeID[i]))
416 : return false;
417 : return true;
418 : }
419 :
420 : };
421 :
422 : /**
423 : * Base class for the CplUserData-interface for the gradient of a grid function.
424 : *
425 : * This base class provides the CplUserData interface for the gradient of a GridFunction
426 : * at integration points. It can be used in particular to specify constant (or explicit)
427 : * import parameters for discretizations. Note that this interface does NOT compute
428 : * derivatives (so, assuming that the data are constant and do not depend on the
429 : * current guess of the solution). For the gradient of the current solution in
430 : * discretizations, so that the derivatives w.r.t. the DoFs are taken into account, use
431 : * the GridFunctionGradientData class.
432 : */
433 : template <typename TGridFunction>
434 : class ExplicitGridFunctionGradient
435 : : public StdExplicitGridFunctionData<ExplicitGridFunctionGradient<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction >
436 : {
437 : typedef StdExplicitGridFunctionData<ExplicitGridFunctionGradient<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction > base_type;
438 :
439 : using base_type::m_spGridFct;
440 :
441 : size_t m_fct; ///< component of the function
442 : LFEID m_lfeID; ///< type of the shape functions
443 : std::map<std::string, double> m_diffCoeffMap;
444 :
445 : public:
446 : // world dimension of grid function
447 : static const int dim = TGridFunction::dim;
448 :
449 : /// constructor
450 0 : ExplicitGridFunctionGradient
451 : (
452 : SmartPtr<TGridFunction> gf, ///< grid function
453 : const char* cmp ///< its component
454 : )
455 0 : : base_type(gf), m_diffCoeffMap()
456 : {
457 : // get function id of name
458 0 : m_fct = m_spGridFct->fct_id_by_name(cmp);
459 :
460 : // check that function exists
461 0 : if(m_fct >= m_spGridFct->num_fct())
462 0 : UG_THROW("ExplicitGridFunctionGradient: Function space does not contain"
463 : " a function with name " << cmp << ".");
464 :
465 : // local finite element id
466 0 : m_lfeID = m_spGridFct->local_finite_element_id(m_fct);
467 0 : }
468 :
469 : // evaluate gradient
470 : template <int refDim>
471 0 : void evaluate(MathVector<dim> vValue[],
472 : const MathVector<dim> vGlobIP[],
473 : number time, int si,
474 : GridObject* elem,
475 : const MathVector<dim> vCornerCoords[],
476 : const MathVector<refDim> vLocIP[],
477 : const size_t nip,
478 : LocalVector* u,
479 : const MathMatrix<refDim, dim>* vJT = NULL) const
480 : {
481 : // reference object id
482 0 : const ReferenceObjectID roid = elem->reference_object_id();
483 :
484 : // get reference element mapping by reference object id
485 0 : std::vector<MathMatrix<refDim, dim> > vJTTmp(nip);
486 0 : if(vJT == NULL){
487 : try{
488 : DimReferenceMapping<refDim, dim>& mapping
489 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
490 :
491 : // compute transformation matrices
492 0 : mapping.jacobian_transposed(&(vJTTmp[0]), vLocIP, nip);
493 :
494 : // store tmp Gradient
495 : vJT = &(vJTTmp[0]);
496 0 : }UG_CATCH_THROW("ExplicitGridFunctionGradient: failed.");
497 : }
498 :
499 : // scale with coeficient
500 0 : const int subsetInd = m_spGridFct->domain()->subset_handler()->get_subset_index(elem);
501 0 : const char* subsetName = m_spGridFct->domain()->subset_handler()->get_subset_name(subsetInd);
502 :
503 0 : double diffCoeff = get_subset_coeff(std::string(subsetName));
504 :
505 : // get trial space
506 : try{
507 : const LocalShapeFunctionSet<refDim>& rTrialSpace =
508 0 : LocalFiniteElementProvider::get<refDim>(roid, m_lfeID);
509 :
510 : // storage for shape function at ip
511 : std::vector<MathVector<refDim> > vLocGrad;
512 : MathVector<refDim> locGrad;
513 :
514 : // Reference Mapping
515 : MathMatrix<dim, refDim> JTInv;
516 :
517 : // loop ips
518 0 : for(size_t ip = 0; ip < nip; ++ip)
519 : {
520 : // evaluate at shapes at ip
521 0 : rTrialSpace.grads(vLocGrad, vLocIP[ip]);
522 :
523 : // get multiindices of element
524 : std::vector<DoFIndex > ind;
525 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
526 :
527 : // compute grad at ip
528 : VecSet(locGrad, 0.0);
529 0 : for(size_t sh = 0; sh < vLocGrad.size(); ++sh)
530 : {
531 0 : const number valSH = DoFRef( *m_spGridFct, ind[sh]);
532 : VecScaleAppend(locGrad, valSH, vLocGrad[sh]);
533 : }
534 :
535 0 : Inverse(JTInv, vJT[ip]);
536 0 : MatVecMult(vValue[ip], JTInv, locGrad);
537 :
538 : // scale with diff coeff (TODO: matrix)
539 : VecScale(vValue[ip], vValue[ip], diffCoeff);
540 : }
541 0 : }
542 0 : UG_CATCH_THROW("ExplicitGridFunctionGradient: Shape Function Set missing for"
543 : " Reference Object: "<<roid<<", Trial Space: "
544 : <<m_lfeID<<", refDim="<<refDim);
545 0 : }
546 :
547 0 : bool continuous() const { return false; }
548 :
549 0 : void add_subset_coeff(const std::string &key, double val)
550 : {
551 0 : m_diffCoeffMap.insert(std::pair<std::string, double>(key, val));
552 0 : }
553 :
554 0 : double get_subset_coeff(const std::string &key) const
555 : {
556 : std::map<std::string, double>::const_iterator it = m_diffCoeffMap.find(key);
557 0 : double val = (it != m_diffCoeffMap.end()) ? it->second : 1.0;
558 0 : return val;
559 : }
560 : };
561 :
562 :
563 : } // end namespace ug
564 :
565 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_USER_DATA__ */
|
