1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file LagrangeTriangle.cpp
19 \authors D. Martin, N. Kielbasiewicz, E. Lunéville
20 \since 17 dec 2002
21 \date 12 feb 2015
22
23 \brief Implementation of xlifepp::LagrangeTriangle class members and related functions
24 */
25
26 #include "LagrangeTriangle.hpp"
27 #include "../Interpolation.hpp"
28 #include "utils.h"
29
30 #include <iostream>
31
32 namespace xlifepp
33 {
34
35 //! LagrangeTriangle constructor for Lagrange reference elements
LagrangeTriangle(const Interpolation * interp_p)36 LagrangeTriangle::LagrangeTriangle(const Interpolation* interp_p)
37 : RefTriangle(interp_p)
38 {
39 name_ += "_Lagrange";
40 trace_p->push("LagrangeTriangle::LagrangeTriangle (" + name_ + ")");
41
42 // build element interpolation data
43 interpolationData();
44 // local numbering of points on sides
45 sideNumbering();
46 // local numbering of points on sides
47 sideOfSideNumbering();
48 // Reference Element on edges
49 sideRefElement();
50 maxDegree = interpolation_p->numtype;
51
52 trace_p->pop();
53 }
54
~LagrangeTriangle()55 LagrangeTriangle::~LagrangeTriangle() {}
56 /*
57 --------------------------------------------------------------------------------
58 interp defines Reference Element interpolation data
59 --------------------------------------------------------------------------------
60 */
interpolationData()61 void LagrangeTriangle::interpolationData()
62 {
63 trace_p->push("LagrangeTriangle::interpolationData");
64 number_t interpNum = interpolation_p->numtype;
65
66 switch(interpNum)
67 {
68 case _P0: // Lagrange P_0 discontinuous
69 nbDofs_ = nbInternalDofs_ = 1;
70 break;
71 case _P1BubbleP3: // Lagrange P_1 + P_3 bubble function
72 nbDofs_ = 4;
73 nbInternalDofs_ = 1;
74 nbDofsOnVertices_ = 3;
75 break;
76 default: // Lagrange standard triangle of degree > 0
77 nbDofs_ = ((interpNum + 1)*(interpNum + 2))/2;
78 nbInternalDofs_ = ((interpNum - 2)*(interpNum - 1))/2;
79 nbDofsOnVertices_ = 3;
80 nbDofsInSides_ = nbDofs_ - nbDofsOnVertices_ - nbInternalDofs_;
81 break;
82 }
83
84 // Creating Degrees Of Freedom of reference element
85 // in Lagrange standard (H1-conforming) elements all D.o.F on edges are sharable
86 refDofs.reserve(nbDofs_);
87 lagrangeRefDofs(2, nbDofsOnVertices_, nbInternalDofs_, 3, nbDofsInSides_);
88
89 // allocate storage for shape functions
90 shapeValues.resize(*this, 2);
91 nbPts_ = nbDofs_;
92
93 trace_p->pop();
94 }
95
96 //! create standard Lagrange numbering on sides
sideNumbering()97 void LagrangeTriangle::sideNumbering()
98 {
99 number_t interpNum = interpolation_p->numtype;
100 number_t nbSides = geomRefElem_p->nbSides();
101 sideDofNumbers_.resize(nbSides);
102 if(interpNum==0) //P0 interpolation
103 {
104 for(number_t side = 0; side < nbSides; side++)
105 {
106 sideDofNumbers_[side].resize(1);
107 sideDofNumbers_[side][0]=1;
108 }
109 return;
110 }
111
112 if(interpNum > 0)
113 {
114 number_t intN = interpNum;
115 if(interpNum == _P1BubbleP3) { intN = 1; }
116 number_t nbVertices = geomRefElem_p->nbVertices();
117 number_t more;
118 number_t nbVerticesPerSide = geomRefElem_p->sideVertexNumbers()[0].size();
119 number_t nbDofsPerSide = intN + 1;
120
121 for(number_t side = 0; side < nbSides; side++)
122 {
123 sideDofNumbers_[side].resize(nbDofsPerSide);
124 for(number_t j = 0; j < nbVerticesPerSide; j++)
125 {
126 sideDofNumbers_[side][j] = geomRefElem_p->sideVertexNumber(j + 1, side + 1);
127 }
128 if(intN > 1) // order > 1
129 {
130 more = nbVertices + 1;
131 for(number_t k = nbVerticesPerSide; k < nbDofsPerSide; k++, more += nbSides)
132 {
133 sideDofNumbers_[side][k] = side + more;
134 }
135 }
136 }
137 }
138 }
139
140 //! create standard Lagrange numbering on side of sides
sideOfSideNumbering()141 void LagrangeTriangle::sideOfSideNumbering()
142 {
143 number_t interpNum = interpolation_p->numtype;
144
145 if(interpNum > 0)
146 {
147 number_t nbSideOfSides = geomRefElem_p->nbSideOfSides();
148 number_t nbVerticesPerSideOfSide = geomRefElem_p->sideOfSideVertexNumbers()[0].size();
149
150 sideOfSideDofNumbers_.resize(nbSideOfSides);
151 for(number_t i = 0; i < nbSideOfSides; i++)
152 {
153 sideOfSideDofNumbers_[i].resize(nbVerticesPerSideOfSide);
154 for(number_t j = 0; j < nbVerticesPerSideOfSide; j++)
155 {
156 sideOfSideDofNumbers_[i][j] = geomRefElem_p->sideOfSideVertexNumber(j + 1, i + 1);
157 }
158 }
159 }
160 }
161
162 /*! compute shape functions coefficients in basis monoms by solving linear systems n=nbdofs=(k+1)*(k+2)/2
163 | 1 y1 y1^2 ... x1 x1y1 x1y1^2 ... x1^k | | A11 A12 ... A1n | | 1 0 ... 0 |
164 | 1 y2 y2^2 ... x2 x2y2 x2y2^2 ... x2^k | | A21 A22 ... A2n | | 0 1 ... 0 |
165 | ... | | ... | | ... |
166 | 1 yi yi^2 ... xi xiyi xiyi^2 ... xi^k | | Ai1 Ai2 ... Ain | = | 0 ...1...0 |
167 | ... | | ... | | ... |
168 | 1 yn yn^2 ... xn xnyn xnyn^2 ... xn^k | | An1 An2 ... Ann | | 0 0 ... 1 |
169 */
170
171 //! compute shape functions as polynomials using general algorithm (node coordinates has to be set)
172 // fill polynomials space Wk and dWk
computeShapeFunctions()173 void LagrangeTriangle::computeShapeFunctions()
174 {
175 //compute Matrix Lij=pj(ni) for all pi in Vk, ni nodes
176 number_t k=interpolation_p->numtype;
177 Matrix<real_t> L(nbDofs_,nbDofs_);
178 PolynomialBasis Pk(_Pk,2,k); //Pk[X,Y] space equipped with canonical basis
179 std::vector<RefDof*>::iterator itd=refDofs.begin();
180 for (number_t i = 1; i <= nbDofs_; ++i, ++itd) //loop on nodes
181 {
182 std::vector<real_t>::const_iterator itpt=(*itd)->coords();
183 PolynomialBasis::iterator itp = Pk.begin();
184 for (number_t j=1; j<=nbDofs_; ++j, ++itp) L(i,j)=itp->eval(*itpt, *(itpt+1)); //loop on monomials
185 }
186 //compute coefficients of shape functions by solving linear systems
187 Matrix<real_t> sf(nbDofs_,_idMatrix);
188 real_t eps=theZeroThreshold;
189 number_t r;
190 bool ok=gaussMultipleSolver(L, sf, nbDofs_, eps, r);
191 if (!ok)
192 {
193 where("LagrangeTrianglePk::computeShapeFunctions()");
194 error("mat_noinvert");
195 }
196 sf.roundToZero();
197 Matrix<real_t>::iterator itm=sf.begin();
198 Vector<real_t> v(nbDofs_);
199 Vector<real_t>::iterator itv;
200 Wk.dimVar=2; Wk.dimVec=1;
201 Wk.name="P"+tostring(k);
202 for (number_t i=0; i<nbDofs_; i++)
203 {
204 for(itv=v.begin(); itv!=v.end(); ++itv, ++itm) *itv=*itm;
205 Wk.add(combine(Pk,v));
206 }
207 // Wk.clean(0.000001); // to be used eventually for printing purpose
208 dWk.resize(2);
209 dWk[0]=dx(Wk); dWk[1]=dy(Wk);
210 }
211
212 //! return nodes numbers of first order triangle elements when splitting current element
splitP1() const213 std::vector<std::vector<number_t> > LagrangeTriangle::splitP1() const {
214 std::vector<std::vector<number_t> > splitnum;
215 std::vector<number_t> tuple(3,0);
216 switch(interpolation_p->numtype) {
217 case 0 :
218 case 1 :
219 tuple[0]=1; tuple[1]=2; tuple[2]=3;
220 splitnum.push_back(tuple);
221 return splitnum;
222 case 2 :
223 tuple[0]=6; tuple[1]=5; tuple[2]=3;
224 splitnum.push_back(tuple);
225 tuple[0]=1; tuple[1]=4; tuple[2]=6;
226 splitnum.push_back(tuple);
227 tuple[0]=5; tuple[1]=6; tuple[2]=4;
228 splitnum.push_back(tuple);
229 tuple[0]=4; tuple[1]=2; tuple[2]=5;
230 splitnum.push_back(tuple);
231 return splitnum;
232 case 3 :
233 tuple[0]=6; tuple[1]=8; tuple[2]=3;
234 splitnum.push_back(tuple);
235 tuple[0]=8; tuple[1]=6; tuple[2]=10;
236 splitnum.push_back(tuple);
237 tuple[0]=9; tuple[1]=10; tuple[2]=6;
238 splitnum.push_back(tuple);
239 tuple[0]=10; tuple[1]=9; tuple[2]=4;
240 splitnum.push_back(tuple);
241 tuple[0]=1; tuple[1]=4; tuple[2]=9;
242 splitnum.push_back(tuple);
243 tuple[0]=10; tuple[1]=5; tuple[2]=8;
244 splitnum.push_back(tuple);
245 tuple[0]=5; tuple[1]=10; tuple[2]=7;
246 splitnum.push_back(tuple);
247 tuple[0]=4; tuple[1]=7; tuple[2]=10;
248 splitnum.push_back(tuple);
249 tuple[0]=7; tuple[1]=2; tuple[2]=5;
250 splitnum.push_back(tuple);
251 return splitnum;
252 case 4 :
253 tuple[0]=6; tuple[1]=11; tuple[2]=3;
254 splitnum.push_back(tuple);
255 tuple[0]=11; tuple[1]=6; tuple[2]=15;
256 splitnum.push_back(tuple);
257 tuple[0]=6; tuple[1]=9; tuple[2]=15;
258 splitnum.push_back(tuple);
259 tuple[0]=15; tuple[1]=9; tuple[2]=13;
260 splitnum.push_back(tuple);
261 tuple[0]=12; tuple[1]=13; tuple[2]=9;
262 splitnum.push_back(tuple);
263 tuple[0]=13; tuple[1]=12; tuple[2]=4;
264 splitnum.push_back(tuple);
265 tuple[0]=1; tuple[1]=4; tuple[2]=12;
266 splitnum.push_back(tuple);
267 tuple[0]=15; tuple[1]=8; tuple[2]=11;
268 splitnum.push_back(tuple);
269 tuple[0]=8; tuple[1]=15; tuple[2]=14;
270 splitnum.push_back(tuple);
271 tuple[0]=13; tuple[1]=14; tuple[2]=15;
272 splitnum.push_back(tuple);
273 tuple[0]=14; tuple[1]=13; tuple[2]=7;
274 splitnum.push_back(tuple);
275 tuple[0]=4; tuple[1]=7; tuple[2]=13;
276 splitnum.push_back(tuple);
277 tuple[0]=14; tuple[1]=5; tuple[2]=8;
278 splitnum.push_back(tuple);
279 tuple[0]=5; tuple[1]=14; tuple[2]=10;
280 splitnum.push_back(tuple);
281 tuple[0]=7; tuple[1]=10; tuple[2]=14;
282 splitnum.push_back(tuple);
283 tuple[0]=10; tuple[1]=2; tuple[2]=5;
284 splitnum.push_back(tuple);
285 return splitnum;
286 case 5 :
287 tuple[0]=6; tuple[1]=14; tuple[2]=3;
288 splitnum.push_back(tuple);
289 tuple[0]=14; tuple[1]=6; tuple[2]=18;
290 splitnum.push_back(tuple);
291 tuple[0]=9; tuple[1]=18; tuple[2]=6;
292 splitnum.push_back(tuple);
293 tuple[0]=18; tuple[1]=9; tuple[2]=21;
294 splitnum.push_back(tuple);
295 tuple[0]=12; tuple[1]=21; tuple[2]=9;
296 splitnum.push_back(tuple);
297 tuple[0]=21; tuple[1]=12; tuple[2]=16;
298 splitnum.push_back(tuple);
299 tuple[0]=15; tuple[1]=16; tuple[2]=12;
300 splitnum.push_back(tuple);
301 tuple[0]=16; tuple[1]=15; tuple[2]=4;
302 splitnum.push_back(tuple);
303 tuple[0]=1; tuple[1]=4; tuple[2]=15;
304 splitnum.push_back(tuple);
305 tuple[0]=18; tuple[1]=11; tuple[2]=14;
306 splitnum.push_back(tuple);
307 tuple[0]=11; tuple[1]=18; tuple[2]=20;
308 splitnum.push_back(tuple);
309 tuple[0]=21; tuple[1]=20; tuple[2]=18;
310 splitnum.push_back(tuple);
311 tuple[0]=20; tuple[1]=21; tuple[2]=19;
312 splitnum.push_back(tuple);
313 tuple[0]=16; tuple[1]=19; tuple[2]=21;
314 splitnum.push_back(tuple);
315 tuple[0]=19; tuple[1]=16; tuple[2]=7;
316 splitnum.push_back(tuple);
317 tuple[0]=4; tuple[1]=7; tuple[2]=16;
318 splitnum.push_back(tuple);
319 tuple[0]=20; tuple[1]=8; tuple[2]=11;
320 splitnum.push_back(tuple);
321 tuple[0]=8; tuple[1]=20; tuple[2]=17;
322 splitnum.push_back(tuple);
323 tuple[0]=19; tuple[1]=17; tuple[2]=20;
324 splitnum.push_back(tuple);
325 tuple[0]=17; tuple[1]=19; tuple[2]=10;
326 splitnum.push_back(tuple);
327 tuple[0]=7; tuple[1]=10; tuple[2]=19;
328 splitnum.push_back(tuple);
329 tuple[0]=17; tuple[1]=5; tuple[2]=8;
330 splitnum.push_back(tuple);
331 tuple[0]=5; tuple[1]=17; tuple[2]=13;
332 splitnum.push_back(tuple);
333 tuple[0]=10; tuple[1]=13; tuple[2]=17;
334 splitnum.push_back(tuple);
335 tuple[0]=13; tuple[1]=2; tuple[2]=5;
336 splitnum.push_back(tuple);
337 return splitnum;
338 case 6 :
339 tuple[0]=6; tuple[1]=17; tuple[2]=3;
340 splitnum.push_back(tuple);
341 tuple[0]=17; tuple[1]=6; tuple[2]=21;
342 splitnum.push_back(tuple);
343 tuple[0]=9; tuple[1]=21; tuple[2]=6;
344 splitnum.push_back(tuple);
345 tuple[0]=21; tuple[1]=9; tuple[2]=24;
346 splitnum.push_back(tuple);
347 tuple[0]=12; tuple[1]=24; tuple[2]=9;
348 splitnum.push_back(tuple);
349 tuple[0]=24; tuple[1]=12; tuple[2]=27;
350 splitnum.push_back(tuple);
351 tuple[0]=15; tuple[1]=27; tuple[2]=12;
352 splitnum.push_back(tuple);
353 tuple[0]=27; tuple[1]=15; tuple[2]=19;
354 splitnum.push_back(tuple);
355 tuple[0]=18; tuple[1]=19; tuple[2]=15;
356 splitnum.push_back(tuple);
357 tuple[0]=19; tuple[1]=18; tuple[2]=4;
358 splitnum.push_back(tuple);
359 tuple[0]=1; tuple[1]=4; tuple[2]=18;
360 splitnum.push_back(tuple);
361 tuple[0]=21; tuple[1]=14; tuple[2]=17;
362 splitnum.push_back(tuple);
363 tuple[0]=14; tuple[1]=21; tuple[2]=26;
364 splitnum.push_back(tuple);
365 tuple[0]=24; tuple[1]=26; tuple[2]=21;
366 splitnum.push_back(tuple);
367 tuple[0]=26; tuple[1]=24; tuple[2]=28;
368 splitnum.push_back(tuple);
369 tuple[0]=27; tuple[1]=28; tuple[2]=24;
370 splitnum.push_back(tuple);
371 tuple[0]=28; tuple[1]=27; tuple[2]=22;
372 splitnum.push_back(tuple);
373 tuple[0]=19; tuple[1]=22; tuple[2]=27;
374 splitnum.push_back(tuple);
375 tuple[0]=22; tuple[1]=19; tuple[2]=7;
376 splitnum.push_back(tuple);
377 tuple[0]=4; tuple[1]=7; tuple[2]=19;
378 splitnum.push_back(tuple);
379 tuple[0]=26; tuple[1]=11; tuple[2]=14;
380 splitnum.push_back(tuple);
381 tuple[0]=11; tuple[1]=26; tuple[2]=23;
382 splitnum.push_back(tuple);
383 tuple[0]=28; tuple[1]=23; tuple[2]=26;
384 splitnum.push_back(tuple);
385 tuple[0]=23; tuple[1]=28; tuple[2]=25;
386 splitnum.push_back(tuple);
387 tuple[0]=22; tuple[1]=25; tuple[2]=28;
388 splitnum.push_back(tuple);
389 tuple[0]=25; tuple[1]=22; tuple[2]=10;
390 splitnum.push_back(tuple);
391 tuple[0]=7; tuple[1]=10; tuple[2]=22;
392 splitnum.push_back(tuple);
393 tuple[0]=23; tuple[1]=8; tuple[2]=11;
394 splitnum.push_back(tuple);
395 tuple[0]=8; tuple[1]=23; tuple[2]=20;
396 splitnum.push_back(tuple);
397 tuple[0]=25; tuple[1]=20; tuple[2]=23;
398 splitnum.push_back(tuple);
399 tuple[0]=20; tuple[1]=25; tuple[2]=13;
400 splitnum.push_back(tuple);
401 tuple[0]=10; tuple[1]=13; tuple[2]=25;
402 splitnum.push_back(tuple);
403 tuple[0]=20; tuple[1]=5; tuple[2]=8;
404 splitnum.push_back(tuple);
405 tuple[0]=5; tuple[1]=20; tuple[2]=16;
406 splitnum.push_back(tuple);
407 tuple[0]=13; tuple[1]=16; tuple[2]=20;
408 splitnum.push_back(tuple);
409 tuple[0]=16; tuple[1]=2; tuple[2]=5;
410 splitnum.push_back(tuple);
411 return splitnum;
412 default: return splitLagrange2DToP1();
413 }
414 return splitnum;
415 }
416
417 //! return nodes numbers of first order triangle elements when splitting current element
splitO1() const418 splitvec_t LagrangeTriangle::splitO1() const
419 {
420 splitvec_t splitnum;
421 std::vector<std::vector<number_t> > buf;
422 buf=this->splitP1();
423 for (number_t i=0; i< buf.size(); i++) { splitnum.push_back(std::make_pair(_triangle, buf[i])); }
424 return splitnum;
425 }
426
427 /*! for a side of current element, return the list of (first order element number, side number) of the first order elements (P1)
428 as a list of pair<P1 element number, P1 side number>
429 example : P2 Lagrange
430 \verbatim
431
432 2 2
433 | \ | \
434 S2 5 4 S1 ==> 5---4
435 | \ | \ | \
436 3---6---1 3---6---1
437 S3
438 \endverbatim
439 Side S1 :
440 - Elt 2 (1 4 6), side (1,4) --> (2,1)
441 - Elt 4 (4 2 5), side (4,2) --> (4,1)
442
443 Side S2 :
444 - Elt 1 (6 5 3), side (5 3) --> (1,2)
445 - Elt 4 (4 2 5), side (2,5) --> (4,2)
446
447 Side S3 :
448 - Elt 1 (6 5 3), side (3 6) --> (1,3)
449 - Elt 2 (1 4 6), side (6,1) --> (4,3)
450 */
451
splitP1Side(number_t side) const452 std::vector<std::pair<number_t,number_t> > LagrangeTriangle::splitP1Side(number_t side) const {
453 if (side<1 || side>3) error("split_bad_side_num",side);
454 std::vector<std::pair<number_t,number_t> > sidenum;
455 switch(interpolation_p->numtype) {
456 case 0 :
457 case 1 :
458 if (side==1) { sidenum.push_back(std::pair<number_t,number_t>(1,1)); }
459 else if (side==2) { sidenum.push_back(std::pair<number_t,number_t>(1,2)); }
460 else { sidenum.push_back(std::pair<number_t,number_t>(1,3)); }
461 return sidenum;
462 case 2 :
463 if (side==1) {
464 sidenum.push_back(std::pair<number_t,number_t>(2,1));
465 sidenum.push_back(std::pair<number_t,number_t>(4,1));
466 }
467 else if (side==2) {
468 sidenum.push_back(std::pair<number_t,number_t>(1,2));
469 sidenum.push_back(std::pair<number_t,number_t>(4,2));
470 }
471 else {
472 sidenum.push_back(std::pair<number_t,number_t>(1,3));
473 sidenum.push_back(std::pair<number_t,number_t>(4,3));
474 }
475 return sidenum;
476 default : error("split_low_deg_Lagrange_elt", words("shape",_triangle) , 2, words("shape",_triangle));
477 }
478 return sidenum;
479 }
480
481 //! pointCoordinates defines Lagrange Reference Element point coordinates
vertexCoordinates()482 std::vector<RefDof*>::iterator LagrangeTriangle::vertexCoordinates()
483 {
484 std::vector<RefDof*>::iterator it_rd(refDofs.begin());
485 (*it_rd++)->coords(1., 0.);
486 (*it_rd++)->coords(0., 1.);
487 (*it_rd++)->coords(0., 0.);
488 return it_rd;
489 }
490
pointCoordinates()491 void LagrangeTriangle::pointCoordinates()
492 {
493 // trace_p->push("LagrangeTriangle::pointCoordinates");
494
495 /*
496 Local numbering of Lagrange Pk triangles
497 2 2 2 2 2
498 | \ | \ | \ | \ | \
499 3---1 5 4 5 7 5 10 5 13
500 k=1 | \ | \ | \ | \
501 .... 3---6---1 8 10 4 8 14 7 8 17 10
502 ....... k=2 | \ | \ | \
503 .............. 3---6---9---1 11 15 13 4 11 20 19 7
504 ................... k=3 | \ | \
505 ............................. 3---6---9--12---1 14 18 21 16 4
506 ..................................... k=4 | \
507 ................................................ 3---6---9--12--15---1
508 ......................................................... k=5
509 */
510 int interpNum = interpolation_p->numtype;
511 std::vector<RefDof*>::iterator it_rd=refDofs.begin();
512 switch (interpNum)
513 {
514 case 0: (*it_rd)->coords(over3_, over3_); break;
515 case _P1BubbleP3:
516 it_rd = vertexCoordinates();
517 (*it_rd++)->coords(over3_, over3_);
518 break;
519 default: // Lagrange triangle of order > 0
520 it_rd = vertexCoordinates();
521 if (interpNum > 1)
522 {
523 // correspondence between segment and triangle local numbering
524 std::vector<number_t> s2t(2*nbPts_);
525 tensorNumberingTriangle(interpNum, s2t);
526
527 RefElement* sideRef = sideRefElems_[0];
528 // start after triangle last vertex
529 number_t pt = 2*geomRefElem_p->nbVertices();
530 while (it_rd != refDofs.end())
531 {
532 real_t x = *(sideRef->refDofs[s2t[pt++]]->coords());
533 real_t y = *(sideRef->refDofs[s2t[pt++]]->coords());
534 (*it_rd++)->coords(x, y);
535 }
536 }
537 break;
538 }
539 // trace_p->pop();
540 }
541
542 //! output of Lagrange D.o.f's numbers on underlying P1 mesh D.o.f's
outputAsP1(std::ofstream & os,const int refNum,const std::vector<number_t> & ranks) const543 void LagrangeTriangle::outputAsP1(std::ofstream& os, const int refNum, const std::vector<number_t>& ranks) const
544 {
545 number_t interpNum = interpolation_p->numtype;
546
547 /*
548 Numbering of d.o.f is made starting with the outtermost layer of elements
549 First three triangles are {1,4,12}, {13,12,4} and {4,7,13} for interpNum=3
550 (edges of which are part of the edges of the original triangle)
551 then recursively to the next inner layer
552 */
553 int pr(0), pr1(-1), corner_1(1), corner, p1, p2(1), p3, p4;
554 while (interpNum > 3)
555 {
556 // first bottom right external "corner" bears number corner
557 corner = corner_1 - pr;
558 // first bottom right internal "corner" bears number corner_1
559 // 13 for interpNum = 4, 16 for interpNum = 5, 19 for interpNum = 6
560 corner_1 += int(3*interpNum);
561
562 // layers of elements "parallel" to { x_1 + x_2 = 1 }
563 p1 = corner++; p2 = p1 + 3; p4 = corner_1 - pr; p3 = p4 - 1;
564 for (number_t j = 1; j < interpNum - 1; j++)
565 {
566 simplexVertexOutput(os, refNum, ranks[p1 + pr1], ranks[p2 + pr1], ranks[p3 + pr1]);
567 simplexVertexOutput(os, refNum, ranks[p4 + pr1], ranks[p3 + pr1], ranks[p2 + pr1]);
568 p1 = p2; p2 += 3; p3 = p4; p4 += 3;
569 if (j == interpNum - 3) { corner_1 += 1; p4 = corner_1 - pr; }
570 }
571
572 simplexVertexOutput(os, refNum, ranks[p1 + pr1], ranks[p2 + pr1], ranks[corner_1 - 1]);
573
574 // layers of elements "parallel" to { x_1 = 0 }
575 p1 = p2; p2 = corner++; p3 = p2 + 3; p4 = corner_1 - pr;
576 for (number_t j = 1; j < interpNum - 1; j++)
577 {
578 simplexVertexOutput(os, refNum, ranks[p1 + pr1], ranks[p2 + pr1], ranks[p3 + pr1]);
579 simplexVertexOutput(os, refNum, ranks[p3 + pr1], ranks[p4 + pr1], ranks[p1 + pr1]);
580 p1 = p4; p2 = p3; p3 += 3; p4 += 3;
581 if (j == interpNum - 3) { corner_1 += 1; p4 = corner_1 - pr; }
582 }
583 simplexVertexOutput(os, refNum, ranks[p1 + pr1], ranks[p2 + pr1], ranks[p3 + pr1]);
584
585 // layers of elements "parallel" to { x_2 = 0 }
586 p4 = p1; p2 = p3; p3 = corner++; p1 = p3 + 3;
587 for (number_t j = 1; j < interpNum - 1; j++)
588 {
589 simplexVertexOutput(os, refNum, ranks[p1 + pr1], ranks[p2 + pr1], ranks[p3 + pr1]);
590 simplexVertexOutput(os, refNum, ranks[p2 + pr1], ranks[p1 + pr1], ranks[p4 + pr1]);
591 p3 = p1; p1 += 3; p2 = p4; p4 += 3;
592 if (j == interpNum - 3) { corner_1 -= 2; p4 = corner_1 - pr; }
593 }
594
595 simplexVertexOutput(os, refNum, ranks[p1 + pr1], ranks[p2 + pr1], ranks[p3 + pr1]);
596
597 // initialize for next inner layer
598 pr = corner_1 - 1; pr1 = pr - 1; interpNum -= 3;
599 }
600
601 pr = pr1;
602
603 // finally consider the three cases for innermost triangles(s)
604 if (interpNum == 1) { simplexVertexOutput(os, refNum, ranks[1 + pr], ranks[2 + pr], ranks[3 + pr]); }
605 else if (interpNum == 2)
606 {
607 /*
608 2 3
609 | \ | \
610 | \ | \
611 | \ | #2 \
612 5 4 ------> 5-------4
613 | \ | \ #4 | \
614 | \ | \ | \
615 | \ | #3 \ | #1 \
616 3-------6-------1 3-------6-------1
617 */
618 simplexVertexOutput(os, refNum, ranks[1 + pr], ranks[4 + pr], ranks[6 + pr]);
619 simplexVertexOutput(os, refNum, ranks[4 + pr], ranks[2 + pr], ranks[5 + pr]);
620 simplexVertexOutput(os, refNum, ranks[6 + pr], ranks[5 + pr], ranks[3 + pr]);
621 simplexVertexOutput(os, refNum, ranks[5 + pr], ranks[6 + pr], ranks[4 + pr]);
622 }
623 else if (interpNum == 3)
624 {
625 /*
626 2 2
627 | \ | \
628 | \ | \
629 | \ | #4 \
630 5 7 5------7
631 | \ | \ #5 | \
632 | \ -------> | \ | \
633 | \ | #6 \| #3 \
634 8 10 4 8------10-----4
635 | \ | \ #8 | \ #2 | \
636 | \ | \ | \ | \
637 | \ | #7 \| #9 \| #1 \
638 3------6------9------1 3------6------9------1
639 */
640
641 simplexVertexOutput(os, refNum, ranks[1 + pr], ranks[4 + pr], ranks[9 + pr]);
642 simplexVertexOutput(os, refNum, ranks[10 + pr], ranks[9 + pr], ranks[4 + pr]);
643 simplexVertexOutput(os, refNum, ranks[4 + pr], ranks[7 + pr], ranks[10 + pr]);
644 simplexVertexOutput(os, refNum, ranks[7 + pr], ranks[2 + pr], ranks[5 + pr]);
645 simplexVertexOutput(os, refNum, ranks[5 + pr], ranks[10 + pr], ranks[7 + pr]);
646 simplexVertexOutput(os, refNum, ranks[10 + pr], ranks[5 + pr], ranks[8 + pr]);
647 simplexVertexOutput(os, refNum, ranks[6 + pr], ranks[8 + pr], ranks[3 + pr]);
648 simplexVertexOutput(os, refNum, ranks[6 + pr], ranks[8 + pr], ranks[10 + pr]);
649 simplexVertexOutput(os, refNum, ranks[9 + pr], ranks[10 + pr], ranks[6 + pr]);
650 }
651 else { noSuchFunction("outputAsP1"); }
652 } // end of outputAsP1
653
654 //! Lagrange standard P0 triangle Reference Element
655 template<>
LagrangeStdTriangle(const Interpolation * interp_p)656 LagrangeStdTriangle<_P0>::LagrangeStdTriangle(const Interpolation* interp_p)
657 : LagrangeTriangle(interp_p)
658 {
659 name_ += "_0";
660 // local coordinates of Points supporting D.o.F
661 pointCoordinates();
662 // build O1 splitting scheme
663 splitO1Scheme = splitO1();
664 }
665
666 template<>
~LagrangeStdTriangle()667 LagrangeStdTriangle<_P0>::~LagrangeStdTriangle() {}
668
669 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const670 void LagrangeStdTriangle<_P0>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
671 {
672 std::vector<real_t>::iterator it_w=shv.w.begin();
673 *it_w = 1.;
674
675 if (withDeriv)
676 {
677 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
678 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
679 *it_dwdx1i = 0.; *it_dwdx2i = 0.;
680 }
681 }
682
683 //! Lagrange standard P1 triangle Reference Element
684 template<>
LagrangeStdTriangle(const Interpolation * interp_p)685 LagrangeStdTriangle<_P1>::LagrangeStdTriangle(const Interpolation* interp_p)
686 : LagrangeTriangle(interp_p)
687 {
688 name_ += "_1";
689 // local coordinates of Points supporting D.o.F
690 pointCoordinates();
691 // build O1 splitting scheme
692 splitO1Scheme = splitO1();
693
694 //#ifdef FIG4TEX_OUTPUT
695 // Draw a LagrangeStdTriangle using TeX and fig4TeX
696 // fig4TeX();
697 //#endif
698 }
699
700 template<>
~LagrangeStdTriangle()701 LagrangeStdTriangle<_P1>::~LagrangeStdTriangle() {}
702
703 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const704 void LagrangeStdTriangle<_P1>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
705 {
706 // In the (x1, x2) reference plane shape functions of the Lagrange P1 triangle
707 // are given by (see local numbering)
708 // w_1(x1,x2) = x1
709 // w_2(x1,x2) = x2
710 // w_3(x1,x2) = x3 = 1 - x1 - x2
711 //
712 real_t x1=*it_pt, x2=*(it_pt+1);
713
714 std::vector<real_t>::iterator it_w=shv.w.begin();
715 *it_w++ = x1;
716 *it_w++ = x2;
717 *it_w = 1.-x1-x2;
718
719 if (withDeriv)
720 {
721 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
722 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
723 *it_dwdx1i++ = 1.; *it_dwdx2i++ = 0.;
724 *it_dwdx1i++ = 0.; *it_dwdx2i++ = 1.;
725 *it_dwdx1i = -1.; *it_dwdx2i = -1.;
726 }
727 }
728
729 //! Lagrange standard Pk triangle Reference Element
730 template<number_t Pk>
LagrangeStdTriangle(const Interpolation * interp_p)731 LagrangeStdTriangle<Pk>::LagrangeStdTriangle(const Interpolation* interp_p)
732 : LagrangeTriangle(interp_p)
733 {
734 name_ += "_" + tostring(Pk);
735 pointCoordinates(); // local coordinates of Points supporting D.o.F
736 // build O1 splitting scheme
737 splitO1Scheme = splitO1();
738
739 //#ifdef FIG4TEX_OUTPUT
740 // Draw a LagrangeStdTriangle using TeX and fig4TeX
741 // fig4TeX();
742 //#endif
743 }
744
745 template<number_t Pk>
~LagrangeStdTriangle()746 LagrangeStdTriangle<Pk>::~LagrangeStdTriangle() {}
747
748 //! computeShapeValues defines Lagrange Pk Reference Element shape functions
749 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const750 void LagrangeStdTriangle<_P2>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
751 {
752 // In the (x1, x2) reference plane shape functions of the Lagrange P2 triangle
753 // are given, with x3 = 1 - x1 - x2, by (see local numbering)
754 // w_i(x1,x2) = xi(2xi-1) (i = 1,2,3)
755 // w_4(x1,x2) = 4 x1 x2
756 // w_5(x1,x2) = 4 x2 x3
757 // w_6(x1,x2) = 4 x3 x1
758 //
759 real_t x1=*it_pt, x2=*(it_pt+1), x3=1.-x1-x2, qx1=4.*x1, qx2=4.*x2, qx3=4.*x3;
760 std::vector<real_t>::iterator it_w=shv.w.begin();
761 *it_w++ = x1*(2.*x1-1.);
762 *it_w++ = x2*(2.*x2-1.);
763 *it_w++ = x3*(2.*x3-1.);
764 *it_w++ = qx1*x2;
765 *it_w++ = qx2*x3;
766 *it_w = qx3*x1;
767
768 if (withDeriv)
769 {
770 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
771 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
772 *it_dwdx1i++ = qx1 - 1.; *it_dwdx2i++ = 0.;
773 *it_dwdx1i++ = 0.; *it_dwdx2i++ = qx2 - 1.;
774 *it_dwdx1i++ = -qx3 + 1.; *it_dwdx2i++ = -qx3 + 1.;
775 *it_dwdx1i++ = qx2; *it_dwdx2i++ = qx1;
776 *it_dwdx1i++ = -qx2; *it_dwdx2i++ = qx3 - qx2;
777 *it_dwdx1i = qx3 - qx1; *it_dwdx2i = -qx1;
778 }
779 }
780
781 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const782 void LagrangeStdTriangle<_P3>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
783 {
784 // In the (x1, x2) reference plane shape functions of the Lagrange P_3 triangle
785 // are given (with x3 = 1 - x1 - x2) by
786 // w_i(x1,x2) = (1/2) xi(3xi-1)(3xi-2) (i = 1,2,3)
787 // w_4(x1,x2) = (9/2) x1(2x1-1)x2
788 // w_5(x1,x2) = (9/2) x2(2x2-1)x3
789 // w_6(x1,x2) = (9/2) x3(2x3-1)x1
790 // w_7(x1,x2) = (9/2) x2(2x2-1)x1
791 // w_8(x1,x2) = (9/2) x3(2x3-1)x2
792 // w_9(x1,x2) = (9/2) x1(2x1-1)x3
793 // w_10(x1,x2) = 27 x1 x2 x3
794
795 real_t x1=*it_pt, x2=*(it_pt+1), x3=1.-x1-x2, x1x3=x1*x3, x2x3=x2*x3,
796 kx1m1=3.*x1-1., kx2m1=3.*x2-1., kx3m1=3.*x3-1.,
797 x1kx1m=x1*kx1m1, x2kx2m=x2*kx2m1, x3kx3m=x3*kx3m1;
798 real_t f9_2=9./2., f27_2=27./2.;
799
800 std::vector<real_t>::iterator it_w=shv.w.begin();
801 *it_w++ = x1kx1m*(kx1m1-1.)/2.;
802 *it_w++ = x2kx2m*(kx2m1-1.)/2.;
803 *it_w++ = x3kx3m*(kx3m1-1.)/2.;
804 *it_w++ = f9_2*x1kx1m*x2;
805 *it_w++ = f9_2*x2kx2m*x3;
806 *it_w++ = f9_2*x3kx3m*x1;
807 *it_w++ = f9_2*x2kx2m*x1;
808 *it_w++ = f9_2*x3kx3m*x2;
809 *it_w++ = f9_2*x1kx1m*x3;
810 *it_w = 27.*x1*x2x3;
811
812 if (withDeriv)
813 {
814 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
815 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
816 *it_dwdx1i++ = 1. + (f27_2*x1-9.)*x1; *it_dwdx2i++ = 0.;
817 *it_dwdx1i++ = 0.; *it_dwdx2i++ = 1. + (f27_2*x2-9.)*x2;
818 *it_dwdx1i++ = -1. - (f27_2*x3-9.)*x3; *it_dwdx2i++ = -1. - (f27_2*x3-9.)*x3;
819 *it_dwdx1i++ = -(f9_2-27.*x1)*x2; *it_dwdx2i++ = f9_2*x1kx1m;
820 *it_dwdx1i++ = -f9_2*x2kx2m; *it_dwdx2i++ = f9_2*((x3-x2)*kx2m1+3.*x2x3);
821 *it_dwdx1i++ = f9_2*((x3-x1)*kx3m1-3.*x1x3); *it_dwdx2i++ = (f9_2-27.*x3)*x1;
822 *it_dwdx1i++ = f9_2*x2kx2m; *it_dwdx2i++ = -(f9_2-27.*x2)*x1;
823 *it_dwdx1i++ = (f9_2-27.*x3)*x2; *it_dwdx2i++ = f9_2*((x3-x2)*kx3m1-3.*x2x3);
824 *it_dwdx1i++ = f9_2*((x3-x1)*kx1m1+3.*x1x3); *it_dwdx2i++ = -f9_2*x1kx1m;
825 *it_dwdx1i = 27.*x2*(x3-x1); *it_dwdx2i = 27.*x1*(x3-x2);
826 }
827 }
828
829 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const830 void LagrangeStdTriangle<_P4>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
831 {
832 // In the (x1, x2) reference plane shape functions of the Lagrange P_4 triangle
833 // are given (with x3 = 1 - x1 - x2) by
834 //
835 // ... well, they are not that given .....
836
837 real_t x1=*it_pt, x2=*(it_pt+1), x3=1.-x1-x2, x1x2=x1*x2,
838 kx1=4.*x1, kx1m1=kx1-1., kx1m2=kx1-2., x1kx1m=x1*kx1m1,
839 kx2=4.*x2, kx2m1=kx2-1., kx2m2=kx2-2., x2kx2m=x2*kx2m1,
840 kx3=4.*x3, kx2x3=kx2*x3;
841 real_t t18=3.-kx1-kx2, t19=kx3*t18, t20=2.-kx1-kx2;
842 real_t f08_3=8./3.;
843
844 std::vector<real_t>::iterator it_w=shv.w.begin();
845 *it_w++ = x1kx1m*kx1m2*(kx1-3.)/6.;
846 *it_w++ = x2kx2m*kx2m2*(-3+kx2)/6.;
847 *it_w++ = t19*t20*(1-kx1-kx2)/24.;
848 *it_w++ = f08_3*x1kx1m*kx1m2*x2;
849 *it_w++ = f08_3*x2kx2m*kx2m2*x3;
850 *it_w++ = f08_3*x1*x3*t18*t20;
851 *it_w++ = 4.*x1kx1m*x2kx2m;
852 *it_w++ = x2kx2m*t19;
853 *it_w++ = x1kx1m*t19;
854 *it_w++ = f08_3*x1x2*kx2m1*kx2m2;
855 *it_w++ = 2./3.*kx2x3*t18*t20;
856 *it_w++ = f08_3*x1kx1m*kx1m2*x3;
857 *it_w++ = 8.*x1kx1m*kx2x3;
858 *it_w++ = 32.*x1x2*kx2m1*x3;
859 *it_w = 8.*x1x2*t19;
860
861 if (withDeriv)
862 {
863 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
864 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
865 real_t t4_1=128./3.*x1, t4_2=128./3.*x2;
866 real_t t4_12=128.*x2, t4_20=16./3.*x2, t4_21=64.*x2;
867 real_t t4_23=128.*x1x2;
868 real_t t4_32=f08_3*x2kx2m*kx2m2, t4_33=1.-x2, t4_34=3.-kx2, t4_38=256.*x2, t4_39=384.-t4_38;
869 real_t t4_41=384.*x2, t4_42=512./3.*x1, t4_49=x1x2*kx2m1;
870 real_t t4_50=32.*t4_49;
871 real_t t4_67=(192.-t4_12)*x2;
872 real_t t4_80=4.*x2*t4_33, t4_84=384.*x1x2;
873 real_t t4_94=512.*x2;
874 real_t t4_109=-25./3.+(140./3.+(-80.+t4_2)*x2)*x2+(140./3+(-160.+t4_12)*x2+(-80.+t4_12+t4_1)*x1)*x1;
875 real_t t4_110=f08_3*x1*(4.*x1-1.)*(4.*x1-2.), t4_111=512./3.*x2;
876 real_t t4_116=64.-t4_12, t4_121=128.*x1, t4_126=32.*x2, t4_128=(t4_12-16.)*x1;
877
878 *it_dwdx1i++ = -1. + (44./3.+(-48.+t4_1)*x1)*x1; *it_dwdx2i++ = 0.;
879 *it_dwdx1i++ = 0.; *it_dwdx2i++ = -1. + (44./3.+(-48.+t4_2)*x2)*x2;
880 *it_dwdx1i++ = t4_109; *it_dwdx2i++ = t4_109;
881 *it_dwdx1i++ = t4_20 + (-t4_21+t4_23)*x1; *it_dwdx2i++ = t4_110;
882 *it_dwdx1i++ = -t4_32; *it_dwdx2i++ = 16./3. + (-224./3.+(224.-t4_111)*x2)*x2 + (-16./3.+t4_116*x2)*x1;
883 *it_dwdx1i++ = -f08_3*t4_33*t4_34*kx2m2 + (-416./3.+t4_39*x2+(288.-384.*x2-t4_42)*x1)*x1; *it_dwdx2i++ = (-208./3.+t4_67+(192.-t4_38-t4_121)*x1)*x1;
884 *it_dwdx1i++ = -4.*x2kx2m + t4_50; *it_dwdx2i++ = (4.-t4_126+t4_128)*x1;
885 *it_dwdx1i++ = (28.+(-144.+t4_12)*x2)*x2 + t4_50; *it_dwdx2i++ = -12. + (152.-t4_39*x2)*x2 + (28.+(-288.+384.*x2)*x2+t4_128)*x1;
886 *it_dwdx1i++ = -4.*t4_33*t4_34 + (152.+(-288.+t4_12)*x2+(-384.+384.*x2+256.*x1)*x1)*x1; *it_dwdx2i++ = (28.-t4_126+(-144.+t4_12+t4_121)*x1)*x1;
887 *it_dwdx1i++ = t4_32; *it_dwdx2i++ = (16./3.-t4_116*x2)*x1;
888 *it_dwdx1i++ = (-208./3.+t4_67)*x2 + ((192.-256.*x2)*x2-t4_23)*x1; *it_dwdx2i++ = 16. + (-416./3.+(288.-t4_111)*x2)*x2 + (-208./3.+384.*x2*t4_33+(96.-t4_38-t4_1)*x1)*x1;
889 *it_dwdx1i++ = 16./3. - t4_20 + (-224./3.+t4_21+(224.-t4_12-t4_42)*x1)*x1; *it_dwdx2i++ = -t4_110;
890 *it_dwdx1i++ = -8.*t4_80 + ((320.-256.*x2)*x2-t4_84)*x1; *it_dwdx2i++ = (-32.+t4_21+(160.-t4_38-t4_121)*x1)*x1;
891 *it_dwdx1i++ = 32.*x2kx2m*t4_33 - 64.*t4_49; *it_dwdx2i++ = (-32.+(320.-t4_41)*x2+(-t4_38+32.)*x1)*x1;
892 *it_dwdx1i = 8.*t4_80*t4_34 + ((-448.+t4_94)*x2+t4_84)*x1; *it_dwdx2i = (96.+(-448.+t4_41)*x2+(t4_94-224.+t4_121)*x1)*x1;
893 }
894 }
895
896 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const897 void LagrangeStdTriangle<_P5>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
898 {
899 real_t x1=*it_pt, x2=*(it_pt+1), x3=1.-x1-x2, x1x2=x1*x2,
900 kx1=5.*x1, kx1m1=kx1-1., kx1m2=kx1-2., kx1m3=kx1-3., kx1m4=kx1-4., x1kx1m=x1*kx1m1,
901 kx2=5.*x2, kx2m1=kx2-1., kx2m2=kx2-2., kx2m3=kx2-3., kx2m4=kx2-4., x2kx2m=x2*kx2m1, kx2x3=kx2*x3;
902 real_t f125_3=125./3., f25_4=25./4., f125_4=125./4., f25_6=25./6., f125_6=125./6.,
903 f05_12=5./12., f25_12=25./12., f05_24=5./24., f25_24=25./24.;
904
905 std::vector<real_t>::iterator it_w=shv.w.begin();
906 real_t t5_22 = 4.-kx1-kx2;
907 real_t t5_23 = 5.*x3*t5_22;
908 real_t t5_24 = 3.-kx1-kx2;
909 real_t t5_25 = 2.-kx1-kx2;
910 real_t t5_26 = t5_22*t5_24*t5_25;
911 real_t t5_27 = t5_23*t5_24;
912
913 *it_w++ = x1kx1m*kx1m2*kx1m3*kx1m4/24.;
914 *it_w++ = x2kx2m*kx2m2*kx2m3*kx2m4/24.;
915 *it_w++ = t5_23*t5_24*t5_25*(1.-kx1-kx2)/120.;
916 *it_w++ = f25_24*x1kx1m*kx1m2*kx1m3*x2;
917 *it_w++ = f25_24*x2kx2m*kx2m2*kx2m3*x3;
918 *it_w++ = f25_24*x1*x3*t5_26;
919 *it_w++ = f25_12*x1kx1m*kx1m2*x2*kx2m1;
920 *it_w++ = f25_12*x2kx2m*kx2m2*x3*t5_22;
921 *it_w++ = f05_12*x1kx1m*t5_27;
922 *it_w++ = f25_12*x1kx1m*x2kx2m*kx2m2;
923 *it_w++ = f05_12*x2kx2m*t5_27;
924 *it_w++ = f25_12*x1kx1m*kx1m2*x3*t5_22;
925 *it_w++ = f25_24*x1x2*kx2m1*kx2m2*kx2m3;
926 *it_w++ = f05_24*kx2x3*t5_26;
927 *it_w++ = f25_24*x1kx1m*kx1m2*kx1m3*x3;
928 *it_w++ = f125_6*x1kx1m*kx1m2*x2*x3;
929 *it_w++ = f125_6*x1x2*kx2m1*kx2m2*x3;
930 *it_w++ = f25_6*x1x2*t5_27;
931 *it_w++ = f125_4*x1kx1m*x2kx2m*x3;
932 *it_w++ = f125_4*x1x2*kx2m1*x3*t5_22;
933 *it_w++ = f25_4*x1kx1m*kx2x3*t5_22;
934
935 if (withDeriv)
936 {
937 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
938 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
939 real_t T1 = 3125./24.*x1;
940 real_t T2 = 3125./24.*x2;
941 real_t T3 = 3125./6.*x2;
942 real_t T4 = 3125./4.*x2;
943 real_t T30 = -137./12. + (375./4.+(-2125./8.+(312.5-T2)*x2)*x2)*x2 + (375./4.+(-2125./4.+(1875./2.-T3)*x2)*x2+(-2125./8.+(1875./2.-T4)*x2+(312.5-T3-T1)*x1)*x1)*x1;
944 real_t T6 = f25_4*x2;
945 real_t T7 = 1375./12.*x2;
946 real_t T8 = 1875./4.*x2;
947 real_t T9 = 3125./6.*x1x2;
948 real_t T48 = f25_24*x2kx2m*kx2m2*kx2m3;
949 real_t T49 = 1.-x2;
950 real_t T51 = -5.*T49*kx2m4;
951 real_t T54 = 3125./3.*x2;
952 real_t T59 = 9375./4.*x2;
953 real_t T62 = 6250./3.*x2;
954 real_t T63 = 15625./24.*x1;
955 real_t T74 = 625./4.*x1x2*kx2m1;
956 real_t T85 = x1x2*kx2m1*kx2m2;
957 real_t T86 = f125_6*T85;
958 real_t T97 = 3125.*x2;
959 real_t T98 = 15625./12.*x1;
960 real_t T123 = 312.5*x2;
961 real_t T138 = (8875./12.+(-4375./4.+T3)*x2)*x2;
962 real_t T141 = 1562.5*x2;
963 real_t T142 = -4375./2.+T141;
964 real_t T143 = T142*x2;
965 real_t T160 = 5.*T49*x2;
966 real_t T162 = 625.*x2;
967 real_t T167 = 6250./3.*x1x2;
968 real_t T185 = 9375./2.*x2;
969 real_t T198 = 1875./4.*x1x2*kx2m1;
970 real_t T238 = f25_24*x1*kx1m1*kx1m2*kx1m3;
971 real_t T239 = 15625./24.*x2;
972 real_t T246 = 1875./4.-T3;
973 real_t T253 = 3125./6.*x1;
974 real_t T260 = f125_3*x2;
975 real_t T261 = T3-625./12.;
976 real_t T267 = 15625./12.*x2;
977 real_t T281 = (f125_6+(-312.5+T4)*x2)*x1;
978 real_t T285 = 625./4.*x2;
979 real_t T290 = 3125./4.*x1;
980 real_t T340 = -4375./2.+T59;
981 real_t T381 = -T141+625./4.;
982 real_t T382 = f05_24*T51*kx2m2*kx2m3 + (-1925./6.+(8875./6.+(-4375./2.+T54)*x2)*x2+(8875./8.+(-13125./4.+T59)*x2+(-4375./3.+T62+T63)*x1)*x1)*x1;
983 real_t T383 = -f25_12*T49*kx2m4*kx2m3 + (2675./6.+(-8875./6.+(1562.5-T3)*x2)*x2+(-7375./4.+(16875./4.-T59)*x2+(8125./3.-T97-T98)*x1)*x1)*x1;
984 real_t T384 = (1175./12.+(-8875./12.+(5625./4.-T4)*x2)*x2)*x2 + ((-250.+1562.5*T49*x2)*x2-T74)*x1;
985 real_t T385 = -f25_6*T49*kx2m4+(-325.+(3875./6.-T123)*x2+(6125./4.+(-9375./4.+T4)*x2+(-2500.+T62+T98)*x1)*x1)*x1;
986
987 *it_dwdx1i++ = 1. + (-f125_6+(875./8.+(-625./3.+T1)*x1)*x1)*x1; *it_dwdx2i++ = 0.;
988 *it_dwdx1i++ = 0.; *it_dwdx2i++ = 1. + (-f125_6+(875./8.+(-625./3.+T2)*x2)*x2)*x2;
989 *it_dwdx1i++ = T30; *it_dwdx2i++ = T30;
990 *it_dwdx1i++ = -T6 + (T7 +(-T8+T9)*x1)*x1; *it_dwdx2i++ = T238;
991 *it_dwdx1i++ = -T48; *it_dwdx2i++ = -f25_4 + (1525./12.+(-5125./8.+(6875./6.-T239)*x2)*x2)*x2 + (f25_4+(-1375./12.+T246*x2)*x2)*x1;
992 *it_dwdx1i++ = T382; *it_dwdx2i++ = (-1925./12.+T138+(8875./12.+T143+(-4375./4.+T141+T253)*x1)*x1)*x1;
993 *it_dwdx1i++ = f25_6*x2kx2m + (-62.5*x2kx2m+T74)*x1; *it_dwdx2i++ = (T260-f25_6+(-T123+f125_4+T261*x1)*x1)*x1;
994 *it_dwdx1i++ = (-37.5+(3875./12.+(-3125./4.+T3)*x2)*x2)*x2 + T86; *it_dwdx2i++ = 50./3. + (-325.+(6125./4.+(-2500.+T267)*x2)*x2)*x2 + (-37.5+(3875./6.+(-9375./4.+T62)*x2)*x2+T281)*x1;
995 *it_dwdx1i++ = T383; *it_dwdx2i++ = (1175./12.+(-250.+T285)*x2+(-8875./12.+(1562.5-T4)*x2+(5625./4.-T141-T290)*x1)*x1)*x1;
996 *it_dwdx1i++ = -f25_12*x2kx2m*kx2m2 + T86; *it_dwdx2i++ = (-f25_6+(62.5-T285)*x2+T281)*x1;
997 *it_dwdx1i++ = T384; *it_dwdx2i++ = -25. + (2675./6.+(-7375./4.+(8125./3.-T267)*x2)*x2)*x2 + (1175./12.+(-8875./6.+(16875./4.-T97)*x2)*x2 + (-125.+(1562.5-T59)*x2-T261*x1)*x1)*x1;
998 *it_dwdx1i++ = T385; *it_dwdx2i++ = (-37.5+T260+(3875./12.-T123+(-3125./4.+T3+T253)*x1)*x1)*x1;
999 *it_dwdx1i++ = T48; *it_dwdx2i++ = (-f25_4+(1375./12.-T246*x2)*x2)*x1;
1000 *it_dwdx1i++ = (-1925./12.+T138)*x2 + ((8875./12.+T143)*x2+((-4375./4.+T141)*x2+T9)*x1)*x1; *it_dwdx2i++ = 25. + (-1925./6.+(8875./8.+(-4375./3.+T239)*x2)*x2)*x2 + (-1925./12.+(8875./6.+(-13125./4.+T62)*x2)*x2+(8875./24.+T340*x2+(T54-4375./12.+T1)*x1)*x1)*x1;
1001 *it_dwdx1i++ = -f25_4 + T6 + (1525./12.-T7+(-5125./8.+T8+(6875./6.-T3-T63)*x1)*x1)*x1; *it_dwdx2i++ = -T238;
1002 *it_dwdx1i++ = 25./3.*T160 + ((-2125./3.+T162)*x2+((2500.-T141)*x2-T167)*x1)*x1; *it_dwdx2i++ = (f125_3-250./3.*x2+(-2125./6.+T162+(2500./3.-T54-T253)*x1)*x1)*x1;
1003 *it_dwdx1i++ = f125_6*x2kx2m*kx2m2*T49 - f125_3*T85; *it_dwdx2i++ = (f125_3+(-2125./3.+(2500.-T62)*x2)*x2+(-f125_3+(625.-T141)*x2)*x1)*x1;
1004 *it_dwdx1i++ = f25_6*T160*kx2m4*kx2m3 + ((-5875./3.+(5000.-T97)*x2)*x2+((3750.-T185)*x2-T167)*x1)*x1; *it_dwdx2i++ = (250.+(-5875./3.+(3750.-T62)*x2)*x2+(-5875./6.+(5000.-T185)*x2+(1250.-T97-T253)*x1)*x1)*x1;
1005 *it_dwdx1i++ = -f125_4*x2kx2m*T49 + ((-375.-T142*x2)*x2-T198)*x1; *it_dwdx2i++ = (f125_4+(-375.+T8)*x2+(-375./2.-T340*x2+T381*x1)*x1)*x1;
1006 *it_dwdx1i++ = f25_4*x2kx2m*T51 + ((1125./2.+(-6875./2.+T97)*x2)*x2+T198)*x1; *it_dwdx2i++ = (-125.+(3625./2.+(-9375./2.+T97)*x2)*x2+(1125./4.+(-6875./2.+T185)*x2-T381*x1)*x1)*x1;
1007 *it_dwdx1i = f25_4*T160*kx2m4 + ((3625./2.+(-6875./2.+T141)*x2)*x2+(-9375./2.*T49*x2+3125.*x1x2)*x1)*x1; *it_dwdx2i = (-125.+(1125./2.-T8)*x2+(3625./4.+(-6875./2.+T59)*x2+(-1562.5+T97+T290)*x1)*x1)*x1;
1008 }
1009 }
1010
1011 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const1012 void LagrangeStdTriangle<_P6>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
1013 {
1014 real_t x1=*it_pt, x2=*(it_pt+1);
1015
1016 std::vector<real_t>::iterator it_w=shv.w.begin();
1017 real_t f3_10 = 0.3, f3_4 = 0.75;
1018
1019 real_t T1 = 6.*x1;
1020 real_t T3 = x1*(T1-1.);
1021 real_t T5 = T3*(T1-2.);
1022 real_t T6 = T1-3.;
1023 real_t T8 = T6*(T1-4.);
1024 real_t T13 = 6.*x2;
1025 real_t T14 = T13-1.;
1026 real_t T15 = x2*T14;
1027 real_t T16 = T13-2.;
1028 real_t T17 = T15*T16;
1029 real_t T18 = T13-3.;
1030 real_t T19 = T13-4.;
1031 real_t T20 = T18*T19;
1032 real_t T25 = 1.-x1-x2;
1033 real_t T26 = 5.-T1-T13;
1034 real_t T28 = 4.-T1-T13;
1035 real_t T29 = 6.*T25*T26*T28;
1036 real_t T30 = 3.-T1-T13;
1037 real_t T31 = 2.-T1-T13;
1038 real_t T46 = T28*T30*T31;
1039 real_t T49 = T6*x2;
1040 real_t T59 = T26*T28*T30;
1041 real_t T67 = T3*x2;
1042 real_t T68 = T14*T16;
1043 real_t T72 = 6.*T15*T25;
1044 real_t T79 = x1*x2;
1045 real_t T80 = T79*T14;
1046 real_t T81 = T16*T18;
1047 real_t T86 = 6.*x2*T25*T26;
1048
1049 *it_w++ = T5*T8*(T1-5.)/120.;
1050 *it_w++ = T17*T20*(T13-5.)/120.;
1051 *it_w++ = T29*T30*T31*(1.-T1-T13)/720.;
1052 *it_w++ = f3_10*T5*T8*x2;
1053 *it_w++ = f3_10*T17*T20*T25;
1054 *it_w++ = f3_10*x1*T25*T26*T46;
1055 *it_w++ = f3_4*T5*T49*T14;
1056 *it_w++ = f3_4*T17*T18*T25*T26;
1057 *it_w++ = f3_4*T3*T25*T59;
1058 *it_w++ = T5*T17;
1059 *it_w++ = T17*T29/6.;
1060 *it_w++ = T5*T29/6.;
1061 *it_w++ = f3_4*T67*T68*T18;
1062 *it_w++ = T72*T59/8.;
1063 *it_w++ = f3_4*T5*T6*T25*T26;
1064 *it_w++ = f3_10*T80*T81*T19;
1065 *it_w++ = T86*T46/20.;
1066 *it_w++ = f3_10*T5*T8*T25;
1067 *it_w++ = 9.*T5*T49*T25;
1068 *it_w++ = 9.*T80*T81*T25;
1069 *it_w++ = 9.*T79*T25*T59;
1070 *it_w++ = 3.*T5*T72;
1071 *it_w++ = 18.*T80*T16*T25*T26;
1072 *it_w++ = 3.*T67*T29;
1073 *it_w++ = 18.*T67*T68*T25;
1074 *it_w++ = 3.*T80*T29;
1075 *it_w++ = 3.*T5*T86;
1076 *it_w = 27.*T67*T14*T25*T26;
1077
1078 if (withDeriv)
1079 {
1080 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1+1;
1081 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
1082 real_t f1944_5=1944./5., f812_5=812./5., fm1323_2=-1323./2., f36_5=36./5.;
1083 real_t f11664_5 = 11664./5.;
1084
1085 real_t T1 = f1944_5*x1;
1086 real_t T11 = f1944_5*x2;
1087 real_t T20 = 1944.*x2;
1088 real_t T27 = 3888.*x2;
1089 real_t T42 = -14.7 + (f812_5+(fm1323_2+(1260.+(-1134.+T11)*x2)*x2)*x2)*x2
1090 + (f812_5+(-1323.+(3780.+(-4536.+T20)*x2)*x2)*x2+(fm1323_2+(3780.+(-6804.+T27)*x2)*x2+(1260.+(-4536.+T27)*x2+(-1134.+T20+T1)*x1)*x1)*x1)*x1;
1091 real_t T43 = f36_5*x2;
1092 real_t T44 = 180.*x2;
1093 real_t T45 = 1134.*x2;
1094 real_t T46 = 2592.*x2;
1095 real_t T47 = x1*x2;
1096 real_t T48 = 1944.*T47;
1097 real_t T56 = 6.*x2;
1098 real_t T57 = T56 - 1.;
1099 real_t T58 = x2*T57;
1100 real_t T59 = T56 - 2.;
1101 real_t T60 = T56 - 3.;
1102 real_t T61 = T59*T60;
1103 real_t T62 = T56 - 4.;
1104 real_t T65 = f3_10*T58*T61*T62;
1105 real_t T66 = 1. - x2;
1106 real_t T67 = 5. - T56;
1107 real_t T68 = 6.*T66*T67;
1108 real_t T69 = T62*T60;
1109 real_t T74 = (10368.-T27)*x2;
1110 real_t T79 = 11664.*x2;
1111 real_t T84 = 15552.*x2;
1112 real_t T87 = 9720.*x2;
1113 real_t T88 = f11664_5*x1;
1114 real_t T101 = T47*T57;
1115 real_t T102 = 648.*T101;
1116 real_t T116 = T57*T59;
1117 real_t T118 = T47*T116*T60;
1118 real_t T119 = 9.*T118;
1119 real_t T134 = 23328.*x2;
1120 real_t T137 = 19440.*x2;
1121 real_t T138 = 5832.*x1;
1122 real_t T148 = T58*T59;
1123 real_t T151 = T47*T116;
1124 real_t T152 = 108.*T151;
1125 real_t T164 = 7776.*x2;
1126 real_t T209 = -19440. + T79;
1127 real_t T210 = T209*x2;
1128 real_t T226 = 594.*x2;
1129 real_t T248 = (2088.+(-5022.+(5184.-T20)*x2)*x2)*x2;
1130 real_t T252 = (15552.-T164)*x2;
1131 real_t T254 = (-10044.+T252)*x2;
1132 real_t T258 = (15552.-T79)*x2;
1133 real_t T279 = 6.*x2*T66;
1134 real_t T281 = 1188.*x2;
1135 real_t T284 = 5832.*x2;
1136 real_t T287 = 9720.*T47;
1137 real_t T300 = -T67*T62;
1138 real_t T304 = -34992. + T84;
1139 real_t T310 = 34992.*x2;
1140 real_t T315 = 31104.*x2;
1141 real_t T334 = 2592.*T101;
1142 real_t T340 = 6.*T59*T66;
1143 real_t T350 = 324.*T151;
1144 real_t T366 = 46656.*x2;
1145 real_t T369 = 19440.*T47;
1146 real_t T451 = 6.*x1;
1147 real_t T460 = f3_10*x1*(T451-1.)*(T451-2.)*(T451-3.)*(T451-4.);
1148 real_t T461 = f11664_5*x2;
1149 real_t T470 = 2592. - T20;
1150 real_t T479 = 1944.*x1;
1151 real_t T488 = 54.*x2;
1152 real_t T490 = (T20-162.)*x1;
1153 real_t T516 = (-27.+(594.+(-2916.+T27)*x2)*x2)*x1;
1154 real_t T520 = 648.*x2;
1155 real_t T529 = 3888.*x1;
1156 real_t T538 = 216.*x2;
1157 real_t T543 = -1296. + T27;
1158 real_t T714 = (-T164+648)*x1;
1159 real_t T731 = -3888. + T79;
1160
1161 real_t f137_5=137./5., fm405_2=-405./2., fm3132_5=-3132./5., f9_2=9./2.,
1162 f99_2=99./2., fm8289_2=-8289./2., fm1197_2=-1197./2., fm12447_2=-12447./2.,
1163 f513_2=513./2., fm1566_5=-1566./5., fm972_5=-972./5., fm3_2=-3./2.,
1164 fm45_2=-45./2., fm1071_2=-1071./2.;
1165
1166 *it_dwdx1i++ = -1. + (f137_5+(fm405_2+(612.+(-810.+T1)*x1)*x1)*x1)*x1;
1167 *it_dwdx1i++ = 0.;
1168 *it_dwdx1i++ = T42;
1169 *it_dwdx1i++ = T43 + (-T44+(T45+(-T46+T48)*x1)*x1)*x1;
1170 *it_dwdx1i++ = -T65;
1171 *it_dwdx1i++ = -T68*T69*T59/20. + (fm3132_5+(4176.+(-10044.+T74)*x2)*x2+(3132.+(-15066.+(23328.-T79)*x2)*x2+(-6696.+(20736.-T84)*x2+(6480.-T87-T88)*x1)*x1)*x1)*x1;
1172 *it_dwdx1i++ = -f9_2*T58 + (99.*T58+(-486.*T58+T102)*x1)*x1;
1173 *it_dwdx1i++ = (f99_2+(fm1197_2+(2376.+(-3726.+T20)*x2)*x2)*x2)*x2+T119;
1174 *it_dwdx1i++ = -f3_4*T66*T67*T69 + (1053.+(-5220.+(9342.+(-7128.+T20)*x2)*x2)*x2+(fm12447_2+(23652.+(-29160.+T79)*x2)*x2+(14796.+(-37584.+T134)*x2+(-15390.+T137+T138)*x1)*x1)*x1)*x1;
1175 *it_dwdx1i++ = 2.*T148 + (-36.*T148+T152)*x1;
1176 *it_dwdx1i++ = (-148.+(1692.+(-6120.+(8424.-T27)*x2)*x2)*x2)*x2 + ((360.+(-3672.+(10368.-T164)*x2)*x2)*x2-T152)*x1;
1177 *it_dwdx1i++ = -2.*T66*T67*T62 + (-1016.+(3384.+(-3672.+1296*x2)*x2)*x2+(6696.+(-18360.+(15552.-T27)*x2)*x2+(-17424.+(33696.-T84)*x2+(19440.-T137-7776.*x1)*x1)*x1)*x1)*x1;
1178 *it_dwdx1i++ = -f3_4*T58*T61 + T119;
1179 *it_dwdx1i++ = (f513_2+(-2610.+(7884.+(-9396.+T27)*x2)*x2)*x2)*x2 + ((-1071.+(9342.+T210)*x2)*x2+((1458.+(-10692.+T79)*x2)*x2+T102)*x1)*x1;
1180 *it_dwdx1i++ = -f9_2*T66*T67 + (594.+(-1197.+T226)*x2+(fm8289_2+(7128.-2916.*x2)*x2+(11556.+(-14904.+T27)*x2+(-13770.+T87+T138)*x1)*x1)*x1)*x1;
1181 *it_dwdx1i++ = T65;
1182 *it_dwdx1i++ = (fm1566_5+T248)*x2 + ((2088.+T254)*x2+((-5022.+T258)*x2+((5184.-T164)*x2-T48)*x1)*x1)*x1;
1183 *it_dwdx1i++ = f36_5 - T43 + (fm972_5+T44+(1404.-T45+(-4104.+T46+(5184.-T20-T88)*x1)*x1)*x1)*x1;
1184 *it_dwdx1i++ = -9.*T279 + ((1296.-T281)*x2+((-7614.+T284)*x2+(T252-T287)*x1)*x1)*x1;
1185 *it_dwdx1i++ = 9.*T58*T61*T66 - 18.*T118;
1186 *it_dwdx1i++ = fm3_2*T279*T300*T60 + ((-6156.+(25704.+T304*x2)*x2)*x2+((19278.+(-52488.+T310)*x2)*x2+((-23328.+T315)*x2+T287)*x1)*x1)*x1;
1187 *it_dwdx1i++ = 36.*T58*T66 + ((720.+(-4968.+T27)*x2)*x2+((-2916.-T209*x2)*x2-T334)*x1)*x1;
1188 *it_dwdx1i++ = 3.*T58*T340*T67 + ((-792.+(7992.+(-22032.+T84)*x2)*x2)*x2+T350)*x1;
1189 *it_dwdx1i++ = -3.*T279*T300 + ((6984.+(-22464.+(23328.-T164)*x2)*x2)*x2+((-28836.+(64152.-T310)*x2)*x2+((41472.-T366)*x2-T369)*x1)*x1)*x1;
1190 *it_dwdx1i++ = -3.*T58*T340 + ((504.+(-4968.+(12960.-T164)*x2)*x2)*x2-T350)*x1;
1191 *it_dwdx1i++ = -3.*T58*T68*T62 + ((2664.+(-22464.+(42768.-T134)*x2)*x2)*x2 + ((-4860.+34992.*x2*T66)*x2-T334)*x1)*x1;
1192 *it_dwdx1i++ = 6.*T279*T67 + ((-4032.+(7992.-T27)*x2)*x2+((21060.+(-33048.+T79)*x2)*x2+((-36288.+T315)*x2+T369)*x1)*x1)*x1;
1193 *it_dwdx1i = -f9_2*T58*T68 + ((-2214.+(17496.+(-27216.+T79)*x2)*x2)*x2+((5832.+(-40824.+T310)*x2)*x2+3888.*T101)*x1)*x1;
1194
1195 *it_dwdx2i++ = 0.;
1196 *it_dwdx2i++ = -1. + (f137_5+(fm405_2+(612.+(-810.+T11)*x2)*x2)*x2)*x2;
1197 *it_dwdx2i++ = T42;
1198 *it_dwdx2i++ = T460;
1199 *it_dwdx2i++ = f36_5 + (fm972_5+(1404.+(-4104.+(5184.-T461)*x2)*x2)*x2)*x2 + (-f36_5+(180.+(-1134.+T470*x2)*x2)*x2)*x1;
1200 *it_dwdx2i++ = (fm1566_5+T248+(2088.+T254+(-5022.+T258+(5184.-T164-T479)*x1)*x1)*x1)*x1;
1201 *it_dwdx2i++ = (-T488+f9_2+(T226-f99_2+(-T20+162.+T490)*x1)*x1)*x1;
1202 *it_dwdx2i++ = fm45_2 + (594.+(fm8289_2+(11556.+(-13770.+T284)*x2)*x2)*x2)*x2 + (f99_2+(-1197.+(7128.+(-14904.+T87)*x2)*x2)*x2+T516)*x1;
1203 *it_dwdx2i++ = (f513_2+(-1071.+(1458.-T520)*x2)*x2+(-2610.+(9342.+(-10692.+T27)*x2)*x2+(7884.+T210+(-9396.+T79+T529)*x1)*x1)*x1)*x1;
1204 *it_dwdx2i++ = (4.+(-72.+T538)*x2+(-36.+(648.-T20)*x2+(72.+T543*x2)*x1)*x1)*x1;
1205 *it_dwdx2i++ = 40. + (-1016.+(6696+(-17424.+(19440.-T164)*x2)*x2)*x2)*x2 + (-148.+(3384.+(-18360.+(33696.-T137)*x2)*x2)*x2+(180.+(-3672.+15552.*x2*T66)*x2+(-72.-T543*x2)*x1)*x1)*x1;
1206 *it_dwdx2i++ = (-148.+ (360.-T538)*x2+(1692.+(-3672.+T20)*x2+(-6120.+T74+(8424.-T164-T529)*x1)*x1)*x1)*x1;
1207 *it_dwdx2i++ = (f9_2+(-99.+(486.-T520)*x2)*x2+T516)*x1;
1208 *it_dwdx2i++ = -45. + (1053.+(fm12447_2+(14796.+(-15390.+T284)*x2)*x2)*x2)*x2 + (f513_2+(-5220.+(23652.+(-37584.+T137)*x2)*x2)*x2+(fm1071_2+(9342.+(-29160.+T134)*x2)*x2+(486.+(-7128.+T79)*x2+T490)*x1)*x1)*x1;
1209 *it_dwdx2i++ = (f99_2-T488+(fm1197_2+T226+(2376.-T20+(-3726.+T20+T479)*x1)*x1)*x1)*x1;
1210 *it_dwdx2i++ = (f36_5+(-180.+(1134.-T470*x2)*x2)*x2)*x1;
1211 *it_dwdx2i++ = 36. + (fm3132_5+(3132.+(-6696.+(6480.-T461)*x2)*x2)*x2)*x2 + (fm1566_5+(4176.+(-15066.+(20736.-T87)*x2)*x2)*x2+(1044.+(-10044.+(23328.-T84)*x2)*x2+(-1674.+(10368.-T79)*x2+(1296.-T27-T1)*x1)*x1)*x1)*x1;
1212 *it_dwdx2i++ = -T460;
1213 *it_dwdx2i++ = (-54.+108.*x2+(-T281+648.+(T27-2538.+(-T27+3888.-T479)*x1)*x1)*x1)*x1;
1214 *it_dwdx2i++ = (-54.+(1296.+(-7614.+(15552.-T87)*x2)*x2)*x2+(54.+(-1188.+(5832.-T164)*x2)*x2)*x1)*x1;
1215 *it_dwdx2i++ = (540.+(-6156.+(19278.+(-23328.+T87)*x2)*x2)*x2+(-3078.+(25704.+(-52488.+T315)*x2)*x2+(6426.-34992*x2*T66+(T84-5832.+T479)*x1)*x1)*x1)*x1;
1216 *it_dwdx2i++ = (-36.+(504.-T520)*x2+(360.+(-4968.+T284)*x2+(-972.+(12960.-T79)*x2+T714)*x1)*x1)*x1;
1217 *it_dwdx2i++ = (180.+(-4032.+(21060.+(-36288.+T137)*x2)*x2)*x2+(-396.+(7992.+(-33048.+T315)*x2)*x2+(216.+T731*x2)*x1)*x1)*x1;
1218 *it_dwdx2i++ = (-360.+(2664.+(-4860.+T46)*x2)*x2+(3492.+(-22464.-T304*x2)*x2+(-9612.+(42768.-T310)*x2+(-T134+10368.-T529)*x1)*x1)*x1)*x1;
1219 *it_dwdx2i++ = (-36.+(720.+(-2916.+T46)*x2)*x2+(252.+(-4968.+(19440.-T84)*x2)*x2+(-216.-T731*x2)*x1)*x1)*x1;
1220 *it_dwdx2i++ = (-360.+(6984.+(-28836.+(41472.-T137)*x2)*x2)*x2+(1332.+(-22464.+(64152.-T366)*x2)*x2+(-1620.+(23328.-T310)*x2+T714)*x1)*x1)*x1;
1221 *it_dwdx2i++ = (180.+(-792.+T520)*x2+(-2016.+(7992.-T284)*x2+(7020.+(-22032.+T79)*x2+(T84-9072.+T529)*x1)*x1)*x1)*x1;
1222 *it_dwdx2i = (135.+(-2214.+(5832.-T27)*x2)*x2+(-1107.+(17496.+(-40824.+T134)*x2)*x2+(1944.+(-27216.+T310)*x2+(T79-972.)*x1)*x1)*x1)*x1;
1223 }
1224 }
1225
1226 template class LagrangeStdTriangle<_P0>;
1227 template class LagrangeStdTriangle<_P1>;
1228 template class LagrangeStdTriangle<_P2>;
1229 template class LagrangeStdTriangle<_P3>;
1230 template class LagrangeStdTriangle<_P4>;
1231 template class LagrangeStdTriangle<_P5>;
1232 template class LagrangeStdTriangle<_P6>;
1233
1234 /*
1235 --------------------------------------------------------------------------------
1236 Lagrange standard P1+bubbleP3 triangle Reference Element
1237 --------------------------------------------------------------------------------
1238 */
1239 template<>
LagrangeStdTriangle(const Interpolation * interp_p)1240 LagrangeStdTriangle<_P1BubbleP3>::LagrangeStdTriangle(const Interpolation* interp_p)
1241 : LagrangeTriangle(interp_p)
1242 {
1243 name_ += "_1+BubbleP3";
1244 // local coordinates of Points supporting D.o.F
1245 pointCoordinates();
1246 // build O1 splitting scheme
1247 splitO1Scheme = splitO1();
1248 }
1249
1250 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const1251 void LagrangeStdTriangle<_P1BubbleP3>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
1252 {
1253 real_t x1 = *it_pt, x2 = *(it_pt + 1), x3 = 1. - x1 - x2;
1254
1255 std::vector<real_t>::iterator it_w(shapeValues.w.begin());
1256 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2(it_dwdx1 + 1);
1257 std::vector<real_t>::iterator it_dwdx1i((*it_dwdx1).begin()), it_dwdx2i((*it_dwdx2).begin());
1258 real_t Ninex1x2x3 = 9.*x1*x2*x3;
1259 *it_w++ = x1 - Ninex1x2x3;
1260 *it_w++ = x2 - Ninex1x2x3;
1261 *it_w++ = x3 - Ninex1x2x3;
1262 *it_w++ = 3.*Ninex1x2x3;
1263
1264 if(withDeriv)
1265 {
1266 real_t d1Ninex1x2x3 = 9.*x2*(x3 - x1), d2Ninex1x2x3 = 9.*x1*(x3 - x2);
1267 *it_dwdx1i++ = 1. - d1Ninex1x2x3; *it_dwdx2i++ = -d2Ninex1x2x3;
1268 *it_dwdx1i++ = -d1Ninex1x2x3; *it_dwdx2i++ = 1. - d2Ninex1x2x3;
1269 *it_dwdx1i++ = -1. - d1Ninex1x2x3; *it_dwdx2i++ = -1. - d2Ninex1x2x3;
1270 *it_dwdx1i = 3.*d1Ninex1x2x3; *it_dwdx2i++ = 3.*d2Ninex1x2x3;
1271 }
1272 }
1273
1274
1275 /*!
1276 point number p has coordinates defined by ( x = coords[s2q[2*p]], y = coords[s2q[2*p+1]] ) where coords are coordinates of the associated 1D Lagrange reference element
1277
1278 Examples of local numbering of Lagrange Pk 1D elements
1279 2---1 2---3---1 2---4---3---1 2---5---4---3---1 2---6---5---4---3---1 2---7---6---5---4---3---1
1280 . k=1 k=2 k=3 k=4 k=5 k=6
1281 Examples of local numbering of Lagrange Pk triangles
1282
1283 2 2 2 2 2 2
1284 | \ | \ | \ | \ | \ | \
1285 3---1 5 4 5 7 5 10 5 13 5 16
1286 . k=1 | \ | \ | \ | \ | \
1287 ....... 3---6---1 8 10 4 8 14 7 8 17 10 8 20 13
1288 .......... k=2 | \ | | \ \ | | \ \ | | \ \
1289 .................. 3---6---9---1 11 15--13 4 11 20 19 7 11 23 25 10
1290 ........................ k=3 | \ | | \ \ | | \ \
1291 ................................. 3---6---9--12---1 14 18--21--16 4 14 26 28 22 7
1292 .......................................... k=4 | \ | | \ \
1293 .................................................... 3---6---9--12--15---1 17 21--24--27--19 4
1294 ............................................................. k=5 | \
1295 ........................................................................... 3---6---9--12--15--18---1
1296 ...................................................................................... k=6
1297 */
tensorNumberingTriangle(const int interpNum,std::vector<number_t> & s2t)1298 void tensorNumberingTriangle(const int interpNum, std::vector<number_t>& s2t)
1299 {
1300 std::vector<number_t>::iterator p=s2t.begin();
1301 int nk=interpNum;
1302 number_t q1 = 0, q2 = 1, q3 = 2, q4 = interpNum, q3_p, q4_m;
1303
1304 while(nk > 0)
1305 {
1306 // 'current perimeter' vertices
1307 (*p++) = q1; (*p++) = q2;
1308 (*p++) = q2; (*p++) = q1;
1309 (*p++) = q2; (*p++) = q2;
1310
1311 // non-vertex points on edges of current 'perimeter'
1312 q3_p = q3; q4_m = q4;
1313 for(dimen_t j = 2; j <= nk; j++)
1314 {
1315 (*p++) = q3_p; (*p++) = q4_m;
1316 (*p++) = q2; (*p++) = q3_p++;
1317 (*p++) = q4_m--; (*p++) = q2;
1318 }
1319
1320 q1 += 2; q2--; q3 += 2 ; q4--;
1321 if(nk == interpNum) { q1++; q2 = interpNum; }
1322 nk -= 3;
1323 }
1324 // central point if any
1325 if(nk == 0) { (*p++) = q1; (*p) = q1; }
1326 }
1327
1328 /*================================================================================================
1329 Lagrange standard Pk triangle (any order)
1330 ================================================================================================*/
1331
LagrangeStdTrianglePk(const Interpolation * interp_p)1332 LagrangeStdTrianglePk::LagrangeStdTrianglePk(const Interpolation* interp_p)
1333 : LagrangeTriangle(interp_p)
1334 {
1335 name_ += "_" + tostring(interp_p->numtype);
1336 pointCoordinates(); // local coordinates of Points supporting D.o.F
1337 computeShapeFunctions(); //construct shape functions as polynomials
1338 buildPolynomialTree(); //construct tree representation
1339 // build O1 splitting scheme
1340 splitO1Scheme = splitO1();
1341 }
1342
~LagrangeStdTrianglePk()1343 LagrangeStdTrianglePk::~LagrangeStdTrianglePk() {}
1344
1345
1346 // // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
1347 // #ifdef FIG4TEX_OUTPUT
1348 // /*
1349 // --------------------------------------------------------------------------------
1350 // fig4TeX builds a tex file to display local numbering and coordinates
1351 // --------------------------------------------------------------------------------
1352 // */
1353 // void LagrangeTriangle::fig4TeX()
1354 // {
1355 // Number nbVertices = geomRefElem_p->nbVertices();
1356 // char b = char(92);
1357
1358 // std::ofstream os;
1359 // os.open((fig4TexFileName()).c_str(), ios::out);
1360 // os << "%" << b << "input fig4tex.tex" << b << "input TeXnicolor.tex" << b << "input Fonts.tex";
1361 // os << std::endl << b << "fourteenpoint" << b << "nopagenumbers";
1362 // os << std::endl << b << "newbox" << b << "mybox";
1363 // if ( nbDofs(1, 1) < 6 ) { os << std::endl << b << "figinit{150pt}"; }
1364 // else { os << std::endl << b << "figinit{250pt}"; }
1365 // string longName = "Lagrange Triangle of degree ";
1366
1367 // // Draw element and display vertices
1368 // number_t pts(0);
1369
1370 // // define vertices
1371 // fig4TeXVertexPt(os, geomRefElem_p->vertices(), pts);
1372
1373 // // create postscript file
1374 // os << std::endl << b << "psbeginfig{}" << std::endl << b << "pssetwidth{1}" << b << "psline[";
1375 // for ( size_t p = 1; p <= nbVertices; p++ ) { os << p << ","; }
1376 // os << "1]";
1377 // os << std::endl << b << "psendfig"; // end of postcript output
1378
1379 // os << std::endl << b << "figvisu{" << b << "mybox}{" << longName << interpolation_p->number() << "}{%"
1380 // << std::endl << b << "figsetmark{$" << b << "bullet$}" << b << "figsetptname{{" << b << "bf#1}}";
1381 // // vertex nodes
1382 // os << std::endl << b << "figwritese1:(5pt)"
1383 // << std::endl << b << "figwritenw2:(5pt)"
1384 // << std::endl << b << "figwritesw3:(5pt)";
1385
1386 // // edge vertices
1387 // os << std::endl << b << "figsetmark{.}"
1388 // << std::endl << b << "color{" << b << "Red}"
1389 // << std::endl << b << "figwritegw1:${}_2$(7pt,7pt)" << std::endl << b << "figwritegw1:${}_1$(12pt,1pt)"
1390 // << std::endl << b << "figwritege2:${}_2$(1pt,-12pt)" << std::endl << b << "figwritege2:${}_1$(7pt,-7pt)"
1391 // << std::endl << b << "figwritege3:${}_1$(1pt,7pt)" << std::endl << b << "figwritege3:${}_2$(7pt,1pt)";
1392 // os << std::endl << b << "endcolor";
1393 // // nodes on edge 1 (north-east), edge 2 (west), edge 3 (south)
1394 // number_t ed(0);
1395 // fig4TeXEdgePt(os , ed++, pts, "ne", "sw");
1396 // fig4TeXEdgePt(os , ed++, pts, "w", "e");
1397 // fig4TeXEdgePt(os , ed++, pts, "s", "n");
1398
1399 // // internal nodes
1400 // if ( nb_internal_dof_ > 0 ) { fig4TeXInternalPoints2d(os, pts); }
1401
1402 // // close figvisu and TeX file
1403 // os << std::endl << "}" << b << "par" << b << "centerline{" << b << "box" << b << "mybox}";
1404 // os << std::endl << b << "bye";
1405 // os.close();
1406 // }
1407 // #endif /* FIG4TEX_OUTPUT */
1408 // // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
1409
1410 } // end of namespace xlifepp
1411