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__FINITE_ELEMENT_GEOMETRY__
34 : #define __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_ELEMENT_GEOMETRY__
35 :
36 : #include "lib_grid/tools/subset_handler_interface.h"
37 : #include "lib_disc/quadrature/quadrature.h"
38 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
39 : #include "lib_disc/reference_element/reference_mapping_provider.h"
40 : #include "lib_disc/reference_element/reference_mapping.h"
41 : #include "common/util/provider.h"
42 :
43 : #include <cmath>
44 :
45 : namespace ug{
46 :
47 : template < typename TElem, int TWorldDim,
48 : typename TTrialSpace, typename TQuadratureRule>
49 : class FEGeometry
50 : {
51 : public:
52 : /// type of reference element
53 : typedef typename reference_element_traits<TElem>::reference_element_type ref_elem_type;
54 :
55 : /// reference element dimension
56 : static const int dim = ref_elem_type::dim;
57 :
58 : /// world dimension
59 : static const int worldDim = TWorldDim;
60 :
61 : /// type of trial space
62 : typedef TTrialSpace trial_space_type;
63 :
64 : /// type of quadrature rule
65 : typedef TQuadratureRule quad_rule_type;
66 :
67 : /// number of shape functions
68 : static const size_t nsh = trial_space_type::nsh;
69 :
70 : /// number of integration points
71 : static const size_t nip = quad_rule_type::nip;
72 :
73 : /// flag indicating if local data may change
74 : static const bool staticLocalData = true;
75 :
76 : public:
77 : /// Constructor
78 : FEGeometry();
79 :
80 : /// number of integration points
81 : size_t num_ip() const {return nip;}
82 :
83 : /// number of shape functions
84 : size_t num_sh() const {return nsh;}
85 :
86 : /// weight for integration point
87 0 : number weight(size_t ip) const {return m_vDetJ[ip] * m_rQuadRule.weight(ip);}
88 :
89 : /// local integration point
90 : const MathVector<dim>& local_ip(size_t ip) const {return m_rQuadRule.point(ip);}
91 :
92 : /// global integration point
93 : const MathVector<worldDim>& global_ip(size_t ip) const
94 : {
95 : UG_ASSERT(ip < nip, "Wrong ip.");
96 : return m_vIPGlobal[ip];
97 : }
98 :
99 : /// local integration point
100 : const MathVector<dim>* local_ips() const {return m_rQuadRule.points();}
101 :
102 : /// global integration point
103 0 : const MathVector<worldDim>* global_ips() const{return &m_vIPGlobal[0];}
104 :
105 : /// shape function at ip
106 : number shape(size_t ip, size_t sh) const
107 : {
108 : UG_ASSERT(ip < nip, "Wrong index"); UG_ASSERT(sh < nsh, "Wrong index");
109 0 : return m_vvShape[ip][sh];
110 : }
111 :
112 : /// local gradient at ip
113 : const MathVector<dim>& local_grad(size_t ip, size_t sh) const
114 : {
115 : UG_ASSERT(ip < nip, "Wrong index"); UG_ASSERT(sh < nsh, "Wrong index");
116 : return m_vvGradLocal[ip][sh];
117 : }
118 :
119 : /// global gradient at ip
120 : const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
121 : {
122 : UG_ASSERT(ip < nip, "Wrong index"); UG_ASSERT(sh < nsh, "Wrong index");
123 : return m_vvGradGlobal[ip][sh];
124 : }
125 :
126 : public:
127 : /// update Geometry for roid
128 : void update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad);
129 : void update_local(ReferenceObjectID roid, const LFEID& lfeID){
130 : update_local(roid, lfeID, 2*lfeID.order() + 1);
131 : }
132 :
133 : /// update Geometry for corners
134 : void update(GridObject* pElem, const MathVector<worldDim>* vCorner)
135 : {
136 0 : update(pElem, vCorner, LFEID(), m_rQuadRule.order());
137 : }
138 :
139 : /// update Geometry for corners
140 : void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
141 : const LFEID& lfeID, size_t orderQuad);
142 : void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
143 : const LFEID& lfeID){
144 : update(pElem, vCorner, lfeID, 2*lfeID.order() + 1);
145 : }
146 :
147 : protected:
148 : /// current element
149 : TElem* m_pElem;
150 :
151 : /// Quadrature rule
152 : const quad_rule_type& m_rQuadRule;
153 :
154 : /// Quadrature rule
155 : const trial_space_type& m_rTrialSpace;
156 :
157 : /// reference mapping
158 : ReferenceMapping<ref_elem_type, worldDim> m_mapping;
159 :
160 : /// global integration points
161 : MathVector<worldDim> m_vIPGlobal[nip];
162 :
163 : /// shape functions evaluated at ip
164 : number m_vvShape[nip][nsh];
165 :
166 : /// global gradient evaluated at ip
167 : MathVector<dim> m_vvGradLocal[nip][nsh];
168 :
169 : /// local gradient evaluated at ip
170 : MathVector<worldDim> m_vvGradGlobal[nip][nsh];
171 :
172 : /// jacobian of transformation at ip
173 : MathMatrix<worldDim,dim> m_vJTInv[nip];
174 :
175 : /// determinate of transformation at ip
176 : number m_vDetJ[nip];
177 : };
178 :
179 :
180 :
181 : template <int TWorldDim, int TRefDim = TWorldDim>
182 : class DimFEGeometry
183 : {
184 : public:
185 : /// reference element dimension
186 : static const int dim = TRefDim;
187 :
188 : /// world dimension
189 : static const int worldDim = TWorldDim;
190 :
191 : /// flag indicating if local data may change
192 : static const bool staticLocalData = false;
193 :
194 : public:
195 : /// default Constructor
196 : DimFEGeometry();
197 :
198 : /// Constructor
199 : DimFEGeometry(size_t order, LFEID lfeid);
200 :
201 : /// Constructor
202 : DimFEGeometry(ReferenceObjectID roid, size_t order, LFEID lfeid);
203 :
204 : /// number of integration points
205 0 : size_t num_ip() const {return m_nip;}
206 :
207 : /// number of shape functions
208 0 : size_t num_sh() const {return m_nsh;}
209 :
210 : /// weight for integration point
211 0 : number weight(size_t ip) const
212 : {
213 : UG_ASSERT(ip < m_nip, "Wrong index");
214 0 : return m_vDetJ[ip] * m_vQuadWeight[ip];
215 : }
216 :
217 : /// local integration point
218 0 : const MathVector<dim>& local_ip(size_t ip) const
219 : {
220 : UG_ASSERT(ip < m_nip, "Wrong index");
221 0 : return m_vIPLocal[ip];
222 : }
223 :
224 : /// global integration point
225 0 : const MathVector<worldDim>& global_ip(size_t ip) const
226 : {
227 : UG_ASSERT(ip < m_vIPGlobal.size(), "Wrong ip.");
228 0 : return m_vIPGlobal[ip];
229 : }
230 :
231 : /// local integration point
232 0 : const MathVector<dim>* local_ips() const {return m_vIPLocal;}
233 :
234 : /// global integration point
235 0 : const MathVector<worldDim>* global_ips() const{return &m_vIPGlobal[0];}
236 :
237 : /// shape function at ip
238 0 : number shape(size_t ip, size_t sh) const
239 : {
240 : UG_ASSERT(ip < m_vvShape.size(), "Wrong index");
241 : UG_ASSERT(sh < m_vvShape[ip].size(), "Wrong index");
242 0 : return m_vvShape[ip][sh];
243 : }
244 :
245 : /// local gradient at ip
246 0 : const MathVector<dim>& local_grad(size_t ip, size_t sh) const
247 : {
248 : UG_ASSERT(ip < m_vvGradLocal.size(), "Wrong index");
249 : UG_ASSERT(sh < m_vvGradLocal[ip].size(), "Wrong index");
250 0 : return m_vvGradLocal[ip][sh];
251 : }
252 :
253 : /// global gradient at ip
254 0 : const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
255 : {
256 : UG_ASSERT(ip < m_vvGradGlobal.size(), "Wrong index");
257 : UG_ASSERT(sh < m_vvGradGlobal[ip].size(), "Wrong index");
258 0 : return m_vvGradGlobal[ip][sh];
259 : }
260 :
261 : /// update Geometry for roid
262 : void update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad);
263 0 : void update_local(ReferenceObjectID roid, const LFEID& lfeID){
264 0 : update_local(roid, lfeID, 2*lfeID.order() + 1);
265 0 : }
266 :
267 : /// update Geometry for corners
268 0 : void update(GridObject* pElem, const MathVector<worldDim>* vCorner)
269 : {
270 0 : update(pElem, vCorner, m_lfeID, m_quadOrder);
271 0 : }
272 :
273 : /// update Geometry for corners
274 : void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
275 : const LFEID& lfeID, size_t orderQuad);
276 0 : void update(GridObject* pElem, const MathVector<worldDim>* vCorner,
277 : const LFEID& lfeID){
278 0 : update(pElem, vCorner, lfeID, 2*lfeID.order() + 1);
279 0 : }
280 :
281 : public:
282 : /// boundary face
283 : class BF
284 : {
285 : public:
286 0 : BF() {}
287 :
288 : /// outer normal on bf. Norm is equal to area
289 0 : inline const MathVector<worldDim>& normal() const {return Normal;} // includes area
290 :
291 : /// volume of bf
292 0 : inline number volume() const {return Vol;}
293 :
294 : /// number of integration points on scvf
295 0 : inline size_t num_ip() const {return nip;}
296 :
297 : /// integration weight
298 0 : inline number weight(size_t ip) const
299 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip] * vWeight[ip];}
300 :
301 : /// local integration point of scvf
302 0 : inline const MathVector<dim>& local_ip(size_t ip) const
303 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vLocalIP[ip];}
304 :
305 : /// global integration point of scvf
306 0 : inline const MathVector<worldDim>& global_ip(size_t ip) const
307 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vGlobalIP[ip];}
308 :
309 : /// Transposed Inverse of Jacobian in integration point
310 0 : inline const MathMatrix<worldDim,dim>& JTInv(size_t ip) const
311 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vJtInv[ip];}
312 :
313 : /// Determinant of Jacobian in integration point
314 0 : inline number detJ(size_t ip) const
315 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vDetJ[ip];}
316 :
317 : /// number of shape functions
318 0 : inline size_t num_sh() const {return nsh;}
319 :
320 : /// value of shape function i in integration point
321 0 : inline number shape(size_t ip, size_t sh) const
322 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return vvShape[ip][sh];}
323 :
324 : /// vector of shape functions in ip point
325 0 : inline const number* shape_vector(size_t ip) const
326 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvShape[ip][0];}
327 :
328 : /// value of local gradient of shape function i in integration point
329 0 : inline const MathVector<dim>& local_grad(size_t ip, size_t sh) const
330 : {UG_ASSERT(sh < num_sh(), "Invalid index");
331 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvLocalGrad[ip][sh];}
332 :
333 : /// vector of local gradients in ip point
334 0 : inline const MathVector<dim>* local_grad_vector(size_t ip) const
335 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvLocalGrad[ip][0];}
336 :
337 : /// value of global gradient of shape function i in integration point
338 0 : inline const MathVector<worldDim>& global_grad(size_t ip, size_t sh) const
339 : {UG_ASSERT(sh < num_sh(), "Invalid index");
340 0 : UG_ASSERT(ip<num_ip(), "Wrong index"); return vvGlobalGrad[ip][sh];}
341 :
342 : /// vector of global gradients in ip point
343 0 : inline const MathVector<worldDim>* global_grad_vector(size_t ip) const
344 0 : {UG_ASSERT(ip<num_ip(), "Wrong index"); return &vvGlobalGrad[ip][0];}
345 :
346 : protected:
347 : /// let outer class access private members
348 : friend class DimFEGeometry<worldDim, dim>;
349 :
350 : size_t nip;
351 : std::vector<MathVector<dim> > vLocalIP; // local integration point (size: nip)
352 : std::vector<MathVector<worldDim> > vGlobalIP; // global integration point (size: nip)
353 : const number* vWeight; // weight at ip
354 : MathVector<worldDim> Normal; // normal (incl. area)
355 : number Vol; // volume of bf
356 : std::vector<MathMatrix<worldDim,dim> > vJtInv; // Jacobian transposed at ip (size: nip)
357 : std::vector<number> vDetJ; // Jacobian det at ip (size: nip)
358 :
359 : size_t nsh;
360 : std::vector<std::vector<number> > vvShape; // shapes at ip (size: nip x nsh)
361 : std::vector<std::vector<MathVector<dim> > > vvLocalGrad; // local grad at ip (size: nip x nsh)
362 : std::vector<std::vector<MathVector<worldDim> > > vvGlobalGrad; // global grad at ip (size: nip x nsh)
363 : };
364 :
365 : /// update boundary data for given element
366 : void update_boundary_faces(GridObject* pElem,
367 : const MathVector<worldDim>* vCornerCoords,
368 : size_t orderQuad,
369 : const ISubsetHandler* ish);
370 :
371 : public:
372 : /// add subset that is interpreted as boundary subset.
373 0 : inline void add_boundary_subset(int subsetIndex) {m_mapVectorBF[subsetIndex];}
374 :
375 : /// removes subset that is interpreted as boundary subset.
376 0 : inline void remove_boundary_subset(int subsetIndex) {m_mapVectorBF.erase(subsetIndex);}
377 :
378 : /// reset all boundary subsets
379 0 : inline void clear_boundary_subsets() {m_mapVectorBF.clear();}
380 :
381 : /// number of registered boundary subsets
382 0 : inline size_t num_boundary_subsets() {return m_mapVectorBF.size();}
383 :
384 : /// number of all boundary faces
385 0 : inline size_t num_bf() const
386 : {
387 : typename std::map<int, std::vector<BF> >::const_iterator it;
388 : size_t num = 0;
389 0 : for ( it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); it++ )
390 0 : num += (*it).second.size();
391 0 : return num;
392 : }
393 :
394 : /// number of boundary faces on subset 'subsetIndex'
395 0 : inline size_t num_bf(int si) const
396 : {
397 : typename std::map<int, std::vector<BF> >::const_iterator it;
398 : it = m_mapVectorBF.find(si);
399 0 : if(it == m_mapVectorBF.end()) return 0;
400 0 : else return (*it).second.size();
401 : }
402 :
403 : /// returns the boundary face i for subsetIndex
404 0 : inline const BF& bf(int si, size_t i) const
405 : {
406 : UG_ASSERT(i < num_bf(si), "Invalid index.");
407 : typename std::map<int, std::vector<BF> >::const_iterator it;
408 : it = m_mapVectorBF.find(si);
409 0 : if(it == m_mapVectorBF.end()) UG_THROW("DimFVGeom: No BndSubset: "<<si);
410 0 : return (*it).second[i];
411 : }
412 :
413 : /// returns reference to vector of boundary faces for subsetIndex
414 0 : inline const std::vector<BF>& bf(int si) const
415 : {
416 : typename std::map<int, std::vector<BF> >::const_iterator it;
417 : it = m_mapVectorBF.find(si);
418 0 : if(it == m_mapVectorBF.end()) return m_vEmptyVectorBF;
419 0 : return (*it).second;
420 : }
421 :
422 : protected:
423 : std::map<int, std::vector<BF> > m_mapVectorBF;
424 : std::vector<BF> m_vEmptyVectorBF;
425 :
426 : protected:
427 : /// current reference object id the local values are prepared for
428 : ReferenceObjectID m_roid;
429 :
430 : /// current element
431 : GridObject* m_pElem;
432 :
433 : /// current integration order
434 : int m_quadOrder;
435 :
436 : /// current local finite element id
437 : LFEID m_lfeID;
438 :
439 : /// number of integration point
440 : size_t m_nip;
441 :
442 : /// local quadrature points
443 : const MathVector<dim>* m_vIPLocal;
444 :
445 : /// local quadrature weights
446 : const number* m_vQuadWeight;
447 :
448 : /// global integration points (size = nip)
449 : std::vector<MathVector<worldDim> > m_vIPGlobal;
450 :
451 : /// jacobian of transformation at ip (size = nip)
452 : std::vector<MathMatrix<worldDim,dim> > m_vJTInv;
453 :
454 : /// determinant of transformation at ip (size = nip)
455 : std::vector<number> m_vDetJ;
456 :
457 : /// number of shape functions
458 : size_t m_nsh;
459 :
460 : /// shape functions evaluated at ip (size = nip x nsh)
461 : std::vector<std::vector<number> > m_vvShape;
462 :
463 : /// global gradient evaluated at ip (size = nip x nsh)
464 : std::vector<std::vector<MathVector<dim> > > m_vvGradLocal;
465 :
466 : /// local gradient evaluated at ip (size = nip x nsh)
467 : std::vector<std::vector<MathVector<worldDim> > > m_vvGradGlobal;
468 : };
469 :
470 : } // end namespace ug
471 :
472 : #include "fe_geom_impl.h"
473 :
474 : #endif /* __H__UG__LIB_DISC__SPATIAL_DISC__DISC_HELPER__FINITE_ELEMENT_GEOMETRY__ */
|
