001/** 002 * CSTCurveUtils -- A collection of methods for working with CSTCurve objects. 003 * 004 * Copyright (C) 2015-2019, by Joseph A. Huwaldt. All rights reserved. 005 * 006 * This library is free software; you can redistribute it and/or modify it under the terms 007 * of the GNU Lesser General Public License as published by the Free Software Foundation; 008 * either version 2.1 of the License, or (at your option) any later version. 009 * 010 * This library is distributed in the hope that it will be useful, but WITHOUT ANY 011 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 012 * PARTICULAR PURPOSE. See the GNU Library General Public License for more details. 013 * 014 * You should have received a copy of the GNU Lesser General Public License along with 015 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - 016 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html 017 */ 018package geomss.geom.cst; 019 020import geomss.geom.*; 021import geomss.geom.nurbs.BasicNurbsCurve; 022import geomss.geom.nurbs.CurveFactory; 023import geomss.geom.nurbs.KnotVector; 024import jahuwaldt.js.param.Parameter; 025import jahuwaldt.tools.math.BasisFunction; 026import jahuwaldt.tools.math.MathTools; 027import jahuwaldt.tools.math.ModelData; 028import jahuwaldt.tools.math.RootException; 029import java.text.MessageFormat; 030import java.util.Arrays; 031import static java.util.Objects.requireNonNull; 032import java.util.ResourceBundle; 033import javax.measure.quantity.Dimensionless; 034import javax.measure.quantity.Length; 035import javolution.context.ArrayFactory; 036import javolution.context.StackContext; 037 038/** 039 * A collection of methods for working with Class-Shape-Transform or {@link CSTCurve} 040 * curves. 041 * 042 * <p> Modified by: Joseph A. Huwaldt </p> 043 * 044 * @author Joseph A. Huwaldt, Date: October 9, 2015 045 * @version December 23, 2019 046 */ 047public final class CSTCurveUtils { 048 /** 049 * The resource bundle for this package. 050 */ 051 private static final ResourceBundle RESOURCES = CSTCurve.RESOURCES; 052 053 /** 054 * Method that returns a {@link BasicCSTCurve} that approximates the input list of 055 * points on either the upper or lower surface of a blunt nosed airfoil. This is 056 * useful to approximate a specific airfoil or blunt LE body shape. This method offers 057 * control over the leading edge radius and the trailing edge slope. The input points 058 * must be ordered from leading edge to trailing edge and must be co-planar. In 059 * addition, the input xhat (chord- wise direction vector) must be parallel to or 060 * contained in the plane containing the points. The origin of the resulting CSTCurve 061 * will be taken as the first point in the input list of points and the trailing edge 062 * offset will be taken as the distance from the last point in the input list to the 063 * line formed by passing xhat through the origin point. Finally, the points must 064 * represent a shape that can be modeled by a blunt-nosed CST curve with the given 065 * class function or the fit will be very poor. 066 * 067 * @param order The order for the Bezier curve used for the shape function. This 068 * controls the accuracy that the input points are approximated. 069 * @param xhat The chord-wise direction to assume for the input points. This vector 070 * must lie parallel to the plane containing the points being 071 * approximated. May not be null. 072 * @param rLE The LE radius to use for the airfoil. May not be null. 073 * @param tanTE The trailing edge tangent vector on the airfoil curve. May not be null. 074 * @param points The list of airfoil points ordered from leading edge to trailing edge 075 * to be approximated by a CSTCurve. May not be null. 076 * @param tol The tolerance to use when determining if the input points are 077 * co-planar. May not be null. 078 * @return A new BasicCSTCurve that approximates the input list of points. The user 079 * data with the key "CHISQ" contains the chi-squared measure of the fit. A 080 * value near zero indicates a perfect fit. 081 * @throws RootException if a fit to the points could not be made. 082 */ 083 public static BasicCSTCurve approxAirfoilPoints(int order, GeomVector<Dimensionless> xhat, 084 Parameter<Length> rLE, GeomVector<Dimensionless> tanTE, PointString<?> points, 085 Parameter<Length> tol) throws RootException { 086 requireNonNull(rLE); 087 requireNonNull(tanTE); 088 089 StackContext.enter(); 090 try { 091 // Is the input curve planar? 092 BasicNurbsCurve ncrv = CurveFactory.fitPoints(2, requireNonNull(points)); 093 if (!ncrv.isPlanar(requireNonNull(tol))) 094 throw new IllegalArgumentException(RESOURCES.getString("approxPntsNotCoplanarErr")); 095 096 // Find the largest dimension among the inputs. 097 int numDims = points.getPhyDimension(); 098 if (numDims < xhat.getPhyDimension()) 099 numDims = xhat.getPhyDimension(); 100 if (numDims < 3) 101 numDims = 3; 102 103 // Convert to common dimensions. 104 points = points.toDimension(numDims); 105 xhat = xhat.toDimension(numDims); 106 ncrv = ncrv.toDimension(numDims); 107 108 // Since curve is planar, any binormal will work as the plane normal vector. 109 Vector<Dimensionless> nhat = ncrv.getBinormal(0.5); 110 Point pMid = nhat.getOrigin(); 111 112 // Determine the yhat direction for the curve. 113 Vector<Dimensionless> yhat = xhat.cross(nhat).toUnitVector(); 114 115 // The origin is the first point in the input list. 116 Point origin = points.get(0).immutable(); 117 118 // Determine the distance from the end of the points curve to the chord-line. 119 Point pTE = GeomUtil.pointLineClosest(points.get(-1), origin, xhat); 120 Parameter<Length> chord = pTE.distance(origin); 121 Parameter<Length> offsetTE = GeomUtil.pointPlaneDistance(points.get(-1), pTE, yhat); 122 double offsetTE_c = offsetTE.divide(chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT); 123 if (MathTools.isApproxZero(offsetTE_c)) 124 offsetTE_c = 0; 125 126 // Check to make sure that xhat is in the plane of the curve. 127 double ftol = tol.divide(chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT); 128 IntersectType itype = GeomUtil.linePlaneIntersect(origin, xhat, Plane.valueOf(nhat, pMid), null, ftol); 129 if (!itype.equals(IntersectType.COINCIDENT)) 130 throw new IllegalArgumentException(RESOURCES.getString("xhatNotCoincidentWithPntsPlane")); 131 132 // Create the class function we need. 133 CSTClassFunction Cf = CSTClassFunction.BLUNTNOSED_AIRFOIL; 134 135 // Create the shape function values we are trying to approximate. 136 int size = points.size(); 137 double[] Sf = new double[size]; 138 double[] x_cArr = new double[size]; 139 int sizem1 = size - 1; 140 for (int i = 1; i < sizem1; ++i) { 141 Vector<Length> pv = points.get(i).toGeomVector(); 142 double x_c = xhat.times(pv).divide(chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT); 143 double y_c = yhat.times(pv).divide(chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT); 144 if (x_c < 0.) 145 throw new IllegalArgumentException("X-location of point #" + (i+1) + " is < 0. All X-coordinates must be positive."); 146 x_cArr[i] = x_c; 147 Sf[i] = (y_c - offsetTE_c * x_c) / Cf.getValue(x_c); 148 } 149 150 // Use the input leading edge radius to compute Sf[0] = sqrt(2*rLE/c). 151 x_cArr[0] = 0; 152 double rLE_c = rLE.divide(chord).asType(Dimensionless.class).getValue(Dimensionless.UNIT); 153 Sf[0] = Math.sqrt(2*rLE_c); 154 155 // Use the input trailing edge tangent vector to compute Sf[1] = tan(beta) + offsetTE/c. 156 x_cArr[sizem1] = 1; 157 Parameter tanTEx = xhat.times(tanTE); 158 Parameter tanTEy = yhat.times(tanTE).opposite(); 159 double tanBeta = tanTEy.divide(tanTEx).getValue(); 160 Sf[sizem1] = tanBeta + offsetTE_c; 161 162 // Fit a set of basis function coefficients to the Sf data in a least-square error sense. 163 double[] coefs = new double[order]; 164 double[] sig = new double[size]; 165 Arrays.fill(sig, 1.0); 166 167 // Force the curve to match the input LE radius and TE slope more accurately than other points. 168 sig[0] = 0.5; 169 sig[sizem1] = 0.5; 170 171 double chisq = ModelData.fit(x_cArr, Sf, sig, new BezierFitBasisFunction(order), coefs); 172 173 // Create the Shape Function from the basis function fit coefficients. 174 CSTShapeFunction Sfunc = CSTShapeFunction.newInstance(order, coefs); 175 176 // Create the CST curve for the airfoil. 177 BasicCSTCurve crv = BasicCSTCurve.newInstance(Cf, Sfunc, offsetTE, origin, chord, xhat, yhat); 178 crv.putUserData("CHISQ", chisq); 179 180 return StackContext.outerCopy(crv); 181 } finally { 182 StackContext.exit(); 183 } 184 } 185 186 /** 187 * BasisFunction used to fit a Bezier curve to a set of data where the parameters of 188 * the function are the coefficients of the Bernstein polynomial of the Bezier curve. 189 */ 190 private static class BezierFitBasisFunction implements BasisFunction { 191 private final KnotVector kv; 192 193 public BezierFitBasisFunction(int order) { 194 kv = KnotVector.bezierKnotVector(order-1); 195 } 196 197 @Override 198 public int getMinNumCoef() { 199 return kv.getDegree()+1; 200 } 201 202 @Override 203 public void func(double x, double[] p) { 204 double[] Nik = kv.basisFunctions(x); 205 int size = p.length; 206 for (int i=0; i < size; ++i) 207 p[i] = Nik[i]; 208 ArrayFactory.DOUBLES_FACTORY.recycle(Nik); 209 } 210 } 211 212 /** 213 * Return a {@link BasicCSTCurve} that represents the thickness distribution for an 214 * airfoil formed by the input upper and lower airfoil curves. The input curves must 215 * be of the same order. 216 * 217 * @param upper The {@link CSTCurve} that represents the upper surface of the airfoil. 218 * May not be null. 219 * @param lower The {@link CSTCurve} that represents the lower surface of the airfoil. 220 * May not be null. 221 * @return A curve that represents the thickness distribution for the airfoil. The 222 * returned curve will have the same units, dimensions, origin, chord and 223 * trailing edge thickness as the input upper surface curve. 224 */ 225 public static BasicCSTCurve getThicknessDist(CSTCurve upper, CSTCurve lower) { 226 // Reference: Kulfan, B.M., Bussoletti, J.E., "'Fundamental' Parametric Geometry 227 // Representations for Aircraft Component Shapes", AIAA-2006-6948, pg. 32. 228 229 CSTShapeFunction Supper = upper.getSFunction(); 230 CSTShapeFunction Slower = lower.getSFunction(); 231 int order = Supper.getOrder(); 232 if (order != Slower.getOrder()) 233 throw new IllegalArgumentException(MessageFormat.format( 234 RESOURCES.getString("crvsDifferentOrder"), "upper", "lower airfoil")); 235 236 // Get the Bernstein Polynomial coefficients for each curve. 237 double[] Aupper = Supper.getCoefficients(); 238 double[] Alower = Slower.getCoefficients(); 239 240 // Thickness distribution is formed form the average of the 241 // upper and lower surface coefficients. 242 double[] Ath = Aupper; 243 for (int i=0; i < order; ++i) { 244 Ath[i] = (Aupper[i] + Alower[i])/2; 245 } 246 247 // Create the Shape Function from the thickness profile coefficients. 248 CSTShapeFunction Sfunc = CSTShapeFunction.newInstance(order, Ath); 249 250 // Create the CST curve for the airfoil. 251 BasicCSTCurve crv = BasicCSTCurve.newInstance(upper.getCFunction(), Sfunc, upper.getOffsetTE(), 252 upper.getOrigin(), upper.getChord(), upper.getXhat(), upper.getYhat()); 253 254 return crv; 255 } 256 257 /** 258 * Return a {@link BasicCSTCurve} that represents the camber distribution for an 259 * airfoil formed by the input upper and lower airfoil curves. The input curves must 260 * be of the same order. 261 * 262 * @param upper The {@link CSTCurve} that represents the upper surface of the airfoil. 263 * May not be null. 264 * @param lower The {@link CSTCurve} that represents the lower surface of the airfoil. 265 * May not be null. 266 * @return A curve that represents the camber distribution for the airfoil. The 267 * returned curve will have the same units, dimensions, origin, chord and 268 * trailing edge thickness as the input upper surface curve. 269 */ 270 public static BasicCSTCurve getCamberDist(CSTCurve upper, CSTCurve lower) { 271 // Reference: Kulfan, B.M., Bussoletti, J.E., "'Fundamental' Parametric Geometry 272 // Representations for Aircraft Component Shapes", AIAA-2006-6948, pg. 32. 273 274 CSTShapeFunction Supper = upper.getSFunction(); 275 CSTShapeFunction Slower = lower.getSFunction(); 276 int order = Supper.getOrder(); 277 if (order != Slower.getOrder()) 278 throw new IllegalArgumentException(MessageFormat.format( 279 RESOURCES.getString("crvsDifferentOrder"), "upper", "lower airfoil")); 280 281 // Get the Bernstein Polynomial coefficients for each curve. 282 double[] Aupper = Supper.getCoefficients(); 283 double[] Alower = Slower.getCoefficients(); 284 285 // Camber distribution is formed form the difference of the 286 // upper and lower surface coefficients. 287 double[] Acamb = Aupper; 288 for (int i=0; i < order; ++i) { 289 Acamb[i] = (Aupper[i] - Alower[i])/2; 290 } 291 292 // Create the Shape Function from the thickness profile coefficients. 293 CSTShapeFunction Sfunc = CSTShapeFunction.newInstance(order, Acamb); 294 295 // Create the CST curve for the airfoil. 296 BasicCSTCurve crv = BasicCSTCurve.newInstance(upper.getCFunction(), Sfunc, upper.getOffsetTE(), 297 upper.getOrigin(), upper.getChord(), upper.getXhat(), upper.getYhat()); 298 299 return crv; 300 } 301}