Line data Source code
1 : /*
2 : * Copyright (c) 2013-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__HANGING_CR_FINITE_VOLUME_GEOMETRY__
38 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__HANGING_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 :
62 : namespace ug{
63 :
64 : /* hcrfv traits */
65 :
66 : template <int TWorldDim,int nrfaceco> struct hcrfv_traits
67 : {
68 : typedef void scv_type;
69 : typedef void face_type;
70 : static const size_t maxNumSCVF;
71 : static const size_t maxNumSCV;
72 : static const size_t maxNSH;
73 : static const size_t maxNumCo;
74 : };
75 :
76 : template <> struct hcrfv_traits<2, 2>
77 : {
78 : typedef ReferenceTriangle scv_type;
79 : typedef ReferenceEdge face_type;
80 : static const size_t maxNumSCVF = 8;
81 : static const size_t maxNumSCV = 8;
82 : static const size_t maxNSH = maxNumSCV;
83 : static const size_t maxNumCo = 4;
84 : };
85 :
86 : template <> struct hcrfv_traits<2, 3>
87 : {
88 : typedef ReferenceTriangle scv_type;
89 : typedef ReferenceEdge face_type;
90 : static const size_t maxNumSCVF = 8;
91 : static const size_t maxNumSCV = 8;
92 : static const size_t maxNSH = maxNumSCV;
93 : static const size_t maxNumCo = 4;
94 : };
95 :
96 : template <> struct hcrfv_traits<3, 3>
97 : {
98 : typedef ReferenceTetrahedron scv_type;
99 : typedef ReferenceTriangle face_type;
100 : static const size_t maxNumSCVF = 40;
101 : static const size_t maxNumSCV = 24;
102 : static const size_t maxNSH = maxNumSCV;
103 : static const size_t maxNumCo = 8;
104 : };
105 :
106 : template <> struct hcrfv_traits<3, 4>
107 : {
108 : typedef ReferencePyramid scv_type;
109 : typedef ReferenceQuadrilateral face_type;
110 : static const size_t maxNumSCVF = 40;
111 : static const size_t maxNumSCV = 24;
112 : static const size_t maxNSH = maxNumSCV;
113 : static const size_t maxNumCo = 8;
114 : };
115 :
116 : ////////////////////////////////////////////////////////////////////////////////
117 : ////////////////////////////////////////////////////////////////////////////////
118 : // Dimension-indipendent Finite Volume Geometry
119 : ////////////////////////////////////////////////////////////////////////////////
120 : ////////////////////////////////////////////////////////////////////////////////
121 :
122 : template < typename TElem, int TWorldDim>
123 : class HCRFVGeometry : public FVGeometryBase
124 : {
125 : public:
126 : /// type of element
127 : typedef TElem elem_type;
128 :
129 : /// type of reference element
130 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
131 :
132 : /// dimension of reference element
133 : static const int dim = ref_elem_type::dim;
134 :
135 : /// dimension of world
136 : static const int worldDim = TWorldDim;
137 :
138 : /// Hanging node flag: this geometry supports hanging nodes
139 : static const bool usesHangingNodes = true;
140 :
141 : /// flag indicating if local data may change
142 : static const bool staticLocalData = true;
143 :
144 : /// type of Shape function used
145 : typedef CrouzeixRaviartLSFS<ref_elem_type> local_shape_fct_set_type;
146 :
147 : /// number of shape functions
148 : static const size_t nsh = local_shape_fct_set_type::nsh;
149 :
150 : /// number of SubControlVolumes
151 : static const size_t numNaturalSCV = nsh;
152 :
153 : /// number of SubControlVolumeFaces
154 : static const size_t numNaturalSCVF = ref_elem_type::numEdges;
155 :
156 : /// used traits
157 : typedef hcrfv_traits<dim,worldDim> traits;
158 :
159 : static const size_t maxNumSCV = traits::maxNumSCV;
160 :
161 : static const size_t maxNumSCVF = traits::maxNumSCVF;
162 :
163 : public:
164 : /// order
165 : static const int order = 1;
166 :
167 : /// traits
168 : typedef typename hcrfv_traits<TWorldDim,TWorldDim>::scv_type scv_type0;
169 : typedef typename hcrfv_traits<TWorldDim,TWorldDim>::face_type face_type0;
170 : typedef typename hcrfv_traits<TWorldDim,TWorldDim+1>::scv_type scv_type1;
171 : typedef typename hcrfv_traits<TWorldDim,TWorldDim+1>::face_type face_type1;
172 :
173 : /// number of integration points
174 : static const size_t nip = 1;
175 :
176 : // local/global element barycenter
177 : MathVector<dim> localBary;
178 : MathVector<worldDim> globalBary;
179 :
180 : public:
181 : /// Sub-Control Volume Face structure
182 : /**
183 : * Each finite element is cut by several sub-control volume faces. The idea
184 : * is that the "dual" skeleton formed by the sub control volume faces of
185 : * all elements again gives rise to a regular mesh with closed
186 : * (lipschitz-bounded) control volumes. The SCVF are the boundary of the
187 : * control volume. In computation the flux over each SCVF must be the same
188 : * in both directions over the face in order to guarantee the conservation
189 : * property.
190 : */
191 : class SCVF
192 : {
193 : public:
194 : /// Number of corners of scvf
195 : static const size_t numCo = dim;
196 :
197 : public:
198 0 : SCVF() {}
199 :
200 : /// index of SubControlVolume on one side of the scvf
201 0 : inline size_t from() const {return From;}
202 :
203 : /// index of SubControlVolume on one side of the scvf
204 0 : inline size_t to() const {return To;}
205 :
206 : /// number of integration points on scvf
207 0 : inline size_t num_ip() const {return nip;}
208 :
209 : /// local integration point of scvf
210 0 : inline const MathVector<dim>& local_ip() const {return localIP;}
211 :
212 : /// global integration point of scvf
213 0 : inline const MathVector<worldDim>& global_ip() const {return globalIP;}
214 :
215 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
216 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
217 :
218 : /// Transposed Inverse of Jacobian in integration point
219 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
220 :
221 : /// Determinant of Jacobian in integration point
222 0 : inline number detJ() const {return detj;}
223 :
224 : /// number of shape functions
225 0 : inline size_t num_sh() const {return nsh;}
226 :
227 : /// value of shape function i in integration point
228 0 : inline number shape(size_t sh) const {return vShape[sh];}
229 :
230 : /// vector of shape functions in ip point
231 0 : inline const number* shape_vector() const {return vShape;}
232 :
233 : /// value of local gradient of shape function i in integration point
234 0 : inline const MathVector<dim>& local_grad(size_t sh) const
235 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
236 :
237 : /// vector of local gradients in ip point
238 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
239 :
240 : /// value of global gradient of shape function i in integration point
241 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
242 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
243 :
244 : /// vector of global gradients in ip point
245 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
246 :
247 : /// number of corners, that bound the scvf
248 0 : inline size_t num_corners() const {return numCo;}
249 :
250 : /// return local corner number i
251 0 : inline const MathVector<dim>& local_corner(size_t co) const
252 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
253 :
254 : /// return glbal corner number i
255 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
256 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
257 :
258 : private:
259 : // let outer class access private members
260 : friend class HCRFVGeometry<TElem, TWorldDim>;
261 :
262 : // This scvf separates the scv with the ids given in "from" and "to"
263 : // The computed normal points in direction from->to
264 : size_t From, To;
265 :
266 : // The normal on the SCVF pointing (from -> to)
267 : MathVector<worldDim> Normal; // normal (incl. area)
268 :
269 : MathVector<dim> vLocPos[numCo]; // local corners of scvf
270 : MathVector<worldDim> vGloPos[numCo]; // global corners of scvf
271 :
272 : // scvf part
273 : MathVector<dim> localIP; // local integration point
274 : MathVector<worldDim> globalIP; // global integration point
275 :
276 : // shapes and derivatives
277 : number vShape[nsh]; // shapes at ip
278 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
279 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
280 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
281 : number detj; // Jacobian det at ip
282 : };
283 :
284 : /// sub control volume structure
285 : class SCV
286 : {
287 : public:
288 : /// Number of corners of scv
289 : static const size_t maxNumCo = 5;
290 :
291 : public:
292 0 : SCV() : numCorners(maxNumCo) {};
293 :
294 : /// volume of scv
295 0 : inline number volume() const {return Vol;}
296 :
297 : /// number of corners, that bound the scvf
298 0 : inline size_t num_corners() const {return numCorners;}
299 :
300 : /// return local corner number i
301 0 : inline const MathVector<dim>& local_corner(size_t co) const
302 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vLocPos[co];}
303 :
304 : /// return global corner number i
305 0 : inline const MathVector<worldDim>& global_corner(size_t co) const
306 0 : {UG_ASSERT(co < num_corners(), "Invalid index."); return vGloPos[co];}
307 :
308 : /// node id that this scv is associated to
309 0 : inline size_t node_id() const {return nodeID;}
310 :
311 : /// number of integration points
312 0 : inline size_t num_ip() const {return nip;}
313 :
314 : /// local integration point of scv
315 0 : inline const MathVector<dim>& local_ip() const {return vLocIP;}
316 :
317 : /// global integration point
318 0 : inline const MathVector<worldDim>& global_ip() const {return vGlobIP;}
319 :
320 : /// normal on scvf (points direction "from"->"to"). Norm is equal to area
321 0 : inline const MathVector<worldDim>& normal() const {return Normal;}
322 :
323 : /// Transposed Inverse of Jacobian in integration point
324 0 : inline const MathMatrix<worldDim,dim>& JTInv() const {return JtInv;}
325 :
326 : /// Determinant of Jacobian in integration point
327 0 : inline number detJ() const {return detj;}
328 :
329 : /// number of shape functions
330 0 : inline size_t num_sh() const {return nsh;}
331 :
332 : /// value of shape function i in integration point
333 0 : inline number shape(size_t sh) const {return vShape[sh];}
334 :
335 : /// vector of shape functions in ip point
336 0 : inline const number* shape_vector() const {return vShape;}
337 :
338 : /// value of local gradient of shape function i in integration point
339 0 : inline const MathVector<dim>& local_grad(size_t sh) const
340 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vLocalGrad[sh];}
341 :
342 : /// vector of local gradients in ip point
343 0 : inline const MathVector<dim>* local_grad_vector() const {return vLocalGrad;}
344 :
345 : /// value of global gradient of shape function i in integration point
346 0 : inline const MathVector<worldDim>& global_grad(size_t sh) const
347 0 : {UG_ASSERT(sh < num_sh(), "Invalid index"); return vGlobalGrad[sh];}
348 :
349 : /// vector of global gradients in ip point
350 0 : inline const MathVector<worldDim>* global_grad_vector() const {return vGlobalGrad;}
351 :
352 : private:
353 : // let outer class access private members
354 : friend class HCRFVGeometry<TElem, TWorldDim>;
355 :
356 : // The normal on the associated face pointing outward
357 : MathVector<worldDim> Normal; // normal (incl. area)
358 :
359 : // node id of associated node
360 : size_t nodeID;
361 :
362 : // volume of scv
363 : number Vol;
364 :
365 : // number of corners of this element
366 : int numCorners;
367 :
368 : // local and global positions of this element
369 : MathVector<dim> vLocPos[maxNumCo]; // local position of node
370 : MathVector<worldDim> vGloPos[maxNumCo]; // global position of node
371 :
372 : MathVector<dim> vLocIP; // local position of node
373 : MathVector<worldDim> vGlobIP; // global position of node
374 :
375 : // shapes and derivatives
376 : number vShape[nsh]; // shapes at ip
377 : MathVector<dim> vLocalGrad[nsh]; // local grad at ip
378 : MathVector<worldDim> vGlobalGrad[nsh]; // global grad at ip
379 : MathMatrix<worldDim,dim> JtInv; // Jacobian transposed at ip
380 : number detj; // Jacobian det at ip
381 : };
382 :
383 : class CONSTRAINED_DOF
384 : {
385 : public:
386 : static const size_t maxNumConstrainingDofs = 4;
387 0 : inline size_t constraining_dofs_index(size_t i) const{
388 0 : return cDofInd[i];
389 : }
390 0 : inline number constraining_dofs_weight(size_t i) const{
391 0 : return cDofWeights[i];
392 : }
393 0 : inline size_t index() const{
394 0 : return i;
395 : }
396 0 : inline size_t num_constraining_dofs() const{
397 0 : return numConstrainingDofs;
398 : }
399 : private:
400 : // let outer class access private members
401 : friend class HCRFVGeometry<TElem, TWorldDim>;
402 :
403 : // constraining dofs indices
404 : size_t cDofInd[maxNumConstrainingDofs];
405 : // weights
406 : number cDofWeights[maxNumConstrainingDofs];
407 : // local index of dof in element
408 : size_t i;
409 : // nr of constraining dofs
410 : size_t numConstrainingDofs;
411 : };
412 :
413 : class HandledEdge
414 : {
415 : public:
416 : size_t index;
417 : size_t associatedSCV[2];
418 : size_t scvfIndex;
419 : // indicates if the handled side is from or to
420 : bool from;
421 : };
422 :
423 : public:
424 : /// construct object and initialize local values and sizes
425 : HCRFVGeometry();
426 :
427 : /// update local data
428 : void update_local_data();
429 :
430 : /// update data for given element
431 : void update(GridObject* elem, const MathVector<worldDim>* vCornerCoords,
432 : const ISubsetHandler* ish = NULL);
433 :
434 : /// debug output
435 : void print();
436 :
437 0 : const MathVector<worldDim>* corners() const {return m_vCo;}
438 :
439 : /// number of SubControlVolumeFaces
440 0 : inline size_t num_scvf() const {return numSCVF;};
441 :
442 : /// const access to SubControlVolumeFace number i
443 0 : inline const SCVF& scvf(size_t i) const
444 0 : {UG_ASSERT(i < maxNumSCVF, "Invalid Index."); return m_vSCVF[i];}
445 :
446 : /// number of SubControlVolumes
447 0 : inline size_t num_scv() const {return numSCV;}
448 :
449 : /// const access to SubControlVolume number i
450 0 : inline const SCV& scv(size_t i) const
451 0 : {UG_ASSERT(i < maxNumSCV, "Invalid Index."); return m_vSCV[i];}
452 :
453 : /// number of constrained dofs
454 0 : inline size_t num_constrained_dofs() const {return numConstrainedDofs;}
455 :
456 : /// const access to constrained dof i
457 0 : inline const CONSTRAINED_DOF& constrained_dof(size_t i) const
458 0 : {UG_ASSERT(i < numConstrainedDofs, "Invalid Index."); return m_vCD[i];}
459 :
460 : /// number of shape functions
461 0 : inline size_t num_sh() const {return nsh;};
462 :
463 : public:
464 : /// returns number of all scvf ips
465 0 : size_t num_scvf_ips() const {return numSCVF;}
466 :
467 : /// returns all ips of scvf as they appear in scv loop
468 0 : const MathVector<dim>* scvf_local_ips() const {return m_vLocSCVF_IP;}
469 :
470 : /// returns all ips of scvf as they appear in scv loop
471 0 : const MathVector<worldDim>* scvf_global_ips() const {return m_vGlobSCVF_IP;}
472 :
473 : /// returns number of all scv ips
474 0 : size_t num_scv_ips() const {return numSCV;}
475 :
476 : /// returns all ips of scv as they appear in scv loop
477 0 : const MathVector<dim>* scv_local_ips() const {return m_vLocUnkCoords;}
478 :
479 : /// returns all ips of scv as they appear in scv loop
480 0 : const MathVector<worldDim>* scv_global_ips() const {return m_vGlobUnkCoords;}
481 :
482 : /// returns local barycenter
483 0 : const MathVector<dim> local_bary() const {return localBary;}
484 :
485 : /// returns global barycenter
486 0 : const MathVector<worldDim> global_bary() const {return globalBary;}
487 :
488 : protected:
489 : // global and local ips on SCVF
490 : MathVector<worldDim> m_vGlobSCVF_IP[maxNumSCVF];
491 : MathVector<dim> m_vLocSCVF_IP[maxNumSCVF];
492 : // coord of location for unknowns in faces (edge/face barycenter)
493 : MathVector<worldDim> m_vGlobUnkCoords[maxNumSCV];
494 : MathVector<dim> m_vLocUnkCoords[maxNumSCV];
495 :
496 : static const size_t numMaxCo = 8;
497 : // corner coordinates
498 : MathVector<worldDim> m_vCo[numMaxCo];
499 :
500 : private:
501 : MathVector<dim> m_ipCoord[maxNumSCVF];
502 :
503 : /// SubControlVolumeFaces
504 : SCVF m_vSCVF[maxNumSCVF];
505 :
506 : /// SubControlVolumes
507 : SCV m_vSCV[maxNumSCV];
508 :
509 : /// constrained Dofs
510 : CONSTRAINED_DOF m_vCD[maxNumSCV];
511 :
512 : std::vector<HandledEdge> handledEdges;
513 :
514 : /// pointer to current element
515 : TElem* m_pElem;
516 :
517 : /// Reference Mapping
518 : ReferenceMapping<ref_elem_type, worldDim> m_mapping;
519 :
520 : /// Reference Element
521 : const ref_elem_type& m_rRefElem;
522 :
523 : /// Shape function set
524 : const local_shape_fct_set_type& m_rTrialSpace;
525 :
526 : size_t numSCV;
527 : size_t numSCVF;
528 : size_t numConstrainedDofs;
529 : // numDofs number of all dofs including constraining and constrained dofs
530 : size_t numDofs;
531 :
532 : bool localUpdateNecessary;
533 :
534 : static const size_t deleted = 117;
535 :
536 : };
537 :
538 :
539 : }
540 :
541 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__HANGING_CR_FINITE_VOLUME_GEOMETRY__ */
|
