44 namespace NekMeshUtils
47 CADSurf::CADSurf(
int i, TopoDS_Shape in, vector<EdgeLoop> ein) : m_edges(ein)
52 gp_Pnt ori(0.0, 0.0, 0.0);
53 transform.SetScale(ori, 1.0 / 1000.0);
54 TopLoc_Location mv(transform);
55 m_s = BRep_Tool::Surface(TopoDS::Face(in));
67 gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
69 GeomAPI_ProjectPointOnSurf projection(
72 Extrema_ExtAlgo_Tree);
75 if (projection.NbPoints() == 0)
78 ShapeAnalysis_Surface sas(m_s);
83 gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
87 gp_Pnt p3 = sas.Value(p2);
88 ASSERTL0(p3.Distance(loc) < 1e-3,
"large locuv distance sas");
92 Quantity_Parameter ui;
93 Quantity_Parameter vi;
95 projection.Parameters(1, ui, vi);
100 if (projection.Distance(1) > 1.0)
103 cerr <<
"large locuv distance " << projection.Distance(1) / 1000.0
137 ASSERTL0(
false,
"Cannot correct locuv");
152 if (n[0] == 0 && n[1] == 0 && n[2] == 0)
160 NekDouble E = r[3] * r[3] + r[4] * r[4] + r[5] * r[5];
161 NekDouble F = r[3] * r[6] + r[4] * r[7] + r[5] * r[8];
162 NekDouble G = r[6] * r[6] + r[7] * r[7] + r[8] * r[8];
163 NekDouble e = n[0] * r[9] + n[1] * r[10] + n[2] * r[11];
164 NekDouble f = n[0] * r[15] + n[1] * r[16] + n[2] * r[17];
165 NekDouble g = n[0] * r[12] + n[1] * r[13] + n[2] * r[14];
168 if (E * G - F * F < 1E-30)
175 K = (e * g - f * f) / (E * G - F * F);
176 H = 0.5 * (e * G - 2 * f * F + g * E) / (E * G - F * F);
179 kv[0] = abs(H + sqrt(H * H - K));
180 kv[1] = abs(H - sqrt(H * H - K));
182 return kv[0] > kv[1] ? kv[0] : kv[1];
187 gp_Pnt loc(p[0] * 1000.0, p[1] * 1000.0, p[2] * 1000.0);
190 ShapeAnalysis_Surface sas(m_s);
195 gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
197 gp_Pnt p3 = sas.Value(p2);
199 return p3.Distance(loc);
204 gp_Pnt loc(tp[0] * 1000.0, tp[1] * 1000.0, tp[2] * 1000.0);
207 ShapeAnalysis_Surface sas(m_s);
212 gp_Pnt2d p2 = sas.ValueOfUV(loc, 1e-7);
214 gp_Pnt p3 = sas.Value(p2);
216 tp[0] = p3.X() / 1000.0;
217 tp[1] = p3.Y() / 1000.0;
218 tp[2] = p3.Z() / 1000.0;
243 location[0] = loc.X();
244 location[1] = loc.Y();
245 location[2] = loc.Z();
257 gp_Vec n = D1U.Crossed(D1V);
264 if (n.X() == 0 && n.Y() == 0 && n.Z() == 0)
310 gp_Vec D1U, D1V, D2U, D2V, D2UV;
311 m_occSurface.D2(uv[0], uv[1], Loc, D1U, D1V, D2U, D2V, D2UV);
339 error <<
"Point not within parameter plane: ";
345 if (fabs(uv[0] -
m_occSurface.FirstUParameter()) > 1E-8)
347 error <<
"U(" << uv[0] <<
") is less than Umin("
356 error <<
"U(" << uv[0] <<
") is greater than Umax("
363 if (fabs(uv[1] -
m_occSurface.FirstVParameter()) > 1E-8)
365 error <<
"V(" << uv[1] <<
") is less than Vmin("
374 error <<
"V(" << uv[1] <<
") is greater than Vmax("
380 error <<
" On Surface: " <<
GetId();
Array< OneD, NekDouble > D1(Array< OneD, NekDouble > uv)
Get the set of first derivatives at parametric point u,v.
#define ASSERTL0(condition, msg)
Array< OneD, NekDouble > P(Array< OneD, NekDouble > uv)
Get the x,y,z at parametric point u,v.
bool m_hasTwoCurves
flag to alert the mesh generation to a potential problem is both curves have only two points in the m...
Array< OneD, NekDouble > locuv(Array< OneD, NekDouble > p)
Performs a reverse look up to find u,v and x,y,z.
BRepAdaptor_Surface m_occSurface
OpenCascade object for surface.
Array< OneD, NekDouble > N(Array< OneD, NekDouble > uv)
Get the normal vector at parametric point u,v.
NekDouble DistanceTo(Array< OneD, NekDouble > p)
does unconstrained locuv to project point from anywhere
bool m_correctNormal
normal
int GetId()
Return ID of the vertex.
bool IsPlane()
returns true if the surface is flat (2D)
void ProjectTo(Array< OneD, NekDouble > &tp, Array< OneD, NekDouble > &uv)
void Test(Array< OneD, NekDouble > uv)
Function which tests the the value of uv used is within the surface.
cadType m_type
type of the cad object
Array< OneD, NekDouble > D2(Array< OneD, NekDouble > uv)
Get the set of second derivatives at parametric point u,v.
NekDouble Curvature(Array< OneD, NekDouble > uv)
returns curvature at point uv