Write mesh to output file.
129 cout <<
"ProcessBL: Refining prismatic boundary layer..." << endl;
133 set<LibUtilities::ShapeType> validElTypes;
137 int nodeId =
m_mesh->m_vertexSet.size();
138 int nl =
m_config[
"layers"].as<
int>();
142 LibUtilities::AnalyticExpressionEvaluator rEval;
145 bool ratioIsString =
false;
147 if (
m_config[
"r"].isType<NekDouble>())
153 std::string rstr =
m_config[
"r"].as<
string>();
154 rExprId = rEval.DefineFunction(
"x y z", rstr);
155 ratioIsString =
true;
159 int prismFaceNodes[5][4] = {
160 {0, 1, 2, 3}, {0, 1, 4, -1}, {1, 2, 5, 4}, {3, 2, 5, -1}, {0, 3, 5, 4}};
161 int hexFaceNodes[6][4] = {{0, 1, 2, 3},
167 map<LibUtilities::ShapeType, int **> faceNodeMap;
176 boost::unordered_map<int, int> splitEls;
184 map<LibUtilities::ShapeType, map<int, SplitMapHelper> > splitMap;
185 int po = nq * (nl + 1);
187 SplitMapHelper splitPrism;
188 int splitMapEdgePrism[6] = {0, 2, 4, 5, 6, 7};
189 int splitMapOffsetPrism[6] = {0, nq, 0, nq - 1, nq + nq - 1, nq};
190 int splitMapIncPrism[6] = {1, 1, po, po, po, po};
191 int splitMapBFacesPrism[3] = {0, 2, 4};
192 int splitMapConnPrism[6][2] = {
193 {0, 0}, {1, 0}, {1, 1}, {0, 1}, {2, 0}, {2, 1}};
195 splitPrism.layerOff = nq;
196 splitPrism.edge = splitMapEdgePrism;
197 splitPrism.offset = splitMapOffsetPrism;
198 splitPrism.inc = splitMapIncPrism;
199 splitPrism.conn =
helper2d(6, splitMapConnPrism);
200 splitPrism.bfacesSize = 3;
201 splitPrism.bfaces = splitMapBFacesPrism;
205 int ho = nq * (nq - 1);
207 SplitMapHelper splitHex0;
208 int splitMapEdgeHex0[8] = {0, 1, 2, 3, 8, 9, 10, 11};
209 int splitMapOffsetHex0[8] = {
210 0, nq - 1, tl - 1, ho, tl, tl + nq - 1, 2 * tl - 1, tl + ho};
211 int splitMapIncHex0[8] = {1, nq, -1, -nq, 1, nq, -1, -nq};
212 int splitMapBFacesHex0[4] = {1, 2, 3, 4};
213 int splitMapConnHex0[8][2] = {
214 {0, 0}, {1, 0}, {2, 0}, {3, 0}, {0, 1}, {1, 1}, {2, 1}, {3, 1}};
216 splitHex0.layerOff = nq * nq;
217 splitHex0.edge = splitMapEdgeHex0;
218 splitHex0.offset = splitMapOffsetHex0;
219 splitHex0.inc = splitMapIncHex0;
220 splitHex0.conn =
helper2d(8, splitMapConnHex0);
221 splitHex0.bfacesSize = 4;
222 splitHex0.bfaces = splitMapBFacesHex0;
236 map<LibUtilities::ShapeType, map<int, SplitEdgeHelper> > splitEdge;
238 int splitPrismEdges[3] = {3, 1, 8};
239 int splitPrismEdgeVert[3][2] = {{0, 3}, {1, 2}, {4, 5}};
240 int splitPrismOffset[3] = {0, nq - 1, nq * (nl + 1) * (nq - 1)};
241 int splitPrismInc[3] = {nq, nq, nq};
242 SplitEdgeHelper splitPrismEdge;
243 splitPrismEdge.size = 3;
244 splitPrismEdge.edge = splitPrismEdges;
245 splitPrismEdge.edgeVert =
helper2d(3, splitPrismEdgeVert);
246 splitPrismEdge.offset = splitPrismOffset;
247 splitPrismEdge.inc = splitPrismInc;
251 int splitHex0Edges[4] = {4, 5, 6, 7};
252 int splitHex0EdgeVert[4][2] = {{0, 4}, {1, 5}, {2, 6}, {3, 7}};
253 int splitHex0Offset[4] = {0, nq - 1, nq * nq - 1, nq * (nq - 1)};
254 int splitHex0Inc[4] = {nq * nq, nq * nq, nq * nq, nq * nq};
255 SplitEdgeHelper splitHex0Edge;
256 splitHex0Edge.size = 4;
257 splitHex0Edge.edge = splitHex0Edges;
258 splitHex0Edge.edgeVert =
helper2d(4, splitHex0EdgeVert);
259 splitHex0Edge.offset = splitHex0Offset;
260 splitHex0Edge.inc = splitHex0Inc;
264 map<LibUtilities::ShapeType, map<int, bool> > revPoints;
275 boost::unordered_map<int, vector<NodeSharedPtr> > edgeMap;
276 boost::unordered_map<int, vector<NodeSharedPtr> >
::iterator eIt;
281 vector<unsigned int> surfs;
283 sort(surfs.begin(), surfs.end());
287 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
290 int nSurf = el->GetFaceCount();
292 for (
int j = 0; j < nSurf; ++j)
294 int bl = el->GetBoundaryLink(j);
302 vector<int> tags = bEl->GetTagList();
305 sort(tags.begin(), tags.end());
306 set_intersection(surfs.begin(),
310 back_inserter(inter));
311 ASSERTL0(inter.size() <= 1,
"Intersection of surfaces wrong");
313 if (inter.size() == 1)
319 cerr <<
"WARNING: Found quadrilateral face " << j
320 <<
" on surface " << surf
321 <<
" connected to prism; ignoring." << endl;
325 if (splitEls.count(el->GetId()) > 0)
327 cerr <<
"WARNING: prism already found; "
328 <<
"ignoring" << endl;
331 splitEls[el->GetId()] = j;
333 else if (validElTypes.count(el->GetConf().m_e) == 0)
335 cerr <<
"WARNING: Unsupported element type "
336 <<
"found in surface " << j <<
"; "
337 <<
"ignoring" << endl;
348 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
354 splitEls[el->GetId()] = 1;
356 else if (validElTypes.count(el->GetConf().m_e) > 0)
358 splitEls[el->GetId()] = 0;
367 if (splitEls.size() == 0)
369 cerr <<
"WARNING: No elements detected to split." << endl;
375 vector<ElementSharedPtr> el =
m_mesh->m_element[
m_mesh->m_expDim];
378 map<int, SpatialDomains::Geometry3DSharedPtr> geomMap;
379 for (
int i = 0; i < el.size(); ++i)
381 const int elId = el[i]->GetId();
382 sIt = splitEls.find(elId);
383 if (sIt == splitEls.end())
389 geomMap[elId] = boost::dynamic_pointer_cast<SpatialDomains::Geometry3D>(
390 el[i]->GetGeom(
m_mesh->m_spaceDim));
394 for (
int i = 0; i < el.size(); ++i)
396 const int elId = el[i]->GetId();
397 sIt = splitEls.find(elId);
399 if (sIt == splitEls.end())
407 const int faceNum = sIt->second;
410 SplitMapHelper &sMap = splitMap[elType][faceNum];
411 SplitEdgeHelper &sEdge = splitEdge[elType][faceNum];
414 std::map<int, int> bLink;
415 for (
int j = 0; j < sMap.bfacesSize; ++j)
417 int bl = el[i]->GetBoundaryLink(sMap.bfaces[j]);
420 bLink[sMap.bfaces[j]] = bl;
434 int nverts = geom->GetNumVerts();
438 for (
int i = 0; i < nverts; ++i)
440 geom->GetVertex(i)->GetCoords(x1, y1, z1);
448 r = rEval.Evaluate(rExprId, x, y, z, 0.0);
456 LibUtilities::BasisKey B0(
460 LibUtilities::PointsKey(nl + 1, t, r));
461 LibUtilities::BasisKey B2(
466 boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>(geom);
473 LibUtilities::BasisKey B0(
477 LibUtilities::PointsKey(nl + 1, t, r));
481 boost::dynamic_pointer_cast<SpatialDomains::HexGeom>(geom);
487 Array<OneD, NekDouble> x(nq * nq * (nl + 1));
488 Array<OneD, NekDouble> y(nq * nq * (nl + 1));
489 Array<OneD, NekDouble> z(nq * nq * (nl + 1));
490 q->GetCoords(x, y, z);
492 int nSplitEdge = sEdge.size;
493 vector<vector<NodeSharedPtr> > edgeNodes(nSplitEdge);
496 for (
int j = 0; j < nSplitEdge; ++j)
498 int locEdge = sEdge.edge[j];
499 int edgeId = el[i]->GetEdge(locEdge)->m_id;
503 eIt = edgeMap.find(edgeId);
505 if (eIt == edgeMap.end())
508 edgeNodes[j].resize(nl + 1);
512 edgeNodes[j][0] = el[i]->GetVertex(sEdge.edgeVert[j][0]);
513 edgeNodes[j][nl] = el[i]->GetVertex(sEdge.edgeVert[j][1]);
523 x0 = x[sEdge.offset[j]];
524 y0 = y[sEdge.offset[j]];
525 z0 = z[sEdge.offset[j]];
527 x1 = x[sEdge.offset[j] + nl * nq];
528 y1 = y[sEdge.offset[j] + nl * nq];
529 z1 = z[sEdge.offset[j] + nl * nq];
531 xm = 0.5 * (x0 + x1);
532 ym = 0.5 * (y0 + y1);
533 zm = 0.5 * (z0 + z1);
537 rnew = rEval.Evaluate(rExprId, xm, ym, zm, 0.0);
540 LibUtilities::PointsKey Pkey(nl + 1, t, rnew);
544 const Array<OneD, const NekDouble> z = newP->GetZ();
547 for (
int k = 1; k < nl; ++k)
549 xm = 0.5 * (1 + z[k]) * (x1 - x0) + x0;
550 ym = 0.5 * (1 + z[k]) * (y1 - y0) + y0;
551 zm = 0.5 * (1 + z[k]) * (z1 - z0) + z0;
559 for (
int k = 1; k < nl; ++k)
561 int pos = sEdge.offset[j] + k * sEdge.inc[j];
563 new Node(nodeId++, x[pos], y[pos], z[pos]));
568 edgeMap[edgeId] = edgeNodes[j];
572 edgeNodes[j] = eIt->second;
577 for (
int j = 0; j < nl; ++j)
581 int offset = j * sMap.layerOff;
584 vector<NodeSharedPtr> nodeList(sMap.size);
585 for (
int k = 0; k < sMap.size; ++k)
587 nodeList[k] = edgeNodes[sMap.conn[k][0]][j + sMap.conn[k][1]];
591 ElmtConfig conf(elType, 1,
true,
true,
false);
593 elType, conf, nodeList, el[i]->GetTagList());
596 for (
int l = 0; l < sMap.size; ++l)
599 for (
int k = 1; k < nq - 1; ++k)
601 int pos = offset + sMap.offset[l] + k * sMap.inc[l];
603 new Node(nodeId++, x[pos], y[pos], z[pos])));
605 HOedge->m_curveType = pt;
611 for (it = bLink.begin(); it != bLink.end(); ++it)
616 vector<NodeSharedPtr> qNodeList(4);
617 for (
int k = 0; k < 4; ++k)
619 qNodeList[k] = nodeList[faceNodeMap[elType][fid][k]];
623 m_mesh->m_element[
m_mesh->m_expDim - 1][bl]->GetTagList();
635 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.
map< string, ConfigOption > m_config
List of configuration values.
Principle Modified Functions .
MeshSharedPtr m_mesh
Mesh object.
boost::shared_ptr< Points< NekDouble > > PointsSharedPtr
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)
Represents a point in the domain.
virtual void ProcessEdges(bool ReprocessEdges=true)
Extract element edges.
virtual void ProcessVertices()
Extract element vertices.
virtual void ProcessElements()
Generate element IDs.
Principle Modified Functions .
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
boost::shared_ptr< Expansion > ExpansionSharedPtr
virtual void ProcessComposites()
Generate composites.
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
virtual void ProcessFaces(bool ReprocessFaces=true)
Extract element faces.
1D power law distribution for boundary layer points
boost::shared_ptr< Element > ElementSharedPtr
1D Gauss-Lobatto-Legendre quadrature points
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr