587 set<LibUtilities::ShapeType> validElTypes;
591 int nodeId =
m_mesh->m_vertexSet.size();
592 int nl =
m_config[
"layers"].as<
int>();
596 LibUtilities::AnalyticExpressionEvaluator rEval;
599 bool ratioIsString =
false;
601 if (
m_config[
"r"].isType<NekDouble>())
607 std::string rstr =
m_config[
"r"].as<
string>();
608 rExprId = rEval.DefineFunction(
"x y z", rstr);
609 ratioIsString =
true;
613 int prismFaceNodes[5][4] = {
614 {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
615 int hexFaceNodes[6][4] = {{0, 1, 2, 3},
621 map<LibUtilities::ShapeType, int **> faceNodeMap;
630 boost::unordered_map<int, int> splitEls;
638 map<LibUtilities::ShapeType, map<int, SplitMapHelper> > splitMap;
639 int po = nq * (nl + 1);
641 SplitMapHelper splitPrism;
642 int splitMapEdgePrism[6] = {0, 2, 4, 5, 6, 7};
643 int splitMapOffsetPrism[6] = {0, nq, 0, nq - 1, nq + nq - 1, nq};
644 int splitMapIncPrism[6] = {1, 1, po, po, po, po};
645 int splitMapBFacesPrism[3] = {0, 2, 4};
646 int splitMapConnPrism[6][2] = {
647 {0, 0}, {1, 0}, {1, 1}, {0, 1}, {2, 0}, {2, 1}};
649 splitPrism.layerOff = nq;
650 splitPrism.edge = splitMapEdgePrism;
651 splitPrism.offset = splitMapOffsetPrism;
652 splitPrism.inc = splitMapIncPrism;
653 splitPrism.conn =
helper2d(6, splitMapConnPrism);
654 splitPrism.bfacesSize = 3;
655 splitPrism.bfaces = splitMapBFacesPrism;
659 int ho = nq * (nq - 1);
661 SplitMapHelper splitHex0;
662 int splitMapEdgeHex0[8] = {0, 1, 2, 3, 8, 9, 10, 11};
663 int splitMapOffsetHex0[8] = {
664 0, nq - 1, tl - 1, ho, tl, tl + nq - 1, 2 * tl - 1, tl + ho};
665 int splitMapIncHex0[8] = {1, nq, -1, -nq, 1, nq, -1, -nq};
666 int splitMapBFacesHex0[4] = {1, 2, 3, 4};
667 int splitMapConnHex0[8][2] = {
668 {0, 0}, {1, 0}, {2, 0}, {3, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}};
670 splitHex0.layerOff = nq * nq;
671 splitHex0.edge = splitMapEdgeHex0;
672 splitHex0.offset = splitMapOffsetHex0;
673 splitHex0.inc = splitMapIncHex0;
674 splitHex0.conn =
helper2d(8, splitMapConnHex0);
675 splitHex0.bfacesSize = 4;
676 splitHex0.bfaces = splitMapBFacesHex0;
690 map<LibUtilities::ShapeType, map<int, SplitEdgeHelper> > splitEdge;
692 int splitPrismEdges[3] = {3, 1, 8};
693 int splitPrismEdgeVert[3][2] = {{0, 3}, {1, 2}, {4, 5}};
694 int splitPrismOffset[3] = {0, nq - 1, nq * (nl + 1) * (nq - 1)};
695 int splitPrismInc[3] = {nq, nq, nq};
696 SplitEdgeHelper splitPrismEdge;
697 splitPrismEdge.size = 3;
698 splitPrismEdge.edge = splitPrismEdges;
699 splitPrismEdge.edgeVert =
helper2d(3, splitPrismEdgeVert);
700 splitPrismEdge.offset = splitPrismOffset;
701 splitPrismEdge.inc = splitPrismInc;
705 int splitHex0Edges[4] = {4, 5, 6, 7};
706 int splitHex0EdgeVert[4][2] = {{0, 4}, {1, 5}, {2, 6}, {3, 7}};
707 int splitHex0Offset[4] = {0, nq - 1, nq * nq - 1, nq * (nq - 1)};
708 int splitHex0Inc[4] = {nq * nq, nq * nq, nq * nq, nq * nq};
709 SplitEdgeHelper splitHex0Edge;
710 splitHex0Edge.size = 4;
711 splitHex0Edge.edge = splitHex0Edges;
712 splitHex0Edge.edgeVert =
helper2d(4, splitHex0EdgeVert);
713 splitHex0Edge.offset = splitHex0Offset;
714 splitHex0Edge.inc = splitHex0Inc;
718 map<LibUtilities::ShapeType, map<int, bool> > revPoints;
729 boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
730 boost::unordered_map<int, vector<NodeSharedPtr> >
::iterator eIt;
732 string surf =
m_config[
"surf"].as<
string>();
735 vector<unsigned int> surfs;
737 sort(surfs.begin(), surfs.end());
741 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
744 int nSurf = el->GetFaceCount();
746 for (
int j = 0; j < nSurf; ++j)
748 int bl = el->GetBoundaryLink(j);
756 vector<int> tags = bEl->GetTagList();
759 sort(tags.begin(), tags.end());
760 set_intersection(surfs.begin(),
764 back_inserter(inter));
765 ASSERTL0(inter.size() <= 1,
"Intersection of surfaces wrong");
767 if (inter.size() == 1)
773 cerr <<
"WARNING: Found quadrilateral face " << j
774 <<
" on surface " << surf
775 <<
" connected to prism; ignoring." << endl;
779 if (splitEls.count(el->GetId()) > 0)
781 cerr <<
"WARNING: prism already found; "
782 <<
"ignoring" << endl;
785 splitEls[el->GetId()] = j;
787 else if (validElTypes.count(el->GetConf().m_e) == 0)
789 cerr <<
"WARNING: Unsupported element type "
790 <<
"found in surface " << j <<
"; "
791 <<
"ignoring" << endl;
802 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
808 splitEls[el->GetId()] = 1;
810 else if (validElTypes.count(el->GetConf().m_e) > 0)
812 splitEls[el->GetId()] = 0;
821 if (splitEls.size() == 0)
823 cerr <<
"WARNING: No elements detected to split." << endl;
829 vector<ElementSharedPtr> el =
m_mesh->m_element[
m_mesh->m_expDim];
832 map<int, SpatialDomains::Geometry3DSharedPtr> geomMap;
833 for (
int i = 0; i < el.size(); ++i)
835 const int elId = el[i]->GetId();
836 sIt = splitEls.find(elId);
837 if (sIt == splitEls.end())
843 geomMap[elId] = boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
844 el[i]->GetGeom(
m_mesh->m_spaceDim));
848 for (
int i = 0; i < el.size(); ++i)
850 const int elId = el[i]->GetId();
851 sIt = splitEls.find(elId);
853 if (sIt == splitEls.end())
861 const int faceNum = sIt->second;
864 SplitMapHelper &sMap = splitMap[elType][faceNum];
865 SplitEdgeHelper &sEdge = splitEdge[elType][faceNum];
868 std::map<int, int> bLink;
869 for (
int j = 0; j < sMap.bfacesSize; ++j)
871 int bl = el[i]->GetBoundaryLink(sMap.bfaces[j]);
874 bLink[sMap.bfaces[j]] = bl;
888 int nverts = geom->GetNumVerts();
892 for (
int i = 0; i < nverts; ++i)
894 geom->GetVertex(i)->GetCoords(x1, y1, z1);
902 r = rEval.Evaluate(rExprId, x, y, z, 0.0);
910 LibUtilities::BasisKey B0(
914 LibUtilities::PointsKey(nl + 1, t, r));
915 LibUtilities::BasisKey B2(
920 boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(geom);
927 LibUtilities::BasisKey B0(
931 LibUtilities::PointsKey(nl + 1, t, r));
935 boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(geom);
941 Array<OneD, NekDouble> x(nq * nq * (nl + 1));
942 Array<OneD, NekDouble> y(nq * nq * (nl + 1));
943 Array<OneD, NekDouble> z(nq * nq * (nl + 1));
944 q->GetCoords(x, y, z);
946 int nSplitEdge = sEdge.size;
947 vector<vector<NodeSharedPtr> > edgeNodes(nSplitEdge);
950 for (
int j = 0; j < nSplitEdge; ++j)
952 int locEdge = sEdge.edge[j];
953 int edgeId = el[i]->GetEdge(locEdge)->m_id;
957 eIt = edgeMap.find(edgeId);
959 if (eIt == edgeMap.end())
962 edgeNodes[j].resize(nl + 1);
966 edgeNodes[j][0] = el[i]->GetVertex(sEdge.edgeVert[j][0]);
967 edgeNodes[j][nl] = el[i]->GetVertex(sEdge.edgeVert[j][1]);
977 x0 = x[sEdge.offset[j]];
978 y0 = y[sEdge.offset[j]];
979 z0 = z[sEdge.offset[j]];
981 x1 = x[sEdge.offset[j] + nl * nq];
982 y1 = y[sEdge.offset[j] + nl * nq];
983 z1 = z[sEdge.offset[j] + nl * nq];
985 xm = 0.5 * (x0 + x1);
986 ym = 0.5 * (y0 + y1);
987 zm = 0.5 * (z0 + z1);
991 rnew = rEval.Evaluate(rExprId, xm, ym, zm, 0.0);
994 LibUtilities::PointsKey Pkey(nl + 1, t, rnew);
998 const Array<OneD, const NekDouble> z = newP->GetZ();
1001 for (
int k = 1; k < nl; ++k)
1003 xm = 0.5 * (1 + z[k]) * (x1 - x0) + x0;
1004 ym = 0.5 * (1 + z[k]) * (y1 - y0) + y0;
1005 zm = 0.5 * (1 + z[k]) * (z1 - z0) + z0;
1013 for (
int k = 1; k < nl; ++k)
1015 int pos = sEdge.offset[j] + k * sEdge.inc[j];
1017 new Node(nodeId++, x[pos], y[pos], z[pos]));
1022 edgeMap[edgeId] = edgeNodes[j];
1026 edgeNodes[j] = eIt->second;
1031 for (
int j = 0; j < nl; ++j)
1035 int offset = j * sMap.layerOff;
1038 vector<NodeSharedPtr> nodeList(sMap.size);
1039 for (
int k = 0; k < sMap.size; ++k)
1041 nodeList[k] = edgeNodes[sMap.conn[k][0]][j + sMap.conn[k][1]];
1045 ElmtConfig conf(elType, 1,
true,
true,
false);
1047 elType, conf, nodeList, el[i]->GetTagList());
1050 for (
int l = 0; l < sMap.size; ++l)
1053 for (
int k = 1; k < nq - 1; ++k)
1055 int pos = offset + sMap.offset[l] + k * sMap.inc[l];
1057 new Node(nodeId++, x[pos], y[pos], z[pos])));
1059 HOedge->m_curveType = pt;
1065 for (it = bLink.begin(); it != bLink.end(); ++it)
1067 int fid = it->first;
1068 int bl = it->second;
1070 vector<NodeSharedPtr> qNodeList(4);
1071 for (
int k = 0; k < 4; ++k)
1073 qNodeList[k] = nodeList[faceNodeMap[elType][fid][k]];
1077 m_mesh->m_element[
m_mesh->m_expDim - 1][bl]->GetTagList();
1089 m_mesh->m_element[
m_mesh->m_expDim - 1][bl] = boundaryElmt;
#define ASSERTL0(condition, msg)
Basic information about an element.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Principle Modified Functions .
boost::shared_ptr< Points< NekDouble > > PointsSharedPtr
MeshSharedPtr m_mesh
Mesh object.
boost::shared_ptr< HexGeom > HexGeomSharedPtr
int ** helper2d(int lda, int arr[][2])
ElementFactory & GetElementFactory()
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Principle Modified Functions .
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Expansion > ExpansionSharedPtr
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
1D power law distribution for boundary layer points
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
1D power law distribution for boundary layer points
boost::shared_ptr< Element > ElementSharedPtr
1D Gauss-Lobatto-Legendre quadrature points
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr