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 :
205 : /// evaluate value on all procs
206 0 : inline void evaluate_global(number& value, const MathVector<dim>& x) const
207 : {
208 : // evaluate at this proc
209 0 : bool bFound = this->evaluate(value, x);
210 :
211 : #ifdef UG_PARALLEL
212 : // share value between all procs
213 : pcl::ProcessCommunicator com;
214 : int numFound = (bFound ? 1 : 0);
215 : numFound = com.allreduce(numFound, PCL_RO_SUM);
216 :
217 : // get overall value
218 : if(!bFound) value = 0.0;
219 : number globValue = com.allreduce(value, PCL_RO_SUM) / numFound;
220 :
221 : if(numFound == 0)
222 : UG_THROW("Point "<<x<<" not found on all "<<pcl::NumProcs()<<" procs.");
223 :
224 : // check correctness for continuous spaces
225 : // note: if the point is found more than one it is located on the
226 : // boundary of some element. thus, if the space is continuous, those
227 : // values should match on all procs.
228 : if(bFound)
229 : if(LocalFiniteElementProvider::continuous(m_lfeID)){
230 : if( fabs(value) > 1e-10 && fabs((globValue - value) / value) > 1e-8)
231 : UG_THROW("Global mean "<<globValue<<" != local value "<<value);
232 : }
233 :
234 : // set as global value
235 : value = globValue;
236 : bFound = true;
237 : #endif
238 :
239 0 : if(!bFound)
240 0 : UG_THROW("Couldn't find an element containing the specified point: " << x);
241 0 : }
242 :
243 : // evaluates at given position
244 0 : number evaluate_global(std::vector<number> vPos)
245 : {
246 0 : if((int)vPos.size() != dim)
247 0 : UG_THROW("Expected "<<dim<<" components, but given "<<vPos.size());
248 :
249 : MathVector<dim> x;
250 0 : for(int i = 0; i < dim; i++) x[i] = vPos[i];
251 :
252 : number value;
253 0 : evaluate_global(value, x);
254 :
255 0 : return value;
256 : }
257 : };
258 :
259 :
260 :
261 : template <typename TGridFunction>
262 : class GlobalGridFunctionGradientData
263 : : public StdGlobPosData<GlobalGridFunctionGradientData<TGridFunction>, MathVector<TGridFunction::dim>, TGridFunction::dim>
264 : {
265 : public:
266 : /// world dimension of grid function
267 : static const int dim = TGridFunction::dim;
268 : typedef typename TGridFunction::element_type element_t;
269 :
270 : private:
271 : /// grid function
272 : SmartPtr<TGridFunction> m_spGridFct;
273 :
274 : /// component of function
275 : size_t m_fct;
276 :
277 : /// local finite element id
278 : LFEID m_lfeID;
279 :
280 : typedef lg_ntree<dim, dim, element_t> tree_t;
281 : tree_t m_tree;
282 :
283 : public:
284 : /// constructor
285 0 : GlobalGridFunctionGradientData(SmartPtr<TGridFunction> spGridFct, const char* cmp)
286 : : m_spGridFct(spGridFct),
287 0 : m_tree(*spGridFct->domain()->grid(), spGridFct->domain()->position_attachment())
288 : {
289 : //this->set_functions(cmp);
290 :
291 : // get function id of name
292 0 : m_fct = spGridFct->fct_id_by_name(cmp);
293 :
294 : // check that function exists
295 0 : if(m_fct >= spGridFct->num_fct())
296 0 : UG_THROW("GridFunctionNumberData: Function space does not contain"
297 : " a function with name " << cmp << ".");
298 :
299 : // local finite element id
300 0 : m_lfeID = spGridFct->local_finite_element_id(m_fct);
301 :
302 :
303 :
304 : // const MGSubsetHandler& ssh = *m_spGridFct->domain()->subset_handler();
305 :
306 0 : SubsetGroup ssGrp(m_spGridFct->domain()->subset_handler());
307 0 : ssGrp.add_all();
308 :
309 : std::vector<size_t> subsetsOfGridFunction;
310 : std::vector<element_t*> elemsWithGridFunctions;
311 :
312 : typename TGridFunction::const_element_iterator iterEnd, iter;
313 :
314 :
315 0 : for(size_t si = 0; si < ssGrp.size(); si++){
316 0 : if( spGridFct->is_def_in_subset(m_fct, si) )
317 0 : subsetsOfGridFunction.push_back(si);
318 : }
319 :
320 :
321 0 : for(size_t i = 0; i<subsetsOfGridFunction.size(); i++){
322 0 : size_t si = subsetsOfGridFunction[i];
323 0 : iter = spGridFct->template begin<element_t>(si);
324 0 : iterEnd = spGridFct->template end<element_t>(si);
325 :
326 0 : for(;iter!=iterEnd; ++iter){
327 0 : element_t *elem = *iter;
328 0 : elemsWithGridFunctions.push_back(elem);
329 : }
330 : }
331 :
332 : // m_tree.create_tree(spGridFct->template begin<element_t>(si),
333 : // spGridFct->template end<element_t>(si));
334 :
335 0 : m_tree.create_tree(elemsWithGridFunctions.begin(), elemsWithGridFunctions.end());
336 :
337 0 : };
338 :
339 0 : virtual ~GlobalGridFunctionGradientData() {}
340 :
341 0 : virtual bool continuous() const
342 : {
343 0 : return false;
344 : }
345 :
346 : /// to full-fill UserData-Interface
347 0 : inline void evaluate(MathVector<dim>& value, const MathVector<dim>& x, number time, int si) const
348 : {
349 0 : if(!evaluate(value, x))
350 0 : UG_THROW("For function "<<m_fct<<" couldn't find an element containing the specified point: " << x);
351 0 : }
352 :
353 : /// evaluates the data at a given point, returns false if point not found
354 0 : inline bool evaluate(MathVector<dim>& value, const MathVector<dim>& x) const
355 : {
356 : static const int refDim = dim;
357 :
358 0 : element_t* elem = NULL;
359 : try{
360 :
361 0 : if(!FindContainingElement(elem, m_tree, x)){
362 : return false;
363 : }
364 :
365 : // get corners of element
366 : std::vector<MathVector<dim> > vCornerCoords;
367 0 : CollectCornerCoordinates(vCornerCoords, *elem, *m_spGridFct->domain());
368 :
369 : // reference object id
370 0 : const ReferenceObjectID roid = elem->reference_object_id();
371 :
372 : // get local position of DoF
373 : DimReferenceMapping<dim, dim>& map
374 : = ReferenceMappingProvider::get<dim, dim>(roid, vCornerCoords);
375 : MathVector<dim> locPos;
376 : VecSet(locPos, 0.5);
377 0 : map.global_to_local(locPos, x);
378 :
379 :
380 : MathMatrix<refDim, dim> JT;
381 : try{
382 : DimReferenceMapping<refDim, dim>& mapping
383 : = ReferenceMappingProvider::get<refDim, dim>(roid, vCornerCoords);
384 :
385 : // compute transformation matrices
386 0 : mapping.jacobian_transposed(JT, x);
387 :
388 0 : }UG_CATCH_THROW("GlobalGridFunctionGradientData: failed.");
389 :
390 : // evaluate at shapes at ip
391 : const LocalShapeFunctionSet<dim>& rTrialSpace =
392 0 : LocalFiniteElementProvider::get<dim>(roid, m_lfeID);
393 : std::vector<MathVector<refDim> > vLocGrad;
394 0 : rTrialSpace.grads(vLocGrad, locPos);
395 :
396 :
397 : // Reference Mapping
398 : MathMatrix<dim, refDim> JTInv;
399 0 : RightInverse (JTInv, JT);
400 :
401 : // get multiindices of element
402 : std::vector<DoFIndex> ind;
403 0 : m_spGridFct->dof_indices(elem, m_fct, ind);
404 :
405 : // storage for shape function at ip
406 : MathVector<refDim> locGrad;
407 :
408 : // compute grad at ip
409 : VecSet(locGrad, 0.0);
410 0 : for(size_t sh = 0; sh < vLocGrad.size(); ++sh)
411 : {
412 0 : const number valSH = DoFRef( *m_spGridFct, ind[sh]);
413 : VecScaleAppend(locGrad, valSH, vLocGrad[sh]);
414 : }
415 :
416 : // transform to global space
417 : MatVecMult(value, JTInv, locGrad);
418 :
419 : // point is found
420 : return true;
421 0 : }
422 0 : UG_CATCH_THROW("GlobalGridFunctionGradientData: Evaluation failed."
423 : << "Point: " << x << ", Element: "
424 : << ElementDebugInfo(*m_spGridFct->domain()->grid(), elem));
425 : }
426 :
427 : /// evaluate value on all procs
428 0 : inline void evaluate_global(MathVector<dim>& value, const MathVector<dim>& x) const
429 : {
430 : // evaluate at this proc
431 0 : bool bFound = this->evaluate(value, x);
432 :
433 : // \todo: (optinal) check for evaluation on other procs
434 :
435 0 : if(!bFound)
436 0 : UG_THROW("Couldn't find an element containing the specified point: " << x);
437 0 : }
438 :
439 : // evaluates at given position
440 0 : std::vector<number> evaluate_global(std::vector<number> vPos)
441 : {
442 0 : if((int)vPos.size() != dim)
443 0 : UG_THROW("Expected "<<dim<<" components, but given "<<vPos.size());
444 :
445 : MathVector<dim> x;
446 0 : for(int i = 0; i < dim; i++) x[i] = vPos[i];
447 :
448 : MathVector<dim> value;
449 0 : evaluate_global(value, x);
450 :
451 0 : for(int i = 0; i < dim; i++) vPos[i] = value[i];
452 0 : return vPos;
453 : }
454 : };
455 :
456 : } // end namespace ug
457 :
458 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_GLOBAL_USER_DATA__ */
|
