Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExpList3D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList3D.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Expansion list 3D definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iomanip>
37 #include <MultiRegions/ExpList3D.h>
38 
39 #include <LocalRegions/HexExp.h>
40 #include <LocalRegions/PrismExp.h>
41 #include <LocalRegions/PyrExp.h>
42 #include <LocalRegions/TetExp.h>
43 
46 
47 namespace Nektar
48 {
49  namespace MultiRegions
50  {
51 
53  {
54  SetExpType(e3D);
55  }
56 
58  {
59  SetExpType(e3D);
60  }
61 
63  {
64  }
65 
66 
68  const LibUtilities::BasisKey &TBa,
69  const LibUtilities::BasisKey &TBb,
70  const LibUtilities::BasisKey &TBc,
71  const LibUtilities::BasisKey &HBa,
72  const LibUtilities::BasisKey &HBb,
73  const LibUtilities::BasisKey &HBc,
75  const LibUtilities::PointsType TetNb):
76  ExpList(pSession,graph3D)
77  {
78  SetExpType(e3D);
79 
84 
85  const SpatialDomains::ExpansionMap &expansions = graph3D->GetExpansions();
86 
87  SpatialDomains::ExpansionMap::const_iterator expIt;
88  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
89  {
94 
95  if((TetGeom = boost::dynamic_pointer_cast<SpatialDomains::TetGeom>(expIt->second->m_geomShPtr)))
96  {
98  {
99 // Ntet = MemoryManager<LocalRegions::NodalTetExp>::AllocateSharedPtr(TetBa,TetBb,TetBc,TetNb,TetGeom);
100 // (*m_exp).push_back(Ntet);
101  }
102  else
103  {
105  (*m_exp).push_back(tet);
106  }
107 
109 
110  m_npoints += TBa.GetNumPoints()*TBb.GetNumPoints()*TBc.GetNumPoints();
111  }
112 /*
113  else if((PrismGeom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(expansions[i]->m_geomShPtr)))
114  {
115  prism = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(Ba,Bb,Bc,PrismGeom);
116  (*m_exp).push_back(prism);
117 
118  m_ncoeffs += StdRegions::StdPrismData::getNumberOfCoefficients(Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes());
119  m_npoints += Ba.GetNumPoints()*Bb.GetNumPoints()*Bc.GetNumPoints();
120 
121  }
122  else if((PyrGeom = boost::dynamic_pointer_cast<SpatialDomains::PyrGeom>(expansions[i]->m_geomShPtr)))
123  {
124  pyramid = MemoryManager<LocalRegions::PyrExp>::AllocateSharedPtr(Ba,Bb,Bc,PyrGeom);
125  (*m_exp).push_back(pyramid);
126 
127  m_ncoeffs += StdRegions::StdPyrData::getNumberOfCoefficients(Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes());
128  m_npoints += Ba.GetNumPoints()*Bb.GetNumPoints()*Bc.GetNumPoints();
129 
130  }
131 */
132  else if((HexGeom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(expIt->second->m_geomShPtr)))
133  {
134  hex = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(HBa,HBb,HBc, HexGeom);
135  (*m_exp).push_back(hex);
136 
137  m_ncoeffs += HBa.GetNumModes()*HBb.GetNumModes()*HBc.GetNumModes();
138  m_npoints += HBa.GetNumPoints()*HBb.GetNumPoints()*HBc.GetNumPoints();
139  }
140  else
141  {
142  ASSERTL0(false,"dynamic cast to a proper Geometry3D failed");
143  }
144 
145  }
146 
147  // Setup Default optimisation information.
148  int nel = GetExpSize();
151 
152  SetCoeffPhys();
153 
155  }
156 
157  /**
158  * Given a mesh \a graph3D, containing information about the domain and
159  * the spectral/hp element expansion, this constructor fills the list
160  * of local expansions \texttt{m_exp} with the proper expansions,
161  * calculates the total number of quadrature points
162  * \f$\boldsymbol{x}_i\f$ and the local expansion coefficients
163  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays #m_coeffs and
164  * #m_phys.
165  *
166  * @param graph3D A mesh, containing information about the domain
167  * and the spectral/hp element expansion.
168  */
170  const SpatialDomains::MeshGraphSharedPtr &graph3D,
171  const std::string &variable) :
172  ExpList(pSession,graph3D)
173  {
174  SetExpType(e3D);
175 
180 
181  const SpatialDomains::ExpansionMap &expansions
182  = graph3D->GetExpansions(variable);
183 
184  SpatialDomains::ExpansionMap::const_iterator expIt;
185  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
186  {
191 
192  if((TetGeom = boost::dynamic_pointer_cast<
193  SpatialDomains::TetGeom>(expIt->second->m_geomShPtr)))
194  {
196  = expIt->second->m_basisKeyVector[0];
198  = expIt->second->m_basisKeyVector[1];
200  = expIt->second->m_basisKeyVector[2];
201 
204  {
205  ASSERTL0(false,"LocalRegions::NodalTetExp is not "
206  "implemented yet");
207  }
208  else
209  {
211  ::AllocateSharedPtr(TetBa,TetBb,TetBc,
212  TetGeom);
213  (*m_exp).push_back(tet);
214  }
215  }
216  else if((PrismGeom = boost::dynamic_pointer_cast<SpatialDomains
217  ::PrismGeom>(expIt->second->m_geomShPtr)))
218  {
219  LibUtilities::BasisKey PrismBa
220  = expIt->second->m_basisKeyVector[0];
221  LibUtilities::BasisKey PrismBb
222  = expIt->second->m_basisKeyVector[1];
223  LibUtilities::BasisKey PrismBc
224  = expIt->second->m_basisKeyVector[2];
225 
227  ::AllocateSharedPtr(PrismBa,PrismBb,
228  PrismBc,PrismGeom);
229  (*m_exp).push_back(prism);
230  }
231  else if((PyrGeom = boost::dynamic_pointer_cast<
232  SpatialDomains::PyrGeom>(expIt->second->m_geomShPtr)))
233  {
235  = expIt->second->m_basisKeyVector[0];
237  = expIt->second->m_basisKeyVector[1];
239  = expIt->second->m_basisKeyVector[2];
240 
242  ::AllocateSharedPtr(PyrBa,PyrBb,PyrBc,
243  PyrGeom);
244  (*m_exp).push_back(pyramid);
245  }
246  else if((HexGeom = boost::dynamic_pointer_cast<
247  SpatialDomains::HexGeom>(expIt->second->m_geomShPtr)))
248  {
250  = expIt->second->m_basisKeyVector[0];
252  = expIt->second->m_basisKeyVector[1];
254  = expIt->second->m_basisKeyVector[2];
255 
257  ::AllocateSharedPtr(HexBa,HexBb,HexBc,
258  HexGeom);
259  (*m_exp).push_back(hex);
260  }
261  else
262  {
263  ASSERTL0(false,"dynamic cast to a proper Geometry3D "
264  "failed");
265  }
266 
267  }
268 
269  // Setup Default optimisation information.
270  int nel = GetExpSize();
273 
274  SetCoeffPhys();
276  }
277 
278  /**
279  * Given an expansion vector \a expansions, containing
280  * information about the domain and the spectral/hp element
281  * expansion, this constructor fills the list of local
282  * expansions \texttt{m_exp} with the proper expansions,
283  * calculates the total number of quadrature points
284  * \f$\boldsymbol{x}_i\f$ and the local expansion coefficients
285  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays
286  * #m_coeffs and #m_phys.
287  *
288  * @param expansions An expansion vector, containing
289  * information about the domain and the
290  * spectral/hp element expansion.
291  */
293  ExpList()
294  {
295  SetExpType(e3D);
296 
301 
302 
303  for(int i = 0; i < expansions.size(); ++i)
304  {
309 
310  SpatialDomains::ExpansionMap::const_iterator expmap = expansions.find(i);
311  ASSERTL1(expmap != expansions.end(), "Unable to find expansion.");
312  const SpatialDomains::ExpansionShPtr exp = expmap->second;
313 
314  if((TetGeom = boost::dynamic_pointer_cast<
315  SpatialDomains::TetGeom>(exp->m_geomShPtr)))
316  {
318  = exp->m_basisKeyVector[0];
320  = exp->m_basisKeyVector[1];
322  = exp->m_basisKeyVector[2];
323 
326  {
327  ASSERTL0(false,"LocalRegions::NodalTetExp is not "
328  "implemented yet");
329  }
330  else
331  {
333  ::AllocateSharedPtr(TetBa,TetBb,TetBc,
334  TetGeom);
335  (*m_exp).push_back(tet);
336  }
337  }
338  else if((PrismGeom = boost::dynamic_pointer_cast<
339  SpatialDomains::PrismGeom>(exp->m_geomShPtr)))
340  {
341  LibUtilities::BasisKey PrismBa
342  = exp->m_basisKeyVector[0];
343  LibUtilities::BasisKey PrismBb
344  = exp->m_basisKeyVector[1];
345  LibUtilities::BasisKey PrismBc
346  = exp->m_basisKeyVector[2];
347 
349  ::AllocateSharedPtr(PrismBa,PrismBb,
350  PrismBc,PrismGeom);
351  (*m_exp).push_back(prism);
352  }
353  else if((PyrGeom = boost::dynamic_pointer_cast<
354  SpatialDomains::PyrGeom>(exp->m_geomShPtr)))
355  {
357  = exp->m_basisKeyVector[0];
359  = exp->m_basisKeyVector[1];
361  = exp->m_basisKeyVector[2];
362 
364  ::AllocateSharedPtr(PyrBa,PyrBb,PyrBc,
365  PyrGeom);
366  (*m_exp).push_back(pyramid);
367  }
368  else if((HexGeom = boost::dynamic_pointer_cast<
369  SpatialDomains::HexGeom>(exp->m_geomShPtr)))
370  {
372  = exp->m_basisKeyVector[0];
374  = exp->m_basisKeyVector[1];
376  = exp->m_basisKeyVector[2];
377 
379  ::AllocateSharedPtr(HexBa,HexBb,HexBc,
380  HexGeom);
381  (*m_exp).push_back(hex);
382  }
383  else
384  {
385  ASSERTL0(false,"dynamic cast to a proper Geometry3D "
386  "failed");
387  }
388 
389  }
390 
391  // Setup Default optimisation information.
392  int nel = GetExpSize();
395 
396  SetCoeffPhys();
397  }
398 
399  /**
400  * Set up the storage for the concatenated list of
401  * coefficients and physical evaluations at the quadrature
402  * points. Each expansion (local element) is processed in turn
403  * to determine the number of coefficients and physical data
404  * points it contributes to the domain. Three arrays,
405  * #m_coeff_offset, #m_phys_offset and #m_offset_elmt_id, are
406  * also initialised and updated to store the data offsets of
407  * each element in the #m_coeffs and #m_phys arrays, and the
408  * element id that each consecutive block is associated
409  * respectively.
410  */
412  {
413  int i;
414 
415  // Set up offset information and array sizes
416  m_coeff_offset = Array<OneD,int>(m_exp->size());
417  m_phys_offset = Array<OneD,int>(m_exp->size());
418  m_offset_elmt_id = Array<OneD,int>(m_exp->size());
419 
420  m_ncoeffs = m_npoints = 0;
421 
422  for(i = 0; i < m_exp->size(); ++i)
423  {
425  m_phys_offset [i] = m_npoints;
426  m_offset_elmt_id[i] = i;
427  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
428  m_npoints += (*m_exp)[i]->GetTotPoints();
429  }
430 
431  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
432  m_phys = Array<OneD, NekDouble>(m_npoints);
433  }
434 
436  {
437  Array<OneD, int> NumShape(4,0);
438 
439  for(int i = 0; i < GetExpSize(); ++i)
440  {
441  switch ((*m_exp)[i]->DetShapeType())
442  {
443  case LibUtilities::eTetrahedron: NumShape[0]++; break;
444  case LibUtilities::ePyramid: NumShape[1]++; break;
445  case LibUtilities::ePrism: NumShape[2]++; break;
446  case LibUtilities::eHexahedron: NumShape[3]++; break;
447  default:
448  ASSERTL0(false, "Unknown expansion type.");
449  break;
450  }
451  }
452 
453  int three = 3;
455  ::AllocateSharedPtr(m_session,three,NumShape);
456  }
457 
458  void ExpList3D::v_WriteVtkPieceHeader(std::ofstream &outfile, int expansion)
459  {
460  int i,j,k;
461  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
462  int nquad1 = (*m_exp)[expansion]->GetNumPoints(1);
463  int nquad2 = (*m_exp)[expansion]->GetNumPoints(2);
464  int ntot = nquad0*nquad1*nquad2;
465  int ntotminus = (nquad0-1)*(nquad1-1)*(nquad2-1);
466 
467  Array<OneD,NekDouble> coords[3];
468  coords[0] = Array<OneD,NekDouble>(ntot);
469  coords[1] = Array<OneD,NekDouble>(ntot);
470  coords[2] = Array<OneD,NekDouble>(ntot);
471  (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
472 
473  outfile << " <Piece NumberOfPoints=\""
474  << ntot << "\" NumberOfCells=\""
475  << ntotminus << "\">" << endl;
476  outfile << " <Points>" << endl;
477  outfile << " <DataArray type=\"Float64\" "
478  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
479  outfile << " ";
480  for (i = 0; i < ntot; ++i)
481  {
482  for (j = 0; j < 3; ++j)
483  {
484  outfile << setprecision(8) << scientific
485  << (float)coords[j][i] << " ";
486  }
487  outfile << endl;
488  }
489 
490  outfile << endl;
491  outfile << " </DataArray>" << endl;
492  outfile << " </Points>" << endl;
493  outfile << " <Cells>" << endl;
494  outfile << " <DataArray type=\"Int32\" "
495  << "Name=\"connectivity\" format=\"ascii\">" << endl;
496  for (i = 0; i < nquad0-1; ++i)
497  {
498  for (j = 0; j < nquad1-1; ++j)
499  {
500  for (k = 0; k < nquad2-1; ++k)
501  {
502  outfile << k*nquad0*nquad1 + j*nquad0 + i << " "
503  << k*nquad0*nquad1 + j*nquad0 + i + 1 << " "
504  << k*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
505  << k*nquad0*nquad1 + (j+1)*nquad0 + i << " "
506  << (k+1)*nquad0*nquad1 + j*nquad0 + i << " "
507  << (k+1)*nquad0*nquad1 + j*nquad0 + i + 1 << " "
508  << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
509  << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i << " " << endl;
510  }
511  }
512  }
513  outfile << endl;
514  outfile << " </DataArray>" << endl;
515  outfile << " <DataArray type=\"Int32\" "
516  << "Name=\"offsets\" format=\"ascii\">" << endl;
517  for (i = 0; i < ntotminus; ++i)
518  {
519  outfile << i*8+8 << " ";
520  }
521  outfile << endl;
522  outfile << " </DataArray>" << endl;
523  outfile << " <DataArray type=\"UInt8\" "
524  << "Name=\"types\" format=\"ascii\">" << endl;
525  for (i = 0; i < ntotminus; ++i)
526  {
527  outfile << "12 ";
528  }
529  outfile << endl;
530  outfile << " </DataArray>" << endl;
531  outfile << " </Cells>" << endl;
532  outfile << " <PointData>" << endl;
533  }
534 
536  {
537  int i, j;
538  for (i = 0; i < m_exp->size(); ++i)
539  {
540  for (j = 0; j < (*m_exp)[i]->GetNfaces(); ++j)
541  {
542  (*m_exp)[i]->ComputeFaceNormal(j);
543  }
544  }
545  }
546 
548  const Array<OneD, NekDouble> &inarray,
549  Array<OneD, NekDouble> &outarray)
550  {
551  int cnt,cnt1;
552 
553  cnt = cnt1 = 0;
554  for(int i = 0; i < GetExpSize(); ++i)
555  {
556  // get new points key
557  int pt0 = (*m_exp)[i]->GetNumPoints(0);
558  int pt1 = (*m_exp)[i]->GetNumPoints(1);
559  int pt2 = (*m_exp)[i]->GetNumPoints(2);
560  int npt0 = (int) pt0*scale;
561  int npt1 = (int) pt1*scale;
562  int npt2 = (int) pt2*scale;
563 
564  LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0));
565  LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1));
566  LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2));
567 
568  // Interpolate points;
569  LibUtilities::Interp3D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
570  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
571  (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
572  &inarray[cnt], newPointsKey0,
573  newPointsKey1, newPointsKey2,
574  &outarray[cnt1]);
575 
576  cnt += pt0*pt1*pt2;
577  cnt1 += npt0*npt1*npt2;
578  }
579  }
580 
582  const Array<OneD, NekDouble> &inarray,
583  Array<OneD, NekDouble> &outarray)
584  {
585  int cnt,cnt1;
586 
587  cnt = cnt1 = 0;
588  for(int i = 0; i < GetExpSize(); ++i)
589  {
590  // get new points key
591  int pt0 = (*m_exp)[i]->GetNumPoints(0);
592  int pt1 = (*m_exp)[i]->GetNumPoints(1);
593  int pt2 = (*m_exp)[i]->GetNumPoints(2);
594  int npt0 = (int) pt0*scale;
595  int npt1 = (int) pt1*scale;
596  int npt2 = (int) pt2*scale;
597 
598  LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0));
599  LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1));
600  LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2));
601 
602  // Project points;
604  newPointsKey1,
605  newPointsKey2,
606  &inarray[cnt],
607  (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
608  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
609  (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
610  &outarray[cnt1]);
611 
612  cnt += npt0*npt1*npt2;
613  cnt1 += pt0*pt1*pt2;
614  }
615 
616  }
617  } //end of namespace
618 } //end of namespace
619