ExpList3D.cpp 23.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
///////////////////////////////////////////////////////////////////////////////
//
// File ExpList3D.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Expansion list 3D definition
//
///////////////////////////////////////////////////////////////////////////////
Mike Kirby's avatar
Mike Kirby committed
35

36
#include <iomanip>
Mike Kirby's avatar
Mike Kirby committed
37
#include <MultiRegions/ExpList3D.h>
38
39
40
41
42

#include <LocalRegions/HexExp.h>
#include <LocalRegions/PrismExp.h>
#include <LocalRegions/PyrExp.h>
#include <LocalRegions/TetExp.h>
Mike Kirby's avatar
Mike Kirby committed
43

44
45
namespace Nektar
{
46
    namespace MultiRegions
47
    {
Chris Cantwell's avatar
Chris Cantwell committed
48

49
        ExpList3D::ExpList3D(): ExpList()
50
51
        {
        }
Chris Cantwell's avatar
Chris Cantwell committed
52

53
54
55
        ExpList3D::ExpList3D(const ExpList3D &In): ExpList(In)
        {
        }
Chris Cantwell's avatar
Chris Cantwell committed
56

57
58
59
        ExpList3D::~ExpList3D()
        {
        }
Chris Cantwell's avatar
Chris Cantwell committed
60

61

Chris Cantwell's avatar
Chris Cantwell committed
62
        ExpList3D::ExpList3D(const LibUtilities::SessionReaderSharedPtr &pSession,
63
                             const LibUtilities::BasisKey &TBa,
64
65
66
67
68
                             const LibUtilities::BasisKey &TBb,
                             const LibUtilities::BasisKey &TBc,
                             const LibUtilities::BasisKey &HBa,
                             const LibUtilities::BasisKey &HBb,
                             const LibUtilities::BasisKey &HBc,
Chris Cantwell's avatar
Chris Cantwell committed
69
                             const SpatialDomains::MeshGraphSharedPtr &graph3D,
70
                             const LibUtilities::PointsType TetNb):
71
            ExpList(pSession,graph3D)
72
73
74
75
76
77
78
        {

            LocalRegions::TetExpSharedPtr tet;
            LocalRegions::HexExpSharedPtr hex;
            LocalRegions::PrismExpSharedPtr prism;
            LocalRegions::PyrExpSharedPtr pyramid;

79
            const SpatialDomains::ExpansionMap &expansions = graph3D->GetExpansions();
Chris Cantwell's avatar
Chris Cantwell committed
80

81
82
            SpatialDomains::ExpansionMap::const_iterator expIt;
            for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
83
84
85
86
87
            {
                SpatialDomains::TetGeomSharedPtr TetGeom;
                SpatialDomains::HexGeomSharedPtr HexGeom;
                SpatialDomains::PrismGeomSharedPtr PrismGeom;
                SpatialDomains::PyrGeomSharedPtr PyrGeom;
Chris Cantwell's avatar
Chris Cantwell committed
88

89
                if(TetGeom = boost::dynamic_pointer_cast<SpatialDomains::TetGeom>(expIt->second->m_geomShPtr)) // Tetrahedron
90
91
92
93
94
95
96
97
                {
                    if(TetNb < LibUtilities::SIZE_PointsType)
                    {
//                         Ntet = MemoryManager<LocalRegions::NodalTetExp>::AllocateSharedPtr(TetBa,TetBb,TetBc,TetNb,TetGeom);
//                         (*m_exp).push_back(Ntet);
                    }
                    else
                    {
98
                        tet = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(TBa,TBb,TBc,TetGeom);
99
100
101
                        (*m_exp).push_back(tet);
                    }

102
                    m_ncoeffs += StdRegions::StdTetData::getNumberOfCoefficients(TBa.GetNumModes(), TBb.GetNumModes(), TBc.GetNumModes());
Chris Cantwell's avatar
Chris Cantwell committed
103

104
                       m_npoints += TBa.GetNumPoints()*TBb.GetNumPoints()*TBc.GetNumPoints();
105
                }
106
/*                else if(PrismGeom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(expansions[i]->m_geomShPtr)) // Prism
107
108
109
110
                {
                      prism = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(Ba,Bb,Bc,PrismGeom);
                      (*m_exp).push_back(prism);

111
                      m_ncoeffs += StdRegions::StdPrismData::getNumberOfCoefficients(Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes());
112
                      m_npoints +=  Ba.GetNumPoints()*Bb.GetNumPoints()*Bc.GetNumPoints();
Chris Cantwell's avatar
Chris Cantwell committed
113

114
                }
115
                else if(PyrGeom = boost::dynamic_pointer_cast<SpatialDomains::PyrGeom>(expansions[i]->m_geomShPtr)) // Pyramid
116
117
118
119
                {
                     pyramid = MemoryManager<LocalRegions::PyrExp>::AllocateSharedPtr(Ba,Bb,Bc,PyrGeom);
                     (*m_exp).push_back(pyramid);

120
                     m_ncoeffs += StdRegions::StdPyrData::getNumberOfCoefficients(Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes());
121
122
123
                      m_npoints +=  Ba.GetNumPoints()*Bb.GetNumPoints()*Bc.GetNumPoints();

                }
124
*/                else if(HexGeom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(expIt->second->m_geomShPtr)) // Hexahedron
125
                {
126
                    hex = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(HBa,HBb,HBc, HexGeom);
127
                    (*m_exp).push_back(hex);
Chris Cantwell's avatar
Chris Cantwell committed
128

129
130
                    m_ncoeffs += HBa.GetNumModes()*HBb.GetNumModes()*HBc.GetNumModes();
                    m_npoints += HBa.GetNumPoints()*HBb.GetNumPoints()*HBc.GetNumPoints();
131
132
133
134
                }
                else
                {
                    ASSERTL0(false,"dynamic cast to a proper Geometry3D failed");
Chris Cantwell's avatar
Chris Cantwell committed
135
136
                }

137
            }
Chris Cantwell's avatar
Chris Cantwell committed
138

139
140
141
142
143
            // Setup Default optimisation information.
            int nel = GetExpSize();
            m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
                ::AllocateSharedPtr(nel);

144
            SetCoeffPhys();
145
146

            ReadGlobalOptimizationParameters();
147
        }
148

Chris Cantwell's avatar
Chris Cantwell committed
149
150
151
152
153
154
155
156
157
158
159
160
        /**
         * Given a mesh \a graph3D, containing information about the domain and
         * the spectral/hp element expansion, this constructor fills the list
         * of local expansions \texttt{m_exp} with the proper expansions,
         * calculates the total number of quadrature points
         * \f$\boldsymbol{x}_i\f$ and the local expansion coefficients
         * \f$\hat{u}^e_n\f$ and allocates memory for the arrays #m_coeffs and
         * #m_phys.
         *
         * @param   graph3D     A mesh, containing information about the domain
         *                      and the spectral/hp element expansion.
         */
Chris Cantwell's avatar
Chris Cantwell committed
161
162
163
        ExpList3D::ExpList3D(const LibUtilities::SessionReaderSharedPtr &pSession,
                const SpatialDomains::MeshGraphSharedPtr &graph3D) :
            ExpList(pSession,graph3D)
164
165
166
167
168
169
        {
            LocalRegions::TetExpSharedPtr tet;
            LocalRegions::HexExpSharedPtr hex;
            LocalRegions::PrismExpSharedPtr prism;
            LocalRegions::PyrExpSharedPtr pyramid;

170
            const SpatialDomains::ExpansionMap &expansions
171
                                        = graph3D->GetExpansions();
Chris Cantwell's avatar
Chris Cantwell committed
172

173
174
            SpatialDomains::ExpansionMap::const_iterator expIt;
            for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
175
176
177
178
179
180
181
            {
                SpatialDomains::TetGeomSharedPtr TetGeom;
                SpatialDomains::HexGeomSharedPtr HexGeom;
                SpatialDomains::PrismGeomSharedPtr PrismGeom;
                SpatialDomains::PyrGeomSharedPtr PyrGeom;

                if(TetGeom = boost::dynamic_pointer_cast<
182
                        SpatialDomains::TetGeom>(expIt->second->m_geomShPtr))
183
184
                {
                    LibUtilities::BasisKey TetBa
185
                                        = expIt->second->m_basisKeyVector[0];
186
                    LibUtilities::BasisKey TetBb
187
                                        = expIt->second->m_basisKeyVector[1];
188
                    LibUtilities::BasisKey TetBc
189
                                        = expIt->second->m_basisKeyVector[2];
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204

                    if(TetBa.GetBasisType() == LibUtilities::eGLL_Lagrange)
                    {
                      ASSERTL0(false,"LocalRegions::NodalTetExp is not "
                                     "implemented yet");
                    }
                    else
                    {
                        tet = MemoryManager<LocalRegions::TetExp>
                                        ::AllocateSharedPtr(TetBa,TetBb,TetBc,
                                                            TetGeom);
                        (*m_exp).push_back(tet);
                    }
                }
                else if(PrismGeom = boost::dynamic_pointer_cast<
205
                        SpatialDomains::PrismGeom>(expIt->second->m_geomShPtr))
206
207
                {
                    LibUtilities::BasisKey PrismBa
208
                                        = expIt->second->m_basisKeyVector[0];
209
                    LibUtilities::BasisKey PrismBb
210
                                        = expIt->second->m_basisKeyVector[1];
211
                    LibUtilities::BasisKey PrismBc
212
                                        = expIt->second->m_basisKeyVector[2];
213
214
215
216
217
218
219

                    prism = MemoryManager<LocalRegions::PrismExp>
                                        ::AllocateSharedPtr(PrismBa,PrismBb,
                                                            PrismBc,PrismGeom);
                    (*m_exp).push_back(prism);
                }
                else if(PyrGeom = boost::dynamic_pointer_cast<
220
                        SpatialDomains::PyrGeom>(expIt->second->m_geomShPtr))
221
222
                {
                    LibUtilities::BasisKey PyrBa
223
                                        = expIt->second->m_basisKeyVector[0];
224
                    LibUtilities::BasisKey PyrBb
225
                                        = expIt->second->m_basisKeyVector[1];
226
                    LibUtilities::BasisKey PyrBc
227
                                        = expIt->second->m_basisKeyVector[2];
228
229
230
231
232
233
234

                    pyramid = MemoryManager<LocalRegions::PyrExp>
                                        ::AllocateSharedPtr(PyrBa,PyrBb,PyrBc,
                                                            PyrGeom);
                    (*m_exp).push_back(pyramid);
                }
                else if(HexGeom = boost::dynamic_pointer_cast<
235
                        SpatialDomains::HexGeom>(expIt->second->m_geomShPtr))
236
237
                {
                    LibUtilities::BasisKey HexBa
238
                                        = expIt->second->m_basisKeyVector[0];
239
                    LibUtilities::BasisKey HexBb
240
                                        = expIt->second->m_basisKeyVector[1];
241
                    LibUtilities::BasisKey HexBc
242
                                        = expIt->second->m_basisKeyVector[2];
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262

                    hex = MemoryManager<LocalRegions::HexExp>
                                        ::AllocateSharedPtr(HexBa,HexBb,HexBc,
                                                            HexGeom);
                    (*m_exp).push_back(hex);
                }
                else
                {
                    ASSERTL0(false,"dynamic cast to a proper Geometry3D "
                                   "failed");
                }

            }

            // Setup Default optimisation information.
            int nel = GetExpSize();
            m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
                ::AllocateSharedPtr(nel);

            SetCoeffPhys();
263
            ReadGlobalOptimizationParameters();
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
        }

        /**
         * Given an expansion vector \a expansions, containing
         * information about the domain and the spectral/hp element
         * expansion, this constructor fills the list of local
         * expansions \texttt{m_exp} with the proper expansions,
         * calculates the total number of quadrature points
         * \f$\boldsymbol{x}_i\f$ and the local expansion coefficients
         * \f$\hat{u}^e_n\f$ and allocates memory for the arrays
         * #m_coeffs and #m_phys.
         *
         * @param expansions An expansion vector, containing
         *                   information about the domain and the
         *                   spectral/hp element expansion.
         */
Chris Cantwell's avatar
Chris Cantwell committed
280
        ExpList3D::ExpList3D(const SpatialDomains::ExpansionMap &expansions):
281
282
283
284
285
286
287
288
            ExpList()
        {
            LocalRegions::TetExpSharedPtr tet;
            LocalRegions::HexExpSharedPtr hex;
            LocalRegions::PrismExpSharedPtr prism;
            LocalRegions::PyrExpSharedPtr pyramid;


289
290
291
292
293
294
            for(int i = 0; i < expansions.size(); ++i)
            {
                SpatialDomains::TetGeomSharedPtr TetGeom;
                SpatialDomains::HexGeomSharedPtr HexGeom;
                SpatialDomains::PrismGeomSharedPtr PrismGeom;
                SpatialDomains::PyrGeomSharedPtr PyrGeom;
Chris Cantwell's avatar
Chris Cantwell committed
295

Chris Cantwell's avatar
Chris Cantwell committed
296
297
298
299
                SpatialDomains::ExpansionMap::const_iterator expmap = expansions.find(i);
                ASSERTL1(expmap != expansions.end(), "Unable to find expansion.");
                const SpatialDomains::ExpansionShPtr exp = expmap->second;

Chris Cantwell's avatar
Chris Cantwell committed
300
                if(TetGeom = boost::dynamic_pointer_cast<
Chris Cantwell's avatar
Chris Cantwell committed
301
                        SpatialDomains::TetGeom>(exp->m_geomShPtr))
302
                {
Chris Cantwell's avatar
Chris Cantwell committed
303
                    LibUtilities::BasisKey TetBa
Chris Cantwell's avatar
Chris Cantwell committed
304
                                        = exp->m_basisKeyVector[0];
Chris Cantwell's avatar
Chris Cantwell committed
305
                    LibUtilities::BasisKey TetBb
Chris Cantwell's avatar
Chris Cantwell committed
306
                                        = exp->m_basisKeyVector[1];
Chris Cantwell's avatar
Chris Cantwell committed
307
                    LibUtilities::BasisKey TetBc
Chris Cantwell's avatar
Chris Cantwell committed
308
                                        = exp->m_basisKeyVector[2];
Chris Cantwell's avatar
Chris Cantwell committed
309

310
                    if(TetBa.GetBasisType() == LibUtilities::eGLL_Lagrange)
311
                    {
Chris Cantwell's avatar
Chris Cantwell committed
312
313
                      ASSERTL0(false,"LocalRegions::NodalTetExp is not "
                                     "implemented yet");
314
315
316
                    }
                    else
                    {
Chris Cantwell's avatar
Chris Cantwell committed
317
318
319
                        tet = MemoryManager<LocalRegions::TetExp>
                                        ::AllocateSharedPtr(TetBa,TetBb,TetBc,
                                                            TetGeom);
320
321
322
                        (*m_exp).push_back(tet);
                    }
                }
Chris Cantwell's avatar
Chris Cantwell committed
323
                else if(PrismGeom = boost::dynamic_pointer_cast<
Chris Cantwell's avatar
Chris Cantwell committed
324
                        SpatialDomains::PrismGeom>(exp->m_geomShPtr))
325
                {
Chris Cantwell's avatar
Chris Cantwell committed
326
                    LibUtilities::BasisKey PrismBa
Chris Cantwell's avatar
Chris Cantwell committed
327
                                        = exp->m_basisKeyVector[0];
Chris Cantwell's avatar
Chris Cantwell committed
328
                    LibUtilities::BasisKey PrismBb
Chris Cantwell's avatar
Chris Cantwell committed
329
                                        = exp->m_basisKeyVector[1];
Chris Cantwell's avatar
Chris Cantwell committed
330
                    LibUtilities::BasisKey PrismBc
Chris Cantwell's avatar
Chris Cantwell committed
331
                                        = exp->m_basisKeyVector[2];
Chris Cantwell's avatar
Chris Cantwell committed
332
333
334
335

                    prism = MemoryManager<LocalRegions::PrismExp>
                                        ::AllocateSharedPtr(PrismBa,PrismBb,
                                                            PrismBc,PrismGeom);
336
337
                    (*m_exp).push_back(prism);
                }
Chris Cantwell's avatar
Chris Cantwell committed
338
                else if(PyrGeom = boost::dynamic_pointer_cast<
Chris Cantwell's avatar
Chris Cantwell committed
339
                        SpatialDomains::PyrGeom>(exp->m_geomShPtr))
340
                {
Chris Cantwell's avatar
Chris Cantwell committed
341
                    LibUtilities::BasisKey PyrBa
Chris Cantwell's avatar
Chris Cantwell committed
342
                                        = exp->m_basisKeyVector[0];
Chris Cantwell's avatar
Chris Cantwell committed
343
                    LibUtilities::BasisKey PyrBb
Chris Cantwell's avatar
Chris Cantwell committed
344
                                        = exp->m_basisKeyVector[1];
Chris Cantwell's avatar
Chris Cantwell committed
345
                    LibUtilities::BasisKey PyrBc
Chris Cantwell's avatar
Chris Cantwell committed
346
                                        = exp->m_basisKeyVector[2];
Chris Cantwell's avatar
Chris Cantwell committed
347
348
349
350

                    pyramid = MemoryManager<LocalRegions::PyrExp>
                                        ::AllocateSharedPtr(PyrBa,PyrBb,PyrBc,
                                                            PyrGeom);
351
352
                    (*m_exp).push_back(pyramid);
                }
Chris Cantwell's avatar
Chris Cantwell committed
353
                else if(HexGeom = boost::dynamic_pointer_cast<
Chris Cantwell's avatar
Chris Cantwell committed
354
                        SpatialDomains::HexGeom>(exp->m_geomShPtr))
355
                {
Chris Cantwell's avatar
Chris Cantwell committed
356
                    LibUtilities::BasisKey HexBa
Chris Cantwell's avatar
Chris Cantwell committed
357
                                        = exp->m_basisKeyVector[0];
Chris Cantwell's avatar
Chris Cantwell committed
358
                    LibUtilities::BasisKey HexBb
Chris Cantwell's avatar
Chris Cantwell committed
359
                                        = exp->m_basisKeyVector[1];
Chris Cantwell's avatar
Chris Cantwell committed
360
                    LibUtilities::BasisKey HexBc
Chris Cantwell's avatar
Chris Cantwell committed
361
                                        = exp->m_basisKeyVector[2];
Chris Cantwell's avatar
Chris Cantwell committed
362
363
364
365

                    hex = MemoryManager<LocalRegions::HexExp>
                                        ::AllocateSharedPtr(HexBa,HexBb,HexBc,
                                                            HexGeom);
366
367
368
369
                    (*m_exp).push_back(hex);
                }
                else
                {
Chris Cantwell's avatar
Chris Cantwell committed
370
371
372
373
374
                    ASSERTL0(false,"dynamic cast to a proper Geometry3D "
                                   "failed");
                }

            }
375

376
            // Setup Default optimisation information.
377
378
379
380
            int nel = GetExpSize();
            m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
                ::AllocateSharedPtr(nel);

381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
            SetCoeffPhys();
        }

        /**
         * Set up the storage for the concatenated list of
         * coefficients and physical evaluations at the quadrature
         * points. Each expansion (local element) is processed in turn
         * to determine the number of coefficients and physical data
         * points it contributes to the domain. Three arrays,
         * #m_coeff_offset, #m_phys_offset and #m_offset_elmt_id, are
         * also initialised and updated to store the data offsets of
         * each element in the #m_coeffs and #m_phys arrays, and the
         * element id that each consecutive block is associated
         * respectively.
         */
        void ExpList3D::SetCoeffPhys()
        {
            int i;

            // Set up offset information and array sizes
            m_coeff_offset   = Array<OneD,int>(m_exp->size());
            m_phys_offset    = Array<OneD,int>(m_exp->size());
            m_offset_elmt_id = Array<OneD,int>(m_exp->size());

            m_ncoeffs = m_npoints = 0;

            for(i = 0; i < m_exp->size(); ++i)
            {
                m_coeff_offset[i]   = m_ncoeffs;
                m_phys_offset [i]   = m_npoints;
                m_offset_elmt_id[i] = i;
                m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
                m_npoints += (*m_exp)[i]->GetTotPoints();
            }

416
417
418
419
            m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
            m_phys   = Array<OneD, NekDouble>(m_npoints);
        }

420
        void ExpList3D::v_ReadGlobalOptimizationParameters()
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
        {
            Array<OneD, int> NumShape(4,0);

            for(int i = 0; i < GetExpSize(); ++i)
            {
                switch ((*m_exp)[i]->DetExpansionType())
                {
                case StdRegions::eTetrahedron:  NumShape[0]++; break;
                case StdRegions::ePyramid:      NumShape[1]++; break;
                case StdRegions::ePrism:        NumShape[2]++; break;
                case StdRegions::eHexahedron:   NumShape[3]++; break;
                }
            }

            int three = 3;
            m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
437
                ::AllocateSharedPtr(m_session,three,NumShape);
438
        }
439

440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
        void ExpList3D::v_WriteVtkPieceHeader(std::ofstream &outfile, int expansion)
        {
            int i,j,k;
            int coordim  = (*m_exp)[expansion]->GetCoordim();
            int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
            int nquad1 = (*m_exp)[expansion]->GetNumPoints(1);
            int nquad2 = (*m_exp)[expansion]->GetNumPoints(2);
            int ntot = nquad0*nquad1*nquad2;
            int ntotminus = (nquad0-1)*(nquad1-1)*(nquad2-1);

            Array<OneD,NekDouble> coords[3];
            coords[0] = Array<OneD,NekDouble>(ntot);
            coords[1] = Array<OneD,NekDouble>(ntot);
            coords[2] = Array<OneD,NekDouble>(ntot);
            (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);

            outfile << "    <Piece NumberOfPoints=\""
                    << ntot << "\" NumberOfCells=\""
                    << ntotminus << "\">" << endl;
            outfile << "      <Points>" << endl;
            outfile << "        <DataArray type=\"Float32\" "
                    << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
            outfile << "          ";
            for (i = 0; i < ntot; ++i)
            {
                for (j = 0; j < 3; ++j)
                {
467
468
                    outfile << setprecision(8) << scientific 
                            << (float)coords[j][i] << " ";
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
                }
                outfile << endl;
            }
            outfile << endl;
            outfile << "        </DataArray>" << endl;
            outfile << "      </Points>" << endl;
            outfile << "      <Cells>" << endl;
            outfile << "        <DataArray type=\"Int32\" "
                    << "Name=\"connectivity\" format=\"ascii\">" << endl;
            for (i = 0; i < nquad0-1; ++i)
            {
                for (j = 0; j < nquad1-1; ++j)
                {
                    for (k = 0; k < nquad2-1; ++k)
                    {
                        outfile << k*nquad0*nquad1 + j*nquad0 + i << " "
                                << k*nquad0*nquad1 + j*nquad0 + i + 1 << " "
                                << k*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
                                << k*nquad0*nquad1 + (j+1)*nquad0 + i << " "
                                << (k+1)*nquad0*nquad1 + j*nquad0 + i << " "
                                << (k+1)*nquad0*nquad1 + j*nquad0 + i + 1 << " "
                                << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
                                << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i << " " << endl;
                    }
                }
            }
            outfile << endl;
            outfile << "        </DataArray>" << endl;
            outfile << "        <DataArray type=\"Int32\" "
                    << "Name=\"offsets\" format=\"ascii\">" << endl;
            for (i = 0; i < ntotminus; ++i)
            {
                outfile << i*8+8 << " ";
            }
            outfile << endl;
            outfile << "        </DataArray>" << endl;
            outfile << "        <DataArray type=\"UInt8\" "
                    << "Name=\"types\" format=\"ascii\">" << endl;
            for (i = 0; i < ntotminus; ++i)
            {
                outfile << "12 ";
            }
            outfile << endl;
            outfile << "        </DataArray>" << endl;
            outfile << "      </Cells>" << endl;
            outfile << "      <PointData>" << endl;
        }

517
518
519
520
521
522
523
524
525
526
527
528
        void ExpList3D::v_SetUpPhysNormals()
        {
            int i, j;
            for (i = 0; i < m_exp->size(); ++i)
            {
                for (j = 0; j < (*m_exp)[i]->GetNfaces(); ++j)
                {
                    (*m_exp)[i]->ComputeFaceNormal(j);
                }
            }
        }

529
  } //end of namespace
Mike Kirby's avatar
Mike Kirby committed
530
531
} //end of namespace