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 LagrangeTetrahedron.cpp
19 \authors D. Martin, N. Kielbasiewicz
20 \since 17 dec 2002
21 \date 6 aug 2012
22
23 \brief Implementation of xlifepp::LagrangeTetrahedron class members and related functions
24
25 xlifepp::LagrangeTetrahedron defines Lagrange Reference Element interpolation data
26 xlifepp::LagrangeTetrahedron is child to xlifepp::RefTetrahedron parent to template xlifepp::LagrangeStdTetrahedron<Pk>
27 */
28
29 #include "LagrangeTetrahedron.hpp"
30 #include "../Interpolation.hpp"
31 #include "../splitUtils.hpp"
32 #include "utils.h"
33 #include <iostream>
34
35 namespace xlifepp
36 {
37
38 //! LagrangeTetrahedron constructor for Lagrange reference elements
LagrangeTetrahedron(const Interpolation * interp_p)39 LagrangeTetrahedron::LagrangeTetrahedron(const Interpolation* interp_p)
40 : RefTetrahedron(interp_p)
41 {
42 name_ += "_Lagrange";
43 trace_p->push("LagrangeTetrahedron::LagrangeTetrahedron (" + name_ + ")");
44 // build element interpolation data
45 interpolationData();
46 // define local numbering of points on edges
47 sideOfSideNumbering();
48 // find or create Reference Element for element edges
49 sideOfSideRefElement();
50 // define local numbering of points on sides
51 sideNumbering();
52 // find or create Reference Element for element sides
53 sideRefElement();
54 //build barycentricSideDofMap : internal dof on side (face) ordered by barycentric coordinates
55 buildBarycentricSideDof(); //no effect if order < 3
56 maxDegree = interpolation_p->numtype;
57
58 trace_p->pop();
59 }
60
~LagrangeTetrahedron()61 LagrangeTetrahedron::~LagrangeTetrahedron() {}
62
63 //! interpolationData defines Reference Element interpolation data
interpolationData()64 void LagrangeTetrahedron::interpolationData()
65 {
66 trace_p->push("LagrangeTetrahedron::interpolationData");
67
68 number_t interpNum = interpolation_p->numtype;
69 switch (interpNum)
70 {
71 case _P0:
72 nbDofs_ = nbInternalDofs_ = 1;
73 break;
74 case _P1:
75 nbDofs_ = nbDofsOnVertices_ = 4;
76 break;
77 default: // Lagrange standard tetrahedron of degree > 1
78 nbDofs_ = ((interpNum + 1)*(interpNum + 2)*(interpNum + 3)) / 6;
79 nbDofsOnVertices_ = 4;
80 nbInternalDofs_ = ((interpNum - 3)*(interpNum - 2)*(interpNum - 1)) / 6;
81 nbDofsInSideOfSides_= 6*(interpNum - 1);
82 nbDofsInSides_ = nbDofs_ - nbInternalDofs_ - nbDofsOnVertices_ - nbDofsInSideOfSides_;
83 break;
84 }
85 // Creating Degrees Of Freedom of reference element
86 // in Lagrange standard elements all D.o.F on sides are sharable
87 // number of internal nodes (or D.o.F)
88 refDofs.reserve(nbDofs_);
89 lagrangeRefDofs(3, nbDofsOnVertices_, nbInternalDofs_, 6, nbDofsInSideOfSides_, 4, nbDofsInSides_);
90 //lagrangeRefDofs(1, nbDofsOnVertices_, nbDofs_ - nbInternalDofs_, 3);
91 // allocate storage for shape functions
92 shapeValues.resize(*this, 3);
93 nbPts_ = nbDofs_;
94 trace_p->pop();
95 }
96
97 /*!internal side dofs mapping when vertices of side are permuted
98 this function returns the number of the n-th internal dof of a face where vertices are permuted
99 if no permutation, it returns n
100
101 Numbering of internal side dofs of Pk tetrahedron, p= k-3
102
103 1 2 2 2 2 2 2
104 p=0 | \ | \ | \ | \ | \ | \
105 3---1 5 4 5 7 5 10 5 13 5 16
106 p=1 | \ | \ | \ | \ | \
107 3---6---1 8 10 4 8 14 7 8 17 10 8 20 13
108 p=2 | \ | | \ \ | | \ \ | | \ \
109 3---6---9---1 11 15--13 4 11 20 19 7 11 23 25 10
110 p=3 | \ | | \ \ | | \ \
111 3---6---9--12---1 14 18--21--16 4 14 26 28 22 7
112 p=4 | \ | | \ \
113 3---6---9--12--15---1 17 21--24--27--19 4
114 p=5 | \
115 3---6---9--12--15--18---1
116 p=6
117 Map exemple for P4 tetrahedron
118
119 2 1
120 | \ | \
121 | \ | \
122 | \ | \ sideDofsMap(1,3,1,2)=2
123 | 2 \ | 1 \ sideDofsMap(2,3,1,2)=3
124 | | \ \ | | \ \ sideDofsMap(3,3,1,2)=1
125 | 3--1 \ | 2--3 \
126 | \ | \
127 3---------------1 2--------------3
128
129 standard face permuted face
130
131 */
sideDofsMap(const number_t & n,const number_t & i,const number_t & j,const number_t & k) const132 number_t LagrangeTetrahedron::sideDofsMap(const number_t& n, const number_t& i, const number_t& j, const number_t& k) const
133 {
134 if (i==1 && j==2 && k==3) return n; //case (1,2,3)
135 //permute barycentric coordinate of dof n
136 Triplet<number_t> t=barycentricSideDofs[n-1], tp=t;
137 switch (i)
138 {
139 case 1 :
140 tp.second=t.third;
141 tp.third=t.second; // case (1,3,2)
142 break;
143 case 2:
144 tp.first=t.second;
145 if (j==1) tp.second=t.first; // case (2,1,3)
146 else
147 {
148 tp.second=t.third; // case (2,3,1)
149 tp.third=t.first;
150 }
151 break;
152 case 3:
153 tp.first=t.third;
154 if (j==1)
155 {
156 tp.second=t.first; // case (3,1,2)
157 tp.third=t.second;
158 }
159 else tp.third=t.first; // case (3,2,1)
160 break;
161 default :
162 where("LagrangeTetrahedron::sideDofsMap()");
163 error("index_out_of_range","i", 1, 3);
164 }
165 //find new dof number after permutation
166 std::map<Triplet<number_t>,number_t>::const_iterator it=barycentricSideDofMap.find(tp);
167 if (it==barycentricSideDofMap.end())
168 {
169 where("LagrangeTetrahedron::sideDofsMap()");
170 error("triplet_not_found");
171 }
172 return it->second;
173 }
174
175 //create the map of internal dofs on side, ordered by barycentric coordinates
buildBarycentricSideDof()176 void LagrangeTetrahedron::buildBarycentricSideDof()
177 {
178 number_t k=interpolation_p->numtype;
179 if (k<3) return; //no internal side dofs
180 number_t n=k-3; //number of layers
181 number_t i=0; //dof index
182 number_t A=n+1, B=1,C=1;
183 for (int p=n; p>=0; p-=3, A-=2, B++, C++) //layer loop
184 {
185 number_t a=A, b=B,c=C;
186 for (number_t j=0; int_t(j)<=p; j++, a--, b++) //dof loop
187 {
188 Triplet<number_t> t(a,b,c);
189 if (barycentricSideDofMap.find(t)==barycentricSideDofMap.end()) barycentricSideDofMap.insert(std::make_pair(t,++i));
190 t=Triplet<number_t>(c,a,b);
191 if (barycentricSideDofMap.find(t)==barycentricSideDofMap.end()) barycentricSideDofMap.insert(std::make_pair(t,++i));
192 t=Triplet<number_t>(b,c,a);
193 if (barycentricSideDofMap.find(t)==barycentricSideDofMap.end()) barycentricSideDofMap.insert(std::make_pair(t,++i));
194 }
195 }
196 if (i>0)
197 {
198 barycentricSideDofs.resize(i);
199 std::vector<Triplet<number_t> >::iterator itvb=barycentricSideDofs.begin();
200 std::map<Triplet<number_t>,number_t>::iterator itm=barycentricSideDofMap.begin(), itme=barycentricSideDofMap.end();
201 for (; itm!=itme; ++itm) *(itvb + (itm->second-1)) = itm->first;
202 }
203 }
204
205 //! output of Lagrange D.o.f's numbers on underlying tetrahedron P1 D.o.f's
outputAsP1(std::ofstream & os,const int refNum,const std::vector<number_t> & rks) const206 void LagrangeTetrahedron::outputAsP1(std::ofstream& os, const int refNum, const std::vector<number_t>& rks) const
207 {
208 //const number_t w(9);
209 number_t interpNum = interpolation_p->numtype;
210 switch (interpNum)
211 {
212 case _P0:
213 case _P1:
214 simplexVertexOutput(os, refNum, rks[0], rks[1], rks[2], rks[3]);
215 break;
216 case _P2:
217 // 4 half-size tetrahedra with nodes ( 1,10, 9, 5), (10, 2, 8, 6), ( 9, 8, 3, 7), ( 5, 6, 7, 4)
218 simplexVertexOutput(os, refNum, rks[0], rks[9], rks[8], rks[4]);
219 simplexVertexOutput(os, refNum, rks[9], rks[1], rks[7], rks[5]);
220 simplexVertexOutput(os, refNum, rks[8], rks[7], rks[2], rks[6]);
221 simplexVertexOutput(os, refNum, rks[4], rks[5], rks[6], rks[3]);
222 // 4 tetrahedra forming 2 facing pyramids : ( 5, 8,10, 9), ( 5, 8, 7, 6), ( 8, 5, 7, 9), ( 8, 5,10, 6);
223 simplexVertexOutput(os, refNum, rks[4], rks[7], rks[9], rks[8]);
224 simplexVertexOutput(os, refNum, rks[4], rks[7], rks[6], rks[5]);
225 simplexVertexOutput(os, refNum, rks[7], rks[4], rks[6], rks[8]);
226 simplexVertexOutput(os, refNum, rks[7], rks[4], rks[9], rks[5]);
227 break;
228 default:
229 noSuchFunction("outputAsP1");
230 break;
231 }
232 }
233
234 //! pointCoordinates defines Lagrange Reference Element point coordinates
vertexCoordinates()235 std::vector<RefDof*>::iterator LagrangeTetrahedron::vertexCoordinates()
236 {
237 std::vector<RefDof*>::iterator it_rd = refDofs.begin();
238 (*it_rd++)->coords(1., 0., 0.);
239 (*it_rd++)->coords(0., 1., 0.);
240 (*it_rd++)->coords(0., 0., 1.);
241 (*it_rd++)->coords(0., 0., 0.);
242 return it_rd;
243 }
244 //NEW VERSION of point coordinates calculation - compliant with face dofs ordering -
245 //does not used the tensorNumberingTetrahedron function
246 //checked up to degree 7
pointCoordinates()247 void LagrangeTetrahedron::pointCoordinates()
248 {
249 trace_p->push("LagrangeTetrahedron::pointCoordinates");
250 // Compute D.o.F's from the correspondence between segment and tetrahedron local numbering
251 number_t interpNum = interpolation_p->numtype;
252 std::vector<RefDof*>::iterator it_rd = refDofs.begin();
253
254 switch (interpNum)
255 {
256 case _P0:
257 (*it_rd)->coords(.25, .25, .25);
258 break;
259 case _P1BubbleP3:
260 it_rd = vertexCoordinates();
261 (*it_rd++)->coords(.25, .25, .25);
262 break;
263 default: // Lagrange standard tetrahedron of degree > 0
264 std::vector<Point> vs(4,Point(0.,0.,0.));
265 vs[0]=Point(1.,0.,0.); vs[1]=Point(0.,1.,0.); vs[2]=Point(0.,0.,1.);
266 for (number_t i=0;i<4;i++) (*it_rd++)->coords(vs[i]);
267 if (interpNum > 1)
268 {
269 for (number_t k=1; k<=interpNum-1; k++)
270 {
271 real_t a=real_t(k)/interpNum;
272 for (number_t e=1; e<=6; e++) //edge dofs
273 (*it_rd++)->coords((1-a)*vs[geomRefElem_p->sideOfSideVertexNumber(1,e)-1]+a*vs[geomRefElem_p->sideOfSideVertexNumber(2,e)-1]);
274 }
275 }
276 if (interpNum < 3) { trace_p->pop(); return; } //no face nodes
277 for (number_t f=1; f<=4; f++) //face dofs
278 {
279 Point &v1=vs[geomRefElem_p->sideVertexNumber(1,f)-1],
280 &v2=vs[geomRefElem_p->sideVertexNumber(2,f)-1],
281 &v3=vs[geomRefElem_p->sideVertexNumber(3,f)-1];
282 number_t n=interpNum-3; //number of layers
283 number_t A=n+1, B=1, C=1;
284 for (int p=n; p>=0; p-=3, A-=2, B++, C++) //layer loop
285 {
286 number_t a=A, b=B,c=C;
287 if (p==0)(*it_rd++)->coords((v1+v2+v3)/3.);
288 else
289 for (number_t j=1; int_t(j)<=p; j++, a--, b++) //dof loop
290 {
291 (*it_rd++)->coords((a*v1+b*v2+c*v3)/interpNum);
292 (*it_rd++)->coords((c*v1+a*v2+b*v3)/interpNum);
293 (*it_rd++)->coords((b*v1+c*v2+a*v3)/interpNum);
294 }
295 }
296 }
297 if (interpNum < 4) { trace_p->pop(); return;} //no internal nodes
298 if (interpNum==4) //internal node for P4
299 {
300 (*it_rd)->coords(.25, .25, .25);
301 trace_p->pop();
302 return;
303 }
304 //internal nodes : Pk-4 nodes on tetrahedron (1-3/k, 1/k, 1/k ), (1/k, 1-3/k, 1/k), (1/k, 1/k, 1-3/k) , (1/k, 1/k, 1/k)
305 //when using general Pk Lagrange, ordering of internal dofs is indifferent because shape functions are computed from coords
306 Interpolation* intm1=findInterpolation(_Lagrange,_standard,interpNum-4,H1);
307 RefElement* pkm1 = findRefElement(_tetrahedron, intm1);
308 std::vector<RefDof*>::iterator itdm1=pkm1->refDofs.begin();
309 real_t b=1./interpNum, a=1-4*b;
310 for(;itdm1!=pkm1->refDofs.end(); itdm1++, it_rd++)
311 (*it_rd)->coords(a*(*itdm1)->point()+b); //x -> (1-4/k)x+1/k maps [0,1] to [1/k, 1-3/k]
312 break;
313 }
314 trace_p->pop();
315 }
316
317 //! sideNumbering defines Lagrange Reference Element numbering on sides
sideNumbering()318 void LagrangeTetrahedron::sideNumbering()
319 {
320 trace_p->push("LagrangeTetrahedron::sideNumbering");
321
322 number_t interpNum = interpolation_p->numtype;
323 number_t interpNum1 = interpNum - 1;
324 number_t nbSides = geomRefElem_p->nbSides();
325 number_t nbSideOfSides = geomRefElem_p->nbSideOfSides();
326 number_t nbVertices = geomRefElem_p->nbVertices();
327 number_t nbVerticesPerSide = geomRefElem_p->sideVertexNumbers()[0].size();
328 number_t nbDofsPerSideOfSide = interpNum + 1;
329 number_t nbDofsPerSide = ((interpNum + 1)*(interpNum + 2)) / 2;
330 number_t nbInternalDofsPerSide = ((interpNum - 2)*(interpNum - 1)) / 2;
331 number_t more = nbVertices + interpNum1*nbSideOfSides;
332
333 sideDofNumbers_.resize(nbSides);
334 if(interpNum==0) //Q0 interpolation
335 {
336 for(number_t side = 0; side < nbSides; side++)
337 {
338 sideDofNumbers_[side].resize(1);
339 sideDofNumbers_[side][0]=1;
340 }
341 trace_p->pop();
342 return;
343 }
344 if (interpNum > 0)
345 {
346 number_t sideNum = 0;
347 for (number_t side = 0; side < nbSides; side++, sideNum++)
348 {
349 number_t nbSideOfSidesPerSide = 3;
350 sideDofNumbers_[side].resize(nbDofsPerSide);
351 for (number_t i = 0; i < nbVerticesPerSide; i++)
352 sideDofNumbers_[side][i] = geomRefElem_p->sideVertexNumber(i + 1, side + 1);
353 number_t i = 0;
354 for (number_t sideOfSideNum = 0; sideOfSideNum < nbSideOfSidesPerSide ; sideOfSideNum++)
355 {
356 i = sideOfSideNum;
357 int_t sideOfSideInSideNum = geomRefElem_p->sideOfSideNumber(sideOfSideNum + 1, side + 1);
358 if (sideOfSideInSideNum > 0)
359 {
360 for (number_t e = 2; e < nbDofsPerSideOfSide; e++)
361 {
362 i += nbSideOfSidesPerSide;
363 sideDofNumbers_[side][i] = sideOfSideDofNumbers_[sideOfSideInSideNum - 1][e];
364 }
365 }
366 else
367 {
368 sideOfSideInSideNum = -sideOfSideInSideNum;
369 for (number_t e = nbDofsPerSideOfSide - 1; e > 1; e--)
370 {
371 i += nbSideOfSidesPerSide;
372 sideDofNumbers_[side][i] = sideOfSideDofNumbers_[sideOfSideInSideNum - 1][e];
373 }
374 }
375 }
376 for (number_t d = 0; d < nbInternalDofsPerSide; d++) // internal D.o.F in side
377 sideDofNumbers_[side][++i] = ++more;
378 }
379 }
380 trace_p->pop();
381 }
382
383 //! sideOfSideNumbering defines Lagrange Reference Element numbering on side of sides
sideOfSideNumbering()384 void LagrangeTetrahedron::sideOfSideNumbering()
385 {
386 trace_p->push("LagrangeTetrahedron::sideOfSideNumbering");
387 number_t interpNum = interpolation_p->numtype;
388
389 if (interpNum > 0)
390 {
391 number_t intN = interpNum;
392 if (interpNum == _P1BubbleP3) intN = 1;
393 number_t nbVertices = geomRefElem_p->nbVertices();
394 number_t nbSideOfSides = geomRefElem_p->nbSideOfSides();
395 number_t more;
396 number_t nbVerticesPerSideOfSide = geomRefElem_p->sideOfSideVertexNumbers()[0].size();
397 number_t nbDofsPerSideOfSide = intN + 1;
398 sideOfSideDofNumbers_.resize(nbSideOfSides);
399 for (number_t sideOfSide = 0; sideOfSide < nbSideOfSides; sideOfSide++)
400 {
401 sideOfSideDofNumbers_[sideOfSide].resize(nbDofsPerSideOfSide);
402 for (number_t j = 0; j < nbVerticesPerSideOfSide; j++)
403 sideOfSideDofNumbers_[sideOfSide][j] = geomRefElem_p->sideOfSideVertexNumber(j + 1, sideOfSide + 1);
404 if (intN > 1)
405 {
406 more = nbVertices + 1;
407 for (number_t k = nbVerticesPerSideOfSide; k < nbDofsPerSideOfSide; k++, more += nbSideOfSides)
408 sideOfSideDofNumbers_[sideOfSide][k] = sideOfSide + more;
409 }
410 }
411 }
412 trace_p->pop();
413 }
414
415 //! return nodes numbers of first order tetrahedron elements when splitting current element
splitP1() const416 std::vector<std::vector<number_t> > LagrangeTetrahedron::splitP1() const {
417 std::vector<std::vector<number_t> > splitnum, buf, buf2;
418 std::vector<number_t> tuple(4, 0);
419 std::vector<number_t> num;
420 switch(interpolation_p->numtype) {
421 case 0:
422 case 1:
423 num.resize(4);
424 for (number_t i=0; i<num.size(); i++) num[i]=i+1;
425 splitnum.push_back(num);
426 return splitnum;
427 case 2:
428 num.resize(10);
429 for (number_t i=0; i<num.size(); i++) num[i]=i+1;
430 return splitTetrahedronP2ToTetrahedraP1(num);
431 case 3:
432 num.resize(20);
433 for (number_t i=0; i<num.size(); i++) num[i]=i+1;
434 return splitTetrahedronP3ToTetrahedraP1(num);
435 default:
436 return splitLagrange3DToP1();
437 }
438 return splitnum;
439 }
440
441 //! return nodes numbers of first order tetrahedron elements when splitting current element
splitO1() const442 splitvec_t LagrangeTetrahedron::splitO1() const
443 {
444 splitvec_t splitnum;
445 std::vector<std::vector<number_t> > buf;
446 buf=this->splitP1();
447 for(number_t i=0; i< buf.size(); i++)
448 splitnum.push_back(std::pair<ShapeType, std::vector<number_t> >(_tetrahedron, buf[i]));
449 return splitnum;
450 }
451
452 //! Lagrange standard P0 tetrahedron Reference Element
453 template<>
LagrangeStdTetrahedron(const Interpolation * interp_p)454 LagrangeStdTetrahedron<_P0>::LagrangeStdTetrahedron(const Interpolation* interp_p)
455 : LagrangeTetrahedron(interp_p)
456 {
457 name_ += "_0";
458 pointCoordinates();
459 // build O1 splitting scheme
460 splitO1Scheme = splitO1();
461 }
462
463 //! Lagrange standard P1 tetrahedron Reference Element
464 template<>
LagrangeStdTetrahedron(const Interpolation * interp_p)465 LagrangeStdTetrahedron<_P1>::LagrangeStdTetrahedron(const Interpolation* interp_p)
466 : LagrangeTetrahedron(interp_p)
467 {
468 name_ += "_1";
469 pointCoordinates();
470 // build O1 splitting scheme
471 splitO1Scheme = splitO1();
472 }
473
474 //! Lagrange standard Pk tetrahedron Reference Element
475 template<number_t Pk>
LagrangeStdTetrahedron(const Interpolation * interp_p)476 LagrangeStdTetrahedron<Pk>::LagrangeStdTetrahedron(const Interpolation* interp_p)
477 : LagrangeTetrahedron(interp_p)
478 {
479 name_ += "_" + tostring(Pk);
480 pointCoordinates();
481 // build O1 splitting scheme
482 splitO1Scheme = splitO1();
483 //#ifdef FIG4TEX_OUTPUT
484 // // Draw a LagrangeStdTetrahedron using TeX and fig4TeX
485 // fig4TeX();
486 //#endif
487 }
488
489 //! computeShapeValues defines Lagrange Reference Element shape functions
490 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const491 void LagrangeStdTetrahedron<_P0>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
492 {
493 std::vector<real_t>::iterator sh0w_i=shv.w.begin();
494 *sh0w_i = 1.;
495 if (withDeriv)
496 {
497 std::vector<std::vector<real_t> >::iterator sh1w=shv.dw.begin(), sh2w=sh1w+1, sh3w=sh1w+2;
498 std::vector<real_t>::iterator sh1w_i=(*sh1w).begin(), sh2w_i=(*sh2w).begin(), sh3w_i=(*sh3w).begin();
499 *sh1w_i = 0.;
500 *sh2w_i = 0.;
501 *sh3w_i = 0.;
502 }
503 }
504
505 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const506 void LagrangeStdTetrahedron<_P1>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
507 {
508 std::vector<real_t>::iterator sh0w_i(shv.w.begin());
509 *sh0w_i++ = *it_pt;
510 *sh0w_i++ = *(it_pt+1);
511 *sh0w_i++ = *(it_pt+2);
512 *sh0w_i = 1-*it_pt-*(it_pt+1)-*(it_pt+2);
513
514 if (withDeriv)
515 {
516 std::vector< std::vector<real_t> >::iterator sh1w(shv.dw.begin()), sh2w(sh1w + 1), sh3w(sh1w + 2);
517 std::vector<real_t>::iterator sh1w_i((*sh1w).begin()), sh2w_i((*sh2w).begin()), sh3w_i((*sh3w).begin());
518 *sh1w_i++ = 1.;
519 *sh2w_i++ = 0.;
520 *sh3w_i++ = 0.;
521 *sh1w_i++ = 0.;
522 *sh2w_i++ = 1.;
523 *sh3w_i++ = 0.;
524 *sh1w_i++ = 0.;
525 *sh2w_i++ = 0.;
526 *sh3w_i++ = 1.;
527 *sh1w_i = -1.;
528 *sh2w_i = -1.;
529 *sh3w_i = -1.;
530 }
531 }
532
533 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const534 void LagrangeStdTetrahedron<_P2>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
535 {
536 real_t x1=*it_pt, x2=*(it_pt+1), x3=*(it_pt+2), x4=1.-x1-x2-x3;
537 real_t qx1=4.*x1, qx2=4.*x2, qx3=4.*x3, qx4=4.*x4, qx4m1=qx4-1.;
538 std::vector<real_t>::iterator sh0w_i = shv.w.begin();
539 *sh0w_i++ = x1*(x1+x1-1.);
540 *sh0w_i++ = x2*(x2+x2-1.);
541 *sh0w_i++ = x3*(x3+x3-1.);
542 *sh0w_i++ = x4*(x4+x4-1.);
543 *sh0w_i++ = qx1*x4;
544 *sh0w_i++ = qx2*x4;
545 *sh0w_i++ = qx3*x4;
546 *sh0w_i++ = qx2*x3;
547 *sh0w_i++ = qx1*x3;
548 *sh0w_i = qx1*x2;
549
550 if (withDeriv)
551 {
552 std::vector< std::vector<real_t> >::iterator sh1w(shv.dw.begin()), sh2w(sh1w + 1), sh3w(sh1w + 2);
553 std::vector<real_t>::iterator sh1w_i((*sh1w).begin()), sh2w_i((*sh2w).begin()), sh3w_i((*sh3w).begin());
554 *sh1w_i++ = qx1-1.; *sh2w_i++ = 0.; *sh3w_i++ = 0.;
555 *sh1w_i++ = 0.; *sh2w_i++ = qx2-1.; *sh3w_i++ = 0.;
556 *sh1w_i++ = 0.; *sh2w_i++ = 0.; *sh3w_i++ = qx3-1.;
557 *sh1w_i++ = -qx4m1; *sh2w_i++ = -qx4m1; *sh3w_i++ = -qx4m1;
558 *sh1w_i++ = qx4-qx1; *sh2w_i++ = -qx1; *sh3w_i++ = -qx1;
559 *sh1w_i++ = -qx2; *sh2w_i++ = qx4-qx2; *sh3w_i++ = -qx2;
560 *sh1w_i++ = -qx3; *sh2w_i++ = -qx3; *sh3w_i++ = qx4-qx3;
561 *sh1w_i++ = 0.; *sh2w_i++ = qx3; *sh3w_i++ = qx2;
562 *sh1w_i++ = qx3; *sh2w_i++ = 0.; *sh3w_i++ = qx1;
563 *sh1w_i = qx2; *sh2w_i = qx1; *sh3w_i = 0.;
564 }
565 }
566
567 template<>
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const568 void LagrangeStdTetrahedron<_P3>::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
569 {
570 real_t x1 = *it_pt, x2 = *(it_pt + 1), x3 = *(it_pt + 2), x4 = 1. - x1 - x2 - x3;
571 real_t tx1m1 = 3.*x1 - 1., tx2m1 = 3.*x2 - 1., tx3m1 = 3.*x3 - 1., tx4m1 = 3.*x4 - 1.;
572 std::vector<real_t>::iterator sh0w_i = shv.w.begin();
573 *sh0w_i++ = 0.5*x1*tx1m1*(tx1m1-1.);
574 *sh0w_i++ = 0.5*x2*tx2m1*(tx2m1-1.);
575 *sh0w_i++ = 0.5*x3*tx3m1*(tx3m1-1.);
576 *sh0w_i++ = 0.5*x4*tx4m1*(tx4m1-1.);
577 *sh0w_i++ = 4.5*x1*tx1m1*x4;
578 *sh0w_i++ = 4.5*x2*tx2m1*x4;
579 *sh0w_i++ = 4.5*x3*tx3m1*x4;
580 *sh0w_i++ = 4.5*x2*tx2m1*x3;
581 *sh0w_i++ = 4.5*x1*tx1m1*x3;
582 *sh0w_i++ = 4.5*x1*tx1m1*x2;
583 *sh0w_i++ = 4.5*x4*tx4m1*x1;
584 *sh0w_i++ = 4.5*x4*tx4m1*x2;
585 *sh0w_i++ = 4.5*x4*tx4m1*x3;
586 *sh0w_i++ = 4.5*x3*tx3m1*x2;
587 *sh0w_i++ = 4.5*x3*tx3m1*x1;
588 *sh0w_i++ = 4.5*x2*tx2m1*x1;
589 *sh0w_i++ = 27.*x2*x3*x4;
590 *sh0w_i++ = 27.*x1*x3*x4;
591 *sh0w_i++ = 27.*x1*x2*x4;
592 *sh0w_i = 27.*x1*x2*x3;
593
594 if (withDeriv)
595 {
596 std::vector< std::vector<real_t> >::iterator sh1w(shv.dw.begin()), sh2w(sh1w + 1), sh3w(sh1w + 2);
597 std::vector<real_t>::iterator sh1w_i((*sh1w).begin()), sh2w_i((*sh2w).begin()), sh3w_i((*sh3w).begin());
598 *sh1w_i++ = 13.5*x1*x1 - 9.*x1 + 1.; *sh2w_i++ = 0.; *sh3w_i++ = 0.;
599 *sh1w_i++ = 0.; *sh2w_i++ = 13.5*x2*x2 - 9.*x2 + 1.; *sh3w_i++ = 0.;
600 *sh1w_i++ = 0.; *sh2w_i++ = 0.; *sh3w_i++ = 13.5*x3*x3 - 9.*x3 + 1.;
601 *sh1w_i++ = -13.5*x4*x4 + 9.*x4 - 1.; *sh2w_i++ = -13.5*x4*x4 + 9.*x4 - 1.; *sh3w_i++ = -13.5*x4*x4 + 9.*x4 - 1.;
602 *sh1w_i++ = -4.5*(3.*x1*x1-(1.+6.*x4)*x1+x4); *sh2w_i++ = -4.5*x1*tx1m1; *sh3w_i++ = -4.5*x1*tx1m1;
603 *sh1w_i++ = -4.5*x2*tx2m1; *sh2w_i++ = -4.5*(3.*x2*x2-(1.+6.*x4)*x2+x4); *sh3w_i++ = -4.5*x2*tx2m1;
604 *sh1w_i++ = -4.5*x3*tx3m1; *sh2w_i++ = -4.5*x3*tx3m1; *sh3w_i++ = -4.5*(3.*x3*x3-(1.+6.*x4)*x3+x4);
605 *sh1w_i++ = 0.; *sh2w_i++ = 4.5*x3*(6.*x2-1.); *sh3w_i++ = 4.5*x2*tx2m1;
606 *sh1w_i++ = 4.5*x3*(6.*x1-1.); *sh2w_i++ = 0.; *sh3w_i++ = 4.5*x1*tx1m1;
607 *sh1w_i++ = 4.5*x2*(6.*x1-1.); *sh2w_i++ = 4.5*x1*tx1m1; *sh3w_i++ = 0.;
608 *sh1w_i++ = 4.5*((1.-6.*x4)*x1+3.*x4*x4-x4); *sh2w_i++ = 4.5*x1*(1.-6.*x4); *sh3w_i++ = 4.5*x1*(1.-6.*x4);
609 *sh1w_i++ = 4.5*x2*(1.-6.*x4); *sh2w_i++ = 4.5*((1.-6.*x4)*x2+3.*x4*x4-x4); *sh3w_i++ = 4.5*x2*(1.-6.*x4);
610 *sh1w_i++ = 4.5*x3*(1.-6.*x4); *sh2w_i++ = 4.5*x3*(1.-6.*x4); *sh3w_i++ = 4.5*((1.-6.*x4)*x3+3.*x4*x4-x4);
611 *sh1w_i++ = 0.; *sh2w_i++ = 4.5*x3*tx3m1; *sh3w_i++ = 4.5*x2*(6.*x3-1.);
612 *sh1w_i++ = 4.5*x3*tx3m1; *sh2w_i++ = 0.; *sh3w_i++ = 4.5*x1*(6.*x3-1.);
613 *sh1w_i++ = 4.5*x2*tx2m1; *sh2w_i++ = 4.5*x1*(6.*x2-1.); *sh3w_i++ = 0.;
614 *sh1w_i++ = -27.*x2*x3; *sh2w_i++ = 27.*x3*(x4-x2); *sh3w_i++ = 27.*x2*(x4-x3);
615 *sh1w_i++ = 27.*x3*(x4-x1); *sh2w_i++ = -27.*x3*x1; *sh3w_i++ = 27.*x1*(x4-x3);
616 *sh1w_i++ = 27.*x2*(x4-x1); *sh2w_i++ = 27.*x1*(x4-x2); *sh3w_i++ = -27.*x1*x2;
617 *sh1w_i = 27.*x2*x3; *sh2w_i = 27.*x1*x3; *sh3w_i = 27.*x1*x2;
618 }
619 }
620
621 template class LagrangeStdTetrahedron<_P2>;
622 template class LagrangeStdTetrahedron<_P3>;
623
624 /*!
625 correspondance between segment and tetrahedron local numbering defined such that
626 - point number p has coordinates defined by
627 - ( x = coords[s2h[0][p]], y = coords[s2h[1][p]], z = coords[s2h[3][p]] )
628 where coords are coordinates of the associated 1D Lagrange reference element
629 ---- No longer used -----
630 */
tensorNumberingTetrahedron(const int interpNum,number_t ** & s2t)631 void tensorNumberingTetrahedron(const int interpNum, number_t**& s2t)
632 {
633 // CHECKED ONLY UP TO degree 7 !!!
634 int nk(interpNum), nkk, nk2;
635 number_t p=0;//, pa=0;
636 number_t q1=0, q2=1, q3=2, q4=interpNum, z1=1;
637 number_t q3_p, q4_m;
638
639 while (nk > 0)
640 {
641 // vertices of 'current perimeter'
642 s2t[0][p] = q1;
643 s2t[1][p] = q2;
644 s2t[2][p++] = q2; // 1
645 s2t[0][p] = q2;
646 s2t[1][p] = q1;
647 s2t[2][p++] = q2; // 2
648 s2t[0][p] = q2;
649 s2t[1][p] = q2;
650 s2t[2][p++] = q1; // 3
651 s2t[0][p] = q2;
652 s2t[1][p] = q2;
653 s2t[2][p++] = q2; // 4
654 // non-vertex points on edges of 'current perimeter'
655 q3_p = q3;
656 q4_m = q4;
657 for (int j = 2; j <= nk; j++)
658 {
659 s2t[0][p] = q3_p;
660 s2t[1][p] = q2;
661 s2t[2][p++] = q2; //Edge# 1->4, points # 5, 11, ...
662 s2t[0][p] = q2;
663 s2t[1][p] = q3_p;
664 s2t[2][p++] = q2; //Edge# 2->4, points # 6, 12, ...
665 s2t[0][p] = q2;
666 s2t[1][p] = q2;
667 s2t[2][p++] = q3_p; //Edge# 3->4, points # 7, 13, ...
668 s2t[0][p] = q2;
669 s2t[1][p] = q3_p;
670 s2t[2][p++] = q4_m; //Edge# 2->3, points # 8, 14, ...
671 s2t[0][p] = q3_p;
672 s2t[1][p] = q2;
673 s2t[2][p++] = q4_m; //Edge# 1->3, points # 9, 15, ...
674 s2t[0][p] = q3_p;
675 s2t[1][p] = q4_m;
676 s2t[2][p++] = q2; //Edge# 1->2, points # 10, 16, ...
677 q3_p++;
678 q4_m--;
679 }
680
681 nkk = nk - 3;
682
683 // for ( int i = pa; i < p; i++ ) { thePrintStream_ <<std::endl<<" ep "<< i+1 <<" ("<< s2t[0][i] << ", "<<s2t[1][i]<<", "<<s2t[2][i]<<")"; }
684 // pa = p;
685 // 'current perimeter' internal points on sides
686 number_t x1=q4, x2=q3 + 1, x1m, x2p;
687 if (nkk >= 0)
688 {
689 // Face #1 x = "0"
690 tensorTetrahedronSideNumbering(nkk, 0, 1, 2, z1, x1, x2, s2t, p);
691 // Face #2 y = "0"
692 tensorTetrahedronSideNumbering(nkk, 1, 2, 0, z1, x1, x2, s2t, p);
693 // Face #3 z = "0"
694 tensorTetrahedronSideNumbering(nkk, 2, 0, 1, z1, x1, x2, s2t, p);
695 // Face #4 x+y+z = "1"
696 // pa = p;
697 nk2 = nkk;
698 while (nk2 > 0)
699 {
700 x1m = x1;
701 x2p = x2;
702 for (dimen_t j = 1; j <= nk2; j++)
703 {
704 s2t[0][p] = x2p;
705 s2t[1][p] = x1 ;
706 s2t[2][p++] = x1m;
707 s2t[0][p] = x1 ;
708 s2t[1][p] = x1m;
709 s2t[2][p++] = x2p;
710 s2t[0][p] = x1m;
711 s2t[1][p] = x2p;
712 s2t[2][p++] = x1;
713 x1m--;
714 x2p++;
715 }
716 x1--;
717 x2 += 2 ;
718 nk2 -= 3;
719 }
720 if (nk2 == 0)
721 {
722 s2t[0][p] = x1;
723 s2t[1][p] = x1;
724 s2t[2][p++] = x1;
725 }
726 }
727 q1 += 3;
728 q2--;
729 q3 += 3;
730 q4--;
731 z1--;
732 if (nk == interpNum)
733 {
734 q1 = 4;
735 q2 = nk;
736 z1 = nk;
737 }
738 nk -= 4;
739 }
740
741 // central point in element if any
742 if (nk == 0)
743 {
744 s2t[0][p] = q1;
745 s2t[1][p] = q1;
746 s2t[2][p] = q1;
747 }
748 }
749
tensorTetrahedronSideNumbering(const int nk,const dimen_t i0,const dimen_t i1,const dimen_t i2,const number_t z1,const number_t xx1,const number_t xx2,number_t ** & s2t,number_t & p)750 void tensorTetrahedronSideNumbering(const int nk, const dimen_t i0, const dimen_t i1, const dimen_t i2, const number_t z1, const number_t xx1, const number_t xx2, number_t**& s2t, number_t& p)
751 {
752 // Define correspondence between segment and local numbering of internal D.o.F's
753 // on side of the Lagrange Pk tetrahedron in the same order as the order of numbering
754 // of internal D.o.F's of the Lagrange Pk triangle
755 int nk2(nk), x1(xx1), x2(xx2), x1m, x2p;//, pa(p);
756 while (nk2 > 0)
757 {
758 // points on vertices and edges of current 'perimeter' on side
759 x1m = x1;
760 x2p = x2;
761 for (int j = 1; j <= nk2; j++)
762 {
763 s2t[i0][p] = z1;
764 s2t[i1][p] = x1;
765 s2t[i2][p++] = x1m;
766 s2t[i0][p] = z1;
767 s2t[i1][p] = x1m;
768 s2t[i2][p++] = x2p;
769 s2t[i0][p] = z1;
770 s2t[i1][p] = x2p;
771 s2t[i2][p++] = x1;
772 x1m--;
773 x2p++;
774 }
775 x1--;
776 x2 += 2;
777 nk2 -= 3;
778 }
779 // central point on side if any
780 if (nk2 == 0)
781 {
782 s2t[i0][p] = z1;
783 s2t[i1][p] = x1;
784 s2t[i2][p++] = x1;
785 }
786 }
787
788
789 /*================================================================================================
790 Lagrange standard Pk tetrahedron (any order) using formal polynomials to compute shape functions
791 ================================================================================================*/
792
LagrangeStdTetrahedronPk(const Interpolation * interp_p)793 LagrangeStdTetrahedronPk::LagrangeStdTetrahedronPk(const Interpolation* interp_p)
794 : LagrangeTetrahedron(interp_p)
795 {
796 name_ += "_" + tostring(interp_p->numtype);
797 pointCoordinates(); // local coordinates of Points supporting D.o.F
798 initShapeFunCoeffs();
799 // build O1 splitting scheme
800 splitO1Scheme = splitO1();
801 }
802
~LagrangeStdTetrahedronPk()803 LagrangeStdTetrahedronPk::~LagrangeStdTetrahedronPk() {}
804
805 /*! init shape functions coefficients in basis monoms by solving linear systems n=nbdofs=(k+1)*(k+2)*(k+3)/6
806 */
initShapeFunCoeffs()807 void LagrangeStdTetrahedronPk::initShapeFunCoeffs()
808 {
809 //build matrix
810 Matrix<real_t> M(nbDofs_,nbDofs_);
811 std::vector<real_t>::iterator itM=M.begin();
812 std::vector<RefDof*>::iterator itd=refDofs.begin();
813 std::vector<real_t>::const_iterator itp;
814 for (number_t i=0; i<nbDofs_; i++, itM+=nbDofs_, itd++)
815 {
816 itp=(*itd)->coords();
817 monomes(*itp, *(itp+1),*(itp+2), itM, 0);
818 }
819 //solving system
820 shapeFunCoeffs=Matrix<real_t>(nbDofs_,nbDofs_);
821 for (number_t k=1; k<=nbDofs_; k++) shapeFunCoeffs(k,k)=1.;
822 number_t row;
823 real_t piv;
824 bool solved=gaussMultipleSolver(M, shapeFunCoeffs, nbDofs_, piv, row);
825 if (!solved)
826 {
827 where("LagrangeStdTetrahedronPk::initShapeFunCoeffs");
828 error("mat_noinvert");
829 }
830 }
831
832 /*! computeShapeValues defines Lagrange Pk Reference Element shape functions
833 general method using matrix approach :
834 the nth basis function is represented as An0 An1*y An2*y^2 + ...+Ank*y^k + ... + AnNx^k
835 where An. is the nth row of a matrix stored in RefElement class as shapeFunCoeffs
836 this matrix has to be computed before
837 the monomes function returns values or derivatives of monoms at a given point
838 */
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const839 void LagrangeStdTetrahedronPk::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
840 {
841 real_t x=*it_pt, y=*(it_pt+1), z=*(it_pt+2);
842 Vector<real_t> mons(nbDofs_);
843 monomes(x,y,z,mons.begin(),0);
844 shv.w = shapeFunCoeffs*mons;
845 if (withDeriv)
846 {
847 monomes(x,y,z,mons.begin(),1);
848 shv.dw[0] = shapeFunCoeffs*mons;
849 monomes(x,y,z,mons.begin(),2);
850 shv.dw[1] = shapeFunCoeffs*mons;
851 monomes(x,y,z,mons.begin(),3);
852 shv.dw[2] = shapeFunCoeffs*mons;
853 }
854 }
855
856 /*! compute value or derivatives of monomes (1, x, y, z, x2, y2, z2, xy, xz, yz ...) at a point (x,y,z)
857 monomes ordering is given by nested loops i=1,n ; j=1,n-i; k=1,n-i-j
858 der = 0 : values
859 der = 1 : dx values
860 der = 2 : dy values
861 der = 3 : dz values
862 */
monomes(real_t x,real_t y,real_t z,std::vector<real_t>::iterator itm,number_t der) const863 void LagrangeStdTetrahedronPk::monomes(real_t x, real_t y, real_t z, std::vector<real_t>::iterator itm, number_t der) const
864 {
865 number_t n=interpolation_p->numtype;
866 real_t vi=1.;
867 switch (der)
868 {
869 case 0:
870 {
871 for (number_t i=0; i<=n; i++, vi*=x)
872 {
873 real_t vj=1.;
874 for (number_t j=0; j<=n-i; j++, vj*=y)
875 {
876 real_t vk=1.;
877 for (number_t k=0; k<=n-i-j; k++,itm++, vk*=z) *itm = vi*vj*vk;
878 }
879 }
880 return;
881 }
882 case 1:
883 {
884 for (number_t j=0; j<(n+1)*(n+2)/2; j++,itm++) *itm = 0; //does not depends on x
885 for (number_t i=1; i<=n; i++, vi*=x)
886 {
887 real_t vj=1.;
888 for (number_t j=0; j<=n-i; j++,vj*=y)
889 {
890 real_t vk=1.;
891 for(number_t k=0; k<=n-i-j; k++,itm++, vk*=z) *itm = i*vi*vj*vk;
892 }
893 }
894 return;
895 }
896 case 2:
897 {
898 for (number_t i=0; i<=n; i++, vi*=x)
899 {
900 for (number_t j=0; j<=n-i; j++ ,itm++) *itm=0.; //does not depends on y
901 real_t vj=1.;
902 for (number_t j=1; j<=n-i; j++,vj*=y)
903 {
904 real_t vk=1.;
905 for(number_t k=0; k<=n-i-j; k++,itm++, vk*=z) *itm = j*vi*vj*vk;
906 }
907 }
908 return;
909 }
910 case 3 :
911 {
912 for (number_t i=0; i<=n; i++, vi*=x)
913 {
914 real_t vj=1.;
915 for (number_t j=0; j<=n-i; j++,vj*=y)
916 {
917 *itm=0.;
918 itm++;
919 real_t vk=1.;
920 for (number_t k=1; k<=n-i-j; k++,itm++, vk*=z) *itm = k*vi*vj*vk;
921 }
922 }
923 return;
924 }
925 default :
926 where("LagrangeStdTetrahedron::monomes");
927 error("bad_derivative",der);
928 }
929 }
930
931 // // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
932 // #ifdef FIG4TEX_OUTPUT
933 // //! fig4TeX builds a tex file to display local numbering and coordinates
934 // void LagrangeTetrahedron::fig4TeX() {
935 // // number_t of points on each edge
936 // number_t nbPtsPerSide(nbDofs(1, 1));
937 // // number_t of points on all edges
938 // number_t nbPtsOnSides = 3*(nbPtsPerSide-1);
939 // char b=char(92);
940
941 // std::ofstream os;
942 // os.open((fig4TexFileName()).c_str(), ios::out);
943 // os << "%" << b << "input fig4tex.tex"<<b<<"input TeXnicolor.tex"<<b<<"input Fonts.tex";
944 // os <<std::endl<<b<<"fourteenpoint"<<b<<"nopagenumbers";
945 // os <<std::endl<<b<<"newbox"<<b<<"mybox";
946 // if ( interpolation_p->numtype < 4 )
947 // {
948 // os <<std::endl<<b<<"figinit{200pt,realistic}";
949 // os <<b<<"def"<<b<<"Obdist{17}"<<b<<"def"<<b<<"Valpsi{40}"<<b<<"def"<<b<<"Valtheta{20}";
950 // }
951 // else
952 // {
953 // os <<std::endl<<b<<"figinit{300pt,realistic}";
954 // os <<b<<"def"<<b<<"Obdist{17}"<<b<<"def"<<b<<"Valpsi{30}"<<b<<"def"<<b<<"Valtheta{30}";
955 // }
956 // // 3D settings
957 // os <<std::endl<<b<<"figsetview("<<b<<"Valpsi,"<<b<<"Valtheta)"<<b<<"figsetobdist("<<b<<"Obdist)";
958 // string longName = "Lagrange Tetrahedron of degree ";
959
960 // // Draw element and display vertices
961 // number_t pts(0);
962 // // define vertices
963 // fig4TeXVertexPt(os, geom_ref_elt_p->vertices(), pts);
964
965 // // start postscript file
966 // os <<std::endl<<b<<"psbeginfig{}"
967 // <<std::endl<<b<<"pssetwidth{1.2}"<<b<<"psline[1,2,3,1]"
968 // <<std::endl<<b<<"pssetdash{7}"<<b<<"psline[1,4,2,4,3]";
969 // if ( nbPtsPerSide > 4 )
970 // {
971 // // Draw reduced side on each side
972 // for ( number_t fa = 1; fa <= 4 ; fa++ ) fig4TeXSmallFace(os, fa, 3, nbPtsOnSides);
973 // }
974 // real_t transl(0.);
975 // if ( nbInternalDofs_ > 1 )
976 // {
977 // os <<std::endl<<"% internal element";
978 // // Draw reduced tetrahedron inside
979 // fig4TeXInternalElement(os, transl);
980 // if ( interpolation_p->numtype > 5)
981 // {
982 // transl = -0.2*interpolation_p->numtype;
983 // fig4TeXInternalElement(os, transl);
984 // }
985 // }
986 // os <<std::endl<<b<<"psendfig"; // end of postcript output
987
988 // os <<std::endl<<b<<"figvisu{"<<b<<"mybox}{"<<b<<"bf "<< longName << interpolation_p->numtype <<"}{%"
989 // <<std::endl<<b<<"figsetmark{$"<<b<<"figBullet$}"<<b<<"figsetptname{{"<<b<<"bf#1}}";
990 // // vertex nodes
991 // os <<std::endl<<b<<"figwritesw1:(5pt)"<<std::endl<<b<<"figwritee2:(5pt)"
992 // <<std::endl<<b<<"figwriten3:(5pt)"<<std::endl<<b<<"figwritene4:(5pt)";
993
994 // if ( nbPtsPerSide > 2 )
995 // {
996 // // nodes on edges
997 // number_t ed(0);
998 // fig4TeXEdgePt(os , ed++, pts, "w");
999 // fig4TeXEdgePt(os , ed++, pts, "ne");
1000 // fig4TeXEdgePt(os , ed++, pts, "w");
1001 // fig4TeXEdgePt(os , ed++, pts, "ne");
1002 // fig4TeXEdgePt(os , ed++, pts, "nw");
1003 // fig4TeXEdgePt(os , ed++, pts, "s");
1004 // // nodes on sides
1005 // ed = 0;
1006 // fig4TeXFacePt(os , ed++, nbPtsOnSides, pts, "ne");
1007 // fig4TeXFacePt(os , ed++, nbPtsOnSides, pts, "nw");
1008 // fig4TeXFacePt(os , ed++, nbPtsOnSides, pts, "s");
1009 // fig4TeXFacePt(os , ed++, nbPtsOnSides, pts, "s");
1010 // }
1011 // // internal nodes
1012 // if ( nbInternalDofs_ > 0 ) fig4TeXInternalPoints(os, transl, pts);
1013 // // end of output
1014 // os <<std::endl<< "}" <<b<<"par"<<b<<"leftline{"<<b<<"box"<<b<<"mybox"<<b<<"hfill}" <<std::endl<<b<<"bye"<<std::endl;
1015 // os.close();
1016 // }
1017
1018 // void LagrangeTetrahedron::fig4TeXInternalElement(std::ofstream& os, const real_t transl)
1019 // {
1020 // // Draw reduced tetrahedron inside
1021 // const char b=char(92);
1022 // number_t pp(11);
1023 // for ( number_t i_dof = nbDofs_-nbInternalDofs_+1; i_dof < nbDofs_-nbInternalDofs_+5; i_dof++, pp++ )
1024 // {
1025 // std::vector<real_t>::const_iterator coor = refDofs[i_dof-1]->coords();
1026 // os <<std::endl<<b<<"figpt"<<pp<<":(";
1027 // os << setprecision(5) << transl + *coor++ <<"," << setprecision(5) << *coor++ <<",";
1028 // os << setprecision(5) << *coor++ <<")";
1029 // }
1030 // os<<std::endl<<b<<"psset (color="<<b<<"Redrgb)"
1031 // <<std::endl<<b<<"pssetdash{8}"<<b<<"pssetwidth{0.4}"<<b<<"psline[11,12,13,11]"
1032 // <<std::endl<<b<<"pssetdash{8}"<<b<<"psline[11,14,12,14,13]";
1033 // os<<std::endl<<b<<"psset (color="<<b<<"Blackrgb)";
1034 // }
1035
1036 // void LagrangeTetrahedron::fig4TeXInternalPoints(std::ofstream& os, const real_t transl, number_t& pts)
1037 // {
1038 // const char b=char(92);
1039 // os <<std::endl<<b<<"tenpoint"<<b<<"it"<<b<<"color{"<<b<<"Red}"<<std::endl<<b<<"figsetmark{$"<<b<<"diamond$}";
1040 // std::vector<RefDof*>::const_iterator refDofs_b(refDofs.begin()+pts), it_rd;
1041 // for ( it_rd = refDofs_b; it_rd != refDofs.end(); it_rd++, pts++ )
1042 // {
1043 // std::vector<real_t>::const_iterator coor = (*it_rd)->coords();
1044 // os <<std::endl<<b<<"figpt0:" << pts+1 <<"(";
1045 // os << setprecision(5) << transl + *coor++ <<","<< setprecision(5) << *coor++ <<",";
1046 // os << setprecision(5) << *coor++ <<")";
1047 // os <<b<<"figwrites0:"<<pts+1<<"(3pt)";
1048 // }
1049 // os <<std::endl<<b<<"endcolor{}";
1050 // }
1051
1052 // #endif /* FIG4TEX_OUTPUT */
1053 // // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
1054
1055 } // end of namespace xlifepp
1056