/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Contents of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3698 - (show annotations)
Tue Nov 29 00:47:29 2011 UTC (7 years, 4 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Brick.cpp
File size: 18928 byte(s)
[x] Rectangle...
[x] Brick node id's and weipa compatibility


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2011 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 #include <ripley/Brick.h>
15 extern "C" {
16 #include "paso/SystemMatrixPattern.h"
17 }
18
19 #if USE_SILO
20 #include <silo.h>
21 #ifdef ESYS_MPI
22 #include <pmpio.h>
23 #endif
24 #endif
25
26 #include <iomanip>
27
28 using namespace std;
29
30 namespace ripley {
31
32 Brick::Brick(int n0, int n1, int n2, double l0, double l1, double l2, int d0,
33 int d1, int d2) :
34 RipleyDomain(3),
35 m_gNE0(n0),
36 m_gNE1(n1),
37 m_gNE2(n2),
38 m_l0(l0),
39 m_l1(l1),
40 m_l2(l2),
41 m_NX(d0),
42 m_NY(d1),
43 m_NZ(d2)
44 {
45 // ensure number of subdivisions is valid and nodes can be distributed
46 // among number of ranks
47 if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
48 throw RipleyException("Invalid number of spatial subdivisions");
49
50 if (n0%m_NX > 0 || n1%m_NY > 0 || n2%m_NZ > 0)
51 throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
52
53 // local number of elements
54 m_NE0 = n0/m_NX;
55 m_NE1 = n1/m_NY;
56 m_NE2 = n2/m_NZ;
57 // local number of nodes (not necessarily owned)
58 m_N0 = m_NE0+1;
59 m_N1 = m_NE1+1;
60 m_N2 = m_NE2+1;
61 // bottom-left-front node is at (offset0,offset1,offset2) in global mesh
62 m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
63 m_offset1 = m_NE1*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);
64 m_offset2 = m_NE2*(m_mpiInfo->rank/(m_NX*m_NY));
65 populateSampleIds();
66 }
67
68
69 Brick::~Brick()
70 {
71 }
72
73 string Brick::getDescription() const
74 {
75 return "ripley::Brick";
76 }
77
78 bool Brick::operator==(const AbstractDomain& other) const
79 {
80 if (dynamic_cast<const Brick*>(&other))
81 return this==&other;
82
83 return false;
84 }
85
86 void Brick::dump(const string& fileName) const
87 {
88 #if USE_SILO
89 string fn(fileName);
90 if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
91 fn+=".silo";
92 }
93
94 const int NUM_SILO_FILES = 1;
95 const char* blockDirFmt = "/block%04d";
96 int driver=DB_HDF5;
97 string siloPath;
98 DBfile* dbfile = NULL;
99
100 #ifdef ESYS_MPI
101 PMPIO_baton_t* baton = NULL;
102 #endif
103
104 if (m_mpiInfo->size > 1) {
105 #ifdef ESYS_MPI
106 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
107 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
108 PMPIO_DefaultClose, (void*)&driver);
109 // try the fallback driver in case of error
110 if (!baton && driver != DB_PDB) {
111 driver = DB_PDB;
112 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
113 0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
114 PMPIO_DefaultClose, (void*)&driver);
115 }
116 if (baton) {
117 char str[64];
118 snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
119 siloPath = str;
120 dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
121 }
122 #endif
123 } else {
124 dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
125 getDescription().c_str(), driver);
126 // try the fallback driver in case of error
127 if (!dbfile && driver != DB_PDB) {
128 driver = DB_PDB;
129 dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
130 getDescription().c_str(), driver);
131 }
132 }
133
134 if (!dbfile)
135 throw RipleyException("dump: Could not create Silo file");
136
137 /*
138 if (driver==DB_HDF5) {
139 // gzip level 1 already provides good compression with minimal
140 // performance penalty. Some tests showed that gzip levels >3 performed
141 // rather badly on escript data both in terms of time and space
142 DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
143 }
144 */
145
146 boost::scoped_ptr<double> x(new double[m_N0]);
147 boost::scoped_ptr<double> y(new double[m_N1]);
148 boost::scoped_ptr<double> z(new double[m_N2]);
149 double* coords[3] = { x.get(), y.get(), z.get() };
150 pair<double,double> xdx = getFirstCoordAndSpacing(0);
151 pair<double,double> ydy = getFirstCoordAndSpacing(1);
152 pair<double,double> zdz = getFirstCoordAndSpacing(2);
153 #pragma omp parallel
154 {
155 #pragma omp for
156 for (dim_t i0 = 0; i0 < m_N0; i0++) {
157 coords[0][i0]=xdx.first+i0*xdx.second;
158 }
159 #pragma omp for
160 for (dim_t i1 = 0; i1 < m_N1; i1++) {
161 coords[1][i1]=ydy.first+i1*ydy.second;
162 }
163 #pragma omp for
164 for (dim_t i2 = 0; i2 < m_N2; i2++) {
165 coords[2][i2]=zdz.first+i2*zdz.second;
166 }
167 }
168 IndexVector dims = getNumNodesPerDim();
169 DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 3, DB_DOUBLE,
170 DB_COLLINEAR, NULL);
171
172 DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,
173 NULL, 0, DB_INT, DB_NODECENT, NULL);
174
175 // write element ids
176 dims = getNumElementsPerDim();
177 DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
178 &dims[0], 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
179
180 // rank 0 writes multimesh and multivar
181 if (m_mpiInfo->rank == 0) {
182 vector<string> tempstrings;
183 vector<char*> names;
184 for (dim_t i=0; i<m_mpiInfo->size; i++) {
185 stringstream path;
186 path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
187 tempstrings.push_back(path.str());
188 names.push_back((char*)tempstrings.back().c_str());
189 }
190 vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
191 DBSetDir(dbfile, "/");
192 DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
193 &types[0], NULL);
194 tempstrings.clear();
195 names.clear();
196 for (dim_t i=0; i<m_mpiInfo->size; i++) {
197 stringstream path;
198 path << "/block" << setw(4) << setfill('0') << right << i << "/nodeId";
199 tempstrings.push_back(path.str());
200 names.push_back((char*)tempstrings.back().c_str());
201 }
202 types.assign(m_mpiInfo->size, DB_QUADVAR);
203 DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
204 &types[0], NULL);
205 tempstrings.clear();
206 names.clear();
207 for (dim_t i=0; i<m_mpiInfo->size; i++) {
208 stringstream path;
209 path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
210 tempstrings.push_back(path.str());
211 names.push_back((char*)tempstrings.back().c_str());
212 }
213 DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
214 &types[0], NULL);
215 }
216
217 if (m_mpiInfo->size > 1) {
218 #ifdef ESYS_MPI
219 PMPIO_HandOffBaton(baton, dbfile);
220 PMPIO_Finish(baton);
221 #endif
222 } else {
223 DBClose(dbfile);
224 }
225
226 #else // USE_SILO
227 throw RipleyException("dump(): no Silo support");
228 #endif
229 }
230
231 const int* Brick::borrowSampleReferenceIDs(int fsType) const
232 {
233 switch (fsType) {
234 case Nodes:
235 return &m_nodeId[0];
236 case Elements:
237 return &m_elementId[0];
238 case FaceElements:
239 return &m_faceId[0];
240 default:
241 break;
242 }
243
244 stringstream msg;
245 msg << "borrowSampleReferenceIDs() not implemented for function space type "
246 << fsType;
247 throw RipleyException(msg.str());
248 }
249
250 bool Brick::ownSample(int fsCode, index_t id) const
251 {
252 #ifdef ESYS_MPI
253 if (fsCode == Nodes) {
254 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
255 const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
256 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
257 } else
258 throw RipleyException("ownSample() only implemented for Nodes");
259 #else
260 return true;
261 #endif
262 }
263
264 void Brick::interpolateOnDomain(escript::Data& target,
265 const escript::Data& in) const
266 {
267 const Brick& inDomain=dynamic_cast<const Brick&>(*(in.getFunctionSpace().getDomain()));
268 const Brick& targetDomain=dynamic_cast<const Brick&>(*(target.getFunctionSpace().getDomain()));
269 if (inDomain != *this)
270 throw RipleyException("Illegal domain of interpolant");
271 if (targetDomain != *this)
272 throw RipleyException("Illegal domain of interpolation target");
273
274 throw RipleyException("interpolateOnDomain() not implemented");
275 }
276
277 Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
278 bool reducedColOrder) const
279 {
280 if (reducedRowOrder || reducedColOrder)
281 throw RipleyException("getPattern() not implemented for reduced order");
282
283 throw RipleyException("getPattern() not implemented");
284 }
285
286 void Brick::Print_Mesh_Info(const bool full) const
287 {
288 RipleyDomain::Print_Mesh_Info(full);
289 if (full) {
290 cout << " Id Coordinates" << endl;
291 cout.precision(15);
292 cout.setf(ios::scientific, ios::floatfield);
293 pair<double,double> xdx = getFirstCoordAndSpacing(0);
294 pair<double,double> ydy = getFirstCoordAndSpacing(1);
295 pair<double,double> zdz = getFirstCoordAndSpacing(2);
296 for (index_t i=0; i < getNumNodes(); i++) {
297 cout << " " << setw(5) << m_nodeId[i]
298 << " " << xdx.first+(i%m_N0)*xdx.second
299 << " " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
300 << " " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
301 }
302 }
303 }
304
305 IndexVector Brick::getNumNodesPerDim() const
306 {
307 IndexVector ret;
308 ret.push_back(m_N0);
309 ret.push_back(m_N1);
310 ret.push_back(m_N2);
311 return ret;
312 }
313
314 IndexVector Brick::getNumElementsPerDim() const
315 {
316 IndexVector ret;
317 ret.push_back(m_NE0);
318 ret.push_back(m_NE1);
319 ret.push_back(m_NE2);
320 return ret;
321 }
322
323 IndexVector Brick::getNumFacesPerBoundary() const
324 {
325 IndexVector ret(6, 0);
326 //left
327 if (m_offset0==0)
328 ret[0]=m_NE1*m_NE2;
329 //right
330 if (m_mpiInfo->rank%m_NX==m_NX-1)
331 ret[1]=m_NE1*m_NE2;
332 //bottom
333 if (m_offset1==0)
334 ret[2]=m_NE0*m_NE2;
335 //top
336 if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
337 ret[3]=m_NE0*m_NE2;
338 //front
339 if (m_offset2==0)
340 ret[4]=m_NE0*m_NE1;
341 //back
342 if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
343 ret[5]=m_NE0*m_NE1;
344 return ret;
345 }
346
347 pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
348 {
349 if (dim==0)
350 return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
351 else if (dim==1)
352 return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
353 else if (dim==2)
354 return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
355
356 throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
357 }
358
359
360 //protected
361 dim_t Brick::getNumFaceElements() const
362 {
363 dim_t n=0;
364 //left
365 if (m_offset0==0)
366 n+=m_NE1*m_NE2;
367 //right
368 if (m_mpiInfo->rank%m_NX==m_NX-1)
369 n+=m_NE1*m_NE2;
370 //bottom
371 if (m_offset1==0)
372 n+=m_NE0*m_NE2;
373 //top
374 if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
375 n+=m_NE0*m_NE2;
376 //front
377 if (m_offset2==0)
378 n+=m_NE0*m_NE1;
379 //back
380 if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
381 n+=m_NE0*m_NE1;
382
383 return n;
384 }
385
386 //protected
387 void Brick::assembleCoordinates(escript::Data& arg) const
388 {
389 escriptDataC x = arg.getDataC();
390 int numDim = m_numDim;
391 if (!isDataPointShapeEqual(&x, 1, &numDim))
392 throw RipleyException("setToX: Invalid Data object shape");
393 if (!numSamplesEqual(&x, 1, getNumNodes()))
394 throw RipleyException("setToX: Illegal number of samples in Data object");
395
396 pair<double,double> xdx = getFirstCoordAndSpacing(0);
397 pair<double,double> ydy = getFirstCoordAndSpacing(1);
398 pair<double,double> zdz = getFirstCoordAndSpacing(2);
399 arg.requireWrite();
400 #pragma omp parallel for
401 for (dim_t i2 = 0; i2 < m_N2; i2++) {
402 for (dim_t i1 = 0; i1 < m_N1; i1++) {
403 for (dim_t i0 = 0; i0 < m_N0; i0++) {
404 double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
405 point[0] = xdx.first+i0*xdx.second;
406 point[1] = ydy.first+i1*ydy.second;
407 point[2] = zdz.first+i2*zdz.second;
408 }
409 }
410 }
411 }
412
413 //private
414 void Brick::populateSampleIds()
415 {
416 // identifiers are ordered from left to right, bottom to top, front to back
417 // on each rank, except for the shared nodes which are owned by the rank
418 // below / to the left / to the front of the current rank
419
420 // build node distribution vector first.
421 // m_nodeDistribution[i] is the first node id on rank i, that is
422 // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
423 m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
424 m_nodeDistribution[1]=getNumNodes();
425 for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
426 const index_t x = k%m_NX;
427 const index_t y = k%(m_NX*m_NY)/m_NX;
428 const index_t z = k/(m_NX*m_NY);
429 index_t numNodes=getNumNodes();
430 if (x>0)
431 numNodes-=m_N1*m_N2;
432 if (y>0)
433 numNodes-=m_N0*m_N2;
434 if (z>0)
435 numNodes-=m_N0*m_N1;
436 // if an edge was subtracted twice add it back
437 if (x>0 && y>0)
438 numNodes+=m_N2;
439 if (x>0 && z>0)
440 numNodes+=m_N1;
441 if (y>0 && z>0)
442 numNodes+=m_N0;
443 // the corner node was removed 3x and added back 3x, so subtract it
444 if (x>0 && y>0 && z>0)
445 numNodes--;
446 m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
447 }
448 m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
449
450 m_nodeId.resize(getNumNodes());
451
452 // the bottom, left and front planes are not owned by this rank so the
453 // identifiers need to be computed accordingly
454 const index_t left = (m_offset0==0 ? 0 : 1);
455 const index_t bottom = (m_offset1==0 ? 0 : 1);
456 const index_t front = (m_offset2==0 ? 0 : 1);
457
458 // case 1: all nodes on left plane are owned by rank on the left
459 if (left>0) {
460 const int neighbour=m_mpiInfo->rank-1;
461 const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
462 const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
463 #pragma omp parallel for
464 for (dim_t i2=front; i2<m_N2; i2++) {
465 for (dim_t i1=bottom; i1<m_N1; i1++) {
466 m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
467 + (i1-bottom+1)*leftN0
468 + (i2-front)*leftN0*leftN1 - 1;
469 }
470 }
471 }
472 // case 2: all nodes on bottom plane are owned by rank below
473 if (bottom>0) {
474 const int neighbour=m_mpiInfo->rank-m_NX;
475 const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
476 const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
477 #pragma omp parallel for
478 for (dim_t i2=front; i2<m_N2; i2++) {
479 for (dim_t i0=left; i0<m_N0; i0++) {
480 m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
481 + bottomN0*(bottomN1-1)
482 + (i2-front)*bottomN0*bottomN1 + i0-left;
483 }
484 }
485 }
486 // case 3: all nodes on front plane are owned by rank in front
487 if (front>0) {
488 const int neighbour=m_mpiInfo->rank-m_NX*m_NY;
489 const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
490 const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
491 const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
492 #pragma omp parallel for
493 for (dim_t i1=bottom; i1<m_N1; i1++) {
494 for (dim_t i0=left; i0<m_N0; i0++) {
495 m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]
496 + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;
497 }
498 }
499 }
500 // case 4: nodes on front bottom edge are owned by the corresponding rank
501 if (front>0 && bottom>0) {
502 const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);
503 const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
504 const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
505 const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
506 #pragma omp parallel for
507 for (dim_t i0=left; i0<m_N0; i0++) {
508 m_nodeId[i0]=m_nodeDistribution[neighbour]
509 + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;
510 }
511 }
512 // case 5: nodes on left bottom edge are owned by the corresponding rank
513 if (left>0 && bottom>0) {
514 const int neighbour=m_mpiInfo->rank-m_NX-1;
515 const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
516 const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
517 #pragma omp parallel for
518 for (dim_t i2=front; i2<m_N2; i2++) {
519 m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
520 + (1+i2-front)*N0*N1-1;
521 }
522 }
523 // case 6: nodes on left front edge are owned by the corresponding rank
524 if (left>0 && front>0) {
525 const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;
526 const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
527 const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
528 const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
529 #pragma omp parallel for
530 for (dim_t i1=bottom; i1<m_N1; i1++) {
531 m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
532 + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;
533 }
534 }
535 // case 7: bottom-left-front corner node owned by corresponding rank
536 if (left>0 && bottom>0 && front>0) {
537 const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;
538 const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
539 const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
540 const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);
541 m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;
542 }
543
544 // the rest of the id's are contiguous
545 const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
546 #pragma omp parallel for
547 for (dim_t i2=front; i2<m_N2; i2++) {
548 for (dim_t i1=bottom; i1<m_N1; i1++) {
549 for (dim_t i0=left; i0<m_N0; i0++) {
550 m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left
551 +(i1-bottom)*(m_N0-left)
552 +(i2-front)*(m_N0-left)*(m_N1-bottom);
553 }
554 }
555 }
556
557 // elements
558 m_elementId.resize(getNumElements());
559 #pragma omp parallel for
560 for (dim_t k=0; k<getNumElements(); k++) {
561 m_elementId[k]=k;
562 }
563
564 // face elements
565 m_faceId.resize(getNumFaceElements());
566 #pragma omp parallel for
567 for (dim_t k=0; k<getNumFaceElements(); k++) {
568 m_faceId[k]=k;
569 }
570 }
571
572 } // end of namespace ripley
573

  ViewVC Help
Powered by ViewVC 1.1.26