Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Christian Wehner
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /*
34 : * Node centered finite volume geometry for Crouzeix-Raviart-Elements
35 : */
36 :
37 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__CR_FINITE_VOLUME_GEOMETRY__
38 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__CR_FINITE_VOLUME_GEOMETRY__
39 :
40 : // extern libraries
41 : #include <cmath>
42 : #include <map>
43 : #include <vector>
44 :
45 : // other ug4 modules
46 : #include "common/common.h"
47 :
48 : // library intern includes
49 : #include "lib_disc/reference_element/reference_element.h"
50 : #include "lib_disc/reference_element/reference_mapping.h"
51 : #include "lib_disc/reference_element/reference_mapping_provider.h"
52 : #include "lib_disc/reference_element/reference_element_traits.h"
53 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
54 : #include "lib_disc/local_finite_element/crouzeix-raviart/crouzeix_raviart.h"
55 : #include "lib_disc/quadrature/gauss/gauss_quad.h"
56 : #include "lib_disc/common/geometry_util.h"
57 : #include "lib_disc/domain_util_impl.h"
58 :
59 : #include "fv_geom_base.h"
60 : #include "fv_util.h"
61 : #include "hfvcr_geom.h"
62 :
63 : namespace ug{
64 :
65 : /* crfv traits */
66 :
67 : template <int TWorldDim,int nrfaceco> struct crfv_traits
68 : {
69 : typedef void scv_type;
70 : typedef void face_type;
71 : static const size_t maxNumSCVF;
72 : static const size_t maxNumSCV;
73 : static const size_t maxNSH;
74 : static const size_t maxNumCo;
75 : };
76 :
77 : template <> struct crfv_traits<1, 1>
78 : {
79 : typedef ReferenceEdge scv_type;
80 : typedef ReferenceVertex face_type;
81 : static const size_t maxNumSCVF = 1;
82 : static const size_t maxNumSCV = 2;
83 : static const size_t maxNSH = maxNumSCV;
84 : static const size_t maxNumCo = 2;
85 : };
86 :
87 : template <> struct crfv_traits<1, 2>
88 : {
89 : typedef ReferenceEdge scv_type;
90 : typedef ReferenceVertex face_type;
91 : static const size_t maxNumSCVF = 1;
92 : static const size_t maxNumSCV = 2;
93 : static const size_t maxNSH = maxNumSCV;
94 : static const size_t maxNumCo = 2;
95 : };
96 :
97 : template <> struct crfv_traits<2, 2>
98 : {
99 : typedef ReferenceTriangle scv_type;
100 : typedef ReferenceEdge face_type;
101 : static const size_t maxNumSCVF = 4;
102 : static const size_t maxNumSCV = 4;
103 : static const size_t maxNSH = maxNumSCV;
104 : static const size_t maxNumCo = 4;
105 : };
106 :
107 : template <> struct crfv_traits<2, 3>
108 : {
109 : typedef ReferenceTriangle scv_type;
110 : typedef ReferenceEdge face_type;
111 : static const size_t maxNumSCVF = 4;
112 : static const size_t maxNumSCV = 4;
113 : static const size_t maxNSH = maxNumSCV;
114 : static const size_t maxNumCo = 4;
115 : };
116 :
117 : template <> struct crfv_traits<3, 3>
118 : {
119 : typedef ReferenceTetrahedron scv_type;
120 : typedef ReferenceTriangle face_type;
121 : static const size_t maxNumSCVF = 10;
122 : static const size_t maxNumSCV = 6;
123 : static const size_t maxNSH = maxNumSCV;
124 : static const size_t maxNumCo = 8;
125 : };
126 :
127 : template <> struct crfv_traits<3, 4>
128 : {
129 : typedef ReferencePyramid scv_type;
130 : typedef ReferenceQuadrilateral face_type;
131 : static const size_t maxNumSCVF = 10;
132 : static const size_t maxNumSCV = 6;
133 : static const size_t maxNSH = maxNumSCV;
134 : static const size_t maxNumCo = 8;
135 : };
136 :
137 : ////////////////////////////////////////////////////////////////////////////////
138 : ////////////////////////////////////////////////////////////////////////////////
139 : // Dim-dependent Finite Volume Geometry
140 : ////////////////////////////////////////////////////////////////////////////////
141 : ////////////////////////////////////////////////////////////////////////////////
142 :
143 : ///
144 : /// Geometry and shape functions for 1st order Vertex-Centered Finite Volume with Crouzeix-Raviart elements
145 : //
146 : // basics:
147 : // Unknowns are centered in element vertices(1d)/edges(2d)/faces(3d)
148 : // scv i corners are given by face i corners + element barycenter
149 : // scvf i corners are given by edge i + element barycenter
150 : //
151 : // important differences to standard linear Lagrange FV class:
152 : //
153 : // new scv class parameter "Normal" is normal on associated face (this is computed for
154 : // Navier-Stokes-Crouzeix-Raviart-staggered grid continuity equation discretization)
155 : // midpoint computations and related parameters eliminated (not necessary for scv computation)
156 : // parameter nodeID refers to node located in edge/face for the unknown, not to corner node
157 : // ip of scv is not corner but associated edge/face midpoint
158 : // location of unknowns (edge/face barycenter) is stored in m_vLocUnkCoords, m_vGlobUnkCoords
159 : //
160 : /**
161 : * \tparam TDim reference element dim
162 : * \tparam TWorldDim (physical) world dimension
163 : */
164 : template <int TDim, int TWorldDim = TDim>
165 : class DimCRFVGeometry : public FVGeometryBase
166 : {
167 : public:
168 : /// used traits
169 : typedef crfv_traits<TDim, TDim> traits;
170 : typedef hcrfv_traits<TDim, TDim> htraits;
171 :
172 : public:
173 : /// dimension of reference element
174 : static const int dim = TDim;
175 :
176 : /// dimension of world
177 : static const int worldDim = TWorldDim;
178 :
179 : /// Hanging node flag: this Geometry does not support hanging nodes
180 : static const bool usesHangingNodes = false;
181 :
182 : /// flag indicating if local data may change
183 : static const bool staticLocalData = false;
184 :
185 : public:
186 : /// order
187 : static const int order = 1;
188 :
189 : /// number of SubControlVolumes
190 : static const size_t maxNumSCV = htraits::maxNumSCV;
191 :
192 : /// traits
193 : typedef typename crfv_traits<TWorldDim,TWorldDim>::scv_type scv_type0;
194 : typedef typename crfv_traits<TWorldDim,TWorldDim>::face_type face_type0;
195 : typedef typename crfv_traits<TWorldDim,TWorldDim+1>::scv_type scv_type1;
196 : typedef typename crfv_traits<TWorldDim,TWorldDim+1>::face_type face_type1;
197 :
198 : /// max number of SubControlVolumeFaces
199 : static const size_t maxNumSCVF = htraits::maxNumSCVF;
200 :
201 : /// max number of shape functions
202 : static const size_t maxNSH = htraits::maxNSH;
203 :
204 : static const size_t maxNumCo = traits::maxNumCo;
205 :
206 : /// number of integration points
207 : static const size_t nip = 1;
208 :
209 : // local/global element barycenter
210 : MathVector<dim> localBary;
211 : MathVector<worldDim> globalBary;
212 :
213 : public:
214 : /// Sub-Control Volume Face structure
215 : /**
216 : * Each finite element is cut by several sub-control volume faces. The idea
217 : * is that the "dual" skeleton formed by the sub control volume faces of
218 : * all elements again gives rise to a regular mesh with closed
219 : * (lipschitz-bounded) control volumes. The SCVF are the boundary of the
220 : * control volume. In computation the flux over each SCVF must be the same
221 : * in both directions over the face in order to guarantee the conservation
222 : * property.
223 : */
224 : class SCVF
225 : {
226 : public:
227 : /// Number of corners of scvf
228 : static const size_t numCo = dim;
229 :
230 : public:
231 0 : SCVF() {}
232 :
233 : /// index of SubControlVolume on one side of the scvf
234 0 : inline size_t from() const {return From;}
235 :
236 : /// index of SubControlVolume on one side of the scvf
237 0 : inline size_t to() const {return To;}
238 :
239 : /// number of integration points on scvf
240 0 : inline size_t num_ip() const {return nip;}
241 :
242 : /// local integration point of scvf
243 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
244 :
245 : /// global integration point of scvf
246 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
247 :
248 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
249 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
250 :
251 : /// Transposed Inverse of Jacobian in integration point
252 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
253 :
254 : /// Determinant of Jacobian in integration point
255 0 : inline number detJ() const {return detj;}
256 :
257 : /// number of shape functions
258 0 : inline size_t num_sh() const {return numSH;}
259 :
260 : /// value of shape function i in integration point
261 0 : inline number shape(size_t sh) const {return vShape[sh];}
262 :
263 : /// vector of shape functions in ip point
264 0 : inline const number* shape_vector() const {return vShape;}
265 :
266 : /// value of local gradient of shape function i in integration point
267 0 : inline const MathVector<dim>& local_grad(size_t sh) const
268 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
269 :
270 : /// vector of local gradients in ip point
271 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
272 :
273 : /// value of global gradient of shape function i in integration point
274 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
275 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
276 :
277 : /// vector of gloabl gradients in ip point
278 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
279 :
280 : /// number of corners, that bound the scvf
281 0 : inline size_t num_corners() const {return numCo;}
282 :
283 : /// return local corner number i
284 0 : inline const MathVector<dim>& local_corner(size_t co) const
285 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
286 :
287 : /// return glbal corner number i
288 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
289 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
290 :
291 : public:
292 0 : inline void operator = (const SCVF &s)
293 : {
294 0 : From = s.from();
295 0 : To = s.to();
296 : localIP = s.local_ip();
297 : globalIP = s.global_ip();
298 : Normal = s.normal();
299 : JtInv = s.JTInv();
300 0 : detj = s.detJ();
301 0 : numSH = s.num_sh();
302 0 : for (size_t i=0;i<numSH;i++){
303 : vLocalGrad[i] = s.global_grad(i);
304 : vGlobalGrad[i] = s.global_grad(i);
305 0 : vShape[i] = s.shape(i);
306 : }
307 0 : for (size_t i=0;i<numCo;i++){
308 : vLocPos[i] = s.local_corner(i);
309 : vGloPos[i] = s.global_corner(i);
310 : }
311 0 : }
312 :
313 : private:
314 : // let outer class access private members
315 : friend class DimCRFVGeometry<dim, worldDim>;
316 :
317 : // This scvf separates the scv with the ids given in "from" and "to"
318 : // The computed normal points in direction from->to
319 : size_t From, To;
320 :
321 : // The normal on the SCVF pointing (from -> to)
322 : MathVector<worldDim> Normal; // normal (incl. area)
323 :
324 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
325 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
326 :
327 : // scvf part
328 : MathVector<dim> localIP; // local integration point
329 : MathVector<worldDim> globalIP; // global integration point
330 :
331 : // shapes and derivatives
332 : size_t numSH;
333 : number vShape[maxNSH]; // shapes at ip
334 : MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
335 : MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
336 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
337 : number detj; // Jacobian det at ip
338 : };
339 :
340 : /// sub control volume structure
341 : class SCV
342 : {
343 : public:
344 : /// Number of corners of scv
345 : static const size_t maxNumCo = 5;
346 :
347 : public:
348 0 : SCV() : numCorners(maxNumCo) {};
349 :
350 : /// volume of scv
351 0 : inline number volume() const {return Vol;}
352 :
353 : /// number of corners, that bound the scvf
354 0 : inline size_t num_corners() const {return numCorners;}
355 :
356 : /// return local corner number i
357 0 : inline const MathVector<dim>& local_corner(size_t co) const
358 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
359 :
360 : /// return global corner number i
361 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
362 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
363 :
364 : /// node id that this scv is associated to
365 0 : inline size_t node_id() const {return nodeID;}
366 :
367 : /// number of integration points
368 0 : inline size_t num_ip() const {return nip;}
369 :
370 : /// local integration point of scv
371 0 : inline const MathVector<dim>& local_ip() const {return vLocIP;}
372 :
373 : /// global integration point
374 0 : inline const MathVector<worldDim>& global_ip() const {return vGlobIP;}
375 :
376 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
377 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
378 :
379 : /// Transposed Inverse of Jacobian in integration point
380 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
381 :
382 : /// Determinant of Jacobian in integration point
383 0 : inline number detJ() const {return detj;}
384 :
385 : /// number of shape functions
386 0 : inline size_t num_sh() const {return numSH;}
387 :
388 : /// value of shape function i in integration point
389 0 : inline number shape(size_t sh) const {return vShape[sh];}
390 :
391 : /// vector of shape functions in ip point
392 0 : inline const number* shape_vector() const {return vShape;}
393 :
394 : /// value of local gradient of shape function i in integration point
395 0 : inline const MathVector<dim>& local_grad(size_t sh) const
396 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
397 :
398 : /// vector of local gradients in ip point
399 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
400 :
401 : /// value of global gradient of shape function i in integration point
402 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
403 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
404 :
405 : /// vector of global gradients in ip point
406 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
407 :
408 : private:
409 : // let outer class access private members
410 : friend class DimCRFVGeometry<dim, worldDim>;
411 :
412 : // The normal on the associated face pointing outward
413 : MathVector<worldDim> Normal; // normal (incl. area)
414 :
415 : // node id of associated node
416 : size_t nodeID;
417 :
418 : // volume of scv
419 : number Vol;
420 :
421 : // number of corners of this element
422 : int numCorners;
423 :
424 : // local and global positions of this element
425 : MathVector<dim> vLocPos[maxNumCo]; // local position of node
426 : MathVector<worldDim> vGloPos[maxNumCo]; // global position of node
427 : MathVector<dim> vLocIP;
428 : MathVector<dim> vGlobIP;
429 :
430 : // shapes and derivatives
431 : size_t numSH;
432 : number vShape[maxNSH]; // shapes at ip
433 : MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
434 : MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
435 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
436 : number detj; // Jacobian det at ip
437 : };
438 :
439 : /// boundary face
440 : class BF
441 : {
442 : public:
443 : /// max number of corners of bf
444 : static const size_t maxNumCo=4;
445 :
446 : public:
447 0 : BF() {}
448 :
449 : /// index of SubControlVolume of the bf
450 0 : inline size_t node_id() const {return nodeID;}
451 :
452 : /// number of integration points on bf
453 0 : inline size_t num_ip() const {return nip;}
454 :
455 : /// local integration point of bf
456 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
457 :
458 : /// global integration point of bf
459 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
460 :
461 : /// outer normal on bf. Norm is equal to area
462 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
463 :
464 : /// volume of bf
465 0 : inline number volume() const {return Vol;}
466 :
467 : /// Transposed Inverse of Jacobian in integration point
468 0 : inline const MathMatrix<worldDim, dim>& JTInv() const {return JtInv;}
469 :
470 : /// Determinant of Jacobian in integration point
471 0 : inline number detJ() const {return detj;}
472 :
473 : /// number of shape functions
474 0 : inline size_t num_sh() const {return numSH;}
475 :
476 : /// value of shape function i in integration point
477 0 : inline number shape(size_t sh) const
478 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vShape[sh];}
479 :
480 : /// vector of local gradients in ip point
481 0 : inline const number* shape_vector() const {return vShape;}
482 :
483 : /// value of local gradient of shape function i in integration point
484 0 : inline const MathVector<dim>& local_grad(size_t sh) const
485 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
486 :
487 : /// vector of local gradients in ip point
488 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
489 :
490 : /// value of global gradient of shape function i in integration point
491 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
492 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
493 :
494 : /// vector of global gradients in ip point
495 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
496 :
497 : /// number of corners, that bound the scvf
498 0 : inline size_t num_corners() const {return numCo;}
499 :
500 : /// return local corner number i
501 0 : inline const MathVector<dim>& local_corner(size_t co) const
502 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
503 :
504 : /// return global corner number i
505 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
506 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
507 :
508 : private:
509 : /// let outer class access private members
510 : friend class DimCRFVGeometry<dim, worldDim>;
511 :
512 : // id of scv this bf belongs to
513 : size_t nodeID;
514 :
515 : size_t numCo;
516 :
517 : MathVector<dim> vLocPos[maxNumCo]; // local corners of bf
518 : MathVector<worldDim> vGloPos[maxNumCo]; // global corners of bf
519 :
520 : // scvf part
521 : MathVector<dim> localIP; // local integration point
522 : MathVector<worldDim> globalIP; // global integration point
523 : MathVector<worldDim> Normal; // normal (incl. area)
524 : number Vol; // volume of bf
525 :
526 : // shapes and derivatives
527 : size_t numSH;
528 : number vShape[maxNSH]; // shapes at ip
529 : MathVector<dim> vLocalGrad[maxNSH]; // local grad at ip
530 : MathVector<worldDim> vGlobalGrad[maxNSH]; // global grad at ip
531 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
532 : number detj; // Jacobian det at ip
533 : };
534 :
535 : class CONSTRAINED_DOF
536 : {
537 : public:
538 : static const size_t maxNumConstrainingDofs = 4;
539 0 : inline size_t constraining_dofs_index(size_t i) const{
540 0 : return cDofInd[i];
541 : }
542 0 : inline number constraining_dofs_weight(size_t i) const{
543 0 : return cDofWeights[i];
544 : }
545 0 : inline size_t index() const{
546 0 : return i;
547 : }
548 0 : inline size_t num_constraining_dofs() const{
549 0 : return numConstrainingDofs;
550 : }
551 : private:
552 : // let outer class access private members
553 : friend class DimCRFVGeometry<dim, worldDim>;
554 :
555 : // constraining dofs indices
556 : size_t cDofInd[maxNumConstrainingDofs];
557 : // weights
558 : number cDofWeights[maxNumConstrainingDofs];
559 : // local index of dof in element
560 : size_t i;
561 : // nr of constraining dofs
562 : size_t numConstrainingDofs;
563 : };
564 :
565 : class HandledEdge
566 : {
567 : public:
568 : size_t index;
569 : size_t associatedSCV[2];
570 : size_t scvfIndex;
571 : // indicates if the handled side is from or to
572 : bool from;
573 : };
574 :
575 : public:
576 : /// construct object and initialize local values and sizes
577 0 : DimCRFVGeometry() : m_pElem(NULL), m_roid(ROID_UNKNOWN) {};
578 :
579 : /// update local data
580 : void update_local_data();
581 :
582 : /// update data for given element
583 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
584 : const ISubsetHandler* ish = NULL);
585 :
586 : /// update data for given element
587 : void update_hanging(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
588 : const ISubsetHandler* ish = NULL,bool keepSCV=false,bool keepSCVF=false);
589 :
590 : /// update data for given element
591 : void update_geometric_data(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
592 : const ISubsetHandler* ish = NULL);
593 :
594 : /// update boundary data for given element
595 : void update_boundary_faces(GridObject* elem,
596 : const MathVector<worldDim>* vCornerCoords,
597 : const ISubsetHandler* ish = NULL);
598 :
599 : /// number of SubControlVolumeFaces
600 0 : inline size_t num_scvf() const {return m_numSCVF;};
601 :
602 : /// const access to SubControlVolumeFace number i
603 0 : inline const SCVF& scvf(size_t i) const
604 0 : {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
605 :
606 : /// number of SubControlVolumes
607 0 : inline size_t num_scv() const {return m_numSCV;}
608 :
609 : /// const access to SubControlVolume number i
610 0 : inline const SCV& scv(size_t i) const
611 0 : {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
612 :
613 : /// number of constrained dofs
614 0 : inline size_t num_constrained_dofs() const {return m_numConstrainedDofs;}
615 :
616 : /// const access to constrained dof i
617 0 : inline const CONSTRAINED_DOF& constrained_dof(size_t i) const
618 0 : {UG_ASSERT(i < m_numConstrainedDofs, "Invalid Index."); return m_vCD[i];}
619 :
620 : /// number of constrained dof scvfs
621 0 : inline size_t num_constrained_scvf(){return m_numConstrainedSCVF;}
622 :
623 : /// number of shape functions
624 0 : inline size_t num_sh() const {return m_nsh;};
625 :
626 : public:
627 : /// returns number of all scvf ips
628 0 : size_t num_scvf_ips() const {return m_numSCVF;}
629 :
630 : /// returns all ips of scvf as they appear in scv loop
631 0 : const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
632 :
633 : /// returns all ips of scvf as they appear in scv loop
634 0 : const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
635 :
636 : /// returns number of all scv ips
637 0 : size_t num_scv_ips() const {return m_numSCV;}
638 :
639 : /// returns all ips of scv as they appear in scv loop
640 0 : const MathVector<dim>* scv_local_ips() const {return m_vLocUnkCoords;}
641 :
642 : /// returns all ips of scv as they appear in scv loop
643 0 : const MathVector<worldDim>* scv_global_ips() const {return m_vGlobUnkCoords;}
644 :
645 : /// returns local barycenter
646 0 : const MathVector<dim> local_bary() const {return localBary;}
647 :
648 : /// returns global barycenter
649 0 : const MathVector<worldDim> global_bary() const {return globalBary;}
650 :
651 : protected:
652 : // global and local ips on SCVF
653 : MathVector<worldDim> m_vGlobSCVF_IP[maxNumSCVF];
654 : MathVector<dim> m_vLocSCVF_IP[maxNumSCVF];
655 : // coord of location for unknowns in faces (edge/face barycenter)
656 : MathVector<worldDim> m_vGlobUnkCoords[maxNumSCV];
657 : MathVector<dim> m_vLocUnkCoords[maxNumSCV];
658 :
659 : public:
660 : /// add subset that is interpreted as boundary subset.
661 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
662 :
663 : /// removes subset that is interpreted as boundary subset.
664 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
665 :
666 : /// reset all boundary subsets
667 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
668 :
669 : /// number of registered boundary subsets
670 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
671 :
672 : /// number of all boundary faces
673 0 : inline size_t num_bf() const
674 : {
675 : typename std::map<int, std::vector<BF> >::const_iterator it;
676 : size_t num = 0;
677 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
678 0 : num += (*it).second.size();
679 0 : return num;
680 : }
681 :
682 : /// number of boundary faces on subset 'subsetIndex'
683 0 : inline size_t num_bf(int subsetIndex) const
684 : {
685 : typename std::map<int, std::vector<BF> >::const_iterator it;
686 : it = m_mapVectorBF.find(subsetIndex);
687 0 : if(it == m_mapVectorBF.end()) return 0;
688 0 : else return (*it).second.size();
689 : }
690 :
691 : /// returns the boundary face i for subsetIndex
692 0 : inline const BF& bf(int subsetIndex, size_t i) const
693 : {
694 : UG_ASSERT(i < num_bf(subsetIndex), "Invalid index.");
695 : typename std::map<int, std::vector<BF> >::const_iterator it;
696 : it = m_mapVectorBF.find(subsetIndex);
697 : UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
698 0 : return (*it).second[i];
699 : }
700 :
701 : /// returns reference to vector of boundary faces for subsetIndex
702 0 : inline const std::vector<BF>& bf(int subsetIndex) const
703 : {
704 : typename std::map<int, std::vector<BF> >::const_iterator it;
705 : it = m_mapVectorBF.find(subsetIndex);
706 : UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
707 0 : return (*it).second;
708 : }
709 :
710 : protected:
711 : std::map<int, std::vector<BF> > m_mapVectorBF;
712 :
713 : /// constrained Dofs
714 : CONSTRAINED_DOF m_vCD[maxNumSCV];
715 :
716 : std::vector<HandledEdge> handledEdges;
717 :
718 : private:
719 : /// pointer to current element
720 : GridObject* m_pElem;
721 :
722 : /// current reference object id
723 : ReferenceObjectID m_roid;
724 :
725 : /// current reference element
726 : DimReferenceElement<dim> m_rRefElem;
727 :
728 : /// current number of scv
729 : size_t m_numSCV;
730 :
731 : /// current number of scvf
732 : size_t m_numSCVF;
733 :
734 : /// current number of shape functions
735 : size_t m_nsh;
736 :
737 : MathVector<dim> m_ipCoord[maxNumSCVF];
738 : MathVector<dim> cornerCoord[maxNumCo];
739 :
740 : /// SubControlVolumeFaces
741 : SCVF m_vSCVF[maxNumSCVF];
742 :
743 : /// SubControlVolumeFaces for constrained dofs
744 : SCVF m_vConstrainedSCVF[maxNumSCVF];
745 :
746 : /// SubControlVolumes
747 : SCV m_vSCV[maxNumSCV];
748 :
749 : /// numDofs number of all dofs including constraining and constrained dofs
750 : size_t m_numDofs;
751 :
752 : size_t m_numConstrainedDofs;
753 :
754 : size_t m_numConstrainedSCVF;
755 :
756 : static const size_t deleted = 117;
757 : };
758 :
759 : ////////////////////////////////////////////////////////////////////////////////
760 : ////////////////////////////////////////////////////////////////////////////////
761 : // Dimension-indipendent Finite Volume Geometry
762 : ////////////////////////////////////////////////////////////////////////////////
763 : ////////////////////////////////////////////////////////////////////////////////
764 :
765 : template < typename TElem, int TWorldDim>
766 : class CRFVGeometry : public FVGeometryBase
767 : {
768 : public:
769 : /// type of element
770 : typedef TElem elem_type;
771 :
772 : /// type of reference element
773 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
774 :
775 : /// dimension of reference element
776 : static const int dim = ref_elem_type::dim;
777 :
778 : /// dimension of world
779 : static const int worldDim = TWorldDim;
780 :
781 : /// Hanging node flag: this Geometry does not support hanging nodes
782 : static const bool usesHangingNodes = false;
783 :
784 : /// flag indicating if local data may change
785 : static const bool staticLocalData = true;
786 :
787 : /// type of Shape function used
788 : typedef CrouzeixRaviartLSFS<ref_elem_type> local_shape_fct_set_type;
789 :
790 : /// number of shape functions
791 : static const size_t nsh = local_shape_fct_set_type::nsh;
792 :
793 : /// number of SubControlVolumes
794 : static const size_t numSCV = nsh;
795 :
796 : /// number of SubControlVolumeFaces
797 : static const size_t numSCVF = ref_elem_type::numEdges;
798 :
799 : static const size_t maxNumSCVF = ref_elem_type::numEdges;
800 :
801 : public:
802 : /// order
803 : static const int order = 1;
804 :
805 : /// traits
806 : typedef typename crfv_traits<TWorldDim,TWorldDim>::scv_type scv_type0;
807 : typedef typename crfv_traits<TWorldDim,TWorldDim>::face_type face_type0;
808 : typedef typename crfv_traits<TWorldDim,TWorldDim+1>::scv_type scv_type1;
809 : typedef typename crfv_traits<TWorldDim,TWorldDim+1>::face_type face_type1;
810 :
811 : /// number of integration points
812 : static const size_t nip = 1;
813 :
814 : // local/global element barycenter
815 : MathVector<dim> localBary;
816 : MathVector<worldDim> globalBary;
817 :
818 : public:
819 : /// Sub-Control Volume Face structure
820 : /**
821 : * Each finite element is cut by several sub-control volume faces. The idea
822 : * is that the "dual" skeleton formed by the sub control volume faces of
823 : * all elements again gives rise to a regular mesh with closed
824 : * (lipschitz-bounded) control volumes. The SCVF are the boundary of the
825 : * control volume. In computation the flux over each SCVF must be the same
826 : * in both directions over the face in order to guarantee the conservation
827 : * property.
828 : */
829 : class SCVF
830 : {
831 : public:
832 : /// Number of corners of scvf
833 : static const size_t numCo = dim;
834 :
835 : public:
836 0 : SCVF() {}
837 :
838 : /// index of SubControlVolume on one side of the scvf
839 0 : inline size_t from() const {return From;}
840 :
841 : /// index of SubControlVolume on one side of the scvf
842 0 : inline size_t to() const {return To;}
843 :
844 : /// number of integration points on scvf
845 0 : inline size_t num_ip() const {return nip;}
846 :
847 : /// local integration point of scvf
848 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
849 :
850 : /// global integration point of scvf
851 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
852 :
853 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
854 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
855 :
856 : /// Transposed Inverse of Jacobian in integration point
857 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
858 :
859 : /// Determinant of Jacobian in integration point
860 0 : inline number detJ() const {return detj;}
861 :
862 : /// number of shape functions
863 0 : inline size_t num_sh() const {return nsh;}
864 :
865 : /// value of shape function i in integration point
866 0 : inline number shape(size_t sh) const {return vShape[sh];}
867 :
868 : /// vector of shape functions in ip point
869 0 : inline const number* shape_vector() const {return vShape;}
870 :
871 : /// value of local gradient of shape function i in integration point
872 0 : inline const MathVector<dim>& local_grad(size_t sh) const
873 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
874 :
875 : /// vector of local gradients in ip point
876 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
877 :
878 : /// value of global gradient of shape function i in integration point
879 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
880 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
881 :
882 : /// vector of global gradients in ip point
883 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
884 :
885 : /// number of corners, that bound the scvf
886 0 : inline size_t num_corners() const {return numCo;}
887 :
888 : /// return local corner number i
889 0 : inline const MathVector<dim>& local_corner(size_t co) const
890 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
891 :
892 : /// return glbal corner number i
893 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
894 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
895 :
896 : private:
897 : // let outer class access private members
898 : friend class CRFVGeometry<TElem, TWorldDim>;
899 :
900 : // This scvf separates the scv with the ids given in "from" and "to"
901 : // The computed normal points in direction from->to
902 : size_t From, To;
903 :
904 : // The normal on the SCVF pointing (from -> to)
905 : MathVector<worldDim> Normal; // normal (incl. area)
906 :
907 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
908 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
909 :
910 : // scvf part
911 : MathVector<dim> localIP; // local integration point
912 : MathVector<worldDim> globalIP; // global integration point
913 :
914 : // shapes and derivatives
915 : number vShape[nsh]; // shapes at ip
916 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
917 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
918 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
919 : number detj; // Jacobian det at ip
920 : };
921 :
922 : /// sub control volume structure
923 : class SCV
924 : {
925 : public:
926 : /// Number of corners of scv
927 : static const size_t maxNumCo = 5;
928 :
929 : public:
930 0 : SCV() : numCorners(maxNumCo) {};
931 :
932 : /// volume of scv
933 0 : inline number volume() const {return Vol;}
934 :
935 : /// number of corners, that bound the scvf
936 0 : inline size_t num_corners() const {return numCorners;}
937 :
938 : /// return local corner number i
939 0 : inline const MathVector<dim>& local_corner(size_t co) const
940 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
941 :
942 : /// return global corner number i
943 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
944 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
945 :
946 : /// node id that this scv is associated to
947 0 : inline size_t node_id() const {return nodeID;}
948 :
949 : /// number of integration points
950 0 : inline size_t num_ip() const {return nip;}
951 :
952 : /// local integration point of scv
953 0 : inline const MathVector<dim>& local_ip() const {return vLocIP;}
954 :
955 : /// global integration point
956 0 : inline const MathVector<worldDim>& global_ip() const {return vGlobIP;}
957 :
958 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
959 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
960 :
961 : /// Transposed Inverse of Jacobian in integration point
962 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
963 :
964 : /// Determinant of Jacobian in integration point
965 0 : inline number detJ() const {return detj;}
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 sh) const {return vShape[sh];}
972 :
973 : /// vector of shape functions in ip point
974 0 : inline const number* shape_vector() const {return vShape;}
975 :
976 : /// value of local gradient of shape function i in integration point
977 0 : inline const MathVector<dim>& local_grad(size_t sh) const
978 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
979 :
980 : /// vector of local gradients in ip point
981 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
982 :
983 : /// value of global gradient of shape function i in integration point
984 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
985 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
986 :
987 : /// vector of global gradients in ip point
988 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
989 :
990 : private:
991 : // let outer class access private members
992 : friend class CRFVGeometry<TElem, TWorldDim>;
993 :
994 : // The normal on the associated face pointing outward
995 : MathVector<worldDim> Normal; // normal (incl. area)
996 :
997 : // node id of associated node
998 : size_t nodeID;
999 :
1000 : // volume of scv
1001 : number Vol;
1002 :
1003 : // number of corners of this element
1004 : int numCorners;
1005 :
1006 : // local and global positions of this element
1007 : MathVector<dim> vLocPos[maxNumCo]; // local position of node
1008 : MathVector<worldDim> vGloPos[maxNumCo]; // global position of node
1009 :
1010 : MathVector<dim> vLocIP; // local position of node
1011 : MathVector<worldDim> vGlobIP; // global position of node
1012 :
1013 : // shapes and derivatives
1014 : number vShape[nsh]; // shapes at ip
1015 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
1016 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
1017 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
1018 : number detj; // Jacobian det at ip
1019 : };
1020 :
1021 : /// boundary face
1022 : class BF
1023 : {
1024 : public:
1025 : /// max number of corners of bf
1026 : static const size_t maxNumCo=4;
1027 :
1028 : public:
1029 0 : BF() {}
1030 :
1031 : /// index of SubControlVolume of the bf
1032 0 : inline size_t node_id() const {return nodeID;}
1033 :
1034 : /// number of integration points on bf
1035 0 : inline size_t num_ip() const {return nip;}
1036 :
1037 : /// local integration point of bf
1038 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
1039 :
1040 : /// global integration point of bf
1041 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
1042 :
1043 : /// outer normal on bf. Norm is equal to area
1044 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
1045 :
1046 : /// volume of bf
1047 0 : inline number volume() const {return Vol;}
1048 :
1049 : /// Transposed Inverse of Jacobian in integration point
1050 0 : inline const MathMatrix<worldDim, dim>& JTInv() const {return JtInv;}
1051 :
1052 : /// Determinant of Jacobian in integration point
1053 0 : inline number detJ() const {return detj;}
1054 :
1055 : /// number of shape functions
1056 0 : inline size_t num_sh() const {return nsh;}
1057 :
1058 : /// value of shape function i in integration point
1059 0 : inline number shape(size_t sh) const
1060 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vShape[sh];}
1061 :
1062 : /// vector of local gradients in ip point
1063 0 : inline const number* shape_vector() const {return vShape;}
1064 :
1065 : /// value of local gradient of shape function i in integration point
1066 0 : inline const MathVector<dim>& local_grad(size_t sh) const
1067 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
1068 :
1069 : /// vector of local gradients in ip point
1070 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
1071 :
1072 : /// value of global gradient of shape function i in integration point
1073 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
1074 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
1075 :
1076 : /// vector of global gradients in ip point
1077 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
1078 :
1079 : /// number of corners, that bound the scvf
1080 0 : inline size_t num_corners() const {return numCo;}
1081 :
1082 : /// return local corner number i
1083 0 : inline const MathVector<dim>& local_corner(size_t co) const
1084 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
1085 :
1086 : /// return global corner number i
1087 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
1088 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
1089 :
1090 : private:
1091 : /// let outer class access private members
1092 : friend class CRFVGeometry<TElem, TWorldDim>;
1093 :
1094 : // id of scv this bf belongs to
1095 : size_t nodeID;
1096 :
1097 : size_t numCo;
1098 :
1099 : MathVector<dim> vLocPos[maxNumCo]; // local corners of bf
1100 : MathVector<worldDim> vGloPos[maxNumCo]; // global corners of bf
1101 :
1102 : // scvf part
1103 : MathVector<dim> localIP; // local integration point
1104 : MathVector<worldDim> globalIP; // global integration point
1105 : MathVector<worldDim> Normal; // normal (incl. area)
1106 : number Vol; // volume of bf
1107 :
1108 : // shapes and derivatives
1109 : number vShape[nsh]; // shapes at ip
1110 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
1111 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
1112 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
1113 : number detj; // Jacobian det at ip
1114 : };
1115 :
1116 : public:
1117 : /// construct object and initialize local values and sizes
1118 : CRFVGeometry();
1119 :
1120 : /// update local data
1121 : void update_local_data();
1122 :
1123 : /// update data for given element
1124 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
1125 : const ISubsetHandler* ish = NULL);
1126 :
1127 : /// update boundary data for given element
1128 : void update_boundary_faces(GridObject* elem,
1129 : const MathVector<worldDim>* vCornerCoords,
1130 : const ISubsetHandler* ish = NULL);
1131 :
1132 0 : const MathVector<worldDim>* corners() const {return m_vCo;}
1133 :
1134 : /// number of SubControlVolumeFaces
1135 0 : inline size_t num_scvf() const {return numSCVF;};
1136 :
1137 : /// const access to SubControlVolumeFace number i
1138 0 : inline const SCVF& scvf(size_t i) const
1139 0 : {UG_ASSERT(i < num_scvf(), "Invalid Index."); return m_vSCVF[i];}
1140 :
1141 : /// number of SubControlVolumes
1142 0 : inline size_t num_scv() const {return numSCV;}
1143 :
1144 : /// const access to SubControlVolume number i
1145 0 : inline const SCV& scv(size_t i) const
1146 0 : {UG_ASSERT(i < num_scv(), "Invalid Index."); return m_vSCV[i];}
1147 :
1148 : /// number of shape functions
1149 0 : inline size_t num_sh() const {return nsh;};
1150 :
1151 : public:
1152 : /// returns number of all scvf ips
1153 0 : size_t num_scvf_ips() const {return numSCVF;}
1154 :
1155 : /// returns all ips of scvf as they appear in scv loop
1156 0 : const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
1157 :
1158 : /// returns all ips of scvf as they appear in scv loop
1159 0 : const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
1160 :
1161 : /// returns number of all scv ips
1162 0 : size_t num_scv_ips() const {return numSCV;}
1163 :
1164 : /// returns all ips of scv as they appear in scv loop
1165 0 : const MathVector<dim>* scv_local_ips() const {return m_vLocUnkCoords;}
1166 :
1167 : /// returns all ips of scv as they appear in scv loop
1168 0 : const MathVector<worldDim>* scv_global_ips() const {return m_vGlobUnkCoords;}
1169 :
1170 : /// returns local barycenter
1171 0 : const MathVector<dim> local_bary() const {return localBary;}
1172 :
1173 : /// returns global barycenter
1174 0 : const MathVector<worldDim> global_bary() const {return globalBary;}
1175 :
1176 : protected:
1177 : // global and local ips on SCVF
1178 : MathVector<worldDim> m_vGlobSCVF_IP[numSCVF];
1179 : MathVector<dim> m_vLocSCVF_IP[numSCVF];
1180 : // coord of location for unknowns in faces (edge/face barycenter)
1181 : MathVector<worldDim> m_vGlobUnkCoords[numSCV];
1182 : MathVector<dim> m_vLocUnkCoords[numSCV];
1183 :
1184 : static const size_t numMaxCo = 8;
1185 : // corner coordinates
1186 : MathVector<worldDim> m_vCo[numMaxCo];
1187 :
1188 : public:
1189 : /// add subset that is interpreted as boundary subset.
1190 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
1191 :
1192 : /// removes subset that is interpreted as boundary subset.
1193 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
1194 :
1195 : /// reset all boundary subsets
1196 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
1197 :
1198 : /// number of registered boundary subsets
1199 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
1200 :
1201 : /// number of all boundary faces
1202 0 : inline size_t num_bf() const
1203 : {
1204 : typename std::map<int, std::vector<BF> >::const_iterator it;
1205 : size_t num = 0;
1206 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
1207 0 : num += (*it).second.size();
1208 0 : return num;
1209 : }
1210 :
1211 : /// number of boundary faces on subset 'subsetIndex'
1212 0 : inline size_t num_bf(int subsetIndex) const
1213 : {
1214 : typename std::map<int, std::vector<BF> >::const_iterator it;
1215 : it = m_mapVectorBF.find(subsetIndex);
1216 0 : if(it == m_mapVectorBF.end()) return 0;
1217 0 : else return (*it).second.size();
1218 : }
1219 :
1220 : /// returns the boundary face i for subsetIndex
1221 0 : inline const BF& bf(int subsetIndex, size_t i) const
1222 : {
1223 : UG_ASSERT(i < num_bf(subsetIndex), "Invalid index.");
1224 : typename std::map<int, std::vector<BF> >::const_iterator it;
1225 : it = m_mapVectorBF.find(subsetIndex);
1226 : UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
1227 0 : return (*it).second[i];
1228 : }
1229 :
1230 : /// returns reference to vector of boundary faces for subsetIndex
1231 0 : inline const std::vector<BF>& bf(int subsetIndex) const
1232 : {
1233 : typename std::map<int, std::vector<BF> >::const_iterator it;
1234 : it = m_mapVectorBF.find(subsetIndex);
1235 : UG_ASSERT(it != m_mapVectorBF.end(), "Bnd Subset Index not requested.")
1236 0 : return (*it).second;
1237 : }
1238 :
1239 : protected:
1240 : std::map<int, std::vector<BF> > m_mapVectorBF;
1241 :
1242 : private:
1243 : MathVector<dim> m_ipCoord[numSCVF];
1244 :
1245 : /// SubControlVolumeFaces
1246 : SCVF m_vSCVF[numSCVF];
1247 :
1248 : /// SubControlVolumes
1249 : SCV m_vSCV[numSCV];
1250 :
1251 : /// pointer to current element
1252 : TElem* m_pElem;
1253 :
1254 : /// Reference Mapping
1255 : ReferenceMapping<ref_elem_type, worldDim> m_mapping;
1256 :
1257 : /// Reference Element
1258 : const ref_elem_type& m_rRefElem;
1259 :
1260 : /// Shape function set
1261 : const local_shape_fct_set_type& m_rTrialSpace;
1262 :
1263 : /// public hanging node interface for compatibility reasons
1264 : public:
1265 : class CONSTRAINED_DOF
1266 : {
1267 : public:
1268 0 : inline size_t constraining_dofs_index(size_t i) const {
1269 0 : return 0;
1270 : }
1271 0 : inline size_t constraining_dofs_weight(size_t i) const{
1272 0 : return 0;
1273 : }
1274 0 : inline size_t index() const{
1275 0 : return 0;
1276 : }
1277 0 : inline size_t num_constraining_dofs() const{
1278 0 : return 0;
1279 : }
1280 : };
1281 0 : inline size_t num_constrained_dofs() const {return 0;}
1282 0 : inline const CONSTRAINED_DOF& constrained_dof(size_t i) const
1283 0 : {return cd;}
1284 : private:
1285 : CONSTRAINED_DOF cd;
1286 : };
1287 :
1288 :
1289 : }
1290 :
1291 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__CR_FINITE_VOLUME_GEOMETRY__ */
|
