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__LOCAL_SHAPE_FUNCTION_SET__LAGRANGE__LAGRANGE__
34 : #define __H__UG__LIB_DISC__LOCAL_SHAPE_FUNCTION_SET__LAGRANGE__LAGRANGE__
35 :
36 : #include "../common/lagrange1d.h"
37 : #include "../local_finite_element_provider.h"
38 : #include "../local_dof_set.h"
39 : #include "lagrange_local_dof.h"
40 : #include "lib_disc/common/multi_index.h"
41 : #include "common/util/metaprogramming_util.h"
42 : #include "lib_grid/grid/grid_base_objects.h"
43 :
44 : namespace ug{
45 :
46 : /// Lagrange Shape Function Set without virtual functions and fixed order
47 : template <typename TRefElem, int TOrder>
48 : class LagrangeLSFS;
49 :
50 : /// Lagrange Shape Function Set without virtual functions and flexible order
51 : template <typename TRefElem>
52 : class FlexLagrangeLSFS;
53 :
54 :
55 : ///////////////////////////////////////////////////////////////////////////////
56 : // Vertex
57 : ///////////////////////////////////////////////////////////////////////////////
58 :
59 : /// specialization for Vertex
60 : /**
61 : * Lagrange shape function of any order for the Reference Vertex
62 : * \tparam TOrder requested order
63 : */
64 : template <int TOrder>
65 : class LagrangeLSFS<ReferenceVertex, TOrder>
66 : : public LagrangeLDS<ReferenceVertex>,
67 : public BaseLSFS<LagrangeLSFS<ReferenceVertex, TOrder>, 0>
68 : {
69 : private:
70 : /// abbreviation for order
71 : static const size_t p = TOrder;
72 :
73 : /// base class
74 : typedef BaseLSFS<LagrangeLSFS<ReferenceVertex, TOrder>, 0> base_type;
75 :
76 : public:
77 : /// Shape type
78 : typedef typename base_type::shape_type shape_type;
79 :
80 : /// Gradient type
81 : typedef typename base_type::grad_type grad_type;
82 :
83 : public:
84 : /// Order of Shape functions
85 : static const size_t order = TOrder;
86 :
87 : /// Dimension, where shape functions are defined
88 : static const int dim = ReferenceVertex::dim;
89 :
90 : /// Number of shape functions
91 : static const size_t nsh = 1;
92 :
93 : public:
94 : /// Constructor
95 0 : LagrangeLSFS() {}
96 :
97 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
98 : inline bool continuous() const {return true;}
99 :
100 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
101 0 : inline size_t num_sh() const {return nsh;}
102 :
103 : /// \copydoc ug::LocalShapeFunctionSet::position()
104 : inline bool position(size_t i, MathVector<dim>& pos) const
105 : {
106 : return true;
107 : }
108 :
109 : /// \copydoc ug::LocalShapeFunctionSet::shape()
110 : inline shape_type shape(size_t i, const MathVector<dim>& x) const
111 : {
112 : return 1.0;
113 : }
114 :
115 : /// \copydoc ug::LocalShapeFunctionSet::grad()
116 : inline void grad(grad_type& g, size_t i, const MathVector<dim>& x) const
117 : {
118 : }
119 : };
120 :
121 : /// specialization for Edges
122 : /**
123 : * Lagrange shape function of any order for the Reference Edge
124 : */
125 : template <>
126 : class FlexLagrangeLSFS<ReferenceVertex>
127 : : public LagrangeLDS<ReferenceVertex>,
128 : public BaseLSFS<FlexLagrangeLSFS<ReferenceVertex>, 0>
129 : {
130 : public:
131 : /// Dimension, where shape functions are defined
132 : static const int dim = ReferenceVertex::dim;
133 :
134 : public:
135 : /// default Constructor
136 0 : FlexLagrangeLSFS() {}
137 :
138 : /// Constructor
139 : FlexLagrangeLSFS(size_t order) {p = order;}
140 :
141 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
142 : inline bool continuous() const {return true;}
143 :
144 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
145 0 : inline size_t num_sh() const {return nsh;}
146 :
147 : /// \copydoc ug::LocalShapeFunctionSet::position()
148 : inline bool position(size_t i, MathVector<dim>& pos) const
149 : {
150 : return true;
151 : }
152 :
153 : /// \copydoc ug::LocalShapeFunctionSet::shape()
154 : inline shape_type shape(size_t i, const MathVector<dim>& x) const
155 : {
156 : return 1.0;
157 : }
158 :
159 : /// \copydoc ug::LocalShapeFunctionSet::grad()
160 : inline void grad(grad_type& g, size_t i, const MathVector<dim>& x) const
161 : {
162 : }
163 :
164 : protected:
165 : /// order
166 : size_t p;
167 :
168 : /// Number of shape functions
169 : static const size_t nsh = 1;
170 : };
171 :
172 : ///////////////////////////////////////////////////////////////////////////////
173 : // Edge
174 : ///////////////////////////////////////////////////////////////////////////////
175 :
176 : /// specialization for Edges
177 : /**
178 : * Lagrange shape function of any order for the Reference Edge
179 : * \tparam TOrder requested order
180 : */
181 : template <int TOrder>
182 : class LagrangeLSFS<ReferenceEdge, TOrder>
183 : : public LagrangeLDS<ReferenceEdge>,
184 : public BaseLSFS<LagrangeLSFS<ReferenceEdge, TOrder>, 1>
185 : {
186 : private:
187 : /// abbreviation for order
188 : static const size_t p = TOrder;
189 :
190 : /// base class
191 : typedef BaseLSFS<LagrangeLSFS<ReferenceEdge, TOrder>, 1> base_type;
192 :
193 : public:
194 : /// Shape type
195 : typedef typename base_type::shape_type shape_type;
196 :
197 : /// Gradient type
198 : typedef typename base_type::grad_type grad_type;
199 :
200 : public:
201 : /// Order of Shape functions
202 : static const size_t order = TOrder;
203 :
204 : /// Dimension, where shape functions are defined
205 : static const int dim = ReferenceEdge::dim;
206 :
207 : /// Number of shape functions
208 : static const size_t nsh = p+1;
209 :
210 : public:
211 : /// Constructor
212 : LagrangeLSFS();
213 :
214 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
215 0 : inline bool continuous() const {return true;}
216 :
217 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
218 720 : inline size_t num_sh() const {return nsh;}
219 :
220 : /// \copydoc ug::LocalShapeFunctionSet::position()
221 0 : inline bool position(size_t i, MathVector<dim>& pos) const
222 : {
223 27 : pos = EquidistantLagrange1D::position(multi_index(i)[0], p);
224 0 : return true;
225 : }
226 :
227 : /// \copydoc ug::LocalShapeFunctionSet::shape()
228 0 : inline shape_type shape(size_t i, const MathVector<dim>& x) const
229 : {
230 778 : return m_vPolynom[multi_index(i)[0]].value(x[0]);
231 : }
232 :
233 : /// \copydoc ug::LocalShapeFunctionSet::grad()
234 0 : inline void grad(grad_type& g, size_t i, const MathVector<dim>& x) const
235 : {
236 540 : g[0] = m_vDPolynom[multi_index(i)[0]].value(x[0]);
237 0 : }
238 :
239 : /// return Multi index for index i
240 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
241 : {
242 : check_index(i);
243 28 : return m_vMultiIndex[i];
244 : }
245 :
246 : /// return the index for a multi_index
247 9 : inline size_t index(const MathVector<dim,int>& ind) const
248 : {
249 : check_multi_index(ind);
250 19 : for(size_t i=0; i<nsh; ++i)
251 19 : if(multi_index(i) == ind) return i;
252 0 : UG_THROW("Index not found in LagrangeLSFS");
253 : }
254 :
255 : /// return Multi index for index i
256 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
257 : {
258 : check_index(i);
259 0 : return MathVector<1,int>(i);
260 : }
261 :
262 : /// return the index for a multi_index
263 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
264 : {
265 : check_multi_index(ind);
266 0 : return ind[0];
267 : }
268 :
269 : /// checks in debug mode that index is valid
270 0 : inline void check_index(size_t i) const
271 : {
272 : UG_ASSERT(i < nsh, "Wrong index.");
273 0 : }
274 :
275 : /// checks in debug mode that multi-index is valid
276 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
277 : {
278 : UG_ASSERT(ind[0] < (int)nsh && ind[0] >= 0, "Wrong MultiIndex");
279 0 : }
280 :
281 : protected:
282 : Polynomial1D m_vPolynom[p+1]; ///< Shape Polynomials
283 : Polynomial1D m_vDPolynom[p+1]; ///< Derivative of Shape Polynomial
284 :
285 : MathVector<dim,int> m_vMultiIndex[nsh];
286 : };
287 :
288 : /// specialization for Edges
289 : /**
290 : * Lagrange shape function of any order for the Reference Edge
291 : */
292 : template <>
293 : class FlexLagrangeLSFS<ReferenceEdge>
294 : : public LagrangeLDS<ReferenceEdge>,
295 : public BaseLSFS<FlexLagrangeLSFS<ReferenceEdge>, 1>
296 : {
297 : public:
298 : /// Dimension, where shape functions are defined
299 : static const int dim = ReferenceEdge::dim;
300 :
301 : public:
302 : /// default Constructor
303 0 : FlexLagrangeLSFS() {set_order(1);}
304 :
305 : /// Constructor
306 0 : FlexLagrangeLSFS(size_t order) {set_order(order);}
307 :
308 : /// sets the order
309 : void set_order(size_t order);
310 :
311 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
312 : inline bool continuous() const {return true;}
313 :
314 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
315 0 : inline size_t num_sh() const {return nsh;}
316 :
317 : /// \copydoc ug::LocalShapeFunctionSet::position()
318 : inline bool position(size_t i, MathVector<dim>& pos) const
319 : {
320 0 : pos = EquidistantLagrange1D::position(multi_index(i)[0], p);
321 : return true;
322 : }
323 :
324 : /// \copydoc ug::LocalShapeFunctionSet::shape()
325 0 : inline shape_type shape(size_t i, const MathVector<dim>& x) const
326 : {
327 0 : return m_vPolynom[multi_index(i)[0]].value(x[0]);
328 : }
329 :
330 : /// \copydoc ug::LocalShapeFunctionSet::grad()
331 0 : inline void grad(grad_type& g, size_t i, const MathVector<dim>& x) const
332 : {
333 0 : g[0] = m_vDPolynom[multi_index(i)[0]].value(x[0]);
334 0 : }
335 :
336 : /// return Multi index for index i
337 : inline const MathVector<dim,int>& multi_index(size_t i) const
338 : {
339 : check_index(i);
340 : return m_vMultiIndex[i];
341 : }
342 :
343 : /// return the index for a multi_index
344 0 : inline size_t index(const MathVector<dim,int>& ind) const
345 : {
346 : check_multi_index(ind);
347 0 : for(size_t i=0; i<nsh; ++i)
348 0 : if(multi_index(i) == ind) return i;
349 0 : UG_THROW("Index not found in LagrangeLSFS");
350 : }
351 :
352 : /// return Multi index for index i
353 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
354 : {
355 : check_index(i);
356 : return MathVector<1,int>(i);
357 : }
358 :
359 : /// return the index for a multi_index
360 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
361 : {
362 : check_multi_index(ind);
363 : return ind[0];
364 : }
365 :
366 : /// checks in debug mode that index is valid
367 : inline void check_index(size_t i) const
368 : {
369 : UG_ASSERT(i < nsh, "Wrong index.");
370 : }
371 :
372 : /// checks in debug mode that multi-index is valid
373 : inline void check_multi_index(const MathVector<dim,int>& ind) const
374 : {
375 : UG_ASSERT(ind[0] < (int)nsh && ind[0] >= 0, "Wrong MultiIndex");
376 : }
377 :
378 : protected:
379 : /// order
380 : size_t p;
381 :
382 : /// Number of shape functions
383 : size_t nsh;
384 :
385 : std::vector<Polynomial1D> m_vPolynom; ///< Shape Polynomials
386 : std::vector<Polynomial1D> m_vDPolynom; ///< Derivative of Shape Polynomial
387 :
388 : std::vector<MathVector<dim,int> > m_vMultiIndex;
389 : };
390 :
391 : ///////////////////////////////////////////////////////////////////////////////
392 : // Triangle
393 : ///////////////////////////////////////////////////////////////////////////////
394 :
395 : template <int TOrder>
396 : class LagrangeLSFS<ReferenceTriangle, TOrder>
397 : : public LagrangeLDS<ReferenceTriangle>,
398 : public BaseLSFS<LagrangeLSFS<ReferenceTriangle, TOrder>, 2>
399 : {
400 : private:
401 : /// abbreviation for order
402 : static const size_t p = TOrder;
403 :
404 : /// base class
405 : typedef BaseLSFS<LagrangeLSFS<ReferenceTriangle, TOrder>, 2> base_type;
406 :
407 : public:
408 : /// Shape type
409 : typedef typename base_type::shape_type shape_type;
410 :
411 : /// Gradient type
412 : typedef typename base_type::grad_type grad_type;
413 :
414 : public:
415 : /// Order of Shape functions
416 : static const size_t order = TOrder;
417 :
418 : /// Dimension, where shape functions are defined
419 : static const int dim = ReferenceTriangle::dim;
420 :
421 : /// Number of shape functions
422 : static const size_t nsh = BinomialCoefficient<dim + p, p>::value;
423 :
424 : public:
425 : /// Constructor
426 : LagrangeLSFS();
427 :
428 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
429 0 : inline bool continuous() const {return true;}
430 :
431 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
432 1320 : inline size_t num_sh() const {return nsh;}
433 :
434 : /// \copydoc ug::LocalShapeFunctionSet::position()
435 0 : inline bool position(size_t i, MathVector<dim>& pos) const
436 : {
437 : // get Multi Index
438 : MathVector<dim,int> ind = multi_index(i);
439 :
440 : // set position
441 141 : for(int d = 0; d < dim; ++d)
442 94 : pos[d] = TruncatedEquidistantLagrange1D::position(ind[d], p);
443 :
444 0 : return true;
445 : }
446 :
447 : /// \copydoc ug::LocalShapeFunctionSet::shape()
448 0 : inline number shape(const size_t i, const MathVector<dim>& x) const
449 : {
450 : // forward
451 1810 : return shape(multi_index(i), x);
452 : }
453 :
454 : /// shape value for a Multi Index
455 1810 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
456 : {
457 : check_multi_index(ind);
458 : //ReferenceTriangle::check_position(x);
459 :
460 : // get adjoint barycentric index
461 1810 : const size_t i0 = p - ind[0] - ind[1];
462 1810 : const number x0 = 1.0 - x[0] - x[1];
463 :
464 1810 : return m_vPolynom[ ind[0] ].value(x[0])
465 1810 : * m_vPolynom[ ind[1] ].value(x[1])
466 1810 : * m_vPolynom[ i0 ].value(x0);
467 : }
468 :
469 :
470 : /// \copydoc ug::LocalShapeFunctionSet::shape()
471 0 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
472 : {
473 760 : grad(g, multi_index(i), x);
474 0 : }
475 :
476 : /// evaluates the gradient
477 760 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
478 : const MathVector<dim>& x) const
479 : {
480 : check_multi_index(ind);
481 : //ReferenceTriangle::check_position(x);
482 :
483 : // get adjoint barycentric index and position
484 760 : const size_t i0 = p - ind[0] - ind[1];
485 760 : const number x0 = 1.0 - x[0] - x[1];
486 :
487 : UG_ASSERT(i0 <= (int)p, "Wrong Multiindex.");
488 : UG_ASSERT(x0 <= 1.0+SMALL && x0 >= -SMALL, "Wrong Position.");
489 :
490 : // loop dimensions
491 2280 : for(int d = 0; d < dim; ++d)
492 : {
493 1520 : g[d] = m_vDPolynom[ind[d]].value(x[d])
494 1520 : * m_vPolynom[i0].value(x0);
495 1520 : g[d] += (-1) * m_vDPolynom[i0].value(x0)
496 1520 : * m_vPolynom[ind[d]].value(x[d]);
497 :
498 : // multiply by all functions not depending on x[d]
499 4560 : for(int d2 = 0; d2 < dim; ++d2)
500 : {
501 : // skip own value
502 3040 : if(d2 == d) continue;
503 :
504 3040 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
505 : }
506 : }
507 760 : }
508 :
509 : /// return Multi index for index i
510 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
511 : {
512 : check_index(i);
513 1819 : return m_vMultiIndex[i];
514 : }
515 :
516 : /// return the index for a multi_index
517 19 : inline size_t index(const MathVector<dim,int>& ind) const
518 : {
519 : check_multi_index(ind);
520 82 : for(size_t i=0; i<nsh; ++i)
521 82 : if(multi_index(i) == ind) return i;
522 0 : UG_THROW("Index not found in LagrangeLSFS");
523 : }
524 :
525 : /// return the index for a multi_index
526 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
527 : {
528 : check_multi_index(ind);
529 :
530 0 : size_t res = ind[0];
531 0 : for(int i = 0; i < ind[1]; ++i)
532 0 : res += (p+1-i);
533 :
534 : check_index(res);
535 0 : return res;
536 : }
537 :
538 : /// return the multi_index for an index
539 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
540 : {
541 : check_index(i);
542 :
543 0 : int i0 = i, i1;
544 0 : for(i1 = 0; i1 < (int)p; ++i1)
545 : {
546 0 : const int diff = i0 - (p+1-i1);
547 0 : if(diff < 0) break;
548 : i0 = diff;
549 : }
550 :
551 : UG_ASSERT(i0 >= 0, "i0 is negative ("<<i0<<")");
552 : UG_ASSERT(i1 >= 0, "i1 is negative ("<<i1<<")");
553 0 : return MathVector<dim,int>( i0, i1 );
554 : }
555 :
556 : /// checks in debug mode that index is valid
557 0 : inline void check_index(size_t i) const
558 : {
559 : UG_ASSERT(i < nsh, "Wrong index.");
560 0 : }
561 :
562 : /// checks in debug mode that multi-index is valid
563 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
564 : {
565 : UG_ASSERT(ind[0] <= (int)p && ind[0]>=0, "Wrong Multiindex.");
566 : UG_ASSERT(ind[1] <= (int)p && ind[0]>=0, "Wrong Multiindex.");
567 : UG_ASSERT(ind[0] + ind[1] <= (int)p && ind[0]>=0, "Wrong Multiindex.");
568 0 : }
569 :
570 : private:
571 : Polynomial1D m_vPolynom[p+1];
572 : Polynomial1D m_vDPolynom[p+1];
573 :
574 : MathVector<dim,int> m_vMultiIndex[nsh];
575 : };
576 :
577 :
578 : template <>
579 : class FlexLagrangeLSFS<ReferenceTriangle>
580 : : public LagrangeLDS<ReferenceTriangle>,
581 : public BaseLSFS<FlexLagrangeLSFS<ReferenceTriangle>, 2>
582 : {
583 : public:
584 : /// Dimension, where shape functions are defined
585 : static const int dim = ReferenceTriangle::dim;
586 :
587 : public:
588 : /// default Constructor
589 0 : FlexLagrangeLSFS() {set_order(1);}
590 :
591 : /// Constructor
592 0 : FlexLagrangeLSFS(size_t order) {set_order(order);}
593 :
594 : /// sets the order
595 : void set_order(size_t order);
596 :
597 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
598 : inline bool continuous() const {return true;}
599 :
600 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
601 0 : inline size_t num_sh() const {return nsh;}
602 :
603 : /// \copydoc ug::LocalShapeFunctionSet::position()
604 : inline bool position(size_t i, MathVector<dim>& pos) const
605 : {
606 : // get Multi Index
607 : MathVector<dim,int> ind = multi_index(i);
608 :
609 : // set position
610 0 : for(int d = 0; d < dim; ++d)
611 0 : pos[d] = TruncatedEquidistantLagrange1D::position(ind[d], p);
612 :
613 : return true;
614 : }
615 :
616 : /// \copydoc ug::LocalShapeFunctionSet::shape()
617 : inline number shape(const size_t i, const MathVector<dim>& x) const
618 : {
619 : // forward
620 0 : return shape(multi_index(i), x);
621 : }
622 :
623 : /// shape value for a Multi Index
624 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
625 : {
626 : check_multi_index(ind);
627 : //ReferenceTriangle::check_position(x);
628 :
629 : // get adjoint barycentric index
630 0 : const size_t i0 = p - ind[0] - ind[1];
631 0 : const number x0 = 1.0 - x[0] - x[1];
632 :
633 : return m_vPolynom[ ind[0] ].value(x[0])
634 0 : * m_vPolynom[ ind[1] ].value(x[1])
635 0 : * m_vPolynom[ i0 ].value(x0);
636 : }
637 :
638 :
639 : /// \copydoc ug::LocalShapeFunctionSet::shape()
640 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
641 : {
642 0 : grad(g, multi_index(i), x);
643 : }
644 :
645 : /// evaluates the gradient
646 0 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
647 : const MathVector<dim>& x) const
648 : {
649 : check_multi_index(ind);
650 : //ReferenceTriangle::check_position(x);
651 :
652 : // get adjoint barycentric index and position
653 0 : const int i0 = p - ind[0] - ind[1];
654 0 : const number x0 = 1.0 - x[0] - x[1];
655 :
656 : UG_ASSERT(i0 <= (int)p && i0 >= 0, "Wrong Multiindex.");
657 : UG_ASSERT(x0 <= 1.0+SMALL && x0 >= -SMALL, "Wrong Position.");
658 :
659 : // loop dimensions
660 0 : for(int d = 0; d < dim; ++d)
661 : {
662 0 : g[d] = m_vDPolynom[ind[d]].value(x[d])
663 0 : * m_vPolynom[i0].value(x0);
664 0 : g[d] += (-1) * m_vDPolynom[i0].value(x0)
665 0 : * m_vPolynom[ind[d]].value(x[d]);
666 :
667 : // multiply by all functions not depending on x[d]
668 0 : for(int d2 = 0; d2 < dim; ++d2)
669 : {
670 : // skip own value
671 0 : if(d2 == d) continue;
672 :
673 0 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
674 : }
675 : }
676 0 : }
677 :
678 : /// return Multi index for index i
679 : inline const MathVector<dim,int>& multi_index(size_t i) const
680 : {
681 : check_index(i);
682 : return m_vMultiIndex[i];
683 : }
684 :
685 : /// return the index for a multi_index
686 0 : inline size_t index(const MathVector<dim,int>& ind) const
687 : {
688 : check_multi_index(ind);
689 0 : for(size_t i=0; i<nsh; ++i)
690 0 : if(multi_index(i) == ind) return i;
691 0 : UG_THROW("Index not found in LagrangeLSFS");
692 : }
693 :
694 : /// return the index for a multi_index
695 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
696 : {
697 : check_multi_index(ind);
698 :
699 : size_t res = ind[0];
700 : for(int i = 0; i < ind[1]; ++i)
701 : res += (p+1-i);
702 :
703 : check_index(res);
704 : return res;
705 : }
706 :
707 : /// return the multi_index for an index
708 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
709 : {
710 : check_index(i);
711 :
712 : int i0 = i, i1;
713 : for(i1 = 0; i1 < (int)p; ++i1)
714 : {
715 : const int diff = i0 - (p+1-i1);
716 : if(diff < 0) break;
717 : i0 = diff;
718 : }
719 :
720 : UG_ASSERT(i0 >= 0, "i0 is negative ("<<i0<<")");
721 : UG_ASSERT(i1 >= 0, "i1 is negative ("<<i1<<")");
722 : return MathVector<dim,int>( i0, i1 );
723 : }
724 :
725 : /// checks in debug mode that index is valid
726 : inline void check_index(size_t i) const
727 : {
728 : UG_ASSERT(i < nsh, "Wrong index.");
729 : }
730 :
731 : /// checks in debug mode that multi-index is valid
732 : inline void check_multi_index(const MathVector<dim,int>& ind) const
733 : {
734 : UG_ASSERT(ind[0] <= (int)p && ind[0]>=0, "Wrong Multiindex.");
735 : UG_ASSERT(ind[1] <= (int)p && ind[0]>=0, "Wrong Multiindex.");
736 : UG_ASSERT(ind[0] + ind[1] <= (int)p && ind[0]>=0, "Wrong Multiindex.");
737 : }
738 :
739 : private:
740 : /// order
741 : size_t p;
742 :
743 : /// Number of shape functions
744 : size_t nsh;
745 :
746 : std::vector<Polynomial1D> m_vPolynom; ///< Shape Polynomials
747 : std::vector<Polynomial1D> m_vDPolynom; ///< Derivative of Shape Polynomial
748 :
749 : std::vector<MathVector<dim,int> > m_vMultiIndex;
750 : };
751 :
752 : ///////////////////////////////////////////////////////////////////////////////
753 : // Quadrilateral
754 : ///////////////////////////////////////////////////////////////////////////////
755 :
756 : template <int TOrder>
757 : class LagrangeLSFS<ReferenceQuadrilateral, TOrder>
758 : : public LagrangeLDS<ReferenceQuadrilateral>,
759 : public BaseLSFS<LagrangeLSFS<ReferenceQuadrilateral, TOrder>, 2>
760 : {
761 : private:
762 : /// abbreviation for order
763 : static const size_t p = TOrder;
764 :
765 : /// base class
766 : typedef BaseLSFS<LagrangeLSFS<ReferenceQuadrilateral, TOrder>, 2> base_type;
767 :
768 : public:
769 : /// Shape type
770 : typedef typename base_type::shape_type shape_type;
771 :
772 : /// Gradient type
773 : typedef typename base_type::grad_type grad_type;
774 :
775 : public:
776 : /// Order of Shape functions
777 : static const size_t order = TOrder;
778 :
779 : /// Dimension, where shape functions are defined
780 : static const int dim = ReferenceQuadrilateral::dim;
781 :
782 : /// Number of shape functions
783 : static const size_t nsh = (p+1)*(p+1);
784 :
785 : public:
786 : /// Constructor
787 : LagrangeLSFS();
788 :
789 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
790 0 : inline bool continuous() const {return true;}
791 :
792 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
793 1920 : inline size_t num_sh() const {return nsh;}
794 :
795 : /// \copydoc ug::LocalShapeFunctionSet::position()
796 0 : inline bool position(size_t i, MathVector<dim>& pos) const
797 : {
798 : // get Multi Index
799 : MathVector<dim,int> ind = multi_index(i);
800 :
801 : // set position
802 213 : for(int d = 0; d < dim; ++d)
803 142 : pos[d] = EquidistantLagrange1D::position(ind[d], p);
804 :
805 0 : return true;
806 : }
807 :
808 : /// \copydoc ug::LocalShapeFunctionSet::shape()
809 0 : inline number shape(const size_t i, const MathVector<dim>& x) const
810 : {
811 : // forward
812 3026 : return shape(multi_index(i), x);
813 : }
814 :
815 : /// shape value for a Multi Index
816 3026 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
817 : {
818 : check_multi_index(ind);
819 : //ReferenceQuadrilateral::check_position(x);
820 :
821 3026 : return m_vPolynom[ ind[0] ].value(x[0])
822 3026 : * m_vPolynom[ ind[1] ].value(x[1]);
823 : }
824 :
825 : /// \copydoc ug::LocalShapeFunctionSet::grad()
826 0 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
827 : {
828 1160 : grad(g, multi_index(i), x);
829 0 : }
830 :
831 : /// evaluates the gradient
832 1160 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
833 : const MathVector<dim>& x) const
834 : {
835 : check_multi_index(ind);
836 : //ReferenceQuadrilateral::check_position(x);
837 :
838 : // loop dimensions
839 3480 : for(int d = 0; d < dim; ++d)
840 : {
841 2320 : g[d] = m_vDPolynom[ind[d]].value(x[d]);
842 :
843 : // multiply by all functions not depending on x[d]
844 6960 : for(int d2 = 0; d2 < dim; ++d2)
845 : {
846 : // skip own value
847 4640 : if(d2 == d) continue;
848 :
849 4640 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
850 : }
851 : }
852 1160 : }
853 :
854 : /// return Multi index for index i
855 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
856 : {
857 : check_index(i);
858 3039 : return m_vMultiIndex[i];
859 : }
860 :
861 : /// return the index for a multi_index
862 29 : inline size_t index(const MathVector<dim,int>& ind) const
863 : {
864 : check_multi_index(ind);
865 191 : for(size_t i=0; i<nsh; ++i)
866 191 : if(multi_index(i) == ind) return i;
867 0 : UG_THROW("Index not found in LagrangeLSFS");
868 : }
869 :
870 : /// return the index for a multi_index
871 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
872 : {
873 : check_multi_index(ind);
874 :
875 0 : return ind[1] * (p+1) + ind[0];
876 : }
877 :
878 : /// return the multi_index for an index
879 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
880 : {
881 : check_index(i);
882 :
883 0 : return MathVector<dim,int>( i%(p+1), i/(p+1) );
884 : }
885 :
886 : /// checks in debug mode that index is valid
887 0 : inline void check_index(size_t i) const
888 : {
889 : UG_ASSERT(i < nsh, "Wrong index.");
890 0 : }
891 :
892 : /// checks in debug mode that multi-index is valid
893 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
894 : {
895 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
896 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
897 0 : }
898 :
899 : private:
900 : Polynomial1D m_vPolynom[p+1];
901 : Polynomial1D m_vDPolynom[p+1];
902 :
903 : MathVector<dim,int> m_vMultiIndex[nsh];
904 : };
905 :
906 :
907 : template <>
908 : class FlexLagrangeLSFS<ReferenceQuadrilateral>
909 : : public LagrangeLDS<ReferenceQuadrilateral>,
910 : public BaseLSFS<FlexLagrangeLSFS<ReferenceQuadrilateral>, 2>
911 : {
912 : public:
913 : /// Dimension, where shape functions are defined
914 : static const int dim = ReferenceQuadrilateral::dim;
915 :
916 : public:
917 : /// default Constructor
918 0 : FlexLagrangeLSFS() {set_order(1);}
919 :
920 : /// Constructor
921 0 : FlexLagrangeLSFS(size_t order) {set_order(order);}
922 :
923 : /// sets the order
924 : void set_order(size_t order);
925 :
926 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
927 : inline bool continuous() const {return true;}
928 :
929 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
930 0 : inline size_t num_sh() const {return nsh;}
931 :
932 : /// \copydoc ug::LocalShapeFunctionSet::position()
933 : inline bool position(size_t i, MathVector<dim>& pos) const
934 : {
935 : // get Multi Index
936 : MathVector<dim,int> ind = multi_index(i);
937 :
938 : // set position
939 0 : for(int d = 0; d < dim; ++d)
940 0 : pos[d] = EquidistantLagrange1D::position(ind[d], p);
941 :
942 : return true;
943 : }
944 :
945 : /// \copydoc ug::LocalShapeFunctionSet::shape()
946 : inline number shape(const size_t i, const MathVector<dim>& x) const
947 : {
948 : // forward
949 0 : return shape(multi_index(i), x);
950 : }
951 :
952 : /// shape value for a Multi Index
953 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
954 : {
955 : check_multi_index(ind);
956 : //ReferenceQuadrilateral::check_position(x);
957 :
958 0 : return m_vPolynom[ ind[0] ].value(x[0])
959 0 : * m_vPolynom[ ind[1] ].value(x[1]);
960 : }
961 :
962 : /// \copydoc ug::LocalShapeFunctionSet::grad()
963 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
964 : {
965 0 : grad(g, multi_index(i), x);
966 : }
967 :
968 : /// evaluates the gradient
969 0 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
970 : const MathVector<dim>& x) const
971 : {
972 : check_multi_index(ind);
973 : //ReferenceQuadrilateral::check_position(x);
974 :
975 : // loop dimensions
976 0 : for(int d = 0; d < dim; ++d)
977 : {
978 0 : g[d] = m_vDPolynom[ind[d]].value(x[d]);
979 :
980 : // multiply by all functions not depending on x[d]
981 0 : for(int d2 = 0; d2 < dim; ++d2)
982 : {
983 : // skip own value
984 0 : if(d2 == d) continue;
985 :
986 0 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
987 : }
988 : }
989 0 : }
990 :
991 : /// return Multi index for index i
992 : inline const MathVector<dim,int>& multi_index(size_t i) const
993 : {
994 : check_index(i);
995 : return m_vMultiIndex[i];
996 : }
997 :
998 : /// return the index for a multi_index
999 0 : inline size_t index(const MathVector<dim,int>& ind) const
1000 : {
1001 : check_multi_index(ind);
1002 0 : for(size_t i=0; i<nsh; ++i)
1003 0 : if(multi_index(i) == ind) return i;
1004 0 : UG_THROW("Index not found in LagrangeLSFS");
1005 : }
1006 :
1007 : /// return the index for a multi_index
1008 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
1009 : {
1010 : check_multi_index(ind);
1011 :
1012 : return ind[1] * (p+1) + ind[0];
1013 : }
1014 :
1015 : /// return the multi_index for an index
1016 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
1017 : {
1018 : check_index(i);
1019 :
1020 : return MathVector<dim,int>( i%(p+1), i/(p+1) );
1021 : }
1022 :
1023 : /// checks in debug mode that index is valid
1024 : inline void check_index(size_t i) const
1025 : {
1026 : UG_ASSERT(i < nsh, "Wrong index.");
1027 : }
1028 :
1029 : /// checks in debug mode that multi-index is valid
1030 : inline void check_multi_index(const MathVector<dim,int>& ind) const
1031 : {
1032 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
1033 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
1034 : }
1035 :
1036 : private:
1037 : /// order
1038 : size_t p;
1039 :
1040 : /// Number of shape functions
1041 : size_t nsh;
1042 :
1043 : std::vector<Polynomial1D> m_vPolynom; ///< Shape Polynomials
1044 : std::vector<Polynomial1D> m_vDPolynom; ///< Derivative of Shape Polynomial
1045 :
1046 : std::vector<MathVector<dim,int> > m_vMultiIndex;
1047 : };
1048 :
1049 : ///////////////////////////////////////////////////////////////////////////////
1050 : // Tetrahedron
1051 : ///////////////////////////////////////////////////////////////////////////////
1052 :
1053 : template <int TOrder>
1054 : class LagrangeLSFS<ReferenceTetrahedron, TOrder>
1055 : : public LagrangeLDS<ReferenceTetrahedron>,
1056 : public BaseLSFS<LagrangeLSFS<ReferenceTetrahedron, TOrder>, 3>
1057 : {
1058 : private:
1059 : /// abbreviation for order
1060 : static const size_t p = TOrder;
1061 :
1062 : /// base class
1063 : typedef BaseLSFS<LagrangeLSFS<ReferenceTetrahedron, TOrder>, 3> base_type;
1064 :
1065 : public:
1066 : /// Shape type
1067 : typedef typename base_type::shape_type shape_type;
1068 :
1069 : /// Gradient type
1070 : typedef typename base_type::grad_type grad_type;
1071 :
1072 : public:
1073 : /// Order of Shape functions
1074 : static const size_t order = TOrder;
1075 :
1076 : /// Dimension, where shape functions are defined
1077 : static const int dim = ReferenceTetrahedron::dim;
1078 :
1079 : /// Number of shape functions
1080 : static const size_t nsh = BinomialCoefficient<dim + p, p>::value;
1081 :
1082 : public:
1083 : /// Constructor
1084 : LagrangeLSFS();
1085 :
1086 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
1087 0 : inline bool continuous() const {return true;}
1088 :
1089 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
1090 2220 : inline size_t num_sh() const {return nsh;}
1091 :
1092 : /// \copydoc ug::LocalShapeFunctionSet::position()
1093 0 : inline bool position(size_t i, MathVector<dim>& pos) const
1094 : {
1095 : // get Multi Index
1096 : MathVector<dim,int> ind = multi_index(i);
1097 :
1098 : // set position
1099 288 : for(int d = 0; d < dim; ++d)
1100 216 : pos[d] = TruncatedEquidistantLagrange1D::position(ind[d], p);
1101 :
1102 0 : return true;
1103 : }
1104 :
1105 : /// \copydoc ug::LocalShapeFunctionSet::shape()
1106 0 : inline number shape(const size_t i, const MathVector<dim>& x) const
1107 : {
1108 : // forward
1109 3752 : return shape(multi_index(i), x);
1110 : }
1111 :
1112 : /// shape value for a Multi Index
1113 3752 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
1114 : {
1115 : check_multi_index(ind);
1116 : //ReferenceTetrahedron::check_position(x);
1117 :
1118 : // get adjoint barycentric index
1119 3752 : const size_t i0 = p - ind[0] - ind[1] - ind[2];
1120 3752 : const number x0 = 1.0 - x[0] - x[1] - x[2];
1121 :
1122 3752 : return m_vPolynom[ ind[0] ].value(x[0])
1123 3752 : * m_vPolynom[ ind[1] ].value(x[1])
1124 3752 : * m_vPolynom[ ind[2] ].value(x[2])
1125 3752 : * m_vPolynom[ i0 ].value(x0);
1126 : }
1127 :
1128 : /// \copydoc ug::LocalShapeFunctionSet::grad()
1129 1360 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
1130 : {
1131 1360 : grad(g, multi_index(i), x);
1132 1360 : }
1133 :
1134 : /// evaluates the gradient
1135 1360 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
1136 : const MathVector<dim>& x) const
1137 : {
1138 : check_multi_index(ind);
1139 : //ReferenceTetrahedron::check_position(x);
1140 :
1141 : // get adjoint barycentric index and position
1142 1360 : const size_t i0 = p - ind[0] - ind[1] - ind[2];
1143 1360 : const number x0 = 1 - x[0] - x[1] - x[2];
1144 :
1145 : UG_ASSERT(i0 <= p, "Wrong Multiindex.");
1146 : UG_ASSERT(x0 <= 1.0+SMALL && x0 >= -SMALL, "Wrong Position.");
1147 :
1148 : // loop dimensions
1149 5440 : for(int d = 0; d < dim; ++d)
1150 : {
1151 4080 : g[d] = m_vDPolynom[ind[d]].value(x[d])
1152 4080 : * m_vPolynom[i0].value(x0);
1153 4080 : g[d] += (-1) * m_vDPolynom[i0].value(x0)
1154 4080 : * m_vPolynom[ind[d]].value(x[d]);
1155 :
1156 : // multiply by all functions not depending on x[d]
1157 16320 : for(int d2 = 0; d2 < dim; ++d2)
1158 : {
1159 : // skip own value
1160 12240 : if(d2 == d) continue;
1161 :
1162 16320 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
1163 : }
1164 : }
1165 1360 : }
1166 :
1167 : /// return Multi index for index i
1168 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
1169 : {
1170 : check_index(i);
1171 3756 : return m_vMultiIndex[i];
1172 : }
1173 :
1174 : /// return the index for a multi_index
1175 34 : inline size_t index(const MathVector<dim,int>& ind) const
1176 : {
1177 : check_multi_index(ind);
1178 275 : for(size_t i=0; i<nsh; ++i)
1179 275 : if(multi_index(i) == ind) return i;
1180 0 : UG_THROW("Index not found in LagrangeLSFS");
1181 : }
1182 :
1183 : /// return the index for a multi_index
1184 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
1185 : {
1186 : check_multi_index(ind);
1187 :
1188 : // add x range
1189 0 : size_t res = ind[0];
1190 :
1191 : // add y range
1192 0 : for(int i = 0; i < ind[1]; ++i)
1193 0 : res += (p+1-ind[2]-i);
1194 :
1195 : // add z range
1196 0 : for(int i2 = 0; i2 < ind[2]; ++i2)
1197 0 : res += BinomCoeff(p+2-i2, p-i2);
1198 :
1199 : check_index(res);
1200 0 : return res;
1201 : }
1202 :
1203 : /// return the multi_index for an index
1204 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
1205 : {
1206 : check_index(i);
1207 :
1208 0 : int i0 = i, i1 = 0, i2 = 0;
1209 0 : for(i2 = 0; i2 <= (int)p; ++i2)
1210 : {
1211 0 : const int binom = BinomCoeff(p+2-i2, p-i2);
1212 :
1213 : // if i2 is correct
1214 0 : const int diff = i0 - binom;
1215 0 : if(diff < 0)
1216 : {
1217 0 : for(i1 = 0; i1 <= (int)p; ++i1)
1218 : {
1219 : // if i1 is correct return values
1220 0 : const int diff = i0 - (p+1-i2-i1);
1221 0 : if(diff < 0)
1222 : return MathVector<dim,int>( i0, i1, i2);
1223 :
1224 : // else decrease i1
1225 : i0 = diff;
1226 : }
1227 : }
1228 : // else go one level lower
1229 : else
1230 : i0 = diff;
1231 : }
1232 :
1233 : UG_ASSERT(i0 >= 0, "i0 is negative ("<<i0<<")");
1234 : UG_ASSERT(i1 >= 0, "i1 is negative ("<<i1<<")");
1235 : UG_ASSERT(i2 >= 0, "i1 is negative ("<<i2<<")");
1236 :
1237 : UG_ASSERT(0, "Should not reach this line.");
1238 : return MathVector<dim,int>( i0, i1, i2);
1239 : }
1240 :
1241 : /// checks in debug mode that index is valid
1242 0 : inline void check_index(size_t i) const
1243 : {
1244 : UG_ASSERT(i < nsh, "Wrong index.");
1245 0 : }
1246 :
1247 : /// checks in debug mode that multi-index is valid
1248 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
1249 : {
1250 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
1251 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
1252 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
1253 : UG_ASSERT(ind[0] + ind[1] + ind[2] <= (int)p, "Wrong Multiindex.");
1254 0 : }
1255 :
1256 : private:
1257 : Polynomial1D m_vPolynom[p+1];
1258 : Polynomial1D m_vDPolynom[p+1];
1259 :
1260 : MathVector<dim,int> m_vMultiIndex[nsh];
1261 : };
1262 :
1263 :
1264 : template <>
1265 : class FlexLagrangeLSFS<ReferenceTetrahedron>
1266 : : public LagrangeLDS<ReferenceTetrahedron>,
1267 : public BaseLSFS<FlexLagrangeLSFS<ReferenceTetrahedron>, 3>
1268 : {
1269 : public:
1270 : /// Dimension, where shape functions are defined
1271 : static const int dim = ReferenceTetrahedron::dim;
1272 :
1273 : public:
1274 : /// default Constructor
1275 0 : FlexLagrangeLSFS() {set_order(1);}
1276 :
1277 : /// Constructor
1278 0 : FlexLagrangeLSFS(size_t order) {set_order(order);}
1279 :
1280 : /// sets the order
1281 : void set_order(size_t order);
1282 :
1283 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
1284 : inline bool continuous() const {return true;}
1285 :
1286 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
1287 0 : inline size_t num_sh() const {return nsh;}
1288 :
1289 : /// \copydoc ug::LocalShapeFunctionSet::position()
1290 : inline bool position(size_t i, MathVector<dim>& pos) const
1291 : {
1292 : // get Multi Index
1293 : MathVector<dim,int> ind = multi_index(i);
1294 :
1295 : // set position
1296 0 : for(int d = 0; d < dim; ++d)
1297 0 : pos[d] = TruncatedEquidistantLagrange1D::position(ind[d], p);
1298 :
1299 : return true;
1300 : }
1301 :
1302 : /// \copydoc ug::LocalShapeFunctionSet::shape()
1303 : inline number shape(const size_t i, const MathVector<dim>& x) const
1304 : {
1305 : // forward
1306 0 : return shape(multi_index(i), x);
1307 : }
1308 :
1309 : /// shape value for a Multi Index
1310 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
1311 : {
1312 : check_multi_index(ind);
1313 : //ReferenceTetrahedron::check_position(x);
1314 :
1315 : // get adjoint barycentric index
1316 0 : const size_t i0 = p - ind[0] - ind[1] - ind[2];
1317 0 : const number x0 = 1.0 - x[0] - x[1] - x[2];
1318 :
1319 : return m_vPolynom[ ind[0] ].value(x[0])
1320 0 : * m_vPolynom[ ind[1] ].value(x[1])
1321 0 : * m_vPolynom[ ind[2] ].value(x[2])
1322 0 : * m_vPolynom[ i0 ].value(x0);
1323 : }
1324 :
1325 : /// \copydoc ug::LocalShapeFunctionSet::grad()
1326 0 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
1327 : {
1328 0 : grad(g, multi_index(i), x);
1329 0 : }
1330 :
1331 : /// evaluates the gradient
1332 0 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
1333 : const MathVector<dim>& x) const
1334 : {
1335 : check_multi_index(ind);
1336 : //ReferenceTetrahedron::check_position(x);
1337 :
1338 : // get adjoint barycentric index and position
1339 0 : const size_t i0 = p - ind[0] - ind[1] - ind[2];
1340 0 : const number x0 = 1 - x[0] - x[1] - x[2];
1341 :
1342 : UG_ASSERT(i0 <= p, "Wrong Multiindex.");
1343 : UG_ASSERT(x0 <= 1.0+SMALL && x0 >= -SMALL, "Wrong Position.");
1344 :
1345 : // loop dimensions
1346 0 : for(int d = 0; d < dim; ++d)
1347 : {
1348 0 : g[d] = m_vDPolynom[ind[d]].value(x[d])
1349 0 : * m_vPolynom[i0].value(x0);
1350 0 : g[d] += (-1) * m_vDPolynom[i0].value(x0)
1351 0 : * m_vPolynom[ind[d]].value(x[d]);
1352 :
1353 : // multiply by all functions not depending on x[d]
1354 0 : for(int d2 = 0; d2 < dim; ++d2)
1355 : {
1356 : // skip own value
1357 0 : if(d2 == d) continue;
1358 :
1359 0 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
1360 : }
1361 : }
1362 0 : }
1363 :
1364 : /// return Multi index for index i
1365 : inline const MathVector<dim,int>& multi_index(size_t i) const
1366 : {
1367 : check_index(i);
1368 : return m_vMultiIndex[i];
1369 : }
1370 :
1371 : /// return the index for a multi_index
1372 0 : inline size_t index(const MathVector<dim,int>& ind) const
1373 : {
1374 : check_multi_index(ind);
1375 0 : for(size_t i=0; i<nsh; ++i)
1376 0 : if(multi_index(i) == ind) return i;
1377 0 : UG_THROW("Index not found in LagrangeLSFS");
1378 : }
1379 :
1380 : /// return the index for a multi_index
1381 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
1382 : {
1383 : check_multi_index(ind);
1384 :
1385 : // add x range
1386 : size_t res = ind[0];
1387 :
1388 : // add y range
1389 : for(int i = 0; i < ind[1]; ++i)
1390 : res += (p+1-ind[2]-i);
1391 :
1392 : // add z range
1393 : for(int i2 = 0; i2 < ind[2]; ++i2)
1394 : res += BinomCoeff(p+2-i2, p-i2);
1395 :
1396 : check_index(res);
1397 : return res;
1398 : }
1399 :
1400 : /// return the multi_index for an index
1401 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
1402 : {
1403 : check_index(i);
1404 :
1405 : int i0 = i, i1 = 0, i2 = 0;
1406 : for(i2 = 0; i2 <= (int)p; ++i2)
1407 : {
1408 : const int binom = BinomCoeff(p+2-i2, p-i2);
1409 :
1410 : // if i2 is correct
1411 : const int diff = i0 - binom;
1412 : if(diff < 0)
1413 : {
1414 : for(i1 = 0; i1 <= (int)p; ++i1)
1415 : {
1416 : // if i1 is correct return values
1417 : const int diff = i0 - (p+1-i2-i1);
1418 : if(diff < 0)
1419 : return MathVector<dim,int>( i0, i1, i2);
1420 :
1421 : // else decrease i1
1422 : i0 = diff;
1423 : }
1424 : }
1425 : // else go one level lower
1426 : else
1427 : i0 = diff;
1428 : }
1429 :
1430 : UG_ASSERT(i0 >= 0, "i0 is negative ("<<i0<<")");
1431 : UG_ASSERT(i1 >= 0, "i1 is negative ("<<i1<<")");
1432 : UG_ASSERT(i2 >= 0, "i1 is negative ("<<i2<<")");
1433 :
1434 : UG_ASSERT(0, "Should not reach this line.");
1435 : return MathVector<dim,int>( i0, i1, i2);
1436 : }
1437 :
1438 : /// checks in debug mode that index is valid
1439 : inline void check_index(size_t i) const
1440 : {
1441 : UG_ASSERT(i < nsh, "Wrong index.");
1442 : }
1443 :
1444 : /// checks in debug mode that multi-index is valid
1445 : inline void check_multi_index(const MathVector<dim,int>& ind) const
1446 : {
1447 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
1448 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
1449 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
1450 : UG_ASSERT(ind[0] + ind[1] + ind[2] <= (int)p, "Wrong Multiindex.");
1451 : }
1452 :
1453 : private:
1454 : /// order
1455 : size_t p;
1456 :
1457 : /// Number of shape functions
1458 : size_t nsh;
1459 :
1460 : std::vector<Polynomial1D> m_vPolynom; ///< Shape Polynomials
1461 : std::vector<Polynomial1D> m_vDPolynom; ///< Derivative of Shape Polynomial
1462 :
1463 : std::vector<MathVector<dim,int> > m_vMultiIndex;
1464 : };
1465 :
1466 : ///////////////////////////////////////////////////////////////////////////////
1467 : // Prism
1468 : ///////////////////////////////////////////////////////////////////////////////
1469 :
1470 : template <int TOrder>
1471 : class LagrangeLSFS<ReferencePrism, TOrder>
1472 : : public LagrangeLDS<ReferencePrism>,
1473 : public BaseLSFS<LagrangeLSFS<ReferencePrism, TOrder>, 3>
1474 : {
1475 : private:
1476 : /// abbreviation for order
1477 : static const size_t p = TOrder;
1478 :
1479 : /// dofs per layer
1480 : static const size_t dofPerLayer = BinomialCoefficient<2 + p, p>::value;
1481 :
1482 : /// base class
1483 : typedef BaseLSFS<LagrangeLSFS<ReferencePrism, TOrder>, 3> base_type;
1484 :
1485 : public:
1486 : /// Shape type
1487 : typedef typename base_type::shape_type shape_type;
1488 :
1489 : /// Gradient type
1490 : typedef typename base_type::grad_type grad_type;
1491 :
1492 : public:
1493 : /// Order of Shape functions
1494 : static const size_t order = TOrder;
1495 :
1496 : /// Dimension, where shape functions are defined
1497 : static const int dim = ReferencePrism::dim;
1498 :
1499 : /// Number of shape functions
1500 : static const size_t nsh = dofPerLayer*(p+1);
1501 :
1502 : public:
1503 : /// Constructor
1504 : LagrangeLSFS();
1505 :
1506 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
1507 0 : inline bool continuous() const {return true;}
1508 :
1509 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
1510 4020 : inline size_t num_sh() const {return nsh;}
1511 :
1512 : /// \copydoc ug::LocalShapeFunctionSet::position()
1513 0 : inline bool position(size_t i, MathVector<dim>& pos) const
1514 : {
1515 : // get Multi Index
1516 : MathVector<dim,int> ind = multi_index(i);
1517 :
1518 : // set position
1519 402 : for(int d = 0; d < 2; ++d)
1520 268 : pos[d] = TruncatedEquidistantLagrange1D::position(ind[d], p);
1521 :
1522 134 : pos[2] = EquidistantLagrange1D::position(ind[2], p);
1523 :
1524 0 : return true;
1525 : }
1526 :
1527 : /// \copydoc ug::LocalShapeFunctionSet::shape()
1528 0 : inline number shape(const size_t i, const MathVector<dim>& x) const
1529 : {
1530 : // forward
1531 9040 : return shape(multi_index(i), x);
1532 : }
1533 :
1534 : /// shape value for a Multi Index
1535 9040 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
1536 : {
1537 : check_multi_index(ind);
1538 : //ReferencePrism::check_position(x);
1539 :
1540 : // get adjoint barycentric index
1541 9040 : const size_t i0 = p - ind[0] - ind[1];
1542 9040 : const number x0 = 1.0 - x[0] - x[1];
1543 :
1544 : // x-y direction
1545 9040 : return m_vTruncPolynom[ ind[0] ].value(x[0])
1546 9040 : * m_vTruncPolynom[ ind[1] ].value(x[1])
1547 9040 : * m_vTruncPolynom[ i0 ].value( x0 )
1548 : // z direction
1549 9040 : * m_vPolynom[ ind[2] ].value(x[2]);
1550 : }
1551 :
1552 : /// \copydoc ug::LocalShapeFunctionSet::grad()
1553 2560 : void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
1554 : {
1555 2560 : grad(g, multi_index(i), x);
1556 2560 : }
1557 :
1558 : /// evaluates the gradient
1559 2560 : void grad(grad_type& g, const MathVector<dim,int> ind,
1560 : const MathVector<dim>& x) const
1561 : {
1562 : check_multi_index(ind);
1563 : //ReferencePrism::check_position(x);
1564 :
1565 : // get adjoint barycentric index and position
1566 2560 : const size_t i0 = p - ind[0] - ind[1];
1567 2560 : const number x0 = 1 - x[0] - x[1];
1568 :
1569 : UG_ASSERT(i0 <= p, "Wrong Multiindex.");
1570 : UG_ASSERT(x0 <= 1.0+SMALL && x0 >= -SMALL, "Wrong Position.");
1571 :
1572 : // x-y gradient
1573 7680 : for(size_t d = 0; d < 2; ++d)
1574 : {
1575 5120 : g[d] = m_vDTruncPolynom[ind[d]].value(x[d])
1576 5120 : * m_vTruncPolynom[i0].value(x0);
1577 5120 : g[d] += (-1) * m_vDTruncPolynom[i0].value(x0)
1578 5120 : * m_vTruncPolynom[ind[d]].value(x[d]);
1579 :
1580 : // multiply by all functions not depending on x[d]
1581 15360 : for(size_t d2 = 0; d2 < 2; ++d2)
1582 : {
1583 : // skip own value
1584 10240 : if(d2 == d) continue;
1585 :
1586 10240 : g[d] *= m_vTruncPolynom[ind[d2]].value(x[d2]);
1587 : }
1588 :
1589 : // multiply by z coordinate
1590 10240 : g[d] *= m_vPolynom[ind[2]].value(x[2]);
1591 : }
1592 :
1593 : // z gradient
1594 2560 : g[2] = m_vDPolynom[ind[2]].value(x[2]);
1595 2560 : g[2] *= m_vTruncPolynom[ ind[0] ].value(x[0])
1596 2560 : * m_vTruncPolynom[ ind[1] ].value(x[1])
1597 2560 : * m_vTruncPolynom[ i0 ].value( x0 );
1598 2560 : }
1599 :
1600 : /// return Multi index for index i
1601 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
1602 : {
1603 : check_index(i);
1604 9046 : return m_vMultiIndex[i];
1605 : }
1606 :
1607 : /// return the index for a multi_index
1608 64 : inline size_t index(const MathVector<dim,int>& ind) const
1609 : {
1610 : check_multi_index(ind);
1611 1012 : for(size_t i=0; i<nsh; ++i)
1612 1012 : if(multi_index(i) == ind) return i;
1613 0 : UG_THROW("Index not found in LagrangeLSFS");
1614 : }
1615 :
1616 : /// return the index for a multi_index
1617 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
1618 : {
1619 : check_multi_index(ind);
1620 :
1621 0 : size_t res = ind[0];
1622 0 : for(int i = 0; i < ind[1]; ++i)
1623 0 : res += (p+1-i);
1624 :
1625 : // add height
1626 0 : res += ind[2] * dofPerLayer;
1627 :
1628 0 : return res;
1629 : }
1630 :
1631 : /// return the multi_index for an index
1632 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
1633 : {
1634 : check_index(i);
1635 :
1636 0 : const size_t i2 = i / dofPerLayer;
1637 :
1638 0 : int i0 = i - i2*dofPerLayer, i1;
1639 0 : for(i1 = 0; i1 < (int)p; ++i1)
1640 : {
1641 0 : const int diff = i0 - (p+1-i1);
1642 0 : if(diff < 0)
1643 : break;
1644 : i0 = diff;
1645 : }
1646 :
1647 0 : return MathVector<dim,int>( i0, i1, i2);
1648 : }
1649 :
1650 : /// checks in debug mode that index is valid
1651 0 : inline void check_index(size_t i) const
1652 : {
1653 : UG_ASSERT(i < nsh, "Wrong index.");
1654 0 : }
1655 :
1656 : /// checks in debug mode that multi-index is valid
1657 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
1658 : {
1659 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
1660 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
1661 : UG_ASSERT(ind[0] + ind[1] <= (int)p, "Wrong Multiindex.");
1662 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
1663 0 : }
1664 :
1665 : private:
1666 : Polynomial1D m_vPolynom[p+1];
1667 : Polynomial1D m_vDPolynom[p+1];
1668 : Polynomial1D m_vTruncPolynom[p+1];
1669 : Polynomial1D m_vDTruncPolynom[p+1];
1670 :
1671 : MathVector<dim,int> m_vMultiIndex[nsh];
1672 : };
1673 :
1674 :
1675 : template <>
1676 : class FlexLagrangeLSFS<ReferencePrism>
1677 : : public LagrangeLDS<ReferencePrism>,
1678 : public BaseLSFS<FlexLagrangeLSFS<ReferencePrism>, 3>
1679 : {
1680 : public:
1681 : /// Dimension, where shape functions are defined
1682 : static const int dim = ReferencePrism::dim;
1683 :
1684 : public:
1685 : /// default Constructor
1686 0 : FlexLagrangeLSFS() {set_order(1);}
1687 :
1688 : /// Constructor
1689 0 : FlexLagrangeLSFS(size_t order) {set_order(order);}
1690 :
1691 : /// sets the order
1692 : void set_order(size_t order);
1693 :
1694 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
1695 : inline bool continuous() const {return true;}
1696 :
1697 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
1698 0 : inline size_t num_sh() const {return nsh;}
1699 :
1700 : /// \copydoc ug::LocalShapeFunctionSet::position()
1701 0 : inline bool position(size_t i, MathVector<dim>& pos) const
1702 : {
1703 : // get Multi Index
1704 : MathVector<dim,int> ind = multi_index(i);
1705 :
1706 : // set position
1707 0 : for(int d = 0; d < 2; ++d)
1708 0 : pos[d] = TruncatedEquidistantLagrange1D::position(ind[d], p);
1709 :
1710 0 : pos[2] = EquidistantLagrange1D::position(ind[2], p);
1711 :
1712 0 : return true;
1713 : }
1714 :
1715 : /// \copydoc ug::LocalShapeFunctionSet::shape()
1716 : inline number shape(const size_t i, const MathVector<dim>& x) const
1717 : {
1718 : // forward
1719 0 : return shape(multi_index(i), x);
1720 : }
1721 :
1722 : /// shape value for a Multi Index
1723 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
1724 : {
1725 : check_multi_index(ind);
1726 : //ReferencePrism::check_position(x);
1727 :
1728 : // get adjoint barycentric index
1729 0 : const size_t i0 = p - ind[0] - ind[1];
1730 0 : const number x0 = 1.0 - x[0] - x[1];
1731 :
1732 : // x-y direction
1733 : return m_vTruncPolynom[ ind[0] ].value(x[0])
1734 0 : * m_vTruncPolynom[ ind[1] ].value(x[1])
1735 0 : * m_vTruncPolynom[ i0 ].value( x0 )
1736 : // z direction
1737 0 : * m_vPolynom[ ind[2] ].value(x[2]);
1738 : }
1739 :
1740 : /// \copydoc ug::LocalShapeFunctionSet::grad()
1741 0 : void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
1742 : {
1743 0 : grad(g, multi_index(i), x);
1744 0 : }
1745 :
1746 : /// evaluates the gradient
1747 0 : void grad(grad_type& g, const MathVector<dim,int> ind,
1748 : const MathVector<dim>& x) const
1749 : {
1750 : check_multi_index(ind);
1751 : //ReferencePrism::check_position(x);
1752 :
1753 : // get adjoint barycentric index and position
1754 0 : const size_t i0 = p - ind[0] - ind[1];
1755 0 : const number x0 = 1 - x[0] - x[1];
1756 :
1757 : UG_ASSERT(i0 <= p, "Wrong Multiindex.");
1758 : UG_ASSERT(x0 <= 1.0+SMALL && x0 >= -SMALL, "Wrong Position.");
1759 :
1760 : // x-y gradient
1761 0 : for(size_t d = 0; d < 2; ++d)
1762 : {
1763 0 : g[d] = m_vDTruncPolynom[ind[d]].value(x[d])
1764 0 : * m_vTruncPolynom[i0].value(x0);
1765 0 : g[d] += (-1) * m_vDTruncPolynom[i0].value(x0)
1766 0 : * m_vTruncPolynom[ind[d]].value(x[d]);
1767 :
1768 : // multiply by all functions not depending on x[d]
1769 0 : for(size_t d2 = 0; d2 < 2; ++d2)
1770 : {
1771 : // skip own value
1772 0 : if(d2 == d) continue;
1773 :
1774 0 : g[d] *= m_vTruncPolynom[ind[d2]].value(x[d2]);
1775 : }
1776 :
1777 : // multiply by z coordinate
1778 0 : g[d] *= m_vPolynom[ind[2]].value(x[2]);
1779 : }
1780 :
1781 : // z gradient
1782 0 : g[2] = m_vDPolynom[ind[2]].value(x[2]);
1783 0 : g[2] *= m_vTruncPolynom[ ind[0] ].value(x[0])
1784 0 : * m_vTruncPolynom[ ind[1] ].value(x[1])
1785 0 : * m_vTruncPolynom[ i0 ].value( x0 );
1786 0 : }
1787 :
1788 : /// return Multi index for index i
1789 : inline const MathVector<dim,int>& multi_index(size_t i) const
1790 : {
1791 : check_index(i);
1792 : return m_vMultiIndex[i];
1793 : }
1794 :
1795 : /// return the index for a multi_index
1796 0 : inline size_t index(const MathVector<dim,int>& ind) const
1797 : {
1798 : check_multi_index(ind);
1799 0 : for(size_t i=0; i<nsh; ++i)
1800 0 : if(multi_index(i) == ind) return i;
1801 0 : UG_THROW("Index not found in LagrangeLSFS");
1802 : }
1803 :
1804 : /// return the index for a multi_index
1805 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
1806 : {
1807 : check_multi_index(ind);
1808 :
1809 : size_t res = ind[0];
1810 : for(int i = 0; i < ind[1]; ++i)
1811 : res += (p+1-i);
1812 :
1813 : // add height
1814 : res += ind[2] * dofPerLayer;
1815 :
1816 : return res;
1817 : }
1818 :
1819 : /// return the multi_index for an index
1820 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
1821 : {
1822 : check_index(i);
1823 :
1824 : const size_t i2 = i / dofPerLayer;
1825 :
1826 : int i0 = i - i2*dofPerLayer, i1;
1827 : for(i1 = 0; i1 < (int)p; ++i1)
1828 : {
1829 : const int diff = i0 - (p+1-i1);
1830 : if(diff < 0)
1831 : break;
1832 : i0 = diff;
1833 : }
1834 :
1835 : return MathVector<dim,int>( i0, i1, i2);
1836 : }
1837 :
1838 : /// checks in debug mode that index is valid
1839 : inline void check_index(size_t i) const
1840 : {
1841 : UG_ASSERT(i < nsh, "Wrong index.");
1842 : }
1843 :
1844 : /// checks in debug mode that multi-index is valid
1845 : inline void check_multi_index(const MathVector<dim,int>& ind) const
1846 : {
1847 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
1848 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
1849 : UG_ASSERT(ind[0] + ind[1] <= (int)p, "Wrong Multiindex.");
1850 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
1851 : }
1852 :
1853 : private:
1854 : /// order
1855 : size_t p;
1856 :
1857 : /// Number of shape functions
1858 : size_t nsh;
1859 :
1860 : /// number of dofs per layer
1861 : size_t dofPerLayer;
1862 :
1863 : std::vector<Polynomial1D> m_vPolynom;
1864 : std::vector<Polynomial1D> m_vDPolynom;
1865 : std::vector<Polynomial1D> m_vTruncPolynom;
1866 : std::vector<Polynomial1D> m_vDTruncPolynom;
1867 :
1868 : std::vector<MathVector<dim,int> > m_vMultiIndex;
1869 : };
1870 :
1871 : ///////////////////////////////////////////////////////////////////////////////
1872 : // Pyramid
1873 : ///////////////////////////////////////////////////////////////////////////////
1874 :
1875 : namespace {
1876 : template <int p>
1877 : struct NumberOfDoFsOfPyramid{
1878 : enum{value = NumberOfDoFsOfPyramid<p-1>::value + (p+1)*(p+1)};
1879 : };
1880 :
1881 : // specialization to end recursion
1882 : template <> struct NumberOfDoFsOfPyramid<1>{enum {value = 5};};
1883 : template <> struct NumberOfDoFsOfPyramid<0>{enum {value = 0};};
1884 : template <> struct NumberOfDoFsOfPyramid<-1>{enum {value = 0};};
1885 : } // end empty namespace
1886 :
1887 : // todo: Implement higher order (impossible ?)
1888 : // NOTE: Currently only 1st order is implemented. There is no shape function
1889 : // set for pyramids, that is continuous and allows a continuous
1890 : // derivative in the inner of the pyramid. This is basically, since
1891 : // one regards the pyramid as two tetrahedrons, glued together.
1892 : template <int TOrder>
1893 : class LagrangeLSFS<ReferencePyramid, TOrder>
1894 : : public LagrangeLDS<ReferencePyramid>,
1895 : public BaseLSFS<LagrangeLSFS<ReferencePyramid, TOrder>, 3>
1896 : {
1897 : private:
1898 : /// abbreviation for order
1899 : static const size_t p = TOrder;
1900 :
1901 : /// base class
1902 : typedef BaseLSFS<LagrangeLSFS<ReferencePyramid, TOrder>, 3> base_type;
1903 :
1904 : public:
1905 : /// Shape type
1906 : typedef typename base_type::shape_type shape_type;
1907 :
1908 : /// Gradient type
1909 : typedef typename base_type::grad_type grad_type;
1910 :
1911 : public:
1912 : /// Order of Shape functions
1913 : static const size_t order = TOrder;
1914 :
1915 : /// Dimension, where shape functions are defined
1916 : static const int dim = 3; //reference_element_type::dim; (compile error on OSX 10.5)
1917 :
1918 : /// Number of shape functions
1919 : static const size_t nsh = NumberOfDoFsOfPyramid<p>::value;
1920 :
1921 : public:
1922 : /// Constructor
1923 : LagrangeLSFS();
1924 :
1925 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
1926 0 : inline bool continuous() const {return true;}
1927 :
1928 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
1929 360 : inline size_t num_sh() const {return nsh;}
1930 :
1931 : /// \copydoc ug::LocalShapeFunctionSet::position()
1932 0 : inline bool position(size_t i, MathVector<dim>& pos) const
1933 : {
1934 : // get Multi Index
1935 : MathVector<dim,int> ind = multi_index(i);
1936 :
1937 : // set position
1938 60 : for(int d = 0; d < dim; ++d)
1939 45 : pos[d] = EquidistantLagrange1D::position(ind[d], p);
1940 :
1941 0 : return true;
1942 : }
1943 :
1944 : /// \copydoc ug::LocalShapeFunctionSet::shape()
1945 450 : inline number shape(const size_t i, const MathVector<dim>& x) const
1946 : {
1947 : // only first order
1948 0 : if(p != 1) UG_THROW("Only 1. order Lagrange Pyramid implemented.");
1949 :
1950 : // smaller value of x and y
1951 450 : number m = x[0];
1952 450 : if (x[0] > x[1]) m = x[1];
1953 :
1954 450 : switch(i)
1955 : {
1956 90 : case 0 : return((1.0-x[0])*(1.0-x[1]) + x[2]*(m-1.0));
1957 90 : case 1 : return(x[0]*(1.0-x[1]) - x[2]*m);
1958 90 : case 2 : return(x[0]*x[1] + x[2]*m);
1959 90 : case 3 : return((1.0-x[0])*x[1] - x[2]*m);
1960 90 : case 4 : return(x[2]);
1961 0 : default: UG_THROW("Wrong index "<< i<<" in Pyramid");
1962 : }
1963 : }
1964 :
1965 : /// shape value for a Multi Index
1966 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
1967 : {
1968 : check_multi_index(ind);
1969 : //ReferencePyramid::check_position(x);
1970 :
1971 : // forward
1972 0 : return shape(index(ind), x);
1973 : }
1974 :
1975 : /// \copydoc ug::LocalShapeFunctionSet::grad()
1976 200 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
1977 : {
1978 : // only first order
1979 0 : if(p != 1) UG_THROW("Only 1. order Lagrange Pyramid implemented.");
1980 :
1981 : int m = 0;
1982 200 : if (x[0] > x[1]) m = 1;
1983 :
1984 200 : switch(i)
1985 : {
1986 : case 0:
1987 40 : g[0] = -(1.0-x[1]);
1988 40 : g[1] = -(1.0-x[0]) + x[2];
1989 40 : g[2] = -(1.0-x[1]);
1990 40 : g[m] += x[2];
1991 40 : break;
1992 : case 1:
1993 40 : g[0] = (1.0-x[1]);
1994 40 : g[1] = -x[0];
1995 40 : g[2] = -x[1];
1996 40 : g[m] -= x[2];
1997 40 : break;
1998 : case 2:
1999 40 : g[0] = x[1];
2000 40 : g[1] = x[0];
2001 40 : g[2] = x[1];
2002 40 : g[m] += x[2];
2003 40 : break;
2004 : case 3:
2005 40 : g[0] = -x[1];
2006 40 : g[1] = 1.0-x[0];
2007 40 : g[2] = -x[1];
2008 40 : g[m] -= x[2];
2009 40 : break;
2010 : case 4:
2011 40 : g[0] = 0.0;
2012 40 : g[1] = 0.0;
2013 40 : g[2] = 1.0;
2014 40 : break;
2015 0 : default: UG_THROW("Wrong index "<< i<<" in Pyramid");
2016 : }
2017 :
2018 200 : }
2019 :
2020 : /// evaluates the gradient
2021 0 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
2022 : const MathVector<dim>& x) const
2023 : {
2024 0 : grad(g, index(ind), x);
2025 0 : }
2026 :
2027 : /// return Multi index for index i
2028 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
2029 : {
2030 : check_index(i);
2031 5 : return m_vMultiIndex[i];
2032 : }
2033 :
2034 : /// return the index for a multi_index
2035 5 : inline size_t index(const MathVector<dim,int>& ind) const
2036 : {
2037 : check_multi_index(ind);
2038 15 : for(size_t i=0; i<nsh; ++i)
2039 15 : if(multi_index(i) == ind) return i;
2040 0 : UG_THROW("Index not found in LagrangeLSFS");
2041 : }
2042 :
2043 : /// return the index for a multi_index
2044 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
2045 : {
2046 : check_multi_index(ind);
2047 :
2048 : size_t res = 0;
2049 :
2050 : // add layers that are completely filled
2051 0 : for(int i2 = 0; i2 < ind[2]; ++i2)
2052 0 : res += (p+1-i2)*(p+1-i2);
2053 :
2054 : // add dofs of top z-layer
2055 0 : res += ind[1] * (p+1-ind[2]) + ind[0];
2056 :
2057 : check_index(res);
2058 0 : return res;
2059 : }
2060 :
2061 : /// return the multi_index for an index
2062 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
2063 : {
2064 : check_index(i);
2065 :
2066 : // get z layer
2067 0 : int iTmp = i;
2068 : int i2;
2069 0 : for(i2 = 0; i2 < (int)p; ++i2)
2070 : {
2071 0 : const int diff = iTmp - (p+1-i2)*(p+1-i2);
2072 0 : if(diff < 0) break;
2073 :
2074 : iTmp = diff;
2075 : }
2076 :
2077 0 : return MathVector<dim,int>( iTmp%(p+1-i2), iTmp/(p+1-i2), i2);
2078 : }
2079 :
2080 : /// checks in debug mode that index is valid
2081 0 : inline void check_index(size_t i) const
2082 : {
2083 : UG_ASSERT(i < nsh, "Wrong index.");
2084 0 : }
2085 :
2086 : /// checks in debug mode that multi-index is valid
2087 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
2088 : {
2089 : UG_ASSERT(ind[0] <= (int)p-ind[2] && ind[0] >= 0, "Wrong Multiindex.");
2090 : UG_ASSERT(ind[1] <= (int)p-ind[2] && ind[1] >= 0, "Wrong Multiindex.");
2091 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
2092 0 : }
2093 :
2094 : private:
2095 : std::vector<std::vector<Polynomial1D> > m_vvPolynom;
2096 : std::vector<std::vector<Polynomial1D> > m_vvDPolynom;
2097 :
2098 : MathVector<dim,int> m_vMultiIndex[nsh];
2099 : };
2100 :
2101 : ///////////////////////////////////////////////////////////////////////////////
2102 : // Hexahedron
2103 : ///////////////////////////////////////////////////////////////////////////////
2104 :
2105 : template <int TOrder>
2106 : class LagrangeLSFS<ReferenceHexahedron, TOrder>
2107 : : public LagrangeLDS<ReferenceHexahedron>,
2108 : public BaseLSFS<LagrangeLSFS<ReferenceHexahedron, TOrder>, 3>
2109 : {
2110 : private:
2111 : /// abbreviation for order
2112 : static const size_t p = TOrder;
2113 :
2114 : /// base class
2115 : typedef BaseLSFS<LagrangeLSFS<ReferenceHexahedron, TOrder>, 3> base_type;
2116 :
2117 : public:
2118 : /// Shape type
2119 : typedef typename base_type::shape_type shape_type;
2120 :
2121 : /// Gradient type
2122 : typedef typename base_type::grad_type grad_type;
2123 :
2124 : public:
2125 : /// Order of Shape functions
2126 : static const size_t order = TOrder;
2127 :
2128 : /// Dimension, where shape functions are defined
2129 : static const int dim = ReferenceHexahedron::dim;
2130 :
2131 : /// Number of shape functions
2132 : static const size_t nsh = (p+1)*(p+1)*(p+1);
2133 :
2134 : public:
2135 : /// Constructor
2136 : LagrangeLSFS();
2137 :
2138 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
2139 0 : inline bool continuous() const {return true;}
2140 :
2141 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
2142 6120 : inline size_t num_sh() const {return nsh;}
2143 :
2144 : /// \copydoc ug::LocalShapeFunctionSet::position()
2145 0 : inline bool position(size_t i, MathVector<dim>& pos) const
2146 : {
2147 : // get Multi Index
2148 : MathVector<dim,int> ind = multi_index(i);
2149 :
2150 : // set position
2151 824 : for(int d = 0; d < dim; ++d)
2152 618 : pos[d] = EquidistantLagrange1D::position(ind[d], p);
2153 :
2154 0 : return true;
2155 : }
2156 :
2157 : /// \copydoc ug::LocalShapeFunctionSet::shape()
2158 0 : inline number shape(const size_t i, const MathVector<dim>& x) const
2159 : {
2160 : // forward
2161 17698 : return shape(multi_index(i), x);
2162 : }
2163 :
2164 : /// shape value for a Multi Index
2165 17698 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
2166 : {
2167 : check_multi_index(ind);
2168 : //ReferenceHexahedron::check_position(x);
2169 :
2170 17698 : return m_vPolynom[ ind[0] ].value(x[0])
2171 17698 : * m_vPolynom[ ind[1] ].value(x[1])
2172 17698 : * m_vPolynom[ ind[2] ].value(x[2]);
2173 : }
2174 :
2175 : /// \copydoc ug::LocalShapeFunctionSet::grad()
2176 3960 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
2177 : {
2178 3960 : grad(g, multi_index(i), x);
2179 3960 : }
2180 :
2181 : /// evaluates the gradient
2182 3960 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
2183 : const MathVector<dim>& x) const
2184 : {
2185 : check_multi_index(ind);
2186 : //ReferenceHexahedron::check_position(x);
2187 :
2188 : // loop dimensions
2189 15840 : for(int d = 0; d < dim; ++d)
2190 : {
2191 11880 : g[d] = m_vDPolynom[ind[d]].value(x[d]);
2192 :
2193 : // multiply by all functions not depending on x[d]
2194 47520 : for(int d2 = 0; d2 < dim; ++d2)
2195 : {
2196 : // skip own value
2197 35640 : if(d2 == d) continue;
2198 :
2199 47520 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
2200 : }
2201 : }
2202 3960 : }
2203 :
2204 : /// return Multi index for index i
2205 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
2206 : {
2207 : check_index(i);
2208 17706 : return m_vMultiIndex[i];
2209 : }
2210 :
2211 : /// return the index for a multi_index
2212 99 : inline size_t index(const MathVector<dim,int>& ind) const
2213 : {
2214 : check_multi_index(ind);
2215 2494 : for(size_t i=0; i<nsh; ++i)
2216 2494 : if(multi_index(i) == ind) return i;
2217 0 : UG_THROW("Index not found in LagrangeLSFS");
2218 : }
2219 :
2220 : /// return the index for a multi_index
2221 0 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
2222 : {
2223 : check_multi_index(ind);
2224 :
2225 0 : return ind[2] * (p+1)*(p+1) + ind[1] * (p+1) + ind[0];
2226 : }
2227 :
2228 : /// return the multi_index for an index
2229 0 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
2230 : {
2231 : check_index(i);
2232 :
2233 0 : return MathVector<dim,int>( i%(p+1), i/(p+1)%(p+1), i/((p+1)*(p+1)));
2234 : }
2235 :
2236 : /// checks in debug mode that index is valid
2237 0 : inline void check_index(size_t i) const
2238 : {
2239 : UG_ASSERT(i < nsh, "Wrong index.");
2240 0 : }
2241 :
2242 : /// checks in debug mode that multi-index is valid
2243 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
2244 : {
2245 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
2246 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
2247 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
2248 0 : }
2249 :
2250 : private:
2251 : Polynomial1D m_vPolynom[p+1];
2252 : Polynomial1D m_vDPolynom[p+1];
2253 :
2254 : MathVector<dim,int> m_vMultiIndex[nsh];
2255 : };
2256 :
2257 : template <>
2258 : class FlexLagrangeLSFS<ReferenceHexahedron>
2259 : : public LagrangeLDS<ReferenceHexahedron>,
2260 : public BaseLSFS<FlexLagrangeLSFS<ReferenceHexahedron>, 3>
2261 : {
2262 : public:
2263 : /// Dimension, where shape functions are defined
2264 : static const int dim = ReferenceHexahedron::dim;
2265 :
2266 : public:
2267 : /// default Constructor
2268 0 : FlexLagrangeLSFS() {set_order(1);}
2269 :
2270 : /// Constructor
2271 0 : FlexLagrangeLSFS(size_t order) {set_order(order);}
2272 :
2273 : /// sets the order
2274 : void set_order(size_t order);
2275 :
2276 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
2277 : inline bool continuous() const {return true;}
2278 :
2279 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
2280 0 : inline size_t num_sh() const {return nsh;}
2281 :
2282 : /// \copydoc ug::LocalShapeFunctionSet::position()
2283 : inline bool position(size_t i, MathVector<dim>& pos) const
2284 : {
2285 : // get Multi Index
2286 : MathVector<dim,int> ind = multi_index(i);
2287 :
2288 : // set position
2289 0 : for(int d = 0; d < dim; ++d)
2290 0 : pos[d] = EquidistantLagrange1D::position(ind[d], p);
2291 :
2292 : return true;
2293 : }
2294 :
2295 : /// \copydoc ug::LocalShapeFunctionSet::shape()
2296 : inline number shape(const size_t i, const MathVector<dim>& x) const
2297 : {
2298 : // forward
2299 0 : return shape(multi_index(i), x);
2300 : }
2301 :
2302 : /// shape value for a Multi Index
2303 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
2304 : {
2305 : check_multi_index(ind);
2306 : //ReferenceHexahedron::check_position(x);
2307 :
2308 0 : return m_vPolynom[ ind[0] ].value(x[0])
2309 0 : * m_vPolynom[ ind[1] ].value(x[1])
2310 0 : * m_vPolynom[ ind[2] ].value(x[2]);
2311 : }
2312 :
2313 : /// \copydoc ug::LocalShapeFunctionSet::grad()
2314 0 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
2315 : {
2316 0 : grad(g, multi_index(i), x);
2317 0 : }
2318 :
2319 : /// evaluates the gradient
2320 0 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
2321 : const MathVector<dim>& x) const
2322 : {
2323 : check_multi_index(ind);
2324 : //ReferenceHexahedron::check_position(x);
2325 :
2326 : // loop dimensions
2327 0 : for(int d = 0; d < dim; ++d)
2328 : {
2329 0 : g[d] = m_vDPolynom[ind[d]].value(x[d]);
2330 :
2331 : // multiply by all functions not depending on x[d]
2332 0 : for(int d2 = 0; d2 < dim; ++d2)
2333 : {
2334 : // skip own value
2335 0 : if(d2 == d) continue;
2336 :
2337 0 : g[d] *= m_vPolynom[ind[d2]].value(x[d2]);
2338 : }
2339 : }
2340 0 : }
2341 :
2342 : /// return Multi index for index i
2343 : inline const MathVector<dim,int>& multi_index(size_t i) const
2344 : {
2345 : check_index(i);
2346 : return m_vMultiIndex[i];
2347 : }
2348 :
2349 : /// return the index for a multi_index
2350 0 : inline size_t index(const MathVector<dim,int>& ind) const
2351 : {
2352 : check_multi_index(ind);
2353 0 : for(size_t i=0; i<nsh; ++i)
2354 0 : if(multi_index(i) == ind) return i;
2355 0 : UG_THROW("Index not found in LagrangeLSFS");
2356 : }
2357 :
2358 : /// return the index for a multi_index
2359 : inline size_t mapped_index(const MathVector<dim,int>& ind) const
2360 : {
2361 : check_multi_index(ind);
2362 :
2363 : return ind[2] * (p+1)*(p+1) + ind[1] * (p+1) + ind[0];
2364 : }
2365 :
2366 : /// return the multi_index for an index
2367 : inline MathVector<dim,int> mapped_multi_index(size_t i) const
2368 : {
2369 : check_index(i);
2370 :
2371 : return MathVector<dim,int>( i%(p+1), i/(p+1)%(p+1), i/((p+1)*(p+1)));
2372 : }
2373 :
2374 : /// checks in debug mode that index is valid
2375 : inline void check_index(size_t i) const
2376 : {
2377 : UG_ASSERT(i < nsh, "Wrong index.");
2378 : }
2379 :
2380 : /// checks in debug mode that multi-index is valid
2381 : inline void check_multi_index(const MathVector<dim,int>& ind) const
2382 : {
2383 : UG_ASSERT(ind[0] <= (int)p && ind[0] >= 0, "Wrong Multiindex.");
2384 : UG_ASSERT(ind[1] <= (int)p && ind[1] >= 0, "Wrong Multiindex.");
2385 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
2386 : }
2387 :
2388 : private:
2389 : /// order
2390 : size_t p;
2391 :
2392 : /// Number of shape functions
2393 : size_t nsh;
2394 :
2395 : std::vector<Polynomial1D> m_vPolynom; ///< Shape Polynomials
2396 : std::vector<Polynomial1D> m_vDPolynom; ///< Derivative of Shape Polynomial
2397 :
2398 : std::vector<MathVector<dim,int> > m_vMultiIndex;
2399 : };
2400 :
2401 : ///////////////////////////////////////////////////////////////////////////////
2402 : // Octahedron
2403 : ///////////////////////////////////////////////////////////////////////////////
2404 :
2405 : // NOTE: Currently only 1st order is implemented. There is no shape function
2406 : // set for octahedrons, that is continuous and allows a continuous
2407 : // derivative in the inner of the pyramid. This is basically, since
2408 : // one regards the octahedron as 4 tetrahedrons, glued together.
2409 : template <int TOrder>
2410 : class LagrangeLSFS<ReferenceOctahedron, TOrder>
2411 : : public LagrangeLDS<ReferenceOctahedron>,
2412 : public BaseLSFS<LagrangeLSFS<ReferenceOctahedron, TOrder>, 3>
2413 : {
2414 : private:
2415 : /// abbreviation for order
2416 : static const size_t p = TOrder;
2417 :
2418 : /// base class
2419 : typedef BaseLSFS<LagrangeLSFS<ReferenceOctahedron, TOrder>, 3> base_type;
2420 :
2421 : public:
2422 : /// Shape type
2423 : typedef typename base_type::shape_type shape_type;
2424 :
2425 : /// Gradient type
2426 : typedef typename base_type::grad_type grad_type;
2427 :
2428 : public:
2429 : /// Order of Shape functions
2430 : static const size_t order = TOrder;
2431 :
2432 : /// Dimension, where shape functions are defined
2433 : static const int dim = 3; //reference_element_type::dim; (compile error on OSX 10.5)
2434 :
2435 : /// Number of shape functions
2436 : static const size_t nsh = 6;
2437 :
2438 : public:
2439 : /// Constructor
2440 : LagrangeLSFS();
2441 :
2442 : /// \copydoc ug::LocalShapeFunctionSet::continuous()
2443 0 : inline bool continuous() const {return true;}
2444 :
2445 : /// \copydoc ug::LocalShapeFunctionSet::num_sh()
2446 0 : inline size_t num_sh() const {return nsh;}
2447 :
2448 : /// \copydoc ug::LocalShapeFunctionSet::position()
2449 0 : inline bool position(size_t i, MathVector<dim>& pos) const
2450 : {
2451 : // get Multi Index
2452 : MathVector<dim,int> ind = multi_index(i);
2453 :
2454 : // set position
2455 0 : for(int d = 0; d < dim; ++d)
2456 0 : pos[d] = EquidistantLagrange1D::position(ind[d], p);
2457 :
2458 0 : return true;
2459 : }
2460 :
2461 : /// \copydoc ug::LocalShapeFunctionSet::shape()
2462 0 : inline number shape(const size_t i, const MathVector<dim>& x) const
2463 : {
2464 : // only first order
2465 : if(p != 1) UG_THROW("Only 1. order Lagrange Octahedron implemented.");
2466 :
2467 : // shape analogously to pyramidal case introducing additional distinction of cases
2468 : // z >= 0 and z < 0
2469 0 : const number z_sgn = (x[2] < 0) ? -1.0 : 1.0;
2470 :
2471 0 : switch(i)
2472 : {
2473 : case 0 :
2474 0 : if (x[2] < 0)
2475 0 : return(-1.0*x[2]);
2476 : else
2477 : return(0.0);
2478 : case 1 :
2479 0 : if (x[0] > x[1])
2480 0 : return((1.0-x[0])*(1.0-x[1]) + z_sgn*x[2]*(x[1]-1.0));
2481 : else
2482 0 : return((1.0-x[0])*(1.0-x[1]) + z_sgn*x[2]*(x[0]-1.0));
2483 : case 2 :
2484 0 : if (x[0] > x[1])
2485 0 : return(x[0]*(1.0-x[1]) - z_sgn*x[2]*x[1]);
2486 : else
2487 0 : return(x[0]*(1.0-x[1]) - z_sgn*x[2]*x[0]);
2488 : case 3 :
2489 0 : if (x[0] > x[1])
2490 0 : return(x[0]*x[1] + z_sgn*x[2]*x[1]);
2491 : else
2492 0 : return(x[0]*x[1] + z_sgn*x[2]*x[0]);
2493 : case 4 :
2494 0 : if (x[0] > x[1])
2495 0 : return((1.0-x[0])*x[1] - z_sgn*x[2]*x[1]);
2496 : else
2497 0 : return((1.0-x[0])*x[1] - z_sgn*x[2]*x[0]);
2498 : case 5 :
2499 0 : if (x[2] < 0)
2500 0 : return(0.0);
2501 : else
2502 : return(x[2]);
2503 0 : default: UG_THROW("Wrong index "<< i<<" in Octahedron.");
2504 : }
2505 : }
2506 :
2507 : /// shape value for a Multi Index
2508 0 : inline number shape(const MathVector<dim,int>& ind, const MathVector<dim>& x) const
2509 : {
2510 : check_multi_index(ind);
2511 :
2512 : // forward
2513 0 : return shape(index(ind), x);
2514 : }
2515 :
2516 : /// \copydoc ug::LocalShapeFunctionSet::grad()
2517 0 : inline void grad(grad_type& g, const size_t i, const MathVector<dim>& x) const
2518 : {
2519 : // only first order
2520 : if(p != 1) UG_THROW("Only 1. order Lagrange Octahedron implemented.");
2521 :
2522 : // shape analogously to pyramidal case introducing additional distinction of cases
2523 : // z >= 0 and z < 0
2524 0 : const number z_sgn = (x[2] < 0) ? -1.0 : 1.0;
2525 :
2526 0 : switch(i)
2527 : {
2528 : case 0:
2529 0 : if (x[2] < 0.0)
2530 : {
2531 0 : g[0] = 0.0;
2532 0 : g[1] = 0.0;
2533 0 : g[2] = -1.0;
2534 0 : break;
2535 : }
2536 : else
2537 : {
2538 0 : g[0] = 0.0;
2539 0 : g[1] = 0.0;
2540 0 : g[2] = 0.0;
2541 0 : break;
2542 : }
2543 : case 1:
2544 0 : if (x[0] > x[1])
2545 : {
2546 0 : g[0] = -(1.0-x[1]);
2547 0 : g[1] = -(1.0-x[0]) + z_sgn*x[2];
2548 0 : g[2] = -z_sgn*(1.0-x[1]);
2549 0 : break;
2550 : }
2551 : else
2552 : {
2553 0 : g[0] = -(1.0-x[1]) + z_sgn*x[2];
2554 0 : g[1] = -(1.0-x[0]);
2555 0 : g[2] = -z_sgn*(1.0-x[0]);
2556 0 : break;
2557 : }
2558 : case 2:
2559 0 : if (x[0] > x[1])
2560 : {
2561 0 : g[0] = (1.0-x[1]);
2562 0 : g[1] = -x[0] - z_sgn*x[2];
2563 0 : g[2] = -z_sgn*x[1];
2564 0 : break;
2565 : }
2566 : else
2567 : {
2568 0 : g[0] = (1.0-x[1]) - z_sgn*x[2];
2569 0 : g[1] = -x[0];
2570 0 : g[2] = -z_sgn*x[0];
2571 0 : break;
2572 : }
2573 : case 3:
2574 0 : if (x[0] > x[1])
2575 : {
2576 0 : g[0] = x[1];
2577 0 : g[1] = x[0] + z_sgn*x[2];
2578 0 : g[2] = z_sgn*x[1];
2579 0 : break;
2580 : }
2581 : else
2582 : {
2583 0 : g[0] = x[1] + z_sgn*x[2];
2584 0 : g[1] = x[0];
2585 0 : g[2] = z_sgn*x[0];
2586 0 : break;
2587 : }
2588 : case 4:
2589 0 : if (x[0] > x[1])
2590 : {
2591 0 : g[0] = -x[1];
2592 0 : g[1] = 1.0-x[0] - z_sgn*x[2];
2593 0 : g[2] = -z_sgn*x[1];
2594 0 : break;
2595 : }
2596 : else
2597 : {
2598 0 : g[0] = -x[1] - z_sgn*x[2];
2599 0 : g[1] = 1.0-x[0];
2600 0 : g[2] = -z_sgn*x[0];
2601 0 : break;
2602 : }
2603 : case 5:
2604 0 : if (x[2] < 0.0)
2605 : {
2606 0 : g[0] = 0.0;
2607 0 : g[1] = 0.0;
2608 0 : g[2] = 0.0;
2609 0 : break;
2610 : }
2611 : else
2612 : {
2613 0 : g[0] = 0.0;
2614 0 : g[1] = 0.0;
2615 0 : g[2] = 1.0;
2616 0 : break;
2617 : }
2618 0 : default: UG_THROW("Wrong index "<< i<<" in Octahedron.");
2619 : }
2620 0 : }
2621 :
2622 : /// evaluates the gradient
2623 0 : inline void grad(grad_type& g, const MathVector<dim,int> ind,
2624 : const MathVector<dim>& x) const
2625 : {
2626 0 : grad(g, index(ind), x);
2627 0 : }
2628 :
2629 : /// return Multi index for index i
2630 0 : inline const MathVector<dim,int>& multi_index(size_t i) const
2631 : {
2632 : check_index(i);
2633 0 : return m_vMultiIndex[i];
2634 : }
2635 :
2636 : /// return the index for a multi_index
2637 0 : inline size_t index(const MathVector<dim,int>& ind) const
2638 : {
2639 : check_multi_index(ind);
2640 0 : for(size_t i=0; i<nsh; ++i)
2641 0 : if(multi_index(i) == ind) return i;
2642 0 : UG_THROW("Index not found in LagrangeLSFS");
2643 : }
2644 :
2645 : /// checks in debug mode that index is valid
2646 0 : inline void check_index(size_t i) const
2647 : {
2648 : UG_ASSERT(i < nsh, "Wrong index.");
2649 0 : }
2650 :
2651 : /// checks in debug mode that multi-index is valid
2652 0 : inline void check_multi_index(const MathVector<dim,int>& ind) const
2653 : {
2654 : UG_ASSERT(ind[0] <= (int)p-ind[2] && ind[0] >= 0, "Wrong Multiindex.");
2655 : UG_ASSERT(ind[1] <= (int)p-ind[2] && ind[1] >= 0, "Wrong Multiindex.");
2656 : UG_ASSERT(ind[2] <= (int)p && ind[2] >= 0, "Wrong Multiindex.");
2657 0 : }
2658 :
2659 : private:
2660 :
2661 : MathVector<dim,int> m_vMultiIndex[nsh];
2662 : };
2663 :
2664 :
2665 : } //namespace ug
2666 :
2667 : #endif /* __H__UG__LIB_DISC__LOCAL_SHAPE_FUNCTION_SET__LAGRANGE__LAGRANGE__ */
2668 :
|
