Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
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 : /*
34 : * Classes for user data that compute normals or project vectors to low-dimensional
35 : * subsets. Note that these classes are mainly implemented for visualization purposes
36 : * and cannot be used for coupling.
37 : */
38 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__ONC_USER_DATA__
39 : #define __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__ONC_USER_DATA__
40 :
41 : #include <vector>
42 :
43 : // ug4 headers
44 : #include "common/common.h"
45 : #include "common/math/ugmath.h"
46 : #include "lib_grid/grid_objects/grid_dim_traits.h"
47 : #include "lib_disc/domain_traits.h"
48 : #include "lib_disc/domain_util.h"
49 : #include "lib_disc/common/geometry_util.h"
50 : #include "lib_disc/reference_element/reference_mapping_provider.h"
51 : #include "lib_disc/spatial_disc/user_data/std_user_data.h"
52 :
53 : namespace ug {
54 :
55 : /// Projection to the outer normal for a given vector user data
56 : /**
57 : * This class implements the UserData for evaluation of the normal component of vectors
58 : * computed by the other user data object. The normal component is evaluation only for
59 : * WDim-1-dimensional subsets - sides of WDim-dimensional elements.
60 : *
61 : * ToDo: Note that the LocalVector provided for the evaluate function is not properly
62 : * used. Thus, the passed UserData objects should not use it.
63 : *
64 : */
65 : template <typename TDomain>
66 : class OutNormCmp
67 : : public StdUserData<OutNormCmp<TDomain>, MathVector<TDomain::dim>, TDomain::dim, void, UserData<MathVector<TDomain::dim>, TDomain::dim, void> >
68 : {
69 : /// the world dimension
70 : static const int dim = TDomain::dim;
71 :
72 : /// the domain type
73 : typedef TDomain domain_type;
74 :
75 : /// the grid type
76 : typedef typename TDomain::grid_type grid_type;
77 :
78 : /// the original data type
79 : typedef MathVector<dim> vec_type;
80 :
81 : /// the original data
82 : SmartPtr<UserData<vec_type, dim, void> > m_spData;
83 :
84 : /// the domain
85 : SmartPtr<domain_type> m_spDomain;
86 :
87 : /// the subset group for the inner part
88 : SubsetGroup m_ssGrp;
89 :
90 : public:
91 :
92 : /// constructor
93 0 : OutNormCmp
94 : (
95 : SmartPtr<domain_type> spDomain, ///< the domain
96 : SmartPtr<UserData<vec_type, dim, void> > spData, ///< the original data
97 : const char * ss_names ///< the subsets
98 : )
99 0 : : m_spData (spData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
100 : {
101 : // Parse the subset names:
102 : std::vector<std::string> vssNames;
103 : try
104 : {
105 0 : TokenizeString (ss_names, vssNames);
106 0 : for (size_t k = 0; k < vssNames.size (); k++)
107 0 : RemoveWhitespaceFromString (vssNames [k]);
108 : m_ssGrp.clear ();
109 0 : m_ssGrp.add (vssNames);
110 0 : } UG_CATCH_THROW ("SubsetIndicatorUserData: Failed to parse subset names.");
111 0 : }
112 :
113 : /// constructor
114 0 : OutNormCmp
115 : (
116 : SmartPtr<domain_type> spDomain, ///< the domain
117 : SmartPtr<UserData<vec_type, dim, void> > spData ///< the original data
118 : )
119 0 : : m_spData (spData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
120 0 : {}
121 :
122 : /// Indicator functions are discontinuous
123 0 : virtual bool continuous () const {return false;}
124 :
125 : /// Returns true to get the grid element in the evaluation routine
126 0 : virtual bool requires_grid_fct () const {return m_spData->requires_grid_fct ();}
127 : // Remark: Note that actually we have not enought data for that. This dependence is neglected.
128 :
129 : /// Evaluator
130 : template <int refDim>
131 0 : inline void evaluate
132 : (
133 : vec_type vValue[],
134 : const MathVector<dim> vGlobIP [],
135 : number time,
136 : int si,
137 : GridObject * elem,
138 : const MathVector<dim> vCornerCoords [],
139 : const MathVector<refDim> vLocIP [],
140 : const size_t nip,
141 : LocalVector * u,
142 : const MathMatrix<refDim, dim> * vJT = NULL
143 : ) const
144 : {
145 : if (refDim == dim)
146 : {
147 0 : (* m_spData) (vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
148 0 : return;
149 : }
150 : if (refDim != dim - 1)
151 : {
152 0 : UG_THROW ("OutNormCmp: Wrong dimensionality of the subset.");
153 : }
154 :
155 : // Get the full-dim. element associated with the given side
156 : typedef typename grid_dim_traits<dim>::grid_base_object fd_elem_type;
157 : fd_elem_type * fd_elem = NULL;
158 : int fd_si = -1;
159 :
160 : typedef typename Grid::traits<fd_elem_type>::secure_container fd_elem_secure_container_type;
161 : fd_elem_secure_container_type fd_elem_list;
162 0 : grid_type& grid = * (grid_type*) (m_spDomain->grid().get ());
163 0 : grid.associated_elements (fd_elem_list, elem);
164 :
165 0 : if (m_ssGrp.empty ()) // no subsets specified; we assume, elem should be the bnd of the whole domain
166 : {
167 0 : if (fd_elem_list.size () == 1)
168 : {
169 : fd_elem = fd_elem_list[0];
170 0 : fd_si = m_spDomain->subset_handler()->get_subset_index (fd_elem);
171 : }
172 : // otherwise this is no boundary of the domain and we return 0
173 : }
174 : else
175 : {
176 0 : for (size_t k = 0; k < fd_elem_list.size (); k++)
177 : {
178 : fd_elem_type * e = fd_elem_list[k];
179 0 : int e_si = m_ssGrp.subset_handler()->get_subset_index (e);
180 0 : if (m_ssGrp.contains (e_si))
181 : {
182 0 : if (fd_elem != NULL) // i.e. e is already the second one
183 : { // this is no boundary of the subset; return 0
184 : fd_elem = NULL;
185 : break;
186 : }
187 : fd_elem = e;
188 : fd_si = e_si;
189 : }
190 : }
191 : }
192 0 : if (fd_elem == NULL)
193 : { // this is no bondary, return 0
194 0 : for (size_t ip = 0; ip < nip; ip++) vValue[ip] = 0;
195 : return;
196 : }
197 :
198 : // Get the coordinates of the corners of the full-dim. element
199 0 : std::vector<MathVector<dim> > fd_elem_corner_coords (domain_traits<dim>::MaxNumVerticesOfElem);
200 0 : CollectCornerCoordinates (fd_elem->base_object_id (), fd_elem_corner_coords,
201 : *fd_elem, *m_spDomain, true);
202 :
203 : // Convert the local ip coordinates (use the global coordinates as they are the same)
204 0 : const ReferenceObjectID fd_roid = fd_elem->reference_object_id ();
205 : DimReferenceMapping<dim, dim> & fd_map
206 : = ReferenceMappingProvider::get<dim, dim> (fd_roid, fd_elem_corner_coords);
207 0 : std::vector<MathVector<dim> > fd_loc_ip (nip);
208 0 : for (size_t ip = 0; ip < nip; ip++)
209 0 : fd_map.global_to_local(fd_loc_ip[ip], vGlobIP[ip]);
210 :
211 : // Call the original UserData, get the vectors
212 0 : (* m_spData) (vValue, vGlobIP, time, fd_si, fd_elem, &(fd_elem_corner_coords[0]),
213 : &(fd_loc_ip[0]), nip, u);
214 : //TODO: Note that u is here a dummy parameter: We do not have enough data for it.
215 :
216 : // Project the vectors
217 0 : const ReferenceObjectID roid = elem->reference_object_id ();
218 : DimReferenceMapping<refDim, dim> & map
219 : = ReferenceMappingProvider::get<refDim, dim> (roid, vCornerCoords);
220 0 : for (size_t ip = 0; ip < nip; ip++)
221 : {
222 : MathMatrix<dim, refDim> J;
223 0 : MathVector<dim> p (vValue[ip]);
224 :
225 0 : map.jacobian (J, vLocIP[ip]);
226 0 : OrthogProjectVec (p, J);
227 : vValue[ip] -= p;
228 : }
229 : };
230 :
231 : /// This function should not be used
232 0 : void operator()
233 : (
234 : vec_type & vValue,
235 : const MathVector<dim> & globIP,
236 : number time,
237 : int si
238 : )
239 : const
240 : {
241 0 : UG_THROW("OutNormCmp: Element required for evaluation, but not passed. Cannot evaluate.");
242 : }
243 :
244 : /// This function should not be used
245 0 : void operator()
246 : (
247 : vec_type vValue [],
248 : const MathVector<dim> vGlobIP [],
249 : number time,
250 : int si,
251 : const size_t nip
252 : ) const
253 : {
254 0 : UG_THROW("OutNormCmp: Element required for evaluation, but not passed. Cannot evaluate.");
255 : }
256 :
257 : };
258 :
259 : /// Scaled projection to the outer normal for a given vector user data
260 : /**
261 : * This class implements the UserData for evaluation of the normal component of vectors
262 : * computed by the other user data object. The normal component is evaluation only for
263 : * WDim-1-dimensional subsets - sides of WDim-dimensional elements. It is scaled by
264 : * a furthe (scalar) UserData.
265 : *
266 : * ToDo: Note that the LocalVector provided for the evaluate function is not properly
267 : * used. Thus, the passed UserData objects should not use it.
268 : *
269 : */
270 : template <typename TDomain>
271 : class ScaledOutNormCmp
272 : : public StdUserData<ScaledOutNormCmp<TDomain>, MathVector<TDomain::dim>, TDomain::dim, void, UserData<MathVector<TDomain::dim>, TDomain::dim, void> >
273 : {
274 : /// the world dimension
275 : static const int dim = TDomain::dim;
276 :
277 : /// the domain type
278 : typedef TDomain domain_type;
279 :
280 : /// the grid type
281 : typedef typename TDomain::grid_type grid_type;
282 :
283 : /// the original data type
284 : typedef MathVector<dim> vec_type;
285 :
286 : /// the original vector field data
287 : SmartPtr<UserData<vec_type, dim, void> > m_spVecData;
288 :
289 : /// the scaling data
290 : SmartPtr<UserData<number, dim, void> > m_spScalData;
291 :
292 : /// the domain
293 : SmartPtr<domain_type> m_spDomain;
294 :
295 : /// the subset group for the inner part
296 : SubsetGroup m_ssGrp;
297 :
298 : public:
299 :
300 : /// constructor
301 0 : ScaledOutNormCmp
302 : (
303 : SmartPtr<domain_type> spDomain, ///< the domain
304 : SmartPtr<UserData<number, dim, void> > spScalData, ///< the scaling data
305 : SmartPtr<UserData<vec_type, dim, void> > spVecData, ///< the original vector field data
306 : const char * ss_names ///< the subsets
307 : )
308 0 : : m_spVecData (spVecData), m_spScalData (spScalData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
309 : {
310 : // Parse the subset names:
311 : std::vector<std::string> vssNames;
312 : try
313 : {
314 0 : TokenizeString (ss_names, vssNames);
315 0 : for (size_t k = 0; k < vssNames.size (); k++)
316 0 : RemoveWhitespaceFromString (vssNames [k]);
317 : m_ssGrp.clear ();
318 0 : m_ssGrp.add (vssNames);
319 0 : } UG_CATCH_THROW ("SubsetIndicatorUserData: Failed to parse subset names.");
320 0 : }
321 :
322 : /// constructor
323 0 : ScaledOutNormCmp
324 : (
325 : SmartPtr<domain_type> spDomain, ///< the domain
326 : SmartPtr<UserData<number, dim, void> > spScalData, ///< the scaling data
327 : SmartPtr<UserData<vec_type, dim, void> > spVecData ///< the original vector field data
328 : )
329 0 : : m_spVecData (spVecData), m_spScalData (spScalData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
330 0 : {}
331 :
332 : /// Indicator functions are discontinuous
333 0 : virtual bool continuous () const {return false;}
334 :
335 : /// Returns true to get the grid element in the evaluation routine
336 0 : virtual bool requires_grid_fct () const {return m_spVecData->requires_grid_fct () || m_spScalData->requires_grid_fct ();}
337 : // Remark: Note that actually we have not enought data for that. This dependence is neglected.
338 :
339 : /// Evaluator
340 : template <int refDim>
341 0 : inline void evaluate
342 : (
343 : vec_type vValue[],
344 : const MathVector<dim> vGlobIP [],
345 : number time,
346 : int si,
347 : GridObject * elem,
348 : const MathVector<dim> vCornerCoords [],
349 : const MathVector<refDim> vLocIP [],
350 : const size_t nip,
351 : LocalVector * u,
352 : const MathMatrix<refDim, dim> * vJT = NULL
353 : ) const
354 : {
355 : if (refDim == dim)
356 : {
357 0 : std::vector<number> scaling (nip);
358 0 : (* m_spVecData) (vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
359 0 : (* m_spScalData) (&(scaling[0]), vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
360 0 : for (size_t ip = 0; ip < nip; ip++)
361 0 : vValue[ip] *= scaling[ip];
362 : return;
363 0 : }
364 : if (refDim != dim - 1)
365 : {
366 0 : UG_THROW ("ScaledOutNormCmp: Wrong dimensionality of the subset.");
367 : }
368 :
369 : // Get the full-dim. element associated with the given side
370 : typedef typename grid_dim_traits<dim>::grid_base_object fd_elem_type;
371 : fd_elem_type * fd_elem = NULL;
372 : int fd_si = -1;
373 :
374 : typedef typename Grid::traits<fd_elem_type>::secure_container fd_elem_secure_container_type;
375 : fd_elem_secure_container_type fd_elem_list;
376 0 : grid_type& grid = * (grid_type*) (m_spDomain->grid().get ());
377 0 : grid.associated_elements (fd_elem_list, elem);
378 :
379 0 : if (m_ssGrp.empty ()) // no subsets specified; we assume, elem should be the bnd of the whole domain
380 : {
381 0 : if (fd_elem_list.size () == 1)
382 : {
383 : fd_elem = fd_elem_list[0];
384 0 : fd_si = m_spDomain->subset_handler()->get_subset_index (fd_elem);
385 : }
386 : // otherwise this is no boundary of the domain and we return 0
387 : }
388 : else
389 : {
390 0 : for (size_t k = 0; k < fd_elem_list.size (); k++)
391 : {
392 : fd_elem_type * e = fd_elem_list[k];
393 0 : int e_si = m_ssGrp.subset_handler()->get_subset_index (e);
394 0 : if (m_ssGrp.contains (e_si))
395 : {
396 0 : if (fd_elem != NULL) // i.e. e is already the second one
397 : { // this is no boundary of the subset; return 0
398 : fd_elem = NULL;
399 : break;
400 : }
401 : fd_elem = e;
402 : fd_si = e_si;
403 : }
404 : }
405 : }
406 0 : if (fd_elem == NULL)
407 : { // this is no bondary, return 0
408 0 : for (size_t ip = 0; ip < nip; ip++) vValue[ip] = 0;
409 : return;
410 : }
411 :
412 : // Get the coordinates of the corners of the full-dim. element
413 0 : std::vector<MathVector<dim> > fd_elem_corner_coords (domain_traits<dim>::MaxNumVerticesOfElem);
414 0 : CollectCornerCoordinates (fd_elem->base_object_id (), fd_elem_corner_coords,
415 : *fd_elem, *m_spDomain, true);
416 :
417 : // Convert the local ip coordinates (use the global coordinates as they are the same)
418 0 : const ReferenceObjectID fd_roid = fd_elem->reference_object_id ();
419 : DimReferenceMapping<dim, dim> & fd_map
420 : = ReferenceMappingProvider::get<dim, dim> (fd_roid, fd_elem_corner_coords);
421 0 : std::vector<MathVector<dim> > fd_loc_ip (nip);
422 0 : for (size_t ip = 0; ip < nip; ip++)
423 0 : fd_map.global_to_local(fd_loc_ip[ip], vGlobIP[ip]);
424 :
425 : // Call the original UserData, get the vectors
426 0 : std::vector<number> scaling (nip);
427 0 : (* m_spVecData) (vValue, vGlobIP, time, fd_si, fd_elem, &(fd_elem_corner_coords[0]),
428 : &(fd_loc_ip[0]), nip, u);
429 0 : (* m_spScalData) (&(scaling[0]), vGlobIP, time, fd_si, fd_elem, &(fd_elem_corner_coords[0]),
430 : &(fd_loc_ip[0]), nip, u);
431 : //TODO: Note that u is here a dummy parameter: We do not have enough data for it.
432 :
433 : // Scale and project the vectors
434 0 : const ReferenceObjectID roid = elem->reference_object_id ();
435 : DimReferenceMapping<refDim, dim> & map
436 : = ReferenceMappingProvider::get<refDim, dim> (roid, vCornerCoords);
437 0 : for (size_t ip = 0; ip < nip; ip++)
438 : {
439 0 : vValue[ip] *= scaling[ip];
440 :
441 : MathMatrix<dim, refDim> J;
442 : MathVector<dim> p (vValue[ip]);
443 :
444 0 : map.jacobian (J, vLocIP[ip]);
445 0 : OrthogProjectVec (p, J);
446 : vValue[ip] -= p;
447 : }
448 0 : };
449 :
450 : /// This function should not be used
451 0 : void operator()
452 : (
453 : vec_type & vValue,
454 : const MathVector<dim> & globIP,
455 : number time,
456 : int si
457 : )
458 : const
459 : {
460 0 : UG_THROW("ScaledOutNormCmp: Element required for evaluation, but not passed. Cannot evaluate.");
461 : }
462 :
463 : /// This function should not be used
464 0 : void operator()
465 : (
466 : vec_type vValue [],
467 : const MathVector<dim> vGlobIP [],
468 : number time,
469 : int si,
470 : const size_t nip
471 : ) const
472 : {
473 0 : UG_THROW("ScaledOutNormCmp: Element required for evaluation, but not passed. Cannot evaluate.");
474 : }
475 :
476 : };
477 :
478 : /// Scaled flux for a given vector user data
479 : /**
480 : * This class implements the UserData for evaluation of the flux of a vector field
481 : * computed by the other user data object. The flux is evaluation only for
482 : * WDim-1-dimensional subsets - sides of WDim-dimensional elements. It is scaled by
483 : * a furthe (scalar) UserData.
484 : *
485 : * ToDo: Note that the LocalVector provided for the evaluate function is not properly
486 : * used. Thus, the passed UserData objects should not use it.
487 : *
488 : */
489 : template <typename TDomain>
490 : class ScaledFluxData
491 : : public StdUserData<ScaledFluxData<TDomain>, number, TDomain::dim, void, UserData<number, TDomain::dim, void> >
492 : {
493 : /// the world dimension
494 : static const int dim = TDomain::dim;
495 :
496 : /// the domain type
497 : typedef TDomain domain_type;
498 :
499 : /// the grid type
500 : typedef typename TDomain::grid_type grid_type;
501 :
502 : /// the original data type
503 : typedef MathVector<dim> vec_type;
504 :
505 : /// the original vector field data
506 : SmartPtr<UserData<vec_type, dim, void> > m_spVecData;
507 :
508 : /// the scaling data
509 : SmartPtr<UserData<number, dim, void> > m_spScalData;
510 :
511 : /// the domain
512 : SmartPtr<domain_type> m_spDomain;
513 :
514 : /// the subset group for the inner part
515 : SubsetGroup m_ssGrp;
516 :
517 : public:
518 :
519 : /// constructor
520 0 : ScaledFluxData
521 : (
522 : SmartPtr<domain_type> spDomain, ///< the domain
523 : SmartPtr<UserData<number, dim, void> > spScalData, ///< the scaling data
524 : SmartPtr<UserData<vec_type, dim, void> > spVecData, ///< the original vector field data
525 : const char * ss_names ///< the subsets
526 : )
527 0 : : m_spVecData (spVecData), m_spScalData (spScalData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
528 : {
529 : // Parse the subset names:
530 : std::vector<std::string> vssNames;
531 : try
532 : {
533 0 : TokenizeString (ss_names, vssNames);
534 0 : for (size_t k = 0; k < vssNames.size (); k++)
535 0 : RemoveWhitespaceFromString (vssNames [k]);
536 : m_ssGrp.clear ();
537 0 : m_ssGrp.add (vssNames);
538 0 : } UG_CATCH_THROW ("SubsetIndicatorUserData: Failed to parse subset names.");
539 0 : }
540 :
541 : /// constructor
542 0 : ScaledFluxData
543 : (
544 : SmartPtr<domain_type> spDomain, ///< the domain
545 : SmartPtr<UserData<number, dim, void> > spScalData, ///< the scaling data
546 : SmartPtr<UserData<vec_type, dim, void> > spVecData ///< the original vector field data
547 : )
548 0 : : m_spVecData (spVecData), m_spScalData (spScalData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
549 0 : {}
550 :
551 : /// Indicator functions are discontinuous
552 0 : virtual bool continuous () const {return false;}
553 :
554 : /// Returns true to get the grid element in the evaluation routine
555 0 : virtual bool requires_grid_fct () const {return m_spVecData->requires_grid_fct () || m_spScalData->requires_grid_fct ();}
556 : // Remark: Note that actually we have not enought data for that. This dependence is neglected.
557 :
558 : /// Evaluator
559 : template <int refDim>
560 0 : inline void evaluate
561 : (
562 : number vValue[],
563 : const MathVector<dim> vGlobIP [],
564 : number time,
565 : int si,
566 : GridObject * elem,
567 : const MathVector<dim> vCornerCoords [],
568 : const MathVector<refDim> vLocIP [],
569 : const size_t nip,
570 : LocalVector * u,
571 : const MathMatrix<refDim, dim> * vJT = NULL
572 : ) const
573 : {
574 : if (refDim != dim - 1)
575 : {
576 0 : UG_THROW ("ScaledFluxData: Wrong dimensionality of the subset.");
577 : }
578 :
579 : // Get the full-dim. element associated with the given side
580 : typedef typename grid_dim_traits<dim>::grid_base_object fd_elem_type;
581 : fd_elem_type * fd_elem = NULL;
582 : int fd_si = -1;
583 :
584 : typedef typename Grid::traits<fd_elem_type>::secure_container fd_elem_secure_container_type;
585 : fd_elem_secure_container_type fd_elem_list;
586 0 : grid_type& grid = * (grid_type*) (m_spDomain->grid().get ());
587 0 : grid.associated_elements (fd_elem_list, elem);
588 :
589 0 : if (m_ssGrp.empty ()) // no subsets specified; we assume, elem should be the bnd of the whole domain
590 : {
591 0 : if (fd_elem_list.size () == 1)
592 : {
593 : fd_elem = fd_elem_list[0];
594 0 : fd_si = m_spDomain->subset_handler()->get_subset_index (fd_elem);
595 : }
596 : // otherwise this is no boundary of the domain and we return 0
597 : }
598 : else
599 : {
600 0 : for (size_t k = 0; k < fd_elem_list.size (); k++)
601 : {
602 : fd_elem_type * e = fd_elem_list[k];
603 0 : int e_si = m_ssGrp.subset_handler()->get_subset_index (e);
604 0 : if (m_ssGrp.contains (e_si))
605 : {
606 0 : if (fd_elem != NULL) // i.e. e is already the second one
607 : { // this is no boundary of the subset; return 0
608 : fd_elem = NULL;
609 : break;
610 : }
611 : fd_elem = e;
612 : fd_si = e_si;
613 : }
614 : }
615 : }
616 0 : if (fd_elem == NULL)
617 : { // this is no bondary, return 0
618 0 : for (size_t ip = 0; ip < nip; ip++) vValue[ip] = 0;
619 : return;
620 : }
621 :
622 : // Get the coordinates of the corners of the full-dim. element
623 0 : std::vector<MathVector<dim> > fd_elem_corner_coords (domain_traits<dim>::MaxNumVerticesOfElem);
624 0 : CollectCornerCoordinates (fd_elem->base_object_id (), fd_elem_corner_coords,
625 : *fd_elem, *m_spDomain, true);
626 :
627 : // Convert the local ip coordinates (use the global coordinates as they are the same)
628 0 : const ReferenceObjectID fd_roid = fd_elem->reference_object_id ();
629 : DimReferenceMapping<dim, dim> & fd_map
630 : = ReferenceMappingProvider::get<dim, dim> (fd_roid, fd_elem_corner_coords);
631 0 : std::vector<MathVector<dim> > fd_loc_ip (nip);
632 0 : for (size_t ip = 0; ip < nip; ip++)
633 0 : fd_map.global_to_local(fd_loc_ip[ip], vGlobIP[ip]);
634 :
635 : // Get the normal to the face:
636 : MathVector<dim> normal;
637 0 : ElementNormal<dim> (elem->reference_object_id (), normal, vCornerCoords);
638 0 : VecScale (normal, normal, 1 / VecLength (normal));
639 0 : for (size_t co = 0; co < fd_elem->num_vertices (); co++)
640 : {
641 : MathVector<dim> s_vec;
642 : VecSubtract (s_vec, fd_elem_corner_coords[co], vCornerCoords[0]);
643 0 : if (VecDot (s_vec, normal) > 1e-12 * VecLength (s_vec))
644 : { // this is the inner normal, invert its direction
645 : VecScale (normal, normal, -1);
646 0 : break;
647 : }
648 : }
649 :
650 : // Call the original UserData, get the vectors
651 0 : std::vector<MathVector<dim> > vecfield (nip);
652 0 : std::vector<number> scaling (nip);
653 0 : (* m_spVecData) (&(vecfield[0]), vGlobIP, time, fd_si, fd_elem, &(fd_elem_corner_coords[0]),
654 : &(fd_loc_ip[0]), nip, u);
655 0 : (* m_spScalData) (&(scaling[0]), vGlobIP, time, fd_si, fd_elem, &(fd_elem_corner_coords[0]),
656 : &(fd_loc_ip[0]), nip, u);
657 : //TODO: Note that u is here a dummy parameter: We do not have enough data for it.
658 :
659 : // Scale and project the vectors
660 0 : for (size_t ip = 0; ip < nip; ip++)
661 0 : vValue[ip] = scaling[ip] * VecDot (vecfield[ip], normal);
662 0 : };
663 :
664 : /// This function should not be used
665 0 : void operator()
666 : (
667 : number & vValue,
668 : const MathVector<dim> & globIP,
669 : number time,
670 : int si
671 : )
672 : const
673 : {
674 0 : UG_THROW("ScaledFluxData: Element required for evaluation, but not passed. Cannot evaluate.");
675 : }
676 :
677 : /// This function should not be used
678 0 : void operator()
679 : (
680 : number vValue [],
681 : const MathVector<dim> vGlobIP [],
682 : number time,
683 : int si,
684 : const size_t nip
685 : ) const
686 : {
687 0 : UG_THROW("ScaledFluxData: Element required for evaluation, but not passed. Cannot evaluate.");
688 : }
689 :
690 : };
691 :
692 : } // end namespace ug
693 :
694 : #endif // __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__ONC_USER_DATA__
695 :
696 : /* End of File */
|
