Branch data Line data Source code
1 : : #include "meshkit/AF2PlaneProjection.hpp"
2 : :
3 : : // C++
4 : : #include <cmath>
5 : :
6 : : // CGM
7 : : #include "CubitDefines.h"
8 : : #include "CubitVector.hpp"
9 : : #include "RefFace.hpp"
10 : : #include "Surface.hpp"
11 : :
12 : : // MeshKit
13 : : #include "meshkit/Error.hpp"
14 : : #include "meshkit/iGeom.hpp"
15 : : #include "meshkit/VecUtil.hpp"
16 : :
17 : 10056 : AF2PlaneProjection::AF2PlaneProjection(iGeom* iGeomPtrArg,
18 : : iGeom::EntityHandle srfcHandle,
19 : : MeshKit::Vector<3> const & planeOrigin,
20 : : MeshKit::Vector<3> const & planeNormal,
21 : : MeshKit::Vector<3> const & planeXDir,
22 : : double scaleFactor) :
23 : : iGeomPtr(iGeomPtrArg),
24 : : surface(srfcHandle),
25 : : pOrigin(planeOrigin),
26 : : pNormal(planeNormal),
27 [ + - ]: 10056 : pXDir(planeXDir)
28 : : {
29 : : // exercise the iGeom pointer and check that the entity is a face
30 : : iBase_EntityType eTypeResult;
31 [ + - ][ + - ]: 10056 : MeshKit::IBERRCHK(iGeomPtr->getEntType(surface, eTypeResult), *iGeomPtr);
32 [ - + ]: 10056 : if (eTypeResult != iBase_FACE)
33 : : {
34 [ # # ]: 0 : throw MeshKit::Error(MeshKit::MK_BAD_INPUT, "Geometry entity handle does not refer to a surface in AF2PlaneProjection constructor.");
35 : : }
36 : :
37 : : // compute the squared length of the normal vector
38 [ + - ][ + - ]: 10056 : double nvNormSqrd = pNormal[0] * pNormal[0] +
39 [ + - ][ + - ]: 10056 : pNormal[1] * pNormal[1] + pNormal[2] * pNormal[2];
[ + - ][ + - ]
40 : :
41 : : // check that the length is sufficiently non-zero
42 [ - + ]: 10056 : if (std::fabs(nvNormSqrd) < 1.0e-16)
43 : : {
44 [ # # ]: 0 : throw MeshKit::Error(MeshKit::MK_BAD_INPUT, "Normal vector length is near zero in AF2PlaneProjection constructor.");
45 : : }
46 : :
47 : : // normalize the normal vector to have length one, if necessary
48 [ - + ]: 10056 : if (std::fabs(nvNormSqrd - 1.0) > 1.0e-8)
49 : : {
50 [ # # ]: 0 : pNormal *= 1.0 / std::sqrt(nvNormSqrd);
51 : : }
52 : :
53 : : // compute the squared length of the x-direction vector
54 [ + - ][ + - ]: 10056 : double xdvNormSqrd = pXDir[0] * pXDir[0] +
55 [ + - ][ + - ]: 10056 : pXDir[1] * pXDir[1] + pXDir[2] * pXDir[2];
[ + - ][ + - ]
56 : :
57 : : // check that the length is sufficiently non-zero
58 [ - + ]: 10056 : if (std::fabs(xdvNormSqrd) < 1.0e-16)
59 : : {
60 [ # # ]: 0 : throw MeshKit::Error(MeshKit::MK_BAD_INPUT, "X-direction vector length is near zero in AF2PlaneProjection constructor.");
61 : : }
62 : :
63 : : // normalize the x-direction vector to have length one, if necessary
64 [ + + ]: 10056 : if (std::fabs(xdvNormSqrd - 1.0) > 1.0e-8)
65 : : {
66 [ + - ]: 2 : pXDir *= 1.0 / std::sqrt(xdvNormSqrd);
67 : : }
68 : :
69 [ + - ]: 10056 : double nxDot = MeshKit::VecUtil::dot((double*)pNormal.data(),
70 [ + - ][ + - ]: 20112 : (double*)pXDir.data());
71 [ - + ]: 10056 : if (std::fabs(nxDot) > 1.0e-12)
72 : : {
73 [ # # ]: 0 : throw MeshKit::Error(MeshKit::MK_BAD_INPUT, "X-direction vector is not normal to the normal vector in AF2PlaneProjection constructor.");
74 : : }
75 : :
76 : : // compute the y-direction vector
77 [ + - ][ + - ]: 10056 : MeshKit::VecUtil::cross((double*)pNormal.data(), (double*)pXDir.data(),
78 [ + - ][ + - ]: 20112 : (double*)pYDir.data());
79 : :
80 [ - + ]: 10056 : if (scaleFactor <= 0.0)
81 : : {
82 : : throw MeshKit::Error(MeshKit::MK_BAD_INPUT,
83 [ # # ]: 0 : "Scale factor is not positive.");
84 : : }
85 : 10056 : scale = 1.0 / scaleFactor;
86 : 10056 : }
87 : :
88 : 288951 : AF2Point2D* AF2PlaneProjection::transformFromSurface(
89 : : AF2Point3D const & srfcPnt, bool & legal) const
90 : : {
91 [ + - ][ + - ]: 288951 : double diffX = srfcPnt.getX() - pOrigin[0];
92 [ + - ][ + - ]: 288951 : double diffY = srfcPnt.getY() - pOrigin[1];
93 [ + - ][ + - ]: 288951 : double diffZ = srfcPnt.getZ() - pOrigin[2];
94 : :
95 [ + - ]: 288951 : double diffDotNormal = diffX * pNormal[0] +
96 [ + - ][ + - ]: 288951 : diffY * pNormal[1] + diffZ * pNormal[2];
97 [ + - ]: 288951 : double normalCmpX = diffDotNormal * pNormal[0];
98 [ + - ]: 288951 : double normalCmpY = diffDotNormal * pNormal[1];
99 [ + - ]: 288951 : double normalCmpZ = diffDotNormal * pNormal[2];
100 : :
101 : 288951 : double planarCmpX = diffX - normalCmpX;
102 : 288951 : double planarCmpY = diffY - normalCmpY;
103 : 288951 : double planarCmpZ = diffZ - normalCmpZ;
104 : :
105 [ + - ]: 288951 : double planePntX = scale * (planarCmpX * pXDir[0] +
106 [ + - ][ + - ]: 288951 : planarCmpY * pXDir[1] + planarCmpZ * pXDir[2]);
107 [ + - ]: 288951 : double planePntY = scale * (planarCmpX * pYDir[0] +
108 [ + - ][ + - ]: 288951 : planarCmpY * pYDir[1] + planarCmpZ * pYDir[2]);
109 : :
110 : : // compute normal at the surface point; if too different orientation, we should not
111 : : // consider the point legal
112 : : // maybe we should keep the normal to the points on the front available always, we repeat this too much
113 : : double normalX, normalY, normalZ;
114 : : iGeomPtr->getEntNrmlXYZ(surface, srfcPnt.getX(), srfcPnt.getY(),
115 [ + - ][ + - ]: 288951 : srfcPnt.getZ(), normalX, normalY, normalZ);
[ + - ][ + - ]
116 : : // normalize?
117 : 288951 : double sqrN= sqrt(normalX*normalX + normalY*normalY + normalZ*normalZ);
118 [ - + ]: 288951 : if (sqrN<=1.e-12)
119 : 0 : legal= false;
120 : : else
121 : : {
122 : 288951 : normalX/=sqrN; normalY/=sqrN; normalZ/=sqrN;
123 [ + - ]: 288951 : double nrmlDot = normalX * pNormal[0] +
124 [ + - ][ + - ]: 288951 : normalY * pNormal[1] + normalZ * pNormal[2];
125 : :
126 [ + + ]: 288951 : if (nrmlDot < 0.0009765625) // 1.0/1024 == 0.0009765625)
127 : 6986 : legal= false;
128 : : }
129 : :
130 [ + - ][ + - ]: 288951 : return new AF2Point2D(planePntX, planePntY);
131 : : }
132 : :
133 : 4056 : AF2Point3D* AF2PlaneProjection::transformToSurface(
134 : : AF2Point2D const & planePnt, unsigned long const & pntId) const
135 : : {
136 : 4056 : MeshKit::Vector<3> rayOrigin(pOrigin);
137 [ + - ][ + - ]: 4056 : rayOrigin += (planePnt.getX() / scale) * pXDir +
[ + - ]
138 [ + - ][ + - ]: 8112 : (planePnt.getY() / scale) * pYDir;
[ + - ]
139 : : /* // we are looking for the point on the surface, within some distance away
140 : : // from our plane;
141 : : const double MYTOL = 1.e-7;
142 : : // maxDistanceFromPlane is a good measure of how far is the surface from the projection plane
143 : : double alfa = 1000*MYTOL + 3.0*maxDistanceFromPlane;
144 : : if (alfa >2.0)
145 : : alfa =2.0; // do not go overboard
146 : : if (alfa < 0.5)
147 : : alfa = 0.5;*/
148 : :
149 [ + - ][ + - ]: 4056 : CubitVector cvRayOrigin(rayOrigin[0], rayOrigin[1], rayOrigin[2]);
[ + - ][ + - ]
150 [ + - ][ + - ]: 4056 : CubitVector cvRayDir(pNormal[0], pNormal[1], pNormal[2]);
[ + - ][ + - ]
151 [ + - ]: 4056 : CubitVector cvPointOnSrfc;
152 [ + - ]: 4056 : CubitVector closestPointOnSrfc;
153 : :
154 : : // closest point in positive normal direction
155 : : RefFace* cgmFacePtr =
156 [ - + ]: 4056 : dynamic_cast<RefFace*>(reinterpret_cast<RefEntity*>(surface));
157 [ + - ]: 4056 : Surface* cgmSrfcPtr = cgmFacePtr->get_surface_ptr();
158 : : CubitStatus posDirResult = cgmSrfcPtr->closest_point_along_vector(
159 [ + - ]: 4056 : cvRayOrigin, cvRayDir, cvPointOnSrfc);
160 [ + + ]: 4056 : if (posDirResult == CUBIT_SUCCESS)
161 : : {
162 [ + - ]: 2213 : closestPointOnSrfc = cvPointOnSrfc;
163 : : }
164 : :
165 [ + - ]: 4056 : cvRayDir *= -1;
166 : : CubitStatus negDirResult = cgmSrfcPtr->closest_point_along_vector(
167 [ + - ]: 4056 : cvRayOrigin, cvRayDir, cvPointOnSrfc);
168 [ + + ]: 4056 : if (negDirResult == CUBIT_SUCCESS)
169 : : {
170 [ + + ]: 3903 : if (posDirResult == CUBIT_SUCCESS)
171 : : {
172 : : // REMARK: It might be better to measure the distance in
173 : : // parametric coordinates from some source point . . .
174 : :
175 : : // compute square distance to point in positive direction
176 [ + - ][ + - ]: 2060 : double posXDiff = closestPointOnSrfc.x() - rayOrigin[0];
177 [ + - ][ + - ]: 2060 : double posYDiff = closestPointOnSrfc.y() - rayOrigin[1];
178 [ + - ][ + - ]: 2060 : double posZDiff = closestPointOnSrfc.z() - rayOrigin[2];
179 : 4120 : double posSqDist = posXDiff * posXDiff +
180 : 4120 : posYDiff * posYDiff + posZDiff * posZDiff;
181 : :
182 : : // compute square distance to point in negative direction
183 [ + - ][ + - ]: 2060 : double negXDiff = cvPointOnSrfc.x() - rayOrigin[0];
184 [ + - ][ + - ]: 2060 : double negYDiff = cvPointOnSrfc.y() - rayOrigin[1];
185 [ + - ][ + - ]: 2060 : double negZDiff = cvPointOnSrfc.z() - rayOrigin[2];
186 : 4120 : double negSqDist = negXDiff * negXDiff +
187 : 4120 : negYDiff * negYDiff + negZDiff * negZDiff;
188 : :
189 [ + + ]: 2060 : if (negSqDist < posSqDist)
190 : : {
191 [ + - ]: 2060 : closestPointOnSrfc = cvPointOnSrfc;
192 : : }
193 : : }
194 : : else
195 : : {
196 [ + - ]: 3903 : closestPointOnSrfc = cvPointOnSrfc;
197 : : }
198 : : }
199 : :
200 [ + - ]: 4056 : return new AF2Point3D(pntId, closestPointOnSrfc.x(),
201 [ + - ][ + - ]: 4056 : closestPointOnSrfc.y(), closestPointOnSrfc.z());
[ + - ][ + - ]
202 [ + - ][ + - ]: 156 : }
|
