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 LagrangeQuadrangle.cpp
19 \authors D. Martin, N. Kielbasiewicz
20 \since 17 dec 2002
21 \date 6 aug 2012
22
23 \brief Implementation of xlifepp::LagrangeQuadrangle class members and related functions
24 */
25
26 #include "LagrangeQuadrangle.hpp"
27 #include "../Interpolation.hpp"
28 #include "utils.h"
29
30 #include <iostream>
31
32 namespace xlifepp
33 {
34
35 //! LagrangeQuadrangle constructor for Lagrange reference elements
LagrangeQuadrangle(const Interpolation * interp_p)36 LagrangeQuadrangle::LagrangeQuadrangle(const Interpolation* interp_p)
37 : RefQuadrangle(interp_p)
38 {
39 name_ += "_Lagrange";
40 trace_p->push("LagrangeQuadrangle::LagrangeQuadrangle (" + name_ + ")");
41
42 // build element interpolation data
43 interpolationData();
44 // local numbering of points on side of sides
45 sideOfSideNumbering();
46 // local numbering of points on sides
47 sideNumbering();
48 // Reference Element on edges
49 sideRefElement();
50 maxDegree = 2*interpolation_p->numtype;
51
52 trace_p->pop();
53 }
54
~LagrangeQuadrangle()55 LagrangeQuadrangle::~LagrangeQuadrangle() {}
56
57 //! interpolationData defines Reference Element interpolation data
interpolationData()58 void LagrangeQuadrangle::interpolationData()
59 {
60 number_t interpNum = interpolation_p->numtype;
61
62 switch(interpNum)
63 {
64 case 0:
65 nbDofs_ = nbInternalDofs_ = 1;
66 break;
67 // Serendipity case ??: // Lagrange quadrangle Q2 serendipity element
68 // Serendipity nbDofs_= 8
69 // Serendipity nbDofsOnVertices_ = 4;
70 // Serendipity break;
71 default: // Lagrange quadrangle of order > 0
72 nbDofs_ = (interpNum + 1) * (interpNum + 1);
73 nbInternalDofs_ = (interpNum - 1) * (interpNum - 1);
74 nbDofsOnVertices_ = 4;
75 nbDofsInSides_ = nbDofs_ - nbInternalDofs_ - nbDofsOnVertices_;
76 break;
77 }
78
79 // Creating Degrees Of Freedom of reference element
80 // in Lagrange standard (H1-conforming) elements all D.o.F on edges are sharable
81 refDofs.reserve(nbDofs_);
82 lagrangeRefDofs(2, nbDofsOnVertices_, nbInternalDofs_, 4, nbDofsInSides_);
83 // lagrangeRefDofs(1, nbDofsOnVertices_, nbDofs_ - nbInternalDofs_, 2);
84
85 // allocate storage for shape functions
86 shapeValues.resize(*this, 2);
87 // Lagrange Quadrangle Q2 serendipity element
88 // {
89 // nbDofs_=8;
90 // }
91 nbPts_ = nbDofs_;
92 }
93
94
95 //! create standard Lagrange numbering on sides
sideNumbering()96 void LagrangeQuadrangle::sideNumbering()
97 {
98 number_t interpNum = interpolation_p->numtype;
99 if(interpNum==0) //Q0 interpolation
100 {
101 for(number_t side = 0; side < geomRefElem_p->nbSides(); side++)
102 {
103 sideDofNumbers_[side].resize(1);
104 sideDofNumbers_[side][0]=1;
105 }
106 return;
107 }
108
109 if(interpNum > 0)
110 {
111 number_t intN(interpNum);
112 // if(interpNum == _P1BubbleP3) { intN = 1; }
113 number_t nbVertices = geomRefElem_p->nbVertices();
114 number_t nbSides = geomRefElem_p->nbSides();
115 number_t more;
116 number_t nbVerticesPerSide = geomRefElem_p->sideVertexNumbers()[0].size();
117 number_t nbDofsPerSide = intN + 1;
118
119 sideDofNumbers_.resize(nbSides);
120 for(number_t side = 0; side < nbSides; side++)
121 {
122 sideDofNumbers_[side].resize(nbDofsPerSide);
123 for(number_t j = 0; j < nbVerticesPerSide; j++)
124 {
125 sideDofNumbers_[side][j] = geomRefElem_p->sideVertexNumber(j + 1, side + 1);
126 }
127 if(intN > 1)
128 {
129 more = nbVertices + 1;
130 for(number_t k = nbVerticesPerSide; k < nbDofsPerSide; k++, more += nbSides)
131 {
132 sideDofNumbers_[side][k] = side + more;
133 }
134 }
135 }
136 }
137 }
138
139 //! create standard Lagrange numbering on side of sides
sideOfSideNumbering()140 void LagrangeQuadrangle::sideOfSideNumbering()
141 {
142 number_t interpNum = interpolation_p->numtype;
143
144 if(interpNum > 0)
145 {
146 number_t nbSideOfSides = geomRefElem_p->nbSideOfSides();
147 number_t nbVerticesPerSideOfSide = geomRefElem_p->sideOfSideVertexNumbers()[0].size();
148
149 sideOfSideDofNumbers_.resize(nbSideOfSides);
150 for(number_t i = 0; i < nbSideOfSides; i++)
151 {
152 sideOfSideDofNumbers_[i].resize(nbVerticesPerSideOfSide);
153 for(number_t j = 0; j < nbVerticesPerSideOfSide; j++)
154 {
155 sideOfSideDofNumbers_[i][j] = geomRefElem_p->sideOfSideVertexNumber(j + 1, i + 1);
156 }
157 }
158 }
159 }
160
161 //! pointCoordinates defines Lagrange Reference Element point coordinates
vertexCoordinates()162 std::vector<RefDof*>::iterator LagrangeQuadrangle::vertexCoordinates()
163 {
164 std::vector<RefDof*>::iterator it_rd = refDofs.begin();
165 (*it_rd++)->coords(1., 0.);
166 (*it_rd++)->coords(1., 1.);
167 (*it_rd++)->coords(0., 1.);
168 (*it_rd++)->coords(0., 0.);
169 return it_rd;
170 }
171
pointCoordinates()172 void LagrangeQuadrangle::pointCoordinates()
173 {
174 /*
175 Examples of local numbering of Lagrange Qk quadrangles, k=1,2,3,4,5
176
177 3------2 3-----6-----2 3--10---6---2 3--14--10---6---2 3--18--14--10---6---2
178 | | | | | | | | | |
179 | | | | 7 15 14 9 7 19 22 18 13 7 23 30 26 22 17
180 4------1 7 9 5 | | | | | |
181 interp=1 | | 11 16 13 5 11 23 25 21 9 11 27 35 34 29 13
182 | | | | | | | |
183 4-----8-----1 4---8--12---1 15 20 24 17 5 15 31 36 33 25 9
184 interpNum = 2 interpNum = 3 | | | |
185 4---8--12--16---1 19 24 28 32 21 5
186 interpNum = 4 | |
187 4---8--12--16--20---1
188 interpNum = 5
189 */
190
191 int interpNum = interpolation_p->numtype;
192 std::vector<RefDof*>::iterator it_rd = refDofs.begin();
193 switch(interpNum)
194 {
195 case 0:
196 (*it_rd)->coords(.5, .5);
197 break;
198 default: // Lagrange quadrangle of degree > 0
199 it_rd = vertexCoordinates();
200 if(interpNum > 1)
201 {
202 // correspondence between segment and quadrangle local numbering
203 std::vector<number_t> s2q(2 * nbPts_);
204 tensorNumberingQuadrangle(interpNum, s2q);
205
206 number_t pt = 2 * geomRefElem_p->nbVertices();
207 // find 1d reference element associated with any edge
208 RefElement* sideRef = sideRefElems_[0];
209 while(it_rd != refDofs.end())
210 {
211 real_t x = *(sideRef->refDofs[s2q[pt++]]->coords());
212 real_t y = *(sideRef->refDofs[s2q[pt++]]->coords());
213 (*it_rd++)->coords(x, y);
214 }
215 }
216 }
217 }
218
219 /*
220 Examples of local numbering of Lagrange Qk quadrangles, k=1,2,3,4,5
221
222 3------2 3-----6-----2 3--10---6---2 3--14--10---6---2 3--18--14--10---6---2
223 | | | | | | | | | |
224 | | | | 7 15 14 9 7 19 22 18 13 7 23 30 26 22 17
225 4------1 7 9 5 | | | | | |
226 interp=1 | | 11 16 13 5 11 23 25 21 9 11 27 35 34 29 13
227 | | | | | | | |
228 4-----8-----1 4---8--12---1 15 20 24 17 5 15 31 36 33 25 9
229 interpNum = 2 interpNum = 3 | | | |
230 4---8--12--16---1 19 24 28 32 21 5
231 interpNum = 4 | |
232 4---8--12--16--20---1
233 interpNum = 5
234 */
235 //! return nodes numbers of first order triangle elements when splitting current element
splitP1() const236 std::vector<std::vector<number_t> > LagrangeQuadrangle::splitP1() const
237 {
238 std::vector<std::vector<number_t> > splitnum;
239 std::vector<number_t> tuple(3,0);
240 switch(interpolation_p->numtype)
241 {
242 case 0 :
243 case 1 :
244 tuple[0]=1; tuple[1]=3; tuple[2]=4;
245 splitnum.push_back(tuple);
246 tuple[0]=3; tuple[1]=1; tuple[2]=2;
247 splitnum.push_back(tuple);
248 return splitnum;
249 case 2 :
250 tuple[0]=8; tuple[1]=7; tuple[2]=4;
251 splitnum.push_back(tuple);
252 tuple[0]=7; tuple[1]=8; tuple[2]=9;
253 splitnum.push_back(tuple);
254 tuple[0]=1; tuple[1]=9; tuple[2]=8;
255 splitnum.push_back(tuple);
256 tuple[0]=9; tuple[1]=1; tuple[2]=5;
257 splitnum.push_back(tuple);
258 tuple[0]=9; tuple[1]=3; tuple[2]=7;
259 splitnum.push_back(tuple);
260 tuple[0]=3; tuple[1]=9; tuple[2]=6;
261 splitnum.push_back(tuple);
262 tuple[0]=5; tuple[1]=6; tuple[2]=9;
263 splitnum.push_back(tuple);
264 tuple[0]=6; tuple[1]=5; tuple[2]=2;
265 splitnum.push_back(tuple);
266 return splitnum;
267 default :
268 return splitLagrange2DToP1();
269 }
270 return splitnum;
271 }
272
273 // creates the list of nodes numbers of first order quadrangle when splitting 2D Lagrange element
274 // Recursive algorithm based on the definition of the numbering rule over a quadrangle.
splitQuad(dimen_t k,dimen_t ori)275 splitvec_t LagrangeQuadrangle::splitQuad(dimen_t k, dimen_t ori)
276 {
277 splitvec_t splitnum;
278 switch(k)
279 {
280 case 1:
281 {
282 int quad[] = {ori+1, ori+2, ori+3, ori+4};
283 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
284 } break;
285 case 2:
286 {
287 int quad[] = {ori+1, ori+5, ori+9, ori+8};
288 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
289 }
290 {
291 int quad[] = {ori+2, ori+6, ori+9, ori+5};
292 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
293 }
294 {
295 int quad[] = {ori+3, ori+7, ori+9, ori+6};
296 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
297 }
298 {
299 int quad[] = {ori+4, ori+8, ori+9, ori+7};
300 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
301 }
302 break;
303 default:
304 // treat inside of the quadrangle first
305 splitnum = splitQuad(k-2, ori+4*k);
306 // then proceed with the boundary quadrangles
307 {
308 number_t n = ori+4*k, cn[] = {n+2, n+3, n+4, n+1}, s2 = n, s3 = s2+1;
309 for(number_t num=1; num<5; num++) // k-1 quadrangles along each side
310 {
311 number_t s1 = ori+num, s4 = s2; s2 = s1+4;
312 number_t quad[] = {s1, s2, s3, s4}; // first one
313 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
314 for(number_t j=0; j<k-3; j++) // k-3 others
315 {
316 s2 += 4, s3 += 4;
317 number_t quad[] = {s2-4, s2, s3, s3-4};
318 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
319 }
320 s2 += 4, s4 = s3, s3 = cn[num-1];
321 {
322 number_t quad[] = {s2-4, s2, s3, s4}; // last one
323 splitnum.push_back(make_pair(_quadrangle, std::vector<number_t>(quad, quad+4)));
324 }
325 }
326 } break;
327 }
328 return splitnum;
329 }
330
331 //! return nodes numbers of first order quadrangle elements of same shape when splitting current element
splitO1() const332 splitvec_t LagrangeQuadrangle::splitO1() const
333 {
334 return splitQuad(interpolation_p->numtype, 0);
335 }
336 /*
337 splitvec_t splitnum;
338 std::vector<number_t> tuple(4,0);
339 switch (interpolation_p->numtype)
340 {
341 case 1 :
342 tuple[0] = 1; tuple[1]=2; tuple[2]=3; tuple[3]=4;
343 splitnum.push_back(std::make_pair(_quadrangle,tuple));
344 return splitnum;
345 case 2 :
346 tuple[0]=8; tuple[1]=9; tuple[2]=7; tuple[3]=4;
347 splitnum.push_back(std::make_pair(_quadrangle,tuple));
348 tuple[0]=1; tuple[1]=5; tuple[2]=9; tuple[3]=8;
349 splitnum.push_back(std::make_pair(_quadrangle,tuple));
350 tuple[0]=9; tuple[1]=6; tuple[2]=3; tuple[3]=7;
351 splitnum.push_back(std::make_pair(_quadrangle,tuple));
352 tuple[0]=5; tuple[1]=2; tuple[2]=6; tuple[3]=9;
353 splitnum.push_back(std::make_pair(_quadrangle,tuple));
354 return splitnum;
355 case 3 :
356 tuple[0]=8; tuple[1]=16; tuple[2]=11; tuple[3]=4;
357 splitnum.push_back(std::make_pair(_quadrangle,tuple));
358 tuple[0]=12; tuple[1]=13; tuple[2]=16; tuple[3]=8;
359 splitnum.push_back(std::make_pair(_quadrangle,tuple));
360 tuple[0]=1; tuple[1]=5; tuple[2]=13; tuple[3]=12;
361 splitnum.push_back(std::make_pair(_quadrangle,tuple));
362 tuple[0]=16; tuple[1]=15; tuple[2]=7; tuple[3]=11;
363 splitnum.push_back(std::make_pair(_quadrangle,tuple));
364 tuple[0]=13; tuple[1]=14; tuple[2]=15; tuple[3]=16;
365 splitnum.push_back(std::make_pair(_quadrangle,tuple));
366 tuple[0]=5; tuple[1]=9; tuple[2]=14; tuple[3]=13;
367 splitnum.push_back(std::make_pair(_quadrangle,tuple));
368 tuple[0]=15; tuple[1]=10; tuple[2]=3; tuple[3]=7;
369 splitnum.push_back(std::make_pair(_quadrangle,tuple));
370 tuple[0]=14; tuple[1]=6; tuple[2]=10; tuple[3]=15;
371 splitnum.push_back(std::make_pair(_quadrangle,tuple));
372 tuple[0]=9; tuple[1]=2; tuple[2]=6; tuple[3]=14;
373 splitnum.push_back(std::make_pair(_quadrangle,tuple));
374 return splitnum;
375 case 4 :
376 tuple[0]=8; tuple[1]=20; tuple[2]=15; tuple[3]=4;
377 splitnum.push_back(std::make_pair(_quadrangle,tuple));
378 tuple[0]=12; tuple[1]=24; tuple[2]=20; tuple[3]=8;
379 splitnum.push_back(std::make_pair(_quadrangle,tuple));
380 tuple[0]=16; tuple[1]=17; tuple[2]=24; tuple[3]=12;
381 splitnum.push_back(std::make_pair(_quadrangle,tuple));
382 tuple[0]=1; tuple[1]=5; tuple[2]=17; tuple[3]=16;
383 splitnum.push_back(std::make_pair(_quadrangle,tuple));
384 tuple[0]=20; tuple[1]=23; tuple[2]=11; tuple[3]=15;
385 splitnum.push_back(std::make_pair(_quadrangle,tuple));
386 tuple[0]=24; tuple[1]=25; tuple[2]=23; tuple[3]=20;
387 splitnum.push_back(std::make_pair(_quadrangle,tuple));
388 tuple[0]=17; tuple[1]=21; tuple[2]=25; tuple[3]=24;
389 splitnum.push_back(std::make_pair(_quadrangle,tuple));
390 tuple[0]=5; tuple[1]=9; tuple[2]=21; tuple[3]=17;
391 splitnum.push_back(std::make_pair(_quadrangle,tuple));
392 tuple[0]=23; tuple[1]=19; tuple[2]=7; tuple[3]=11;
393 splitnum.push_back(std::make_pair(_quadrangle,tuple));
394 tuple[0]=25; tuple[1]=22; tuple[2]=19; tuple[3]=23;
395 splitnum.push_back(std::make_pair(_quadrangle,tuple));
396 tuple[0]=21; tuple[1]=18; tuple[2]=22; tuple[3]=25;
397 splitnum.push_back(std::make_pair(_quadrangle,tuple));
398 tuple[0]=9; tuple[1]=13; tuple[2]=18; tuple[3]=21;
399 splitnum.push_back(std::make_pair(_quadrangle,tuple));
400 tuple[0]=19; tuple[1]=14; tuple[2]=3; tuple[3]=7;
401 splitnum.push_back(std::make_pair(_quadrangle,tuple));
402 tuple[0]=22; tuple[1]=10; tuple[2]=14; tuple[3]=19;
403 splitnum.push_back(std::make_pair(_quadrangle,tuple));
404 tuple[0]=18; tuple[1]=6; tuple[2]=10; tuple[3]=22;
405 splitnum.push_back(std::make_pair(_quadrangle,tuple));
406 tuple[0]=13; tuple[1]=2; tuple[2]=6; tuple[3]=18;
407 splitnum.push_back(std::make_pair(_quadrangle,tuple));
408 return splitnum;
409 case 5 :
410 tuple[0]=8; tuple[1]=24; tuple[2]=19; tuple[3]=4;
411 splitnum.push_back(std::make_pair(_quadrangle,tuple));
412 tuple[0]=12; tuple[1]=28; tuple[2]=24; tuple[3]=8;
413 splitnum.push_back(std::make_pair(_quadrangle,tuple));
414 tuple[0]=16; tuple[1]=32; tuple[2]=28; tuple[3]=12;
415 splitnum.push_back(std::make_pair(_quadrangle,tuple));
416 tuple[0]=20; tuple[1]=21; tuple[2]=32; tuple[3]=16;
417 splitnum.push_back(std::make_pair(_quadrangle,tuple));
418 tuple[0]=1; tuple[1]=5; tuple[2]=21; tuple[3]=20;
419 splitnum.push_back(std::make_pair(_quadrangle,tuple));
420 tuple[0]=24; tuple[1]=31; tuple[2]=15; tuple[3]=19;
421 splitnum.push_back(std::make_pair(_quadrangle,tuple));
422 tuple[0]=28; tuple[1]=36; tuple[2]=31; tuple[3]=24;
423 splitnum.push_back(std::make_pair(_quadrangle,tuple));
424 tuple[0]=32; tuple[1]=33; tuple[2]=36; tuple[3]=28;
425 splitnum.push_back(std::make_pair(_quadrangle,tuple));
426 tuple[0]=21; tuple[1]=25; tuple[2]=33; tuple[3]=32;
427 splitnum.push_back(std::make_pair(_quadrangle,tuple));
428 tuple[0]=5; tuple[1]=9; tuple[2]=25; tuple[3]=21;
429 splitnum.push_back(std::make_pair(_quadrangle,tuple));
430 tuple[0]=31; tuple[1]=27; tuple[2]=11; tuple[3]=15;
431 splitnum.push_back(std::make_pair(_quadrangle,tuple));
432 tuple[0]=36; tuple[1]=35; tuple[2]=27; tuple[3]=31;
433 splitnum.push_back(std::make_pair(_quadrangle,tuple));
434 tuple[0]=33; tuple[1]=34; tuple[2]=35; tuple[3]=36;
435 splitnum.push_back(std::make_pair(_quadrangle,tuple));
436 tuple[0]=25; tuple[1]=29; tuple[2]=34; tuple[3]=33;
437 splitnum.push_back(std::make_pair(_quadrangle,tuple));
438 tuple[0]=9; tuple[1]=13; tuple[2]=29; tuple[3]=25;
439 splitnum.push_back(std::make_pair(_quadrangle,tuple));
440 tuple[0]=27; tuple[1]=23; tuple[2]=7; tuple[3]=11;
441 splitnum.push_back(std::make_pair(_quadrangle,tuple));
442 tuple[0]=35; tuple[1]=30; tuple[2]=23; tuple[3]=27;
443 splitnum.push_back(std::make_pair(_quadrangle,tuple));
444 tuple[0]=34; tuple[1]=26; tuple[2]=30; tuple[3]=35;
445 splitnum.push_back(std::make_pair(_quadrangle,tuple));
446 tuple[0]=29; tuple[1]=22; tuple[2]=26; tuple[3]=34;
447 splitnum.push_back(std::make_pair(_quadrangle,tuple));
448 tuple[0]=13; tuple[1]=17; tuple[2]=22; tuple[3]=29;
449 splitnum.push_back(std::make_pair(_quadrangle,tuple));
450 tuple[0]=23; tuple[1]=18; tuple[2]=3; tuple[3]=7;
451 splitnum.push_back(std::make_pair(_quadrangle,tuple));
452 tuple[0]=30; tuple[1]=14; tuple[2]=18; tuple[3]=23;
453 splitnum.push_back(std::make_pair(_quadrangle,tuple));
454 tuple[0]=26; tuple[1]=10; tuple[2]=14; tuple[3]=30;
455 splitnum.push_back(std::make_pair(_quadrangle,tuple));
456 tuple[0]=22; tuple[1]=6; tuple[2]=10; tuple[3]=26;
457 splitnum.push_back(std::make_pair(_quadrangle,tuple));
458 tuple[0]=17; tuple[1]=2; tuple[2]=6; tuple[3]=22;
459 splitnum.push_back(std::make_pair(_quadrangle,tuple));
460 return splitnum;
461 case 6 :
462 tuple[0]=8; tuple[1]=28; tuple[2]=23; tuple[3]=4;
463 splitnum.push_back(std::make_pair(_quadrangle,tuple));
464 tuple[0]=12; tuple[1]=32; tuple[2]=28; tuple[3]=8;
465 splitnum.push_back(std::make_pair(_quadrangle,tuple));
466 tuple[0]=16; tuple[1]=36; tuple[2]=32; tuple[3]=12;
467 splitnum.push_back(std::make_pair(_quadrangle,tuple));
468 tuple[0]=20; tuple[1]=40; tuple[2]=36; tuple[3]=16;
469 splitnum.push_back(std::make_pair(_quadrangle,tuple));
470 tuple[0]=24; tuple[1]=25; tuple[2]=40; tuple[3]=20;
471 splitnum.push_back(std::make_pair(_quadrangle,tuple));
472 tuple[0]=1; tuple[1]=5; tuple[2]=25; tuple[3]=24;
473 splitnum.push_back(std::make_pair(_quadrangle,tuple));
474 tuple[0]=28; tuple[1]=39; tuple[2]=19; tuple[3]=23;
475 splitnum.push_back(std::make_pair(_quadrangle,tuple));
476 tuple[0]=32; tuple[1]=44; tuple[2]=39; tuple[3]=28;
477 splitnum.push_back(std::make_pair(_quadrangle,tuple));
478 tuple[0]=36; tuple[1]=48; tuple[2]=44; tuple[3]=32;
479 splitnum.push_back(std::make_pair(_quadrangle,tuple));
480 tuple[0]=40; tuple[1]=41; tuple[2]=48; tuple[3]=36;
481 splitnum.push_back(std::make_pair(_quadrangle,tuple));
482 tuple[0]=25; tuple[1]=29; tuple[2]=41; tuple[3]=40;
483 splitnum.push_back(std::make_pair(_quadrangle,tuple));
484 tuple[0]=5; tuple[1]=9; tuple[2]=29; tuple[3]=25;
485 splitnum.push_back(std::make_pair(_quadrangle,tuple));
486 tuple[0]=39; tuple[1]=35; tuple[2]=15; tuple[3]=19;
487 splitnum.push_back(std::make_pair(_quadrangle,tuple));
488 tuple[0]=44; tuple[1]=47; tuple[2]=35; tuple[3]=39;
489 splitnum.push_back(std::make_pair(_quadrangle,tuple));
490 tuple[0]=48; tuple[1]=49; tuple[2]=47; tuple[3]=44;
491 splitnum.push_back(std::make_pair(_quadrangle,tuple));
492 tuple[0]=41; tuple[1]=45; tuple[2]=49; tuple[3]=48;
493 splitnum.push_back(std::make_pair(_quadrangle,tuple));
494 tuple[0]=29; tuple[1]=33; tuple[2]=45; tuple[3]=41;
495 splitnum.push_back(std::make_pair(_quadrangle,tuple));
496 tuple[0]=9; tuple[1]=13; tuple[2]=33; tuple[3]=29;
497 splitnum.push_back(std::make_pair(_quadrangle,tuple));
498 tuple[0]=35; tuple[1]=31; tuple[2]=11; tuple[3]=15;
499 splitnum.push_back(std::make_pair(_quadrangle,tuple));
500 tuple[0]=47; tuple[1]=43; tuple[2]=31; tuple[3]=35;
501 splitnum.push_back(std::make_pair(_quadrangle,tuple));
502 tuple[0]=49; tuple[1]=46; tuple[2]=43; tuple[3]=47;
503 splitnum.push_back(std::make_pair(_quadrangle,tuple));
504 tuple[0]=45; tuple[1]=42; tuple[2]=46; tuple[3]=49;
505 splitnum.push_back(std::make_pair(_quadrangle,tuple));
506 tuple[0]=33; tuple[1]=37; tuple[2]=42; tuple[3]=45;
507 splitnum.push_back(std::make_pair(_quadrangle,tuple));
508 tuple[0]=13; tuple[1]=17; tuple[2]=37; tuple[3]=33;
509 splitnum.push_back(std::make_pair(_quadrangle,tuple));
510 tuple[0]=31; tuple[1]=27; tuple[2]=7; tuple[3]=11;
511 splitnum.push_back(std::make_pair(_quadrangle,tuple));
512 tuple[0]=43; tuple[1]=38; tuple[2]=27; tuple[3]=31;
513 splitnum.push_back(std::make_pair(_quadrangle,tuple));
514 tuple[0]=46; tuple[1]=34; tuple[2]=38; tuple[3]=43;
515 splitnum.push_back(std::make_pair(_quadrangle,tuple));
516 tuple[0]=42; tuple[1]=30; tuple[2]=34; tuple[3]=46;
517 splitnum.push_back(std::make_pair(_quadrangle,tuple));
518 tuple[0]=37; tuple[1]=26; tuple[2]=30; tuple[3]=42;
519 splitnum.push_back(std::make_pair(_quadrangle,tuple));
520 tuple[0]=17; tuple[1]=21; tuple[2]=26; tuple[3]=37;
521 splitnum.push_back(std::make_pair(_quadrangle,tuple));
522 tuple[0]=27; tuple[1]=22; tuple[2]=3; tuple[3]=7;
523 splitnum.push_back(std::make_pair(_quadrangle,tuple));
524 tuple[0]=38; tuple[1]=18; tuple[2]=22; tuple[3]=27;
525 splitnum.push_back(std::make_pair(_quadrangle,tuple));
526 tuple[0]=34; tuple[1]=14; tuple[2]=18; tuple[3]=38;
527 splitnum.push_back(std::make_pair(_quadrangle,tuple));
528 tuple[0]=30; tuple[1]=10; tuple[2]=14; tuple[3]=34;
529 splitnum.push_back(std::make_pair(_quadrangle,tuple));
530 tuple[0]=26; tuple[1]=6; tuple[2]=10; tuple[3]=30;
531 splitnum.push_back(std::make_pair(_quadrangle,tuple));
532 tuple[0]=21; tuple[1]=2; tuple[2]=6; tuple[3]=26;
533 splitnum.push_back(std::make_pair(_quadrangle,tuple));
534 return splitnum;
535 default :
536 return splitLagrange2DToQ1();
537 }
538 return splitnum;
539 }*/
540
541 /*! for a side of current element, return the list of (first order element number, side number) of the first order elements (P1)
542 as a list of pair<P1 element number, P1 side number>
543
544 \verbatim
545 s2
546 3----2 Side 1 : Elt 2 (3 1 2), side 2 (1,2) --> (2,2)
547 |\ | Side 2 : Elt 2 (3 1 2), side 3 (2 3) --> (2,3)
548 Q1 s3| \ |s1 Side 3 : Elt 1 (1 3 4), side 2 (3 4) --> (1,2)
549 4----1 Side 4 : Elt 1 (1 3 4), side 3 (4 1) --> (1,3)
550 s4
551 s2
552 3-----6-----2 Side 1 -> (4,2) (8,2)
553 | \ 6 | \ 8|
554 | 5 \ | 7 \ | Side 2 -> (6,3) (8,3)
555 Q2 s3 7-----9-----5 s1
556 | \ 2| \ 4| Side 3 -> (1,2) (5,2)
557 | 1 \ | 3 \ |
558 4-----8-----1 Side 4 -> (1,3) (3,3)
559 s4
560 \endverbatim
561 */
splitP1Side(number_t side) const562 std::vector<std::pair<number_t,number_t> > LagrangeQuadrangle::splitP1Side(number_t side) const
563 {
564 if(side<1 || side>4) { error("split_bad_side_num",side); }
565 std::vector<std::pair<number_t,number_t> > sidenum;
566 switch(interpolation_p->numtype)
567 {
568 case 0 :
569 case 1 :
570 if(side==1) { sidenum.push_back(std::pair<number_t,number_t>(2,2)); }
571 else if(side==2) { sidenum.push_back(std::pair<number_t,number_t>(2,3)); }
572 else if(side==3) { sidenum.push_back(std::pair<number_t,number_t>(1,2)); }
573 else { sidenum.push_back(std::pair<number_t,number_t>(1,3)); }
574 return sidenum;
575 case 2 :
576 if(side==1)
577 {
578 sidenum.push_back(std::pair<number_t,number_t>(4,2));
579 sidenum.push_back(std::pair<number_t,number_t>(8,2));
580 }
581 else if(side==2)
582 {
583 sidenum.push_back(std::pair<number_t,number_t>(6,3));
584 sidenum.push_back(std::pair<number_t,number_t>(8,3));
585 }
586 else if(side==3)
587 {
588 sidenum.push_back(std::pair<number_t,number_t>(1,2));
589 sidenum.push_back(std::pair<number_t,number_t>(5,2));
590 }
591 else
592 {
593 sidenum.push_back(std::pair<number_t,number_t>(1,3));
594 sidenum.push_back(std::pair<number_t,number_t>(3,3));
595 }
596 return sidenum;
597 default : error("split_low_deg_Lagrange_elt", words("shape",_triangle), 2, words("shape",_triangle));
598 }
599 return sidenum;
600 }
601
602 //! 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) const603 void LagrangeQuadrangle::outputAsP1(std::ofstream& os, const int refNum, const std::vector<number_t>& ranks) const
604 {
605 int interpNum = interpolation_p->numtype;
606 /*
607 Numbering of d.o.f is made starting with the outtermost layer of elements
608 For instance first and second triangles are {13,12,1} and {5,13,1} for interpNum=3
609 (edges of which are part of the edges of the original quadrangle)
610 then recursively to the next inner layer
611 */
612 int pr(0), pr_1(-1), corner_1(1), corner, p1, p2, p3, p4;
613 while(interpNum > 2)
614 {
615 // first outtermost lower right corner: 13 for interpNum=3, 17 for interpNum=4, 21 for interpNum=5
616 corner_1 += 4 * interpNum;
617 corner = corner_1 - pr;
618 // part vertical layer parallel and close to { x_1 = 1 }
619 p1 = corner_1 - 1; p2 = 1; p3 = p2 + 4; p4 = corner;
620 for(int j = 1; j < interpNum; j++)
621 {
622 simplexVertexOutput(os, refNum, ranks[p1 + pr_1], ranks[p2 + pr_1], ranks[p4 + pr_1]);
623 simplexVertexOutput(os, refNum, ranks[p3 + pr_1], ranks[p4 + pr_1], ranks[p2 + pr_1]);
624 p1 = p4; p2 += 4; p3 = p2 + 4;
625 if(j == interpNum - 2) { corner += 1; p4 = corner; }
626 else { p4 += 4; }
627 }
628
629 // part vertical layer parallel and close to { x_2 = 1 }
630 p1 = corner; p2 = p1 - 5; p3 = 2; p4 = p3 + 4;
631 for(int j = 1; j < interpNum; j++)
632 {
633 simplexVertexOutput(os, refNum, ranks[p1 + pr_1], ranks[p2 + pr_1], ranks[p4 + pr_1]);
634 simplexVertexOutput(os, refNum, ranks[p3 + pr_1], ranks[p4 + pr_1], ranks[p2 + pr_1]);
635 p3 = p4; p2 = p1; p4 = p3 + 4;
636 if(j == interpNum - 2) { corner += 1; p1 = corner; }
637 else { p1 += 4; }
638 }
639
640 // part vertical layer parallel and close to { x_1 = 0 }
641 p2 = corner; p3 = p2 - 5; p4 = 3; p1 = p4 + 4;
642 for(int j = 1; j < interpNum; j++)
643 {
644 simplexVertexOutput(os, refNum, ranks[p1 + pr_1], ranks[p2 + pr_1], ranks[p4 + pr_1]);
645 simplexVertexOutput(os, refNum, ranks[p3 + pr_1], ranks[p4 + pr_1], ranks[p2 + pr_1]);
646 p3 = p2; p4 = p1; p1 = p4 + 4;
647 if(j == interpNum - 2) { corner += 1; p2 = corner; }
648 else { p2 += 4; }
649 }
650
651 // part vertical layer parallel and close to { x_2 = 0 }
652 p1 = 4; p2 = p1 + 4; p3 = corner; p4 = p3 - 5;
653 for(int j = 1; j < interpNum; j++)
654 {
655 simplexVertexOutput(os, refNum, ranks[p1 + pr_1], ranks[p2 + pr_1], ranks[p4 + pr_1]);
656 simplexVertexOutput(os, refNum, ranks[p3 + pr_1], ranks[p4 + pr_1], ranks[p2 + pr_1]);
657 p1 = p2; p4 = p3; p2 = p1 + 4;
658 if(j == interpNum - 2) { corner = corner_1 - pr; p3 = corner; }
659 else { p3 += 4; }
660 }
661 pr = corner_1 - 1; pr_1 = pr - 1; interpNum -= 2;
662 }
663
664 // finally consider the two cases for innermost quadrangle(s)
665 if(interpNum == 1)
666 {
667 /*
668 3-------2 3-------2
669 | | | \ #2 |
670 | | --> | \ |
671 | | | #1 \ |
672 4-------1 4-------1
673 */
674 simplexVertexOutput(os, refNum, ranks[4 + pr_1], ranks[1 + pr_1], ranks[3 + pr_1]);
675 simplexVertexOutput(os, refNum, ranks[2 + pr_1], ranks[3 + pr_1], ranks[1 + pr_1]);
676 }
677 else if(interpNum == 2)
678 {
679 /*
680 3-------6-------2 3-------6-------2
681 | | | \ #5 | #4 / |
682 | | | \ | / |
683 | | | #6 \ | / #3 |
684 7 9 5 --> 7-------9-------5
685 | | | #7 / | \ #2 |
686 | | | / | \ |
687 | | | / #8 | #1 \ |
688 4-------8-------1 4-------8-------1
689 */
690 simplexVertexOutput(os, refNum, ranks[8 + pr_1], ranks[1 + pr_1], ranks[9 + pr_1]);
691 simplexVertexOutput(os, refNum, ranks[1 + pr_1], ranks[5 + pr_1], ranks[9 + pr_1]);
692 simplexVertexOutput(os, refNum, ranks[5 + pr_1], ranks[2 + pr_1], ranks[9 + pr_1]);
693 simplexVertexOutput(os, refNum, ranks[2 + pr_1], ranks[6 + pr_1], ranks[9 + pr_1]);
694 simplexVertexOutput(os, refNum, ranks[6 + pr_1], ranks[3 + pr_1], ranks[9 + pr_1]);
695 simplexVertexOutput(os, refNum, ranks[3 + pr_1], ranks[7 + pr_1], ranks[9 + pr_1]);
696 simplexVertexOutput(os, refNum, ranks[7 + pr_1], ranks[4 + pr_1], ranks[9 + pr_1]);
697 simplexVertexOutput(os, refNum, ranks[4 + pr_1], ranks[8 + pr_1], ranks[9 + pr_1]);
698 }
699 else
700 {
701 noSuchFunction("outputAsP1");
702 }
703 }
704
705 /*!
706 standard Lagrange Q_k quadrangle Reference Element shape functions ((k+1)^2 nodes)
707 build using standard Lagrange P_k 1D Reference Element shape functions (k+1 nodes)
708 at standard (equidistant on edges) points or Gauss-Lobatto abscissae.
709 */
computeShapeValues(std::vector<real_t>::const_iterator it_pt,ShapeValues & shv,const bool withDeriv) const710 void LagrangeQuadrangle::computeShapeValues(std::vector<real_t>::const_iterator it_pt, ShapeValues& shv, const bool withDeriv) const
711 {
712 int interpNum = interpolation_p->numtype;
713 std::vector<real_t>::iterator it_w=shv.w.begin();
714 std::vector< std::vector<real_t> >::iterator it_dwdx1=shv.dw.begin(), it_dwdx2=it_dwdx1 + 1;
715 std::vector<real_t>::iterator it_dwdx1i=(*it_dwdx1).begin(), it_dwdx2i=(*it_dwdx2).begin();
716 real_t x1 = *it_pt;
717 real_t x2 = *(it_pt + 1);
718 real_t x1m1 = x1 - 1;
719 real_t x2m1 = x2 - 1;
720
721 switch(interpNum)
722 {
723 case _P0:
724 *it_w = 1.;
725 if(withDeriv) { *it_dwdx1i = 0.; *it_dwdx2i = 0.; }
726 break;
727 case _P1:
728 // In the (x1, x2) reference plane, shape functions of the Lagrange
729 // standard P1 quadrangle are given by
730 // w_1(x1,x2) = x1 (1-x2) | w_2(x1,x2) = x1 x2
731 // w_3(x1,x2) = (1-x1) x2 | w_4(x1,x2) = (1-x1)(1-x2)
732 *it_w++ = -x1 * x2m1; // 0
733 *it_w++ = x1 * x2; // 1
734 *it_w++ = -x1m1 * x2; // 2
735 *it_w = x1m1 * x2m1; // 3
736 if(withDeriv)
737 {
738 *it_dwdx1i++ = -x2m1; *it_dwdx2i++ = -x1; // 0
739 *it_dwdx1i++ = x2; *it_dwdx2i++ = x1; // 1
740 *it_dwdx1i++ = -x2; *it_dwdx2i++ = -x1m1; // 2
741 *it_dwdx1i = x2m1; *it_dwdx2i = x1m1; // 3
742 }
743 break;
744 case _P2:
745 {
746 real_t dx1m1=2*x1-1, dx2m1=2*x2-1, ax1=x1*dx1m1, ax2=x2*dx2m1,
747 bx1=x1m1*dx1m1, bx2=x2m1*dx2m1, cx1=-x1*x1m1, cx2=-x2*x2m1;
748 *it_w++ = ax1*bx2; *it_w++ = ax1*ax2; *it_w++ = bx1*ax2;
749 *it_w++ = bx1*bx2; *it_w++ = 4*ax1*cx2; *it_w++ = 4*cx1*ax2;
750 *it_w++ = 4*bx1*cx2; *it_w++ = 4*cx1*bx2; *it_w++ = 16*cx1*cx2;
751 if(withDeriv)
752 {
753 real_t dax1 = 4*x1-1, dax2 = 4*x2-1, dbx1 = 4*x1-3,
754 dbx2 = 4*x2-3, dcx1 = 1-2*x1, dcx2 = 1-2*x2;
755 *it_dwdx1i++ = dax1*bx2; *it_dwdx2i++ =ax1*dbx2;
756 *it_dwdx1i++ = dax1*ax2; *it_dwdx2i++ =ax1*dax2;
757 *it_dwdx1i++ = dbx1*ax2; *it_dwdx2i++ =bx1*dax2;
758 *it_dwdx1i++ = dbx1*bx2; *it_dwdx2i++ =bx1*dbx2;
759 *it_dwdx1i++ = 4*dax1*cx2; *it_dwdx2i++ =4*ax1*dcx2;
760 *it_dwdx1i++ = 4*dcx1*ax2; *it_dwdx2i++ =4*cx1*dax2;
761 *it_dwdx1i++ = 4*dbx1*cx2; *it_dwdx2i++ =4*bx1*dcx2;
762 *it_dwdx1i++ = 4*dcx1*bx2; *it_dwdx2i++ =4*cx1*dbx2;
763 *it_dwdx1i++ = 16*dcx1*cx2; *it_dwdx2i++ =16*cx1*dcx2;
764 }
765 }
766 break;
767 default:
768 // 1D shape functions
769 sideRefElems_[0]->computeShapeValues(it_pt++, withDeriv);
770 std::vector<real_t> shape_0 = sideRefElems_[0]->shapeValues.w;
771 std::vector<real_t> dshape_0 = sideRefElems_[0]->shapeValues.dw[0];
772 sideRefElems_[0]->computeShapeValues(it_pt, withDeriv);
773 std::vector<real_t> shape_1 = sideRefElems_[0]->shapeValues.w;
774 std::vector<real_t> dshape_1 = sideRefElems_[0]->shapeValues.dw[0];
775
776 // correspondence between segment and quadrangle local numbering
777 std::vector<number_t> s2q(2 * (interpNum + 1) * (interpNum + 1));
778 tensorNumberingQuadrangle(interpNum, s2q);
779
780 // tensor product (square) of 1D shape functions
781 std::vector<number_t>::iterator r_i(s2q.begin()), r_i_b;
782 for(it_w = shv.w.begin(); it_w != shv.w.end(); it_w++)
783 {
784 *it_w = shape_0[*r_i++] * shape_1[*r_i++];
785 }
786
787 if(withDeriv)
788 {
789 r_i = s2q.begin();
790 for(it_dwdx1i = (*it_dwdx1).begin(); it_dwdx1i != (*it_dwdx1).end(); it_dwdx1i++, it_dwdx2i++)
791 {
792 r_i_b = r_i++;
793 *it_dwdx1i = dshape_0[*r_i_b] * shape_1[*r_i];
794 *it_dwdx2i = shape_0[*r_i_b] * dshape_1[*r_i++];
795 }
796 }
797 break;
798 }
799 }
800
801 //! Lagrange standard Qk quadrangle Reference Element
LagrangeStdQuadrangle(const Interpolation * interp_p)802 LagrangeStdQuadrangle::LagrangeStdQuadrangle(const Interpolation* interp_p)
803 : LagrangeQuadrangle(interp_p)
804 {
805 number_t kd = interp_p->numtype;
806 name_ += "_" + tostring(kd);
807
808 // local coordinates of Points supporting D.o.F
809 pointCoordinates();
810 // build O1 splitting scheme
811 splitO1Scheme = splitO1();
812
813 #ifdef FIG4TEX_OUTPUT
814 // Draw a LagrangeStdQuadrangle using TeX and fig4TeX
815 if(kd > 0) { fig4TeX(); }
816 #endif
817 }
818
819 //! Lagrange Gauss-Lobatto Qk quadrangle Reference Element
LagrangeGLQuadrangle(const Interpolation * interp_p)820 LagrangeGLQuadrangle::LagrangeGLQuadrangle(const Interpolation* interp_p)
821 : LagrangeQuadrangle(interp_p)
822 {
823 number_t kd = interp_p->numtype;
824 name_ += "Gauss-Lobatto_" + tostring(kd);
825
826 // local coordinates of Points supporting D.o.F
827 pointCoordinates();
828 // build O1 splitting scheme
829 splitO1Scheme = splitO1();
830
831 #ifdef FIG4TEX_OUTPUT
832 // Draw a LagrangeGLQuadrangle using TeX and fig4TeX
833 if(kd > 0) { fig4TeX(); }
834 #endif
835 }
836
837
838 //! tensorNumberingQuadrangle : correspondence between segment and quadrangle local numbering
tensorNumberingQuadrangle(const int interpNum,std::vector<number_t> & s2q)839 void tensorNumberingQuadrangle(const int interpNum, std::vector<number_t>& s2q)
840 {
841 /*
842 Examples of local numbering of Lagrange Pk 1D elements
843
844 2----1 2---3---1 2---4---3---1 2---5---4---3---1 2---6---5---4---3---1
845 k=1 k=2 k=3 k=4 k=5
846
847 Examples of local numbering of Lagrange Qk quadrangles
848
849 3----2 3---6---2 3--10---6---2 3--14--10---6---2 3--18--14--10---6---2
850 | | | | | | | | | |
851 4----1 7 9 5 7 15 14 9 7 19 22 18 13 7 23 30 26 22 17
852 k=1 | | | | | | | |
853 ...... 4---8---1 11 16 13 5 11 23 25 21 9 11 27 35 34 29 13
854 k=2 | | | | | |
855 ................. 4---8--12---1 15 20 24 17 5 15 31 36 33 25 9
856 k=3 | | | |
857 ................................ 4---8--12--16---1 19 24 28 32 21 5
858 k=4 | |
859 ................................................... 4---8--12--16--20---1
860 k=5
861 */
862
863 std::vector<number_t>::iterator p = s2q.begin();
864 int nk = interpNum;
865 number_t q1 = 0, q2 = 1, q3 = 2, q4 = interpNum, q3_p, q4_m;
866
867 while(nk > 0)
868 {
869 // 'current perimeter' vertices
870 (*p++) = q1; (*p++) = q2;
871 (*p++) = q1; (*p++) = q1;
872 (*p++) = q2; (*p++) = q1;
873 (*p++) = q2; (*p++) = q2;
874
875 // non-vertex points on edges of current 'perimeter'
876 q3_p = q3; q4_m = q4;
877 for(dimen_t j = 2; j <= nk; j++)
878 {
879 (*p++) = q1; (*p++) = q4_m;
880 (*p++) = q3_p; (*p++) = q1;
881 (*p++) = q2; (*p++) = q3_p++;
882 (*p++) = q4_m--; (*p++) = q2;
883 }
884
885 q1 += 1; q2--; q3 += 1; q4--;
886 if(nk == interpNum) { q1++ ; q2 = interpNum; }
887 nk -= 2;
888 }
889
890 // central point if any
891 if(nk == 0) { (*p++) = q1; *p = q1; }
892
893 }
894
895 // // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
896 // #ifdef FIG4TEX_OUTPUT
897 // //! fig4TeX builds a tex file to display local numbering and coordinates
898 // void LagrangeQuadrangle::fig4TeX() {
899 // number_t nbVertices = geomRefElem_p->nbVertices();
900 // char b=char(92);
901
902 // std::ofstream os;
903 // os.open((fig4TexFileName()).c_str(), ios::out);
904 // os <<"%"<<b<< "input fig4tex.tex"<<b<<"input TeXnicolor.tex"<<b<<"input Fonts.tex";
905 // os <<std::endl<<b<<"fourteenpoint"<<b<<"nopagenumbers";
906 // os <<std::endl<<b<<"newbox"<<b<<"mybox";
907 // if ( nbDofs(1, 1) < 6 ) os <<std::endl<<b<<"figinit{150pt}";
908 // else os <<std::endl<<b<<"figinit{250pt}";
909 // string longName = "Lagrange Quadrangle of degree ";
910
911 // // Draw element and display vertices
912 // number_t pts(0);
913
914 // // define vertices
915 // fig4TeXVertexPt(os, geomRefElem_p->vertices(), pts);
916
917 // // start postscript file
918 // os <<std::endl<<b<<"psbeginfig{}"<<std::endl<<b<<"pssetwidth{1}"<<b<<"psline[";
919 // for ( size_t p = 1; p <= nbVertices; p++ ) os << p <<",";
920 // os <<"1]";
921 // os <<std::endl<<b<<"psendfig"; // end of postcript output
922
923 // os <<std::endl<<b<<"figvisu{"<<b<<"mybox}{"<< longName << interpolation_p->number() <<"}{%"
924 // <<std::endl<<b<<"figsetmark{$"<<b<<"bullet$}"<<b<<"figsetptname{{"<<b<<"bf#1}}";
925 // // vertex nodes
926 // os <<std::endl<<b<<"figwritese1:(5pt)"
927 // <<std::endl<<b<<"figwritene2:(5pt)"
928 // <<std::endl<<b<<"figwritenw3:(5pt)"
929 // <<std::endl<<b<<"figwritesw4:(5pt)";
930
931 // if ( nbDofs(1, 1) <= 7 ) {
932 // // edge vertices
933 // os <<std::endl<<b<<"figsetmark{}"
934 // <<std::endl<<b<<"color{"<<b<<"Red}"
935 // <<std::endl<<b<<"figwritegw1:${}_2$(1pt,7pt)"<<std::endl<<b<<"figwritegw1:${}_1$(7pt,1pt)"
936 // <<std::endl<<b<<"figwritegw2:${}_1$(1pt,-7pt)"<<std::endl<<b<<"figwritegw2:${}_2$(7pt,-1pt)"
937 // <<std::endl<<b<<"figwritege3:${}_2$(1pt,-7pt)"<<std::endl<<b<<"figwritege3:${}_1$(7pt,-1pt)"
938 // <<std::endl<<b<<"figwritege4:${}_1$(1pt,7pt)"<<std::endl<<b<<"figwritege4:${}_2$(7pt,1pt)";
939 // os <<std::endl<<b<<"endcolor";
940 // }
941
942 // // nodes on edge 1 (east), edge 2 (north), edge 3 (west), edge 4 (south)
943 // number_t ed(0);
944 // fig4TeXEdgePt(os , ed++, pts, "e", "w");
945 // fig4TeXEdgePt(os , ed++, pts, "n", "s");
946 // fig4TeXEdgePt(os , ed++, pts, "w", "e");
947 // fig4TeXEdgePt(os , ed++, pts, "s", "n");
948
949 // // internal nodes
950 // if ( nbInternalDofs_ > 0 ) fig4TeXInternalPoints2d(os, pts);
951
952 // // close figvisu and TeX file
953 // os <<std::endl<<"}" <<b<<"par"<<b<<"centerline{"<<b<<"box"<<b<<"mybox}";
954 // os <<std::endl<<b<<"bye";
955 // os.close();
956 // }
957 // #endif /* FIG4TEX_OUTPUT */
958 // // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST // FOR TEST //
959
960 } // end of namespace xlifepp
961