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