Line data Source code
1 : /*
2 : * Copyright (c) 2011-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 : #include "lagrange.h"
34 : #include <sstream>
35 : #include "common/util/provider.h"
36 :
37 : namespace ug{
38 :
39 : template <typename TRefElem>
40 : void SetLagrangeVertexMultiIndex( MathVector<TRefElem::dim,int>* vMultiIndex,
41 : const TRefElem& rRef,
42 : size_t p,
43 : size_t& index)
44 : {
45 : // dimension of Reference element
46 : static const int dim = TRefElem::dim;
47 :
48 : // get corner position integer
49 : const MathVector<dim,int>* vCo = rRef.corner();
50 :
51 : // loop all vertices
52 105 : for(size_t i = 0; i< rRef.num(0); ++i)
53 : {
54 : // set multi index
55 299 : for(int d = 0; d<dim; ++d)
56 : {
57 : UG_ASSERT(((long)p)*vCo[i][d] >= 0,
58 : "Wrong multi index m["<<d<<"]="<<p*vCo[i][d]);
59 219 : vMultiIndex[index][d] = p*vCo[i][d];
60 : }
61 :
62 : // next index
63 86 : ++index;
64 : }
65 : }
66 :
67 : template <typename TRefElem>
68 19 : void SetLagrangeEdgeMultiIndex( MathVector<TRefElem::dim,int>* vMultiIndex,
69 : const TRefElem& rRef,
70 : size_t p,
71 : size_t& index)
72 : {
73 : // dimension of Reference element
74 : static const int dim = TRefElem::dim;
75 :
76 : // only for 2d,3d elems we do something
77 : if(dim < 1) return;
78 :
79 : // get corner position integer
80 : const MathVector<dim,int>* vCo = rRef.corner();
81 :
82 : // loop all edges
83 132 : for(size_t e = 0; e< rRef.num(1); ++e)
84 : {
85 : // get ids of corners of edge
86 113 : const size_t co0 = rRef.id(1,e, 0,0);
87 113 : const size_t co1 = rRef.id(1,e, 0,1);
88 :
89 : // add dofs on the edge
90 218 : for(size_t i = 1; i < p; ++i)
91 : {
92 : // compute multi index
93 : MathVector<dim,int> m;
94 105 : VecScaleAdd(m, i, vCo[co1], (p-i), vCo[co0]);
95 :
96 : // set multi index
97 393 : for(int d = 0; d<dim; ++d)
98 : {
99 : UG_ASSERT(m[d] >= 0, "Wrong multi index m["<<d<<"]="<<m[d]);
100 288 : vMultiIndex[index][d] = m[d];
101 : }
102 :
103 : // next index
104 105 : index++;
105 : }
106 : }
107 : }
108 :
109 : template <typename TRefElem>
110 16 : void SetLagrangeFaceMultiIndex( MathVector<TRefElem::dim,int>* vMultiIndex,
111 : const TRefElem& rRef,
112 : size_t p,
113 : size_t& index)
114 : {
115 : // dimension of Reference element
116 : static const int dim = TRefElem::dim;
117 :
118 : // only for 2d,3d elems we do something
119 : if(dim < 2) return;
120 :
121 : // get corner position integer
122 : const MathVector<dim,int>* vCo = rRef.corner();
123 :
124 : // add dof on quadrilateral
125 72 : for(size_t f = 0; f< rRef.num(2); ++f)
126 : {
127 : // get ids of corners of edge
128 : const size_t co0 = rRef.id(2,f, 0,0);
129 56 : const size_t co1 = rRef.id(2,f, 0,1);
130 : size_t co2 = rRef.id(2,f, 0,2);
131 56 : if(rRef.num(2,f,0)==4) co2 = rRef.id(2,f, 0,3);
132 :
133 : // directions of counting
134 : MathVector<dim,int> d1, d2;
135 56 : VecSubtract(d1, vCo[co1], vCo[co0]);
136 56 : VecSubtract(d2, vCo[co2], vCo[co0]);
137 :
138 : // loop 'y'-direction
139 107 : for(size_t j = 1; j < p; ++j)
140 : {
141 : // for a quadrilateral we have a quadratic loop, but for a
142 : // triangle we need to stop at the diagonal
143 51 : const size_t off = ((rRef.num(2,f,0)==3) ? j : 0);
144 :
145 : // loop 'x'-direction
146 108 : for(size_t i = 1; i < p-off; ++i)
147 : {
148 : // compute multi index
149 57 : MathVector<dim,int> m = vCo[co0]; m *= p;
150 57 : VecScaleAppend(m, i, d1, j, d2);
151 :
152 : // set multi index
153 222 : for(int d = 0; d<dim; ++d)
154 : {
155 : UG_ASSERT(m[d] >= 0, "Wrong multi index m["<<d<<"]="<<m[d]);
156 165 : vMultiIndex[index][d] = m[d];
157 : }
158 :
159 : // next index
160 57 : ++index;
161 : }
162 : }
163 : }
164 : }
165 :
166 : template <typename TRefElem>
167 10 : void SetLagrangeVolumeMultiIndex( MathVector<TRefElem::dim,int>* vMultiIndex,
168 : const TRefElem& rRef,
169 : size_t p,
170 : size_t& index)
171 : {
172 : // dimension of Reference element
173 : static const int dim = TRefElem::dim;
174 :
175 : // only for 3d elems we do something
176 : if(dim < 3) return;
177 :
178 : // get corner position integer
179 : // const MathVector<dim,int>* vCo = rRef.corner();
180 :
181 : // get type of reference element
182 : ReferenceObjectID type = rRef.roid(dim, 0);
183 :
184 : // handle elems
185 10 : switch(type)
186 : {
187 : case ROID_TETRAHEDRON:
188 6 : for(size_t m2 = 1; m2 < p; ++m2)
189 4 : for(size_t m1 = 1; m1 < p-m2; ++m1)
190 1 : for(size_t m0 = 1; m0 < p-m2-m1; ++m0)
191 : {
192 : // use regular mapping for inner DoFs
193 0 : vMultiIndex[index][0] = m0;
194 0 : vMultiIndex[index][1] = m1;
195 0 : vMultiIndex[index++][2] = m2;
196 : }
197 : break;
198 :
199 1 : case ROID_PYRAMID:
200 1 : if(p>1) UG_THROW("LagrangeLSFS: Higher order Pyramid not implemented.");
201 : break;
202 :
203 0 : case ROID_OCTAHEDRON:
204 0 : if(p>1) UG_THROW("LagrangeLSFS: Higher order Octahedron not implemented.");
205 : break;
206 :
207 : case ROID_PRISM:
208 6 : for(size_t m2 = 1; m2 < p; ++m2)
209 8 : for(size_t m1 = 1; m1 < p; ++m1)
210 7 : for(size_t m0 = 1; m0 < p-m1; ++m0)
211 : {
212 : // use regular mapping for inner DoFs
213 2 : vMultiIndex[index][0] = m0;
214 2 : vMultiIndex[index][1] = m1;
215 2 : vMultiIndex[index++][2] = m2;
216 : }
217 : break;
218 :
219 : case ROID_HEXAHEDRON:
220 6 : for(size_t m2 = 1; m2 < p; ++m2)
221 8 : for(size_t m1 = 1; m1 < p; ++m1)
222 14 : for(size_t m0 = 1; m0 < p; ++m0)
223 : {
224 : // use regular mapping for inner DoFs
225 9 : vMultiIndex[index][0] = m0;
226 9 : vMultiIndex[index][1] = m1;
227 9 : vMultiIndex[index++][2] = m2;
228 : }
229 : break;
230 0 : default: UG_THROW("LagrangeLSFS: Missing 3d mapping for type '"<<type<<"'."
231 : " roid="<<rRef.reference_object_id());
232 : }
233 : }
234 :
235 : template <typename TRefElem>
236 19 : void SetLagrangeMultiIndex( MathVector<TRefElem::dim,int>* vMultiIndex,
237 : const TRefElem& rRef,
238 : size_t p)
239 : {
240 : // init shape -> multi-index mapping
241 19 : size_t index = 0;
242 :
243 : // vertices
244 : SetLagrangeVertexMultiIndex(vMultiIndex, rRef, p, index);
245 :
246 : // edge
247 19 : SetLagrangeEdgeMultiIndex(vMultiIndex, rRef, p, index);
248 :
249 : // faces
250 16 : SetLagrangeFaceMultiIndex(vMultiIndex, rRef, p, index);
251 :
252 : // volumes
253 10 : SetLagrangeVolumeMultiIndex(vMultiIndex, rRef, p, index);
254 19 : }
255 :
256 : ///////////////////////////////////////////////////////////////////////////////
257 : // Edge
258 : ///////////////////////////////////////////////////////////////////////////////
259 :
260 : template <int TOrder>
261 3 : LagrangeLSFS<ReferenceEdge, TOrder>::LagrangeLSFS()
262 30 : : LagrangeLDS<ReferenceEdge>(p)
263 : {
264 : // init polynomials
265 12 : for(size_t i = 0; i < nsh; ++i)
266 : {
267 : // create equidistant polynomials and its derivatives
268 9 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
269 18 : m_vDPolynom[i] = m_vPolynom[i].derivative();
270 : }
271 :
272 : // reference element
273 3 : const ReferenceEdge& rRef = Provider<ReferenceEdge>::get();
274 :
275 : // init shape -> multi-index mapping
276 3 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
277 3 : }
278 :
279 0 : void FlexLagrangeLSFS<ReferenceEdge>::set_order(size_t order)
280 : {
281 0 : LagrangeLDS<ReferenceEdge>::set_order(order);
282 :
283 : // resize
284 0 : p = order;
285 0 : nsh = p+1;
286 0 : m_vPolynom.resize(nsh);
287 0 : m_vDPolynom.resize(nsh);
288 0 : m_vMultiIndex.resize(nsh);
289 :
290 : // init polynomials
291 0 : for(size_t i = 0; i < nsh; ++i)
292 : {
293 : // create equidistant polynomials and its derivatives
294 0 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
295 0 : m_vDPolynom[i] = m_vPolynom[i].derivative();
296 : }
297 :
298 : // reference element
299 0 : const ReferenceEdge& rRef = Provider<ReferenceEdge>::get();
300 :
301 : // init shape -> multi-index mapping
302 0 : SetLagrangeMultiIndex(&m_vMultiIndex[0], rRef, p);
303 0 : }
304 :
305 : template class LagrangeLSFS<ReferenceEdge, 1>;
306 : template class LagrangeLSFS<ReferenceEdge, 2>;
307 : template class LagrangeLSFS<ReferenceEdge, 3>;
308 : template class LagrangeLSFS<ReferenceEdge, 4>;
309 : template class LagrangeLSFS<ReferenceEdge, 5>;
310 :
311 : //template class FlexLagrangeLSFS<ReferenceEdge>;
312 :
313 : ///////////////////////////////////////////////////////////////////////////////
314 : // Triangle
315 : ///////////////////////////////////////////////////////////////////////////////
316 :
317 : template <int TOrder>
318 3 : LagrangeLSFS<ReferenceTriangle, TOrder>::LagrangeLSFS()
319 40 : : LagrangeLDS<ReferenceTriangle>(p)
320 : {
321 : // init polynomials
322 12 : for(size_t i = 0; i <= p; ++i)
323 : {
324 : // create trancated equidistant polynomials and its derivatives
325 9 : m_vPolynom[i] = TruncatedEquidistantLagrange1D(i, p);
326 18 : m_vDPolynom[i] = m_vPolynom[i].derivative();
327 : }
328 :
329 : // reference element
330 3 : const ReferenceTriangle& rRef = Provider<ReferenceTriangle>::get();
331 :
332 : // init shape -> multi-index mapping
333 3 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
334 3 : }
335 :
336 0 : void FlexLagrangeLSFS<ReferenceTriangle>::set_order(size_t order)
337 : {
338 0 : LagrangeLDS<ReferenceTriangle>::set_order(order);
339 :
340 : // resize
341 0 : p = order;
342 0 : nsh = BinomCoeff(dim+p, p);
343 :
344 0 : const size_t polSize = p+1;
345 0 : m_vPolynom.resize(polSize);
346 0 : m_vDPolynom.resize(polSize);
347 0 : m_vMultiIndex.resize(nsh);
348 :
349 : // init polynomials
350 0 : for(size_t i = 0; i <= p; ++i)
351 : {
352 : // create trancated equidistant polynomials and its derivatives
353 0 : m_vPolynom[i] = TruncatedEquidistantLagrange1D(i, p);
354 0 : m_vDPolynom[i] = m_vPolynom[i].derivative();
355 : }
356 :
357 : // reference element
358 0 : const ReferenceTriangle& rRef = Provider<ReferenceTriangle>::get();
359 :
360 : // init shape -> multi-index mapping
361 0 : SetLagrangeMultiIndex(&m_vMultiIndex[0], rRef, p);
362 0 : }
363 :
364 : template class LagrangeLSFS<ReferenceTriangle, 1>;
365 : template class LagrangeLSFS<ReferenceTriangle, 2>;
366 : template class LagrangeLSFS<ReferenceTriangle, 3>;
367 : template class LagrangeLSFS<ReferenceTriangle, 4>;
368 : template class LagrangeLSFS<ReferenceTriangle, 5>;
369 :
370 : //template class FlexLagrangeLSFS<ReferenceTriangle>;
371 :
372 : ///////////////////////////////////////////////////////////////////////////////
373 : // Quadrilateral
374 : ///////////////////////////////////////////////////////////////////////////////
375 :
376 : template <int TOrder>
377 3 : LagrangeLSFS<ReferenceQuadrilateral, TOrder>::LagrangeLSFS()
378 50 : : LagrangeLDS<ReferenceQuadrilateral>(p)
379 : {
380 : // init polynomials
381 12 : for(size_t i = 0; i <= p; ++i)
382 : {
383 : // create truncated equidistant polynomials and its derivatives
384 9 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
385 18 : m_vDPolynom[i] = m_vPolynom[i].derivative();
386 : }
387 :
388 : // reference element
389 3 : const ReferenceQuadrilateral& rRef = Provider<ReferenceQuadrilateral>::get();
390 :
391 : // init shape -> multi-index mapping
392 3 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
393 3 : }
394 :
395 0 : void FlexLagrangeLSFS<ReferenceQuadrilateral>::set_order(size_t order)
396 : {
397 0 : LagrangeLDS<ReferenceQuadrilateral>::set_order(order);
398 :
399 : // resize
400 0 : p = order;
401 0 : nsh = (p+1)*(p+1);
402 :
403 : const size_t polSize = p+1;
404 0 : m_vPolynom.resize(polSize);
405 0 : m_vDPolynom.resize(polSize);
406 0 : m_vMultiIndex.resize(nsh);
407 :
408 : // init polynomials
409 0 : for(size_t i = 0; i <= p; ++i)
410 : {
411 : // create truncated equidistant polynomials and its derivatives
412 0 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
413 0 : m_vDPolynom[i] = m_vPolynom[i].derivative();
414 : }
415 :
416 : // reference element
417 0 : const ReferenceQuadrilateral& rRef = Provider<ReferenceQuadrilateral>::get();
418 :
419 : // init shape -> multi-index mapping
420 0 : SetLagrangeMultiIndex(&m_vMultiIndex[0], rRef, p);
421 0 : }
422 :
423 : template class LagrangeLSFS<ReferenceQuadrilateral, 1>;
424 : template class LagrangeLSFS<ReferenceQuadrilateral, 2>;
425 : template class LagrangeLSFS<ReferenceQuadrilateral, 3>;
426 : template class LagrangeLSFS<ReferenceQuadrilateral, 4>;
427 : template class LagrangeLSFS<ReferenceQuadrilateral, 5>;
428 :
429 : //template class FlexLagrangeLSFS<ReferenceQuadrilateral>;
430 :
431 : ///////////////////////////////////////////////////////////////////////////////
432 : // Tetrahedron
433 : ///////////////////////////////////////////////////////////////////////////////
434 :
435 : template <int TOrder>
436 3 : LagrangeLSFS<ReferenceTetrahedron, TOrder>::LagrangeLSFS()
437 55 : : LagrangeLDS<ReferenceTetrahedron>(p)
438 : {
439 12 : for(size_t i = 0; i <= p; ++i)
440 : {
441 : // create trancated equidistant polynomials and its derivatives
442 9 : m_vPolynom[i] = TruncatedEquidistantLagrange1D(i, p);
443 18 : m_vDPolynom[i] = m_vPolynom[i].derivative();
444 : }
445 :
446 : // reference element
447 3 : const ReferenceTetrahedron& rRef = Provider<ReferenceTetrahedron>::get();
448 :
449 : // init shape -> multi-index mapping
450 3 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
451 3 : }
452 :
453 0 : void FlexLagrangeLSFS<ReferenceTetrahedron>::set_order(size_t order)
454 : {
455 0 : LagrangeLDS<ReferenceTetrahedron>::set_order(order);
456 :
457 : // resize
458 0 : p = order;
459 0 : nsh = BinomCoeff(dim + p, p);
460 :
461 0 : const size_t polSize = p+1;
462 0 : m_vPolynom.resize(polSize);
463 0 : m_vDPolynom.resize(polSize);
464 0 : m_vMultiIndex.resize(nsh);
465 :
466 : // init polynomials
467 0 : for(size_t i = 0; i <= p; ++i)
468 : {
469 : // create trancated equidistant polynomials and its derivatives
470 0 : m_vPolynom[i] = TruncatedEquidistantLagrange1D(i, p);
471 0 : m_vDPolynom[i] = m_vPolynom[i].derivative();
472 : }
473 :
474 : // reference element
475 0 : const ReferenceTetrahedron& rRef = Provider<ReferenceTetrahedron>::get();
476 :
477 : // init shape -> multi-index mapping
478 0 : SetLagrangeMultiIndex(&m_vMultiIndex[0], rRef, p);
479 0 : }
480 :
481 : template class LagrangeLSFS<ReferenceTetrahedron, 1>;
482 : template class LagrangeLSFS<ReferenceTetrahedron, 2>;
483 : template class LagrangeLSFS<ReferenceTetrahedron, 3>;
484 : template class LagrangeLSFS<ReferenceTetrahedron, 4>;
485 : template class LagrangeLSFS<ReferenceTetrahedron, 5>;
486 :
487 : //template class FlexLagrangeLSFS<ReferenceTetrahedron>;
488 :
489 : ///////////////////////////////////////////////////////////////////////////////
490 : // Prism
491 : ///////////////////////////////////////////////////////////////////////////////
492 :
493 : template <int TOrder>
494 3 : LagrangeLSFS<ReferencePrism, TOrder>::LagrangeLSFS()
495 103 : : LagrangeLDS<ReferencePrism>(p)
496 : {
497 12 : for(size_t i = 0; i <= p; ++i)
498 : {
499 : // create truncated equidistant polynomials and its derivatives
500 9 : m_vTruncPolynom[i] = TruncatedEquidistantLagrange1D(i, p);
501 9 : m_vDTruncPolynom[i] = m_vTruncPolynom[i].derivative();
502 :
503 : // create equidistant polynomials and its derivatives
504 9 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
505 18 : m_vDPolynom[i] = m_vPolynom[i].derivative();
506 : }
507 :
508 : // reference element
509 3 : const ReferencePrism& rRef = Provider<ReferencePrism>::get();
510 :
511 : // init shape -> multi-index mapping
512 3 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
513 3 : }
514 :
515 0 : void FlexLagrangeLSFS<ReferencePrism>::set_order(size_t order)
516 : {
517 0 : LagrangeLDS<ReferencePrism>::set_order(order);
518 :
519 : // resize
520 0 : p = order;
521 0 : dofPerLayer = BinomCoeff(2 + p, p);
522 0 : nsh = dofPerLayer * (p+1);
523 :
524 : const size_t polSize = p+1;
525 0 : m_vPolynom.resize(polSize);
526 0 : m_vDPolynom.resize(polSize);
527 0 : m_vTruncPolynom.resize(polSize);
528 0 : m_vDTruncPolynom.resize(polSize);
529 0 : m_vMultiIndex.resize(nsh);
530 :
531 : // init polynomials
532 0 : for(size_t i = 0; i <= p; ++i)
533 : {
534 : // create truncated equidistant polynomials and its derivatives
535 0 : m_vTruncPolynom[i] = TruncatedEquidistantLagrange1D(i, p);
536 0 : m_vDTruncPolynom[i] = m_vTruncPolynom[i].derivative();
537 :
538 : // create equidistant polynomials and its derivatives
539 0 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
540 0 : m_vDPolynom[i] = m_vPolynom[i].derivative();
541 : }
542 :
543 : // reference element
544 0 : const ReferencePrism& rRef = Provider<ReferencePrism>::get();
545 :
546 : // init shape -> multi-index mapping
547 0 : SetLagrangeMultiIndex(&m_vMultiIndex[0], rRef, p);
548 0 : }
549 :
550 : template class LagrangeLSFS<ReferencePrism, 1>;
551 : template class LagrangeLSFS<ReferencePrism, 2>;
552 : template class LagrangeLSFS<ReferencePrism, 3>;
553 : template class LagrangeLSFS<ReferencePrism, 4>;
554 : template class LagrangeLSFS<ReferencePrism, 5>;
555 :
556 : //template class FlexLagrangeLSFS<ReferencePrism>;
557 :
558 : ///////////////////////////////////////////////////////////////////////////////
559 : // Pyramid
560 : ///////////////////////////////////////////////////////////////////////////////
561 :
562 : template <int TOrder>
563 1 : LagrangeLSFS<ReferencePyramid, TOrder>::LagrangeLSFS()
564 6 : : LagrangeLDS<ReferencePyramid>(p)
565 : {
566 1 : m_vvPolynom.resize(p+1);
567 1 : m_vvDPolynom.resize(p+1);
568 :
569 3 : for(size_t i2 = 0; i2 <= p; ++i2)
570 : {
571 2 : m_vvPolynom[i2].resize(p+1);
572 2 : m_vvDPolynom[i2].resize(p+1);
573 :
574 5 : for(size_t i = 0; i <= p-i2; ++i)
575 : {
576 6 : m_vvPolynom[i2][i] = BoundedEquidistantLagrange1D(i, p, p-i2);
577 6 : m_vvDPolynom[i2][i] = m_vvPolynom[i2][i].derivative();
578 : }
579 : }
580 :
581 : // reference element
582 : const ReferencePyramid& rRef =
583 1 : Provider<ReferencePyramid>::get();
584 :
585 : // init shape -> multi-index mapping
586 1 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
587 1 : }
588 :
589 : template class LagrangeLSFS<ReferencePyramid, 1>;
590 : template class LagrangeLSFS<ReferencePyramid, 2>;
591 :
592 : ///////////////////////////////////////////////////////////////////////////////
593 : // Octahedron
594 : ///////////////////////////////////////////////////////////////////////////////
595 :
596 : template <int TOrder>
597 0 : LagrangeLSFS<ReferenceOctahedron, TOrder>::LagrangeLSFS()
598 0 : : LagrangeLDS<ReferenceOctahedron>(p)
599 : {
600 : //UG_THROW("LagrangeLSFS<ReferenceOctahedron, TOrder>::LagrangeLSFS(): Octahedral elements currently not implemented. Use LagrangeP1 implementation instead.");
601 :
602 : if(p != 1)
603 : UG_THROW("LagrangeLSFS<ReferenceOctahedron, TOrder>::LagrangeLSFS(): Octahedral elements only implemented for order p = 1.");
604 :
605 : // reference element
606 : const ReferenceOctahedron& rRef =
607 0 : Provider<ReferenceOctahedron>::get();
608 :
609 : // init shape -> multi-index mapping
610 0 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
611 0 : }
612 :
613 : template class LagrangeLSFS<ReferenceOctahedron, 1>;
614 :
615 : ///////////////////////////////////////////////////////////////////////////////
616 : // Hexahedron
617 : ///////////////////////////////////////////////////////////////////////////////
618 :
619 : template <int TOrder>
620 3 : LagrangeLSFS<ReferenceHexahedron, TOrder>::LagrangeLSFS()
621 120 : : LagrangeLDS<ReferenceHexahedron>(p)
622 : {
623 : // init polynomials
624 12 : for(size_t i = 0; i <= p; ++i)
625 : {
626 : // create trancated equidistant polynomials and its derivatives
627 9 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
628 18 : m_vDPolynom[i] = m_vPolynom[i].derivative();
629 : }
630 :
631 : // reference element
632 3 : const ReferenceHexahedron& rRef = Provider<ReferenceHexahedron>::get();
633 :
634 : // init shape -> multi-index mapping
635 3 : SetLagrangeMultiIndex(m_vMultiIndex, rRef, p);
636 3 : }
637 :
638 0 : void FlexLagrangeLSFS<ReferenceHexahedron>::set_order(size_t order)
639 : {
640 0 : LagrangeLDS<ReferenceHexahedron>::set_order(order);
641 :
642 : // resize
643 0 : p = order;
644 0 : nsh = (p+1) * (p+1) * (p+1);
645 :
646 : const size_t polSize = p+1;
647 0 : m_vPolynom.resize(polSize);
648 0 : m_vDPolynom.resize(polSize);
649 0 : m_vMultiIndex.resize(nsh);
650 :
651 : // init polynomials
652 0 : for(size_t i = 0; i <= p; ++i)
653 : {
654 : // create trancated equidistant polynomials and its derivatives
655 0 : m_vPolynom[i] = EquidistantLagrange1D(i, p);
656 0 : m_vDPolynom[i] = m_vPolynom[i].derivative();
657 : }
658 :
659 : // reference element
660 0 : const ReferenceHexahedron& rRef = Provider<ReferenceHexahedron>::get();
661 :
662 : // init shape -> multi-index mapping
663 0 : SetLagrangeMultiIndex(&m_vMultiIndex[0], rRef, p);
664 0 : }
665 :
666 : template class LagrangeLSFS<ReferenceHexahedron, 1>;
667 : template class LagrangeLSFS<ReferenceHexahedron, 2>;
668 : template class LagrangeLSFS<ReferenceHexahedron, 3>;
669 : template class LagrangeLSFS<ReferenceHexahedron, 4>;
670 : template class LagrangeLSFS<ReferenceHexahedron, 5>;
671 :
672 : } // end namespace ug
|
