diff --git a/doc/mesh-grading-function.pdf b/doc/mesh-grading-function.pdf new file mode 100644 index 000000000..78903c59b Binary files /dev/null and b/doc/mesh-grading-function.pdf differ diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index 2b226af9e..07521629f 100644 --- a/src/ASM/ASMbase.C +++ b/src/ASM/ASMbase.C @@ -23,6 +23,7 @@ bool ASMbase::fixHomogeneousDirichlet = true; +int ASMbase::dbgElm = 0; /*! @@ -1008,8 +1009,8 @@ bool ASMbase::extractNodalVec (const Vector& globRes, Vector& nodeVec, bool ASMbase::injectNodalVec (const Vector& nodeVec, Vector& globVec, const std::vector& madof, int basis) const { - if (madof.empty()) - { + if (madof.empty()) + { std::cerr <<" *** ASMbase::injectNodalVec: Empty MADOF array."<< std::endl; return false; } diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 07736746d..abe5ae5a7 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -654,6 +654,8 @@ class ASMbase public: static bool fixHomogeneousDirichlet; //!< If \e true, pre-eliminate fixed DOFs + static int dbgElm; //!< One-based element index to print debugging info for + size_t idx; //!< Index of this patch in the multi-patch model protected: diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 1ac5dfdb7..3bc795201 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -688,7 +688,11 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers, for (i = 0; i < n1; i++, node += i1) { int slave = slaveNodes[revers ? n1-i-1 : i]; - if (coordCheck && !neighbor.getCoord(node).equal(this->getCoord(slave),xtol)) + if (!coordCheck) + ASMbase::collapseNodes(neighbor,node,*this,slave); + else if (neighbor.getCoord(node).equal(this->getCoord(slave),xtol)) + ASMbase::collapseNodes(neighbor,node,*this,slave); + else { std::cerr <<" *** ASMs2D::connectPatch: Non-matching nodes " << node <<": "<< neighbor.getCoord(node) @@ -696,8 +700,6 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers, << slave <<": "<< this->getCoord(slave) << std::endl; return false; } - else - ASMbase::collapseNodes(neighbor,node,*this,slave); } return true; @@ -1017,7 +1019,7 @@ size_t ASMs2D::constrainEdgeLocal (int dir, bool open, int dof, int code, void ASMs2D::constrainCorner (int I, int J, int dof, int code, char basis) { int node = this->getCorner(I, J, basis); - if (node > -1) + if (node > 0) this->prescribe(node,dof,code); } @@ -1535,6 +1537,11 @@ bool ASMs2D::integrate (Integrand& integrand, fe.iel = MLGE[iel]; if (fe.iel < 1) continue; // zero-area element +#ifdef SP_DEBUG + if (dbgElm < 0 && 1+iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + int i1 = p1 + iel % nel1; int i2 = p2 + iel / nel1; @@ -1686,10 +1693,13 @@ bool ASMs2D::integrate (Integrand& integrand, utl::getGmat(Jac,dXidu,fe.G); #if SP_DEBUG > 4 - std::cout <<"\niel, ip = "<< iel <<" "<< ip - <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; - if (!fe.d2NdX2.empty()) - std::cout <<"d2NdX2 ="<< fe.d2NdX2; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + { + std::cout <<"\niel, ip = "<< iel <<" "<< ip + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (!fe.d2NdX2.empty()) + std::cout <<"d2NdX2 ="<< fe.d2NdX2; + } #endif // Cartesian coordinates of current integration point @@ -1714,6 +1724,10 @@ bool ASMs2D::integrate (Integrand& integrand, ok = false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } } } @@ -1788,6 +1802,11 @@ bool ASMs2D::integrate (Integrand& integrand, fe.iel = MLGE[iel]; if (fe.iel < 1) continue; // zero-area element +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + int i1 = p1 + iel % nel1; int i2 = p2 + iel / nel1; @@ -1864,8 +1883,9 @@ bool ASMs2D::integrate (Integrand& integrand, utl::getGmat(Jac,dXidu,fe.G); #if SP_DEBUG > 4 - std::cout <<"\niel, jp = "<< iel <<" "<< jp - <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + std::cout <<"\niel, jp = "<< iel <<" "<< jp + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; #endif // Cartesian coordinates of current integration point @@ -1890,6 +1910,10 @@ bool ASMs2D::integrate (Integrand& integrand, ok = false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } } } @@ -2134,6 +2158,11 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, fe.iel = abs(MLGE[doXelms+iel-1]); if (fe.iel < 1) continue; // zero-area element +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + // Skip elements that are not on current boundary edge bool skipMe = false; switch (edgeDir) @@ -2220,6 +2249,11 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex, A->destruct(); if (!ok) return false; + +#ifdef SP_DEBUG + if (dbgElm < 0 && iel == -dbgElm) + break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -2728,7 +2762,7 @@ bool ASMs2D::evaluate (const RealFunc* func, RealArray& vec, } -int ASMs2D::getCorner(int I, int J, int basis) const +int ASMs2D::getCorner (int I, int J, int basis) const { int n1, n2, node = 1; for (char i = 1; i <= basis; i++) diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index 8534e89c9..c24f89ad0 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -756,7 +756,11 @@ bool ASMs3D::connectBasis (int face, ASMs3D& neighbor, int nface, int norient, } int slave = slaveNodes[k][l]; - if (coordCheck && !neighbor.getCoord(node).equal(this->getCoord(slave),xtol)) + if (!coordCheck) + ASMbase::collapseNodes(neighbor,node,*this,slave); + else if (neighbor.getCoord(node).equal(this->getCoord(slave),xtol)) + ASMbase::collapseNodes(neighbor,node,*this,slave); + else { std::cerr <<" *** ASMs3D::connectPatch: Non-matching nodes " << node <<": "<< neighbor.getCoord(node) @@ -764,8 +768,6 @@ bool ASMs3D::connectBasis (int face, ASMs3D& neighbor, int nface, int norient, << slave <<": "<< this->getCoord(slave) << std::endl; return false; } - else - ASMbase::collapseNodes(neighbor,node,*this,slave); } return true; @@ -1110,7 +1112,7 @@ size_t ASMs3D::constrainFaceLocal (int dir, bool open, int dof, int code, } -std::vector ASMs3D::getEdge(int lEdge, bool open, int basis) const +std::vector ASMs3D::getEdge (int lEdge, bool open, int basis) const { std::vector result; int n1, n2, n3, n, inc = 1; @@ -1254,10 +1256,8 @@ void ASMs3D::constrainCorner (int I, int J, int K, int dof, int code, char basis) { int node = this->getCorner(I, J, K, basis); - if (node < 1) - return; - - this->prescribe(node,dof,code); + if (node > 0) + this->prescribe(node,dof,code); } @@ -1785,6 +1785,11 @@ bool ASMs3D::integrate (Integrand& integrand, fe.iel = MLGE[iel]; if (fe.iel < 1) continue; // zero-volume element +#ifdef SP_DEBUG + if (dbgElm < 0 && 1+iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + int i1 = p1 + iel % nel1; int i2 = p2 + (iel / nel1) % nel2; int i3 = p3 + iel / (nel1*nel2); @@ -1946,8 +1951,13 @@ bool ASMs3D::integrate (Integrand& integrand, utl::getGmat(Jac,dXidu,fe.G); #if SP_DEBUG > 4 - std::cout <<"\niel, ip = "<< iel <<" "<< ip - <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + { + std::cout <<"\niel, ip = "<< iel <<" "<< ip + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (!fe.d2NdX2.empty()) + std::cout <<"d2NdX2 ="<< fe.d2NdX2; + } #endif // Cartesian coordinates of current integration point @@ -1972,6 +1982,10 @@ bool ASMs3D::integrate (Integrand& integrand, ok = false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } } } @@ -2050,6 +2064,11 @@ bool ASMs3D::integrate (Integrand& integrand, fe.iel = MLGE[iel]; if (fe.iel < 1) continue; // zero-volume element +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + int i1 = p1 + iel % nel1; int i2 = p2 + (iel / nel1) % nel2; int i3 = p3 + iel / (nel1*nel2); @@ -2129,8 +2148,9 @@ bool ASMs3D::integrate (Integrand& integrand, utl::getGmat(Jac,dXidu,fe.G); #if SP_DEBUG > 4 - std::cout <<"\niel, jp = "<< iel <<" "<< jp - <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + std::cout <<"\niel, jp = "<< iel <<" "<< jp + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; #endif // Cartesian coordinates of current integration point @@ -2155,6 +2175,10 @@ bool ASMs3D::integrate (Integrand& integrand, ok = false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } } } @@ -2278,6 +2302,11 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex, fe.iel = abs(MLGE[doXelms+iel]); if (fe.iel < 1) continue; // zero-volume element +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + int i1 = p1 + iel % nel1; int i2 = p2 + (iel / nel1) % nel2; int i3 = p3 + iel / (nel1*nel2); @@ -2395,6 +2424,11 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex, ok = false; A->destruct(); + +#ifdef SP_DEBUG + if (dbgElm < 0 && iel == -dbgElm) + break; // Skipping all elements, except for -dbgElm +#endif } } } @@ -3190,7 +3224,7 @@ bool ASMs3D::evaluate (const RealFunc* func, RealArray& vec, } -int ASMs3D::getCorner(int I, int J, int K, int basis) const +int ASMs3D::getCorner (int I, int J, int K, int basis) const { int n1, n2, n3; int node = this->findStartNode(n1,n2,n3,basis); diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index 465c87c08..97365fec1 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -448,15 +448,14 @@ bool ASMs3Dmx::integrate (Integrand& integrand, // Evaluate basis function derivatives at all integration points std::vector> splinex(m_basis.size()); std::vector> splinex2(m_basis.size()); - if (use2ndDer) { + if (use2ndDer) #pragma omp parallel for schedule(static) for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],gpar[2],splinex2[i]); - } else { + else #pragma omp parallel for schedule(static) for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]); - } std::vector elem_sizes; for (auto& it : m_basis) @@ -513,7 +512,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, if (useElmVtx) this->getElementCorners(i1-1,i2-1,i3-1,fe.XC); - + if (integrand.getIntegrandType() & Integrand::G_MATRIX) { // Element size in parametric space @@ -553,13 +552,12 @@ bool ASMs3Dmx::integrate (Integrand& integrand, fe.w = gpar[2](k+1,i3-p3+1); // Fetch basis function derivatives at current integration point - if (use2ndDer) { + if (use2ndDer) for (size_t b = 0; b < m_basis.size(); ++b) SplineUtils::extractBasis(splinex2[b][ip],fe.basis(b+1),dNxdu[b], d2Nxdu2[b]); - } else { + else for (size_t b = 0; b < m_basis.size(); ++b) SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),dNxdu[b]); - } // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 510ad9ee0..d5aebc0bd 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -637,13 +637,13 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis) // get all elements and functions on this edge std::vector edgeFunctions; std::vector edgeElements; - lr->getEdgeFunctions(edgeFunctions, edge); - lr->getEdgeElements( edgeElements, edge); + lr->getEdgeFunctions(edgeFunctions,edge); + lr->getEdgeElements (edgeElements ,edge); // find the corners since these are not to be included in the L2-fitting // of the inhomogenuous dirichlet boundaries; corners are interpolatory. - // Optimization note: loop over the "edge"-container to manually pick up the - // end nodes. LRspine::getEdgeFunctions() does a global search. + // Optimization note: loop over the "edge"-container to manually pick up + // the end nodes. LRspine::getEdgeFunctions() does a global search. std::vector c1, c2; switch (edge) { @@ -676,7 +676,7 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis) int bcode = abs(code); int j = 0; - for (auto b : edgeFunctions ) + for (auto b : edgeFunctions) { de.MLGN[j++] = b->getId(); // skip corners for open boundaries @@ -993,6 +993,11 @@ bool ASMu2D::integrate (Integrand& integrand, std::vector::iterator el = lrspline->elementBegin(); for (int iel = 1; el != lrspline->elementEnd(); el++, iel++) { +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + FiniteElement fe(MNPC[iel-1].size()); fe.iel = MLGE[iel-1]; @@ -1091,13 +1096,18 @@ bool ASMu2D::integrate (Integrand& integrand, lrspline->computeBasis(fe.u,fe.v,spline, iel-1); SplineUtils::extractBasis(spline,fe.N,dNdu); #if SP_DEBUG > 4 - std::cout <<"\nBasis functions at a integration point " - <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1] - <<" left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1]; - for (size_t ii = 0; ii < spline.basisValues.size(); ii++) - std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t' - << spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii]; - std::cout << std::endl; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + { + std::cout <<"\nBasis functions at a integration point " + <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1] + <<" left_idx = "<< spline.left_idx[0] + <<" "<< spline.left_idx[1]; + for (size_t ii = 0; ii < spline.basisValues.size(); ii++) + std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t' + << spline.basisDerivs_u[ii] <<'\t' + << spline.basisDerivs_v[ii]; + std::cout << std::endl; + } #endif } @@ -1111,7 +1121,8 @@ bool ASMu2D::integrate (Integrand& integrand, return false; #if SP_DEBUG > 4 - std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; #endif // Cartesian coordinates of current integration point @@ -1134,6 +1145,10 @@ bool ASMu2D::integrate (Integrand& integrand, return false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -1178,6 +1193,11 @@ bool ASMu2D::integrate (Integrand& integrand, double dA = this->getParametricArea(iel); if (dA < 0.0) return false; // topology error (probably logic error) +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + // Set up control point (nodal) coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; @@ -1230,8 +1250,9 @@ bool ASMu2D::integrate (Integrand& integrand, return false; #if SP_DEBUG > 4 - std::cout <<"\niel, ip = "<< iel <<" "<< ip - <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + std::cout <<"\niel, ip = "<< iel <<" "<< ip + <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; #endif // Cartesian coordinates of current integration point @@ -1254,6 +1275,10 @@ bool ASMu2D::integrate (Integrand& integrand, return false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -1306,6 +1331,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, std::vector::iterator el = lrspline->elementBegin(); for (int iel = 1; el != lrspline->elementEnd(); el++, iel++) { +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + // Skip elements that are not on current boundary edge bool skipMe = false; switch (edgeDir) @@ -1383,6 +1413,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, return false; A->destruct(); + +#ifdef SP_DEBUG + if (dbgElm < 0 && iel == -dbgElm) + break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -1870,9 +1905,9 @@ bool ASMu2D::updateDirichlet (const std::map& func, (*mit)->setSlaveCoeff(edgeControlmatrix[dof-1][nit->first]); else //scalar condition (*mit)->setSlaveCoeff(edgeControlpoints[nit->first]); - #if SP_DEBUG > 1 +#if SP_DEBUG > 1 std::cout <<"Updated constraint: "<< **mit; - #endif +#endif } } } diff --git a/src/SIM/SIMdummy.h b/src/SIM/SIMdummy.h index 0a31e6b19..33d3bb025 100644 --- a/src/SIM/SIMdummy.h +++ b/src/SIM/SIMdummy.h @@ -14,7 +14,13 @@ #ifndef _SIM_DUMMY_H_ #define _SIM_DUMMY_H_ -#include "SIMbase.h" +#include +#include + +class ASMbase; +class IntegrandBase; +class ModelGenerator; +class TiXmlElement; /*! @@ -34,11 +40,16 @@ template class SIMdummy : public Base //! \brief Returns the number of parameter dimensions in the model. virtual unsigned short int getNoParamDim() const { return 0; } //! \brief Reads a patch from given input stream. - virtual ASMbase* readPatch(std::istream&,int, const SIMbase::CharVec&) const + virtual ASMbase* readPatch(std::istream&,int, + const std::vector&) const { return nullptr; } //! \brief Reads patches from given input stream. - virtual bool readPatches(std::istream&,SIMdependency::PatchVec&,const char*) const + virtual bool readPatches(std::istream&, + std::vector&,const char*) const { return false; } + //! \brief Creates the computational FEM model from the spline patches. + virtual bool createFEMmodel(char) { return false; } + protected: //! \brief Preprocesses a user-defined Dirichlet boundary property. virtual bool addConstraint(int,int,int,int,int,int&,char) diff --git a/src/Utility/Test/TestUtilities.C b/src/Utility/Test/TestUtilities.C index c427e3c5b..7ba33b771 100644 --- a/src/Utility/Test/TestUtilities.C +++ b/src/Utility/Test/TestUtilities.C @@ -34,7 +34,7 @@ TEST(TestUtilities, ParseIntegers) TEST(TestUtilities, ParseKnots) { - std::string simple("0 0 0.5 0.9 1.0"); + std::string simple("xi 0 0.5 0.9 1.0"); strtok(const_cast(simple.c_str())," "); std::vector xi; utl::parseKnots(xi); @@ -43,6 +43,17 @@ TEST(TestUtilities, ParseKnots) ASSERT_FLOAT_EQ(xi[1], 0.5); ASSERT_FLOAT_EQ(xi[2], 0.9); ASSERT_FLOAT_EQ(xi[3], 1.0); + + std::string graded5("xi C5 79 0.01 2.0"); + strtok(const_cast(graded5.c_str())," "); + xi.clear(); + utl::parseKnots(xi); + EXPECT_EQ(xi.size(), 79U); + EXPECT_FLOAT_EQ(xi[39], 0.5); + EXPECT_FLOAT_EQ(xi.front()+xi.back(), 1.0); + + for (size_t i = 0; i < xi.size(); i++) + std::cout <<"xi["<< i <<"] = "<< xi[i] << std::endl; } TEST(TestUtilities, ReadLine) diff --git a/src/Utility/Utilities.C b/src/Utility/Utilities.C index 3437a0b78..5cf4cff4f 100644 --- a/src/Utility/Utilities.C +++ b/src/Utility/Utilities.C @@ -32,11 +32,53 @@ void utl::parseIntegers (std::vector& values, const char* argv) } +/*! + If the first character of the first token is a digit, it is assumed that + the knot values are listed explicitly by the remaining tokens. Otherwise, + various mesh grading schemes are assumed with the following syntax: +
    +
  • G n alpha xi1 xi2 - Geometric grading
    • +
  • B n alpha xi1 xi2 - Biased geometric grading
    • +
  • C n alpha xi1 xi2 - Centered geometric grading
    • +
  • C5 n X1 X2 xi1 xi2 - Centred grading based on a 5th order function
    • +
+ In each case the starting and ending knot values, xi1 and xi2, are optional. + If not specified, 0.0 and 1.0 is assumed. + The 4th scheme was proposed by Tymofiy Gerasimov 14.12.2016, and is described + in the file mesh-grading-function.pdf located in the doc folder. +*/ + bool utl::parseKnots (std::vector& xi) { char* cline = strtok(nullptr," "); if (!cline) return false; + else if (toupper(cline[0]) == 'C' && cline[1] == '5') + { + // Centred grading using a 5'th order polynomial + int nX = atoi(strtok(nullptr," ")); + double X1 = (cline = strtok(nullptr," ")) ? atof(cline) : 1.0; + double X2 = (cline = strtok(nullptr," ")) ? atof(cline) : 1.0; + double y0 = (cline = strtok(nullptr," ")) ? atof(cline) : 0.0; + double y1 = (cline = strtok(nullptr," ")) ? atof(cline) : 1.0; + double dy = y1 - y0; + if (nX < 2 || nX%2 == 0 || + X1 <= 0.0 || X1 >= 1.0 || X2 <= 1.0 || + y0 < 0.0 || dy <= 0.0 || y1 > 1.0) + return false; + + double A = ( 16.0*X1 + 8.0*X2 - 24.0)*dy; + double B = (-40.0*X1 - 20.0*X2 + 60.0)*dy; + double C = ( 32.0*X1 + 18.0*X2 - 50.0)*dy; + double E = ( -8.0*X1 - 7.0*X2 + 15.0)*dy; + double F = X2*dy; + dy /= (1+nX); + double y = y0 + dy; + if (y0 > 0.0) xi.push_back(y0); + for (int i = 0; i < nX; i++, y += dy) + xi.push_back(((((A*y+B)*y+C)*y+E)*y+F)*y+y0); + if (y1 < 1.0) xi.push_back(y1); + } else if (isalpha(cline[0])) { // Geometric grading @@ -48,6 +90,7 @@ bool utl::parseKnots (std::vector& xi) double xi2 = (cline = strtok(nullptr," ")) ? atof(cline) : 1.0; if (xi1 < 0.0 || xi2 <= xi1 || xi2 > 1.0 || ru < 1) return false; + if (biased && ru > 1 && alpha != 1.0) alpha = pow(alpha,1.0/double(ru)); else if (centred && ru > 1)