Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter, Andreas Vogel
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__GRID_FUNCTION_GLOBAL_USER_DATA__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_GLOBAL_USER_DATA__
35 :
36 : #include "common/common.h"
37 :
38 : #include "lib_grid/tools/subset_group.h"
39 :
40 : #include "lib_disc/common/function_group.h"
41 : #include "lib_disc/common/groups_util.h"
42 : #include "lib_disc/quadrature/quadrature.h"
43 : #include "lib_disc/domain_util.h"
44 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
45 : #include "lib_disc/spatial_disc/user_data/std_glob_pos_data.h"
46 : #include "lib_disc/reference_element/reference_mapping_provider.h"
47 : #include "lib_grid/algorithms/space_partitioning/lg_ntree.h"
48 :
49 : #include <math.h> /* fabs */
50 :
51 : namespace ug{
52 :
53 :
54 :
55 : template <typename TGridFunction, int elemDim = TGridFunction::dim>
56 : class GlobalGridFunctionNumberData
57 : : public StdGlobPosData<GlobalGridFunctionNumberData<TGridFunction, elemDim>, number, TGridFunction::dim>
58 : {
59 : public:
60 : /// world dimension of grid function
61 : static const int dim = TGridFunction::dim;
62 : typedef typename TGridFunction::template dim_traits<elemDim>::grid_base_object element_t;
63 :
64 : private:
65 : /// grid function
66 : SmartPtr<TGridFunction> m_spGridFct;
67 :
68 : /// component of function
69 : size_t m_fct;
70 :
71 : /// local finite element id
72 : LFEID m_lfeID;
73 :
74 : typedef lg_ntree<dim, dim, element_t> tree_t;
75 : tree_t m_tree;
76 :
77 : public:
78 : /// constructor
79 0 : GlobalGridFunctionNumberData(SmartPtr<TGridFunction> spGridFct, const char* cmp)
80 : : m_spGridFct(spGridFct),
81 0 : m_tree(*spGridFct->domain()->grid(), spGridFct->domain()->position_attachment())
82 : {
83 : //this->set_functions(cmp);
84 :
85 : // get function id of name
86 0 : m_fct = spGridFct->fct_id_by_name(cmp);
87 :
88 : // check that function exists
89 0 : if(m_fct >= spGridFct->num_fct())
90 0 : UG_THROW("GridFunctionNumberData: Function space does not contain"
91 : " a function with name " << cmp << ".");
92 :
93 : // local finite element id
94 0 : m_lfeID = spGridFct->local_finite_element_id(m_fct);
95 :
96 :
97 :
98 : // const MGSubsetHandler& ssh = *m_spGridFct->domain()->subset_handler();
99 :
100 0 : SubsetGroup ssGrp(m_spGridFct->domain()->subset_handler());
101 0 : ssGrp.add_all();
102 :
103 : std::vector<size_t> subsetsOfGridFunction;
104 : std::vector<element_t*> elemsWithGridFunctions;
105 :
106 : typename TGridFunction::template dim_traits<elemDim>::const_iterator iterEnd, iter;
107 :
108 :
109 0 : for(size_t si = 0; si < ssGrp.size(); si++){
110 0 : if( spGridFct->is_def_in_subset(m_fct, si) )
111 0 : subsetsOfGridFunction.push_back(si);
112 : }
113 :
114 :
115 0 : for(size_t i = 0; i<subsetsOfGridFunction.size(); i++){
116 0 : size_t si = subsetsOfGridFunction[i];
117 0 : iter = spGridFct->template begin<element_t>(si);
118 0 : iterEnd = spGridFct->template end<element_t>(si);
119 :
120 0 : for(;iter!=iterEnd; ++iter){
121 0 : element_t *elem = *iter;
122 0 : elemsWithGridFunctions.push_back(elem);
123 : }
124 : }
125 :
126 : // m_tree.create_tree(spGridFct->template begin<element_t>(si),
127 : // spGridFct->template end<element_t>(si));
128 :
129 0 : m_tree.create_tree(elemsWithGridFunctions.begin(), elemsWithGridFunctions.end());
130 :
131 0 : };
132 :
133 0 : virtual ~GlobalGridFunctionNumberData() {}
134 :
135 0 : virtual bool continuous() const
136 : {
137 0 : return LocalFiniteElementProvider::continuous(m_lfeID);
138 : }
139 :
140 : /// to full-fill UserData-Interface
141 0 : inline void evaluate(number& value, const MathVector<dim>& x, number time, int si) const
142 : {
143 0 : if(!evaluate(value, x))
144 0 : UG_THROW("For function "<<m_fct<<" couldn't find an element containing the specified point: " << x);
145 0 : }
146 :
147 0 : inline number evaluate(const MathVector<dim>& x) const
148 : {
149 : number value;
150 0 : if(!evaluate(value, x))
151 0 : UG_THROW("For function "<<m_fct<<" couldn't find an element containing the specified point: " << x);
152 0 : return value;
153 : }
154 :
155 : /// evaluates the data at a given point, returns false if point not found
156 0 : inline bool evaluate(number& value, const MathVector<dim>& x) const
157 : {
158 0 : element_t* elem = NULL;
159 : //try{
160 :
161 0 : if(!FindContainingElement(elem, m_tree, x)){
162 : return false;
163 : }
164 :
165 : // get corners of element
166 : std::vector<MathVector<dim> > vCornerCoords;
167 0 : CollectCornerCoordinates(vCornerCoords, *elem, *m_spGridFct->domain());
168 :
169 : // reference object id
170 0 : const ReferenceObjectID roid = elem->reference_object_id();
171 :
172 : // get local position of DoF
173 : DimReferenceMapping<elemDim, dim>& map
174 : = ReferenceMappingProvider::get<elemDim, dim>(roid, vCornerCoords);
175 : MathVector<elemDim> locPos;
176 : VecSet(locPos, 0.5);
177 0 : map.global_to_local(locPos, x);
178 :
179 : // evaluate at shapes at ip
180 : const LocalShapeFunctionSet<elemDim>& rTrialSpace =
181 0 : LocalFiniteElementProvider::get<elemDim>(roid, m_lfeID);
182 : std::vector<number> vShape;
183 0 : rTrialSpace.shapes(vShape, locPos);
184 :
185 : // get multiindices of element
186 : std::vector<DoFIndex> ind;
187 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
188 :
189 : // compute solution at integration point
190 0 : value = 0.0;
191 0 : for(size_t sh = 0; sh < vShape.size(); ++sh)
192 : {
193 0 : const number valSH = DoFRef(*m_spGridFct, ind[sh]);
194 0 : value += valSH * vShape[sh];
195 : }
196 :
197 : // point is found
198 : return true;
199 : //}
200 : //UG_CATCH_THROW("GlobalGridFunctionNumberData: Evaluation failed."
201 : // << "Point: " << x << ", Element: "
202 : // << ElementDebugInfo(*m_spGridFct->domain()->grid(), elem));
203 0 : }
204 : /// evaluate value on all procs, throws when no containing element is found
205 0 : inline void evaluate_global(number& value, const MathVector<dim>& x) const
206 : {
207 : // evaluate at this proc
208 : bool bFound = this->try_evaluate_global(value, x);
209 :
210 0 : if(!bFound)
211 0 : UG_THROW("Couldn't find an element containing the specified point: " << x);
212 0 : }
213 :
214 :
215 : /// evaluate value on all procs, does not throw an error if element is not found.
216 : /// returns truth value whether element was found or not
217 : inline bool try_evaluate_global(number& value, const MathVector<dim>& x) const
218 : {
219 : // evaluate at this proc
220 0 : bool bFound = this->evaluate(value, x);
221 :
222 : #ifdef UG_PARALLEL
223 : // share value between all procs
224 : pcl::ProcessCommunicator com;
225 : int numFound = (bFound ? 1 : 0);
226 : numFound = com.allreduce(numFound, PCL_RO_SUM);
227 :
228 : // get overall value
229 : if(!bFound) value = 0.0;
230 : number globValue = com.allreduce(value, PCL_RO_SUM) / numFound;
231 :
232 : // return false if no process contains an element with coordinates x
233 : if(numFound == 0)
234 : return false;
235 :
236 :
237 :
238 : // check correctness for continuous spaces
239 : // note: if the point is found more than one it is located on the
240 : // boundary of some element. thus, if the space is continuous, those
241 : // values should match on all procs.
242 : if(bFound)
243 : if(LocalFiniteElementProvider::continuous(m_lfeID)){
244 : if( fabs(value) > 1e-10 && fabs((globValue - value) / value) > 1e-8)
245 : UG_THROW("Global mean "<<globValue<<" != local value "<<value);
246 : }
247 :
248 : // set as global value
249 : value = globValue;
250 : // if we get here, we found an element on some process
251 : bFound = true;
252 :
253 : #endif
254 : return bFound;
255 : }
256 :
257 : // evaluates at given position
258 0 : number evaluate_global(std::vector<number> vPos)
259 : {
260 0 : if((int)vPos.size() != dim)
261 0 : UG_THROW("Expected "<<dim<<" components, but given "<<vPos.size());
262 :
263 : MathVector<dim> x;
264 0 : for(int i = 0; i < dim; i++) x[i] = vPos[i];
265 :
266 : number value;
267 0 : evaluate_global(value, x);
268 :
269 0 : return value;
270 : }
271 : };
272 :
273 :
274 :
275 : template <typename TGridFunction>
276 : class GlobalGridFunctionGradientData
277 : : public StdGlobPosData<GlobalGridFunctionGradientData<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction::dim>
278 : {
279 : public:
280 : /// world dimension of grid function
281 : static const int dim = TGridFunction::dim;
282 : typedef typename TGridFunction::element_type element_t;
283 :
284 : private:
285 : /// grid function
286 : SmartPtr<TGridFunction> m_spGridFct;
287 :
288 : /// component of function
289 : size_t m_fct;
290 :
291 : /// local finite element id
292 : LFEID m_lfeID;
293 :
294 : typedef lg_ntree<dim, dim, element_t> tree_t;
295 : tree_t m_tree;
296 :
297 : public:
298 : /// constructor
299 0 : GlobalGridFunctionGradientData(SmartPtr<TGridFunction> spGridFct, const char* cmp)
300 : : m_spGridFct(spGridFct),
301 0 : m_tree(*spGridFct->domain()->grid(), spGridFct->domain()->position_attachment())
302 : {
303 : //this->set_functions(cmp);
304 :
305 : // get function id of name
306 0 : m_fct = spGridFct->fct_id_by_name(cmp);
307 :
308 : // check that function exists
309 0 : if(m_fct >= spGridFct->num_fct())
310 0 : UG_THROW("GridFunctionNumberData: Function space does not contain"
311 : " a function with name " << cmp << ".");
312 :
313 : // local finite element id
314 0 : m_lfeID = spGridFct->local_finite_element_id(m_fct);
315 :
316 :
317 :
318 : // const MGSubsetHandler& ssh = *m_spGridFct->domain()->subset_handler();
319 :
320 0 : SubsetGroup ssGrp(m_spGridFct->domain()->subset_handler());
321 0 : ssGrp.add_all();
322 :
323 : std::vector<size_t> subsetsOfGridFunction;
324 : std::vector<element_t*> elemsWithGridFunctions;
325 :
326 : typename TGridFunction::const_element_iterator iterEnd, iter;
327 :
328 :
329 0 : for(size_t si = 0; si < ssGrp.size(); si++){
330 0 : if( spGridFct->is_def_in_subset(m_fct, si) )
331 0 : subsetsOfGridFunction.push_back(si);
332 : }
333 :
334 :
335 0 : for(size_t i = 0; i<subsetsOfGridFunction.size(); i++){
336 0 : size_t si = subsetsOfGridFunction[i];
337 0 : iter = spGridFct->template begin<element_t>(si);
338 0 : iterEnd = spGridFct->template end<element_t>(si);
339 :
340 0 : for(;iter!=iterEnd; ++iter){
341 0 : element_t *elem = *iter;
342 0 : elemsWithGridFunctions.push_back(elem);
343 : }
344 : }
345 :
346 : // m_tree.create_tree(spGridFct->template begin<element_t>(si),
347 : // spGridFct->template end<element_t>(si));
348 :
349 0 : m_tree.create_tree(elemsWithGridFunctions.begin(), elemsWithGridFunctions.end());
350 :
351 0 : };
352 :
353 0 : virtual ~GlobalGridFunctionGradientData() {}
354 :
355 0 : virtual bool continuous() const
356 : {
357 0 : return false;
358 : }
359 :
360 : /// to full-fill UserData-Interface
361 0 : inline void evaluate(MathVector<dim>& value, const MathVector<dim>& x, number time, int si) const
362 : {
363 0 : if(!evaluate(value, x))
364 0 : UG_THROW("For function "<<m_fct<<" couldn't find an element containing the specified point: " << x);
365 0 : }
366 :
367 : /// evaluates the data at a given point, returns false if point not found
368 0 : inline bool evaluate(MathVector<dim>& value, const MathVector<dim>& x) const
369 : {
370 : static const int refDim = dim;
371 :
372 0 : element_t* elem = NULL;
373 : try{
374 :
375 0 : if(!FindContainingElement(elem, m_tree, x)){
376 : return false;
377 : }
378 :
379 : // get corners of element
380 : std::vector<MathVector<dim> > vCornerCoords;
381 0 : CollectCornerCoordinates(vCornerCoords, *elem, *m_spGridFct->domain());
382 :
383 : // reference object id
384 0 : const ReferenceObjectID roid = elem->reference_object_id();
385 :
386 : // get local position of DoF
387 : DimReferenceMapping<dim, dim>& map
388 : = ReferenceMappingProvider::get<dim, dim>(roid, vCornerCoords);
389 : MathVector<dim> locPos;
390 : VecSet(locPos, 0.5);
391 0 : map.global_to_local(locPos, x);
392 :
393 :
394 : MathMatrix<refDim, dim> JT;
395 : try{
396 : DimReferenceMapping<refDim, dim>& mapping
397 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
398 :
399 : // compute transformation matrices
400 0 : mapping.jacobian_transposed(JT, x);
401 :
402 0 : }UG_CATCH_THROW("GlobalGridFunctionGradientData: failed.");
403 :
404 : // evaluate at shapes at ip
405 : const LocalShapeFunctionSet<dim>& rTrialSpace =
406 0 : LocalFiniteElementProvider::get<dim>(roid, m_lfeID);
407 : std::vector<MathVector<refDim> > vLocGrad;
408 0 : rTrialSpace.grads(vLocGrad, locPos);
409 :
410 :
411 : // Reference Mapping
412 : MathMatrix<dim, refDim> JTInv;
413 0 : RightInverse (JTInv, JT);
414 :
415 : // get multiindices of element
416 : std::vector<DoFIndex> ind;
417 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
418 :
419 : // storage for shape function at ip
420 : MathVector<refDim> locGrad;
421 :
422 : // compute grad at ip
423 : VecSet(locGrad, 0.0);
424 0 : for(size_t sh = 0; sh < vLocGrad.size(); ++sh)
425 : {
426 0 : const number valSH = DoFRef( *m_spGridFct, ind[sh]);
427 : VecScaleAppend(locGrad, valSH, vLocGrad[sh]);
428 : }
429 :
430 : // transform to global space
431 : MatVecMult(value, JTInv, locGrad);
432 :
433 : // point is found
434 : return true;
435 0 : }
436 0 : UG_CATCH_THROW("GlobalGridFunctionGradientData: Evaluation failed."
437 : << "Point: " << x << ", Element: "
438 : << ElementDebugInfo(*m_spGridFct->domain()->grid(), elem));
439 : }
440 :
441 : /// evaluate value on all procs
442 0 : inline void evaluate_global(MathVector<dim>& value, const MathVector<dim>& x) const
443 : {
444 : // evaluate at this proc
445 0 : bool bFound = this->evaluate(value, x);
446 :
447 : // \todo: (optinal) check for evaluation on other procs
448 :
449 0 : if(!bFound)
450 0 : UG_THROW("Couldn't find an element containing the specified point: " << x);
451 0 : }
452 :
453 : // evaluates at given position
454 0 : std::vector<number> evaluate_global(std::vector<number> vPos)
455 : {
456 0 : if((int)vPos.size() != dim)
457 0 : UG_THROW("Expected "<<dim<<" components, but given "<<vPos.size());
458 :
459 : MathVector<dim> x;
460 0 : for(int i = 0; i < dim; i++) x[i] = vPos[i];
461 :
462 : MathVector<dim> value;
463 0 : evaluate_global(value, x);
464 :
465 0 : for(int i = 0; i < dim; i++) vPos[i] = value[i];
466 0 : return vPos;
467 : }
468 : };
469 :
470 : } // end namespace ug
471 :
472 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_GLOBAL_USER_DATA__ */
|
