36 #include <boost/geometry.hpp>
37 #include <boost/geometry/geometries/box.hpp>
38 #include <boost/geometry/geometries/point.hpp>
39 #include <boost/geometry/index/rtree.hpp>
51 namespace bg = boost::geometry;
52 namespace bgi = boost::geometry::index;
54 typedef bg::model::point<double, 3, bg::cs::cartesian>
point;
55 typedef bg::model::box<point>
box;
56 typedef std::pair<box, unsigned int>
boxI;
61 namespace NekMeshUtils
66 NekDouble xmin = numeric_limits<double>::max(),
67 xmax = -1.0 * numeric_limits<double>::max(),
68 ymin = numeric_limits<double>::max(),
69 ymax = -1.0 * numeric_limits<double>::max(),
70 zmin = numeric_limits<double>::max(),
71 zmax = -1.0 * numeric_limits<double>::max();
73 vector<NodeSharedPtr> ns = el->GetVertexList();
74 for (
int i = 0; i < ns.size(); i++)
76 xmin = min(xmin, ns[i]->m_x);
77 xmax = max(xmax, ns[i]->m_x);
78 ymin = min(ymin, ns[i]->m_y);
79 ymax = max(ymax, ns[i]->m_y);
80 zmin = min(zmin, ns[i]->m_z);
81 zmax = max(zmax, ns[i]->m_z);
84 return box(
point(xmin - ov, ymin - ov, zmin - ov),
85 point(xmax + ov, ymax + ov, zmax + ov));
90 NekDouble xmin = numeric_limits<double>::max(),
91 xmax = -1.0 * numeric_limits<double>::max(),
92 ymin = numeric_limits<double>::max(),
93 ymax = -1.0 * numeric_limits<double>::max(),
94 zmin = numeric_limits<double>::max(),
95 zmax = -1.0 * numeric_limits<double>::max();
97 for (
int j = 0; j < els.size(); j++)
99 vector<NodeSharedPtr> ns = els[j]->GetVertexList();
100 for (
int i = 0; i < ns.size(); i++)
102 xmin = min(xmin, ns[i]->m_x);
103 xmax = max(xmax, ns[i]->m_x);
104 ymin = min(ymin, ns[i]->m_y);
105 ymax = max(ymax, ns[i]->m_y);
106 zmin = min(zmin, ns[i]->m_z);
107 zmax = max(zmax, ns[i]->m_z);
111 return box(
point(xmin - ov, ymin - ov, zmin - ov),
112 point(xmax + ov, ymax + ov, zmax + ov));
117 return box(
point(n->m_x - ov, n->m_y - ov, n->m_z - ov),
118 point(n->m_x + ov, n->m_y + ov, n->m_z + ov));
132 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
134 vector<blInfoSharedPtr> infos = m_nToNInfo[bit->first];
135 for (
int i = 0; i < infos.size(); i++)
137 if (bit->second->bl > infos[i]->bl + 1)
139 cout <<
"non smooth error " << bit->second->bl <<
" "
140 << infos[i]->bl << endl;
145 for (
int i = 0; i < m_mesh->m_element[3].size(); i++)
148 if (!IsPrismValid(el))
150 cout <<
"validity error " << el->GetId() << endl;
159 map<NodeSharedPtr, NodeSharedPtr> BLMesh::GetSymNodes()
161 map<NodeSharedPtr, NodeSharedPtr> ret;
164 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
166 if (!bit->second->onSym)
174 bit->second->pNode->SetCADSurf(bit->second->symsurf, s, uv);
175 ret[bit->first] = bit->second->pNode;
182 vector<NodeSharedPtr> ns1 = el->GetVertexList();
186 V[0] = n->m_x - ns1[0]->m_x;
187 V[1] = n->m_y - ns1[0]->m_y;
188 V[2] = n->m_z - ns1[0]->m_z;
190 NekDouble Vmag = sqrt(V[0] * V[0] + V[1] * V[1] + V[2] * V[2]);
192 NekDouble ang = (N1[0] * V[0] + N1[1] * V[1] + N1[2] * V[2]) / Vmag;
197 void BLMesh::GrowLayers()
212 map<int, vector<ElementSharedPtr> > psElements;
213 for (
int i = 0; i < m_mesh->m_element[2].size(); i++)
217 find(m_blsurfs.begin(), m_blsurfs.end(), el->m_parentCAD->GetId());
220 m_symSurfs.begin(), m_symSurfs.end(), el->m_parentCAD->GetId());
222 if (f == m_blsurfs.end() && s == m_symSurfs.end())
224 psElements[el->m_parentCAD->GetId()].push_back(el);
227 for (
int i = 0; i < m_psuedoSurface.size(); i++)
229 psElements[m_psuedoSurface[i]->m_parentCAD->GetId()].push_back(
233 bgi::rtree<boxI, bgi::quadratic<16> > TopTree;
234 map<int, bgi::rtree<boxI, bgi::quadratic<16> > > SubTrees;
240 for (
int l = 1; l < m_layer; l++)
242 NekDouble delta = (m_layerT[l] - m_layerT[l - 1]);
245 map<int, vector<ElementSharedPtr> >
::iterator it;
246 for (it = psElements.begin(); it != psElements.end(); it++)
248 TopTree.insert(make_pair(
GetBox(it->second, m_bl), it->first));
249 vector<boxI> toInsert;
250 for (
int i = 0; i < it->second.size(); i++)
252 toInsert.push_back(make_pair(
GetBox(it->second[i], m_bl), i));
254 SubTrees[it->first].insert(toInsert.begin(), toInsert.end());
257 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
259 if (bit->second->stopped)
264 vector<boxI> results;
265 TopTree.query(bgi::intersects(
point(bit->second->pNode->m_x,
266 bit->second->pNode->m_y,
267 bit->second->pNode->m_z)),
268 back_inserter(results));
270 for (
int i = 0; i < results.size(); i++)
273 bit->second->surfs.find(results[i].second);
274 if (f == bit->second->surfs.end())
277 surfs.insert(results[i].second);
283 for (iit = surfs.begin(); iit != surfs.end(); iit++)
286 SubTrees[*iit].query(
287 bgi::intersects(
GetBox(bit->second->pNode, m_bl)),
288 back_inserter(results));
289 for (
int i = 0; i < results.size(); i++)
291 if (
Infont(bit->second->pNode,
292 psElements[*iit][results[i].second]))
295 Proximity(bit->second->pNode,
296 psElements[*iit][results[i].second]);
297 if (prox < delta * 2.5)
301 bit->second->stopped =
true;
321 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
323 if (bit->second->stopped)
329 bool shouldStop =
false;
330 vector<blInfoSharedPtr> ne = m_nToNInfo[bit->first];
331 for (
int i = 0; i < ne.size(); i++)
333 if (ne[i]->bl < bit->second->bl)
341 bit->second->stopped =
true;
345 bit->second->AlignNode(m_layerT[l]);
354 return (a * b > 0.0);
359 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
364 vector<NodeSharedPtr> ns = el->GetVertexList();
367 E0[0] = ns[1]->m_x - ns[0]->m_x;
368 E0[1] = ns[1]->m_y - ns[0]->m_y;
369 E0[2] = ns[1]->m_z - ns[0]->m_z;
371 E1[0] = ns[2]->m_x - ns[0]->m_x;
372 E1[1] = ns[2]->m_y - ns[0]->m_y;
373 E1[2] = ns[2]->m_z - ns[0]->m_z;
387 NekDouble s = b * e - c * d, t = b * d - a * e;
401 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e / c));
407 s = (d >= 0 ? 0 : (-d >= a ? 1 : -d / a));
438 s = (numer >= denom ? 1 : numer / denom);
445 new Node(0, B[0] + s * E0[0] + t * E1[0], B[1] + s * E0[1] + t * E1[1],
446 B[2] + s * E0[2] + t * E1[2]));
448 return n->Distance(point);
453 vector<NodeSharedPtr> ns1 = e1->GetVertexList();
454 vector<NodeSharedPtr> ns2 = e2->GetVertexList();
455 if (ns1[0] == ns2[0] || ns1[0] == ns2[1] || ns1[0] == ns2[2] ||
456 ns1[1] == ns2[0] || ns1[1] == ns2[1] || ns1[1] == ns2[2] ||
457 ns1[2] == ns2[0] || ns1[2] == ns2[1] || ns1[2] == ns2[2])
465 N1[0] = (ns1[1]->m_y - ns1[0]->m_y) * (ns1[2]->m_z - ns1[0]->m_z) -
466 (ns1[2]->m_y - ns1[0]->m_y) * (ns1[1]->m_z - ns1[0]->m_z);
467 N1[1] = -1.0 * ((ns1[1]->m_x - ns1[0]->m_x) * (ns1[2]->m_z - ns1[0]->m_z) -
468 (ns1[2]->m_x - ns1[0]->m_x) * (ns1[1]->m_z - ns1[0]->m_z));
469 N1[2] = (ns1[1]->m_x - ns1[0]->m_x) * (ns1[2]->m_y - ns1[0]->m_y) -
470 (ns1[2]->m_x - ns1[0]->m_x) * (ns1[1]->m_y - ns1[0]->m_y);
472 N2[0] = (ns2[1]->m_y - ns2[0]->m_y) * (ns2[2]->m_z - ns2[0]->m_z) -
473 (ns2[2]->m_y - ns2[0]->m_y) * (ns2[1]->m_z - ns2[0]->m_z);
474 N2[1] = -1.0 * ((ns2[1]->m_x - ns2[0]->m_x) * (ns2[2]->m_z - ns2[0]->m_z) -
475 (ns2[2]->m_x - ns2[0]->m_x) * (ns2[1]->m_z - ns2[0]->m_z));
476 N2[2] = (ns2[1]->m_x - ns2[0]->m_x) * (ns2[2]->m_y - ns2[0]->m_y) -
477 (ns2[2]->m_x - ns2[0]->m_x) * (ns2[1]->m_y - ns2[0]->m_y);
480 (N1[0] * ns1[0]->m_x + N1[1] * ns1[0]->m_y + N1[2] * ns1[0]->m_z);
482 (N2[0] * ns2[0]->m_x + N2[1] * ns2[0]->m_y + N2[2] * ns2[0]->m_z);
487 N2[0] * ns1[0]->m_x + N2[1] * ns1[0]->m_y + N2[2] * ns1[0]->m_z + d2;
489 N2[0] * ns1[1]->m_x + N2[1] * ns1[1]->m_y + N2[2] * ns1[1]->m_z + d2;
491 N2[0] * ns1[2]->m_x + N2[1] * ns1[2]->m_y + N2[2] * ns1[2]->m_z + d2;
494 N1[0] * ns2[0]->m_x + N1[1] * ns2[0]->m_y + N1[2] * ns2[0]->m_z + d1;
496 N1[0] * ns2[1]->m_x + N1[1] * ns2[1]->m_y + N1[2] * ns2[1]->m_z + d1;
498 N1[0] * ns2[2]->m_x + N1[1] * ns2[2]->m_y + N1[2] * ns2[2]->m_z + d1;
500 if (
sign(dv1[0], dv1[1]) &&
sign(dv1[1], dv1[2]))
504 if (
sign(dv2[0], dv2[1]) &&
sign(dv2[1], dv2[2]))
510 D[0] = N1[1] * N2[2] - N1[2] * N2[1];
511 D[1] = -1.0 * (N1[0] * N2[2] - N1[2] * N2[0]);
512 D[2] = N1[0] * N2[1] - N1[0] * N2[1];
514 int base1 = 0, base2 = 0;
515 if (!
sign(dv2[0], dv2[1]) &&
sign(dv2[1], dv2[2]))
519 else if (!
sign(dv2[1], dv2[2]) &&
sign(dv2[2], dv2[0]))
523 else if (!
sign(dv2[2], dv2[0]) &&
sign(dv2[0], dv2[1]))
529 cout <<
"base not set" << endl;
532 if (!
sign(dv1[0], dv1[1]) &&
sign(dv1[1], dv1[2]))
536 else if (!
sign(dv1[1], dv1[2]) &&
sign(dv1[2], dv1[0]))
540 else if (!
sign(dv1[2], dv1[0]) &&
sign(dv1[0], dv1[1]))
546 cout <<
"base not set" << endl;
551 p1[0] = D[0] * ns1[0]->m_x + D[1] * ns1[0]->m_y + D[2] * ns1[0]->m_z;
552 p1[1] = D[0] * ns1[1]->m_x + D[1] * ns1[1]->m_y + D[2] * ns1[1]->m_z;
553 p1[2] = D[0] * ns1[2]->m_x + D[1] * ns1[2]->m_y + D[2] * ns1[2]->m_z;
555 p2[0] = D[0] * ns2[0]->m_x + D[1] * ns2[0]->m_y + D[2] * ns2[0]->m_z;
556 p2[1] = D[0] * ns2[1]->m_x + D[1] * ns2[1]->m_y + D[2] * ns2[1]->m_z;
557 p2[2] = D[0] * ns2[2]->m_x + D[1] * ns2[2]->m_y + D[2] * ns2[2]->m_z;
577 t11 = p1[o1] + (p1[base1] - p1[o1]) * dv1[o1] / (dv1[o1] - dv1[base1]);
578 t12 = p1[o2] + (p1[base1] - p1[o2]) * dv1[o2] / (dv1[o2] - dv1[base1]);
596 t21 = p2[o1] + (p2[base2] - p2[o1]) * dv2[o1] / (dv2[o1] - dv2[base2]);
597 t22 = p2[o2] + (p2[base2] - p2[o2]) * dv2[o2] / (dv2[o2] - dv2[base2]);
614 if (!
sign(t21 - t11, t22 - t11) || !
sign(t21 - t12, t22 - t12))
622 void BLMesh::Shrink()
630 vector<ElementSharedPtr> inv;
631 for (
int i = 0; i < m_mesh->m_element[3].size(); i++)
634 if (!IsPrismValid(el))
640 smsh = (inv.size() > 0);
642 for (
int i = 0; i < inv.size(); i++)
645 vector<blInfoSharedPtr> bls;
646 vector<NodeSharedPtr> ns = t->GetVertexList();
647 for (
int j = 0; j < ns.size(); j++)
649 bls.push_back(m_blData[ns[j]]);
656 for (
int j = 0; j < 3; j++)
658 mx = max(mx, bls[j]->bl);
660 ASSERTL0(mx > 0,
"shrinking to nothing");
661 for (
int j = 0; j < 3; j++)
668 bls[j]->AlignNode(m_layerT[bls[j]->bl]);
670 if (!IsPrismValid(inv[i]))
681 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
683 vector<blInfoSharedPtr> infos = m_nToNInfo[bit->first];
684 for (
int i = 0; i < infos.size(); i++)
686 if (bit->second->bl > infos[i]->bl + 1)
689 bit->second->AlignNode(m_layerT[bit->second->bl]);
700 NekDouble mn = numeric_limits<double>::max();
701 NekDouble mx = -1.0 * numeric_limits<double>::max();
702 vector<NodeSharedPtr> ns = el->GetVertexList();
704 for (
int j = 0; j < ns.size(); j++)
723 for (
int j = 0; j < 6; j++)
737 dxdz(0, 0) * (dxdz(1, 1) * dxdz(2, 2) - dxdz(2, 1) * dxdz(1, 2)) -
738 dxdz(0, 1) * (dxdz(1, 0) * dxdz(2, 2) - dxdz(2, 0) * dxdz(1, 2)) +
739 dxdz(0, 2) * (dxdz(1, 0) * dxdz(2, 1) - dxdz(2, 0) * dxdz(1, 1));
740 mn = min(mn, jacDet);
741 mx = max(mx, jacDet);
753 void BLMesh::BuildElements()
756 map<CADOrientation::Orientation, vector<int> > baseTri;
757 map<CADOrientation::Orientation, vector<int> > topTri;
787 for (
int i = 0; i < m_mesh->m_element[2].size(); i++)
791 find(m_blsurfs.begin(), m_blsurfs.end(), el->m_parentCAD->GetId());
793 if (f == m_blsurfs.end())
799 vector<NodeSharedPtr> tn(3);
800 vector<NodeSharedPtr> pn(6);
801 vector<NodeSharedPtr> n = el->GetVertexList();
804 for (
int j = 0; j < 3; j++)
806 pn[baseTri[o][j]] = n[j];
807 pn[topTri[o][j]] = m_blData[n[j]]->pNode;
808 tn[j] = m_blData[n[j]]->pNode;
812 tags.push_back(m_id);
817 m_mesh->m_element[3].push_back(E);
822 m_psuedoSurface.push_back(T);
824 T->m_parentCAD = el->m_parentCAD;
830 NekDouble BLMesh::Visability(vector<ElementSharedPtr> tris,
833 NekDouble mn = numeric_limits<double>::max();
835 for (
int i = 0; i < tris.size(); i++)
838 NekDouble dt = tmp[0] * N[0] + tmp[1] * N[1] + tmp[2] * N[2];
847 vector<Array<OneD, NekDouble> > N;
848 for (
int i = 0; i < tris.size(); i++)
850 N.push_back(tris[i]->Normal(
true));
853 vector<NekDouble> w(N.size());
856 for (
int i = 0; i < N.size(); i++)
858 w[i] = 1.0 / N.size();
861 for (
int i = 0; i < N.size(); i++)
863 Np[0] += w[i] * N[i][0];
864 Np[1] += w[i] * N[i][1];
865 Np[2] += w[i] * N[i][2];
867 NekDouble mag = sqrt(Np[0] * Np[0] + Np[1] * Np[1] + Np[2] * Np[2]);
876 vector<NekDouble> a(N.size());
877 while (fabs(dot - 1) > 1e-6)
886 for (
int i = 0; i < N.size(); i++)
889 Np[0] * N[i][0] + Np[1] * N[i][1] + Np[2] * N[i][2];
890 if (fabs(dot2 - 1) < 1e-9)
892 a[i] = dot2 / fabs(dot2) * 1e-9;
903 for (
int i = 0; i < N.size(); i++)
905 w[i] = w[i] * a[i] / aSum;
909 for (
int i = 0; i < N.size(); i++)
915 for (
int i = 0; i < N.size(); i++)
917 NpN[0] += w[i] * N[i][0];
918 NpN[1] += w[i] * N[i][1];
919 NpN[2] += w[i] * N[i][2];
921 mag = sqrt(NpN[0] * NpN[0] + NpN[1] * NpN[1] + NpN[2] * NpN[2]);
926 Np[0] = 0.8 * NpN[0] + (1.0 - 0.8) * Np[0];
927 Np[1] = 0.8 * NpN[1] + (1.0 - 0.8) * Np[1];
928 Np[2] = 0.8 * NpN[2] + (1.0 - 0.8) * Np[2];
929 mag = sqrt(Np[0] * Np[0] + Np[1] * Np[1] + Np[2] * Np[2]);
934 dot = Np[0] * Nplast[0] + Np[1] * Nplast[1] + Np[2] * Nplast[2];
938 cout <<
"run out of iterations" << endl;
949 NekDouble a = 2.0 * (1.0 - m_prog) / (1.0 - pow(m_prog, m_layer + 1));
951 m_layerT[0] = a * m_prog * m_bl;
952 for (
int i = 1; i < m_layer; i++)
954 m_layerT[i] = m_layerT[i - 1] + a * pow(m_prog, i) * m_bl;
957 if(m_mesh->m_verbose)
959 cout <<
"First layer height " << m_layerT[0] << endl;
971 for (it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end();
974 vector<pair<int, CADSurfSharedPtr> > ss = (*it)->GetCADSurfs();
975 vector<unsigned int> surfs;
976 for (
int i = 0; i < ss.size(); i++)
978 surfs.push_back(ss[i].first);
980 sort(surfs.begin(), surfs.end());
981 vector<unsigned int> inter, diff;
983 set_intersection(m_blsurfs.begin(), m_blsurfs.end(), surfs.begin(),
984 surfs.end(), back_inserter(inter));
985 set_symmetric_difference(inter.begin(), inter.end(), surfs.begin(),
986 surfs.end(), back_inserter(diff));
989 if (inter.size() > 0)
994 bln->stopped =
false;
1003 ASSERTL0(diff.size() <= 1,
"not setup for curve bl refinement");
1004 symSurfs.insert(diff[0]);
1005 bln->symsurf = diff[0];
1013 m_blData[(*it)] = bln;
1020 for (
int i = 0; i < m_mesh->m_element[2].size(); i++)
1023 find(m_blsurfs.begin(), m_blsurfs.end(),
1024 m_mesh->m_element[2][i]->m_parentCAD->GetId());
1026 if (f == m_blsurfs.end())
1032 vector<NodeSharedPtr> ns = m_mesh->m_element[2][i]->GetVertexList();
1033 for (
int j = 0; j < ns.size(); j++)
1035 m_blData[ns[j]]->els.push_back(m_mesh->m_element[2][i]);
1036 m_blData[ns[j]]->surfs.insert(
1037 m_mesh->m_element[2][i]->m_parentCAD->GetId());
1042 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1045 bit->second->N = GetNormal(bit->second->els);
1047 if (Visability(bit->second->els, bit->second->N) < 0.0)
1049 cerr <<
"failed " << bit->first->m_x <<
" " << bit->first->m_y
1050 <<
" " << bit->first->m_z <<
" "
1051 << Visability(bit->second->els, bit->second->N) << endl;
1056 for (
int k = 0; k < 3; k++)
1058 loc[k] += bit->second->N[k] * m_layerT[0];
1061 bit->second->pNode = boost::shared_ptr<Node>(
1062 new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
1063 bit->second->bl = 0;
1066 m_symSurfs = vector<unsigned int>(symSurfs.begin(), symSurfs.end());
1070 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1072 if (!bit->second->onSym)
1079 m_mesh->m_cad->GetSurf(bit->second->symsurf)->ProjectTo(loc, uv);
1082 m_mesh->m_cad->GetSurf(bit->second->symsurf)->P(uv);
1085 N[0] = nl[0] - bit->first->m_x;
1086 N[1] = nl[1] - bit->first->m_y;
1087 N[2] = nl[2] - bit->first->m_z;
1089 NekDouble mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1095 bit->second->AlignNode(m_layerT[0]);
1100 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1103 added.insert(bit->first->m_id);
1104 for (
int i = 0; i < bit->second->els.size(); i++)
1106 vector<NodeSharedPtr> ns = bit->second->els[i]->GetVertexList();
1107 for (
int j = 0; j < ns.size(); j++)
1110 if (t == added.end())
1112 m_nToNInfo[bit->first].push_back(m_blData[ns[j]]);
1118 for (
int l = 0; l < 10; l++)
1120 for (bit = m_blData.begin(); bit != m_blData.end(); bit++)
1122 if (bit->first->GetNumCADSurf() > 1)
1128 vector<blInfoSharedPtr> data = m_nToNInfo[bit->first];
1130 for (
int i = 0; i < data.size(); i++)
1132 NekDouble d = bit->first->Distance(data[i]->oNode);
1134 sumV[0] += data[i]->N[0] / d;
1135 sumV[1] += data[i]->N[1] / d;
1136 sumV[2] += data[i]->N[2] / d;
1142 sqrt(sumV[0] * sumV[0] + sumV[1] * sumV[1] + sumV[2] * sumV[2]);
1149 N[0] = (1.0 - 0.8) * bit->second->N[0] + 0.8 * sumV[0];
1150 N[1] = (1.0 - 0.8) * bit->second->N[1] + 0.8 * sumV[1];
1151 N[2] = (1.0 - 0.8) * bit->second->N[2] + 0.8 * sumV[2];
1153 mag = sqrt(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1159 bit->second->AlignNode(m_layerT[0]);
1177 ASSERTL0(failed == 0,
"some normals failed to generate");
1188 VandermondeI.Invert();
static NekDouble ang(NekDouble x, NekDouble y)
#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.
#define sign(a, b)
return the sign(b)*a
ElementFactory & GetElementFactory()
SharedMatrix GetVandermondeForDeriv(int dir)
Return the Vandermonde matrix of the derivative of the basis functions for the nodal distribution...
bg::model::box< point > box
bool Infont(NodeSharedPtr n, ElementSharedPtr el)
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
box GetBox(NodeSharedPtr n, NekDouble ov)
Specialisation of the NodalUtil class to support nodal prismatic elements.
bg::model::point< double, 3, bg::cs::cartesian > point
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
boost::shared_ptr< CADSurf > CADSurfSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::pair< box, unsigned int > boxI
SharedMatrix GetVandermonde()
Return the Vandermonde matrix for the nodal distribution.
boost::shared_ptr< blInfo > blInfoSharedPtr
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
boost::shared_ptr< Element > ElementSharedPtr
3D electrostatically spaced points on a Prism