Line data Source code
1 : /*
2 : * Copyright (c) 2010-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__FV1_GEOMETRY__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FV1_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/local_finite_element/lagrange/lagrangep1.h"
52 : #include "lib_disc/quadrature/gauss/gauss_quad.h"
53 : #include "fv_util.h"
54 : #include "fv_geom_base.h"
55 :
56 : namespace ug{
57 :
58 : ////////////////////////////////////////////////////////////////////////////////
59 : // FV1 Geometry for Reference Element Type
60 : ////////////////////////////////////////////////////////////////////////////////
61 :
62 : // forward declaration
63 : template < typename TElem, int TWorldDim, bool TCondensed> class FV1Geometry_gen;
64 :
65 : /// Geometry and shape functions for 1st order Vertex-Centered Finite Volume
66 : /**
67 : * The class provides the geometry and shape functions for 1st order Vertex-Centered
68 : * Finite Element Finite Volume method based on the Donald Diagrams.
69 : *
70 : * Cf. class FV1Geometry_gen for the implementation.
71 : *
72 : * \tparam TElem Element type
73 : * \tparam TWorldDim (physical) world dimension
74 : */
75 : template <typename TElem, int TWorldDim>
76 0 : class FV1Geometry : public FV1Geometry_gen<TElem, TWorldDim, false> {};
77 :
78 : /// Geometry and shape functions for 1st order Vertex-Centered Finite Volume
79 : /**
80 : * The class provides the geometry and shape functions for 1st order Vertex-Centered
81 : * Finite Element Finite Volume method based on the Donald Diagrams.
82 : *
83 : * This class shifts the subcontrol volume face integration
84 : * points to the edges. This allows to reduce the matrix pattern and to avoid positive
85 : * off-diagonal entries in some cases. (For ex., the discretization of the Laplacian on
86 : * a grid of rectangles retains results in the 5-point stencil.) However note that, in many
87 : * cases, this leads to the discretization order reduction.
88 : *
89 : * Cf. class FV1Geometry_gen for the implementation.
90 : *
91 : * \tparam TElem Element type
92 : * \tparam TWorldDim (physical) world dimension
93 : */
94 : template <typename TElem, int TWorldDim>
95 0 : class FV1CondensedGeometry : public FV1Geometry_gen<TElem, TWorldDim, true> {};
96 :
97 : /// Geometry and shape functions for 1st order Vertex-Centered Finite Volume
98 : /**
99 : * The class provides the geometry and shape functions for 1st order Vertex-Centered
100 : * Finite Element Finite Volume method based on the Donald Diagrams.
101 : *
102 : * The class provides an option (TCondensed) to shift the subcontrol volume face integration
103 : * points to the edges. This allows to reduce the matrix pattern and to avoid positive
104 : * off-diagonal entries in some cases. (For ex., the discretization of the Laplacian on
105 : * a grid of rectangles retains results in the 5-point stencil.) However note that, in many
106 : * cases, this leads to the discretization order reduction.
107 : *
108 : * \tparam TElem Element type
109 : * \tparam TWorldDim (physical) world dimension
110 : * \tparam TCondensed if to shift the scvf ip's to midpoints of the edges
111 : */
112 : template <typename TElem, int TWorldDim, bool TCondensed>
113 0 : class FV1Geometry_gen : public FVGeometryBase
114 : {
115 : public:
116 : /// type of element
117 : typedef TElem elem_type;
118 :
119 : /// type of reference element
120 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
121 :
122 : /// used traits
123 : typedef fv1_traits<ref_elem_type, TWorldDim> traits;
124 :
125 : public:
126 : /// dimension of reference element
127 : static const int dim = ref_elem_type::dim;
128 :
129 : /// dimension of world
130 : static const int worldDim = TWorldDim;
131 :
132 : /// Hanging node flag: this Geometry does not support hanging nodes
133 : static const bool usesHangingNodes = false;
134 :
135 : /// flag indicating if local data may change
136 : static const bool staticLocalData = true;
137 :
138 : /// whether the scheme shifts the scvf ip's to midpoints of the edges
139 : static const bool condensed_scvf_ips = TCondensed;
140 :
141 : public:
142 : /// order
143 : static const int order = 1;
144 :
145 : /// number of SubControlVolumes
146 : static const size_t numSCV = traits::numSCV;
147 :
148 : /// type of SubControlVolume
149 : typedef typename traits::scv_type scv_type;
150 :
151 : /// number of SubControlVolumeFaces
152 : static const size_t numSCVF = traits::numSCVF;
153 :
154 : /// type of Shape function used
155 : typedef LagrangeP1<ref_elem_type> local_shape_fct_set_type;
156 :
157 : /// number of shape functions
158 : static const size_t nsh = local_shape_fct_set_type::nsh;
159 :
160 : /// number of integration points
161 : static const size_t nip = 1;
162 :
163 : public:
164 : /// Sub-Control Volume Face structure
165 : /**
166 : * Each finite element is cut by several sub-control volume faces. The idea
167 : * is that the "dual" skeleton formed by the sub control volume faces of
168 : * all elements again gives rise to a regular mesh with closed
169 : * (lipschitz-bounded) control volumes. The SCVF are the boundary of the
170 : * control volume. In computation the flux over each SCVF must be the same
171 : * in both directions over the face in order to guarantee the conservation
172 : * property.
173 : */
174 : class SCVF
175 : {
176 : public:
177 : /// Number of corners of scvf
178 : static const size_t numCo = traits::NumCornersOfSCVF;
179 :
180 : public:
181 0 : SCVF() {}
182 :
183 : /// index of SubControlVolume on one side of the scvf
184 0 : inline size_t from() const {return From;}
185 :
186 : /// index of SubControlVolume on one side of the scvf
187 0 : inline size_t to() const {return To;}
188 :
189 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
190 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
191 :
192 : /// number of integration points on scvf
193 0 : inline size_t num_ip() const {return nip;}
194 :
195 : /// local integration point of scvf
196 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
197 :
198 : /// global integration point of scvf
199 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
200 :
201 : /// Transposed Inverse of Jacobian in integration point
202 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
203 :
204 : /// Determinant of Jacobian in integration point
205 0 : inline number detJ() const {return detj;}
206 :
207 : /// number of shape functions
208 0 : inline size_t num_sh() const {return nsh;}
209 :
210 : /// value of shape function i in integration point
211 0 : inline number shape(size_t sh) const {return vShape[sh];}
212 :
213 : /// vector of shape functions in ip point
214 0 : inline const number* shape_vector() const {return vShape;}
215 :
216 : /// value of local gradient of shape function i in integration point
217 0 : inline const MathVector<dim>& local_grad(size_t sh) const
218 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
219 :
220 : /// vector of local gradients in ip point
221 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
222 :
223 : /// value of global gradient of shape function i in integration point
224 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
225 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
226 :
227 : /// vector of global gradients in ip point
228 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
229 :
230 : /// number of corners, that bound the scvf
231 0 : inline size_t num_corners() const {return numCo;}
232 :
233 : /// return local corner number i
234 0 : inline const MathVector<dim>& local_corner(size_t co) const
235 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
236 :
237 : /// return global corner number i
238 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
239 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
240 :
241 : private:
242 : // let outer class access private members
243 : friend class FV1Geometry_gen<TElem, TWorldDim, TCondensed>;
244 :
245 : // This scvf separates the scv with the ids given in "from" and "to"
246 : // The computed normal points in direction from->to
247 : size_t From, To;
248 :
249 : // The normal on the SCVF pointing (from -> to)
250 : MathVector<worldDim> Normal; // normal (incl. area)
251 :
252 : // ordering is:
253 : // 1D: edgeMidPoint
254 : // 2D: edgeMidPoint, CenterOfElement
255 : // 3D: edgeMidPoint, Side #1, CenterOfElement, Side #2
256 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
257 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
258 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scvf
259 :
260 : // scvf part
261 : MathVector<dim> localIP; // local integration point
262 : MathVector<worldDim> globalIP; // global integration point
263 :
264 : // shapes and derivatives
265 : number vShape[nsh]; // shapes at ip
266 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
267 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
268 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
269 : number detj; // Jacobian det at ip
270 : };
271 :
272 : /// sub control volume structure
273 : class SCV
274 : {
275 : public:
276 : /// Number of corners of scvf
277 : static const size_t numCo = traits::NumCornersOfSCV;
278 :
279 : public:
280 0 : SCV() {};
281 :
282 : /// volume of scv
283 0 : inline number volume() const {return Vol;}
284 :
285 : /// number of corners, that bound the scvf
286 0 : inline size_t num_corners() const {return numCo;}
287 :
288 : /// return local corner number i
289 0 : inline const MathVector<dim>& local_corner(size_t co) const
290 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
291 :
292 : /// return global corner number i
293 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
294 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
295 :
296 : /// return local corners
297 0 : inline const MathVector<dim>* local_corners() const
298 0 : {return &vLocPos[0];}
299 :
300 : /// return global corners
301 0 : inline const MathVector<worldDim>* global_corners() const
302 0 : {return &vGloPos[0];}
303 :
304 : /// node id that this scv is associated to
305 0 : inline size_t node_id() const {return nodeId;}
306 :
307 : /// number of integration points
308 0 : inline size_t num_ip() const {return nip;}
309 :
310 : /// local integration point of scv
311 0 : inline const MathVector<dim>& local_ip() const {return vLocPos[0];}
312 :
313 : /// global integration point
314 0 : inline const MathVector<worldDim>& global_ip() const {return vGloPos[0];}
315 :
316 : /// Transposed Inverse of Jacobian in integration point
317 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
318 :
319 : /// Determinant of Jacobian in integration point
320 0 : inline number detJ() const {return detj;}
321 :
322 : /// number of shape functions
323 0 : inline size_t num_sh() const {return nsh;}
324 :
325 : /// value of shape function i in integration point
326 0 : inline number shape(size_t sh) const {return vShape[sh];}
327 :
328 : /// vector of shape functions in ip point
329 0 : inline const number* shape_vector() const {return vShape;}
330 :
331 : /// value of local gradient of shape function i in integration point
332 0 : inline const MathVector<dim>& local_grad(size_t sh) const
333 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
334 :
335 : /// vector of local gradients in ip point
336 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
337 :
338 : /// value of global gradient of shape function i in integration point
339 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
340 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
341 :
342 : /// vector of global gradients in ip point
343 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
344 :
345 : private:
346 : // let outer class access private members
347 : friend class FV1Geometry_gen<TElem, TWorldDim, TCondensed>;
348 :
349 : // node id of associated node
350 : size_t nodeId;
351 :
352 : // volume of scv
353 : number Vol;
354 :
355 : // local and global positions of this element
356 : MathVector<dim> vLocPos[numCo]; // local position of node
357 : MathVector<worldDim> vGloPos[numCo]; // global position of node
358 : MidID midId[numCo]; // dimension and id of object, that's midpoint bounds the scv
359 :
360 : // shapes and derivatives
361 : number vShape[nsh]; // shapes at ip
362 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
363 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
364 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
365 : number detj; // Jacobian det at ip
366 : };
367 :
368 : /// boundary face
369 : class BF
370 : {
371 : public:
372 : /// Number of corners of bf
373 : static const size_t numCo = traits::NumCornersOfBF;
374 :
375 : public:
376 0 : BF() {}
377 :
378 : /// index of SubControlVolume of the bf
379 0 : inline size_t node_id() const {return nodeId;}
380 :
381 : /// number of integration points on bf
382 0 : inline size_t num_ip() const {return nip;}
383 :
384 : /// local integration point of bf
385 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
386 :
387 : /// global integration point of bf
388 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
389 :
390 : /// outer normal on bf. Norm is equal to area
391 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
392 :
393 : /// volume of bf
394 0 : inline number volume() const {return Vol;}
395 :
396 : /// Transposed Inverse of Jacobian in integration point
397 0 : inline const MathMatrix<worldDim, dim>& JTInv() const {return JtInv;}
398 :
399 : /// Determinant of Jacobian in integration point
400 0 : inline number detJ() const {return detj;}
401 :
402 : /// number of shape functions
403 0 : inline size_t num_sh() const {return nsh;}
404 :
405 : /// value of shape function i in integration point
406 0 : inline number shape(size_t sh) const
407 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vShape[sh];}
408 :
409 : /// vector of local gradients in ip point
410 0 : inline const number* shape_vector() const {return vShape;}
411 :
412 : /// value of local gradient of shape function i in integration point
413 0 : inline const MathVector<dim>& local_grad(size_t sh) const
414 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
415 :
416 : /// vector of local gradients in ip point
417 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
418 :
419 : /// value of global gradient of shape function i in integration point
420 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
421 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
422 :
423 : /// vector of global gradients in ip point
424 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
425 :
426 : /// number of corners, that bound the scvf
427 0 : inline size_t num_corners() const {return numCo;}
428 :
429 : /// return local corner number i
430 0 : inline const MathVector<dim>& local_corner(size_t co) const
431 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
432 :
433 : /// return global corner number i
434 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
435 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
436 :
437 : private:
438 : /// let outer class access private members
439 : friend class FV1Geometry_gen<TElem, TWorldDim, TCondensed>;
440 :
441 : // id of scv this bf belongs to
442 : size_t nodeId;
443 :
444 : // ordering is:
445 : // 1D: edgeMidPoint
446 : // 2D: edgeMidPoint, CenterOfElement
447 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
448 : MathVector<dim> vLocPos[numCo]; // local corners of bf
449 : MathVector<worldDim> vGloPos[numCo]; // global corners of bf
450 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the bf
451 :
452 : // scvf part
453 : MathVector<dim> localIP; // local integration point
454 : MathVector<worldDim> globalIP; // global integration point
455 : MathVector<worldDim> Normal; // normal (incl. area)
456 : number Vol; // volume of bf
457 :
458 : // shapes and derivatives
459 : number vShape[nsh]; // shapes at ip
460 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
461 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
462 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
463 : number detj; // Jacobian det at ip
464 : };
465 :
466 : public:
467 : /// construct object and initialize local values and sizes
468 : FV1Geometry_gen();
469 :
470 : /// update local data
471 : void update_local_data();
472 :
473 : /// update data for given element
474 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
475 : const ISubsetHandler* ish = NULL);
476 :
477 : /// update boundary data for given element
478 : void update_boundary_faces(GridObject* elem,
479 : const MathVector<worldDim>* vCornerCoords,
480 : const ISubsetHandler* ish = NULL);
481 :
482 : /// get the element
483 0 : TElem* elem() const {return m_pElem;}
484 :
485 : /// get vector of the global coordinates of corners for current element
486 0 : const MathVector<worldDim>* corners() const {return m_vvGloMid[0];}
487 :
488 : /// number of SubControlVolumeFaces
489 0 : size_t num_scvf() const {return numSCVF;};
490 :
491 : /// const access to SubControlVolumeFace number i
492 0 : const SCVF& scvf(size_t i) const
493 0 : {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
494 :
495 : /// number of SubControlVolumes
496 : // do not use this method to obtain the number of shape functions,
497 : // since this is NOT the same for pyramids; use num_sh() instead.
498 0 : size_t num_scv() const {return numSCV;}
499 :
500 : /// const access to SubControlVolume number i
501 0 : const SCV& scv(size_t i) const
502 0 : {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
503 :
504 : /// number of shape functions
505 0 : size_t num_sh() const {return nsh;};
506 :
507 : /// returns reference object id
508 0 : ReferenceObjectID roid() const {return ref_elem_type::REFERENCE_OBJECT_ID;}
509 :
510 :
511 : public:
512 : /// returns number of all scvf ips
513 0 : size_t num_scvf_ips() const {return numSCVF;}
514 :
515 : /// returns all ips of scvf as they appear in scv loop
516 0 : const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
517 :
518 : /// returns all ips of scvf as they appear in scv loop
519 0 : const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
520 :
521 : /// returns number of all scv ips
522 0 : size_t num_scv_ips() const {return numSCV;}
523 :
524 : /// returns all ips of scv as they appear in scv loop
525 0 : const MathVector<dim>* scv_local_ips() const {
526 : if(ref_elem_type::REFERENCE_OBJECT_ID == ROID_PYRAMID || ref_elem_type::REFERENCE_OBJECT_ID == ROID_OCTAHEDRON)
527 0 : return &(m_vLocSCV_IP[0]);
528 : else
529 0 : return &(m_vvLocMid[0][0]);
530 : }
531 :
532 : /// returns all ips of scv as they appear in scv loop
533 0 : const MathVector<worldDim>* scv_global_ips() const {
534 : if(ref_elem_type::REFERENCE_OBJECT_ID == ROID_PYRAMID || ref_elem_type::REFERENCE_OBJECT_ID == ROID_OCTAHEDRON)
535 0 : return &(m_vGlobSCV_IP[0]);
536 : else
537 0 : return &(m_vvGloMid[0][0]);
538 : }
539 :
540 : /// return local coords for node ID
541 0 : const MathVector<dim>& local_node_position(size_t nodeID) const
542 : {
543 : UG_ASSERT(nodeID < (size_t) ref_elem_type::numCorners, "Invalid node id.");
544 0 : return m_vvLocMid[0][nodeID];
545 : }
546 :
547 : /// return global coords for node ID
548 0 : const MathVector<worldDim>& global_node_position(size_t nodeID) const
549 : {
550 : UG_ASSERT(nodeID < (size_t) ref_elem_type::numCorners, "Invalid node id.");
551 0 : return m_vvGloMid[0][nodeID];
552 : }
553 :
554 : /// returns the local coordinates of the center of mass of the element
555 0 : const MathVector<dim>* coe_local() const {return &(m_vvLocMid[dim][0]);}
556 :
557 : /// returns the global coordinates of the center of mass of the element
558 0 : const MathVector<worldDim>* coe_global() const {return &(m_vvGloMid[dim][0]);}
559 :
560 : protected:
561 : // global and local ips on SCVF
562 : MathVector<worldDim> m_vGlobSCVF_IP[numSCVF];
563 : MathVector<dim> m_vLocSCVF_IP[numSCVF];
564 :
565 : // global and local ips on SCV (only needed for Pyramid and Octahedron)
566 : MathVector<worldDim> m_vGlobSCV_IP[numSCV];
567 : MathVector<dim> m_vLocSCV_IP[numSCV];
568 :
569 : public:
570 : /// add subset that is interpreted as boundary subset.
571 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
572 :
573 : /// removes subset that is interpreted as boundary subset.
574 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
575 :
576 : /// reset all boundary subsets
577 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
578 :
579 : /// number of registered boundary subsets
580 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
581 :
582 : /// number of all boundary faces
583 0 : inline size_t num_bf() const
584 : {
585 : typename std::map<int, std::vector<BF> >::const_iterator it;
586 : size_t num = 0;
587 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
588 0 : num += (*it).second.size();
589 0 : return num;
590 : }
591 :
592 : /// number of boundary faces on subset 'subsetIndex'
593 0 : inline size_t num_bf(int si) const
594 : {
595 : typename std::map<int, std::vector<BF> >::const_iterator it;
596 : it = m_mapVectorBF.find(si);
597 0 : if(it == m_mapVectorBF.end()) return 0;
598 0 : else return (*it).second.size();
599 : }
600 :
601 : /// returns the boundary face i for subsetIndex
602 0 : inline const BF& bf(int si, size_t i) const
603 : {
604 : UG_ASSERT(i < num_bf(si), "Invalid index.");
605 : typename std::map<int, std::vector<BF> >::const_iterator it;
606 : it = m_mapVectorBF.find(si);
607 0 : if(it == m_mapVectorBF.end()) UG_THROW("FVGeom: No bnd face for subset"<<si);
608 0 : return (*it).second[i];
609 : }
610 :
611 : /// returns reference to vector of boundary faces for subsetIndex
612 0 : inline const std::vector<BF>& bf(int si) const
613 : {
614 : typename std::map<int, std::vector<BF> >::const_iterator it;
615 : it = m_mapVectorBF.find(si);
616 0 : if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
617 0 : return (*it).second;
618 : }
619 :
620 0 : void reset_curr_elem() {m_pElem = NULL;}
621 :
622 : protected:
623 : std::map<int, std::vector<BF> > m_mapVectorBF;
624 : std::vector<BF> m_vEmptyVectorBF;
625 :
626 : private:
627 : /// pointer to current element
628 : TElem* m_pElem;
629 :
630 : /// max number of geom objects in all dimensions
631 : // (most objects in 1 dim, i.e. number of edges, but +1 for 1D)
632 : static const int maxMid = numSCVF + 1;
633 :
634 : /// local and global geom object midpoints for each dimension
635 : MathVector<dim> m_vvLocMid[dim+1][maxMid];
636 : MathVector<worldDim> m_vvGloMid[dim+1][maxMid];
637 :
638 : /// SubControlVolumeFaces
639 : SCVF m_vSCVF[numSCVF];
640 :
641 : /// SubControlVolumes
642 : SCV m_vSCV[numSCV];
643 :
644 : /// Reference Mapping
645 : ReferenceMapping<ref_elem_type, worldDim> m_mapping;
646 :
647 : /// Reference Element
648 : const ref_elem_type& m_rRefElem;
649 :
650 : /// Shape function set
651 : const local_shape_fct_set_type& m_rTrialSpace;
652 : };
653 :
654 : ////////////////////////////////////////////////////////////////////////////////
655 : // Dim-dependent Finite Volume Geometry
656 : ////////////////////////////////////////////////////////////////////////////////
657 :
658 : /// Geometry and shape functions for 1st order Vertex-Centered Finite Volume
659 : /**
660 : * \tparam TDim reference element dim
661 : * \tparam TWorldDim (physical) world dimension
662 : */
663 : template <int TDim, int TWorldDim = TDim>
664 0 : class DimFV1Geometry : public FVGeometryBase
665 : {
666 : public:
667 : /// used traits
668 : typedef fv1_dim_traits<TDim, TWorldDim> traits;
669 :
670 : public:
671 : /// dimension of reference element
672 : static const int dim = TDim;
673 :
674 : /// dimension of world
675 : static const int worldDim = TWorldDim;
676 :
677 : /// Hanging node flag: this Geometry does not support hanging nodes
678 : static const bool usesHangingNodes = false;
679 :
680 : /// flag indicating if local data may change
681 : static const bool staticLocalData = false;
682 :
683 : public:
684 : /// order
685 : static const int order = 1;
686 :
687 : /// number of SubControlVolumes
688 : static const size_t maxNumSCV = traits::maxNumSCV;
689 :
690 : /// type of SubControlVolume
691 : typedef typename traits::scv_type scv_type;
692 :
693 : /// max number of SubControlVolumeFaces
694 : static const size_t maxNumSCVF = traits::maxNumSCVF;
695 :
696 : /// max number of shape functions
697 : static const size_t maxNSH = traits::maxNSH;
698 :
699 : /// number of integration points
700 : static const size_t nip = 1;
701 :
702 : public:
703 : /// Sub-Control Volume Face structure
704 : /**
705 : * Each finite element is cut by several sub-control volume faces. The idea
706 : * is that the "dual" skeleton formed by the sub control volume faces of
707 : * all elements again gives rise to a regular mesh with closed
708 : * (lipschitz-bounded) control volumes. The SCVF are the boundary of the
709 : * control volume. In computation the flux over each SCVF must be the same
710 : * in both directions over the face in order to guarantee the conservation
711 : * property.
712 : */
713 : class SCVF
714 : {
715 : public:
716 : /// Number of corners of scvf
717 : static const size_t numCo = traits::NumCornersOfSCVF;
718 :
719 : public:
720 0 : SCVF() {}
721 :
722 : /// index of SubControlVolume on one side of the scvf
723 0 : inline size_t from() const {return From;}
724 :
725 : /// index of SubControlVolume on one side of the scvf
726 0 : inline size_t to() const {return To;}
727 :
728 : /// number of integration points on scvf
729 0 : inline size_t num_ip() const {return nip;}
730 :
731 : /// local integration point of scvf
732 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
733 :
734 : /// global integration point of scvf
735 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
736 :
737 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
738 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
739 :
740 : /// Transposed Inverse of Jacobian in integration point
741 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
742 :
743 : /// Determinant of Jacobian in integration point
744 0 : inline number detJ() const {return detj;}
745 :
746 : /// number of shape functions
747 0 : inline size_t num_sh() const {return numSH;}
748 :
749 : /// value of shape function i in integration point
750 0 : inline number shape(size_t sh) const {return vShape[sh];}
751 :
752 : /// vector of shape functions in ip point
753 0 : inline const number* shape_vector() const {return vShape;}
754 :
755 : /// value of local gradient of shape function i in integration point
756 0 : inline const MathVector<dim>& local_grad(size_t sh) const
757 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
758 :
759 : /// vector of local gradients in ip point
760 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
761 :
762 : /// value of global gradient of shape function i in integration point
763 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
764 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
765 :
766 : /// vector of gloabl gradients in ip point
767 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
768 :
769 : /// number of corners, that bound the scvf
770 0 : inline size_t num_corners() const {return numCo;}
771 :
772 : /// return local corner number i
773 0 : inline const MathVector<dim>& local_corner(size_t co) const
774 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
775 :
776 : /// return glbal corner number i
777 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
778 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
779 :
780 : private:
781 : // let outer class access private members
782 : friend class DimFV1Geometry<dim, worldDim>;
783 :
784 : // This scvf separates the scv with the ids given in "from" and "to"
785 : // The computed normal points in direction from->to
786 : size_t From, To;
787 :
788 : // The normal on the SCVF pointing (from -> to)
789 : MathVector<worldDim> Normal; // normal (incl. area)
790 :
791 : // ordering is:
792 : // 1D: edgeMidPoint
793 : // 2D: edgeMidPoint, CenterOfElement
794 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
795 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
796 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
797 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scvf
798 :
799 : // scvf part
800 : MathVector<dim> localIP; // local integration point
801 : MathVector<worldDim> globalIP; // global integration point
802 :
803 : // shapes and derivatives
804 : size_t numSH;
805 : number vShape[maxNSH]; // shapes at ip
806 : MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
807 : MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
808 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
809 : number detj; // Jacobian det at ip
810 : };
811 :
812 : /// sub control volume structure
813 : class SCV
814 : {
815 : public:
816 : /// Number of corners of scv
817 : static const size_t numCo = traits::NumCornersOfSCV;
818 :
819 : public:
820 0 : SCV() {};
821 :
822 : /// volume of scv
823 0 : inline number volume() const {return Vol;}
824 :
825 : /// number of corners, that bound the scv
826 0 : inline size_t num_corners() const {return numCo;}
827 :
828 : /// return local corner number i
829 0 : inline const MathVector<dim>& local_corner(size_t co) const
830 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
831 :
832 : /// return global corner number i
833 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
834 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
835 :
836 : /// node id that this scv is associated to
837 0 : inline size_t node_id() const {return nodeId;}
838 :
839 : /// number of integration points
840 0 : inline size_t num_ip() const {return nip;}
841 :
842 : /// local integration point of scv
843 0 : inline const MathVector<dim>& local_ip() const {return vLocPos[0];}
844 :
845 : /// global integration point
846 0 : inline const MathVector<worldDim>& global_ip() const {return vGloPos[0];}
847 :
848 : /// Transposed Inverse of Jacobian in integration point
849 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
850 :
851 : /// Determinant of Jacobian in integration point
852 0 : inline number detJ() const {return detj;}
853 :
854 : /// number of shape functions
855 0 : inline size_t num_sh() const {return numSH;}
856 :
857 : /// value of shape function i in integration point
858 0 : inline number shape(size_t sh) const {return vShape[sh];}
859 :
860 : /// vector of shape functions in ip point
861 0 : inline const number* shape_vector() const {return vShape;}
862 :
863 : /// value of local gradient of shape function i in integration point
864 0 : inline const MathVector<dim>& local_grad(size_t sh) const
865 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
866 :
867 : /// vector of local gradients in ip point
868 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
869 :
870 : /// value of global gradient of shape function i in integration point
871 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
872 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
873 :
874 : /// vector of global gradients in ip point
875 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
876 :
877 : private:
878 : // let outer class access private members
879 : friend class DimFV1Geometry<dim, worldDim>;
880 :
881 : // node id of associated node
882 : size_t nodeId;
883 :
884 : // volume of scv
885 : number Vol;
886 :
887 : // local and global positions of this element
888 : MathVector<dim> vLocPos[numCo]; // local position of node
889 : MathVector<worldDim> vGloPos[numCo]; // global position of node
890 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the scv
891 :
892 : // shapes and derivatives
893 : size_t numSH;
894 : number vShape[maxNSH]; // shapes at ip
895 : MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
896 : MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
897 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
898 : number detj; // Jacobian det at ip
899 : };
900 :
901 : /// boundary face
902 : class BF
903 : {
904 : public:
905 : /// Number of corners of bf
906 : static const size_t numCo = traits::NumCornersOfSCVF;
907 :
908 : public:
909 0 : BF() {}
910 :
911 : /// index of SubControlVolume of the bf
912 0 : inline size_t node_id() const {return nodeId;}
913 :
914 : /// number of integration points on bf
915 0 : inline size_t num_ip() const {return nip;}
916 :
917 : /// local integration point of bf
918 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
919 :
920 : /// global integration point of bf
921 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
922 :
923 : /// outer normal on bf. Norm is equal to area
924 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
925 :
926 : /// volume of bf
927 0 : inline number volume() const {return Vol;}
928 :
929 : /// Transposed Inverse of Jacobian in integration point
930 0 : inline const MathMatrix<worldDim, dim>& JTInv() const {return JtInv;}
931 :
932 : /// Determinant of Jacobian in integration point
933 0 : inline number detJ() const {return detj;}
934 :
935 : /// number of shape functions
936 0 : inline size_t num_sh() const {return numSH;}
937 :
938 : /// value of shape function i in integration point
939 0 : inline number shape(size_t sh) const
940 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vShape[sh];}
941 :
942 : /// vector of local gradients in ip point
943 0 : inline const number* shape_vector() const {return vShape;}
944 :
945 : /// value of local gradient of shape function i in integration point
946 0 : inline const MathVector<dim>& local_grad(size_t sh) const
947 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
948 :
949 : /// vector of local gradients in ip point
950 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
951 :
952 : /// value of global gradient of shape function i in integration point
953 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
954 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
955 :
956 : /// vector of global gradients in ip point
957 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
958 :
959 : /// number of corners, that bound the scvf
960 0 : inline size_t num_corners() const {return numCo;}
961 :
962 : /// return local corner number i
963 0 : inline const MathVector<dim>& local_corner(size_t co) const
964 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
965 :
966 : /// return global corner number i
967 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
968 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
969 :
970 : private:
971 : /// let outer class access private members
972 : friend class DimFV1Geometry<dim, worldDim>;
973 :
974 : // id of scv this bf belongs to
975 : size_t nodeId;
976 :
977 : // ordering is:
978 : // 1D: edgeMidPoint
979 : // 2D: edgeMidPoint, CenterOfElement
980 : // 3D: edgeMidPoint, Side one, CenterOfElement, Side two
981 : MathVector<dim> vLocPos[numCo]; // local corners of bf
982 : MathVector<worldDim> vGloPos[numCo]; // global corners of bf
983 : MidID vMidID[numCo]; // dimension and id of object, that's midpoint bounds the bf
984 :
985 : // scvf part
986 : MathVector<dim> localIP; // local integration point
987 : MathVector<worldDim> globalIP; // global integration point
988 : MathVector<worldDim> Normal; // normal (incl. area)
989 : number Vol; // volume of bf
990 :
991 : // shapes and derivatives
992 : size_t numSH;
993 : number vShape[maxNSH]; // shapes at ip
994 : MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
995 : MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
996 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
997 : number detj; // Jacobian det at ip
998 : };
999 :
1000 : public:
1001 : /// construct object and initialize local values and sizes
1002 0 : DimFV1Geometry() : m_pElem(NULL), m_roid(ROID_UNKNOWN) {};
1003 :
1004 : /// update data for given element
1005 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
1006 : const ISubsetHandler* ish = NULL);
1007 :
1008 : /// update boundary data for given element
1009 : void update_boundary_faces(GridObject* elem,
1010 : const MathVector<worldDim>* vCornerCoords,
1011 : const ISubsetHandler* ish = NULL);
1012 :
1013 : /// get the element
1014 0 : GridObject* elem() const {return m_pElem;}
1015 :
1016 : /// get vector of corners for current element
1017 0 : const MathVector<worldDim>* corners() const {return m_vvGloMid[0];}
1018 :
1019 : /// number of SubControlVolumeFaces
1020 0 : size_t num_scvf() const {return m_numSCVF;};
1021 :
1022 : /// const access to SubControlVolumeFace number i
1023 0 : const SCVF& scvf(size_t i) const
1024 0 : {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
1025 :
1026 : /// number of SubControlVolumes
1027 : // do not use this method to obtain the number of shape functions,
1028 : // since this is NOT the same for pyramids; use num_sh() instead.
1029 0 : size_t num_scv() const {return m_numSCV;}
1030 :
1031 : /// const access to SubControlVolume number i
1032 0 : const SCV& scv(size_t i) const
1033 0 : {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
1034 :
1035 : /// number of shape functions
1036 0 : size_t num_sh() const {return m_nsh;};
1037 :
1038 : public:
1039 : /// returns number of all scvf ips
1040 0 : size_t num_scvf_ips() const {return m_numSCVF;}
1041 :
1042 : /// returns all ips of scvf as they appear in scv loop
1043 0 : const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
1044 :
1045 : /// returns all ips of scvf as they appear in scv loop
1046 0 : const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
1047 :
1048 : /// returns number of all scv ips
1049 0 : size_t num_scv_ips() const {return m_numSCV;}
1050 :
1051 : /// returns all ips of scv as they appear in scv loop
1052 0 : const MathVector<dim>* scv_local_ips() const {
1053 0 : if(m_roid == ROID_PYRAMID || m_roid == ROID_OCTAHEDRON)
1054 0 : return &(m_vLocSCV_IP[0]);
1055 : else
1056 0 : return &(m_vvLocMid[0][0]);
1057 : }
1058 :
1059 : /// returns all ips of scv as they appear in scv loop
1060 0 : const MathVector<worldDim>* scv_global_ips() const {
1061 0 : if(m_roid == ROID_PYRAMID || m_roid == ROID_OCTAHEDRON)
1062 0 : return &(m_vGlobSCV_IP[0]);
1063 : else
1064 0 : return &(m_vvGloMid[0][0]);
1065 : }
1066 :
1067 : /// return local coords for node ID
1068 0 : const MathVector<dim>& local_node_position(size_t nodeID) const
1069 : {
1070 : UG_ASSERT(nodeID < (size_t) maxMid, "Invalid node id.");
1071 0 : return m_vvLocMid[0][nodeID];
1072 : }
1073 :
1074 : /// return global coords for node ID
1075 0 : const MathVector<worldDim>& global_node_position(size_t nodeID) const
1076 : {
1077 : UG_ASSERT(nodeID < (size_t) maxMid, "Invalid node id.");
1078 0 : return m_vvGloMid[0][nodeID];
1079 : }
1080 :
1081 : /// returns the local coordinates of the center of mass of the element
1082 0 : const MathVector<dim>* coe_local() const {return &(m_vvLocMid[dim][0]);}
1083 :
1084 : /// returns the global coordinates of the center of mass of the element
1085 0 : const MathVector<worldDim>* coe_global() const {return &(m_vvGloMid[dim][0]);}
1086 :
1087 : /// returns reference object id
1088 0 : ReferenceObjectID roid() const {return m_roid;}
1089 :
1090 : /// update local data
1091 : void update_local(ReferenceObjectID roid);
1092 :
1093 : protected:
1094 : // global and local ips on SCVF
1095 : MathVector<worldDim> m_vGlobSCVF_IP[maxNumSCVF];
1096 : MathVector<dim> m_vLocSCVF_IP[maxNumSCVF];
1097 :
1098 : // global and local ips on SCV (only needed for Pyramid and Octahedron)
1099 : MathVector<worldDim> m_vGlobSCV_IP[maxNumSCV];
1100 : MathVector<dim> m_vLocSCV_IP[maxNumSCV];
1101 :
1102 : public:
1103 : /// add subset that is interpreted as boundary subset.
1104 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
1105 :
1106 : /// removes subset that is interpreted as boundary subset.
1107 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
1108 :
1109 : /// reset all boundary subsets
1110 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
1111 :
1112 : /// number of registered boundary subsets
1113 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
1114 :
1115 : /// number of all boundary faces
1116 0 : inline size_t num_bf() const
1117 : {
1118 : typename std::map<int, std::vector<BF> >::const_iterator it;
1119 : size_t num = 0;
1120 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
1121 0 : num += (*it).second.size();
1122 0 : return num;
1123 : }
1124 :
1125 : /// number of boundary faces on subset 'subsetIndex'
1126 0 : inline size_t num_bf(int si) const
1127 : {
1128 : typename std::map<int, std::vector<BF> >::const_iterator it;
1129 : it = m_mapVectorBF.find(si);
1130 0 : if(it == m_mapVectorBF.end()) return 0;
1131 0 : else return (*it).second.size();
1132 : }
1133 :
1134 : /// returns the boundary face i for subsetIndex
1135 0 : inline const BF& bf(int si, size_t i) const
1136 : {
1137 : typename std::map<int, std::vector<BF> >::const_iterator it;
1138 : it = m_mapVectorBF.find(si);
1139 0 : if(it == m_mapVectorBF.end()) UG_THROW("DimFV1Geom: No BndSubset "<<si);
1140 0 : return (*it).second[i];
1141 : }
1142 :
1143 : /// returns reference to vector of boundary faces for subsetIndex
1144 0 : inline const std::vector<BF>& bf(int si) const
1145 : {
1146 : typename std::map<int, std::vector<BF> >::const_iterator it;
1147 : it = m_mapVectorBF.find(si);
1148 0 : if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
1149 0 : return (*it).second;
1150 : }
1151 :
1152 : protected:
1153 : std::map<int, std::vector<BF> > m_mapVectorBF;
1154 : std::vector<BF> m_vEmptyVectorBF;
1155 :
1156 : private:
1157 : /// pointer to current element
1158 : GridObject* m_pElem;
1159 :
1160 : /// current reference object id
1161 : ReferenceObjectID m_roid;
1162 :
1163 : /// current number of scv
1164 : size_t m_numSCV;
1165 :
1166 : /// current number of scvf
1167 : size_t m_numSCVF;
1168 :
1169 : /// current number of shape functions
1170 : size_t m_nsh;
1171 :
1172 : /// max number of geometric objects in a dimension
1173 : // (most objects in 1 dim, i.e. number of edges, but +1 for 1D)
1174 : static const int maxMid = maxNumSCVF + 1;
1175 :
1176 : /// local and global geom object midpoints for each dimension
1177 : MathVector<dim> m_vvLocMid[dim+1][maxMid];
1178 : MathVector<worldDim> m_vvGloMid[dim+1][maxMid];
1179 :
1180 : /// SubControlVolumeFaces
1181 : SCVF m_vSCVF[maxNumSCVF];
1182 :
1183 : /// SubControlVolumes
1184 : SCV m_vSCV[maxNumSCV];
1185 : };
1186 :
1187 :
1188 : ////////////////////////////////////////////////////////////////////////////////
1189 : // FV1 Manifold Geometry
1190 : ////////////////////////////////////////////////////////////////////////////////
1191 :
1192 : template <typename TElem, int TWorldDim>
1193 : class FV1ManifoldGeometry
1194 : {
1195 : public:
1196 : // type of element
1197 : typedef TElem elem_type;
1198 :
1199 : // type of reference element
1200 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
1201 :
1202 : public:
1203 : // order
1204 : static const int order = 1;
1205 :
1206 : // number of BoundaryFaces
1207 : static const size_t m_numBF = ref_elem_type::numCorners;
1208 :
1209 : // dimension of reference element
1210 : static const int dim = ref_elem_type::dim;
1211 :
1212 : // dimension of world
1213 : static const int worldDim = TWorldDim;
1214 :
1215 : // type of BoundaryFaces
1216 : typedef typename fv1_traits<ref_elem_type, dim>::scv_type bf_type;
1217 :
1218 : // Hanging node flag: this Geometry does not support hanging nodes
1219 : static const bool usesHangingNodes = false;
1220 :
1221 : /// flag indicating if local data may change
1222 : static const bool staticLocalData = true;
1223 :
1224 : protected:
1225 : struct MidID
1226 : {
1227 0 : MidID() : dim(0), id(0) {};
1228 0 : MidID(size_t dim_, size_t id_) : dim(dim_), id(id_) {};
1229 : size_t dim;
1230 : size_t id;
1231 : };
1232 :
1233 : public:
1234 0 : class BF
1235 : {
1236 : private:
1237 : // let outer class access private members
1238 : friend class FV1ManifoldGeometry<TElem, TWorldDim>;
1239 :
1240 : // number of integration points
1241 : static const size_t m_numIP = 1;
1242 :
1243 : // max number of corners of bf
1244 : static const size_t numCorners = fv1_traits<ref_elem_type, dim>::NumCornersOfSCV;
1245 :
1246 : public:
1247 0 : BF() {};
1248 :
1249 : /// node id that this bf is associated to
1250 0 : inline size_t node_id() const {return nodeId;}
1251 :
1252 : /// number of integration points
1253 0 : inline size_t num_ip() const {return m_numIP;}
1254 :
1255 : /// local integration point of bf
1256 0 : inline const MathVector<dim>& local_ip() const
1257 0 : {return vLocPos[0];} // <-- always the vertex
1258 :
1259 : /// global integration point
1260 0 : inline const MathVector<worldDim>& global_ip() const
1261 0 : {return vGloPos[0];} // <-- here too
1262 :
1263 : /// volume of bf
1264 0 : inline number volume() const {return vol;}
1265 :
1266 : /// number of corners, that bound the bf
1267 0 : inline size_t num_corners() const {return numCorners;}
1268 :
1269 : /// return local position of corner number i
1270 0 : inline const MathVector<dim>& local_corner(size_t i) const
1271 0 : {UG_ASSERT(i < num_corners(), "Invalid index."); return vLocPos[i];}
1272 :
1273 : /// return global position of corner number i
1274 0 : inline const MathVector<worldDim>& global_corner(size_t i) const
1275 0 : {UG_ASSERT(i < num_corners(), "Invalid index."); return vGloPos[i];}
1276 :
1277 : /// number of shape functions
1278 0 : inline size_t num_sh() const {return vShape.size();}
1279 :
1280 : /// value of shape function i in integration point
1281 0 : inline number shape(size_t i, size_t ip) const
1282 0 : {UG_ASSERT(ip < num_ip(), "Invalid index"); return vShape[i];}
1283 :
1284 : private:
1285 : size_t nodeId; // id of associated node
1286 :
1287 : // CORNERS: ordering is:
1288 : // 1D: edgeMidPoint, CenterOfElement
1289 : // 2D: edgeMidPoint, Side one, CenterOfElement, Side two
1290 : MathVector<dim> vLocPos[numCorners]; // local position of node
1291 : MathVector<worldDim> vGloPos[numCorners]; // global position of node
1292 : MidID midId[numCorners]; // dimension and id of object, whose midpoint bounds the scv
1293 :
1294 : // IPs & shapes
1295 : std::vector<number> vShape; // shapes at ip
1296 :
1297 : number vol;
1298 : };
1299 :
1300 : protected:
1301 0 : void copy_local_corners(BF& bf)
1302 : {
1303 0 : for (size_t i = 0; i < bf.num_corners(); ++i)
1304 : {
1305 0 : const size_t dim = bf.midId[i].dim;
1306 0 : const size_t id = bf.midId[i].id;
1307 : bf.vLocPos[i] = m_locMid[dim][id];
1308 : }
1309 0 : }
1310 :
1311 0 : void copy_global_corners(BF& bf)
1312 : {
1313 0 : for (size_t i = 0; i < bf.num_corners(); ++i)
1314 : {
1315 0 : const size_t dim = bf.midId[i].dim;
1316 0 : const size_t id = bf.midId[i].id;
1317 : bf.vGloPos[i] = m_gloMid[dim][id];
1318 : }
1319 0 : }
1320 :
1321 : std::vector<MathVector<dim> > m_vLocBFIP;
1322 : std::vector<MathVector<worldDim> > m_vGlobBFIP;
1323 :
1324 :
1325 : public:
1326 : /// constructor
1327 : FV1ManifoldGeometry();
1328 :
1329 : /// update data for given element
1330 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
1331 : const ISubsetHandler* ish = NULL);
1332 :
1333 : /// get the element
1334 0 : GridObject* elem() const {return m_pElem;}
1335 :
1336 : /// get vector of corners for current element
1337 0 : const MathVector<worldDim>* corners() const {return m_gloMid[0];}
1338 :
1339 : /// number of BoundaryFaces
1340 0 : inline size_t num_bf() const {return m_numBF;}
1341 :
1342 : /// const access to Boundary Face number i
1343 0 : inline const BF& bf(size_t i) const
1344 0 : {UG_ASSERT(i < num_bf(), "Invalid Index."); return m_vBF[i];}
1345 :
1346 : /// returns all ips of scvf as they appear in scv loop
1347 0 : const MathVector<worldDim>* bf_global_ips() const {return &m_vGlobBFIP[0];}
1348 :
1349 : /// returns number of all BF ips
1350 0 : size_t num_bf_global_ips() const {return m_vGlobBFIP.size();}
1351 :
1352 : /// returns all ips of BF as they appear in scv loop
1353 0 : const MathVector<dim>* bf_local_ips() const {return &m_vLocBFIP[0];}
1354 :
1355 : /// returns number of all BF ips
1356 0 : size_t num_bf_local_ips() const {return m_vLocBFIP.size();}
1357 :
1358 : private:
1359 : // pointer to current element
1360 : GridObject* m_pElem;
1361 :
1362 : // local and global geom object midpoints for each dimension
1363 : MathVector<dim> m_locMid[dim+1][m_numBF];
1364 : MathVector<worldDim> m_gloMid[dim+1][m_numBF];
1365 :
1366 : // BndFaces
1367 : BF m_vBF[m_numBF];
1368 :
1369 : // Reference Mapping
1370 : ReferenceMapping<ref_elem_type, worldDim> m_rMapping;
1371 :
1372 : // Reference Element
1373 : const ref_elem_type& m_rRefElem;
1374 : };
1375 :
1376 : }
1377 :
1378 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FV1_GEOMETRY__ */
|
