60 "Refines a prismatic or quadrilateral boundary layer.");
64 int **ret =
new int *[lda];
65 for (
int i = 0; i < lda; ++i)
68 ret[i][0] = arr[i][0];
69 ret[i][1] = arr[i][1];
76 int **ret =
new int *[lda];
77 for (
int i = 0; i < lda; ++i)
80 ret[i][0] = arr[i][0];
81 ret[i][1] = arr[i][1];
82 ret[i][2] = arr[i][2];
83 ret[i][3] = arr[i][3];
113 ConfigOption(
false,
"2",
"Number of layers to refine.");
115 ConfigOption(
false,
"5",
"Number of points in high order elements.");
117 ConfigOption(
false,
"",
"Tag identifying surface connected to prism.");
119 ConfigOption(
false,
"2.0",
"Ratio to use in geometry progression.");
130 cout <<
"ProcessBL: Refining boundary layer..." << endl;
132 int dim =
m_mesh->m_expDim;
148 ASSERTL0(0,
"Dimension not supported");
162 int nodeId =
m_mesh->m_vertexSet.size();
163 int nl =
m_config[
"layers"].as<
int>();
170 bool ratioIsString =
false;
172 if (
m_config[
"r"].isType<NekDouble>())
178 std::string rstr =
m_config[
"r"].as<
string>();
179 rExprId = rEval.DefineFunction(
"x y z", rstr);
180 ratioIsString =
true;
188 map<int, int> splitEls;
194 boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
195 boost::unordered_map<int, vector<NodeSharedPtr> >
::iterator eIt;
197 string surf =
m_config[
"surf"].as<
string>();
200 vector<unsigned int> surfs;
202 sort(surfs.begin(), surfs.end());
206 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
209 int nSurf = el->GetEdgeCount();
211 for (
int j = 0; j < nSurf; ++j)
213 int bl = el->GetBoundaryLink(j);
220 vector<int> tags = bEl->GetTagList();
223 sort(tags.begin(), tags.end());
224 set_intersection(surfs.begin(), surfs.end(),
225 tags .begin(), tags .end(),
226 back_inserter(inter));
228 "Intersection of surfaces wrong");
230 if (inter.size() == 1)
234 cerr <<
"WARNING: Found non-quad element "
235 <<
"to split in surface " << surf
236 <<
"; ignoring" << endl;
240 if (splitEls.count(el->GetId()) > 0)
242 cerr <<
"WARNING: quad already found; "
243 <<
"ignoring" << endl;
247 splitEls[el->GetId()] = j;
254 ASSERTL0(
false,
"Surface must be specified.");
257 if (splitEls.size() == 0)
259 cerr <<
"WARNING: No elements detected to split." << endl;
265 vector<ElementSharedPtr> el =
m_mesh->m_element[
m_mesh->m_expDim];
269 for (
int i = 0; i < el.size(); ++i)
273 if (splitEls.count(el[i]->GetId()) == 0)
280 std::map<int, int> bLink;
281 for (
int j = 0; j < 4; j++)
283 int bl = el[i]->GetBoundaryLink(j);
284 if ( (bl != -1) && ( j != splitEls[el[i]->GetId()]) )
293 el[i]->GetGeom(
m_mesh->m_spaceDim));
306 int nverts = geom->GetNumVerts();
310 for(
int i = 0; i < nverts; ++i)
312 geom->GetVertex(i)->GetCoords(x1,y1,z1);
313 x += x1; y += y1; z += z1;
318 r = rEval.Evaluate(rExprId,x,y,z,0.0);
331 if (splitEls[el[i]->GetId()] % 2)
348 vector<vector<NodeSharedPtr> > edgeNodes(2);
351 for (
int j = 0; j < 2; ++j)
353 int locEdge = (splitEls[el[i]->GetId()]+1+2*j)%4;
354 int edgeId = el[i]->GetEdge(locEdge)->m_id;
358 eIt = edgeMap.find(edgeId);
360 if (eIt == edgeMap.end())
363 edgeNodes[j].resize(nl+1);
367 edgeNodes[j][0] = el[i]->GetVertex(locEdge);
368 edgeNodes[j][nl] = el[i]->GetVertex((locEdge+1)%4);
378 x0 = edgeNodes[j][0]->m_x;
379 y0 = edgeNodes[j][0]->m_y;
381 x1 = edgeNodes[j][nl]->m_x;
382 y1 = edgeNodes[j][nl]->m_y;
389 rnew = rEval.Evaluate(rExprId,xm,ym,zm,0.0);
401 for (
int k = 1; k < nl; ++k)
403 xm = 0.5*(1+z[k])*(x1-x0) + x0;
404 ym = 0.5*(1+z[k])*(y1-y0) + y0;
407 new Node(nodeId++, xm,ym,zm));
414 for (
int k = 1; k < nl; ++k)
425 pos = nq*(nl+1) -1 - k;
431 ASSERTL0(0,
"Quad edge should be < 4.");
435 new Node(nodeId++, x[pos], y[pos], z[pos]));
440 edgeMap[edgeId] = edgeNodes[j];
445 if (eIt->second[0] == el[i]->GetVertex(locEdge))
448 edgeNodes[j] = eIt->second;
453 edgeNodes[j].resize(nl+1);
454 for (
int k = 0; k < nl+1; ++k)
456 edgeNodes[j][k] = eIt->second[nl-k];
464 for (
int j = 0; j < nl; ++j)
467 vector<NodeSharedPtr> nodeList(4);
468 switch (splitEls[el[i]->GetId()])
472 nodeList[0] = edgeNodes[1][nl-j ];
473 nodeList[1] = edgeNodes[0][j ];
474 nodeList[2] = edgeNodes[0][j+1];
475 nodeList[3] = edgeNodes[1][nl-j-1];
480 nodeList[0] = edgeNodes[1][j];
481 nodeList[1] = edgeNodes[1][j+1];
482 nodeList[2] = edgeNodes[0][nl-j-1];
483 nodeList[3] = edgeNodes[0][nl-j];
488 nodeList[0] = edgeNodes[0][nl-j ];
489 nodeList[1] = edgeNodes[1][j ];
490 nodeList[2] = edgeNodes[1][j+1];
491 nodeList[3] = edgeNodes[0][nl-j-1];
496 nodeList[0] = edgeNodes[0][j];
497 nodeList[1] = edgeNodes[0][j+1];
498 nodeList[2] = edgeNodes[1][nl-j-1];
499 nodeList[3] = edgeNodes[1][nl-j];
510 for (
int l = 0; l < 2; ++l)
512 int locEdge = (splitEls[el[i]->GetId()]+2*l)%4;
516 for (
int k = 1; k < nq-1; ++k)
524 pos = j+1 + k*(nl+1);
527 pos = (j+1)*nq + (nq-1) - k;
530 pos = (nl+1)*(nq-1) + j - k*(nl+1);
533 ASSERTL0(0,
"Quad edge should be < 4.");
536 HOedge->m_edgeNodes.push_back(
538 new Node(nodeId++,x[pos],y[pos],0.0)));
540 HOedge->m_curveType = pt;
546 for (it = bLink.begin(); it != bLink.end(); ++it)
555 for (
int k = 0; k < 2; ++k)
558 k, nodeList[(eid+k)%4]);
564 vector<NodeSharedPtr> qNodeList(2);
565 for (
int k = 0; k < 2; ++k)
567 qNodeList[k] = nodeList[(eid+k)%4];
570 tagBE =
m_mesh->m_element[
m_mesh->m_expDim-1][bl]->GetTagList();
575 m_mesh->m_element[
m_mesh->m_expDim-1].push_back(boundaryElmt);
587 set<LibUtilities::ShapeType> validElTypes;
591 int nodeId =
m_mesh->m_vertexSet.size();
592 int nl =
m_config[
"layers"].as<
int>();
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);
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);
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}};
671 splitHex0.
edge = splitMapEdgeHex0;
672 splitHex0.
offset = splitMapOffsetHex0;
673 splitHex0.
inc = splitMapIncHex0;
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};
697 splitPrismEdge.
size = 3;
698 splitPrismEdge.
edge = splitPrismEdges;
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};
710 splitHex0Edge.
size = 4;
711 splitHex0Edge.
edge = splitHex0Edges;
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())
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;
868 std::map<int, int> bLink;
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);
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]);
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);
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)
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.
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Principle Modified Functions .
pair< ModuleType, string > ModuleKey
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)
virtual NEKMESHUTILS_EXPORT void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
virtual void Process()
Write mesh to output file.
Principle Modified Functions .
virtual NEKMESHUTILS_EXPORT void ProcessElements()
Generate element IDs.
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
boost::shared_ptr< QuadExp > QuadExpSharedPtr
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Expansion > ExpansionSharedPtr
This class defines evaluator of analytic (symbolic) mathematical expressions. Expressions are allowed...
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< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Abstract base class for processing modules.
virtual NEKMESHUTILS_EXPORT void ProcessVertices()
Extract element vertices.
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
1D power law distribution for boundary layer points
virtual NEKMESHUTILS_EXPORT void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
boost::shared_ptr< Element > ElementSharedPtr
std::pair< ModuleType, std::string > ModuleKey
Describes the specification for a Basis.
virtual NEKMESHUTILS_EXPORT void ProcessComposites()
Generate composites.
int ** helper2d(int lda, int arr[][4])
1D Gauss-Lobatto-Legendre quadrature points
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr