49 namespace MultiRegions
87 SpatialDomains::ExpansionMap::const_iterator expIt;
88 for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
95 if((TetGeom = boost::dynamic_pointer_cast<SpatialDomains::TetGeom>(expIt->second->m_geomShPtr)))
105 (*m_exp).push_back(tet);
132 else if((HexGeom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(expIt->second->m_geomShPtr)))
135 (*m_exp).push_back(hex);
142 ASSERTL0(
false,
"dynamic cast to a proper Geometry3D failed");
171 const std::string &variable) :
182 = graph3D->GetExpansions(variable);
184 SpatialDomains::ExpansionMap::const_iterator expIt;
185 for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
192 if((TetGeom = boost::dynamic_pointer_cast<
196 = expIt->second->m_basisKeyVector[0];
198 = expIt->second->m_basisKeyVector[1];
200 = expIt->second->m_basisKeyVector[2];
205 ASSERTL0(
false,
"LocalRegions::NodalTetExp is not "
213 (*m_exp).push_back(tet);
216 else if((PrismGeom = boost::dynamic_pointer_cast<SpatialDomains
217 ::PrismGeom>(expIt->second->m_geomShPtr)))
220 = expIt->second->m_basisKeyVector[0];
222 = expIt->second->m_basisKeyVector[1];
224 = expIt->second->m_basisKeyVector[2];
229 (*m_exp).push_back(prism);
231 else if((PyrGeom = boost::dynamic_pointer_cast<
235 = expIt->second->m_basisKeyVector[0];
237 = expIt->second->m_basisKeyVector[1];
239 = expIt->second->m_basisKeyVector[2];
244 (*m_exp).push_back(pyramid);
246 else if((HexGeom = boost::dynamic_pointer_cast<
250 = expIt->second->m_basisKeyVector[0];
252 = expIt->second->m_basisKeyVector[1];
254 = expIt->second->m_basisKeyVector[2];
259 (*m_exp).push_back(hex);
263 ASSERTL0(
false,
"dynamic cast to a proper Geometry3D "
303 for(
int i = 0; i < expansions.size(); ++i)
310 SpatialDomains::ExpansionMap::const_iterator expmap = expansions.find(i);
311 ASSERTL1(expmap != expansions.end(),
"Unable to find expansion.");
314 if((TetGeom = boost::dynamic_pointer_cast<
318 = exp->m_basisKeyVector[0];
320 = exp->m_basisKeyVector[1];
322 = exp->m_basisKeyVector[2];
327 ASSERTL0(
false,
"LocalRegions::NodalTetExp is not "
335 (*m_exp).push_back(tet);
338 else if((PrismGeom = boost::dynamic_pointer_cast<
342 = exp->m_basisKeyVector[0];
344 = exp->m_basisKeyVector[1];
346 = exp->m_basisKeyVector[2];
351 (*m_exp).push_back(prism);
353 else if((PyrGeom = boost::dynamic_pointer_cast<
357 = exp->m_basisKeyVector[0];
359 = exp->m_basisKeyVector[1];
361 = exp->m_basisKeyVector[2];
366 (*m_exp).push_back(pyramid);
368 else if((HexGeom = boost::dynamic_pointer_cast<
372 = exp->m_basisKeyVector[0];
374 = exp->m_basisKeyVector[1];
376 = exp->m_basisKeyVector[2];
381 (*m_exp).push_back(hex);
385 ASSERTL0(
false,
"dynamic cast to a proper Geometry3D "
422 for(i = 0; i <
m_exp->size(); ++i)
426 m_offset_elmt_id[i] = i;
427 m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
428 m_npoints += (*m_exp)[i]->GetTotPoints();
437 Array<OneD, int> NumShape(4,0);
441 switch ((*
m_exp)[i]->DetShapeType())
448 ASSERTL0(
false,
"Unknown expansion type.");
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);
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]);
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;
480 for (i = 0; i < ntot; ++i)
482 for (j = 0; j < 3; ++j)
484 outfile << setprecision(8) << scientific
485 << (float)coords[j][i] <<
" ";
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)
498 for (j = 0; j < nquad1-1; ++j)
500 for (k = 0; k < nquad2-1; ++k)
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;
514 outfile <<
" </DataArray>" << endl;
515 outfile <<
" <DataArray type=\"Int32\" "
516 <<
"Name=\"offsets\" format=\"ascii\">" << endl;
517 for (i = 0; i < ntotminus; ++i)
519 outfile << i*8+8 <<
" ";
522 outfile <<
" </DataArray>" << endl;
523 outfile <<
" <DataArray type=\"UInt8\" "
524 <<
"Name=\"types\" format=\"ascii\">" << endl;
525 for (i = 0; i < ntotminus; ++i)
530 outfile <<
" </DataArray>" << endl;
531 outfile <<
" </Cells>" << endl;
532 outfile <<
" <PointData>" << endl;
538 for (i = 0; i <
m_exp->size(); ++i)
540 for (j = 0; j < (*m_exp)[i]->GetNfaces(); ++j)
542 (*m_exp)[i]->ComputeFaceNormal(j);
548 const Array<OneD, NekDouble> &inarray,
549 Array<OneD, NekDouble> &outarray)
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;
570 (*
m_exp)[i]->GetBasis(1)->GetPointsKey(),
571 (*
m_exp)[i]->GetBasis(2)->GetPointsKey(),
572 &inarray[cnt], newPointsKey0,
573 newPointsKey1, newPointsKey2,
577 cnt1 += npt0*npt1*npt2;
582 const Array<OneD, NekDouble> &inarray,
583 Array<OneD, NekDouble> &outarray)
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;
607 (*
m_exp)[i]->GetBasis(0)->GetPointsKey(),
608 (*
m_exp)[i]->GetBasis(1)->GetPointsKey(),
609 (*
m_exp)[i]->GetBasis(2)->GetPointsKey(),
612 cnt += npt0*npt1*npt2;