Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: 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__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_HIGHER_ORDER_GEOMETRY__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_HIGHER_ORDER_GEOMETRY__
35 :
36 : // extern libraries
37 : #include <cmath>
38 : #include <map>
39 : #include <vector>
40 :
41 : // other ug4 modules
42 : #include "common/common.h"
43 :
44 : // library intern includes
45 : #include "lib_grid/tools/subset_handler_interface.h"
46 : #include "lib_disc/reference_element/reference_element.h"
47 : #include "lib_disc/reference_element/reference_element_traits.h"
48 : #include "lib_disc/reference_element/reference_mapping.h"
49 : #include "lib_disc/reference_element/reference_mapping_provider.h"
50 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
51 : #include "lib_disc/quadrature/gauss/gauss_quad.h"
52 : #include "fv_util.h"
53 : #include "fv_geom_base.h"
54 :
55 : namespace ug{
56 :
57 : ////////////////////////////////////////////////////////////////////////////////
58 : // FV Geometry for Reference Element Type (all orders)
59 : ////////////////////////////////////////////////////////////////////////////////
60 :
61 : /// Geometry and shape functions for any order Vertex-Centered Finite Volume
62 : /**
63 : * \tparam TOrder order
64 : * \tparam TElem Element type
65 : * \tparam TWorldDim (physical) world dimension
66 : * \tparam TQuadOrderSCVF integration order for scvf
67 : * \tparam TQuadOrderSCV integration order for scv
68 : */
69 : template < int TOrder, typename TElem, int TWorldDim,
70 : int TQuadOrder = TOrder + 1>
71 : class FVGeometry : public FVGeometryBase
72 : {
73 : private:
74 : /// small abbreviation for order
75 : static const int p = TOrder;
76 :
77 : public:
78 : /// type of element
79 : typedef TElem elem_type;
80 :
81 : /// type of reference element
82 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
83 :
84 : /// dimension of reference element
85 : static const int dim = ref_elem_type::dim;
86 :
87 : /// dimension of world
88 : static const int worldDim = TWorldDim;
89 :
90 : public:
91 : /// order
92 : static const int order = TOrder;
93 :
94 : /// number of subelements
95 : static const size_t numSubElem = Pow<p, dim>::value;
96 :
97 :
98 : /// traits used
99 : typedef fv1_traits<ref_elem_type, worldDim> traits;
100 :
101 : /// type of SubControlVolumeFace
102 : typedef typename traits::scvf_type scvf_type;
103 :
104 : /// number of SCVF per SubElement
105 : static const size_t numSCVFPerSubElem = ref_elem_type::numEdges;
106 :
107 : /// number of SubControlVolumeFaces
108 : static const size_t numSCVF = numSubElem * numSCVFPerSubElem;
109 :
110 : /// quadrature order
111 : static const int quadOrderSCVF = TQuadOrder;
112 :
113 : /// type of quadrature rule
114 : typedef GaussQuadrature<scvf_type, quadOrderSCVF> scvf_quad_rule_type;
115 :
116 : /// number of scvf ip
117 : static const size_t numSCVFIP = scvf_quad_rule_type::nip * numSCVF;
118 :
119 :
120 : /// type of SubControlVolume
121 : typedef typename traits::scv_type scv_type;
122 :
123 : /// number of SCV per SubElement
124 : static const size_t numSCVPerSubElem = ref_elem_type::numCorners;
125 :
126 : /// number of SubControlVolumes
127 : static const size_t numSCV = numSubElem * numSCVPerSubElem;
128 :
129 : /// quadrature order
130 : static const int quadOrderSCV = TQuadOrder;
131 :
132 : /// type of quadrature rule
133 : typedef GaussQuadrature<scv_type, quadOrderSCV> scv_quad_rule_type;
134 :
135 : /// number of scv ip
136 : static const size_t numSCVIP = scv_quad_rule_type::nip * numSCV;
137 :
138 :
139 : /// type of Shape function used
140 : typedef LagrangeLSFS<ref_elem_type, p> local_shape_fct_set_type;
141 :
142 : /// number of shape functions
143 : static const size_t nsh = local_shape_fct_set_type::nsh;
144 :
145 : /// Hanging node flag: this Geometry does not support hanging nodes
146 : static const bool usesHangingNodes = false;
147 :
148 : /// flag indicating if local data may change
149 : static const bool staticLocalData = true;
150 :
151 : public:
152 : /// Sub-Control Volume Face structure
153 : class SCVF
154 : {
155 : public:
156 : /// Number of integration points
157 : static const size_t nip = scvf_quad_rule_type::nip;
158 :
159 : /// Number of corners of scvf
160 : static const size_t numCo = traits::NumCornersOfSCVF;
161 :
162 : public:
163 0 : SCVF() {}
164 :
165 : /// index of SubControlVolume on one side of the scvf
166 0 : inline size_t from() const {return From;}
167 :
168 : /// index of SubControlVolume on one side of the scvf
169 0 : inline size_t to() const {return To;}
170 :
171 : /// number of integration points on scvf
172 0 : inline size_t num_ip() const {return nip;}
173 :
174 : /// integration weight
175 0 : inline number weight(size_t ip) const
176 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];;}
177 :
178 : /// local integration point of scvf
179 0 : inline const MathVector<dim>& local_ip(size_t ip) const
180 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
181 :
182 : /// global integration point of scvf
183 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
184 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
185 :
186 : /// normal on scvf (points direction "from"->"to"), normalized
187 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
188 :
189 : /// Transposed Inverse of Jacobian in integration point
190 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
191 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
192 :
193 : /// Determinant of Jacobian in integration point
194 0 : inline number detJ(size_t ip) const
195 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
196 :
197 : /// number of shape functions
198 0 : inline size_t num_sh() const {return nsh;}
199 :
200 : /// value of shape function i in integration point
201 0 : inline number shape(size_t ip, size_t sh) const
202 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
203 :
204 : /// vector of shape functions in ip point
205 0 : inline const number* shape_vector(size_t ip) const
206 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip];}
207 :
208 : /// value of local gradient of shape function i in integration point
209 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
210 : {UG_ASSERT(sh < num_sh(), "Invalid index");
211 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
212 :
213 : /// vector of local gradients in ip point
214 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
215 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip];}
216 :
217 : /// value of global gradient of shape function i in integration point
218 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
219 : {UG_ASSERT(sh < num_sh(), "Invalid index");
220 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
221 :
222 : /// vector of global gradients in ip point
223 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
224 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip];}
225 :
226 : /// number of corners, that bound the scvf
227 0 : inline size_t num_corners() const {return numCo;}
228 :
229 : /// return local corner number i
230 0 : inline const MathVector<dim>& local_corner(size_t co) const
231 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
232 :
233 : /// return global corner number i
234 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
235 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
236 :
237 : private:
238 : // let outer class access private members
239 : friend class FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>;
240 :
241 : // This scvf separates the scv with the ids given in "from" and "to"
242 : // The computed normal points in direction from->to
243 : size_t From, To;
244 :
245 : // The normal on the SCVF pointing (from -> to)
246 : MathVector<worldDim> Normal; // normal (incl. area)
247 :
248 : // ordering is:
249 : // 1D: edgeMidPoint
250 : // 2D: edgeMidPoint, CenterOfElement
251 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
252 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
253 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
254 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scvf
255 :
256 : // scvf part
257 : MathVector<dim> vLocalIP[nip]; // local integration point
258 : MathVector<worldDim> vGlobalIP[nip]; // global integration point
259 : const number* vWeight; // weight at ip
260 :
261 : // shapes and derivatives
262 : number vvShape[nip][nsh]; // shapes at ip
263 : MathVector<dim> vvLocalGrad[nip][nsh]; // local grad at ip
264 : MathVector<worldDim> vvGlobalGrad[nip][nsh]; // global grad at ip
265 : MathMatrix<worldDim,dim> vJtInv[nip]; // Jacobian transposed at ip
266 : number vDetJ[nip]; // Jacobian det at ip
267 : number vDetJMap[nip]; // Jacobian det at ip
268 : };
269 :
270 : /// sub control volume structure
271 : class SCV
272 : {
273 : public:
274 : /// Number of integration points
275 : static const size_t nip = scv_quad_rule_type::nip;
276 :
277 : /// Number of corners of scvf
278 : static const size_t numCo = traits::NumCornersOfSCV;
279 :
280 : public:
281 0 : SCV() {};
282 :
283 : /// number of corners, that bound the scvf
284 0 : inline size_t num_corners() const {return numCo;}
285 :
286 : /// return local corner number i
287 0 : inline const MathVector<dim>& local_corner(size_t co) const
288 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
289 :
290 : /// return glbal corner number i
291 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
292 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
293 :
294 : /// node id that this scv is associated to
295 0 : inline size_t node_id() const {return nodeId;}
296 :
297 : /// number of integration points
298 0 : inline size_t num_ip() const {return nip;}
299 :
300 : /// weigth of integration point
301 0 : inline number weight(size_t ip) const
302 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];}
303 :
304 : /// local integration point of scv
305 0 : inline const MathVector<dim>& local_ip(size_t ip) const
306 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
307 :
308 : /// global integration point
309 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
310 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
311 :
312 : /// Transposed Inverse of Jacobian in integration point
313 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
314 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vJtInv[ip];}
315 :
316 : /// Determinant of Jacobian in integration point
317 0 : inline number detJ(size_t ip) const
318 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vDetJ[ip];}
319 :
320 : /// number of shape functions
321 0 : inline size_t num_sh() const {return nsh;}
322 :
323 : /// value of shape function i in integration point
324 0 : inline number shape(size_t ip, size_t sh) const
325 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vvShape[ip][sh];}
326 :
327 : /// vector of shape functions in ip point
328 0 : inline const number* shape_vector(size_t ip) const
329 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vvShape[ip];}
330 :
331 : /// value of local gradient of shape function i in integration point
332 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
333 : {UG_ASSERT(sh < num_sh(), "Invalid index");
334 0 : UG_ASSERT(ip<num_ip(),"Wring index."); return vvLocalGrad[ip][sh];}
335 :
336 : /// vector of local gradients in ip point
337 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
338 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vvLocalGrad[ip];}
339 :
340 : /// value of global gradient of shape function i in integration point
341 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
342 : {UG_ASSERT(sh < num_sh(), "Invalid index");
343 0 : UG_ASSERT(ip<num_ip(),"Wring index."); return vvGlobalGrad[ip][sh];}
344 :
345 : /// vector of global gradients in ip point
346 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
347 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vvGlobalGrad[ip];}
348 :
349 : private:
350 : // let outer class access private members
351 : friend class FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>;
352 :
353 : // node id of associated node
354 : size_t nodeId;
355 :
356 : // scv part
357 : MathVector<dim> vLocalIP[nip]; // local integration point
358 : MathVector<worldDim> vGlobalIP[nip]; // global intergration point
359 : const number* vWeight; // weight at ip
360 :
361 : // local and global positions of this element
362 : MathVector<dim> vLocPos[numCo]; // local position of node
363 : MathVector<worldDim> vGloPos[numCo]; // global position of node
364 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scv
365 :
366 : // shapes and derivatives
367 : number vvShape[nip][nsh]; // shapes at ip
368 : MathVector<dim> vvLocalGrad[nip][nsh]; // local grad at ip
369 : MathVector<worldDim> vvGlobalGrad[nip][nsh]; // global grad at ip
370 : MathMatrix<worldDim,dim> vJtInv[nip]; // Jacobian transposed at ip
371 : number vDetJ[nip]; // Jacobian det at ip
372 : number vDetJMap[nip]; // Jacobian det at ip for scv integral map
373 : };
374 :
375 : /// boundary face
376 : class BF
377 : {
378 : public:
379 : /// number of integration points
380 : static const size_t nip = scvf_quad_rule_type::nip;
381 :
382 : /// Number of corners of bf
383 : static const size_t numCo = traits::NumCornersOfSCVF;
384 :
385 : public:
386 0 : BF() {}
387 :
388 : /// index of SubControlVolume of the bf
389 0 : inline size_t node_id() const {return nodeId;}
390 :
391 : /// outer normal on bf. Norm is equal to area
392 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
393 :
394 : /// volume of bf
395 0 : inline number volume() const {return Vol;}
396 :
397 : /// number of integration points on scvf
398 0 : inline size_t num_ip() const {return nip;}
399 :
400 : /// integration weight
401 0 : inline number weight(size_t ip) const
402 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip];}
403 :
404 : /// local integration point of scvf
405 0 : inline const MathVector<dim>& local_ip(size_t ip) const
406 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
407 :
408 : /// global integration point of scvf
409 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
410 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
411 :
412 : /// Transposed Inverse of Jacobian in integration point
413 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
414 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
415 :
416 : /// Determinant of Jacobian in integration point
417 0 : inline number detJ(size_t ip) const
418 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
419 :
420 : /// number of shape functions
421 0 : inline size_t num_sh() const {return nsh;}
422 :
423 : /// value of shape function i in integration point
424 0 : inline number shape(size_t ip, size_t sh) const
425 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
426 :
427 : /// vector of shape functions in ip point
428 0 : inline const number* shape_vector(size_t ip) const
429 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip];}
430 :
431 : /// value of local gradient of shape function i in integration point
432 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
433 : {UG_ASSERT(sh < num_sh(), "Invalid index");
434 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
435 :
436 : /// vector of local gradients in ip point
437 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
438 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip];}
439 :
440 : /// value of global gradient of shape function i in integration point
441 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
442 : {UG_ASSERT(sh < num_sh(), "Invalid index");
443 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
444 :
445 : /// vector of global gradients in ip point
446 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
447 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip];}
448 :
449 : /// number of corners, that bound the scvf
450 0 : inline size_t num_corners() const {return numCo;}
451 :
452 : /// return local corner number i
453 0 : inline const MathVector<dim>& local_corner(size_t co) const
454 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
455 :
456 : /// return glbal corner number i
457 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
458 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
459 :
460 : private:
461 : /// let outer class access private members
462 : friend class FVGeometry<TOrder, TElem, TWorldDim, TQuadOrder>;
463 :
464 : // id of scv this bf belongs to
465 : size_t nodeId;
466 :
467 : // ordering is:
468 : // 1D: edgeMidPoint
469 : // 2D: edgeMidPoint, CenterOfElement
470 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
471 : MathVector<dim> vLocPos[numCo]; // local corners of bf
472 : MathVector<worldDim> vGloPos[numCo]; // global corners of bf
473 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the bf
474 :
475 : // scvf part
476 : MathVector<dim> vLocalIP[nip]; // local integration point
477 : MathVector<worldDim> vGlobalIP[nip]; // global integration point
478 : const number* vWeight; // weight at ip
479 : MathVector<worldDim> Normal; // normal (incl. area)
480 : number Vol; // volume of bf
481 :
482 : // shapes and derivatives
483 : number vvShape[nip][nsh]; // shapes at ip
484 : MathVector<dim> vvLocalGrad[nip][nsh]; // local grad at ip
485 : MathVector<worldDim> vvGlobalGrad[nip][nsh]; // global grad at ip
486 : MathMatrix<worldDim,dim> vJtInv[nip]; // Jacobian transposed at ip
487 : number vDetJ[nip]; // Jacobian det at ip
488 : };
489 :
490 :
491 : public:
492 : /// construct object and initialize local values and sizes
493 : FVGeometry();
494 :
495 : /// update local data
496 : void update_local_data();
497 :
498 : /// update Geometry for roid
499 : void update_local(ReferenceObjectID roid,
500 : const LFEID& lfeID = LFEID(LFEID::LAGRANGE, worldDim, 1),
501 : size_t orderQuad = TQuadOrder);
502 :
503 : /// update data for given element
504 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
505 : const ISubsetHandler* ish = NULL);
506 :
507 : /// update boundary data for given element
508 : void update_boundary_faces(GridObject* elem,
509 : const MathVector<worldDim>* vCornerCoords,
510 : const ISubsetHandler* ish = NULL);
511 :
512 : /// number of SubControlVolumeFaces
513 0 : inline size_t num_scvf() const {return numSCVF;};
514 :
515 : /// const access to SubControlVolumeFace number i
516 0 : inline const SCVF& scvf(size_t i) const
517 0 : {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
518 :
519 : /// number of SubControlVolumes
520 0 : inline size_t num_scv() const {return numSCV;}
521 :
522 : /// const access to SubControlVolume number i
523 0 : inline const SCV& scv(size_t i) const
524 0 : {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
525 :
526 : /// number of shape functions
527 0 : inline size_t num_sh() const {return nsh;}
528 :
529 : public:
530 : /// returns number of all scvf ips
531 0 : size_t num_scvf_ips() const {return numSCVFIP;}
532 :
533 : /// returns all ips of scvf as they appear in scv loop
534 0 : const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
535 :
536 : /// returns all ips of scvf as they appear in scv loop
537 0 : const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
538 :
539 : /// returns number of all scv ips
540 0 : size_t num_scv_ips() const {return numSCVIP;}
541 :
542 : /// returns all ips of scv as they appear in scv loop
543 0 : const MathVector<dim>* scv_local_ips() const {return m_vLocSCV_IP;}
544 :
545 : /// returns all ips of scv as they appear in scv loop
546 0 : const MathVector<worldDim>* scv_global_ips() const {return m_vGlobSCV_IP;}
547 :
548 :
549 : protected:
550 : // global and local ips on SCVF
551 : MathVector<worldDim> m_vGlobSCVF_IP[numSCVFIP];
552 : MathVector<dim> m_vLocSCVF_IP[numSCVFIP];
553 :
554 : // global and local ips on SCVF
555 : MathVector<worldDim> m_vGlobSCV_IP[numSCVIP];
556 : MathVector<dim> m_vLocSCV_IP[numSCVIP];
557 :
558 : protected:
559 : // maximum number of geom objects in a dimension
560 : static const int maxMid = numSCVF + 1;
561 :
562 : // subelement
563 0 : struct SubElement
564 : {
565 : // shape number of corners of subelement
566 : size_t vDoFID[ref_elem_type::numCorners];
567 :
568 : // local and global geom object midpoints for each dimension
569 : // (most objects in 1 dim, i.e. number of edges, but +1 for 1D)
570 : MathVector<dim> vvLocMid[dim+1][maxMid];
571 : MathVector<worldDim> vvGloMid[dim+1][maxMid];
572 :
573 : // flag is subelement has boundary sides
574 : bool isBndElem;
575 :
576 : // -1 is no bnd side, >= 0 corresponding side of whole element
577 : std::vector<int> vElemBndSide;
578 : };
579 :
580 : /// subelements
581 : SubElement m_vSubElem[numSubElem];
582 :
583 : public:
584 : /// add subset that is interpreted as boundary subset.
585 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
586 :
587 : /// removes subset that is interpreted as boundary subset.
588 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
589 :
590 : /// reset all boundary subsets
591 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
592 :
593 : /// number of registered boundary subsets
594 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
595 :
596 : /// number of all boundary faces
597 0 : inline size_t num_bf() const
598 : {
599 : typename std::map<int, std::vector<BF> >::const_iterator it;
600 : size_t num = 0;
601 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
602 0 : num += (*it).second.size();
603 0 : return num;
604 : }
605 :
606 : /// number of boundary faces on subset 'subsetIndex'
607 0 : inline size_t num_bf(int si) const
608 : {
609 : typename std::map<int, std::vector<BF> >::const_iterator it;
610 : it = m_mapVectorBF.find(si);
611 0 : if(it == m_mapVectorBF.end()) return 0;
612 0 : else return (*it).second.size();
613 : }
614 :
615 : /// returns the boundary face i for subsetIndex
616 0 : inline const BF& bf(int si, size_t i) const
617 : {
618 : typename std::map<int, std::vector<BF> >::const_iterator it;
619 : it = m_mapVectorBF.find(si);
620 0 : if(it == m_mapVectorBF.end()) UG_THROW("FVGeom: No BndSubset"<<si);
621 0 : return (*it).second[i];
622 : }
623 :
624 : /// returns reference to vector of boundary faces for subsetIndex
625 0 : inline const std::vector<BF>& bf(int si) const
626 : {
627 : typename std::map<int, std::vector<BF> >::const_iterator it;
628 : it = m_mapVectorBF.find(si);
629 0 : if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
630 0 : return (*it).second;
631 : }
632 :
633 0 : void reset_curr_elem() {m_pElem = NULL;}
634 :
635 : protected:
636 : std::map<int, std::vector<BF> > m_mapVectorBF;
637 : std::vector<BF> m_vEmptyVectorBF;
638 :
639 : private:
640 : /// pointer to current element
641 : TElem* m_pElem;
642 :
643 : /// corners of reference element
644 : MathVector<dim> m_vLocCorner[ref_elem_type::numCorners];
645 :
646 : /// SubControlVolumeFaces
647 : SCVF m_vSCVF[numSCVF];
648 :
649 : /// SubControlVolumes
650 : SCV m_vSCV[numSCV];
651 :
652 : /// Reference Mapping
653 : ReferenceMapping<ref_elem_type, worldDim> m_rMapping;
654 :
655 : /// Reference Element
656 : const ref_elem_type& m_rRefElem;
657 :
658 : /// Shape function set
659 : const local_shape_fct_set_type& m_rTrialSpace;
660 :
661 : /// Quad Rule scvf
662 : const scvf_quad_rule_type& m_rSCVFQuadRule;
663 :
664 : /// Quad Rule scv
665 : const scv_quad_rule_type& m_rSCVQuadRule;
666 : };
667 :
668 : ////////////////////////////////////////////////////////////////////////////////
669 : // Dimension FV Geometry (all orders) DIM FV
670 : ////////////////////////////////////////////////////////////////////////////////
671 :
672 : /// Geometry and shape functions for any order Vertex-Centered Finite Volume
673 : /**
674 : * \tparam TDim reference element dim
675 : * \tparam TWorldDim (physical) world dimension
676 : */
677 : template <int TWorldDim, int TDim = TWorldDim>
678 : class DimFVGeometry : public FVGeometryBase
679 : {
680 : public:
681 : /// traits used
682 : typedef fv1_dim_traits<TDim, TWorldDim> traits;
683 :
684 : /// dimension of reference element
685 : static const int dim = TDim;
686 :
687 : /// dimension of world
688 : static const int worldDim = TWorldDim;
689 :
690 : public:
691 : /// type of SubControlVolumeFace
692 : typedef typename traits::scvf_type scvf_type;
693 :
694 : /// type of SubControlVolume
695 : typedef typename traits::scv_type scv_type;
696 :
697 : /// Hanging node flag: this Geometry does not support hanging nodes
698 : static const bool usesHangingNodes = false;
699 :
700 : /// flag indicating if local data may change
701 : static const bool staticLocalData = false;
702 :
703 : public:
704 : /// Sub-Control Volume Face structure
705 : class SCVF
706 : {
707 : public:
708 : /// Number of corners of scvf
709 : static const size_t numCo = traits::NumCornersOfSCVF;
710 :
711 : public:
712 0 : SCVF() {}
713 :
714 : /// index of SubControlVolume on one side of the scvf
715 0 : inline size_t from() const {return From;}
716 :
717 : /// index of SubControlVolume on one side of the scvf
718 0 : inline size_t to() const {return To;}
719 :
720 : /// number of integration points on scvf
721 0 : inline size_t num_ip() const {return nip;}
722 :
723 : /// integration weight
724 0 : inline number weight(size_t ip) const
725 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];}
726 :
727 : /// local integration point of scvf
728 0 : inline const MathVector<dim>& local_ip(size_t ip) const
729 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
730 :
731 : /// global integration point of scvf
732 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
733 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
734 :
735 : /// normal on scvf (points direction "from"->"to"), normalized
736 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
737 :
738 : /// Transposed Inverse of Jacobian in integration point
739 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
740 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
741 :
742 : /// Determinant of Jacobian in integration point
743 0 : inline number detJ(size_t ip) const
744 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
745 :
746 : /// number of shape functions
747 0 : inline size_t num_sh() const {return nsh;}
748 :
749 : /// value of shape function i in integration point
750 0 : inline number shape(size_t ip, size_t sh) const
751 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
752 :
753 : /// vector of shape functions in ip point
754 0 : inline const number* shape_vector(size_t ip) const
755 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvShape[ip][0];}
756 :
757 : /// value of local gradient of shape function i in integration point
758 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
759 : {UG_ASSERT(sh < num_sh(), "Invalid index");
760 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
761 :
762 : /// vector of local gradients in ip point
763 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
764 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvLocalGrad[ip][0];}
765 :
766 : /// value of global gradient of shape function i in integration point
767 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
768 : {UG_ASSERT(sh < num_sh(), "Invalid index");
769 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
770 :
771 : /// vector of gloabl gradients in ip point
772 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
773 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvGlobalGrad[ip][0];}
774 :
775 : /// number of corners, that bound the scvf
776 0 : inline size_t num_corners() const {return numCo;}
777 :
778 : /// return local corner number i
779 0 : inline const MathVector<dim>& local_corner(size_t co) const
780 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
781 :
782 : /// return glbal corner number i
783 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
784 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
785 :
786 : private:
787 : // let outer class access private members
788 : friend class DimFVGeometry<worldDim, dim>;
789 :
790 : // This scvf separates the scv with the ids given in "from" and "to"
791 : // The computed normal points in direction from->to
792 : size_t From, To;
793 :
794 : // The normal on the SCVF pointing (from -> to)
795 : MathVector<worldDim> Normal; // normal (incl. area)
796 :
797 : // ordering is:
798 : // 1D: edgeMidPoint
799 : // 2D: edgeMidPoint, CenterOfElement
800 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
801 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
802 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
803 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scvf
804 :
805 : // scvf part
806 : size_t nip;
807 : std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
808 : std::vector<MathVector<worldDim> > vGlobalIP; // global integration point (size: nip)
809 : const number* vWeight; // weight at ip
810 :
811 : // shapes and derivatives
812 : size_t nsh;
813 : std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
814 : std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
815 : std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
816 : std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
817 : std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
818 : std::vector<number> vDetJMap; // Jacobian det at ip (size: nip)
819 : };
820 :
821 : /// sub control volume structure
822 : class SCV
823 : {
824 : public:
825 : /// Number of corners of scvf
826 : static const size_t numCo = traits::NumCornersOfSCV;
827 :
828 : public:
829 0 : SCV() {};
830 :
831 : /// number of corners, that bound the scvf
832 0 : inline size_t num_corners() const {return numCo;}
833 :
834 : /// return local corner number i
835 0 : inline const MathVector<dim>& local_corner(size_t co) const
836 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
837 :
838 : /// return glbal corner number i
839 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
840 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
841 :
842 : /// node id that this scv is associated to
843 0 : inline size_t node_id() const {return nodeId;}
844 :
845 : /// number of integration points
846 0 : inline size_t num_ip() const {return nip;}
847 :
848 : /// weigth of integration point
849 0 : inline number weight(size_t ip) const
850 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip]*vDetJMap[ip];}
851 :
852 : /// local integration point of scv
853 0 : inline const MathVector<dim>& local_ip(size_t ip) const
854 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
855 :
856 : /// global integration point
857 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
858 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
859 :
860 : /// Transposed Inverse of Jacobian in integration point
861 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
862 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vJtInv[ip];}
863 :
864 : /// Determinant of Jacobian in integration point
865 0 : inline number detJ(size_t ip) const
866 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vDetJ[ip];}
867 :
868 : /// number of shape functions
869 0 : inline size_t num_sh() const {return nsh;}
870 :
871 : /// value of shape function i in integration point
872 0 : inline number shape(size_t ip, size_t sh) const
873 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return vvShape[ip][sh];}
874 :
875 : /// vector of shape functions in ip point
876 0 : inline const number* shape_vector(size_t ip) const
877 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return &vvShape[ip][0];}
878 :
879 : /// value of local gradient of shape function i in integration point
880 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
881 : {UG_ASSERT(sh < num_sh(), "Invalid index");
882 0 : UG_ASSERT(ip<num_ip(),"Wring index."); return vvLocalGrad[ip][sh];}
883 :
884 : /// vector of local gradients in ip point
885 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
886 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return &vvLocalGrad[ip][0];}
887 :
888 : /// value of global gradient of shape function i in integration point
889 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
890 : {UG_ASSERT(sh < num_sh(), "Invalid index");
891 0 : UG_ASSERT(ip<num_ip(),"Wring index."); return vvGlobalGrad[ip][sh];}
892 :
893 : /// vector of global gradients in ip point
894 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
895 0 : {UG_ASSERT(ip<num_ip(),"Wring index."); return &vvGlobalGrad[ip][0];}
896 :
897 : private:
898 : // let outer class access private members
899 : friend class DimFVGeometry<worldDim, dim>;
900 :
901 : // node id of associated node
902 : size_t nodeId;
903 :
904 : // local and global positions of this element
905 : MathVector<dim> vLocPos[numCo]; // local position of node
906 : MathVector<worldDim> vGloPos[numCo]; // global position of node
907 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scv
908 :
909 : // scv part
910 : size_t nip;
911 : std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
912 : std::vector<MathVector<worldDim> > vGlobalIP; // global intergration point (size: nip)
913 : const number* vWeight; // weight at ip
914 :
915 : // shapes and derivatives
916 : size_t nsh;
917 : std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
918 : std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
919 : std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
920 : std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
921 : std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
922 : std::vector<number> vDetJMap; // Jacobian det at ip (size: nip)
923 : };
924 :
925 : /// boundary face
926 : class BF
927 : {
928 : public:
929 : /// Number of corners of bf
930 : static const size_t numCo = traits::NumCornersOfSCVF;
931 :
932 : public:
933 0 : BF() {}
934 :
935 : /// index of SubControlVolume of the bf
936 0 : inline size_t node_id() const {return nodeId;}
937 :
938 : /// outer normal on bf. Norm is equal to area
939 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
940 :
941 : /// volume of bf
942 0 : inline number volume() const {return Vol;}
943 :
944 : /// number of integration points on scvf
945 0 : inline size_t num_ip() const {return nip;}
946 :
947 : /// integration weight
948 0 : inline number weight(size_t ip) const
949 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vWeight[ip];}
950 :
951 : /// local integration point of scvf
952 0 : inline const MathVector<dim>& local_ip(size_t ip) const
953 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
954 :
955 : /// global integration point of scvf
956 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
957 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
958 :
959 : /// Transposed Inverse of Jacobian in integration point
960 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
961 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
962 :
963 : /// Determinant of Jacobian in integration point
964 0 : inline number detJ(size_t ip) const
965 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
966 :
967 : /// number of shape functions
968 0 : inline size_t num_sh() const {return nsh;}
969 :
970 : /// value of shape function i in integration point
971 0 : inline number shape(size_t ip, size_t sh) const
972 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
973 :
974 : /// vector of shape functions in ip point
975 0 : inline const number* shape_vector(size_t ip) const
976 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvShape[ip][0];}
977 :
978 : /// value of local gradient of shape function i in integration point
979 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
980 : {UG_ASSERT(sh < num_sh(), "Invalid index");
981 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
982 :
983 : /// vector of local gradients in ip point
984 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
985 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvLocalGrad[ip][0];}
986 :
987 : /// value of global gradient of shape function i in integration point
988 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
989 : {UG_ASSERT(sh < num_sh(), "Invalid index");
990 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
991 :
992 : /// vector of global gradients in ip point
993 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
994 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvGlobalGrad[ip][0];}
995 :
996 : /// number of corners, that bound the scvf
997 0 : inline size_t num_corners() const {return numCo;}
998 :
999 : /// return local corner number i
1000 0 : inline const MathVector<dim>& local_corner(size_t co) const
1001 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
1002 :
1003 : /// return glbal corner number i
1004 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
1005 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
1006 :
1007 : private:
1008 : /// let outer class access private members
1009 : friend class DimFVGeometry<worldDim, dim>;
1010 :
1011 : // id of scv this bf belongs to
1012 : size_t nodeId;
1013 :
1014 : // ordering is:
1015 : // 1D: edgeMidPoint
1016 : // 2D: edgeMidPoint, CenterOfElement
1017 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
1018 : MathVector<dim> vLocPos[numCo]; // local corners of bf
1019 : MathVector<worldDim> vGloPos[numCo]; // global corners of bf
1020 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the bf
1021 :
1022 : // scvf part
1023 : size_t nip;
1024 : std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
1025 : std::vector<MathVector<worldDim> > vGlobalIP; // global integration point (size: nip)
1026 : const number* vWeight; // weight at ip
1027 : MathVector<worldDim> Normal; // normal (incl. area)
1028 : number Vol; // volume of bf
1029 :
1030 : // shapes and derivatives
1031 : size_t nsh;
1032 : std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
1033 : std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
1034 : std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
1035 : std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
1036 : std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
1037 : };
1038 :
1039 :
1040 : public:
1041 : /// construct object and initialize local values and sizes
1042 : DimFVGeometry();
1043 :
1044 : /// update local data
1045 : void update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad);
1046 0 : void update_local(ReferenceObjectID roid, const LFEID& lfeID)
1047 : {
1048 0 : update_local(roid, lfeID, lfeID.order() + 1);
1049 0 : }
1050 :
1051 : /// update data for given element
1052 0 : void update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords,
1053 : const ISubsetHandler* ish = NULL)
1054 : {
1055 0 : update(pElem, vCornerCoords, m_lfeID, m_quadOrderSCV, ish);
1056 0 : }
1057 :
1058 : /// update data for given element
1059 : void update(GridObject* pElem, const MathVector<worldDim>* vCornerCoords,
1060 : const LFEID& lfeID, size_t orderQuad,
1061 : const ISubsetHandler* ish = NULL);
1062 :
1063 : /// update boundary data for given element
1064 : void update_boundary_faces(GridObject* pElem,
1065 : const MathVector<worldDim>* vCornerCoords,
1066 : const ISubsetHandler* ish = NULL);
1067 :
1068 : /// number of SubControlVolumeFaces
1069 0 : inline size_t num_scvf() const {return m_numSCVF;};
1070 :
1071 : /// const access to SubControlVolumeFace number i
1072 0 : inline const SCVF& scvf(size_t i) const
1073 0 : {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
1074 :
1075 : /// number of SubControlVolumes
1076 0 : inline size_t num_scv() const {return m_numSCV;}
1077 :
1078 : /// const access to SubControlVolume number i
1079 0 : inline const SCV& scv(size_t i) const
1080 0 : {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
1081 :
1082 : /// number of shape functions
1083 0 : inline size_t num_sh() const {return m_nsh;};
1084 :
1085 : public:
1086 : /// returns number of all scvf ips
1087 0 : size_t num_scvf_ips() const {return m_numSCVFIP;}
1088 :
1089 : /// returns all ips of scvf as they appear in scv loop
1090 0 : const MathVector<dim>* scvf_local_ips() const {return &(m_vLocSCVF_IP[0]);}
1091 :
1092 : /// returns all ips of scvf as they appear in scv loop
1093 0 : const MathVector<worldDim>* scvf_global_ips() const {return &(m_vGlobSCVF_IP[0]);}
1094 :
1095 : /// returns number of all scv ips
1096 0 : size_t num_scv_ips() const {return m_numSCVIP;}
1097 :
1098 : /// returns all ips of scv as they appear in scv loop
1099 0 : const MathVector<dim>* scv_local_ips() const {return &(m_vLocSCV_IP[0]);}
1100 :
1101 : /// returns all ips of scv as they appear in scv loop
1102 0 : const MathVector<worldDim>* scv_global_ips() const {return &(m_vGlobSCV_IP[0]);}
1103 :
1104 :
1105 : protected:
1106 : // global and local ips on SCVF (size: numSCVFIP)
1107 : std::vector<MathVector<worldDim> > m_vGlobSCVF_IP;
1108 : std::vector<MathVector<dim> > m_vLocSCVF_IP;
1109 :
1110 : // global and local ips on SCVF
1111 : std::vector<MathVector<worldDim> > m_vGlobSCV_IP;
1112 : std::vector<MathVector<dim> > m_vLocSCV_IP;
1113 :
1114 : protected:
1115 : // miximum number of geom objects in a dimension
1116 : static const int maxMid = traits::maxNumSCVF +1;
1117 :
1118 : // subelement
1119 0 : struct SubElement
1120 : {
1121 : // shape number of corners of subelement
1122 : size_t vDoFID[traits::maxNumSCV];
1123 :
1124 : // local and global geom object midpoints for each dimension
1125 : // (most objects in 1 dim, i.e. number of edges, but +1 for 1D)
1126 : MathVector<dim> vvLocMid[dim+1][maxMid];
1127 : MathVector<worldDim> vvGloMid[dim+1][maxMid];
1128 :
1129 : // flag is subelement has boundary sides
1130 : bool isBndElem;
1131 :
1132 : // -1 is no bnd side, >= 0 corresponding side of whole element
1133 : std::vector<int> vElemBndSide;
1134 : };
1135 :
1136 : /// subelements (size: numSubElem)
1137 : std::vector<SubElement> m_vSubElem;
1138 :
1139 : public:
1140 : /// add subset that is interpreted as boundary subset.
1141 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
1142 :
1143 : /// removes subset that is interpreted as boundary subset.
1144 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
1145 :
1146 : /// reset all boundary subsets
1147 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
1148 :
1149 : /// number of registered boundary subsets
1150 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
1151 :
1152 : /// number of all boundary faces
1153 0 : inline size_t num_bf() const
1154 : {
1155 : typename std::map<int, std::vector<BF> >::const_iterator it;
1156 : size_t num = 0;
1157 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
1158 0 : num += (*it).second.size();
1159 0 : return num;
1160 : }
1161 :
1162 : /// number of boundary faces on subset 'subsetIndex'
1163 0 : inline size_t num_bf(int si) const
1164 : {
1165 : typename std::map<int, std::vector<BF> >::const_iterator it;
1166 : it = m_mapVectorBF.find(si);
1167 0 : if(it == m_mapVectorBF.end()) return 0;
1168 0 : else return (*it).second.size();
1169 : }
1170 :
1171 : /// returns the boundary face i for subsetIndex
1172 0 : inline const BF& bf(int si, size_t i) const
1173 : {
1174 : UG_ASSERT(i < num_bf(si), "Invalid index.");
1175 : typename std::map<int, std::vector<BF> >::const_iterator it;
1176 : it = m_mapVectorBF.find(si);
1177 0 : if(it == m_mapVectorBF.end()) UG_THROW("DimFVGeom: No BndSubset: "<<si);
1178 0 : return (*it).second[i];
1179 : }
1180 :
1181 : /// returns reference to vector of boundary faces for subsetIndex
1182 0 : inline const std::vector<BF>& bf(int si) const
1183 : {
1184 : typename std::map<int, std::vector<BF> >::const_iterator it;
1185 : it = m_mapVectorBF.find(si);
1186 0 : if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
1187 0 : return (*it).second;
1188 : }
1189 :
1190 0 : void reset_curr_elem() {m_pElem = NULL;}
1191 :
1192 : protected:
1193 : std::map<int, std::vector<BF> > m_mapVectorBF;
1194 : std::vector<BF> m_vEmptyVectorBF;
1195 :
1196 : private:
1197 : /// pointer to current element
1198 : GridObject* m_pElem;
1199 :
1200 : /// current reference object id
1201 : ReferenceObjectID m_roid;
1202 :
1203 : /// current order
1204 : int m_orderShape;
1205 :
1206 : /// current trial space
1207 : LFEID m_lfeID;
1208 :
1209 : /// current number of subelements
1210 : size_t m_numSubElem;
1211 :
1212 : /// number of shape functions
1213 : size_t m_nsh;
1214 :
1215 : /// current number of scvf
1216 : size_t m_numSCVF;
1217 :
1218 : /// current number of SCVF per SubElement
1219 : size_t m_numSCVFPerSubElem;
1220 :
1221 : /// SubControlVolumeFaces (size: numSCVF)
1222 : std::vector<SCVF> m_vSCVF;
1223 :
1224 : /// quadrature order
1225 : int m_quadOrderSCVF;
1226 :
1227 : /// number of scvf ip
1228 : size_t m_numSCVFIP;
1229 :
1230 :
1231 : /// current number of scv
1232 : size_t m_numSCV;
1233 :
1234 : /// number of SCV per SubElement
1235 : size_t m_numSCVPerSubElem;
1236 :
1237 : /// SubControlVolumes (size: numSCV)
1238 : std::vector<SCV> m_vSCV;
1239 :
1240 : /// current quadrature order scv
1241 : int m_quadOrderSCV;
1242 :
1243 : /// number of scv ip
1244 : size_t m_numSCVIP;
1245 : };
1246 :
1247 : }
1248 :
1249 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_VOLUME_HIGHER_ORDER_GEOMETRY__ */
|
